# Waveform Relaxation with asynchronous time-integration

Peter Meisrimel${}^{\mbox{\tiny\rm*},\mbox{\tiny\rm 1}}$, Philipp
Birken${}^{\mbox{\tiny\rm 1}}$

###### Abstract

We consider Waveform Relaxation (WR) methods for partitioned time-integration
of surface-coupled multiphysics problems. WR allows independent time-
discretizations on independent and adaptive time-grids, while maintaining high
time-integration orders. Classical WR methods such as Jacobi or Gauss-Seidel
WR are typically either parallel or converge quickly.

We present a novel parallel WR method utilizing asynchronous communication
techniques to get both properties. Classical WR methods exchange discrete
functions after time-integration of a subproblem. We instead asynchronously
exchange time-point solutions during time-integration and directly incorporate
all new information in the interpolants. We show both continuous and time-
discrete convergence in a framework that generalizes existing linear WR
convergence theory. An algorithm for choosing optimal relaxation in our new WR
method is presented.

Convergence is demonstrated in two conjugate heat transfer examples. Our new
method shows an improved performance over classical WR methods. In one example
we show a partitioned coupling of the compressible Euler equations with a
nonlinear heat equation, with subproblems implemented using the open source
libraries DUNE and FEniCS.

${}^{\mbox{\tiny\rm 1}}$Centre for the mathematical sciences, Numerical
Analysis, Lund University, Lund, Sweden
email<EMAIL_ADDRESS><EMAIL_ADDRESS>${}^{\mbox{\tiny\rm*}}$Corresponding author

Keywords: Asynchronous iteration, Waveform Relaxation, Dynamic Iteration,
Coupled Problems, Thermal Fluid-Structure Interaction
Mathematics Subject Classification (2000): 65B99, 65F99, 65L05, 65Y05, 80M10,
80M25
The authors gratefully acknowledge support from the Swedish e-science
collaboration eSSENCE.

##  1 Introduction

We consider multiphysics problems, which are comprised of coupled systems with
different physics. In particular, we consider problems with a bidirectional
surface coupling. I.e., the subproblems interact via a lower dimensional
interface. Examples are within fluid structure interaction in the simulation
of blood flow in large arteries [4], cooling of rocket engines [12, 13] or gas
quenching [33].

We follow the partitioned approach, which allows re-use of existing codes and
solving the subproblems with different computational methods on individual
grids. Our focus is time-integration. We want to solve subproblems using
independent and higher-order time-discretizations on adaptive time-grids.
Additionally, we want to perform time-integration of the subproblems in
parallel, on top of a parallelization in space.

A technique that promises to meet all these requirements is the so called
Waveform relaxation (WR). An iteration requires solving the subproblems on a
time window. Thereby, continuous interface functions, obtained via suitable
interpolation, are provided from the respective other problem. WR methods were
originally introduced in [15] for systems of ordinary differential equations
(ODEs), and used for the first time to solve time dependent PDEs in [6, 7]. WR
appears in the literature under a variety of names: Waveform
relaxation/iteration, dynamic iteration/relaxation and Picard(-Lindelöf)
iteration.

The most common type of WR methods are Gauss-Seidel (GS) WR, solving all
subproblems in sequence, and Jacobi WR, which solves all subproblems in
parallel. However, the parallelism of Jacobi WR typically comes at the cost of
slower convergence rates, due to less information exchange.

In this article, we construct a novel and inherently parallel WR method with
more information exchange than Jacobi WR. Our ansatz to increase communication
is to exchange the results of each timestep directly after computation. Any
new information is directly incorporated by updating the interpolants,
affecting their subsequent evaluations. This increases the information
exchange and thus enhances convergence rates. We use asynchronous One-sided-
communication that allows solving subproblems in parallel and does not require
function calls on the receiving processor.

WR methods require convergence acceleration to achieve fast convergence rates.
We consider classical convergence acceleration by weighting updates using
relaxation parameters, which are highly problem specific [11, 21]. Other
acceleration techniques involve using an additional convolution relaxation
term [28, 11] or Krylov-subspace acceleration [17], see [17] for a wider
overview of different acceleration techniques. However, many of these are not
applicable in the partitioned approach. Black-box convergence acceleration
techniques such as quasi-Newton methods can also be applied to WR and have
been shown to work well [29].

With classical WR methods, data dependencies between the subproblems are
fixed. In our new method, dependencies can vary in time and differ in each
iteration, since One-sided-communication is not deterministic. We present an
analytical description of our new method and convergence proofs in the
continuous and time-discrete setting for linear problems. This generalizes
existing WR theory [11].

We present an algorithm for optimal relaxation in our new method. Here,
optimal relaxation critically depends on the realized communication, which is
not deterministic. Thus, we deduce the realized communication between the
subsolvers in every timestep and choose suitable relaxation for each time-
point solution.

We demonstrate our method using two conjugate heat transfer test cases,
showing convergence in both. In the first test case, two coupled heterogeneous
linear heat equations, performance results show a runtime speed-up of our new
method compared to classical Jacobi and GS WR methods. The second experiment
is a gas quenching test case, which consists of the compressible Euler
equations coupled to a nonlinear heat equation. Here, we demonstrate a black-
box coupling of heterogeneous space discretizations and subsolver codes. The
fluid is solved using a finite volume discretization implemented in DUNE [1]
and the solid is solved via a finite element discretization implemented in
FEniCS [16].

The paper is structured as follows: We first introduce general continuous and
time-discrete WR methods in Sections 2 and 3. In Section 4 we provide a brief
overview over the principles of one-sided asynchronous communication. We
formalize our new approach in Section 5 and present a first algorithm. We
discuss convergence in the linear case in Section 6, showing time-discrete and
continuous convergence. Our algorithm for choosing optimal relaxation for two
coupled problems is shown in Section 7. Finally, we show numerical results,
followed by summary and conclusions.

##  2 Continuous Waveform Relaxation

Consider the following coupled system of initial value problems

$\displaystyle\begin{split}\dot{\bm{v}}(t)&=\bm{g}(t,\bm{v}(t),\bm{w}(t)),\qquad\bm{v}(0)=\bm{v}_{0}\in\mathbb{R}^{d_{v}},\\\
\dot{\bm{w}}(t)&=\bm{h}(t,\bm{v}(t),\bm{w}(t)),\qquad\bm{w}(0)=\bm{w}_{0}\in\mathbb{R}^{d_{w}},\end{split}\quad
t\in[0,T_{f}<\infty].$ (1)

We now define a general continuous Waveform Relaxation method. Given
$\bm{v}^{(k)}$ and $\bm{w}^{(k)}$, a single iteration consists of solving two
differential equations and performing two relaxation steps as follows:

$\displaystyle\dot{\hat{\bm{v}}}^{(k+1)}(t)$
$\displaystyle=\bm{g}\left(t,\hat{\bm{v}}^{(k+1)}(t),\bm{w}_{\,*}^{(k)}(t)\right),\quad\hat{\bm{v}}^{(k+1)}(0)=\bm{v}_{0},\quad
t\in[0,T_{f}],$ (2a) $\displaystyle\bm{v}^{(k+1)}$
$\displaystyle=\bm{v}^{(k)}+\bm{\Theta}_{v}(\hat{\bm{v}}^{(k+1)}-\bm{v}^{(k)}),$
(2b) $\displaystyle\dot{\hat{\bm{w}}}^{(k+1)}(t)$
$\displaystyle=\bm{h}\left(t,\bm{v}_{\,*}^{(k)}(t),\hat{\bm{w}}^{(k+1)}(t)\right),\quad\hat{\bm{w}}^{(k+1)}(0)=\bm{w}_{0},\quad
t\in[0,T_{f}],$ (2c) $\displaystyle\bm{w}^{(k+1)}$
$\displaystyle=\bm{w}^{(k)}+\bm{\Theta}_{w}(\hat{\bm{w}}^{(k+1)}-\bm{w}^{(k)}),$
(2d)

with nonsingular diagonal matrices $\bm{\Theta}_{v}\in\mathbb{R}^{d_{v}\times
d_{v}}$, $\bm{\Theta}_{w}\in\mathbb{R}^{d_{w}\times d_{w}}$ for relaxation.
Extensions to more than two systems are straight-forward, c.f. [30]. The
specific WR method is defined by the choices for $\bm{v}_{\,*}^{(k)}$ and
$\bm{w}_{\,*}^{(k)}$. The trivial initial guesses for $\bm{v}^{(0)}$ and
$\bm{w}^{(0)}$ are to extrapolate the initial value.

The most common WR methods are Gauss-Seidel (GS) and Jacobi WR, c.f., [31,
30]. Continuous GS WR is given by

$\bm{v}_{\,*}^{(k)}=\bm{v}^{(k+1)}\quad\text{and}\quad\bm{w}_{\,*}^{(k)}=\bm{w}^{(k)}.$
(3)

GS WR is sequential, which makes it sensitive to the order of the systems in
(1).

Jacobi WR is given by

$\bm{v}_{\,*}^{(k)}=\bm{v}^{(k)}\quad\text{and}\quad\bm{w}_{\,*}^{(k)}=\bm{w}^{(k)},$
(4)

which allows parallel computation of (2a), (2b) and (2c), (2d).

The iteration is commonly terminated if

$\|\bm{y}^{(k+1)}(T_{f})-\bm{y}^{(k)}(T_{f})\|<\|\bm{y}^{(k+1)}(0)\|\,TOL_{WR},$
(5)

i.e., via the relative update measured at $t=T_{f}$, where updates tend to be
the largest. Here, $\bm{y}$ is a subset of the unknowns of $\bm{v}$ resp.
$\bm{w}$ that $\bm{h}$ resp. $\bm{g}$ in (1) depends on.

##  3 Time-discrete Waveform Relaxation

We enable the use of independent time-grids and time-integration schemes by
using interpolants of the respective discrete solutions in the right-hand
sides of (2a) and (2c).

We denote discrete solutions by

$\underline{\bm{\bm{v}}}^{(k)}:=\\{\bm{v}^{(k)}_{n}\\}_{n=0,\ldots,N_{v}^{(k)}}\quad\text{and}\quad\underline{\bm{\bm{w}}}^{(k)}:=\\{\bm{w}^{(k)}_{n}\\}_{n=0,\ldots,N_{w}^{(k)}},$

on time-grids $0=t_{0}^{(v),(k)}<\ldots<t_{N_{v}^{(k)}}^{(v),(k)}=T_{f}$ and
$0=t_{0}^{(w),(k)}<\ldots<t_{N_{w}^{(k)}}^{(w),(k)}=T_{f}$. The interpolants
are as follows:

$\mathcal{I}(\underline{\bm{\bm{v}}}^{(k)})\in\mathcal{C}\left([0,T_{f}];\mathbb{R}^{d_{v}}\right)\quad\text{and}\quad\mathcal{I}(\underline{\bm{\bm{w}}}^{(k)})\in\mathcal{C}\left([0,T_{f}];\mathbb{R}^{d_{w}}\right).$

Here, we omit the time-grids as input to the interpolants for ease of
notation. We obtain a time-discrete WR method by using discrete time-
integration to solve

