A Parareal algorithm without Coarse Propagator?
Martin J. Gander Mario Ohlberger Stephan Rave








A Parareal algorithm without Coarse Propagator?
    Martin J. GanderUniversité de Genève, martin.gander@unige.ch
Mario OhlbergerUniversity of Münster, mario.ohlberger@uni-muenster.de
Stephan RaveUniversity of Münster, stephan.rave@uni-muenster.de
    September 4, 2024
======================================================================================================================================================================================================



The Parareal algorithm was invented in 2001 in order to
  parallelize the solution of evolution problems in the time
  direction. It is based on parallel fine time propagators called F
  and sequential coarse time propagators called G, which alternatingly
  solve the evolution problem and iteratively converge to the fine
  solution. The coarse propagator G is a very important component of
  Parareal, as one sees in the convergence analyses. We present here
  for the first time a Parareal algorithm without coarse propagator,
  and explain why this can work very well for parabolic problems. We
  give a new convergence proof for coarse propagators approximating in
  space, in contrast to the more classical coarse propagators which
  are approximations in time, and our proof also applies in the
  absence of the coarse propagator. We illustrate our theoretical
  results with numerical experiments, and also explain why this
  approach can not work for hyperbolic problems.




§ INTRODUCTION


The Parareal algorithm was introduced in <cit.> as
a non-intrusive way to parallelize time stepping methods in the time
direction for solving partial differential equations. Its convergence
is now well understood: for linear problems of parabolic type,
Parareal converges superlinearly on bounded time intervals, and
satisfies a linear convergence bound on unbounded time intervals
<cit.>, which means it is a contraction even if
computations are performed on arbitrary long time intervals. It was
shown in the same reference also that for hyperbolic problems, the
linear convergence estimate over long time does not predict
contraction, and the superlinear estimate only indicates contraction
when there are already too many iterations performed and one is
approaching the finite step convergence property of Parareal. A
non-linear convergence analysis for Parareal can be found in
<cit.>, and also the performance of Parareal for
Hamiltonian problems is well understood, see
<cit.>, where also a derivative variant of Parareal
was proposed.

An essential ingredient in the Parareal algorithm is the coarse
propagator, and it is its accuracy that has a decisive influence on
the convergence of Parareal, as one can see in the estimates from
<cit.>. We consider here a
Parareal algorithm in the standard form for solving an evolution
problem ∂_t u= F(u,t), and time partition
0=T_0<T_1<…<T_N=T,

    U_n+1^k+1=F(U_n^k,T_n,T_n+1)+G(U_n^k+1,T_n,T_n+1)-G(U_n^k,T_n,T_n+1),

where the fine solver F(U_n^k,T_n,T_n+1) and the coarse solver
G(U_n^k,T_n,T_n+1) solve the underlying evolution problem on time
intervals Δ T:=T_n+1-T_n,

    ∂_t u= F(u,t),   u(T_n)=U_n^k,

with different accuracy. We are interested in understanding what
happens when we remove the coarse solver in the Parareal algorithm
(<ref>), i.e.  we run instead the iteration

    U_n+1^k+1=F(U_n^k,T_n,T_n+1).

Note that this is quite different from the Identity Parareal algorithm
<cit.>, where the coarse propagator was replaced by
the identity, a very coarse approximation, while here we remove the
coarse propagator altogether.




§ A NUMERICAL EXPERIMENT


We start with a numerical experiment for the one dimensional heat
equation on the spatial domain Ω:=(0,1) and the time interval (0,T],

    ∂_t u=∂_xxu+f,   u(x,0)=u_0(x),   u(0,t)=u(1,t)=0.

We show in Figure <ref>

the results obtained with zero initial conditions, u_0(x)=0, and the
source term

    f(x,t)=10e^-100(x-0.5)^2∑_j=1^4e^-100(t-t_j)^2,  
      t_j=0.1,0.6,1.35,1.85,

like a heating device which is turned on and of at specific time
instances. We discretized the heat equation (<ref>)
using a centered finite difference discretization in space with mesh
size Δ x=1/128, and Backward Euler in time with time
step Δ t=1/96, and run the Parareal algorithm up to
T=3 with various numbers N of coarse time intervals, where the
coarse propagator G just does one big Backward Euler step.  We see
two really interesting results: first, the Parareal algorithm without
coarse propagator is actually also converging when applied to the heat
equation, both methods are a contraction. And second, if there are not
many coarse time intervals, i.e. if the coarse time interval length
Δ T is becoming large, the Parareal algorithm without coarse
correction is converging even faster than with coarse correction!

In order to test this surprising property further, we repeat the
above experiment, but use now homogeneous Neumann boundary conditions
in (<ref>), ∂_xu(0,t)=∂_xu(1,t)=0,
instead of the Dirichlet conditions. We show the corresponding results
in Figure <ref>.

