# Coupling Parareal with Non-Overlapping Domain Decomposition Method

Rim Guetat

Faculty of Sciences of Tunis, University of Tunis El Manar
2060 Tunis, Tunisia
rymguetat@yahoo.fr

**ABSTRACT.** In this paper, we present a new parallel algorithm for time dependent problems based on coupling parareal with non-overlapping domain decomposition method in order to increase parallelism in time and in space. For this we focus on the iterative methods of parallelization in space to solve the interface problem like Neumann-Neumann method. In the new algorithm, the coarse temporel propagator is defined on the global domain and the Neumann-Neumann method is chosen as a fine propagator with a few iterations. We present the rigorous convergence analysis of the new coupled algorithm on bounded time interval. Numerical experiments illustrate the performance of this new algorithm and confirm our analysis.

**RÉSUMÉ.** Dans ce papier, nous présentons un nouvel algorithme parallèle pour les problèmes dépendant du temps basé sur le couplage du pararéel avec les méthodes de décomposition de domaine sans recouvrement afin d'augmenter le parallélisme dans le temps et l'espace. Nous nous concentrons sur les méthodes itératives de parallélisation en espace pour résoudre le problème d'interface par la méthode de Neumann-Neumann. Dans ce nouvel algorithme, le propagateur grossier est défini sur le domaine global et la méthode de Neumann-Neumann est choisie pour le propagateur fin avec quelques itérations. Nous présentons l'analyse rigoureuse de convergence du nouvel algorithme couplé sur un intervalle de temps borné. Des expériences numériques illustrent les performances de ce nouvel algorithme et confirment notre analyse.

**KEYWORDS:** Parallel Algorithm, Parareal algorithm, Non-overlapping domain decomposition method, Schur complement, Neumann-Neumann method.

**MOTS-CLÉS:** Algorithme parallèle, Algorithme du pararéel, Méthodes de décomposition de domaine sans recouvrement, Complément de Schur, Méthode de Neumann-Neumann.
# 1. Introduction

In this paper we will focus on domain decomposition methods (DDM)[16, 14] applied to Finite Element Method (FEM). Though non-overlapping domain decomposition methods are powerful techniques for solving elliptic partial differential equations (PDEs). These methods consist of splitting the global domain into several non overlapping subdomains and using a parallel solver to solve the local problems. This procedure leads to iterative methods to find interface unknowns connecting domains together, namely the interface displacements in the Schur complement methods.

The performance of the method depends on the choice of efficient preconditioner for the interface operators. In this context three methods have been developed in 91's such as the Neumann-Neumann method, Finite Element Tearing and Interconnect (FETI) [6] and Balanced Domain Decomposition (BDD)[13]. In the present paper, the method of non-overlapping domain decomposition used is Neumann-Neumann method [9]. This method consists of solving the subdomain problems using Dirichlet interface conditions, followed by a correction step involving Neumann interface conditions.

Non-overlapping domain decomposition method can be extended to time dependent problems. The classical approach consists in discretizing in time with an implicit scheme and applying a domain decomposition method for the steady problems obtained at each time step [5]. The numerical temporal evolution schemes are sequential in nature and can be costly on a parallel computer. To reduce the computational cost over the whole time interval, recent interest arose in order to exploit modern engineering and modern hardware architectures, increased parallelism in numerical algorithms are of great importance. This development have lead to an increased focus on temporal parallelization in research, for this purpose we adopt the parareal algorithm introduced for the first time by Lions and al [11] which need to be solved in real-time. The parareal algorithm can be interpreted as different approaches, such as mutltigrid-in-time algorithm [7] or predictor-corrector method [3]. To be valuable to solve stiff models in real-time simulations, parareal use different models for the coarse and the fine solvers. The parareal algorithm has been applied to many fields by many authors [2, 12].

For problems with a larger size, or larger complexity, the parareal method suffers from the size of the spatial subproblem to solve, in this case we propose to combine the parareal algorithm with the domain decomposition method in order to reduce the complexity and performs the parallelism both in space and in time subdomains. There are at least three variants of the parareal algorithm which all contain some spatial parallelism: the first variant, which is the subject of this contribution, uses for the fine propagator a domain decomposition method; the second variant uses for the coarse propagator a domain decomposition method, and finally the third variant employs for both the coarse and fine propagators a domain decomposition method.

We develop in this paper a new parallel algorithm for advection-dominated diffusion problem by coupling the parareal algorithm and the non-overlapping domain decomposition like Neumann-Neumann method.

In this new algorithm, the Neumann-Neumann method is chosen as the fine propagator on each sub-interval with only a few iterations, where the interface problem is solved by Krylov type such as GMRES. In particular there is no coarse mesh needed in the case of many subdomains, if the time step is not too big. The remainder of the paper proceeds as follows. The model problem is presented in Section 2. In Section 3, we present the Neumann-Neumann method. Section 4 describes the parareal in time algorithm. We
propose our new coupled-algorithm in Section 5. Some experimental results are reported in Section 6. Our conclusion is given finally in Section 7.

## 2. Model problem

Let $\Omega$ be a bounded domain in $\mathbb{R}^d$, $d = 1, 2, 3$ with a Lipschitz boundary $\partial\Omega$. Our model problem is the advection-dominated diffusion problem defined by

$$ egin{cases} 
\frac{\partial u}{\partial t} - \nu \Delta u + \alpha \cdot \nabla u = f(t) \quad \text{in } \Omega \times ]0, T], \\ 
u(\mathbf{x}, t = 0) = u_0 \quad \text{in } \Omega, \\ 
u(\mathbf{x}, t) = 0 \quad \text{on } \partial \Omega \times ]0, T], 
\end{cases} 
\quad (1) $$

where $\alpha = (\alpha_1, \alpha_2)$ is the advection velocity and $f$ is the second term. We note that for the coefficient of viscosity $\nu$ is important and the order of the Peclet number $Pe$ is small, the advection term is significantly smaller than the diffusion term. Physically, diffusion dominates and advection is negligible.

If $u_0 \in H_0^1(\Omega)$, the unique weak solution $u$ is in $L^2(0, T; H^2(\Omega)) \cap H^1(0, T; L^2(\Omega))$. Once equation (1) discretized in time by an implicit scheme of first order, with a uniform time step equal to $\Delta t$, the model problem reduces at every time step $T_{n+1} = (n+1)\Delta t$ to the solution of

$$ \bar{\mathcal{L}}(u^{n+1}) := \left(\frac{1}{\Delta t} - \nu\Delta + \alpha\nabla\right)u^{n+1} = f^{n+1} + \frac{u^n}{\Delta t} = \underline{f}^{n+1}. \quad (2) $$

After discretization of the problem by finite element we have

$$ A\mathbf{u}^{n+1} := (I + \Delta t K)\mathbf{u}^{n+1} = b^{n+1} \quad (3) $$

where $\mathbf{u}^n(\mathbf{x})$ is the vector of nodal values of the approximation of the solution $u(\mathbf{x}, T_n)$, $I$ is a identity matrix, $K$ is a symmetric positive definite matrix arising from the discretization of the advection-dominated diffusion problem and $b^{n+1}$ the vector of nodal values of the right hand side depending on $\mathbf{u}^n(\mathbf{x})$. At a fixed time level $T_{n+1}$, we have to solve the problem (2) depending on the time step parameter by applying domain decomposition methods.

## 3. The Neumann-Neumann preconditioner

We decompose the computational domain $\Omega$ into finitely many non-overlapping subdomains $\Omega_i$, $i = 1...N_s$ and denote the subdomain boundaries by $\partial\Omega_i$. The local interfaces are defined by $\Gamma_{ij} = (\partial\Omega_i \cap \Omega_j) \setminus \partial\Omega$. Furthermore, at a fixed time level $T_{n+1}$, we define the global interface of this decomposition $\Gamma = \bigcup_{i \neq j} \Gamma_{ij}$. Let $A^{(i)}$ be the local stiffness matrix corresponding to subdomain $\Omega_i$, $\mathbf{u}_i^{n+1}$ be the vector of degrees of freedom corresponding to all elements in $\Omega_i$, and let $R^{(i)}$ denote the restriction matrix with entries 0 or 1 mapping the degrees of freedom $\mathbf{u}_i$ into global degrees of freedom. We introduce the interface degrees of freedom, $\gamma^{n+1} := R^{(i)}\mathbf{u}_i^{n+1}|_{\Gamma}$ and $b^{(i,n+1)}$ the vector of nodal values of the right hand side in $\Omega_i$. The stiffness matrix $A$ is obtained by the standard subassembly process

$$ A = \sum_{i=1}^{N_s} R^{(i)T} A^{(i)} R^{(i)}, \quad (4) $$
the system $Au^{n+1} = b^{n+1}$ can now be rewritten

