# Schwarz Waveform Relaxation Physics-Informed Neural Networks for Solving
Advection-Diffusion-Reaction Equations

Emmanuel LORIN<EMAIL_ADDRESS>Xu YANG<EMAIL_ADDRESS>School
of Mathematics and Statistics, Carleton University, Ottawa, Canada, K1S 5B6
Centre de Recherches Mathématiques, Université de Montréal, Montréal, Canada,
H3T 1J4  Department of Mathematics, University of California, Santa Barbara,
CA 93106, USA

###### Abstract

This paper develops a physics-informed neural network (PINN) based on the
Schwarz waveform relaxation (SWR) method for solving local and nonlocal
advection-diffusion-reaction equations. Specifically, we derive the
formulation by constructing subdomain-dependent local solutions by minimizing
local loss functions, allowing the decomposition of the training process in
different domains in an embarrassingly parallel procedure. Provided the
convergence of PINN, the overall proposed algorithm is convergent. By
constructing local solutions, one can, in particular, adapt the depth of the
deep neural networks, depending on the solution’s spectral space and time
complexity in each subdomain. We present some numerical experiments based on
classical and Robin-SWR to illustrate the performance and comment on the
convergence of the proposed method.

###### keywords:

Physics-informed neural network, Schwarz waveform relaxation, domain
decomposition, advection-diffusion-reaction equations.

##  1 Introduction

This paper focuses on the derivation and analysis of a Neural-Network (NN)
based Schwarz Waveform Relaxation (SWR) Domain Decomposition Method (DDM) for
solving partial differential equations (PDE) in parallel. We will focus in
this paper on a simple diffusion-advection-reaction equation. Still, the
proposed strategy applies to any other evolution (in particular wave-like)
equations, for which the convergence of SWR is proven [1, 2, 3, 4, 5, 6, 7, 8,
9, 10, 11]. We derive a combined DDM-SWR and Physics-Informed Neural Network
(PINN) method for solving local and nonlocal diffusion-advection-reaction
equations. The latter was developed by Karniadakis et al. [12, 13, 14] and is
a general strategy in scientific machine learning for solving PDE using deep
neural networks via the minimization of well-designed loss functions. Notice
that in [15], was also proposed a more direct DDM for solving PDE.
Interestingly, both ([15] and the one presented here) methods could actually
be combined; this was however not tested in this paper. Let us also mention a
recent paper [16], where a combination of Schwarz DDM with NN-based solvers is
proposed for stationary PDE. Beyond the derivation of the SWR-NN method, this
paper’s objective is to exhibit some fundamental properties that make this
methodology very promising. The general principle is to solve Initial Boundary
Value Problems (IBVPs) by constructing local solutions (subdomain-dependent)
obtained by minimizing local loss functions. The overall strategy is
convergent (provided that the PINN method is convergent) and allows, in
particular, to locally decompose the training process in different subdomains
within an embarrassingly parallel procedure. The construction of local
solutions also allows to locally adapt the depth of the deep neural network,
depending on the solution’s spectral space and time complexity in each
subdomain.

In this paper, we will primarily focus on the derivation aspects and will not
necessarily detail all the computational aspects, particularly regarding the
selection of the training points. This will, however, be specified in the
Numerical Experiment section. For convenience, we shall recall some basic
aspects about PINNs, neural networks and SWR method for evolution equations,
which shall be used in the paper later.

###  1.1 Basics on PINNs

Let us recall the principle of PINNs for solving, e.g., an evolution PDE over
$\Omega\times[0,T]$,

$\displaystyle\left\\{\begin{array}[]{lcl}\partial_{t}u+Pu&=&f\,\,\hbox{ in
}\Omega\times[0,T]\\\ Mu&=&0\,\,\hbox{ in }\Gamma\times[0,T]\\\
u(\cdot,0)&=&u_{0}\,\,\hbox{ in }\Omega\end{array}\right.\,,$ (4)

where i) $P$ is a differential operator in space, and $M$ a differential or
algebraic boundary operator over the domain boundary $\Gamma$; ii) are $f$ and
$u_{0}$ are imposed functions. The PINN approach, which generalizes the DE-
solver from Lagaris [17] consists in parameterizing (by say
${\boldsymbol{w}}$) a NN, $N({\boldsymbol{w}},{\boldsymbol{x}},t)$
approximating the solution to (4), by minimizing (a discrete version of) the
following loss function

$\displaystyle\mathcal{L}({\boldsymbol{w}})$ $\displaystyle=$
$\displaystyle\|(\partial_{t}+P)N({\boldsymbol{w}},\cdot,\cdot)-f\|_{L^{2}(\Omega\times(0,T))}+\lambda_{1}\|MN({\boldsymbol{w}},\cdot,\cdot)\|_{L^{2}(\Gamma\times(0,T))}$
$\displaystyle+\lambda_{2}\|N({\boldsymbol{w}},\cdot,0)-u_{0}\|_{L^{2}(\Omega)}\,,$

where $\lambda_{1,2}$ are some free positive parameters and
${\boldsymbol{w}}\in W\in\mathbb{R}^{P}$ for some large $P$, and where
$\|\cdot\|_{L^{2}(\Omega\times(0,T))}$ (resp.
$\|\cdot\|_{L^{2}(\Gamma\times(0,T))}$) denotes the $L^{2}-$norm over
$\Omega\times(0,T)$ (resp. $\Gamma\times(0,T)$). Practically, the loss
functions are constructed by estimating the values at a very large number of
space$\&$time-training points $\\{({\boldsymbol{x}}_{j},t_{n})\\}_{j,n}$.
Hence the $L^{2}$-norm are not exactly computed, but only approximated.
Karniadakis and collaborators have developed numerous techniques to improve
the efficiency and rate of convergence of PINN-algorithms for different types
of PDE. We refer for instance to [12, 13, 14] for details. Ultimately, the
PINN strategy is to provide more efficient solvers than standard methods
(finite difference, finite-volume, finite-elements, spectral methods,
pseudospectral methods, etc) for high dimensional (stochastic or
deterministic) PDEs. As far as we know, this is not clearly yet established,
which justifies the active research in this field and the developed of new
methods.

###  1.2 Basics of Neural Networks

We here recall the basics of neural networks. We denote the neural network
$N({\boldsymbol{w}},{\boldsymbol{x}})$ with
${\boldsymbol{x}}=(x_{1},\cdots,x_{d})\in\Omega\subseteq\mathbb{R}^{d}$ and
where we denote by ${\boldsymbol{w}}$ the unknown parameters. Neural networks
usually read (for 1 hidden layer, machine learning)

$\displaystyle\left.\begin{array}[]{lcl}N({\boldsymbol{w}},{\boldsymbol{x}})&=&\sum_{i=1}^{H}v_{i}\sigma_{i}\big{(}\sum_{j=1}^{d}w_{ij}x_{j}+u_{i}\big{)}\,,\end{array}\right.$
(6)

where $\\{\sigma_{i}\\}_{1\leqslant i\leqslant H}$ are the sigmoid transfer
functions, and $H$ is the number of sigmoid units, $\\{w_{ij}\\}_{ij}$ are the
weights and $\\{u_{i}\\}_{i}$ the bias. When considering several hidden layers
(deep learning), we have then to compose functions of the form (6), [18]. That
is

$\displaystyle N$ $\displaystyle=$
$\displaystyle\mathcal{N}_{p}\circ\mathcal{N}_{p-1}\circ\cdots\mathcal{N}_{1}\circ\mathcal{N}_{0}\,,$

where for $0\leqslant r\leqslant p$, $\mathcal{N}_{r}$ is defined from
$\mathbb{R}^{a_{r}}$ (with $a_{r}\in\mathbb{N}$) to $\mathbb{R}^{a_{r+1}}$ by
$\sigma_{r}(W_{r}X_{r}+b_{r})$, $\sigma_{r}$ is an activation function,
$X_{r}\in\mathbb{R}^{a_{r}}$ and where $(a_{0},\cdots,a_{p+1})$ where $p+1$
layers are considered. The layer $r=0$ is the input layer and $r=p+1$ is the
output layer, such that $a_{0}=a_{p+1}=m$. In fine, $N$ from
$X\in\mathbb{R}^{m}$ to $\mathbb{R}^{m}$.

###  1.3 Basics on SWR methods for evolution equations

