Parareal

Parareal is an algorithm from numerical analysis and used for the solution of initial value problems. In contrast to e.g. Runge-Kutta or multi-step methods, some of the computations in Parareal can be performed in parallel and Parareal is therefore one example of a parallel-in-time integration method. While historically most efforts to parallelize the numerical solution of partial differential equations focussed on the spatial discretization, in view of the challenges from exascale computing, parallel methods for temporal discretization have been identified as a possible way to increase concurrency in numerical software. Because Parareal computes the numerical solution for multiple time steps in parallel, it is categorized as a parallel across the steps method. This is in contrast to approaches using parallelism across the method like parallel Runge-Kutta method where independent stages can be computed in parallel or parallel across the system methods like waveform relaxation.

History

Parareal can be derived as both a multigrid method in time method or as multiple shooting along the time axis. Both ideas, multigrid in time as well as adopting multiple shooting for time integration, go back to the 1980s and 1990s. Parareal is a widely studied method and has been used and modified for a range of different applications. Ideas to parallelize the solution of initial value problems go back even further: the first paper proposing a parallel-in-time integration method appeared in 1964.

Algorithm

Parareal solves an initial value problem of the form

${\dot {y}}(t)=f(y(t),t),\quad y(t_{0})=y_{0}\quad {\text{with}}\quad t_{0}\leq t\leq T.$

Here, the right hand side $f$ can correspond to the spatial discretization of a partial differential equation in a method of lines approach.

Parareal now requires a decomposition of the time interval $[t_{0},T]$ into $P$ so-called time slices $[t_{j},t_{j+1}]$ such that

$[t_{0},T]=[t_{0},t_{1}]\cup [t_{1},t_{2}]\cup \ldots \cup [t_{P-1},t_{P}].$

Each time slice is assigned to one processing unit when parallelizing the algorithm, so that $P$ is equal to the number of processing units used for Parareal: in an MPI based code for example, this would be the number of processes, while in an OpenMP based code, $P$ would be equal to the number of threads.

Parareal is based on the iterative application of two methods for integration of ordinary differential equations. One, commonly labeled ${\mathcal {F}}$, should be of high accuracy and computational cost while the other, typically labelled ${\mathcal {G}}$, must be computationally cheap but can be much less accurate. Typically, some form of Runge-Kutta method is chosen for both coarse and fine integrator, where ${\mathcal {G}}$ might be of lower order and use a larger time step than ${\mathcal {F}}$. If the initial value problem stems from the discretization of a PDE, ${\mathcal {G}}$ can also use a coarser spatial discretization. The result of numerical integration with one of these methods over a time slice $[t_{j},t_{j+1}]$ for some starting value $y_{j}$ given at $t_{j}$ is then written as

$y={\mathcal {F}}(y_{j},t_{j},t_{j+1})$ or $y={\mathcal {G}}(y_{j},t_{j},t_{j+1})$.

Serial time integration with the fine method would then correspond to a step by step computation of

$y_{j+1}={\mathcal {F}}(y_{j},t_{j},t_{j+1}),\quad j=0,\ldots ,P-1.$

Parareal instead uses the following iteration

$y_{j+1}^{k+1}={\mathcal {G}}(y_{j}^{k+1},t_{j},t_{j+1})+{\mathcal {F}}(y_{j}^{k},t_{j},t_{j+1})-{\mathcal {G}}(y_{j}^{k},t_{j},t_{j+1}),\quad j=0,\ldots ,P-1,\quad k=0,\ldots ,K-1,$

where $k$ is the iteration counter. As the iteration converges and $y_{j}^{k+1}-y_{j}^{k}\to 0$, the terms from the coarse method cancel out and Parareal reproduces the solution that is obtained by the serial execution of the fine method only. It can be shown that Parareal converges after a maximum of $P$ iterations.. For Parareal to provide speedup, however, it has to convergence in a number of iterations significantly smaller than the number of time slices, that is $K\ll P$.

In the Parareal iteration, the computationally expensive evaluation of ${\mathcal {F}}(y_{j}^{k},t_{j},t_{j+1})$ can be performed in parallel on $P$ processing units. By contrast, the dependency of $y_{j+1}^{k+1}$ on ${\mathcal {G}}(y_{j}^{k+1},t_{j},t_{j+1})$ means that the coarse correction has to be computed in serial order.

Speedup

Under some assumptions, a simple theoretical model for the speedup of Parareal can be derived. Although in applications these assumptions can be too restrictive, the model still is useful to illustrate the trade offs that are involved in obtaining speedup with Parareal.

First, assume that every time slice $[t_{j},t_{j+1}]$ consists of exactly $N_{f}$ steps of the fine integrator and of $N_{c}$ steps of the coarse integrator. This includes in particular the assumption that all time slices are of identical length and that both coarse and fine integrator use a constant step size over the full simulation. Second, denote by $\tau _{f}$ and $\tau _{c}$ the computing time required for a single step of the coarse or fine method and assume that both are constant. This is typically not exactly true when an implicit method is used, because then runtimes vary depending on the number of iterations required by the iterative solver.

Under these two assumptions, the runtime for the fine method integrating over $P$ time slices can be modelled as

$c_{\text{fine}}=PN_{f}\tau _{f}.$

The runtime of Parareal using $P$ processing units and performing $K$ iterations is

$c_{\text{parareal}}=(K+1)PN_{c}\tau _{c}+KN_{f}\tau _{f}.$

Speedup of Parareal then is

$S_{p}={\frac {c_{\text{fine}}}{c_{\text{parareal}}}}={\frac {1}{(K+1){\frac {N_{c}}{N_{f}}}{\frac {\tau _{c}}{\tau _{f}}}+{\frac {K}{P}}}}\leq \min \left\{{\frac {N_{f}\tau _{f}}{N_{c}\tau _{c}}},{\frac {P}{K}}\right\}.$

These two bounds illustrate the trade off that has to be made in choosing the coarse method: one the one hand, it has to be cheap and/or use a much larger time step to make the first bound as large as possible, on the other hand the number of iterations $K$ has to be kept low to keep the second bound large. In particular, Parareal's parallel efficiency is bounded by

$E_{p}={\frac {S_{p}}{P}}\leq {\frac {1}{K}},$

that is by the inverse of the number of required iterations.