Modify FEniCS tutorials to enable parallel runs (#120)

* Modifying all FEniCS based tutorials to be compatible with parallel design of FEniCS-Adapter
Co-authored-by: BenjaminRueth <benjamin.rueth@tum.de>

Author: IshaanDesai

CHT/flow-over-plate/buoyantPimpleFoam-fenics/Solid/heat.py

@@ -6,7 +6,7 @@
 from fenics import Function, SubDomain, RectangleMesh, BoxMesh, FunctionSpace, VectorFunctionSpace, Point, \
     Expression, Constant, DirichletBC, \
     TrialFunction, TestFunction, File, solve, plot, lhs, rhs, grad, inner, dot, dx, ds, interpolate, project, \
-    near
+    near, MeshFunction, MPI
 from fenicsprecice import Adapter
 import numpy as np
 
@@ -140,8 +140,17 @@ def determine_heat_flux(V_g, u, k, flux):
 t = 0
 u_D.t = t + dt
 
+# mark mesh w.r.t ranks
+ranks = File("Solid/VTK/ranks%s.pvd.pvd" % precice.get_participant_name())
+mesh_rank = MeshFunction("size_t", mesh, mesh.topology().dim())
+mesh_rank.set_all(MPI.rank(MPI.comm_world))
+mesh_rank.rename("myRank", "")
+ranks << mesh_rank
+
+# Create output file
 file_out = File("Solid/VTK/%s.pvd" % precice.get_participant_name())
 file_out << u_n
+
 print("output vtk for time = {}".format(float(t)))
 n = 0
 

CHT/flow-over-plate/buoyantPimpleFoam-fenics/precice-config.xml

@@ -6,7 +6,7 @@
         <sink filter="%Severity% > debug and %Rank% = 0" format="---[precice] %ColorizedSeverity% %Message%" enabled="true"/>
     </log>
 
-    <solver-interface dimensions="3">
+    <solver-interface dimensions="2">
 
     <data:scalar name="Temperature"/>
     <data:scalar name="Heat-Flux"/>
@@ -23,32 +23,32 @@
 
     <participant name="Fluid">
         <use-mesh name="Fluid-Mesh" provide="yes"/>
+        <use-mesh name="Solid-Mesh" from="Solid"/>
         <read-data name="Heat-Flux" mesh="Fluid-Mesh"/>
         <write-data name="Temperature" mesh="Fluid-Mesh"/>
+        <mapping:nearest-neighbor direction="read" from="Solid-Mesh" to="Fluid-Mesh" constraint="consistent"/>
     </participant>
 
     <participant name="Solid">
         <use-mesh name="Fluid-Mesh" from="Fluid"/>
         <use-mesh name="Solid-Mesh" provide="yes"/>
         <read-data name="Temperature" mesh="Solid-Mesh"/>
         <write-data name="Heat-Flux" mesh="Solid-Mesh"/>
-        <mapping:nearest-neighbor direction="read" from="Fluid-Mesh" to="Solid-Mesh" constraint="consistent" timing="initial"/>
-        <mapping:nearest-neighbor direction="write" from="Solid-Mesh" to="Fluid-Mesh" constraint="consistent" timing="initial"/>
+        <mapping:nearest-neighbor direction="read" from="Fluid-Mesh" to="Solid-Mesh" constraint="consistent"/>
     </participant>
 
     <m2n:sockets from="Fluid" to="Solid"/>
 
     <coupling-scheme:serial-implicit>
         <time-window-size value="0.01"/>
         <max-time value="1"/>
-        <max-iterations value="200"/>
         <participants first="Fluid" second="Solid"/>
         <exchange data="Temperature" mesh="Fluid-Mesh" from="Fluid" to="Solid"/>
-        <exchange data="Heat-Flux" mesh="Fluid-Mesh" from="Solid" to="Fluid"/>
+        <exchange data="Heat-Flux" mesh="Solid-Mesh" from="Solid" to="Fluid"/>
+        <max-iterations value="200"/>
         <relative-convergence-measure limit="1.0e-6" data="Temperature" mesh="Fluid-Mesh"/>
-        <extrapolation-order value="0"/>
         <acceleration:IQN-ILS>
-            <data mesh="Fluid-Mesh" name="Heat-Flux" />
+            <data mesh="Solid-Mesh" name="Heat-Flux" />
             <initial-relaxation value="0.01" />
             <max-used-iterations value="80" />
             <time-windows-reused value="10" />

FSI/cylinderFlap/OpenFOAM-FEniCS/precice-config.xml

@@ -6,7 +6,7 @@
         <sink filter="%Severity% > debug and %Rank% = 0" format="---[precice] %ColorizedSeverity% %Message%" enabled="true"/>
     </log>
 
-    <solver-interface dimensions="3">
+    <solver-interface dimensions="2">
 
     <data:vector name="Forces0"/>
     <data:vector name="Displacements0"/>
@@ -30,15 +30,15 @@
         <use-mesh name="Solid" from="fenics"/>
         <write-data name="Forces0" mesh="Fluid-Mesh-Faces"/>
         <read-data name="Displacements0" mesh="Fluid-Mesh-Nodes"/>
-        <mapping:rbf-thin-plate-splines direction="write" from="Fluid-Mesh-Faces" to="Solid" constraint="conservative" z-dead="true"/>
-        <mapping:rbf-thin-plate-splines direction="read" from="Solid" to="Fluid-Mesh-Nodes" constraint="consistent" z-dead="true"/>
+        <mapping:rbf-thin-plate-splines direction="write" from="Fluid-Mesh-Faces" to="Solid" constraint="conservative"/>
+        <mapping:rbf-thin-plate-splines direction="read" from="Solid" to="Fluid-Mesh-Nodes" constraint="consistent"/>
     </participant>
 
     <participant name="fenics">
         <use-mesh name="Solid" provide="yes"/>
         <read-data name="Forces0" mesh="Solid"/>
         <write-data name="Displacements0" mesh="Solid"/>
-        <watch-point mesh="Solid" name="point1" coordinate="0.6;0.2;0." />
+        <watch-point mesh="Solid" name="point1" coordinate="0.6;0." />
     </participant>
 
     <m2n:sockets from="Fluid" to="fenics"/>

FSI/flap_perp/OpenFOAM-FEniCS/Solid/perp-flap.py

@@ -1,7 +1,7 @@
 # Import required libs
 from fenics import Constant, Function, AutoSubDomain, RectangleMesh, VectorFunctionSpace, interpolate, \
-    TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, \
-    Identity, inner, dx, ds, sym, grad, lhs, rhs, dot, File, solve, PointSource, assemble_system
+    TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, SubDomain, \
+    Identity, inner, dx, ds, sym, grad, lhs, rhs, dot, File, solve, PointSource, assemble_system, MPI, MeshFunction
 import dolfin
 
 from ufl import nabla_div
@@ -11,17 +11,22 @@
 from enum import Enum
 
 
-# define the two kinds of boundary: clamped and coupling Neumann Boundary
-def clamped_boundary(x, on_boundary):
-    return on_boundary and abs(x[1]) < tol
+class clampedBoundary(SubDomain):
+    def inside(self, x, on_boundary):
+        tol = 1E-14
+        if on_boundary and abs(x[1]) < tol:
+            return True
+        else:
+            return False
 
 
-def neumann_boundary(x, on_boundary):
-    """
-    determines whether a node is on the coupling boundary
-    
-    """
-    return on_boundary and ((abs(x[1] - 1) < tol) or abs(abs(x[0]) - W / 2) < tol)
+class neumannBoundary(SubDomain):
+    def inside(self, x, on_boundary):
+        tol = 1E-14
+        if on_boundary and ((abs(x[1] - 1) < tol) or abs(abs(x[0]) - W / 2) < tol):
+            return True
+        else:
+            return False
 
 
 # Geometry and material properties
@@ -46,9 +51,6 @@ def neumann_boundary(x, on_boundary):
 # create Function Space
 V = VectorFunctionSpace(mesh, 'P', 2)
 
-# BCs
-tol = 1E-14
-
 # Trial and Test Functions
 du = TrialFunction(V)
 v = TestFunction(V)
@@ -64,24 +66,19 @@ def neumann_boundary(x, on_boundary):
 f_N_function = interpolate(Expression(("1", "0"), degree=1), V)
 u_function = interpolate(Expression(("0", "0"), degree=1), V)
 
-coupling_boundary = AutoSubDomain(neumann_boundary)
-
 precice = Adapter(adapter_config_filename="precice-adapter-config-fsi-s.json")
 
-clamped_boundary_domain = AutoSubDomain(clamped_boundary)
-force_boundary = AutoSubDomain(neumann_boundary)
-
 # Initialize the coupling interface
 # Function space V is passed twice as both read and write functions are defined using the same space
-precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V,
-                                fixed_boundary=clamped_boundary_domain)
+precice_dt = precice.initialize(neumannBoundary(), read_function_space=V, write_object=V,
+                                fixed_boundary=clampedBoundary())
 
 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]))
 
 # clamp the beam at the bottom