In this subsection, we recall the principle of SWR-DDM for solving evolution
PDE. Consider a $d$-dimensional first order in time evolution partial
differential equation $\partial_{t}u+Pu=f$ in the spatial domain
$\Omega\subseteq\mathbb{R}^{d}$, and time domain $(0,T)$, where $P$ is a
linear differential operator in space. The initial data is denoted by $u_{0}$,
and we impose, say, null Dirichlet boundary conditions on
$\Gamma_{\textrm{ext}}=\partial\Omega$. We present the method for 2
subdomains, although in practice an arbitrary number of subdomains can be
employed. We first split $\Omega$ into two open subdomains
$\Omega_{\varepsilon}^{\pm}$, with or without overlap
($\Omega_{\varepsilon}^{+}\cap\Omega_{\varepsilon}^{-}=\emptyset$ or
$\Omega_{\varepsilon}^{+}\cap\Omega_{\varepsilon}^{-}\neq\emptyset$), and
$\varepsilon\geqslant 0$. The SWR algorithm consists in iteratively solving
IBVPs in $\Omega_{\varepsilon}^{\pm}\times(0,T)$, using transmission
conditions at the subdomain interfaces
$\Gamma_{\varepsilon}^{\pm}:=\partial\Omega_{\varepsilon}^{\pm}\cap\Omega_{\varepsilon}^{\mp}$.
The imposed transmission conditions are established using the preceding
Schwarz iteration data in the adjacent subdomain. That is, for $k\geqslant 1$,
and denoting $u^{\pm}$ the solution in $\Omega_{\varepsilon}^{\pm}$, we
consider

