Split partitioned-heat-conduction case into complex and simple one (#148)

* Strip case of everything besides simple heat conduction setup.

* Adds precice-config.xml supporting parallel runs.

* Add documentation and nutils case.

* Add complex case. Based on ae77ae5.

* Cleanup and Update to newest version of python bindings.

* Use scalar-valued flux.

* Apply new structure to complex case.

* Remove subcycling.

* Simplify heat nutils to one mesh, make config consistent

* Fix names.

* BCs working.

* Minor restructuring

* Fix initial condition and output.

* Mark nutils case as under construction. See #152.

* Add link to Nutils

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

Co-authored-by: Benjamin Uekermann <[email protected]>

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

Co-authored-by: Benjamin Uekermann <[email protected]>

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

Co-authored-by: Benjamin Uekermann <[email protected]>

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

Co-authored-by: Benjamin Uekermann <[email protected]>

* Add note for website.

* Provide more details and use x_c = 1

* Cleanup for consistency with non-complex case.

* Fix formatting

* Cleanup w.r.t visualization

* Add summary.

* Remove FEniCS from setup.

Co-authored-by: uekerman <[email protected]>

Author: BenjaminRodenberg

partitioned-heat-conduction-complex/README.md

@@ -0,0 +1,27 @@
+---
+title: Partitioned heat conduction (complex setup)
+permalink: tutorials-partitioned-heat-conduction-complex.html
+keywords: FEniCS, Heat conduction
+summary: This tutorial is the advanced version of the `partitioned-heat-conduction` tutorial. More advanced features and geometries may be used.
+---
+
+{% include important.html content="We have not yet ported the documentation of the preCICE tutorials from the preCICE wiki to here. Please go to the [preCICE wiki](https://github.com/precice/precice/wiki#2-getting-started---tutorials)" %}
+
+## Setup
+
+This case is an advanced version of `partitioned-heat-conduction`. Some advanced features offered by this case:
+
+* Geometries may be chosen arbitrarily. One possibility is to use a circle and a rectangular plate with a hole, but you can also provide your own geometry, if you want.
+* You may combine arbitrary mesh resolutions at the coupling interface.
+* Nearest projection mapping is used.
+* The Dirichlet and Neumann participants may be swapped arbitrarily.
+* The exchanged temperature is still scalar valued, but the heat flux is vector valued.
+* You can decide to use a time dependent heat flux and right-hand side to make the problem more challenging.
+
+## Available solvers and dependencies
+
+See `partitioned-heat-conduction`, only `fenics` is provided as a solver.
+
+## Running the simulation
+
+See `partitioned-heat-conduction`. The additional featured mentioned above can be activated via command line arguments. Please run `python3 fenics/heat.py --help` for a full list of provided arguments.

partitioned-heat-conduction/fenics-fenics/.gitignore renamed to partitioned-heat-conduction-complex/fenics/.gitignore

@@ -29,7 +29,7 @@
     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, DomainPart
+from my_enums import ProblemType, DomainPart
 import argparse
 import numpy as np
 from problem_setup import get_geometry, get_problem_setup
@@ -53,9 +53,12 @@ def determine_gradient(V_g, u, flux):
     solve(a == L, flux)
 
 
-parser = argparse.ArgumentParser(description='Solving heat equation for simple or complex interface case')
-parser.add_argument("-d", "--dirichlet", help="create a dirichlet problem", dest='dirichlet', action='store_true')
-parser.add_argument("-n", "--neumann", help="create a neumann problem", dest='neumann', action='store_true')
+parser = argparse.ArgumentParser(
+    description='Solving heat equation for simple or complex interface case')
+parser.add_argument("-d", "--dirichlet", help="create a dirichlet problem",
+                    dest='dirichlet', action='store_true')
+parser.add_argument("-n", "--neumann", help="create a neumann problem",
+                    dest='neumann', action='store_true')
 parser.add_argument("-g", "--gamma", help="parameter gamma to set temporal dependence of heat flux", default=0.0,
                     type=float)
 parser.add_argument("-a", "--arbitrary-coupling-interface",
@@ -69,32 +72,9 @@ def determine_gradient(V_g, u, flux):
 
 args = parser.parse_args()
 
-subcycle = Subcycling.NONE
-
-fenics_dt = None
-error_tol = None
-
-# for all scenarios, we assume precice_dt == .1
-if subcycle is Subcycling.NONE and not args.arbitrary_coupling_interface:
-    fenics_dt = .1  # time step size
-    error_tol = 10 ** -6  # Error is bounded by coupling accuracy. In theory we would obtain the analytical solution.
-    print("Subcycling = NO, Arbitrary coupling interface = NO, error tolerance = {}".format(error_tol))
-elif subcycle is Subcycling.NONE and args.arbitrary_coupling_interface:
-    fenics_dt = .1  # time step size
-    error_tol = 10 ** -3  # error low, if we do not subcycle. In theory we would obtain the analytical solution.
-    # TODO For reasons, why we currently still have a relatively high error, see milestone https://github.com/precice/fenics-adapter/milestone/1
-    print("Subcycling = NO, Arbitrary coupling interface = YES, error tolerance = {}".format(error_tol))
-elif subcycle is Subcycling.MATCHING:
-    fenics_dt = .01  # time step size
-    error_tol = 10 ** -2  # error increases. If we use subcycling, we cannot assume that we still get the exact solution.
-    # TODO Using waveform relaxation, we should be able to obtain the exact solution here, as well.
-    print("Subcycling = YES, Matching. error tolerance = {}".format(error_tol))
-elif subcycle is Subcycling.NONMATCHING:
-    fenics_dt = .03  # time step size
-    error_tol = 10 ** -1  # error increases. If we use subcycling, we cannot assume that we still get the exact solution.
-    # TODO Using waveform relaxation, we should be able to obtain the exact solution here, as well.
-    print("Subcycling = YES, Non-matching. error tolerance = {}".format(error_tol))
-
+fenics_dt = .1
+# Error is bounded by coupling accuracy. In theory we can obtain the analytical solution.
+error_tol = 10 ** -6
 alpha = 3  # parameter alpha
 beta = 1.3  # parameter beta
 gamma = args.gamma  # parameter gamma, dependence of heat flux on time
@@ -112,7 +92,8 @@ def determine_gradient(V_g, u, flux):
                  beta=beta, gamma=gamma, t=0)
 u_D_function = interpolate(u_D, V)
 # Define flux in x direction on coupling interface (grad(u_D) in normal direction)
-f_N = Expression(("2 * gamma*t*x[0] + 2 * (1-gamma)*x[0]", "2 * alpha*x[1]"), degree=1, gamma=gamma, alpha=alpha, t=0)
+f_N = Expression(("2 * gamma*t*x[0] + 2 * (1-gamma)*x[0]",
+                  "2 * alpha*x[1]"), degree=1, gamma=gamma, alpha=alpha, t=0)
 f_N_function = interpolate(f_N, V_g)
 
 # Define initial value
@@ -124,10 +105,12 @@ 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_dt = 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_dt = precice.initialize(
+        coupling_boundary, read_function_space=V_g, write_object=u_D_function)
 
 boundary_marker = False
 
@@ -152,7 +135,8 @@ def determine_gradient(V_g, u, flux):
     # modify Neumann boundary condition on coupling interface, modify weak form correspondingly
     if not boundary_marker:  # there is only 1 Neumann-BC which is at the coupling boundary -> integration over whole boundary
         if coupling_expression.is_scalar_valued():
-            F += v * coupling_expression * dolfin.ds  # this term has to be added to weak form to add a Neumann BC (see e.g. p. 83ff Langtangen, Hans Petter, and Anders Logg. "Solving PDEs in Python The FEniCS Tutorial Volume I." (2016).)
+            # this term has to be added to weak form to add a Neumann BC (see e.g. p. 83ff Langtangen, Hans Petter, and Anders Logg. "Solving PDEs in Python The FEniCS Tutorial Volume I." (2016).)
+            F += v * coupling_expression * dolfin.ds
         elif coupling_expression.is_vector_valued():
             normal = FacetNormal(mesh)
             F += -v * dot(normal, coupling_expression) * dolfin.ds
@@ -199,15 +183,17 @@ def determine_gradient(V_g, u, flux):
 error_out << error_pointwise
 
 # set t_1 = t_0 + dt, this gives u_D^1
-u_D.t = t + dt(0)  # call dt(0) to evaluate FEniCS Constant. Todo: is there a better way?
+# call dt(0) to evaluate FEniCS Constant. Todo: is there a better way?
+u_D.t = t + dt(0)
 f.t = t + dt(0)
 
 flux = Function(V_g)
 flux.rename("Flux", "")
 
 while precice.is_coupling_ongoing():
 
-    if precice.is_action_required(precice.action_write_iteration_checkpoint()):  # write checkpoint
+    # write checkpoint
+    if precice.is_action_required(precice.action_write_iteration_checkpoint()):
         precice.store_checkpoint(u_n, t, n)
 
     read_data = precice.read_data()
@@ -235,7 +221,8 @@ def determine_gradient(V_g, u, flux):
 
     precice_dt = precice.advance(dt(0))
 
-    if precice.is_action_required(precice.action_read_iteration_checkpoint()):  # roll back to checkpoint
+    # roll back to checkpoint
+    if precice.is_action_required(precice.action_read_iteration_checkpoint()):
         u_cp, t_cp, n_cp = precice.retrieve_checkpoint()
         u_n.assign(u_cp)
         t = t_cp
@@ -248,7 +235,8 @@ def determine_gradient(V_g, u, flux):
     if precice.is_time_window_complete():
         u_ref = interpolate(u_D, V)
         u_ref.rename("reference", " ")
-        error, error_pointwise = compute_errors(u_n, u_ref, V, total_error_tol=error_tol)
+        error, error_pointwise = compute_errors(
+            u_n, u_ref, V, total_error_tol=error_tol)
         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))

