use tuple valued state for proper checkpointing

Author: BenjaminRodenberg

FSI/cylinderFlap/OpenFOAM-FEniCS/Solid/cyl-flap.py

@@ -10,6 +10,12 @@
 import matplotlib.pyplot as plt
 from fenicsadapter import Adapter
 from enum import Enum
+import argparse
+
+parser = argparse.ArgumentParser()
+parser.add_argument("--wr-tag", "-wr", help="wr tag", choices=["11", "52"], default="11", type=str)
+parser.add_argument("--waveform-interpolation-strategy", "-wri", help="interpolation", choices=["linear", "quadratic"], default="linear", type=str)
+args = parser.parse_args()
 
 # Beam geometry
 dim = 2  # number of dimensions
@@ -71,7 +77,8 @@ def remaining_boundary(x, on_boundary):
 
 # displacement fields
 u_np1 = Function(V)
-saved_u_old = Function(V)
+v_np1 = Function(V)
+a_np1 = Function(V)
 
 # function known from previous timestep
 u_n = Function(V)
@@ -90,11 +97,11 @@ def remaining_boundary(x, on_boundary):
 # get the adapter ready
 
 # read fenics-adapter json-config-file)
-adapter_config_filename = "../tools/precice-adapter-config-fsi-s_wr.json"
-other_adapter_config_filename = "../tools/precice-adapter-config-excite.json"
+adapter_config_filename = "../tools/precice-adapter-config-fsi-s_wr{wr_tag}.json".format(wr_tag=args.wr_tag)
+other_adapter_config_filename = "../tools/precice-adapter-config-excite_wr{wr_tag}.json".format(wr_tag=args.wr_tag)
 
 # create Adapter
-precice = Adapter(adapter_config_filename, other_adapter_config_filename)
+precice = Adapter(adapter_config_filename, other_adapter_config_filename, wr_interpolation_strategy=args.waveform_interpolation_strategy)
 
 # create subdomains used by the adapter
 clamped_boundary_domain = AutoSubDomain(left_boundary)
@@ -104,7 +111,7 @@ def remaining_boundary(x, on_boundary):
                                 mesh=mesh,
                                 read_field=f_N_function,
                                 write_field=u_function,
-                                u_n=u_n,
+                                u_n=(u_n, v_n, a_n),
                                 dimension=dim,
                                 dirichlet_boundary=clamped_boundary_domain)
 
@@ -115,8 +122,8 @@ def remaining_boundary(x, on_boundary):
 
 
 # generalized alpha method (time stepping) parameters
-alpha_m = Constant(0.2)
-alpha_f = Constant(0.4)
+alpha_m = Constant(0)
+alpha_f = Constant(0)
 gamma = Constant(0.5 + alpha_f - alpha_m)
 beta = Constant((gamma + 0.5)**2 * 0.25)
 
@@ -187,10 +194,10 @@ def avg(x_old, x_new, alpha):
 
 
 # residual
-a_np1 = update_acceleration(du, u_n, v_n, a_n, ufl=True)
-v_np1 = update_velocity(a_np1, u_n, v_n, a_n, ufl=True)
+da = update_acceleration(du, u_n, v_n, a_n, ufl=True)
+#v_np1 = update_velocity(a_np1, u_n, v_n, a_n, ufl=True)
 
-res = m(avg(a_n, a_np1, alpha_m), v) + k(avg(u_n, du, alpha_f), v)
+res = m(avg(a_n, da, alpha_m), v) + k(avg(u_n, du, alpha_f), v)
 
 Forces_x, Forces_y = precice.create_force_boundary_condition(V)
 
@@ -214,7 +221,8 @@ def avg(x_old, x_new, alpha):
 u_np1.rename("Displacement", "")
 displacement_out << u_n
 
-
+u_tip_buffer = []
+t_buffer = []
 # time loop for coupling
 while precice.is_coupling_ongoing():
     A, b = assemble_system(a_form, L_form, bc)
@@ -229,19 +237,25 @@ def avg(x_old, x_new, alpha):
     assert(b is not b_forces)
     solve(A, u_np1.vector(), b_forces)
 
-    t, n, precice_timestep_complete, precice_dt, Forces_x, Forces_y = precice.advance(u_np1, u_np1, u_n, t, float(dt), n)
+    a_np1 = project(update_acceleration(u_np1, u_n, v_n, a_n, ufl=False), V)
+    v_np1 = project(update_velocity(a_np1, u_n, v_n, a_n, ufl=False), V)
+    
+    t, n, precice_timestep_complete, precice_dt, Forces_x, Forces_y = precice.advance(u_np1, (u_np1, v_np1, a_np1), (u_n, v_n, a_n), t, float(dt), n)
     dt = Constant(np.min([precice_dt, fenics_dt]))
 
