An adaptive parareal algorithm.

In this paper, we consider the problem of accelerating the numerical simulation of time dependent problems by time domain decomposition. The available algorithms enabling such decompositions present severe efficiency limitations and are an obstacle for the solution of large scale and high dimensional problems. Our main contribution is the improvement of the parallel efficiency of the parareal in time method. The parareal method is based on combining predictions made by a numerically inexpensive solver (with coarse physics and/or coarse resolution) with corrections coming from an expensive solver (with high-fidelity physics and high resolution). At convergence, the algorithm provides a solution that has the fine solver's high-fidelity physics and high resolution. In the classical version, the fine solver has a fixed high accuracy which is the major obstacle to achieve a competitive parallel efficiency. In this paper, we develop an adaptive variant that overcomes this obstacle by dynamically increasing the accuracy of the fine solver across the parareal iterations. We theoretically show that the parallel efficiency becomes very competitive in the ideal case where the cost of the coarse solver is small, thus proving that the only remaining factors impeding full scalability become the cost of the coarse solver and communication time. The developed theory has also the merit of setting a general framework to understand the success of several extensions of parareal based on iteratively improving the quality of the fine solver and re-using information from previous parareal steps. We illustrate the actual performance of the method in stiff ODEs, which are a challenging family of problems since the only mechanism for adaptivity is time and efficiency is affected by the cost of the coarse solver.

Introduction

Solving complex models with high accuracy and within a reasonable computing time has motivated the search for numerical schemes that exploit efficiently parallel computing architectures. In this paper, the model consists of a Partial Differential Equation (PDE) set on a domain . In this context, one of the main ideas to parallelize a simulation is to break the problem into subproblems defined over subdomains of a partition of . The domain can potentially have high dimensionality and be composed of different variables like space, time, velocity or even more specific variables for some problems. While there exist algorithms with very good scalability properties for the decomposition of the spatial variable in elliptic and saddle-point problems (see [1] or [2] for an overview), the same cannot be said for the decomposition of time of even simple systems of ODEs. This is despite the fact that research on time domain decomposition is currently very active and has by now a history of at least 50 years (back to at least [3]) during which several algorithms have been explored (see [4] for an overview). As a consequence, time domain decomposition is to date only a secondary option when it comes to deciding what algorithm/method distributes the tasks in a parallel cluster.

The main goal of this work is to address this efficiency limitation in the framework of one particular scheme: the parareal in time algorithm. The method was first introduced in [5] and has been well accepted by the community because it is easily applicable to a relatively large spectrum of problems. (Some specific difficulties are nevertheless encountered on certain types of PDEs as reported in, e.g., [6], [7] for hyperbolic systems or [8], [9] for hamiltonian problems). Another ingredient for its success is that, even though its scalability properties are limited, they are in general competitive in comparison with other methods. Without entering into very specific details of the algorithm at this stage, we can summarize the procedure by saying that we build iteratively a sequence to approximate the exact solution of the problem by a predictor–corrector algorithm. At every iteration, predictions are made by a solver which has to be as numerically inexpensive as possible since it is run on the full time interval. It usually involves coarse physics and/or coarse resolution. Corrections involve an expensive solver with high-fidelity physics and high resolution which is propagated in parallel over small time subdomains. In the classical version of parareal, the fine solver has a fixed high accuracy across all iterations. It is set to the one that we would use to solve the dynamics at the desired accuracy with a purely sequential solver. It is well-known that this point is the major obstacle to achieve better parallel efficiency. In this paper, we propose an adaptive variant where the accuracy of the fine solver is increased across the iterations. Our main goal is to show that this new point of view overcomes the obstacle of the cost of the fine solver and that the only remaining factors limiting high performance become the cost of the coarse solver and communication time. We refer to, e.g., [10] for contributions on the lowering of the cost of that coarse solver.

