Table Of Content

NonamemanuscriptNo. (willbeinsertedbytheeditor) Wave propagation characteristics of Parareal DanielRuprecht 7 1 0 Received:date/Accepted:date 2 n Abstract Thepaperderivesandanalysesthe(semi-)discrete ForthePararealparallel-in-timealgorithmthereissome a J dispersion relation of the Parareal parallel-in-time integra- theory available illustrating its limitations in this respect. 5 tionmethod.ItinvestigatesParareal’swavepropagationchar- Bal shows that Parareal with a sufficiently damping coarse acteristics with the aim to better understand what causes methodisunconditionallystableforparabolicproblemsbut ] A thewelldocumentedstabilityproblemsforhyperbolicequa- not for hyperbolic equations [2]. An early investigation of tions. The analysis shows that the instability is caused by Parareal’sstabilitypropertiesshowedinstabilitiesforimag- N convergence of the amplification factor to the exact value inary eigenvalues [29]. Gander and Vandewalle [16] give a . h from above for medium to high wave numbers. Phase er- detailed analysis of Parareal’s convergence and show that t a rors in the coarse propagator are identified as the culprit, evenforthesimpleadvectionequationut+ux=0,Parareal m which suggests that specifically tailored coarse level meth- is either unstable or inefficient. Numerical experiments re- [ odscouldprovidearemedy. veal that the instability emerges in the nonlinear case as a degradationofconvergencewithincreasingReynoldsnum- 1 Keywords Parareal · convergence · dispersion relation · v ber [30]. Approaches exist to stabilise Parareal for hyper- hyperbolicproblem·advection-dominatedproblem 9 bolicequations[3,4,7,13,15,22,27],buttypicallywithsig- 5 nificantoverhead,leadingtofurtherdegradationofefficiency, 3 orlimitedapplicability. 1 1 Introduction 0 Sinceakeycharacteristicofhyperbolicproblemsisthe . existence of waves propagating with finite speeds, under- 1 Parallel computing has become ubiquitous in science and 0 standing Parareal’s wave propagation characteristics is im- engineering, but requires suitable numerical algorithms to 7 portanttounderstandand,hopefully,resolvetheseproblems. beefficient.Parallel-in-timeintegrationmethodshavebeen 1 However,nosuchanalysisexiststhatgivesinsightintohow : identified as a promising direction to increase the level of v theinstabilityemerges.Abetterunderstandingoftheinsta- concurrency in computer simulations that involve the nu- i bility could show the way to novel methods that allow the X mericalsolutionoftimedependentpartialdifferentialequa- efficient and robust parallel-in-time solution of flows gov- r tions [5]. A variety of methods has been proposed [8,10, a erned by advection. Additionally, just like for “classical” 12,21,24],theearliestgoingbackto1964[25].Whileeven timesteppingmethods,detailedknowledgeofParareal’sthe- complexdiffusiveproblemscanbetackledsuccessfully[11, oreticalpropertiesfortestcaseswillhelpunderstandingits 14,23,31]–althoughparallelefficienciesremainlow–hy- performanceforcomplextestcaseswheremathematicalthe- perbolic or advection-dominated problems have proved to oryisnotavailable. bemuchhardertoparalleliseintime.Thiscurrentlyprevents Tothisend,thepaperderivesadiscretedispersionrela- theuseofparallel-in-timeintegrationformostproblemsin tionforPararealtostudyhowplanewavesolutionsu(x,t)= computational fluid dynamics, even though many applica- exp(−iωt)exp(iκx) are propagated in time. It studies the tionsstrugglewithexcessivesolutiontimesandcouldbene- discretephasespeedandamplificationfactorandhowthey fitgreatlyfromnewparallelisationstrategies. depend on e.g. the number of processors, choice of coarse DanielRuprecht propagator and other parameters. The analysis reveals that SchoolofMechanicalEngineering,UniversityofLeeds,LS29JT,UK thesourceoftheinstabilityisaconvergencefromabovein 2 DanielRuprecht the amplification factor in higher wave number modes. In 2.1 ThePararealiterationinmatrixform diffusiveproblems,wherehighwavenumbersarenaturally strongly damped, this does not cause any problems, but in AsafirststeptowardderivingParareal’sdispersionrelation hyperbolicproblemswithlittleornodiffusionitcausesthe we will need to derive its stability function which will re- amplificationfactorstoexceedavalueofoneandthustrig- quirewritingitinmatrixform.Considernowthecasewhere gers instability. Furthermore, the paper identifies phase er- boththecoarseandthefineintegratorareone-stepmethods rors in the coarse propagator as the source of these issues. with stability functions Rf and Rc. Then, G and F can be Thissuggeststhatcontrollingcoarselevelphaseerrorscould expressedasmatrices bekeytodevisingefficientparallel-in-timemethodsforhy- u =F(u )=Fu , u =G(u )=Gu (5) perbolicequations. p p−1 p−1 p p−1 p−1 All results presented in this paper have been produced withF :=Rf(Aδt)Nf andG:=Rc(A∆t)Nc.Denoteasuk= using pyParareal, a simple open-source Python code. It (u ,...,u )∈R(P+1)navectorthatcontainstheapproximate 0 P is freely available [26] to maximise the usefulness of the solutions at all time points T , j =1,...,P and the initial j here presented analysis, allowing other researchers to test valueu .Simplealgebrashowsthatonestepofiteration(4) 0 different equations or to explore sets of parameters which isequivalenttotheblockmatrixformulation arenotanalysedinthispaper. M uk=(cid:0)M −M (cid:1)uk−1+b (6) g g f 2 Pararealforlinearproblems withmatrices   Parareal[24]isaparallel-in-timeintegrationmethodforan I initialvalueproblem −F I  u˙(t)=Au(t), u(t)=u ∈Cn, 0≤t≤T. (1) Mf := ... ... ∈R(P+1)n,(P+1)n (7) 0 −F I For the sake of simplicity we consider only the linear case wheretherighthandsideisgivenbyamatrixA∈Cn,n.To and parallelise integration of (1) in time, Parareal decomposes   I thetimeinterval[0,T]intoPtimeslices −G I  [0,T]=[0,T1)∪[T1,T2)∪...∪[TP−1,T], (2) Mg:= ... ... ∈R(P+1)n,(P+1)n (8) with P indicating the number of processors. Denote as F −GI δt andG two“classical”timeintegrationmethodswithtime ∆t and a vector b=(u ,0,...,0)∈R(P+1)n. Formulation (6) steps of length δt and ∆t (e.g. Runge-Kutta methods). For 0 interpretesPararealasapreconditionedlineariteration[1]. the sake of simplicity, assume that all slice [T ,T ) have j−1 j thesamelength∆T andthatthislengthisanintegermulti- ple of both δt and ∆t so that ∆T =N ∆t and ∆T =Nδt. c f 2.2 StabilityfunctionofParareal Below, δt will always denote the time step size of the fine method and ∆t the time step size of the coarse method, so FromthematrixformulationofasinglePararealiteration(6), thatweomittheindicesandjustwriteGandFtoavoidclut- we can now derive its stability function, that is we can ex- ter.Standardserialtimemarchingusingthemethoddenoted press the update from the initial value u to an approxima- asFwouldcorrespondtoevaluating 0 tion u at time T =P using Parareal with K iterations as P u =F(u ), p=1,...,P, (3) multiplicationbyasinglematrix.Thefinepropagatorsolu- p p−1 tionsatisfies whereu ≈u(T ).Instead,afteraninitialisationprocedure p p toprovidevaluesu0p –typicallyrunningthecoarsemethod Mfuf =b (9) once–Pararealcomputestheiteration and is a fixed point of iteration (6). Therefore, propagation uk =G(uk )+F(uk−1)−G(uk−1), p=1,...,P (4) p p−1 p−1 p−1 oftheerror fork=1,...,K wherethecomputationallyexpensiveeval- ek:=uk−u (10) f uationofthefinemethodcanbeparallelisedoverPproces- sors. If the number of iterations K is small enough and the isgovernedbythematrix coarsemethodmuchcheaperthanthefine,iteration(4)can runinlesswallclocktimethanseriallycomputing(3). E:=(cid:0)I−M−1M (cid:1). (11) g f WavepropagationcharacteristicsofParareal 3 Using this notation and applying (6) recursively, it is now 2.4 Maximumsingularvalue easytoshowthat ThematrixEdefinedin(11)determineshowquicklyPara- k−1 uk=Eu0+∑EjM−1b. (12) real converges. Note that E is nil-potent with EP =0, ow- g j=0 ing to the fact that after P iterations Parareal will always reproducethefinesolutionexactly.Therefore,alleigenval- If,asistypicallydone,thestartvalueu0 fortheiterationis ues of E are zero and the spectral radius is not useful to producedbyasinglerunofthecoarsemethod,thatisif analyseconvergence.Below,toinvestigateconvergence,we Mgu0=b, (13) willthereforecomputethemaximumsingularvalueσ ofE instead.Since Equation(12)furthersimplifiesto (cid:32) (cid:33) k σ =(cid:107)E(cid:107) , uk= ∑Ej M−1b. (14) 2 g j=0 Pararealwillconvergeifσ <1.However,incontrasttothe Therighthandsidevectorcanbegeneratedfromtheinitial spectral radius being smaller than one, this is only a suffi- valueu by cientcondition,butnotanecessaryone. 0 b=e ⊗u (15) 1 0 withe =(1,0,...,0)∈RP+1.Finally,denoteas 3 Discretedispersionrelationforthe 1 advection-diffusionequation C=[0,I], 0∈Rn,Pn, I∈Rn,n, (16) the matrix that selects the last n entries out of uk. Now, a Starting from (18) we now derive the (semi-)discrete dispersion relation of Parareal for the one dimensional linear fullPararealupdatefromsomeinitialvalueu toanapprox- 0 advectiondiffusionproblem imation u using K iterations can be written compactly as P u +Uu =νu . (19) t x xx (cid:32) (cid:33) K u =C ∑Ej M−1(e⊗u ). (17) First,assumeaplanewavesolutioninspace P g 0 j=0 u(x,t)=uˆ(t)eiκx (20) By inserting the unit basis vectors for u into (17), we can 0 nowconstructamatrixRsuchthat withwavenumberκ sothat(19)reducestotheinitialvalue problem u =Ru . (18) P 0 The stability matrix R depends on K, T, ∆T, P, ∆t, δt, F, uˆt(t)=−(cid:0)Uiκ+νκ2(cid:1)uˆ(t)=:δ(κ,U,ν)uˆ(t) (21) G and A. Note that Staff and Rønquist derived the stability with initial value uˆ(0)=1. Integrating (21) from t =0 to functionforthescalarcaseusingPascal’stree[29]. t=T inonestepgives uˆ =R(δ,T)uˆ (22) T 0 2.3 Strongversusweakscaling where R is the stability function of the employed method. There are two different scenarios that we can study when Nowassumethattheapproximationofuîsadiscreteplane increasing the number of processors P. If we fix the final wave so that the solution at the end of the nth time slice is timeT,increasingPwillleadtobetterresolutionsincethe givenby coarse time step ∆t cannot be larger than the length of a uˆ =e−iωn∆T. (23) time slice ∆T – the coarse method has to perform at least n onesteppertimeslice.Inthisscenario,moreprocessorsare Insertingthisin(22)gives usedtoabsorbthecostofhighertemporalresolution(“weak scaling”). e−iωTuˆ =Ruˆ ⇒ω =ilog(R). (24) 0 0 Here,however,weareinterestedinthecasewhereT and T Pincreasetogether.WethereforealwaysassumethatT =P, ForRinpolarcoordinates,thatisR=|R|exp(iθ)withθ = thatiseachtimeslicehaslengthoneandincreasingPmeans angle(R),weget parallelising over more time slices covering a longer time interval. Since dispersion properties of numerical methods ω =i(log(|R|)+iθ)T−1. (25) are typically analysed for a unit interval, this causes some Theexactintegrator,forexample,wouldread issuesthatweresolveby“normalising”thePararealstability function,see§3.1. R =eδ(κ,U,ν)T. (26) exact 4 DanielRuprecht Using (24) to compute the resulting frequency yields ω = where θ is the angle function and θ some target angle, targ iδ(k,U,ν)andretrievesthecontinuousplanewavesolution whichwestillneedtodefine. Wecomputeω andtheresultingphasespeedandampli- ficationfactorforarangeofwavenumbers0≤κ ≤κ ≤ u(x,T)=uˆ eiκx=e−νκ2Teiκ(x−UT) (27) 1 2 T ...≤κ ≤π.Forκ ,θ issettotheangleofthefrequency N 1 targ ωcomputedfromtheanalyticdispersionrelation.Afterthat, of(19).Italsoreproducesthedispersionrelationofthecon- θ issettotheangleoftherootselectedfortheprevious tinuoussystem targ value of κ. The rationale is that small changes in κ should log(R) onlyresultinsmallchangesoffrequencyandphasesothat ω =i =iδ(κ,U,ν)=Uκ−iνκ2. (28) T θ(ω )≈θ(ω)iftheincrementbetweenwavenumbersis i−1 i √ smallenough.Fromtheselectedroot PRwethencompute However, if we use an approximate stability function R in- ω using(24),theresultingdiscretephasespeedandampli- stead,wegetsomeapproximateω =ω +iω withω ,ω ∈ r i r i R.Theresultingsemi-discretesolutionbecomes fication factor and the target angle θtarg for the next wave number. unj =e−iωtneiκx=eωitneiκ(x−ωκrtn). (29) Therefore, ω /κ governs the propagation speed of the so- 4 Analysisofthedispersionrelation r lution while ω governs the growth or decay in amplitude. i Consequently, ω /κ is referred to as phase velocity while After showing how to derive Parareal’s dispersion relation r exp(ω)iscalledamplificationfactor.Inthecontinuouscase, andnormalisingittotheunitinterval,thissectionnowpro- i the phase speed is equal toU and the amplification factor vides a detailed analysis of different factors influencing its equal to e−νκ2. Note that for (19) the exact phase speed discretephasespeedandamplificationfactor. shouldbeindependentofthewavenumberκ.However,the discretephasespeedofanumericalmethodoftenwillchange 4.1 Influenceofdiffusion with κ, thus introducing numerical dispersion. Also note that for ν >0 higher wave numbers decay faster because Figure 1 shows the discrete phase speed and amplification theamplificationfactordecreasesrapidlyasκ increases. factorofPararealforP=16,backwardEulerwith∆t=1.0 as coarse and the exact integrator as fine propagator. Both 3.1 Normalisation levelsuseδ =−(cid:0)Uiκ+νκ2(cid:1),thatisthesymbolofthecon- tinuousspatialoperator.ThelefttwofiguresareforU=1.0 TheupdatefunctionRforPararealinEquation(18)denotes andν=0.0(nodiffusion)whilethetwoontherightarefor notanupdateover[0,1]butover[0,P]wherePisthenum- U =1.0andν =0.1(diffusive). ber of processors. A phase speed of ω /κ = 1.0, for ex- Inbothcases,thediscretephasespeedofPararealcon- r ample, indicates a wave that travels across a whole inter- verges almost monotonically toward the continuous phase val [0,1] during the step. If scaled up to an interval [0,P], speed. Even for ten iterations, Parareal still causes signifi- thecorrespondingphasespeedwouldbecomeω /κ=Pin- cant slowing of medium to large wave number modes. Pa- r stead. rarealrequiresalmostthefullnumberofiterations,k=15, Thisenlargedrangeofvaluescausesproblemswiththe beforeitfaithfullyreproducesthecorrectphasespeedacross complex logarithm in (24) because of the discontinuity in mostofthespectrum.However,foranynumberofiterations the atan2 function for R=−i. Therefore, we “normalise” where speed up might still be possible, Parareal will intro- theupdateRforPararealto[0,1].Tothisend,decompose ducesignificantnumericaldispersion.Slightartificialaccel- √ √ eration is also observed for high wave number modes for R= PR·...· PR (30) k=15inthenon-diffusiveandk=10inthediffusivecase, √ butgenerallyphasespeedsarequitesimilarinthediffusive where PRcorrespondstothepropagationover[0,1]instead √ andnon-diffusivecase. of[0,P].SincetherearePmanyroots PR,wehavetoselect Theamplificationfactorinthenon-diffusivecase(centre- therightone.First,weusethezerosfunctionofnumPyto leftfigure)illustratesParareal’sinstabilityforhyperbolicequa- findallPcomplexrootsz of i tions:fork=10andk=15,itislargerthanoneforasig- zn−R=0. (31) nificantpartofthespectrum,indicatingthatthesemodesare √ unstable and will be artificially amplified. For k = 5, the Then,weselectasroot PRthevaluezithatsatisfies iteration has not yet corrected for the strong diffusivity of (cid:12)(cid:12)(cid:12)θ(√PR)−θtarg(cid:12)(cid:12)(cid:12)= min (cid:12)(cid:12)θ(zp)−θtarg(cid:12)(cid:12) (32) twhiethcosiagrnseifipcraonptanguamtoerriacnadldraemmpaiinnsgsotfambleedfiourmatlol mhiogdhewsabvuet p=1,...,P WavepropagationcharacteristicsofParareal 5 1.0 1.0 1.0 1.0 0.8 or0.8 0.8 or0.8 Phase speed00..46 Exact plification fact00..46 Exact Phase speed00..46 Exact plification fact00..46 Exact Fine m Fine Fine m Fine Coarse A Coarse Coarse A Coarse 0.2 Parareal k=5 0.2 Parareal k=5 0.2 Parareal k=5 0.2 Parareal k=5 Parareal k=10 Parareal k=10 Parareal k=10 Parareal k=10 Parareal k=15 Parareal k=15 Parareal k=15 Parareal k=15 0.0 0.0 0.0 0.0 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 Wave number Wave number Wave number Wave number Fig.1 DiscretephasespeedandamplificationfactorforPararealwithbackwardEulerasGandtheexactintegratorforF.Thesymbolforthe spatialdiscretisationisδ=−(cid:0)iκ+νκ2(cid:1).Thediffusioncoefficientisν=0.0(left)andν=0.1(right). numbers.Thereasonforthestabilityproblemsisdiscernible dispersionleadstoasignificanttroughatthesidesofthedo- from the amplification factor for the diffusive case (right- main.Afterfifteeniterations,thesolutionapproximatesthe most figure): from k =0 (blue circles) to k =5, Parareal mainpartoftheGausspeakreasonablywell,butdispersion reproduces the correct amplification factor for small wave stillleadstovisibleerrorsalongtheflanks.Therightfigure numbermodesbutsignificantlyoverestimatestheamplitude shows a part of the resulting spectrum. For k=5, only the ofmediumtolargewavenumbers.Itthencontinuestocon- lowestwavenumbermodesarepresent,leadingtothesine vergetothecorrectvaluefromabove.Forthediffusivecase shaped solution. After k=10 iterations, most of the spec- where the exact values are smaller than one this does not trumisstillbeingtruncatedbutasmallrangeofwavenum- cause instabilities. In the non-diffusive case, however, any bers around κ =0.05 is being artificially amplified which overestimation of the analytical amplification factor imme- createsthe“bulge”seenintheleftfigure.Finally,fork=15 diatelycausesinstability.Thereis,inasense,“noroom”for iterations, Parareal start to correctly capture the spectrum theamplificationfactortoconvergetothecorrectvaluefrom butthestillsignificantoverestimationoflowwavenumber above.Thisalsomeansthatusinganon-diffusivemethodas amplitudesandunderestimationofhighermodesstillcauses coarsepropagator,forexampletrapezoidalrule,leadstodis- visibleerrors. astrous consequences (not shown) where most parts of the Observation1 TheamplificationfactorinPararealforhigh- spectrumareunstableforalmostanyvalueofk. erwavenumbersconverges“fromabove”.Indiffusiveprob- lemsthesewavenumbersaredamped,sotheexactamplifi- 1.4 100 cationfactorissignificantlysmallerthanone,leavingroom Parareal k=5 Exact 1.2 Parareal k=10 Parareal k=5 for Parareal to overestimate it without crossing the thresh- Parareal k=15 Parareal k=10 1.0 Exact 10-3 Parareal k=15 oldtoinstability.Fornon-diffusiveproblemswheretheexact 0.8 ) amplificationfactorisone,everyoverestimationcausesthe û u 0.6 bs( 10-6 modetobecomeunstable. 0.4 a 0.2 10-9 0.0 4.2 Loworderfinitedifferenceincoarsemethod 0.2 10-12 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 0.0 0.1 0.3 x Wave number In a realistic scenario, some approximation of the spatial Fig. 2 Gauss peak in physical space (left) and corresponding spec- derivativeswouldhavetobeusedinsteadoftheexactsym- trum(right)forU=1.0andν=0.0integratedusingapseudo-spectral bolδ.Forsimplefinitedifferences,wecanstudytheeffect methodwith64modesinspaceandPararealwithP=16processors intime. this has on the dispersion relation. Consider the first order upwindfinitedifference u −u Figure2illustratehowthephasespeedandamplitudeer- j j−1 u (x )≈ x j rorsdiscussedabovemanifestthemselves.Itshowsasingle ∆x GausspeakadvectedwithavelocityofU=1.0withν=0.0 as approximationfor ux in (19).Assuming adiscrete plane on a spatial domain [0,4] over a time interval [0,16] dis- wave tributedoverP=16processors.Apseudo-spectraldiscreti- u =uˆ(t)eiκj∆x j sationisusedinspace,allowingtorepresentthe derivative exactly.Fork=5iterations,mosthigherwavenumbersare onauniformspatialmeshxj = j∆xinsteadofthecontinu- dampedoutandtheresultlooksessentiallylikealowwave ousplanewave(20),thisleadsto number sine function. The artificially amplified medium to eiκxj−eiκ(xj−∆x) 1−e−iκ∆x highwavenumbermodescreatea“bulge”fork=10while ux(xj)≈ =eiκxj ∆x ∆x 6 DanielRuprecht Forν =0thisresultsintheinitialvalueproblem 1.0 1.0 uˆt(t)=−U1−∆e−xiκ∆xuˆ(t)=:δ˜(k,U,δx)uˆ(t) eed0.8 eed0.8 p0.6 p0.6 e s e s withinitialvalueuˆ(0)=1andadiscretesymbolδ˜ instead Phas0.4 Phas0.4 ofδ asin(21).Notethat 0.2 Continuous 0.2 Continuous Implicit Euler Implicit Euler 1−e−iκ∆x Implicit Euler + FD Implicit Euler + FD lim =iκ 0.00 1 2 3 0.00 1 2 3 ∆x→0 ∆x Wave number Wave number Fig. 3 Phase speed (left) and amplification factor (right) of the im- sothatδ˜ →δ as∆x→0.Durrangivesdetailsfordifferent plicitEulermethodusingtheexactsymbolδ (blackcircles)ortheap- stencils[6]. proximatesymbolδõfthesecondordercentredfinitedifference(blue The dispersion properties of the implicit Euler method squares). togetherwiththefirstorderupwindfinitedifferencearequal- itativelysimilartotheonesforimplicitEulerwiththeexact the instability. Whereas for Parareal using δ on the coarse symbol(notshown).1Usingtheupwindfinitedifferencein- leveltheiterationk=5waswillstable(seeFigure1),itis stead of the exact symbol gives qualitatively similar wave nowunstable.Foriterationsk=10andk=15,largepartsof propagation characteristics for the coarse propagator. Nu- thespectrumremainunstable.Also,thestrongernumerical merical slowdown increases (up to the point where modes slow down of the coarse method makes it harder for Para- at the higher end of the spectrum almost do not propagate realtocorrectforphasespeederrors.WherebeforeParareal atall)andnumericaldiffusionbecomessomewhatstronger. withk=15iterationcapturedtheexactphasespeedreason- Asaresult,Parareal’sdispersionproperties(alsonotshown) ablywell,inFigure4westillseesignificantnumericalslow are also relatively similar even, except for too small phase downofthehigherwavenumbermodes. speedsevenfork=15. However,ifweusethecentredfinitedifference 1.2 uj+1−uj−1 1.0 ux(xj)≈ 1.0 2∆x 0.8 or iδñs=te−adU,tehiiκs∆lxe−ades−tioκ∆anx =app−rUoxiismina(tκe∆syxm).bol (33) Phase speed00..46 EFixnaect mplification fact000...468 EFixnaect 2∆x ∆x Coarse A Coarse 0.2 Parareal k=5 0.2 Parareal k=5 Parareal k=10 Parareal k=10 In this case, it turns out that the dispersion properties of Parareal k=15 Parareal k=15 0.0 0.0 implicit Euler with δ and δ˜ are quite different. Figure 3 0 1 2 3 0 1 2 3 Wave number Wave number showsthediscretephasespeed(left)andamplificationfac- Fig.4 Phasespeed(left)andamplificationfactor(right)forsamecon- tor(right)forbothconfigurations.Forthephasespeed,both figurationasinFigure1butwithafirstorderupwindfinitedifference version agree qualitatively, even though using δ˜ leads to inthecoarsepropagatorinsteadoftheexactsymbol. noticeable stronger slow down, particularly of higher wave numbers. For the amplification factor, however, there is a significant difference between the semi-discrete and fully Observation2 The choice of finite difference stencil used discrete method. While the former damps high wave num- in the coarse propagator can have a significant effect on bers strongly, the combination of implicit Euler and cen- Parareal.Itseemsthatevenorderedstencilsthatfailtore- tred finite differences strongly damps medium wave num- move high wave number modes cause similar problems as berswhiledampingofhighwavenumbersisweak. non-diffusive time stepping methods, suggesting that sten- In Parareal, this causes a situation similar to what hap- cilswithupwind-biasareamuchbetterchoice. pens when using the trapezoidal rule as coarse propagator, albeit less drastic. Figure 4 shows again the phase speed (left) and amplification factor (right) for the same configu- 4.3 Influenceofphaseandamplitudeerrors rationasbeforebutimplicitEulerwithcentredfinitediffer- enceforG.Thefailureofthecoarsemethodtoremovehigh To investigate whether phase errors or amplitude errors in wave number modes again leads to an earlier triggering of thecoarsemethodtriggertheinstability,weconstructcoarse propagators where either phase or the amplitude is exact. 1 Thescriptplot ieuler dispersion.pysuppliedtogetherwith DenoteasR thestabilityfunctionofbackwardEulerand thePythoncodecanbeusedtovisualizethedispersionpropertiesof euler thecoarsepropagatoralone. asRexactthestabilityfunctionoftheexactintegrator.Then,a WavepropagationcharacteristicsofParareal 7 1.2 1.4 100 1.0 Parareal k=5 Exact 1.0 1.2 Parareal k=10 Parareal k=5 Parareal k=15 Parareal k=10 or0.8 or 1.0 Exact 10-3 Parareal k=15 plification fact00..46 Exact plification fact000...468 Exact u 000...468 âbs()u 10-6 m Fine m Fine A Coarse A Coarse 0.2 10-9 0.2 Parareal k=5 0.2 Parareal k=5 Parareal k=10 Parareal k=10 0.0 Parareal k=15 Parareal k=15 0.0 0.0 0.2 10-12 0 1 2 3 0 1 2 3 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 0.0 0.1 0.3 Wave number Wave number x Wave number Fig.5 AmplificationfactorofPararealfortheadvectionequationfor Fig.6 ThesameGausspeak(left)andcorrespondingspectrum(right) anartificiallyconstructedcoarsemethodwithexactphasespeed(left) asinFigure2butnowcomputedwiththeR coarsepropagatorwith 1 orexactamplificationfactor(right). exactphasespeed. methodwiththeamplificationfactorofbackwardEulerand Theeffectofeliminatingphaseerrorsinthecoarsemethod nophasespeederrorcanbeconstructedas can also be illustrated by analysing the maximum singular R :=|R |eiθ(Rexact) (34) value σ of the error propagation matrix. Figure 7 shows σ 1 euler dependingonthewavenumberκ forthreedifferentcoarse while a method with no amplification error and the phase propagators:thebackwardEuler,theartificiallyconstructed speedofbackwardEulercanbeconstructedas propagator R with no phase error and the artificially con- 1 structedpropagatorR withnoamplitudeerror.Fortheback- R2:=|Rexact|eiθ(Reuler). (35) 2 wardEulermethod,σ islargerthanoneforsignificantparts Theseartificiallyconstructedpropagatorsarenowusedwithin of the spectrum, indicating possible divergence of Parareal Parareal. for these modes. The situation is even worse for R2, mir- Figure 5 shows the resulting amplification factor when roringtheproblemswithanon-diffusivecoarsemethodlike using R1 (left) or R2 (right) as coarse propagator. For R1, the trapezoidal rule mentioned above. For R1, however, σ where there is no phase speed error in the coarse propaga- remainsbelowoneacrossthewholespectrum,sothatPara- tor,thereisnoinstability.Alreadyfork=10itproducesa realwillconvergeforeverymode.Sinceσ approachesone goodapproximationoftheexactamplificationfactoracross formediumtohighwavenumbers,convergencethereispo- the whole spectrum. In contrast, for R where there is no tentiallyveryslow,inlinewiththeerrorsseenintheupper 2 amplification error produced by G, the instability is clearly partofthespectrumoftheGausspeak.However,incontrast presentfork=5,k=10andk=15. totheothertwocases,thesewavenumberswillnottrigger Figure6showsthesolutionforthesamesetupthatwas instabilities. usedforFigure2,exceptusingtheR artificialcoarseprop- 1 agator without phase errors instead of the backward Euler. 1.4 For k=5 iterations, the peak is strongly damped but, be- causeGhasnophaseerrors,inthecorrectplace.Afterk= ue σ1.2 al1.0 10 iterations, Parareal has corrected for most the numeri- ular v0.8 caldampingandalreadyprovidesareasonableapproxima- g n tion,eventhoughtheamplitudeofmostwavenumbersinthe m si0.6 u m0.4 spectrumisstillseverelyunderestimated.However,thelack axi Backward Euler ofphaseerrorsandresultingnumericaldispersionavoidsthe M0.2 R1 R2 “bulge” and distortions that were present in Figure 2. Fi- 0.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Wave number nally,fork=15iterations,thesolutionprovidedbyParareal Fig.7 Maximumsingularvalueσ oferrorpropagationmatrixEde- isindistinguishablefromtheexactsolution.Smallunderes- pendingonthewavenumberofthreechoicesofcoarsepropagator.R 1 timationoftheamplitudesoflargerwavenumberscanstill hasexactphasespeedwhileR hasexactamplificationfactor. 2 beseeninthespectrum,buttheeffectisminimal.Notethat thisdoesnotmeanthatPararealwillprovidespeedup-ina realistic scenario, where F is not exact but a time stepping Insummary,theseresultsstronglysuggestthatphaseer- method,too,itwoulddependonhowmanyiterationsarere- rors in the coarse method are responsible for the instabil- quiredforPararealtobeasaccurateasthefinemethodrun ity,whichisinlinewithpreviousfindingsthatPararealcan seriallyandtheactualruntimesofbothpropagators.Allthat quickly correct even for very strong numerical diffusion as can be said so far is that avoiding coarse propagator phase long as a wave is placed at roughly the correct position by errorsavoidstheinstabilityandleadstofasterconvergence. thecoarsepredictor[27]. 8 DanielRuprecht Observation3 TheinstabilityinPararealseemstobecaused 1.2 1.0 by phase errors in the coarse propagator while amplitude 1.0 eRrerloartsioanretoqausicykmlyptcootircrePcaterdarbeyaltheiteration. Phase speed000...468 Exact plification factor000...468 Exact Fine m Fine Coarse A Coarse ItisinterestingtopointouthowtheR1 propagatorwithex- 0.2 PPaarraarreeaall kk==510 0.2 PPaarraarreeaall kk==510 actphasespeedisrelatedtotheasymptoticPararealmethod Parareal k=15 Parareal k=15 0.0 0.0 0 1 2 3 0 1 2 3 developedbyHautetal.[17].Theexactpropagatorfor(19) Wave number Wave number isgivenby Fig.8 Phasespeed(left)andamplificationfactor(right)foranarti- ficially constructed fine propagator with the same phase error as the R =eδ(U,κ,ν)=e−νκ2te−Uiκ. (36) implicitEulercoarsepropagator. exact Therefore,wehave 1.4 100 Parareal k=5 Exact 1.2 Parareal k=10 Parareal k=5 |Rexact|=e−νκ2 (37) 1.0 PEaxraacrteal k=15 10-3 PPaarraarreeaall kk==1105 0.8 ) and u 0.6 ˆs(u 10-6 b 0.4 a θ(R )=−Uκ. (38) exact 0.2 10-9 Equivalent to the use of a coarse propagator R with exact 0.0 1 0.2 10-12 phase propagation would be solving a transformed coarse 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 0.0 0.1 0.3 x Wave number levelprobleminsteadbysetting Fig.9 GausspeakcomputedwiththeR finepropagator.Theincor- 3 rectphasespeedofthefinemethodputsthepeakatacompletelywrong u˜(t):=eUiκtuˆ(t). (39) position(left),butthereisnoinstabilityandthespectrum(right)converges as quickly as for the exact phase speed coarse propagator in Thisleadstothepurelydiffusivecoarselevelproblem Figure6. u˜(t)=−νκ2u˜(t) (40) t thecoarsemethod.Whilesuchanintegratorwouldnotmake withrestrictionoperatoreUiκt andinterpolatione−Uiκt tak- for a very useful method in practice it is valuable for illus- ingcareofthepropagationpart.Thisispreciselythestrat- trative purposes. The coarse method is again the standard egypursuedinthenonlinearcaseby“asymptoticParareal” implicitEuler. wheretheyfactoroutafasttermwithpurelyimaginaryeigen- Figure8showsthephasespeed(left)andamplification values,relatedtoacousticwaves.Inasense,theirapproach factor(right)ofPararealforthisconfiguration.Sincethefine canbeunderstoodasanattempttoconstructacoarsemethod andcoarsemethodhavethesame(highlyinaccurate)phase withminimalphasespeederror.Ofcourse,evaluationofthe speed, Parareal matches the fine method’s phase speed ex- transformationisnottrivialformorecomplexproblemsand actlyfromthefirstiterationandalllinescoincide.Theam- requiresasophisticatedapproach[28],incontrasttothehere plificationfactorconvergesquicklytothecorrectvalueand studiedlinearadvection-diffusionequationwherethetrans- looks almost identical to the case where a coarse propaga- formationissimplymultiplicationbye−Uiκt andeUiκt. torwithexactphasespeedwasused,compareforFigure 5 (left).Noinstabilityoccursandamplificationfactorsarebe- Phaseerrorormismatch lowoneacrossthewholespectrumforalliterations. Figure9showshowPararealconvergesforthisconfigu- Sofar,wehavealwaysassumedthatthefinemethodisex- rationinphysicalandspectralspace.Becausebothfineand act. This leaves the question whether the instability is trig- coarse method now have substantial phase error, the Gauss gered by phase errors in the coarse method or simply by a peakisatacompletelywrongposition.However,fork=10, mismatchinphasespeedsbetweenfineandcoarselevel.In Pararealalreadyapproximatesitreasonablywellandshows order to see if the instability arises if both fine and coarse no sign of instability. Convergence looks again very simi- level have the same large phase error, we replace the fine lar to the results shown in Figure 6 except for the wrong propagatorstabilityfunctionby position of the Gauss peak. While making the fine method R =|R |eiθ(Rcoarse). (41) as inaccurate as the coarse method is clearly not a useful 3 fine strategy to stabilise Parareal, this experiment nevertheless Now,thefinepropagatorisamethodwithexactamplifica- demonstratesthattheinstabilityistriggeredbydifferentdis- tionfactorbutadiscretephasespeedthatisasinaccurateas cretephasespeedsinthecoarseandfinemethod. WavepropagationcharacteristicsofParareal 9 Observation4 Analysing further the issue of phase errors Observation5 Sincephaseerrorsofthecoarsemethodob- shows that the instability seems to arise from mismatches viouslydependonitstimestepsize,reducingthecoarsetime betweenthephasespeedofcoarseandfinepropagator. stephelpstoreducetherangeofunstablewavenumbersand theseverityoftheinstability. It is interesting to note that a very similar observation was made by Ernst and Gander for multi-grid methods for the Helmholtz equation. There, the “coarse grid correction 4.5 Numberoftimeslices fails because of the incorrect dispersion relation (phase error) on coarser and coarser grids [...]” [9]. They find that AllexamplessofaronlyconsideredP=16timeslicesand adjusting the wave number of the coarse level problem in processors.ToillustratetheeffectofincreasingP,Figure11 relationtothemeshsizeleadstorapidconvergenceoftheir showsthediscrete dispersionrelationforstandard Parareal multi-grid solver. Investigating if and how their approach forP=64timeslicesorprocessors(sameconfigurationas mightbeappliedtoParareal(whichcanalsobeconsidered in Figure 1 except for P). Even for k =15 iterations, Pa- asamulti-gridintimemethod[16])wouldbeaveryinter- rareal reproduces the correct phase speed (left figure) very estingdirectionforfutureresearch. poorly–wavesacrosslargepartsofthespectrumsufferfrom substantialnumericalslowdown.Also,convergeisslowand thereisonlymarginalimprovementfromk=5tok=15it- 4.4 Coarsetimestep erations.Convergenceissomewhatfasterfortheamplifica- tionfactor(rightfigure)withmoresubstantialimprovement Usingafinertimestepforthecoarsemethodwillobviously fork=15overthecoarsemethod.However,therealsore- reduce its phase error and can thus be expected to benefit mainssignificantnumericalattenuationoftheupperhalfof Parareal convergence. Figure 10 shows that this is indeed thespectrum.IfintegratingtheGausspeakwiththisconfig- true. It shows discrete phase speed (left) and amplification uration, the result at T =64 after k=5 iterations is essen- factor (right) for the same configuration as used for Fig- tiallyastraightline(notshown)asalmostallmodesbeyond ure 1, except now using two coarse step per time slice in- κ = 0 are strongly damped. A small overshoot at around stead of one. Since the coarse propagator alone is now al- κ=1isnoticeablefork=15iterationsandthiswillworsen ready significantly more accurate, Parareal with k=5 and askincreases.Ingeneral,asPincreases,ittakesmoreiter- k=10iterationsprovidesmoreaccuratephasespeedsand, ations to trigger the instability since the slow convergence fork=15,reproducestheexactvalueexactly,Thereduced requires longer to correct for the strong diffusivity of the phase errors translate into a milder instability. For k=10, coarsemethod. some wave numbers have amplification factors above one, butboththerangeofunstablewavenumbersandthesever- ity of the instability are much smaller than if only a single 1.0 1.0 coarse time step is used. This explains why configurations canbequitesuccessfulwherebothFandGusenearlyiden- 0.8 or0.8 totiitcohanelrftoimmreteahnsetsefi,pnes.egaa.nnaddnthaexecphdeeinfafspeivrleeonwhciegoihrndoerrurdndetirismscpreeatitisisaaaltcdihoiisnecvoreendtitsbhaye- Phase speed00..46 EFixnaect mplification fact00..46 EFixnaect Coarse A Coarse coarselevel. 0.2 Parareal k=5 0.2 Parareal k=5 Parareal k=10 Parareal k=10 Parareal k=15 Parareal k=15 0.0 0.0 0 1 2 3 0 1 2 3 1.2 Wave number Wave number 1.0 Fig.11 PhaseerrorandamplificationfactorforP=64incontrastto 1.0 P=16asinFigure1. 0.8 or Phase speed00..46 Exact plification fact000...468 Exact These results suggest that the high wave numbers are FCionaerse Am FCionaerse the slowest to converge and that convergence deteriorates 0.2 PPaarraarreeaall kk==510 0.2 PPaarraarreeaall kk==510 as P increases. This is confirmed by Figure 12, showing Parareal k=15 Parareal k=15 0.0 0.0 themaximumsingularvalueforthreewavenumbersplotted 0 1 2 3 0 1 2 3 Wave number Wave number againstP.ConvergencegenerallygetsworseasPincreases, Fig.10 Phasespeed(left)andamplificationfactor(right)forstandard but note that even for P=64 the low wave number mode PararealwiththesameconfigurationasforFigure1exceptusingtwo (blue circles) still converges while the high wave number coarsetimestepspertimeslice. mode (green crosses) might already be divergent for only P=4processors.Therealsoseemstobealimitforσ asP 10 DanielRuprecht increases,withhigherwavenumberslevellingoffathigher in and for k ≥8 wave number κ =2.69 is again causing valuesofσ. the largest defect. As ν increases, however, the defects for Therefore,Pararealcouldprovidesomespeedupforlin- κ=2.69willreducefurther,thecross-overpointwillmove ear hyperbolic problems if the initial data consists only of to later iterations and eventually lower wave numbers will verylowwavenumbermodes.Thisis,ofcourse,mainlyof determine convergence for all iterations. In a sense, in line academicinterest,sincealmosteveryrealisticproblemwill withtheanalysisabove,Pararealpropagateshighwavenum- havenonlinearitiesgeneratinghighwavenumbers.Thisex- bermodesverywronglyinbothcases,butsincehighwave plains, however, why divergence damping in the fine prop- number modes are quickly attenuated if ν >0, it does not agatorcanaccelerateconvergenceofParareal[27],asitre- matterverymuchinthediffusivecase. moves exactly the high wave number modes that converge theslowest. 102 102 100 100 1.6 10-2 10-2 1.4 ect 10-4 ect 10-4 ular value σ 11..02 Parareal def 111000-1--086 Parareal def 111000-1--086 m sing 00..68 1100--1142 ==00..4950 1100--1142 ==00..4950 Maximu 00..42 ==00..4950 1 3P=a2.r56a9rea7l Ite9rati1o1n k13 15 1 3P=a2.r56a9rea7l Ite9rati1o1n k13 15 =2.69 0.0 Fig.13 Pararealdefectversusnumberofiterationsforν=0.0(left) 10 20 30 40 50 60 andν=0.1(right). Number of processors Fig. 12 Maximum singular value of E depending on the number of processorsPforthreedifferentwavenumbersκ. Figure14confirmsthisforawiderrangeofwavenum- bersκ.Itshowsthemaximumsingularvalueσ oftheerror propagationmatrixEoverthewholespectrumforthreedif- Observation6 While convergence becomes slower as the ferentvaluesofν.Forν =0(hyperboliccase),σ increases number of processors P increases, low wave numbers re- monotonically with κ and the highest wave number con- mainconvergentevenforlargenumbersofPbuthighwave vergestheslowest.Afteraroundκ≥0.8,thesingularvalues numbersmightnotconvergealreadyforP=4. arelargerthanoneandmodesbecomepotentiallyunstable. Forν =0.1,σ increasesuntilaroundκ =1.8andthende- creasesagain.Therefore,theslowestconvergingmodeisno 4.6 Wavenumber longer at the end but in the middle of the spectrum. Also, we now have σ <1 for all κ so that all modes will con- Theanalysisaboveshowedthathigherwavenumberscon- verge,eventhoughsomepotentiallyveryslowly.Increasing verge slower and are more susceptible to instabilities. This diffusion further to ν =0.5 greatly improves convergence is confirmed in Figure 13 showing the difference between for all modes, the largest σ across the whole spectrum is thePararealandfineintegratorstabilityfunction nowbelow0.5.Theworstconvergingmodehasalsomoved (cid:12) (cid:12) “furtherdown”thespectrumandisnowataroundκ =1.0. d(k):=(cid:12)Rparareal(k)−Rfine(cid:12). (42) This shows how the strong natural damping of high wave numbers from diffusion counteracts Parareal’s tendency to Thesmallestwavenumber,κ =0.45,convergesquicklyin amplifythemandthusstabilisesit. the hyperbolic and diffusive case. For κ =0.9, both cases show some initial stalling before the mode converges. Fi- Observation7 Sincediffusionnaturallydampshigherwave nally, κ =2.69 shows first a significant increase in the de- numbers,itremovestheissueofslowornoconvergenceat fectbeforeconvergencesetsinaskcomesclosetoP=16. theendofthespectrum.Therefore,asthediffusivityparame- Interestingly, this is the case for both ν =0 and ν =0.1. terν increases,thewavenumberthatconvergestheslowest Whilethe“bulge”isevenmorepronouncedinthediffusive anddeterminesperformanceofPararealbecomessmaller. case,sincemodesdecayproportionaltoe−νκ2tinamplitude, the absolute values of the defect are orders of magnitude smaller. Therefore, at least until k = 7 iterations, it is no 5 Conclusions longer the high wave number mode κ =2.69 that will re- strict performance, but rather the lower wave number κ = Efficient parallel-in-time integration of hyperbolic and ad- 0.9. Then, the instability for the high wave number kicks vection-dominatedproblemshasbeenshowntobeproblem-