$\displaystyle\dot{\hat{\bm{v}}}^{(k+1)}(t)$
$\displaystyle=\bm{g}\left(t,\hat{\bm{v}}^{(k+1)}(t),\mathcal{I}(\underline{\bm{\bm{w}}}_{\,*}^{(k)})(t)\right),\quad\hat{\bm{v}}^{(k+1)}(0)=\bm{v}_{0},\quad
t\in[0,T_{f}],$ (6a) $\displaystyle\dot{\hat{\bm{w}}}^{(k+1)}(t)$
$\displaystyle=\bm{h}\left(t,\mathcal{I}(\underline{\bm{\bm{v}}}_{\,*}^{(k)})(t),\hat{\bm{w}}^{(k+1)}(t)\right),\quad\hat{\bm{w}}^{(k+1)}(0)=\bm{w}_{0},\quad
t\in[0,T_{f}].$ (6b)

Here, one chooses $\underline{\bm{\bm{v}}}_{\,*}^{(k)}$,
$\underline{\bm{\bm{w}}}_{\,*}^{(k)}$ in accordance with e.g., (3) or (4). We
consider polynomial interpolation.

Relaxation is performed in the discrete data-points as follows

$\displaystyle\bm{v}^{(k+1)}_{n}$
$\displaystyle=\left(\bm{I}-\bm{\Theta_{v}}\right)\mathcal{I}(\underline{\bm{\bm{v}}}^{(k)})(t_{n}^{(v),(k+1)})+\bm{\Theta_{v}}\hat{\bm{v}}^{(k+1)}_{n},\quad
n=1,\ldots,N_{v}^{(k+1)},$ (7a) $\displaystyle\bm{w}^{(k+1)}_{n}$
$\displaystyle=\left(\bm{I}-\bm{\Theta_{w}}\right)\mathcal{I}(\underline{\bm{\bm{w}}}^{(k)})(t_{n}^{(w),(k+1)})+\bm{\Theta_{w}}\hat{\bm{w}}^{(k+1)}_{n},\quad
n=1,\ldots,N_{w}^{(k+1)}.$ (7b)

See Algorithm 1 for a pseudocode of time-discrete Jacobi WR. We use
interpolation with evaluation at run-time. I.e., we define the interpolants
once using fixed data-structures for $\underline{\bm{\bm{v}}}_{\,*}$ and
$\underline{\bm{\bm{w}}}_{\,*}$. We then update
$\underline{\bm{\bm{v}}}_{\,*}$ and $\underline{\bm{\bm{w}}}_{\,*}$ during the
iteration, affecting the results of subsequent evaluations. Our discrete
initial guesses for $\bm{v}^{(0)}$ and $\bm{w}^{(0)}$ are to extrapolate the
initial value.

Pseudocode: Time-discrete Jacobi WR

1:Process $0$ (p$0$)  2:$\underline{\bm{\bm{w}}}^{(0)}$ discrete initial guess  3:Initialize $\underline{\bm{\bm{w}}}_{\,*}$ and $\mathcal{I}(\underline{\bm{\bm{w}}}_{\,*})$ 4:for $k=0,\ldots,k_{\max}-1$ do 5: $\underline{\bm{\bm{w}}}_{\,*}\leftarrow\underline{\bm{\bm{w}}}^{(k)}$ Update interpolant  6: $\underline{\bm{\hat{\bm{v}}}}^{(k+1)}\leftarrow$ Discr. solve (6a)  7: $\underline{\bm{\bm{v}}}^{(k+1)}\leftarrow$ Relaxation (7a)  8: $\underline{\bm{\bm{v}}}^{(k+1)}\rightarrow$ Send to p$1$ 9: $\underline{\bm{\bm{w}}}^{(k+1)}\leftarrow$ Recv. from p$1$ 10: Check (5), break if true 11:end for |  Process $1$ (p$0$)  $\underline{\bm{\bm{v}}}^{(0)}$ discrete initial guess  Initialize $\underline{\bm{\bm{v}}}_{\,*}$ and $\mathcal{I}(\underline{\bm{\bm{v}}}_{\,*})$ for $k=0,\ldots,k_{\max}-1$ do $\underline{\bm{\bm{v}}}_{\,*}\leftarrow\underline{\bm{\bm{v}}}^{(k)}$ Update interpolant  $\underline{\bm{\hat{\bm{w}}}}^{(k+1)}\leftarrow$ Discr. solve (6b)  $\underline{\bm{\bm{w}}}^{(k+1)}\leftarrow$ Relaxation (7b)  $\underline{\bm{\bm{w}}}^{(k+1)}\rightarrow$ Send to p$1$ $\underline{\bm{\bm{v}}}^{(k+1)}\leftarrow$ Recv. from p$1$ Check (5), break if true end for
---|---

List of algfloats 1 Pseudocode of the time-discrete Jacobi WR method.
Obtaining the initial guesses may involve communication. Here,
$\underline{\bm{\bm{v}}}_{\,*}^{(k)}$ and
$\underline{\bm{\bm{w}}}_{\,*}^{(k)}$ in (6) are defined by (4).

##  4 One-sided communication

The standard in parallel computations is Point-to-Point communication,
primarily using MPI_Send and MPI_Recv. Here, every MPI_Send requires a
matching MPI_Recv function call. This works well for algorithms with fixed
synchronization points, e.g., the termination check in Algorithm 1. Here, we
want to perform time-integration of the subproblems in parallel on independent
grids, exchanging information after each timestep.

We solve this by using One-sided communication, also called remote memory
access (RMA) [20, Chpt.11], which is asynchronous. Since it is non-standard,
we give a brief overview over the associated concepts and methods. It is also
worth noting that RMA via MPI is not available in mpi4py, requiring an
implementation in C++ or Fortran.

MPI_Window objects and their allocated memory facilitate RMA. One accesses the
memory associated with a given window via MPI_Get (read) or MPI_Put (write)
operations. The memory of a target window is only accessible during access
periods. We use passive target synchronization, in which a processor creates
an access period on a target window (including windows on its own memory) by
locking and unlocking the target window using MPI_Win_(un)lock. This does not
require active participation in terms of function calls on the processor of
the targeted window. An additional parameter in the MPI_Win_(un)lock call can
specify the type of lock to ensure exclusive access for writing or shared
access for reading, preventing access conflicts. The order in which locks are
obtained is determined at runtime and can differ with each execution, which
means RMA is not deterministic.

The RMA memory model uses public and private copies of windows. All MPI_Put
and MPI_Get operations act on the public copy. However, the present variable
values are given by the private window copies. As such, a process needs to
synchronize its public and private window copies to obtain any updates
received via MPI_Put operations. This is done explicitly using MPI_Win_sync.
Some RMA functions perform this synchronization implicitly, note that
MPI_Barrier does not.

##  5 Waveform relaxation with asynchronous time-integration

While Jacobi WR is parallel, convergence rates are typically slower than those
of GS WR, due to less information exchange. The goal is to develop a parallel
WR method with more information exchange than Jacobi WR and thus a faster
convergence rate. Our ansatz is to use asynchronous communication already
during time-integration to increase the information exchange.

We start with the parallel Algorithm 1 and modify it to increase
communication. First, we expose the interpolant data
$\underline{\bm{\bm{v}}}_{\,*}$ and $\underline{\bm{\bm{w}}}_{\,*}$ via
MPI_Window objects. Thus, remote updates are directly incorporated in
subsequent interpolant evaluations. Next, we move relaxation and communication
to the time-step level, i.e., into the time-stepping loop. We asynchronously
communicate new time-point solutions, remotely updating the corresponding
values in $\underline{\bm{\bm{v}}}_{\,*}$ resp.
$\underline{\bm{\bm{w}}}_{\,*}$ on the other process, using MPI_Put.

Algorithm 2 shows the pseudocode for two coupled problems, with a different
number of timesteps $N_{v}\neq N_{w}$ for each subproblem and with constant
relaxation. In Section 7 we present an algorithm with variable relaxation,
based on the realized communication.

Pseudocode: WR with asynchronous time-integration

1:Process $0$ (p$0$)  2:$\underline{\bm{\bm{w}}}^{(0)}$ discrete initial guess  3:Initialize $\underline{\bm{\bm{w}}}_{\,*}\leftarrow\underline{\bm{\bm{w}}}^{(0)}$ and $\mathcal{I}(\underline{\bm{\bm{w}}}_{\,*})$ 4:Expose $\underline{\bm{\bm{w}}}_{\,*}$ via MPI_Window 5:for $k=0,\ldots,k_{\max}-1$ do 6: for $n=1,\ldots,N_{v}$ do 7: $\hat{\bm{v}}^{(k+1)}_{n}\leftarrow$ Solve* (6a)  8: $\bm{v}^{(k+1)}_{n}\leftarrow$ Relaxation (7a)  9: $\bm{v}^{(k+1)}_{n}\rightarrow$ MPI_Put to $\underline{\bm{\bm{v}}}_{\,*}$ on p$1$ 10: end for 11: Sync. + Termination check  12:end for |  Process $1$ (p$1$)  $\underline{\bm{\bm{v}}}^{(0)}$ discrete initial guess  Initialize $\underline{\bm{\bm{v}}}_{\,*}\leftarrow\underline{\bm{\bm{v}}}^{(0)}$ and $\mathcal{I}(\underline{\bm{\bm{v}}}_{\,*})$ Expose $\underline{\bm{\bm{v}}}_{\,*}$ via MPI_Window for $k=0,\ldots,k_{\max}-1$ do for $n=1,\ldots,N_{w}$ do $\hat{\bm{w}}^{(k+1)}_{n}\leftarrow$ Solve* (6b)  $\bm{w}^{(k+1)}_{n}\leftarrow$ Relaxation (7b)  $\bm{w}^{(k+1)}_{n}\rightarrow$ MPI_Put to $\underline{\bm{\bm{w}}}_{\,*}$ on p$0$ end for Sync. + Termination check  end for
---|---

List of algfloats 2 New proposed method using asynchronous communication
during time-integration. Obtaining the initial guesses may involve
communication. Solve* denotes a single discrete timestep in solving (6a) resp.
(6b), with the interpolants defined in Line 3 in the right-hand sides.

Our new method is defined by

$(\underline{\bm{\bm{v}}}_{\,*}^{(k)})_{n}=\begin{cases}\bm{v}_{n}^{(k+1)},&\text{if
available}\\\ \bm{v}_{n}^{(k)},&\text{else}\end{cases},$ (8)

$\underline{\bm{\bm{w}}}_{\,*}^{(k)}$ analogous. Availability is determined at
run-time, by the present data when evaluating the interpolant. Since
asynchronous communication is not deterministic,
$\underline{\bm{\bm{v}}}_{\,*}^{(k)}$ and
$\underline{\bm{\bm{w}}}_{\,*}^{(k)}$ can vary for different timesteps and
with $k$. Due to remote updates, evaluations of the interpolants for the same
$t$ and $k$, but at different real-life times can differ.

The corresponding continuous WR method is (2), with $\bm{v}_{\,*}^{(k)}$ and
$\bm{w}_{\,*}^{(k)}$ varying in both time and with $k$.

###  5.1 Variable relaxation

With our new method, $\bm{v}_{\,*}^{(k)}$ and $\bm{w}_{\,*}^{(k)}$ vary with
$t$ and $k$. It is possible that we obtain Jacobi or GS WR. Since optimal
relaxation matrices can notably differ for Jacobi and GS WR, constant
relaxation is unlike to achieve optimal convergence acceleration in our new
method. We instead consider relaxation varying with $t$ and $k$:

