# Possible quantum gravity effects on the gravitational deflection of light

Xin Li1,3<EMAIL_ADDRESS>Zhe Chang2,3<EMAIL_ADDRESS>1Institute of
Theoretical Physics, Chinese Academy of Sciences, 100190 Beijing, China
2Institute of High Energy Physics, Chinese Academy of Sciences, 100049
Beijing, China
3Theoretical Physics Center for Science Facilities, Chinese Academy of
Sciences

###### Abstract

We investigate possible quantum gravity (QG) effects on the gravitational
deflection of light. Two forms of deformation of the Schwarzschild spacetime
are proposed. The first ansatz is a given Finslerian line element, it could be
regarded as a weak QG effect on the Schwarzschild spacetime. Starting from
this ansatz, we deduce the deflection angle of the light ray which passes a
weak gravitational source. The second ansatz could be regarded as a strong QG
effect on the Schwarzschild spacetime. The deflection angle of the light ray
which passes a weak gravitational source is deduced in this Riemannian
spacetime. This QG effect may distinguish the mixed light rays in the absence
of gravitational source by a “spectroscope” (the gravitational source). The
solutions of gravitational field equation in this Riemannian spacetime
indicate that the QG effect could be regarded as the vacuum energy and the
energy density of vacuum is related to the spacetime deformation parameter.

###### pacs:

04.60.Bc

##  I Introduction

The quantum theory of gravity (QG) has been studied for more than 70 years
Stachel . The most famous theories of quantum gravity, such as string theory
string  and loop quantum gravity loop , provide us some key features of QG.
Amelino-Camelia Amelino  listed possible candidates of QG effects: the
violation of Lorentz symmetry of Standard Model Garay  as well as discrete
symmetry (CPT symmetry) Ellis ; Kostelecky , Planck-scale fuzziness of
spacetime Amelino1 . Other possible QG effects include: deviations from Newton
s law at very short distances Hoyle , possible production of mini-black holes
Dimopoulos . The researches of such possible consequence of QG are refered as
quantum gravity phenomenology (QGP). The rapid progress in technology makes
experiments have the opportunities to test the sub-Planckian consequences of
QG scenarios.

The QGP covers a wide range of subjects. One of the most important QG effects
is the violation of the Lorentz invariance (LI)Mattingly . A feature of QG,
which is the most debated possibility for a quantum spacetime, manifests that
the spacetime in Planck-scale may be noncommutative Connes ; Majid . A huge
numbers of investigations of noncommutative spacetime manifest that the Lie-
algebra Poincare symmetries are either broken to a smaller symmetry (Lie
algebra or deformed into Hopf-algebra Majid1  symmetries). The LI violation
implies that the dispersion relations for elementary particles should be
modified. Moreover, studying on the dispersion relation is a convenient way
for physicists to test the departure from LI. In the past few years, Amelino-
Camelia and Smolin as well as their collaborators have developed the Doubly
Special Relativity (DSR) Amelino2 ; Amelino3 ; Amelino4 ; Smolin1 ; Smolin2
to take Planck-scale effects into account by introducing an invariant Planckin
parameter in the theory of Special Relativity. The general form of dispersion
relation for free particles in the DSR is of the form

$E^{2}=m^{2}+p^{2}+\sum_{n=1}^{\infty}\alpha_{n}(\mu,M_{p})p^{n}~{},$ (1)

where $p=\sqrt{\parallel\vec{p}\parallel^{2}}$, $\mu$ denotes a parameter of
the theory with mass scale and $M_{p}$ is the Planck mass. The modified
dispersion relations (MDR) have been tested through observations on gamma-ray
bursts and ultra-high energy cosmic raysJacobson . Girelli et al.Girelli
showed that the MDR can be incorporated into the framework of Finsler
geometry. The symmetry of the MDR was described in the Hamiltonian formalism.
The generators of symmetry commute with ${\cal M}(p)$ (here ${\cal
M}(p)=m^{2}$ gives the mass shell condition). The mass shell condition is
invariant under the deformed Lorentz transformations.

The research of Girelli et al. Girelli  gives a possible origin of MDR, which
means the quantum spacetime may have a Finslerian form. Randers space, as a
special kind of Finsler space, was first proposed by G. Randers Randers .
Within the framework of Randers space, modified dispersion relation has been
discussed RF , and the threshold of ultra high energy cosmic rays was
investigated UHC . The investigation of the isometric group of Finsler space
indicates that the Lorentz symmetry is broken in Finsler space Li .

Regardless of any particular theories of QG, we generally agree that the
spacetime in QG scenario should be deformed. In this paper, we start from a
giving deformed spacetime to investigate possible QG effects on the
gravitational deflection of light. This paper is organized as follows. In Sec.
II, we introduce a Finslerian line element which could be regarded as a weak
deformation from the Schwarzschild spacetime. Starting from such line element,
we investigate the trajectory of the light ray. In Sec. III, we introduce a
Riemannian line element which could be regarded as a strong deformation from
the Schwarzschild spacetime. The trajectory of the light ray in this
Riemannian spacetime is studied. The Einstein’s gravitational field equation
is presented in this Riemannian spacetime. The light ray deflected by a weak
gravitational source is studied. In Sec. IV, we give conclusions on the
possible QG effects on the gravitational deflection of light.

##  II Weak quantum gravity effect on the gravitational deflection of light

Instead of defining an inner product structure over the tangent bundle in
Riemann geometry, Finsler geometry is based on the so called Finsler structure
$F$ with the property $F(x,\lambda y)=\lambda F(x,y)$ for all $\lambda>0$,
where $x$ represents position and $y\equiv\frac{dx}{d\tau}$ represents
velocity. The Finsler metric is given asBook by Bao

$g_{\mu\nu}\equiv\frac{\partial}{\partial y^{\mu}}\frac{\partial}{\partial
y^{\nu}}\left(\frac{1}{2}F^{2}\right).$ (2)

Finsler geometry has its genesis in integrals of the form

$\int^{r}_{s}F(x^{1},\cdots,x^{n};\frac{dx^{1}}{d\tau},\cdots,\frac{dx^{n}}{d\tau})d\tau~{}.$
(3)

So that the Finsler structure represents the length element of Finsler space.

Following the calculus of variations, one get the geodesic equation of Finsler
spaceBook by Bao

$\frac{d^{2}\sigma^{\lambda}}{d\tau^{2}}+\gamma^{\lambda}_{\mu\nu}\frac{d\sigma^{\mu}}{d\tau}\frac{d\sigma^{\nu}}{d\tau}=\frac{d\sigma^{\mu}}{d\tau}\frac{d}{d\tau}\left(\log
F(\sigma,\frac{d\sigma}{d\tau})\right),$ (4)

where

$\gamma^{\lambda}_{\mu\nu}=\frac{g^{\lambda\alpha}}{2}\left(\frac{\partial
g_{\mu\alpha}}{\partial x^{\nu}}+\frac{\partial g_{\nu\alpha}}{\partial
x^{\mu}}-\frac{\partial g_{\mu\nu}}{\partial x^{\alpha}}\right)$ (5)

is the formal Christoffel symbols of the second kind with the same form of
Riemannian connection. The parallel transport has been studied in the
framework of Cartan connectionMatsumoto ; Antonelli ; Szabo . The notation of
parallel transport in Finsler manifold means that the length
$F\left(\frac{d\sigma}{d\tau}\right)$ is constant. Thus, the autoparallel
equation can be got from the equation (4)

$\frac{d^{2}\sigma^{\lambda}}{d\tau^{2}}+\gamma^{\lambda}_{\mu\nu}\frac{d\sigma^{\mu}}{d\tau}\frac{d\sigma^{\nu}}{d\tau}=0.$
(6)

Since the geodesic equation (4) is directly derived from the integral length
of $\sigma$

$L=\int F\left(\frac{d\sigma}{d\tau}\right)d\tau,$ (7)

the inner product
$\left(\sqrt{g_{\mu\nu}\frac{d\sigma^{\mu}}{d\tau}\frac{d\sigma^{\nu}}{d\tau}}=F\left(\frac{d\sigma}{d\tau}\right)\right)$
of two parallel transported vectors is preserved.

In a phenomenological approach, it is reasonable to suppose that the line
element of spacetime is given as

$\displaystyle F^{2}d\tau^{2}=$ $\displaystyle
adt^{2}-a^{-1}dr^{2}-r^{2}(\sin^{2}\theta d\theta^{2}+d\phi^{2})$ (8)
$\displaystyle+$
$\displaystyle\kappa\sqrt{a}\frac{GM}{r}\left(1-a\frac{J^{2}}{E^{2}r^{2}}\right)^{\frac{3}{4}}dt\sqrt{dtdr},$

where $a=1-\frac{2GM}{r}$, $\frac{GM}{r}\ll 1$, $M$ is the mass of
gravitational source, $E$ is the energy per unit mass of the particle, $J$ is
the angular momentum per unit mass of the particle and $\kappa$ is
dimensionless and constancy parameter which is function of the new physics
scale. The dimensionless parameter $\kappa$ just plays the role of the
measurement of QG effect. While $\kappa$ vanishes, the line element (8)
returns to the famous Schwarzschild metric. The line element (8) could be
regarded as the deformation of the Schwarzschild metric. Since the coefficient
of the non Riemannian term in (8) is isotropic, it is convenient to consider
that the motion of particle is confined in the plane of
$\theta=\frac{\pi}{2}$. Therefore, the trajectory of the particle is described
by the $r$ and $\phi$ coordinates, and the line element reduces to

$\displaystyle F^{2}d\tau^{2}=$ $\displaystyle
adt^{2}-a^{-1}dr^{2}-r^{2}d\phi^{2}$ (9) $\displaystyle+$
$\displaystyle\kappa\sqrt{a}\frac{GM}{r}\left(1-a\frac{J^{2}}{E^{2}r^{2}}\right)^{\frac{3}{4}}dt\sqrt{dtdr}.$

It should be noticed that the coefficients of the line element (9) do not
depend either on $r$ or $\phi$. By making use of the autoparallel geodesic
equation (6), we obtain two integrals of motion

$\displaystyle
a\dot{t}+\kappa\sqrt{a}\frac{3GM}{4r}\left(1-a\frac{J^{2}}{E^{2}r^{2}}\right)^{\frac{3}{4}}\sqrt{\dot{t}\dot{r}}$
$\displaystyle=$ $\displaystyle E,$ (10) $\displaystyle r^{2}\dot{\phi}$
$\displaystyle=$ $\displaystyle J,$ (11)

where a dot denotes the derivatives with respect to trajectory parameter
$\tau$. The deformation term in (9) is very small, here we can approximately
choose the two integral constants to be $E$ and $J$, respectively. In this
paper, we mainly consider the motion of light ray. There is an additional
constraint arising from $F=0$ for the motion of light ray

$a\dot{t}^{2}-a^{-1}\dot{r}^{2}-r^{2}\dot{\phi}^{2}+\kappa\sqrt{a}\frac{GM}{r}\left(1-a\frac{J^{2}}{E^{2}r^{2}}\right)^{\frac{3}{4}}\dot{t}^{\frac{3}{2}}\dot{r}^{\frac{1}{2}}=0.$
(12)

By making use of the equations (10), (11) and (12), we obtain an approximate
equation for $\dot{r}$, which is valid for $\frac{GM}{r}\ll 1$ (it was assumed
in the beginning of this section)

$\dot{r}=\sqrt{E^{2}-a\frac{J^{2}}{r^{2}}}\left(1-\kappa\frac{GM}{4r}\right).$
(13)

Combining (13) with equation (11), we obtain

$\left(\frac{1}{r^{2}}\frac{dr}{d\phi}\right)^{2}=\left(\frac{E^{2}}{J^{2}}-\frac{a}{r^{2}}\right)\left(1-\kappa\frac{GM}{2r}\right).$
(14)

In terms of the variable $u=\frac{GM}{r}$, the equation (14) changes as

$\left(\frac{du}{d\phi}\right)^{2}=\left(\left(\frac{EGM}{J}\right)^{2}-u^{2}(1-2u)\right)\left(1-\frac{\kappa
u}{2}\right).$ (15)

The derivatives $\frac{d}{d\phi}$ of the equation (15) gives

$\frac{d^{2}u}{d\phi}^{2}+u=3u^{2}-\frac{\kappa}{4}\left(\left(\frac{EGM}{J}\right)^{2}-3u^{2}\right)+O(u^{3}).$
(16)

Solving the equation (16) to order $u$, we get that

$u=u_{0}\cos\phi.$ (17)

Substituting the above solution into (15), we obtain the approximate solution
for $u_{0}$

$u_{0}=\frac{EGM}{J}=\frac{GM}{\xi},$ (18)

where $\xi$ is the minimum distance of the light ray to the gravitational
source with mass $M$. The closest approach of the light ray to the
gravitational source implies $\frac{dr}{d\phi}=0$. Then, deducing from the
equation (14), we obtain approximate solution for $\xi$$(=\frac{J}{E}$). This
is the reason for the second equation in (18). By making use of the first
order approximation (17), one may get the second order approximation of (16)