$$ 
\begin{bmatrix}
A_{II}^{(1)} & \dots & 0 & A_{\Gamma\Gamma}^{(1)} R^{(1)} \\
\vdots & \ddots & \vdots & \vdots \\
0 & \dots & A_{II}^{(N_s)} & A_{\Gamma\Gamma}^{(N_s)} R^{(N_s)} \\
R^{(1)T} A_{\Gamma\Gamma}^{(1)T} & \dots & R^{(N_s)T} A_{\Gamma\Gamma}^{(N_s)T} & \sum_{i=1}^{N_s} R^{(i)T} A_{\Gamma\Gamma}^{(i)} R^{(i)}
\end{bmatrix}
\begin{bmatrix}
\mathbf{u}_{1,I}^{n+1} \\
\vdots \\
\mathbf{u}_{N_s,I}^{n+1} \\
\gamma^{\mathbf{n}+\mathbf{1}}
\end{bmatrix}
=
\begin{bmatrix}
b_I^{(1,n+1)} \\
\vdots \\
b_I^{(N_s,n+1)} \\
\sum_{i=1}^{N_s} R^{(i)T} b_{\Gamma}^{(i,n+1)}
\end{bmatrix}
\tag{5}
$$

Eliminating the unknowns corresponding to the subdomain interior i.e. $\mathbf{u}_{i,I}^{n+1}$ and in (3) we arrive at

$$ \left( \sum_{i=1}^{N_s} R^{(i)T} \left( A_{\Gamma\Gamma}^{(i)} - A_{\Gamma\Gamma}^{(i)T} A_{II}^{(i)T} A_{\Gamma\Gamma}^{(i)} \right) R^{(i)} \right) \gamma_{\mathbf{n}+1} = \sum_{i=1}^{N_s} R^{(i)T} \left( b_{\Gamma,n+1}^{(i)} - A_{\Gamma\Gamma}^{(i)T} A_{II}^{(i)T} b_I^{(i,n+1)} \right). \tag{6} $$

We introduce the following notation

$$ 
\begin{align}
S_i &= \left( A_{\Gamma\Gamma}^{(i)} - A_{\Gamma\Gamma}^{(i)T} A_{II}^{(i)T} A_{\Gamma\Gamma}^{(i)} \right), & S &= \sum_{i=1}^{M} R^{(i)T} S_i R^{(i)}, \tag{7} \\
g_i^{n+1} &= \left( b_{\Gamma}^{(i,,n+1)} - A_{\Gamma\Gamma}^{(i)T} A_{II}^{(i)T} b_{I}^{(i,n+1)} \right), & g_{n+1} &= \sum_{i=1}^{M} R^{(i)T} g_i^{n+1}, \nonumber
\end{align}
$$

where $S_i$ is the local primal Schur complement. We obtain the following Schur complement problem on interface

$$ S\gamma^{\mathbf{n}+1} = g^{n+1}. \tag{8} $$

The parallel preconditioner "Neumann-Neumann" $P_{NN}$ is then defined as follows

$$ P_{NN}^{-1} = \sum_{i=1}^{N_s} R^{(i)T} D^{(i)} S_i^{-1} D^{(i)} R^{(i)}, \tag{9} $$

with $D_i$ is a scaling diagonal subdomain matrix and form a decomposition of unity on the interface.

$$ \sum_{i=1}^{N_s} R^{(i)T} D^{(i)} R^{(i)} = I_{\Gamma}, \tag{10} $$

the preconditioner Neumann-Neumann remains efficient for an arbitrary number of subdomains as long as the advective term is not too strong. For solving the preconditioned linear system $P_{NN}^{-1} S \gamma_{n+1} = P_{NN}^{-1} g$, we use the Generalized Minimal Residual (GMRES) [15, 17], an additional advantage in using GMRES is that it gives faster convergence and it does not require a relaxation parameter.

# 3.1. Convergence analysis

## 3.1.1. Schur complement algorithm in the continuous case