-bc = DirichletBC(V, Constant((0, 0)), clamped_boundary)
+bc = DirichletBC(V, Constant((0, 0)), clampedBoundary())
 
 # alpha method parameters
 alpha_m = Constant(0.2)
@@ -179,11 +176,18 @@ def avg(x_old, x_new, alpha):
 u_tip.append(0.0)
 E_ext = 0
 
+# mark mesh w.r.t ranks
+mesh_rank = MeshFunction("size_t", mesh, mesh.topology().dim())
+mesh_rank.set_all(MPI.rank(MPI.comm_world) + 0)
+mesh_rank.rename("myRank", "")
+
 displacement_out = File("Solid/FSI-S/u_fsi.pvd")
+ranks = File("Solid/FSI-S/ranks%s.pvd" % precice.get_participant_name())
 
 u_n.rename("Displacement", "")
 u_np1.rename("Displacement", "")
 displacement_out << u_n
+ranks << mesh_rank
 
 while precice.is_coupling_ongoing():
 

FSI/flap_perp/OpenFOAM-FEniCS/precice-config.xml

@@ -6,7 +6,7 @@
         <sink filter="%Severity% > debug and %Rank% = 0" format="---[precice] %ColorizedSeverity% %Message%" enabled="true"/>
     </log>
 