$\displaystyle\left\\{\begin{array}[]{lcl}(\partial_{t}+P)u^{\pm,(k)}&=&f,\,\hbox{
in }\Omega_{\varepsilon}^{\pm}\times(0,T),\\\
u^{\pm,(k)}(\cdot,0)&=&u_{0}^{\pm},\,\hbox{ in }\Omega_{\varepsilon}^{\pm},\\\
\mathcal{T}_{\pm}u^{\pm,(k)}&=&\mathcal{T}_{\pm}u^{\mp,(k-1)},\,\hbox{ on
}\Gamma_{\varepsilon}^{\pm}\times(0,T),\\\ u^{\pm,(k)}&=&0,\,\hbox{ on
}\Lambda_{\varepsilon}^{\pm}\times(0,T)\,,\end{array}\right.$ (11)

with a given initial guess $u^{\pm,(0)}$, where $\mathcal{T}^{\pm}$ denotes a
boundary$/$transmission operator and where
$\Lambda_{\varepsilon}^{\pm}:=\partial\Omega_{\varepsilon}^{\pm}\backslash\Gamma_{\varepsilon}^{\pm}$
are internal boundaries. Classical SWR (CSWR) method consists in taking
$\mathcal{T}^{\pm}$ as the identity operator while Optimized-SWR (OSWR) method
consists in taking
$\mathcal{T}^{\pm}=\nabla_{{\boldsymbol{n}}^{\pm}}\pm\lambda^{\pm}_{\Gamma_{\varepsilon}}$
for some well chosen (optimized from the convergence rate point of view)
$\lambda^{\pm}_{\Gamma^{\pm}_{\varepsilon}}\in\mathbb{R}^{*}$, and outward
normal vector ${\boldsymbol{n}}^{\pm}$ to $\Gamma^{\pm}_{\varepsilon}$. The
OSWR method is then a special case of Robin-SWR methods. In addition, in order
to provide a faster convergence than CSWR, the OSWR method is often convergent
even for non-overlapping DDM. The latter is hence of crucial interest from the
computational complexity point. We refer to [6, 9, 10, 11, 7] for details. The
convergence criterion for the Schwarz DDM is typically given for any
$0<t\leqslant T$, by

$\displaystyle\big{\|}\hskip
5.69046pt\|u^{+,(k)}(\cdot,t)-u^{-,(k)}(\cdot,t)\|_{\infty;\overline{\Omega}_{\varepsilon}^{+}\cap\overline{\Omega}_{\varepsilon}^{-}}\big{\|}_{L^{2}(0,t)}\leqslant\delta^{\textrm{Sc}},$
(12)

with $\delta^{\textrm{Sc}}$ small enough. When the convergence of the full
iterative algorithm is obtained at Schwarz iteration $k^{\textrm{cvg}}$, one
gets the converged global solution $u^{\textrm{cvg}}:=u^{(k^{\textrm{cvg}})}$
in $\Omega$. The reconstructed solution $u$, is finally defined as
$u_{|\Omega_{\varepsilon}^{\pm}}=u^{\pm,(k^{\textrm{cvg}})}$.

###  1.4 Advection-diffusion-reaction equation

Rather than considering a general situation, for which the rapid convergence
of the SWR method and efficiency are not necessarily proven, we propose to
focus on the advection-diffusion-reaction equation, for which both properties
are established in [6] (see also [5, 1, 2, 3, 4] for the Schrödinger
equation). Let us consider the following initial boundary-value problem: find
the real function $u({\boldsymbol{x}},t)$ solution to the advection-diffusion-
reaction equation on $\mathbb{R}^{d}$, $d\geqslant 1$,

$\displaystyle\left\\{\begin{array}[]{l}\partial_{t}u=\nu({\boldsymbol{x}})\triangle
u+{\boldsymbol{a}}({\boldsymbol{x}})\cdot\nabla
u+r({\boldsymbol{x}})u,\,{\boldsymbol{x}}\in\mathbb{R}^{d},\,t>0,\\\
u({\boldsymbol{x}},0)=u_{0}({\boldsymbol{x}}),\,{\boldsymbol{x}}\in\mathbb{R}^{d},\end{array}\right.$
(15)

with initial condition $u_{0}$, and the real-valued space-dependent smooth
reaction term $r$, advection vector ${\boldsymbol{a}}$ and diffusion $\nu$.
We recall from [6], that for considering $\Omega=\mathbb{R}$ with constant
coefficients in (15), and for $u_{0}\in L^{2}(\Omega)$, $f\in
L^{2}(0,T;L^{2}(\Omega))$, there exists a unique weak solution in
$C(0,T;L^{2}(\Omega))\cap L^{2}(0,T;H^{1}(\Omega))$. Moreover, if $u_{0}\in
H^{2}(\Omega)$ and $f\in L^{2}(0,T;L^{2}(\Omega))$, there exists a unique weak
solution in $L^{2}(0,T;H^{3}(\Omega))\cap H^{3/2}(0,T;L^{2}(\Omega))$.

###  1.5 Organization of the paper

The rest of the paper is organized as follows. In Section 2, we derive the
combined SWR-PINN method, and some properties are proposed in Sections 2.2 and
2.4. Section 3 is devoted to some numerical experiments illustrating the
convergence of the overall SWR-PINN method. We make conclusive remarks in
Section 4.

##  2 SWR-PINN method

In this section, we propose to combine PINN-based solvers with SWR-DDM to
solve the advection-diffusion-reaction equation on a bounded domain
$\Omega\subset\mathbb{R}^{d}$, imposing null Dirichlet boundary conditions at
$\Gamma$. For the sake of simplicity of the presentation, the derivation is
proposed for two subdomains; the extension to an arbitrary number of
subdomains is straightforward.

###  2.1 Derivation of the SWR-PINN method

The standard SWR method for two subdomains consists in solving the IVBP using
the following algorithm

$\displaystyle\left\\{\begin{array}[]{lcl}\partial_{t}u^{\pm,(k)}&=&\nu({\boldsymbol{x}})\triangle
u^{\pm,(k)}+{\boldsymbol{a}}({\boldsymbol{x}})\cdot\nabla
u^{\pm,(k)}+r({\boldsymbol{x}})u^{\pm,(k)},\,\hbox{ in
}\Omega_{\varepsilon}^{\pm}\times(0,T),\\\
u^{\pm,(k)}(\cdot,0)&=&u_{0}^{\pm},\,\hbox{ in }\Omega_{\varepsilon}^{\pm},\\\
\mathcal{T}_{\pm}u^{\pm,(k)}&=&\mathcal{T}_{\pm}u^{\mp,(k-1)},\,\hbox{ on
}\Gamma_{\varepsilon}^{\pm}\times(0,T),\\\ u^{\pm,(k)}&=&0,\,\hbox{ on
}\Lambda_{\varepsilon}^{\pm}\times(0,T)\,,\end{array}\right.$ (20)

where $\mathcal{T}_{\pm}$ is a boundary operator, and where we recall that
$\Lambda_{\varepsilon}^{\pm}=\partial\Omega_{\varepsilon}^{\pm}\backslash\Gamma_{\varepsilon}^{\pm}$.
The well-posedness and the convergence of this method and its rate of
convergence were established in [6] for different types of transmission
conditions. SWR algorithms can actually be reformulated as a fixed point
methods (FPM), and their rate of convergence is hence determined by a
contraction factor of the FPM. More specifically, it is proven in [5] that for

$\displaystyle\partial_{t}u+a\partial_{x}u-\nu\partial_{xx}u+ru=f\,\,\,(x,t)\in\mathbb{R}^{2}\times(0,T),$

where $f\in L^{2}(\Omega)$, $\nu>0$, $a,b$ in $\mathbb{R}$, the CSWR method is
convergent and has a convergence rate $C_{\textrm{CSWR}}$ (contract factor),
at least, given by

$\displaystyle C_{\textrm{CSWR}}$ $\displaystyle=$
$\displaystyle\exp\Big{(}-\cfrac{\varepsilon}{\nu}(\sqrt{a^{2}+4\nu
r})\Big{)}\,.$

In fact, this can be refined to superlinear convergence rate
$2/\sqrt{\pi}\int^{\infty}_{\varepsilon/\sqrt{\nu T}}e^{-s^{2}}ds$. For Robin-
SWR methods, with transmission conditions $\partial_{x}\pm\lambda$, the
convergence rate is actually improved

$\displaystyle C_{\textrm{Robin}}$ $\displaystyle=$
$\displaystyle\sup_{\omega\in\mathbb{R}}\Big{|}\cfrac{\sqrt{a^{2}+4\nu(r+{\tt
i}\omega)}-\lambda}{\sqrt{a^{2}+4\nu(r+{\tt
i}\omega)}+\lambda}\Big{|}\exp\Big{(}-\cfrac{\varepsilon}{\nu}(\sqrt{a^{2}+4\nu
r})\Big{)}\,,$

as $|\sqrt{a^{2}+4\nu(r+{\tt i}\omega)}-\lambda|/|\sqrt{a^{2}+4\nu(r+{\tt
i}\omega)}+\lambda|<1$. We notice in particular, that a crucial element for
the rate of convergence of SWR methods, is the size of the overlapping zone.
However, overlapping is not required for Robin-SWR methods to converge.
Rather than a standard approximation of (20) using a finite
elements$/$difference or pseudospectral methods [4, 5, 3], we then propose to
solve this system using a PINN method. We denote by
$N({\boldsymbol{w}},{\boldsymbol{x}},t)$ the generic NN to optimize, where
${\boldsymbol{w}}$ denotes the unknown parameters. The SWR-NN hence consists
in searching for an approximate solution to the SWR method by applying local
PINN algorithms. That is, we now consider

$\displaystyle\left\\{\begin{array}[]{lcl}\partial_{t}N^{\pm,(k)}&=&\nu({\boldsymbol{x}})\triangle
N^{\pm,(k)}+{\boldsymbol{a}}({\boldsymbol{x}})\cdot\nabla
N^{\pm,(k)}+r({\boldsymbol{x}})N^{\pm,(k)},\,\hbox{ in
}\Omega_{\varepsilon}^{\pm}\times(0,T),\\\
N^{\pm,(k)}(\cdot,0)&=&u_{0}^{\pm},\,\hbox{ in }\Omega_{\varepsilon}^{\pm},\\\
\mathcal{T}_{\pm}N^{\pm,(k)}&=&\mathcal{T}_{\pm}N^{\mp,(k-1)},\,\hbox{ on
}\Gamma_{\varepsilon}^{\pm}\times(0,T),\\\ N^{\pm,(k)}&=&0,\,\hbox{ on
}\Lambda_{\varepsilon}^{\pm}\times(0,T)\,.\end{array}\right.$ (25)

###### Remark 2.1

For the CSWR method ($\mathcal{T}_{\pm}$ is the identity operator), it is
proven in [6], among many other well-posed results, that for $f\in
L^{2}(0,T;H^{1}(\Omega_{\varepsilon}^{\pm}))$ and $u_{0}\in H^{1}(\Omega)$
with $u^{\pm,(0)}\in H^{3/4}(0,T)$ and some compatibility conditions, the
algorithms (20) is well-posed in $L^{2}(0,T;H^{2}(\Omega)\cap
H^{1}(0,T;L^{2}(\Omega))$.
Let us now denote $e^{\pm,(k)}:=u^{\pm,(k)}-u$ and
$\widetilde{e}^{\pm,(k)}:=N^{\pm,(k)}-u$. In Theorem 3.3 from [6] it is stated
that for any $k\geqslant 0$ and for some positive constant $C>0$

$\displaystyle\|(e^{+,(2k+1)},e^{-,(2k+1)})\|_{\mathcal{H}_{\varepsilon}}$
$\displaystyle\leqslant$ $\displaystyle
C\exp\Big{(}-k\cfrac{\varepsilon(\sqrt{a^{2}+4\nu
r})}{\nu}\Big{)}\|(u(\varepsilon,\cdot),u(-\varepsilon/2,\cdot))\|_{(_{0}H^{3/4}(0,T))^{2}}\,,$

where we have denoted

$\displaystyle\mathcal{H}_{\varepsilon}$ $\displaystyle:=$ $\displaystyle
L^{2}(0,T;H^{2}(\Omega_{\varepsilon}^{+})\cap
H^{1}(0,T;L^{2}(\Omega_{\varepsilon}^{+}))\times
L^{2}(0,T;H^{2}(\Omega_{\varepsilon}^{+})\cap
H^{1}(0,T;L^{2}(\Omega_{\varepsilon}^{-}))\,.$

Now if we assume (convergence of the PINN-method) that the NN solution to (25)
is such that there exist
$\eta(\overline{\boldsymbol{w}}^{\pm;(k)};\varepsilon)$

$\displaystyle\|N^{\pm,(k)}-u^{\pm,(k)}\|_{\mathcal{H}_{\varepsilon}}$
$\displaystyle\leqslant$
$\displaystyle\eta(\overline{\boldsymbol{w}}^{\pm;(k)};\varepsilon)\,,$

then

$\displaystyle\|\widetilde{e}^{\pm,(2k+1)}\|_{\mathcal{H}_{\varepsilon}}$
$\displaystyle\leqslant$
$\displaystyle\|e^{\pm,(2k+1)})\|_{\mathcal{H}_{\varepsilon}}+\|N^{\pm,(2k+1)}-u^{\pm,(2k+1)}\|_{\mathcal{H}_{\varepsilon}}\,.$

We then trivially deduce the convergence of the overall PINN-CSWR method.
Similar conclusions can be reached for the PINN-OSWR algorithms. In particular
in [6], the OSWR method is shown to be convergent in
$L^{2}(0,T;H^{3}(\Omega))\cap H^{3/2}(0,T;L^{2}(\Omega))$.

First, we use the standard technique to include the initial condition [17], by
searching for a trial network in the form

$\displaystyle T({\boldsymbol{w}},{\boldsymbol{x}},t)$ $\displaystyle=$
$\displaystyle
u_{0}({\boldsymbol{x}})+tN({\boldsymbol{w}},{\boldsymbol{x}},t)\,.$

For the sake of simplicity of the notations, we will yet denote by $N$ the
neural networks, which includes the contribution of the initial data. Notice
that this step is not essential, but allows to simplify the loss functions.
Hence at each Schwarz iteration we solve (25), by minimization the following
two local “independent” loss functions $\mathcal{L}^{\pm}$, for some positive
parameters $\lambda_{\textrm{Int}}^{\pm}$, $\lambda_{\textrm{Ext}}^{\pm}$ and
where we have denoted
${\boldsymbol{w}}=({\boldsymbol{w}}^{-},{\boldsymbol{w}}^{+})$. In particular,
we benefit from local training processes (subdomain-dependent), which allows
us to potentially avoid using the stochastic gradient method or$/$and improve
its convergence. Typically, the mini-batches would actually correspond to
training points for the local loss functions under consideration. At Schwarz
iteration $k$, we hence minimize

$\displaystyle\left.\begin{array}[]{lcl}\mathcal{L}^{\pm,(k)}({\boldsymbol{w}})&=&\big{\|}\partial_{t}N^{\pm,(k)}({\boldsymbol{w}}^{\pm},\cdot,\cdot)-\nu({\boldsymbol{x}})\triangle
N^{\pm,(k)}({\boldsymbol{w}}^{\pm},\cdot,\cdot)-{\boldsymbol{a}}({\boldsymbol{x}})\cdot\nabla
N^{\pm,(k)}({\boldsymbol{w}}^{\pm},\cdot,\cdot)\\\
&&-r({\boldsymbol{x}})N^{\pm,(k)}({\boldsymbol{w}}^{\pm},\cdot,\cdot)\big{\|}_{L^{2}(\Omega_{\varepsilon}^{\pm}\times(0,T))}+\lambda_{\textrm{Ext}}^{\pm}\big{\|}N^{\pm,(k)}({\boldsymbol{w}}^{\pm},\cdot,\cdot)\big{\|}_{L^{2}(\Lambda_{\varepsilon}^{\pm}\times(0,T))}\\\
&&+\lambda_{\textrm{Int}}^{\pm}\big{\|}\mathcal{T}_{\pm}N^{\pm,(k)}({\boldsymbol{w}}^{\pm},\cdot,\cdot)-\mathcal{T}_{\pm}N^{\mp,(k-1)}(\overline{{\boldsymbol{w}}}^{\pm},\cdot,\cdot)\big{\|}_{L^{2}(\Gamma_{\varepsilon}^{\pm}\times(0,T))}\,,\end{array}\right.$

where $\overline{{\boldsymbol{w}}}^{\pm}$ were computed at the Schwarz
iteration $k-1$. Recall that practically the loss functions are numerically
evaluated by approximating the norm using training points, typically randomly
chosen in $\Omega^{\pm}_{\varepsilon}$. This method allows for a complete
spatial decoupling of the problem over 2 (or arbitrary number of) subdomains.
Finally, the reconstructed solution is hence defined as
$N_{|\Omega_{\varepsilon}^{\pm}}(\cdot,t)=N^{\pm}(\cdot,t)$ for all
$t\geqslant 0$. More specifically

$\displaystyle\lim_{k\rightarrow+\infty}\big{\|}\hskip
5.69046pt\|N^{+,(k)}(\overline{{\boldsymbol{w}}}^{+},\cdot,\cdot)-N^{-,(k)}(\overline{{\boldsymbol{w}}}^{+},\cdot,\cdot)\|_{\infty,\overline{\Omega}_{\varepsilon}^{+}\cap\overline{\Omega}_{\varepsilon}^{-}}\big{\|}_{L^{2}(0,T)}=0\,,$
(27)

and we define the solution to the advection-diffusion-reaction equation as

$\displaystyle N$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}N^{+,(k^{\textrm{cvg}})}(\overline{{\boldsymbol{w}}}^{+},\cdot,\cdot),&\textrm{
in }\Omega_{\varepsilon}^{+}\times(0,T)\,,\\\
N^{-,(k^{\textrm{cvg}})}(\overline{{\boldsymbol{w}}}^{-},\cdot,\cdot),&\textrm{
in }\Omega_{\varepsilon}^{-}\times(0,T)\,.\end{array}\right.$

Practically, in order to evaluate the loss functions, it is necessary to
compute the equation at some very large $N_{T}N_{\boldsymbol{x}}^{\pm}$
randomly chosen training points $\\{({\boldsymbol{x}}^{\pm}_{j},t)\\}_{n;j}$
in $\Omega_{\varepsilon}^{\pm}\times(0,T)$, as the $L^{2}-$norms are not
exactly performed. From the point of view of the computation of the loss
function (requiring the evaluation of the PDE at the training points), the
algorithm is hence trivially embarrassingly parallel. From the optimization
point of view, the method now requires minimizing two loss functions.
Naturally, the computation of the minimization problems is embarrassingly
parallel as the two IBVPs are totally decoupled. As we are now considering two
IBVPs on smaller spatial domains, we can locally adapt the depth of the local
networks.
It is important to mention that, unlike SWR methods combined with standard
numerical (finite-difference, -volume, -elements, pseudospectral) methods, for
which convergence can be proven, the combination of SWR and PINN methods will
not necessarily ensure convergence to zero of the residual history. This is
due to the fact that from one Schwarz iteration to the next, the reconstructed
solutions may slightly differ as the minima obtained by minimization of the
local loss functions will a priori slightly differ. This fact is actually
inherent to the NN-based method. However, we expect the residual history to be
small from a practical point of view and for loss functions sufficiently
small. In addition to this argument, let us mention that the transmission
condition is naturally not exactly satisfied if it is included in the loss
function. A large enough weight can, for instance, be imposed on the
transmission constraint to ensure that it is accurately satisfied.

###  2.2 About the interest of using SWR DDM for NN-based algorithms

The estimation of the loss function using the direct PINN-method for solving
local PDE is trivially embarrassingly parallel, as the estimation is
independently performed at any given training point. However, this associated
minimization problem (batch-size related) is not locally specific, and
Stochastic Gradient Method (SGM) is hence an essential tool. In the proposed
approach, the local loss functions which are evaluated have specific meanings;
and allows to get accurate approximations of the solution in any given
subdomain.

The SWR method is a domain decomposition method in space for solving PDE.
Using standard advection-diffusion-reaction equation solvers, the main
algorithmic costs are the loss function estimations and the computation of
solutions to linear systems at each time iteration, involved in implicit or
semi-implicit stable schemes [6]. The latter has a polynomial complexity
$O(N^{\alpha})$, where $1<\alpha\leqslant 3$ is typically dependent on the
structure of the matrix involved in the linear system. Using a PINN approach,
there are naturally no more linear systems to solve “to estimate” the
solution. Instead, an optimization algorithm is necessary to parameterize the
NN-based solution. Denoting by $N_{{\boldsymbol{x}};t}$ the (a priori very
large) number of space-time training points to construct the loss function,
and $N_{W}$ the total number of parameters. The computation of solution using
the direct PINN method is decomposed into two parts:

  * 1.

Estimation of the loss function, with a complexity $O(N_{{\boldsymbol{x}};t})$
with $N_{{\boldsymbol{x}};t}\gg 1$. This step is hence embarrassingly parallel
for local PDE (with or without combination with the SWR method.)

  * 2.

Minimization of the loss function with a complexity $O(N_{W}^{p})$ for
$p\in(1,3)$. Typically stochastic gradient methods [19, 20, 21] are used to
deal with possibly high dimensionality (for very accurate solutions) of this
minimization problem and allows for a relatively efficient parallelization.

Within the framework of DDM and for two subdomains, the SWR-NN indeed requires
the embarrassingly parallel minimization of two independent loss functions
constructed using local training. The empirical argument which justifies the
proposed methodology is as follows. The structure and complexity of the
solution is thought to be “divided” in the two (much more in practice of
course) spatial subdomains. As a consequence, based on the local structure of
the solutions in $\Omega^{\pm}_{\varepsilon}$, the depth of the local neural
networks $N^{\pm}$ can then be adapted$/$reduced compared to the one-domain
PINN approach with one unique neural network. The extreme case in that matter,
would be a domain decomposition into small finite volumes, where the solution
would be approximated by a constant (cell-center finite volume method) that is
0-depth NNs, even if the overall solution has a complex spatial structure.
Naturally, the larger the subdomain size, the deeper the depth of the searched
local neural network associated to this subdomain. For two subdomains, the
minimization step within the SWR-NN consists in solving in parallel

$\displaystyle\overline{{\boldsymbol{w}}}^{\pm,(k)}$ $\displaystyle=$
$\displaystyle\textrm{argmin}_{{\boldsymbol{w}}\in
W^{\pm}}\mathcal{L}^{\pm,(k)}({\boldsymbol{w}})\,,$

rather than (for direct method)

$\displaystyle\overline{{\boldsymbol{w}}}$ $\displaystyle=$
$\displaystyle\textrm{argmin}_{{\boldsymbol{w}}\in
W}\mathcal{L}({\boldsymbol{w}})\,,$

where $N_{W^{\pm}}:=\\#W^{\pm}\leqslant N_{W}:=\\#W$. That is, it is possible
to decompose the minimization problem in several spatial subregions, where the
spectral structure of the solution can be very different from on subdomain to
the next, requiring locally smaller depths than using a unique global deep NN.
Hence, for SWR-NN method we have to perform the following tasks.

  * 1.

Estimation of the loss functions with a complexity
$O(N^{\pm}_{{\boldsymbol{x}};T})$. This step is hence embarrassingly parallel
within the SWR method and allows to deal with local (subdomain-dependent)
training points.

  * 2.

Minimization of the local loss functions with a complexity
$O(N_{W_{+}}^{\alpha})+O(N^{\alpha}_{W_{-}})$ for $\alpha\in(1,3)$, where in
principle $N_{W_{\pm}}\leqslant N_{W}$.

The downsize of SWR-NN methods is that it requires to repeat
$k^{\textrm{cvg}}$ times (that until convergence) the computation of the
uncoupled systems. Unlike standard SWR-DDM, where the gain is on the
computation of local linear systems of smaller size, the main interest is that
we locally solve local (and less complex) minimization problems, where we
expect the size of the search space to be smaller.

Notice that the SWR-DDM allows for an embarrassing parallelization of the
overall PINN PDE solver. Indeed, unlike the standard computation of the
(local) minima of the loss function, which requires non-trivial non-
embarrassingly parallelization, the proposed approach allows for the
embarrassingly parallel computation of minima of local loss functions. Three
levels of parallelization are then possible

  * 1.

Trivial parallelization of the estimation of the local loss functions.

  * 2.

Embarrassingly parallel computation of the minima of the local loss functions.

  * 3.

In addition, the minimum of a local loss function can also be performed in
parallel using the domain decomposition method for PINN, as proposed in see
[15].

From the computational point of view, the SWR-PINN algorithm allows i) to
adapt the depth of (most of) local NNs compared to using one unique (global)
NN, and ii) to estimate the local loss functions using local subdomain-
dependent training points and potentially allows for using direct (none-
stochastic) gradient methods for a sufficiently large number of subdomains.
This step is the analog of the reduction of the size of the linear systems to
be solved (scaling effect) within standard SWR when are applied as real space
solvers [1, 2, 3].

We here summarize the overall computational complexity of SWR-PINN and direct
PINN methods.

  * 1.

Direct approach: $O(N_{{\boldsymbol{x};t}})+O(N^{P}_{W})$. In this case, we
expect $N_{W}$ to be large, and $P$ depends on the used optimization
algorithm.

  * 2.

SWR approach:
$O\big{(}k^{\textrm{cvg}}(N^{P}_{W_{+}}N^{P}_{W_{-}})\big{)}+O(N^{-}_{{\boldsymbol{x};t}})+O(N^{+}_{{\boldsymbol{x};t}})$.
In this case, we expect $N_{W^{\pm}}<N_{W}$. As $P$ is strictly greater than
$1$ there is a scaling effect which makes this approach potential more
efficient. Moreover the prefactor is also though to be much smaller using SWR
methods. Practically it is required for $k^{\textrm{cvg}}$ to be small enough.
As it is well-known the choice of the transmission conditions is a crucial
element to minimize $k^{\textrm{cvg}}$. Dirichlet transmission conditions is
known to provide very slow convergence. At the opposite of the spectrum and
for wave-like equations, Dirichlet-to-Neumann like transmission conditions are
known to provide extremely fast convergence, but can be computationally
complex to approximate. Another way to accelerate the convergence of the SWR
algorithm consists in increasing the subdomain overlap (that is increase
$\varepsilon$) For the advection-diffusion-reaction equation, optimized SWR
method, based on optimized Robin-transmission conditions is a good compromise
between convergence rate and computational complexity [5]. As specified above,
the computation of the loss function is embarrassingly parallel unlike the
direct approach.

###  2.3 Nonlocal operator

We have argued above that the use of SWR methods allows for an efficient
parallel computation of the overall loss functions through the efficient
estimation (using local training points) of local loss functions. We show
below that whenever nonlocal terms are present in the equation, the efficiency
of the SWR-PINN method is not deteriorated by those terms. In the following,
we assume that the equation contains a nonlocal operator $F$, typically
defined as a convolution product:

  * 1.

$F(u)=\partial_{x}^{\alpha}u$, with $\alpha>0$, fractional derivative in
space, modeling nonlocal effect [22]. The latter is actually defined as a
convolution. We refer to [23] for details.

  * 2.

$F(u)=\rho*_{\boldsymbol{x}}u$ where $*_{\boldsymbol{x}}$ denotes the spatial
convolution product, a nonlocal potential for some given function $\rho$.

We consider the equation on a truncated domain $\Omega\subset\mathbb{R}^{d}$
with boundary $\Gamma$, as follows

$\displaystyle\left\\{\begin{array}[]{l}\partial_{t}u=\nu({\boldsymbol{x}})\triangle
u+F(u),\,{\boldsymbol{x}}\in\Omega,\,t>0,\\\
u({\boldsymbol{x}},0)=u_{0}({\boldsymbol{x}}),\,{\boldsymbol{x}}\in\Omega,\\\
u=0,\,{\boldsymbol{x}}\in\Gamma\,,\end{array}\right.$ (32)

and such that $F$ is defined as a convolution product in space

$\displaystyle
F(u)=u*_{\boldsymbol{x}}\rho=\int_{\Omega}u({\boldsymbol{x}}-{\boldsymbol{y}},t)\rho({\boldsymbol{y}})d{\boldsymbol{y}}\,,$

Then the SWR-PINN scheme reads

$\displaystyle\left\\{\begin{array}[]{lcl}\partial_{t}N^{\pm,(k)}&=&\nu({\boldsymbol{x}})\triangle
N^{\pm,(k)}+\int_{\Omega_{\varepsilon}^{\pm}}N^{\pm,(k)}({\boldsymbol{w}}^{\pm},{\boldsymbol{x}}-{\boldsymbol{y}},t)\big{)}\rho({\boldsymbol{y}})d{\boldsymbol{y}}\\\
&&+\int_{\Omega_{\varepsilon}^{\mp}}N^{\mp,(k-1)}(\overline{{\boldsymbol{w}}}^{\mp},{\boldsymbol{x}}-{\boldsymbol{y}},t)\big{)}\rho({\boldsymbol{y}})d{\boldsymbol{y}},\,\hbox{
in }\Omega_{\varepsilon}^{\pm}\times(0,T),\\\
N^{\pm,(k)}(\cdot,0)&=&u_{0}^{\pm},\,\hbox{ in }\Omega_{\varepsilon}^{\pm},\\\
\mathcal{T}_{\pm}N^{\pm,(k)}&=&\mathcal{T}_{\pm}N^{\mp,(k-1)},\,\hbox{ on
}\Gamma_{\varepsilon}^{\pm}\times(0,T),\\\ N^{\pm,(k)}&=&0,\,\hbox{ on
}\Lambda\backslash\Gamma_{\varepsilon}^{\pm}\times(0,T)\,,\end{array}\right.$
(38)

where $\overline{\boldsymbol{w}}^{\pm}\in W_{\pm}$ was computed at the
previous Schwarz iteration, with some transmission operator
$\mathcal{T}_{\pm}$. Hence in this case, we still have to minimize local loss
functions

$\displaystyle\left.\begin{array}[]{lcl}\mathcal{L}^{\pm,(k)}({\boldsymbol{w}})&=&\big{\|}\partial_{t}N^{\pm,(k)}({\boldsymbol{w}}^{\pm},\cdot,\cdot)-\nu({\boldsymbol{x}})\triangle
N^{\pm,(k)}({\boldsymbol{w}}^{\pm},\cdot,\cdot)\\\
&&-\int_{\Omega_{\varepsilon}^{\pm}}N^{\pm,(k)}({\boldsymbol{w}}^{\pm},{\boldsymbol{x}}-{\boldsymbol{y}},t)\big{)}\rho({\boldsymbol{y}})d{\boldsymbol{y}}\\\
&&-\int_{\Omega_{\varepsilon}^{\mp}}N^{\mp,(k-1)}(\overline{{\boldsymbol{w}}}^{\mp},{\boldsymbol{x}}-{\boldsymbol{y}},t)\big{)}\rho({\boldsymbol{y}})d{\boldsymbol{y}}\big{\|}_{L^{2}(\Omega_{\varepsilon}^{\pm}\times(0,T))}\\\
&&+\lambda_{\textrm{Int}}^{\pm}\big{\|}\mathcal{T}_{\pm}N^{\pm,(k)}({\boldsymbol{w}}^{\pm},\cdot,\cdot)-\mathcal{T}_{\pm}N^{\mp,(k-1)}(\overline{{\boldsymbol{w}}}^{\mp},\cdot,\cdot)\big{\|}_{L^{2}(\Gamma_{\varepsilon}^{\pm}\times(0,T))}\\\
&&+\lambda_{\textrm{Ext}}^{\pm}\big{\|}N^{\pm,(k)}({\boldsymbol{w}}^{\pm},\cdot,\cdot)\big{\|}_{L^{2}(\Lambda_{\varepsilon}^{\pm}\times(0,T))}\,.\end{array}\right.$

Practically, we can approximate the convolution product as follows. Denoting
by $\\{{\boldsymbol{x}}_{j}\\}_{j\in\mathcal{J}^{\pm}}\in\Omega^{\pm}$ the
local spatial training points, for
${\boldsymbol{x}}\in\Omega_{\varepsilon}^{\pm}$, we approximate

$\displaystyle\int_{\Omega_{\varepsilon}^{\pm}}N^{\pm,(k)}({\boldsymbol{w}}^{\pm},{\boldsymbol{x}}-{\boldsymbol{y}},t)\big{)}\rho({\boldsymbol{y}})d{\boldsymbol{y}}+\int_{\Omega_{\varepsilon}^{\mp}}N^{\mp,(k-1)}(\overline{{\boldsymbol{w}}}^{\mp},{\boldsymbol{x}}-{\boldsymbol{y}},t)\big{)}\rho({\boldsymbol{y}})d{\boldsymbol{y}}$