$u=u_{0}^{2}\left((1+\kappa/4)(1+\sin^{2}\phi)-\kappa/4\right).$ (19)

Thus, the solution of equation (16) is

$u=u_{0}\cos\phi+u_{0}^{2}\left((1+\kappa/4)(1+\sin^{2}\phi)-\kappa/4\right).$
(20)

At infinity ($u=0$), the solution (20) shows that the angle
$\phi=\pm(\frac{\pi}{2}+\alpha)$, and the small angle $\alpha$ satisfies the
constraint

$-u_{0}\sin\alpha+u_{0}^{2}\left((1+\kappa/4)(1+\cos^{2}\alpha)-\kappa/4\right)=0.$
(21)

Since the angle $\alpha$ is very small, the solution of the above equation
(21) is

$\alpha=u_{0}(2+\kappa/4)=(2+\kappa/4)\frac{GM}{\xi}.$ (22)

The two asymptotic directions differs from $\pi$ by the deflection angle

$\hat{\alpha}=2\alpha=(4+\kappa/2)\frac{GM}{\xi}$ (23)

In general relativity, the light passing a massive object $M$ at a minimum
distance $\xi$ suffers deflection, and the deflection angle (“Einstein angle”)
is $\hat{\alpha}_{E}=\frac{4GM}{\xi}$ Weinberg . The spacetime (8) deformed
from Schwarzschild spacetime implies a deformed deflection angle. The formula
for the deformed deflection angle (23) differs from $\hat{\alpha}_{E}$ by
$\kappa\frac{GM}{2\xi}$, which is proportional to the deformed parameter
$\kappa$.

In astronomical observations, the gravitational lensing surveys is used to
calculate the mass distribution that projected onto the sky. A large amount of
observations manifest that the expected gravitational lensing effects deducing
by the Einstein angle are not in accord with the experimental data. An example
is the full-sky data product for the Bullet Cluster 1E0657-558 Clowe . We wish
the deformed deflection angle (23) may account to these observations.

##  III Strong quantum gravity effect on the gravitational deflection of light

In Sec. II, we got a deformed deflection angle for the light ray. The
Finslerian metric (8) could be regarded as a weak deformation from the
Schwarzschild metric. In this section, we discuss the strong quantum gravity
effect on the gravitational deflection of light. Again, in a phenomenological
approach, we propose that the line element of spacetime is given as

$ds^{2}=adt^{2}-(4n+1)^{2}a^{-1}dr^{2}-r^{2}d\phi^{2},$ (24)

where the deformation parameter $n=0,1,2\cdots$. In ansatz (24), we already
confined the motion of particle in the plane of $\theta=\frac{\pi}{2}$. The
reason is the same within Sec. II. The ansatz (24) manifests a large
deformation from the Schwarzschild metric, it returns to the Schwarzschild
metric while the deformation parameter vanishes.

By making use of the autoparallel geodesic equation of Riemannian space
Weinberg , we obtain two integrals of motion

$\displaystyle a\dot{t}$ $\displaystyle=$ $\displaystyle E,$ (25)
$\displaystyle r^{2}\dot{\phi}$ $\displaystyle=$ $\displaystyle J.$ (26)

The additional constraint arising from $\frac{ds}{d\tau}=0$ for the motion of
light is

$a\dot{t}^{2}-(4n+1)^{2}a^{-1}\dot{r}^{2}-r^{2}\dot{\phi}^{2}=0$ (27)

By making use of the equations (25), (26) and (27), we get an approximate
equation for $\dot{r}$

$(4n+1)^{2}\dot{r}^{2}=E^{2}-a\frac{J^{2}}{r^{2}}.$ (28)

Combining this equation (28) with (26), we obtain that

$\left(\frac{4n+1}{r^{2}}\frac{dr}{d\phi}\right)^{2}=\frac{E^{2}}{J^{2}}-\frac{a}{r^{2}}.$
(29)

The closest approach of the light ray to the gravitational source $M$ implies
$\frac{dr}{d\phi}=0$. Then, deducing from the equation (29), we obtain
approximate solution for the distance of closest approach
$\xi$($=\frac{J}{E}$). Changing variable to $u=\frac{GM}{r}$, we obtain from
equation (14) that

$(4n+1)^{2}\left(\frac{du}{d\phi}\right)^{2}=\left(\frac{EGM}{J}\right)^{2}-u^{2}(1-2u).$
(30)

Calculating the derivatives $\frac{d}{d\phi}$ of the equation (30) gives

$\frac{d^{2}u}{d\phi^{2}}+\frac{u}{(4n+1)^{2}}=\frac{3u^{2}}{(4n+1)^{2}}.$
(31)

Noticing that $u=\frac{GM}{r}\ll 1$, the equation (31) has solution

$u=u_{0}\cos\frac{\phi}{4n+1}+u_{0}^{2}\left(1+\sin^{2}\frac{\phi}{4n+1}\right).$
(32)

Substituting the above solution into (30), we get the approximate solution for
$u_{0}$

$u_{0}=\frac{EGM}{J}=\frac{GM}{\xi}.$ (33)

At infinity ($u=0$), to first order in $u_{0}$, the solution (32) shows

$\phi=\pm(4n+1)\frac{\pi}{2}.$ (34)

The angle $\phi$ in (34) for different $n$ only differ from an integer times
of $2\pi$. The formula (34) means each light ray which corresponds to a given
parameter $n$ is mixed at infinity in the absence of the gravitational source
$M$, and the light rays all move with the same direction. When they suffer
from a week gravitational source, to second order in $u_{0}$, the solution
(32) implies

$\phi=\pm\left((4n+1)\frac{\pi}{2}+\alpha\right),$ (35)

and the small angle $\alpha$ satisfies the constraint

$-u_{0}\sin\frac{\alpha}{4n+1}+u_{0}^{2}\left(1+\cos^{2}\frac{\alpha}{4n+1}\right)=0.$
(36)

Since the angle $\alpha$ is very small, the solution of the above equation
(36) is

$\alpha=2(4n+1)u_{0}=2(4n+1)\frac{GM}{\xi}.$ (37)

The difference of two asymptotic directions
$\pm\left((4n+1)\frac{\pi}{2}+\alpha\right)$ differs from $(4n+1)\pi$ by the
deflection angle

$\hat{\alpha}=2\alpha=4(4n+1)\frac{GM}{\xi}.$ (38)

This formula (38) means the deflection angles for different $n$ are different.
Therefore, the weak gravitational source $M$ just plays the role of the
“spectroscope”. The mixed light rays come from infinity, refracted by the
“spectroscope” $M$ at the same distance of closest approach $\xi$, go into
different directions. And the deflection angle is in proportion to the
spacetime deformation parameter $n$.

One should notice from the formula (28) that the radial momentum of light rays
is modified. The ansatz (24) is a Riemannian line element. Therefore, the
motion of particle in this spacetime satisfies the Einstein’s gravitational
field equation. The Einstein’s gravitational field equation is taken the from

$R_{\mu\nu}=-8\pi
G\left(T_{\mu\nu}-\frac{1}{2}g_{\mu\nu}T^{\lambda}_{~{}\lambda}\right),$ (39)

where $R_{\mu\nu}$ is the Ricci tensor. The energy-momentum tensor of a
perfect fluid is taken the from as

$T^{\mu\nu}=-pg^{\mu\nu}+(p+\rho)U^{\mu}U^{\nu},$ (40)

here $U^{\mu}$ is the fluid four-velocity and satisfies
$g_{\mu\nu}U^{\mu}U^{\nu}=1$ , $p$ and $\rho$ are pressure and energy density
respectively. By making use of the line element (24), we show that the
combination of the $tt$ and $rr$ component of the field equation (39) gives

$p=-\rho,$ (41)

and the $\phi\phi$ component of the field equation (39) gives

$-1+\frac{1}{(4n+1)^{2}}=-4\pi Gr^{2}(\rho-p).$ (42)

The equation (41) implies that there is a quantum vacuum with energy density
$\rho$ outside the gravitational source. Combining the equation (41) with
(42), we obtain

$\rho=\frac{1}{8\pi Gr^{2}}\left(1-\frac{1}{(4n+1)^{2}}\right).$ (43)

The equation (43) indicates that the energy density of vacuum is the function
of spacetime deformation parameter $n$. From the point of view of particle
physics, the cosmological constant, a popular candidate of the dark energy,
naturally arises as an energy density of the vacuum Copeland . It implies that
this kind of QG effect (9) may arise from the influence of the dark energy.

##  IV Conclusion

In this paper, we investigated the possible quantum gravity effects on the
gravitational deflection of light. One of the most expected QG effect is the
deformation of spacetime geometry. In a phenomenological approach, we proposed
two deformations of the Schwarzschild spacetime. The first ansatz (8) is a
given Finslerian line element, it could be regarded as a weak QG effect on the
Schwarzschild spacetime. The deformation term (the non Riemannian term) is
very small. Starting from the ansatz (8), we deduced the deflection angle
(23). This deformed deflection angle may account to the observations of
gravitational lensing which can not be explained in the framework of general
relativity.

The second ansatz (24) could be regarded as a strong QG effect on the
Schwarzschild spacetime, for the deformation term $(4n+1)^{2}\geq 1$. Starting
from the ansatz (24), we deduced the deflection angle (38). This QG effect may
distinguish the mixed light rays in the absence of gravitational source by a
“spectroscope” $M$. The solutions of gravitational field equation in spacetime
(24) indicate that the QG effect could be regarded as the vacuum effect and
the energy density of vacuum is related to the spacetime deformation parameter
$n$. While $n$ vanishes, the spacetime with quantum effect (24) returns to the
Schwarzschild spacetime and the vacuum energy vanishes, it means that there is
nothing exist outside the gravitational source $M$.

Acknowledgements

We would like to thank Prof. C. J. Zhu for useful discussions. The work was
supported by the NSF of China under Grant No. 10525522 and 10875129.

## References

  * (1) J. Stachel, Early History of Quantum Gravity. In Black Holes, Gravitational Radiation and the Universe, ed by B.R. Iyer, B. Bhawal eds. (Kluwer Academic Publisher, Netherlands, 1999).
  * (2) M. B. Green, J. H. Schwarz and E. Witten, Superstring theory (Cambridge Univ. Press, Cambridge, 1987); J. Polchinski, Superstring Theory and Beyond, (Cambridge University Press, Cambridge, 1998).
  * (3) C. Rovelli, arXiv:gr-qc/9710008, Living Rev. Rel. 1, 1 (1998); A. Ashtekar, arXiv:gr-qc/0112038; L. Smolin, arXiv:hep-th/0303185; T. Thiemann, arXiv:gr-qc/0210094, Lect. Notes Phys. 631, 41 (2003).
  * (4) G. Amelino-Camelia, arXiv:gr-qc/0412136.
  * (5) L. J. Garay, Int. J. Mod. Phys. A 10, 145 (1995).
  * (6) J. R. Ellis, J. L. Lopez, N. E. Mavromatos and D. V. Nanopoulos, Phys. Rev. D 53, 3846 (1996).
  * (7) V. A. Kostelecky, CPT and Lorentz symmetry . Proceedings of the 3rd Meeting on CPT and Lorentz Symmetry (CPT 04), Bloomington, Indiana, 4-7 Aug 2004. Published in Hackensack, USA: World Scientific (2005).
  * (8) G. Amelino-Camelia, Phys. Rev. D 62, 024015 (2000).
  * (9) C. D. Hoyle, et al., Phys. Rev. D 70, 042004 (2004); Phys. Rev. Lett. 86, 1418 (2001).
  * (10) S. Dimopoulos and G. Landsberg, Phys. Rev. Lett. 87, 161602 (2001); L. A. Anchordoqui, J. L. Feng, H. Goldberg and A. D. Shapere, Phys. Rev. D 65, 124027 (2002).
  * (11) D. Mattingly, Living Rev. Rel. 8, 5 (2005).
  * (12) A. Connes: Noncommutative Geometry, (Academic Press, 1995).
  * (13) S. Majid: Foundations of Quantum Group Theory, (Cambridge University Press, 1995).
  * (14) S. Majid and H. Ruegg: Phys. Lett. B 334, 348 (1994).
  * (15) G. Amelino-Camelia, Phys. Lett. B510, 255 (2001).
  * (16) G. Amelino-Camelia, Int. J. Mod. Phys. D11, 35 (2002).
  * (17) G. Amelino-Camelia, Nature 418, 34 (2002).
  * (18) J. Magueijo and L. Smolin, Phys. Rev. Lett. 88, 190403 (2002).
  * (19) J. Magueijo and L. Smolin, Phys. Rev. D67, 044017 (2003).
  * (20) T. Jacobson, S. Liberati and D. Mattingly, Annals Phys. 321, 150 (2006).
  * (21) F. Girelli, S. Liberati and L. Sindoni, Phys. Rev. D75, 064015 (2007).
  * (22) G. Randers, Phys. Rev. 59, 195 (1941).
  * (23) Z. Chang and X. Li, Phys. Lett. B 663, 103 (2008).
  * (24) Z. Chang and X. Li, Chinese Physics C 33, 626 (2009).
  * (25) X. Li, Z. Chang and X. H. Mo, arXiv:hep-th/1001.2667.
  * (26) D. Bao, S. S. Chern and Z. Shen, An Introduction to Riemann–Finsler Geometry, Graduate Texts in Mathematics 200, Springer, New York, 2000.
  * (27) M. Matsumoto, Foundations of Finsler Geometry and Special Finsler Spaces, Kaiseisha Press, Saikawa Shigaken, Japan 1986.
  * (28) P. L. Antonelli and S. F. Rutz, “Finsler Geometry” Advanced studies in Pure Mathematics 48, Sapporo (2005) p. 210 -In memory of M.Matsumoto.
  * (29) Z. Szabo, Ann. Glob. Anal. Geom 34, 381 (2008).
  * (30) S. Weinberg, Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity, Wiley, New York, 1972.
  * (31) D. Clowe, S. W. Randall and M. Markevitch, http://flamingos.astro.ufl.edu/1e0657/index.html (2006); D. Clowe, S. W. Randall and M. Markevitch, arXiv:astro-ph/0611496.
  * (32) E. J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753 (2006).


