AxialGeoMultiscale Tutorial

partitioned-pipe-multiscale/README.md

@@ -0,0 +1,40 @@
+---
+title: Partitioned Pipe Multiscale
+permalink: tutorials-partitioned-pipe-multiscale.html
+keywords: OpenFOAM, python
+summary: The 1D-3D Partitioned Pipe is a simple geometric multiscale case, that consists of flow between two pipes with heterogeneous dimensionality. The flow is incompressible and laminar. This tutorial contains an OpenFOAM case for the 3D participant and Python solver for the 1D participant.  
+---
+
+{% note %}
+Get the [case files of this tutorial](https://github.com/ezonta/tutorials/tree/GeoMultiScaleTutorial/1D-3D-partitioned-pipe). Read how in the [tutorials introduction](https://www.precice.org/tutorials.html).
+{% endnote %}
+
+## Prerequisites
+
+- preCICE
+- Python bindings of preCICE
+- OpenFOAM together with the OpenFOAM adapter of preCICE
+
+## Setup
+
+We exchange velocity data from the 1D to the 3D participant and for the pressure data vice versa. The config looks as follows:
+
+![Config Visualization](images/tutorials-partitioned-pipe-multiscale-config.png)
+
+## How to run
+
+In two different terminals execute
+
+```bash
+cd fluid1d-python && ./run.sh
+```
+
+```bash
+cd fluid3d-openfoam && ./run.sh
+```
+
+## Results
+
+You should be able to see an established laminar profile at the inlet of the 3D participant.
+
+![Expected Result](images/tutorials-partitioned-pipe-multiscale-profile.png)

partitioned-pipe-multiscale/fluid1d-python/Fluid1D.py

@@ -0,0 +1,169 @@
+#!/usr/bin/env python3
+
+from numpy.core.fromnumeric import reshape
+# TODO: uncomment for vtk output
+# 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 = "Fluid1D"
+    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 = "Fluid1D-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()

partitioned-pipe-multiscale/fluid1d-python/clean.sh

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

partitioned-pipe-multiscale/fluid1d-python/run.sh

@@ -0,0 +1,4 @@
+#!/bin/sh
+
+cd ..
+./fluid1d-python/Fluid1D.py

partitioned-pipe-multiscale/fluid3d-openfoam/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;
+    }
+}

partitioned-pipe-multiscale/fluid3d-openfoam/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;
+    }
+}

partitioned-pipe-multiscale/fluid3d-openfoam/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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+
+// ************************************************************************* //

partitioned-pipe-multiscale/fluid3d-openfoam/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="Fluid3D"
+
+# 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!"
+#------------------------------------------------------------------------------

partitioned-pipe-multiscale/fluid3d-openfoam/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;

partitioned-pipe-multiscale/fluid3d-openfoam/constant/turbulenceProperties

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