We see that while the Parareal algorithm with coarse propagator is
converging as before with Dirichlet conditions, the Parareal algorithm
without coarse propagator is not converging any more, except in the
last case when only 6 coarse time intervals are used, and here it is
the finite step convergence property of the Parareal algorithm when
iterating 6 times that leads to convergence, there is no contraction
without the coarse propagator with Neumann boundary conditions.

It is very helpful at this point to take a closer look at the
solutions we were computing in these two experiments for zero
Dirichlet and zero Neumann conditions, shown in Figure <ref>
in the two leftmost panels.

We see that the solution with Dirichlet boundary conditions is very
much local in time: in order to compute the solution around t=2, we
only need to know if the right hand side function was on just before,
but not at the much earlier times t<<2, since the solution does not
contain this information any more. This is very different in the case
of Neumann conditions, when using the same source term: now the
solution at time t=2 strongly depends on the fact that the right
hand side function was on at earlier times, since the heat injected into
the system is kept, as a constant value in space. In analogy, it is
very easy to predict the temperature in the room you are currently in
a year from now, you only need to know if the heater or air-conditioner
is on and the windows and doors are open or closed then, but not
what happened with the room over the entire year before.



§ CONVERGENCE ANALYSIS


In order to study our observation mathematically, it is easiest to
consider a Parareal algorithm which uses a spectral method in space
and is performing exact integrations in time for our heat equation
model problem (<ref>). The fine propagator
F(U_n^k,T_n,T_n+1) is thus using a highly accurate spectral method
in space with m_F basis functions to solve (<ref>),
and the coarse propagator G(U_n^k,T_n,T_n+1) solves the same
problem using a very cheap spectral method in space using only
m_G<<m_F basis functions. Note that we allow explicitly that
m_G=0, which means the coarse solver is not present in the Parareal
algorithm (<ref>)! 

A spectral representation of the solution of the heat equation
(<ref>) can be obtained using separation of variables,
and considering a spatial domain Ω:=(0,π) to avoid to have to
carry the constant π in all the computations (see also Remark
<ref> later), we get

    u(x,t):=∑_m=1^∞û_m(t)sin(m x).

Expanding also the right hand side in a sine series,

    f(x,t):=∑_m=1^∞f̂_m(t)sin(m x),

and the initial condition,

    u_0(x):=∑_m=1^∞û_0,m(t)sin(m x),

we find that the Fourier coefficients in the solution satisfy

    ∂_tû_m=-m^2û_m+f̂_m,  û_m(0)=û_0,m.

The solution of this equation can be written in closed form using an
integrating factor,

    û_m(t)=û_0,me^-m^2t
        +∫_0^tf̂_m(τ)e^-m^2(t-τ)dτ.

If we use in the Parareal algorithm (<ref>) a spectral
approximation of the solution for the fine propagator using m_F
modes, we obtain 

    F(U_n^k,T_n,T_n+1)=∑_m=1^m_F(Û_m,n^k
      e^-m^2Δ T+∫_T_n^T_n+1f̂_m(τ)
      e^-m^2(t-τ)dτ)sin(m x),

and similarly for the coarse approximation when using m_G modes

    G(U_n^k,T_n,T_n+1)=∑_m=1^m_G(Û_m,n^k
      e^-m^2Δ T+∫_T_n^T_n+1f̂_m(τ)
      e^-m^2(t-τ)dτ)sin(m x).

Because of the orthogonality of the sine functions, the Parareal
algorithm (<ref>) is diagonalized in this representation,
which is the main reason that we chose a spectral method, to simplify
the present analysis for this short note. For Fourier modes with m≤
m_G, the first and the third term in (<ref>) coincide, as we
see from (<ref>) and (<ref>),
and thus cancel, and hence the Parareal update formula simplifies to

    [c][0.85]
            Û_m,n+1^k+1=Û_m,n^k+1
            e^-m^2Δ T+∫_T_n^T_n+1f̂_m(τ)
            e^-m^2(t-τ)dτ,

For m>m_G, the coarse correction term is not present in the Parareal update
formula (<ref>), and thus only the contribution from the fine propagator remains,

    [c][0.85]
            Û_m,n+1^k+1=Û_m,n^k
            e^-m^2Δ T+∫_T_n^T_n+1f̂_m(τ)
            e^-m^2(t-τ)dτ,

