Channel transport with particles

channel-transport-particles/README.md

@@ -0,0 +1,55 @@
+---
+title: Channel transport with particles
+permalink: tutorials-channel-transport-particles.html
+keywords: volume coupling, particles, OpenFOAM, MercuryDPM, transport, just-in-time mapping
+summary: A CFD problem is coupled to a particles in a uni-directional way for particles tracing.
+---
+
+{% note %}
+Get the [case files of this tutorial](https://github.com/precice/tutorials/tree/master/channel-transport-particles). Read how in the [tutorials introduction](https://precice.org/tutorials.html).
+{% endnote %}
+
+## Setup
+
+We model a two-dimensional incompressible fluid flowing through a channel with an obstacle. The fluid problem is coupled to a particle participant for particle tracing. Similar to the transport problem (see the [channel-transport tutorial](tutorials-channel-transport.html)), particles are arranged in a circular blob close to the inflow.
+
+## Configuration
+
+preCICE configuration (image generated using the [precice-config-visualizer](https://precice.org/tooling-config-visualization.html)):
+
+![preCICE configuration visualization](images/tutorials-channel-transport-precice-config.png
+)
+
+## Available solvers
+
+Fluid participant:
+
+* Nutils. For more information, have a look at the [Nutils adapter documentation](https://precice.org/adapter-nutils.html). This Nutils solver requires at least Nutils v7.0.
+
+* OpenFOAM (pimpleFoam). For more information, have a look at the [OpenFOAM adapter documentation](https://precice.org/adapter-openfoam-overview.html).
+
+Particle participant:
+
+* MercuryDPM
+
+## Running the simulation
+
+For the fluid solver, use Nutils for ease of installation and OpenFOAM for speed.
+
+Open two separate terminals and start one fluid and one transport participant by calling the respective run scripts `run.sh` located in each of the participants' directory. For example:
+
+```bash
+cd fluid-nutils
+./run.sh
+```
+
+and either the non-adaptive mesh transport solver
+
+```bash
+cd particles-mercurydpm
+./run.sh
+```
+
+## Post-processing
+
+All solvers generate `vtk` files which can be visualized using, e.g., ParaView.

channel-transport-particles/clean-tutorial.sh

@@ -0,0 +1 @@
+../tools/clean-tutorial-base.sh

channel-transport-particles/fluid-nutils/clean.sh

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

channel-transport-particles/fluid-nutils/fluid.py

@@ -0,0 +1,111 @@
+#! /usr/bin/env python3
+
+#
+# Incompressible NSE solved within a channel geometry with parabolic inflow profile and an obstacle attached to the bottom towards the middle of the domain. The fluid field is initialized with a Stokes solution. The resulting velocity field is written to preCICE on the complete volume.
+#
+
+from nutils import function, mesh, cli, solver, export
+import treelog as log
+import numpy as np
+import precice
+from mpi4py import MPI
+
+
+def main():
+
+    print("Running Nutils")
+
+    # define the Nutils mesh
+    nx = 48
+    ny = 16
+    step_start = nx // 3
+    step_end = nx // 2
+    step_hight = ny // 2
+
+    grid = np.linspace(0, 6, nx + 1), np.linspace(0, 2, ny + 1)
+    domain, geom = mesh.rectilinear(grid)
+    domain = domain.withboundary(inflow="left", outflow="right", wall="top,bottom") - domain[
+        step_start:step_end, :step_hight
+    ].withboundary(wall="left,top,right")
+
+    # cloud of Gauss points
+    gauss = domain.sample("gauss", degree=4)
+
+    # Nutils namespace
+    ns = function.Namespace()
+    ns.x = geom
+
+    ns.ubasis = domain.basis("std", degree=2).vector(2)
+    ns.pbasis = domain.basis("std", degree=1)
+    ns.u_i = "ubasis_ni ?u_n"  # solution
+    ns.p = "pbasis_n ?p_n"  # solution
+    ns.dudt_i = "ubasis_ni (?u_n - ?u0_n) / ?dt"  # time derivative
+    ns.μ = 0.5  # viscosity
+    ns.σ_ij = "μ (u_i,j + u_j,i) - p δ_ij"
+    ns.uin = "10 x_1 (2 - x_1)"  # inflow profile
+
+    # define the weak form, Stokes problem
+    ures = gauss.integral("ubasis_ni,j σ_ij d:x" @ ns)
+    pres = gauss.integral("pbasis_n u_k,k d:x" @ ns)
+
+    # define Dirichlet boundary condition
+    sqr = domain.boundary["inflow"].integral("(u_0 - uin)^2 d:x" @ ns, degree=2)
+    sqr += domain.boundary["inflow,outflow"].integral("u_1^2 d:x" @ ns, degree=2)
+    sqr += domain.boundary["wall"].integral("u_k u_k d:x" @ ns, degree=2)
+    cons = solver.optimize(["u"], sqr, droptol=1e-15)
+
+    # preCICE setup
+    participant = precice.Participant("Fluid", "../precice-config.xml", 0, 1)
+
+    # define coupling mesh
+    mesh_name = "Fluid-Mesh"
+    vertices = gauss.eval(ns.x)
+    vertex_ids = participant.set_mesh_vertices(mesh_name, vertices)
+
+    # coupling data
+    data_name = "Velocity"
+
+    participant.initialize()
+
+    timestep = 0
+    solver_dt = 0.005
+    precice_dt = participant.get_max_time_step_size()
+    dt = min(precice_dt, solver_dt)
+
+    state = solver.solve_linear(("u", "p"), (ures, pres), constrain=cons)  # initial condition
+
+    # add convective term and time derivative for Navier-Stokes
+    ures += gauss.integral("ubasis_ni (dudt_i + μ (u_i u_j)_,j) d:x" @ ns)
+
+    while participant.is_coupling_ongoing():
+
+        if timestep % 1 == 0:  # visualize
+            bezier = domain.sample("bezier", 2)
+            x, u, p = bezier.eval(["x_i", "u_i", "p"] @ ns, **state)
+            with log.add(log.DataLog()):
+                export.vtk("Fluid_" + str(timestep), bezier.tri, x, u=u, p=p)
+
+        precice_dt = participant.get_max_time_step_size()
+
+        # potentially adjust non-matching timestep sizes
+        dt = min(solver_dt, precice_dt)
+
+        # solve Nutils timestep
+        state["u0"] = state["u"]
+        state["dt"] = dt
+        state = solver.newton(("u", "p"), (ures, pres), constrain=cons, arguments=state).solve(1e-10)
+
+        velocity_values = gauss.eval(ns.u, **state)
+        participant.write_data(mesh_name, data_name, vertex_ids, velocity_values)
+
+        # do the coupling
+        participant.advance(dt)
+
+        # advance variables
+        timestep += 1
+
+    participant.finalize()
+
+
+if __name__ == "__main__":
+    cli.run(main)

channel-transport-particles/fluid-nutils/requirements.txt

@@ -0,0 +1,5 @@
+setuptools # required by nutils
+nutils==7
+numpy >1, <2
+pyprecice~=3.0
+setuptools

channel-transport-particles/fluid-nutils/run.sh

@@ -0,0 +1,16 @@
+#!/usr/bin/env bash
+set -e -u
+
+. ../../tools/log.sh
+exec > >(tee --append "$LOGFILE") 2>&1
+
+if [ ! -v PRECICE_TUTORIALS_NO_VENV ]
+then
+    python3 -m venv .venv
+    . .venv/bin/activate
+    pip install -r requirements.txt && pip freeze > pip-installed-packages.log
+fi
+
+python3 fluid.py
+
+close_log

channel-transport-particles/fluid-openfoam/0/U

@@ -0,0 +1,39 @@
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+
+dimensions      [0 1 -1 0 0 0 0];
+internalField   uniform (10 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    obstacle
+    {
+        type            noSlip;
+    }
+    upperWall
+    {
+        type            noSlip;
+    }
+    lowerWall
+    {
+        type            noSlip;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}

channel-transport-particles/fluid-openfoam/0/p

@@ -0,0 +1,45 @@
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    obstacle
+    {
+        type            zeroGradient;
+    }
+
+    upperWall
+    {
+        type            zeroGradient;
+    }
+
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}

channel-transport-particles/fluid-openfoam/clean.sh

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

channel-transport-particles/fluid-openfoam/constant/transportProperties

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

channel-transport-particles/fluid-openfoam/constant/turbulenceProperties

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

channel-transport-particles/fluid-openfoam/run.sh

@@ -0,0 +1,12 @@
+#!/usr/bin/env bash
+set -e -u
+
+. ../../tools/log.sh
+exec > >(tee --append "$LOGFILE") 2>&1
+
+blockMesh
+
+../../tools/run-openfoam.sh "$@"
+. ../../tools/openfoam-remove-empty-dirs.sh && openfoam_remove_empty_dirs
+
+close_log