by

$\displaystyle\left.\begin{array}[]{l}\sum_{j\in\mathcal{J}^{\pm}}c_{j}F\big{(}N^{\pm,(k)}({\boldsymbol{w}}^{\pm},{\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{j},t)\big{)}+\sum_{j\in\mathcal{J}^{\mp}}c_{j}F\big{(}N^{\mp,(k-1)}(\overline{{\boldsymbol{w}}}^{\mp},{\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{j},t)\big{)}\,,\end{array}\right.$

for some weights $\\{c_{j}\\}_{j\in\mathcal{J}^{\pm}}$.

As it was discussed above the interest of using a DDM is to decompose the
training and search of the local solution over smaller set of parameters.
However, whenever the equation is nonlocal, it is necessary to extend the
search of the parameters in the global computational domains. More
specifically, for local equations, the local NN-solution in
$\Omega_{\varepsilon}^{+}$ (resp. $\Omega_{\varepsilon}^{-}$) only requires
parameters in $W_{+}$ (resp. $W_{-}$). However, if the equation is nonlocal,
in order to construct the solution in $\Omega_{\varepsilon}^{\pm}$, we have in
principle to search the NN parameters in all $W$, containing both $W_{+}$ and
$W_{-}$, as the solution values in $\Omega_{\varepsilon}^{\pm}$ depend on
values of the solution in $\Omega^{\mp}_{\varepsilon}$. This problem would
also occur to construct the loss function, using the direct PINNs method,
within the term

$\displaystyle\int_{\Omega}N({\boldsymbol{w}},{\boldsymbol{x}}-{\boldsymbol{y}},t)\big{)}\rho({\boldsymbol{y}})d{\boldsymbol{y}}\,.$

The SWR-PINN method allows to deal with this issue, as at Schwarz iteration
$k$, the loss function in $\Omega^{\pm}_{\varepsilon}$ is evaluated through
the solution in $\Omega_{\varepsilon}^{\pm}$ at the previous Schwarz iteration
$(k-1)$ from $N^{\mp}(\overline{{\boldsymbol{w}}}^{\mp},{\boldsymbol{x}},t)$
thanks to the previously evaluated parameter
$\overline{{\boldsymbol{w}}}^{\mp}$.

###  2.4 How about a non-iterative domain decomposition in space?

The domain decomposition method which is derived in this paper is a SWR-in
space method which is an iterative method allowing for the convergence of the
decomposed solution towards the exact solution of the PDE under consideration.
The main weakness of this DDM is the fact that the decoupled system has to be
solved several times (iterative method). It is hence natural to ask if a
“similar spatial domain decomposition”, but non-iterative, is possible.
In this goal, we decompose the domain $\Omega$ as above:
$\Omega_{\varepsilon}^{\pm}$, with or without overlap
($\Omega_{\varepsilon}^{+}\cap\Omega_{\varepsilon}^{-}=\emptyset$ or
$\Omega_{\varepsilon}^{+}\cap\Omega_{\varepsilon}^{-}\neq\emptyset$), with
$\varepsilon\geqslant 0$ and consider (32). That is we search for a solution
of the form
$N_{|\Omega_{\varepsilon}^{\pm}}(\cdot,t)=N^{\pm}(\overline{{\boldsymbol{w}}}^{\pm}\cdot,t)$
for all $t\geqslant 0$ such that

$\displaystyle\left\\{\begin{array}[]{lcl}\partial_{t}N^{\pm}&=&\nu({\boldsymbol{x}})\triangle
N^{\pm}+\int_{\Omega_{\varepsilon}^{\pm}}N^{\pm}({\boldsymbol{w}}^{\pm},{\boldsymbol{x}}-{\boldsymbol{y}},t)\big{)}\rho({\boldsymbol{y}})d{\boldsymbol{y}}\\\
&&+\int_{\Omega_{\varepsilon}^{\mp}}N^{\mp}({\boldsymbol{w}}^{\mp},{\boldsymbol{x}}-{\boldsymbol{y}},t)\big{)}\rho({\boldsymbol{y}})d{\boldsymbol{y}},\,\hbox{
in }\Omega_{\varepsilon}^{\pm}\times(0,T),\\\ N^{\pm}&=&0,\,\hbox{ on
}\Gamma\times(0,T)\,,\end{array}\right.$ (44)