$\displaystyle\begin{split}\bm{v}^{(k+1)}(t)&=\bm{v}^{(k)}(t)+\bm{\Theta}_{v}^{(k+1)}(t)\left(\hat{\bm{v}}^{(k+1)}(t)-\bm{v}^{(k)}(t)\right),\\\
\bm{w}^{(k+1)}(t)&=\bm{w}^{(k)}(t)+\bm{\Theta}_{w}^{(k+1)}(t)\left(\hat{\bm{w}}^{(k+1)}(t)-\bm{w}^{(k)}(t)\right),\end{split}$
(9)

where $\bm{\Theta}_{v}^{(k+1)}(t)$ and $\bm{\Theta}_{w}^{(k+1)}(t)$ are non-
singular diagonal matrices. We first discuss convergence results, since these
provide us the necessary insight on how to choose relaxation. Following that,
we extend Algorithm 2 to variable relaxation in Section 7. There, we present a
concrete approach to determining $\bm{v}_{\,*}^{(k)}$ and $\bm{w}_{\,*}^{(k)}$
at runtime, to then choose appropriate relaxation for each timestep.

##  6 Convergence Analysis

Similar to [24, 25, 9, 10] we analyze convergence in the linear setting. We
first present the established WR theory, before extending it to include our
new method. Consider the following monolithic system

$\bm{B}\dot{\bm{u}}(t)+\bm{A}\bm{u}(t)=\bm{f}(t),\quad\bm{u}(0)=\bm{u}_{0}\in\mathbb{R}^{d},\quad
t\in[0,T_{f}<\infty],\quad\bm{B}\text{ nonsingular},$ (10)

with $\bm{B}$, $\bm{A}$ $\in\mathbb{R}^{d\times d}$, and $\bm{f}$ Lipschitz-
continuous. Here, one can express classical WR methods such as Jacobi and GS
WR via constant splittings [24, 9]

$\bm{B}=\bm{M}_{B}-\bm{N}_{B}\quad\text{and}\quad\bm{A}=\bm{M}_{A}-\bm{N}_{A},\quad\bm{M}_{B}\text{
nonsingular},$ (11)

and the iteration

$\displaystyle\begin{split}\bm{M}_{B}\dot{\bm{u}}^{(k+1)}(t)+\bm{M}_{A}\bm{u}^{(k+1)}(t)=\bm{N}_{B}\dot{\bm{u}}^{(k)}(t)+\bm{N}_{A}\bm{u}^{(k)}(t)+\bm{f}(t),\\\
\bm{u}^{(k+1)}(0)=\bm{u}_{0},\quad t\in[0,T_{f}].\end{split}$ (12)

The particular splitting (11) depends on the WR method, e.g., Jacobi or GS WR,
and includes constant relaxation (2b), (2d). We omit dependencies of the
splitting matrices on the relaxation matrices in (11) for readability.

Consider for example the following system of ODEs:

$\underbrace{\begin{pmatrix}\bm{B}_{1}&\bm{B}_{2}\\\
\bm{B}_{3}&\bm{B}_{4}\end{pmatrix}}_{\bm{B}}\underbrace{\begin{pmatrix}\dot{\bm{v}}(t)\\\
\dot{\bm{w}}(t)\end{pmatrix}}_{\dot{\bm{u}}(t)}+\underbrace{\begin{pmatrix}\bm{A}_{1}&\bm{A}_{2}\\\
\bm{A}_{3}&\bm{A}_{4}\end{pmatrix}}_{\bm{A}}\underbrace{\begin{pmatrix}\bm{v}(t)\\\
\bm{w}(t)\end{pmatrix}}_{\bm{u}(t)}=\underbrace{\begin{pmatrix}\bm{f}_{1}(t)\\\
\bm{f}_{2}(t)\end{pmatrix}}_{\bm{f}(t)},\quad\underbrace{\begin{pmatrix}\bm{v}(0)\\\
\bm{w}(0)\end{pmatrix}}_{\bm{u}(0)}=\underbrace{\begin{pmatrix}\bm{v}_{0}\\\
\bm{w}_{0}\end{pmatrix}}_{\bm{u}_{0}},$

with $t\in[0,T_{f}]$ and $\bm{B}_{1}$, $\bm{B}_{4}$ nonsingular. Then, Jacobi
WR, without relaxation, is given by

$\displaystyle\begin{split}&\begin{pmatrix}\bm{B}_{1}&\bm{0}\\\
\bm{0}&\bm{B}_{4}\end{pmatrix}\dot{\bm{u}}^{(k+1)}(t)+\begin{pmatrix}\bm{A}_{1}&\bm{0}\\\
\bm{0}&\bm{A}_{4}\end{pmatrix}\bm{u}^{(k+1)}(t)\\\
=&\begin{pmatrix}\bm{0}&-\bm{B}_{2}\\\
-\bm{B}_{3}&\bm{0}\end{pmatrix}\dot{\bm{u}}^{(k)}(t)+\begin{pmatrix}\bm{0}&-\bm{A}_{2}\\\
-\bm{A}_{3}&\bm{0}\end{pmatrix}\bm{u}^{(k)}(t)+\bm{f}(t),\end{split}$

with $t\in[0,T_{f}]$ and $\bm{u}^{(k+1)}(0)=\bm{u}_{0}$. The inherent
parallelism of this method is reflected by the block-diagonal structure of the
matrices on the left-hand side.

One determines the convergence properties of continuous WR methods by
analyzing the iteration (12). Similarly, time-discrete WR methods are
described via time-discretizations of (12).

###  6.1 Time-discrete WR with asynchronous communication

We now consider the time-discrete case. Similar to [10, 25], we use convergent
and zero-stable $m$-step linear multistep methods (LMM), see e.g. [8,
Chap.3.2], on matching time-grids with the constant step-size $\Delta
t=T_{f}/N$. A $m$-step LMM applied to (10) is

$\sum_{\ell=0}^{m}\left(a_{\ell}\bm{B}+b_{\ell}\Delta
t\bm{A}\right)\bm{u}_{n+\ell}=\Delta
t\sum_{\ell=0}^{m}b_{\ell}\,\bm{f}(t_{n+\ell}),\quad n=m\ldots,N,$ (13)

which defines the discrete monolithic solution. A particular LMM is defined by
its coefficients $a_{\ell},b_{\ell}$ and requires starting values
$\bm{u}_{\ell}$, $\ell=0,\ldots,m-1$.

Classical WR methods such as GS and Jacobi WR can be described by time-
discretizations of (12). Using the same LMM, this is

$\displaystyle\begin{split}\sum_{\ell=0}^{m}\left(a_{\ell}\bm{M}_{B}+b_{\ell}\Delta
t\,\bm{M}_{A}\right)\bm{u}_{n+\ell}^{(k+1)}=\sum_{\ell=0}^{m}\left(a_{\ell}\bm{N}_{B}+b_{\ell}\Delta
t\,\bm{N}_{A}\right)\bm{u}^{(k)}_{n+\ell}+b_{\ell}\Delta
t\,\bm{f}(t_{n+\ell}),\quad n=m\ldots,N,\end{split}$

with starting values $\bm{u}_{\ell}^{(k)}=\bm{u}_{\ell}$, $\ell=0,\ldots,m-1$,
$\forall k\geq 0$.

We can describe Algorithm 2 by an analogous iteration, where the splitting
matrices (11) can differ for each $\ell$, $n$ and $k$. That is, we consider

$\displaystyle\begin{split}&\sum_{\ell=0}^{m}\left(a_{\ell}\bm{M}_{B,n,\ell}^{(k+1)}+b_{\ell}\Delta
t\,\bm{M}^{(k+1)}_{A,n,\ell}\right)\bm{u}_{n+\ell}^{(k+1)}\\\
=&\sum_{\ell=0}^{m}\left(a_{\ell}\bm{N}_{B,n,\ell}^{(k+1)}+b_{\ell}\Delta
t\,\bm{N}^{(k+1)}_{A,n,\ell}\right)\bm{u}^{(k)}_{n+\ell}+b_{\ell}\Delta
t\,\bm{f}(t_{n+\ell}),\quad n=m\ldots,N.\end{split}$ (14)

Here, the concrete matrices $\bm{M}_{B,n,\ell}^{(k+1)}$ and
$\bm{M}_{A,n,\ell}^{(k+1)}$ are, for each $\ell$, $n$ and $k$, determined by
$\underline{\bm{\bm{v}}}_{*}^{(k)}$ and $\underline{\bm{\bm{w}}}_{*}^{(k)}$,
as emerging from the realized communication in e.g., Algorithm 2, including
relaxation. These matrices fulfill the splitting property
$\bm{B}=\bm{M}_{B,n,\ell}^{(k+1)}-\bm{N}_{B,n,\ell}^{(k+1)}$, $\bm{A}$
analogous, by which the discrete monolithic solution defined by (13) is a
fixed point of the discrete WR method defined by (14).

We define the discrete WR error as

$\bm{e}_{n}^{(k)}:=\bm{u}_{n}-\bm{u}_{n}^{(k)}.$

Taking the difference between (13) and (14) shows that it fulfills

$\sum_{\ell=0}^{m}\bm{C}^{(k+1)}_{n,\ell}\bm{e}_{n+\ell}^{(k+1)}=\sum_{\ell=0}^{m}\bm{D}^{(k+1)}_{n,\ell}\bm{e}^{(k)}_{n+\ell},$
(15)

with

$\bm{C}^{(k)}_{n,\ell}:=a_{\ell}\bm{M}_{B,n,\ell}^{(k)}+b_{\ell}\,\Delta
t\bm{M}^{(k)}_{A,n,\ell},\quad\bm{D}^{(k)}_{n,\ell}:=a_{\ell}\bm{N}_{B,n,\ell}^{(k)}+b_{\ell}\,\Delta
t\bm{N}^{(k)}_{A,n,\ell}.$

The starting values $\bm{u}^{(k+1)}_{\ell}$ define the starting errors
$\bm{e}^{(k+1)}_{\ell}$, $\ell=0,\ldots,m-1$.

In the following theorem we show convergence in the form of
$\|\bm{e}^{(k)}_{n}\|\rightarrow 0$, for $k\rightarrow\infty$, for all
$n=m,\ldots,N$. It is an extension of the convergence result from [10]. There,
we have constant splittings, whereas in our method the splittings vary with
$n$, $\ell$ and $k$. Here $\|\cdot\|:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a
norm and we similar use $\|\cdot\|$ to denote the induced matrix norm.

######  Theorem 6.1.

Let the splittings

$\bm{B}=\bm{M}_{B,n,\ell}^{(k)}-\bm{N}_{B,n,\ell}^{(k)}\quad\text{and}\quad\bm{A}=\bm{M}_{A,n,\ell}^{(k)}-\bm{N}_{A,n,\ell}^{(k)},\quad\bm{M}_{B,n,\ell}^{(k)}\text{
nonsingular},$

with

$\bm{C}_{n,m}^{(k)}\text{
nonsingular},\quad\|{\bm{C}_{n,m}^{(k)}}^{-1}\bm{D}_{n,m}^{(k)}\|<1,\quad
n=m,\ldots,N$ (16)

