Add Nutils participants for turek-hron-fsi3 (#648)

* Add Nutils participants for turek-hron-fsi3

* Add python requirements

* Change Nutils cases from FSI2 to FSI3

* Remove dummy mock from solid.py

* Adjust README and plotting script

* Update turek-hron-fsi3/README.md

Co-authored-by: Gerasimos Chourdakis <gerasimos.chourdakis@ipvs.uni-stuttgart.de>

* Restructure time loop

This patch conforms the code to the standard control flow order of
requires_writing_checkpoint, requires_reading_checkpoint and
is_time_window_complete.

* Prepare for dynamic timesteps

This patch removes the Dynamic.set_timestep method, which is not a good fit in
a situation where time steps may vary. Instead the add_and_plot method is
rewritten to use dimensional times rather than a fixed-length deque, and the
newmark_defo and newmark_velo methods receive a timestep argument to replace
the attribute. The latter will be further modified in a subsequent commit to
allow for truly dynamic time steps.

* Rename t->traction to make place for time variable

* Update precice timestep at every iteration

* Acknowledgements

---------

Co-authored-by: Benjamin Uekermann <benjamin.uekermann@ipvs.uni-stuttgart.de>
Co-authored-by: Benjamin Uekermann <benjamin.uekermann@gmail.com>
Co-authored-by: Gerasimos Chourdakis <gerasimos.chourdakis@ipvs.uni-stuttgart.de>

Author: 3 people

turek-hron-fsi3/README.md

@@ -1,8 +1,8 @@
 ---
 title: Turek Hron FSI3
 permalink: tutorials-turek-hron-fsi3.html
-keywords: OpenFOAM, deal.II, verification
-summary: The Turek-Hron FSI cases are well-established numerical benchmarks and, therefore, well suited for verification of preCICE itself and the used adapters. In this tutorial, we focus on the FSI3 case, which presents the most challenging case in terms of added mass. Please note that the meshes of this case are significantly finer than for other tutorials. Running the simulation might take a few hours. We do not recommend to run this tutorials as your first preCICE tutorial.  
+keywords: OpenFOAM, deal.II, Nutils, verification
+summary: The Turek-Hron FSI cases are well-established numerical benchmarks and, therefore, well suited for verification of preCICE itself and the used adapters. In this tutorial, we focus on the FSI3 case, which presents the most challenging case in terms of added mass. Please note that the meshes of this case are significantly finer than for other tutorials. Running the simulation might take a few hours. We do not recommend to run this tutorials as your first preCICE tutorial.
 ---
 
 {% note %}
@@ -28,14 +28,16 @@ preCICE configuration (image generated using the [precice-config-visualizer](htt
 Fluid participant:
 
 * OpenFOAM (pimpleFoam). In case you are using a very old OpenFOAM version, you will need to adjust the solver to `pimpleDyMFoam` in the `Fluid/system/controlDict` file. For more information, have a look at the [OpenFOAM adapter documentation](https://precice.org/adapter-openfoam-overview.html).
+* 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 v9.0. This case takes significantly longer to run than the OpenFOAM case, see [related issue](https://github.com/precice/tutorials/issues/506).
 
 Solid participant:
 
 * deal.II. For more information, have a look at the [deal.II adapter documentation](https://precice.org/adapter-dealii-overview.html). This tutorial requires the nonlinear solid solver. Please copy the nonlinear solver executable to the `solid-dealii` folder or make it discoverable at runtime and update the `solid-dealii/run.sh` script.
+* 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 v9.0.
 
 ## Running the Simulation
 
-Open two separate terminals and start each participant by calling the respective run script.
+Open two separate terminals and start each participant by calling the respective run script. For example:
 
 ```bash
 cd fluid-openfoam
@@ -57,13 +59,16 @@ In the first few timesteps, many coupling iterations are required for convergenc
 
 ## Post-processing
 
-You can visualize the results of the coupled simulation using e.g. ParaView. Fluid results are in the OpenFOAM format and you may load the `fluid-openfoam.foam` file. Solid results are in VTK format.
+You can visualize the results of the coupled simulation using e.g. ParaView. OpenFOAM uses an OpenFOAM-specific format, and you can directly load the (empty) file `fluid-openfoam.foam` in Paraview or convert the results to VTK with `foamToVTK`. deal.II writes VTK files. Both Nutils solvers currently do not write VTK files (but use their own in-situ visualization), but this can be easily added similarly to the [perpendicular flap solvers](https://github.com/precice/tutorials/blob/98a78fe2dc2f6c5d84b2b30d35d00352782236f8/perpendicular-flap/fluid-nutils/fluid.py#L227).
 
-If you want to visualize both domains with ParaView, keep in mind that the deal.II solver writes results every few timesteps, while the OpenFOAM solver writes in reference to simulated time. For this reason, make sure that you use compatible write intervals. You may also need to convert the OpenFOAM results to VTK (with the command `foamToVTK`).
+If you want to visualize both domains with ParaView, keep in mind that the different solvers may write results with different output frequencies, which you might want to [synchronize](https://precice.org/tutorials-visualization.html#synchronizing-results).
 
-There is an [known issue](https://github.com/precice/openfoam-adapter/issues/26) that leads to additional "empty" result directories when running with some OpenFOAM versions, leading to inconveniences during post-processing. At the end of `run.sh`, we call `openfoam_remove_empty_dirs` (provided by `tools/openfoam-remove-empty-dirs`) to delete the additional files before importing the results in ParaView.
+There is a [known issue](https://github.com/precice/openfoam-adapter/issues/26) that leads to additional "empty" result directories when running with some OpenFOAM versions, leading to inconveniences during post-processing. At the end of `run.sh`, we call `openfoam_remove_empty_dirs` (provided by `tools/openfoam-remove-empty-dirs`) to delete the additional files before importing the results in ParaView.
 
-Moreover, as we defined a watchpoint at the flap tip (see `precice-config.xml`), we can plot it with gnuplot using the script `plot-displacement.sh`.  The resulting graph shows the vertical (y) displacement of the tip of the flap.
+Moreover, as we defined a watchpoint at the flap tip (see `precice-config.xml`), we can plot it with gnuplot using the script `plot-displacement.sh`, which expects the directory of the selected solid participant as a command line argument. For example:
+
+ ```shell
+ plot-displacement.sh solid-dealii
 
 ![FSI3 watchpoint](images/tutorials-turek-hron-fsi3-tip-plot.png)
 
@@ -84,14 +89,10 @@ mv blockMeshDict blockMeshDict_original
 mv blockMeshDict_refined blockMeshDict
 ```
 
-For the double-refined mesh, it is wisely to use local basis functions in the RBF data mapping method instead of global ones. You can use:
-
-```xml
-<mapping:rbf-compact-tps-c2 direction="read" from="Fluid-Mesh-Centers" to="Solid-Mesh"
-                            support-radius="0.011" constraint="consistent" />
-```
+## Acknowledgements
 
-You can find more information on RBF data mapping in the [documentation](https://precice.org/configuration-mapping.html#radial-basis-function-mapping).
+Thanks to the Technical University of Eindhoven for funding the development of
+the Nutils participants for this tutorial.
 
 ## References
 

turek-hron-fsi3/fluid-nutils/clean.sh

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

turek-hron-fsi3/fluid-nutils/fluid.geo

@@ -0,0 +1,79 @@
+// Geometry file for the turek.py example.
+//
+// This is a generalized description of the setup defined by Turek and Hron,
+// which requires the following lenghts to be supplied externally, using the
+// numbers argument of mesh.gmsh or the -setnumber switch when invoking the
+// gmsh application directly:
+//
+// - channel_length: length of the fluid domain
+// - channel_height: height of the fluid domain
+// - x_center: horizontal position of the cylinder measured from the left edge
+// - y_center: vertical position of the cylinder measured from the bottom edge
+// - cylinder_radius: radius of the cylinder
+// - structure_length: length of the elastic structure measured from the cylinder wall
+// - structure_thickness: thickness of the elastic structure
+// - min_elemsize: mesh element size at the solid/fluid interface
+// - max_elemsize: mesh element size at the channel wall and far field
+//
+// The parameterization matches largely that of Table 1 of Turek and Hron 2006,
+// with the main difference that reference points A and B cannot be
+// independently placed but are always located at the tip of the elastic
+// structure and the leading edge of the cylinder, respectively.
+
+SetFactory("OpenCASCADE");
+
+Rectangle(1) = {0, 0, 0, channel_length, channel_height};
+Rectangle(2) = {x_center, y_center - structure_thickness/2, 0, cylinder_radius + structure_length, structure_thickness, 0};
+Disk(3) = {x_center, y_center, 0, cylinder_radius};
+BooleanDifference(4) = { Surface{2}; }{ Surface{3}; };
+BooleanDifference(5) = { Surface{1}; }{ Surface{2,3}; };
+A = newp; Point(A) = {x_center + cylinder_radius + structure_length, y_center, 0};
+B = newp; Point(B) = {x_center - cylinder_radius, y_center, 0};
+
+// At this point surface 3 (cylinder), 4 (solid domain) and 5 (fluid domain) are
+// non-overlapping. Gmsh promises that the boolean fragments operation with
+// deletion will reuse the surface IDs for the new objects.
+
+_() = BooleanFragments{ Surface{3,4,5}; Point{A,B}; Delete; }{};
+
+// Fragments deduplicates boundary segments, which means that we can now
+// perform boolean operations on the index sets.
+
+bnd_cylinder() = Abs(Boundary{ Surface{3}; });
+bnd_structure() = Abs(Boundary{ Surface{4}; });
+tmp = bnd_structure();
+bnd_structure -= bnd_cylinder();
+bnd_cylinder -= tmp();
+bnd_fluid() = Abs(Boundary{ Surface{5}; });
+bnd_fluid -= bnd_structure();
+bnd_fluid -= bnd_cylinder();
+
+// After subtracting the inner boundaries, only the four boundary segments of
+// rectangle 1 remain in bnd_fluid, and we are going to assume that they are
+// ordered bottom, right, top, left.
+
+Physical Surface("fluid") = {5};
+Physical Line("inlet") = {bnd_fluid(3)};
+Physical Line("outlet") = {bnd_fluid(1)};
+Physical Line("wall") = {bnd_fluid(0), bnd_fluid(2)};
+Physical Line("cylinder") = {bnd_cylinder()};
+Physical Line("structure") = {bnd_structure()};
+Physical Point("A") = {A};
+Physical Point("B") = {B};
+
+// The element size is set to be uniformly min_elemsize inside the elastic
+// structure, and grow linearly to max_elemsize in the fluid domain over a
+// distance of half the channel height.
+
+Mesh.MeshSizeFromPoints = 0;
+Mesh.MeshSizeFromCurvature = 0;
+Mesh.MeshSizeExtendFromBoundary = 0;
+Field[1] = Distance;
+Field[1].SurfacesList = {3,4};
+Field[2] = Threshold;
+Field[2].InField = 1;
+Field[2].DistMin = 0;
+Field[2].DistMax = channel_height/2;
+Field[2].SizeMin = min_elemsize;
+Field[2].SizeMax = max_elemsize;
+Background Field = 2;