The only difference is the iteration index k in the first term on
the right, but this term makes all the difference: for m≤ m_G,
the update formula (<ref>) represent the exact
integration of the corresponding mode sequentially going through the
entire time domain (0,T). For m> m_G, the update formula
(<ref>) represents simply a block Jacobi iteration
solving all time subintervals (T_n,T_n+1) in parallel starting
from the current approximation at iteration k. Using this insight,
we obtain the following convergence result.


  The Parareal algorithm (<ref>) for the heat equation
  (<ref>) using the spectral coarse propagator
  (<ref>) and the spectral fine propagator
  (<ref>) satisfies for any initial guess U_n^0
  the convergence estimate
  
    sup_n||U_n^k(·)-u(·,T_n)||_2≤
        e^-(m_g+1)^2kΔ Tsup_n||U_n^0(·)-u(·,T_n)||_2,

  and this estimate also holds if the coarse propagator does not
  contain any modes, m_G=0, which means it is not present.
  The Parareal algorithm therefore converges also without coarse
  propagator.

 We introduce the error in Fourier space at the
interfaces, Ê_m,n^k:=Û_m,n-Û_m,n^k,
where the converged solution Û_m,n satisfies

    Û_m,n+1=Û_m,n
      e^-m^2Δ T+∫_T_n^T_n+1f̂_m(τ)
      e^-m^2(t-τ)dτ,

Taking the difference with the Parareal update formula for small m
in (<ref>) and using the fact that
U_m,0^k=U_m,0=u^0, we thus obtain a vanishing error in the
first Fourier coefficients after having performed one Parareal
iteration,

    Ê_m,n^k=0,   m=1,2,…,m_G,   k=1,2,….

For m>m_G, taking the difference between (<ref>)
and the update formula for larger m in (<ref>), the
part with the source term f cancels, and we are left with

    Ê_m,n^k+1=Ê_m,n^ke^-m^2Δ T,  
       m=m_G+1,m_G+2,….

Now using the Parseval-Plancherel identity, and taking the largest
convergence factor e^-(m_G+1)^2Δ T out of the sum, we obtain
the convergence estimate (<ref>) taking the sup
over all time intervals n.

We show in Figure <ref> a graphical comparison of the
convergence factor of each Fourier mode without coarse propagator from
(<ref>), i.e. e^-m^2Δ T, and with coarse
propagator |e^-m^2Δ T-R_G(-m^2Δ
  T)|/1-|R_G(-m^2Δ T)| from <cit.> for one
Backward Euler step, R_G(z)=1/1-z, as we used for Figures
<ref> and
<ref>.

We see that indeed there are situations where the coarse correction
with one Backward Euler step is detrimental, as observed in Figures
<ref> and
<ref>, and there is also quantitative
agreement.


  Theorem <ref> shows that the Parareal algorithm (<ref>)
  for the heat equation is scalable, provided the time interval length
  Δ T is held constant: no matter how many such time intervals
  are solved simultaneously, the convergence rate remains the same,
  and this even in the case without coarse propagator! This is like
  the scalability of one level Schwarz methods discovered in
  <cit.> for molecular simulations in solvent models,
  and rigorously proved in
  <cit.>
  using three different techniques: Fourier analysis, maximum
  principle, and projection arguments in Hilbert spaces.



  If we had zero Neumann boundary conditions, a similar result could
  also be derived using a cosine expansion, but then no coarse
  correction would not lead to a convergent method, since the zero
  mode would have the contraction factor one. One must therefore in
  this case have at least the constant mode in the coarse correction
  for the Parareal algorithm to converge. Otherwise however all
  remains the same.



  A precise dependence on the domain size can be obtained by
  introducing the spatial domain of length L, Ω=(0,L). In
  this case the convergence factor simply becomes
  e^-((m_G+1)π/L)^2Δ T. The same
  estimate can also be obtained for a more general parabolic equation,
  the convergence factor then becomes e^-λ_m_G+1Δ T,
  where λ_m_G+1 is the m_G plus first eigenvalue of the
  corresponding spatial operator.




§ CONCLUSIONS


We have shown that Parareal without a coarse propagator can converge
very well for the prototype parabolic problem of the heat equation,
because such problems “forget” all fine information over time, and
only very coarse information, if at all, remains, depending on the
boundary conditions. This is very different for hyperbolic problems:
while for the advection equation u_t+au_x=f with Dirichlet boundary
condition, whose solution is shown in Figure <ref> (middle)
the information is also “forgotten” by leaving the domain on the
right and Parareal converges very well <cit.>, also
without coarse propagator, for the advection equation with periodic
boundary conditions in the fourth plot of Figure <ref> this is
not the case any more: nothing is forgotten, and the coarse propagator
needs to have fine accuracy quality for the Parareal algorithm not to
fail. The situation is similar for the second order wave equation
u_tt=u_xx, whose solution is shown on the right in Figure
<ref>.



  
Acknowledgement. Funded by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) under Germany's Excellence Strategy
EXC 2044 – 390685587, Mathematics Münster: Dynamics–Geometry–Structure
and the Swiss National Science Foundation.


plain