and an initial guess

$\underline{\bm{e}}^{(0)}:=\left({\bm{e}_{m}^{(0)}}^{T},\ldots,{\bm{e}_{N}^{(0)}}^{T}\right)^{T}\in\mathbb{R}^{d(N-m+1)}$

be given. Then, the solution of the discrete WR method defined by (14)
converges to the solution of (13).

###### Proof.

We consider the so called ”all-at-once system” system, which is the system for
all timesteps of a given iteration. This is

$\underline{\bm{C}}^{(k+1)}\underline{\bm{e}}^{(k+1)}=\underline{\bm{D}}^{(k+1)}\underline{\bm{e}}^{(k)},$

with

$\displaystyle\underline{\bm{C}}^{(k)}$
$\displaystyle:=\begin{pmatrix}\bm{C}^{(k)}_{m,m}&\bm{0}&\ldots&\ldots&\ldots&\bm{0}\\\
\vdots&\bm{C}^{(k)}_{m+1,m}&\ddots&&&\vdots\\\
\bm{C}^{(k)}_{m^{2},0}&&\ddots&\ddots&&\vdots\\\
\bm{0}&\ddots&&\ddots&\ddots&\vdots\\\ \vdots&\ddots&\ddots&&\ddots&\bm{0}\\\
\bm{0}&\ldots&\bm{0}&\bm{C}^{(k)}_{N,0}&\ldots&\bm{C}^{(k)}_{N,m}\end{pmatrix},$
$\displaystyle\underline{\bm{D}}^{(k)}$
$\displaystyle:=\begin{pmatrix}\bm{D}^{(k)}_{m,m}&\bm{0}&\ldots&\ldots&\ldots&\bm{0}\\\
\vdots&\bm{D}^{(k)}_{m+1,m}&\ddots&&&\vdots\\\
\bm{D}^{(k)}_{m^{2},0}&&\ddots&\ddots&&\vdots\\\
\bm{0}&\ddots&&\ddots&\ddots&\vdots\\\ \vdots&\ddots&\ddots&&\ddots&\bm{0}\\\
\bm{0}&\ldots&\bm{0}&\bm{D}^{(k)}_{N,0}&\ldots&\bm{D}^{(k)}_{N,m}\end{pmatrix},$

with $\underline{\bm{C}}$,
$\underline{\bm{D}}\in\mathbb{R}^{(d(N-m+1))\times(d(N-m+1))}$. Using standard
index notation for $\underline{\bm{C}}^{(k)}$ to reference the above blocks,
its inverse is

$\displaystyle{\underline{\bm{C}}^{(k)}}^{-1}=\begin{pmatrix}\bm{G}^{(k)}_{1,1}&\bm{0}&\ldots&\bm{0}\\\
\vdots&\ddots&\ddots&\vdots\\\ \vdots&&\ddots&\bm{0}\\\
\bm{G}^{(k)}_{N-m+1,1}&\ldots&\ldots&\bm{G}^{(k)}_{N-m+1,N-m+1}\end{pmatrix},$
$\displaystyle\bm{G}^{(k)}_{i,j}=\begin{cases}\left(\underline{\bm{C}}_{i,j}^{(k)}\right)^{-1},&\quad
i=j\\\
-\left(\underline{\bm{C}}_{i,i}^{(k)}\right)^{-1}\sum_{\ell=1}^{i-1}\underline{\bm{C}}^{(k)}_{i,\ell}\bm{G}_{\ell,j},&\quad
i\neq j\end{cases}.$

The inverse ${\underline{\bm{C}}^{(k)}}^{-1}$ only requires inverses of its
diagonal blocks $\bm{C}^{(k)}_{n,m}$, for which we assume existence, see (16).
Since the resulting iteration matrix
${\underline{\bm{C}}^{(k)}}^{-1}\underline{\bm{D}}^{(k)}$ is block lower-
triangular, we get the following forward elimination:

$\displaystyle\begin{split}\bm{e}_{m}^{(k+1)}&={\bm{C}_{m,m}^{(k+1)}}^{-1}\bm{D}_{m,m}^{(k+1)}\bm{e}_{m}^{(k)},\\\
\bm{e}_{m+1}^{(k+1)}&={\bm{C}_{m+1,m}^{(k+1)}}^{-1}\bm{D}_{m+1,m}^{(k+1)}\bm{e}_{m+1}^{(k)}\\\
&\quad+\left(\bm{G}_{2,1}^{(k+1)}\bm{D}_{m,m}^{(k+1)}+{\bm{C}_{m+1,m}^{(k+1)}}^{-1}\bm{D}_{m+1,m-1}^{(k+1)}\right)\bm{e}_{m}^{(k)},\\\
\bm{e}_{m+2}^{(k+1)}&={\bm{C}_{m+2,m}^{(k+1)}}^{-1}\bm{D}_{m+2,m}^{(k+1)}\bm{e}_{m+2}^{(k)}\\\
&\quad+\left(\bm{G}_{3,2}^{(k+1)}\bm{D}_{m+1,m-1}^{(k+1)}+{\bm{C}_{m+2,m}^{(k+1)}}^{-1}\bm{D}_{m+2,m-1}^{(k+1)}\right)\bm{e}_{m+1}^{(k)}\\\
&\quad+\left(\bm{G}_{3,1}^{(k+1)}\bm{D}_{m,m}^{(k+1)}+\bm{G}_{3,2}^{(k+1)}\bm{D}_{m+1,m-1}^{(k+1)}+{\bm{C}_{m+2,m}^{(k+1)}}^{-1}\bm{D}_{m+2,m-1}^{(k+1)}\right)\bm{e}_{m}^{(k)},\\\
\bm{e}_{m+3}^{(k+1)}&=\ldots.\end{split}$ (17)

We get $\|\bm{e}_{m}^{(k+1)}\|=0$, $k\rightarrow\infty$ from
$\|{\bm{C}_{m,m}^{(k+1)}}^{-1}\bm{D}_{m,m}^{(k+1)}\|<1$. With
$\bm{e}_{m}^{(k+1)}$ vanishing and
$\|{\bm{C}_{m+1,m}^{(k+1)}}^{-1}\bm{D}_{m+1,m}^{(k+1)}\|<1$, we get
$\|\bm{e}_{m+1}^{(k+1)}\|=0$, $k\rightarrow\infty$. By induction we get
$\|\bm{e}_{n}^{(k+1)}\|=0$, $k\rightarrow\infty$, for all $n=m,\ldots,N$,
which implies $\|\underline{\bm{e}}^{(k+1)}\|=0$, $k\rightarrow\infty$. ∎

The assumption of $\bm{C}_{n,m}^{(k)}$ nonsingular in (16) is a solvability
assumption on the occurring linear systems in (15). The growing number of
terms on the right hand sides in (17) shows the potential for initial growth
of $\|\underline{\bm{e}}^{(k)}\|$, even if the iteration matrix may be normal.

######  Remark 6.2.

Consider Algorithm 2 with matching time-grids, i.e., $N=N_{v}=N_{w}$ and
constant relaxation. Then, the splitting matrices ${\bm{M}_{B,n,\ell}^{(k)}}$
resp. $\bm{M}_{A,n,\ell}^{(k)}$ match either those of Jacobi or GS WR for each
$\ell,n,k$, due to matching $(\underline{\bm{\bm{v}}}_{*}^{(k)})_{n}$ and
$(\underline{\bm{\bm{w}}}_{*}^{(k)})_{n}$, for each $n$, see (8). Thus, one
only needs to consider three distinct cases for the matrices (16). Each one of
these is a convergence requirement for either Jacobi or GS WR. Consequently,
time-discrete convergence of Jacobi and GS WR (in either ordering of (1))
means (16) is met.

###  6.2 Continuous WR with asynchronous communication

We now consider continuous WR methods with $\bm{w}_{\,*}^{(k)}$ and
$\bm{v}_{\,*}^{(k)}$ in (2a), (2c) varying with $t$ and $k$. By straight-
forward substitutions of $\hat{\bm{v}}^{(k+1)}$ and $\hat{\bm{w}}^{(k+1)}$ in
the variable relaxation steps (9) into (2), we get the iteration

$\displaystyle\begin{split}&\bm{M}_{B}^{(k+1)}(t)\dot{\bm{u}}^{(k+1)}(t)+\bm{M}_{A}^{(k+1)}(t)\bm{u}^{(k+1)}(t)\\\
=&\,\,\bm{N}_{B}^{(k+1)}(t)\dot{\bm{u}}^{(k)}(t)+\bm{N}_{A}^{(k+1)}(t)\bm{u}^{(k)}(t)+\bm{f}(t),\quad\bm{u}^{(k+1)}(0)=\bm{u}_{0},\quad
t\in[0,T_{f}],\end{split}$ (18)

with splittings

$\bm{B}=\bm{M}_{B}^{(k)}(t)-\bm{N}_{B}^{(k)}(t),\quad\bm{A}=\bm{M}_{A}^{(k)}(t)-\bm{N}_{A}^{(k)}(t),\quad\bm{M}_{B}^{(k)}(t)\text{
nonsingular},$ (19)

$t\in[0,T_{f}]$, $k>0$. Here, we omit the dependencies on the relaxation
matrices for readability.

In the previous section, we considered the convergence of the fully discrete
WR iteration (14) for $\Delta t$ fixed for $k\rightarrow\infty$, which gives
(13). Here, we instead discuss the convergence of the continuous iteration
(18) for $k\rightarrow\infty$. Before doing so, we would like to point out
that we cannot guarantee that we obtain (18) from (14) in the limit $\Delta
t\rightarrow 0$. The reason is that the splittings chosen and thus the
matrices $\bm{M}_{B,n,\ell}^{(k+1)}$ resp. $\bm{M}_{A,n,\ell}^{(k+1)}$ in (14)
can change with every time step. This would yield $\bm{M}_{B}^{(k)}$,
$\bm{M}_{A}^{(k)}$ discontinuous everywhere in the limit and (18) would not be
well defined.

The typical scenario for (14), as implemented via Algorithm 2, is that
splittings match those of Jacobi WR until one subsolver is at least one
timestep ahead of another subsolver. From then on, the splitting matrices
match those of GS WR. I.e., for a given $k$, $\bm{M}_{B}^{(k)}$ and
$\bm{M}_{A}^{(k)}$ are piece-wise constant, with a single discontinuity.

We thus assume that the limit has only a finite number of jumps and consider
(18) in a piece-wise sense with piece-wise Lipschitz-continuous data. This
guarantees existence of a piece-wise solution of (18) for all $k>0$.
Additionally, we assume that splitting matrices corresponding to the same
time-point, in the same iteration, are identical. E.g., in (14) both
$\bm{M}_{B,n,\ell}^{(k+1)}$ and $\bm{M}_{B,n+1,\ell-1}^{(k+1)}$ correspond to
$t_{n+\ell}$. This can be guaranteed by implementation, storing interpolant
evaluations. Now we analyze the convergence properties of (18) under these
assumptions.

Consider (10) with $\bm{A}$, $\bm{B}$, $\bm{f}$ time-dependent and piece-wise
Lipschitz-continuous. Then the solution is

$\bm{u}(t)=\text{e}^{-C(t)}\left(\bm{u}_{0}+\int_{0}^{t}\text{e}^{C(s)}\bm{B}^{-1}(s)\bm{f}(s)\text{d}s\right),$
(20)