To define the Schur complement in the continuous case for the model problem (1) at each time step $T_{n+1}$ on the space domain $\Omega$, where $\Omega = [0, 1] \times \mathbb{R}$, we assume that the spatial domain is partitioned into two non-overlapping subdomains $\Omega_1 = [0, c] \times \mathbb{R}$ and
$\Omega_2 = [c, 1] \times \mathbb{R}$ where $c \in ]0, 1[$. We denote by $\Gamma = \{(x, y) \in [0, 1] \times \mathbb{R} | x = c\}$, $\mathbf{n}_i$ the unit outward normal for $\Omega$ on the interface $\Gamma$ and $\mathbf{n}$ the unit outward normal vector on $\Gamma$ directed from $\Omega_1$ to $\Omega_2$. We split the original problem (2) into two locals problems as follows

$$ \left\{ \begin{array}{ll} \bar{\mathcal{L}}u_i^{n+1}(x,y) & = \displaystyle\frac{f_i^{n+1}(x,y)}{\partial_i} \quad \text{in } \Omega_i, \quad i=1,2, \\ u_1^{n+1}(x,y) & = u_2^{n+1}(x,y) \quad \text{on } \Gamma, \\ \displaystyle\frac{\partial u_1^{n+1}}{\partial \mathbf{n}}(x,y) & = \displaystyle\frac{\partial u_2^{n+1}}{\partial \mathbf{n}}(x,y) \quad \text{on } \Gamma, \\ u_i^{n+1}(x,y) & = 0 \quad \text{on } \partial\Omega_i \cap \partial\Omega, \quad i=1,2. \end{array} \right. \qquad (11) $$

The split problem (11) with transmission conditions (11)$_2$, (11)$_3$ is equivalent to the original one, in the sense we that introduce the Dirichlet data on the interface $\gamma^{n+1} = u^{n+1}|_{\Gamma}$. Introducing the extension operators $\mathcal{H}_i$ and $\mathcal{G}_i$ as follows

$$ \left\{ \begin{array}{ll} \bar{\mathcal{L}}\mathcal{H}_i\gamma^{n+1} & = 0 & \text{in } \Omega_i, \quad i=1,2, \\ \mathcal{H}_i\gamma^{n+1} & = \gamma^{n+1} & \text{on } \Gamma, \\ \mathcal{H}_i\gamma^{n+1} & = 0 & \text{on } \partial\Omega_i \cap \partial\Omega, \end{array} \right. \qquad (12) $$

and

$$ \left\{ \begin{array}{ll} \bar{\mathcal{L}}\mathcal{G}_i f_i^{n+1} & = f_i^{n+1} \quad \text{in } \Omega_i, \quad i=1,2, \\ \mathcal{G}_i f_i^{n+1} & = 0 \quad \text{on } \Gamma, \\ \mathcal{G}_i f_i^{n+1} & = 0 \quad \text{on } \partial\Omega_i \cap \partial\Omega. \end{array} \right. \qquad (13) $$

The solution is constructed by defining

$$ u_i^{n+1} = \mathcal{H}_i \gamma^{n+1} + \mathcal{G}_i f_i^{n+1}, \quad i=1,2. \qquad (14) $$

The Steklov-Poincaré equation is given by

$$ S\gamma_{n+1} = \chi_{n+1} \qquad (15) $$

where $\chi_{n+1}$

$$ \chi_{n+1} := \frac{\partial \mathcal{G}_2 f_i^{n+1}}{\partial \mathbf{n}} - \frac{\partial \mathcal{G}_1 f_i^{n+1}}{\partial \mathbf{n}} \\ = -\sum_{i=1}^{2} \frac{\partial \mathcal{G}_i f_i^{n+1}}{\partial \mathbf{n}_i}, $$

and

$$ S\gamma := \frac{\partial \mathcal{H}_1 \gamma}{\partial \mathbf{n}} - \frac{\partial \mathcal{H}_2 \gamma}{\partial \mathbf{n}} \\ = \sum_{i=1}^{2} \frac{\partial \mathcal{H}_i \gamma}{\partial \mathbf{n}_i}. $$

The local Steklov-Poincaré operators $S_i$ are then defined as follows

$$ S_i \gamma := \frac{\partial \mathcal{H}_i \gamma}{\partial \mathbf{n}_i}. \qquad (16) $$

For accelerating the convergence we initialize the interface condition with the previous solution at convergence (computed in time $T_n$). $\gamma^{n+1,0} = \gamma^{n,P_{max}}$. Now the Schur
complement algorithm consists of the following steps: at every time step $T_{n+1}$, given an initial guess $\gamma^{0,0} = u_{0|\Gamma}$ along the interface $\Gamma$, compute for $p = 1, 2, \dots$

$$ \left\{ \begin{array}{ll} \bar{\mathcal{L}}(u_i^{n+1,p+1})(x,y) & = f_{i,n+1}(x,y) \quad \text{in } \Omega_i, i = 1, 2, \\ u_i^{n+1,p+1}(x,y) & = 0 \quad \text{on } \partial\Omega_i \cap \partial\Omega, \\ u_i^{n+1,p+1}(x,y) & = \gamma^{n+1,p} \quad \text{on } \Gamma, \end{array} \right. \qquad (17) $$

and then update the value along the interface using

$$ \gamma^{n+1,p+1} = \gamma^{n+1,p} - \theta(S\gamma^{n+1,p} - g^{n+1}), \quad \text{on } \Gamma, \qquad (18) $$

$\theta \in ]0, 1]$ being a relaxation parameter. To investigate the convergence of the iterative solutions of problem (17) we define the error $e_i^{n,p}$ between the exact solution $u^n$ and the iterates $u_i^{n,p}$ at a fixed time level $T_n$. By linearity, the error satisfies a homogeneous advection-dominated diffusion. Let $\tilde{\gamma}^{n,p} = \gamma^{n,p} - u_{|\Gamma}^n$ with $\tilde{\gamma}^p(0) = 0$ to satisfy the compatibility condition. The error equations are given by

$$ \left\{ \begin{array}{ll} \bar{\mathcal{L}}(e_1^{n+1,p+1}) & = 0 \quad \text{in } \Omega_1, \\ e_1^{n+1,p+1}(\cdot) & = 0 \quad \text{on } \partial\Omega_i \cap \partial\Omega, \\ e_1^{n+1,p+1}(\cdot) & = \tilde{\gamma}^{n+1,p}, \quad \text{on } \Gamma, \end{array} \right. \qquad (19) $$

$$ \left\{ \begin{array}{ll} \bar{\mathcal{L}}(e_2^{n+1,p+1}) & = 0 \quad \text{in } \Omega_2, \\ e_2^{n+1,p+1}(\cdot) & = 0 \quad \text{on } \partial\Omega_i \cap \partial\Omega, \\ e_2^{n+1,p+1}(\cdot) & = \tilde{\gamma}^{n+1,p} \quad \text{on } \Gamma, \end{array} \right. \qquad (20) $$

$$ \tilde{\gamma}^{n+1,p+1} = \tilde{\gamma}^{n+1,p} - \theta(S_1 + S_2)\tilde{\gamma}^{n+1,p} \quad \text{on } \Gamma, \\ = (Id - \theta(S_1 + S_2))\tilde{\gamma}^{n+1,p}, \quad \text{on } \Gamma. $$

We define the Fourier transform of $e_i^{n,p}$ in the $y$ direction by

$$ \hat{e}_i^{n,p}(x, w) := \int_{\mathbb{R}} \exp(-iwy) e_i^{n,p}(x, y) dy. \qquad (21) $$

Applying formally the Fourier transform to the equations (19) and (35), we obtain

$$ \left\{ \begin{array}{ll} \hat{\mathcal{L}}\hat{e}_1^{n+1,p+1}(x, w) & = 0 \quad \text{for } (x, w) \in ]0, c[ \times \mathbb{R}, \\ \hat{e}_1^{n+1,p+1}(0, w) & = 0 \quad \text{on } \mathbb{R}, \\ \hat{e}_1^{n+1,p+1}(c, w) & = \hat{\gamma}^{n+1,p}(w) \quad \text{on } \mathbb{R}, \end{array} \right. \qquad (22) $$

$$ \left\{ \begin{array}{ll} \hat{\mathcal{L}}\hat{e}_2^{n+1,p+1}(x, w) & = 0 \quad \text{pour } (x, w) \in ]c, 1[ \times \mathbb{R}, \\ \hat{e}_2^{n+1,p+1}(1, w) & = 0 \quad \text{on } \mathbb{R}, \\ \hat{e}_2^{n+1,p+1}(c, w) & = \hat{\gamma}^{n+1,p}(w) \quad \text{on } \mathbb{R}, \end{array} \right. \qquad (23) $$

where

$$ \hat{\mathcal{L}} = -\nu \frac{\partial^2}{\partial^2 x} + \alpha_1 \frac{\partial}{\partial x} + (i\alpha_2 \nu w + \nu w^2). \qquad (24) $$

The characteristic equation of the second order ordinary differential equation in the $x$-variable gives

$$ -\nu x^2 + \alpha_1 x + \left( \frac{1}{\Delta t} + i\alpha_2 w + \nu w^2 \right) = 0, \qquad (25) $$
where the solutions have the form

$$
\begin{align}
\hat{e}_1^{n+1, p+1}(x, w) &= C_{1,1}^{+} \exp(\sigma^{+}x) + C_{1,2}^{-} \exp(\sigma^{-}x), \quad x \in [0, c], \tag{26} \\
\hat{e}_2^{n+1, p+1}(x, w) &= C_{2,1}^{+} \exp(\sigma^{+}(x-1)) + C_{2,2}^{-} \exp(\sigma^{-}(x-1)), \quad x \in [c, 1], \nonumber
\end{align}
$$

with the characteristic roots

$$
\sigma^{+} = \frac{\alpha_{1} + \sqrt{r}}{2\nu} \quad \text{and} \quad \sigma^{-} = \frac{\alpha_{1} - \sqrt{r}}{2\nu},
$$

where $r = Re + iIm$ is a complex number with a real part $Re = \alpha_1^2 + 4\nu(\frac{1}{\Delta t} + \nu w^2)$
and a imaginary part $Im = 4\nu\alpha_2 w$. The complex square root is defined by

$$
\sqrt{r} = \frac{1}{\sqrt{2}} (Re + \sqrt{Re^2 + Im^2})^{1/2} + i \frac{\operatorname{Im}}{|\operatorname{Im}| \sqrt{2}} (-Re + \sqrt{Re^2 + Im^2})^{1/2} \quad (27)
$$

The coefficients $C_{i,1}$ and $C_{i,2}$ are computed from the boundary conditions in the error equations $(26)_1$, $(26)_2$, we get

$$
\left\{
\begin{array}{ll}
\hat{e}_1^{n+1,p+1}(0,w) = 0 & \Rightarrow C_{1,1} = -C_{1,2} \\
\hat{e}_1^{n+1,p+1}(c,w) = \hat{\gamma}^{n+1,p} & \Rightarrow C_{1,1} = \frac{\hat{\gamma}^{n+1,p}}{\exp(\frac{\alpha_1}{2\nu}c) \sinh(\frac{\sqrt{r}}{2\nu}c)} \\
\hat{e}_2^{n+1,p+1}(1,w) = 0 & \Rightarrow C_{2,1} = -C_{2,2} \\
\hat{e}_2^{n+1,p+1}(c,w) = \hat{\gamma}^{n+1,p} & \Rightarrow C_2^+ = \frac{\hat{\gamma}^{n+1,p}}{\exp(\frac{\alpha_1}{2\nu}(c-1)) \sinh(\frac{\sqrt{r}}{2\nu}(c-1))},
\end{array}
\right.
\quad (28)
$$

and inserting them into the error equations (26)$_1$, (26)$_2$ we deduce the subdomain solu-
tions

$$
\left\{
\begin{array}{ll}
\hat{e}_1^{n+1,p+1}(x, \Delta t, w) & = \exp\left(\frac{\alpha_1}{2\nu}(x-c)\right) \frac{\sinh\left(\frac{\sqrt{r}}{2\nu}x\right)}{\sinh\left(\frac{\sqrt{r}}{2\nu}c\right)} \hat{\gamma}^{n+1,p}, & x \in [0, c], \\
\\
\hat{e}_2^{n+1,p+1}(x, \Delta t, w) & = \exp\left(\frac{\alpha_1}{2\nu}(x-c)\right) \frac{\sinh\left(\frac{\sqrt{r}}{2\nu}(x-1)\right)}{\sinh\left(\frac{\sqrt{r}}{2\nu}(c-1)\right)} \hat{\gamma}^{n+1,p}, & x \in [c, 1],
\end{array}
\right.
\qquad (29)
$$

We obtain for the transformed local Steklov-Poincaré operators

$$
\left\{
\begin{aligned}
\hat{S}_1 \hat{\tilde{\gamma}}^{n+1,p} &= \frac{\partial \hat{e}_1^{n+1,p+1}}{\partial x} \Big|_{x=c} = \left[ \frac{\alpha_1}{2\nu} + \frac{\sqrt{r}}{2\nu} \coth\left(\frac{\sqrt{r}}{2\nu}c\right) \right] \hat{\tilde{\gamma}}^{n+1,p}, \\
\hat{S}_2 \hat{\tilde{\gamma}}^{n+1,p} &= -\frac{\partial \hat{e}_2^{n+1,p+1}}{\partial x} \Big|_{x=c} = -\left[ \frac{\alpha_1}{2\nu} + \frac{\sqrt{r}}{2\nu} \coth\left(\frac{\sqrt{r}}{2\nu}(c-1)\right) \right] \hat{\tilde{\gamma}}^{n+1,p},
\end{aligned}
\right.
\qquad (30)
$$

hence the global Steklov-Poincaré operator can be written as

$$
\hat{S} = \hat{S}_1 + \hat{S}_2 = \frac{\sqrt{r}}{2\nu} \left[ \coth\left(\frac{\sqrt{r}}{2\nu}c\right) - \coth\left(\frac{\sqrt{r}}{2\nu}(c-1)\right) \right]. \quad (31)
$$

**Remark** For studying the error solutions of the subproblems (17) over bounded subdom-
mains Ω₁ = [0, c] × [a, b] and Ω₂ = [c, 1] × [a, b], where 0 < a < b. We define for any
bounded functions h the norm

$$
\|h(\cdot, \cdot)\|_{\infty} := \sup_{(x,y) \in \Omega} |h(x, y)|.
$$
The error-estimates of the error solutions (19) and (35) are presented in the following lemma :

**Lemma 1.** (Super-Solutions) At each time step $T_{n+1}$, for all $p \in \mathbb{N}$ and $r = \alpha_1^2 + 4\nu(\frac{1}{\Delta t} + \nu b^2)$, the error solutions of the subproblems (17) with $\alpha_2 = 0$ satisfy

$$0 \le e_i^{n+1,p}(x, y) \le \tilde{e}_i^{n+1,p}, i = 1, 2,$$

where the super-solutions $\tilde{e}_i^{n+1,p}$ are the solutions of the linear problems

$$\left\{
\begin{array}{ll}
\bar{\mathcal{L}}(\tilde{e}_1^{n+1,p+1}) & =0 \quad \text{in } \Omega_1, \\
\tilde{e}_1^{n+1,p+1} & =0 \quad \text{on } \partial\Omega_i \cap \partial\Omega, \\
\tilde{e}_1^{n+1,p+1}(x) & =\exp\left(\frac{\alpha_1}{2\nu}(x-c)\right) \frac{\sinh\left(\frac{\sqrt{r}}{2\nu}x\right)}{\sinh\left(\frac{\sqrt{r}}{2\nu}c\right)} \|\tilde{\gamma}^{n+1,p}\|_{\infty}, x \in [0, c]
\end{array}
\right.
\quad (32)$$

$$\left\{
\begin{array}{ll}
\bar{\mathcal{L}}(\tilde{e}_2^{n+1,p+1}) & = 0 \quad \text{in } \Omega_2, \\
\tilde{e}_2^{n+1,p+1} & = 0 \quad \text{on } \partial\Omega_i \cap \partial\Omega, \\
\tilde{e}_2^{n+1,p+1}(x) & = \exp\left(\frac{\alpha_1}{2\nu}(x-c)\right) \frac{\sinh\left(\frac{\sqrt{r}}{2\nu}(x-1)\right)}{\sinh\left(\frac{\sqrt{r}}{2\nu}(c-1)\right)} \|\tilde{\gamma}^{n+1,p}\|_{\infty}, x \in [c, 1].
\end{array}
\right.
\quad (33)$$

The solutions to (32) and (33) state non-negative solutions

$$\left\{
\begin{array}{ll}
\tilde{e}_1^{n+1,p+1}(x, \Delta t) & = \exp\left(\frac{\alpha_1}{2\nu}(x-c)\right) \frac{\sinh\left(\frac{\sqrt{r}}{2\nu}x\right)}{\sinh\left(\frac{\sqrt{r}}{2\nu}c\right)} \|\tilde{\gamma}^{n+1,p}\|_{\infty}, \quad x \in [0, c], \\
\\
\tilde{e}_2^{n+1,p+1}(x, \Delta t, w) & = \exp\left(\frac{\alpha_1}{2\nu}(x-c)\right) \frac{\sinh\left(\frac{\sqrt{r}}{2\nu}(x-1)\right)}{\sinh\left(\frac{\sqrt{r}}{2\nu}(c-1)\right)} \|\tilde{\gamma}^{n+1,p}\|_{\infty}, \quad x \in [c, 1],
\end{array}
\right.
\quad (34)$$

*Proof.* The difference $E_i := \tilde{e}_i^{p+1} - e_i^{p+1}$, $i = 1, 2$ satisfy

$$\left\{
\begin{array}{ll}
\bar{\mathcal{L}}(E_i^{n+1,p+1}) & = 0 \quad \text{in } \Omega_i, \\
E_i^{n+1,p+1}(\cdot) & = 0 \quad \text{on } \partial\Omega_i \cap \partial\Omega, \\
E_i^{n+1,p+1}(\cdot) & \ge 0 \quad \text{on } \Gamma,
\end{array}
\right.
\quad (35)$$

Now we apply the maximum principle lemma in [4] therefore $E_i^{n+1,p+1} \ge 0$ and similarly for $\tilde{E}_i^{n+1,p+1} = \tilde{e}_i^{n+1,p+1} + e_i^{n+1,p+1} \ge 0$ the error solutions (19) and (35) can be bounded by

$$\left\{
\begin{array}{l}
| e_1^{n+1,p+1} | \le \tilde{e}_1^{n+1,p+1}(x, \Delta t) = \exp\left(\frac{\alpha_1}{2\nu}(x-c)\right) \frac{\sinh\left(\frac{\sqrt{r}}{2\nu}x\right)}{\sinh\left(\frac{\sqrt{r}}{2\nu}c\right)} \|\tilde{\gamma}^{n+1,p}\|_{\infty}, \quad x \in [0, c], \\
\\
| e_2^{n+1,p+1} | \le \tilde{e}_2^{n+1,p+1}(x, \Delta t) = \exp\left(\frac{\alpha_1}{2\nu}(x-c)\right) \frac{\sinh\left(\frac{\sqrt{r}}{2\nu}(x-1)\right)}{\sinh\left(\frac{\sqrt{r}}{2\nu}(c-1)\right)} \|\tilde{\gamma}^{n+1,p}\|_{\infty}, \quad x \in [c, 1],
\end{array}
\right.
\quad (36)$$

Evaluating the derivatives of these equations (36) at $x = c$ taking the supremum over $\Gamma$ the symbol of the global Steklov-Poincaré operator can be bounded by the quantity

$$\frac{\sqrt{r}}{2\nu} \left[ \coth\left(\frac{\sqrt{r}}{2\nu}c\right) - \coth\left(\frac{\sqrt{r}}{2\nu}(c-1)\right) \right]$$

□
### 3.1.2. Neumann-Neumann algorithm

The Neumann-Neumann algorithm starts with an initial guess $\psi_{i,n+1}^0$ along the interface $\Gamma \times [0, T]$ and then computes for $\theta \in ]0, 1]$ simultaneously for $i = 1, 2$ with $p = 1, 2, \dots$

$$ \left\{ \begin{array}{ll} \bar{\mathcal{L}} u_i^{n+1, p+1} & = \displaystyle \frac{f}{\Gamma}^{n+1} \quad \text{in } \Omega_i, \quad i = 1, 2, \\ u_i^{n+1, p+1} & = \gamma^{n+1, p} \quad \text{on } \Gamma, \\ u_i^{n+1, p+1} & = 0 \quad \text{on } \partial\Omega_i \cap \partial\Omega, \quad i = 1, 2, \end{array} \right. \tag{37} $$

$$ \left\{ \begin{array}{ll} \bar{\mathcal{L}}\psi_i^{n+1,p+1} & = 0 & \text{in } \Omega_i, \quad i = 1, 2, \\ \frac{\partial\psi_i^{n+1,p+1}}{\partial n_i} & = \frac{1}{2} \left( \frac{\partial u_1^{n+1,p+1}}{\partial n_1} + \frac{\partial u_2^{n+1,p+1}}{\partial n_2} \right) & \text{on } \Gamma, \\ \psi_i^{n+1,p+1} & = 0 & \text{on } \partial\Omega_i \cap \partial\Omega, \quad i = 1, 2, \end{array} \right. \tag{38} $$

$$ \gamma^{n+1,p+1} = \gamma^{n+1,p} - \theta \left( \frac{1}{2} \psi_1^{n+1,p+1} + \frac{1}{2} \psi_2^{n+1,p+1} \right) \text{ on } \Gamma. \tag{39} $$

The Neumann-Neumann preconditioner is a weighted sum of inverse of transform local Steklov-Poinaré. Thus the Fourier transform preconditioner for two sub-domains is defined by:

$$ \widehat{P_{NN}}^{-1} = \frac{1}{4}(\hat{S}_1^{-1} + \hat{S}_2^{-1}), $$

where from (30)

$$ \begin{aligned} \hat{S}_1^{-1} &= \frac{2\nu}{\alpha_1 + \sqrt{r} \coth\left(\frac{\sqrt{r}}{2\nu}c\right)}, \\ \hat{S}_2^{-1} &= \frac{2\nu}{\sqrt{r} \coth\left(\frac{\sqrt{r}}{2\nu}(c-1)\right) - \alpha_1}. \end{aligned} $$

However, the Fourier transform of the interface equation (39) becomes

$$ \hat{\gamma}^{n+1,p+1} = (Id - \theta \widehat{P_{NN}}^{-1} \hat{S}) \hat{\gamma}^{n+1,p} \text{ on } \Gamma, \tag{40} $$

where the symbol of the preconditioned operator is easily determined as

$$ \widehat{P}_{NN}^{-1} \hat{S} = \frac{r \left( \coth\left(\frac{\sqrt{r}}{2\nu}c\right) + \coth\left(\frac{\sqrt{r}}{2\nu}(c-1)\right) \right)^2}{\left( \sqrt{r} \coth\left(\frac{\sqrt{r}}{2\nu}c\right) + \alpha_1 \right) \left( \sqrt{r} \coth\left(\frac{\sqrt{r}}{2\nu}(c-1)\right) - \alpha_1 \right)}. \tag{41} $$

### 3.1.3. Convergence analysis

The main goal of the analysis is to estimate the efficiency of the proposed preconditioner. For this purpose we study the convergence factor. In numerical calculation, the frequency parameter $w$ can not be high and we have $|w| \in [w_{min}, w_{max}]$, where $w_{min} = \pi$ and $w_{max} = \pi/h$ are the frequency minimum and maximum values. The updating boundary condition in (40) satisfies the following inequality

$$ \max_{|w| \in [w_{min}, w_{max}]} |\hat{\gamma}^{n+1,p}| = \rho^p \max_{|w| \in [w_{min}, w_{max}]} |\hat{\gamma}^{n+1,0}|, \tag{42} $$
where the optimal convergence factor $\rho^p$ is defined by

$$ \rho := \max_{|w| \in [w_{\min}, w_{\max}]} \left| 1 - \theta \frac{r \left( \coth\left(\frac{\sqrt{r}}{2\nu}c\right) + \coth\left(\frac{\sqrt{r}}{2\nu}(c-1)\right) \right)^2}{\left( \sqrt{r} \coth\left(\frac{\sqrt{r}}{2\nu}c\right) + \alpha_1 \right) \left( \sqrt{r} \coth\left(\frac{\sqrt{r}}{2\nu}(c-1)\right) - \alpha_1 \right)} \right| \quad (43) $$

**Remark** The convergence factor depends on the relaxation parameter $\theta$ but in practical computation Krylov subspace method is more robust and gives better results. For this purpose, in the discrete case, we replace the Richardson iteration by GMRES iteration.

**Lemma 2.** For a symmetric decomposition $c = 0.5$. The convergence factor $\rho_{GMRES}$ (see [1]) of Neumann-Neumann preconditioning in a GMRES iteration is bounded by

$$ \rho_{GMRES} \le \left( 1 - \frac{\left( \min_{w \in [\pi, \frac{\pi}{h}]} (\hat{P}_{NN}^{-1} \hat{S}) \right)^2}{\max_{w \in [\pi, \frac{\pi}{h}]} |\hat{P}_{NN}^{-1} \hat{S}|^2} \right)^{1/2}, \quad (44) $$

where $\hat{P}_{NN}^{-1} \hat{S}$ is given by (41) for $c = 0.5$.

*Proof.* The reduction factor in a GMRES iteration, for a positive real matrix $A$ with symmetric part $M = \frac{1}{2}(A + A^T)$, is bounded from above [15, 1]

$$ \rho_{GMRES} = \left( 1 - \frac{(\lambda_{min}(M))^2}{\lambda_{max}(A^T A)} \right)^{1/2}. \quad (45) $$

where $A = \hat{P}_{NN}^{-1} \hat{S}$, $\lambda_{min}$ and $\lambda_{max}$ are the minimum and maximum eigenvalues. After a short calculation, we deduce from (41) for $c = 0.5$ the quantity for the preconditioned operator $\hat{P}_{NN}^{-1} \hat{S}$

$$ \hat{P}_{NN}^{-1} \hat{S}(w) = \frac{r \coth^2\left(\frac{\sqrt{r}c}{2\nu}\right)}{r \coth^2\left(\frac{\sqrt{r}c}{2\nu}\right) - \alpha_1^2}. \quad (46) $$

□

# 4. The parareal in time algorithm

The exact flow of the model problem (1) is denoted $u(t) = \mathcal{E}_t(u_0)$. The strategy of parareal algorithm is to do a time decomposition in the spirit of domain decomposition. We divide the total simulation time into $N$ subintervals $[T_n, T_{n+1}]$, $n = 0, \dots, N-1$ of size $\Delta T$. Then $u(T_n) = \mathcal{E}_{\Delta T}(u(T_n)) = \mathcal{E}_{n\Delta T}(u_0)$, from the semigroup property of $\mathcal{E}$. In general $\mathcal{E}$ is not realizable and can only be approximated, we thus introduce a fine and precise propagator $\mathcal{F}$, it is based on an appropriate classical discretization scheme based on a small enough time step $\delta t$, then we have the approximation $u(T_n) \simeq U_n = \mathcal{F}_{n\Delta T}(u_0)$. The original parareal algorithm makes possible to define iteratively a sequence $U_n^k$ that converges toward $U_n$ as $k$ goes to infinity. The definition of this sequence requires a coarse propagator $\mathcal{G}$, which represents a cheap and less accurate
approximation to $\mathcal{F}_{\Delta T}$ and based on a similar scheme in time with the larger time step $\Delta T$, $\delta t \ll \Delta T$. The propagator $\mathcal{G}$ has to be run sequentially and $\mathcal{F}$ in parallel. In this way, the algorithm starts with an initial approximation $U_n^0$ given for example by the sequential computation:

$$U_0^0 = u_0, \quad U_n^0 = \mathcal{G}_{\Delta T}(U_{n-1}^0),$$

the numerical scheme of the parareal method is given by:

$$U_n^{k+1} = \mathcal{G}_{\Delta T}(U_{n-1}^{k+1}) + \mathcal{F}_{\Delta T}(U_{n-1}^k) - \mathcal{G}_{\Delta T}(U_{n-1}^k), \quad (47)$$

with the initial condition $U_0^k = u_0$. Note that after $k$ iterations of the parareal method, the solution $U_n^k$ for $m \le k$ is exactly equal to the numerical solution given by using the $\mathcal{F}$ propagator. In another hand, it is clear that we have $U_n = \mathcal{F}_{n\Delta T}(0, u_0)$ for any $n$. The parareal algorithm remains efficient for a large number $N$ of time-subdomains, the parareal must converge with $k \ll N$.

## 4.1. Convergence analysis

We review here the stability and the convergence of the standard parareal algorithm in two dimensions. The convergence have been first presented in [3] on the entire domain $\mathbb{R}$. We define the error derived from the standard algorithm (47) for $n = 1, \dots, N$ and $k = 0, \dots, k_{conv}$ by $e_n^k := U_n^k - U_n$. After adding and subtracting $\mathcal{G}_{\Delta T}(U_{n-1})$ and $\tilde{\mathcal{F}}_{\Delta T}(U_{n-1}^{k-1})$ this error $e_n^k$ satisfies for $n \ge 1$ and $k \ge 1$

$$\begin{aligned} e_n^k &= \mathcal{G}_{\Delta T}(U_{n-1}^k) - \mathcal{G}_{\Delta T}(U_{n-1}^{k-1}) \\ &\quad + \mathcal{F}_{\Delta T}(U_{n-1}^{k-1}) - \mathcal{F}_{\Delta T}(U_{n-1}^k) \\ &= [\mathcal{G}_{\Delta T}(U_{n-1}^k) - \mathcal{G}_{\Delta T}(U_{n-1})] \\ &\quad + [\mathcal{G}_{\Delta T}(\mathcal{F}_{T_{n-1}}(u_0)) - \mathcal{F}_{\Delta T}(U_{n-1})] \end{aligned} \quad (48)$$

We now denote by $\Delta\mathcal{G}$ the difference $\mathcal{F} - \mathcal{G}$. Then the relation (48) becomes over $\Omega$; for every $n$ and $k$

$$e_n^k = [ \mathcal{G}_{\Delta T}(U_{n-1}^k) - \mathcal{G}_{\Delta T}(U_{n-1}) ] + [ \Delta\mathcal{G}_{\Delta T}(U_{n-1}) - \Delta\mathcal{G}_{\Delta T}(U_{n-1}^{k-1}) ] \quad (49)$$

Our purpose now is to have an estimation of the error $e_n^k(x)$ based on the discrete Fourier-transform of the exact and discrete solutions which is given for any $h$ in the bounded domain in $\mathbb{R}^2$ by the Fourier series

$$\hat{h}(t, \xi_1, \xi_2) = \sum_{\xi_1 = -\infty}^{+\infty} \sum_{\xi_2 = -\infty}^{+\infty} h(t, x, y) \exp(-i(\xi_1 x + \xi_2 y)).$$

We apply the Fourier transform in $x$ and $y$ directions to the model problem (1) with $u_0$ a periodic function, we obtain

$$\left\{ \begin{array}{l} \displaystyle \frac{\partial \hat{u}}{\partial t}(t, \xi_1, \xi_2) + (\nu(\xi_1^2 + \xi_2^2) + i(\alpha_1 \xi_1 + \alpha_2 \xi_2))\hat{u}(t, \xi_1, \xi_2) = 0, \quad (\xi_1, \xi_2) \in \mathbb{Z} \times \mathbb{Z} \\[1em] \displaystyle \hat{u}(t, \xi_1, \xi_2) = \hat{u}_0(\xi_1, \xi_2) \end{array} \right. \quad (50)$$

We note $\delta(\xi_1, \xi_2) := ((\nu(\xi_1^2 + \xi_2^2) + i(\alpha_1 \xi_1 + \alpha_2 \xi_2))\Delta T$ and $\xi = (\xi_1, \xi_2)$. The exact solution of problem (50) is given at time $T_n$ by

$$\hat{u}(\xi, T_n) = \exp(-n\delta(\xi))\hat{u}_0(\xi) = \hat{\mathcal{E}}_{T_n}(\delta(\xi), T_0, \hat{u}_0) \quad (51)$$
In a similar manner to (47), we introduce two approximations $\hat{\mathcal{F}}$ and $\hat{\mathcal{G}}$ of $\hat{\mathcal{E}}_{T_n}$ and we define a sequence $\{\hat{U}_n\}_{n \ge 1}$ by

$$ \hat{U}_n(\xi) = \hat{\mathcal{F}}_{\Delta T}(\delta(\xi), \hat{U}_{n-1}(\xi)). \quad \forall n, \hat{U}_0 = \hat{u}_0, $$

we construct a sequence $\{\hat{U}_n^k\}_k$ by

$$ \hat{U}_0^0(\xi) = \hat{u}_0(\xi), \quad \hat{U}_n^0(\xi) = \hat{\mathcal{G}}_{\Delta T}(\delta(\xi), U_{n-1}^0(\xi)) \quad \text{for } n = 1, \dots, N, $$

and for $k = 0, \dots, k_{conv}$

$$ \hat{U}_n^{k+1}(\xi) = \hat{\mathcal{G}}_{\Delta T}(\delta(\xi), \hat{U}_{n-1}^{k+1}(\xi) - \hat{\mathcal{G}}_{\Delta T}(\delta(\xi), \hat{U}_{n-1}^k(\xi)) \quad \text{for } n = 1, \dots, N. $$

Here specially, $\hat{\mathcal{F}}$ is taken exact so that $\hat{U}_n(\xi) = \exp(-n\delta(\xi))\hat{U}_0(\xi)$ and for implicit first order Euler scheme $\hat{\mathcal{G}}_{\Delta T}(\xi, \cdot) = (1 + \delta(\xi))^{-1}$. Let $\hat{e}_n^k := \hat{U}_n^k - \hat{U}_n$ and $\Delta\hat{\mathcal{G}} = \hat{\mathcal{F}} - \hat{\mathcal{G}}$. This error satisfies for all $\xi \in \mathbb{Z}^2$

$$ \hat{e}_n^k(\xi) = [ \hat{\mathcal{G}}_{\Delta T}(\delta(\xi), \hat{U}_{n-1}^k(\xi)) - \hat{\mathcal{G}}_{\Delta T}(\delta(\xi), \hat{U}_{n-1}(\xi))] \\ + [ \Delta\hat{\mathcal{G}}_{\Delta T}(\delta(\xi), \hat{U}_{n-1}(\xi)) - \Delta\hat{\mathcal{G}}_{\Delta T}(\delta(\xi), \hat{U}_{n-1}^{k-1}(\xi))] \quad (52) $$

We introduce the stability function $R$

$$ R_{k,n}(\delta) = |\hat{\mathcal{G}}_{\Delta T}(\delta, \hat{U}_{n-1}^k(\xi))|^{n-k} |\hat{\mathcal{F}}_{\Delta T}(\delta, \hat{U}_{n-1}^k(\xi)) - \hat{\mathcal{G}}_{\Delta T}(\delta, \hat{U}_{n-1}^k(\xi))|^{k} \begin{pmatrix} n \\ k-1 \end{pmatrix}. \quad (53) $$

The parareal algorithm replaces an algorithm of order $m$ by an algorithm of order $km$ and requires to solve a fine scale problem $k$ times to obtain this accuracy. The stability is then stated in the following theorem (see Theorem 2 in [3]).

**Theorem 1.** Let us assume that the coarse scheme is an approximation of first order of the exact propagator and that it is dissipative in the sense that there exist two constants $C$ and $\gamma > 0$ such that for all $\xi$

$$ \begin{aligned} |\mathcal{F}_{\Delta T}(\delta(\xi)) - \mathcal{G}_{\Delta T}(\delta(\xi))| &\le C(|\delta(\xi)|^2 \wedge 1), \\ |\mathcal{G}_{\Delta T}(\delta(\xi))| &\le (1 + C\Delta T) \exp(-\gamma(|\delta(\xi)| \wedge 1)). \end{aligned} \quad (54) $$

Then the parallel algorithm is stable in the sense that $R_{k,n}(\delta)$ given by (53) is bounded uniformly in $k$ and $\delta = \delta(\xi)$ for $\xi \in \mathbb{Z}^2$ where $a \wedge b = \min(a, b)$.

The one dimensional results (see Theorem 3 in [3]) have been generalized to two dimensions and we have the following result of convergence at final time $T = 1$ by using the norm $H^1(\Omega)$ in [10]

$$ \| u \|_{H^1(\Omega)} = \sum_{\xi_1 = -\infty}^{+\infty} \sum_{\xi_2 = -\infty}^{+\infty} |\hat{u}|^2 (1 + \xi_1^2 + \xi_2^2). $$

**Theorem 2.** Under the assumptions of theorem 1 and if $u_0$ is periodic, the error $e_n^k$ given by (52) is bounded by

$$ \| e_n^k \|_{H^1(\Omega)} \le C \Delta T^k \| u_0 \|_{H^1(\Omega)}. \quad (55) $$
## 5. Parareal-Non-Overlapping Domain Decomposition Method

We use the Neumann-Neumann algorithm for the approximate resolution of the fine propagator independently over each block $\Omega_i \times [T_{n-1}, T_n]$, $i = 1, \dots, N_s$, $n = 0, \dots, N$. The coupled parareal is based on the dual iterative definition of the snapshots $U_n^k$ at each time $T_n$ and on the interface condition $\gamma_{n-1}^{k,0} := u_{n-1,\Gamma}^{k,0}$. The algorithm is thus expressed as follows

$$U_n^{k+1} = \mathcal{G}_{\Delta T}(U_{n-1}^{k+1}) + \tilde{\mathcal{F}}_{\Delta T}(U_{n-1}^k, \gamma_{n-1}^{k,0}) - \mathcal{G}_{\Delta T}(U_{n-1}^k), \quad (56)$$

where $\tilde{\mathcal{F}}_{\Delta T}(U_{n-1}^k, \gamma_{n-1}^{k,0})$ is the approximate propagator with only some P domain decomposition iterations.

More precisely we construct $\tilde{\mathcal{F}}_{\Delta T}$ as follows: we denote by $u_{i,n-1}^{k,P}(x,t)$ the local solutions on the space-time slices $\Omega_i \times [T_{n-1}, T_n]$ at iteration $k$. The fine solution approximated by the algorithm is built, for every fixed $k$ and in each subdomain $\Omega_i$, in several steps by a propagation over $\Omega_i \times [T_{n-1}, T_n]$, an update of the boundary conditions and a new propagation over $\Omega_i \times [T_{n-1}, T_n]$ using the new boundary conditions. The definition of the fine solver is then completed by the computation of the new global solution over $\Omega \times [T_{n-1}, T_n]$. More precisely, for $p = 1, \dots, P$ we solve

$$\left\{ \begin{array}{ll} \bar{\mathcal{L}}(u_{i,n-1}^{k,p}) & = 0 \quad \text{in } \Omega_i \times [T_{n-1}, T_n[ \\ u_{i,n-1}^{k,p}(x,t) & = \gamma_{n-1}^{k,p-1}(t) \quad \text{on } \Gamma \times [T_{n-1}, T_n[ \\ u_{i,n-1}^{k,p}(0,t) & = 0 \quad \text{in } ]T_{n-1}, T_n[ \end{array} \right. $$

$$\left\{ \begin{array}{ll} \bar{\mathcal{L}}(\psi_{i,n-1}^{k,p}) & = 0 \quad \text{in } \Omega_i \times [T_{n-1}, T_n[ \\ \psi_{i,n-1}^{k,p}(x,t) & = \frac{1}{2} \left( \frac{\partial u_{1,n-1}^{k,p}}{\partial \mathbf{n}_1} + \frac{\partial u_{2,n+1}^{k,p}}{\partial \mathbf{n}_2} \right) \quad \text{on } \Gamma \times [T_{n-1}, T_n[ \\ \psi_{i,n-1}^{k,p}(x,t) & = 0 \quad \text{in } ]T_{n-1}, T_n[. \end{array} \right. $$

Then we define the global fine solution by

$$\tilde{\mathcal{F}}_{\Delta T}(U_{n-1}^k, u_{n-1,\Gamma}^{k+1,0}) = \sum_{i=1}^{N_s} \chi^i(x) u_{i,n-1}^{k,p}(x, T_n),$$

where $(\chi^i)_{1 \le i \le n}$ is a regular partition of unity defined by

$$\forall x \in \Omega, \sum_{i=1}^{N_s} \chi^i(x) = 1, \chi^i \ge 0, i = 1, \dots, N_s.$$

### Remark

To accelerate the convergence, we initialized the interface conditions by the previous iterate parareal after a few iterations of domain decomposition methods

$$u_{n-1,\Gamma}^{k,0} = u_{n-1,\Gamma}^{k-1,P},$$

where the boundary condition at the interface is compatible with the initial condition provided at time $T_{n-1}$ for each local problem.
For $P \to \infty$ we have $u_{n-1,\Gamma}^{k-1,\infty}$ verify the exact interface conditions associated with the initial condition $U_{n-1}^{k-1}$ which then $\tilde{\mathcal{F}}_{\Delta T}(U_{n-1}^k, u_{n-1,\Gamma}^{k,\infty}) = \mathcal{F}_{\Delta T}(U_{n-1}^k)$ we obtain the standard parareal algorithm (47).

The parareal-Neumann-Neumann algorithm for the problem (1) can be written as follows (see Algorithm 1):

<caption>**Algorithm 1:** Parareal-Neumann-Neumann algorithm</caption>

```
U_0^0 := u_0;

epsilon a tolerance of the algorithm;

**for** n ← 1 **to** N **do**
    $\gamma_{n-1}^{0,0} := u_0$;
    Step 1. solve $\mathcal{G}$ to obtain the initial guess $U_n^0$
    $Un^0 := \mathcal{G}_{\Delta T}(U_{n-1}^0)$
    **end**
    **repeat**
        Step 2. solve in parallel $\tilde{\mathcal{F}}$ by Neumann-Neumann method
        **parallel for** n ← 1 **to** N **do**
            **parallel for** i ← 1 **to** Ns **do**
                $u_{i,n}^k := \tilde{\mathcal{F}}_{\Delta T}(U_{n-1}^0, \gamma_{n-1}^{k,P})$
            **end for**
            **end for**
        Step 3. initialize the interface condition by the previous iterate parareal after a few iterations of Neumann-Neumann method
        **for** n ← 1 **to** N **do**
            $\gamma_{n-1}^{k+1,0} = \gamma_{n-1}^{k,P}$
        **end**
        Step 4. compute the new iterate $U_n^{k+1}$, n = 1...N by a series of correction iterations
        **for** n ← 1 **to** N **do**
            $U_n^{k+1} = \mathcal{G}_{\Delta T}(U_{n-1}^{k+1}) + \tilde{\mathcal{F}}_{\Delta T}(U_{n-1}^k, \gamma_{n-1}^{k,0}) - \mathcal{G}_{\Delta T}(U_{n-1}^k)$;
        **end**
        $k := k + 1$
    **until** $\|U_n^{k+1} - U_n^k\| < \epsilon$;
```

## 5.1. Convergence analysis

We present now a convergence result of the coupled parareal algorithm.

**Theorem 3.** Under the hypothesis (54) and assuming that

$$ \rho^P \approx \frac{1}{N^{\frac{5}{2}}}, \qquad (57) $$

The parareal algorithm (56) applied to the one dimensional convection diffusion problem (1) with $u_0$ in $H^1(\Omega)$ satisfies the estimate

$$ \| U_n^k(x) - u_n(x) \|_{H^1(\Omega)} \le C \Delta T^{k+1} + \rho^P \Delta T^{k-1} \| u_0 \|_{H^1(\Omega)}. \quad (58) $$
*Proof.* The proof of this theorem will appear in [8]. The estimations of the error $v(x) = U_n^k(x) - u_n(x)$ provided from the convergence rate $\Delta T^{k+1}$ of standard parareal algorithm and from the convergence factor $\rho$ of the domain decomposition method. We present however a detailed numerical study of the convergence in the following section. $\square$

# 6. Numerical Results

We present now numerical experiments to study the theoretical convergence factor $\rho$ according to Lemma (2) of the discretized Neumann-Neumann algorithm for the model problem on the unit square in two dimensions with the initial condition $u_0(x) = \sin(\pi x) \sin(\pi y)$, advection velocity $\alpha = (\alpha_1, \alpha_2)$ and $f(t) = 1 + t$. We choose $\nu = 0.5$, $\alpha_1 = 1$ and $\alpha_2 = 0$. The prototype implementation was done in FreeFem++, The library MPI is used for the parallelization of the code. We discretize in space using finite elements $P^1$-Lagrange on a uniform mesh with meshsize $h$. We use the implicit time discretization of the advection-dominated diffusion equation (1). We decompose the domain into two non-overlapping subdomains $\Omega_1 = [0, 0.5] \times [0, 1]$ and $\Omega_2 = [0.5, ] \times [0, 1]$. For determining the values of interface conditions for domain decomposition algorithms the GMRES method is applied. We test the Neumann-Neumann algorithm by choosing the interface condition $\gamma^{0,0} = u_{0,\Gamma}$.

## 6.1. Convergence of the Neumann-Neumann algorithm

We wish to study the convergence factor. From Figure 1 to Figure 4, we show the convergence factor behavior for different values of mesh size $h$ for Schur complement algorithm and for the best preconditionner Neumann-Neumann for different time-steps $\Delta t$. One can clearly see that the behavior of the algorithms depend on time step. We observe in Figure 2 and Figure 4 that the Neumann-Neumann algorithm converges as fast as the Schur complement algorithm for bigger time step. We deduce that the Schur complement algorithm is much less sensitive with respect the mesh size. The use of parallel preconditioner Neumann-Neumann propagator is necessary to accelerate the convergence. From Figure 1 and Figure 3 we observe convergence results for no preconditioner Schur complement which shows a strong mesh size $h$-dependence.

## 6.2. Convergence of the coupled parareal-Neumann-Neumann algorithm

We now focus on the coupled parareal algorithm. We choose for the spatial discretization step $h = \frac{1}{128}$, for the coarse time step $\Delta T = \frac{1}{10}$ and for the fine time step $\delta t = \frac{\Delta T}{50}$. For the Neumann-Neumann algorithm we make the number GMRES iterations constant (2 iterations per subdomain) for fine time step $\Delta t = \delta t$. We initialize the interface condi-
<figure><img src="image_1.png" /><figcaption>**Figure 2.** The theoretical convergence factor $\rho$ of Neumann-Neumann preconditioner for different values of mesh size $h = 1/8, 1/16, 1/32$ and $\Delta t$ varying through $[10^{-3}, 10^{-2}]$.</figcaption></figure>

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 4.** The theoretical convergence factor $\rho$ of Neumann-Neumann preconditioner for different values of mesh size $h=1/64, 1/128$ and $\Delta t$ varying through $[10^{-3}, 10^{-2}]$.</figcaption></figure>

<figure><img src="image_4.png" /></figure>
tion by $\gamma_{n}^{k,0} = \gamma_{n}^{k-1,2}$. To test the convergence of the solution obtained from the coupled parareal algorithm to the monodomain fine solution $\mathcal{F}$, we compare the relative errors

$$ \frac{\|U_{\text{exact}}(T_n) - \mathcal{F}_{n\Delta T}(u_0)\|_{H^1(\Omega)}}{\|U_{\text{exact}}(T_n)\|_{H^1(\Omega)}} \qquad (59) $$

with

$$ \frac{\|U_{\text{exact}}(T_n) - U_n^k\|_{H^1(\Omega)}}{\|U_{\text{exact}}(T_n)\|_{H^1(\Omega)}} \qquad (60) $$

where $U_{\text{exact}}$ is supposed to be the exact solution (reference solution) computed in the global domain with a very fine time step $\delta t/10$. Table 1 illustrates the errors between the exact solution and the pure sequential fine solver $\mathcal{F}$ at the times $T_n, 1 \le n \le 10$.

<table><caption>**Table 1.** Error between a fine solution and the exact solution obtained with a very fine time step $\delta t/10$ at the times $T_n, 1 \le n \le 10$.</caption><thead><tr><th>T<sub>n</sub></th><th>||U<sub>exact</sub>(T<sub>n</sub>) - F<sub>nΔT</sub>(u<sub>0</sub>)||<sub>H<sup>1</sup>(Ω)</sub> / ||U<sub>exact</sub>(T<sub>n</sub>)||<sub>H<sup>1</sup>(Ω)</sub></th></tr></thead><tbody><tr><td>0.1</td><td>1.85 10<sup>-4</sup></td></tr><tr><td>0.2</td><td>8.80 10<sup>-5</sup></td></tr><tr><td>0.3</td><td>1.16 10<sup>-4</sup></td></tr><tr><td>0.4</td><td>1.63 10<sup>-4</sup></td></tr><tr><td>0.5</td><td>1.91 10<sup>-4</sup></td></tr><tr><td>0.6</td><td>2.11 10<sup>-4</sup></td></tr><tr><td>0.7</td><td>2.28 10<sup>-4</sup></td></tr><tr><td>0.8</td><td>2.43 10<sup>-4</sup></td></tr><tr><td>0.9</td><td>2.58 10<sup>-4</sup></td></tr><tr><td>1</td><td>2.73 10<sup>-4</sup></td></tr></tbody></table>

<table><caption>**Table 2.** Coupled parareal: Historic of the errors between the exact solution and solution of coupled parareal for different values of iteration of parareal $k = 1, 2, 3, 4$ for $P = 2$.</caption><thead><tr><th>T<sub>n</sub></th><th>k=1</th><th>k=2</th><th>k=3</th><th>k=4</th></tr></thead><tbody><tr><td>0.1</td><td>4.43 10<sup>-3</sup></td><td>1.36 10<sup>-3</sup></td><td>1.92 10<sup>-4</sup></td><td>1.72 10<sup>-4</sup></td></tr><tr><td>0.2</td><td>2.37 10<sup>-3</sup></td><td>1.11 10<sup>-3</sup></td><td>1.13 10<sup>-4</sup></td><td>9.14 10<sup>-5</sup></td></tr><tr><td>0.3</td><td>1.21 10<sup>-3</sup></td><td>1.17 10<sup>-3</sup></td><td>1.19 10<sup>-3</sup></td><td>1.18 10<sup>-3</sup></td></tr><tr><td>0.4</td><td>1.42 10<sup>-3</sup></td><td>1.29 10<sup>-3</sup></td><td>1.28 10<sup>-3</sup></td><td>1.26 10<sup>-3</sup></td></tr><tr><td>0.5</td><td>1.76 10<sup>-3</sup></td><td>1.39 10<sup>-3</sup></td><td>1.37 10<sup>-3</sup></td><td>1.32 10<sup>-3</sup></td></tr><tr><td>0.6</td><td>1.97 10<sup>-3</sup></td><td>1.50 10<sup>-3</sup></td><td>1.47 10<sup>-3</sup></td><td>2.02 10<sup>-4</sup></td></tr><tr><td>0.7</td><td>2.11 10<sup>-3</sup></td><td>1.60 10<sup>-3</sup></td><td>1.57 10<sup>-3</sup></td><td>2.07 10<sup>-4</sup></td></tr><tr><td>0.8</td><td>2.22 10<sup>-3</sup></td><td>1.69 10<sup>-3</sup></td><td>3.66 10<sup>-4</sup></td><td>2.12 10<sup>-4</sup></td></tr><tr><td>0.9</td><td>2.31 10<sup>-3</sup></td><td>1.79 10<sup>-3</sup></td><td>3.76 10<sup>-4</sup></td><td>2.18 10<sup>-4</sup></td></tr><tr><td>1</td><td>2.40 10<sup>-3</sup></td><td>1.89 10<sup>-3</sup></td><td>3.86 10<sup>-4</sup></td><td>2.25 10<sup>-4</sup></td></tr></tbody></table>
The Table 2 represent the error between the coupled parareal obtained after two iterations of Neumann-Neumann algorithm and the numerical solution computed in the global domain. Table 2 shows the convergence of the coupled algorithm after $k = 4$. The conclusion reached from these experiment is that with few iterations in each time step are sufficient to obtain convergence. The measured convergence factor of GMRES is $\rho^P = 1.078 \ 10^{-2}$, we notice that at convergence for modified algorithm it is sufficient that $\rho^P \approx C\Delta T^2$. The results confirm the theoretical hypothesis (57) derived in Theorem (6.3).

### 6.3. Experiments with more than two subdomains

These experiments are of interest to generalize the convergence results for two subdomains to many subdomains. We consider the same model problem with the same parameters. We decompose the global domain $\Omega = [0, 1] \times [0, 1]$ into $N_s = 4$ non-overlapping subdomains with cross-points, namely interface nodes that belong to more than two domains, $\Omega_1 = [0, 0.5] \times [0, 0.5]$, $\Omega_2 = [0.5, 1] \times [0, 0.5]$, $\Omega_3 = [0, 0.5] \times [0.5, 1]$ and $\Omega_4 = [0.5, 1] \times [0.5, 1]$. We fixed two iterations of GMRES. We compare the maximum error in time of (59) with the error (60), the maximum error defined in the global domain is $2.73 \ 10^{-4}$. We report in Table 3 the maximum error in time (60).

We note that convergence can be achieved in a finite number of a GMRES in two iterations only after 4 iterations of parareal. The measured convergence factor of GMRES $\rho^2 = 1.257 \ 10^{-2}$. The results confirm the theoretical hypothesis in Theorem. Furthermore the convergence is independent of the number of subdomains.

<table><caption>**Table 3.** Historic of the error using the norm $L^\infty(0, T; H^1(\Omega))$ between parareal-domain decomposition solutions and the reference one, for different values of $k$</caption><thead><tr><th>$k=1$</th><th>$k=2$</th><th>$k=3$</th><th>$k=4$</th><th>$k=5$</th></tr></thead><tbody><tr><td>$5.39 \ 10^{-2}$</td><td>$4.42 \ 10^{-3}$</td><td>$3.05 \ 10^{-4}$</td><td>$2.74 \ 10^{-4}$</td><td>$2.74 \ 10^{-4}$</td></tr></tbody></table>

**Numerical scalability of the parallel** Finally we investigate the numerical scalability of the parallel with respect to the coarse time step $\Delta T$. We split the global domain $\Omega$ into $N_s = 8$ from $\Omega_1 = [0, 0.25] \times [0, 0.5]$ to $\Omega_8 = [0.75, 1] \times [0.5, 1]$. We use two iterations of GMRES, next we fix the fine time step $\delta t = \frac{1}{6000}$ and vary the coarse time step $\Delta T$. The $log_{10}$ scale error between the global fine solution and the reference one is -3.38. From Table 4, we observe that the convergence can be achieved after $k = 4$. One can see that when increasing the number of time decompositions subdomains parareal is barely affected.

## 7. Conclusion

The purpose of this paper was to develop an efficient algorithm for advection-dominated diffusion problem by combining the parareal algorithm with Neumann-Neumann domain
<table><caption>**Table 4.** Historic of the error using the norm $L^\infty(0, T; H^1(\Omega))$ in $\log_{10}$ scale between parareal-domain decomposition solutions and the reference one, for different values of $k$, decomposition algorithm, in order to perform the parallelism both in time and space. We have studied more precisely the convergence factor of Neumann-Neumann algorithm, the main properties used in the convergence analysis is a Fourier analysis and a maximum principle for the 2D case and a two subdomain decomposition. Numerical results show that, with only two iterations of GMRES with Neumann-Neumann preconditioner in each parareal iteration, the new algorithm leads to an accurate monodomain solution.</caption><thead><tr><th>k</th><th>ΔT</th><th>ΔT/2</th><th>ΔT/4</th><th>ΔT/8</th></tr></thead><tbody><tr><td>1</td><td>-1.16</td><td>-1.57</td><td>-2.05</td><td>-2.1</td></tr><tr><td>2</td><td>-1.82</td><td>-2.26</td><td>-2.74</td><td>-2.85</td></tr><tr><td>3</td><td>-2.55</td><td>-3.03</td><td>-3.37</td><td>-3.37</td></tr><tr><td>4</td><td>-3.25</td><td>-3.31</td><td>-3.37</td><td>-3.37</td></tr><tr><td>5</td><td>-3.37</td><td>-3.38</td><td>-3.37</td><td>-3.37</td></tr></tbody></table>

## 8. References

[1] Y. ACHDOU, F. NATAF, "A Robin-Robin preconditioner for an advection-diffusion problem, In Domain decomposition methods, 10 (Boulder, CO, 1997), pages 377-383. Amer. Math. Soc., Providence, RI, 1998.

[2] L. BAFFICO, S. BERNARD, MADAY, Y., TURINICI, G. & ZÉRAH, G. "Parallel-in-time molecular-dynamics simulations, Phys. Rev. E, vol. 66, num. 057701, 2002.

[3] G. BAL, Y. MADAY, "On the convergence and the stability of the parareal algorithm to solve partial differential equations", C. R. Acad. Sci. Paris Sér. I Math., vol. 40 pages 425-432, 2005.

[4] H. BREZIS, "Analyse fonctionnelle: Théorie et applications Collection: Sciences Sup, Dunod 2005 - 248 pages.

[5] X.-C. CAI, "Additive Schwarz algorithms for parabolic convection- diffusion equations", Numerische Mathematik, vol. 60 num. 1 pages 41-61, 1991.

[6] C. FARHAT, "A Lagrange multiplier based divide and conquer finite element algorithm", Journal of Computer and Systems Engineering, num. 2 pages 149-156, 1991.

[7] M. J. GANDER, S. VANDEWALLE, "Analysis of the parareal time-parallel time-integration method", SIAM J. Sci. Comput., vol. 29 pages 556-578, 2007.

[8] R. GUETAT "Méthode de parallélisation en temps: Application aux méthodes de décomposition de domaine", PhD thesis, UPMC & EPT, 203 pages, 2011.

[9] P. LE TALLEC, "Domain decomposition methods in computational mechanics", in: Oden (J. T.), Computational Mechanics Advances vol. 1 num. 2, 121-220, 1994.

[10] G. LEBORGNE "Complément: espaces de Sobolev fractionnaires", Notes du cours d'Équations aux Dérivées Partielles de l'ISIMA, deuxième année http://www.isima.fr/leborgne, 2007.

[11] J.-L. LIONS Y. MADAY, G. TURINICI, "Résolution d'EDP par un schéma en temps pararéel", C. R. Acad. Sci. Paris Sér. I Math., vol. 332 num. 7 pages 661-668, 2001.

[12] Y. MADAY, G. TURINICI, "A parareal in time procedure for the control of partial differential equation, C.R. Acad. Sci. Paris, Ser.Ivol. 335 pages 387-392, 2002.
[13] J. MANDEL “Balancing domain decomposition”, *Communications in Numerical Methods in Engineering*, vol. 233 num. 9(3), 1993.

[14] A. QUARTERONI, A. VALLI “Domain decomposition methods for partial differential equations”, *Oxford University Press*. 2000.

[15] Y. SAAD, H. SCHULTZ, “GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems, *SIAM J. Sci. Stat. Comput.* 7 (3) pages 856-869, 1986.

[16] A. TOSELLI, O. WIDLUNG, “Domain Decomposition Methods - Algorithms and Theory”, *Springer Series in Computational Mathematics*. Springer, vol. 20, 2005.

[17] J. XU, X.-C. CAI, “A preconditioned gmres method for nonsymmetric or indefinite problems, *Math. Comp.*59, 311–319, 1992.