We present in Section 2 the new adaptive point of view. This requires to formulate an idealized version of the parareal algorithm in an infinite dimensional function space where the fine propagations are replaced by the exact ones (Section 2.2). Since this scheme is obviously not implementable in practice, we formulate a feasible “perturbed” version that involves approximations of the exact propagations at increasing accuracy across the iterations (Section 2.3). The accuracies are tightened in such a way that the feasible adaptive algorithm converges at the same rate as the ideal one and with a near-minimal numerical cost. The identified tolerances involve quantities that are difficult to estimate in practice. In addition, they may not be optimal because they are derived from a theoretical convergence analysis based on abstract conditions for the coarse and fine solvers. We bridge this gap between theory and actual implementation by proposing practical guidelines to set these tolerances. We next explain in Section 2.4 how the new formulation invites to use adaptive schemes not only in the time variable, but also in other variables that may be involved in the dynamics. The performance of the algorithm could also be enhanced by re-using informations from previous iterations in order to limit the cost of internal solvers. The techniques for this will strongly depend on the nature of the specific problem. We discuss common situations in Appendix. We close Section 2.5 by listing the main advantages of the new framework and how the classical parareal paradigm can be formulated with the optics of the new standpoint.

The parallel performance of the adaptive scheme is difficult to predict a priori but in Section 3 we carry a discussion where we show that it will always be superior to the classical approach. In the idealized situation where the cost of the coarse solver and communication delays are negligible, we show that the algorithm would exhibit a very high parallel efficiency.

We emphasize that our theory is general in the sense that it is applicable to ODEs and also to PDEs involving time, space and possibly other variables. We defer to a future work the presentation of a numerical PDE example since it requires the deployment of space–time adaptive methods which is a challenging topic in itself and the techniques usually depend very specifically on the problem nature. Instead, we illustrate the performance of the algorithm on stiff ODEs in Section 4. They are a challenging family of problems because the only source of adaptivity is time and they do not allow to use a very inexpensive coarse solver. The tested ODEs are the Brusselator, the Van der Pol, the Oregonator equations, and an SEIR model which has very recently been proposed in [11] to model the spread of the COVID-19 virus in the Wuhan city area. The code to reproduce our results (and experiment with other ODEs) is available online.1 Our first two examples are relatively stiff, while the other two are highly stiff, and serve to illustrate the limitations of the approach. We show that in the relatively stiff problems, the adaptive parareal algorithm performs between 2 to 3 times better than the classical one when the solution is approximated at high accuracy. In addition, we confirm in these two examples that the only remaining obstacles to achieve a very competitive performance are the cost of the coarse solver and communication time between processors. Due to the nature of the algorithm, it is not clear how to overcome these limitations, especially the one coming from the coarse solver. This issue will come at the forefront for future research since, as our two highly stiff examples illustrate, the cost of the coarse solver may even prevent from obtaining any speed-up at all. We show that if we could find an inexpensive coarse solver, perhaps based on empirical data or on good back-of the envelope calculations, our adaptive algorithm would yield interesting speeds-ups even in highly stiff cases.

In addition to the improvement in parallel efficiency (except, of course, in the above discussed extreme cases), the adaptive version of parareal brings an important conceptual novelty to the field of time domain decomposition which is the one of error controlled computations. By this we mean that the deviation of the numerical result from the exact continuous solution is certifiably quantified and set to meet a given target accuracy with respect to a problem relevant norm. This requires the formulation of the algorithm at an infinite dimensional level as is done in this paper. This point of view is fundamentally different from the fully discrete setting in which the parareal algorithm has always been thought of in practice. That is, we first fix a discretization (often on a uniform grid) and then speed-up the computation of the discrete evolution with parareal. In this way, we do not have any rigorous control of the error with respect to the actual continuous solution and we do not have any systematic procedure to minimize the number of degrees of freedom.

