WIP: aste Replay Mode tutorial on elastic flap case

perpendicular-flap/.gitignore

@@ -0,0 +1,5 @@
+preCICE-output/
+Solid-fenics*
+*.log
+vtk*
+Solid/FSI-S

perpendicular-flap/fluid-aste/README.md

@@ -0,0 +1,96 @@
+# Tutorial for replay mode
+
+## Overview
+
+This tutorial is an example how the Artificial Solver Testing Environment (aste) can be used. The idea is that you run a coupled simulation with two regular solvers and save the coupling data in every timestep. Then, this data is converted to aste-format. Finally, aste replaces one of the solvers and the simulation can be "replayed" using only one solver and aste with the previously computed results.
+
+This tutorial uses the results from the [OpenFOAM-FEniCS perpendicular flap tutorial](https://github.com/precice/tutorials/tree/master/perpendicular-flap) as a basis. aste then replaces OpenFOAM.
+
+## Requirements
+
+To run this tutorial you need to install the following components:
+
+- [preCICE](https://github.com/precice/precice/wiki/Get-preCICE)
+- [aste](https://github.com/precice/aste/tree/develop)
+- [FEniCS](https://fenicsproject.org/)
+- [FEniCS-Adapter](https://github.com/precice/fenics-adapter)
+- OpenFOAM, e.g. [OpenFOAM 7](https://openfoam.org/version/7/)
+- [OpenFOAM Adapter](https://github.com/precice/openfoam-adapter/wiki/Building) matching the OpenFOAM version.
+
+Make sure to add `aste/build` to the `PATH` such that the python scripts and `preciceMap` can be found from anywhere on your system.
+
+## Step-by-Step explanations
+
+### Generating vtk output during a simulation
+
+The base case for this tutorial is the OpenFOAM-FEniCS perpendicular flap tutorial. So let's start in the the root-directory of the perpendicular flap case: [`tutorials/perpendicular-flap`](https://github.com/precice/tutorials/tree/master/perpendicular-flap).
+
+To generate vtk output in preCICE, add the statement `<export:vtk directory="preCICE-output" />` to the `precice-config.xml` in `tutorials/perpendicular-flap` in the Solid participant. The result should look the following way:
+
+```xml
+<participant name="Solid">
+    <use-mesh name="Solid-Mesh" provide="yes"/>
+    <read-data name="Force" mesh="Solid-Mesh"/>
+    <write-data name="Displacement" mesh="Solid-Mesh"/>
+    <export:vtk directory="preCICE-output" />
+</participant>
+```
+
+Then, run the simulation as explained in the [`README.md`](https://github.com/precice/tutorials/blob/develop/perpendicular-flap/README.md) of the case.
+
+The exports can be found in the `preCICE-output` directory.
+
+### Converting the output to aste format
+
+Copy the `preCICE-output` folder to the root directory of the aste tutorial `tutorials/perpendicular-flap/fluid-aste`.
+To convert the files to the correct format, open the `preCICE-output` folder and run
+
+`precice_to_aste.py Solid-fenics -n 500 -t Forces0 --datadim 3`
+
+### Replay of the simulation with aste
+
+#### The quick way to run
+
+1. Open two terminals and navigate to `solid-fenics` and `fluid-aste`.
+1. Run `.../solid-fenics$ python3 solid.py` and `.../fluid-aste$ preciceMap -v -c ../precice-config-aste.xml -p Fluid --mesh preCICE-output/Solid-fenics` in two terminals.
+
+Read on if you want to know what to change in the configuration files starting from the OpenFOAM-FEniCS tutorial.
+
+#### The way to do it yourself
+
+The aste-FEniCS tutorial provides its own `precice-config.xml` and `precice-adapter-config-fsi-s.json` for getting started quickly. They are called `precice-config-aste.xml` and `precice-adapter-config-fsi-s-aste.json`.
+
+However, these can also be generated by modifiying the files from the original `perpendicular-flap` case. The main purpose of this tutorial is to explain how the replay mode with aste can be used for arbirtrary simulation setups.
+
+1. Copy the files `perpendicular-flap/precice-config.xml` and `perpendicular-flap/solid-fenics/precice-adapter-config-fsi-s.json` from the original `perpendicular-flap` tutorial into the respective directories.
+2. Remove the line `"write_data_name": "Displacement"` from `precice-adapter-config-fsi-s.json`, because we will not write any data to aste.
+3. Change the `read_data_name` from `Force` to `Data` since aste works with the general data name `Data`.
+4. In ```precice-config.xml``` we require some more changes:
+    1. Rename `Force` to `Data` throughout the whole file.
+    2. Delete all lines where `Displacement` occurs, since we are only coupling in one direction (from aste to FEniCS).
+5. Then we process the Fluid meshes to fake the Fluid participant using aste:
+    1. Rename `Fluid-Mesh` to `MeshA` throughout the whole file.
+    2. Delete the `mesh` `"Fluid-Mesh-Nodes"` and all lines where `"Fluid-Mesh-Nodes"` occurs.
+    3. Rename the participant `Fluid` to `A`
+6. Change the coupling scheme to an explicit scheme:
+    1. Change `serial-implicit` to `serial-explicit`.
+    2. Explicit schemes don't iterate such that they don't need convergence measures, extrapolation or acceleration. Delete everything in the coulping scheme except for
+
+    ```xml
+      <participants first="A" second="Solid" />
+      <max-time value="5.0" />
+      <time-window-size value="1.e-2" />
+      <exchange data="Data" mesh="Solid-Mesh" from="A" to="Solid" />
+    ```
+
+#### Run
+
+After applying these preparations, we are now able to run the tutorial in two terminals as described above.
+
+### Play Around
+
+If you want to explore more possibilities of the Replay-Mode here are some ideas:
+
+- Export the output at the Fluid participant. Instead of `preCICE-output/Solid-fenics.*` the files containing the forces will be `preCICE-output/Fluid-Mesh-Faces-Fluid.*`
+- Export the forces in the [Calculix-Openfoam-tutorial](https://github.com/precice/tutorials/tree/master/perpendicular-flap) with the same geometry and feed them to FEniCS with aste.
+- If you are familiar with the OpenFOAM-Adapter you can also replace FEniCS with aste such that aste writes previously exported displacements to OpenFOAM.

perpendicular-flap/precice-config-aste.xml

@@ -0,0 +1,50 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+
+<precice-configuration>
+
+  <log enabled="true">
+    <sink type="stream" output="stdout" filter="%Severity% > debug" enabled="true" />
+    <sink type="file" output="filtered.log"  filter= "(%Severity% > debug) or (%Severity% >= trace and %Module% contains SolverInterfaceImpl)" enabled="true" />
+    <sink type="file" output="debug.log" filter="" enabled="true"/>
+  </log>
+
+  <solver-interface dimensions="2">
+
+    <!-- Data fields that are exchanged between the solvers -->
+    <data:vector name="Force" />
+
+    <!-- A common mesh that uses these data fields -->
+    <mesh name="Fluid-Mesh" flip-normals="no">
+      <use-data name="Force" />
+    </mesh>
+
+    <mesh name="Solid-Mesh" flip-normals="no">
+      <use-data name="Force" />
+    </mesh>
+
+    <participant name="Fluid">
+      <use-mesh name="Fluid-Mesh" provide="yes" />
+      <use-mesh name="Solid-Mesh" provide="no" from="Solid" />
+      <write-data name="Force" mesh="Fluid-Mesh"/>
+      <mapping:nearest-neighbor constraint="conservative" direction="write" from="Fluid-Mesh" to="Solid-Mesh" />
+      <export:vtk every-n-time-windows="1" directory="vtkFluid/" normals="0"/>
+    </participant>
+
+    <participant name="Solid">
+      <use-mesh name="Solid-Mesh" provide="yes" />
+      <read-data name="Force" mesh="Solid-Mesh" />
+      <export:vtk every-n-time-windows="1" directory="vtkSolid/" normals="0"/>
+    </participant>
+
+    <m2n:sockets from="Fluid" to="Solid" exchange-directory=".." />
+
+    <coupling-scheme:serial-explicit>
+      <participants first="Fluid" second="Solid" />
+      <max-time value="5.0" />
+      <time-window-size value="1.e-2" />
+      <exchange data="Force" mesh="Solid-Mesh" from="Fluid" to="Solid" />
+    </coupling-scheme:serial-explicit>
+
+  </solver-interface>
+
+</precice-configuration>

perpendicular-flap/precice-config.xml

@@ -42,10 +42,19 @@
       <use-mesh name="Solid-Mesh" provide="yes" />
       <write-data name="Displacement" mesh="Solid-Mesh" />
       <read-data name="Force" mesh="Solid-Mesh" />
-      <watch-point mesh="Solid-Mesh" name="Flap-Tip" coordinate="0.0;1" />
+      <export:vtk directory="preCICE-output" />
     </participant>
 
     <m2n:sockets from="Fluid" to="Solid" exchange-directory=".." />
+    <!--
+    <coupling-scheme:serial-explicit>
+      <time-window-size value="1e-2"/>
+      <max-time value="5"/>
+      <participants first="Fluid" second="Solid"/>
+      <exchange data="Force" mesh="Solid-Mesh" from="Fluid" to="Solid"/>
+      <exchange data="Displacement" mesh="Solid-Mesh" from="Solid" to="Fluid"/>
+    </coupling-scheme:serial-explicit>
+    -->
 
     <coupling-scheme:parallel-implicit>
       <time-window-size value="0.01" />

perpendicular-flap/solid-fenics/precice-adapter-config-fsi-s-aste.json

@@ -0,0 +1,8 @@
+{
+  "participant_name": "Solid",
+  "config_file_name": "../precice-config-aste.xml",
+  "interface": {
+      "coupling_mesh_name": "Solid-Mesh",
+      "read_data_name": "Force"
+    }
+}

perpendicular-flap/solid-fenics/solid.py

@@ -1,7 +1,7 @@
 # Import required libs
 from fenics import Constant, Function, AutoSubDomain, RectangleMesh, VectorFunctionSpace, interpolate, \
     TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, project, \
-    Identity, inner, dx, ds, sym, grad, lhs, rhs, dot, File, solve, PointSource, assemble_system
+    Identity, inner, dx, ds, sym, grad, lhs, rhs, dot, File, solve, PointSource, assemble_system, TensorFunctionSpace, XDMFFile
 from ufl import nabla_div
 import numpy as np
 import matplotlib.pyplot as plt
@@ -66,12 +66,12 @@ def neumann_boundary(x, on_boundary):
 coupling_boundary = AutoSubDomain(neumann_boundary)
 fixed_boundary = AutoSubDomain(clamped_boundary)
 
-precice = Adapter(adapter_config_filename="precice-adapter-config-fsi-s.json")
+precice = Adapter(adapter_config_filename="precice-adapter-config-fsi-s-aste.json")
 
 # Initialize the coupling interface
-precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V, fixed_boundary=fixed_boundary)
+precice_dt = precice.initialize(coupling_boundary, read_function_space=V, fixed_boundary=fixed_boundary)
 
-fenics_dt = precice_dt  # if fenics_dt == precice_dt, no subcycling is applied
+fenics_dt = precice_dt # if fenics_dt == precice_dt, no subcycling is applied
 # fenics_dt = 0.02  # if fenics_dt < precice_dt, subcycling is applied
 dt = Constant(np.min([precice_dt, fenics_dt]))
 
@@ -168,6 +168,12 @@ def avg(x_old, x_new, alpha):
 # parameters for Time-Stepping
 t = 0.0
 n = 0
+
+time = []
+u_tip = []
+time.append(0.0)
+u_tip.append(0.0)
+
 E_ext = 0
 
 displacement_out = File("Solid/FSI-S/u_fsi.pvd")
@@ -176,10 +182,35 @@ def avg(x_old, x_new, alpha):
 u_np1.rename("Displacement", "")
 displacement_out << u_n
 
-while precice.is_coupling_ongoing():
+# stress computation
+def local_project(v, V, u=None):
+    """Element-wise projection using LocalSolver"""
+    dv = TrialFunction(V)
+    v_ = TestFunction(V)
+    a_proj = inner(dv, v_)*dx
+    b_proj = inner(v, v_)*dx
+    solver = LocalSolver(a_proj, b_proj)
+    solver.factorize()
+    if u is None:
+        u = Function(V)
+        solver.solve_local_rhs(u)
+        return u
+    else:
+        solver.solve_local_rhs(u)
+        return
+
+Vsig = TensorFunctionSpace(mesh, "CG", 1)
+sig = Function(Vsig, name="sigma")
+
+xdmf_file = XDMFFile("elastodynamics-results.xdmf")
+xdmf_file.parameters["flush_output"] = True
+xdmf_file.parameters["functions_share_mesh"] = True
+xdmf_file.parameters["rewrite_function_mesh"] = False
 
-    if precice.is_action_required(precice.action_write_iteration_checkpoint()):  # write checkpoint
-        precice.store_checkpoint(u_n, t, n)
+# time loop for coupling
+
+
+while precice.is_coupling_ongoing():
 
     # read data from preCICE and get a new coupling expression
     read_data = precice.read_data()
@@ -201,29 +232,39 @@ def avg(x_old, x_new, alpha):
 
     dt = Constant(np.min([precice_dt, fenics_dt]))
 
-    # Write new displacements to preCICE
-    precice.write_data(u_np1)
-
-    # Call to advance coupling, also returns the optimum time step value
     precice_dt = precice.advance(dt(0))
 
-    # Either revert to old step if timestep has not converged or move to next timestep
-    if precice.is_action_required(precice.action_read_iteration_checkpoint()):  # roll back to checkpoint
-        u_cp, t_cp, n_cp = precice.retrieve_checkpoint()
-        u_n.assign(u_cp)
-        t = t_cp
-        n = n_cp
-    else:
-        u_n.assign(u_np1)
-        t += float(dt)
-        n += 1
+    u_n.assign(u_np1)
+    t += float(dt)
+    n += 1
 
     if precice.is_time_window_complete():
+        local_project(sigma(u_n), Vsig, sig)
+
+        # Plot tip displacement evolution
+        displacement_out << u_n
+
+        xdmf_file.write(u_n,t)
+        xdmf_file.write(sig,t)           
         update_fields(u_np1, saved_u_old, v_n, a_n)
-        if n % 10 == 0:
-            displacement_out << (u_n, t)
 
-# Plot tip displacement evolution
-displacement_out << u_n
+        if n % 10==0:
+            local_project(sigma(u_n), Vsig, sig)
+
+            displacement_out << (u_n,t)
+            xdmf_file.write(u_n,t)
+            xdmf_file.write(sig,t)
+
+        u_tip.append(u_n(0.,1.)[0])
+        time.append(t)
+
 
 precice.finalize()
+
+# Plot tip displacement evolution
+displacement_out << u_n
+plt.figure()
+plt.plot(time, u_tip)
+plt.xlabel("Time")
+plt.ylabel("Tip displacement")
+plt.show()