# Optimized Schwarz waveform relaxation method for the incompressible Stokes problem

Duc-Quang Bui, Caroline Japhet, Pascal Omnes

► To cite this version:

Duc-Quang Bui, Caroline Japhet, Pascal Omnes. Optimized Schwarz waveform relaxation method for the incompressible Stokes problem. ESAIM: Mathematical Modelling and Numerical Analysis, 2024, 58 (4), pp.1229-1261. 10.1051/m2an/2024020. hal-04649051

HAL Id: hal-04649051

https://hal.science/hal-04649051v1

Submitted on 16 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# OPTIMIZED SCHWARZ WAVEFORM RELAXATION METHOD FOR THE INCOMPRESSIBLE STOKES PROBLEM

DUC-QUANG BUI<fnref n="1" />, CAROLINE JAPHET<fnref n="1" /> AND PASCAL OMNES<fnref n="1" />,<fnref n="2" />,*

**Abstract.** We propose and analyse the optimized Schwarz waveform relaxation (OSWR) method for the unsteady incompressible Stokes equations. Well-posedness of the local subdomain problems with Robin boundary conditions is proved. Convergence of the velocity is shown through energy estimates; however, pressure converges only up to constant values in the subdomains, and an astute correction technique is proposed to recover these constants from the velocity. The convergence factor of the OSWR algorithm is obtained through a Fourier analysis, and allows to efficiently optimize the space-time Robin transmission conditions involved in the OSWR method. Then, numerical illustrations for the two-dimensional unsteady incompressible Stokes system are presented to illustrate the performance of the OSWR algorithm.

**Mathematics Subject Classification.** 65M55, 35K45, 76D07, 65M12, 65M22, 65B99.

Received May 24, 2023. Accepted March 13, 2024.

## 1. INTRODUCTION

The study of physical phenomena, whether natural or industrial, is frequently based on numerical simulations involving an increasing number of degrees of freedom. This growing complexity may require the use of resolution techniques which on the one hand are suitable for parallel computing architectures, and on the other hand allow local space and time stepping adapted to the physics, such as space-time domain decomposition (DD) methods. In this article we are concerned with such methods, with Robin transmission conditions at the interfaces between subdomains, for solving applications related to incompressible fluid mechanics, that are modelled by the unsteady (Navier)-Stokes system.

The well-posedness of such systems with Robin conditions (without domain decomposition) has been the subject of several works in the steady case, see e.g. [48] for the Stokes problem (where the Robin condition is expressed with the symmetric part of the velocity gradient, instead of the gradient), references [39, 46] for the Oseen and Navier-Stokes systems, and [17] for the Stokes-Darcy Coupling. On the other hand, there are few works in the unsteady case; in [40] existence and uniqueness of a solution with a time-dependent Robin boundary condition of the type curl **u** × **n** = β(t) **u** is addressed. In [29] the Stokes problem with Robin conditions is

<footnote>*Keywords and phrases.* Unsteady incompressible Stokes system, space-time domain decomposition, optimized Schwarz waveform relaxation, Robin transmission conditions, correction technique for the pressure.</footnote>

<footnote n="1">¹ Université Sorbonne Paris Nord, LAGA, CNRS, UMR 7539, Institut Galilée, 99 Av. J.-B. Clément, 93430 Villetaneuse, France.</footnote>

<footnote n="2">² Université Paris-Saclay, CEA, Service de Génie Logiciel pour la Simulation, 91191 Gif-sur-Yvette, France.</footnote>

<footnote>*Corresponding author: pascal.omnes@cea.fr</footnote>
studied, in the context of a global-in-time DD method applied the coupled nonlinear Stokes and Darcy Flows. The well-posedness is not shown.

In this article we study the well-posedness of the unsteady incompressible Stokes system with Robin boundary conditions of type $\alpha(\nu\partial_n\mathbf{u} \cdot \mathbf{n} - p) + \mathbf{u} \cdot \mathbf{n} = g(t)$ and $\beta\nu\partial_n\mathbf{u} \times \mathbf{n} + \mathbf{u} \times \mathbf{n} = \xi(t)$, in the context of space-time DD methods.

Concerning the DD approaches with Robin conditions, several studies have been carried out for the incompressible (Navier)-Stokes equations: in [35, 41–44] the steady Oseen equation (and its application to the non-stationary Navier-Stokes equations, using a spatial DD at each time step) is considered. More precisely, in [35, 43, 44] a stabilized finite element approximation is proposed (with non-standard Robin conditions due to the stabilization). The convergence of the DD method is proven for the velocity. For the pressure, the convergence is proven when the original monodomain problem involves Robin boundary conditions on a part of the physical boundary. However, the authors point out that for an Oseen problem with Dirichlet conditions on the whole physical boundary, the pressure of the Robin-Robin DD algorithm will converge up to a constant which can differ for different subdomains. This important observation is also mentioned in [12] for the steady Stokes problem, where the DD method is based on a penalty term on the interface (in that case the Robin conditions are not equivalent to the physical ones). The convergence is shown for a modified pressure in the two-subdomains case. This issue of pressure converging up to a constant that depends on the subdomains is also raised in [24, 34] for the discrete Schwarz algorithm with a DDFV scheme applied to the semi-discrete in time Navier-Stokes system. In [7, 13], an optimized Schwarz DD method is studied, and applied at each time step to the semi-discrete in time Navier-Stokes equations. Other transmission conditions (Dirichlet/Neumann) are considered e.g. in [22, 45, 47, 50] for Stokes and Navier-Stokes equations.

In this article we consider global-in-time Schwarz methods which use waveform relaxation techniques, i.e. Schwarz waveform relaxation (SWR). Such iterative methods use computations in the subdomains over the whole time interval, exchanging space-time boundary data through transmission conditions on the space-time interfaces. The main advantage is that space-time discretizations can be chosen independently on each subdomain, and, at the end of each iteration, only a small amount of information is exchanged, which makes the parallelization (in space and time) very efficient.

The space-time boundary data play an important role in the convergence process and can be of Dirichlet [19, 23], absorbing, Robin (or Ventcell) type [4, 21, 25, 26, 36]. The value of the Robin (or Ventcell) parameters can be optimized to improve convergence rates (see [21, 31, 33, 36]), and the corresponding method is called optimized Schwarz waveform relaxation (OSWR). This method is widely used and analyzed for fluid dynamics, see references above, and e.g. [1, 3, 6, 20, 30, 36, 37, 49].

For the application of the SWR method on the Navier-Stokes equations, we are aware of the article [3] where an OSWR method is proposed for the rotating 3D incompressible hydrostatic Navier-Stokes equations with free surface. However, the hydrostatic nature of the model modifies the structure of the continuity equation which now involves a transport term for the free surface (which plays the same role as the pressure in the momentum equation of the standard Navier-Stokes system), so that the results in [3] cannot apply to the problem considered in the present work. In [13], an SWR method for the Oseen equations is studied; optimal transparent boundary conditions are derived, and local approximations for these nonlocal conditions are proposed. No general convergence analysis of the resulting algorithm (e.g. via energy estimates) is given. A convergence factor is obtained in the idealized case of two half-space subdomains and unbounded time interval, via Laplace-Fourier transforms.

Concerning the compressible Euler and Navier-Stokes equations, in [14, 15] an SWR method is proposed and various numerical experiments are shown.

However, until now, there exists no convergence proof (for SWR or OSWR) for the incompressible Navier-Stokes equations. We contribute to the understanding of the behaviour of the OSWR method by attacking representative, though simpler, model problems. To begin with, we analyze the method on the evolutionary Stokes equations, a simplified version of the evolutionary Navier-Stokes system in which the convection is simply discarded. The convergence analysis of the velocity iterates involved in the OSWR method, for the
Stokes equations, can be performed in a similar manner as for parabolic equations. An extension of this analysis to the evolutionary Oseen equations (a linearization of the Navier-Stokes equations in which the convective velocity field is considered as a given datum) is given in [10]. However, the convergence analysis of the OSWR method has its own obstacle related to the pressure converging only up to constants in the various subdomains, as discussed above. A second purpose of this article is to propose a new technique, in the multidomain case, to recover the pressure from the velocity (at any iteration).

A third purpose of this article is to discuss the choice of the Robin parameters, which play a crucial role in the optimization of the convergence rate. Until recently, the common practice was to derive and optimize a convergence rate in the idealized case of two half-space subdomains and unbounded time interval, via Laplace-Fourier transforms performed on the continuous model (i.e. without taking into account the actual discretization method). We first follow this standard approach in this work, but in a second step modify it to also include the effect of the discretization in the time direction; the Robin parameters obtained with such a modification improve the convergence rate over the standard choice in our numerical tests. Note that studying the influence of the numerical scheme over the OSWR convergence rate is a recent approach, pursued for example in [2,16,27].

The remainder of this article is organized as follows. In Section 2, we present the model problem and its multidomain form. Since the multi-domain formulation involves local Stokes problems with Robin boundary conditions, we prove the well-posedness of such problems in Section 3. Next, Section 4 is dedicated to the algorithm. In Section 5 we show that, in general, the pressure calculated by the OSWR algorithm will not converge to the monodomain solution. In Section 6, we obtain a convergence result on the velocity through an energy estimate, and in Section 7, we propose an astute technique to recover the pressure from the velocity. In Section 8, a Fourier analysis is done to get a formulation for the convergence factor of the OSWR algorithm. In Section 9, an optimization procedure (based on the convergence factor of the method), that allows to obtain efficient Robin parameters, is given. Then, numerical illustrations for the unsteady Stokes system follow in Section 10.

## 2. PRESENTATION OF THE MODEL AND MULTIDOMAIN FORMULATION

For a bounded domain $\Omega \subseteq \mathbb{R}^2$, and for a given viscosity coefficient $\nu > 0$ that we suppose constant and uniform, for given initial condition $u_0$ and source term $f$, we denote respectively by $u, p$ the velocity and pressure unknowns in the incompressible non-stationary Stokes system:

$$ 
\begin{array}{lll}
\partial_t u - \nu \Delta u + \nabla p = \mathbf{f} & \text{in } \Omega \times (0, T), \\
\nabla \cdot \mathbf{u} = 0 & \text{in } \Omega \times (0, T), \\
u(\cdot, t=0) = u_0 & \text{in } \Omega, \\
u = 0 & \text{on } \partial\Omega \times (0, T).
\end{array} 
\quad (1) $$

This system does not have a unique solution: if $(u, p)$ is a solution, then $(u, p + c)$ is also a solution, for any constant $c$. Then, for uniqueness, one needs, for example, the zero-mean condition on the pressure

$$ \int_{\Omega} p = 0. \quad (2) $$

Thus, we introduce the notation $L_0^2(\Omega) = \{p \in L^2(\Omega), \int_{\Omega} p = 0\}$.

Next, we shall introduce the following spaces, which are the completions, in $H^1(\Omega)$ and in $L^2(\Omega)$, respectively, of the set of compactly supported $C^\infty$ functions with vanishing divergence:

$$ V = \{ \mathbf{u} \in [H_0^1(\Omega)]^2, \nabla \cdot \mathbf{u} = 0 \}, $$

$$ H = \{ \mathbf{u} \in [L^2(\Omega)]^2, \nabla \cdot \mathbf{u} = 0, \mathbf{u} \cdot \mathbf{n}_{\partial\Omega} = 0 \text{ on } \partial\Omega \}, $$
where $\mathbf{n}_{\partial\Omega}$ is the unit, outward pointing, normal vector field on $\partial\Omega$. We denote by $V'$ the dual space of $V$ and denote by $\langle \cdot, \cdot \rangle_{V',V}$ the duality bracket between the two spaces. We recall ([8, Proposition IV.5.13]) that, if $\Omega$, $\mathbf{f}$ and $\mathbf{u}_0$ are regular enough, problem (1)–(2) has a unique solution $(\mathbf{u}, p)$ such that