where

$C(t)=\int_{0}^{t}\bm{B}^{-1}(s)\bm{A}(s)\text{d}s.$

We can apply this solution formula to (18). Replacing $\dot{\bm{u}}^{(k)}$ via
integration by parts and performing lengthy, but straight-forward
rearrangements, yield the solution:

$\bm{u}^{(k+1)}(t)=\bm{K}^{(k+1)}(t)\bm{u}^{(k)}(t)+\int_{0}^{t}\bm{\mathcal{K}}_{c}^{(k+1)}(s)\bm{u}^{(k)}(s)\text{d}s+\bm{\varphi}^{(k+1)}(t),$
(21)

with

$\displaystyle\begin{split}\bm{K}^{(k)}(t)&={\bm{M}_{B}^{(k)}}^{-1}(t)\bm{N}_{B}^{(k)}(t),\\\
\bm{C}^{(k)}(t)&=\int_{0}^{t}{\bm{M}_{B}^{(k)}}^{-1}(s)\bm{M}_{A}^{(k)}(s)\text{d}s,\\\
\bm{\mathcal{K}}_{c}^{(k)}(t)&=\text{e}^{\bm{C}^{(k)}(s)-\bm{C}^{(k)}(t)}\left({\bm{M}_{B}^{(k)}}^{-1}(s)\bm{N}_{A}(s)-\frac{\text{d}}{\text{d}s}\left(\text{e}^{\bm{C}^{(k)}(s)}\right)\bm{K}^{(k)}(s)-\frac{\text{d}}{\text{d}s}\bm{K}^{(k)}(s)\right),\\\
\bm{\varphi}^{(k)}(t)&=\text{e}^{-\bm{C}^{(k)}(t)}\left(\left(\bm{I}-\bm{K}^{(k)}(0)\right)\bm{u}_{0}+\int_{0}^{t}\text{e}^{\bm{C}^{(k)}(s)}{\bm{M}_{B}^{(k)}}^{-1}(s)\bm{f}(s)\text{d}s\right),\end{split}$
(22)

where

$\frac{\text{d}}{\text{d}s}\left(\text{e}^{\bm{C}^{(k)}(s)}\right)=\int_{0}^{1}\text{e}^{\alpha\bm{C}^{(k)}(s)}\frac{\text{d}\bm{C}^{(k)}(s)}{\text{d}s}\text{e}^{(1-\alpha)\bm{C}^{(k)}(s)}\text{d}\alpha,$

c.f. [32].

Consider the continuous WR error

$\bm{e}^{(k)}:=\bm{u}^{(k)}-\bm{u},$ (23)

where $\bm{u}$ is the solution to (10). Taking the difference between (21) and
(20) gives:

$\bm{e}^{(k+1)}(t)=\bm{K}^{(k+1)}(t)\bm{e}^{(k)}(t)+\int_{0}^{t}\bm{\mathcal{K}}_{c}^{(k+1)}(s)\bm{e}^{(k)}(s)\text{d}s,\quad\bm{e}^{(k+1)}(0)=\bm{0},\quad
t\in[0,T_{f}].$

For the following theorem, we define the (vector) function norm
$\|\bm{e}\|_{[0,t]}:=\sup_{\tau\in[0,t]}\|\bm{e}(\tau)\|$ and (matrix)
function norm $\|\bm{A}\|_{[0,t]}:=\sup_{\tau\in[0,t]}\|\bm{A}(\tau)\|$, based
on the induced matrix norm. This result is an extension of a convergence
result in [9], to the situation where splittings vary with $t$ and $k$.

######  Theorem 6.3.

Let splittings (19) with $\bm{M}_{B}^{(k)}$, $\bm{M}_{A}^{(k)}$ and
$\bm{e}^{(0)}$ piece-wise Lipschitz-continuous for all $k>0$ be given. Then,
the error (23) fulfills

$\|\bm{e}^{(k)}\|_{[0,t]}\leq\left(\sum_{j=0}^{k}{k\choose
j}{K^{\,\max}(t)}^{k-j}{\mathcal{K}_{c}^{\,\max}(t)}^{j}\frac{t^{\,j}}{j!}\right)\|\bm{e}^{(0)}\|_{[0,t]},$

where $K^{\,\max}(t):=\sup_{k\in\mathbb{N}}\|\bm{K}^{(k)}\|_{[0,t]}$ and
$\mathcal{K}_{c}^{\,\max}(t):=\sup_{k\in\mathbb{N}}\|\bm{\mathcal{K}}_{c}^{(k)}\|_{[0,t]}$,
c.f. (22).

###### Proof.

To avoid ambiguity for the function norm, we denote the relevant variable by
$\tau$. Following the same principles as in the constant splitting case [14,
9], applications of the triangle inequality, submultiplicativity and straight-
forward upper bounds yield

$\displaystyle\|\bm{e}^{(k)}\|_{[0,t]}$
$\displaystyle=\left\|\bm{K}^{(k)}\bm{e}^{(k-1)}+\int_{0}^{\tau}\bm{\mathcal{K}}_{c}^{(k)}(s)\bm{e}^{(k-1)}(s)\text{d}s\right\|_{[0,t]}$
$\displaystyle\leq\big{\|}\bm{K}^{(k)}\bm{e}^{(k-1)}\big{\|}_{[0,t]}+\int_{0}^{t}\big{\|}\bm{\mathcal{K}}_{c}^{(k)}(s)\bm{e}^{(k-1)}(s)\big{\|}\text{d}s$
$\displaystyle\leq\big{\|}\bm{K}^{(k)}\big{\|}_{[0,t]}\big{\|}\bm{e}^{(k-1)}\big{\|}_{[0,t]}+\big{\|}\bm{\mathcal{K}}_{c}^{(k)}\big{\|}_{[0,t]}\int_{0}^{t}\big{\|}\bm{e}^{(k-1)}\big{\|}_{[0,t]}\text{d}s.$

Repeated application and taking the supremum over $k$ then gives

$\displaystyle\|\bm{e}^{(k)}\|_{[0,t]}$
$\displaystyle\leq\left(\sum_{j=0}^{k}{k\choose
j}K^{\,\max}(t)^{k-j}\mathcal{K}_{c}^{\,\max}(t)^{j}\int\ldots\int
1\,\,\text{d}s_{1}\ldots\text{d}s_{j}\right)\|\bm{e}^{(0)}\|_{[0,t]}$
$\displaystyle\leq\left(\sum_{j=0}^{k}{k\choose
j}K^{\,\max}(t)^{k-j}\mathcal{K}_{c}^{\,\max}(t)^{\,j}\frac{t^{\,j}}{j\,!}\right)\|\bm{e}^{(0)}\|_{[0,t]}.$

∎

Our results is consistent with the time-discrete result of Theorem 6.1 for
$\Delta t\rightarrow 0$, under the aforementioned assumptions. The
$\mathcal{K}_{c}^{\,\max}(t)^{\,j}\frac{t^{\,j}}{j\,!}$ term converges super-
linearly for $j\rightarrow\infty$, but can lead to large error bounds for
small $j$ and large $t$. The asymptotic convergence rate for
$k\rightarrow\infty$ is bounded from above by $\|K^{\,\max}\|_{[0,t]}$.

##  7 Variable relaxation algorithm for two coupled problems

We now provide an algorithm and an implementation for variable relaxation when
using asynchronous communication. Theorem 6.1, which includes variable
relaxation, shows that the discrete asymptotic convergence rate is bounded by

$\displaystyle\begin{split}\max_{n=m,\ldots,N,\,\,k>0}\left\|{\bm{C}_{n,m}^{(k)}}^{-1}\bm{D}_{n,m}^{(k)}\right\|.\end{split}$
(24)

These ${\bm{C}_{n,m}^{(k)}}^{-1}\bm{D}_{n,m}^{(k)}$ are the diagonal blocks of
the iteration matrix and they depend on the chosen relaxation. Thus, we choose
relaxation to minimize the spectral radii, resp. norms of all diagonal blocks,
which minimizes (24). The optimal relaxation depends on the problem and the
splitting. We discuss how to determine the specific values for our numerical
experiments in Section 8.1.2.

Here, we present an algorithm for the separate processes to determine at
runtime which splitting occurs in each timestep. With two coupled problems,
there are exactly three cases, corresponding to Jacobi and GS WR, see Remark
6.2. This algorithm does not have a straight-forward extension to more than
two coupled problems. We consider relaxation of the subset of unknowns
exchanged between the processors. Furthermore, we consider the more general
case of non-matching time-grids with constant stepsizes.

In the following we reference the subsolver for a given subproblem as process,
not excluding usage of multiple processors to solve a subproblem. The basic
structure for the WR iteration and time-integration are analogous to Algorithm
2. The differences will be within interpolation, communication and relaxation.

######  Definition 7.1.

In Algorithm 2, we say a process is ahead of another process, if all
interpolant evaluations in the $n$-th timestep and $k$-th iteration depend on
data-points from the $k$-th iteration, rather than the $(k-1)$-st iteration.

######  Definition 7.2.

We say a timestep has a local Jacobi shape, if no process is ahead of another
process, and a local GS shape, if one process is ahead of the other.

Our new algorithm for variable relaxation consists of first determining the
local shape and then updating the interpolant, using appropriate relaxation.
Thus, we first communicate $\hat{\bm{v}}^{(k+1)}_{n}$ and
$\hat{\bm{w}}^{(k+1)}_{n}$ to a buffer on the respective other process. Then,
that process determines the local shape and updates the interpolant data
$\underline{\bm{\bm{v}}}_{\,*}$ resp. $\underline{\bm{\bm{w}}}_{\,*}$ using
appropriate relaxation. With relaxation independent of the local shapes, e.g.,
Algorithm 2 or constant relaxation, buffers are not required and one directly
updates $\underline{\bm{\bm{v}}}_{\,*}$ and $\underline{\bm{\bm{w}}}_{\,*}$.

By comparison with the constant splittings we see the following: Jacobi WR
requires no relaxation during time-integration, see (4), since there is no
dependency on the new iterate. With GS WR, see (3), one only requires
relaxation of the incoming data for the process that is not ahead. Thus, we
only update the interpolant during time-integration, if the other process is
ahead. We determine if a process is ahead, based on additionally communicated
update indicators. For each data-point, these indicate if it has been remotely
updated. We explain the method to determine if the other process is ahead
later on.

The principle data-structures involved are shown in Figure 1. We denote the
time-grids by $\mathcal{T}_{v}=\\{t^{(v)}_{n}\\}_{n}$ and
$\mathcal{T}_{w}=\\{t^{(w)}_{n}\\}_{n}$, with
$\mathcal{T}=\mathcal{T}_{v}\cup\mathcal{T}_{w}$.

The local shape and thus resulting relaxation is determined by a single
process and thus defined on the time-points of a single time-grid. To ensure
consistent relaxation for both interpolants, we define them on the shared
time-grid $\mathcal{T}$, rather than the time-grid of the respective other
process, as in Algorithm 2. The base interpolant data is determined by the
time-grid of the respective other process, which we interpolate to the shared
grid before relaxation.