# Asynchronous parareal time discretization for partial differential equations

Frédéric Magoulès CentraleSupélec, Université Paris-Saclay, France
(correspondence, frederic.magoules@hotmail.com). Guillaume Gbikpi-Benissan IRT
SystemX, France.

###### Abstract

Asynchronous iterations are more and more investigated for both scaling and
fault-resilience purpose on high performance computing platforms. While so
far, they have been exclusively applied within space domain decomposition
frameworks, this paper advocates a novel application direction targeting time-
decomposed time-parallel approaches. Specifically, an asynchronous iterative
model is derived from the Parareal scheme, for which convergence and speedup
analysis are then conducted. It turned out that Parareal and async-Parareal
feature very close convergence conditions, asymptotically equivalent,
including the finite-time termination property. Based on a computational cost
model aware of unsteady communication delays, our speedup analysis shows the
potential performance gain from asynchronous iterations, which is confirmed by
some experimental case of heat evolution on a homogeneous supercomputer. This
primary work clearly suggests possible further benefits from asynchronous
iterations.

Keywords: time domain decomposition; Parareal; asynchronous iterations; time-
parallel algorithms; time-dependent problems

##  1 Introduction

While high-performance computing (HPC) architectures are more and more
parallel, the development of efficient and robust algorithms that can fully
exploit high degrees of parallelism is still an actual challenge in many
scientific computing fields. For solving partial differential equations (PDEs)
defined on the interior of arbitrary geometrical shapes, strong advances
resulted in highly scalable domain decomposition methods (DDMs) (see, e.g.,
[24, 38]), which feature high rates of speedup and efficiency. Further,
although it first appears unnatural to break the time line causality in time-
dependent differential equations, thorough applications of DDMs still led
researchers to time-parallel methods for allowing actual parallel time-
integration iterations. Waveform relaxation methods [21], for instance, solve
space-time PDEs by means of independent time-integration processes performed
on subregions of the spatial domain, while ensuring global consistency through
time-dependent interface conditions. Multiple shooting [9] and time-multigrid
[20] methods led to the more recent Parareal time discretization scheme [23],
and a general analysis of parallel implicit time-integration algorithms [14],
where the time domain itself is decomposed. Basically, two levels of time
discretization grid are considered. An approximated solution of the PDE on a
coarse grid provides initial values on each time subinterval, when more
accurate solutions are thus independently computed on the basis of a
sufficiently fine time grid. These fine solutions are then used to enhance new
another coarse solution. In this sense, the Parareal computational model can
also be considered as an iterative predictor-corrector scheme. We refer to
[18] for an introductory insight into time-decomposed time-integration
approaches.

Despite such great advances with DDMs, there remain scaling limits upon these
parallel algorithms, due to the fact that they require, at some point,
sequential global iterations, even if each of these iterations is performed in
a parallel fashion. This sequence requirement implies iterative global
synchronization, which causes both bottleneck and fault-resilience issues at
HPC scales. For such methods then, the only way to completely break sequential
computation is to resort to asynchronous iterations [7], which therefore
constitute today one of the most interesting algorithmic approaches toward HPC
applications. Arisen from linear relaxation schemes (see [12]), asynchronous
iterations are basically described by a non-deterministic computational model,
such that there is no more a single sequence of parallel iterations for
approximating a global solution. Instead, one can independently run a proper
sequence of iterations on each subregion of a domain, therefore communication
can be completely overlapped by meaningful computation. So far, only fixed-
point iterative models have been proved to be asynchronously convergent, under
various general contraction conditions (see, e.g., [31, 4, 13, 5]). Still,
despite the lower convergence rate of asynchronous iterations, their practical
efficiency in execution time over classical fixed-point methods is undoubtedly
established, thanks to extensive studies addressing various computational
problems (see, e.g., [22, 19, 8, 10, 11]). The asynchronous iterations theory
has been largely applied within the framework of overlapping Schwarz domain
decomposition methods (see, e.g., [15, 36, 37]), according to their tight
relation with block-Jacobi and Gauss-Seidel relaxation schemes (see, e.g.,
[17]). This included waveform relaxation methods as well (see, e.g., [2, 16]),
which constitutes so far the only application targeting time-parallel methods.
Very recently, non-overlapping Schwarz methods have been investigated (see
[29]), just as some spatial domain substructuring approach (see [30]). The
goal of this paper is to introduce a novel asynchronous iterations application
direction by addressing time-decomposed time-parallel algorithms through the
well-known Parareal iterative scheme. Besides a successful experimental
investigation in [28], where we only focused on implementation aspects, we
present here a theoretical analysis to fully assess both convergence
conditions and speedup potential.

Section 2 provides a short general introduction to asynchronous iterations,
with most important analysis tools for assessing their application to the
Parareal iterative scheme. In section 3 then, we deal with its asynchronous
convergence analysis, after deriving the corresponding iterative model.
Speedup analysis is conducted in section 4, where we first propose a
computational cost model not neglecting communication effects. Section 5
finally presents practical performance experimented on a supercomputer, for a
heat evolution problem, and our conclusions follow in section 6.

##  2 General asynchronous iterations

###  2.1 Computational model

Asynchronous iterations arose in the framework of linear relaxation methods,
from the successive experimental and theoretical works of Rosenfeld [33] and
Chazan & Miranker [12], respectively. Given a linear system

$Ax=b,$

where $x$ is a vector of $n$ unknowns, let $M$ be a nonsingular matrix and
consider the mapping $f$ defined by

$f(x)=(I-M^{-1}A)x+M^{-1}b,$

with $I$ denoting the identity matrix. Given then that

$A=M-(M-A),$

this linear system can be seen as a fixed-point problem, as we have

$Ax=b\quad\iff\quad x=f(x).$

For instance, the Jacobi relaxation scheme consists of taking $M$ as the
diagonal part of $A$. A classical iterative fixed point search generates a
sequence $\\{x^{k}\\}_{k\in\mathbb{N}}$ of vectors such that

$x^{k+1}=f(x^{k}),$ (1)

and such that, for any given initial vector $x^{0}$,

$\lim_{k\to+\infty}x^{k}=x^{*},$

with $x^{*}$ being the unique fixed point of $f$. In a parallel computation
context with $p$ processors, $p\leq n$, one may assume a decomposition of the
form

$x=(x_{1},\ldots,x_{p})^{\mathsf{T}},\qquad
f(x)=(f_{1}(x),\ldots,f_{p}(x))^{\mathsf{T}},$

which gives an equivalent sequence of parallel iterations

$x_{i}^{k+1}=f_{i}(x_{1}^{k},\ldots,x_{p}^{k}),\qquad\forall
i\in\\{1,\ldots,p\\}.$

The first step toward asynchronous iterations is the introduction of the free
steering (see, e.g., [34]), where, at each iteration, one is allowed to relax
an arbitrary subset of unknowns. Note however that parallel computing were not
yet targeted at this time, and that, except for free steering, the choice of
the unknown(s) to be relaxed was based on deterministic criteria (like Gauss-
Seidel relaxation, for instance, which follows a cyclic pattern). Let us then
consider, at each iteration, an arbitrary subset

$P^{(k)}\subseteq\\{1,\ldots,p\\},$

defining the relaxed subvectors of unknowns. The corresponding free-steering
iterative model is given by