We conclude this introduction by some bibliographical remarks. To the best of our knowledge, the current abstract and broad formulation of an adaptive version of parareal has never been proposed in the literature. However, previous works have instantiated in a variety of particular applications the idea of re-using information from previous parareal iterations, either with the purpose of improving the quality of the fine solver or to build an initial guess of internal iterative routines. Among the most relevant ones stand the coupling of the parareal algorithm with spatial domain decomposition (see[12], [13], [14]), the combination of the parareal algorithm with iterative high order methods in time like spectral deferred corrections (see [15], [16], [17]). Parareal has also been combined with multigrid iterative techniques which, in addition, involve a hierarchy of space–time meshes (see [18], [19]). In a similar spirit, we can also cite applications of the parareal algorithm to solve optimal control problems where information from previous steps is used (see [12], [20]). In Appendix, we briefly explain in what sense the above strategies can be seen as particular instances of the current approach and how our viewpoint could help to give them more solid theoretical foundations. The idea of re-using information as a starting guess for internal iterative solvers is also discussed in the Appendix. It has been explored in several different works, e.g. [17], [18], and [21] provides a convergence analysis in simple situations (a more complete analysis will be proposed in a forthcoming work).

An adaptive parareal algorithm

In this section, after introducing some preliminary notations in Section 2.1, we formulate an ideal parareal scheme on an infinite dimensional functional setting (Section 2.2). We then present feasible realizations involving a fine solver whose accuracy is adaptively increased across the iterations (Section 2.3). We prove that the feasible adaptive algorithm converges at the same rate as the ideal one provided that the tolerances of the fine solver are increased at certain rate which will be discussed. Finally, we discuss how the new paradigm can be realized thanks to adaptive schemes and/or the re-use of information from previous steps (Section 2.4 and Appendix). In Section 2.5, we connect the new adaptive formulation with the classical parareal algorithm and list the main advantages of the new standpoint.

Setting and preliminary notations

Let be a Banach space of functions defined over a domain (), e.g. . Let be a propagator, that is, an operator such that, for any given time , and any function , takes as an initial value at time and propagates it at time . We assume that satisfies the semi group property We further assume that is implicitly defined through the solution of the time-dependent problem where is an operator from into with adequate regularity we shall detail latter. Then, given any , denotes the solution to (1) at time with initial condition at time . In our problem of interest, we study the evolution given by (1) when the initial condition is . Note that could also be associated to a discretized version of the evolution equation or be defined through an operator that is not necessary related to an evolution equation (see [22]).

Since, in general, the problem does not have an explicit solution, we seek to approximate it at a given target accuracy. For any initial value , any , and any we denote by an element of that approximates such that we have where, here and in the following, denotes the norm in . Any realization of involves three main ingredients:

a numerical scheme to discretize the time dependent problem (1) (e.g. an Euler scheme in time),

a certain expected error size associated with the choice of the discretization (e.g. error associated with the time step size of the Euler scheme),

a numerical implementation to solve the resulting discrete systems (e.g. conjugate gradient, Newton method, SSOR, …).

In the following, we will use the term solver to denote a particular choice for (i), (ii) and (iii). Given a solver , we will use the same notation as for the exact propagator to express that is an approximation of with a certain accuracy . In other words, we can write .

An idealized version of the parareal algorithm

We introduce a decomposition of the time interval into subintervals , . Without loss of generality, we will take them of uniform size which means that for . For a given target accuracy , the primary goal of the parareal in time algorithm is to build an approximation of such that The classical way to achieve this is to set where is some sequential solver in with in (2). Since this comes at the cost of solving over the whole time interval , the main goal of the parareal in time algorithm is to speed up the computing time, while maintaining the same target accuracy . This is made possible by first decomposing the computations over the time domain. Instead of solving over , we perform parallel solves over each interval of size . We next introduce an idealized version of it which will not be feasible in practice but will be the starting point of subsequent implementable versions. The algorithm relies on the use of a solver (known as the coarse solver) with the following properties involving the operator

Hypotheses (H): There exists constants such that for any function and for any and , Note that these hypothesis are the classical abstract formulations of the properties of numerical schemes related to stability and accuracy. Hypothesis (4b) is a Lipschitz condition and the quantity is a small constant which, in the case of a Euler scheme, would be equal to the time step size.