where $\overline{\boldsymbol{w}}^{\pm}\in W_{\pm}$. The PINN-method then
consists in solving

$\displaystyle\overline{{\boldsymbol{w}}}^{\pm}$ $\displaystyle=$
$\displaystyle\textrm{argmin}_{{\boldsymbol{w}}\in
W^{\pm}}\mathcal{L}^{\pm}({\boldsymbol{w}})\,,$

where

$\displaystyle\left.\begin{array}[]{lcl}\mathcal{L}^{\pm}({\boldsymbol{w}})&=&\big{\|}\partial_{t}N^{\pm,(k)}({\boldsymbol{w}}^{\pm},\cdot,\cdot)-\nu({\boldsymbol{x}})\triangle
N^{\pm,(k)}({\boldsymbol{w}}^{\pm},\cdot,\cdot)\\\
&&-\int_{\Omega_{\varepsilon}^{\pm}}N^{\pm}({\boldsymbol{w}}^{\pm},{\boldsymbol{x}}-{\boldsymbol{y}},t)\big{)}\rho({\boldsymbol{y}})d{\boldsymbol{y}}\\\
&&-\int_{\Omega_{\varepsilon}^{\mp}}N^{\mp}({\boldsymbol{w}}^{\mp},{\boldsymbol{x}}-{\boldsymbol{y}},t)\big{)}\rho({\boldsymbol{y}})d{\boldsymbol{y}}\big{\|}_{L^{2}(\Omega_{\varepsilon}^{\pm}\times(0,T))}\\\
&&+\lambda_{\textrm{Ext}}^{\pm}\big{\|}N^{\pm,(k)}({\boldsymbol{w}}^{\pm},\cdot,\cdot)\big{\|}_{L^{2}(\Gamma\times(0,T))}\,.\end{array}\right.$