partitioned-heat-conduction/fenics-fenics/Allclean renamed to partitioned-heat-conduction-complex/fenics/Allclean

@@ -1,5 +1,6 @@
 from enum import Enum
 
+
 class ProblemType(Enum):
     """
     Enum defines problem type. Details see above.
@@ -16,15 +17,3 @@ class DomainPart(Enum):
     RIGHT = 2  # right part of domain in simple interface case
     CIRCULAR = 3  # circular part of domain in complex interface case
     RECTANGLE = 4  # domain excluding circular part of complex interface case
-
-
-class Subcycling(Enum):
-    """
-    Enum defines which kind of subcycling is used
-    """
-    NONE = 0  # no subcycling, precice_dt == fenics_dt
-    MATCHING = 1  # subcycling, where fenics_dt fits into precice_dt, mod(precice_dt, fenics_dt) == 0
-    NONMATCHING = 2  # subcycling, where fenics_dt does not fit into precice_dt, mod(precice_dt, fenics_dt) != 0
-
-    # note: the modulo expressions above should be understood in an exact way (no floating point round off problems. For
-    # details, see https://stackoverflow.com/questions/14763722/python-modulo-on-floats)

partitioned-heat-conduction/fenics-fenics/errorcomputation.py renamed to partitioned-heat-conduction-complex/fenics/errorcomputation.py

@@ -0,0 +1,9 @@
+{
+  "participant_name": "Dirichlet",
+  "config_file_name": "../precice-config.xml",
+  "interface": {
+    "coupling_mesh_name": "Dirichlet-Mesh",
+    "write_data_name": "Flux",
+    "read_data_name": "Temperature"
+  }
+}

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

@@ -0,0 +1,9 @@
+{
+  "participant_name": "Neumann",
+  "config_file_name": "../precice-config.xml",
+  "interface": {
+    "coupling_mesh_name": "Neumann-Mesh",
+    "write_data_name": "Temperature",
+    "read_data_name": "Flux"
+  }
+}

partitioned-heat-conduction/fenics-fenics/my_enums.py renamed to partitioned-heat-conduction-complex/fenics/my_enums.py

@@ -0,0 +1,61 @@
+<?xml version="1.0"?>
+
+<precice-configuration>
+
+  <log>
+    <sink filter="%Severity% > debug and %Rank% = 0" format="---[precice] %ColorizedSeverity% %Message%" enabled="true"/>
+  </log>
+
+  <solver-interface dimensions="2">
+
+    <data:scalar name="Temperature"/>
+    <data:vector name="Flux"/>
+
+    <mesh name="Dirichlet-Mesh">
+      <use-data name="Temperature"/>
+      <use-data name="Flux"/>
+    </mesh>
+
+    <mesh name="Neumann-Mesh">
+      <use-data name="Temperature"/>
+      <use-data name="Flux"/>
+    </mesh>
+
+    <participant name="Dirichlet">
+      <use-mesh name="Dirichlet-Mesh" provide="yes"/>
+      <use-mesh name="Neumann-Mesh" from="Neumann"/>
+      <write-data name="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>
+
+    <participant name="Neumann">
+      <use-mesh name="Neumann-Mesh" provide="yes"/>
+      <use-mesh name="Dirichlet-Mesh" from="Dirichlet"/>
+      <write-data name="Temperature" mesh="Neumann-Mesh"/>
+      <read-data name="Flux" mesh="Neumann-Mesh"/>
+      <mapping:nearest-projection direction="read" from="Dirichlet-Mesh" to="Neumann-Mesh" constraint="consistent" />
+    </participant>
+
+    <m2n:sockets from="Dirichlet" to="Neumann" exchange-directory=".."/>
+
+    <coupling-scheme:serial-implicit>
+      <participants first="Dirichlet" second="Neumann"/>
+      <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="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="Temperature" mesh="Neumann-Mesh" limit="1e-5"/>
+      <acceleration:IQN-ILS>
+        <data name="Temperature" mesh="Neumann-Mesh"/>
+        <initial-relaxation value="0.1"/>
+        <max-used-iterations value="10"/>
+        <time-windows-reused value="5"/>
+        <filter type="QR2" limit="1e-3"/>
+      </acceleration:IQN-ILS>
+    </coupling-scheme:serial-implicit>
+
+  </solver-interface>
+</precice-configuration>