The idealized version of the algorithm consists in building iteratively a series of approximations of for following the recursive formula At this point, several comments are in order. The first one is that the computation of only requires propagations with over intervals of size . As follows from (5), for a given iteration , propagations of this size are required, each of them over distinct intervals of size , each of them with independent initial conditions. Since they are independent from each other, they can be computed over parallel processors and the original computation over is decomposed into parallel computations over subintervals of size . The second observation is that the algorithm may not be implementable in practice because it involves the exact propagator . Feasible instantiations consist of replacing by some approximation with a certain accuracy which has to be carefully chosen. We will come to this point in the next section. The third observation is to note that, in the current version of the algorithm, for all , the exact solution is obtained after exactly parareal iterations. This number can be reduced when we only look for an approximate solution with accuracy . Depending on the problem, the final number of iterations can actually be much smaller than . The convergence result of Theorem 2.1 and its proof are helpful to understand the main mechanisms driving the convergence of the algorithm and explaining its behavior. To present it, we introduce the shorthand notation for the error norm and the quantities

If and satisfy Hypothesis (4) , then,

The proof is in the spirit of existing results from the literature (see [5], [23], [24], [25]) but it is instructive to give it for subsequent developments in the paper. We introduce the following quantities as shorthand notations for the proof.

If , using definition (5) for , we have for , where we have used (4a), (4b) to derive the second to last inequality.

For , starting from (5), we have Taking norms and using (4b), (4c), we derive Following [25], we consider the sequence defined recursively as follows. For , and for , Since for and , we analyze the behavior of to derive a bound for . For this, we consider the generating function From (8), (9) we get from which we derive Since, , we can bound the term in the denominator by . Next, using the binomial expansion and identifying the term in in the expansion, we derive the bound Hence, using definition (7) for and , which ends the proof of the theorem. □

Note that at least one step is not sharp in the above proof: it is the step where is replaced by . Note also that is the quantity driving convergence and its speed.

Introducing the quantity we can write and we note that a sufficient condition to converge is that In other words, is the minimal accuracy that the coarse solver has to satisfy in order to guarantee convergence of the ideal parareal algorithm. In the following, we will work under the assumption that satisfies (11).

As we will see in the next section, plays also a critical role in certain convergence properties of the perturbed algorithm so we finish this section by discussing the behavior of depending on several scenarios. First, and are Lipschitz constants (fixed by the properties of the evolution problem) so they could be potentially large numbers. As a result, could be a large number and condition (11) would not be very stringent. The value of can be small for very long time simulations where becomes large or if becomes small compared to (that is, if the number of processors becomes large).

Feasible realizations of the parareal algorithm

Feasible versions of algorithm (5) involve approximations of with a certain accuracy . This leads to consider algorithms of the form Since no feasible version will converge at a better rate than (6), we analyze here what is the minimal accuracy that preserves it. A result in this direction is given in the following theorem. It requires to introduce the quantity which tends to as .

Let and satisfy Hypothesis (4) . Let be any given positive integer. If for all and all , the approximation has accuracy then the of the feasible parareal scheme (12) satisfy with

Let us make a couple of remarks before giving the proof of the theorem. First, the sufficient condition to converge is now so the minimal accuracy required for the coarse solver is stronger than in (11) for the ideal case. Note however that when is small (roughly, ), the condition on is similar in the ideal and perturbed case.

Second, comparing (6), (14), the rate of convergence of the feasible parareal algorithm deviates from , the ideal one, by a factor The parameter plays again a critical role in the convergence properties and determines whether convergence is close to the ideal rate , or deviates from it by a potentially important factor.

The proof follows the same lines as the one for Theorem 2.1 and , are defined exactly as before. In addition, it will be useful to introduce the sequence We concentrate on the case since the case is identical as in Theorem 2.1. For , using (12), we have Taking norms, using (4b), (4c) and the definition (2) applied to , we derive Similarly to Theorem 2.1, we introduce the sequence defined for as for all and for , The associated generating function satisfies Hence By replacing again at the denominator the factor by and using the binomial expansion (10), we derive the bound where we have used that . The coefficient associated to the term above gives the inequality From the definition of , we have that Therefore, recalling the definition (7) of and , we derive where we have used the definition of and in the last line. This inequality ends the proof since for . □