Therefore, in this case, we still have to minimize local loss functions.
However, there are 2 main issues:

  * 1.

Even if $\varepsilon$ is taken equal to zero, the decoupling of the solution
in the two subdomains, naturally induces a discontinuity at the subdomain
interfaces. It is possible to impose additional compatibility conditions (to
be included in the loss function), in the form of continuity condition
$N^{+}=N^{-}$ at $\Gamma_{\varepsilon}^{\pm}$, differentiability, but the
reconstructed global solution $N$ (such that
$N_{|\Omega^{\pm}_{\varepsilon}}=N^{\pm}$), will obviously not be an
approximate solution to the equation under consideration. Moreover, the
compatibility conditions will induce a re-coupling of the two systems in the
spirit of the following item.

  * 2.

The two systems, in $N^{+}$ and $N^{-}$ are actually also coupled through the
nonlocal term. This effect is similar to the addition of a compatibility
condition described above. Hence, the computation of the loss functions
$\mathcal{L}^{\pm}$ would not be embarrassingly parallel anymore. This is not
an issue in the SWR framework; as in the latter case, say at Schwarz iteration
$k$, the nonlocal term uses the approximate solution at the Schwarz iteration
$k-1$, which is a known quantity.

Hence unlike the SWR-PINN method, for which (27) occurs

$\displaystyle
N^{+}_{\overline{\Omega}_{\varepsilon}^{+}\cap\overline{\Omega}_{\varepsilon}^{-}}(\overline{{\boldsymbol{w}}}^{+},\cdot,\cdot)\neq
N_{\overline{\Omega}_{\varepsilon}^{+}\cap\overline{\Omega}_{\varepsilon}^{-}}^{-}(\overline{{\boldsymbol{w}}}^{-},\cdot,\cdot)\,.$
(46)

##  3 Numerical Experiments

In this section, we propose basic experiments in order to numerically
illustrate the convergence of the overall method. The PINN-algorithm was
implemented using deep-learning and optimization toolboxes from matlab,
DeepXDE [24] and tensorflow [25]. Although relatively simple, these
experiments illustrate the proof of concept of the proposed strategy and not
to provide the best convergence possible (which will be the purpose of a
future work).
Experiment 1. We consider the standard advection-diffusion-reaction equation

$\displaystyle\partial_{t}u$ $\displaystyle=$ $\displaystyle
a\partial_{x}u+\nu\partial_{xx}u+ru\,,$ (47)

