Wave propagation characteristics of Parareal

Wave propagation characteristics of Parareal The paper derives and analyses the (semi-)discrete dispersion relation of the Parareal parallel-in-time integration method. It investigates Parareal’s wave propagation characteristics with the aim to better understand what causes the well documented stability problems for hyperbolic equations. The analysis shows that the instability is caused by convergence of the ampli- fication factor to the exact value from above for medium to high wave numbers. Phase errors in the coarse propagator are identified as the culprit, which suggests that specifically tailored coarse level methods could provide a remedy. Keywords Parareal · Convergence · Dispersion relation · Hyperbolic problem · Advection-dominated problem 1 Introduction even for the simple advection equation u + u = 0, Parareal t x is either unstable or inefficient. Numerical experiments reveal Parallel computing has become ubiquitous in science and that the instability emerges in the nonlinear case as a engineering, but requires suitable numerical algorithms to degradation of convergence with increasing Reynolds num- be efficient. Parallel-in-time integration methods have been ber [32]. Approaches exist to stabilise Parareal for hyperbolic identified as a promising direction to increase the level of con- equations [3,4,7,13,15,23,29], but typically with significant currency in computer simulations that involve the numerical overhead, leading to further degradation of efficiency, or lim- solution of time dependent partial differential equations [5]. ited applicability. A variety of methods has been proposed [8,10,12,22,25], Since a key characteristic of hyperbolic problems is the earliest going back to 1964 [27]. While even complex the existence of waves propagating with finite speeds, diffusive problems can be tackled successfully [11,14,24, understanding Parareal’s wave propagation characteristics is 33]—although parallel efficiencies remain low—hyperbolic important to understand and, hopefully, resolve these prob- or advection-dominated problems have proved to be much lems. However, no such analysis exists that gives insight into harder to parallelise in time. This currently prevents the use how the instability emerges. A better understanding of the of parallel-in-time integration for most problems in com- instability could show the way to novel methods that allow putational fluid dynamics, even though many applications the efficient and robust parallel-in-time solution of flows gov- struggle with excessive solution times and could benefit erned by advection. Additionally, just like for “classical” greatly from new parallelisation strategies. time stepping methods, detailed knowledge of Parareal’s the- For the Parareal parallel-in-time algorithm there is some oretical properties for test cases will help understanding its theory available illustrating its limitations in this respect. performance for complex test problems where mathematical Bal shows that Parareal with a sufficiently damping coarse theory is not available. method is unconditionally stable for parabolic problems but To this end, the paper derives a discrete dispersion rela- not for hyperbolic equations [2]. An early investigation of tion for Parareal to study how plane wave solutions u(x , t ) = Parareal’s stability properties showed instabilities for imag- exp(−i ωt ) exp(i κx ) are propagated in time. It studies the inary eigenvalues [31]. Gander and Vandewalle [17]givea discrete phase speed and amplification factor and how they detailed analysis of Parareal’s convergence and show that depend on e.g. the number of processors, choice of coarse propagator and other parameters. The analysis reveals that B Daniel Ruprecht the source of the instability is a convergence from above in d.ruprecht@leeds.ac.uk the amplification factor in higher wave number modes. In diffusive problems, where high wave numbers are naturally School of Mechanical Engineering, University of Leeds, Leeds LS2 9JT, UK 123 D. Ruprecht strongly damped, this does not cause any problems, but in cessors. If the number of iterations K is small enough and hyperbolic problems with little or no diffusion it causes the the coarse method much cheaper than the fine, iteration (4) amplification factors to exceed a value of one and thus trig- can run in less wall clock time than serially computing (3). gers instability. Furthermore, the paper identifies phase errors in the coarse propagator as the source of these issues. This 2.1 The Parareal iteration in matrix form suggests that controlling coarse level phase errors could be key to devising efficient parallel-in-time methods for hyper- As a first step toward deriving Parareal’s dispersion relation bolic equations. we will need to derive its stability function which will require All results presented in this paper have been produced writing it in matrix form. Consider now the case where both using pyParareal, a simple open-source Python code. It the coarse and the fine integrator are one-step methods with is freely available [28] to maximise the usefulness of the stability functions R and R . Then, G and F can be expressed f c here presented analysis, allowing other researchers to test as matrices different equations or to explore sets of parameters which are not analysed in this paper. u = F(u ) = Fu , u = G(u ) = Gu (5) p p−1 p−1 p p−1 p−1 f N 2 Parareal for linear problems with F := R (Aδt ) and G := (R (At )) . Denote f c k (P+1)n as u = (u ,..., u ) ∈ R a vector that contains the 0 P Parareal [25] is a parallel-in-time integration method for an approximate solutions at all time points T , j = 1,..., P initial value problem and the initial value u . Simple algebra shows that one step of iteration (4) is equivalent to the block matrix formulation u ˙(t ) = Au(t ), u(0) = u ∈ C , 0 ≤ t ≤ T . (1) k k−1 M u = M − M u + b (6) g g f For the sake of simplicity we consider only the linear case n,n where the right hand side is given by a matrix A ∈ C .To with matrices parallelise integration of (1) in time, Parareal decomposes the time interval [0, T ] into P time slices ⎡ ⎤ ⎢ ⎥ −FI [0, T]=[0, T ) ∪[T , T ) ∪ ··· ∪ [T , T ], (2) 1 1 2 P−1 ⎢ ⎥ (P+1)n,(P+1)n M := ∈ R (7) ⎢ ⎥ . . . . ⎣ ⎦ . . with P indicating the number of processors. Denote as F δt −FI and G two “classical” time integration methods with time steps of length δt and t (e.g. Runge–Kutta methods). For the and sake of simplicity, assume that all slice [T , T ) have the j −1 j same length T and that this length is an integer multiple of ⎡ ⎤ both δt and t so that T = N t and T = N δt.Below, c f ⎢ ⎥ −GI δt will always denote the time step size of the fine method and ⎢ ⎥ (P+1)n,(P+1)n M := ∈ R (8) ⎢ ⎥ . . . . ⎣ ⎦ t the time step size of the coarse method, so that we omit . . the indices and just write G and F to avoid clutter. Standard −GI serial time marching using the method denoted as F would correspond to evaluating (P+1)n and a vector b = (u , 0,..., 0) ∈ R . Formulation (6) interpretes Parareal as a preconditioned linear iteration [1]. u = F(u ), p = 1,..., P, (3) p p−1 where u ≈ u(T ). Instead, after an initialisation procedure 2.2 Stability function of Parareal p p to provide values u —typically running the coarse method once—Parareal computes the iteration From the matrix formulation of a single Parareal iteration (6), we can now derive its stability function, that is we can express k k k−1 k−1 the update from the initial value u to an approximation u at 0 P u = G u + F u − G u , p = 1,..., P p p−1 p−1 p−1 time T = T using Parareal with K iterations as multiplica- (4) tion by a single matrix. The fine propagator solution satisfies for k = 1,..., K where the computationally expensive evaluation of the fine method can be parallelised over P pro- M u = b (9) f f 123 Wave propagation characteristics of Parareal and is a fixed point of iteration (6). Therefore, propagation 2.3 Weak scaling versus longer simulation times of the error There are two different application scenarios for Parareal that k k e := u − u (10) f we can study when increasing the number of processors P. If we fix the final time T , increasing P will lead to better is governed by the matrix resolution since the coarse time step t cannot be larger than the length of a time slice T —the coarse method has to −1 E := I − M M (11) f perform at least one step per time slice. In this scenario, more processors are used to absorb the cost of higher temporal in the sense that resolution (“weak scaling”). Alternatively, we can use additional processors to compute k k−1 e = Ee . (12) until a later final time T and this is the scenario investigated here. Consequently, we study here the case where T and Using this notation and applying (6) recursively, it is now P increase together and always assume that T = P, that easy to show that is each time slice has length one and increasing P means parallelising over more time slices covering a longer time k−1 k 0 j −1 interval. Since dispersion properties of numerical methods u = Eu + E M b. (13) are typically analysed for a unit interval, this causes some j =0 issues that we resolve by “normalising” the Parareal stability 0 function, see Sect. 3.1. If, as is typically done, the start value u for the iteration is produced by a single run of the coarse method, that is if 2.4 Maximum singular value M u = b, (14) The matrix E defined in (11) determines how quickly Parareal Equation (13) further simplifies to converges. Note that E is nil-potent with E = 0, owing to the fact that after P iterations Parareal will always reproduce ⎛ ⎞ the fine solution exactly. Therefore, all eigenvalues of E are k j −1 ⎝ ⎠ u = E M b. (15) g zero and the spectral radius is not useful to analyse conver- j =0 gence. Below, to investigate convergence, we will therefore compute the maximum singular value σ of E instead. Since The right hand side vector can be generated from the initial value u via σ = E , b = C u (16) 1 0 it follows from (12) that by defining k k−1 k−1 k 0 e  ≤ E e  = σ e  ≤ σ e  (20) (P+1)n,n n,n Pn,n 2 2 2 2 C = [I; 0] ∈ R , I ∈ R , 0 ∈ R . (17) so that if σ< 1 Parareal will converge monotonically with- Finally, denote as out stalling. In particular, this rules out behaviour as found by Gander and Vandewalle for hyperbolic problems, where the n,Pn n,n C = [0, I] , 0 ∈ R , I ∈ R , (18) error would first increase substantially over the first P/2 iter- ations before beginning to decrease [16]. However, achieving the matrix that selects the last n entries out of u .Now,a fast convergence and good efficiency will typically require full Parareal update from some initial value u to an approx- σ  1. Note that if the coarse method is used to generate u , imation u using K iterations can be written compactly as it follows from (9) and (14) that the initial error is ⎛ ⎞ j −1 ⎝ ⎠ 0 0 −1 −1 u = C E M C u =: M u . (19) P 2 1 0 Parareal 0 e = u − u = M − M b. (21) g f g f j =0 n,n The stability matrix M ∈ R depends on K , T , The size of σ depends on the accuracy “gap” between the Parareal T , P, t, δt, F, G and A. Note that Staff and Rønquist coarse and fine integrator and the wave number. Figure 1 derived the stability function for the scalar case using Pas- shows σ for varying values of t when backward Euler is cal’s tree [31]. used for both coarse and fine method. Clearly, as the coarse 123 D. Ruprecht Fig. 1 Maximum singular value σ for decreasing coarse time step t Fig. 2 Projected speedup for pipelined Parareal [26]with P = 16 for ν = 0. Backward Euler is used for both F and G and for t = δt = processors for the same parameters as in Fig. 1 and the number of 0.05, both methods coincide so that σ = 0 iterations k fixed such that σ ≤ tol with tol = 0.01 time step approaches to fine time step of δt = 0.05, the max- 2.5 Convergence and (in)stability of Parareal imum singular value approaches zero. However, in this limit the coarse and fine propagator are identical and no speedup Two different but connected issues with Parareal are dis- is possible. The larger t compared to δt, the cheaper the cussed throughout the paper, convergence and (in)stability. coarse method becomes but since σ also grows, more iter- Here, convergence refers to how fast Parareal approaches the ations are likely required. Note that higher wave numbers fine solution within a single instance of Parareal, that is lead to higher values of σ while lower wave numbers tend to have values of σ  1 even for large coarse-to-fine time step M → F as k → P. (23) Parareal ratios. Looking at σ also provides a way to refine performance models for Parareal. Typically, in models projecting speedup, As discussed above, after k = P iterations we always have the number of iterations has to be fixed in addition to t, δt M = F, but particularly for hyperbolic problems Parareal and P. Instead, at least for linear problems, we can fix k such this can happen only at the final iteration k = P while at that k = P − 1 there is still a substantial difference [16]. Clearly, speedup is not obtainable in such a situation. The maximum singular value σ of E gives an upper bound or worst-case σ ≤ tol (22) scenario of how fast Parareal converges to the fine solution, cf. Equation (20). While σ< 1 does not necessarily guaran- for some fixed tolerance tol. The resulting projected speedup tee converge quick enough to generate meaningful speedup, for P = 16 processors and a tolerance of tol = 1e−2is it guarantees monotonous convergence and rules out an error shown in Fig. 2. First, as the coarse time step increases, the that increases first before decreasing only in later iterations. reduced cost of the coarse propagator improves achievable The other issue investigated in the paper is that of stabil- speedup. Simultaneously, the decreasing accuracy of G man- ity of repeated application of Parareal (“restarting”). Below, ifests itself in an increasing number of iterations required to stability is normally assessed for Parareal with a fixed num- match the selected tolerance. These two counteracting effects ber of iterations k. A configuration of Parareal is referred create a “sweet spot“ where G is accurate enough to still to as unstable if it leads to an amplification factor of more enable relatively fast convergence but cheap enough to allow than unity. This corresponds to an artificial increase in wave for speedup. It is noteworthy, however, that this sweet spot is amplitudes and, just as for classical time stepping methods, different for lower and higher wave numbers. Therefore, the would result in a diverging numerical approximation if the potential for speedup from Parareal does not solely depend method is used recursively on the solved equations and discretization parameters but also on the solution—the more prominent high wave num- ber modes are, the more restricted achievable speedup. M →∞ as n →∞. (24) ( ) Parareal 123 Wave propagation characteristics of Parareal While for classical methods this recursive application simply Using (30) to compute the resulting frequency yields ω = means stepping through time steps, for Parareal with P pro- i δ(k, U,ν) and retrieves the continuous plane wave solution cessors it would mean computing one window [0, T ], then i κx −νκ T i κ(x −UT ) restarting it with the final approximation as initial value for u(x , T ) =ˆ u e = e e (33) the next window [T , 2T ] and so on. P P of (25). It also reproduces the dispersion relation of the con- tinuous system 3 Discrete dispersion relation for the log(R) advection–diffusion equation 2 ω = i = i δ(κ, U,ν) = U κ − iνκ . (34) Starting from (19) we now derive the (semi-)discrete dis- However, if we use an approximate stability function R persion relation of Parareal for the one dimensional linear instead, we get some approximate ω = ω + i ω with r i advection diffusion problem ω ,ω ∈ R. The resulting semi-discrete solution becomes r i u + Uu = νu . (25) n −i ωt i κx ω t i κ x − t t x xx n i n ( n) u = e e = e e . (35) Therefore, ω /κ governs the propagation speed of the solu- First, assume a plane wave solution in space tion while ω governs the growth or decay in amplitude. Consequently, ω /κ is referred to as phase velocity while i κx u(x , t ) =ˆ u(t )e (26) exp(ω ) is called amplification factor. In the continuous case, the phase speed is equal to U and the amplification factor with wave number κ so that (25) reduces to the initial value −νκ equal to e . Note that for (25) the exact phase speed should problem be independent of the wave number κ. However, the discrete phase speed of a numerical method often will change with u ˆ (t ) =− Ui κ + νκ u ˆ(t ) =: δ(κ, U,ν)u ˆ(t ) (27) t κ, thus introducing numerical dispersion. Also note that for ν> 0 higher wave numbers decay faster because the ampli- fication factor decreases rapidly as κ increases. with initial value u ˆ(0) = 1. Integrating (27) from t = 0to t = T in one step gives 3.1 Normalisation u ˆ = R(δ, T )u ˆ (28) T 0 The update function R for Parareal in Eq. (19) denotes not ] where T = P is the an update over [0, 1] but over [0, T P P where R is the stability function of the employed method. number of processors. A phase speed of ω /κ = 1.0, for Now assume that the approximation of u ˆ is a discrete plane example, indicates a wave that travels across a whole inter- wave so that the solution at the end of the nth time slice is val [0, 1] during the step. If scaled up to an interval [0, P], given by the corresponding phase speed would become ω /κ = P instead. −i ωnT u ˆ = e . (29) n This enlarged range of values causes problems with the complex logarithm in (30). As an example, take the stabil- ity function of the exact propagator (32). Analytically, the Inserting this in (28)gives identity log(R) −i ωT e u ˆ = Ru ˆ ⇒ ω = i . (30) 0 0 R T ω = i = i δ(κ, U,ν) = i δ(κ, U,ν) (36) T T For R in polar coordinates, that is R = |R| exp(iθ) with holds, resulting in the correct dispersion relation (34)ofthe θ = angle(R), we get continuous system. However, depending on the values of κ, U, T and ν, this identity is not necessarily satisfied when −1 ω = i (log(|R|) + i θ ) T . (31) computing the complex logarithm using np.log. For exam- ple, for κ = 2, U = 1, ν = 0 and T > 0, the exact stability −2iT The exact integrator, for example, would read function is R = e . In Python, we obtain for T = 1 δ(κ,U ,ν)T −2i R = e . (32) 1j ∗ np.log(e )/1 = 2 = κ (37) exact 123 D. Ruprecht but for T = 2 we obtain continuous spatial operator. The two upper figures are for U = 1.0 and ν = 0.0 (no diffusion) while the two lower −4i 1j ∗ np.log(e )/2 ≈−1.1416 = κ (38) figures are for U = 1.0 and ν = 0.1 (diffusive). In both cases, the discrete phase speed of Parareal con- and so identity (36) is not fulfilled. The reason is that the verges almost monotonically toward the continuous phase logarithm of a complex number is not unique and so np.log speed. Even for ten iterations, Parareal still causes significant returns only a complex logarithm but not necessarily the right slowing of medium to large wave number modes. Parareal one in terms of the dispersion relation. requires almost the full number of iterations, k = 15, before To circumvent this problem, we “normalise” the update R it faithfully reproduces the correct phase speed across most for Parareal to [0, 1]. To this end, decompose of the spectrum. However, for any number of iterations where speedup might still be possible, Parareal will introduce sig- √ √ P P R = R · ... · R (39) nificant numerical dispersion. Slight artificial acceleration is also observed for high wave number modes for k = 15 where R corresponds to the propagation over [0, 1] instead in the non-diffusive and k = 10 in the diffusive case, but of [0, P]. Since there are P many roots R,wehavetoselect generally phase speeds are quite similar in the diffusive and the right one. First, we use the zeros function of numPy non-diffusive case. to find all P complex roots z of The amplification factor in the non-diffusive case (upper right figure) illustrates Parareal’s instability for hyperbolic z − R = 0. (40) equations: for k = 10 and k = 15 it is larger than one for a √ significant part of the spectrum, indicating that these modes Then, we select as root R the value z that satisfies are unstable and will be artificially amplified. For k = 5, the √ iteration has not yet corrected for the strong diffusivity of the θ( R) − θ = min θ(z ) − θ (41) targ p targ coarse propagator and remains stable for all modes but with p=1,...,P significant numerical damping of medium to high wave num- where θ is the angle function and θ some target angle, bers. The reason for the stability problems is discernible from targ which we still need to define. the amplification factor for the diffusive case (lower-right): from k = 0 (blue circles) to k = 5, Parareal reproduces the We compute ω and the resulting phase speed and ampli- fication factor for a range of wave numbers 0 ≤ κ ≤ κ ≤ correct amplification factor for small wave number modes 1 2 but significantly overestimates the amplitude of medium to ··· ≤ κ ≤ π.For κ , θ is set to the angle of the frequency N 1 targ ω computed from the analytic dispersion relation. After that, large wave numbers. It then continues to converge to the cor- θ is set to the angle of the root selected for the previous rect value from above. For the diffusive case where the exact targ value of κ. The rationale is that small changes in κ should values are smaller than one this does not cause instabilities. only result in small changes of frequency and phase so that In the non-diffusive case, however, any overestimation of the θ(ω ) ≈ θ(ω ) if the increment between wave numbers is analytical amplification factor immediately causes instabil- i −1 i small enough. From the selected root R we then compute ity. There is, in a sense, “no room” for the amplification factor ω using (30), the resulting discrete phase speed and ampli- to converge to the correct value from above. This also means fication factor and the target angle θ for the next wave that using a non-diffusive method as coarse propagator, for targ example trapezoidal rule, leads to disastrous consequences number. (not shown) where most parts of the spectrum are unstable for almost any value of k. 4 Analysis of the dispersion relation Figure 4 illustrate how the phase speed and amplitude errors discussed above manifest themselves. It shows a sin- After showing how to derive Parareal’s dispersion relation gle Gauss peak advected with a velocity of U = 1.0 with and normalising it to the unit interval, this section now pro- ν = 0.0 on a spatial domain [0, 4] over a time interval [0, 16] vides a detailed analysis of different factors influencing its distributed over P = 16 processors and N = 2 coarse time discrete phase speed and amplification factor. steps per slice. A spectral discretisation is used in space, allowing to represent the derivative exactly. For k = 5 iter- 4.1 Influence of diffusion ations, most higher wave numbers are damped out and the result looks essentially like a low wave number sine func- Figure 3 shows the discrete phase speed and amplification tion. The artificially amplified medium to high wave number modes create a “bulge” for k = 10 while dispersion leads factor of Parareal for P = 16, backward Euler with t = 1.0 as coarse and the exact integrator as fine propagator. Both to a significant trough at the sides of the domain. After fif- levels use δ =− Ui κ + νκ , that is the symbol of the teen iterations, the solution approximates the main part of 123 Wave propagation characteristics of Parareal Fig. 3 Discrete phase speed and amplification factor for Parareal with backward Euler as G and the exact integrator for F. The symbol for the spatial discretisation is δ =− i κ + νκ . The diffusion coefficient is ν = 0.0 (upper) and ν = 0.1(lower) the Gauss peak reasonably well, but dispersion still leads to ing room for Parareal to overestimate it without crossing the visible errors along the flanks. The right figure shows a part threshold to instability. For non-diffusive problems where the of the resulting spectrum. For k = 5, only the lowest wave exact amplification factor is one, every overestimation causes number modes are present, leading to the sine shaped solu- the mode to become unstable. tion. After k = 10 iterations, most of the spectrum is still being truncated but a small range of wave numbers around 4.2 Low order finite difference in coarse method κ = 0.05 is being artificially amplified which creates the “bulge” seen in the left figure. Finally, for k = 15 iterations, In a realistic scenario, some approximation of the spatial Parareal starts to correctly capture the spectrum but the still derivatives would have to be used instead of the exact symbol significant overestimation of low wave number amplitudes δ. For simple finite differences, we can study the effect this and underestimation of higher modes causes visible errors. has on the dispersion relation. Consider the first order upwind Observation 1 The amplification factor in Parareal for high- finite difference er wave numbers converges “from above”. In diffusive problems these wave numbers are damped, so the exact u − u j j −1 u (x ) ≈ amplification factor is significantly smaller than one, leav- x j 123 D. Ruprecht Fig. 4 Gauss peak in physical space (left) and corresponding spectrum (right) for U = 1.0and ν = 0.0 integrated using a pseudo-spectral method with 64 modes in space and Parareal with P = 16 processors in time Fig. 5 Phase speed (left) and amplification factor (right) of the implicit Euler method using the exact symbol δ (black circles) or the approximate symbol δ of the second order centred finite difference (blue squares) as approximation for u in (25). Assuming a discrete plane with initial value u ˆ(0) = 1 and a discrete symbol δ instead wave of δ as in (27). Note that i κ j x −iκx 1 − e u =ˆ u(t )e lim = i κ x →0 on a uniform spatial mesh x = j x instead of the continu- so that δ → δ as x → 0. Durran gives details for different ous plane wave (26), this leads to stencils [6]. i κx iκ(x −x ) −iκx j j e − e 1 − e The dispersion properties of the implicit Euler method i κx u (x ) ≈ = e . x j together with the first order upwind finite difference are x x qualitatively similar to the ones for implicit Euler with the For ν = 0 this results in the initial value problem −iκx 1 − e u ˆ (t ) =−U u ˆ(t ) =: δ(k, U,δx )u ˆ(t ) 123 Wave propagation characteristics of Parareal Fig. 6 Phase speed (left) and amplification factor (right) for same configuration as in Fig. 3 but with a second order centred finite difference in the coarse propagator instead of the exact symbol exact symbol (not shown). Using the upwind finite differ- strongly damps medium wave numbers while damping of ence instead of the exact symbol gives qualitatively similar high wave numbers is weak. wave propagation characteristics for the coarse propagator. In Parareal, this causes a situation similar to what hap- Numerical slowdown increases (up to the point where modes pens when using the trapezoidal rule as coarse propagator, at the higher end of the spectrum almost do not propagate at albeit less drastic. Figure 6 shows again the phase speed (left) all) and numerical diffusion becomes somewhat stronger. As and amplification factor (right) for the same configuration as a result, Parareal’s dispersion properties (also not shown) are before but implicit Euler with centred finite difference for G. also relatively similar, except for too small phase speeds even The failure of the coarse method to remove high wave number for k = 15. modes again leads to an earlier triggering of the instability. However, if we use the centred finite difference Whereas for Parareal using δ on the coarse level the iteration k = 5 was will stable (see Fig. 3), it is now unstable. For u − u j +1 j −1 iterations k = 10 and k = 15, large parts of the spectrum u (x ) ≈ x j 2x remain unstable. Also, the stronger numerical slow down of the coarse method makes it harder for Parareal to correct instead, this leads to an approximate symbol for phase speed errors. Where before Parareal with k = 15 iteration captured the exact phase speed reasonably well, in iκx −iκx Fig. 6 we still see significant numerical slow down of the e − e sin(κx ) δ =−U =−Ui . (42) higher wave number modes. 2x x Observation 2 The choice of finite difference stencil used In this case, it turns out that the dispersion properties of ˜ in the coarse propagator can have a significant effect on implicit Euler with δ and δ are quite different. Figure 5 shows Parareal. It seems that centred stencils that fail to remove the discrete phase speed (left) and amplification factor (right) high wave number modes cause similar problems as non- for both configurations. For the phase speed, both version ˜ diffusive time stepping methods, suggesting that stencils with agree qualitatively, even though using δ leads to noticeable upwind-bias are a much better choice. stronger slow down, particularly of higher wave numbers. For the amplification factor, however, there is a significant dif- ference between the semi-discrete and fully discrete method. 4.3 Influence of phase and amplitude errors While the former damps high wave numbers strongly, the combination of implicit Euler and centred finite differences To investigate whether phase errors or amplitude errors in the coarse method trigger the instability, we construct coarse propagators where either the phase or the amplitude is exact. The script plot_ieuler_dispersion.py supplied together Denote as R the stability function of backward Euler and euler with the Python code can be used to visualize the dispersion proper- ties of the coarse propagator alone. as R the stability function of the exact integrator. Then, a exact 123 D. Ruprecht method with the amplification factor of backward Euler and no phase speed error can be constructed as iθ(R ) exact R := |R | e (43) 1 euler while a method with no amplification error and the phase speed of backward Euler can be constructed as iθ(R ) euler R := |R | e . (44) 2 exact These artificially constructed propagators are now used within Parareal. Figure 7 shows the resulting amplification factor when using R (upper) or R (lower) as coarse propagator. For R , 1 2 1 where there is no phase speed error in the coarse propagator, there is no instability. Already for k = 10 it produces a good approximation of the exact amplification factor across the whole spectrum. In contrast, for R where there is no amplification error produced by G, the instability is clearly present for k = 5, k = 10 and k = 15. Figure 8 shows the solution for the same setup that was used for Fig. 4, except using the R artificial coarse propaga- tor without phase errors instead of the backward Euler. For k = 5 iterations, the peak is strongly damped but, because G has no phase errors, in the correct place. After k = 10 iterations, Parareal has corrected for most the numerical damping and already provides a reasonable approximation, even though the amplitude of most wave numbers in the spec- trum is still severely underestimated. However, the lack of phase errors and resulting numerical dispersion avoids the “bulge” and distortions that were present in Fig. 4. Finally, for k = 15 iterations, the solution provided by Parareal is indistinguishable from the exact solution. Small underesti- mation of the amplitudes of larger wave numbers can still Fig. 7 Amplification factor of Parareal for the advection equation for an artificially constructed coarse method with exact phase speed (upper) be seen in the spectrum, but the effect is minimal. Note that or exact amplification factor (lower) this does not mean that Parareal will provide speedup—in a realistic scenario, where F is not exact but a time stepping method, too, it would depend on how many iterations are method like the trapezoidal rule mentioned above. For R , required for Parareal to be as accurate as the fine method run however, σ remains below one across the whole spectrum, serially and the actual runtimes of both propagators. All that so that Parareal will converge monotonically for every mode. can be said so far is that avoiding coarse propagator phase Since σ approaches one for medium to high wave numbers, errors avoids the instability and leads to faster convergence. convergence there is potentially very slow, in line with the The effect of eliminating phase errors in the coarse method errors seen in the upper part of the spectrum of the Gauss can also be illustrated by analysing the maximum singular peak. However, in contrast to the other two cases, these wave value σ of the error propagation matrix. Figure 9 shows σ numbers will not trigger instabilities. depending on the wave number κ for three different coarse In summary, these results strongly suggest that phase propagators: the backward Euler, the artificially constructed errors in the coarse method are responsible for the insta- propagator R with no phase error and the artificially con- bility, which is in line with previous findings that Parareal structed propagator R with no amplitude error. For the can quickly correct even for very strong numerical diffusion backward Euler method, σ is larger than one for significant as long as a wave is placed at roughly the correct position by parts of the spectrum, indicating possible non-monotonous the coarse predictor [29]. convergence for these modes. The situation is even worse for R , mirroring the problems with a non-diffusive coarse 123 Wave propagation characteristics of Parareal Fig. 8 The same Gauss peak (left) and corresponding spectrum (right) as in Fig. 4 but now computed with the R coarse propagator with exact phase speed is given by δ(U ,κ,ν) −νκ t −Ui κ R = e = e e . (45) exact Therefore, we have −νκ |R | = e (46) exact and θ(R ) =−U κ. (47) exact Equivalent to the use of a coarse propagator R with exact phase propagation would be solving a transformed coarse level problem instead by setting Ui κt u ˜(t ) := e u ˆ(t ). (48) This leads to the purely diffusive coarse level problem Fig. 9 Maximum singular value σ of error propagation matrix E depending on the wave number of three choices of coarse propaga- u ˜ (t ) =−νκ u ˜(t ) (49) tor. R has exact phase speed while R has exact amplification factor 1 2 Ui κt −Ui κt with restriction operator e and interpolation e tak- ing care of the propagation part. This is precisely the strategy Observation 3 The instability in Parareal seems to be caused pursued in the nonlinear case by “asymptotic Parareal” where by phase errors in the coarse propagator while amplitude they factor out a fast term with purely imaginary eigenval- errors are quickly corrected by the iteration. ues, related to acoustic waves. In a sense, their approach can be understood as an attempt to construct a coarse method 4.3.1 Relation to asymptotic Parareal with minimal phase speed error. Of course, evaluation of the transformation is not trivial for more complex problems and It is interesting to point out how the R propagator with exact requires a sophisticated approach [30], in contrast to the here phase speed is related to the asymptotic Parareal method studied linear advection–diffusion equation where the trans- −Ui κt Ui κt developed by Haut et al. [18]. The exact propagator for (25) formation is simply multiplication by e and e . 123 D. Ruprecht 4.3.2 Phase error or mismatch So far, we have always assumed that the fine method is exact. This leaves the question whether the instability is triggered by phase errors in the coarse method or simply by a mismatch in phase speeds between fine and coarse level. In order to see if the instability arises if both fine and coarse level have the same large phase error, we replace the fine propagator stability function by iθ(R ) coarse R = |R | e . (50) 3 fine Now, the fine propagator is a method with exact amplification factor but a discrete phase speed that is as inaccurate as the coarse method. While such an integrator would not make for a very useful method in practice it is valuable for illustrative purposes. The coarse method is again the standard implicit Euler. Figure 10 shows the phase speed (upper) and amplifica- tion factor (lower) of Parareal for this configuration. Since the fine and coarse method have the same (highly inaccurate) phase speed, Parareal matches the fine method’s phase speed exactly from the first iteration and all lines coincide. The amplification factor converges quickly to the correct value and looks almost identical to the case where a coarse prop- agator with exact phase speed was used, compare for Fig. 7 (upper). No instability occurs and amplification factors are below one across the whole spectrum for all iterations. Figure 11 shows how Parareal converges for this configu- ration in physical and spectral space. Because both fine and coarse method now have substantial phase error, the Gauss peak is at a completely wrong position. However, for k = 10, Parareal already approximates it reasonably well and shows Fig. 10 Phase speed (upper) and amplification factor (lower) for an no sign of instability. Convergence looks again very similar artificially constructed fine propagator with the same phase error as the to the results shown in Fig. 8 except for the wrong position implicit Euler coarse propagator of the Gauss peak. While making the fine method as inaccu- rate as the coarse method is clearly not a useful strategy to multi-grid in time method [17]) would be a very interesting stabilise Parareal, this experiment nevertheless demonstrates direction for future research. Furthermore, it seems possi- that the instability is triggered by different discrete phase ble that parallel-in-time methods with more than two levels speeds in the coarse and fine method. like MGRIT [10] or PFASST [8] could yield some improve- Observation 4 Analysing further the issue of phase errors ment, because they would allow for less drastic changes in shows that the instability seems to arise from mismatches resolution compared to two-level Parareal. between the phase speed of coarse and fine propagator. It is interesting to note that a very similar observation was 4.4 Coarse time step made by Ernst and Gander for multi-grid methods for the Helmholtz equation. There, the “coarse grid correction fails Using a smaller time step for the coarse method will obvi- because of the incorrect dispersion relation (phase error) on ously reduce its phase error and can thus be expected to coarser and coarser grids […]” [9]. They find that adjust- benefit Parareal convergence. Figure 12 shows that this is ing the wave number of the coarse level problem in relation indeed true. It shows discrete phase speed (left) and ampli- to the mesh size leads to rapid convergence of their multi- fication factor (right) for the same configuration as used grid solver. Investigating if and how their approach might for Fig. 3, except now using two coarse step per time slice be applied to Parareal (which can also be considered as a instead of one. Since the coarse propagator alone is now 123 Wave propagation characteristics of Parareal Fig. 11 Gauss peak computed with the R fine propagator. The incorrect phase speed of the fine method puts the peak at a completely wrong position (left), but there is no instability and the spectrum (right) converges as quickly as for the exact phase speed coarse propagator in Fig. 8 Fig. 12 Phase speed (left) and amplification factor (right) for standard Parareal with the same configuration as for Fig. 3 except using two coarse time steps per time slice already significantly more accurate, Parareal with k = 5 and Observation 5 Since phase errors of the coarse method obvi- k = 10 iterations provides more accurate phase speeds and, ously depend on its time step size, reducing the coarse time for k = 15, reproduces the exact value exactly, The reduced step helps to reduce the range of unstable wave numbers and phase errors translate into a milder instability. For k = 10, the severity of the instability. some wave numbers have amplification factors above one, but both the range of unstable wave numbers and the sever- 4.5 Number of time slices ity of the instability are much smaller than if only a single coarse time step is used. This explains why configurations All examples so far only considered P = 16 time slices and can be quite successful where both F and G use nearly iden- processors. To illustrate the effect of increasing P,Fig. 13 tical time steps and the difference in runtime is achieved by shows the discrete dispersion relation for standard Parareal other means, e.g. an expensive high order spatial discretisa- for P = 64 time slices or processors (same configuration tion for the fine and a cheap low order discretisation on the as in Fig. 3 except for P). Even for k = 15 iterations, coarse level. Parareal reproduces the correct phase speed (left figure) very 123 D. Ruprecht Fig. 13 Phase error (left) and amplification factor for (right) P = 64 in contrast to P = 16 as in Fig. 3 poorly—waves across large parts of the spectrum suffer from substantial numerical slowdown. Also, converge is slow and there is only marginal improvement from k = 5to k = 15 iterations. Convergence is somewhat faster for the amplifica- tion factor (right figure) with more substantial improvement for k = 15 over the coarse method. However, there also remains significant numerical attenuation of the upper half of the spectrum. If integrating the Gauss peak with this con- figuration, the result at T = 64 after k = 5 iterations is essentially a straight line (not shown) as almost all modes beyond κ = 0 are strongly damped. A small overshoot at around κ = 1 is noticeable for k = 15 iterations and this will worsen as k increases. In general, as P increases, it takes more iterations to trigger the instability since the slow con- vergence requires longer to correct for the strong diffusivity of the coarse method. These results suggest that the high wave numbers are the Fig. 14 Maximum singular value of E depending on the number of processors P for three different wave numbers κ slowest to converge and that convergence deteriorates as P increases. This is confirmed by Fig. 14, showing the maxi- mum singular value for three wave numbers plotted against Parareal [29], as it removes exactly the high wave number P. Convergence generally gets worse as P increases, but modes that converge the slowest. note that even for P = 64 the low wave number mode (blue circles) still converges monotonically while the high wave number mode (green crosses) might already converge non- Observation 6 While convergence becomes slower as the monotonically for only P = 4 processors. There also seems number of processors P increases, low wave numbers con- to be a limit for σ as P increases, with higher wave numbers verge monotonically even for large numbers of P but high levelling off at higher values of σ . wave numbers might not do so already for P = 4. Therefore, Parareal could provide some speedup for linear hyperbolic problems if the solution consists mainly of very low wave number modes and/or numerical diffusion in the 4.6 Wave number fine propagator is sufficiently strong to remove higher wave number modes. This also explains why divergence damp- The analysis above showed that higher wave numbers con- verge slower and are more susceptible to instabilities. This ing in the fine propagator can accelerate convergence of is confirmed in Fig. 15 showing the difference between the 123 Wave propagation characteristics of Parareal Fig. 16 Maximum singular value of E depending on the wave number κ for three values of diffusivity ν the largest defect. As ν increases, however, the defects for κ = 2.69 will reduce further, the cross-over point will move to later iterations and eventually lower wave numbers will determine convergence for all iterations. In a sense, in line with the analysis above, Parareal propagates high wave num- ber modes very wrongly in both cases, but since high wave number modes are quickly attenuated if ν> 0, it does not matter very much in the diffusive case. Figure 16 confirms this for a wider range of wave numbers Fig. 15 Parareal defect versus number of iterations for ν = 0.0 (upper) κ. It shows the maximum singular value σ of the error prop- and ν = 0.1(lower) agation matrix E over the whole spectrum for three different values of ν.For ν = 0 (hyperbolic case), σ increases mono- tonically with κ and the highest wave number converges Parareal and fine integrator stability function the slowest. After around κ ≥ 0.8, the singular values are larger than one and convergence becomes potentially non- d(k) := |R (k) − R | . (51) Parareal fine monotone. For ν = 0.1, σ increases until around κ = 1.8 and then decreases again. Therefore, the slowest converging The smallest wave number, κ = 0.45, converges quickly mode is no longer at the end but in the middle of the spectrum. in the hyperbolic (upper) and diffusive case (lower). For Also, we now have σ< 1 for all κ so that all modes will con- κ = 0.9, both cases show some initial stalling before the verge, even though some potentially very slowly. Increasing mode converges. Finally, κ = 2.69 shows first a significant diffusion further to ν = 0.5 greatly improves convergence increase in the defect before convergence sets in as k comes for all modes, the largest σ across the whole spectrum is close to P = 16. Interestingly, this is the case for both ν = 0 now below 0.5. The worst converging mode has also moved and ν = 0.1. While the “bulge” is even more pronounced in “further down” the spectrum and is now at around κ = 1.0. −νκ t the diffusive case, since modes decay proportional to e This shows how the strong natural damping of high wave in amplitude, the absolute values of the defect are orders of numbers from diffusion counteracts Parareal’s tendency to magnitude smaller. Therefore, at least until k = 7 iterations, amplify them and thus stabilises it. it is no longer the high wave number mode κ = 2.69 that will restrict performance, but rather the lower wave number Observation 7 Since diffusion naturally damps higher wave κ = 0.9. Then, the instability for the high wave number kicks numbers, it removes the issue of slow or no convergence at the in and for k ≥ 8 wave number κ = 2.69 is again causing end of the spectrum. Therefore, as the diffusivity parameter 123 D. Ruprecht ν increases, the wave number that converges the slowest and analysis presented here to systems with multiple waves, e.g. determines performance of Parareal becomes smaller. the shallow water equations, or to nonlinear problem where wave numbers interact would be another interesting line of inquiry. Furthermore, the framework used here to analyse 5 Conclusions Parareal is straightforward to adopt for other parallel-in-time integration methods as long as a matrix representation for them is available. Efficient parallel-in-time integration of hyperbolic and ad- vection-dominated problems has been shown to be prob- Acknowledgements I would like to thank Martin Gander, Martin lematic. This prevents application of a promising new Schreiber and Beth Wingate for the very interesting discussions at the parallelisation strategy to many problems in computational 5th Workshop on Parallel-in-time integration at the Banff International fluid dynamics, despite the urgent need for better utilisation Research Station (BIRS) in November 2016, which led to the comments about asymptotic Parareal and multi-grid for 1D Helmholtz equation in of massively parallel computers. For the Parareal parallel- the paper. Parts of the pseudo-spectral code used to illustrate the effect in-time method, mathematical theory has shown that the of numerical dispersion come from David Ketcheson’s short course algorithm is either unstable or inefficient when applied to PseudoSpectralPython [21]. The pyParareal code written for this paper hyperbolic equations, but so far no detailed analysis exists of relies heavily on the open Python packages NumPy [34], SciPy [20] and Matplotlib [19]. how exactly the instability emerges. The paper presents the first detailed analysis of how Open Access This article is distributed under the terms of the Creative Parareal propagates waves and the ways in which the insta- Commons Attribution 4.0 International License (http://creativecomm bility is triggered. It uses a formulation of Parareal as a ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit preconditioned fixed point iteration for linear problems to to the original author(s) and the source, provide a link to the Creative derive its stability function. From there, a discrete disper- Commons license, and indicate if changes were made. sion relation is obtained that allows to study the phase speed and amplitude errors from Parareal when computing wave- like solutions. To deal with issues arising from increasing the References time interval together with the number of processors, a simple procedure is introduced to normalise the stability function to 1. Amodio, P., Brugnano, L.: Parallel solution in time of ODEs: some the unit interval. achievements and perspectives. Appl. Numer. Math. 59(3–4), 424– Analysis of the discrete dispersion relation and the max- 435 (2009). https://doi.org/10.1016/j.apnum.2008.03.024 imum singular value of the error propagation matrix then 2. Bal, G.: On the convergence and the stability of the Parareal algo- rithm to solve partial differential equations. In: Kornhuber, R., et al. allows to make a range of observations, illustrating where (eds.) Domain Decomposition Methods in Science and Engineer- the issues of Parareal for wave problems originate. A key ing, Lecture Notes in Computational Science and Engineering, vol. finding is that the source of the instability are different dis- 40, pp. 426–432. Springer, Berlin (2005). https://doi.org/10.1007/ crete phase speeds on the coarse and fine level, which cause 3-540-26825-1_43 3. Chen, F., Hesthaven, J.S., Zhu, X.: On the use of reduced basis instability of higher wave number modes. Interestingly, the methods to accelerate and stabilize the Parareal method. In: Quar- overestimation of high wave number amplitudes is present in teroni, A., Rozza, G. (eds.) Reduced Order Methods for Modeling diffusive problems, too, but since there these amplitudes are and Computational Reduction, MS&A—Modeling, Simulation naturally strongly damped, it does not trigger instabilities. and Applications, vol. 9, pp. 187–214. Springer, Berlin (2014). https://doi.org/10.1007/978-3-319-02090-7_7 Further analysis addresses the role of the number of proces- 4. Dai, X., Maday, Y.: Stable Parareal in time method for first- and sors, the coarse time step size and comments on possible second-order hyperbolic systems. SIAM J. Sci. Comput. 35(1), connections to asymptotic Parareal and multi-grid methods A52–A78 (2013). https://doi.org/10.1137/110861002 for the Helmholtz equation. 5. Dongarra, J., et al.: Applied mathematics research for exascale computing. Technical Report LLNL-TR-651000, Lawrence Liver- The analysis presented here will be useful to interpret and more National Laboratory (2014). http://science.energy.gov/%7E/ understand performance of Parareal for more complex prob- media/ascr/pdf/research/am/docs/EMWGreport.pdf lems in computational fluid dynamics. A natural line of future 6. Durran, D.R.: Numerical Methods for Fluid Dynamics, Texts in research would be to attempt to develop a new, more stable, Applied Mathematics, vol. 32. Springer, New York (2010). https:// doi.org/10.1007/978-1-4419-6412-0 parallel-in-time method for hyperbolic problems based on the 7. Eghbal, A., Gerber, A.G., Aubanel, E.: Acceleration of unsteady provided observations. For example, the update in Parareal hydrodynamic simulations using the Parareal algorithm. J. Comput. proceeds component wise. That means that if the coarse prop- Sci. (2016). https://doi.org/10.1016/j.jocs.2016.12.006 agator moves a wave at the wrong speed, the update will not 8. Emmett, M., Minion, M.L.: Toward an efficient parallel in time method for partial differential equations. Commun. Appl. Math. know that a simple shift of entries could provide a good cor- Comput. Sci. 7, 105–132 (2012). https://doi.org/10.2140/camcos. rection. Attempting to somehow modify the Parareal update 2012.7.105 to take into account this type of information seems promising, 9. Ernst, O.G., Gander, M.J.: Multigrid methods for Helmholtz prob- even though probably challenging to do in 3D. Extending the lems: a convergent scheme in 1D using standard components. In: 123 Wave propagation characteristics of Parareal Direct and Inverse Problems in Wave Propagation and Applica- 23. Kooij, G., Botchev, M., Geurts, B.: A block Krylov subspace imple- tions. De Gruyer (2013). https://doi.org/10.1515/9783110282283. mentation of the time-parallel Paraexp method and its extension 135 for nonlinear partial differential equations. J. Comput. Appl. Math. 10. Falgout, R.D., Friedhoff, S., Kolev, T.V., MacLachlan, S.P., 316, 229–246 (2017). https://doi.org/10.1016/j.cam.2016.09.036.. Schroder, J.B.: Parallel time integration with multigrid. SIAM (Selected papers from NUMDIFF-14) J. Sci. Comput. 36, C635–C661 (2014). https://doi.org/10.1137/ 24. Kreienbuehl, A., Naegel, A., Ruprecht, D., Speck, R., Wittum, 130944230 G., Krause, R.: Numerical simulation of skin transport using 11. Falgout, R.D., Manteuffel, T.A., O’Neill, B., Schroder, J.B.: Multi- Parareal. Comput. Vis. Sci. 17, 99–108(2015). https://doi.org/10. grid reduction in time for nonlinear parabolic problems: A case 1007/s00791-015-0246-y study. SIAM J. Sci. Comput. 39(5), S298–S322. (2016). https:// 25. Lions, J.L., Maday, Y., Turinici, G.: A “Parareal” in time discretiza- doi.org/10.1137/16M1082330 tion of PDE’s. Comptes Rendus l’Acad. Sci. Ser. I Math. 332, 12. Farhat, C., Chandesris, M.: Time-decomposed parallel time- 661–668 (2001). https://doi.org/10.1016/S0764-4442(00)01793-6 integrators: theory and feasibility studies for fluid, structure, and 26. Minion, M.L.: A hybrid Parareal spectral deferred corrections fluid–structure applications. Int. J. Numer. Methods Eng. 58(9), method. Commun. Appl. Math. Comput. Sci. 5(2), 265–301 (2010). 1397–1434 (2003). https://doi.org/10.1002/nme.860 https://doi.org/10.2140/camcos.2010.5.265 13. Farhat, C., Cortial, J., Dastillung, C., Bavestrello, H.: Time-parallel 27. Nievergelt, J.: Parallel methods for integrating ordinary differential implicit integrators for the near-real-time prediction of linear struc- equations. Commun. ACM 7(12), 731–733 (1964). https://doi.org/ tural dynamic responses. Int. J. Numer. Methods Eng. 67, 697–724 10.1145/355588.365137 (2006). https://doi.org/10.1002/nme.1653 28. Ruprecht, D.: v2.0 Parallel-in-time/pyparareal: wave propagation 14. Fischer, P.F., Hecht, F., Maday, Y.: A Parareal in time semi-implicit characteristics of parareal (2017). https://doi.org/10.5281/zenodo. approximation of the Navier–Stokes equations. In: Kornhuber, R., 1012274 et al. (eds.) Domain Decomposition Methods in Science and Engi- 29. Ruprecht, D., Krause, R.: Explicit parallel-in-time integration of neering, Lecture Notes in Computational Science and Engineering, a linear acoustic–advection system. Comput. Fluids 59, 72–83 vol. 40, pp. 433–440. Springer, Berlin (2005). https://doi.org/10. (2012). https://doi.org/10.1016/j.compfluid.2012.02.015 1007/3-540-26825-1_44 30. Schreiber, M., Peixoto, P.S., Haut, T.,Wingate, B.: Beyond spatial 15. Gander, M.J., Petcu, M.: Analysis of a Krylov subspace enhanced scalability limitations with a massively parallel method for linear Parareal algorithm for linear problem. ESAIM: Proc. 25, 114–129 oscillatory problems. Int. J. High Perform. Comput. Appl. (2017). (2008). https://doi.org/10.1051/proc:082508 https://doi.org/10.1177/1094342016687625 16. Gander, M.J., Vandewalle, S.: Analysis of the Parareal time-parallel 31. Staff, G.A., Rønquist, E.M.: Stability of the Parareal algorithm. time-integration method. SIAM J. Sci. Comput. 29(2), 556–578 In: Kornhuber, R., et al. (eds.) Domain Decomposition Methods in (2007). https://doi.org/10.1137/05064607X Science and Engineering, Lecture Notes in Computational Science 17. Gander, M.J., Vandewalle, S.: On the superlinear and linear con- and Engineering, vol. 40, pp. 449–456. Springer, Berlin (2005). vergence of the parareal algorithm. In: Widlund, O.B., Keyes, D.E. https://doi.org/10.1007/3-540-26825-1_46 (eds.) Domain Decomposition Methods in Science and Engineer- 32. Steiner, J., Ruprecht, D., Speck, R., Krause, R.: Convergence ing, Lecture Notes in Computational Science and Engineering, vol. of Parareal for the Navier–Stokes equations depending on the 55, pp. 291–298. Springer, Berlin (2007). https://doi.org/10.1007/ Reynolds number. In: Abdulle, A., Deparis, S., Kressner, D., 978-3-540-34469-8_34 Nobile, F., Picasso, M. (eds.) Numerical Mathematics and 18. Haut, T., Wingate, B.: An asymptotic parallel-in-time method for Advanced Applications—ENUMATH 2013, Lecture Notes in highly oscillatory PDEs. SIAM J. Sci. Comput. 36(2), A693–A713 Computational Science and Engineering, vol. 103, pp. 195– (2014). https://doi.org/10.1137/130914577 202. Springer, Berlin (2015). https://doi.org/10.1007/978-3-319- 19. Hunter, J.D.: Matplotlib: a 2D graphics environment. Comput. Sci. 10705-9_19 Eng. 9(3), 90–95 (2007). https://doi.org/10.1109/MCSE.2007.55 33. Trindade, J.M.F., Pereira, J.C.F.: Parallel-in-time simulation of the 20. Jones, E., Oliphant, T., Peterson, P., et al.: SciPy: open source sci- unsteady Navier–Stokes equations for incompressible flow. Int. J. entific tools for python (2001–). http://www.scipy.org/. Accessed Numer. Methods Fluids 45(10), 1123–1136 (2004). https://doi.org/ 28 May 2018 10.1002/fld.732 21. Ketcheson, D.I.: Pseudo spectral python (2015). https://github. 34. van der Walt, S., Colbert, S.C., Varoquaux, G.: The numpy array: com/ketch/PseudoSpectralPython. Accessed 28 May 2018 a structure for efficient numerical computation. Comput. Sci. Eng. 22. Kiehl, M.: Parallel multiple shooting for the solution of initial value 13(2), 22–30 (2011). https://doi.org/10.1109/MCSE.2011.37 problems. Parallel Comput. 20(3), 275–295 (1994). https://doi.org/ 10.1016/S0167-8191(06)80013-X http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Computing and Visualization in Science Springer Journals

Wave propagation characteristics of Parareal

Free
17 pages
Loading next page...
 
/lp/springer_journal/wave-propagation-characteristics-of-parareal-dRvdP5QWXs
Publisher
Springer Berlin Heidelberg
Copyright
Copyright © 2018 by The Author(s)
Subject
Mathematics; Computational Mathematics and Numerical Analysis; Computer Applications in Chemistry; Algorithms; Visualization; Numerical Analysis; Calculus of Variations and Optimal Control; Optimization
ISSN
1432-9360
eISSN
1433-0369
D.O.I.
10.1007/s00791-018-0296-z
Publisher site
See Article on Publisher Site

Abstract

The paper derives and analyses the (semi-)discrete dispersion relation of the Parareal parallel-in-time integration method. It investigates Parareal’s wave propagation characteristics with the aim to better understand what causes the well documented stability problems for hyperbolic equations. The analysis shows that the instability is caused by convergence of the ampli- fication factor to the exact value from above for medium to high wave numbers. Phase errors in the coarse propagator are identified as the culprit, which suggests that specifically tailored coarse level methods could provide a remedy. Keywords Parareal · Convergence · Dispersion relation · Hyperbolic problem · Advection-dominated problem 1 Introduction even for the simple advection equation u + u = 0, Parareal t x is either unstable or inefficient. Numerical experiments reveal Parallel computing has become ubiquitous in science and that the instability emerges in the nonlinear case as a engineering, but requires suitable numerical algorithms to degradation of convergence with increasing Reynolds num- be efficient. Parallel-in-time integration methods have been ber [32]. Approaches exist to stabilise Parareal for hyperbolic identified as a promising direction to increase the level of con- equations [3,4,7,13,15,23,29], but typically with significant currency in computer simulations that involve the numerical overhead, leading to further degradation of efficiency, or lim- solution of time dependent partial differential equations [5]. ited applicability. A variety of methods has been proposed [8,10,12,22,25], Since a key characteristic of hyperbolic problems is the earliest going back to 1964 [27]. While even complex the existence of waves propagating with finite speeds, diffusive problems can be tackled successfully [11,14,24, understanding Parareal’s wave propagation characteristics is 33]—although parallel efficiencies remain low—hyperbolic important to understand and, hopefully, resolve these prob- or advection-dominated problems have proved to be much lems. However, no such analysis exists that gives insight into harder to parallelise in time. This currently prevents the use how the instability emerges. A better understanding of the of parallel-in-time integration for most problems in com- instability could show the way to novel methods that allow putational fluid dynamics, even though many applications the efficient and robust parallel-in-time solution of flows gov- struggle with excessive solution times and could benefit erned by advection. Additionally, just like for “classical” greatly from new parallelisation strategies. time stepping methods, detailed knowledge of Parareal’s the- For the Parareal parallel-in-time algorithm there is some oretical properties for test cases will help understanding its theory available illustrating its limitations in this respect. performance for complex test problems where mathematical Bal shows that Parareal with a sufficiently damping coarse theory is not available. method is unconditionally stable for parabolic problems but To this end, the paper derives a discrete dispersion rela- not for hyperbolic equations [2]. An early investigation of tion for Parareal to study how plane wave solutions u(x , t ) = Parareal’s stability properties showed instabilities for imag- exp(−i ωt ) exp(i κx ) are propagated in time. It studies the inary eigenvalues [31]. Gander and Vandewalle [17]givea discrete phase speed and amplification factor and how they detailed analysis of Parareal’s convergence and show that depend on e.g. the number of processors, choice of coarse propagator and other parameters. The analysis reveals that B Daniel Ruprecht the source of the instability is a convergence from above in d.ruprecht@leeds.ac.uk the amplification factor in higher wave number modes. In diffusive problems, where high wave numbers are naturally School of Mechanical Engineering, University of Leeds, Leeds LS2 9JT, UK 123 D. Ruprecht strongly damped, this does not cause any problems, but in cessors. If the number of iterations K is small enough and hyperbolic problems with little or no diffusion it causes the the coarse method much cheaper than the fine, iteration (4) amplification factors to exceed a value of one and thus trig- can run in less wall clock time than serially computing (3). gers instability. Furthermore, the paper identifies phase errors in the coarse propagator as the source of these issues. This 2.1 The Parareal iteration in matrix form suggests that controlling coarse level phase errors could be key to devising efficient parallel-in-time methods for hyper- As a first step toward deriving Parareal’s dispersion relation bolic equations. we will need to derive its stability function which will require All results presented in this paper have been produced writing it in matrix form. Consider now the case where both using pyParareal, a simple open-source Python code. It the coarse and the fine integrator are one-step methods with is freely available [28] to maximise the usefulness of the stability functions R and R . Then, G and F can be expressed f c here presented analysis, allowing other researchers to test as matrices different equations or to explore sets of parameters which are not analysed in this paper. u = F(u ) = Fu , u = G(u ) = Gu (5) p p−1 p−1 p p−1 p−1 f N 2 Parareal for linear problems with F := R (Aδt ) and G := (R (At )) . Denote f c k (P+1)n as u = (u ,..., u ) ∈ R a vector that contains the 0 P Parareal [25] is a parallel-in-time integration method for an approximate solutions at all time points T , j = 1,..., P initial value problem and the initial value u . Simple algebra shows that one step of iteration (4) is equivalent to the block matrix formulation u ˙(t ) = Au(t ), u(0) = u ∈ C , 0 ≤ t ≤ T . (1) k k−1 M u = M − M u + b (6) g g f For the sake of simplicity we consider only the linear case n,n where the right hand side is given by a matrix A ∈ C .To with matrices parallelise integration of (1) in time, Parareal decomposes the time interval [0, T ] into P time slices ⎡ ⎤ ⎢ ⎥ −FI [0, T]=[0, T ) ∪[T , T ) ∪ ··· ∪ [T , T ], (2) 1 1 2 P−1 ⎢ ⎥ (P+1)n,(P+1)n M := ∈ R (7) ⎢ ⎥ . . . . ⎣ ⎦ . . with P indicating the number of processors. Denote as F δt −FI and G two “classical” time integration methods with time steps of length δt and t (e.g. Runge–Kutta methods). For the and sake of simplicity, assume that all slice [T , T ) have the j −1 j same length T and that this length is an integer multiple of ⎡ ⎤ both δt and t so that T = N t and T = N δt.Below, c f ⎢ ⎥ −GI δt will always denote the time step size of the fine method and ⎢ ⎥ (P+1)n,(P+1)n M := ∈ R (8) ⎢ ⎥ . . . . ⎣ ⎦ t the time step size of the coarse method, so that we omit . . the indices and just write G and F to avoid clutter. Standard −GI serial time marching using the method denoted as F would correspond to evaluating (P+1)n and a vector b = (u , 0,..., 0) ∈ R . Formulation (6) interpretes Parareal as a preconditioned linear iteration [1]. u = F(u ), p = 1,..., P, (3) p p−1 where u ≈ u(T ). Instead, after an initialisation procedure 2.2 Stability function of Parareal p p to provide values u —typically running the coarse method once—Parareal computes the iteration From the matrix formulation of a single Parareal iteration (6), we can now derive its stability function, that is we can express k k k−1 k−1 the update from the initial value u to an approximation u at 0 P u = G u + F u − G u , p = 1,..., P p p−1 p−1 p−1 time T = T using Parareal with K iterations as multiplica- (4) tion by a single matrix. The fine propagator solution satisfies for k = 1,..., K where the computationally expensive evaluation of the fine method can be parallelised over P pro- M u = b (9) f f 123 Wave propagation characteristics of Parareal and is a fixed point of iteration (6). Therefore, propagation 2.3 Weak scaling versus longer simulation times of the error There are two different application scenarios for Parareal that k k e := u − u (10) f we can study when increasing the number of processors P. If we fix the final time T , increasing P will lead to better is governed by the matrix resolution since the coarse time step t cannot be larger than the length of a time slice T —the coarse method has to −1 E := I − M M (11) f perform at least one step per time slice. In this scenario, more processors are used to absorb the cost of higher temporal in the sense that resolution (“weak scaling”). Alternatively, we can use additional processors to compute k k−1 e = Ee . (12) until a later final time T and this is the scenario investigated here. Consequently, we study here the case where T and Using this notation and applying (6) recursively, it is now P increase together and always assume that T = P, that easy to show that is each time slice has length one and increasing P means parallelising over more time slices covering a longer time k−1 k 0 j −1 interval. Since dispersion properties of numerical methods u = Eu + E M b. (13) are typically analysed for a unit interval, this causes some j =0 issues that we resolve by “normalising” the Parareal stability 0 function, see Sect. 3.1. If, as is typically done, the start value u for the iteration is produced by a single run of the coarse method, that is if 2.4 Maximum singular value M u = b, (14) The matrix E defined in (11) determines how quickly Parareal Equation (13) further simplifies to converges. Note that E is nil-potent with E = 0, owing to the fact that after P iterations Parareal will always reproduce ⎛ ⎞ the fine solution exactly. Therefore, all eigenvalues of E are k j −1 ⎝ ⎠ u = E M b. (15) g zero and the spectral radius is not useful to analyse conver- j =0 gence. Below, to investigate convergence, we will therefore compute the maximum singular value σ of E instead. Since The right hand side vector can be generated from the initial value u via σ = E , b = C u (16) 1 0 it follows from (12) that by defining k k−1 k−1 k 0 e  ≤ E e  = σ e  ≤ σ e  (20) (P+1)n,n n,n Pn,n 2 2 2 2 C = [I; 0] ∈ R , I ∈ R , 0 ∈ R . (17) so that if σ< 1 Parareal will converge monotonically with- Finally, denote as out stalling. In particular, this rules out behaviour as found by Gander and Vandewalle for hyperbolic problems, where the n,Pn n,n C = [0, I] , 0 ∈ R , I ∈ R , (18) error would first increase substantially over the first P/2 iter- ations before beginning to decrease [16]. However, achieving the matrix that selects the last n entries out of u .Now,a fast convergence and good efficiency will typically require full Parareal update from some initial value u to an approx- σ  1. Note that if the coarse method is used to generate u , imation u using K iterations can be written compactly as it follows from (9) and (14) that the initial error is ⎛ ⎞ j −1 ⎝ ⎠ 0 0 −1 −1 u = C E M C u =: M u . (19) P 2 1 0 Parareal 0 e = u − u = M − M b. (21) g f g f j =0 n,n The stability matrix M ∈ R depends on K , T , The size of σ depends on the accuracy “gap” between the Parareal T , P, t, δt, F, G and A. Note that Staff and Rønquist coarse and fine integrator and the wave number. Figure 1 derived the stability function for the scalar case using Pas- shows σ for varying values of t when backward Euler is cal’s tree [31]. used for both coarse and fine method. Clearly, as the coarse 123 D. Ruprecht Fig. 1 Maximum singular value σ for decreasing coarse time step t Fig. 2 Projected speedup for pipelined Parareal [26]with P = 16 for ν = 0. Backward Euler is used for both F and G and for t = δt = processors for the same parameters as in Fig. 1 and the number of 0.05, both methods coincide so that σ = 0 iterations k fixed such that σ ≤ tol with tol = 0.01 time step approaches to fine time step of δt = 0.05, the max- 2.5 Convergence and (in)stability of Parareal imum singular value approaches zero. However, in this limit the coarse and fine propagator are identical and no speedup Two different but connected issues with Parareal are dis- is possible. The larger t compared to δt, the cheaper the cussed throughout the paper, convergence and (in)stability. coarse method becomes but since σ also grows, more iter- Here, convergence refers to how fast Parareal approaches the ations are likely required. Note that higher wave numbers fine solution within a single instance of Parareal, that is lead to higher values of σ while lower wave numbers tend to have values of σ  1 even for large coarse-to-fine time step M → F as k → P. (23) Parareal ratios. Looking at σ also provides a way to refine performance models for Parareal. Typically, in models projecting speedup, As discussed above, after k = P iterations we always have the number of iterations has to be fixed in addition to t, δt M = F, but particularly for hyperbolic problems Parareal and P. Instead, at least for linear problems, we can fix k such this can happen only at the final iteration k = P while at that k = P − 1 there is still a substantial difference [16]. Clearly, speedup is not obtainable in such a situation. The maximum singular value σ of E gives an upper bound or worst-case σ ≤ tol (22) scenario of how fast Parareal converges to the fine solution, cf. Equation (20). While σ< 1 does not necessarily guaran- for some fixed tolerance tol. The resulting projected speedup tee converge quick enough to generate meaningful speedup, for P = 16 processors and a tolerance of tol = 1e−2is it guarantees monotonous convergence and rules out an error shown in Fig. 2. First, as the coarse time step increases, the that increases first before decreasing only in later iterations. reduced cost of the coarse propagator improves achievable The other issue investigated in the paper is that of stabil- speedup. Simultaneously, the decreasing accuracy of G man- ity of repeated application of Parareal (“restarting”). Below, ifests itself in an increasing number of iterations required to stability is normally assessed for Parareal with a fixed num- match the selected tolerance. These two counteracting effects ber of iterations k. A configuration of Parareal is referred create a “sweet spot“ where G is accurate enough to still to as unstable if it leads to an amplification factor of more enable relatively fast convergence but cheap enough to allow than unity. This corresponds to an artificial increase in wave for speedup. It is noteworthy, however, that this sweet spot is amplitudes and, just as for classical time stepping methods, different for lower and higher wave numbers. Therefore, the would result in a diverging numerical approximation if the potential for speedup from Parareal does not solely depend method is used recursively on the solved equations and discretization parameters but also on the solution—the more prominent high wave num- ber modes are, the more restricted achievable speedup. M →∞ as n →∞. (24) ( ) Parareal 123 Wave propagation characteristics of Parareal While for classical methods this recursive application simply Using (30) to compute the resulting frequency yields ω = means stepping through time steps, for Parareal with P pro- i δ(k, U,ν) and retrieves the continuous plane wave solution cessors it would mean computing one window [0, T ], then i κx −νκ T i κ(x −UT ) restarting it with the final approximation as initial value for u(x , T ) =ˆ u e = e e (33) the next window [T , 2T ] and so on. P P of (25). It also reproduces the dispersion relation of the con- tinuous system 3 Discrete dispersion relation for the log(R) advection–diffusion equation 2 ω = i = i δ(κ, U,ν) = U κ − iνκ . (34) Starting from (19) we now derive the (semi-)discrete dis- However, if we use an approximate stability function R persion relation of Parareal for the one dimensional linear instead, we get some approximate ω = ω + i ω with r i advection diffusion problem ω ,ω ∈ R. The resulting semi-discrete solution becomes r i u + Uu = νu . (25) n −i ωt i κx ω t i κ x − t t x xx n i n ( n) u = e e = e e . (35) Therefore, ω /κ governs the propagation speed of the solu- First, assume a plane wave solution in space tion while ω governs the growth or decay in amplitude. Consequently, ω /κ is referred to as phase velocity while i κx u(x , t ) =ˆ u(t )e (26) exp(ω ) is called amplification factor. In the continuous case, the phase speed is equal to U and the amplification factor with wave number κ so that (25) reduces to the initial value −νκ equal to e . Note that for (25) the exact phase speed should problem be independent of the wave number κ. However, the discrete phase speed of a numerical method often will change with u ˆ (t ) =− Ui κ + νκ u ˆ(t ) =: δ(κ, U,ν)u ˆ(t ) (27) t κ, thus introducing numerical dispersion. Also note that for ν> 0 higher wave numbers decay faster because the ampli- fication factor decreases rapidly as κ increases. with initial value u ˆ(0) = 1. Integrating (27) from t = 0to t = T in one step gives 3.1 Normalisation u ˆ = R(δ, T )u ˆ (28) T 0 The update function R for Parareal in Eq. (19) denotes not ] where T = P is the an update over [0, 1] but over [0, T P P where R is the stability function of the employed method. number of processors. A phase speed of ω /κ = 1.0, for Now assume that the approximation of u ˆ is a discrete plane example, indicates a wave that travels across a whole inter- wave so that the solution at the end of the nth time slice is val [0, 1] during the step. If scaled up to an interval [0, P], given by the corresponding phase speed would become ω /κ = P instead. −i ωnT u ˆ = e . (29) n This enlarged range of values causes problems with the complex logarithm in (30). As an example, take the stabil- ity function of the exact propagator (32). Analytically, the Inserting this in (28)gives identity log(R) −i ωT e u ˆ = Ru ˆ ⇒ ω = i . (30) 0 0 R T ω = i = i δ(κ, U,ν) = i δ(κ, U,ν) (36) T T For R in polar coordinates, that is R = |R| exp(iθ) with holds, resulting in the correct dispersion relation (34)ofthe θ = angle(R), we get continuous system. However, depending on the values of κ, U, T and ν, this identity is not necessarily satisfied when −1 ω = i (log(|R|) + i θ ) T . (31) computing the complex logarithm using np.log. For exam- ple, for κ = 2, U = 1, ν = 0 and T > 0, the exact stability −2iT The exact integrator, for example, would read function is R = e . In Python, we obtain for T = 1 δ(κ,U ,ν)T −2i R = e . (32) 1j ∗ np.log(e )/1 = 2 = κ (37) exact 123 D. Ruprecht but for T = 2 we obtain continuous spatial operator. The two upper figures are for U = 1.0 and ν = 0.0 (no diffusion) while the two lower −4i 1j ∗ np.log(e )/2 ≈−1.1416 = κ (38) figures are for U = 1.0 and ν = 0.1 (diffusive). In both cases, the discrete phase speed of Parareal con- and so identity (36) is not fulfilled. The reason is that the verges almost monotonically toward the continuous phase logarithm of a complex number is not unique and so np.log speed. Even for ten iterations, Parareal still causes significant returns only a complex logarithm but not necessarily the right slowing of medium to large wave number modes. Parareal one in terms of the dispersion relation. requires almost the full number of iterations, k = 15, before To circumvent this problem, we “normalise” the update R it faithfully reproduces the correct phase speed across most for Parareal to [0, 1]. To this end, decompose of the spectrum. However, for any number of iterations where speedup might still be possible, Parareal will introduce sig- √ √ P P R = R · ... · R (39) nificant numerical dispersion. Slight artificial acceleration is also observed for high wave number modes for k = 15 where R corresponds to the propagation over [0, 1] instead in the non-diffusive and k = 10 in the diffusive case, but of [0, P]. Since there are P many roots R,wehavetoselect generally phase speeds are quite similar in the diffusive and the right one. First, we use the zeros function of numPy non-diffusive case. to find all P complex roots z of The amplification factor in the non-diffusive case (upper right figure) illustrates Parareal’s instability for hyperbolic z − R = 0. (40) equations: for k = 10 and k = 15 it is larger than one for a √ significant part of the spectrum, indicating that these modes Then, we select as root R the value z that satisfies are unstable and will be artificially amplified. For k = 5, the √ iteration has not yet corrected for the strong diffusivity of the θ( R) − θ = min θ(z ) − θ (41) targ p targ coarse propagator and remains stable for all modes but with p=1,...,P significant numerical damping of medium to high wave num- where θ is the angle function and θ some target angle, bers. The reason for the stability problems is discernible from targ which we still need to define. the amplification factor for the diffusive case (lower-right): from k = 0 (blue circles) to k = 5, Parareal reproduces the We compute ω and the resulting phase speed and ampli- fication factor for a range of wave numbers 0 ≤ κ ≤ κ ≤ correct amplification factor for small wave number modes 1 2 but significantly overestimates the amplitude of medium to ··· ≤ κ ≤ π.For κ , θ is set to the angle of the frequency N 1 targ ω computed from the analytic dispersion relation. After that, large wave numbers. It then continues to converge to the cor- θ is set to the angle of the root selected for the previous rect value from above. For the diffusive case where the exact targ value of κ. The rationale is that small changes in κ should values are smaller than one this does not cause instabilities. only result in small changes of frequency and phase so that In the non-diffusive case, however, any overestimation of the θ(ω ) ≈ θ(ω ) if the increment between wave numbers is analytical amplification factor immediately causes instabil- i −1 i small enough. From the selected root R we then compute ity. There is, in a sense, “no room” for the amplification factor ω using (30), the resulting discrete phase speed and ampli- to converge to the correct value from above. This also means fication factor and the target angle θ for the next wave that using a non-diffusive method as coarse propagator, for targ example trapezoidal rule, leads to disastrous consequences number. (not shown) where most parts of the spectrum are unstable for almost any value of k. 4 Analysis of the dispersion relation Figure 4 illustrate how the phase speed and amplitude errors discussed above manifest themselves. It shows a sin- After showing how to derive Parareal’s dispersion relation gle Gauss peak advected with a velocity of U = 1.0 with and normalising it to the unit interval, this section now pro- ν = 0.0 on a spatial domain [0, 4] over a time interval [0, 16] vides a detailed analysis of different factors influencing its distributed over P = 16 processors and N = 2 coarse time discrete phase speed and amplification factor. steps per slice. A spectral discretisation is used in space, allowing to represent the derivative exactly. For k = 5 iter- 4.1 Influence of diffusion ations, most higher wave numbers are damped out and the result looks essentially like a low wave number sine func- Figure 3 shows the discrete phase speed and amplification tion. The artificially amplified medium to high wave number modes create a “bulge” for k = 10 while dispersion leads factor of Parareal for P = 16, backward Euler with t = 1.0 as coarse and the exact integrator as fine propagator. Both to a significant trough at the sides of the domain. After fif- levels use δ =− Ui κ + νκ , that is the symbol of the teen iterations, the solution approximates the main part of 123 Wave propagation characteristics of Parareal Fig. 3 Discrete phase speed and amplification factor for Parareal with backward Euler as G and the exact integrator for F. The symbol for the spatial discretisation is δ =− i κ + νκ . The diffusion coefficient is ν = 0.0 (upper) and ν = 0.1(lower) the Gauss peak reasonably well, but dispersion still leads to ing room for Parareal to overestimate it without crossing the visible errors along the flanks. The right figure shows a part threshold to instability. For non-diffusive problems where the of the resulting spectrum. For k = 5, only the lowest wave exact amplification factor is one, every overestimation causes number modes are present, leading to the sine shaped solu- the mode to become unstable. tion. After k = 10 iterations, most of the spectrum is still being truncated but a small range of wave numbers around 4.2 Low order finite difference in coarse method κ = 0.05 is being artificially amplified which creates the “bulge” seen in the left figure. Finally, for k = 15 iterations, In a realistic scenario, some approximation of the spatial Parareal starts to correctly capture the spectrum but the still derivatives would have to be used instead of the exact symbol significant overestimation of low wave number amplitudes δ. For simple finite differences, we can study the effect this and underestimation of higher modes causes visible errors. has on the dispersion relation. Consider the first order upwind Observation 1 The amplification factor in Parareal for high- finite difference er wave numbers converges “from above”. In diffusive problems these wave numbers are damped, so the exact u − u j j −1 u (x ) ≈ amplification factor is significantly smaller than one, leav- x j 123 D. Ruprecht Fig. 4 Gauss peak in physical space (left) and corresponding spectrum (right) for U = 1.0and ν = 0.0 integrated using a pseudo-spectral method with 64 modes in space and Parareal with P = 16 processors in time Fig. 5 Phase speed (left) and amplification factor (right) of the implicit Euler method using the exact symbol δ (black circles) or the approximate symbol δ of the second order centred finite difference (blue squares) as approximation for u in (25). Assuming a discrete plane with initial value u ˆ(0) = 1 and a discrete symbol δ instead wave of δ as in (27). Note that i κ j x −iκx 1 − e u =ˆ u(t )e lim = i κ x →0 on a uniform spatial mesh x = j x instead of the continu- so that δ → δ as x → 0. Durran gives details for different ous plane wave (26), this leads to stencils [6]. i κx iκ(x −x ) −iκx j j e − e 1 − e The dispersion properties of the implicit Euler method i κx u (x ) ≈ = e . x j together with the first order upwind finite difference are x x qualitatively similar to the ones for implicit Euler with the For ν = 0 this results in the initial value problem −iκx 1 − e u ˆ (t ) =−U u ˆ(t ) =: δ(k, U,δx )u ˆ(t ) 123 Wave propagation characteristics of Parareal Fig. 6 Phase speed (left) and amplification factor (right) for same configuration as in Fig. 3 but with a second order centred finite difference in the coarse propagator instead of the exact symbol exact symbol (not shown). Using the upwind finite differ- strongly damps medium wave numbers while damping of ence instead of the exact symbol gives qualitatively similar high wave numbers is weak. wave propagation characteristics for the coarse propagator. In Parareal, this causes a situation similar to what hap- Numerical slowdown increases (up to the point where modes pens when using the trapezoidal rule as coarse propagator, at the higher end of the spectrum almost do not propagate at albeit less drastic. Figure 6 shows again the phase speed (left) all) and numerical diffusion becomes somewhat stronger. As and amplification factor (right) for the same configuration as a result, Parareal’s dispersion properties (also not shown) are before but implicit Euler with centred finite difference for G. also relatively similar, except for too small phase speeds even The failure of the coarse method to remove high wave number for k = 15. modes again leads to an earlier triggering of the instability. However, if we use the centred finite difference Whereas for Parareal using δ on the coarse level the iteration k = 5 was will stable (see Fig. 3), it is now unstable. For u − u j +1 j −1 iterations k = 10 and k = 15, large parts of the spectrum u (x ) ≈ x j 2x remain unstable. Also, the stronger numerical slow down of the coarse method makes it harder for Parareal to correct instead, this leads to an approximate symbol for phase speed errors. Where before Parareal with k = 15 iteration captured the exact phase speed reasonably well, in iκx −iκx Fig. 6 we still see significant numerical slow down of the e − e sin(κx ) δ =−U =−Ui . (42) higher wave number modes. 2x x Observation 2 The choice of finite difference stencil used In this case, it turns out that the dispersion properties of ˜ in the coarse propagator can have a significant effect on implicit Euler with δ and δ are quite different. Figure 5 shows Parareal. It seems that centred stencils that fail to remove the discrete phase speed (left) and amplification factor (right) high wave number modes cause similar problems as non- for both configurations. For the phase speed, both version ˜ diffusive time stepping methods, suggesting that stencils with agree qualitatively, even though using δ leads to noticeable upwind-bias are a much better choice. stronger slow down, particularly of higher wave numbers. For the amplification factor, however, there is a significant dif- ference between the semi-discrete and fully discrete method. 4.3 Influence of phase and amplitude errors While the former damps high wave numbers strongly, the combination of implicit Euler and centred finite differences To investigate whether phase errors or amplitude errors in the coarse method trigger the instability, we construct coarse propagators where either the phase or the amplitude is exact. The script plot_ieuler_dispersion.py supplied together Denote as R the stability function of backward Euler and euler with the Python code can be used to visualize the dispersion proper- ties of the coarse propagator alone. as R the stability function of the exact integrator. Then, a exact 123 D. Ruprecht method with the amplification factor of backward Euler and no phase speed error can be constructed as iθ(R ) exact R := |R | e (43) 1 euler while a method with no amplification error and the phase speed of backward Euler can be constructed as iθ(R ) euler R := |R | e . (44) 2 exact These artificially constructed propagators are now used within Parareal. Figure 7 shows the resulting amplification factor when using R (upper) or R (lower) as coarse propagator. For R , 1 2 1 where there is no phase speed error in the coarse propagator, there is no instability. Already for k = 10 it produces a good approximation of the exact amplification factor across the whole spectrum. In contrast, for R where there is no amplification error produced by G, the instability is clearly present for k = 5, k = 10 and k = 15. Figure 8 shows the solution for the same setup that was used for Fig. 4, except using the R artificial coarse propaga- tor without phase errors instead of the backward Euler. For k = 5 iterations, the peak is strongly damped but, because G has no phase errors, in the correct place. After k = 10 iterations, Parareal has corrected for most the numerical damping and already provides a reasonable approximation, even though the amplitude of most wave numbers in the spec- trum is still severely underestimated. However, the lack of phase errors and resulting numerical dispersion avoids the “bulge” and distortions that were present in Fig. 4. Finally, for k = 15 iterations, the solution provided by Parareal is indistinguishable from the exact solution. Small underesti- mation of the amplitudes of larger wave numbers can still Fig. 7 Amplification factor of Parareal for the advection equation for an artificially constructed coarse method with exact phase speed (upper) be seen in the spectrum, but the effect is minimal. Note that or exact amplification factor (lower) this does not mean that Parareal will provide speedup—in a realistic scenario, where F is not exact but a time stepping method, too, it would depend on how many iterations are method like the trapezoidal rule mentioned above. For R , required for Parareal to be as accurate as the fine method run however, σ remains below one across the whole spectrum, serially and the actual runtimes of both propagators. All that so that Parareal will converge monotonically for every mode. can be said so far is that avoiding coarse propagator phase Since σ approaches one for medium to high wave numbers, errors avoids the instability and leads to faster convergence. convergence there is potentially very slow, in line with the The effect of eliminating phase errors in the coarse method errors seen in the upper part of the spectrum of the Gauss can also be illustrated by analysing the maximum singular peak. However, in contrast to the other two cases, these wave value σ of the error propagation matrix. Figure 9 shows σ numbers will not trigger instabilities. depending on the wave number κ for three different coarse In summary, these results strongly suggest that phase propagators: the backward Euler, the artificially constructed errors in the coarse method are responsible for the insta- propagator R with no phase error and the artificially con- bility, which is in line with previous findings that Parareal structed propagator R with no amplitude error. For the can quickly correct even for very strong numerical diffusion backward Euler method, σ is larger than one for significant as long as a wave is placed at roughly the correct position by parts of the spectrum, indicating possible non-monotonous the coarse predictor [29]. convergence for these modes. The situation is even worse for R , mirroring the problems with a non-diffusive coarse 123 Wave propagation characteristics of Parareal Fig. 8 The same Gauss peak (left) and corresponding spectrum (right) as in Fig. 4 but now computed with the R coarse propagator with exact phase speed is given by δ(U ,κ,ν) −νκ t −Ui κ R = e = e e . (45) exact Therefore, we have −νκ |R | = e (46) exact and θ(R ) =−U κ. (47) exact Equivalent to the use of a coarse propagator R with exact phase propagation would be solving a transformed coarse level problem instead by setting Ui κt u ˜(t ) := e u ˆ(t ). (48) This leads to the purely diffusive coarse level problem Fig. 9 Maximum singular value σ of error propagation matrix E depending on the wave number of three choices of coarse propaga- u ˜ (t ) =−νκ u ˜(t ) (49) tor. R has exact phase speed while R has exact amplification factor 1 2 Ui κt −Ui κt with restriction operator e and interpolation e tak- ing care of the propagation part. This is precisely the strategy Observation 3 The instability in Parareal seems to be caused pursued in the nonlinear case by “asymptotic Parareal” where by phase errors in the coarse propagator while amplitude they factor out a fast term with purely imaginary eigenval- errors are quickly corrected by the iteration. ues, related to acoustic waves. In a sense, their approach can be understood as an attempt to construct a coarse method 4.3.1 Relation to asymptotic Parareal with minimal phase speed error. Of course, evaluation of the transformation is not trivial for more complex problems and It is interesting to point out how the R propagator with exact requires a sophisticated approach [30], in contrast to the here phase speed is related to the asymptotic Parareal method studied linear advection–diffusion equation where the trans- −Ui κt Ui κt developed by Haut et al. [18]. The exact propagator for (25) formation is simply multiplication by e and e . 123 D. Ruprecht 4.3.2 Phase error or mismatch So far, we have always assumed that the fine method is exact. This leaves the question whether the instability is triggered by phase errors in the coarse method or simply by a mismatch in phase speeds between fine and coarse level. In order to see if the instability arises if both fine and coarse level have the same large phase error, we replace the fine propagator stability function by iθ(R ) coarse R = |R | e . (50) 3 fine Now, the fine propagator is a method with exact amplification factor but a discrete phase speed that is as inaccurate as the coarse method. While such an integrator would not make for a very useful method in practice it is valuable for illustrative purposes. The coarse method is again the standard implicit Euler. Figure 10 shows the phase speed (upper) and amplifica- tion factor (lower) of Parareal for this configuration. Since the fine and coarse method have the same (highly inaccurate) phase speed, Parareal matches the fine method’s phase speed exactly from the first iteration and all lines coincide. The amplification factor converges quickly to the correct value and looks almost identical to the case where a coarse prop- agator with exact phase speed was used, compare for Fig. 7 (upper). No instability occurs and amplification factors are below one across the whole spectrum for all iterations. Figure 11 shows how Parareal converges for this configu- ration in physical and spectral space. Because both fine and coarse method now have substantial phase error, the Gauss peak is at a completely wrong position. However, for k = 10, Parareal already approximates it reasonably well and shows Fig. 10 Phase speed (upper) and amplification factor (lower) for an no sign of instability. Convergence looks again very similar artificially constructed fine propagator with the same phase error as the to the results shown in Fig. 8 except for the wrong position implicit Euler coarse propagator of the Gauss peak. While making the fine method as inaccu- rate as the coarse method is clearly not a useful strategy to multi-grid in time method [17]) would be a very interesting stabilise Parareal, this experiment nevertheless demonstrates direction for future research. Furthermore, it seems possi- that the instability is triggered by different discrete phase ble that parallel-in-time methods with more than two levels speeds in the coarse and fine method. like MGRIT [10] or PFASST [8] could yield some improve- Observation 4 Analysing further the issue of phase errors ment, because they would allow for less drastic changes in shows that the instability seems to arise from mismatches resolution compared to two-level Parareal. between the phase speed of coarse and fine propagator. It is interesting to note that a very similar observation was 4.4 Coarse time step made by Ernst and Gander for multi-grid methods for the Helmholtz equation. There, the “coarse grid correction fails Using a smaller time step for the coarse method will obvi- because of the incorrect dispersion relation (phase error) on ously reduce its phase error and can thus be expected to coarser and coarser grids […]” [9]. They find that adjust- benefit Parareal convergence. Figure 12 shows that this is ing the wave number of the coarse level problem in relation indeed true. It shows discrete phase speed (left) and ampli- to the mesh size leads to rapid convergence of their multi- fication factor (right) for the same configuration as used grid solver. Investigating if and how their approach might for Fig. 3, except now using two coarse step per time slice be applied to Parareal (which can also be considered as a instead of one. Since the coarse propagator alone is now 123 Wave propagation characteristics of Parareal Fig. 11 Gauss peak computed with the R fine propagator. The incorrect phase speed of the fine method puts the peak at a completely wrong position (left), but there is no instability and the spectrum (right) converges as quickly as for the exact phase speed coarse propagator in Fig. 8 Fig. 12 Phase speed (left) and amplification factor (right) for standard Parareal with the same configuration as for Fig. 3 except using two coarse time steps per time slice already significantly more accurate, Parareal with k = 5 and Observation 5 Since phase errors of the coarse method obvi- k = 10 iterations provides more accurate phase speeds and, ously depend on its time step size, reducing the coarse time for k = 15, reproduces the exact value exactly, The reduced step helps to reduce the range of unstable wave numbers and phase errors translate into a milder instability. For k = 10, the severity of the instability. some wave numbers have amplification factors above one, but both the range of unstable wave numbers and the sever- 4.5 Number of time slices ity of the instability are much smaller than if only a single coarse time step is used. This explains why configurations All examples so far only considered P = 16 time slices and can be quite successful where both F and G use nearly iden- processors. To illustrate the effect of increasing P,Fig. 13 tical time steps and the difference in runtime is achieved by shows the discrete dispersion relation for standard Parareal other means, e.g. an expensive high order spatial discretisa- for P = 64 time slices or processors (same configuration tion for the fine and a cheap low order discretisation on the as in Fig. 3 except for P). Even for k = 15 iterations, coarse level. Parareal reproduces the correct phase speed (left figure) very 123 D. Ruprecht Fig. 13 Phase error (left) and amplification factor for (right) P = 64 in contrast to P = 16 as in Fig. 3 poorly—waves across large parts of the spectrum suffer from substantial numerical slowdown. Also, converge is slow and there is only marginal improvement from k = 5to k = 15 iterations. Convergence is somewhat faster for the amplifica- tion factor (right figure) with more substantial improvement for k = 15 over the coarse method. However, there also remains significant numerical attenuation of the upper half of the spectrum. If integrating the Gauss peak with this con- figuration, the result at T = 64 after k = 5 iterations is essentially a straight line (not shown) as almost all modes beyond κ = 0 are strongly damped. A small overshoot at around κ = 1 is noticeable for k = 15 iterations and this will worsen as k increases. In general, as P increases, it takes more iterations to trigger the instability since the slow con- vergence requires longer to correct for the strong diffusivity of the coarse method. These results suggest that the high wave numbers are the Fig. 14 Maximum singular value of E depending on the number of processors P for three different wave numbers κ slowest to converge and that convergence deteriorates as P increases. This is confirmed by Fig. 14, showing the maxi- mum singular value for three wave numbers plotted against Parareal [29], as it removes exactly the high wave number P. Convergence generally gets worse as P increases, but modes that converge the slowest. note that even for P = 64 the low wave number mode (blue circles) still converges monotonically while the high wave number mode (green crosses) might already converge non- Observation 6 While convergence becomes slower as the monotonically for only P = 4 processors. There also seems number of processors P increases, low wave numbers con- to be a limit for σ as P increases, with higher wave numbers verge monotonically even for large numbers of P but high levelling off at higher values of σ . wave numbers might not do so already for P = 4. Therefore, Parareal could provide some speedup for linear hyperbolic problems if the solution consists mainly of very low wave number modes and/or numerical diffusion in the 4.6 Wave number fine propagator is sufficiently strong to remove higher wave number modes. This also explains why divergence damp- The analysis above showed that higher wave numbers con- verge slower and are more susceptible to instabilities. This ing in the fine propagator can accelerate convergence of is confirmed in Fig. 15 showing the difference between the 123 Wave propagation characteristics of Parareal Fig. 16 Maximum singular value of E depending on the wave number κ for three values of diffusivity ν the largest defect. As ν increases, however, the defects for κ = 2.69 will reduce further, the cross-over point will move to later iterations and eventually lower wave numbers will determine convergence for all iterations. In a sense, in line with the analysis above, Parareal propagates high wave num- ber modes very wrongly in both cases, but since high wave number modes are quickly attenuated if ν> 0, it does not matter very much in the diffusive case. Figure 16 confirms this for a wider range of wave numbers Fig. 15 Parareal defect versus number of iterations for ν = 0.0 (upper) κ. It shows the maximum singular value σ of the error prop- and ν = 0.1(lower) agation matrix E over the whole spectrum for three different values of ν.For ν = 0 (hyperbolic case), σ increases mono- tonically with κ and the highest wave number converges Parareal and fine integrator stability function the slowest. After around κ ≥ 0.8, the singular values are larger than one and convergence becomes potentially non- d(k) := |R (k) − R | . (51) Parareal fine monotone. For ν = 0.1, σ increases until around κ = 1.8 and then decreases again. Therefore, the slowest converging The smallest wave number, κ = 0.45, converges quickly mode is no longer at the end but in the middle of the spectrum. in the hyperbolic (upper) and diffusive case (lower). For Also, we now have σ< 1 for all κ so that all modes will con- κ = 0.9, both cases show some initial stalling before the verge, even though some potentially very slowly. Increasing mode converges. Finally, κ = 2.69 shows first a significant diffusion further to ν = 0.5 greatly improves convergence increase in the defect before convergence sets in as k comes for all modes, the largest σ across the whole spectrum is close to P = 16. Interestingly, this is the case for both ν = 0 now below 0.5. The worst converging mode has also moved and ν = 0.1. While the “bulge” is even more pronounced in “further down” the spectrum and is now at around κ = 1.0. −νκ t the diffusive case, since modes decay proportional to e This shows how the strong natural damping of high wave in amplitude, the absolute values of the defect are orders of numbers from diffusion counteracts Parareal’s tendency to magnitude smaller. Therefore, at least until k = 7 iterations, amplify them and thus stabilises it. it is no longer the high wave number mode κ = 2.69 that will restrict performance, but rather the lower wave number Observation 7 Since diffusion naturally damps higher wave κ = 0.9. Then, the instability for the high wave number kicks numbers, it removes the issue of slow or no convergence at the in and for k ≥ 8 wave number κ = 2.69 is again causing end of the spectrum. Therefore, as the diffusivity parameter 123 D. Ruprecht ν increases, the wave number that converges the slowest and analysis presented here to systems with multiple waves, e.g. determines performance of Parareal becomes smaller. the shallow water equations, or to nonlinear problem where wave numbers interact would be another interesting line of inquiry. Furthermore, the framework used here to analyse 5 Conclusions Parareal is straightforward to adopt for other parallel-in-time integration methods as long as a matrix representation for them is available. Efficient parallel-in-time integration of hyperbolic and ad- vection-dominated problems has been shown to be prob- Acknowledgements I would like to thank Martin Gander, Martin lematic. This prevents application of a promising new Schreiber and Beth Wingate for the very interesting discussions at the parallelisation strategy to many problems in computational 5th Workshop on Parallel-in-time integration at the Banff International fluid dynamics, despite the urgent need for better utilisation Research Station (BIRS) in November 2016, which led to the comments about asymptotic Parareal and multi-grid for 1D Helmholtz equation in of massively parallel computers. For the Parareal parallel- the paper. Parts of the pseudo-spectral code used to illustrate the effect in-time method, mathematical theory has shown that the of numerical dispersion come from David Ketcheson’s short course algorithm is either unstable or inefficient when applied to PseudoSpectralPython [21]. The pyParareal code written for this paper hyperbolic equations, but so far no detailed analysis exists of relies heavily on the open Python packages NumPy [34], SciPy [20] and Matplotlib [19]. how exactly the instability emerges. The paper presents the first detailed analysis of how Open Access This article is distributed under the terms of the Creative Parareal propagates waves and the ways in which the insta- Commons Attribution 4.0 International License (http://creativecomm bility is triggered. It uses a formulation of Parareal as a ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit preconditioned fixed point iteration for linear problems to to the original author(s) and the source, provide a link to the Creative derive its stability function. From there, a discrete disper- Commons license, and indicate if changes were made. sion relation is obtained that allows to study the phase speed and amplitude errors from Parareal when computing wave- like solutions. To deal with issues arising from increasing the References time interval together with the number of processors, a simple procedure is introduced to normalise the stability function to 1. Amodio, P., Brugnano, L.: Parallel solution in time of ODEs: some the unit interval. achievements and perspectives. Appl. Numer. Math. 59(3–4), 424– Analysis of the discrete dispersion relation and the max- 435 (2009). https://doi.org/10.1016/j.apnum.2008.03.024 imum singular value of the error propagation matrix then 2. Bal, G.: On the convergence and the stability of the Parareal algo- rithm to solve partial differential equations. In: Kornhuber, R., et al. allows to make a range of observations, illustrating where (eds.) Domain Decomposition Methods in Science and Engineer- the issues of Parareal for wave problems originate. A key ing, Lecture Notes in Computational Science and Engineering, vol. finding is that the source of the instability are different dis- 40, pp. 426–432. Springer, Berlin (2005). https://doi.org/10.1007/ crete phase speeds on the coarse and fine level, which cause 3-540-26825-1_43 3. Chen, F., Hesthaven, J.S., Zhu, X.: On the use of reduced basis instability of higher wave number modes. Interestingly, the methods to accelerate and stabilize the Parareal method. In: Quar- overestimation of high wave number amplitudes is present in teroni, A., Rozza, G. (eds.) Reduced Order Methods for Modeling diffusive problems, too, but since there these amplitudes are and Computational Reduction, MS&A—Modeling, Simulation naturally strongly damped, it does not trigger instabilities. and Applications, vol. 9, pp. 187–214. Springer, Berlin (2014). https://doi.org/10.1007/978-3-319-02090-7_7 Further analysis addresses the role of the number of proces- 4. Dai, X., Maday, Y.: Stable Parareal in time method for first- and sors, the coarse time step size and comments on possible second-order hyperbolic systems. SIAM J. Sci. Comput. 35(1), connections to asymptotic Parareal and multi-grid methods A52–A78 (2013). https://doi.org/10.1137/110861002 for the Helmholtz equation. 5. Dongarra, J., et al.: Applied mathematics research for exascale computing. Technical Report LLNL-TR-651000, Lawrence Liver- The analysis presented here will be useful to interpret and more National Laboratory (2014). http://science.energy.gov/%7E/ understand performance of Parareal for more complex prob- media/ascr/pdf/research/am/docs/EMWGreport.pdf lems in computational fluid dynamics. A natural line of future 6. Durran, D.R.: Numerical Methods for Fluid Dynamics, Texts in research would be to attempt to develop a new, more stable, Applied Mathematics, vol. 32. Springer, New York (2010). https:// doi.org/10.1007/978-1-4419-6412-0 parallel-in-time method for hyperbolic problems based on the 7. Eghbal, A., Gerber, A.G., Aubanel, E.: Acceleration of unsteady provided observations. For example, the update in Parareal hydrodynamic simulations using the Parareal algorithm. J. Comput. proceeds component wise. That means that if the coarse prop- Sci. (2016). https://doi.org/10.1016/j.jocs.2016.12.006 agator moves a wave at the wrong speed, the update will not 8. Emmett, M., Minion, M.L.: Toward an efficient parallel in time method for partial differential equations. Commun. Appl. Math. know that a simple shift of entries could provide a good cor- Comput. Sci. 7, 105–132 (2012). https://doi.org/10.2140/camcos. rection. Attempting to somehow modify the Parareal update 2012.7.105 to take into account this type of information seems promising, 9. Ernst, O.G., Gander, M.J.: Multigrid methods for Helmholtz prob- even though probably challenging to do in 3D. Extending the lems: a convergent scheme in 1D using standard components. In: 123 Wave propagation characteristics of Parareal Direct and Inverse Problems in Wave Propagation and Applica- 23. Kooij, G., Botchev, M., Geurts, B.: A block Krylov subspace imple- tions. De Gruyer (2013). https://doi.org/10.1515/9783110282283. mentation of the time-parallel Paraexp method and its extension 135 for nonlinear partial differential equations. J. Comput. Appl. Math. 10. Falgout, R.D., Friedhoff, S., Kolev, T.V., MacLachlan, S.P., 316, 229–246 (2017). https://doi.org/10.1016/j.cam.2016.09.036.. Schroder, J.B.: Parallel time integration with multigrid. SIAM (Selected papers from NUMDIFF-14) J. Sci. Comput. 36, C635–C661 (2014). https://doi.org/10.1137/ 24. Kreienbuehl, A., Naegel, A., Ruprecht, D., Speck, R., Wittum, 130944230 G., Krause, R.: Numerical simulation of skin transport using 11. Falgout, R.D., Manteuffel, T.A., O’Neill, B., Schroder, J.B.: Multi- Parareal. Comput. Vis. Sci. 17, 99–108(2015). https://doi.org/10. grid reduction in time for nonlinear parabolic problems: A case 1007/s00791-015-0246-y study. SIAM J. Sci. Comput. 39(5), S298–S322. (2016). https:// 25. Lions, J.L., Maday, Y., Turinici, G.: A “Parareal” in time discretiza- doi.org/10.1137/16M1082330 tion of PDE’s. Comptes Rendus l’Acad. Sci. Ser. I Math. 332, 12. Farhat, C., Chandesris, M.: Time-decomposed parallel time- 661–668 (2001). https://doi.org/10.1016/S0764-4442(00)01793-6 integrators: theory and feasibility studies for fluid, structure, and 26. Minion, M.L.: A hybrid Parareal spectral deferred corrections fluid–structure applications. Int. J. Numer. Methods Eng. 58(9), method. Commun. Appl. Math. Comput. Sci. 5(2), 265–301 (2010). 1397–1434 (2003). https://doi.org/10.1002/nme.860 https://doi.org/10.2140/camcos.2010.5.265 13. Farhat, C., Cortial, J., Dastillung, C., Bavestrello, H.: Time-parallel 27. Nievergelt, J.: Parallel methods for integrating ordinary differential implicit integrators for the near-real-time prediction of linear struc- equations. Commun. ACM 7(12), 731–733 (1964). https://doi.org/ tural dynamic responses. Int. J. Numer. Methods Eng. 67, 697–724 10.1145/355588.365137 (2006). https://doi.org/10.1002/nme.1653 28. Ruprecht, D.: v2.0 Parallel-in-time/pyparareal: wave propagation 14. Fischer, P.F., Hecht, F., Maday, Y.: A Parareal in time semi-implicit characteristics of parareal (2017). https://doi.org/10.5281/zenodo. approximation of the Navier–Stokes equations. In: Kornhuber, R., 1012274 et al. (eds.) Domain Decomposition Methods in Science and Engi- 29. Ruprecht, D., Krause, R.: Explicit parallel-in-time integration of neering, Lecture Notes in Computational Science and Engineering, a linear acoustic–advection system. Comput. Fluids 59, 72–83 vol. 40, pp. 433–440. Springer, Berlin (2005). https://doi.org/10. (2012). https://doi.org/10.1016/j.compfluid.2012.02.015 1007/3-540-26825-1_44 30. Schreiber, M., Peixoto, P.S., Haut, T.,Wingate, B.: Beyond spatial 15. Gander, M.J., Petcu, M.: Analysis of a Krylov subspace enhanced scalability limitations with a massively parallel method for linear Parareal algorithm for linear problem. ESAIM: Proc. 25, 114–129 oscillatory problems. Int. J. High Perform. Comput. Appl. (2017). (2008). https://doi.org/10.1051/proc:082508 https://doi.org/10.1177/1094342016687625 16. Gander, M.J., Vandewalle, S.: Analysis of the Parareal time-parallel 31. Staff, G.A., Rønquist, E.M.: Stability of the Parareal algorithm. time-integration method. SIAM J. Sci. Comput. 29(2), 556–578 In: Kornhuber, R., et al. (eds.) Domain Decomposition Methods in (2007). https://doi.org/10.1137/05064607X Science and Engineering, Lecture Notes in Computational Science 17. Gander, M.J., Vandewalle, S.: On the superlinear and linear con- and Engineering, vol. 40, pp. 449–456. Springer, Berlin (2005). vergence of the parareal algorithm. In: Widlund, O.B., Keyes, D.E. https://doi.org/10.1007/3-540-26825-1_46 (eds.) Domain Decomposition Methods in Science and Engineer- 32. Steiner, J., Ruprecht, D., Speck, R., Krause, R.: Convergence ing, Lecture Notes in Computational Science and Engineering, vol. of Parareal for the Navier–Stokes equations depending on the 55, pp. 291–298. Springer, Berlin (2007). https://doi.org/10.1007/ Reynolds number. In: Abdulle, A., Deparis, S., Kressner, D., 978-3-540-34469-8_34 Nobile, F., Picasso, M. (eds.) Numerical Mathematics and 18. Haut, T., Wingate, B.: An asymptotic parallel-in-time method for Advanced Applications—ENUMATH 2013, Lecture Notes in highly oscillatory PDEs. SIAM J. Sci. Comput. 36(2), A693–A713 Computational Science and Engineering, vol. 103, pp. 195– (2014). https://doi.org/10.1137/130914577 202. Springer, Berlin (2015). https://doi.org/10.1007/978-3-319- 19. Hunter, J.D.: Matplotlib: a 2D graphics environment. Comput. Sci. 10705-9_19 Eng. 9(3), 90–95 (2007). https://doi.org/10.1109/MCSE.2007.55 33. Trindade, J.M.F., Pereira, J.C.F.: Parallel-in-time simulation of the 20. Jones, E., Oliphant, T., Peterson, P., et al.: SciPy: open source sci- unsteady Navier–Stokes equations for incompressible flow. Int. J. entific tools for python (2001–). http://www.scipy.org/. Accessed Numer. Methods Fluids 45(10), 1123–1136 (2004). https://doi.org/ 28 May 2018 10.1002/fld.732 21. Ketcheson, D.I.: Pseudo spectral python (2015). https://github. 34. van der Walt, S., Colbert, S.C., Varoquaux, G.: The numpy array: com/ketch/PseudoSpectralPython. Accessed 28 May 2018 a structure for efficient numerical computation. Comput. Sci. Eng. 22. Kiehl, M.: Parallel multiple shooting for the solution of initial value 13(2), 22–30 (2011). https://doi.org/10.1109/MCSE.2011.37 problems. Parallel Comput. 20(3), 275–295 (1994). https://doi.org/ 10.1016/S0167-8191(06)80013-X

Journal

Computing and Visualization in ScienceSpringer Journals

Published: May 29, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off