Practical realization of

Since the accuracy needs to improve with , the most natural way to build the approximations is with adaptive techniques and with adaptive refinements at every step . The implementation ultimately rests on the use of a posteriori error estimators. It opens the door to local time step adaptation in the parareal algorithm as well as spatial coarsening or refinement if the problem involves additional spatial variables.

In principle, as decreases with , the numerical cost increases in terms of degrees of freedom and also in terms of computing time. This actually reveals the key idea of this new approach which is that we would like that only the last fine solver is expensive and the cost of the previous ones is a small fraction of the cost of the last one (we refer to the next sub-section for a more precise statement). By re-using information from previous iterations, we can limit the cost of internal solvers required in and enhance the speed-up. This depends of course on the nature of the specific problem. We discuss several common situations in Appendix.

Connection to the classical formulation of the parareal algorithm and advantages of the current view-point

In the original version of the algorithm, is approximated with an accuracy which is kept constant in and across the parareal iterations . This has usually been done by using a solver defined in the same spirit as , but satisfying Hypothesis (4) with a better accuracy . We have in this case and we recover the classical algorithm (see [26] and [27]) Compared to this classical version of the parareal algorithm, the adaptive approach offers the following important advantages:

The algorithm converges to the exact solution and not to the solution achieved by the fixed chosen fine solver (indeed, for any , as ).

For a final target accuracy , the parallel efficiency will always be superior to the classical approach (see Section 3).

We minimize the computational resources (degrees of freedom) because we identify the minimal required accuracies at each iteration (see Eq. (13)). Early iterations use a loose tolerance, thus avoiding unnecessary work due to oversolving, while later iterations use tighter tolerances to deliver accuracy.

The dynamical refinements of the fine solver invite to incorporate adaptive solvers with a posteriori error estimators to the parareal scheme.

Parallel efficiency

It is difficult to give accurate a priori estimations for the speed-up and efficiency of the method due to its adaptive nature so the actual performance can only be established through relevant examples. In Section 4, we give some results for the case of the Brusselator system. Despite this difficulty in estimation, we make some general remarks in this section, which aim primarily at highlighting the relevance of the cost of the coarse solver. The speed-up is defined as the ratio between the cost to run a sequential fine solver achieving a target accuracy with the cost to run an adaptive parareal algorithm providing at the end the same target accuracy . The parallel efficiency of the method is then defined as the ratio of the above speed up with the number of processor which gives a target of 1 to any parallel solver:

Assume that and are the numerical costs to realize and . Since the tolerances decrease with , we have . Neglecting the communication delays, the cost of the adaptive solver is The classical parareal algorithm involves, at every iteration propagations at the highest accuracy . Thus its cost is from which it directly follows that (at least if, with obvious notation, ) Thus the parallel performance of the adaptive algorithm is at least the one of the classical version. Note that this holds even when communications are not negligible since there is the same amount of information exchange in both algorithms.

We next give a more quantitative statement on an admittedly idealized setting. Assume that the cost of the coarse solve is negligible, that there is no communication delay and that the cost to realize is This assumption for the cost is, for instance, reasonable when we use an explicit time-stepping method of order . It would also hold for an implicit method where a direct solver can be used. Note that could actually depend on but we stick to this simple model for clarity of exposition.

If , for some and if the cost of the coarse solver is negligible with respect to for any , then Therefore

The cost of the scalable adaptive parareal scheme after iterations is

Since we are in a range where the scheme converges, the quantity is bounded and thus there exists such that for all . We will account for this with the notation . Note that in fact and are close to one. Let us start with the simple case and denote : we thus derive In the general case (), the same is true with Up to this last factor, the current conclusion is that the global cost of the parareal procedure is equal to the last fine solver on each sub-interval with size (both the coarse and the previous fine propagations are negligible).

