Split partitioned-heat-conduction case into complex and simple one

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/heat.py → partitioned-heat-conduction-complex/fenics/heat.py

@@ -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/my_enums.py → partitioned-heat-conduction-complex/fenics/my_enums.py

@@ -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-complex/fenics/precice-adapter-config-D.json

@@ -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-complex/fenics/precice-adapter-config-N.json

@@ -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-complex/precice-config.xml

@@ -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>

partitioned-heat-conduction/README.md

@@ -1,8 +1,70 @@
 ---
 title: Partitioned heat conduction
 permalink: tutorials-partitioned-heat-conduction.html
-keywords:
-summary:
+keywords: FEniCS, Nutils, Heat conduction
+summary: In this tutorial we solve a simple heat equation. The domain is partitioned and the coupling is established in a Dirichlet-Neumann fashion.
 ---
 
 {% 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
+
+In this tutorial we solve a partitioned heat equation. For information on the non-partitioned case, please refer to [1, p.37ff]. In this tutorial the computational domain is partitioned and coupled via preCICE. The coupling roughly follows the approach described in [2].
+
+![Case setup of partitioned-heat-conduction case](images/setup.png)
+
+Case setup from [3]. `D` denotes the Dirichlet participant and `N` denotes the Neumann participant.
+
+The heat equation is solved on a rectangular domain `Omega = [0,2] x [0,1]` with given Dirichlet boundary conditions. We split the domain at `x_c = 1` using a straight vertical line, the coupling interface. The left part of the domain will be referred to as the Dirichlet partition and the right part as the Neumann partition. To couple the two participant we use Dirichlet-Neumann coupling. Here, the Dirichlet participant receives Dirichlet boundary conditions (`Temperature`) at the coupling interface and solves the heat equation using these boundary conditions on the left part of the domain. Then the Dirichlet participant computes the resulting heat flux (`Flux`) from the solution and sends it to the Neumann participant. The Neumann participant uses the flux as a Neumann boundary condition to solve the heat equation on the right part of the domain. We then extract the temperature from the solution and send it back to the Dirichlet participant. This establishes the coupling between the two participants.
+
+This simple case allows us to compare the solution for the partitioned case to a known analytical solution (method of manufactures solutions, see [1, p.37ff]). For more usage examples and details, please refer to [3, sect. 4.1].
+
+## Available solvers and dependencies
+
+You can either couple a solver with itself or different solvers with each other. In any case you will need to have preCICE and the python bindings installed on your system.
+
+* `fenics`, requires you to install [FEniCS](https://fenicsproject.org/download/) and the [FEniCS-adapter](https://github.com/precice/fenics-adapter). The code is largely based on this [fenics-tutorial](https://github.com/hplgit/fenics-tutorial/blob/master/pub/python/vol1/ft03_heat.py) from [1].
+
+
+* :construction: This case is still under construction. See https://github.com/precice/tutorials/issues/152. :construction: `nutils`, requires you to install [Nutils](http://www.nutils.org/en/latest/).
+
+## Running the simulation
+
+For choosing whether you want to run the Dirichlet-kind and a Neumann-kind participant, please provide the following commandline input:
+
+* `-d` flag will enforce Dirichlet boundary conditions on the coupling interface.
+* `-n` flag will enforce Neumann boundary conditions on the coupling interface.
+
+For running the case, open two terminals and:
+
+```
+cd fenics
+python3 heat.py -d
+```
+
+and
+
+```
+cd fenics
+python3 heat.py -n
+```
+
+If you want to use Nutils for one or both sides of the setup, just `cd nutils`. The FEniCS case also supports parallel runs. Here, simply execute
+
+```
+mpirun -n <N_PROC> python3 heat.py -d
+```
+
+## Visualization
+
+Output is written into the folder `fenics/out`. You can visualize the content with paraview by opening the `*.pvd` files. The files `Dirichlet.pvd` and `Neumann.pvd` correspond to the numerical solution of the Dirichlet, respectively Neumann, problem, while the files with the prefix `ref` correspond to the analytical reference solution, the files with `error` show the error and the files with `ranks` the ranks of the solvers (if executed in parallel).
+
+![Animation of the partitioned heat equation](images/HT_FEniCS_movie.gif)
+
+Visualization in paraview for `x_c = 1.5`.
+
+## References
+
+[1] Hans Petter Langtangen and Anders Logg. "Solving PDEs in Minutes-The FEniCS Tutorial Volume I." (2016). [pdf](https://fenicsproject.org/pub/tutorial/pdf/fenics-tutorial-vol1.pdf)  
+[2] Azahar Monge and Philipp Birken. "Convergence Analysis of the Dirichlet-Neumann Iteration for Finite Element Discretizations." (2016). Proceedings in Applied Mathematics and Mechanics. [doi](https://doi.org/10.1002/pamm.201610355)  
+[3] Benjamin Rüth, Benjamin Uekermann, Miriam Mehl, Philipp Birken, Azahar Monge, and Hans Joachim Bungartz. "Quasi-Newton waveform iteration for partitioned surface-coupled multiphysics applications." (2020). International Journal for Numerical Methods in Engineering. [doi](https://doi.org/10.1002/nme.6443)