Update all FEniCS participants to be compatible with precice:develop (v3) (precice#333)

* Update all FEniCS participants to be compatible with precice:develop and fenics-adapter updates
* Change from thin plate splines mapping to general rbf mapping
* Update how dt is used.
* Store time in vtk files for better visualization.
* Partitioned heat conduction:
  * Fix config and case definition. (precice/precice#1610)
  * Use stricter error tolerance.

Author: IshaanDesai

channel-transport-reaction/chemical-fenics/chemical-reaction-advection-diffusion.py

@@ -27,9 +27,8 @@ def inside(self, x, on_boundary):
 # Initialize preCICE
 precice = fenicsprecice.Adapter(
     adapter_config_filename="chemical-reaction-advection-diffusion.json")
-precice_dt = precice.initialize(
-    coupling_subdomain=CouplingDomain(),
-    read_function_space=W)
+precice.initialize(coupling_subdomain=CouplingDomain(), read_function_space=W)
+precice_dt = precice.get_max_time_step_size()
 
 flow_expr = precice.create_coupling_expression()
 
@@ -95,13 +94,14 @@ def inside(self, x, on_boundary):
   # No implicit coupling
 while precice.is_coupling_ongoing():
 
-    read_data = precice.read_data()
+    precice_dt = precice.get_max_time_step_size()
+    dt = np.min([default_dt, precice_dt])
+    read_data = precice.read_data(dt)
     precice.update_coupling_expression(flow_expr, read_data)
     flow_old.assign(flow)
     flow.interpolate(flow_expr)
     # If we add writing, do it here
 
-    dt = np.min([default_dt, precice_dt])
     k.assign(1. / dt)
 
     t += dt
@@ -128,6 +128,6 @@ def inside(self, x, on_boundary):
     vtkfileC << u_C, t
     vtkfileFlow << flow, t
 
-    precice_dt = precice.advance(dt)
+    precice.advance(dt)
 
 precice.finalize()

channel-transport-reaction/fluid-fenics/fluid.py

@@ -118,10 +118,8 @@ def inside(self, x, on_boundary):
 
 
 precice = fenicsprecice.Adapter(adapter_config_filename="fluid-config.json")
-precice_dt = precice.initialize(
-    coupling_subdomain=CouplingDomain(),
-    write_object=u_)
-
+precice.initialize(coupling_subdomain=CouplingDomain(), write_object=u_)
+precice_dt = precice.get_max_time_step_size()
 
 dt = np.min([default_dt, precice_dt])
 k.assign(dt)
@@ -155,7 +153,8 @@ def inside(self, x, on_boundary):
 
     vtkfile << u_
 
-    precice_dt = precice.advance(dt)
+    precice.advance(dt)
+    precice_dt = precice.get_max_time_step_size()
     dt = np.min([default_dt, precice_dt])
     k.assign(dt)
 

elastic-tube-3d/solid-fenics/solid.py

@@ -65,7 +65,8 @@ def neumann_boundary(x, on_boundary):
 precice = Adapter(adapter_config_filename="precice-adapter-config-fsi-s.json")
 
 # Initialize the coupling interface
-precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V, fixed_boundary=fixed_boundary)
+precice.initialize(coupling_boundary, read_function_space=V, write_object=V, fixed_boundary=fixed_boundary)
+precice_dt = precice.get_max_time_step_size()
 
 fenics_dt = precice_dt  # if fenics_dt == precice_dt, no subcycling is applied
 dt = Constant(np.min([precice_dt, fenics_dt]))
@@ -168,15 +169,18 @@ def avg(x_old, x_new, alpha):
 
 u_n.rename("Displacement", "")
 u_np1.rename("Displacement", "")
-displacement_out << u_n
+displacement_out << (u_n, t)
 
 while precice.is_coupling_ongoing():
 
-    if precice.is_action_required(precice.action_write_iteration_checkpoint()):  # write checkpoint
+    if precice.requires_writing_checkpoint():  # write checkpoint
         precice.store_checkpoint(u_n, t, n)
 
+    precice_dt = precice.get_max_time_step_size()
+    dt = Constant(np.min([precice_dt, fenics_dt]))
+
     # read data from preCICE and get a new coupling expression
-    read_data = precice.read_data()
+    read_data = precice.read_data(dt)
 
     # Update the point sources on the coupling boundary with the new read data
     forces_x, forces_y, forces_z = precice.get_point_sources(read_data)
@@ -195,17 +199,15 @@ def avg(x_old, x_new, alpha):
     assert (b is not b_forces)
     solve(A, u_np1.vector(), b_forces)
 
-    dt = Constant(np.min([precice_dt, fenics_dt]))
-
     # Write relative displacements to preCICE
     u_delta.vector()[:] = u_np1.vector()[:] - u_n.vector()[:]
     precice.write_data(u_delta)
 
     # Call to advance coupling, also returns the optimum time step value
-    precice_dt = precice.advance(dt(0))
+    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
+    if precice.requires_reading_checkpoint():  # roll back to checkpoint
         u_cp, t_cp, n_cp = precice.retrieve_checkpoint()
         u_n.assign(u_cp)
         t = t_cp
@@ -221,6 +223,6 @@ def avg(x_old, x_new, alpha):
             displacement_out << (u_n, t)
 
 # Plot tip displacement evolution
-displacement_out << u_n
+displacement_out << (u_n, t)
 
 precice.finalize()

flow-over-heated-plate/solid-fenics/solid.py

@@ -115,13 +115,14 @@ def determine_heat_flux(V_g, u, k, flux):
 # Adapter definition and initialization
 precice = Adapter(adapter_config_filename="precice-adapter-config.json")
 
-precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V_flux_y)
+precice.initialize(coupling_boundary, read_function_space=V, write_object=V_flux_y)
 
 # Create a FEniCS Expression to define and control the coupling boundary values
 coupling_expression = precice.create_coupling_expression()
 
 # Assigning appropriate dt
 dt = Constant(0)