Since the accuracy that is obtained at the end of the parareal procedure (see (14)) is of the same order as the accuracy provided with classical parareal solver (compare with (6)), it follows that if we now take the last target accuracy of the adaptive algorithm as the accuracy of the fine scheme in the classical parareal algorithm, the cost would be Therefore, In addition, we know that when the cost of the coarse solver is negligible, Dividing (20) by (19) yields and

In the ideal setting of Proposition 3.1:

The parallel efficiency of the adaptive parareal algorithm does not depend on the final number of iterations. This is in contrast to the classical version whose efficiency decreases with the final number of iterations as .

The efficiency behaves like in the adaptive version, and rapidly goes to zero with . As soon as becomes negligible with respect to 1, we will be in the range of full scalability.

We emphasize that, obviously, the above idealized setting will never hold in practice, but the result is interesting in its own right since it highlights that the cost of the fine solver is no longer the main obstacle for full scalability in the adaptive setting: the cost of the coarse solver becomes now the major obstruction towards full efficiency.

Numerical tests

Guidelines for a practical implementation

Practical choice of .

Formula (13) of the convergence analysis of Section 2.3 gives an estimate for that one could in principle use for the implementation. However, these tolerances may not be optimal because they are derived from a theoretical convergence analysis based on abstract conditions for the coarse and fine solvers. This was confirmed during our numerical tests where we observed that using estimates (13) for did not deliver satisfactory enough results. This is the reason why it is necessary to devise a practical rule to set . We have explored the following choice: if is the final target accuracy, the classical parareal algorithm is usually run with a solver that delivers a slightly higher accuracy, say . Assume that the classical algorithm converges in iterations. We propose to build the tolerances of in such a way to target that and such that the cost of the last fine propagation is of the order of the sum of the previous ones. This motivates to set The numerical example of the next section uses these tolerances.

Load balancing.

For simplicity of exposition, the algorithm has so far been discussed for subintervals of uniform size . However, this decomposition may lead to a task imbalance because some time intervals may have more complex dynamics than others, requiring more degrees of freedom, thus more computational time. In order to balance tasks as efficiently as possible, we dynamically adapt the size of the subintervals in a way to have the fine solver propagations as balanced as possible among processors.

Results for several stiff ODEs

We apply our adaptive algorithm to several stiff ODEs where the only mechanism for adaptivity is time. Our results illustrate that our approach improves the speed-up and efficiency with respect to the classical non-adaptive parareal method. We also show that the main element affecting performance is no longer the cost of the fine solver but the cost of the coarse solver. In extreme cases, this cost may even prevent any speed-up at all (see Section 4.2.3) and puts this obstruction at the forefront for future research. The code to reproduce the numerical results is available online at: Other ODEs can easily be tested as indicated in the instructions.

https://plmlab.math.cnrs.fr/mulahernandez/parareal-adaptive

Note that the algorithm could also be applied to PDEs but we defer the presentation of numerical examples to future works since this requires full space–time adaptive techniques which are a topic in itself since they are challenging to formulate and deploy and very specific to each type of problem.

The brusselator system

We consider the brusselator system with initial condition and . This is a stiff ODE that models a chain of chemical reactions. It was already studied in a previous work on the parareal algorithm (see [25]). The system has a fixed point at and which becomes unstable when and leads to oscillations. We place ourselves in this oscillatory regime by setting and . The dynamics present large velocity variations in some time subintervals, making the use of adaptive time-stepping schemes particularly desirable for an appropriate treatment of the transient.

For the coarse solver, we set and use an explicit Runge–Kutta method of order 5 with an adaptive time-stepping (see [28]). For the fine solver, we use the implicit Runge–Kutta method of the Radau IIA family of order 5 with adaptive time-stepping (see [29], [30]). Both integrators are available in the ODE integration library of Scipy2 which we have used in our library.