-    <solver-interface dimensions="3">
+    <solver-interface dimensions="2">
 
     <data:vector name="Forces0"/>
     <data:vector name="Displacements0"/>
@@ -30,15 +30,15 @@
         <use-mesh name="Solid" from="fenics"/>
         <write-data name="Forces0" mesh="Fluid-Mesh-Faces"/>
         <read-data name="Displacements0" mesh="Fluid-Mesh-Nodes"/>
-        <mapping:rbf-thin-plate-splines direction="write" from="Fluid-Mesh-Faces" to="Solid" constraint="conservative" z-dead="true"/>
-        <mapping:rbf-thin-plate-splines direction="read" from="Solid" to="Fluid-Mesh-Nodes" constraint="consistent" z-dead="true"/>
+        <mapping:rbf-thin-plate-splines direction="write" from="Fluid-Mesh-Faces" to="Solid" constraint="conservative"/>
+        <mapping:rbf-thin-plate-splines direction="read" from="Solid" to="Fluid-Mesh-Nodes" constraint="consistent"/>
     </participant>
 
     <participant name="fenics">
         <use-mesh name="Solid" provide="yes"/>
         <read-data name="Forces0" mesh="Solid"/>
         <write-data name="Displacements0" mesh="Solid"/>
-        <watch-point mesh="Solid" name="point1" coordinate="-0.05;0;1" />
+        <watch-point mesh="Solid" name="point1" coordinate="-0.05;1" />
     </participant>
 
     <m2n:sockets from="Fluid" to="fenics"/>