on $\Omega\times[0,T]=(-1,1)\times[0,0.25]$ with Dirichlet boundary conditions
at $x=\pm 1$ and such that $a=0.5$, $\nu=0.0$, $r=0.1$ and the initial
conditions is $u_{0}(x)=\exp(-30(x-0.1)^{2})$. We decompose the domain in two
subdomains $\Omega^{-}_{\varepsilon}=(-1,\varepsilon/2)$ and
$\Omega^{+}_{\varepsilon}=(-\varepsilon/2,1)$ with $\varepsilon=0.1$. We here
use the Classical Schwarz Waveform Relaxation method, based on Dirichlet
transmission conditions
$N^{\pm,(k)}(\cdot,\pm\varepsilon/2,\cdot)=N^{\mp,(k-1)}(\cdot,\mp\varepsilon/2,\cdot)$,
where $N^{\pm}$ are the two local NN defined in $\Omega_{\varepsilon}^{\pm}$.
We consider the following data: the NN have both $5$ layers, with $25$ neurons
each. We select $5000$ internal collocation points. We also use a local SGM
with $5$ epochs and mini-batches size of $500$. In the gradient method the
learning rate with decay rate of $5\times 10^{-3}$ starting at $10^{-3}$. We
reconstruct the overall solution using a total of $5000$ prediction points.
Initially we take $N^{\pm,(0)}=0$. We report the reconstructed solution after
the first SWR iteration (resp. converged SWR algorithm) in Fig. 1 (Left)
(resp. 1 (Right)) from two the local solutions in $\Omega_{\varepsilon}^{\pm}$
at final time $T=0.25$, and overlapping zone of size $1/99$.

Figure 1: Experiment 1. (Left) Reconstructed solutions after the first
Schwartz iteration (Right) Reconstructed solutions at convergence of the SWR
method.

The SWR convergence rate is defined as the slope of the logarithm of the
residual history according to the Schwarz iteration number, that is
$\\{(k,\log(\mathcal{E}^{(k)}))\,:\,k\in\mathbb{N}\\}$, with (for 2
subdomains)

$\displaystyle\mathcal{E}^{(k)}:=\sum_{i=1}^{2}\big{\|}\hskip
5.69046pt\|N^{+,(k)}-N^{-,(k)}\|_{\infty;\overline{\Omega}_{\varepsilon}^{+}\cap\overline{\Omega}_{\varepsilon}^{-}}\big{\|}_{L^{2}(0,T)}\leqslant\delta^{\textrm{Sc}},$
(48)

$\delta^{\textrm{Sc}}$ being a small parameter.
We report in Fig 2 (Left) the graph of convergence of the stochastic gradient
methods applied to each local loss functions. Notice that each “oscillation”
corresponds to a new Schwarz iteration. We report in Fig. 2 (Right) the graph
of convergence of the SWR-method in the form of the residual history in the
overlapping zone. The combination of convergent SWR methods with standard
numerical (finite element, finite-difference) methods for which there is a
uniform convergence to zero of the residual history as a function of the
Schwarz iterations.

Figure 2: Experiment 1. (Left) Converged reconstructed solutions (Right) Graph
of local loss functions.

Notice that, unlike the converged SWR method combined with standard numerical
methods where the residual goes to zero when $k$ goes to infinity, the
residual does not exactly go to zero. This is due to the fact that from one
Schwarz iteration to the next, the (local) solution are obtained by
constructing “new” local minima as the local loss functions are small but not
null, and hence change from one iteration to the next.

Experiment 2. In the following experiment, we implement a Robin-SWR method for
solving (47), which is expected to provide better convergence than CSWR [5].
As it was discussed in [5] and recalled above, the optimized SWR (and more
generally Robin-SWR) methods is convergent, even without overlap, that is when
$\varepsilon$ is null. We consider the same equation as above with $a=1$,
$r=0$, $\nu=5\times 10^{-2}$ and $T=0.5$. The initial conditions is
$u_{0}(x)=10\exp(-4x^{2})+$. We decompose the domain in two subdomains
$\Omega_{\varepsilon}^{-}=(-4,0)$ and $\Omega_{\varepsilon}^{+}=(0,4)$, with
hence $\varepsilon=0$. The Robin transmission conditions
$(\lambda\pm\partial_{x})N^{\pm,(k)}(\cdot,0,\cdot)=(\lambda+\partial_{x})N^{\mp,(k-1)}(\cdot,0,\cdot)$,
where $N^{\pm}$ are the two local NN defined in $\Omega_{\varepsilon}^{\pm}$
and we have taken $\lambda=5$. We consider the following data: the NN have
both $8$ layers, with $30$ neurons each. We select $5000$ internal collocation
points. We also use local SGM with $5$ epochs and mini-batches size of $500$.
In the gradient method the learning rate with decay rate of $10^{-3}$ starting
at $10^{-3}$. We reconstruct the overall solution using a total of $5000$
prediction points. Initially we take $N^{\pm,(0)}=0$. We report the
reconstructed solution after the first SWR iteration (resp. converged SWR
algorithm) in Fig. 3 (Left) (resp. 3 (Right)) from two the local solutions in
$\Omega_{\varepsilon}^{\pm}$ at final time $T=0.5$.

Figure 3: Experiment 2. (Left) Reconstructed solutions after the first
Schwartz iteration (Right) Reconstructed solutions at convergence of the SWR
method.

We next report in Fig 4 (Left) the graph of convergence of the stochastic
gradient methods applied to each local loss functions. We report in Fig. 4
(Right) the graph of convergence of the SWR-method in the form of the residual
history in the overlapping zone.

Importantly, we observe that Robin-SWR-PINN still converges even if the two
subdomains do not overlap.

Figure 4: Experiment 2. (Left) Converged reconstructed solutions (Right) Graph
of local loss functions.

Experiment 2bis. In the following non-overlapping 2-domain Robin-SWR
experiment, we now consider that the diffusion coefficient is space-dependent;
more specifically, $\nu(x)=0.1$ (resp. $2.5\times 10^{-2}$) for
$x\in\Omega_{\varepsilon}^{+}$ (resp. $x\in\Omega_{\varepsilon}^{-}$). The
rest of the data are as follows: $a=1$, $r=0$, $\nu=5\times 10^{-2}$, $T=0.1$,
$\lambda=5$. The initial condition is given by
$u_{0}(x)=10\exp(-3(x-1)^{2})\cos(10(x-1)^{2})$ and is such that the solution
has a very different structure in the two subdomains. We want here to
illustrate the ability of the derived approach to select different depths of
the local neural networks, depending on the structure of the solution: in
$\Omega^{-}_{\varepsilon}$ (resp. $\Omega_{\varepsilon}^{+}$) the solution is
mainly null (resp. oscillatory), except close to the interface. The two
subdomains are $\Omega_{\varepsilon}^{-}=(-2,0)$ and
$\Omega_{\varepsilon}^{+}=(0,2)$, with hence $\varepsilon=0$. The two local NN
$N^{\pm}$, over $\Omega_{\varepsilon}^{\pm}$ have the following structure:
$N^{-}$ (resp. $N^{+}$) possesses $3$ (resp. $10$) layers and $10$ (resp.
$50$) neurons. The minimization process in $N^{-}$ is much more efficiently
performed than in $N^{+}$ with a relatively similar accuracy. As above, we
select $5000$ internal collocation points. We also use local SGM with $5$
epochs and mini-batches size of $5000$. In the gradient method the learning
rate with decay rate of $10^{-3}$ starting at $10^{-3}$. We reconstruct the
overall solution using a total of $20000$ prediction points. Initially we take
$N^{\pm,(0)}=0$. We report the reconstructed solution after the first SWR
iteration (resp. converged SWR algorithm) in Fig. 5 (Top-Left) (resp. 5 (Top-
Right)) from two the local solutions in $\Omega_{\varepsilon}^{\pm}$ at final
time $T=0.1$. We also zoom in (5, (Bottom-Left)), in the interface region to
better observe the SWR convergence. The local loss functions are represented
in Fig. 5 (Bottom-Right). We observe that roughly, the computational time to
perform the solution in $\Omega_{\varepsilon}^{+}$ was $2.5$ times faster than
in $\Omega_{\varepsilon}^{-}$.

Figure 5: Experiment 2bis. (Top-Left) Converged solutions after the first
Schwarz iteration. (Top-Right) Reconstructed converged solution. (Bottom-Left)
Zoom in of the converged solution in the interface region. (Bottom-Right)
Local loss function convergence.

Experiment 3. In this last experiment, we consider a two-dimensional
advection-diffusion equation on a square $[-1,1]^{2}\times[0,T]$.

$\displaystyle\partial_{t}u+{\bf a}\cdot\nabla u-\nu\triangle u$
$\displaystyle=$ $\displaystyle 0\,,$

