AxialGeoMultiscale Tutorial

1D-3D-partitioned-pipe/Fluid1/Fluid1D.py

@@ -0,0 +1,168 @@
+#!/usr/bin/env python3
+
+from numpy.core.fromnumeric import reshape
+from evtk.hl import pointsToVTK
+import numpy as np
+import precice
+
+def main():
+
+    # number of nodes, length of domain and space interval
+
+    length = 10
+    n = 20
+    h = 1 / n
+    t = 0
+    counter = 0
+
+    # generate mesh
+
+    x = np.zeros(n+1)
+    y = np.zeros(n+1)
+    z = np.linspace(0,length,n+1)
+
+    # initial data
+
+    u = np.zeros((n+1,3))
+    p = np.zeros(n+1)
+    rhs = np.zeros(n+1)
+
+    # preCICE setup
+
+    participant_name = "Fluid1"
+    config_file_name = "../precice-config.xml"
+    solver_process_index = 0
+    solver_process_size = 1
+    interface = precice.Interface(participant_name, config_file_name, solver_process_index, solver_process_size)
+
+    mesh_name = "Fluid1-Mesh"
+    mesh_id = interface.get_mesh_id(mesh_name)
+
+    velocity_name = "Velocity"
+    velocity_id = interface.get_data_id(velocity_name, mesh_id)
+
+    pressure_name = "Pressure"
+    pressure_id = interface.get_data_id(pressure_name, mesh_id)
+
+    positions = [0, 0, 0]
+
+    vertex_ids = interface.set_mesh_vertex(mesh_id, positions)
+
+    precice_dt = interface.initialize()
+
+
+    while interface.is_coupling_ongoing():
+        if interface.is_action_required(
+            precice.action_write_iteration_checkpoint()):
+
+            u_iter = u
+            p_iter = p
+            t_iter = t
+            counter_iter = counter
+
+            interface.mark_action_fulfilled(
+            precice.action_write_iteration_checkpoint())
+
+        # determine time step size
+
+        dt = 0.025
+        dt = np.minimum(dt,precice_dt)
+
+        # set boundary conditions
+
+        u[0,2] = 1
+        # u[1,2] = 1  # dirichlet velocity inlet
+
+        u[-1,2] = u[-2,2]   # neumann velocity outlet
+
+        p[0] = p[1] # neumann pressure inlet
+        # p[-1] = 0   # dirichlet pressure outlet
+        if interface.is_read_data_available():  # get dirichlet pressure outlet value from 3D solver
+            p_read_in = interface.read_scalar_data(pressure_id, vertex_ids)
+            p[-1] = p_read_in
+        else:
+            p[-1] = 0
+
+        # compute right-hand side of 1D PPE
+
+        for i in range(n-1):
+
+            rhs[i+1] = (1 / dt) * ((u[i+2,2] - u[i,2]) / 2*h)
+
+        # solve the PPE using a SOR solver
+
+        tolerance = 0.001
+        error = 1
+        omega = 1.7
+        max_iter = 1000
+        iter = 0
+
+        while error >= tolerance:
+
+            p[0] = p[1] # renew neumann pressure inlet
+
+            for i in range(n-1):
+                p[i+1] = (1-omega) * p[i+1] + ((omega * h**2) / 2) * (((p[i] + p[i+2]) / h**2) - rhs[i+1])
+
+            sum = 0
+            for i in range(n-1):
+                val = ((p[i] - 2*p[i+1] + p[i+2]) / h**2) - rhs[i+1]
+                sum += val*val
+
+            error = np.sqrt(sum/n)
+
+            iter += 1
+            if iter >= max_iter:
+                print("SOR solver did not converge.\n")
+                break
+
+
+        # calculate new velocities
+
+        for i in range(n-1):
+            u[i+1,2] = u[i+1,2] - dt * ((p[i+2] - p[i+1]) / h)
+
+
+        # write velocity to 3D solver
+
+        write_vel = u[-2,:]
+
+        if interface.is_write_data_required(dt):    # write new velocities to 3D solver
+            interface.write_vector_data(
+            velocity_id, vertex_ids, write_vel)
+
+        # transform data and write output data to vtk files
+
+        u_print = np.reshape(u[:,2], n+1)
+        p_print = np.reshape(p, n+1)
+
+        u_print = np.ascontiguousarray(u_print, dtype=np.float32)
+        p_print = np.ascontiguousarray(p_print, dtype=np.float32)
+
+        filename = "./results/Fluid1D_" + str(counter)
+
+        # TODO: uncomment if pyEVTK is installed and vtk output of 1D participant is wanted
+        # pointsToVTK(filename, x, y, z, data = {"U" : u_print, "p" : p_print})
+
+        # advance simulation time
+
+        dt = interface.advance(dt)
+        t = t + dt
+        counter += 1
+
+        if interface.is_action_required(
+            precice.action_read_iteration_checkpoint()):
+
+            u = u_iter
+            p = p_iter
+            t = t_iter
+            counter = counter_iter
+
+            interface.mark_action_fulfilled(
+            precice.action_read_iteration_checkpoint())
+
+    interface.finalize()
+
+
+if __name__ == "__main__":
+    main()