As already discussed, the target accuracies should be ensured by rigorous a posteriori error estimators. However, these type of estimators are unfortunately not available in the Scipy library and we are not aware of any mainstream library with this capability. As a surrogate, we have used the above mentioned classical ODE integrators that only guarantee local accuracy between time-steps , but not global accuracy between macro intervals (composed of several time-steps). The local accuracy can be specified in the library routine via the parameters atol and rtol of the function scipy.integrate.solve_ivp. To relate this local accuracy control to the global one, we have built a priori a “chart” mapping accuracies of the solver on macro-intervals against the tolerance parameters atol and rtol of the library. To simplify, these two parameters have been set to be equal () and their value is fixed according to the chart. As an example, we provide a chart for for the scheme of the fine solver in Fig. 1. The dots are computed values: for a given value of the parameter , we examine the accuracy of the solver. We then interpolate the points with a cubic spline interpolation. This way, for a given intermediate accuracy in the parareal algorithm, we can easily adapt the parameter value that is required.

Fig. 1

Mapping of the accuracies against the tolerance parameters () of the library. The dots are computed values: for a given value of the parameter, we examine the accuracy of the solver. We then interpolate the points with a cubic spline interpolation. This way, for a given intermediate accuracy in the algorithm, we can infer the parameter value and . Case , integrator of the fine solver.

We use formula (16) to compare the speed-up of the classical and adaptive parareal algorithm in terms of the number of operations involved in the numerical solution (communication delays have not been taken into account). For the costs and , we take into account:

the number of time steps (which is adaptively increased as we tightened the accuracy),

the number of right-hand side evaluations,

for the fine solver, we additionally count the number of evaluations of the Jacobian matrix and of the number of linear system inversions.

In Fig. 2, we plot the obtained speed-up for different configurations:

Fig. 2

Speed-up in comparison to running a sequential fine solver as a function of the number of processors . Dashed lines: classical parareal. Continuous lines: Adaptive parareal.

the final time varies from to ,

the final target accuracy is or ,

the number of processors varies from 10 to 100.

As anticipated in Section 3, the speed-up of the adaptive parareal is always superior to the one of the classical parareal. We observe that the gain is marginal for a moderate accuracy () but it is about 2.5 times larger for . Note that sometimes the speed-up does not increase monotonically as the number of processors increases. Also, the speed-up generally increases with but the increase is rather moderate.

The values significantly differ from the range of full scalability and we next explain why this is mainly due to the cost of the coarse solver. Since the problem is stiff and we consider relatively long time intervals, it has been necessary to use a sufficiently accurate coarse solver. This explains our choice of an explicit Runge–Kutta scheme of order 5. To illustrate the impact of its cost, let us fix , and (other parameters would yield similar conclusions). We compare the speed-up and efficiency when we count or do not count the cost of the coarse solver in Table 1. Obviously, when we do not count the cost of the coarse solver, the performance of both algorithms improves but it is particularly increased in the case of the adaptive version. If the cost of was negligible, it would deliver a very satisfactory efficiency of . This is five times larger than what the classical parareal would yield. This analysis illustrates that the major obstacle to achieve competitive scalabilities is no longer the cost of the fine solver like in the classical version, but the cost of the coarse propagator.

Table 1

Brusselator: Impact of the cost of the coarse solver. Speed-up and efficiency with , and .

Speed-up	Classical parareal	Adaptive parareal
With cost G	4.06	7.38
Without cost G	7.38	37.76
Efficiency	Classical parareal	Adaptive parareal
With cost G	8%	14.76%
Without cost G	14.76%	75.52%

We next give some insight on the differences in the convergence behavior of both algorithms. We fix , and and plot in Fig. 3 the convergence history of the parareal solution in terms of:

Fig. 3

Brusselator: Convergence history of the errors for , and . Top: classical parareal. Bottom: adaptive parareal. Left: errors of the fine solver at every fine time-step. Right: maximum parareal error at each iteration .

the errors of the fine solver at every fine time-step

the maximum error of the parareal solution at the macro-intervals

Note that the maximum error in the adaptive scheme steadily decreases to the desired accuracy whereas the error in the classical scheme degrades at iteration before converging. This type of behavior has been observed for all other configurations and we conjecture that an important difference in accuracy between the coarse and the fine solver at early stages of the algorithm may be the cause. Finally, an inspection of the error of the fine solver shows that the adaptive algorithm succeeds to reduce the error at every time in a much more uniform way than the classical algorithm.