HT/partitioned-heat/fenics-fenics/heat.py

@@ -26,7 +26,7 @@
 
 from __future__ import print_function, division
 from fenics import Function, FunctionSpace, Expression, Constant, DirichletBC, TrialFunction, TestFunction, \
-    File, solve, lhs, rhs, grad, inner, dot, dx, ds, interpolate, VectorFunctionSpace
+    File, solve, lhs, rhs, grad, inner, dot, dx, ds, interpolate, VectorFunctionSpace, MeshFunction, MPI
 from fenicsprecice import Adapter
 from errorcomputation import compute_errors
 from my_enums import ProblemType, Subcycling
@@ -174,16 +174,26 @@ def determine_gradient(V_g, u, flux):
 u_ref = interpolate(u_D, V)
 u_ref.rename("reference", " ")
 
+# mark mesh w.r.t ranks
+mesh_rank = MeshFunction("size_t", mesh, mesh.topology().dim())
+if problem is ProblemType.NEUMANN:
+    mesh_rank.set_all(MPI.rank(MPI.comm_world) + 4)
+else:
+    mesh_rank.set_all(MPI.rank(MPI.comm_world) + 0)
+mesh_rank.rename("myRank", "")
+
 # Generating output files
 temperature_out = File("out/%s.pvd" % precice.get_participant_name())
 ref_out = File("out/ref%s.pvd" % precice.get_participant_name())
 error_out = File("out/error%s.pvd" % precice.get_participant_name())
+ranks = File("out/ranks%s.pvd" % precice.get_participant_name())
 
 # 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
 ref_out << u_ref
+ranks << mesh_rank
 
 error_total, error_pointwise = compute_errors(u_n, u_ref, V)
 error_out << error_pointwise

HT/partitioned-heat/fenics-fenics/precice-config.xml

@@ -26,14 +26,15 @@
       <use-mesh name="NeumannNodes" from="HeatNeumann"/>
       <write-data name="Flux" mesh="DirichletNodes"/>
       <read-data  name="Temperature" mesh="DirichletNodes"/>
-      <mapping:nearest-projection direction="write" from="DirichletNodes" to="NeumannNodes"   constraint="consistent" timing="initial"/>
       <mapping:nearest-projection direction="read"  from="NeumannNodes"   to="DirichletNodes" constraint="consistent" timing="initial"/>
     </participant>
 
     <participant name="HeatNeumann">
-      <use-mesh name="NeumannNodes" provide="yes"/>
+	    <use-mesh name="NeumannNodes" provide="yes"/>
+	    <use-mesh name="DirichletNodes" from="HeatDirichlet"/>
       <write-data name="Temperature" mesh="NeumannNodes"/>
       <read-data  name="Flux"        mesh="NeumannNodes"/>
+      <mapping:nearest-projection direction="read" from="DirichletNodes" to="NeumannNodes"   constraint="consistent" timing="initial"/>
     </participant>
 
     <m2n:sockets from="HeatDirichlet" to="HeatNeumann"/>
@@ -43,9 +44,9 @@
          <max-time value="1"/>
          <time-window-size value=".1" valid-digits="8"/>
          <max-iterations value="100"/>
-         <exchange data="Flux"        mesh="NeumannNodes" from="HeatDirichlet" to="HeatNeumann" />
+         <exchange data="Flux"        mesh="DirichletNodes" from="HeatDirichlet" to="HeatNeumann" />
          <exchange data="Temperature" mesh="NeumannNodes" from="HeatNeumann"   to="HeatDirichlet" initialize="true"/>
-         <relative-convergence-measure data="Flux"        mesh="NeumannNodes" limit="1e-5"/>
+         <relative-convergence-measure data="Flux"        mesh="DirichletNodes" limit="1e-5"/>
          <relative-convergence-measure data="Temperature" mesh="NeumannNodes" limit="1e-5"/>
          <acceleration:IQN-ILS> 
             <data name="Temperature" mesh="NeumannNodes"/>