1D-3D-partitioned-pipe/Fluid1/clean.sh

@@ -0,0 +1,3 @@
+#!/bin/sh
+
+rm ./results/Fluid1D_*

1D-3D-partitioned-pipe/Fluid1/run.sh

@@ -0,0 +1,3 @@
+#!/bin/sh
+
+./Fluid1D.py

1D-3D-partitioned-pipe/Fluid2/0.orig/U

@@ -0,0 +1,30 @@
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            fixedGradient;
+        gradient        uniform (0 0 0);
+    }
+    fixedWalls
+    {
+        type            noSlip;
+    }
+}

1D-3D-partitioned-pipe/Fluid2/0.orig/p

@@ -0,0 +1,31 @@
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedGradient;
+        gradient        uniform 0;
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    fixedWalls
+    {
+        type            zeroGradient;
+    }
+}

1D-3D-partitioned-pipe/Fluid2/0.orig/uniform/functionObjects/functionObjectProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2106                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       dictionary;
+    location    "0/uniform/functionObjects";
+    object      functionObjectProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+
+// ************************************************************************* //

1D-3D-partitioned-pipe/Fluid2/clean.sh

@@ -0,0 +1,46 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+echo "Cleaning..."
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+Participant2="Fluid2"
+
+# Clean the case
+cleanCase
+rm -rfv 0
+# Create an empty .foam file for ParaView
+# Note: ".foam" triggers the native OpenFOAM reader of ParaView.
+# Change to ".OpenFOAM" to use the OpenFOAM reader provided with OpenFOAM.
+touch ${Participant2}.foam
+cd ..
+# Remove the log files
+rm -fv ${Participant2}_blockMesh.log
+rm -fv ${Participant2}_checkMesh.log
+rm -fv ${Participant2}_potentialFoam.log
+rm -fv ${Participant2}_decomposePar.log
+rm -fv ${Participant2}.log
+rm -fv ${Participant2}_reconstructPar.log
+
+
+# Remove the preCICE-related log files
+echo "Deleting the preCICE log files..."
+rm -fv \
+    precice-*.log \
+    precice-*.json \
+
+# Output files for preCICE versions before 1.2:
+rm -fv \
+    iterations-${Participant1}.txt iterations-${Participant2}.txt \
+    convergence-${Participant1}.txt convergence-${Participant2}.txt \
+    Events-${Participant1}.log Events-${Participant2}.log \
+    EventTimings-${Participant1}.log EventTimings-${Participant2}.log
+
+rm -fv .*.address
+
+rm -fv watchpointLeft.txt
+
+echo "Cleaning complete!"
+#------------------------------------------------------------------------------

1D-3D-partitioned-pipe/Fluid2/constant/transportProperties

@@ -0,0 +1,12 @@
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+
+transportModel  Newtonian;
+
+nu              nu [ 0 2 -1 0 0 0 0 ] 1e1;

1D-3D-partitioned-pipe/Fluid2/constant/turbulenceProperties

@@ -0,0 +1,10 @@
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+
+simulationType laminar;

1D-3D-partitioned-pipe/Fluid2/run.sh

@@ -0,0 +1,21 @@
+#!/bin/bash
+
+cd ${0%/*} || exit 1    		    		# Run from this directory
+. $WM_PROJECT_DIR/bin/tools/RunFunctions    # Tutorial run functions
+
+echo "Preparing and running the Fluid2 participant..."
+
+rm -rfv ./0/
+cp -r ./0.orig/ ./0/
+cd ..
+blockMesh -case Fluid2
+checkMesh -case Fluid2
+cd Fluid2
+# Run
+solver=$(getApplication)
+
+cd ..
+$solver -case Fluid2 | tee Fluid2.log 2>&1
+
+foamToVTK -case Fluid2
+