MPI_Window: update ind. ($\mathcal{T}_{w}$)MPI_PutMPI_Window: buffer
$\hat{\underline{\bm{\bm{w}}}}^{(k+1)}$
($\mathcal{T}_{w}$)MPI_Put$\mathcal{I}$$\mathcal{I}$Interpolant data
$\underline{\bm{\bm{w}}}_{\,*}$ ($\mathcal{T}$)RelaxSolution:
$\underline{\bm{\hat{\bm{v}}}}^{(k+1)}$ ($\mathcal{T}_{v}$)time-int.MPI_Put
Figure 1: Sketch of the data-structures involved for variable relaxation from
the perspective of the process solving (6a). Iteration indicators and the
buffer are MPI_Window objects to receive asynchronous updates via MPI_Put from
the other process. Their size is determined by the time-grid used on the other
process. When performing relaxation, we update those data-points of the
interpolant that are part of $\mathcal{T}_{w}$. Any remaining points are
computed via interpolation. New time-point solutions are asynchronously
communicated to the buffer of the other process, also updating the
corresponding update indicators.

W.l.o.g., consider the timestep from $t_{n}^{(v)}$ to $t_{n+1}^{(v)}$. To
determine if the other process is ahead, we first need to find the smallest
enclosing interval $[t_{-}^{(w)},t_{+}^{(w)}]$ with $t_{-}^{(w)}$,
$t_{+}^{(w)}$ $\in\mathcal{T}_{w}$, such that all evaluations of the
interpolant during this timestep depend on discrete data points
$\bm{w}^{(k)}_{n}$ or $\bm{w}^{(k+1)}_{n}$ corresponding to the time-points
$t_{-}^{(w)},\ldots,t_{+}^{(w)}$. We determine if the other process is ahead
by checking the update indicator at $t_{+}^{(w)}$. An example for determining
the smallest enclosing interval, on non-matching time-grids, is shown in
Figure 2.

$\mathcal{T}_{v}$$t_{n}^{(v)}$$t_{n+1}^{(v)}$$\mathcal{T}_{w}$$t_{-}^{(w)}$$t_{+}^{(w)}$$\times$
Figure 2:  Consider linear interpolation and a Runge-Kutta scheme with
$c_{i}\in[0,1]$, then $t_{-}^{(w)}=\max(\\{t\in\mathcal{T}_{w}|t\leq
t_{n}^{(v)}\\})$ and $t_{+}^{(w)}=\min(\\{t\in\mathcal{T}_{w}|t\geq
t_{n+1}^{(v)}\\})$. The figure visualizes the smallest enclosing interval on
non-matching $\mathcal{T}_{v}$ and $\mathcal{T}_{w}$. The cross marks the
time-point for which we need to check the update indicator, to determine if
the other process is ahead.

We use markings for all discrete timepoints $t\in\mathcal{T}$, tracking
relaxation type and if relaxation has been performed. All timepoints are
unmarked at the beginning of each iteration.

In a given timestep, we mark all unmarked timepoints within the smallest
enclosing interval for appropriate GS relaxation, if the other process is
ahead, and for Jacobi relaxation otherwise. Additionally, if the other process
is ahead, we perform relaxation on all not previously relaxed data-points for
$t\in\mathcal{T}$: $t_{-}^{(w)}\leq t\leq t_{+}^{(w)}$, according to their
respective markings. Finally, the actual timestep is computed. This procedure
is visualized in Algorithm 3.

timestep $n$Sync. buffer + update indicatorsFind enclosing intervalother
process ahead?mark points (GS)mark points (Jacobi)relax data-pointscompute
stepSend dataYesNo$n=n+1$ List of algfloats 3 Timestepping procedure for
variable relaxation algorithm with two processes. The synchronization of
buffer and update indicators is performed using MPI_Win_sync.

After time-integration, the processes exchange information on which
$t\in\mathcal{T}$ received GS relaxation and perform according relaxation on
all non-relaxed points (over-ruling any markings for Jacobi relaxation).
Afterwards, any remaining Jacobi relaxation is performed. Lastly, the
solutions corresponding to $\bm{v}^{(k+1)}(T_{f})$ and $\bm{w}^{(k+1)}(T_{f})$
are exchanged to facilitate a consistent termination check, c.f. (5), on both
processes.

##  8 Numerical results

We consider two conjugate heat transfer examples. The first one is two
heterogeneous coupled linear heat equations, which is conform with the linear
convergence theory presented in Section 6. We use it to demonstrate
convergence of our new method and for a performance comparison with Jacobi and
GS WR.

The second example is a gas quenching test case, inspired by [2, 3],
simulating cooling of a hot steel plate with pressurized air. We model the air
using the compressible Euler equations and the steel plate via the nonlinear
heat equation. We use a partitioned coupling of different spatial
discretizations implemented using the packages DUNE [1] and FEniCS [16]. This
example is not conform with the assumptions on Section 6, since it is
nonlinear and we use non-matching time-grids. Yet, we demonstrate that our new
method is convergent.

###  8.1 Coupled heat equations

The model equations are

$\displaystyle\begin{split}\alpha_{1}\partial_{t}u_{1}(t,\bm{x})-\lambda_{1}\Delta
u_{1}(t,\bm{x})=0,&\quad(t,\bm{x})\in(0,T_{f}]\times\Omega_{1},\\\
u_{1}(t,\bm{x})=0,&\quad(t,\bm{x})\in[0,T_{f}]\times\Omega_{1}\setminus\Gamma,\\\
u_{1}(t,\bm{x})=u_{\Gamma}(t,\bm{x}),&\quad(t,\bm{x})\in[0,T_{f}]\times\Gamma,\\\
u_{1}(0,\bm{x})=u_{0}(\bm{x}),&\quad\bm{x}\in\Omega_{1},\end{split}$ (25)

and

$\displaystyle\begin{split}\alpha_{2}\partial_{t}u_{2}(t,\bm{x})-\lambda_{2}\Delta
u_{2}(t,\bm{x})=0,&\quad(t,\bm{x})\in(0,T_{f}]\times\Omega_{2},\\\
u_{2}(t,\bm{x})=0,&\quad(t,\bm{x})\in[0,T_{f}]\times\Omega_{2}\setminus\Gamma,\\\
\lambda_{2}\nabla u_{2}(t,\bm{x})\cdot\bm{n}_{2}=-\lambda_{1}\nabla
u_{1}(t,\bm{x})\cdot\bm{n}_{1},&\quad(t,\bm{x})\in[0,T_{f}]\times\Gamma,\\\
u_{2}(0,\bm{x})=u_{0}(\bm{x}),&\quad\bm{x}\in\Omega_{2}.\end{split}$ (26)

Here, $\lambda$ is the thermal conductivity and the thermal diffusivity $D$ is
defined by

$\displaystyle D=\lambda/\alpha,\quad\mbox{with}\quad\alpha=\rho c_{p},$

with density $\rho$ and specific heat capacity $c_{p}$.

The corresponding monolithic problem is a linear heat equation defined on
$\Omega_{1}\cup\Omega_{2}$ with space dependent material parameters that have
a jump in $\alpha(\bm{x})$ and $\lambda(\bm{x})$ at the interface. By
enforcing continuity of temperature and heat flux at the interface, the above
partitioned formulation is equivalent to the monolithic problem, in a weak
sense [26, Chap.7].

$0$$1$$-1$$0$$1$$\Omega_{1}$$\Omega_{2}$$\bm{n}_{1}$$\bm{n}_{2}$$\Gamma$
Figure 3: Geometry of the coupled heat problem.

We consider $T_{f}=10^{4}$ and the initial condition
$u_{0}(\bm{x})=500\sin(\pi/2(x_{1}+1))\sin(\pi x_{2})$. See Table 1 for the
materials considered here and Figure 3 for the geometry.

Material | $\alpha=\rho\cdot c_{p}[J/(Km^{3})]$ | $\lambda[W/(mK)]$
---|---|---
Air | $1.293\cdot 1005$ | $0.0243$
Water | $999.7\cdot 4192.1$ | $0.58\par$
Steel | $7836\cdot 443$ | $48.9$
Table 1: Material parameters.

The subproblems exchange information in the form of interface temperature
$u_{\Gamma}=u_{2}\big{|}_{\Gamma}$ and the heat flux $q:=\lambda_{1}\nabla
u_{1}\cdot\bm{n}_{1}$.

####  8.1.1 Discretizations

We discretize (26) in space using linear finite elements implemented using
FEniCS [16], on a triangulation obtained from a cartesian grid. We use the
Crank-Nicolson method to discretize time. This yields

$\displaystyle\begin{split}&\int_{\Omega_{2}}\alpha_{2}(u^{(2)}_{n+1}-u^{(2)}_{n})\varphi+\frac{\Delta
t}{2}\lambda_{2}\nabla\left(u^{(2)}_{n+1}+u^{(2)}_{n}\right)\cdot\nabla\varphi\,\text{d}\bm{x}\\\
&+\frac{\Delta
t}{2}\int_{\Gamma}\left(q_{n+1}+q_{n}\right)\varphi\,\text{d}\bm{S}=0,\quad\forall\varphi\in
V,\end{split}$

with an appropriate finite element space $V$. The weak form of (25) is
obtained by omitting the heat flux $q$ and including the Dirichlet-boundary
condition at $\Gamma$. On $\Omega_{1}$ we use the above weak form to compute
the heat flux $q_{n+1}\approx q(t_{n+1})$, based on $u^{(1)}_{n+1}$,
$u^{(1)}_{n}$ and $q_{n}$. We compute the initial flux $q_{0}$ required for
the interpolant from the initial condition via

$q_{0}=\lambda_{1}\int_{\Gamma}(\nabla
u_{0}\cdot\bm{n}_{1})\varphi\,\text{d}S.$

The WR methods are implemented following the partitioned approach, treating
the space-discretizations of the subsolvers as black-boxes. That is, instead
of exchanging discrete interface unknowns, the subsolvers exchange the
interface temperatures resp. heat fluxes corresponding to points at the
interface $\Gamma$.

####  8.1.2 Relaxation

We use single parameter relaxation on either one or both of the discrete
exchange variables $\bm{u}_{\Gamma}$ and $\bm{q}$, i.e.,

$\displaystyle\bm{u}_{\Gamma}^{(k+1)}(t)$
$\displaystyle=\Theta(t)^{(k+1)}\hat{\bm{u}}_{\Gamma}^{(k+1)}+(1-\Theta(t)^{(k+1)})\bm{u}_{\Gamma}^{(k)},$
$\displaystyle\bm{q}^{(k+1)}(t)$
$\displaystyle=\Theta(t)^{(k+1)}\hat{\bm{q}}^{(k+1)}+(1-\Theta(t)^{(k+1)})\bm{q}^{(k)},$

with $\Theta(t)^{(k+1)}\in\mathbb{R}$. To perform the algorithm from Section
5.1, we require relaxation parameters for Jacobi and GS WR.

Optimal constant relaxation on $\bm{u}_{\Gamma}$ for time-discrete GS WR in
the $\Omega_{1}\rightarrow\Omega_{2}$ order has been determined in [22, 23].
This was done for 1D linear finite elements on a uniform space discretization
and implicit Euler for constant and matching step-sizes. However, results from
[19] show $\Theta_{opt}$ to be robust, working well in 2D and with the second
order in time SDIRK2 scheme.