+precice_dt = precice.get_max_time_step_size()
 dt.assign(np.min([fenics_dt, precice_dt]))
 
 # Define variational problem
@@ -149,7 +150,7 @@ def determine_heat_flux(V_g, u, k, flux):
 
 # Create output file
 file_out = File("output/%s.pvd" % precice.get_participant_name())
-file_out << u_n
+file_out << (u_n, t)
 
 print("output vtk for time = {}".format(float(t)))
 n = 0
@@ -159,16 +160,16 @@ def determine_heat_flux(V_g, u, k, flux):
 
 while precice.is_coupling_ongoing():
 
-    if precice.is_action_required(precice.action_write_iteration_checkpoint()):  # write checkpoint
+    if precice.requires_writing_checkpoint():  # write checkpoint
         precice.store_checkpoint(u_n, t, n)
 
-    read_data = precice.read_data()
+    precice_dt = precice.get_max_time_step_size()
+    dt.assign(np.min([fenics_dt, precice_dt]))
+    read_data = precice.read_data(dt)
 
     # Update the coupling expression with the new read data
     precice.update_coupling_expression(coupling_expression, read_data)
 
-    dt.assign(np.min([fenics_dt, precice_dt]))
-
     # Compute solution
     solve(a == L, u_np1, bcs)
 
@@ -177,9 +178,9 @@ def determine_heat_flux(V_g, u, k, flux):
     fluxes_y = fluxes.sub(1)  # only exchange y component of flux.
     precice.write_data(fluxes_y)
 
-    precice_dt = precice.advance(dt(0))
+    precice.advance(dt(0))
 
-    if precice.is_action_required(precice.action_read_iteration_checkpoint()):  # roll back to checkpoint
+    if precice.requires_reading_checkpoint():  # roll back to checkpoint
         u_cp, t_cp, n_cp = precice.retrieve_checkpoint()
         u_n.assign(u_cp)
         t = t_cp
@@ -193,7 +194,7 @@ def determine_heat_flux(V_g, u, k, flux):
         tol = 10e-5  # we need some tolerance, since otherwise output might be skipped.
         if abs((t + tol) % dt_out) < 2 * tol:  # output if t is a multiple of dt_out
             print("output vtk for time = {}".format(float(t)))
-            file_out << u_n
+            file_out << (u_n, t)
 
     # Update dirichlet BC
     u_D.t = t + float(dt)

partitioned-heat-conduction-complex/fenics/heat.py

@@ -103,14 +103,15 @@ def determine_gradient(V_g, u, flux):
 # Initialize the adapter according to the specific participant
 if problem is ProblemType.DIRICHLET:
     precice = Adapter(adapter_config_filename="precice-adapter-config-D.json")
-    precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=f_N_function)
+    precice.initialize(coupling_boundary, read_function_space=V, write_object=f_N_function)
 elif problem is ProblemType.NEUMANN:
     precice = Adapter(adapter_config_filename="precice-adapter-config-N.json")
-    precice_dt = precice.initialize(coupling_boundary, read_function_space=V_g, write_object=u_D_function)
+    precice.initialize(coupling_boundary, read_function_space=V_g, write_object=u_D_function)
 
 boundary_marker = False
 
 dt = Constant(0)
+precice_dt = precice.get_max_time_step_size()
 dt.assign(np.min([fenics_dt, precice_dt]))
 
 # Define variational problem
@@ -173,7 +174,7 @@ def determine_gradient(V_g, u, flux):
 # output solution and reference solution at t=0, n=0
 n = 0
 print('output u^%d and u_ref^%d' % (n, n))
-temperature_out << u_n
+temperature_out << (u_n, t)
 ref_out << u_ref
 ranks << mesh_rank
 
@@ -191,10 +192,12 @@ def determine_gradient(V_g, u, flux):
 while precice.is_coupling_ongoing():
 
     # write checkpoint
-    if precice.is_action_required(precice.action_write_iteration_checkpoint()):
+    if precice.requires_writing_checkpoint():
         precice.store_checkpoint(u_n, t, n)
 
-    read_data = precice.read_data()
+    precice_dt = precice.get_max_time_step_size()
+    dt.assign(np.min([fenics_dt, precice_dt]))
+    read_data = precice.read_data(dt)
     if problem is ProblemType.DIRICHLET and (domain_part is DomainPart.CIRCULAR or domain_part is DomainPart.RECTANGLE):
         # We have to data for an arbitrary point that is not on the circle, to obtain exact solution.
         # See https://github.com/precice/fenics-adapter/issues/113 for details.