with ${\bf a}=(-0.5,0)^{T}$ and $\nu=0.1$ and $T=2\pi/10$. The two subdomains
are $\Omega_{\varepsilon}^{+}=(-1,\varepsilon)\times(-1,1)$ and
$\Omega_{\varepsilon}^{+}=(-1,\varepsilon)\times(-1,1)$ where
$\varepsilon=0.1$; hence the interfaces are located at $\pm 0.1$. The initial
data is a Gaussian function
$u_{0}({\boldsymbol{x}})=\exp(-5\|{\boldsymbol{x}}\|^{2})$ and the final
computational time is $T=0.1$. A classical SWR algorithm is here combined with
the PINN method. On the other subdomain boundaries, we impose null Dirichlet
boundary conditions. The equation is solved using the library DeepXDE [24]
combined with tensorflow [25]. In each subdomain is used a neural network with
$3$ layers and $10$ neurons; Adam’s optimizer is used (learning rate
$10^{-3}$, epoch=$10^{3}$) along with $\tanh$ activation function. In Fig. 6
(Top), we report the initial data in $\Omega_{\varepsilon}^{\pm}$. In Fig. 6,
we represent the solution at the end of the first Schwarz iteration (Left) and
fully converged solution (Right) at final time $T_{f}$. In future works, we
will propose more advanced simulations. The corresponding code is available on
GitHub, where the interested reader could find all the relevant information
regarding the code.

Figure 6: Experiment 3. (Top) Initial data in $\Omega_{\varepsilon}^{\pm}$.
(Bottom-Left) Solutions after the first Schwarz iteration at time $T_{f}$ .
(Bottom-Right) Reconstructed converged solution at time $T_{f}$.

##  4 Conclusion

In this paper, we have derived a Schwarz Waveform Relaxed Physics-Informed
Neural Networks (SWR-PINN) method for solving advection-diffusion-reaction
equations in parallel. Some preliminary illustrating experiments are presented
to validate the approach.

###  4.1 Pros. and cons. of the SWR-PINN method

We summarize below the pros and cons of the proposed method.
Pros.

  * 1.

Embarrassingly parallelization of the local loss function training.

  * 2.

Parallel construction of local neural networks with adaptive depth and
complexity.

  * 3.

For convergent PINN algorithms, the SWR-PINN is convergent.

  * 4.

Flexible choice of the transmission conditions.

Cons.

  * 1.

As a fixed point method, SWR methods require several iterations.

  * 2.

The transmission conditions must be accurately satisfied through a
penalization term in the loss function in order to accurately implement the
SWR algorithm. Ideally, we should directly include the transmission within the
definition in the NN. This is possible, considering the CSWR (Dirichlet-based
transmission conditions) method and the following NN, $T^{\pm}$

$\displaystyle T^{\pm;(k)}({\boldsymbol{w}},{\boldsymbol{x}},t)$
$\displaystyle=$ $\displaystyle
N_{|\Gamma_{\varepsilon}^{\pm}\times(0,T)}^{\mp;(k-1)}(\overline{{\boldsymbol{w}}}^{\mp},{\boldsymbol{x}},t)+\big{(}{\bf
1}_{\Gamma_{\varepsilon}^{\pm}\times(0,T)}({\boldsymbol{x}},t)-1\big{)}N^{\pm;(k)}({\boldsymbol{w}},{\boldsymbol{x}},t)\,.$

  * 3.

Convergence or high precision of the overall algorithm can be hard to reach if
the PINN algorithm is not used with sufficiently high precision. Instable
numerical behavior can also be observed with the CSWR method.

###  4.2 Concluding remarks and future investigations

As far as we know, this paper is the first attempt to combine the SWR and PINN
methods. Although the theory of SWR-DDM is now well developed in terms of
convergence and convergence rate for different types of evolution PDE and
their approximation with finite difference and finite element methods, the
theory of convergence of PINN is not yet complete. Consequently, the
convergence of the overall SWR-PINN method is still subject to the proof of
convergence of the latter, which is largely empirically established. In future
works, we plan to focus on “real-life” experiments where the main benefits of
the SWR-PINN will be exhibited and illustrated.

## References

  * [1] X. Antoine and E. Lorin.  Multilevel preconditioning technique for Schwarz waveform relaxation domain decomposition method for real- and imaginary-time nonlinear Schrödinger equation.  Appl. Math. Comput., 336:403–417, 2018.
  * [2] X. Antoine and E. Lorin.  An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations.  Numer. Math., 137(4):923–958, 2017.
  * [3] X. Antoine, F. Hou, and E. Lorin.  Asymptotic estimates of the convergence of classical Schwarz waveform relaxation domain decomposition methods for two-dimensional stationary quantum waves.  ESAIM Math. Model. Numer. Anal., 52(4):1569–1596, 2018.
  * [4] X. Antoine and E. Lorin.  On the rate of convergence of Schwarz waveform relaxation methods for the time-dependent Schrödinger equation.  J. Comput. Appl. Math., 354:15–30, 2019.
  * [5] L. Halpern and J. Szeftel.  Optimized and quasi-optimal Schwarz waveform relaxation for the one-dimensional Schrödinger equation.  Math. Models Methods Appl. Sci., 20(12):2167–2199, 2010.
  * [6] M. Gander and L. Halpern.  Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems.  SIAM J. Numer. Anal., 45(2), 2007.
  * [7] M.J. Gander, L. Halpern, and F. Nataf.  Optimal convergence for overlapping and non-overlapping Schwarz waveform relaxation.  In Proceedings of the 11th International Conference on Domain decomposition, pages 27–36, 1999.
  * [8] X. Antoine and E. Lorin.  Asymptotic convergence rates of Schwarz waveform relaxation algorithms for Schrödinger equations with an arbitrary number of subdomains.  Multiscale Science and Engineering, 1(1):34–46, 2019.
  * [9] M.J. Gander.  Optimal Schwarz waveform relaxation methods for the one-dimensional wave equation.  SIAM J. Numer. Anal., 41:1643–1681, 2003.
  * [10] M.J. Gander.  Optimized Schwarz methods.  SIAM J. Numer. Anal., 44:699–731, 2006.
  * [11] V. Dolean, P. Jolivet, and F. Nataf.  An introduction to domain decomposition methods: theory and parallel implementation.  2015\.
  * [12] M. Raissi, P. Perdikaris, and G.E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [13] G. Pang, L. Lu, and G.E. Karniadakis.  fPINNs: fractional physics-informed neural networks.  SIAM J. Sci. Comput., 41(4):A2603–A2626, 2019.
  * [14] L. Yang, D. Zhang, and G. E. Karniadakis.  Physics-informed generative adversarial networks for stochastic differential equations.  SIAM J. Sci. Comput., 42(1):A292–A317, 2020.
  * [15] A. D. Jagtap and G. E. Karniadakis.  Extended physics-informed neural networks (XPINNs): a generalized space-time domain decomposition based deep learning framework for nonlinear partial differential equations.  Commun. Comput. Phys., 28(5):2002–2041, 2020.
  * [16] A. Heinlein, A. Klawonn, M. Lanser, and J. Weber.  Combining machine learning and domain decomposition methods for the solution of partial differential equations—a review.  GAMM Mitteilungen, 44(1), 2021.
  * [17] I.E. Lagaris, A. Likas, and D.I. Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE Transactions on Neural Networks, 9(5):987–1000, 1998.
  * [18] B. Després.  Analyse numérique et neural networks.  Technical report, Université de Paris, https://www.ljll.math.upmc.fr/despres, 2021.
  * [19] H. Robbins and S. Monro.  A stochastic approximation method.  Ann. Math. Statistics, 22:400–407, 1951.
  * [20] L. Bottou.  Large-scale machine learning with stochastic gradient descent.  In Proceedings of COMPSTAT’2010, pages 177–186. Physica-Verlag/Springer, Heidelberg, 2010.
  * [21] S. Sun, Z. Cao, H. Zhu, and J. Zhao.  A survey of optimization methods from a machine learning perspective.  IEEE Transactions on Cybernetics, 50(8):3668–3681, 2020.
  * [22] A. Lischke, G. Pang, M. Gulian, F. Song, C. Glusa, X. Zheng, Z. Mao, W. Cai, M.M. Meerschaert, M. Ainsworth, and G. E. Karniadakis.  What is the fractional Laplacian? A comparative review with new results.  J. Comput. Phys., 404, 2020.
  * [23] E.C. De Oliveira and J.A. Tenreiro Machado.  A review of definitions for fractional derivatives and integral.  Mathematical Problems in Engineering, 2014, 2014.
  * [24] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis.  Deepxde: A deep learning library for solving differential equations, 2020\.
  * [25] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  Software available from tensorflow.org.