-    if precice_timestep_complete:
-        
-        update_fields(u_np1, saved_u_old, v_n, a_n)
-        
+    u_tip_buffer.append(u_n(0.6, 0.2)[1])
+    t_buffer.append(t)
+
+    if precice_timestep_complete:              
         if n % 20==0:
             displacement_out << (u_n, t)
 
-        u_tip.append(u_n(0.6, 0.2)[1])
+        u_tip += u_tip_buffer
+        u_tip_buffer = []
+        time += t_buffer
+        t_buffer = []
         print(u_n(0.6, 0.2))
-        time.append(t)
+
 
 
 # Plot tip displacement evolution
@@ -252,6 +266,8 @@ def avg(x_old, x_new, alpha):
 plt.ylabel("Tip displacement")
 plt.show()
 
-output_file = open("subiteration_out.txt", "a")
+output_file = open("wr11_out.txt", "a")
 output_file.write("{};{}\n".format(u_tip[-1], dt(0)))
 precice.finalize()
+print(np.array2string(np.array(time), separator=','))
+print(np.array2string(np.array(u_tip), separator=','))

FSI/cylinderFlap/OpenFOAM-FEniCS/tools/README.md

@@ -8,16 +8,16 @@ You will have to replace `adapter_config_filename = "precice-adapter-config-fsi-
 Run the following commands from `tutorials/FSI/cylinderFlap/OpenFOAM-FEniCS`:
 
 * `python3 Solid/cyl-flap.py`
-* `python3 tools/excitement.py tools/precice-config.xml`
+* `python3 tools/excitement.py`
 
-A file `out.txt` will be created in this folder. The result of the simulation is written to this folder
+A file `out*.txt` will be created in this folder. The result of the simulation is written to this folder
 
-## Using multirate / subiterations
+## Using multirate
 
 Depending on `precice_dt` and `fenics_dt` the computations are done using multirate timeintegration or not.
 
-* *no multirate:* the input force is provided for every single timestep
-* *multirate:* the input force is only provided for the timestep size `precice_dt`, if `fenics_dt < precice_dt` the same force will be applied for all fenics timesteps inside of a window.
+* *no multirate:* the input force is provided for every single timestep append `-wr 11` to the commands for running
+* *multirate:* the input force is only provided for some timesteps, this grid does not match the grid of the structure solver. We interpolate between force evaluations. Use `-wr 52`, additionally, you can choose the interpolation strategy via `-wri linear` or `-wri quadratic`.
 
 # Convergence Study
 

FSI/cylinderFlap/OpenFOAM-FEniCS/tools/excitement.py

@@ -4,7 +4,7 @@
 import numpy as np
 import precice_future as precice
 
-omega_excitement = 1
+omega_excitement = 4
 
 # beam geometry
 L = 0.35  # length
@@ -24,7 +24,7 @@ def compute_force(time, force_excitement):
 
 
 parser = argparse.ArgumentParser()
-parser.add_argument("configurationFileName", help="Name of the xml config file.", type=str)
+parser.add_argument("--wr-tag", "-wr", help="wr tag", choices=["11", "52"], default="11", type=str)
 
 try:
     args = parser.parse_args()
@@ -33,11 +33,9 @@ def compute_force(time, force_excitement):
     print("Usage: python ./solverdummy precice-config participant-name mesh-name")
     quit()
 
-configuration_file_name = args.configurationFileName
+configuration_file_name = "tools/precice-config_wr{wr_tag}.xml".format(wr_tag=args.wr_tag)
 participant_name = "Excitement"
 mesh_name = "Excitement-Mesh"
-write_data_name = "Forces0"
-read_data_name = "Displacements0"
 
 n_vertices = 1
 
@@ -59,8 +57,6 @@ def compute_force(time, force_excitement):
 vertices[:, 1] = y_attack
 
 vertex_ids = interface.set_mesh_vertices(mesh_id, vertices)
-write_data_id = interface.get_data_id(write_data_name, mesh_id)
-read_data_id = interface.get_data_id(read_data_name, mesh_id)
 
 dt = interface.initialize()
 t = 0
@@ -78,11 +74,20 @@ def compute_force(time, force_excitement):
         print("DUMMY: Writing iteration checkpoint")
         interface.fulfilled_action(precice.action_write_iteration_checkpoint())
 
-    forces = compute_force(t+dt, force_excitement)
-    print("sending forces: {}".format(forces))
-    interface.write_block_vector_data(write_data_id, vertex_ids, forces)
+    if args.wr_tag == "11":
+        forces = compute_force(t+dt, force_excitement)
+        print("sending forces: {}".format(forces))
+        interface.write_block_vector_data(interface.get_data_id("Forces1", mesh_id), vertex_ids, forces)
+    elif args.wr_tag == "52":
+        forces = compute_force(t+.5 * dt, force_excitement)
+        print("sending forces: {}".format(forces))
+        interface.write_block_vector_data(interface.get_data_id("Forces1", mesh_id), vertex_ids, forces)
+
+        forces = compute_force(t+dt, force_excitement)
+        print("sending forces: {}".format(forces))
+        interface.write_block_vector_data(interface.get_data_id("Forces2", mesh_id), vertex_ids, forces)
+
     dt = interface.advance(dt)
-    displacement = interface.read_block_vector_data(read_data_id, vertex_ids)
 
     if interface.is_action_required(precice.action_read_iteration_checkpoint()):
         print("DUMMY: Reading iteration checkpoint")

FSI/cylinderFlap/OpenFOAM-FEniCS/tools/plot_convergence_study.py

@@ -17,8 +17,8 @@ def outfile_to_data(filename):
         uu.append(u)
         dtdt.append(dt)
 
-    u_values = np.array(uu[1:], dtype=float)
-    dt_values = np.array(dtdt[1:], dtype=float)
+    u_values = np.array(uu, dtype=float)
+    dt_values = np.array(dtdt, dtype=float)
 
 
     def unzip(iterable):
@@ -31,7 +31,7 @@ def unzip(iterable):
     return [u for u, _ in paired], [dt for _, dt in paired]
 
 
-u_values, dt_values = outfile_to_data("../out.txt")
+u_values, dt_values = outfile_to_data("wr11_out.txt")
 
 reference_value = u_values[0]
 
@@ -40,18 +40,29 @@ def unzip(iterable):
 u_values = np.array(u_values)
 dt_values = np.array(dt_values)
 
-plt.loglog(dt_values[1:], errors[1:], '.', label='num')
+plt.loglog(dt_values, errors, '.', label='newmark')
 
-u_values, dt_values = outfile_to_data("../subiteration_out.txt")
+u_values, dt_values = outfile_to_data("wr52_lin_out.txt")
 
 errors = np.abs(u_values - reference_value)
 
 u_values = np.array(u_values)
 dt_values = np.array(dt_values)
 
-plt.loglog(dt_values, errors, '.', label='num w subiterations')
-plt.loglog(dt_values[1:], dt_values[1:], '--', label='O(h^1)')
-plt.loglog(dt_values[1:], dt_values[1:]**2, ':', label='O(h^2)')
+plt.loglog(dt_values[1:], errors[1:], '1', label='wr52_lin')
+
+u_values, dt_values = outfile_to_data("wr52_quad_out.txt")
+
+errors = np.abs(u_values - reference_value)
+
+u_values = np.array(u_values)
+dt_values = np.array(dt_values)
+
+plt.loglog(dt_values[1:], errors[1:], '2', label='wr52_quad')
+
+
+plt.loglog(dt_values[1:], errors[1]/dt_values[1]*dt_values[1:], '--', label='O(h^1)')
+plt.loglog(dt_values[1:], errors[1]/dt_values[1]**2*dt_values[1:]**2, ':', label='O(h^2)')
 plt.legend()
 plt.show()
 

FSI/cylinderFlap/OpenFOAM-FEniCS/tools/precice-adapter-config-excite_wr11.json

@@ -0,0 +1,12 @@
+{
+  "solver_name": "fenics",
+  "config_file_name": "../tools/precice-config_wr11.xml",  
+  "interface": {
+      "coupling_mesh_name": "Excitement",
+      "write_data_name": "Forces",
+      "read_data_name": "Displacements"
+    },
+  "waveform": {
+      "n_substeps": 1
+    }
+}

FSI/cylinderFlap/OpenFOAM-FEniCS/tools/precice-adapter-config-excite.json FSI/cylinderFlap/OpenFOAM-FEniCS/tools/precice-adapter-config-excite_wr52.json

@@ -1,12 +1,12 @@
 {
   "solver_name": "fenics",
-  "config_file_name": "../tools/precice-config.xml",  
+  "config_file_name": "../tools/precice-config_wr52.xml",  
   "interface": {
       "coupling_mesh_name": "Excitement",
       "write_data_name": "Forces0",
       "read_data_name": "Displacements0"
     },
   "waveform": {
-      "n_substeps": 1
+      "n_substeps": 2
     }
 }

FSI/cylinderFlap/OpenFOAM-FEniCS/tools/precice-adapter-config-fsi-s_wr.json

@@ -42,7 +42,7 @@
         <coupling-scheme:serial-explicit>
             <participants first="Excitement" second="fenics" />
             <max-time value="1"/>
-            <timestep-length value="0.1" />
+            <timestep-length value="0.001" />
             <exchange data="Forces1" mesh="Solid" from="Excitement" to="fenics"/>
             <exchange data="Displacements1" mesh="Solid" from="fenics" to="Excitement" initialize="0"/>
         </coupling-scheme:serial-explicit>

FSI/cylinderFlap/OpenFOAM-FEniCS/tools/precice-config_wr.xml FSI/cylinderFlap/OpenFOAM-FEniCS/tools/precice-config_wr11.xml

@@ -0,0 +1,81 @@
+<?xml version="1.0"?>
+
+<precice-configuration>
+
+    <log>
+      <sink type="stream" output="stdout"  filter= "%Severity% > debug" format="preCICE:%ColorizedSeverity% %Message%" enabled="true" />
+    </log>
+
+    <solver-interface dimensions="2" >
+
+        <data:vector name="Forces1"/>
+        <data:vector name="Forces2"/>
+        <data:vector name="Displacements1"/>
+        <data:vector name="Displacements2"/>
+        <data:vector name="Displacements3"/>
+        <data:vector name="Displacements4"/>
+        <data:vector name="Displacements5"/>
+
+        <mesh name="Excitement-Mesh">
+            <use-data name="Forces1"/>
+            <use-data name="Forces2"/>
+            <use-data name="Displacements1"/>
+            <use-data name="Displacements2"/>
+            <use-data name="Displacements3"/>
+            <use-data name="Displacements4"/>
+            <use-data name="Displacements5"/>
+        </mesh>
+
+        <mesh name="Solid">
+            <use-data name="Forces1"/>
+            <use-data name="Forces2"/>
+            <use-data name="Displacements1"/>
+            <use-data name="Displacements2"/>
+            <use-data name="Displacements3"/>
+            <use-data name="Displacements4"/>
+            <use-data name="Displacements5"/>
+        </mesh>
+
+        <participant name="Excitement">
+            <use-mesh name="Excitement-Mesh" provide="yes"/>
+            <use-mesh name="Solid" from="fenics"/>
+            <write-data name="Forces1" mesh="Excitement-Mesh"/>
+            <write-data name="Forces2" mesh="Excitement-Mesh"/>
+            <read-data name="Displacements1" mesh="Excitement-Mesh"/>
+            <read-data name="Displacements2" mesh="Excitement-Mesh"/>
+            <read-data name="Displacements3" mesh="Excitement-Mesh"/>
+            <read-data name="Displacements4" mesh="Excitement-Mesh"/>
+            <read-data name="Displacements5" mesh="Excitement-Mesh"/>
+            <mapping:petrbf-thin-plate-splines direction="write" from="Excitement-Mesh" to="Solid" constraint="conservative" />
+            <mapping:petrbf-thin-plate-splines direction="read" from="Solid" to="Excitement-Mesh" constraint="consistent" />
+        </participant>
+
+        <participant name="fenics">
+            <use-mesh name="Solid" provide="yes"/>
+            <read-data name="Forces1" mesh="Solid"/>
+            <read-data name="Forces2" mesh="Solid"/>
+            <write-data name="Displacements1" mesh="Solid"/>
+            <write-data name="Displacements2" mesh="Solid"/>
+            <write-data name="Displacements3" mesh="Solid"/>
+            <write-data name="Displacements4" mesh="Solid"/>
+            <write-data name="Displacements5" mesh="Solid"/>
+            <watch-point mesh="Solid" name="point1" coordinate="0.6;0.2;0." />
+        </participant>
+
+        <m2n:sockets from="Excitement" to="fenics" distribution-type="gather-scatter"/>
+
+        <coupling-scheme:serial-explicit>
+            <participants first="Excitement" second="fenics" />
+            <max-time value="1"/>
+            <timestep-length value="0.5" />
+            <exchange data="Forces1" mesh="Solid" from="Excitement" to="fenics"/>
+            <exchange data="Forces2" mesh="Solid" from="Excitement" to="fenics"/>
+            <exchange data="Displacements1" mesh="Solid" from="fenics" to="Excitement" initialize="0"/>
+            <exchange data="Displacements2" mesh="Solid" from="fenics" to="Excitement" initialize="0"/>
+            <exchange data="Displacements3" mesh="Solid" from="fenics" to="Excitement" initialize="0"/>
+            <exchange data="Displacements4" mesh="Solid" from="fenics" to="Excitement" initialize="0"/>
+            <exchange data="Displacements5" mesh="Solid" from="fenics" to="Excitement" initialize="0"/>
+        </coupling-scheme:serial-explicit>
+    </solver-interface>
+
+</precice-configuration>