$$ \mathbf{u} \in (L^2(0, T; V) \cap C^0([0, T]; H)), \quad \partial_t \mathbf{u} \in L^2(0, T; V'), \\ p \in W^{-1,\infty}(0, T; L_0^2(\Omega)). $$

In order to apply a domain-decomposition strategy for this problem, we decompose $\Omega$ into $M$ non-overlapping subdomains $\Omega_i$, i.e. $\Omega_i \cap \Omega_j = \emptyset$ for $i \neq j$, and $\bar{\Omega} = \bigcup_{i=1}^M \bar{\Omega}_i$. For $i = 1, 2, \dots, M$, we denote by $\mathcal{I}_i$ the set of indices of the neighbouring subdomain(s) of $\Omega_i$: it holds that $j \in \mathcal{I}_i$ if and only if $|\partial\Omega_i \cap \partial\Omega_j| > 0$, where $|\cdot|$ denotes the one dimensional measure. We denote by $\Gamma_{ij}$ the interface (if it exists) between $\Omega_i$ and $\Omega_j$, $\mathbf{n}_{ij}$ the unit normal vector on $\Gamma_{ij}$, directed from $\Omega_i$ to $\Omega_j$. Note that this implies that $\mathbf{n}_{ij} = -\mathbf{n}_{ji}$.

Denoting by $\mathbf{u}_i$, $(\mathbf{u}_0)_i$, $p_i$ and $\mathbf{f}_i$ the respective restrictions of $\mathbf{u}$, $\mathbf{u}_0$, $p$ and $\mathbf{f}$ to $\Omega_i$, the monodomain problem is equivalent to the following multidomain one

$$ \begin{array}{ll} \partial_t \mathbf{u}_i - \nu \Delta \mathbf{u}_i + \nabla p_i = \mathbf{f}_i & \text{in } \Omega_i \times (0, T), \\ \nabla \cdot \mathbf{u}_i = 0 & \text{in } \Omega_i \times (0, T), \\ \mathbf{u}_i(\cdot, t=0) = (\mathbf{u}_0)_i & \text{in } \Omega_i, \\ \mathbf{u}_i = 0 & \text{on } (\partial\Omega_i \cap \partial\Omega_j) \times (0, T), \end{array} \quad (3) $$

for all $i \in [1, M]$, together with the physical transmission conditions on the space-time interfaces $\Gamma_{ij} \times (0, T)$, $j \in \mathcal{I}_i$, $i \in [1, M]$,

$$ \begin{array}{ll} \mathbf{u}_i \cdot \mathbf{n}_{ij} = -\mathbf{u}_j \cdot \mathbf{n}_{ji}, & \\ \mathbf{u}_i \times \mathbf{n}_{ij} = -\mathbf{u}_j \times \mathbf{n}_{ji}, & \\ \nu \partial_{\mathbf{n}_{ij}} \mathbf{u}_i \cdot \mathbf{n}_{ij} - p_i = \nu \partial_{\mathbf{n}_{ji}} \mathbf{u}_j \cdot \mathbf{n}_{ji} - p_j, & (4) \\ \nu \partial_{\mathbf{n}_{ij}} \mathbf{u}_i \times \mathbf{n}_{ij} = \nu \partial_{\mathbf{n}_{ji}} \mathbf{u}_j \times \mathbf{n}_{ji}. & \end{array} $$

For any choice of $(\alpha_{ij}, \alpha_{ji}, \beta_{ij}, \beta_{ji}) \in (\mathbb{R}^{+*})^4$, those conditions are equivalent to the following Robin transmission conditions on $\Gamma_{ij} \times (0, T) = \Gamma_{ji} \times (0, T)$:

$$ \begin{align} \alpha_{ij}(\nu \partial_{\mathbf{n}_{ij}} \mathbf{u}_i \cdot \mathbf{n}_{ij} - p_i) + \mathbf{u}_i \cdot \mathbf{n}_{ij} &= \alpha_{ij}(\nu \partial_{\mathbf{n}_{ij}} \mathbf{u}_j \cdot \mathbf{n}_{ij} - p_j) + \mathbf{u}_j \cdot \mathbf{n}_{ij}, \tag{5} \\ \alpha_{ji}(\nu \partial_{\mathbf{n}_{ji}} \mathbf{u}_j \cdot \mathbf{n}_{ji} - p_j) + \mathbf{u}_j \cdot \mathbf{n}_{ji} &= \alpha_{ji}(\nu \partial_{\mathbf{n}_{ji}} \mathbf{u}_i \cdot \mathbf{n}_{ji} - p_i) + \mathbf{u}_i \cdot \mathbf{n}_{ji}, \nonumber \\ \beta_{ij} \nu \partial_{\mathbf{n}_{ij}} \mathbf{u}_i \times \mathbf{n}_{ij} + \mathbf{u}_i \times \mathbf{n}_{ij} &= \beta_{ij} \nu \partial_{\mathbf{n}_{ij}} \mathbf{u}_j \times \mathbf{n}_{ij} + \mathbf{u}_j \times \mathbf{n}_{ij}, \nonumber \\ \beta_{ji} \nu \partial_{\mathbf{n}_{ji}} \mathbf{u}_j \times \mathbf{n}_{ji} + \mathbf{u}_j \times \mathbf{n}_{ji} &= \beta_{ji} \nu \partial_{\mathbf{n}_{ji}} \mathbf{u}_i \times \mathbf{n}_{ji} + \mathbf{u}_i \times \mathbf{n}_{ji}. \nonumber \end{align} $$

Finally, the zero-mean condition for the pressure is equivalent to

$$ \sum_{i=1}^{M} \int_{\Omega_i} p_i = 0. \quad (6) $$

This setting requires that we should study the Stokes system in a domain where Robin boundary conditions are applied on a part of the boundary. This is what is done in the next section.

### 3. THE STOKES PROBLEM WITH ROBIN BOUNDARY CONDITIONS

We now consider a domain, still denoted by $\Omega$, for which the boundary is decomposed into two parts: $\partial\Omega = \Gamma_D \cup \Gamma_R$, with $|\Gamma_R| > 0$. Let $\mathbf{n}$ be the outgoing normal vector on $\Gamma_R$; we consider the following system,
with $\alpha > 0$ and $\beta > 0$

$$
\begin{align}
& \partial_t \mathbf{u} - \nu \Delta \mathbf{u} + \nabla p = \mathbf{f} \quad \text{in } \Omega \times (0, T), \nonumber \\
& \nabla \cdot \mathbf{u} = 0 \quad \text{in } \Omega \times (0, T), \nonumber \\
& \mathbf{u}(\cdot, t=0) = \mathbf{u}_0 \quad \text{in } \Omega, \nonumber \\
& \mathbf{u} = 0 \quad \text{on } \Gamma_D \times (0, T), \nonumber \\
& \alpha(\nu \partial_n \mathbf{u} \cdot \mathbf{n} - p) + \mathbf{u} \cdot \mathbf{n} = g \quad \text{on } \Gamma_R \times (0, T), \nonumber \\
& \beta \nu \partial_n \mathbf{u} \times \mathbf{n} + \mathbf{u} \times \mathbf{n} = \xi \quad \text{on } \Gamma_R \times (0, T), \tag{7}
\end{align}
$$

where $\mathbf{f}$ is at least in $[L^2(\Omega \times (0, T))]^2$, $g$ and $\xi$ are at least in $[L^2(\Gamma_R \times (0, T))]$.

In order to set this problem under an appropriate (parabolic) variational form, we multiply the first equation by a divergence-free test function $\mathbf{v}$ (smooth enough) that vanishes on $\Gamma_D$ and integrate by parts on $\Omega$. The flux $(-\nu\partial_n\mathbf{u} + p\mathbf{n})$ is then decomposed into normal and tangential parts and boundary conditions of (7) are used. We obtain then the following parabolic variational problem

$$ \langle \partial_t \mathbf{u}, \mathbf{v} \rangle_{V'_D, V_D} + a(t, \mathbf{u}, \mathbf{v}) = c(t, \mathbf{v}), \quad \text{a.e. } t \in (0, T), \forall \mathbf{v} \in V_D, \quad (8) $$

$$ \mathbf{u}(0) = \mathbf{u}_0, \quad (9) $$

where the spaces are defined as

$$
\begin{align*}
V_D &= \left\{ \mathbf{u} \in [H^1(\Omega)]^2, \mathbf{u} = 0 \text{ on } \Gamma_D, \nabla \cdot \mathbf{u} = 0 \right\}, \\
H_D &= \left\{ \mathbf{u} \in [L^2(\Omega)]^2, \mathbf{u} \cdot \mathbf{n} = 0 \text{ on } \Gamma_D, \nabla \cdot \mathbf{u} = 0 \right\},
\end{align*}
$$

together with their linear and bilinear forms

$$ a(\mathbf{u}, \mathbf{v}) = \nu (\nabla \mathbf{u}, \nabla \mathbf{v})_{\Omega} + \frac{1}{\alpha} (\mathbf{u} \cdot \mathbf{n}, \mathbf{v} \cdot \mathbf{n})_{\Gamma_R} + \frac{1}{\beta} (\mathbf{u} \times \mathbf{n}, \mathbf{v} \times \mathbf{n})_{\Gamma_R}, \quad (10) $$

$$ c(t, \mathbf{v}) = (\mathbf{f}(t), \mathbf{v})_{\Omega} + \frac{1}{\alpha} (g(t), \mathbf{v} \cdot \mathbf{n})_{\Gamma_R} + \frac{1}{\beta} (\xi(t), \mathbf{v} \times \mathbf{n})_{\Gamma_R}. \quad (11) $$

Here, $(\cdot, \cdot)_{D}$ denotes, for any set $D$ (whatever the space-dimension of $D$) the standard scalar or the matrix-valued scalar $L^2$ product on $D$. In the same way, we shall use the notation $\| \cdot \|_D$ for the associated $L^2(D)$ norm. All terms in the definition of the forms $a$ and $c$ are well-defined for $(\mathbf{u}, \mathbf{v}) \in V_D \times V_D$.

From these definitions, $V_D$ is dense in $H_D$ and the embedding $V_D \subset H_D$ is continuous. We can identify $H_D$ with its dual space, and we are in the situation where $V_D \subset H_D \equiv H'_D \subset V'_D$, which is the classical setting for parabolic equations (see e.g. [18], Sect. 6.1, [9], P. 218). In this context, we recall the following theorem Theorem 6.6 from [18].

**Theorem 3.1.** *Problem (8)–(9) admits a unique solution*

$$ \mathbf{u} \in (L^2(0, T; V_D) \cap C^0([0, T]; H_D)), $$

with $\partial_t \mathbf{u} \in L^2(0, T; V'_D)$ if the following properties are verified

- $\mathbf{u}_0 \in H_D$ and $c \in L^2(0, T; V'_D)$,

- The function $t \mapsto a(t, \mathbf{u}, \mathbf{v})$ is measurable for all $(\mathbf{u}, \mathbf{v}) \in V_D^2$,

- $\exists M \in \mathbb{R}$ such that $|a(t, \mathbf{u}, \mathbf{v})| \le M \|\mathbf{u}\|_{V_D} \|\mathbf{v}\|_{V_D}$ for almost every $t$ and for all $(\mathbf{u}, \mathbf{v}) \in V_D^2$,

- $\exists m > 0$ such that $a(t, \mathbf{u}, \mathbf{u}) \ge m \|\mathbf{u}\|_{V_D}^2$ for almost every $t$ and for all $\mathbf{u} \in V_D$.

We shall apply this result to our setting, with the simplification that the bilinear form defined by (10) does not depend on time. We obtain the following result:
**Theorem 3.2.** Assume that $\mathbf{f} \in L^2(0, T; [L^2(\Omega)]^2)$, $g, \xi \in L^2(0, T; L^2(\Gamma_R))$, and $\mathbf{u}_0 \in H_D$. Let $a$ and $c$ be defined by (10) and (11), respectively. Then, problem (8)–(9) admits a unique solution $\mathbf{u} \in L^2(0, T; V_D) \cap C^0([0, T]; H_D)$, which is such that $\partial_t \mathbf{u} \in L^2(0, T; V'_D)$.

*Proof.* We shall show that $a$ and $c$ verify the hypothesis of Theorem 3.1. First, it is well-known that, as soon as $|\Gamma_R| > 0$, then

$$ \|\mathbf{u}\|_{V_D} := (\|\nabla \mathbf{u}\|_{\Omega}^2 + \|\mathbf{u}\|_{\Gamma_R}^2)^{\frac{1}{2}} = (\|\nabla \mathbf{u}\|_{\Omega}^2 + \|\mathbf{u} \cdot \mathbf{n}\|_{\Gamma_R}^2 + \|\mathbf{u} \times \mathbf{n}\|_{\Gamma_R}^2)^{\frac{1}{2}} $$

is a norm equivalent to the $H^1$ norm on $V_D$, and we shall therefore work with this norm.

Let $M = \max(\nu, \frac{1}{\alpha}, \frac{1}{\beta})$. From the Cauchy-Schwarz inequality, we get the continuity of $a(\cdot, \cdot)$:

$$ |a(\mathbf{u}, \mathbf{v})| \leq M \|\mathbf{u}\|_{V_D} \|\mathbf{v}\|_{V_D}, \quad \forall \mathbf{u}, \mathbf{v} \in V_D. $$

Let $m = \min(\nu, \frac{1}{\alpha}, \frac{1}{\beta}) > 0$. From the definition of $\|\cdot\|_{V_D}$, we get the coercivity of $a(\cdot, \cdot)$:

$$ a(\mathbf{u}, \mathbf{u}) \geq m \|\mathbf{u}\|_{V_D}^2, \quad \forall \mathbf{u} \in V_D. $$

Then, for a.e. $t \in (0, T)$, the continuity of $c(t, \cdot)$ is deduced from the Cauchy-Schwarz inequality and the equivalence between the $H^1(\Omega)$-norm and $\|\cdot\|_{V_D}$:

$$ |c(t, \mathbf{v})| \leq \left[ C_1 \|\mathbf{f}(t)\|_{\Omega} + \frac{1}{\alpha} \|g(t)\|_{\Gamma_R} + \frac{1}{\beta} \|\xi(t)\|_{\Gamma_R} \right] \|\mathbf{v}\|_{V_D}. $$

Moreover, thanks to the hypothesis on the time dependence of $\mathbf{f}$, $g$ and $\xi$, the quantity

$$ C_1 \|\mathbf{f}(t)\|_{\Omega} + \frac{1}{\alpha} \|g(t)\|_{\Gamma_R} + \frac{1}{\beta} \|\xi(t)\|_{\Gamma_R} $$

is square integrable on $(0, T)$, and we can now apply Theorem 3.1, which finishes the proof. $\square$

**Remark 3.3.** Since $V_D$ is continuously and densely embedded in $H_D$, the fact that $\mathbf{u} \in C^0([0, T]; H_D)$ is a consequence of the fact that the space

$$ \mathcal{W}(V_D, V'_D) := \{ \mathbf{v} : (0, T) \mapsto V_D; \mathbf{v} \in L^2(0, T; V_D); \partial_t \mathbf{v} \in L^2(0, T; V'_D) \} $$

is included in $C^0([0, T]; H_D)$, as stated, for example, by Lemma 6.2 from [18] and Theorem II.5.13 from [8].

This has the important implication that it is legitimate to consider $\mathbf{u}(t) \in H_D$ for all $t \in [0, T]$. Moreover, the following integral equality holds for all $t \in [0, T]$ and for all $(\mathbf{u}, \mathbf{v}) \in [\mathcal{W}(V_D, V'_D)]^2$ (see [18], Lem. 6.3 and [8],Thm. II.5.12):

$$ \int_0^t (\langle \partial_t \mathbf{u}(s), \mathbf{v}(s) \rangle_{V'_D, V_D} + \langle \partial_t \mathbf{v}(s), \mathbf{u}(s) \rangle_{V'_D, V_D}) ds = (\mathbf{u}(t), \mathbf{v}(t))_{\Omega} - (\mathbf{u}(0), \mathbf{v}(0))_{\Omega}. \quad (12) $$

Now, since we have obtained the velocity $\mathbf{u}$ from the constrained variational problem (8)–(9), we shall construct the pressure by relaxing the divergence free condition on the velocity test functions, and we shall therefore consider the space

$$ X_D = \{ \mathbf{v} \in [H^1(\Omega)]^2, \mathbf{v} = 0 \text{ on } \Gamma_D \}, $$

equipped with the norm $\|\cdot\|_{X_D} := \|\cdot\|_{V_D}$ defined above. Like often with the Stokes problem, we shall rely on the surjectivity of the divergence operator, and on general properties of surjective mappings in Hilbert spaces. More precisely, we shall use the following results.
**Lemma 3.4.** The mapping $B$ from $X_D$ into $L^2(\Omega)$ defined by $B(v) = -\nabla \cdot v$ is continuous and surjective.

*Proof.* This is a special case of Lemma 4.9 from [18] (with, using the notations of [18], $\partial\Omega_1 = \Gamma_D$, $\partial\Omega_2 = \emptyset$, $\partial\Omega_3 = \emptyset$ and $\partial\Omega_4 = \Gamma_R$). □

**Lemma 3.5.** Let $L$ be in $\mathcal{L}(E; F)$ and $L^T$ be its adjoint in $\mathcal{L}(F'; E')$. If $L$ is surjective in $F$, then $\mathrm{Im} \, L^T$ is closed in $E'$.

For the proof of Lemma 3.5, see, e.g., Lemma A.40 from [18]. Next, we recall the following definition (see, e.g., [8], Def. IV.2.1)

**Definition 3.6.** Let $E$ be a Banach space with dual space $E'$; then for any subset $A \subset E$, we define $A^{\perp} \subset E'$ as follows:

$$A^{\perp} := \{\phi \in E', \forall x \in A, \langle\phi, x\rangle_{E',E} = 0\}$$

as well as the following properties stated, e.g., by Remark IV.2.1 from [8]

**Lemma 3.7.** If $A \subset C \subset E$, then $C^{\perp} \subset A^{\perp}$.

**Lemma 3.8.** If $A$ is a linear subspace of $E$, then $(A^{\perp})^{\perp} = A$ if and only if $A$ is closed in $E$.

Moreover, we also recall the following general result

**Lemma 3.9.** Let $L$ be in $\mathcal{L}(E; F)$, then $(\mathrm{Im} \, L^T)^{\perp} \subset \mathrm{Ker} \, L$

*Proof.* If $f \in (\mathrm{Im} \, L^T)^{\perp}$, then $\langle L^T q, f \rangle_{E',E} = 0, \forall q \in F'$. Thus $\langle q, Lf \rangle_{F',F} = 0$ for all $q \in F'$, which means that $Lf = 0$, and thus $f \in \mathrm{Ker} \, L$. □

From these results, we obtain the following Lemma, which will be useful in the construction of the pressure field:

**Lemma 3.10.** Let $B^T$ be the adjoint operator of $B$, from $L^2(\Omega)$ into $X'_D$. Then for any $\phi$ in $X'_D$ that vanishes on $V_D$, there exists $P \in L^2(\Omega)$ such that $\phi = B^T P$.

*Proof.* Since $B$ is in $\mathcal{L}(X_D; L^2(\Omega))$ and is surjective (Lem. 3.4), then $(\mathrm{Im} \, B^T)$ is closed in $X'_D$ (Lem. 3.5), and $((\mathrm{Im} \, B^T)^{\perp})^{\perp} = \mathrm{Im} \, B^T$ (Lem. 3.8). Now, using Lemmas 3.9 and 3.7, we get $(\mathrm{Ker} \, B)^{\perp} \subset ((\mathrm{Im} \, B^T)^{\perp})^{\perp} = \mathrm{Im} \, B^T$. So if $\phi$ in $X'_D$ vanishes on $V_D = \mathrm{Ker} \, B$, then $\phi$ is in $(\mathrm{Ker} \, B)^{\perp}$ and so in $\mathrm{Im} \, B^T$, which exactly means that there exists $P \in L^2(\Omega)$ such that $\phi = B^T P$. □

Using this result, we can now state the following theorem.

**Theorem 3.11.** Assume that $f \in L^2(0, T; [L^2(\Omega)]^2)$, $\xi, g \in L^2(0, T; L^2(\Gamma_R))$ and $u_0 \in H_D$, then there exists unique $u \in (L^2(0, T; V_D) \cap C^0([0, T]; H_D))$ and $p \in W^{-1,\infty}(0, T; L^2(\Omega))$, with $\partial_t u \in L^2(0, T; V'_D)$ such that $(u, p)$ verifies problem (7) in the sense that

* $u$ verifies (8)-(9)

* $p = \partial_t P$ with $P \in L^\infty(0, T; L^2(\Omega))$ that satisfies

$$\int_0^t c(s, \mathbf{v})ds - (\mathbf{u}(t), \mathbf{v})_{\Omega} + (\mathbf{u}_0, \mathbf{v})_{\Omega} - \int_0^t a(\mathbf{u}(s), \mathbf{v})ds = - \int_{\Omega} P(t)\nabla \cdot \mathbf{v}, \forall \mathbf{v} \in X_D. \quad (13)$$
*Proof.* Using Theorem 3.2, the existence and uniqueness of **u** being the solution of (8)–(9) and satisfying the regularity claim of the present theorem is obtained. Thus, we only have to show the existence and uniqueness of the pressure field *p* as defined above. For this, we first show that the linear form over $X_D$ defined by the left-hand side of (13) is continuous on $X_D$ and thus an element of $X'_D$. Secondly, we use Lemma 3.10 to show the existence of $P(t) \in L^2(\Omega)$ such that (13) is verified for all $t \ge 0$. The third step is to prove that $\|P(t)\|_{L^2(\Omega)}$ is bounded with respect to time. Finally, we shall prove that *p* is unique.

First step: Consider, for **u** being the solution of (8)–(9), the function $t \mapsto a(\mathbf{u}(t), \mathbf{v})$ and the function $t \mapsto c(t, \mathbf{v})$ where *a* and *c* are defined by (10) and (11). Then their definitions can be straightforwardly extended to consider $\mathbf{v} \in X_D$ and, for any $t \in (0, T)$, the following element of $X'_D$ is well-defined:

$$b(t, \mathbf{v}) := \int_{0}^{t} c(s, \mathbf{v})ds - (\mathbf{u}(t), \mathbf{v})_{\Omega} + (\mathbf{u}_{0}, \mathbf{v})_{\Omega} - \int_{0}^{t} a(\mathbf{u}(s), \mathbf{v})ds, \quad \forall \mathbf{v} \in X_{D}.$$

Indeed, one has that

$$\begin{aligned} \left| \int_{0}^{t} a(\mathbf{u}(s), \mathbf{v})ds \right| &\le \int_{0}^{t} M \|\mathbf{u}(s)\|_{X_D} \|\mathbf{v}\|_{X_D} ds \\ &\le M \sqrt{t} \left[ \int_{0}^{t} \|\mathbf{u}(s)\|_{X_D}^{2} ds \right]^{\frac{1}{2}} \|\mathbf{v}\|_{X_D} \\ &\le M \sqrt{T} \|\mathbf{u}\|_{L^2(0,T;V_D)} \|\mathbf{v}\|_{X_D}, \end{aligned}$$

and

$$\begin{aligned} \left| \int_{0}^{t} c(s, \mathbf{v})ds \right| &\le \int_{0}^{t} \left( C_1 \|\mathbf{f}(s)\|_{\Omega} + \frac{1}{\alpha} \|g(s)\|_{\Gamma_R} + \frac{1}{\beta} \|\xi(s)\|_{\Gamma_R} \right) \|\mathbf{v}\|_{X_D} \\ &\le \gamma_1 \|\mathbf{v}\|_{X_D}, \end{aligned}$$

with

$$\gamma_1 = C_1 \sqrt{T} \|\mathbf{f}\|_{L^2(0,T;[L^2(\Omega)]^2)} + \frac{\sqrt{T}}{\alpha} \|g\|_{L^2(0,T;L^2(\Gamma_R))} + \frac{\sqrt{T}}{\beta} \|\xi\|_{L^2(0,T;L^2(\Gamma_R))}.$$

In addition, since **u** belongs to $C^0([0, T]; H_D)$, then

$$\begin{aligned} &\left| - (\mathbf{u}(t), \mathbf{v})_{\Omega} + (\mathbf{u}_0, \mathbf{v})_{\Omega} \right| \\ &\le 2 \|\mathbf{u}\|_{L^{\infty}([0,T];[L^2(\Omega)]^2)} \|\mathbf{v}\|_{\Omega} \\ &\le 2C_1 \|\mathbf{u}\|_{L^{\infty}([0,T];[L^2(\Omega)]^2)} \|\mathbf{v}\|_{X_D}. \end{aligned}$$

Defining $C_2 = 2C_1 \|\mathbf{u}\|_{L^{\infty}([0,T];[L^2(\Omega)]^2)} + \gamma_1 + M\sqrt{T} \|\mathbf{u}\|_{L^2(0,T;V_D)}$, this leads to the fact that

$$|b(t, \mathbf{v})| \le C_2 \|\mathbf{v}\|_{X_D}, \quad \forall \mathbf{v} \in X_D, \forall t \in (0, T), \quad (14)$$

and thus proves that the left-hand side of (13) defines and element of $X'_D$, which ends the first step of the proof.

Second step: From (8) and (12) (with **v** not depending on time), we obtain that $b(t, \mathbf{v}) = 0$ for all $\mathbf{v} \in V_D$, for all $t \in (0, T)$. Thus, using Lemma 3.10, we conclude that, for all $t \in (0, T)$, there exists $P(t) \in L^2(\Omega)$ satisfying

$$b(t, \mathbf{v}) = \langle B^T P(t), \mathbf{v} \rangle_{X'_D, X_D} = -(P(t), \nabla \cdot \mathbf{v})_{\Omega} = - \int_{\Omega} P(t) \nabla \cdot \mathbf{v}, \quad \forall \mathbf{v} \in X_D, \quad (15)$$

which is exactly (13) and ends the second step of the proof.

Third step: the surjectivity of the divergence mapping leads to the following inf-sup condition: there exists $\gamma_2 > 0$, s.t.

$$\inf_{q \in L^2(\Omega)} \sup_{\mathbf{v} \in X_D} \frac{(B\mathbf{v}, q)_{\Omega}}{\|\mathbf{v}\|_{X_D} \|q\|_{L^2(\Omega)}} = \gamma_2 > 0,$$
which implies, for all $q \in L^2(\Omega)$

$$ \gamma_2 \|q\|_{L^2(\Omega)} \le \sup_{\mathbf{v} \in X_D} \frac{(B\mathbf{v}, q)_\Omega}{\|\mathbf{v}\|_{X_D}}. \quad (16) $$

In order to use $q = P(t)$ in (16), we need to evaluate $(B\mathbf{v}, P(t))_\Omega$. From (15), we obtain that $(B\mathbf{v}, P(t))_\Omega = \langle B^T P(t), \mathbf{v} \rangle_{X'_D, X_D} = b(t, \mathbf{v})$; together with (14), we get

$$ \|P(t)\|_{L^2(\Omega)} \le \frac{1}{\gamma_2} \sup_{\mathbf{v} \in X_D} \frac{b(t, \mathbf{v})}{\|\mathbf{v}\|_{X_D}} \le \frac{C_2}{\gamma_2}. $$

We conclude that $P(t) \in L^\infty(0, T; L^2(\Omega))$. Then, we define the pressure $p = \partial_t P$; thus $p \in H^{-1,\infty}(0, T; L^2(\Omega))$, which concludes the third step of the proof.

The final step of the proof is to show that $p$ is unique. Consider the case $\mathbf{u}_0 = 0$ and $c = 0$. Then, we have $\mathbf{u} = 0$, and (13) leads to $\int_{\Omega} P(t) \nabla \cdot \mathbf{v} = 0, \forall \mathbf{v} \in X_D$. From the surjectivity of the divergence mapping, one gets that $P(t) = 0$ for all $t$, and then $p = 0$. □

## 4. OPTIMIZED SCHWARZ WAVEFORM RELAXATION ALGORITHM

The OSWR algorithm for solving the multidomain problem (3)–(4) is as follows.

**Algorithm 1 (OSWR)**

Choose initial Robin data $g_{ij}^0, \xi_{ij}^0$ on $\Gamma_{ij} \times (0, T)$, $j \in \mathcal{I}_i$, $i = 1, 2, \dots, M$ for $\ell = 1, 2, \dots, \mathbf{d}$ do

1. Solve the local space-time Robin problems, for $i = 1, 2, \dots, M$

$$ \begin{align} \partial_t \mathbf{u}_i^\ell - \nu \Delta \mathbf{u}_i^\ell + \nabla p_i^\ell &= \mathbf{f}_i && \text{in } \Omega_i \times (0, T) \notag \\ \nabla \cdot \mathbf{u}_i^\ell &= 0 && \text{in } \Omega_i \times (0, T) \notag \\ \mathbf{u}_i^\ell(\cdot, t=0) &= \mathbf{u}_{0,i} && \text{in } \Omega_i \notag \\ \alpha_{ij}(\nu \partial_{n_{ij}} \mathbf{u}_i^\ell \cdot \mathbf{n}_{ij} - p_i^\ell) + \mathbf{u}_i^\ell \cdot \mathbf{n}_{ij} &= g_{ij}^{\ell-1} && \text{on } \Gamma_{ij} \times (0, T), \quad j \in \mathcal{I}_i \notag \\ \beta_{ij} \nu \partial_{n_{ij}} \mathbf{u}_i^\ell \times \mathbf{n}_{ij} + \mathbf{u}_i^\ell \times \mathbf{n}_{ij} &= \xi_{ij}^{\ell-1} && \text{on } \Gamma_{ij} \times (0, T), \quad j \in \mathcal{I}_i \notag \\ \mathbf{u}_i^\ell &= 0 && \text{on } (\partial \Omega_i \cap \partial \Omega) \times (0, T) \tag{17} \end{align} $$

2. Update the Robin terms $g_{ij}^\ell, \xi_{ij}^\ell$ on $\Gamma_{ij} \times (0, T)$, for $j \in \mathcal{I}_i$, $i = 1, 2, \dots, M$

$$ g_{ij}^{\ell} = \alpha_{ij}(\nu \partial_{n_{ij}} \mathbf{u}_j^{\ell} \cdot \mathbf{n}_{ij} - p_j^{\ell}) + \mathbf{u}_j^{\ell} \cdot \mathbf{n}_{ij}, \quad (18a) $$

$$ \xi_{ij}^{\ell} = \beta_{ij} \nu \partial_{n_{ij}} \mathbf{u}_j^{\ell} \times \mathbf{n}_{ij} + \mathbf{u}_j^{\ell} \times \mathbf{n}_{ij}. \quad (18b) $$

end for

**Remark 1.** Let $i \in [1, M]$, $j \in \mathcal{I}_i$. Formulas given by (18) can be rewritten as

$$ g_{ij}^{\ell} = \frac{\alpha_{ij}}{\alpha_{ji}} (\alpha_{ji}(\nu \partial_{n_{ji}} \mathbf{u}_j^{\ell} \cdot \mathbf{n}_{ji} - p_j^{\ell}) + \mathbf{u}_j^{\ell} \cdot \mathbf{n}_{ji}) - \frac{\alpha_{ij}}{\alpha_{ji}} \mathbf{u}_j^{\ell} \cdot \mathbf{n}_{ji} + \mathbf{u}_j^{\ell} \cdot \mathbf{n}_{ij} $$

$$ \xi_{ij}^{\ell} = \frac{\beta_{ij}}{\beta_{ji}} (\beta_{ji} \nu \partial_{n_{ji}} \mathbf{u}_j^{\ell} \times \mathbf{n}_{ji} + \mathbf{u}_j^{\ell} \times \mathbf{n}_{ji}) - \frac{\beta_{ij}}{\beta_{ji}} \mathbf{u}_j^{\ell} \times \mathbf{n}_{ji} + \mathbf{u}_j^{\ell} \times \mathbf{n}_{ij}, $$

or equivalently, using the Robin transmission conditions in (17),

$$ g_{ij}^{\ell} = \frac{\alpha_{ij}}{\alpha_{ji}} g_{ji}^{\ell-1} - \frac{\alpha_{ij} + \alpha_{ji}}{\alpha_{ji}} \mathbf{u}_j^{\ell} \cdot \mathbf{n}_{ji}, \quad (19a) $$

$$ \xi_{ij}^{\ell} = \frac{\beta_{ij}}{\beta_{ji}} \xi_{ji}^{\ell-1} - \frac{\beta_{ij} + \beta_{ji}}{\beta_{ji}} \mathbf{u}_j^{\ell} \times \mathbf{n}_{ji}. \quad (19b) $$
One advantage of formula (19) is that, if $g_{ij}^{\ell-1}$ and $\xi_{ij}^{\ell-1}$ have $L^2(\Gamma_{ij})$ regularity, so will $g_{ij}^{\ell}$ and $\xi_{ij}^{\ell}$. Indeed, in (19) the regularities of $g_{ij}^{\ell}$ and $\xi_{ij}^{\ell}$ depend only on those of $g_{ji}^{\ell-1}$, $\xi_{ji}^{\ell-1}$ and $\mathbf{u}_j^{\ell}$, whose trace is in $L^2(0, T; H^{\frac{1}{2}}(\Gamma_{ij}))$ (recall that we have $\mathbf{u}_j^{\ell} \in L^2(0, T; [H^1(\Omega_j)]^2)$, see Section 3). On the other hand, formula (18) will return new Robin boundary data $g_{ij}^{\ell}$ and $\xi_{ij}^{\ell}$ with a lower regularity, which is not satisfying for an iterative algorithm. Another advantage of formula (19) is that it is easier to implement in practice, than formula (18) (as one doesn't need to calculate the discrete counterpart of the fluxes).

Now, we may express the iterative algorithm in the following way. We first define

$$
\begin{align*}
V_i &= \{\mathbf{u} \in [H^1(\Omega_i)]^2, \mathbf{u} = 0 \text{ on } \partial\Omega_i \cap \partial\Omega, \nabla \cdot \mathbf{u} = 0 \text{ in } \Omega_i\}, \\
H_i &= \{\mathbf{u} \in [L^2(\Omega_i)]^2, \mathbf{u} \cdot \mathbf{n}_{\partial\Omega_i} = 0 \text{ on } \partial\Omega_i \cap \partial\Omega, \nabla \cdot \mathbf{u} = 0 \text{ in } \Omega_i\}. \\
X_i &= \{\mathbf{u} \in [H^1(\Omega_i)]^2, \mathbf{u} = 0 \text{ on } \partial\Omega_i \cap \partial\Omega\},
\end{align*}
$$

Then, we set, for all $\mathbf{u}, \mathbf{v} \in X_i$ and $t \in (0, T)$,

$$
\begin{align}
a_i(\mathbf{u}, \mathbf{v}) &:= \nu (\nabla \mathbf{u}, \nabla \mathbf{v})_{\Omega_i} + \sum_{j \in \mathcal{I}_i} \frac{1}{\alpha_{ij}} (\mathbf{u} \cdot \mathbf{n}_{ij}, \mathbf{v} \cdot \mathbf{n}_{ij})_{\Gamma_{ij}} + \frac{1}{\beta_{ij}} (\mathbf{u} \times \mathbf{n}_{ij}, \mathbf{v} \times \mathbf{n}_{ij})_{\Gamma_{ij}}, \tag{20} \\
c_i^\ell(t, \mathbf{v}) &:= (\mathbf{f}(t), \mathbf{v})_{\Omega_i} + \sum_{j \in \mathcal{I}_i} \frac{1}{\alpha_{ij}} (g_{ij}^{\ell-1}(t), \mathbf{v} \cdot \mathbf{n}_{ij})_{\Gamma_{ij}} + \frac{1}{\beta_{ij}} (\xi_{ij}^{\ell-1}(t), \mathbf{v} \times \mathbf{n}_{ij})_{\Gamma_{ij}},
\end{align}
$$

and the algorithm reads: for all $\ell \ge 1$, given $g_{ij}^{\ell-1}$, $\xi_{ij}^{\ell-1}$ on each space-time interface $\Gamma_{ij} \times (0, T)$, solve, for each $i = 1...M$:

$$
\begin{align}
\langle \partial_t \mathbf{u}_i^\ell, \mathbf{v} \rangle_{V_i', V_i} + a_i(\mathbf{u}_i^\ell, \mathbf{v}) &= c_i^\ell(t, \mathbf{v}), \quad \text{a.e. } t \in (0, T), \forall \mathbf{v} \in V_i, \tag{21} \\
\mathbf{u}_i^\ell(0) &= \mathbf{u}_{0,i}. \nonumber
\end{align}
$$

Then we construct $p_i^\ell = \partial_t P_i^\ell$, where $P_i^\ell$ is such that

$$
\begin{equation}
\begin{split}
& (\mathbf{u}_i^\ell(t), \mathbf{v})_{\Omega_i} - (\mathbf{u}_{0,i}, \mathbf{v})_{\Omega_i} + \int_0^t a_i(\mathbf{u}_i^\ell(s), \mathbf{v})ds - (P_i^\ell, \nabla \cdot \mathbf{v})_{\Omega_i} - \int_0^t c_i^\ell(s, \mathbf{v})ds = 0, \\
& \qquad \forall \mathbf{v} \in X_i.
\end{split}
\tag{22}
\end{equation}
$$

Finally, the data are updated by using (19a)–(19b) on the space-time interfaces.

With this formulation, we can state the following result

**Theorem 4.1.** Assume that $g_{ij}^0, \xi_{ij}^0 \in L^2(0, T; L^2(\Gamma_{ij}))$ and $\mathbf{u}_0|_{\Omega_i} \in H_i$. Then, the OSWR algorithm is well-defined and for all $\ell$, $\mathbf{u}_i^\ell \in L^2(0, T; V_i) \cap C^0([0, T]; H_i)$, $\partial_t \mathbf{u}_i^\ell \in L^2(0, T; V_i')$, $p_i^\ell \in W^{-1,\infty}(0, T; L^2(\Omega_i))$ and $g_{ij}^\ell, \xi_{ij}^\ell \in L^2(0, T; L^2(\Gamma_{ij}))$.

*Proof.* By Theorem 3.11, if $g_{ij}^{\ell-1}$, $\xi_{ij}^{\ell-1} \in L^2(0, T; L^2(\Gamma_{ij}))$, then one gets $\mathbf{u}_i^\ell$ verifying (21) with $\mathbf{u}_i^\ell \in L^2(0, T; V_i) \cap C^0([0, T]; H_i)$ and $\partial_t \mathbf{u}_i^\ell \in L^2(0, T; V_i')$. Additionally, Theorem 3.11 tells us that there exists $P_i^\ell$ verifying (22). We take $p_i^\ell = \partial_t P_i^\ell \in W^{-1,\infty}(0, T; L^2(\Omega_i))$.

Using the trace theorem, the normal and tangent traces of $\mathbf{u}_i^\ell$ on $\Gamma_{ij} \times (0, T)$ belong to $L^2(0, T; L^2(\Gamma_{ij}))$. Hence, using the update formula (19), we infer that $g_{ij}^\ell, \xi_{ij}^\ell \in L^2(0, T; L^2(\Gamma_{ij}))$. The proof is then carried out by a simple induction. $\square$

**Remark 4.2.** The OSWR algorithm is constructed without considering the last condition (6), hence it may not converge to the monodomain solution. We shall show in the next section that, indeed, the pressure in each subdomain may not converge to the restriction of the monodomain pressure.
## 5. FIRST OBSERVATIONS ON THE TWO SUBDOMAINS CASE

For the trivial case of a one-dimensional problem and two subdomains, one can show that the velocity iterates converge, while the pressure iterates do not converge in general, see [10].

This result generalizes to higher dimensions as follows: let us consider the two-subdomain case, i.e. $M = 2$. To simplify notation, we set $\Gamma := \Gamma_{12} = \Gamma_{21}$, and for any $\phi$ in $(\alpha, g, \mathbf{u})$, we write $\phi_1$ and $\phi_2$ instead of $\phi_{12}$ and $\phi_{21}$, respectively.

The divergence-free condition of the velocity in each subdomain leads to

$$ \int_{\partial\Omega_i} \mathbf{u}_i^\ell \cdot \mathbf{n}_{\partial\Omega_i} = 0 = \int_{\Gamma} \mathbf{u}_i^\ell \cdot \mathbf{n}_i, \quad i = 1, 2. \qquad (23) $$

The update of Robin terms for the normal components can also be written as

$$ g_i^\ell = \frac{\alpha_i}{\alpha_j} g_j^{\ell-1} - \frac{\alpha_i + \alpha_j}{\alpha_j} \mathbf{u}_j^\ell \cdot \mathbf{n}_j, \quad j = 3-i, i = 1, 2. $$

Integrating over $\Gamma$, and taking (23) into account, we get

$$ \int_{\Gamma} g_i^\ell = \frac{\alpha_i}{\alpha_j} \int_{\Gamma} g_j^{\ell-1} = \int_{\Gamma} g_i^{\ell-2}, \quad j = 3-i, i = 1, 2. $$

Therefore, a necessary condition for the convergence of the algorithm to the monodomain solution is

$$ \int_{\Gamma} g_i^0 = \int_{\Gamma} g_i, \quad i = 1, 2, \qquad (24) $$

with $g_i = \alpha_i(\nu \partial_{\mathbf{n}_i} \mathbf{u} \cdot \mathbf{n}_i - p) + \mathbf{u} \cdot \mathbf{n}_i$, $i = 1, 2$, in which $(\mathbf{u}, p)$ is the monodomain solution of problem (1). Condition (24) cannot be achieved in practice because the quantity $g_i$, $i = 1, 2$, is not known.

More precisely, whereas the convergence of the velocity iterates will be proven in Section 6 below, independently of condition (24), the pressure iterates will converge only if condition (24) is satisfied, and thus will not converge in general. A correction technique to recover the pressure from the velocity will be proposed in Section 7.

## 6. CONVERGENCE OF THE VELOCITY via ENERGY ESTIMATE

In this Section, we suppose additional regularity on $\mathbf{u}_0$, $\mathbf{f}$ and $\Omega$, which leads to regularity properties of the strong solution of problem (1)-(2). Namely, we recall Theorem 1, Page 86 from [32].

**Theorem 6.1.** Let $\Omega$ be a bounded domain of $\mathbb{R}^2$ with twice continuously differentiable boundary. For any $\mathbf{u}_0 \in V$ and $\mathbf{f} \in L^2(0, T; L^2(\Omega))^2$, problem (1)-(2) has a unique solution $(\mathbf{u}, p)$ such that

$$ \mathbf{u} \in C^0([0, T]; V) \cap L^2(0, T; (H^2(\Omega))^2), \quad \partial_t \mathbf{u} \in L^2(0, T; L^2(\Omega))^2, \\ p \in L^2(0, T; H^1(\Omega)). $$

Using Theorem 6.1, we prove that, if its hypotheses are satisfied, then the velocity iterates converge to the monodomain velocity.

**Theorem 6.2.** Assume that the hypotheses of Theorem 6.1 are satisfied. Let $g_{ij}^0$ and $\xi_{ij}^0$ belong to $L^2(0, T; L^2(\Gamma_{ij}))$ and let $\mathbf{u}_i^\ell$ be the velocity component of the solution of Algorithm 1 (OSWR). Then, if $\alpha_{ij} = \alpha_{ji}$ and $\beta_{ij} = \beta_{ji}$, the sequence $\mathbf{u}_i^\ell$ converges to $\mathbf{u}_i = \mathbf{u}|_{\Omega_i}$ in $C^0([0, T]; H_i) \cap L^2(0, T; V_i)$.
*Proof.* Denote by $p_i = p|_{\Omega_i}$. Then, thanks to the extra regularity of $(\mathbf{u}, p)$ given by Theorem 6.1, we can define its Robin trace on any space-time interface $\Gamma_{ij} \times (0, T)$

$$g_{ij} = \alpha_{ij}(\nu \partial_{\mathbf{n}_{ij}} \mathbf{u}_j \cdot \mathbf{n}_{ij} - p_j) + \mathbf{u}_j \cdot \mathbf{n}_{ij}, \quad (25a)$$

$$\xi_{ij} = \beta_{ij} \nu \partial_{\mathbf{n}_{ij}} \mathbf{u}_j \times \mathbf{n}_{ij} + \mathbf{u}_j \times \mathbf{n}_{ij}, \quad (25b)$$

and they both belong to $L^2(0, T; L^2(\Gamma_{ij}))$. Then (5) implies

$$g_{ij} = \frac{\alpha_{ij}}{\alpha_{ji}} g_{ji} - \frac{\alpha_{ij} + \alpha_{ji}}{\alpha_{ji}} \mathbf{u}_j \cdot \mathbf{n}_{ji}, \quad (26a)$$

$$\xi_{ij} = \frac{\beta_{ij}}{\beta_{ji}} \xi_{ji} - \frac{\beta_{ij} + \beta_{ji}}{\beta_{ji}} \mathbf{u}_j \times \mathbf{n}_{ji}. \quad (26b)$$

Moreover, $(\mathbf{u}_i, p_i)$ is the strong solution of each local Robin boundary problem with source term $\mathbf{f}_i$, initial condition $\mathbf{u}_{0,i}$ and Robin terms $g_{ij}$ and $\xi_{ij}$ on $\Gamma_{ij}$. We can write these local problems in variational forms similar to (20)–(21), in which we replace $g_{ij}^\ell$ by $g_{ij}$ and $\xi_{ij}^\ell$ by $\xi_{ij}$.

We define the errors as the differences between the iterates and the restrictions (to each subdomain) of the monodomain solution and denote by

$$\mathbf{e}_i^\ell := \mathbf{u}_i^\ell - \mathbf{u}_i, \quad h_{ij}^\ell = g_{ij}^\ell - g_{ji}^\ell, \quad \zeta_{ij}^\ell = \xi_{ij}^\ell - \xi_{ji}, \quad j \in \mathcal{I}_i, i \in [1, M]. \quad (27)$$

Then, the errors also verify the following variational problems similar to (20)–(21): for a.e. $t \in (0, T)$, $\forall \mathbf{v} \in V_i$,

$$\langle \partial_t \mathbf{e}_i^\ell, \mathbf{v} \rangle_{V_i', V_i} + a_i(\mathbf{e}_i^\ell, \mathbf{v}) = \sum_{j \in \mathcal{I}_i} \frac{1}{\alpha_{ij}} (h_{ij}^{\ell-1}, \mathbf{v} \cdot \mathbf{n}_{ij})_{\Gamma_{ij}} + \sum_{j \in \mathcal{I}_i} \frac{1}{\beta_{ij}} (\zeta_{ij}^{\ell-1}, \mathbf{v} \times \mathbf{n}_{ij})_{\Gamma_{ij}}, \quad (28)$$

with initial condition $\mathbf{e}_i^\ell(0) = 0$. All integrals on $\Gamma_{ij}$ are well defined since $g_{ij}$ and $\xi_{ij}$ are both in $L^2(0, T; L^2(\Gamma_{ij}))$, and since we have proved that this is also the case for $g_{ij}^\ell$ and $\xi_{ij}^\ell$ as soon as it is true for $\ell = 0$.

With $\alpha_{ij} = \alpha_{ji}$ and $\beta_{ij} = \beta_{ji}$, the update formulas (19) and (26) for the Robin terms on $\Gamma_{ij} \times (0, T)$ lead to

$$\mathbf{e}_i^\ell \cdot \mathbf{n}_{ij} = \frac{1}{2} (h_{ij}^{\ell-1} - h_{ji}^\ell) , \quad \mathbf{e}_i^\ell \times \mathbf{n}_{ij} = \frac{1}{2} (\zeta_{ij}^{\ell-1} - \zeta_{ji}^\ell). \quad (29)$$

Choosing $\mathbf{e}_i^\ell$ as test function in (28), one gets

$$\langle \partial_t \mathbf{e}_i^\ell, \mathbf{e}_i^\ell \rangle_{V_i', V_i} + \nu(\nabla \mathbf{e}_i^\ell, \nabla \mathbf{e}_i^\ell)_{\Omega_i} \\ + \sum_{j \in \mathcal{I}_i} \frac{1}{\alpha_{ij}} (\mathbf{e}_i^\ell \cdot \mathbf{n}_{ij}, \mathbf{e}_i^\ell \cdot \mathbf{n}_{ij})_{\Gamma_{ij}} + \sum_{j \in \mathcal{I}_i} \frac{1}{\beta_{ij}} (\mathbf{e}_i^\ell \times \mathbf{n}_{ij}, \mathbf{e}_i^\ell \times \mathbf{n}_{ij})_{\Gamma_{ij}} \\ = \sum_{j \in \mathcal{I}_i} \frac{1}{\alpha_{ij}} (h_{ij}^{\ell-1}, \mathbf{e}_i^\ell \cdot \mathbf{n}_{ij})_{\Gamma_{ij}} + \sum_{j \in \mathcal{I}_i} \frac{1}{\beta_{ij}} (\zeta_{ij}^{\ell-1}, \mathbf{e}_i^\ell \times \mathbf{n}_{ij})_{\Gamma_{ij}}. \quad (30)$$

On the boundary $\Gamma_{ij}$, $j \in \mathcal{I}_i$, replacing (29) into (30), one gets

$$\langle \partial_t \mathbf{e}_i^\ell, \mathbf{e}_i^\ell \rangle_{V_i', V_i} + \nu(\nabla \mathbf{e}_i^\ell, \nabla \mathbf{e}_i^\ell)_{\Omega_i} + \frac{1}{4} \sum_{j \in \mathcal{I}_i} \frac{1}{\alpha_{ij}} (h_{ij}^{\ell-1} - h_{ji}^\ell, h_{ij}^{\ell-1} - h_{ji}^\ell)_{\Gamma_{ij}} \\ + \frac{1}{4} \sum_{j \in \mathcal{I}_i} \frac{1}{\beta_{ij}} (\zeta_{ij}^{\ell-1} - \zeta_{ji}^\ell, \zeta_{ij}^{\ell-1} - \zeta_{ji}^\ell)_{\Gamma_{ij}} \\ = \frac{1}{2} \sum_{j \in \mathcal{I}_i} \frac{1}{\alpha_{ij}} (h_{ij}^{\ell-1}, h_{ij}^{\ell-1} - h_{ji}^\ell)_{\Gamma_{ij}} + \frac{1}{2} \sum_{j \in \mathcal{I}_i} \frac{1}{\beta_{ij}} (\zeta_{ij}^{\ell-1}, \zeta_{ij}^{\ell-1} - \zeta_{ji}^\ell)_{\Gamma_{ij}},$$
or equivalently

$$ \begin{aligned} & \langle \partial_t \mathbf{e}_i^\ell, \mathbf{e}_i^\ell \rangle_{V'_i, V_i} + \nu \|\nabla \mathbf{e}_i^\ell\|_{\Omega_i}^2 + \frac{1}{4} \sum_{j \in \mathcal{I}_i} \frac{1}{\alpha_{ij}} \|h_{ji}^\ell\|_{\Gamma_{ij}}^2 + \frac{1}{4} \sum_{j \in \mathcal{I}_i} \frac{1}{\beta_{ij}} \|\zeta_{ji}^\ell\|_{\Gamma_{ij}}^2 \\ & = \frac{1}{4} \sum_{j \in \mathcal{I}_i} \frac{1}{\alpha_{ij}} \|h_{ij}^{\ell-1}\|_{\Gamma_{ij}}^2 + \frac{1}{4} \sum_{j \in \mathcal{I}_i} \frac{1}{\beta_{ij}} \|\zeta_{ij}^{\ell-1}\|_{\Gamma_{ij}}^2, \end{aligned} \quad (31) $$

(recall that notation $|| \cdot ||_D$ corresponds to the $L^2(D)$-norm for any set $D$).

Adapting (12) to $\Omega_i$, integrating (31) on $(0, T)$, and using that $\mathbf{e}_i^\ell(0) = 0$, we get

$$ \begin{aligned} & \|\mathbf{e}_i^\ell(T)\|_{\Omega_i}^2 + 2\nu \int_0^T \|\nabla \mathbf{e}_i^\ell\|_{\Omega_i}^2 + \sum_{j \in \mathcal{I}_i} \frac{1}{2\alpha_{ij}} \int_0^T \|h_{ji}^\ell\|_{\Gamma_{ij}}^2 + \sum_{j \in \mathcal{I}_i} \frac{1}{2\beta_{ij}} \int_0^T \|\zeta_{ji}^\ell\|_{\Gamma_{ij}}^2 \\ & = \sum_{j \in \mathcal{I}_i} \frac{1}{2\alpha_{ij}} \int_0^T \|h_{ij}^{\ell-1}\|_{\Gamma_{ij}}^2 + \sum_{j \in \mathcal{I}_i} \frac{1}{2\beta_{ij}} \int_0^T \|\zeta_{ij}^{\ell-1}\|_{\Gamma_{ij}}^2. \end{aligned} \quad (32) $$

Then, summing with respect to $i$, from 1 to $M$, we get

$$ \sum_{i=1}^{M} \|\mathbf{e}_i^\ell(\cdot, T)\|_{\Omega_i}^2 + 2\nu \sum_{i=1}^{M} \int_0^T \|\nabla \mathbf{e}_i^\ell\|_{\Omega_i}^2 + E_R^\ell = E_R^{\ell-1}, $$

where $E_R^\ell = \sum_{i=1}^M \sum_{j \in \mathcal{I}_i} \frac{1}{2\beta_{ij}} \int_0^T \|\zeta_{ij}^\ell\|_{\Gamma_{ij}}^2 + \sum_{i=1}^M \sum_{j \in \mathcal{I}_i} \frac{1}{2\alpha_{ij}} \int_0^T \|h_{ij}^\ell\|_{\Gamma_{ij}}^2$.

Summing now with respect to $\ell$, from 1 to $L$, we obtain

$$ \sum_{\ell=1}^{L} \sum_{i=1}^{M} \|\mathbf{e}_i^\ell(\cdot, T)\|_{\Omega_i}^2 + 2\nu \sum_{\ell=1}^{L} \sum_{i=1}^{M} \int_0^T \|\nabla \mathbf{e}_i^\ell\|_{\Omega_i}^2 (t) dt + E_R^L = E_R^0. $$

As $E_R^L \ge 0$ for all $L$, the sums $\sum_{\ell=1}^L \sum_{i=1}^M \|\mathbf{e}_i^\ell(\cdot, T)\|_{\Omega_i}^2$ and $\sum_{\ell=1}^L \sum_{i=1}^M \int_0^T \|\nabla \mathbf{e}_i^\ell\|_{\Omega_i}^2$ are bounded; hence $\|\mathbf{e}_i^\ell(T)\|_{\Omega_i}^2$ and $\int_0^T \|\nabla \mathbf{e}_i^\ell\|_{\Omega_i}^2(t)dt$ tend to 0 when $\ell \to \infty$.

In addition, in (32), we can integrate on $(0, t)$ instead of $(0, T)$, and we get for all $t \in (0, T)$

$$ \sum_{\ell=1}^{L} \sum_{i=1}^{M} \|\mathbf{e}_i^\ell(t)\|_{\Omega_i}^2 \le E_R^0. $$

This first leads to the convergence of $\|\mathbf{e}_i^\ell(t)\|_{\Omega_i}$ to 0 for all $t$ and thus to the convergence of $\mathbf{e}_i^\ell$ to 0 in $C^0([0, T]; H_i)$, but also to the fact that, integrating on $(0, T)$, it holds that

$$ \sum_{\ell=1}^{L} \sum_{i=1}^{M} \int_{0}^{T} \|\mathbf{e}_{i}^{\ell}(t)\|_{\Omega_{i}}^{2} dt \leq T E_{R}^{0}. $$

This implies that $\int_0^T \|\mathbf{e}_i^\ell(t)\|_{\Omega_i}^2 dt$ tends to 0 when $\ell \to +\infty$. Then, summing with $\int_0^T \|\nabla \mathbf{e}_i^\ell(t)\|_{\Omega_i}^2 dt$ that also tends to 0, we have that $\int_0^T \|\mathbf{e}_i^\ell(t)\|_{[H^1(\Omega_i)]^2}^2 dt$ tends to 0, or, in other words, that $\mathbf{e}_i^\ell$ tends to 0 in $L^2(0, T; V_i)$, for $i \in [1, M]$. $\square$

Now, we prove a convergence result for the pressure. We set $P(t) = \int_0^t p(s)ds$ and $P_i = P|_{\Omega_i}$ and denote the error by $D_i^\ell(t) = (P_i^\ell - P_i)(t)$, $i \in [1, M]$. Then we can state the following result.
**Corollary 6.3.** Let all hypotheses of Theorem 6.2 be satisfied. Then for all $t \in [0, T]$ it holds that $\|D_i^\ell(t) - \frac{1}{|\Omega_i|} \int_{\Omega_i} D_i^\ell(t)\|_{\Omega_i} \to 0$ when $\ell \to \infty$.

*Proof.* Let $i \in [1, M]$. As $(\mathbf{u}_i, p_i)$ is the strong solution of the Robin problem with boundary conditions $g_{ij}, \xi_{ij}, j \in \mathcal{I}_i$, then $P_i$ verifies a variational formulation similar to (22): $\forall \mathbf{v} \in X_i$ it holds

$$ (\mathbf{u}_i(t), \mathbf{v})_{\Omega_i} - (\mathbf{u}_{0,i}, \mathbf{v})_{\Omega_i} + \int_0^t a_i(\mathbf{u}_i(s), \mathbf{v})ds - (P_i(t), \nabla \cdot \mathbf{v})_{\Omega_i} - \int_0^t c_i(s, \mathbf{v})ds = 0 \quad (33) $$

Then, from (22) and (33), taking the test function $\mathbf{v} \in [H_0^1(\Omega_i)]^2 \subset X_i$, the boundary terms in $c_i^\ell(s, \mathbf{v})$ and $c_i(s, \mathbf{v})$ vanish and then $c_i^\ell(s, \mathbf{v}) - c_i(s, \mathbf{v})$ also vanishes. Then we get

$$ (D_i^\ell(t), \nabla \cdot \mathbf{v})_{\Omega_i} = (\mathbf{e}_i^\ell(t), \mathbf{v})_{\Omega_i} + \int_0^t a_i(\mathbf{e}_i^\ell(s), \mathbf{v})ds, \forall \mathbf{v} \in [H_0^1(\Omega_i)]^2. $$

As $(c, \nabla \cdot \mathbf{v})_{\Omega_i} = 0$ for all constants $c$ and $\mathbf{v} \in [H_0^1(\Omega_i)]^2$, the above formulation implies that $\forall \mathbf{v} \in [H_0^1(\Omega_i)]^2$

$$ (D_i^\ell(t) - \frac{1}{|\Omega_i|} \int_{\Omega_i} D_i^\ell(t), \nabla \cdot \mathbf{v})_{\Omega_i} = (\mathbf{e}_i^\ell(t), \mathbf{v})_{\Omega_i} + \int_0^t a_i(\mathbf{e}_i^\ell(s), \mathbf{v})ds. $$

Since $(D_i^\ell - \frac{1}{|\Omega_i|} \int_{\Omega_i} D_i^\ell) \in L_0^2(\Omega_i) = \{p \in L^2(\Omega_i), \int_{\Omega_i} p = 0\}, i \in [1, M]$, from the inf-sup condition there exists $\gamma_3$ s.t.

$$ \|D_i^\ell - \frac{1}{|\Omega_i|} \int_{\Omega_i} D_i^\ell\|_{\Omega_i} \le \frac{1}{\gamma_3} \sup_{\mathbf{v} \in [H_0^1(\Omega_i)]^2} \frac{|(\mathbf{e}_i^\ell(t), \mathbf{v})_{\Omega_i} + \int_0^t a_i(\mathbf{e}_i^\ell(s), \mathbf{v})ds|}{\|\mathbf{v}\|_{[H_0^1(\Omega_i)]^2}}. $$

We apply again the continuity of $a_i(.,.)$

$$ \left| \int_{0}^{t} a_{i}(\mathbf{e}_{i}^{\ell}(s), \mathbf{v}) ds \right| \leq M_{i} \int_{0}^{t} \| \mathbf{e}_{i}^{\ell}(s) \|_{X_{i}} \| \mathbf{v} \|_{X_{i}} ds \leq M_{i} \| \mathbf{v} \|_{[H_{0}^{1}(\Omega_{i})]^{2}} \sqrt{T} \| \mathbf{e}_{i}^{\ell} \|_{L^{2}(0,T; X_{i})} $$

as well as the Cauchy-Schwarz and Poincaré inequalities on $(\mathbf{e}_i^\ell(t), \mathbf{v})_{\Omega_i}$, we get

$$ \|D_i^\ell - \frac{1}{|\Omega_i|} \int_{\Omega_i} D_i^\ell\|_{\Omega_i} \le \frac{1}{\gamma_3} \left[ C_{P_i} \|\mathbf{e}_i^\ell(t)\|_{\Omega_i} + M_i \sqrt{T} \|\mathbf{e}_i^\ell\|_{L^2(0,T; X_i)} \right] $$

with $C_{P_i}$ the Poincaré constant of $\Omega_i$. From the convergence of the velocity, we get the corollary. $\square$

**Remark 6.4.** Corollary 6.3 tells us that, when $\ell$ grows, the space fluctuations of (the time primitive of) the pressure error converge to 0 in each subdomain. Thus, the pressure error itself behaves for large $\ell$ like a piecewise constant function, with values that possibly depend on the subdomain $\Omega_i$ and iteration count $\ell$. And, indeed, numerical results given in Section 10 show that pressure iterates do not converge to the monodomain solution, unless a correction is applied, which is the object of the next Section.

## 7. RECOVERING THE PRESSURE

Let us introduce the notation $\langle p \rangle_O = \frac{1}{|O|} \int_O p \, dx$ for the mean value of a function on a domain $O$ (whatever the space dimension of $O$).

We set $d_i^\ell := p_i - p_i^\ell, i \in [1, M]$, and recall that $h_{ij}^\ell$ is defined in (27).
**Hypothesis 7.1.** In this section, we suppose that, for a.e $t \in (0, T)$

$$ - \|d_i^\ell - \langle d_i^\ell \rangle_{\Omega_i} \|_{\Omega_i} \to 0 \text{ for all } i \text{ when } \ell \to +\infty $$

$$ - (\langle d_i^\ell \rangle_{\Gamma_{ij}} - \langle d_i^\ell \rangle_{\Omega_i}) \text{ tends to } 0 \text{ for all } j \in \mathcal{I}_i, \text{ for all } i, \text{ when } \ell \to +\infty $$

$$ - (\langle h_{ij}^{\ell-1} \rangle_{\Gamma_{ij}} + \alpha_{ij} \langle d_i^\ell \rangle_{\Gamma_{ij}}) \to 0 \text{ for all } j \in \mathcal{I}_i, \text{ for all } i, \text{ when } \ell \to +\infty $$

**Remark 7.2.** The above hypothesis can be implied from stronger assumptions on the regularity and convergence of the velocity. Indeed, suppose that $(e_i^\ell, d_i^\ell)$ is the strong solution of the following Robin problem

$$ \begin{array}{ll} \partial_t e_i^\ell - \nu \Delta e_i^\ell + \nabla d_i^\ell = 0 & \text{in } \Omega_i \times (0, T) \\ \nabla \cdot e_i^\ell = 0 & \text{in } \Omega_i \times (0, T) \\ e_i^\ell(\cdot, t=0) = 0 & \text{in } \Omega_i \\ e_i^\ell = 0 & \text{on } (\partial\Omega \cap \partial\Omega_i) \times (0, T) \\ \alpha_{ij}(\nu \partial_{n_{ij}} e_i^\ell \cdot n_{ij} - d_i^\ell) + e_i^\ell \cdot n_{ij} = h_{ij}^{\ell-1} & \text{on } \Gamma_{ij} \times (0, T) \\ \beta_{ij} \nu \partial_{n_{ij}} e_i^\ell \times n + e_i^\ell \times n_{ij} = \zeta_{ij}^{\ell-1} & \text{on } \Gamma_{ij} \times (0, T) \end{array} $$

with the following convergence

$$ \|e_i^\ell\|_{L^\infty(0,T;[H^2(\Omega_i)]^2)} \to 0, \quad \|\partial_t e_i^\ell\|_{L^\infty(0,T;[L^2(\Omega_i)]^2)} \to 0. $$

From this, we get, for a.e. $t \in (0, T)$, $\|\nabla d_i^\ell(t)\|_{\Omega_i} \to 0$, which implies the first and second items in Hypothesis 7.1. This also implies the convergence of trace of the velocity: for a.e. $t \in (0, T)$, we have $\|\alpha_{ij}\nu\partial_{n_{ij}}e_i^\ell(t) \cdot n_{ij} + e_i^\ell(t) \cdot n_{ij}\|_{\Gamma_{ij}} \to 0$ that leads to the third item in Hypothesis 7.1.

One can rewrite the three items in Hypothesis 7.1 on the error as follows : when $\ell \to +\infty$, $\forall i \in [1, M]$,

$$ \|(p_i^\ell - p_i) - (\langle p_i^\ell \rangle_{\Omega_i} - \langle p_i \rangle_{\Omega_i})\|_{\Omega_i} \to 0, \quad (34) $$

$$ (\langle p_i^\ell - p_i \rangle_{\Gamma_{ij}}) - (\langle p_i^\ell - p_i \rangle_{\Omega_i}) \to 0, \quad \forall j \in \mathcal{I}_i, \quad (35) $$

$$ [\langle g_{ij}^{\ell-1} \rangle_{\Gamma_{ij}} - \langle g_{ij} \rangle_{\Gamma_{ij}}] + \alpha_{ij} \langle p_i^\ell - p_i \rangle_{\Gamma_{ij}} \to 0, \quad \forall j \in \mathcal{I}_i. \quad (36) $$

Expression (34) shows that $p_i^\ell(t)$ will tend to $p_i(t)$ if and only if the mean-value of $p_i^\ell(t)$ on $\Omega_i$ tends to the mean value of $p_i(t)$. However, no constraint was imposed on the mean-value of $p_i^\ell(t)$ in the algorithm, since, thanks to the Robin boundary conditions, such constraint is not necessary to obtain local well-posed problems at each iteration. In Section 5, we observed cases in which $p_i^\ell$ does not converge to the monodomain solution $p_i$. In this section, we build a modified pressure $\tilde{p}_i^\ell$ such that $\tilde{p}_i^\ell(t)$ tends to $p_i(t)$ in $L^2(\Omega_i)$, $i=1, \dots, M$. Let us denote $X_i(t) := \langle p_i(t) \rangle_{\Omega_i}$, $\forall i \in [1, M]$. Then, using this notation, (34) reads

$$ \|(p_i^\ell(t) - \langle p_i^\ell(t) \rangle_{\Omega_i} + X_i(t)) - p_i(t)\|_{L^2(\Omega_i)} \to 0 \text{ when } \ell \to \infty. \quad (37) $$

From (37), we see that $(p_i^\ell(t) - \langle p_i^\ell(t) \rangle_{\Omega_i} + X_i(t))$ is the right approximation to calculate at each iteration since it tends to $p_i(t)$. However, we do not know how to calculate it because $X_i$ is not known. A similar question was raised in the thesis of Lissoni Theorem IV.3.9 from [34] at the discrete level, within a Schwarz algorithm applied at each time step of a time marching scheme for the numerical approximation of the incompressible Navier-Stokes equations.

We introduce below a new quantity $Y_i^\ell(t)$, fully computable at any given iteration $\ell$, that tends to $X_i(t)$ when $\ell$ tends to infinity, from which we will define the modified pressure $\tilde{p}_i^\ell$.

To ease the presentation, we shall set $|\Gamma_{ij}| = 0$, $\alpha_{ij} = 0$ and $g_{ij}^{\ell-1} = 0$ if $j \notin \mathcal{I}_i$. Moreover, we introduce the constant matrix

$$ A = (a_{ij})_{1 \le i,j \le M}, \quad \text{with} \quad a_{ii} = \sum_{j=1, j \ne i}^{M} |\Gamma_{ij}|\alpha_{ij}, \quad \text{and} \quad a_{ij} = -|\Gamma_{ji}|\alpha_{ji} \quad \text{if } j \ne i $$
together with the constant vector $C = (|\Omega_1|, |\Omega_2|, \dots, |\Omega_M|)$ and the sequence of vectors $(b^\ell)_\ell$, with $b^\ell = (b_1^\ell, b_2^\ell, \dots, b_M^\ell)^t$ defined as

$$b_i^\ell = \sum_{j=1}^{M} |\Gamma_{ij}| [\langle g_{ij}^{\ell-1} \rangle_{\Gamma_{ij}} + \alpha_{ij} \langle p_i^\ell \rangle_{\Omega_i}] - \sum_{j=1}^{M} |\Gamma_{ji}| [\langle g_{ji}^{\ell-1} \rangle_{\Gamma_{ji}} + \alpha_{ji} \langle p_j^\ell \rangle_{\Omega_j}].$$

**Theorem 7.3.** Assume that $\alpha_{ij} = \alpha_{ji}$, $\forall (i, j)$. We have the following properties

(i) For all $\ell$, the following system

$$\begin{aligned}
&AY^\ell = b^\ell, \\
&CY^\ell = 0,
\end{aligned}
\qquad (38)$$

has a unique solution $Y^\ell \in \mathbb{R}^M$.

(ii) Moreover, we have $Y^\ell \to X := (X_1, X_2, \dots, X_M)$ in $\mathbb{R}^M$, and we get for all $t$: $\|\tilde{p}_i^\ell - p_i\|_{L^2(\Omega_i)} \to 0$ when $\ell \to \infty$, with $\tilde{p}_i^\ell(t) := p_i^\ell(t) - \langle p_i^\ell(t) \rangle_{\Omega_i} + Y_i^\ell(t)$.

*Proof of (i).* The proof of Theorem 7.3-(i) relies on two main steps:

(a) Existence of solutions to the system $AY^\ell = b^\ell$,

(b) Existence and uniqueness of a solution to system (38) thanks to the additional constraint $CY^\ell = 0$.

Let us start with (a). Because $\alpha_{ij} = \alpha_{ji}$, it holds that $A$ is symmetric and then existence of at least one solution to the system $AY^\ell = b^\ell$ is equivalent to proving that $b^\ell \in \mathrm{Im}(A) = (\mathrm{Ker}(A))^\perp$. Thus, we start with the determination of $\mathrm{Ker}(A)$.

Let $Y = (Y_1, Y_2, \dots, Y_M)^t \in \mathrm{Ker}(A)$. Then, we have $\sum_{j=1}^M a_{ij}Y_j = 0, \forall i \in [1, M]$. As $\alpha_{ij} = \alpha_{ji}$, we have $a_{ii} = -\sum_{j=1, j \neq i}^M a_{ij}$, which implies

$$0 = \sum_{j=1}^{M} a_{ij} Y_j Y_i = \left( \sum_{j=1, j \neq i}^{M} a_{ij} Y_j Y_i \right) + a_{ii} Y_i^2 = \sum_{j=1, j \neq i}^{M} a_{ij} (Y_j Y_i - Y_i^2).$$

Summing the above expression in $i$, and using that $a_{ij} = a_{ji}$, we obtain

$$\sum_{i=1}^{M} \sum_{j=1, j \neq i}^{M} a_{ij} (Y_j Y_i - Y_i^2) = - \sum_{i<j} a_{ij} (Y_i - Y_j)^2 = 0.$$

As $a_{ij} \le 0$ for all $(i, j)$ with $i \neq j$, and $a_{ij} < 0$ as soon as subdomains $i$ and $j$ are neighbours, this implies that $Y_i = Y_j$ for any pair of neighbouring subdomains $i$ and $j$. Since $\Omega$ is connected, this finally implies that all $Y_i$ are equal i.e. $\mathrm{Ker}(A) = \mathrm{span}(\mathbf{e})$ with $\mathbf{e} = (1, 1, \dots, 1, 1)$. Then, $b^\ell \in (\mathrm{Ker}(A))^\perp$ is equivalent to $b^\ell \cdot \mathbf{e} = \sum_{i=1}^M b_i^\ell = 0$. The latter property is proved in the following way:

$$\sum_{i=1}^{M} b_i^\ell = \sum_{i=1}^{M} \left[ \sum_{j=1}^{M} |\Gamma_{ij}| (\langle g_{ij}^{\ell-1} \rangle_{\Gamma_{ij}} + \alpha_{ij} \langle p_i^\ell \rangle_{\Omega_i}) - \sum_{j=1}^{M} |\Gamma_{ji}| (\langle g_{ji}^{\ell-1} \rangle_{\Gamma_{ji}} + \alpha_{ji} \langle p_j^\ell \rangle_{\Omega_j}) \right].$$

Denoting $\Delta_{ij} := |\Gamma_{ij}| (\langle g_{ij}^{\ell-1} \rangle_{\Gamma_{ij}} + \alpha_{ij} \langle p_i^\ell \rangle_{\Omega_i})$, we obtain

$$\sum_{i=1}^{M} b_i^\ell = \sum_{i=1}^{M} \sum_{j=1}^{M} \Delta_{ij} - \sum_{i=1}^{M} \sum_{j=1}^{M} \Delta_{ji} = 0,$$
which proves that $b^\ell \cdot e = 0$, which concludes the proof of (a).

Let us now turn to (b). From (a), we know that there exists at least a solution to $AY = b$; we let $Y^* = Y^* + \mu e$, with $\mu \in \mathbb{R}$. Existence of a solution to (38) follows from the fact that $Ce = |\Omega| \neq 0$: Choosing $\mu = -\frac{1}{|\Omega|}CY^*$ leads to $CY = CY^* + \mu Ce = 0$ and then $Y$ solves (38). As far as uniqueness is concerned, let $Y_1$ and $Y_2$ be two solutions of (38); since $(Y_1 - Y_2) \in \text{Ker}(A)$, then $(Y_1 - Y_2) = \tau e$, with $\tau \in \mathbb{R}$. Since $\tau|\Omega| = \tau Ce = C(Y_1 - Y_2) = 0$ it follows that $\tau = 0$ and $Y_1 = Y_2$. This ends the proof of Theorem 7.3-(i). $\square$

*Proof of Theorem 7.3-(ii).* It relies on the two main results:

(c) $b^\ell \to AX$ in $\mathbb{R}^M$,

(d) $CX = 0$.

Let us prove (c): from the divergence-free property of $u_i$, we have

$$0 = \int_{\Omega_i} \nabla \cdot u_i = \int_{\partial \Omega_i} u_i \cdot n_{\partial \Omega_i} = \sum_{j \in \mathcal{I}_i} \int_{\Gamma_{ij}} u_i \cdot n_{ij}. \quad (39)$$

Moreover, from the definition of $g_{ij}$ in (25a) and the physical transmission conditions (4), we have

$$|\Gamma_{ij}|\langle g_{ij} \rangle_{\Gamma_{ij}} - |\Gamma_{ji}|\langle g_{ji} \rangle_{\Gamma_{ji}} = \int_{\Gamma_{ij}} (g_{ij} - g_{ji}) = 2 \int_{\Gamma_{ij}} u_i \cdot n_{ij}. \quad (40)$$

Hence, from (39) and (40) we get

$$\sum_{j \in \mathcal{I}_i} |\Gamma_{ij}|\langle g_{ij} \rangle_{\Gamma_{ij}} = \sum_{j \in \mathcal{I}_i} |\Gamma_{ji}|\langle g_{ji} \rangle_{\Gamma_{ji}}. \quad (41)$$

Expression (36) is equivalent to

$$\langle g_{ij}^{\ell-1} \rangle_{\Gamma_{ij}} + \alpha_{ij} \langle p_i^\ell - p_i \rangle_{\Gamma_{ij}} \to \langle g_{ij} \rangle_{\Gamma_{ij}}. \quad (42)$$

From (35), we may replace $\langle p_i^\ell - p_i \rangle_{\Gamma_{ij}}$ by $\langle p_i^\ell - p_i \rangle_{\Omega_i}$ in (42), then multiply by $|\Gamma_{ij}|$ and sum over $j \in \mathcal{I}_i$ for a given $i$ to obtain

$$\sum_{j \in \mathcal{I}_i} |\Gamma_{ij}| [\langle g_{ij}^{\ell-1} \rangle_{\Gamma_{ij}} + \alpha_{ij} \langle p_i^\ell - p_i \rangle_{\Omega_i}] \to \sum_{j \in \mathcal{I}_i} |\Gamma_{ij}| \langle g_{ij} \rangle_{\Gamma_{ij}}. \quad (43)$$

In exactly the same way, we also obtain

$$\sum_{j \in \mathcal{I}_i} |\Gamma_{ji}| [\langle g_{ji}^{\ell-1} \rangle_{\Gamma_{ji}} + \alpha_{ji} \langle p_j^\ell - p_j \rangle_{\Omega_j}] \to \sum_{j \in \mathcal{I}_i} |\Gamma_{ji}| \langle g_{ji} \rangle_{\Gamma_{ji}}. \quad (44)$$

Using (43), (44) and (41), we obtain

$$\sum_{j \in \mathcal{I}_i} |\Gamma_{ij}| [\langle g_{ij}^{\ell-1} \rangle_{\Gamma_{ij}} + \alpha_{ij} \langle p_i^\ell \rangle_{\Omega_i} - \alpha_{ij} \langle p_i \rangle_{\Omega_i}] \\ - \sum_{j \in \mathcal{I}_i} |\Gamma_{ji}| [\langle g_{ji}^{\ell-1} \rangle_{\Gamma_{ji}} + \alpha_{ji} \langle p_j^\ell \rangle_{\Omega_j} - \alpha_{ji} \langle p_j \rangle_{\Omega_j}] \to 0,$$

or equivalently

$$\sum_{j \in \mathcal{I}_i} |\Gamma_{ij}| [\langle g_{ij}^{\ell-1} \rangle_{\Gamma_{ij}} + \alpha_{ij} \langle p_i^\ell \rangle_{\Omega_i}] - \sum_{j \in \mathcal{I}_i} |\Gamma_{ji}| [\langle g_{ji}^{\ell-1} \rangle_{\Gamma_{ji}} + \alpha_{ji} \langle p_j^\ell \rangle_{\Omega_j}] \\ \to \sum_{j \in \mathcal{I}_i} |\Gamma_{ij}| \alpha_{ij} \langle p_i \rangle_{\Omega_i} - \sum_{j \in \mathcal{I}_i} |\Gamma_{ji}| \alpha_{ji} \langle p_j \rangle_{\Omega_j}.$$
This is exactly $b^\ell \to AX$.

Let us now prove (d): We have

$$ \int_{\Omega} p_i = \sum_{i=1}^{M} \int_{\Omega_i} p_i = \sum_{i=1}^{M} |\Omega_i| \langle p_i \rangle_{\Omega_i} = 0, $$

i.e. $CX = 0$.

We now prove Theorem 7.3-(ii): From the solution $Y^\ell$ of (38) given by Theorem 7.3-(i), and from (c) and (d), we have $A(Y^\ell - X) \to 0$ and $C(Y^\ell - X) = 0$. Uniqueness of a solution to $AZ = B$ and $CZ = 0$ as soon as $B$ is in $\mathrm{Im}(A)$ and finite dimension now imply that $(Y^\ell - X) \to 0$ when $\ell \to \infty$. Then, using (37), with a triangle inequality, ends the proof of Theorem 7.3-(ii). $\square$

**Remark 7.4.** In the general case of $M$ subdomains, the calculation of $\tilde{p}_i^\ell$ can be done only once, at the last OSWR iteration. It involves solving the coarse problem (38) when $M > 2$, and is given by an explicit formula when $M = 2$ (see Cor. 7.6), thus the cost of calculating the modified pressure is negligible. In practice, if the stopping criterion of the OSWR method is based on the jumps of the Robin data at the subdomain interfaces, then additional calculations of the pressure correction may be needed at intermediate iterations; the associated additional computational cost is expected to remain small since problems (38) are of size $M$.

**Remark 7.5.** Recovering the correct pressure could also be performed from the fact that $\nabla(p_i^\ell - p_i)$ tends to zero when $\ell \to \infty$. Indeed, for a given $\Omega_i$, choosing first an arbitrary point $\mathbf{x}_i \in \Omega_i$, then one may write

$$ p_i(\mathbf{x}) = p_i(\mathbf{x}_i) + (\mathbf{x} - \mathbf{x}_i) \cdot \int_0^1 \nabla p_i(\mathbf{x}_i + t(\mathbf{x} - \mathbf{x}_i)) dt, \quad \forall \mathbf{x} \in \Omega_i. $$

Then, one could replace $\nabla p_i$ by $\nabla p_i^\ell$ to obtain approximate values of the pressure at each point $\mathbf{x}$. However, this formula holds on a given subdomain $\Omega_i$. In order to relate values of the pressures in $\Omega_i$ to those in a neighboring subdomain $\Omega_j$ through this kind of formula, one needs to choose a point on the boundary $\Gamma_{ij}$ that will serve as the point $\mathbf{x}_j$ in the subdomain $\Omega_j$, and so on. At the discrete level, there are several drawbacks to that: this requires further communications between subdomains, the pressure gradient at the boundaries may not be easy to define (e.g. when the pressure is defined as a piecewise constant field like in the Crouzeix-Raviart finite element), and finally there are many ways to go from one cell to another in the mesh, and, due to round-off errors, this may lead to different evaluations of the pressure at a given cell in particular in very large scale computations.

In the two-subdomain case, we use the same notation as in Section 5. Then the calculation of $\tilde{p}_i^\ell$ can be done by the following explicit formula.

**Corollary 7.6.** Let $M = 2$, $\alpha = \alpha_1 = \alpha_2$, and define, for $i = 1, 2$ and $j = 3 - i$,

$$ \tilde{p}_i^\ell = p_i^\ell + \frac{|\Omega_j|}{|\Omega|} \left[ \frac{1}{\alpha} (\langle g_i^{\ell-1} \rangle_\Gamma - \langle g_j^{\ell-1} \rangle_\Gamma) \right] - \frac{|\Omega_i|}{|\Omega|} \langle p_i^\ell \rangle_{\Omega_i} - \frac{|\Omega_j|}{|\Omega|} \langle p_j^\ell \rangle_{\Omega_j}. $$

Then $\tilde{p}_i^\ell$ tends to $p_i$ in $L^2(\Omega_i)$, when $\ell$ tends to infinity, for $i = 1, 2$.

*Proof.* For $M = 2$ we have

$$ b_1^\ell = -b_2^\ell = |\Gamma| [\langle g_1^{\ell-1} \rangle_\Gamma + \alpha \langle p_1^\ell \rangle_{\Omega_1}] - |\Gamma| [\langle g_2^{\ell-1} \rangle_\Gamma + \alpha \langle p_2^\ell \rangle_{\Omega_2}], $$

$$ \mathbf{A} = \begin{bmatrix} \alpha|\Gamma| & -\alpha|\Gamma| \\ -\alpha|\Gamma| & \alpha|\Gamma| \end{bmatrix}, $$

$$ C = [|\Omega_1| \quad |\Omega_2|]. $$
System (38) for $M = 2$ has a unique solution given by

$$Y_1^\ell = \frac{|\Omega_2|}{|\Omega|} \left[ \frac{1}{\alpha} \left( \langle g_1^{\ell-1} \rangle_\Gamma - \langle g_2^{\ell-1} \rangle_\Gamma \right) + \langle p_1^\ell \rangle_{\Omega_1} - \langle p_2^\ell \rangle_{\Omega_2} \right],$$

$$Y_2^\ell = \frac{|\Omega_1|}{|\Omega|} \left[ \frac{1}{\alpha} \left( \langle g_2^{\ell-1} \rangle_\Gamma - \langle g_1^{\ell-1} \rangle_\Gamma \right) + \langle p_2^\ell \rangle_{\Omega_2} - \langle p_1^\ell \rangle_{\Omega_1} \right].$$

Thus we directly obtain, for $i = 1, 2$, $j = 3 - i$,

$$p_i^\ell - \langle p_i^\ell \rangle_{\Omega_i} + Y_i^\ell = \tilde{p}_i^\ell.$$

Then, using Theorem 7.3 (ii), we get $\tilde{p}_i^\ell \to p_i$ in $L^2(\Omega_i)$, for $i = 1, 2$, which ends the proof of Corollary 7.6. $\square$

## 8. CONVERGENCE FACTOR via FOURIER TRANSFORM

The aim of this section is to find a way to conveniently choose the parameters $(\alpha, \beta)$ that play an important role in the actual rate of convergence in numerical experiments.

Let $\Omega = \mathbb{R}^2$. We consider two subdomains $\Omega_1 = (-\infty, 0) \times \mathbb{R}$ and $\Omega_2 = (0, +\infty) \times \mathbb{R}$, as commonly done for the analysis of OSWR methods. To simplify notation, we set $\Gamma := \Gamma_{12} = \Gamma_{21} = \{x = 0\} \times \mathbb{R}$, and denote $\alpha_{12}$ and $\alpha_{21}$ by $\alpha_1$ and $\alpha_2$, respectively. We denote $\mathbf{u} = (u, v)$ the two components of the velocity and set $\mathbf{f} = (f_x, f_y)$. Recall here the Stokes problem

$$\begin{array}{l} \partial_t u - \nu \Delta u + \partial_x p = f_x \\ \partial_t v - \nu \Delta v + \partial_y p = f_y, \quad \text{in } \Omega \times (0, T) \\ \partial_x u + \partial_y v = 0 \\ u(\cdot, t = 0) = u_0, \quad \text{in } \Omega \\ v(\cdot, t = 0) = v_0 \\ u, v \to 0, \quad \text{when } |(x, y)| \to +\infty. \end{array}$$

We write the algorithm for the errors using the same notation $(u, v, p)$, which means that, by linearity, we set $f_x = f_y = 0$ and $u_0 = v_0 = 0$. To avoid additional notation for the Robin terms, we write the OSWR algorithm as follows: starting with $u_i^0, v_i^0, p_i^0$, at step $\ell \ge 1$ and provided $u_i^{\ell-1}, v_i^{\ell-1}, p_i^{\ell-1}$ we solve

$$\begin{array}{l} \partial_t u_i^\ell - \nu \Delta u_i^\ell + \partial_x p_i^\ell = 0 \\ \partial_t v_i^\ell - \nu \Delta v_i^\ell + \partial_y p_i^\ell = 0, \quad \text{in } \Omega_i \times (0, T) \\ \partial_x u_i^\ell + \partial_y v_i^\ell = 0 \\ u_i^\ell(\cdot, t = 0) = 0, \quad \text{in } \Omega_i \\ v_i^\ell(\cdot, t = 0) = 0 \\ u_i^\ell, v_i^\ell \to 0 \quad \text{when } |(x, y)| \to +\infty \end{array}$$

together with transmission condition on $\Gamma \times (0, T)$, for $i = 1, 2$ and $j = 3 - i$:

$$\begin{array}{l} \alpha_i(\nu \partial_x u_i^\ell - p_i^\ell) + (-1)^{i+1} u_i^\ell = \alpha_i(\nu \partial_x u_j^{\ell-1} - p_j^{\ell-1}) + (-1)^{i+1} u_j^{\ell-1} \\ \nu \beta_i \partial_x v_i^\ell + (-1)^{i+1} v_i^\ell = \nu \beta_i \partial_x v_j^{\ell-1} + (-1)^{i+1} v_j^{\ell-1} \end{array}$$

Let us consider the system in $\Omega_1$, and let $\ell \ge 1$. Taking the Fourier transform in time and in $y$-direction with time frequency $\omega$ and space frequency $k \ne 0$, and, for the sake of simplicity, keeping notation $u, v$ instead of $\hat{u}, \hat{v}$, we get

$$i\omega u_1^\ell - \nu \partial_{xx} u_1^\ell + \nu k^2 u_1^\ell + \partial_x p_1^\ell = 0, \quad (45a)$$

$$i\omega v_1^\ell - \nu \partial_{xx} v_1^\ell + \nu k^2 v_1^\ell + ikp_1^\ell = 0, \quad (45b)$$

$$\partial_x u_1^\ell + ikv_1^\ell = 0. \quad (45c)$$
By differentiating equation (45b) with respect to $x$, multiplying (45a) by $(-ik)$, and summing the resulting equations, and denoting $w_1^\ell := \partial_x v_1^\ell - iku_1^\ell$ the vorticity, we get the vorticity equation

$$i\omega w_1^\ell - \nu \partial_{xx} w_1^\ell + \nu k^2 w_1^\ell = 0. \quad (46)$$

Denote by $\lambda = \sqrt{k^2 + \frac{i\omega}{\nu}}$ with positive real part. As $w_1$ vanishes at $-\infty$, one gets

$$w_1^\ell = E^\ell \exp(\lambda x). \quad (47)$$

Using the definition of $w_1$ and differentiating (45c), we get, for $u_1$

$$\partial_{xx} u_1^\ell - k^2 u_1^\ell = -ik w_1^\ell. \quad (48)$$

The homogeneous equation associated to (48) has characteristic roots $\pm|k|$. As $u_1$ and $v_1$ vanish at $-\infty$, we only retain the root $|k|$. Given the form (47) of the right-hand side of (48), its solution can be written under the form

$$u_1^\ell = A^\ell \exp(|k|x) + B^\ell \exp(\lambda x),$$

with $A^\ell, B^\ell \in \mathbb{C}$. Then, using (45c) and (45b), we get

$$v_1^\ell = A^\ell \frac{i|k|}{k} \exp(|k|x) + B^\ell \frac{i\lambda}{k} \exp(\lambda x),$$

$$p_1^\ell = -A^\ell \frac{i\omega}{|k|} \exp(|k|x).$$

Similarly, for domain $\Omega_2$, there exist $C^\ell, D^\ell \in \mathbb{C}$ such that

$$u_2^\ell = C^\ell \exp(-|k|x) + D^\ell \exp(-\lambda x)$$

$$v_2^\ell = -C^\ell \frac{i|k|}{k} \exp(-|k|x) - D^\ell \frac{i\lambda}{k} \exp(-\lambda x)$$

$$p_2^\ell = C^\ell \frac{i\omega}{|k|} \exp(-|k|x)$$

Replacing the above expressions in the transmission conditions, we obtain

$$\alpha_1(\nu|k|A^\ell + \nu\lambda B^\ell + \frac{i\omega}{|k|}A^\ell) + A^\ell + B^\ell = \alpha_1(-\nu|k|C^{\ell-1} - \nu\lambda D^{\ell-1} - \frac{i\omega}{|k|}C^{\ell-1}) + C^{\ell-1} + D^{\ell-1},$$

$$\nu\beta_1(ikA^\ell + \frac{i\lambda^2}{k}B^\ell) + \frac{i|k|}{k}A^\ell + \frac{i\lambda}{k}B^\ell = \nu\beta_1(ikC^{\ell-1} + \frac{i\lambda^2}{k}D^{\ell-1}) - \frac{i|k|}{k}C^\ell - \frac{i\lambda}{k}D^{\ell-1}$$

and

$$\alpha_2(-\nu|k|C^\ell - \nu\lambda D^\ell - \frac{i\omega}{|k|}C^\ell) - C^\ell - D^\ell = \alpha_2(\nu|k|A^{\ell-1} + \nu\lambda B^{\ell-1} + \frac{i\omega}{|k|}A^{\ell-1}) - A^{\ell-1} - B^{\ell-1},$$

$$\nu\beta_2(ikC^\ell + \frac{i\lambda^2}{k}D^\ell) + \frac{i|k|}{k}C^\ell + \frac{i\lambda}{k}D^\ell = \nu\beta_2(ikA^{\ell-1} + \frac{i\lambda^2}{k}B^{\ell-1}) - \frac{i|k|}{k}A^\ell - \frac{i\lambda}{k}B^{\ell-1}.$$
These transmission conditions can be written in matrix form as follows :

$$ \mathcal{M}(\alpha_1, \beta_1) \begin{pmatrix} A^\ell \\ B^\ell \end{pmatrix} = \mathcal{N}(\alpha_1, \beta_1) \begin{pmatrix} C^{\ell-1} \\ D^{\ell-1} \end{pmatrix} \quad \text{and} \quad \mathcal{M}(\alpha_2, \beta_2) \begin{pmatrix} C^\ell \\ D^\ell \end{pmatrix} = \mathcal{N}(\alpha_2, \beta_2) \begin{pmatrix} A^{\ell-1} \\ B^{\ell-1} \end{pmatrix} $$

where

$$ \mathcal{M}(\alpha, \beta) := \begin{bmatrix} 1 + \frac{\nu\alpha\lambda^2}{|k|} & 1 + \alpha\nu\lambda \\ \nu\beta k + \frac{|k|}{k} & \frac{\nu\beta\lambda^2}{k} + \frac{\lambda}{k} \end{bmatrix}, \quad \mathcal{N}(\alpha, \beta) := \begin{bmatrix} 1 - \frac{\nu\alpha\lambda^2}{|k|} & 1 - \alpha\nu\lambda \\ \nu\beta k - \frac{|k|}{k} & \frac{\nu\beta\lambda^2}{k} - \frac{\lambda}{k} \end{bmatrix}. \qquad (49) $$

This leads to the following recurrent formulation

$$ \begin{pmatrix} A^\ell \\ B^\ell \end{pmatrix} = \mathcal{R}(\alpha_1, \alpha_2, \beta_1, \beta_2) \begin{pmatrix} A^{\ell-2} \\ B^{\ell-2} \end{pmatrix}, \quad \forall \ell \ge 2, \qquad (50) $$

where

$$ \mathcal{R}(\alpha_1, \alpha_2, \beta_1, \beta_2) = \mathcal{M}^{-1}(\alpha_1, \beta_1)\mathcal{N}(\alpha_1, \beta_1)\mathcal{M}^{-1}(\alpha_2, \beta_2)\mathcal{N}(\alpha_2, \beta_2). \qquad (51) $$

In view of (50), the convergence properties of the OSWR algorithm, and in particular its rate, will depend on the spectral radius of the matrix $\mathcal{R}$ defined in (51).

**Remark 8.1.** If one sets $\tilde{\alpha} := \nu\alpha$ and $\tilde{\beta} := \nu\beta$, as well as $\tilde{\omega} := \frac{\omega}{\nu}$, then matrices $\mathcal{M}$ and $\mathcal{N}$ (defined in (49)), depend only on $\tilde{\alpha}$, $\tilde{\beta}$, on $\tilde{\omega}$ and on $k$.

Thus, when $\nu$ varies, the convergence rate remains unchanged if $\tilde{\alpha}$ and $\tilde{\beta}$ are kept constant and if the range in which $\tilde{\omega}$ is considered does not change. As will be seen in Section 9, this is the case if $\nu\Delta t$ and $\nu T$ are kept unchanged. This observation is in line with the fact that the non-dimensional form of the Stokes equations is not modified when $\nu T$ is kept constant for a fixed domain size<fnref n="1" />.

**Remark 8.2.** When $k$ tends to 0, the spectral radius of the matrix $\mathcal{R}$ tends to 1. This is coherent with what was observed in Section 5 and in Remarks 4.2 and 6.4, which led us to the pressure correction described in Section 7.

**Remark 8.3.** When $k$ and $\omega$ tend to $+\infty$, the spectral radius of the matrix $\mathcal{R}$ tends to 1. This implies that analysing the iteration matrix does not help to prove the general convergence (for all frequencies) of the algorithm, and that one always needs the energy estimate technique of Section 6 (for another example, see [11]).

**Remark 8.4.** In practical experiments, all equations are discretized in space and time. As far as space discretization is concerned, the solution of the discrete version of (46) remains close to (47) if the space discretization parameter is small enough with respect to $\sqrt{\frac{\nu}{\omega}}$; since $\omega$ is in practice bounded by $\frac{\pi}{\Delta t}$, we expect that the above Fourier analysis may remain close to practical experiments if the term $\sqrt{\nu\Delta t}$ is large enough compared to the space discretization parameter. This has indeed recently been observed for the heat equation in [2]. As far as time discretization is concerned, the inclusion of its effect in the convergence analysis of OSWR methods is a current topic of research, and is for example addressed in [16] where a $Z$-transform is used and in [2], where a discrete-time analysis of the OSWR method is proposed. This issue is also addressed in Section 9.2.

## 9. OPTIMIZED ROBIN PARAMETERS

One can choose $\alpha_1, \alpha_2, \beta_1, \beta_2$ to minimize the convergence factor of the continuous OSWR algorithm, defined in the above section. Such parameters are called *continuous optimized parameters*. However, for the incompressible Stokes problem, we will see in the numerical experiments of Section 10 that better results can be obtained by minimizing the discrete-time counterpart of this convergence factor. The corresponding parameters are then called *discrete-time optimized parameters*. Both of these optimization procedures are described below.

<footnote n="1">¹Starting from (1) and performing a simple multiplicative change of variables $t = Tt'$, $\mathbf{x} = L\mathbf{x}'$, $\mathbf{u} = U\mathbf{u}'$, $p = \frac{UL}{T}p'$ with non-dimensional prime variables, then these prime variables verify a system similar to (1) in which $\nu$ is replaced by $\frac{\nu T}{L^2}$, which shows that for a fixed $L$, the Stokes equations are unchanged if $\nu T$ is kept constant.</footnote>
## 9.1. Continuous optimized parameters

From Section 8, the convergence factor is $\rho(\mathcal{R}(\alpha_1, \alpha_2, \beta_1, \beta_2, k, \omega))$, where $\mathcal{R}$ is defined in (51), and $\rho(\mathcal{R})$ denotes the spectral radius of $\mathcal{R}$. While we have $\max_{(k, \omega) \in \mathbb{R}^2} \rho(\mathcal{R}(\alpha_1, \alpha_2, \beta_1, \beta_2, k, \omega)) = 1$ (see Rem. 8.3), we can use this convergence factor to calculate Robin parameters for numerical computations, for which the frequencies $k$ and $\omega$ are bounded (by frequencies relevant to the global space-time domain and the ones supported by the numerical grid). Thus, we set

$$ \tilde{\rho}_c(\alpha_1, \alpha_2, \beta_1, \beta_2) := \max_{\substack{\pi_L \le k \le \frac{\pi}{h_\Gamma}, \\ \pi_\Gamma \le \omega \le \frac{\pi}{\Delta t}}} \rho(\mathcal{R}(\alpha_1, \alpha_2, \beta_1, \beta_2, k, \omega)), $$

where $L$ is a characteristic size of the computational domain and $h_\Gamma$ is a measure of the mesh step size on the interface (typically the mean-value of the segment lengths).

Let us consider the one-sided Robin case $\alpha := \alpha_1 = \alpha_2 = \beta_1 = \beta_2$, and set $\rho_c(\alpha) := \tilde{\rho}_c(\alpha, \alpha, \alpha, \alpha)$. Then, the continuous optimized Robin parameter $\alpha_c$ is defined as a solution of the following minimization problem :

$$ \rho_c(\alpha_c) = \min_{\alpha>0} \rho_c(\alpha). \qquad (52) $$

## 9.2. Discrete-time optimized parameters

One can also consider the semi-discrete in time counterpart of the continuous convergence factor to better capture the discrete-time frequencies, i.e. replace in the expression of $\mathcal{R}$ the term $i\omega$ by its discrete counterpart using the implicit Euler scheme, that is we replace $i\omega$ by $\frac{1-e^{-i\omega\Delta t}}{\Delta t}$. Equivalently, we replace in the expression of $\mathcal{R}$ (in (51)) the term $\omega$ by $\bar{\omega} := -i\left(\frac{1-e^{-i\omega\Delta t}}{\Delta t}\right)$, and set $\mathcal{R}_{\Delta t}(\alpha_1, \alpha_2, \beta_1, \beta_2, k, \omega) := \mathcal{R}(\alpha_1, \alpha_2, \beta_1, \beta_2, k, \bar{\omega})$.

Then, as above, we define

$$ \tilde{\rho}(\alpha_1, \alpha_2, \beta_1, \beta_2) := \max_{\substack{\pi_L \le k \le \frac{\pi}{h_\Gamma}, \\ \pi_\Gamma \le \omega \le \frac{\pi}{\Delta t}}} \rho(\mathcal{R}_{\Delta t}(\alpha_1, \alpha_2, \beta_1, \beta_2, k, \omega)). $$

Let us consider the one-sided Robin case $\alpha := \alpha_1 = \alpha_2 = \beta_1 = \beta_2$, and define $\rho(\alpha) := \tilde{\rho}(\alpha, \alpha, \alpha, \alpha)$. Then, the Discrete-time (DT) optimized Robin parameter $\alpha^*$ is defined as a solution of the following minimization problem :

$$ \rho(\alpha^*) = \min_{\alpha>0} \rho(\alpha). \qquad (53) $$

**Remark 9.1.** On could also consider optimized Robin-2p parameters $(\alpha, \beta)$ with $\alpha := \alpha_1 = \alpha_2$, $\beta := \beta_1 = \beta_2$, or 2-sided parameters $(\gamma, \delta)$ with $\gamma := \alpha_1 = \beta_1$, $\delta := \alpha_2 = \beta_2$, that optimize the continuous or discrete-time convergence factors as done in [10]. Given their additional complexity, these more general cases will not be considered here, and are the subject of a subsequent article.

**Remark 9.2.** In practice, $\alpha_c$ and $\alpha^*$ are computed by solving numerically problems (52) and (53) respectively, using the `fminsearch` function of MATLAB [38].

# 10. NUMERICAL RESULTS

In this section, we present numerical experiments that illustrate the performances of the OSWR method of Section 4, with Freefem++ [28]. For the space discretization we use the nonconforming Crouzeix-Raviart Finite Element method in 2D (i.e. piecewise linear elements continuous only at the midpoints of the edges of the mesh for the velocity $\mathbf{u} = (u_x, u_y)$, and piecewise constant $\mathbb{P}_0$ elements for the pressure $p$), and consider the backward Euler method for the time discretization.

In what follows, the term "monodomain solution" will refer to the fully discrete solution obtained on the global mesh without domain decomposition.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Example 1: mesh 1 (left) and mesh 2 (right).</figcaption></figure>

In Sections 10.1 to 10.3, we set $\Omega = ]0, 1[ \times ]0, 1[$, $T = 1$, and consider the Stokes problem (1), where the value of the diffusion coefficient $\nu$ will be specified in each of the examples below. From Remark 9.1, only one-sided Robin parameter $\alpha := \alpha_1 = \alpha_2 = \beta_1 = \beta_2$ will be considered. In particular, we will use the theoretical optimized values $\alpha_c$ and $\alpha^*$ defined in Section 9, see Remark 9.2. Random initial Robin data on the space-time interfaces will be used, unless specified.

In Section 10.1 some results are shown on the convergence of the OSWR algorithm, without and with modification of the pressure as in Section 7. In Section 10.2 we illustrate the influence of the Robin parameter on the convergence of the algorithm, and then in Section 10.4 we present results on a more realistic test case.

## 10.1. Recovering the pressure: a rotating velocity example

The diffusion coefficient is $\nu = 0.1$ and we choose the right-hand side $f$ and the values of the boundary and initial conditions so that the exact solution is given by

$$ 
\begin{aligned}
\mathbf{u}(\mathbf{x}, t) &= (-\cos(\pi y) \sin(\pi x) \cos(2\pi t), \sin(\pi y) \cos(\pi x) \cos(2\pi t)), \\
p(\mathbf{x}, t) &= \cos(t)(x^2 - y^2), \quad \forall \mathbf{x} \in \Omega, \forall t \in (0, T).
\end{aligned}
$$

The domain $\Omega$ is decomposed into nine subdomains as in Figure 1, and two meshes will be considered (as shown on Fig. 1), with mesh sizes $h = 0.0625$ and $h = 0.0312$ respectively. To each mesh, the associated time step is $\Delta t = h$.

We choose $\alpha_1 = \alpha_2 = \beta_1 = \beta_2 = \alpha^*$, where $\alpha^*$ is the DT-Optimized Robin parameter defined in Section 9.1, whose value here is $\alpha^* \approx 3.0832 \times 10^{-1}$ for mesh 1 and $\alpha^* \approx 2.2719 \times 10^{-1}$ for mesh 2.

On Figure 2 we show the evolution of the relative errors, of $p$, $u_x$ and $u_y$, in the $L^\infty(0, T; L^2(\Omega))$-norm, between the OSWR and monodomain solutions, as functions of the number of OSWR iterations, for mesh 1 (on the left) and mesh 2 (on the right). The top figures are with non-modified pressure, and the bottom figures are with the modified pressure $\tilde{p}_i^\ell$, $i = 1, 2$, at each iteration $\ell$ (defined in Sect. 7). We observe that, with the non-modified pressure, the method converges for the velocity but not for the pressure, as expected from the observations of Section 5 and Theorem 6.2. On the other hand, with the modified pressure, we see that the method now converges both for the velocity and the pressure, accordingly to Theorem 7.3.

**Remark 10.1.** Even if we calculate a modified pressure at each iteration, we do not use it in the transmission conditions of Algorithm 1, thus this does not change the velocity convergence, as shown on Figure 2.

**Remark 10.2.** Here and in what follows, the pressure is modified at each iteration to illustrate the convergence of the multidomain solution to the monodomain one. A consequence of Remark 7.4 is that in practice one needs only to modify the pressure at the last OSWR iteration, which makes the cost of the modification negligible.
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Example 1: relative errors for $u_x$, $u_y$ and $p$ (in logarithmic scale) versus iterations with non-modified pressure (top), and modified pressure (bottom), for mesh 1 (left) and mesh 2 (right).</figcaption></figure>

## 10.2. Optimized Robin parameters

The domain $\Omega$ is decomposed into two subdomains as in Figure 3, and we consider the three uniform meshes of Figure 3, with mesh sizes on the interface and associated time steps equal to $h_{\Gamma} = \Delta t = 1/12$, $h_{\Gamma} = \Delta t = 1/24$, and $h_{\Gamma} = \Delta t = 1/48$, respectively. In order to analyze the convergence behavior of the method, we simulate the error equations (i.e. we take homogeneous initial and boundary conditions, and $f = 0$). Thus, the OSWR solution converges to zero.

### 10.2.1. Case with a fixed mesh and different values of $\nu$

We consider mesh 2 (i.e. $h_{\Gamma} = \Delta t = 1/24$). In Figure 4, we plot the evolution of the continuous convergence factor $\rho_c$ (on the left) and of the discrete-time convergence factor $\rho$ (on the right), as functions of the Robin parameter $\alpha$, for different values of $\nu$: $\nu = 1$ (solid line), $\nu = 0.5$ (dashed line), $\nu = 0.1$ (dash-dotted line), $\nu = 0.05$ (dotted line). The theoretical optimized values $\alpha_c$ (blue circle) and $\alpha^*$ (red star), are also shown. We observe that both $\alpha_c$ and $\alpha^*$ increase when $\nu$ decreases. However, the values of $\alpha_c$ and $\alpha^*$ are very different, and when $\nu$ decreases, $\alpha^*$ increases faster than $\alpha_c$, with an associated $\rho(\alpha^*)$ that increases slower than $\rho_c(\alpha_c)$.

In Figure 5, we plot the evolution of the relative errors, of $p$, $u_x$ and $u_y$, in the $L^\infty(0, T; L^2(\Omega))$-norm, in logarithmic scale, after twenty OSWR iterations, as functions of the Robin parameter $\alpha$. We also show the values of the errors obtained with optimized parameter $\alpha = \alpha_c$ (blue circle) and DT-optimized parameter $\alpha = \alpha^*$ (red star). The figures correspond to $\nu = 1$ (top left), $\nu = 0.5$ (top right), $\nu = 0.1$ (bottom left), $\nu = 0.05$ (bottom
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Example 2: mesh 1 (left), mesh 2 (middle), and mesh 3 (right).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Example 2: continuous (left) and discrete-time (right) convergence factors versus $\alpha$, with $\alpha_c$ (blue circle) and $\alpha^*$ (red star), with $h_{\Gamma} = \Delta t = 1/24$; for $\nu = 1$ (solid line), $\nu = 0.5$ (dashed line), $\nu = 0.1$ (dash-dotted line), $\nu = 0.05$ (dotted line).</figcaption></figure>

right). We see that $\alpha^*$ is close to the numerical Robin value giving the smallest error after the same number of iterations, while $\alpha_c$ gives a larger error.

### 10.2.2. Case with $\nu$ fixed and different space-time meshes

Let us take $\nu = 0.1$. In Figure 6, we plot the evolution of the continuous (on the left) and discrete-time (on the right) convergence factors, versus $\alpha$, for different space-time meshes with $h_{\Gamma} = \Delta t = 1/12$ (solid line), $h_{\Gamma} = \Delta t = 1/24$ (dashed line), and $h_{\Gamma} = \Delta t = 1/48$ (dash-dotted line). The theoretical optimized values $\alpha_c$ (blue circle) and $\alpha^*$ (red star) are also shown. We observe that both $\alpha_c$ and $\alpha^*$ decrease when the space-time mesh is refined. However, the values of $\alpha_c$ and $\alpha^*$ are again very different.

In Figure 7, we plot the relative errors, of $p$, $u_x$ and $u_y$, in the $L^\infty(0,T;L^2(\Omega))$-norm, after twenty OSWR iterations, versus Robin parameter $\alpha$, for mesh 1 (top left), mesh 2 (top right), and mesh 3 (bottom). We also show the values of the errors obtained with $\alpha = \alpha_c$ (blue circle) and $\alpha = \alpha^*$ (red star). We observe that $\alpha^*$ is close to the numerical Robin value giving the smallest error after the same number of iterations, while $\alpha_c$ gives a larger error, for all space-time meshes considered.

## 10.3. Asymptotic behavior with respect to the discretization parameters

In this part, we assess the asymptotic performance of the OSWR method with respect to the mesh size $h$ (with the choice $\Delta t = h$), with DT-optimized parameter $\alpha^*$ by studying the evolution of the number of iterations needed to reduce the initial error by a given factor.
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Example 2: Relative errors after 20 iterations, for $u_x$, $u_y$ and $p$ (in logarithmic scale) versus $\alpha$, with their values at $\alpha_c$ (blue circles) and at $\alpha^*$ (red stars), with $h_Γ = \Delta t = 1/24$; for $\nu = 1$ (top left), $\nu = 0.5$ (top right), $\nu = 0.1$ (bottom left), $\nu = 0.05$ (bottom right).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6. Example 2: continuous (left) and discrete-time (right) convergence factors versus $\alpha$, with $\alpha_c$ (blue circle) and $\alpha^*$ (red star), with $\nu = 0.1$; for $h_Γ = \Delta t = 1/12$ (solid line), $h_Γ = \Delta t = 1/24$ (dashed line), $h_Γ = \Delta t = 1/48$ (dash-dotted line).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIGURE 7. Example 2: Relative errors after 20 iterations, for $u_x$, $u_y$ and $p$ (in logarithmic scale) versus $\alpha$, with their values at $\alpha_c$ (blue circles) and at $\alpha^*$ (red stars), with $\nu = 0.1$; for $h_\Gamma = \Delta t = 1/12$ (top left), $h_\Gamma = \Delta t = 1/24$ (top right), $h_\Gamma = \Delta t = 1/48$ (bottom).</figcaption></figure>

We take $\nu = 0.1$ and consider two subdomains with the interface located at $x = \frac{1}{2}$. We consider different meshes with $\Delta t = h = \frac{1}{8}, \frac{1}{16}, \frac{1}{32}, \frac{1}{64}, \frac{1}{128}$, respectively.

In Figure 8 we plot the number $l_*$ of iterations that it takes to reduce the relative $L^\infty(0, T; L^2(\Omega))$-error by a factor $10^{-3}$, for $u$ (left) and $p$ (right), as a function of $h$, on a loglog plot. We observe that $l_* = O(h^{-\frac{1}{2}})$. We notice that this corresponds to the asymptotic performance of the OSWR algorithm with optimized Robin parameter for advection reaction diffusion equations in two dimensions proved in [5].

## 10.4. A more realistic test case

In this example we take $\nu = \frac{1}{\mathcal{R}e}$ with $\mathcal{R}e = 200$, and $T = 5$. The initial condition is $u_0 = 0$. The computational domain is represented on Figure 9. The yellow part, denoted by $\Omega_f$, corresponds to the location where the source term $f$ in the Stokes equations does not vanish; we set :

$$f = \begin{cases} (-2(\sin(\pi t) + \cos(4\pi t)), 0) & \text{in } \Omega_f, \\ (0,0) & \text{in } \Omega \setminus \Omega_f. \end{cases}$$

The domain is decomposed into nine nonoverlapping subdomains (i.e. $M = 9$), with interfaces that are drawn in magenta on Figure 9 (on the left). The subdomains are represented on Figure 9 (on the right), together with their meshes, with a total of 46026 degrees of freedom. The time step is $\Delta t = 0.05$.
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Example 3: Asymptotic behavior (log-log scale): number of iterations $l_{\star}$ needed to reduce the relative $L^{\infty}(0, T; L^2(\Omega))$-error by a factor $10^{-3}$, for $\mathbf{u}$ (left) and $p$ (right), as a function of $h$ (with $\Delta t = h$), with DT-optimized parameter $\alpha^{*}$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9. Example 4: Computational domain (with $\Omega_f$ in yellow) and domain decomposition with interfaces in magenta (left), and subdomains meshes (right).</figcaption></figure>

In Figure 10, we plot the pressure $p$ (on the left) and the velocity field ($u_x, u_y$) (on the right) at final time $t = T$. The fluid is moved from right to left in the central zone, then interacts with the vertical left boundary and recirculates from left to right above and below the central zone.

In Figure 11, we plot the evolution of the continuous convergence factor $\rho_c$ (on the left) and discrete-time convergence factor $\rho$ (on the right), as functions of the Robin parameter $\alpha$. The theoretical optimized values $\alpha_c$ (blue circle) and $\alpha^*$ (red star) are also shown. Their numerical values, obtained following Remark 9.2 with $h_{\Gamma} = 0.05$ and $L = 4.25^2$, are given by $\alpha_c \approx 3.2283 \times 10^{-2}$ and $\alpha^* \approx 6.6063 \times 10^{-1}$, and differ from about a factor 20.

In order to initialize Algorithm 1, we first need to define initial Robin data $g_{ij}^0$ and $\xi_{ij}^0$. Here, we consider and compare two different ways of choosing these initial Robin data.

First, we choose them to be constant in time, equal to the Robin operator applied to an initial state ($u_0, p_0$). Here the initial condition is $u_0 = 0$, and, if we do not want to calculate an estimation of $p(t = 0)$, we may

<footnote n="2">²Note that $L = \max(4.25, 3)$, where 4.25 and 3 are the horizontal and vertical lengths of $\Omega$, respectively.</footnote>
<figure><img src="image_11.png" /><figcaption>FIGURE 10. Example 4: pressure (left) and velocity field (right) at final time $t = 5$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 11. Example 4: continuous (left) and discrete-time (right) convergence factors versus $\alpha$, with corresponding theoretical optimized values $\alpha_c$ (blue circle) and $\alpha^*$ (red star).</figcaption></figure>

simply set $p_0 = 0$. Thus a first possible choice is

$$g_{ij}^0 = \xi_{ij}^0 = 0, \quad \text{for } i = 1, \dots, M \text{ and } j = 3 - i. \tag{54}$$

Secondly, we also consider another choice for $g_{ij}^0$ and $\xi_{ij}^0$, that are non-constant in time, and computed as follows: we solve, for $k = 0, \dots, 7$, and $t_k = \frac{k}{4}$, the problems:

$$\begin{gathered} (\nabla P_k, \nabla q) = (f(t_k), \nabla q), \quad \forall q \in H^1(\Omega) \\ \int_{\Omega} P_k = 0. \end{gathered}$$

Then using the time-periodicity of $f$, we set $P_k = P_{k-8}$, for $k = 8, \dots, 20$. Finally, using linear interpolation between $P_k$ and $P_{k+1}$ on the subintervals $[t_k, t_{k+1}]$, for all $k = 0, \dots, 19$, we obtain a non constant in time function that we denote $p^*$, and we set:

$$g_{ij}^0 = -\alpha_{ij} p^*|_{\Omega_j}, \quad \xi_{ij}^0 = 0, \quad \text{for } i = 1, \dots, M \text{ and } j = 3 - i. \tag{55}$$

In Figure 12, we show the evolution of the relative errors, between the OSWR and monodomain solutions, of $\mathbf{u}$ in $L^\infty(0, T; H)$ and $L^2(0, T; V)$ norms$^3$, and of $p$ in $L^\infty(0, T; L^2(\Omega))$ norm$^3$, as a function of the number

<footnote>$^3$in the sense of broken norms.</footnote>
<figure><img src="image_13.png" /><figcaption>FIGURE 12. Example 4: Relative errors in $L^\infty(0, T; H)$ and $L^2(0, T; V)$ norms for $u$, and in $L^\infty(0, T; L^2(\Omega))$ norm for $p$ (in logarithmic scale) versus iterations, with optimized Robin parameters $\alpha_c$ (blue circle, blue dashed and green diamond curves) and $\alpha^*$ (solid magenta, dashdot red, and dotted black curves), with initial Robin datum given by (54) (left) and by (55) (right).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 13. Example 4: Relative errors for $u$ and $p$, in logarithmic scale, after 25 iterations (left) and 50 iterations (right), versus $\alpha$, with their values at $\alpha_c$ (blue circles) and at $\alpha^*$ (red stars).</figcaption></figure>

of iterations, with optimized Robin parameters $\alpha_c$ (blue circle, blue dashed and green diamond curves) and $\alpha^*$ (solid magenta, dash-dotted red, and dotted black curves). On the left figure, the initial Robin datum is defined by (54) and on the right figure, the initial Robin datum is given by (55). Clearly, using $\alpha^*$ and starting from a better initial guess for the Robin data significantly improves convergence.

When the initial Robin datum is given by (54) (left figure), for the curve of $p$, the convergence rate is faster with $\alpha = \alpha_c$ on the first iterations, and then with $\alpha = \alpha^*$ after iteration 50. When the initial Robin datum is given by (55) (right figure), the curves for $p$ have almost the same speed of convergence on the first iterations, with $\alpha = \alpha^*$ and $\alpha = \alpha_c$, with a faster convergence rate with $\alpha = \alpha^*$ after iteration 60. For the velocity $u$, in both cases of initial Robin datum, we observe that convergence is much faster with $\alpha = \alpha^*$ than with $\alpha = \alpha_c$.
This illustrates the importance of incorporating the actual numerical time marching scheme in the convergence factor expression.

In Figure 13, we plot the evolution of the relative errors, in logarithmic scale, of $\mathbf{u}$ in the $L^\infty(0, T; H)$ and $L^2(0, T; V)$ norms, and of $p$ in the $L^\infty(0, T; L^2(\Omega))$-norm, after twenty-five OSWR iterations (left), and fifty iterations (right), as functions of the Robin parameter $\alpha$. We also show the values of the errors obtained with optimized parameter $\alpha = \alpha_c$ (blue circles) and DT-optimized parameter $\alpha = \alpha^*$ (red stars). We see that for the velocity, $\alpha^*$ is close to the numerical Robin value giving the smallest error after the same number of iterations (25 or 50), while $\alpha_c$ gives a larger error. For the pressure, the values of the error at $\alpha^*$ and $\alpha_c$ at iterations 25 and 50 are quite close. The fact that the optimum values of $\alpha$ for the pressure and the velocities are different is not yet understood.

## ACKNOWLEDGEMENTS

The work of the authors was supported by the ANR project CINE-PARA under grant ANR-15-CE23-0019.

## REFERENCES

[1] S. Ali Hassan, C. Japhet and M. Vohralík, A posteriori stopping criteria for space-time domain decomposition for the heat equation in mixed formulations. *Electron. Trans. Numer. Anal.* **49** (2018) 151–181.

[2] A. Arnoult, C. Japhet and P. Omnes, Discrete-time analysis of optimized Schwarz waveform relaxation with Robin parameters depending on the targeted iteration count. *ESAIM:M2AN* **57** (2023) 2371–2396.

[3] E. Audusse, P. Dreyfuss and B. Merlet, Optimized Schwarz waveform relaxation for the primitive equations of the ocean. *SIAM J. Sci. Comput.* **32** (2010) 2908–2936.

[4] D. Bennequin, M.J. Gander and L. Halpern, A homographic best approximation problem with application to optimized Schwarz waveform relaxation. *Math. Comput.* **78** (2009) 185–223.

[5] D. Bennequin, M. Gander, L. Gouarin and L. Halpern, Optimized schwarz waveform relaxation for advection reaction diffusion equations in two dimensions. *Numer. Math.* **134** (2016) 513–567.

[6] P.-M. Berthe, C. Japhet and P. Omnes, Space-time domain decomposition with finite volumes for porous media applications. In *Domain decomposition methods in science and engineering XXI. Proceedings of the 21st international conference, Inria Rennes Center, France, June 25–29, 2012*. Cham, Springer (2014) 567–575.

[7] E. Blayo, D. Cherel and A. Rousseau, Towards optimized Schwarz methods for the Navier-Stokes equations. *J. Sci. Comput.* **66** (2016) 275–295.

[8] F. Boyer and P. Fabrie, Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models. Springer, New York (2012).

[9] H. Brezis, Analyse fonctionnelle. Théorie et Applications. Mathématiques appliquées pour la maîtrise, Masson (1987).

[10] D.Q. Bui, *New space-time domain decomposition algorithms combined with the Parareal algorithm*, Ph.D. thesis, Thèse de doctorat, Mathématiques appliquées, Université Sorbonne Paris Nord (2021).

[11] D.Q. Bui, C. Japhet, Y. Maday and P. Omnes, Coupling parareal with optimized Schwarz waveform relaxation for parabolic problems. *SIAM J. Numer. Anal.* **60** (2022) 913–939.

[12] T. Chacón Rebollo and E. Chacón Vera, A non-overlapping domain decomposition method for the Stokes equations via a penalty term on the interface. *C. R. Math. Acad. Sci. Paris* **334** (2002) 221–226.

[13] D. Cherel, *Décomposition de domaine pour des systèmes issus des équations de Navier-Stokes*, Ph.D. thesis, Université Grenoble Alpes (2012).

[14] O.A. Ciobanu, *Méthode de décomposition de domaine avec adaptation de maillage en espace-temps pour les équations d'Euler et de Navier-Stokes*, Ph.D. thesis, Université Paris (2014) 13.

[15] O. Ciobanu, L. Halpern, X. Juvigny and J. Ryan, Overlapping domain decomposition applied to the Navier-Stokes equations, edited by T. Dickopf, M.J. Gander, L. Halpern, R. Krause and L.F. Pavarino. In: *Domain Decomposition Methods in Science and Engineering XXII*. Cham, Springer International Publishing (2016) 461–470.

[16] S. Clement, F. Lemarié and E. Blayo, Discrete analysis of Schwarz waveform relaxation for a diffusion reaction problem with discontinuous coefficients. *SMAI J. Comput. Math.* **8** (2022) 99–124.

[17] M. Discacciati, A. Quarteroni and A. Valli, Robin-Robin domain decomposition methods for the Stokes-Darcy coupling. *SIAM J. Numer. Anal.* **45** (2007) 1246–1268.

[18] A. Ern and J.-L. Guermond, Theory and Practice of Finite Elements. Springer, New York (2004).
[19] M.J. Gander and A.M. Stuart, Space-time continuous analysis of waveform relaxation for the heat equation. *SIAM J. Sci. Comput.* **19** (1998) 2014–2031.

[20] M.J. Gander and L. Halpern, Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems. *SIAM J. Numer. Anal.* **45** (2007) 666–697.

[21] M.J. Gander, L. Halpern and F. Nataf, Optimal Schwarz waveform relaxation for the one dimensional wave equation. *SIAM J. Numer. Anal.* **41** (2003) 1643–1681.

[22] P. Gervasio, A. Quarteroni and F. Saleri, Spectral approximation of Navier-Stokes equations. In: *Fundamental directions in mathematical fluid mechanics*. Basel, Birkhäuser (2000) 71–127.

[23] E. Giladi and H.B. Keller, Space-time domain decomposition for parabolic problems. *Numer. Math.* **93** (2002) 279–313.

[24] T. Goudon, S. Krell and G. Lissoni, Non-overlapping Schwarz algorithms for the incompressible Navier-Stokes equations with DDFV discretizations. *ESAIM:M2AN* **55** (2021) 1271–1321.

[25] L. Halpern and J. Széftel, Nonlinear nonoverlapping Schwarz waveform relaxation for semilinear wave propagation. *Math. Comput.* **78** (2009) 865–889.

[26] L. Halpern, C. Japhet and J. Széftel, Optimized Schwarz waveform relaxation and discontinuous Galerkin time stepping for heterogeneous problems. *SIAM J. Numer. Anal.* **50** (2012) 2588–2611.

[27] R.D. Haynes and K. Mohammad, Fully discrete Schwarz waveform relaxation on two bounded overlapping subdomains, edited by R. Haynes, S. MacLachlan, X.-C. Cai, L. Halpern, H.H. Kim, A. Klawonn and O. Widlund. In: *Domain Decomposition Methods in Science and Engineering XXV*. Cham, Springer International Publishing (2020) 159–166.

[28] F. Hecht, New development in freefem++. *J. Numer. Math.* **20** (2012) 251–265.

[29] T.-T.-P. Hoang and H. Lee, A global-in-time domain decomposition method for the coupled nonlinear Stokes and Darcy flows. *J. Sci. Comput.* **87** (2021) 22.

[30] T.-T.-P. Hoang, C. Japhet, M. Kern and J. E. Roberts, Space-time domain decomposition for advection-diffusion problems in mixed formulations. *Math. Comput. Simul.* **137** (2017) 366–389.

[31] C. Japhet and F. Nataf, The best interface conditions for domain decomposition methods: absorbing boundary conditions. In: *Absorbing boundaries and layers, domain decomposition methods, Nova Sci. Publ.* Huntington, NY (2001) 348–373.

[32] O.A. Ladyzhenskaya, The Mathematical Theory of Viscous Incompressible Flow. Gordon and Breach, Sciences Publishers (1963).

[33] F. Lemarié, L. Debreu and E. Blayo, Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients. Part 2: The variable coefficients case. *Electron. Trans. Numer. Anal.* **40** (2013) 170–186.

[34] G. Lissoni, *DDFV method : applications to fluid mechanics and domain decomposition*, Ph.D. thesis, COMUE Université Côte d'Azur (2019).

[35] G. Lube, L. Müller and H. Müller, A new non-overlapping domain decomposition method for stabilized finite element methods applied to the non-stationary Navier-Stokes equations. *Numer. Linear Algebra Appl.* **7** (2000) 449–472.

[36] V. Martin, An optimized Schwarz waveform relaxation method for the unsteady convection diffusion equation in two dimensions. *Appl. Numer. Math.* **52** (2005) 401–428.

[37] V. Martin, Schwarz waveform relaxation algorithms for the linear viscous equatorial shallow water equations. *SIAM J. Sci. Comput.* **31** (2009) 3595–3625.

[38] The Mathworks, Inc., MathWorks Help Center. Natick, Massachusetts, United States. Available at https://fr.mathworks.com/help/matlab/ref/fminsearch.html.

[39] D. Medková, Weak solutions of the Robin problem for the Oseen system. *J. Elliptic Parabol. Equ.* **5** (2019) 189–213.

[40] S. Monniaux and E.M. Ouhabaz, The incompressible Navier-Stokes system with time-dependent Robin-type boundary conditions. *J. Math. Fluid Mech.* **17** (2015) 707–722.

[41] L. Müller and G. Lube, A nonoverlapping DDM for the nonstationary Navier-Stokes problem. *Z. Angew. Math. Mech.* **81** (2001) 725–726.

[42] F.-C. Otto and G. Lube, Non-overlapping domain decomposition applied to incompressible flow problems. In: *Domain Decomposition Methods 10. The 10th International Conference, Boulder, CO, USA, August 10–14, 1997*. Providence, RI, AMS, American Mathematical Society (1998) 507–514.

[43] F.-C. Otto and G. Lube, A nonoverlapping domain decomposition method for the Oseen equations. *Math. Models Methods Appl. Sci.* **8** (1998) 1091–1117.
[44] F.C. Otto, G. Lube and L. Müller, An iterative substructuring method for div-stable finite element approximations of the Oseen problem. *Computing* **67** (2001) 91–117.

[45] L.F. Pavarino and O.B. Widlund, Balancing Neumann-Neumann methods for incompressible Stokes equations. *Commun. Pure Appl. Math.* **55** (2002) 302–335.

[46] R. Russo and A. Tartaglione, On the Robin problem for Stokes and Navier-Stokes systems. *Math. Models Methods Appl. Sci.* **16** (2006) 701–716.

[47] J.C. Strikwerda and C.D. Scarbnick, A domain decomposition method for incompressible viscous flow. *SIAM J. Sci. Comput.* **14** (1993) 49–67.

[48] A. Tartaglione and G. Starita, A note on the Robin problem for the Stokes system. *Rend. Accad. Sci. Fis. Mat. Napoli* **68** (2001) 129–138.

[49] S. Thery, C. Pelletier, F. Lemarié and E. Blayo, Analysis of schwarz waveform relaxation for the coupled ekman boundary layer problem with continuously variable coefficients. *Numer. Algorithms* **89** (2022) 1145–1181.

[50] X. Xu, C.O. Chow and S. H. Lui, On nonoverlapping domain decomposition methods for the incompressible Navier-Stokes equations. *ESAIM:M2AN* **39** (2005) 1251–1269.

<figure><img src="image_15.png" /></figure>

Please help to maintain this journal in open access!

This journal is currently published in open access under the **Subscribe to Open model (S2O).** We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting subscribers@edpsciences.org.

More information, including a list of supporters and financial transparency reports, is available at https://edpsciences.org/en/subscribe-to-open-s2o.