Skip to content

AxialGeoMultiscale Tutorial #308

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 20 commits into
base: develop
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 34 additions & 0 deletions partitioned-pipe-multiscale/README.md
Original file line number Diff line number Diff line change
@@ -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)
11 changes: 11 additions & 0 deletions partitioned-pipe-multiscale/clean-tutorial.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@
#!/usr/bin/env sh
set -e -u

# shellcheck disable=SC1091
. ../tools/cleaning-tools.sh

clean_tutorial .
clean_precice_logs .
rm -fv ./*.log
rm -fv ./*.vtu

105 changes: 105 additions & 0 deletions partitioned-pipe-multiscale/fluid1d-python/Fluid1D.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,105 @@
#!/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
participant = precice.Participant(participant_name, config_file_name, solver_process_index, solver_process_size)
mesh_name = "Fluid1D-Mesh"
velocity_name = "Velocity"
pressure_name = "Pressure"
positions = numpy.zeros((1, 3)) # one vertex of three dimensions
vertex_ids = participant.set_mesh_vertices(mesh_name, positions)
participant.initialize()
precice_dt = participant.get_max_time_step_size()

# problem definition
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you write a few sentences on what exactly you are trying to solve here?
I have zero experience with nutils, so it could be useful if at least @IshaanDesai could also do a quick check. I think you have discussed it already.

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 participant.is_coupling_ongoing():
# get dirichlet pressure outlet value from 3D solver
p_read = participant.read_data(mesh_name, pressure_name, vertex_ids, precice_dt)
# filter out unphysical negative pressure values
p_read = numpy.maximum(0, p_read)

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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I barely understand what is going on here, but a tolerance of 1e-15 sounds too strict to me.

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 participant.requires_writing_checkpoint():
u_iter = u
p_iter = p

write_vel = [[0, 0, u[-1][0]]]
# write new velocities to 3D solver
participant.write_data(mesh_name, velocity_name, vertex_ids, write_vel)
participant.advance(precice_dt)
precice_dt = participant.get_max_time_step_size()

if participant.requires_reading_checkpoint():
u = u_iter
p = p_iter

participant.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)
8 changes: 8 additions & 0 deletions partitioned-pipe-multiscale/fluid1d-python/clean.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,8 @@
#!/usr/bin/env sh
set -e -u

. ../../tools/cleaning-tools.sh

rm -f ./results/Fluid1D_*
clean_precice_logs .
clean_case_logs .
5 changes: 5 additions & 0 deletions partitioned-pipe-multiscale/fluid1d-python/requirements.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
setuptools # required by nutils
nutils==7
numpy >1, <2
pyprecice~=3.0
matplotlib
2 changes: 2 additions & 0 deletions partitioned-pipe-multiscale/fluid1d-python/results/.gitignore
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
# Keep this directory, but ignore all files
*
13 changes: 13 additions & 0 deletions partitioned-pipe-multiscale/fluid1d-python/run.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
#!/usr/bin/env bash
set -e -u

. ../../tools/log.sh
exec > >(tee --append "$LOGFILE") 2>&1

python3 -m venv .venv
. .venv/bin/activate
pip install -r requirements.txt && pip freeze > pip-installed-packages.log

NUTILS_RICHOUTPUT=no python3 Fluid1D.py

close_log
30 changes: 30 additions & 0 deletions partitioned-pipe-multiscale/fluid3d-openfoam/0/U
Original file line number Diff line number Diff line change
@@ -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;
}
}
31 changes: 31 additions & 0 deletions partitioned-pipe-multiscale/fluid3d-openfoam/0/p
Original file line number Diff line number Diff line change
@@ -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;
}
}
6 changes: 6 additions & 0 deletions partitioned-pipe-multiscale/fluid3d-openfoam/clean.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
#!/usr/bin/env sh
set -e -u

. ../../tools/cleaning-tools.sh

clean_openfoam .
Original file line number Diff line number Diff line change
@@ -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 ] 10;
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}

simulationType laminar;
12 changes: 12 additions & 0 deletions partitioned-pipe-multiscale/fluid3d-openfoam/run.sh
Original file line number Diff line number Diff line change
@@ -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
Loading