Here, we use the results in [22] to determine optimal relaxation for the
remaining local shapes, i.e., Jacobi and GS in the
$\Omega_{2}\rightarrow\Omega_{1}$ order. In the
$\Omega_{1}\rightarrow\Omega_{2}$ order, the iteration, including relaxation,
for $\bm{u}_{\Gamma}$ and $\bm{q}$ is given by the following relations [22]:

$\begin{pmatrix}\bm{u}_{\Gamma}\\\
\bm{q}\end{pmatrix}^{(k+1)}=\begin{pmatrix}(1-\Theta)\bm{I}-\Theta{\bm{S}^{(2)}}^{-1}\bm{S}^{(1)}&\bm{0}\\\
\bm{0}&\bm{S}^{(1)}\end{pmatrix}\begin{pmatrix}\bm{u}_{\Gamma}\\\
\bm{q}\end{pmatrix}^{(k)}+\bm{\psi}^{(k)}.$

Here, $\bm{\psi}^{(k)}$ are additional terms irrelevant to the iteration
matrix. In the 1D case, where $\bm{S}^{(m)}\in\mathbb{R}$, the optimal choice
is $\Theta_{opt}=\frac{1}{\left|1+{\bm{S}^{(2)}}^{-1}\bm{S}^{(1)}\right|}$,
yielding a zero spectral radius for the above iteration matrix. Analytical
expressions for $\bm{S}^{(m)}$ have been computed in [22], which are

$\displaystyle\bm{S}^{(m)}$ $\displaystyle=\frac{6\Delta t\Delta
x(\alpha_{m}\Delta x^{2}+3\lambda_{m}\Delta t)-(\alpha_{m}\Delta
x^{2}-6\lambda\Delta t)^{2}s_{m}}{18\Delta t\Delta x^{3}},$ $\displaystyle
s_{m}$ $\displaystyle=\sum_{i=1}^{N}\frac{3\Delta t\Delta
x^{2}\sin^{2}(i\pi\Delta x)}{2\alpha_{m}\Delta x^{2}+6\lambda_{m}\Delta
t+(\alpha_{m}\Delta x^{2}-6\lambda_{m}\Delta t)\cos(i\pi\Delta x)}.$

Similarly, the relations for GS in the $\Omega_{2}\rightarrow\Omega_{1}$
order, with relaxation on $\bm{q}$, are

$\begin{pmatrix}\bm{u}_{\Gamma}\\\
\bm{q}\end{pmatrix}^{(k+1)}=\begin{pmatrix}\bm{0}&-{\bm{S}^{(2)}}^{-1}\\\
\bm{0}&(1-\Theta)\bm{I}-\Theta\bm{S}^{(1)}{\bm{S}^{(2)}}^{-1}&\end{pmatrix}\begin{pmatrix}\bm{u}_{\Gamma}\\\
\bm{q}\end{pmatrix}^{(k)}+\bm{\psi}^{(k)}.$

Relaxation with the same $\Theta_{opt}$ as in the
$\Omega_{1}\rightarrow\Omega_{2}$ case yields a zero spectral radius in 1D.
Jacobi WR with relaxation on both $\bm{u}_{\Gamma}$ and $\bm{q}$ using the
same $\Theta$ yields

$\begin{pmatrix}\bm{u}_{\Gamma}\\\
\bm{q}\end{pmatrix}^{(k+1)}=\begin{pmatrix}(1-\Theta)\bm{I}&-\Theta{\bm{S}^{(2)}}^{-1}\\\
\Theta\bm{S}^{(1)}&(1-\Theta)\bm{I}\end{pmatrix}\begin{pmatrix}\bm{u}_{\Gamma}\\\
\bm{q}\end{pmatrix}^{(k)}+\bm{\psi}^{(k)}.$

In 1D and with $\bm{S}^{(1)}{\bm{S}^{(2)}}^{-1}>0$, which is the case for all
material combinations here considered, the spectral radius is minimal for
$\Theta=1/(\bm{S}^{(1)}{\bm{S}^{(2)}}^{-1}+1)$. In particular, the spectral
radius with optimal relaxation is

$\rho_{opt}^{\text{Jacobi}}=\sqrt{\frac{\bm{S}^{(1)}{\bm{S}^{(2)}}^{-1}}{\bm{S}^{(1)}{\bm{S}^{(2)}}^{-1}+1}}.$

With $\Delta x=1/513$, $\Delta t=5$ and $\alpha_{m},\lambda_{m}$ for air,
water and steel, see Table 1, the spectral radii are $\rho_{\text{air-
steel}}^{\text{Jacobi}}=0.037$, $\rho_{\text{air-
water}}^{\text{Jacobi}}=0.059$ and $\rho_{\text{water-
steel}}^{\text{Jacobi}}=0.528$.

####  8.1.3 Results

All numerical experiments were run on an Intel i5-2500K 3.30 GHz CPU with
Python 3.6.9, Open MPI 2.1.1, FEniCS 2019.2.0.dev0 [16]. The code is available
at [18].

For GS we denote the different orderings of GS WR by ”GS_DN” and ”GS_ND”. We
use the result for GS in the ”DN” ($\Omega_{1}$ first) order with
$TOL_{WR}=10^{-11}$ as the reference result. Despite the asynchronous method
being non-deterministic, results for $5$ simulations showed no notable
deviations in the number of iterations, we show the mean result.

First, we consider the error of the interface temperature $\bm{u}_{\Gamma}$,
over $k$, for $TOL_{WR}=10^{-10}$, $\Delta x=1/513$ and $N=200$, resulting in
a comparable accuracy in space and time. We use $\bm{u}_{\Gamma}^{(k)}$ in the
discrete interface $\mathcal{L}^{2}$ norm for the termination check (5) and
error computation.

Results in Figure 4 show convergence for all considered WR methods,
numerically verifying the result of Theorem 6.1. The convergence rates of our
new asynchronous method are in between Jacobi and GS in all test cases.
Performance results in Figure 5 show slight performance improvements compared
to the constant splitting WR methods. This is since our new method requires
less than twice the number of iterations of GS WR and is parallel.

Figure 4: Left to right: Air-steel, air-water and water-steel. $k$ marks the
number of iterations until (5), with $TOL_{WR}=10^{-10}$, is met.

Figure 5: Left to right: Air-steel, air-water and water-steel. We use
MPI_Wtime to measure the wall-clock runtime.

###  8.2 Gas quenching

We now consider an example for the gas quenching application [33, 2]; cooling
a hot (steel) work-piece with air. We model the fluid (air) via the
compressible Euler equations

$\displaystyle\begin{split}\partial_{t}\rho+\nabla\cdot\rho\,\bm{v}&=0,\\\
\partial_{t}\rho v_{i}+\sum_{j=1}^{2}\partial_{x_{i}}(\rho
v_{i}v_{j}+\delta_{ij}p)&=0,\quad i=1,2,\\\ \partial_{t}\rho
E+\nabla\cdot(\rho H\bm{v})&=0.\end{split}$ (27)

Here $\rho$, $p$ and $v_{1},\,v_{2}$ are density, pressure and velocities. The
enthalpy is $H=E+p/\rho$ with total energy (per unit mass)
$E=e+|\bm{v}|^{2}/2$ and specific internal energy $e$. The system is completed
with the ideal gas law $p=(\gamma-1)\rho\,e$, where $\gamma$ is the adiabatic
exponent, here $\gamma=1.4$. The temperature is $T=p/(\rho\cdot R)$ with the
specific gas constant $R=287.058$ for dry air.

We model the solid by the nonlinear heat equation

$\displaystyle\begin{split}\alpha(u)\partial_{t}u(t,\bm{x})-\nabla\cdot\left(\lambda(u)\nabla
u(t,\bm{x})\right)=0,&\quad(t,\bm{x})\in(0,T_{f}]\times\Omega_{2},\\\ \nabla
u(t,\bm{x})\cdot\bm{n}=0,&\quad(t,\bm{x})\in[0,T_{f}]\times\partial\Omega_{2}\setminus\Gamma,\\\
\lambda(u)\nabla
u(t,\bm{x})\cdot\bm{n}_{\Gamma}=q(t,\bm{x}),&\quad(t,\bm{x})\in[0,T_{f}]\times\Gamma,\\\
u(0,\bm{x})=u_{0}^{\text{solid}},&\quad\bm{x}\in\Omega_{2}.\end{split}$ (28)

Here, we use the material parameters of 51CrV4 steel from [27], which are
given by

$\displaystyle\lambda(u)$
$\displaystyle=40.1+0.05u-0.0001u^{2}+4.9e-8u^{3},\quad\alpha(u)=7836c_{p}(u)$
$\displaystyle c_{p}(u)$
$\displaystyle=-10\text{ln}\left(\left(\text{e}^{c_{p1}(u)/10}+\text{e}^{c_{p2}(u)/10}\right)/2\right),$
$\displaystyle c_{p1}(u)$ $\displaystyle=34.2\text{e}^{0.0026u}+421.15,\quad
c_{p2}(u)=956.5\text{e}^{-0.012(u-900)}+0.45u.$

The problem geometry, inspired by [3], is shown in Figure 6.

$\Omega_{1}$$\Omega_{2}$$\Gamma$$\bm{n}_{\Gamma}$InletOutlet$x_{1}$$x_{2}$$0$$0.1$$0.3$$0.4$$0$$0.02$$0.07$
Figure 6: Geometry, not to scale, for the gas quenching test case.

In $\Omega_{1}$ we define
$\hat{u}^{0}_{fluid}:=(\rho_{0},\rho_{0}v_{0},0,E_{0})$, with
$\rho_{0}=1.225$, $T_{0}=273.15$ and $v_{0}=0.8\text{Ma}$. We use
$\hat{u}^{0}_{fluid}$ for the inlet boundary and the farfield boundary at the
top, on the outlet we extrapolate interior values to obtain a zero flux. On
the bottom we employ a slip boundary condition. Additionally, at the interface
$\Gamma$ we use $T_{\Gamma}$, on the left we set the wall temperature to
$T_{0}$ and on the right we use $\nabla T\cdot\bm{n}_{\Gamma}=0$.

On $\Omega_{1}$ we discretize (27) in space using a $1$st order finite volume
discretization, implemented in DUNE [1, 5]. The is grid shown in Figure 7. We
discretize (28) in space using linear finite elements, implemented using
FEniCS [16], on a triangulation obtained from splitting the squares of a
cartesian grid with $\Delta x=0.00125$, matching $\Omega_{1}$ at the
interface.

Figure 7: Grid used in gas quenching test case. Starting from a cartesian grid
with $\Delta x=0.005$ we perform two steps of local grid refinement on all
cells whose centres $(x_{1},x_{2})$ fulfill the following criteria: First
step: $x_{1}>0.05$, $x_{2}<0.03$ and second step: $x_{1}>0.08$, $x_{2}<0.02$.
The resulting grid consists of $5132$ cells.

For the time-discretizazions, we use the Crank-Nicolson method with $\Delta
t=10^{-2}$, c.f. Section 8.1, for the solid. In the fluid we use the SDIRK2
method with $\Delta t=2\cdot 10^{-4}$, solving the nonlinear systems using a
Jacobian-free Newton-Krylov method with ILU preconditioning.

