Parareal
Parareal is a parallel algorithm from numerical analysis and used for the solution of initial value problems. It was introduced in 2001 by Lions, Maday and Turinici. Since then, it has become one of the most widely studied parallel-in-time integration methods. Parallel-in-time integration methods. 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 or extrapolation methods, 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 formula_1 Here, the right hand side formula_2 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 formula_3 into formula_4 so-called time slices formula_5 such that formula_6 Each time slice is assigned to one processing unit when parallelizing the algorithm, so that formula_4 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, formula_4 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 labelled formula_9 , should be of high accuracy and computational cost while the other, typically labelled formula_10, 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 formula_10 might be of lower order and use a larger time step than formula_9. If the initial value problem stems from the discretization of a PDE, formula_10 can also use a coarser spatial discretization, but this can negatively impact convergence unless high order interpolation is used. The result of numerical integration with one of these methods over a time slice  formula_14 for some starting value formula_15 given at formula_16 is then written as formula_17 or formula_18. Serial time integration with the fine method would then correspond to a step-by-step computation of formula_19 Parareal instead uses the following iteration formula_20 where formula_21 is the iteration counter. As the iteration converges and formula_22, 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 formula_4 iterations. For Parareal to provide speedup, however, it has to  converge in a number of iterations significantly smaller than the number of time slices, that is formula_24. In the Parareal iteration, the computationally expensive evaluation of formula_25 can be performed in parallel on formula_4 processing units. By contrast, the dependency of formula_27 on formula_28 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 formula_5 consists of exactly formula_30 steps of the fine integrator and of formula_31 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 formula_32 and formula_33 the computing time required for a single step of the fine and coarse methods, respectively, 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 formula_4 time slices can be modelled as formula_35 The runtime of Parareal using formula_4 processing units and performing formula_37 iterations is formula_38 Speedup of Parareal then is formula_39 These two bounds illustrate the trade off that has to be made in choosing the coarse method: on 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 formula_37 has to  be kept low to keep the second bound large. In particular, Parareal's parallel efficiency is bounded by formula_41 that is by the inverse of the number of required iterations. Instability for imaginary eigenvalues. The vanilla version of Parareal has issues for problems with imaginary eigenvalues. It typically only converges toward the very last iterations, that is as formula_21 approaches formula_4, and the speedup formula_44 is always going to be smaller than one. So either the number of iterations is small and Parareal is unstable or, if formula_21 is large enough to make Parareal stable, no speedup is possible. This  also means that Parareal is typically unstable for hyperbolic equations. Even though the formal analysis by Gander and Vandewalle covers only linear problems with constant coefficients, the problem also arises when Parareal is applied to the nonlinear Navier–Stokes equations when the viscosity coefficient becomes too small and the Reynolds number too large. Different approaches exist to stabilise Parareal, one being Krylov-subspace enhanced Parareal. Variants. There are multiple algorithms that are directly based or at least inspired by the original Parareal algorithm. Krylov-subspace enhanced Parareal. Early on it was recognised that for linear problems information generated by the  fine method formula_46 can be used to improve the accuracy of the coarse method formula_47. Originally, the idea was formulated for the parallel implicit time-integrator PITA, a method closely related to Parareal but with small differences in how the correction is done. In every iteration formula_21 the result formula_49 is computed for values formula_50 for formula_51. Based on this information, the subspace formula_52 is defined and updated after every Parareal iteration. Denote as formula_53 the orthogonal projection from formula_54 to formula_55. Then, replace the coarse method with the improved integrator formula_56. As the number of iterations increases, the space  formula_55 will grow and the modified propagator formula_58 will become more accurate. This will lead to faster convergence. This version of Parareal can also stably integrate linear hyperbolic partial differential equations. An extension to nonlinear problems based on the reduced basis method exists as well. Hybrid Parareal spectral deferred corrections. A method with improved parallel efficiency based on a combination of Parareal with spectral deferred corrections (SDC) has been proposed by M. Minion. It limits the choice for coarse and fine integrator to SDC, sacrificing flexibility for improved parallel efficiency. Instead of the limit of formula_59, the bound on parallel  efficiency in the hybrid method becomes formula_60 with formula_61 being the number of iterations of the serial SDC base method and formula_62 the typically greater number of iterations of the parallel hybrid method. The Parareal-SDC hybrid has been further improved by addition of a "full approximation scheme" as used in nonlinear multigrid. This led to the development of the "parallel full approximation scheme in space and time" (PFASST). Performance of PFASST has been studied for PEPC, a Barnes-Hut tree code based particle solver developed at Juelich Supercomputing Centre. Simulations using all 262,144 cores on the IBM BlueGene/P  system JUGENE showed that PFASST could produce additional speedup beyond saturation of the spatial tree parallelisation. Multigrid reduction in time (MGRIT). The multigrid reduction in time method (MGRIT) generalises the interpretation of Parareal as a multigrid-in-time algorithms to multiple levels using different smoothers. It is a more general approach but for a specific choice of parameters it is equivalent to Parareal. The XBraid library implementing MGRIT is being developed by Lawrence Livermore National Laboratory. ParaExp. ParaExp uses exponential integrators within Parareal. While limited to linear problems, it can produce almost optimal parallel speedup.