Partitioned heat conduction with Schwarz-type domain decomposition

partitioned-heat-conduction-overlap/README.md

@@ -0,0 +1,27 @@
+---
+title: Partitioned heat conduction with overlapping Schwarz-type domain decomposition
+permalink: tutorials-partitioned-heat-conduction-overlap.html
+keywords: FEniCS, Nutils, Heat conduction
+summary: We solve a simple heat equation. The domain is partitioned and the coupling is established in an overlapping-Schwarz-type domain decomposition.
+---
+
+{% note %}
+Get the [case files of this tutorial](https://github.com/precice/tutorials/tree/master/partitioned-heat-conduction-overlap). Read how in the [tutorials introduction](https://www.precice.org/tutorials.html).
+{% endnote %}
+
+## Setup
+
+We solve a partitioned heat equation, but apply an overlapping Schwarz-type domain decomposition method in this tutorial.
+
+![Case setup of partitioned-heat-conduction case with overlapping Schwarz-type domain decomposition](images/tutorials-partitioned-heat-conduction-overlap-setup.png)
+
+## Running the simulation
+
+This tutorial is for FEniCS.
+
+For choosing whether you want to run the left or right participant, please provide the following commandline input:
+
+* `python3 heat.py left` flag will run the left participant.
+* `python3 heat.py right` flag will run the right participant.
+
+Like for the case `partitioned-heat-conduction` (using Dirichlet-Neumann coupling), we can also expect for the overlapping domain decomposition applied here to recover the analytical solution. `errorcomputation.py` checks this explicitly, by comparing the numerical to the analytical solution and raising an error, if the approximation error is not within a given tolerance.

partitioned-heat-conduction-overlap/clean-tutorial.sh

@@ -0,0 +1 @@
+../tools/clean-tutorial-base.sh

partitioned-heat-conduction-overlap/fenics/clean.sh

@@ -0,0 +1,6 @@
+#!/bin/sh
+set -e -u
+
+. ../../tools/cleaning-tools.sh
+
+clean_fenics .

partitioned-heat-conduction-overlap/fenics/errorcomputation.py

@@ -0,0 +1,15 @@
+from fenics import inner, assemble, dx, project, sqrt
+
+
+def compute_errors(u_approx, u_ref, v, total_error_tol=10 ** -4):
+    # compute pointwise L2 error
+    error_normalized = (u_ref - u_approx) / u_ref
+    # project onto function space
+    error_pointwise = project(abs(error_normalized), v)
+    # determine L2 norm to estimate total error
+    error_total = sqrt(assemble(inner(error_pointwise, error_pointwise) * dx))
+    error_pointwise.rename("error", " ")
+
+    assert (error_total < total_error_tol)
+
+    return error_total, error_pointwise

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

@@ -0,0 +1,254 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# The basic example is taken from *Langtangen, Hans Petter, and Anders Logg. Solving PDEs in Python: The FEniCS
+# Tutorial I. Springer International Publishing, 2016.*
+#
+# The example code has been extended with preCICE API calls and mixed boundary conditions to allow for a Dirichlet-Neumann
+# coupling of two separate heat equations.
+#
+# The original source code can be found on https://github.com/hplgit/fenics-tutorial/blob/master/pub/python/vol1/ft03_heat.py.
+#
+# Heat equation with Dirichlet conditions. (Dirichlet problem, $f_N = \mathcal{D}\left(u_C\right)$)
+#
+# \begin{align*}
+# \frac{\partial u}{\partial t} &= \Delta u + f && \text{on the unit square $\Omega = \left[0,1\right] \times \left[0,1\right]$}\\
+# u &= u_C             && \text{on the coupling boundary $\Gamma_C$ at $x = 1$}\\
+# u &= u_D             && \text{on the remaining boundary $\Gamma = \partial \Omega \setminus \Gamma_C$}\\
+# u &= u_0             && \text{on $\Omega$ at t = 0}\\
+# u &= 1 + x^2 + \alpha y^2 + \beta t \\
+# f &= \beta - 2 - 2\alpha \\
+# \end{align*}
+#
+# Heat equation with mixed boundary conditions. (Neumann problem, $u_C = \mathcal{N}\left(f_N\right)$)
+#
+# \begin{align*}
+# \frac{\partial u}{\partial t} &= \Delta u + f && \text{on the shifted unit square $\Omega = \left[1,2\right] \times \left[0,1\right]$}\\
+# \frac{\text{d}u}{\text{d}\vec{n}} &= f_N             && \text{on the coupling boundary $\Gamma_C$ at $x = 1$}\\
+# u &= u_D             && \text{on the remaining boundary $\Gamma = \partial \Omega \setminus \Gamma_C$}\\
+# u &= u_0             && \text{on $\Omega$ at t = 0}\\
+# u &= 1 + x^2 + \alpha y^2 + \beta t \\
+# f &= \beta - 2 - 2\alpha \\
+# \end{align*}
+
+# In[1]:
+
+
+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, MeshFunction, MPI
+from fenics import SubDomain, Point, RectangleMesh, near, Function, Expression
+from fenicsprecice import Adapter
+from errorcomputation import compute_errors
+import numpy as np
+import dolfin
+from dolfin import FacetNormal, dot
+from enum import Enum
+import sys
+
+
+class DomainPart(Enum):
+    """
+    Enum defines which part of the domain [x_left, x_right] x [y_bottom, y_top] we compute.
+    """
+    LEFT = 1  # left part of domain in simple interface case
+    RIGHT = 2  # right part of domain in simple interface case
+
+
+fenics_dt = .1  # time step size
+# Error is bounded by coupling accuracy. In theory we would obtain the analytical solution.
+error_tol = 10**-6
+alpha = 3  # parameter alpha
+beta = 1.3  # parameter beta
+
+
+if sys.argv[1] == 'left':
+    domain_part = DomainPart.LEFT
+if sys.argv[1] == 'right':
+    domain_part = DomainPart.RIGHT
+
+y_bottom, y_top = 0, 1
+x_left, x_right = 0, 2
+# x coordinate of coupling interface; for Schwarz Domain Decomposition coupling interface sits between The DoFs.
+x_coupling = 1.0
+
+overlap_cells = 1
+
+nx = 9 + overlap_cells
+ny = 9
+hx = (x_right - x_left) / (2 * nx - overlap_cells)
+
+
+if domain_part is DomainPart.LEFT:
+    p0 = Point(x_left, y_bottom)
+    # rightmost point is the DoF hx/2 right of the coupling interface that belongs to the Right participant
+    x_schwarz_read = x_coupling + hx * (overlap_cells * 0.5)
+    x_schwarz_write = x_coupling - hx * (overlap_cells * 0.5)
+    p1 = Point(x_schwarz_read, y_top)
+elif domain_part is DomainPart.RIGHT:
+    # leftmost point is the DoF hx/2 left of the coupling interface that belongs to the Left participant
+    x_schwarz_read = x_coupling - hx * (overlap_cells * 0.5)
+    x_schwarz_write = x_coupling + hx * (overlap_cells * 0.5)
+    p0 = Point(x_schwarz_read, y_bottom)
+    p1 = Point(x_right, y_top)
+
+
+class ExcludeStraightBoundary(SubDomain):
+    def get_user_input_args(self, args):
+        self._interface = args.interface
+
+    def inside(self, x, on_boundary):
+        tol = 1E-14
+        if on_boundary and not near(x[0], x_schwarz_read, tol) or near(x[1], y_top, tol) or near(x[1], y_bottom, tol):
+            return True
+        else:
+            return False
+
+
+class OverlapDomain(SubDomain):
+    def inside(self, x, on_boundary):
+        tol = 1E-14
+        if (x[0] <= x_coupling + hx * (overlap_cells * 0.5) + tol) and (x[0] >= x_coupling -
+                                                                        hx * (overlap_cells - 0.5) - tol):  # Point lies inside of overlapping domain
+            return True
+        else:
+            return False
+
+
+class ReadBoundary(SubDomain):
+    def inside(self, x, on_boundary):
+        tol = 1E-14
+        if on_boundary and near(x[0], x_schwarz_read, tol):
+            return True
+        else:
+            return False
+
+
+class WriteBoundary(SubDomain):
+    def inside(self, x, on_boundary):
+        tol = 1E-14
+        if near(x[0], x_schwarz_write, tol):  # Point lies inside of the domain!
+            return True
+        else:
+            return False
+
+
+mesh = RectangleMesh(p0, p1, nx, ny, diagonal="left")
+read_boundary = ReadBoundary()
+write_boundary = WriteBoundary()
+remaining_boundary = ExcludeStraightBoundary()
+
+
+# Define function space using mesh
+V = FunctionSpace(mesh, 'P', 2)
+
+# Define boundary conditions
+u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t', degree=2, alpha=alpha, beta=beta, t=0)
+u_D_function = interpolate(u_D, V)
+
+# Define initial value
+u_n = interpolate(u_D, V)
+u_n.rename("Temperature", "")
+
+
+precice, precice_dt, initial_data = None, 0.0, None
+
+# Initialize the adapter according to the specific participant
+if domain_part is DomainPart.LEFT:
+    precice = Adapter(adapter_config_filename="precice-adapter-config-L.json")
+elif domain_part is DomainPart.RIGHT:
+    precice = Adapter(adapter_config_filename="precice-adapter-config-R.json")
+
+precice_dt = precice.initialize(OverlapDomain(), read_function_space=V, write_object=u_D_function)
+
+dt = Constant(0)
+dt.assign(np.min([fenics_dt, precice_dt]))
+
+# Define variational problem
+u = TrialFunction(V)
+v = TestFunction(V)
+f = Expression('beta - 2 - 2*alpha', degree=2, alpha=alpha, beta=beta, t=0)
+F = u * v / dt * dx + dot(grad(u), grad(v)) * dx - (u_n / dt + f) * v * dx
+
+bcs = [DirichletBC(V, u_D, remaining_boundary)]
+coupling_expression = precice.create_coupling_expression()
+bcs.append(DirichletBC(V, coupling_expression, read_boundary))
+
+a, L = lhs(F), rhs(F)
+
+# Time-stepping
+u_np1 = Function(V)
+u_np1.rename("Temperature", "")
+t = 0
+
+# reference solution at t=0
+u_ref = interpolate(u_D, V)
+u_ref.rename("reference", " ")
+
+# Generating output files
+temperature_out = File("output/%s.pvd" % precice.get_participant_name())
+ref_out = File("output/ref%s.pvd" % precice.get_participant_name())
+error_out = File("output/error%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
+
+error_total, error_pointwise = compute_errors(u_n, u_ref, V)
+error_out << error_pointwise
+
+# set t_1 = t_0 + dt, this gives u_D^1
+# call dt(0) to evaluate FEniCS Constant. Todo: is there a better way?
+u_D.t = t + dt(0)
+f.t = t + dt(0)
+
+while precice.is_coupling_ongoing():
+
+    # write checkpoint
+    if precice.is_action_required(precice.action_write_iteration_checkpoint()):
+        precice.store_checkpoint(u_n, t, n)
+
+    read_data = precice.read_data()
+
+    # 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)
+
+    precice.write_data(u_np1)
+
+    precice_dt = precice.advance(dt(0))
+
+    # 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
+        n = n_cp
+    else:  # update solution
+        u_n.assign(u_np1)
+        t += float(dt)
+        n += 1
+
+    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)
+        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_np1
+    ref_out << u_ref
+    error_out << error_pointwise
+
+    # Update Dirichlet BC
+    u_D.t = t + float(dt)
+    f.t = t + float(dt)
+
+# Hold plot
+precice.finalize()

partitioned-heat-conduction-overlap/fenics/precice-adapter-config-L.json

@@ -0,0 +1,9 @@
+{
+  "participant_name": "Left",
+  "config_file_name": "../precice-config.xml",
+  "interface": {
+      "coupling_mesh_name": "Left-Mesh",
+      "write_data_name": "TemperatureLeft",
+      "read_data_name": "TemperatureRight"
+    }
+}

partitioned-heat-conduction-overlap/fenics/precice-adapter-config-R.json

@@ -0,0 +1,9 @@
+{
+  "participant_name": "Right",
+  "config_file_name": "../precice-config.xml",
+  "interface": {
+      "coupling_mesh_name": "Right-Mesh",
+      "write_data_name": "TemperatureRight",
+      "read_data_name": "TemperatureLeft"
+    }
+}

partitioned-heat-conduction-overlap/fenics/run.sh

@@ -0,0 +1,24 @@
+#!/bin/sh
+set -e -u
+
+usage() { echo "Usage: cmd [-d] [-n]" 1>&2; exit 1; }
+
+# Check if no input argument was provided
+if [ -z "$*" ] ; then
+	usage
+fi
+
+# Select appropriate case
+while getopts ":dn" opt; do
+  case ${opt} in
+  d)
+    python3 heat.py -d --error-tol 10e-3
+    ;;
+  n)
+    python3 heat.py -n --error-tol 10e-3
+    ;;
+  *)
+    usage
+    ;;
+  esac
+done