We obtain the initial condition $u^{0}_{fluid}$ on $\Omega_{1}$ for the
coupled simulation as follows: We compute a stationary solution with
$T_{\Gamma}=900$, starting with initial conditions $\hat{u}^{0}_{fluid}$ and
simulating until $T_{f}=0.01$ with a step-size of $\Delta t=10^{-5}$. On
$\Omega_{2}$ we use the constant initial value $u_{0}^{\text{solid}}=900$.

The problems are coupled using a Dirichlet-Neumann approach. That is, on
$\Omega_{1}$ we compute the heat-flux $q$ as the interface boundary condition
on $\Omega_{2}$ and on $\Omega_{1}$, we use the interface temperature
$T_{\Gamma}=u\big{|}_{\Gamma}$, computed in $\Omega_{2}$.

We compute the discrete heat-flux by computing

$q=\kappa(T^{\,*})\nabla T^{\,*}\cdot\bm{n}_{\Gamma},$ (29)

at the interface, where $T^{\,*}$ is the temperature on $\Omega_{1}$, linearly
reconstructed. The heat conductivity $\kappa$ in the fluid is given via the
Sutherland law by

$\frac{\kappa(T\,)}{\kappa_{ref}}=\left(\frac{T}{T_{ref}}\right)^{3/2}\frac{T_{ref}+S_{k}}{T+S_{k}},$

with $\kappa_{ref}=0.0241$, $T_{ref}=273$ and $S_{k}=194$.

For the coupled run, we successively perform WR until $T_{f}=5$, on time-
windows of length $0.1$. We use the same relaxation as Section 8.1, with
relaxation parameters for the air-steel material combination, see Table 1,
$\Delta t=10^{-2}$ and $\Delta x=0.00125$, which is the cell-size at the
interface.

Here, the computational load of the subsolvers is not well balanced, since the
fluid problem has more unknowns and requires smaller step-sizes. As a
consequence, our new method is almost identical to GS WR. The load-balancing
problem can be resolved by allocating more processors in the space
discretiazion of the fluid problem. This is subject to future work and
required for a sensible performance comparison with Jacobi and GS WR.

Figure 8 shows the WR updates for the first $3$ time-windows, which shows the
iteration converges rapidly, as expected [23]. For larger tolerances, chosen
in accordance with the errors in space and time, one requires at most $3-4$
iterations per time-window.

Figure 8 also shows the temperature at the left top tip of the steel plate
over time, showing a steady cooling effect. Figures 9 visualizes the
temperature at $T_{f}=5$.

Figure 8: Left: Updates of the interface temperatures in the $2$-norm in the
first $3$ time-windows. Right: Temperature at the top left tip of the steel
plate, $(x_{1},x_{2})=(0.1,0.02)$, c.f., Figure 6, over time. Figure 9:
Isolines of the temperatures at $T_{f}=5$. Here, the isolines in $\Omega_{1}$
and $\Omega_{2}$ are distinct.

##  9 Summary and conclusions

In this paper we presented a novel parallel WR method utilizing asynchronous
communication during time-integration. Here, MPI One-sided communication is
instrumental in the implementation of this method. The analytical description
and convergence results of our new method extend existing linear WR theory by
including splittings variable in time and iteration. We present an algorithm
for performing optimal variable relaxation for two coupled problems.

Our new method is unconditionally parallel due to the use of asynchronous
communication. If used for coupled problems with poor load-balancing, our
method is equivalent to classical Gauss-Seidel WR.

Numerical results demonstrate convergence of our method, with a convergence
rate faster than Jacobi WR. Additionally, a performance comparison with two
subsolvers with an approximately equal computational workload, show
performance improvements for our new method, in comparison with Gauss-Seidel
and Jacobi WR.

WR methods enable the coupling of subsolvers in a partitioned manner. For our
numerical experiments we implemented our PDE subsolvers using the open source
packages DUNE [1] and FEniCS [16], and developed suitable adapters to
facilitate this coupling. This coupling is shown to work well in a conjugate
heat transfer test case. There, we couple the compressible Euler equations
with a nonlinear heat equation, using finite volume resp. linear finite
element discretizations in space.

#### Acknowledgments

We’d like to thank Joachim Hein for help with MPI, Robert Klöfkorn for help in
implementing the DUNE subsolver and Benjamin Rodenberg for collaboration in
developing the coupling to FEniCS.

## References

  * [1] P. Bastian, M. Blatt, A. Dedner, N. A. Dreier, C. Engwer, R. Fritze, C. Gräser, C. Grüninger, D. Kempf, R. Klöfkorn, M. Ohlberger, and O. Sander, “The DUNE framework: Basic concepts and recent developments,” _Computers and Mathematics with Applications_ , vol. 81, pp. 75–112, 2021.
  * [2] P. Birken, “Numerical Methods for the Unsteady Compressible Navier-Stokes Equations,” Habilitation Thesis, Kassel University, 2012.
  * [3] P. Birken, K. J. Quint, S. Hartmann, and A. Meister, “A time-adaptive fluid-structure interaction method for thermal coupling,” _Computing and Visualization in Science_ , vol. 13, no. 7, pp. 331–340, 2010.
  * [4] P. Crosetto, P. Reymond, S. Deparis, D. Kontaxakis, N. Stergiopulos, and A. Quarteroni, “Fluid-structure interaction simulation of aortic blood flow,” _Computers and Fluids_ , vol. 43, no. 1, pp. 46–57, 2011.
  * [5] A. Dedner and R. Klöfkorn, “Extendible and efficient python framework for solving evolution equations with stabilized discontinuous galerkin methods,” _arXiv_ , 2020.
  * [6] M. J. Gander and A. M. Stuart, “Space-time continuous analysis of Waveform Relaxation for the heat equation,” vol. 19, no. 6, pp. 2014–2031, 1998.
  * [7] E. Giladi and H. B. Keller, “Space-time domain decomposition for parabolic problems,” _Numerische Mathematik_ , vol. 93, no. 2, pp. 279–313, 2002\.
  * [8] E. Hairer, S. P. Nørsett, and G. Wanner, _Solving ordinary differential equations I_. Berlin: Springer, 1993.
  * [9] J. Janssen and S. Vandewalle, “Multigrid Waveform Relaxation on Spatial Finite Element Meshes : The Continuous-Time Case,” _SIAM Journal on Numerical Analysis_ , vol. 33, no. 2, pp. 456–474, 1996.
  * [10] ——, “Multigrid waveform relaxation on spatial finite element meshes: the discrete-time case,” _SIAM Journal on Scientific Computing_ , vol. 17, no. 1, pp. 133–155, 1996.
  * [11] ——, “On SOR waveform relaxation methods,” _SIAM Journal on Numerical Analysis_ , vol. 34, no. 6, pp. 2456–2481, 1997.
  * [12] D. Kowollik, V. Tini, S. Reese, and M. Haupt, “3D fluid–structure interaction analysis of a typical liquid rocket engine cycle based on a novel viscoplastic damage model,” _International journal for numerical methods in engineering_ , vol. 94, no. 13, pp. 1165–1190, 2013.
  * [13] D. S. C. Kowollik, P. Horst, and M. C. Haupt, “Fluid-structure interaction analysis applied to thermal barrier coated cooled rocket thrust chambers with subsequent local investigation of delamination phenomena,” vol. 4, pp. 617–636, 2013.
  * [14] J. W. M. Kranenborg, “Partitioned Multirate Time Integration for Coupled Systems of Ordinary Differential Equations,” Master thesis, Lund University, 2018.
  * [15] E. Lelarasmee, “The Waveform Relaxation Method for time domain analysis of large scale integrated circuits: Theory and Applications,” PHD Thesis, U.C. Berkeley, 1982.
  * [16] A. Logg, K.-A. Mardal, and G. Wells, _Automated solution of differential equations by the finite element method: The FEniCS book_. Springer, 2012, vol. 84.
  * [17] A. Lumsdaine and D. Wu, “Krylov subspace acceleration of waveform relaxation,” _SIAM Journal on Numerical Analysis_ , vol. 41, no. 1, pp. 90–111, 2003.
  * [18] P. Meisrimel, “Waveform relaxation with asynchronous time-integration,” 2021\. [Online]. Available: https://gitlab.maths.lth.se/PeterMeisrimel/asynch-WFR
  * [19] P. Meisrimel, A. Monge, and P. Birken, “A time adaptive multirate Dirichlet–Neumann waveform relaxation method for heterogeneous coupled heat equations,” _arXiv preprint_.
  * [20] MessagePassingInterfaceForum, “MPI: A Message-Passing Interface Standard Version 3.0,” Tech. Rep., 2012.
  * [21] U. Miekkala and O. Nevanlinna, “Convergence of Dynamic Iteration Methods for Initial Value Problems,” _SIAM Journal on Scientific and Statistical Computing_ , vol. 8, no. 4, pp. 459–482, 1987.
  * [22] A. Monge, “Partitioned methods for time-dependent thermal fluid-structure interaction,” Ph.D., Lund University, 2018.
  * [23] A. Monge and P. Birken, “On the convergence rate of the Dirichlet-Neumann iteration for unsteady thermal fluid-structure interaction,” _Computational Mechanics_ , pp. 1–17, 2017.
  * [24] O. Nevanlinna, “Remarks on Picard-Lindelöf Iteration, Part I,” _BIT Numerical Mathematics_ , vol. 29, no. April 1988, pp. 328–346, 1989\.
  * [25] ——, “Remarks on Picard-Lindelöf Iteration, Part II,” _BIT Numerical Mathematics_ , vol. 29, no. April 1988, pp. 535–562, 1989.
  * [26] A. Quarteroni and V. Alberto, _Domain decomposition methods for partial differential equations_. Oxford University Press, 1999.
  * [27] K. J. Quint, S. Hartmann, S. Rothe, N. Saba, and K. Steinhoff, “Experimental validation of high-order time integration for non-linear heat transfer problems,” _Computational Mechanics_ , vol. 48, no. 1, pp. 81–96, 2011\.
  * [28] M. W. Reichelt, J. K. White, and J. Allen, “Optimal convolution SOR acceleration of Waveform Relaxation with application to paralla simulation of semiconductor devices,” _SIAM Journal on Scientific Computing_ , vol. 16, no. 5, pp. 1137–1158, 1995.
  * [29] B. Rüth, B. Uekermann, M. Mehl, P. Birken, A. Monge, and H. J. Bungartz, “Quasi-Newton waveform iteration for partitioned surface-coupled multiphysics applications,” _International Journal for Numerical Methods in Engineering_ , pp. 1–22, 2020.
  * [30] S. Vandewalle, _Parallel multigrid waveform relaxation for parabolic problems_. Springer-Verlag, Berlin, 2013\.
  * [31] J. K. White and A. Sangiovanni-Vincentelli, _Relaxation Techniques for the Simulation of VLSI Circuits_. Boston: Kluwer Academic Publishers, 1987.
  * [32] R. M. Wilcox, “Exponential operators and parameter differentiation in quantum physics,” _Journal of Mathematical Physics_ , vol. 8, no. 4, pp. 962–982, 1967.
  * [33] P. W. Yarrington and E. A. Thornton, “Finite element analysis of low-speed compressible flows within convectively cooled structures,” _Journal of thermophysics and heat transfer_ , vol. 8, no. 4, pp. 678–686, 1994.