$x_{i}^{k+1}=\left\\{\begin{array}[]{ll}f_{i}\left(x_{1}^{k},\ldots,x_{p}^{k}\right),&\forall
i\in P^{(k)},\\\ x_{i}^{k},&\forall i\notin P^{(k)},\end{array}\right.$

where iterations

$x_{i}^{k+1}=x_{i}^{k},\quad\forall i\notin P^{(k)},$

are only implicit. Then, while considering actual parallel free steering, one
would be interested in letting each subvector be continuously relaxed at a
random pace. In this case therefore, there is no more explicit global
iterations, and the sequence $\\{P^{(k)}\\}_{k\in\mathbb{N}}$ instead models
implicit global iterations by accounting each concurrent update of at least
one subvector. In such a parallel context, Chazan & Miranker [12] took into
account the delay of accessing a subvector, in the sense that this subvector
could be relaxed in the meantime, and therefore could become outdated before
being used. Consequently, a subvector $x_{i}^{k+1}$, $i\in\\{1,\ldots,p\\}$,
is no more necessarily computed on the basis of subvectors $x_{j}^{k}$,
$j\in\\{1,\ldots,p\\}$, but on the basis of possibly outdated subvectors
$x_{j}^{\tau_{j}^{i}(k)}$, i.e. with

$\tau_{j}^{i}(k)\leq k.$

This yielded asynchronous iterations

$x_{i}^{k+1}=\left\\{\begin{array}[]{ll}f_{i}\left(x_{1}^{\tau_{1}^{i}(k)},\ldots,x_{p}^{\tau_{p}^{i}(k)}\right),&\forall
i\in P^{(k)},\\\ x_{i}^{k},&\forall i\notin P^{(k)}.\end{array}\right.$ (2)

Two straight assumptions are made to allow for convergence. First, no
component $i\in\\{1,\ldots,p\\}$ should definitively stop being updated, which
implies that $i$ belongs to successive subsets $P^{(k)}$ infinitely often,
i.e.

$\forall i\in\\{1,\dots,p\\},\quad\operatorname{card}\\{k\in\mathbb{N}\ |\
i\in P^{(k)}\\}=+\infty,$ (3)

where $\operatorname{card}$ denotes the cardinality. Secondly, delays to
access updated components are bounded, so that

$\forall
i,j\in\\{1,\dots,p\\},\quad\lim_{k\to+\infty}\tau_{j}^{i}(k)=+\infty.$ (4)

Note finally that classical iterations are thus specialization of asynchronous
iterations where

$P^{(k)}=\\{1,\ldots,p\\},\quad\tau_{j}^{i}(k)=k,\qquad\forall
k\in\mathbb{N}.$

More generally now, we can consider an arbitrary, non necessarily linear,
mapping

$\widetilde{f}:E^{m}\to E,\qquad m\in\mathbb{N},$

with $E^{m}$ being the Cartesian product of $m$ sets
$E=E_{1}\times\cdots\times E_{p}$, and a fixed-point problem

$x=\widetilde{f}(x,x,\ldots,x).$

The study of such a generalization goes back to Baudet [4] for addressing
Newton asynchronous iterations. This generalizes the model (2) to

$x_{i}^{k+1}=\left\\{\begin{array}[]{ll}\widetilde{f}_{i}\left(x_{1}^{\tau_{1,1}^{i}(k)},\ldots,x_{p}^{\tau_{p,1}^{i}(k)},\ldots,x_{1}^{\tau_{1,m}^{i}(k)},\ldots,x_{p}^{\tau_{p,m}^{i}(k)}\right),&\forall
i\in P^{(k)},\\\ x_{i}^{k},&\forall i\notin P^{(k)},\end{array}\right.$ (5)

in the sense that (2) corresponds to the case $m=1$.

###  2.2 Convergence conditions

It is well known (see, e.g., Theorem 2.1 in [3]) that the sequential iterative
model (1) is convergent if and only if

$\rho(I-M^{-1}A)<1,$

where $\rho(.)$ denotes the spectral radius. Chazan & Miranker [12] showed
that the asynchronous iterative model (2) is convergent if and only if

$\rho(|I-M^{-1}A|)<1,$

where $|.|$ denotes the absolute value, which is slightly more restrictive,
since we always have (see, e.g., Corollary 6.3 in [6])

$\rho(I-M^{-1}A)\leq\rho(|I-M^{-1}A|).$

Still, to study the application of asynchronous iterations to the Parareal
scheme, we would rather need the framework of arbitrary mappings, for which we
recall here two most general sufficient contraction conditions guaranteeing
asynchronous convergence.

Let each set $E_{i}$, with $i\in\\{1,\ldots,p\\}$, be a Banach space provided
with some norm $\|.\|_{(i)}$, and consider weighted maximum norms

$\|x\|_{\infty}^{w}=\max_{i\in\\{1,\ldots,p\\}}\frac{\|x_{i}\|_{(i)}}{w_{i}},$
(6)

where $w$ is a positive vector with real entries $w_{i}$,
$i\in\\{1,\ldots,p\\}$. The following result goes back to El Tarazi [13]:

######  Theorem 1 (El Tarazi, 1982).

The asynchronous iterative model (5) is convergent if there exists a positive
vector $w$ and a positive real $\alpha<1$ such that

$\forall X,Y\in
E^{m},\quad\|\widetilde{f}(X)-\widetilde{f}(Y)\|_{\infty}^{w}\leq\alpha\max\\{\|x^{(1)}-y^{(1)}\|_{\infty}^{w},\ldots,\|x^{(m)}-y^{(m)}\|_{\infty}^{w}\\},$

with

$X=(x^{(1)},\ldots,x^{(m)}),\quad Y=(y^{(1)},\ldots,y^{(m)}).$

Bertsekas [5] then provided a more general contraction framework based on the
identification of an implicit sequence of nested sets which converges toward
the solution set $\\{x^{*}\\}$ (actually, this even applies for mappings
admitting several fixed points):

######  Theorem 2 (Bertsekas, 1983).

The asynchronous iterative model (5), with

$m=1,$

is convergent if there exists a sequence $\\{S^{(l)}\\}_{l\in\mathbb{N}}$ of
subsets of $E$, with

$x^{0}\in S^{(0)},\qquad S^{(l)}=S^{(l)}_{1}\times\cdots\times
S^{(l)}_{p},\qquad S^{(l+1)}\subset
S^{(l)},\qquad\lim_{l\to+\infty}S^{(l)}=\\{x^{*}\\},$

and such that

$\forall x\in S^{(l)},\quad\widetilde{f}(x)\in S^{(l+1)}.$

As mentioned above, the sets sequence $\\{S^{(l)}\\}_{l\in\mathbb{N}}$ is only
implicitly generated all along the successive updates of the components
$x_{i}$, and, for more clarity, we do use another sequence variable, $l$,
since several asynchronous iterations $k$ could be performed within a same set
$S^{(l)}$ before entering the next one, $S^{(l+1)}$. The main difficulty in
applying such a result to a specific mapping $\widetilde{f}$ is to be able to
identify the sets $S^{(l)}$ and $S^{(l)}_{i}$. Bertsekas & Tsitsiklis [7]
showed that this was more general than the contraction framework based on
weighted maximum norms. When, for instance, Theorem 1 is applicable, with
$m=1$, then we have (see [6, section 6.3] and [7])

$S^{(l)}=\left\\{x\in E\ |\
\|x-x^{*}\|_{\infty}^{w}\leq\alpha^{l}\|x^{0}-x^{*}\|_{\infty}^{w}\right\\},\qquad
l\in\mathbb{N},$ (7)

and

$S^{(l)}_{i}=\left\\{\dot{x}\in E_{i}\ |\
\frac{\|\dot{x}-x_{i}^{*}\|_{(i)}}{w_{i}}\leq\alpha^{l}\|x^{0}-x^{*}\|_{\infty}^{w}\right\\},\qquad
l\in\mathbb{N},\quad i\in\left\\{1,\ldots,p\right\\}.$

##  3 Parareal asynchronous iterations

###  3.1 Computational model

Let us consider a time-dependent PDE,

$\frac{\partial u}{\partial t}(s,t)=g(u,s,t),\qquad
t\in[0,T]\subset\mathbb{R},\quad s\in\Omega\subset\mathbb{R}^{3},$ (8)

with some given boundary condition,

$h(u,s,t)=0,\qquad t\in[0,T],\quad s\in\partial\Omega,$

and initial values,

$u^{(0)}(s):=u(s,0),\qquad s\in\Omega.$

Let the time domain $[0,T]$ be given by

$[0,T]=[T_{0},T_{1}]\cup[T_{1},T_{2}]\cup\cdots\cup[T_{p-1},T_{p}],\qquad
p\in\mathbb{N},$

and assume initial values $\lambda_{i-1}(s)$, with $s\in\Omega$, on each
subdomain $[T_{i-1},T_{i}]$, with $i\in\\{1,\ldots,p\\}$. Independent
subproblems can then be formulated by

$\left\\{\begin{array}[]{lllll}\dfrac{\partial u_{i}}{\partial
t}(s,t)&=&g(u_{i},s,t),&t\in[T_{i-1},T_{i}],&s\in\Omega,\\\
h(u_{i},s,t)&=&0,&t\in[T_{i-1},T_{i}],&s\in\partial\Omega,\\\
u_{i}(s,T_{i-1})&=&\lambda_{i-1}(s),&s\in\Omega,&\end{array}\right.$ (9)

which is equivalent to the global problem (8) if and only if

$\left\\{\begin{array}[]{llll}\lambda_{0}(s)&=&u^{(0)}(s),&\\\
\lambda_{i-1}(s)&=&u_{i-1}(s,T_{i-1}),&\forall
i\in\\{2,\ldots,p\\}.\end{array}\right.$ (10)

Now let, for $i\in\\{1,\ldots,p\\}$,

$u_{i}^{(t)}(s):=u_{i}(s,t),\qquad s\in\Omega,\quad t\in[T_{i-1},T_{i}],$

and let $F$ be a function which, to initial values $\lambda_{i-1}(s)$,
associates a sufficiently accurate solution of the subproblem (9), at time
$T_{i}$, i.e.,

$F(\lambda_{i-1})\cong u_{i}^{(T_{i})}.$

Similarly, let $G$ be a function which, to initial values $\lambda_{i-1}(s)$,
associates a cheap approximation of the solution of the subproblem (9), at
time $T_{i}$, i.e.,

$G(\lambda_{i-1})\approx u_{i}^{(T_{i})}.$

From initial guess

$\left\\{\begin{array}[]{lcll}\lambda_{0}^{0}&=&u^{(0)},&\\\
\lambda_{i}^{0}&=&G(\lambda_{i-1}^{0}),&\forall
i\in\\{1,\ldots,p\\},\end{array}\right.$

the Parareal iterative scheme [23, 14] defines sequences
$\\{\lambda_{i}^{k}\\}_{k\in\mathbb{N}}$ such that

$\left\\{\begin{array}[]{lcll}\lambda_{0}^{k+1}&=&\lambda_{0}^{k},&\\\
\lambda_{i}^{k+1}&=&G(\lambda_{i-1}^{k+1})+F(\lambda_{i-1}^{k})-G(\lambda_{i-1}^{k}),&\forall
i\in\\{1,\ldots,p\\},\end{array}\right.$ (11)

which should converge to a collection $\\{\lambda_{i}^{*}\\}$ satisfying the
interface conditions (10). This algorithm is thus expected to provide the
solution $u^{(T)}$ of (8) at time $T$, without applying the whole sequence

$u^{(T)}=F\circ F\circ\cdots\circ F(u^{(0)}).$

To easily assess such a computational model within the previous framework of
general asynchronous iterations, let us consider linear operators

$G(\lambda_{i})=\mathcal{G}\lambda_{i}+\delta,\quad
F(\lambda_{i})=\mathcal{F}\lambda_{i}+\zeta,\qquad i\in\\{0,1,\ldots,p\\},$

and a vector

$\lambda=(\lambda_{0},\lambda_{1},\ldots,\lambda_{p}).$

Then, (11) simplifies to

$\left\\{\begin{array}[]{lcll}\lambda_{0}^{k+1}&=&\lambda_{0}^{k},&\\\
\lambda_{i}^{k+1}&=&\mathcal{G}\lambda_{i-1}^{k+1}+(\mathcal{F}-\mathcal{G})\lambda_{i-1}^{k}+\zeta,&\forall
i\in\\{1,\ldots,p\\},\end{array}\right.$

which corresponds to synchronous relaxations (1) of a problem

$A\lambda=b,$

with

$A=\begin{bmatrix}I&0&0&\cdots&0\\\ -\mathcal{F}&I&0&\cdots&0\\\
0&-\mathcal{F}&I&\ddots&\vdots\\\ \vdots&\ddots&\ddots&\ddots&0\\\
0&\cdots&0&-\mathcal{F}&I\end{bmatrix},\quad
M=\begin{bmatrix}I&0&0&\cdots&0\\\ -\mathcal{G}&I&0&\cdots&0\\\
0&-\mathcal{G}&I&\ddots&\vdots\\\ \vdots&\ddots&\ddots&\ddots&0\\\
0&\cdots&0&-\mathcal{G}&I\end{bmatrix},\quad
b=\begin{bmatrix}\lambda_{0}^{*}\\\ \zeta\\\ \zeta\\\ \vdots\\\
\zeta\end{bmatrix},$

recalling that the iterations mapping $f$ is given by

$f(\lambda)=(I-M^{-1}A)\lambda+M^{-1}b.$

######  Remark.

This corresponds to a preconditioned Richardson iterative scheme, as also
formulated, e.g., in [32].

Nonetheless, a straightforward application of asynchronous relaxations (2)
would yield a computation of the form

$\lambda_{i}^{k+1}=\left\\{\begin{array}[]{ll}G(\lambda_{i-1}^{k+1})+F(\lambda_{i-1}^{\tau_{i-1}^{i}(k)})-G(\lambda_{i-1}^{\tau_{i-1}^{i}(k)}),&\forall
i\in P^{(k)},\\\ \lambda_{i}^{k},&\forall i\notin P^{(k)},\end{array}\right.$

which is not of great interest, due to the synchronization on
$G(\lambda_{i-1}^{k+1})$. We would therefore rather be interested in seeing
what happens with a computation, at least, of the form

$\lambda_{i}^{k+1}=\left\\{\begin{array}[]{ll}G(\lambda_{i-1}^{\tau_{i-1,1}^{i}(k)})+F(\lambda_{i-1}^{\tau_{i-1,2}^{i}(k)})-G(\lambda_{i-1}^{\tau_{i-1,2}^{i}(k)}),&\forall
i\in P^{(k)},\\\ \lambda_{i}^{k},&\forall i\notin P^{(k)},\end{array}\right.$
(12)

which corresponds to the more general asynchronous iterations (5), with $m=2$,
and where the iterations mapping $\widetilde{f}$ is given by

$\left\\{\begin{array}[]{lcll}\widetilde{f}_{0}(\lambda^{(1)},\lambda^{(2)})&=&\lambda_{0}^{*},&\\\
\widetilde{f}_{i}(\lambda^{(1)},\lambda^{(2)})&=&G(\lambda_{i-1}^{(1)})+F(\lambda_{i-1}^{(2)})-G(\lambda_{i-1}^{(2)}),&\forall
i\in\\{1,\ldots,p\\}.\end{array}\right.$ (13)

It is easy to verify that

$\lambda=f(\lambda)\quad\iff\quad\lambda=\widetilde{f}(\lambda,\lambda),$

which ensures an equivalent fixed-point problem. We recall that the
asynchronous iterations variable $k$ grows faster than the synchronous
iterations one, therefore, the delays expressed in (12) do not have to be seen
as implying a retarded iterative procedure.

Note, finally, that the mapping $\widetilde{f}$ is explicitly expressed
without linearity assumption on $F$ and $G$. In the following convergence
analysis of the Parareal asynchronous iterations (12), such an assumption will
be sometimes made only to keep the development simple and to use same
notations as with the linear synchronous formulation. Furthermore, to fully
exploit the domain decomposition point of view, one will notice that, in the
asynchronous framework, the extension of these results to a general
formulation using $F_{i}$ and $G_{i}$ is quite straightforward.

###  3.2 Convergence conditions

Let $\|.\|_{\infty}$ denote a maximum norm (6) where $w$ is a vector of $1$.
We first restate preliminary results in our linear synchronous framework.

######  Theorem 3.

The Parareal synchronous iterative model (11) defines a sequence
$\\{\lambda^{k}\\}_{k\in\mathbb{N}}$ such that

$\|\lambda^{k}-\lambda^{*}\|_{\infty}\leq\alpha^{k}\|\lambda^{0}-\lambda^{*}\|_{\infty},$

with

$\alpha=\frac{1-\theta^{p}}{1-\theta}\|\mathcal{F}-\mathcal{G}\|,\qquad\theta\neq
1,\quad\theta\geq\|\mathcal{G}\|,$

and $\|.\|$ denoting any given norm.

###### Proof.

See Theorem 4.5 in [18] for a possible proof, or directly consider the fact
that linear fixed-point formulation implies

$\|\lambda^{k}-\lambda^{*}\|_{\infty}\leq(\|I-M^{-1}A\|_{\infty})^{k}\|\lambda^{0}-\lambda^{*}\|_{\infty}.$

∎

######  Remark.

$\lim_{p\to+\infty}\frac{1-\theta^{p}}{1-\theta}\|\mathcal{F}-\mathcal{G}\|=\left\\{\begin{array}[]{ll}+\infty,&\theta>1,\\\
\dfrac{1}{1-\theta}\|\mathcal{F}-\mathcal{G}\|,&\theta<1.\end{array}\right.$

One is therefore interested by the case $\theta<1$ which ensures a stable
iterative procedure.

######  Corollary 1.

The Parareal synchronous iterative model (11) is convergent if

$\|\mathcal{G}\|<1,$

and

$\|\mathcal{G}\|+\|\mathcal{F}-\mathcal{G}\|<1+\|\mathcal{G}\|^{p}\|\mathcal{F}-\mathcal{G}\|.$

###### Proof.

It directly follows from Theorem 3 that convergence is guaranteed if

$\frac{1-\theta^{p}}{1-\theta}\|\mathcal{F}-\mathcal{G}\|<1.$

If $\theta<1$, this is implied by

$(1-\theta^{p})\|\mathcal{F}-\mathcal{G}\|<1-\theta.$

If then $\|\mathcal{G}\|<1$, we can take $\theta=\|\mathcal{G}\|$, which leads
to the result. ∎

######  Remark.

The Parareal synchronous iterative model (11) is thus asymptotically
convergent ($p\to+\infty$) if

$\|\mathcal{G}\|+\|\mathcal{F}-\mathcal{G}\|<1.$

We need here some words about the terminology. It is well established (see
[14]) that Parareal iterations terminate in finite time $k=p$, due to the
sequential propagation of the accurate solution, from $[T_{0},T_{1}]$ to
$[T_{p-1},T_{p}]$. By “convergence”, we therefore mean the fact that the
iterative procedure actually allows one to reach a given precision by
decreasing the maximum error over the whole time domain $[T_{0},T_{p}]$, which
thus suggests the possibility to converge before getting sequential. To be
rigorous though, “convergence”, here, is a shortcut word for “convergence with
respect to the maximum norm”.

In the following, we would now like to extend such results to the asynchronous
framework.

######  Theorem 4.

The Parareal asynchronous iterative model (12) defines a sequence
$\\{\lambda^{k}\\}_{k\in\mathbb{N}}$ such that

$\|\lambda^{k}-\lambda^{*}\|_{\infty}\leq\widetilde{\alpha}^{\sigma(k)}\|\lambda^{0}-\lambda^{*}\|_{\infty},$

with

$\widetilde{\alpha}=\|\mathcal{G}\|+\|\mathcal{F}-\mathcal{G}\|,$

where $\sigma$ is an integer-valued function on $\mathbb{N}$ satisfying

$\underset{k\to+\infty}{\lim}\sigma(k)=+\infty.$

###### Proof.

According to (12), we have, $\forall i\in P^{(k)}$,

$\displaystyle\lambda_{i}^{k+1}$
$\displaystyle=\mathcal{G}\lambda_{i-1}^{\tau_{i-1,1}^{i}(k)}+(\mathcal{F}-\mathcal{G})\lambda_{i-1}^{\tau_{i-1,2}^{i}(k)}+\zeta,$
$\displaystyle\lambda_{i}^{k+1}-\lambda_{i}^{*}$
$\displaystyle=\mathcal{G}(\lambda_{i-1}^{\tau_{i-1,1}^{i}(k)}-\lambda_{i-1}^{*})+(\mathcal{F}-\mathcal{G})(\lambda_{i-1}^{\tau_{i-1,2}^{i}(k)}-\lambda_{i-1}^{*}),$
$\displaystyle\|\lambda_{i}^{k+1}-\lambda_{i}^{*}\|$
$\displaystyle\leq\|\mathcal{G}\|\|\lambda_{i-1}^{\tau_{i-1,1}^{i}(k)}-\lambda_{i-1}^{*}\|+\|\mathcal{F}-\mathcal{G}\|\|\lambda_{i-1}^{\tau_{i-1,2}^{i}(k)}-\lambda_{i-1}^{*}\|,$
$\displaystyle\|\lambda_{i}^{k+1}-\lambda_{i}^{*}\|$
$\displaystyle\leq(\|\mathcal{G}\|+\|\mathcal{F}-\mathcal{G}\|)\max\\{\|\lambda_{i-1}^{\tau_{i-1,1}^{i}(k)}-\lambda_{i-1}^{*}\|,\|\lambda_{i-1}^{\tau_{i-1,2}^{i}(k)}-\lambda_{i-1}^{*}\|\\},$

which gives

$\left\\{\begin{array}[]{lcll}\|\lambda_{i}^{k+1}-\lambda_{i}^{*}\|&\leq&\widetilde{\alpha}\max\\{\|\lambda_{i-1}^{\tau_{i-1,1}^{i}(k)}-\lambda_{i-1}^{*}\|,\|\lambda_{i-1}^{\tau_{i-1,2}^{i}(k)}-\lambda_{i-1}^{*}\|\\},&\forall
i\in P^{(k)},\\\
\|\lambda_{i}^{k+1}-\lambda_{i}^{*}\|&=&\|\lambda_{i}^{k}-\lambda_{i}^{*}\|,&\forall
i\notin P^{(k)}.\end{array}\right.$

It follows that, $\forall i\in\\{1,\ldots,p\\}$,

$\|\lambda_{i}^{k+1}-\lambda_{i}^{*}\|\leq\max\left\\{\begin{array}[]{l}\widetilde{\alpha}\max\\{\|\lambda_{i-1}^{\tau_{i-1,1}^{i}(k)}-\lambda_{i-1}^{*}\|,\|\lambda_{i-1}^{\tau_{i-1,2}^{i}(k)}-\lambda_{i-1}^{*}\|\\},\\\
\|\lambda_{i}^{k}-\lambda_{i}^{*}\|\end{array}\right\\},$

and then,

$\|\lambda^{k+1}-\lambda^{*}\|_{\infty}\leq\max_{i\in\\{1,\ldots,p\\}}\left\\{\begin{array}[]{l}\widetilde{\alpha}\|\lambda_{i-1}^{\tau_{i-1,1}^{i}(k)}-\lambda_{i-1}^{*}\|,\\\
\widetilde{\alpha}\|\lambda_{i-1}^{\tau_{i-1,2}^{i}(k)}-\lambda_{i-1}^{*}\|,\\\
\|\lambda_{i}^{k}-\lambda_{i}^{*}\|\end{array}\right\\}.$ (14)

Let us define sets

$S^{(k)}=\\{\lambda\ |\
\|\lambda-\lambda^{*}\|_{\infty}\leq\widetilde{\alpha}^{\sigma(k)}\|\lambda^{0}-\lambda^{*}\|_{\infty}\\},$

with $\sigma$ being a function on $\mathbb{N}$ such that

$\sigma(0)=0.$

We obviously verify

$\lambda^{0}\in S^{(0)}.$

Assume then that

$\forall l\leq k,\quad\lambda^{l}\in S^{(l)}.$

Therefore, from (14), it follows that

$\|\lambda^{k+1}-\lambda^{*}\|_{\infty}\leq\max_{i\in\\{1,\ldots,p\\}}\left\\{\begin{array}[]{l}\widetilde{\alpha}\widetilde{\alpha}^{\sigma(\tau_{i-1,1}^{i}(k))}\|\lambda^{0}-\lambda^{*}\|_{\infty},\\\
\widetilde{\alpha}\widetilde{\alpha}^{\sigma(\tau_{i-1,2}^{i}(k))}\|\lambda^{0}-\lambda^{*}\|_{\infty},\\\
\widetilde{\alpha}^{\sigma(k)}\|\lambda^{0}-\lambda^{*}\|_{\infty}\end{array}\right\\},$

and by taking $\sigma(k+1)$ such that

$\widetilde{\alpha}^{\sigma(k+1)}=\max_{i\in\\{1,\ldots,p\\}}\left\\{\widetilde{\alpha}^{\sigma(\tau_{i-1,1}^{i}(k))+1},\widetilde{\alpha}^{\sigma(\tau_{i-1,2}^{i}(k))+1},\widetilde{\alpha}^{\sigma(k)}\right\\},$
(15)

we do have

$\|\lambda^{k+1}-\lambda^{*}\|_{\infty}\leq\widetilde{\alpha}^{\sigma(k+1)}\|\lambda^{0}-\lambda^{*}\|_{\infty},$

which confirms that

$\forall k\in\mathbb{N},\quad\lambda^{k}\in S^{(k)}.$

Finally, according to the classical assumption (3) on $P^{k}$, we eventually
have, $\forall i\notin P^{(k)}$,

$\lambda_{i}^{k+1}=\lambda_{i}^{k}=\lambda_{i}^{k-1}=\cdots=\lambda_{i}^{k_{i}+1}=\mathcal{G}\lambda_{i-1}^{\tau_{i-1,1}^{i}(k_{i})}+(\mathcal{F}-\mathcal{G})\lambda_{i-1}^{\tau_{i-1,2}^{i}(k_{i})}+\zeta,\quad
i\in P^{(k_{i})},$

which implies in (15) that

$\sigma(k+1)\in\left\\{\begin{array}[]{l}\sigma(\tau_{i-1,1}^{i}(k))+1,\ \
\sigma(\tau_{i-1,2}^{i}(k))+1,\\\ \sigma(\tau_{i-1,1}^{i}(k_{i}))+1,\ \
\sigma(\tau_{i-1,2}^{i}(k_{i}))+1\end{array}\right\\}_{i\in\\{1,\ldots,p\\}}.$

Therefore, with the other classical assumption (4) on $\tau_{i-1,1}^{i}$ and
$\tau_{i-1,2}^{i}$, we also eventually have

$\sigma(k+1)=\sigma(k-d_{1})+1=\sigma(k-d_{2})+1+1=\cdots,$

with

$0\leq d_{1}<d_{2}<\cdots\leq k,$

which ensures that

$\underset{k\to+\infty}{\lim}\sigma(k)=+\infty,$

and concludes the proof. ∎

With such an error bound, a convergence condition trivially follows.

######  Corollary 2.

The Parareal asynchronous iterative model (12) is convergent if

$\|\mathcal{G}\|+\|\mathcal{F}-\mathcal{G}\|<1.$ (16)

###### Proof.

This is straightforward, from Theorem 4. ∎

######  Remark.

The proof of Theorem 4 exhibits

$\|\lambda_{i}^{k+1}-\lambda_{i}^{*}\|\leq(\|\mathcal{G}\|+\|\mathcal{F}-\mathcal{G}\|)\max\\{\|\lambda_{i-1}^{\tau_{i-1,1}^{i}(k)}-\lambda_{i-1}^{*}\|,\|\lambda_{i-1}^{\tau_{i-1,2}^{i}(k)}-\lambda_{i-1}^{*}\|\\},$

which clearly indicates that Parareal asynchronous iterations (12) match the
contraction framework of Theorem 1, which also implies convergence under the
same sufficient condition (16).

One of the advantages of asynchronous iterations is that the algorithm can be
studied on the basis of only local behavior $\widetilde{f}_{i}$, independently
of its input which, therefore, can be provided by whatever other local
contracting mapping $\widetilde{f}_{i-1}$. This is an interesting feature
which suggests much more flexibility for further developments. Indeed, for
instance, in the proof of Theorem 4, one can notice that the whole analysis is
based on only

$\lambda_{i}^{k+1}=\mathcal{G}\lambda_{i-1}^{\tau_{i-1,1}^{i}(k)}+(\mathcal{F}-\mathcal{G})\lambda_{i-1}^{\tau_{i-1,2}^{i}(k)}+\zeta.$

Therefore, to fully account for the domain decomposition context, the
resulting error bound can be directly restated using $\mathcal{G}_{i}$ and
$\mathcal{F}_{i}$, with

$\widetilde{\alpha}=\max_{i}\|\mathcal{G}_{i}\|+\|\mathcal{F}_{i}-\mathcal{G}_{i}\|.$

Furthermore, without linearity assumption on functions $G_{i}$ and $F_{i}$,
the same result can be simply obtained, for

$\lambda_{i}^{k+1}=G_{i}(\lambda_{i-1}^{\tau_{i-1,1}^{i}(k)})+F_{i}(\lambda_{i-1}^{\tau_{i-1,2}^{i}(k)})-G_{i}(\lambda_{i-1}^{\tau_{i-1,2}^{i}(k)}),$

by generalizing the arbitrary norm $\|.\|$ to a distance function over the
space of functions $u^{(t)}$, with $u^{(t)}(s):=u(s,t)$.

In conclusion, it is interesting to see that, while asynchronous iterations
generally require more restrictive convergence conditions, they can be
successfully performed in any Parareal-based solver, for large numbers of time
subintervals (${p\to\infty}$), without verification of additional properties,
since convergence conditions are very nearly the same.

When, additionally, it happens that the Parareal iterative scheme itself does
not actually converge, a finite-time termination is still ensured, at $k=p$,
with the solution provided by a sequential application of $F$ from
$[T_{0},T_{1}]$ to $[T_{p-1},T_{p}]$ (see [14]). We use here the idea of
nested sets $\\{S^{(l)}\\}_{l\in\mathbb{N}}$ from Theorem 2 to show the
finite-time termination of Parareal asynchronous iterations. Note that
deriving nested sets of the form (7) applies only when the convergence
condition is satisfied, what is not assumed here for finite-time termination.
We therefore instead define $S^{(l)}$ so as to have the following property.

######  Proposition 1.

Parareal asynchronous iterations (12) generate a sets sequence
$\\{S^{(l)}\\}_{l\in\mathbb{N}}$ such that

$\forall l\geq p,\quad S^{(l)}=\\{\lambda^{*}\\}.$ (17)

###### Proof.

Let us define sets

$S^{(l)}=\left\\{\lambda\ |\ \forall
i\in\\{0,\ldots,\min\\{l,p\\}\\},\quad\lambda_{i}=\lambda_{i}^{*}\right\\},\quad\quad
l\in\mathbb{N},$

so that we have (17) and

$\lambda^{0}\in S^{(0)}.$ (18)

Now let

$\lambda^{(1)},\lambda^{(2)}\in
S^{(l)},\qquad\lambda^{(3)}=\widetilde{f}(\lambda^{(1)},\lambda^{(2)}),$

where $\widetilde{f}$ is the Parareal asynchronous iterations mapping (see
(13)). Then we have

$\left\\{\begin{array}[]{lcll}\lambda^{(3)}_{0}&=&\lambda_{0}^{*},&\\\
\lambda^{(3)}_{i}&=&G(\lambda_{i-1}^{*})+F(\lambda_{i-1}^{*})-G(\lambda_{i-1}^{*})=\lambda_{i}^{*},&\forall
i\in\\{1,\ldots,\min\\{l+1,p\\}\\},\end{array}\right.$

which implies that

$\forall\lambda^{(1)},\lambda^{(2)}\in
S^{(l)},\quad\widetilde{f}(\lambda^{(1)},\lambda^{(2)})\in S^{(l+1)}.$ (19)

At last, let $\Gamma$ denotes the set of all functions $u^{(t)}$, with
$u^{(t)}(s):=u(s,t)$, and let us define sets

$S^{(l)}_{i}=\left\\{\begin{array}[]{lll}\\{\lambda_{i}^{*}\\},&i\in\\{0,\ldots,\min\\{l,p\\}\\},&l\in\mathbb{N},\\\
\Gamma,&i\in\\{l+1,\ldots,p\\},&l\in\\{0,\ldots,p-1\\},\end{array}\right.$

so that we do verify

$S^{(l)}=S^{(l)}_{0}\times\cdots\times S^{(l)}_{p}.$ (20)

(18) to (20) ensure that Parareal asynchronous iterations (12) generate such a
sequence $\\{S^{(l)}\\}_{l\in\mathbb{N}}$, which concludes the proof. ∎

##  4 Performance analysis

###  4.1 Computational cost

Farhat & Chandesris [14] established a bound of $1/2$ on the parallel
efficiency (speedup per processor) of the Parareal scheme, assuming that at
least one corrective iteration would be needed (which implies at least $k=2$
iterations). It would therefore be of some interest to see at which extent
asynchronous iterations could improve the parallel performance, since it is
intuitively expected that, here as well, the finite-time termination causes
the possible time saving to be bounded.

Without questioning this efficiency bound, Aubanel [1] provided additional
insights by taking into account different possible implementations. We
consider here the distributed approach therein proposed, which features better
scaling potential. Basically, it consists of distributing the sequential
application of $G$ over the set of processors, as described by Algorithm 1. We
assume the case where we have $p$ processors, so that $\lambda_{i}$, with
$i\in\\{1,\ldots,p\\}$, is held by the processor $i$.

Algorithm 1 Distributed Parareal scheme (on processor $i\in\\{1,\ldots,p\\}$)

1: $\lambda_{i-1}^{0}:=u^{(0)}$

2: for all $j\in\\{1,\ldots,i-1\\}$ do

3: $\lambda_{i-1}^{0}:=G(\lambda_{i-1}^{0})$

4: end for

5: $w_{i}^{0}:=G(\lambda_{i-1}^{0})$

6: $\lambda_{i}^{0}:=w_{i}^{0}$

7: $k:=0$

8: Allow for message reception from processor $i-1$

9: repeat

10: if $i>k$ then

11: $v_{i}^{k}:=F(\lambda_{i-1}^{k})$

12: if $i-1>k$ then

13: Wait for reception of $\lambda_{i-1}^{k+1}$ from processor $i-1$

14: $w_{i}^{k+1}:=G(\lambda_{i-1}^{k+1})$

15: else

16: $w_{i}^{k+1}:=w_{i}^{k}$

17: end if

18: $\lambda_{i}^{k+1}:=w_{i}^{k+1}+v_{i}^{k}-w_{i}^{k}$

19: if $i<p$ then

20: Request sending of $\lambda_{i}^{k+1}$ to processor $i+1$

21: end if

22: else

23: $\lambda_{i}^{k+1}:=\lambda_{i}^{k}$

24: end if

25: $k:=k+1$

26: until $\lambda^{k}\simeq\lambda^{*}$

Let then $\mathcal{C}_{G}$ and $\mathcal{C}_{F}$ denote the computational
costs of applying $G$ and $F$, respectively. The cost of a sequential
application of $F$ over the whole time interval is thus given by

$\mathcal{C}(p,\mathcal{C}_{F})=p\ \mathcal{C}_{F}.$

What is remarkable with such a distributed implementation is that the
processor $i-1$ starts applying $F$ before the processor $i$. Therefore, in a
perfect load balance context, one may expect $\lambda_{i-1}^{k+1}$ to be
available to the processor $i$ at the moment when it finishes applying $F$,
what would result in a computational cost, after $k$ iterations, given by (see
[1])

$\mathcal{C}^{(k)}(p,\mathcal{C}_{F},\mathcal{C}_{G})=p\
\mathcal{C}_{G}+k(\mathcal{C}_{F}+\mathcal{C}_{G}),$

where $p\ \mathcal{C}_{G}$ corresponds to the cost of the initialization phase
in Algorithm 1.

While accurate enough within the scope of [1], this ideal case however assumes
negligible communication, vector addition and memory access costs, which is
questionable in real world situations. Let then, at least,
$\mathcal{C}_{C,G}(i,k)$ give communication overhead costs which accumulate
with $\mathcal{C}_{G}$ on the processor $i$, at the iteration $k$, without
being overlapped by the application of $F$ on the processor $p$. It yields:

$\mathcal{C}^{(k)}(p,\mathcal{C}_{F},\mathcal{C}_{G},\mathcal{C}_{C,G})=p\
\mathcal{C}_{G}+k(\mathcal{C}_{F}+\mathcal{C}_{G})+\sum_{l=1}^{k}\sum_{j=l+1}^{p-1}\mathcal{C}_{C,G}(j,l).$
(21)

######  Remark.

Actually, not only communication cost is not negligible, but it can even be
higher than $\mathcal{C}_{G}$, since both are closely proportional to the size
of $\lambda_{i}$. Note that $\mathcal{C}_{C,G}(j,l)$ is of the form

$\mathcal{C}_{C,G}(j,l)=\mathcal{C}_{C}(j,l)+\beta_{j,l}\
\mathcal{C}_{G},\qquad\beta_{j,l}\in\\{0,1\\},$

with

$\mathcal{C}_{C}(j,l)>0\quad\iff\quad\beta_{j,l}=1,$

where $\mathcal{C}_{C}(j,l)$ is a non-overlapped communication cost.

Still, let us assume a constant uniform communication cost and, for instance,
a four-steps fine integration

$F=\mathring{F}\circ\mathring{F}\circ\mathring{F}\circ\mathring{F},$

with $\mathring{F}$ being some one-step integration which costs nearly as much
as $G$. Then, by introducing non-negligible communication $C$ into Aubanel’s
model [1], Algorithm 1 would behave like

$\begin{array}[]{cclcccllll}\text{proc
1}&G\mathring{F}\mathring{F}\mathring{F}\mathring{F},&&&&&&&&\\\ &&C&&&&&&&\\\
\text{proc
2}&GG\mathring{F}\mathring{F}\mathring{F}&\mathring{F}G&\mathring{F}&\mathring{F}&\mathring{F}&\mathring{F},&&&\\\
&&&C&&&&C&&\\\ \text{proc
3}&GGG\mathring{F}\mathring{F}&\mathring{F}\mathring{F}&-&G&\mathring{F}&\mathring{F}&\mathring{F}\mathring{F}G&\mathring{F}\mathring{F}\mathring{F}\mathring{F},&\\\
&&&&&C&&&C&C\\\ \text{proc
4}&GGGG\mathring{F}&\mathring{F}\mathring{F}&\mathring{F}&-&-&G&\mathring{F}\mathring{F}\mathring{F}&\mathring{F}G\mathring{F}\mathring{F}&\mathring{F}\mathring{F}G\mathring{F}\mathring{F}\mathring{F}\mathring{F},\end{array}$

which shows that the overhead would lie only in the initialization or the
first iteration phase, implying that

$\sum_{j=l+1}^{p-1}\mathcal{C}_{C,G}(j,l)=0,\qquad\forall l>1.$

Our general cost model including $\mathcal{C}_{C,G}$ therefore becomes handy
when variations can occur in the communication delays, due for instance to a
fault recovery procedure or simple perturbation on the network. We see from
the example that a little cost increase in the second message transfer toward
the processor 4 would result in an overhead cost. At the same time, the impact
of these variations can also be expected to decrease with $k$. Nevertheless,
instead of an exhaustive analysis, we are rather interested in assessing a
global potential effect of communication. Let us then consider a rough
estimation, $\mathcal{\overline{C}}_{C,G}$, of the average overhead cost at
each iteration, on each processor. This leads to:

######  Proposition 2.

The computational cost of $k$ distributed Parareal iterations (Algorithm 1) is
of the form

$\mathcal{C}^{(k)}(p,\mathcal{C}_{F},\mathcal{C}_{G},\mathcal{\overline{C}}_{C,G})=p\
\mathcal{C}_{G}+k\left(\mathcal{C}_{F}+\mathcal{C}_{G}+\left(p-1-\frac{k+1}{2}\right)\mathcal{\overline{C}}_{C,G}\right).$

###### Proof.

From (21), we have

$\displaystyle\mathcal{C}^{(k)}(p,\mathcal{C}_{F},\mathcal{C}_{G},\mathcal{\overline{C}}_{C,G})$
$\displaystyle=p\
\mathcal{C}_{G}+k(\mathcal{C}_{F}+\mathcal{C}_{G})+\sum_{l=1}^{k}(p-1-l)\
\mathcal{\overline{C}}_{C,G}$ $\displaystyle=p\
\mathcal{C}_{G}+k(\mathcal{C}_{F}+\mathcal{C}_{G})+\left(k(p-1)-\frac{k(k+1)}{2}\right)\mathcal{\overline{C}}_{C,G}.$

The proposition thus follows. ∎

######  Remark.

While $k$ might be expected to remain relatively small, the effects of
$\mathcal{\overline{C}}_{C,G}$ however increase with the number $p$ of
processors. It is important to note that accounting such overhead costs does
not contradict the perfect load balance assumption. In this sense, this would
be a minimum relevant ideal model for the distributed Parareal scheme.

Let us finally introduce asynchronous iterations into Algorithm 1, to obtain
Algorithm 2.

Algorithm 2 Distributed async-Parareal scheme (on processor
$i\in\\{1,\ldots,p\\}$)

1: $\lambda_{i-1}^{0}:=u^{(0)}$

2: for all $j\in\\{1,\ldots,i-1\\}$ do

3: $\lambda_{i-1}^{0}:=G(\lambda_{i-1}^{0})$

4: end for

5: $w_{i}^{0}:=G(\lambda_{i-1}^{0})$

6: $\lambda_{i}^{0}:=w_{i}^{0}$

7: $k_{i}:=0$

8: Allow for message reception from processor $i-1$

9: repeat

10: $v_{i}^{k_{i}}:=F(\lambda_{i-1}^{k_{i}})$

11: $w_{i}^{k_{i}+1}:=G(\lambda_{i-1}^{k_{i}+1})$

12: $\lambda_{i}^{k_{i}+1}:=w_{i}^{k_{i}+1}+v_{i}^{k_{i}}-w_{i}^{k_{i}}$

13: if $i<p$ then

14: Request sending of $\lambda_{i}^{k_{i}+1}$ to processor $i+1$

15: end if

16: $k_{i}:=k_{i}+1$

17: until
$(\lambda_{0}^{0},\lambda_{1}^{k_{1}},\ldots,\lambda_{p}^{k_{p}})\simeq\lambda^{*}$

Still assuming perfect load balance, the corresponding computational cost is
that of the processor having performed the most iterations. To make connection
with the implicit global asynchronous iterations variable, say
$\widetilde{k}$, let us define, for $k_{i}$ iterations on processor $i$,

$\kappa_{i}(\widetilde{k}):=k_{i}=\operatorname{card}\\{l\leq\widetilde{k}\ |\
i\in P^{(l)}\\},$

and set

$\kappa(\widetilde{k}):=\max_{i\in\\{1,\ldots,p\\}}\kappa_{i}(\widetilde{k}).$

Then we have:

######  Proposition 3.

The computational cost of $\widetilde{k}$ distributed Parareal asynchronous
iterations (Algorithm 2) is of the form

$\mathcal{\widetilde{C}}^{(\widetilde{k})}(p,\mathcal{C}_{F},\mathcal{C}_{G})=p\
\mathcal{C}_{G}+\kappa(\widetilde{k})(\mathcal{C}_{F}+\mathcal{C}_{G}).$

###### Proof.

This is straightforward, since communication is completely overlapped. ∎

###  4.2 Performance gain

The error bounds from Theorem 3 and Theorem 4 indicate that the convergence
speeds of Parareal synchronous and asynchronous iterations depend on how small
the respective factors

$\alpha=\frac{1-\|\mathcal{G}\|^{p}}{1-\|\mathcal{G}\|}\|\mathcal{F}-\mathcal{G}\|,\qquad\widetilde{\alpha}=\|\mathcal{G}\|+\|\mathcal{F}-\mathcal{G}\|$

are. As intuitively expected, we verify the following comparison.

######  Proposition 4.

$\widetilde{\alpha}<1\quad\implies\quad\alpha<\widetilde{\alpha}.$

###### Proof.

Since $\|\mathcal{G}\|<1$, we have

$\alpha<\frac{1}{1-\|\mathcal{G}\|}\|\mathcal{F}-\mathcal{G}\|.$

Assume that

$\widetilde{\alpha}\leq\frac{1}{1-\|\mathcal{G}\|}\|\mathcal{F}-\mathcal{G}\|.$

Then we have

$\displaystyle(\|\mathcal{G}\|+\|\mathcal{F}-\mathcal{G}\|)(1-\|\mathcal{G}\|)$
$\displaystyle\leq\|\mathcal{F}-\mathcal{G}\|,$
$\displaystyle\|\mathcal{G}\|(1-\|\mathcal{G}\|)+\|\mathcal{F}-\mathcal{G}\|(1-\|\mathcal{G}\|-1)$
$\displaystyle\leq 0,$
$\displaystyle\|\mathcal{G}\|(1-\|\mathcal{G}\|-\|\mathcal{F}-\mathcal{G}\|)$
$\displaystyle\leq 0,$ $\displaystyle
1-\|\mathcal{G}\|-\|\mathcal{F}-\mathcal{G}\|$ $\displaystyle\leq 0,$
$\displaystyle\widetilde{\alpha}$ $\displaystyle\geq 1.$

Consequently,

$\widetilde{\alpha}<1\quad\implies\quad\frac{1}{1-\|\mathcal{G}\|}\|\mathcal{F}-\mathcal{G}\|<\widetilde{\alpha},$

which concludes the proof. ∎

Nonetheless, if we consider synchronous and asynchronous iterations variables
$k$ and $\widetilde{k}$, respectively, we should bear in mind that
$\widetilde{k}$ grows much faster than $k$. Therefore, we can reasonably
expect that $\sigma(\widetilde{k})$ also grows faster than $k$, sufficiently
to make Parareal asynchronous iterations reach the desired precision sooner
than synchronous ones do.

To be more practical, besides convergence rate, we see from Algorithm 2 that,
either $\lambda_{i-1}^{k_{i}+1}$ is always an updated available input for the
processor $i$, and hence, this reduces to a classical Parareal execution, or
it happens that

$\lambda_{i-1}^{k_{i}+1}=\lambda_{i-1}^{k_{i}}$

for some processors $i$, and hence, one has

$\lambda_{i}^{k_{i}+1}=F(\lambda_{i-1}^{k_{i}}).$

The advantage of the asynchronous scheme is to account for the fact that this
could be enough to reach the desired precision, instead of waiting for an
updated input in order to necessarily output an exact Parareal update. Note
that, in any case then, we should have

$\kappa(\widetilde{k})\geq k,$

where $k$ would be the number of synchronous iterations. From our cost model
in Proposition 2 and Proposition 3, we deduce:

######  Corollary 3.

The speedup over classical distributed Parareal is bounded as

$\frac{\mathcal{C}^{(k)}(p,\mathcal{C}_{F},\mathcal{C}_{G},\mathcal{\overline{C}}_{C,G})}{\mathcal{\widetilde{C}}^{(\widetilde{k})}(p,\mathcal{C}_{F},\mathcal{C}_{G})}\
\leq\
1+\frac{\left(p-2\right)\mathcal{\overline{C}}_{C,G}}{\mathcal{C}_{F}+\mathcal{C}_{G}}.$

###### Proof.

If

$\mathcal{C}^{(k)}(p,\mathcal{C}_{F},\mathcal{C}_{G},\mathcal{\overline{C}}_{C,G})<\mathcal{\widetilde{C}}^{(\widetilde{k})}(p,\mathcal{C}_{F},\mathcal{C}_{G}),$

then

$\frac{\mathcal{C}^{(k)}(p,\mathcal{C}_{F},\mathcal{C}_{G},\mathcal{\overline{C}}_{C,G})}{\mathcal{\widetilde{C}}^{(\widetilde{k})}(p,\mathcal{C}_{F},\mathcal{C}_{G})}\
<\ 1\ \leq\
1+\frac{\left(p-2\right)\mathcal{\overline{C}}_{C,G}}{\mathcal{C}_{F}+\mathcal{C}_{G}},\qquad\qquad
p\geq 2.$

Now, considering that

$\mathcal{C}^{(k)}(p,\mathcal{C}_{F},\mathcal{C}_{G},\mathcal{\overline{C}}_{C,G})\geq\mathcal{\widetilde{C}}^{(\widetilde{k})}(p,\mathcal{C}_{F},\mathcal{C}_{G}),$

we have

$\displaystyle\frac{\mathcal{C}^{(k)}(p,\mathcal{C}_{F},\mathcal{C}_{G},\mathcal{\overline{C}}_{C,G})}{\mathcal{\widetilde{C}}^{(\widetilde{k})}(p,\mathcal{C}_{F},\mathcal{C}_{G})}$
$\displaystyle\leq\frac{\mathcal{C}^{(k)}(p,\mathcal{C}_{F},\mathcal{C}_{G},\mathcal{\overline{C}}_{C,G})-p\
\mathcal{C}_{G}}{\mathcal{\widetilde{C}}^{(\widetilde{k})}(p,\mathcal{C}_{F},\mathcal{C}_{G})-p\
\mathcal{C}_{G}}$
$\displaystyle=\frac{k}{\kappa(\widetilde{k})}\left(1+\frac{\left(p-1-\frac{k+1}{2}\right)\mathcal{\overline{C}}_{C,G}}{\mathcal{C}_{F}+\mathcal{C}_{G}}\right),$

where we recall that $p\ \mathcal{C}_{G}$ is the initialization cost. With the
maximum gain expected at $\kappa(\widetilde{k})=k$, this bound reaches its
maximum at $k=1$, which leads to the result. ∎

Still, we should point out the fact that, in HPC context, and particularly for
time-parallel time-integration (see, e.g., introduction and comments in [14,
1, 32]), one is generally more interested in being able to improve performance
by providing more and more processors, in order to keep a given amount of work
per processor (strong scaling). In this sense, the bound of 50% over the
traditional parallel efficiency is not actually limiting the desired scaling
potential. Let us then highlight, for a more complete performance analysis,
that the asymptotic speedup from the distributed Parareal over sequential
application of $F$ is given by

$\lim_{p\to+\infty}\frac{\mathcal{C}(p,\mathcal{C}_{F})}{\mathcal{C}^{(k)}(p,\mathcal{C}_{F},\mathcal{C}_{G},\mathcal{\overline{C}}_{C,G})}\
=\ \frac{\mathcal{C}_{F}}{\mathcal{C}_{G}+k\ \mathcal{\overline{C}}_{C,G}}\
\leq\ \frac{\mathcal{C}_{F}}{\mathcal{C}_{G}+\mathcal{\overline{C}}_{C,G}},$

while that of the distributed async-Parareal reduces to

$\lim_{p\to+\infty}\frac{\mathcal{C}(p,\mathcal{C}_{F})}{\mathcal{\widetilde{C}}^{(\widetilde{k})}(p,\mathcal{C}_{F},\mathcal{C}_{G})}\
=\ \frac{\mathcal{C}_{F}}{\mathcal{C}_{G}},$

which finally implies that

$\lim_{p\to+\infty}\frac{\mathcal{C}^{(k)}(p,\mathcal{C}_{F},\mathcal{C}_{G},\mathcal{\overline{C}}_{C,G})}{\mathcal{\widetilde{C}}^{(\widetilde{k})}(p,\mathcal{C}_{F},\mathcal{C}_{G})}\
=\ 1+k\ \frac{\mathcal{\overline{C}}_{C,G}}{\mathcal{C}_{G}}.$

As expected, the performance gain from async-Parareal should increase with
$\mathcal{\overline{C}}_{C,G}$, even in a perfectly well balanced
configuration, and we saw that $\mathcal{\overline{C}}_{C,G}$ can be of the
same order as $\mathcal{C}_{G}$. It may however seem surprising that the
asymptotic speedup does depend on $k$, and not on $\kappa(\widetilde{k})$ too.
When looking at the classical speedup though (in the proof of Corollary 3), we
may expect that, for given $\mathcal{\overline{C}}_{C,G}$, $F$ and $G$, the
ratio $k/\kappa(\widetilde{k})$ is bounded by a constant. This is very likely
to be true, as the current async-Parareal scheme is also quite close to the
classical Parareal one. Further analysis may be needed to fully assess such a
relation. Still in this case, while the speedup would decrease with $k$, the
asymptotic speedup rather increases. Indeed, in the distributed Parareal
algorithm, the amount of communication decreases with $k$ (as partial
termination is taken into account), which lessens the impact of
$\mathcal{\overline{C}}_{C,G}$. However, for very high numbers of processors,
this diminution is generally negligible, as $p-k$ remains significantly high
for a long period of time.

In the end, while the practical performance gain from applying asynchronous
iterations to an existing distributed scheme is no more questionable for
spatial parallelism on PDEs, the case of time parallelism is quite less
obvious. We therefore present here experimental results only conducted to
confirm the possibility of performance gain here as well. Especially, such a
configuration where each processor depends on only one other processor is
probably the worst for squeezing the most out of asynchronous iterations.
Furthermore, to give full advantage to the classical distributed Parareal, we
experimented using perfect static load balance, on a homogeneous supercomputer
with high-speed communication (which minimizes
$\mathcal{\overline{C}}_{C,G}$). We reported another, more detailed,
experimental investigation in [28], where we mainly focused, throughout the
paper, on implementation aspects within the Message Passing Interface (MPI)
framework. This issue about asynchronous iterations programming has been more
generally discussed in [25, 27], where an MPI-based programming library has
been proposed, including tools to evaluate the loop stopping criterion of
Algorithm 2 (see [26] for a full introduction to exact residual error
evaluation under asynchronous iterations).

##  5 Numerical results

###  5.1 Problem and experimental settings

The experimental case consisted of simulating an apartment gradually cooled by
an air conditioner. We used a simple heat evolution model of the form

$\frac{\partial u}{\partial t}=\Delta u,$

by setting a constant Dirichlet boundary condition on the air conditioner to
simulate a continuous cooling at 73.4 degrees Fahrenheit (23 degrees Celsius),
from an initial homogeneous temperature at 86.0 degrees Fahrenheit (30 degrees
Celsius). Spatial discrete equations were obtained through $P_{1}$-Lagrange
finite-elements integration upon a mesh composed of 171,478 tetrahedrons and
33,796 nodes (see Figure 1), generated by the TetGen [35] free software. The
computer aided design (CAD) was performed using the CATIA V5 proprietary
software tool ($\copyright$Dassault Systèmes).

Figure 1: Finite-elements mesh of the experimental 3D model.

Operators $G$ and $F$ consisted of either backward Euler or trapezoidal rule
time integration, with respective time steps

$\Delta T=0.2,\qquad\delta t=0.002,$

which implied that

$\frac{\mathcal{C}_{F}}{\mathcal{C}_{G}}\simeq 100.$

Each spatial linear system resolution was sequential and based on LU
factorization. Figure 2 shows the evolution of the temperature over time.

Figure 2: Heat distribution at times $t=1$, $t=10$ and $t=20$.

We performed several parallel runs on an SGI Altix ICE supercomputer of 68
nodes interconnected by a QDR Infiniband network (40Gb/s). Each node consisted
of two six-cores Intel Xeon X5650 CPUs (12 cores per node) at 2.66GHz, and
21GB RAM allocated for running jobs. The MPI library SGI-MPT was loaded as
communication middleware. Execution parameters included the number $p$ of MPI
processes, which always equals the number of time intervals, and the loop
stopping criterion

$\|\lambda^{k}-\lambda^{k-1}\|_{\infty}<\varepsilon,\qquad\varepsilon\in\\{10^{-6},10^{-5}\\}.$
(22)

We refer to [27] for the implementation of such a criterion in case of
asynchronous iterations $k_{1}$ to $k_{p}$.

###  5.2 Performance comparison

Table 1 and Table 2 report average results, where $\mathcal{T}$ denotes wall-
clock execution times in seconds. We recall that $T_{p}$ is the simulated end
time, and that $\mathcal{\overline{C}}_{C,G}$ is the average non-overlapped
overhead costs induced by communication at each iteration. This does not
directly corresponds to communication costs, and therefore is hardly
accurately predictable. We rather deduced values
$\mathcal{\overline{C}}_{C,G}$ which make $\mathcal{C}^{(k)}$ match the
measured execution times. We also measured, on average,

$\mathcal{C}_{G}\simeq 0.14\text{ sec}.$

Compared to these $\mathcal{\overline{C}}_{C,G}$ values, it seems clear that,
even in such an ideal HPC computational environment, the actual impact of
communication costs cannot be neglected, despite the optimized distributed
implementation from Aubanel [1].

Table 1: Comparison of synchronous and asynchronous Parareal.
$F$: trapezoidal rule, $G$: backward Euler, $\varepsilon$: $10^{-6}$. | Parareal | Async-Parareal
---|---|---
| $p$ | $T_{p}$
---|---
16 | 3.2
24 | 4.8
32 | 6.4
48 | 9.6
64 | 12.8
90 | 18.0
| $\mathcal{T}$ | $k$ | $\mathcal{\overline{C}}_{C,G}$ | $\|\lambda^{k}-\lambda^{*}\|_{\infty}$
---|---|---|---
280 | 10 | 1.530 | 1.49E-07
420 | 14 | 1.010 | 1.86E-07
691 | 20 | 1.011 | 3.75E-08
1169 | 32 | 0.727 | 4.95E-11
1488 | 44 | 0.486 | 3.75E-08
2676 | 64 | 0.486 | 4.75E-11
| $\mathcal{T}$ | $\kappa(\widetilde{k})$ | $\mathcal{\widetilde{C}}^{(\widetilde{k})}$ | $\|\lambda^{\widetilde{k}}-\lambda^{*}\|_{\infty}$
---|---|---|---
342 | 24 | 337 | 5.01E-08
502 | 34 | 484 | 1.06E-07
644 | 44 | 627 | 1.71E-08
922 | 66 | 940 | 3.33E-11
1255 | 92 | 1310 | 2.79E-10
1938 | 149 | 2119 | 4.78E-11
Table 2: Comparison of synchronous and asynchronous Parareal.
$F$: backward Euler, $G$: backward Euler, $\varepsilon$: $10^{-5}$. | Parareal | Async-Parareal
---|---|---
| $p$ | $T_{p}$
---|---
16 | 3.2
24 | 4.8
32 | 6.4
48 | 9.6
| $\mathcal{T}$ | $k$ | $\mathcal{\overline{C}}_{C,G}$ | $\|\lambda^{k}-\lambda^{*}\|_{\infty}$
---|---|---|---
116 | 6 | 0.233 | 4.43E-07
131 | 6 | 0.307 | 4.99E-07
126 | 6 | 0.211 | 3.41E-07
139 | 6 | 0.196 | 1.32E-07
| $\mathcal{T}$ | $\kappa(\widetilde{k})$ | $\mathcal{\widetilde{C}}^{(\widetilde{k})}$ | $\|\lambda^{\widetilde{k}}-\lambda^{*}\|_{\infty}$
---|---|---|---
348 | 18 | 293 | 4.68E-08
379 | 19 | 297 | 6.44E-08
382 | 19 | 276 | 9.35E-08
360 | 15 | 211 | 1.48E-07

As expected, Table 1 shows that the overhead costs at one iteration tend to
diminish as $k$ grows, which makes the average $\mathcal{\overline{C}}_{C,G}$
diminish as well, while the sum of these costs can however keep increasing,
which explains the increasing better performance of asynchronous execution. On
the other hand, estimated values $\mathcal{\widetilde{C}}^{(\widetilde{k})}$
are quite close to measured times (mostly in Table 1). This is explained by
the fact that asynchronous execution somehow converts overhead communication
costs into additional iterations which are accounted by
$\kappa(\widetilde{k})$. Dealing with the prediction of
$\kappa(\widetilde{k})$ could be a less harassing task, compared to
$\mathcal{\overline{C}}_{C,G}$. One can notice for instance that

$2<\frac{\kappa(\widetilde{k})}{k}<3,\qquad\qquad
2<\frac{\kappa(\widetilde{k})}{k}<4,$

in Table 1 and Table 2 respectively, which confirms, for these test cases, our
expectation of a generally bounded ratio.

Finally, with a stopping criterion (22), it happened that async-Parareal
generally tended to provide a better precision, according to the sequential
solution from $F$. While this was probably due to the particular termination
delay of asynchronous iterations, it did not prevent us, in case of Table 1,
from obtaining a faster solver, for $p\geq 32$. For instance, around 12
minutes are saved over 44, using 90 processors. On the contrary, Table 2 shows
cases where, with sufficiently low values $\mathcal{\overline{C}}_{C,G}$
(regarding numbers of iterations), synchronous iterations keep performing
better than asynchronous ones.

##  6 Conclusion

Efficiently applying asynchronous iterations in time domain decomposition is
rather challenging, and this primary work certainly leaves an amount of
further desirable benchmarks and improvements. While we tackled here the basic
Parareal scheme, it is somehow important to note that introducing asynchronous
iterations is a design approach quite different from, and not concurrent to,
classical scheme improvement, as they could actually be applied to more
advanced Parareal schemes as well, as long as communication is a possible
factor of performance limit in real situations. What is particularly
interesting is that convergence conditions (in maximum error norm) barely
change, and are even asymptotically the same, while also finite-time
termination remains guaranteed. Still, despite the potential theoretical
performance gain and some experimental results confirming it in practice, it
is clear that such a pipeline parallel configuration drastically limits the
benefits from asynchronous iterations, compared to spatial parallelism where
each processor generally depends on several others, such as in undirected grid
configurations. Hopefully, as pointed out by this first theoretical analysis,
much more flexibility is now provided to deal with this issue, which is
currently under consideration.

## References

  * [1] E. Aubanel.  Scheduling of tasks in the parareal algorithm.  Parallel Computing, 37(3):172–182, 2011.
  * [2] J. Bahi, E. Griepentrog, and J. C. Miellou.  Parallel treatment of a class of differential-algebraic systems.  SIAM Journal on Numerical Analysis, 33(5):1969–1980, 1996.
  * [3] J. M. Bahi, S. Contassot-Vivier, and R. Couturier.  Parallel Iterative Algorithms: from sequential to grid computing.  Numerical Analysis and Scientific Computing. Chapman & Hall/CRC, 2007\.
  * [4] G. M. Baudet.  Asynchronous iterative methods for multiprocessors.  J. ACM, 25(2):226–244, 1978.
  * [5] D. P. Bertsekas.  Distributed asynchronous computation of fixed points.  Mathematical Programming, 27(1):107–120, 1983.
  * [6] D. P. Bertsekas and J. N. Tsitsiklis.  Parallel and Distributed Computation: Numerical Methods.  Prentice-Hall, Inc., Upper Saddle River, NJ, USA, 1989.
  * [7] D. P. Bertsekas and J. N. Tsitsiklis.  Some aspects of parallel and distributed iterative algorithms – A survey.  Automatica, 27(1):3–21, 1991.
  * [8] E. D. Chajakis and S. A. Zenios.  Synchronous and asynchronous implementations of relaxation algorithms for nonlinear network optimization.  Parallel Computing, 17(8):873–894, 1991.
  * [9] P. Chartier and B. Philippe.  A parallel shooting technique for solving dissipative ODE’s.  Computing, 51(3):209–236, 1993.
  * [10] M. Chau, D. E. Baz, R. Guivarch, and P. Spiteri.  MPI implementation of parallel subdomain methods for linear and nonlinear convection–diffusion problems.  Journal of Parallel and Distributed Computing, 67(5):581–591, 2007\.
  * [11] M. Chau, R. Couturier, J. Bahi, and P. Spiteri.  Parallel solution of the obstacle problem in grid environments.  International Journal of High Performance Computing Applications, 25(4):488–495, 2011.
  * [12] D. Chazan and W. Miranker.  Chaotic relaxation.  Linear Algebra and its Applications, 2(2):199–222, 1969.
  * [13] M. N. El Tarazi.  Some convergence results for asynchronous algorithms.  Numerische Mathematik, 39(3):325–340, 1982.
  * [14] C. Farhat and M. Chandesris.  Time-decomposed parallel time-integrators: theory and feasibility studies for fluid, structure, and fluid–structure applications.  International Journal for Numerical Methods in Engineering, 58(9):1397–1434, 2003.
  * [15] A. Frommer, H. Schwandt, and D. B. Szyld.  Asynchronous weighted additive Schwarz methods.  Electronic Transactions on Numerical Analysis, 5:48–61, 1997.
  * [16] A. Frommer and D. B. Szyld.  On asynchronous iterations.  Journal of Computational and Applied Mathematics, 123(1–2):201–216, 2000.
  * [17] M. J. Gander.  Schwarz methods over the course of time.  Electronic Transactions on Numerical Analysis, 31:228–255, 2008\.
  * [18] M. J. Gander and S. Vandewalle.  Analysis of the Parareal time-parallel time-integration method.  SIAM J. Sci. Comput., 29(2):556–578, 2007.
  * [19] L. Hart and S. McCormick.  Asynchronous multilevel adaptive methods for solving partial differential equations on multiprocessors: Basic ideas.  Parallel Computing, 12(2):131–144, 1989.
  * [20] B. Leimkuhler.  Timestep acceleration of waveform relaxation.  SIAM Journal on Numerical Analysis, 35(1):31–50, 1998.
  * [21] E. Lelarasmee, A. E. Ruehli, and A. L. Sangiovanni-Vincentelli.  The waveform relaxation method for time-domain analysis of large scale integrated circuits.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 1(3):131–145, 1982.
  * [22] K. Li-Shan, C. Yu-Ping, S. Le-Lin, and Q. Hui-Yun.  The asynchronous parallel algorithms s-cor for solving p.d.e.’s on multiprocessors.  International Journal of Computer Mathematics, 18(2):163–172, 1985\.
  * [23] J.-L. Lions, Y. Maday, and G. Turinici.  Résolution d’EDP par un schéma en temps “pararéel”.  C. R. Acad. Sci. Paris Sér. I Math., 332(7):661 – 668, 2001.
  * [24] P. L. Lions.  On the Schwarz alternating method II.  In T. Chan, R. Glowinski, G. A. Meurant, J. Périaux, and O. Widlund, editors, Domain Decomposition Methods for Partial Differential Equations II, pages 47–70. Society for Industrial and Applied Mathematics, Philadelphia, 1989.
  * [25] F. Magoulès and G. Gbikpi-Benissan.  JACK: An asynchronous communication kernel library for iterative algorithms.  The Journal of Supercomputing, 73(8):3468–3487, 2017.
  * [26] F. Magoulès and G. Gbikpi-Benissan.  Distributed convergence detection based on global residual error under asynchronous iterations.  IEEE Transactions on Parallel and Distributed Systems, 29(4):819–829, 2018.
  * [27] F. Magoulès and G. Gbikpi-Benissan.  JACK2: An MPI-based communication library with non-blocking synchronization for asynchronous iterations.  Advances in Engineering Software, 119:116–133, 2018.
  * [28] F. Magoulès, G. Gbikpi-Benissan, and Q. Zou.  Asynchronous iterations of Parareal algorithm for option pricing models.  Mathematics, 6(4):1–18, 2018.  Paper 45.
  * [29] F. Magoulès, D. B. Szyld, and C. Venet.  Asynchronous optimized Schwarz methods with and without overlap.  Numerische Mathematik, 137(1):199–227, 2017.
  * [30] F. Magoulès and C. Venet.  Asynchronous iterative sub-structuring methods.  Mathematics and Computers in Simulation, 145(Supplement C):34–49, 2018.
  * [31] J. C. Miellou.  Algorithmes de relaxation chaotique à retards.  ESAIM: Mathematical Modelling and Numerical Analysis - Modélisation Mathématique et Analyse Numérique, 9(R1):55–82, 1975.
  * [32] A. S. Nielsen and J. S. Hesthaven.  Fault tolerance in the Parareal method.  In Proceedings of the ACM Workshop on Fault-Tolerance for HPC at Extreme Scale, FTXS ’16, pages 1–8, New York, NY, USA, 2016. ACM.
  * [33] J. L. Rosenfeld.  A case study in programming for parallel-processors.  Communications of the ACM, 12(12):645–655, 1969.
  * [34] S. Schechter.  Relaxation methods for linear equations.  Communications on Pure and Applied Mathematics, 12(2):313–335, 1959\.
  * [35] H. Si.  Tetgen, a delaunay-based quality tetrahedral mesh generator.  ACM Trans. Math. Softw., 41(2):11:1–11:36, 2015.
  * [36] P. Spitéri, J.-C. Miellou, and D. El Baz.  Asynchronous Schwarz alternating method with flexible communication for the obstacle problem.  Réseaux et Systèmes Répartis - Calculateurs Parallèles, 13(1):47–66, 2001.
  * [37] P. Spiteri, J.-C. Miellou, and D. El Baz.  Parallel asynchronous schwarz and multisplitting methods for a nonlinear diffusion problem.  Numerical Algorithms, 33(1):461–474, 2003.
  * [38] P. L. Tallec.  Domain decomposition methods in computational mechanics.  In J. T. Oden, editor, Computational Mechanics Advances, volume 1 (2), pages 121–220. North-Holland, 1994.