Merge pull request precice#124 from precice/fenics-adapter-v1.0.0

Ensure compatibility of tutorials with fenics-adapter v1.0.0

Author: BenjaminRodenberg

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

@@ -3,10 +3,11 @@
 """
 
 from __future__ import print_function, division
-from fenics import Function, SubDomain, RectangleMesh, BoxMesh, FunctionSpace, Point, Expression, Constant, DirichletBC, \
-    TrialFunction, TestFunction, File, solve, plot, lhs, rhs, grad, inner, dot, dx, ds, assemble, interpolate, project, \
-    near
-from fenicsadapter import Adapter
+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, MeshFunction, MPI
+from fenicsprecice import Adapter
 import numpy as np
 
 
@@ -59,27 +60,21 @@ def inside(self, x, on_boundary):
             return False
 
 
-def fluxes_from_temperature_full_domain(f, v_vec, k):
-    """Computes flux from weak form (see p.3 in Toselli, Andrea, and Olof
-    Widlund. Domain decomposition methods-algorithms and theory. Vol. 34.
-    Springer Science & Business Media, 2006.).
-
-    :param f: weak form with known u^{n+1}
-    :param v_vec: vector function space
+def determine_heat_flux(V_g, u, k, flux):
+    """
+    compute flux following http://hplgit.github.io/INF5620/doc/pub/fenics_tutorial1.1/tu2.html#tut-poisson-gradu
+    :param V_g: Vector function space
+    :param u: solution where gradient is to be determined
     :param k: thermal conductivity
-    :return: fluxes function
+    :param flux: returns calculated flux into this value
     """
-    fluxes_vector = assemble(f)  # assemble weak form -> evaluate integral
-    v = TestFunction(v_vec)
-    fluxes = Function(v_vec)  # create function for flux
-    area = assemble(v * ds).get_local()
-    for i in range(area.shape[0]):
-        if area[i] != 0:  # put weight from assemble on function
-            fluxes.vector()[i] = - k * fluxes_vector[i] / area[i]  # scale by surface area
-        else:
-            assert (abs(fluxes_vector[i]) < 1E-9)  # for non surface parts, we expect zero flux
-            fluxes.vector()[i] = - k * fluxes_vector[i]
-    return fluxes
+
+    w = TrialFunction(V_g)
+    v = TestFunction(V_g)
+
+    a = inner(w, v) * dx
+    L = inner(-k * grad(u), v) * dx
+    solve(a == L, flux)
 
 
 # Create mesh and define function space
@@ -88,7 +83,7 @@ def fluxes_from_temperature_full_domain(f, v_vec, k):
 nz = 1
 
 fenics_dt = 0.01  # time step size
-dt_out = 0.2
+dt_out = 0.2  # interval for writing VTK files
 y_top = 0
 y_bottom = y_top - .25
 x_left = 0
@@ -99,16 +94,16 @@ def fluxes_from_temperature_full_domain(f, v_vec, k):
 
 mesh = RectangleMesh(p0, p1, nx, ny)
 V = FunctionSpace(mesh, 'P', 1)
+V_g = VectorFunctionSpace(mesh, 'P', 1)
 
 alpha = 1  # m^2/s, https://en.wikipedia.org/wiki/Thermal_diffusivity
 k = 100  # kg * m / s^3 / K, https://en.wikipedia.org/wiki/Thermal_conductivity
 
 # Define boundary condition
 u_D = Constant('310')
 u_D_function = interpolate(u_D, V)
-# Define flux in x direction on coupling interface (grad(u_D) in normal direction)
-f_N = Constant('0')
-f_N_function = interpolate(f_N, V)
+# We will only exchange flux in y direction on coupling interface. No initialization necessary.
+V_flux_y = V_g.sub(1)
 
 coupling_boundary = TopBoundary()
 bottom_boundary = BottomBoundary()
@@ -120,7 +115,7 @@ def fluxes_from_temperature_full_domain(f, v_vec, k):
 # Adapter definition and initialization
 precice = Adapter(adapter_config_filename="precice-adapter-config.json")
 
-precice_dt = precice.initialize(coupling_boundary, mesh, V)
+precice_dt = 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()
@@ -142,13 +137,26 @@ def fluxes_from_temperature_full_domain(f, v_vec, k):
 
 # Time-stepping
 u_np1 = Function(V)
-F_known_u = u_np1 * v / dt * dx + alpha * dot(grad(u_np1), grad(v)) * dx - u_n * v / dt * dx
 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
 
+fluxes = Function(V_g)
+fluxes.rename("Fluxes", "")
+
 while precice.is_coupling_ongoing():
 
     if precice.is_action_required(precice.action_write_iteration_checkpoint()):  # write checkpoint
@@ -165,8 +173,9 @@ def fluxes_from_temperature_full_domain(f, v_vec, k):
     solve(a == L, u_np1, bcs)
 
     # Dirichlet problem obtains flux from solution and sends flux on boundary to Neumann problem
-    fluxes = fluxes_from_temperature_full_domain(F_known_u, V, k)
-    precice.write_data(fluxes)
+    determine_heat_flux(V_g, u_np1, k, fluxes)
+    fluxes_y = fluxes.sub(1)  # only exchange y component of flux.
+    precice.write_data(fluxes_y)
 
     precice_dt = precice.advance(dt(0))
 
@@ -177,17 +186,17 @@ def fluxes_from_temperature_full_domain(f, v_vec, k):
         n = n_cp
     else:  # update solution
         u_n.assign(u_np1)
-        t += dt
+        t += float(dt)
         n += 1
 
     if precice.is_time_window_complete():
-        # if abs(t % dt_out) < 10e-5:  # output if t is a multiple of dt_out
+        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
 
     # Update dirichlet BC
-    u_D.t = t + dt(0)
-
-file_out << u_n
+    u_D.t = t + float(dt)
 
 # Hold plot
 precice.finalize()

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/.gitignore renamed to FSI/cylinderFlap_2D/OpenFOAM-FEniCS/.gitignore

@@ -1,11 +1,14 @@
-Fluid
+!Fluid
+Fluid/*
 !Fluid/0
 !Fluid/0.orig
 !Fluid/constant
+Fluid/constant/polyMesh
 !Fluid/system
 !Fluid/Fluid.foam
 !Fluid/precice-adapter-config.yml
-Solid
+!Solid
+Solid/*
 !Solid/cyl-flap.py
 !Solid/precice-adapter-config-fsi-s.json
 *.log