diff --git a/perpendicular-flap/.gitignore b/perpendicular-flap/.gitignore new file mode 100644 index 000000000..493b6986c --- /dev/null +++ b/perpendicular-flap/.gitignore @@ -0,0 +1,5 @@ +preCICE-output/ +Solid-fenics* +*.log +vtk* +Solid/FSI-S diff --git a/perpendicular-flap/fluid-aste/README.md b/perpendicular-flap/fluid-aste/README.md new file mode 100644 index 000000000..d2b3c5358 --- /dev/null +++ b/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 `` to the `precice-config.xml` in `tutorials/perpendicular-flap` in the Solid participant. The result should look the following way: + +```xml + + + + + + +``` + +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 + + + + + ``` + +#### 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. diff --git a/perpendicular-flap/precice-config-aste.xml b/perpendicular-flap/precice-config-aste.xml new file mode 100644 index 000000000..68048db4c --- /dev/null +++ b/perpendicular-flap/precice-config-aste.xml @@ -0,0 +1,50 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/perpendicular-flap/precice-config.xml b/perpendicular-flap/precice-config.xml index 97ebd5afc..9469529cd 100644 --- a/perpendicular-flap/precice-config.xml +++ b/perpendicular-flap/precice-config.xml @@ -42,10 +42,19 @@ - + + diff --git a/perpendicular-flap/solid-fenics/precice-adapter-config-fsi-s-aste.json b/perpendicular-flap/solid-fenics/precice-adapter-config-fsi-s-aste.json new file mode 100644 index 000000000..1c478caad --- /dev/null +++ b/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" + } +} diff --git a/perpendicular-flap/solid-fenics/solid.py b/perpendicular-flap/solid-fenics/solid.py index 66bcb3a62..1c89d7421 100644 --- a/perpendicular-flap/solid-fenics/solid.py +++ b/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()