The Van der Pol oscillator

We next consider the Van der Pol oscillator with initial condition and . When , this equation is a simple nonstiff harmonic oscillator. When , the system has a limit cycle and becomes stiffer and stiffer as its value is increased. For our tests, we set which is already a relatively stiff case.

Like in the example of the Brusselator system, we set for the coarse solver and use an explicit Runge–Kutta method of order 5 with an adaptive time-stepping (see [28]). For the fine solver, we use the implicit Runge–Kutta method of the Radau IIA family of order 5 with adaptive time-stepping.

In Fig. 4, we plot the obtained speed-up for different configurations:

Fig. 4

Van der Pol: speed-up in comparison to running a sequential fine solver as a function of the number of processors . Dashed lines: classical parareal. Continuous lines: Adaptive parareal.

the final time is or ,

the final target accuracy is or ,

the number of processors varies from 10 to 100.

Like in the previous example, the adaptive algorithm outperforms the nonadaptive version in terms of speed-up. However, the gain is marginal for moderate accuracies . For high accuracy , the adaptive algorithm improves the speed-up by a factor of about 2 to 3 times with respect to the classical one. The improvement is more significant for large .

In Table 2, we illustrate that the coarse solver is again the main bottleneck to reach high parallel efficiency in the adaptive algorithm: we examine the speed-up and efficiency for , and when we take and do not take into account the cost of the coarse solver.

Table 2

Van der Pol: Impact of the cost of the coarse solver. Speed-up and efficiency with , and .

Speed-up	Classical parareal	Adaptive parareal
With cost G	4.54	11.14
Without cost G	6.61	32.63
Efficiency	Classical parareal	Adaptive parareal
With cost G	11.35%	27.8%
Without cost G	16.5%	81.56%

Discussion on extremely challenging cases of highly stiff ODEs: Oregonator and SEIR epidemic model

In extreme cases, the cost of the coarse solver may prevent any speed-up at all (see Section 4.2.3) and puts this obstruction at the forefront for future research. A prominent case are highly stiff ODEs. In the provided code, we can observe this fact in the case of the Oregonator system of equations. In view of the ongoing pandemic of COVID-19 at the time when this article was written, we have also tested an SEIR epidemic model which has very recently been proposed in [11] as a simple model for the spread of the virus in the Wuhan city area. The model captures the effect of the presence of individual reaction to the risk of infection and governmental action. If we run the parareal algorithm for the latter model, an integrator of type LSODA [31] for both coarse and fine solvers seems well adapted since it alternates between Adams or BDF integration for nonstiff and stiff parts and it has automatic stiffness detection. To reach a target accuracy , it has been necessary to set , making the cost of the coarse solver too expensive to yield any parallel efficiency. However, if we could find an inexpensive coarse solver, perhaps based on empirical data or on good back-of the envelope calculations, our adaptive parareal algorithm would yield interesting speeds-ups. For example, with processors, we would get a speed-up of 3.79 (versus 1.4) in the nonadaptive version.

Conclusions and perspectives

The new adaptive formulation of the parareal algorithm opens the door to improve significantly the parallel efficiency of the method provided that the cost of the coarse solver is moderate. The increasing target tolerances which have to be met at each step allows to use online stopping criteria involving a posteriori estimators. The developed methodology remains theoretical since we have not quantified the impact of communication delays between processors nor potential memory issues (note however that the load balancing is devised with the purpose of equilibrating tasks and memory). In the framework of the ANR project “Ciné-Para (ANR-15-CE23-0019)”, we are working on these issues that are of a different level of theory and involve different collaborators. This will be the topic of another paper. In addition to this, several extensions based on the current findings are subject of ongoing works, in particular, the coupling of the adaptive parareal with adaptive space–time schemes, the coupling of parareal with internal iterative solvers like in the discussion of Appendix A.3 of the Appendix, and also the development of inexpensive coarse solvers.