AxialGeoMultiscale Tutorial

partitioned-pipe-multiscale/README.md

@@ -0,0 +1,34 @@
+---
+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, coupling two pipes with different dimensions. 
+---
+
+{% note %}
+Get the [case files of this tutorial](https://github.com/precice/tutorials/tree/master/partitioned-pipe-multiscale). Read how in the [tutorials introduction](https://www.precice.org/tutorials.html).
+{% endnote %}
+
+## Setup
+
+We exchange velocity data from the upstream 1D to the downstream 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
+
+Visualizing the results in ParaView, we 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,112 @@
+#!/usr/bin/env python3
+
+from nutils import mesh, function, solver, cli
+from nutils.expression_v2 import Namespace
+import numpy
+import treelog
+import matplotlib.pyplot as plt
+import precice
+
+
+def main(nelems: int, etype: str, degree: int, reynolds: float):
+    '''
+    1D channel flow problem.
+    .. arguments::
+       nelems [12]
+         Number of elements along edge.
+       etype [square]
+         Element type (square/triangle/mixed).
+       degree [2]
+         Polynomial degree for velocity; the pressure space is one degree less.
+       reynolds [1000]
+         Reynolds number, taking the domain size as characteristic length.
+    '''
+    # 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()
+
+    # problem definition
+    domain, geom = mesh.rectilinear([numpy.linspace(0, 1, nelems + 1)])
+    ns = Namespace()
+    ns.δ = function.eye(domain.ndims)
+    ns.Σ = function.ones([domain.ndims])
+    ns.Re = reynolds
+    ns.x = geom
+    ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
+    ns.ubasis = domain.basis('std', degree=2).vector(domain.ndims)
+    ns.pbasis = domain.basis('std', degree=1)
+    ns.u = function.dotarg('u', ns.ubasis)
+    ns.p = function.dotarg('p', ns.pbasis)
+    ns.stress_ij = '(∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij'
+    ures = domain.integral('∇_j(ubasis_ni) stress_ij dV' @ ns, degree=4)
+    pres = domain.integral('pbasis_n ∇_k(u_k) dV' @ ns, degree=4)
+
+    while interface.is_coupling_ongoing():
+        if interface.is_read_data_available():  # get dirichlet pressure outlet value from 3D solver
+            p_read = interface.read_scalar_data(pressure_id, vertex_ids)
+            p_read = numpy.maximum(0, p_read)     # filter out unphysical negative pressure values
+        else:
+            p_read = 0
+        usqr = domain.boundary['left'].integral('(u_0 - 1)^2 dS' @ ns, degree=2)
+        diricons = solver.optimize('u', usqr, droptol=1e-15)
+        ucons = diricons
+        stringintegral = '(p - ' + str(numpy.rint(p_read)) + ')^2 dS'
+        psqr = domain.boundary['right'].integral(stringintegral @ ns, degree=2)
+        pcons = solver.optimize('p', psqr, droptol=1e-15)
+        cons = dict(u=ucons, p=pcons)
+        with treelog.context('stokes'):
+            state0 = solver.solve_linear(('u', 'p'), (ures, pres), constrain=cons)
+            x, u, p = postprocess(domain, ns, precice_dt, **state0)
+
+        if interface.is_action_required(
+                precice.action_write_iteration_checkpoint()):
+            u_iter = u
+            p_iter = p
+            interface.mark_action_fulfilled(
+                precice.action_write_iteration_checkpoint())
+
+        write_vel = [0, 0, u[-1]]
+        if interface.is_write_data_required(precice_dt):    # write new velocities to 3D solver
+            interface.write_vector_data(velocity_id, vertex_ids, write_vel)
+        precice_dt = interface.advance(precice_dt)
+
+        if interface.is_action_required(
+                precice.action_read_iteration_checkpoint()):
+            u = u_iter
+            p = p_iter
+            interface.mark_action_fulfilled(
+                precice.action_read_iteration_checkpoint())
+
+    interface.finalize()
+    return state0
+
+
+def postprocess(domain, ns, dt, **arguments):
+
+    bezier = domain.sample('bezier', 9)
+    x, u, p = bezier.eval(['x_i', 'u_i', 'p'] @ ns, **arguments)
+
+    ax1 = plt.subplot(211)
+    ax1.plot(x, u)
+    ax1.set_title("velocity u")
+    ax2 = plt.subplot(212)
+    ax2.plot(x, p)
+    ax2.set_title("pressure p")
+    plt.savefig("./results/Fluid1D_" + str(dt) + ".png")
+    return x, u, p
+
+
+if __name__ == '__main__':
+    cli.run(main)

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

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

partitioned-pipe-multiscale/fluid1d-python/precice-Fluid1D-iterations.log

@@ -0,0 +1,5 @@
+TimeWindow  TotalIterations  Iterations  Convergence  
+     1       2       2       1  
+     2       4       2       1  
+     3       6       2       1  
+     4       8       2       1  

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

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

partitioned-pipe-multiscale/fluid3d-openfoam/0/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/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/clean.sh

@@ -0,0 +1,6 @@
+#!/bin/sh
+set -e -u
+
+. ../../tools/cleaning-tools.sh
+
+clean_openfoam .

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;

partitioned-pipe-multiscale/fluid3d-openfoam/run.sh

@@ -0,0 +1,9 @@
+#!/bin/sh
+set -e -u
+
+blockMesh
+touch Fluid3D.foam
+
+../../tools/run-openfoam.sh "$@"
+. ../../tools/openfoam-remove-empty-dirs.sh && openfoam_remove_empty_dirs
+

partitioned-pipe-multiscale/fluid3d-openfoam/system/blockMeshDict

@@ -0,0 +1,104 @@
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant/polyMesh";
+    object      blockMeshDict;
+}
+
+convertToMeters 1;
+
+
+zmin  20;
+zmax  40;
+
+xcells 10; // per block (5 blocks)
+ycells 10; // per block (5 blocks)
+zcells 20;
+
+vertices
+(
+    (-3.535534 -3.535534 $zmin) // 0
+    ( 3.535534 -3.535534 $zmin) // 1
+    ( 3.535534  3.535534 $zmin) // 2
+    (-3.535534  3.535534 $zmin) // 3
+    (-1.414214 -1.414214 $zmin) // 4
+    ( 1.414214 -1.414214 $zmin) // 5
+    ( 1.414214  1.414214 $zmin) // 6
+    (-1.414214  1.414214 $zmin) // 7
+
+    (-3.535534 -3.535534 $zmax) // 8
+    ( 3.535534 -3.535534 $zmax) // 9
+    ( 3.535534  3.535534 $zmax) // 10
+    (-3.535534  3.535534 $zmax) // 11
+    (-1.414214 -1.414214 $zmax) // 12
+    ( 1.414214 -1.414214 $zmax) // 13
+    ( 1.414214  1.414214 $zmax) // 14
+    (-1.414214  1.414214 $zmax) // 15
+);
+
+blocks
+(
+hex (0 1 5 4 8 9 13 12)   pipe ($xcells $ycells $zcells)  edgeGrading (1 1 1 1 1 1 1 1 1 1 1 1)
+hex (1 2 6 5 9 10 14 13)  pipe ($xcells $ycells $zcells)  edgeGrading (1 1 1 1 1 1 1 1 1 1 1 1)
+hex (2 3 7 6 10 11 15 14) pipe ($xcells $ycells $zcells)  edgeGrading (1 1 1 1 1 1 1 1 1 1 1 1)
+hex (3 0 4 7 11 8 12 15)  pipe ($xcells $ycells $zcells)  edgeGrading (1 1 1 1 1 1 1 1 1 1 1 1)
+hex (4 5 6 7 12 13 14 15) pipe ($xcells $ycells $zcells)  edgeGrading (1 1 1 1 1 1 1 1 1 1 1 1)
+);
+
+edges
+(
+arc 0 1   ( 0   -5   $zmin)
+arc 1 2   ( 5    0   $zmin)
+arc 2 3   ( 0    5   $zmin)
+arc 3 0   (-5    0   $zmin)
+arc 4 5   ( 0   -1.5 $zmin)
+arc 5 6   ( 1.5  0   $zmin)
+arc 6 7   ( 0    1.5 $zmin)
+arc 7 4   (-1.5  0   $zmin)
+
+arc 8 9   ( 0  -5   $zmax)
+arc 9 10  ( 5   0   $zmax)
+arc 10 11 ( 0   5   $zmax)
+arc 11 8  (-5   0   $zmax)
+arc 12 13 ( 0  -1.5 $zmax)
+arc 13 14 (1.5  0   $zmax)
+arc 14 15 ( 0   1.5 $zmax)
+arc 15 12 (-1.5 0   $zmax)
+);
+
+
+patches
+(
+
+    patch fixedWalls
+    (
+        (0 1 9   8)
+        (1 2 10  9)
+        (2 3 11 10)
+        (3 0 8  11)
+    )
+
+    patch inlet
+    (
+        (0 1 5 4)
+        (1 2 6 5)
+        (2 3 7 6)
+        (3 0 4 7)
+        (4 5 6 7)
+    )
+
+    patch outlet
+    (
+        (8 9 13   12)
+        (9 10 14  13)
+        (10 11 15 14)
+        (12 13 14 15)
+        (11 8 12  15)
+    )
+);
+
+mergePatchPairs
+(
+);