@@ -203,8 +206,6 @@ def determine_gradient(V_g, u, flux):
     # Update the coupling expression with the new read data
     precice.update_coupling_expression(coupling_expression, read_data)
 
-    dt.assign(np.min([fenics_dt, precice_dt]))
-
     # Compute solution u^n+1, use bcs u_D^n+1, u^n and coupling bcs
     solve(a == L, u_np1, bcs)
 
@@ -217,10 +218,10 @@ def determine_gradient(V_g, u, flux):
         # Neumann problem reads flux and writes temperature on boundary to Dirichlet problem
         precice.write_data(u_np1)
 
-    precice_dt = precice.advance(dt(0))
+    precice.advance(dt(0))
 
     # roll back to checkpoint
-    if precice.is_action_required(precice.action_read_iteration_checkpoint()):
+    if precice.requires_reading_checkpoint():
         u_cp, t_cp, n_cp = precice.retrieve_checkpoint()
         u_n.assign(u_cp)
         t = t_cp
@@ -237,7 +238,7 @@ def determine_gradient(V_g, u, flux):
         print('n = %d, t = %.2f: L2 error on domain = %.3g' % (n, t, error))
         # output solution and reference solution at t_n+1
         print('output u^%d and u_ref^%d' % (n, n))
-        temperature_out << u_n
+        temperature_out << (u_n, t)
         ref_out << u_ref
         error_out << error_pointwise
 

partitioned-heat-conduction-complex/fenics/precice-adapter-config-D.json

@@ -3,7 +3,7 @@
   "config_file_name": "../precice-config.xml",
   "interface": {
     "coupling_mesh_name": "Dirichlet-Mesh",
-    "write_data_name": "Flux",
+    "write_data_name": "Heat-Flux",
     "read_data_name": "Temperature"
   }
 }

partitioned-heat-conduction-complex/fenics/precice-adapter-config-N.json

@@ -4,6 +4,6 @@
   "interface": {
     "coupling_mesh_name": "Neumann-Mesh",
     "write_data_name": "Temperature",
-    "read_data_name": "Flux"
+    "read_data_name": "Heat-Flux"
   }
 }

partitioned-heat-conduction-complex/precice-config.xml

@@ -6,22 +6,22 @@
   </log>
 
   <data:scalar name="Temperature"/>
-  <data:vector name="Flux"/>
+  <data:vector name="Heat-Flux"/>
 
   <mesh name="Dirichlet-Mesh" dimensions="2">
     <use-data name="Temperature"/>
-    <use-data name="Flux"/>
+    <use-data name="Heat-Flux"/>
   </mesh>
 
   <mesh name="Neumann-Mesh" dimensions="2">
     <use-data name="Temperature"/>
-    <use-data name="Flux"/>
+    <use-data name="Heat-Flux"/>
   </mesh>
 
   <participant name="Dirichlet">
     <provide-mesh name="Dirichlet-Mesh" />
     <receive-mesh name="Neumann-Mesh" from="Neumann" />
-    <write-data name="Flux" mesh="Dirichlet-Mesh" />
+    <write-data name="Heat-Flux" mesh="Dirichlet-Mesh" />
     <read-data name="Temperature" mesh="Dirichlet-Mesh" />
     <mapping:nearest-projection direction="read" from="Neumann-Mesh" to="Dirichlet-Mesh" constraint="consistent" />
   </participant>
@@ -30,7 +30,7 @@
     <provide-mesh name="Neumann-Mesh" />
     <receive-mesh name="Dirichlet-Mesh" from="Dirichlet"/>
     <write-data name="Temperature" mesh="Neumann-Mesh"/>
-    <read-data name="Flux" mesh="Neumann-Mesh"/>
+    <read-data name="Heat-Flux" mesh="Neumann-Mesh"/>
     <mapping:nearest-projection direction="read" from="Dirichlet-Mesh" to="Neumann-Mesh" constraint="consistent" />
   </participant>
 
@@ -41,9 +41,9 @@
     <max-time value="1.0"/>
     <time-window-size value="0.1"/>
     <max-iterations value="100"/>
-    <exchange data="Flux" mesh="Dirichlet-Mesh" from="Dirichlet" to="Neumann" />
+    <exchange data="Heat-Flux" mesh="Dirichlet-Mesh" from="Dirichlet" to="Neumann" />
     <exchange data="Temperature" mesh="Neumann-Mesh" from="Neumann" to="Dirichlet" initialize="true"/>
-    <relative-convergence-measure data="Flux" mesh="Dirichlet-Mesh" limit="1e-5"/>
+    <relative-convergence-measure data="Heat-Flux" mesh="Dirichlet-Mesh" limit="1e-5"/>
     <relative-convergence-measure data="Temperature" mesh="Neumann-Mesh" limit="1e-5"/>
     <acceleration:IQN-ILS>
       <data name="Temperature" mesh="Neumann-Mesh"/>