Semi-Lagrangian Exponential Integration with application to the rotating shallow water equations
SSEMI-LAGRANGIAN EXPONENTIAL INTEGRATION WITHAPPLICATION TO THE ROTATING SHALLOW WATEREQUATIONS ∗ PEDRO S. PEIXOTO † AND
MARTIN SCHREIBER ‡ Abstract.
In this paper we propose a novel way to integrate time-evolving partial differen-tial equations that contain nonlinear advection and stiff linear operators, combining exponentialintegration techniques and semi-Lagrangian methods.The general formulation is built from the solution of an integration factor problem with re-spect to the problem written with a material derivative, so that the exponential integration schemenaturally incorporates the nonlinear advection. Semi-Lagrangian techniques are used to treat thedependence of the exponential integrator on the flow trajectories. The formulation is general, asmany exponential integration techniques could be combined with different semi-Lagrangian meth-ods. This formulation allows an accurate solution of the linear stiff operator, a property inherited bythe exponential integration technique. It also provides a sufficiently accurate representation of thenonlinear advection, even with large time-step sizes, a property inherited by the semi-Lagrangianmethod.Aiming for application in weather and climate modeling, we discuss possible combinations ofwell established exponential integration techniques and state-of-the-art semi-Lagrangian methodsused operationally in the application. We show experiments for the planar rotating shallow waterequations. When compared to traditional exponential integration techniques, the experiments revealthat the coupling with semi-Lagrangian allows stabler integration with larger time-step sizes. Fromthe application perspective, which already uses semi-Lagrangian methods, the exponential treatmentcould improve the solution of wave-dispersion when compared to semi-implicit schemes.
Key words.
Exponential integrator, semi-Lagrangian, nonlinear advection, rotating shallowwater equations, weather and climate modeling.
AMS subject classifications.
1. Introduction.
Consider an autonomous initial value problem of the form(1) ∂u∂t = L ( u ) + N ( u ) , u (0) = u where L is a linear (possibly differential) operator and N is a function (usually non-linear). Exponential integrators are usually derived making use of exponentials of adiscrete form of the linear operator L . Many schemes of this form exists, as one maynotice from the review of [29].Several application models, such as those related to fluid dynamics [1], have animportant advection term in the equations, usually nonlinear ( (cid:126)v · ∇ u ). This can berepresented as(2) DuDt = ∂u∂t + (cid:126)v · ∇ u = L ( u ) + ˜ N ( u ) , u (0) = u , where D/Dt represents a total or material derivative, (cid:126)v = (cid:126)v ( t, (cid:126)x, u ) is the advectionvelocity, u = u ( t, (cid:126)x ), ˜ N represents a general nonlinear term and the gradient ( ∇ ) actsonly on the spatial variables ( ∇ = ( ∂ x , ∂ x , ..., ∂ x n )). ∗ Submitted to the editors 9th of August, 2018.
Funding:
This work was funded by FAPESP grants 2014/10501-0 and 2016/18445-7, USP-Santander Mobility grant 527/2016, Shell Brasil under the ANP R&D levy grant 20714-2 and NCARfunding. † ∼ pedrosp). ‡ a r X i v : . [ phy s i c s . c o m p - ph ] J u l P. S. PEIXOTO, M. SCHREIBER
The treatment of the nonlinear advection in exponential integrators varies andleads to different mathematical properties. It can, for instance, be simply thoughtas a nonlinear term in the exponential integration scheme [9]. Or else, the nonlinearterm can be treated via a linearization procedure [15, 35, 10, 50, 24, 31].A well-established method to solve equations with nonlinear advection is the semi-Lagrangian advection approach [44, 53, 19], sometimes denoted as the characteristicsmethod [41, 8, 7]. The cost-effectiveness of semi-Lagrangian schemes depends on theproblem [5]. They are used in computational fluid dynamics [57, 13], and are verysuccessfully used in weather forecasting [55], hence being adopted by several weatherforecasting centres in operational models [17, 4, 22, 39].Exponential integrators and semi-Lagrangian schemes have an interesting connec-tion. For linear advection, the characteristics (which define a particle trajectory) areprecisely given by the exponential of the linear advection operator [11]. Moreover, fornonlinear advection, it is possible to establish an equivalence between the solution ofa general integration factor problem to a semi-Lagrangian approach [52]. Therefore,it is possible to obtain properties of semi-Lagrangian schemes considering them froman exponential operator point of view. Or, similarly, it is possible to consider thesolution of a semi-Lagrangian problem in place of an operator exponential [12]. Thelatter allows, for example, the development of high order semi-Lagrangian schemes[13].The goal of this work is to further explore a combination of semi-Lagrangian andexponential integrators. The key development in this paper is to consider an expo-nential integration scheme that is built with respect to the total (material) derivative,therefore treats nonlinear advection within the exponentiation framework, which, toour knowledge, has not yet been explored in the literature. With this methodol-ogy, nonlinear advection is calculated accurately with low dispersion error (propertyearned from the semi-Lagrangian approach), in combination with an accurate solutionof the linear problem even for very stiff problems (property earned from the expo-nential integration). In principle, several combinations of exponential integration andsemi-Lagrangian schemes could be explored. We will derive the general principles ofthe method and then illustrate how well-established schemes can be used together.Traditional geophysical fluid dynamics models usually employ either an explicittime stepping scheme, for which the time-step sizes are constrained by faster wavesin the system (e.g inertia-gravity), or implicit time stepping schemes (e.g. Crank-Nicolson), which allow larger time-steps, at the cost of slowing the faster (short wave-length) linear waves. A recent review on the matter of time stepping schemes forweather and climate [39] points out the need of time integration schemes that allowlarge time-steps while preserving wave dispersion properties. Small scale horizontalgravity waves play an important role in the large structure of the middle atmosphere,particularly for climate simulations [38]. Exponential integrators provide a way toobtain large time-steps without affecting these small-scale waves, preserving superiorlinear dispersion properties (see [47, 14]).An important model for the atmosphere and ocean dynamics is formed by thetwo-dimensional nonlinear rotating shallow water equations (SWE), as they providea simple set of equations that already carry many of the complications encounteredin full three-dimensional dynamics. Recent works of [15] and [25] explored the useof exponential integrators in SWE and showed its potential and practical relevanceto weather forecasting. They explored the dynamic linearization procedure of [54] toobtain their exponential integrator, and the nonlinear advection was treated withinthe linearization. Also within this application framework, [24] shows results from
EMI-LAGRANGIAN EXPONENTIAL INTEGRATION
2. Exponential integration.
We start providing a brief review of some existingexponential integration techniques that will be relevant for the semi-Lagrangian expo-nential approach. More details may be found in the review of exponential integratorsof [29] and in references therein.
Numerically, the solution of equation (1), u ( t ), is approximated by ( n ) discrete values that could be, for example, grid pointvalues or spectral coefficients. This defines the discrete solution U ( t ) ∈ R n evolving intime. The linear operator ( L ) can be approximated by a discrete version of it ( L ), witha preferred discretization scheme. Since L may be originated from a partial differentialequation problem, it is prudent to keep in mind that L may carry information fromspatially varying features. However, having derived it for an autonomous system, itis independent of time. So the analogous semi-discrete problem of interest may bewritten as(3) dU ( t ) dt = LU ( t ) + N ( U ( t )) , U (0) = U , where L ∈ R n × R n is the discrete linear operator (an n × n matrix) and N ( U ) is adiscrete version of N ( u ).Now let’s assume that U ( t n ) is given for a current time t n , and that we wishto calculate U ( t n +1 ), for t n +1 = t n + ∆ t . Since L does not depend on time, the P. S. PEIXOTO, M. SCHREIBER integration factor problem,(4) dQ n ( t ) dt = − Q n ( t ) L, Q n ( t n ) = I, where I is the identity matrix, has a unique solution given by(5) Q n ( t ) = e − ( t − t n ) L . Using the integration factor in equation (3) one sees that(6) ddt ( Q n ( t ) U ( t )) = Q n ( t ) N ( U ) . Therefore the problem has an exact solution which may be implicitly represented as,(7) U ( t n +1 ) = e ∆ tL U ( t n ) + e ∆ tL (cid:90) t n +1 t n e − ( s − t n ) L N ( U ( s )) ds, where we note that Q − n ( t ) = e ( t − t n ) L is the inverse of Q n ( t ). This last equation iswell-known as the variation-of-constants formula. Exponential integration makesuse of calculations of the exponentials, and/or exponential related functions, to obtaina time marching scheme along the lines of equation (7).The main difference between variants of exponential integration schemes is in theway the nonlinear term is evaluated. If the equation is purely linear ( N = 0), then theintegral term in equation (7) vanishes and it is possible to solve the problem directlyfrom the matrix exponential calculation for each time-step. For nonlinear problems,there exists several approaches (see [29]). We will use the Runge-Kutta ExponentialTime Differencing (ETDRK) as a representative method, following [16]. However,for the semi-Lagrangian exponential scheme (to be shown), other methods could beconsidered in a similar fashion.As a first order approximation, let the nonlinear term N ( U ) in the integral beconstant in time, for each time-step, with value N ( U ( t n )). Using equation (4) andassuming L − exists, we may then formally derive what is known as the first orderETD1RK method,(8) U n +1ED1 = ϕ (∆ tL ) U n + ∆ t ϕ (∆ tL ) N ( U n ) , where,(9) ϕ ( z ) = e z , ϕ ( z ) = z − ( e z − . with z = ∆ tL .More general (higher order) ETD schemes may be derived using higher order ϕ k functions (see [16]), which may be defined from the recurrence relation(10) ϕ k +1 ( z ) = z − ( ϕ k ( z ) − ϕ k (0)) , where potential singularities may be treated by a series expansions for ϕ n evaluationsclose to their singularity.We will be particularly interested in this paper in the second order ETD2RKscheme, in order to allow a fair comparison to other well-established second order EMI-LAGRANGIAN EXPONENTIAL INTEGRATION U n be the numerical approximationof U ( t n ) at time t , then the ETD2RK scheme may be written, as function of theETD1RK scheme, as(11) U n +1ED2 = U n +1ED1 + ∆ tϕ (∆ tL ) (cid:0) N ( U n +1ED1 ) − N ( U n ) (cid:1) , which is derived substituting the second order approximation for the nonlinear term,(12) N ( U ( s )) = N ( U ( t n )) + ( s − t n )∆ t ( N ( U ED1 ( t n +1 )) − N ( U ( t n ))) + O (∆ t ) , into equation (7).
3. Semi-Lagrangian integration.
Eulerian schemes usually keep a fixed gridand evaluate the movement of the particles that pass through a computational cell.For nonlinear advection, these schemes usually have time-step size limited by theCourant-Friedrichs-Lewy condition (CFL). Lagrangian schemes usually follow parti-cle trajectories (characteristics) through time and may not even rely on a fixed com-putational grid, or else have a grid evolving over time. This can create complicatedgrid structures involving, for example, intersections of trajectories. Semi-Lagrangianschemes compensate for this by keeping a fixed grid, but following the particle tra-jectories for a single time-step (a local version of the classical Lagrangian approach).Since the trajectories may end, or start, in points not on the reference grid, usuallyan interpolation step is required.In the context of atmospheric simulations, this scheme usually allows time-stepsizes significantly larger than CFL-restricted Eulerian schemes [45]. As we will ob-serve, maintaining this property of significantly larger time-step sizes is a non-trivialtask with exponential integration.In this section we introduce classic notations and results about semi-Lagrangianschemes. This will be required as a basis to derive the semi-Lagrangian exponentialschemes in the next section. Further details on semi-Lagrangian methods can befound in [53] and [19].
We start considering Equation (2) on a La-grangian framework, relative to a particle initially positioned at (cid:126)r in space. Thus,the system state is formed by u = u ( t, (cid:126)r ( t )), with advection velocity defined as (cid:126)v = (cid:126)v ( t, (cid:126)r ( t ) , u ( t, (cid:126)r ( t ))). Here, (cid:126)r ( t ) is the Lagrangian trajectory of the particle, there-fore it is the solution of the non-autonomous problem(13) d(cid:126)r ( t ) dt = (cid:126)v ( t, (cid:126)r ( t ) , u ( t, (cid:126)r ( t ))) , (cid:126)r (0) = (cid:126)r . Equation (2) may be written in a Lagrangian framework as(14) du ( t, (cid:126)r ( t )) dt = ∂u ( t, (cid:126)r ( t )) ∂t + (cid:126)v · ∇ u ( t, (cid:126)r ( t )) = L ( u ( t, (cid:126)r ( t ))) + ˜ N ( u ( t, (cid:126)r ( t ))) , with initial condition u (0 , (cid:126)r ) = u , where now L and ˜ N may implicitly also depend onthe position (cid:126)r ( t ). This particular time derivative ( d/dt ) on the Lagrangian frameworkis usually denoted as a total (material) derivative D/Dt , as in equation (2).As in the previous section, we will focus here on a discretized problem, where L may be again directly viewed as a finite dimensional matrix operator, hence linear, P. S. PEIXOTO, M. SCHREIBER and will be denoted by L . In a Lagrangian framework, L may depend on the particleposition (cid:126)r ( t ). Therefore, we may analogously to equation (3) set the general non-autonomous semi-discrete problem to be(15) DU ( t, (cid:126)r ( t )) Dt = L ( U ( t, (cid:126)r ( t ))) + ˜ N ( U ( t, (cid:126)r ( t ))) , U ( t , (cid:126)r ( t )) = U , where ˜ N is a numerical approximation to ˜ N , and now the time differential is a totalderivative and depends on the solution of the problem given in (13).Although there are many forms of semi-Lagrangian schemes [53], these usuallyrely on basically two components with both playing important roles in the accuracyand stability of the schemes [21, 19, 40].Component (i): interpolation of the information to the reference grid. As shownin [21], the interpolation order needs to be chosen in agreement with the accuracyorder of the trajectory calculation.Component (ii): the evaluation of trajectories, which are solutions of the problem(13).We will consider a back-trajectory approach, which is a well-established approach[30] that assumes that the grid is fixed at time t n +1 . The trajectory determines theposition of a departure point at time t n , which is likely not to be a grid point, so aninterpolation of the advected quantity is required. The trajectory evaluation itself canthen be obtained by a direct numerical time integration of differential equation (13),as a sub-cycling procedure, or, which is more common in atmospheric applications,iteratively solve its integral form. In the later, one may obtain the departure point (cid:126)r d = (cid:126)r ( t n ) from the knowledge of the arrival point (cid:126)r a = (cid:126)r ( t n +1 ), which is set to be agrid point, using two-time level schemes [37]. This can be done using the midpointrule integration (for (cid:126)r m = (cid:126)r ( t n +1 / )) and an iterative procedure to solve the nonlinearresulting equation. An example of such procedure, which will be used in the presentwork, uses the following iterative equation,(16) (cid:126)r k +1 m = (cid:126)r a − (cid:126)v ( t n +1 / , (cid:126)r km , u km ) ∆ t , where u km = u ( t n +1 / , (cid:126)r km ) and usually (cid:126)r m = (cid:126)r a is used as initial step for the procedure.The departure point is then obtained by simply considering (cid:126)r d = 2 (cid:126)r m − (cid:126)r a .In case (cid:126)v is not known within [ t n , t n +1 ], for example if (cid:126)v depends on u , its evalua-tion in intermediate times requires an extrapolation from previous time-steps. This ex-trapolation may directly influence the stability of the scheme [19]. A well-establishedapproach is the Stable Extrapolation Two-Time-Level Scheme (SETTLS) of [30], usedat the ECMWF in their global spectral model IFS for operational weather forecasts.All trajectory calculations used in this work follow the SETTLS scheme to obtain thedeparture points. An important scheme foratmospheric modeling is the one used in the IFS-ECMWF model. It uses a semi-Lagrangian scheme coupled with a semi-implicit time stepping of linear terms withspectral horizontal discretization. This scheme, based on [30], will serve as a firstguideline in the development of the semi-Lagrangian exponential schemes. European Centre for Medium-Range Weather ForecastsEMI-LAGRANGIAN EXPONENTIAL INTEGRATION U n +1 − U n ∗ ∆ t = 12 (cid:0) ( LU ) n +1 + ( LU ) n ∗ (cid:1) + ˜ N n +1 / , where the subscript ∗ with superscript n denotes interpolation to departure points ( (cid:126)r d )[6] and the last term represents the non-linearities at the midpoint of the trajectory.This term is computed based on averaging and extrapolation (see [30], Eq. (4.4,4.5))with(18) ˜ N n +1 / = 12 (cid:16) (cid:104) N n − ˜ N n − (cid:105) ∗ (cid:124) (cid:123)(cid:122) (cid:125) Extrapolation to t n +1 + ˜ N n (cid:17) , which is the SETTLS extrapolation, where ˜ N n is the evaluation of the nonlinearterm at time t n . The unknowns in Equation (17) are implicitly given by U n +1 and( LU ) n +1 , requiring a linear solver.To ensure an overall second order accurate scheme (assuming ∆ t ∝ ∆ x ), it issufficient to use cubic interpolations of the advected quantities (with respect to Equa-tions (17) and (18)), and linear interpolations of the velocities in the iterative processof trajectory calculations [40].
4. Semi-Lagrangian exponential integration.
In this section, we discusshow the general exponential integration techniques can be applied in a Lagrangianreference frame, to derive the novel semi-Lagrangian exponential methodology.
The key concept investigated in this paper is to consider,from a numerical perspective, the exponential integration of Equation (15) consideringthe total (material) derivative.As in Section 2, where we built exponential integration schemes from the solutionof an integration factor problem, we would like to be able to define a similar integra-tion factor for the problem with respect to this material derivative. We assume theexistence of an integration factor P n ( t ) that is a solution to the problem(19) D ( P n ( t ) U ( t, (cid:126)r ( t ))) Dt = P n ( t ) ˜ N ( U ( t, (cid:126)r ( t ))) , P n ( t n ) = I. Assuming that U is a solution of (15), P n will also be a solution of(20) DP n ( t ) Dt U ( t, (cid:126)r ( t )) = − P n ( t ) L ( U ( t, (cid:126)r ( t ))) , P n ( t n ) = I. We recall that L may depend on the spatial variations, which are now depending alsoon time due to the Lagrangian framework, so we will explicitly indicate this with asubscript as L = L (cid:126)r ( t ) . If L (cid:126)r ( t ) commutes in time, that is, L (cid:126)r ( t ) L (cid:126)r ( s ) = L (cid:126)r ( s ) L (cid:126)r ( t ) forall times s and t , then the integration factor problem has a solution given by(21) P n ( t ) = e − (cid:82) ttn L (cid:126)r ( s ) ds . For the continuous problem with L based on space-varying coefficients (dependentof the particle position), the commutation assumption of L will most likely not besatisfied. The integration factor may, however, still exist and be well defined, butmight not have the usual matrix exponential form. Assuming that such an integration P. S. PEIXOTO, M. SCHREIBER factor exists, and that it is invertible ( P − n exists for all time), equations (15) and(19) indicate the following implicit relation on U (analogous to (7)),(22) U ( t n +1 , (cid:126)r ( t n +1 )) = P − n ( t n +1 ) U ( t n , (cid:126)r ( t n )) + P − n ( t n +1 ) (cid:90) t n +1 t n P n ( s ) ˜ N ( U ( s, (cid:126)r ( s ))) ds. This is the fundamental equation for the derivation of the semi-Lagrangian exponentialschemes developed in this paper.Numerically, one needs an explicit way of calculating the integration factor. Thiswill depend on the problem of interest. One possibility is to directly numericallyintegrate equation (20), which is the basis of many operator splitting techniques [52].Another possibility, if such integration factor is unknown in its exponential form,is to assume that L does not vary within each time-step for each given local trajectory,since then the problem reduces to a matrix exponential problem. This assumption canbe also thought as the fundamental idea of semi-Langrangian methods which we willapply later: interpolating the Lagrangian solution to a regular Eulerian grid allowsapplying the linear and nonlinear operators on a regular grid. This should provide afirst order approximation to the true integration factor at each time-step. This greatlysimplifies the problem, as in this case P n = Q n , as defined in equation (5), and theproblem reduces to(23) U ( t n +1 , (cid:126)r ( t n +1 )) = e ∆ tL U ( t n , (cid:126)r ( t n )) + e ∆ tL (cid:90) t n +1 t n e − ( s − t n ) L ˜ N ( U ( s, (cid:126)r ( s ))) ds + E L , where E L denotes the potential errors introduced with the assumption of L beingconstant along the trajectories. This is almost identical to what we obtained for theusual exponential integration approach (see equation (7)), but now U is varying alonga particle trajectory in time, resulting in a derivation of what we are calling a semi-Lagrangian exponential integration. Using the semi-Lagrangian notation, we rewritethe numerical method from equation (23) as(24) U n +1 = e ∆ tL U n ∗ + e ∆ tL (cid:90) t n +1 t n e − ( s − t n ) L ˜ N ( U ( s, (cid:126)r ( s ))) ds + E L , where U n +1 is given at grid points and U n ∗ refers to the (interpolated) value at de-parture points.Finally, different semi-Lagrangian exponential schemes can be built depending onhow the integral is approximated, as is the case with the usual exponential integrationtechniques. The integral consists of an application of a linear operator, which we willcompactly write as T ( s ) = e − ( s − t n ) L , and a nonlinear function calculation, compactly w ( s ) = ˜ N ( U ( s, (cid:126)r ( s ))). Interestingly, a method as simple as calculating a midpointrule integration can have different forms, and naive formulations may easily generateinconsistent schemes.We will discuss this issue by first pointing out two very important remarks whichwill play a crucial role for the derivation of a consistent and stable semi-Lagrangianexponential integration scheme:(R1) Integration : The integral term relies on a linear operator, T , acting on anonlinear function, w , integrated along a trajectory. If we wish to evaluatethe term T w at time t n − k , with k ≥
0, at departure points (or trajectorymidpoints), we should first apply the linear operator to the nonlinear func-tion evaluated at grid points at time t n − k , and only then interpolate to the EMI-LAGRANGIAN EXPONENTIAL INTEGRATION at time t n − k , the application of the linear oper-ators should come before the interpolation operation , as, for example, with( T ( t n − k ) w ( t n − k )) ∗ .(R2) Advected quantities : At time t n +1 , interpolated values of quantities comingfrom times t n − k , k ≥
0, are assumed to have already been advected, thereforethe results lay on a regular grid relative to the arrival points. Consequently, attime t n +1 , for advected quantities, the linear operator should be applied afterthe interpolation operation , as, for example, with e ∆ tL U n ∗ term of equation(24).These important remarks are a peculiarity of this kind of semi-Lagrangian expo-nential formulation. They are required since the interpolation operation of advectedquantities in general does not commute with a linear operator. That is, even thought e ∆ tL is a linear operator, possibly independent of time and space, it does not in gen-eral commute with the interpolation operator ( ∗ ), since this interpolation reflects anon-regular grid formed by nonlinear backward trajectories. Therefore, in general, e ∆ tL U n ∗ (cid:54) = (cid:0) e ∆ tL U n (cid:1) ∗ . We provide, in Appendix A, an illustration of this lack ofcommutation, which interestingly happens even in the case of linear advection.Returning to the case of the midpoint rule integration of the nonlinear term, wemay have, for example, these 2 distinct approximations of the integral term,(25) (cid:90) t n +1 t n T ( s ) w ( s ) ds ≈ (cid:40) A = ∆ t T ( t n +1 / ) (cid:2) w ( t n +1 / ) (cid:3) † A = ∆ t (cid:2) T ( t n +1 / ) w ( t n +1 / ) (cid:3) † , where we used the † symbol to indicate the interpolation to the midpoints of thetrajectories (whereas the interpolation to departure points are indicated with ∗ ).In A , first the nonlinear function ( w ) is evaluated at the intermediate time( t n +1 / ) on a regular grid and then this is used to obtain interpolated values atthe trajectories midpoints. Only then the operator T , evaluated at time t n +1 / ,is applied. Since the operator is applied on values given at trajectories midpoints,this application happens on a possibly irregular grid, therefore may results in aninconsistent application of the operator.In A , first both the operator T and the nonlinear function w are evaluated atthe intermediate time ( t n +1 / ), on a regular grid, then the operator is applied on w .Only then interpolated values are obtained for the midpoints of the trajectories. Inthis case the linear operator is always applied on regular grid quantities.As an example, consider the following scheme,(26) U n +1 = e ∆ tL U n ∗ + ∆ t e ∆ t L ˜ N n +1 / , (unstable scheme) . where ˜ N n +1 / is an approximation of the nonlinear term at the trajectory midpoint,for example considering equation (18). Due to the aforementioned remarks and dis-cussion, this scheme applies the linear operator on conceptually irregular grids, andtherefore may be inconsistent with the underlying equations. Experiments with thismethod revealed it is unstable even with small time-step sizes. Fol-lowing the SETTLS scheme [30] for the semi-Lagrangian discretization, but using itwith respect to equation (24), we may derive our first combination of semi-Lagrangian0
P. S. PEIXOTO, M. SCHREIBER exponential scheme, which we will denote as SL-EXP-SETTLS. The scheme is derivedfrom (24) as(27) U n +1SLEX = e ∆ tL U n ∗ + ∆ t e ∆ tL ˜ N n +1 / e , where we use the SETTLS extrapolation to obtain the value of ˜ N at the trajectorymidpoint as(28) ˜ N n +1 / e = 12 (cid:104) N n − e ∆ tL ˜ N n − (cid:105) n ∗ + 12 ˜ N n . We note that ˜ N n +1 / e is an approximation of e − ( s − t n ) L ˜ N ( U ( s, (cid:126)r ( s ))) at the midpointof the trajectory taking into account the remark (R1). One may note that it differsfrom the method presented in Equation (26), as in this case the linear operator andnonlinear calculations are not split, but computed jointly at desired time-steps.It is possible to simplify the above equations in order to require only 2 expo-nential evaluations per time-step. This scheme is a multi-step scheme, that requiresinformation from two previous time-steps. This scheme may also be thought as asemi-Lagrangian version of the Integrating Factor method, proposed in [16] as thesecond order Adams-Bashforth Integrating Factor method (IFAB2), as one can noticefrom their equation (31).As discussed in [16], the concept of stability for Integrating Factor methods isunclear. This is also the case for our semi-Lagrangian version of exponential schemes.Therefore, this is a topic we discuss in this paper purely from a numerical perspective. To constructsemi-Lagrangian Exponential Time Differencing Runge-Kutta schemes (SL-ETDRK)in analogy to usual ETDRK schemes, we need to pay attention to the remarks (R1)and (R2) above. In usual ETD schemes, as shown in Section 2, the exponential infront of the integral in equation (24) would commute with the integral, to be placedwithin the integrand. However, since now the integral is along trajectories, this nolonger results in an equivalent problem in the numerical scheme, due the remarkspointed out above. Therefore, we should first evaluate the integral term, and thenapply the linear operator ( e ∆ tL ).Following this strategy, we may derive the semi-Lagrangian ETD1RK scheme inthe following way. From equation (23), assuming as in ETD1RK that the non-linearityis constant within a time-step, we have(29) U n +1SLED1 = ϕ (∆ tL ) (cid:104) U n + ∆ t ϕ ( − ∆ tL ) ˜ N ( U n ) (cid:105) n ∗ , where we used that ϕ ( z ) = ϕ ( z ) ϕ ( − z ). This scheme can be computed numericallywith only two ϕ function evaluations and one interpolation per time-step.Deriving the second order scheme (SL-ETD2RK) involves a more careful analysisof how the integral in equation (23) is approximated. Let N ( U ( s )) = N ( U ( t n , (cid:126)r ( t n ))) +( s − t n )∆ t ( N ( U SLED1 ( t n +1 , (cid:126)r ( t n +1 ))) − N ( U ( t n , (cid:126)r ( t n )))) + O (∆ t ) , (30) EMI-LAGRANGIAN EXPONENTIAL INTEGRATION U SLED2 ( t n +1 , (cid:126)r ( t n +1 )) = ϕ (∆ tL ) U ( t n , (cid:126)r ( t n ))+ ϕ (∆ tL ) (cid:18)(cid:90) t n +1 t n e − ( s − t n ) L ds (cid:19) N ( U ( t n , (cid:126)r ( t n )))+ ϕ (∆ tL ) (cid:18)(cid:90) t n +1 t n ( s − t n )∆ t e − ( s − t n ) L ds (cid:19) N ( U ( t n +1 , (cid:126)r ( t n +1 ))) − ϕ (∆ tL ) (cid:18)(cid:90) t n +1 t n ( s − t n )∆ t e − ( s − t n ) L ds (cid:19) N ( U ( t n , (cid:126)r ( t n ))) . (31)To be able to preserve e ∆ tL outside of the integral, and still make use of the ϕ functions of ETDRK schemes, we may factor out the ϕ ( z ) = e z function of the ϕ k functions. Let ψ k functions be defined as(32) ψ k ( z ) = ( − k +1 ϕ k ( − z ) + k − (cid:88) l =1 ϕ l ( − z ) . It can be shown that ϕ k ( z ) = ϕ ( z ) ψ k ( z ) by substituting equation (10) into the right-hand-side of the definition of ψ k and using binomial expansions in a similar way asdone in [16].Using the SL-ETD1RK scheme and the properties of the ϕ functions with respectto ψ functions, we may write the SL-ETD2RK scheme as(33) U n +1SLED2 = U n +1SLED1 + ∆ t ϕ (∆ tL ) (cid:2) ψ (∆ tL ) N ( U n +1SLED1 ) − ( ψ (∆ tL ) N ( U n )) n ∗ (cid:3) . After suitably rearranging the equations, the scheme can be coded to require 4evaluations of ϕ (or ψ ) functions and 2 interpolations.
5. Rotating Shallow Water Equations on f-Plane.
In this section we reviewthe basic concepts of the Shallow Water Equations (SWE), which will serve as anapplication for the schemes discussed in the previous sections.Considering a Lagrangian framework, with particle trajectories given by (cid:126)r ( t ) =( x ( t ) , y ( t )) on a plane, we define (cid:126)v = (cid:126)v ( t, (cid:126)r ( t )) = ( u ( t, (cid:126)r ( t )) , v ( t, (cid:126)r ( t ))) to be the flowvelocity, and η = η ( t, (cid:126)r ( t )) a fluid depth perturbation about a constant mean fluiddepth (¯ η ). The rotating SWE on a plane may then be written as(34) DUDt = L U + ˜ N ( U ) , where the time derivative is assumed to be the total (material) derivative, and(35) U = uvη , L = f − g∂ x − f − g∂ y − ¯ η∂ x − ¯ η∂ y , ˜ N ( U ) = − η ∇ · (cid:126)v , where the total fluid depth h is given by h = η + ¯ η . The gravity g and the Coriolisparameter f are assumed to be constant throughout this paper (f-plane approxima-tion). Initial conditions for the prognostic variables ( u, v, η ) are assumed to be given.Bi-periodic boundary conditions will be adopted for ( x, y ) on a rectangular limitedset of R . The bottom topography considered in this work is assumed to be flat.2 P. S. PEIXOTO, M. SCHREIBER
Well-established models adopt semi-implicit schemes [19, 45], with implicit treat-ment of linear terms and explicit treatment of non-linearities. Among the implicitschemes for the linear waves, Crank-Nicolson (trapezoidal differencing) is frequentlyadopted, as done for example in the IFS model of the ECMWF [20, 30], coupled with asemi-Lagrangian approach. Modern models that use non-regular spherical grids, suchas MPAS [51] or DYNAMICO [18], adopt explicit time stepping procedures based onRunge-Kutta time integration. See [39] for an extensive list and description of themain time stepping schemes used for weather and climate models.
We seek to find the exponentialof the linear operator L where we assume the time-step size ∆ t incorporated into L by simple scaling. Assuming a double Fourier expansion of U in space on a [0; 2 π ) periodic domain, we can look at a single mode (single wavenumber) to understandthe action of L in terms of its exponentials. For a fixed time, let U be of the form U (cid:126)k ( (cid:126)x ) = e i(cid:126)k · (cid:126)x ˆ U (cid:126)k , with (cid:126)k = ( k , k ), (cid:126)x = ( x , x ) = ( x, y ), ˆ U (cid:126)k independent of (cid:126)x and i = √−
1. Then(36) L U (cid:126)k = f − gik − f − gik − ¯ ηik − ¯ ηik ˆ U (cid:126)k , where the matrix above is the matrix symbol of L (usually denoted as L ( i(cid:126)k )), withpurely imaginary eigenvalues [36],(37) ω f ( (cid:126)k ) = 0 , ω g ( (cid:126)k ) = ± i (cid:113) f + g ¯ η (cid:126)k · (cid:126)k, where ω f ( (cid:126)k ) is the steady geostrophic mode and ω g defines the 2 inertia-gravity wavemodes ( ω − g ( (cid:126)k ), ω + g ( (cid:126)k )). The eigenvectors can be directly computed from L ( i(cid:126)k ),yielding a matrix of eigenvectors Q . Writing the eigenvalues as a diagonal matrixΛ = [ ω f ( (cid:126)k ) , ω − g ( (cid:126)k ) , ω + g ( (cid:126)k )], and using L ( i(cid:126)k ) = Q Λ Q − , the exponential of L can bedirectly calculated for the shallow water system through its symbol as(38) e L ( i(cid:126)k ) = Qe Λ Q − , where the e Λ is the diagonal matrix with entries given by the exponential of therespective eigenvalues.For the studies conducted in the present work, we exploit features from doubleFourier spectral spatial discretization. This allows us to compute the numerical matrixexponential directly from equation (38), providing an exponential ( ϕ ) of the linearoperator accurate to machine precision. To evaluate ϕ n (∆ tL ) functions (see Eq. (10)),we use that in spectral space(39) ϕ n (∆ t L ( i(cid:126)k )) = Qϕ n (∆ t Λ) Q − , hence computing ϕ n element-wise for each diagonal element in Λ.We would like to emphasize that computing the exponential directly as discussedabove is only possible because we are exploiting the orthogonality of Fourier basis onthe bi-periodic domain. For the SWE on the sphere, we’d like to mention that a similarfeature was recently also discovered by using spherical harmonics to exponentially timeintegrate gravity modes. EMI-LAGRANGIAN EXPONENTIAL INTEGRATION
The linear SWE on an f-plane de-fine a hyperbolic system formed by inertia-gravity (Poincar´e) and geostrophic (steady)waves, with the dispersion relations described in the previous section ( ω g and ω f , re-spectively). Numerical schemes should be able to represent well these two kinds ofwaves. We will adopt in this study spectral spatial discretizations of the linear opera-tor (based on Fourier series), therefore errors in the evaluation of the linear operatorare negligible (of machine precision) for each wavenumber. However, the temporaldiscretization may still be a source of error which can be directly investigated.Linear exponential integration schemes do not introduce any errors in time if thelinear operator and its exponential is calculated analytically. However, state-of-the-artweather forecasting systems do not usually adopt exponential integration schemes, butmostly Runge-Kutta schemes [51] when explicit, or Crank-Nicolson (CN) [30] whenimplicit (see a complete description in [39]). To ensure large time-steps, implicitschemes are preferred, but, in this case, the dispersion relations are usually not veryaccurately attained for the smaller wave-modes (faster gravity waves). Durran [19]shows that the approximate dispersion (˜ ω ) of the CN scheme preserves the steadygeostrophic modes (for ˜ ω f ( (cid:126)k ) = ω f ( (cid:126)k ) = 0). However, the gravity waves will havedispersion of the form,(40) ˜ ω g ( (cid:126)k ) = ω g ( (cid:126)k ) + ∆ t
12 ( ω g ( (cid:126)k )) + O (∆ t ) , which is purely imaginary (the amplitude of the mode is not altered by the scheme),but the phase speed is affected. The odd powers of ω g indicate that the additionalterms (error) will always produce a reduction of the ˜ ω g frequency, and this reductionwill be larger the larger the wavenumber norm ( (cid:126)k · (cid:126)k ), since it depends on ω g ( (cid:126)k ).Therefore, the error in the Crank-Nicolson method slows down the faster (largerwavenumber) inertia-gravity waves, which will be slower when larger time-step sizesare used.For finite difference schemes the spatial errors significantly influence the dispersionrelations. [43] analyzes the effect of different discretizations on the shallow water wavesdispersions. To preserve an adequate representation of the inertia-gravity waves andreduce computational modes arising from spatial discretizations, staggered grids arepreferred. These are usually called C-grids in the geoscientific modeling community,and have the depth variable centered in the cell and the velocities given at the edges4 P. S. PEIXOTO, M. SCHREIBER of cell, normal to the edge [2]. Since many modern atmospheric models that use non-regular C- grids are using finite-differences/volumes with explicit time integration, wewill also consider this approach as reference in our experiments further in the paper.
6. Numerical experiments.
We will consider the following set of schemes tobe analyzed: • RK-FDC: Runge-Kutta second order in time with second order in space en-ergy conserving finite differences discretization on a staggered C-grid [46]. • SL-SI-SETTLS: Semi-Lagrangian, semi-implicit (Crank-Nicolson) scheme us-ing spectral discretization adapted from [30] to the plane. • SL-EXP-SETTLS: Exponential version of SL-SI-SETTLS (see Section 4.2). • ETD2RK: Original ETD2RK scheme with spectral space discretization (seeSection 2). • SL-ETD2RK: Semi-Lagrangian version of ETD2RK (see Section 4.3). • REF: Reference solution. Runge-Kutta forth order in time with small time-step and high resolution Eulerian spectral discretization (pseudo-spectral forall nonlinear terms, such as advection).The schemes are connected in the following way. RK-FDC is a reference explicitscheme well-established for the solution of the SWE of very low cost per time-step,but restricted to smaller time-steps (CFL condition). SL-SI-SETTLS is the state-of-the-art scheme used in many global atmospheric dynamical cores, which we aim tocompare to our semi-Lagrangian exponential schemes (SL-EXP-SETTLS, SL-EXP-ETD2RK). ETD2RK is a well-established exponential integration technique, whichwe aim to compare to our semi-Lagrangian version, SL-ETD2RK, considering thedifferent treatment of the nonlinear advection.
The experiments will be exe-cuted on a scenario that mimics the Earth’s dimensions, and we will follow the stan-dard spherical test case parameters defined in [56]. The domain is set to be [0 , L x ] × [0 , L y ] = [0 , πa ] × [0 , πa ], where a = 6371 .
22 km is the Earth radius, with bi-periodicboundary conditions. The gravity acceleration constant is set to g = 9 . − and the Coriolis frequency constant is f = 2Ω, with Ω = 7 . × − rad · s − . Themean depth is ¯ η = 10 km so that the gravity wave speed is c = √ g ¯ η ≈
313 ms − .The experiments will be performed with a horizontal discretization of 512 spectralmodes in each dimension. This corresponds to 768 physical grid points to avoidaliasing effects, which would result in a grid cell with a length of approximately 52 kmin each coordinate. The exception is the reference solution (REF), for which we willuse 1024 spectral modes per coordinate. Such high horizontal resolution was chosenin order to reduce the errors relative to spatial discretizations and allow a clearercomparison of the different time stepping schemes. The time-step sizes will varyaccording to the analysis to be investigated.We will present results of errors in two metrics: maximum absolute error (Max-Error) and root mean square error (RMSError), always for fixed integration time(timestamp). In case of mismatching resolutions, where point-wise comparison is notwell defined, bi-cubic spline interpolation is used on the highest resolution result toobtain information on the lowest resolution grid. This lack of matching happens as weare using a collocated grid (A-grid in geophysical notation) for the reference solution(REF), with physical representation of the quantities considered in the center of thecell. EMI-LAGRANGIAN EXPONENTIAL INTEGRATION The analysis of the energy spectra is deeplyrelated to the study of turbulence in fluid dynamics models, which is well investigatedfor the atmosphere (e.g. [34, 32]). Here, we do not intend to do turbulence analysis, butrather use spectrum analysis to compare small-scale wave interactions of the differentschemes. Therefore, we will assume a simplified kinetic energy spectrum analysis.The one-dimensional Discrete Power Density Spectra ([42]) is obtained from thetwo-dimensional kinetic energy spectrum using the Fourier transformed velocities(ˆ u ( (cid:126)k ) , ˆ v ( (cid:126)k )), as(41) E n = (cid:88) n ≤(cid:107) (cid:126)k (cid:107)
0. Energy builds up in even wavenumbers dueto nonlinear interactions. Note also that the spectra converges towards the well known − / EMI-LAGRANGIAN EXPONENTIAL INTEGRATION y ( k m ) Vorticity REF t=30.0 days dt=000002 -4e-05-2e-050e+002e-054e-05 / s (b) Horizontal wavelength (km)10 (11 (9 (7 (5 (3 (1 K i n e t i c E n e r g y S e c t r u m ( m s ( ) k (5/3 k (3 Kinetic Energy S ectrum
REF_t1.0REF_t2.0REF_t4.0REF_t8.0REF_t10.0REF_t12.0REF_t20.0
Fig. 3 . (a) Reference solution for vorticity at 30 days using the full nonlinear SWE. (b) Kineticenergy spectrum for reference solution using the full nonlinear SWE for different integration times(from 1 day to 20 days). in small wavelengths stably is usually a major challenge for numerical schemes. Considering ˜ N = 0 simplifies the semi-Lagrangian exponential schemes.In fact, in this case, SL-EXP-SETTLS and SL-ETD2RK are equivalent, since theonly non-linearity left (advection) is treated within the semi-Lagrangian approach.SL-SI-SETTLS also greatly simplifies for similar reasons. RK-FDC, ETD2RK andREF still have to deal with the nonlinear advection as a nonlinear term. The finitedifferences scheme RK-FDC is built about the vector invariant form of the equations,where nonlinear advection is not explicit, therefore it is not clear how to remove thenonlinear divergence and we do not present results of this scheme in this case.(a) Timestep size (sec)10 (cid:1) (cid:1) (cid:1) M a x E rr o r SL-SI-SETTLSSL-ETD2RKSL-EXP-SETTLSETD2RK (b) Timestep size (sec)10 (cid:1) (cid:1) (cid:1) (cid:1) R M SE rr o r SL-SI-SETTLSSL-ETD2RKSL-EXP-SETTLSETD2RK
Fig. 4 . Errors at day 1 of integration of the unstable jet test case without nonlinear divergence.(a) Maximum absolute error and (b) RMS error for day 1 of integration with respect to the referencesolution (REF) for different time-step sizes and schemes. All schemes were tested for all time-stepsizes indicated. If a scheme does not shows a value for a large time-step size it indicates that itbecame unstable for this test. The SWE without nonlinear divergence were adopted in this test,therefore SL-EXP-SETTLS and SL-ETD2RK are identical.
The initial period is dominated by linear gravity waves, so that is where we expectto see benefits of the exponential integration scheme with respect to the semi-implicitscheme. We show in Figure 4 the errors at day 1 of integration for the unstable8
P. S. PEIXOTO, M. SCHREIBER jet test case without nonlinear divergence. It should be noted that for small time-step sizes the dominating error in the semi-Lagrangian schemes becomes the spatialinterpolation errors, not the temporal. These errors for small time-step sizes may bereduced by considering the resolution proportional to the time-step size (∆ x ∝ ∆ t ) orincreasing the accuracy order of the semi-Lagrangian scheme. For smaller time-stepsizes, the ETD2RK is the most accurate one, since ETD2RK has all spatial operatorstreated spectrally, hence no interpolation errors.However, the semi-Lagrangian schemes are stable throughout all time-step sizestested, whereas the ETD2RK scheme is limited by advection CFL time-step size. Ingeneral, the semi-Lagrangian exponential schemes are more accurate than the semi-implicit scheme (SL-SI-SETTLS), due to the accurate treatment of the linear waves.Concluding, the semi-Lagrangian exponential schemes provide a more accurate way,compared to SL-SI-SETLLS, to extend the time-step size allowed by the traditionalexponential scheme (ETD2RK). y ( k m ) Vorticity SL-EXP-SETTLS t= 10.0 days dt= 0225.0 y ( k m ) Vor t icit y SL-EXP-SETTLS t = 1 0 .0 d ay s d t = 0 2 2 5 .0 y ( k m ) Vorticity ETD2RK t=10.0 days dt=0225.0 -4e-05-2e-050e+002e-054e-05 / s y ( k m ) Vorticity ETD2RK t=10.0 days dt=0225.0 -4e-05-2e-050e+002e-054e-05 / s (a) (b) Fig. 5 . Numerical solution of the SWE without nonlinear divergence for the unstable jet testcase at the time of 10 days for the vorticity field using a time-step size of 225 seconds. The boxwithin each figure shows an amplification of the main vortex formation. (a) SL-SI-SETTLS/SL-EXP-SETTLS/SL-ETD2RK (identically looking, only plotted one for sake of brevity), (b) ETD2RK.We can observe more small scale features for the ETD2RK method around the vortex borders.
Due to the dynamically unstable (chaotic) nature of the test case, quantitativeanalysis of errors in longer periods of time is not usually indicated. However, itis interesting to see qualitatively how the schemes behave once the vortices havedeveloped. We show in Figure 5 the vorticity at day 10 for the schemes investigated.All schemes seem to be able to represent well the vortex formation, but we noticethat the ETD2RK has more oscillations at or around the vortices, whereas the semi-Lagrangian schemes show smoother vortices.
In this section, we willanalyze the schemes with respect to the full SWE, including the nonlinear divergence.In this case, the RK-FDC schemes will also be included in the analysis. Also, thedifferent semi-Lagrangian exponential schemes (SL-ETD2RK and SL-EXP-SETTLS)now differ from each other.We start by studying the results of all time integration methods at day 1 of thefull SWE for the unstable jet test case. Convergence plots for varying time-step sizesare presented in Figure 6. As in the previous test, due to the limitation imposed bythe spatial interpolation used in the semi-Lagrangian schemes, the ETD2RK schemeprovides more accurate results for smaller time-step sizes. The ETD2RK scheme isagain restricted in large time-step sizes by CFL condition for advection. The RK-FDC scheme is limited in both time and space: the finite differences scheme limitsthe accuracy, and the gravity wave speed CFL limits the time-step size. With theinclusion of the nonlinear divergence, the SL-EXP-SETTLS scheme turns out to be
EMI-LAGRANGIAN EXPONENTIAL INTEGRATION Timestep size t (sec) L e rr o r o n s u r f a c e h e i g h t h SL-SI-SETTLSETD2RKSL-ETD2RKSL-EXP-SETTLSRK-FDC (b) Timestep size t (sec) L R M S e rr o r o n s u r f a c e h e i g h t h SL-SI-SETTLSETD2RKSL-ETD2RKSL-EXP-SETTLSRK-FDC
Fig. 6 . Errors on surface height at day 1 of integration of the unstable jet test case for thefull nonlinear SWE. (a) Maximum absolute error and (b) RMS error. Same as Figure 4 but nowconsidering the nonlinear divergence. unstable when used with large time-steps. Compared to the SL-SI-SETTLS scheme,the SL-EXP-SETTLS preserves better the high wavenumber gravity waves, whichinteract with each other in the nonlinear divergence and become numerically unstable.Differently, the SL-ETD2RK scheme is stable with large time-steps, and is moreaccurate than the SL-SI-SETTLS scheme, due to the more accurate treatment of thelinear waves. The theoretical stability analysis of the semi-Lagrangian schemes isstill a matter to be investigated and is here considered only in an empirical sense.However, we point out an important difference between them: SL-EXP-SETTLS isa multistep scheme (requires an extrapolation from a previous time-step), whereasthe SL-ETD2RK is a single step method (apart from the extrapolation used in theback trajectory calculation). We also notice that SL-ETD2RK seems to be a viableextension of the ETD2RK scheme to larger time-steps, being more accurate than theSL-SI-SETTLS. Wallclock time (sec) L e rr o r o n s u r f a c e h e i g h t h SL-SI-SETTLSETD2RKSL-ETD2RKSL-EXP-SETTLSRK-FDC
Fig. 7 . Errors vs. wallclock time at day 1 for different time integration methods. Assuming therequest to stay below 100s of wallclock time, we can observe that the SL-SI-SETTLS and SL-EXP-ETDRK methods provide the best results.
So far, we only compared the error with the time-step size of individual timeintegration methods. However, the time-step size does not directly relate to the totalcomputational requirements and therefore also not to the wallclock time. Although0
P. S. PEIXOTO, M. SCHREIBER it is challenging to run representative wallclock times studies which also relate to afully-developed dynamical core running on a super computer, we provide such wall-clock time studies in Figure 7. Small wallclock times indicate larger time-step sizes.Here, we can observe that the exponential semi-Lagrangian scheme SL-ETD2RK hascompetitive wall-clock times compared to the state-of-the-art SL-SI-SETTLS method,particularly for small wall-clock times. There are, nevertheless, additional computingrequirements of the exponentials, but we would like to point out that these are the firsttypes of SL-ETDnRK methods and we expect improved results with further computa-tional, mathematical and modeling optimizations. Again, we also observe a limitationof the accuracy for all semi-Lagrangian methods for small time-step sizes, due to the2nd order accuracy in space. For high accuracy regimes, with small time-step sizes,the ETD2RK scheme provided the best results. y ( k m ) Vorticity SL-ETD2RK t= 10.0 days dt= 0112.5 y ( k m ) Vor t icit y SL-ETD2 RK t = 1 0 .0 d ay s d t = 0 1 1 2 .5 y ( k m ) Vorticity ETD2RK t=10.0 days dt=0112.5 -4e-05-2e-050e+002e-054e-05 / s y ( k m ) Vorticity ETD2RK t=10.0 days dt=0112.5 -4e-05-2e-050e+002e-054e-05 / s (a) (b) Fig. 8 . Vorticity field of the full nonlinear SWE for the unstable jet test case at day 10 using atime-step size of 112.5 seconds. Images refer to (a) SL-ETD2RK and (b) ETD2RK schemes. Theplot for SL-SI-SETTLS is very similar to SL-ETD2RK and therefore not shown for sake of brevity.
Now, we study results for different time integration methods at day 10 of thebenchmark. The vorticity fields for two different schemes (SL-ETD2RK, ETD2RK)are depicted in Figure 8, where, since SL-SI-SETTLS provides results similar to SL-ETD2RK, this scheme was not depicted. They are again qualitatively very similar,although the ETD2RK shows more high wavenumber oscillations around the vortices. y ( k m ) Vorticity SL-SI-SETTLS t= 10.0 days dt= 0225.0 y ( k m ) Vor t icit y SL-SI-SETTLS t = 1 0 .0 d ay s d t = 0 2 2 5 .0 y ( k m ) Vorticity SL-ETD2RK t=10.0 days dt=0225.0 -4e-05-2e-050e+002e-054e-05 / s y ( k m ) Vorticity SL-ETD2RK t=10.0 days dt=0225.0 -4e-05-2e-050e+002e-054e-05 / s (a) (b) Fig. 9 . Numerical solution of the full nonlinear SWE at time 10 days for the vorticity fieldusing a time-step size of 225 seconds. (a) SL-SI-SETTLS, (b) SL-ETD2RK.
For larger time-step sizes, due to the extra energy in the high wavenumber gravitywaves, the SL-ETD2RK triggers small turbulent like features after long runs whencompared to SL-SI-SETTLS which is illustrated in Figure 9.Finally, this motivates to investigate the energy spectrum for selected time inte-gration methods and parameters, see Figure 10. If there is no dissipation of near gridscale energy for the SL-ETD2RK scheme, this energy destabilizes the jet into smallerscale features. This can be clearly observed in the energy spectrum, where we also
EMI-LAGRANGIAN EXPONENTIAL INTEGRATION SL-ETD2RKdt225ETD2RKSL-SI-SETTLSdt225
Fig. 10 . Kinetic energy spectrum for different methods and time-step sizes for the full nonlinearSWE at day 10 of integration. We can observe that the scheme SL-ETD2RK with 225 secondstime-step has developed less energy in large wavelengths and more energy in small wavelengths. TheETD2RK scheme has more energy in the small wavelengths and SL-SI-SETTLS with time-step of225 seconds has the least amount of energy in small scales. notice that the ETD2RK scheme gathers more energy in the smaller scales and theSL-SI-SETTLS with large time-step has the least amount of energy in the smallerscales.
For the pur-pose of weather and climate simulations, a certain amount of small-scale dissipationis usually required, either from a numerical stability perspective or from a physicalpoint of view. The SL-SI-SETTLS scheme, when used in the full IFS dynamical core,adopts a spectral hyper-diffusion filter in the momentum equations in order to bothnumerically stabilize the scheme and physically dissipate energy from the small-scaleenergy tail (see [26] for an analysis of the impacts of the diffusion in a global spectralmodel and [33] for a comprehensive discussion on the use of diffusion in atmosphericmodels). We remark that in full models this energy in high wavenumbers could beused to model physical sub-grid properties, such as convection.With the semi-Lagrangian exponential scheme, it is possible to preserve the pre-cise dispersion of linear waves and apply a term specific dissipation in the nonlineardivergence term. This way, linear waves (long and short) are treated accurately, butonly the longer waves originated from the nonlinear interactions are preserved in themodel. This allows the model to be numerically stable without damping the linearwaves, and also provides dissipation of small-scale features generated by the additionalenergy in high wavenumbers excited by the exponential integration.In the analysis that follows we considered an implicit spectral diffusion filter ( µ ∇ )applied only to the nonlinear divergence term with µ the diffusion coefficient.We start by analyzing different diffusion coefficients with the kinetic energy spec-trum of the SL-ETD2RK scheme at day 10. Figure 11a shows the amount of diffusionrequired to obtain a solution along the lines of the SL-SI-SETTLS with a time-stepsize of 900 seconds, and, following these results, we will adopt µ = 25 . × m s − .This value is similar to what is actually used in weather forecasting systems for the2 P. S. PEIXOTO, M. SCHREIBER Horizontal wavelength (km)10 K i n e t i c E n e ( g y Sp e c t ( u m ( m s ) k k Kinetic Ene(gy Spect(um t=10.0 days dt=0900.0s
SL-SI-SETTLSSL-ETD2RK_v12.8SL-ETD2RK_v25.6SL-ETD2RK_v51.2 Horizontal wavelength (km)10 K i n e t i c E n e ( g y Sp e c t ( u m ( m s ) k k Kinetic Ene(gy Spect(um t=10.0 days
SL-SI-SETTLS_dt0900.0SL-ETD2RK_dt0900.0ETD2RK_dt0225.0SL-EXP-SETTLS_dt0225.0 (a) (b)
Fig. 11 . (a) Kinetic energy spectrum considering an implicit diffusion filter on the nonlineardivergence term various values of µ ( . , . and . × m s − ), for the SL-ETD2RK schemeand no diffusion for the SL-SI-SETTLS scheme. (b) Kinetic energy spectrum considering an im-plicit diffusion on the nonlinear divergence term with µ = 25 . × m s − , for the schemes andparameters shown in Figure 13. (a) Timestep size (sec)10 (cid:1) M a x E rr o r SL-SI-SETTLSSL-ETD2RKSL-EXP-SETTLSETD2RK (b) Timestep size (sec)10 (cid:1) (cid:1) R M SE rr o r SL-SI-SETTLSSL-ETD2RKSL-EXP-SETTLSETD2RK
Fig. 12 . Errors at day 1 of integration of the unstable jet test case for the full nonlinear SWEwith implicit diffusion on the nonlinear divergence term ( µ = 25 . × m s − ). full equations, whereas here, we are only considering it for the nonlinear divergence.Next, we investigate error comparisons at day 1 of integration shown in Figure 12.We can observe that the two semi-Lagrangian exponential schemes deliver more ac-curacy compared to the SL-SI-SETTLS scheme. However, the following results willreveal significant differences in the two semi-Lagrangian exponential schemes for largertime integration ranges.Finally, we time integrate to day 10 and investigate the different time integrationmethods, with the vorticity field depicted in Figure 13. Despite the implicit diffusion,the ETD2RK scheme is still not able to do time-step sizes as large as the semi-Lagrangian schemes, due to the instability originated from the nonlinear advectionterm. In Figure 11 we show the kinetic energy spectra of the results shown in 13,where one may note the effects of the filtering in the small scale features. As before,the SL-EXP-SETTLS scheme is unstable for large time-steps and the diffusion was notable to circumvent this instability. However, the SL-ETD2RK scheme with diffusionnow does not develop near grid scale features (see Fig. 8b for comparison) even with EMI-LAGRANGIAN EXPONENTIAL INTEGRATION y ( k m ) Vorticity SL-SI-SETTLS t=10.0 days dt=0900.0 -4e-05-2e-050e+002e-054e-05 / s y ( k m ) Vorticity SL-SI-SETTLS t=10.0 days dt=0900.0 -4e-05-2e-050e+002e-054e-05 / s y ( k m ) Vorticity SL-ETD2RK t= 10.0 days dt= 0900.0 y ( k m ) Vor t icit y SL-ETD2 RK t = 1 0 .0 d ay s d t = 0 9 0 0 .0 (a) (b) y ( k m ) Vorticity ETD2RK t=10.0 days dt=0225.0 -4e-05-2e-050e+002e-054e-05 / s y ( k m ) Vorticity ETD2RK t=10.0 days dt=0225.0 -4e-05-2e-050e+002e-054e-05 / s y ( k m ) Vorticity SL-EXP-SETTLS t= 10.0 days dt= 0225.0 y ( k m ) Vor t icit y SL-EXP-SETTLS t = 1 0 .0 d ay s d t = 0 2 2 5 .0 (c) (d)
Fig. 13 . Numerical solution of the full SWE at time 10 days for the vorticity field using implicitdiffusion on the nonlinear divergence term with µ = 25 . × m s − . (a) SL-SI-SETTLS with ∆ t = 900 s , (b) SL-ETD2RK with ∆ t = 900 s , (c) ETD2RK with ∆ t = 225 s , (d) SL-EXP-SETTLSwith ∆ t = 225 s . a time-step size of 900 seconds.
7. Concluding remarks.
Semi-Lagrangian schemes and exponential integra-tors both play important roles in different applications. The exploration of a mixedformulation is a challenging problem, partially tackled in this paper. We show anovel approach that combines these methods by exploring the exponential integrationformulation in terms of material derivatives.The approach may be helpful for users of standard exponential integration tech-niques, extending these to larger time-step sizes preserving accurate solutions. Also,from the application of weather and climate modeling perspective, the method showsa way to improve the dispersion properties of a well-established scheme, thereforebetter representing fast linear gravity waves, with competitive computational cost.The results presented in this paper show the potential benefits of the combina-tion of these different methods, semi-Lagrangian and exponential integration, in aspecific scenario, of rotating shallow water equations. Nevertheless, it also points outmeans of using such an approach in other equation sets that present both nonlinearadvection and stiff linear operators, such as the Euler equations in fluid dynamics.We highlight, however, that the development of exponential schemes for fluids thatdevelop shocks, for which the use of conserved variables are usually preferred, stillrepresent a challenge.
Acknowledgements.
The ideas behind this work originated from discussionswith Colin Cotter, Jemma Shipton and Beth Wingate, whom are greatly acknowl-edged. We would also like to thank Saulo Barros for discussions with respect tosemi-Lagrangian spectral schemes.
Appendix A. Lack of commutation between linear operator and in-terpolation at trajectory points.
Consider a general vector (cid:126)w ∈ R n , a linearoperator T ∈ R n × R n , which will represent here, for example, a matrix exponential,4 P. S. PEIXOTO, M. SCHREIBER and I (cid:126)x : R n → R n an interpolation operation with respect to points (cid:126)x ∈ R n . Fol-lowing the semi-Lagrangian notation for interpolation, we may concisely write that I (cid:126)x ( (cid:126)w ) = (cid:126)w ∗ , where the ∗ implicitly indicates the interpolation with respect to (cid:126)x . Thissection is just to point a simple example to illustrate that even in very simple cases( T (cid:126)w ) ∗ (cid:54) = T ( (cid:126)w ∗ ).Consider a 1D periodic grid with uniformly spaced points ( x i ) i =1 ,n , with distance∆ x from each other. In this example we will consider a scalar advection with constantvelocity given by ∆ x/ ∆ t , so that, after one time-step, the departure points will be asimple translation and will match exactly their left neighbours. That is, the trajectorygoes from t n to t n +1 carrying the function value at x i − to the x i point. In this case,the interpolation to departure points will be given by a periodic shift in the indexes,(45) I (cid:126)x ( (cid:126)w ) = I (cid:126)x ([ w , w , w , . . . , w n ]) = [ w n , w , w , . . . , w n − ] = (cid:126)w ∗ . Note that the operator I (cid:126)x is a linear operator.Now consider a simple diagonal linear operator T = ( α ii ) i =1 ,n , with α ii (cid:54) = α jj ,for j (cid:54) = i . In this case,( T (cid:126)w ) ∗ = ([ α w , α w , w , . . . , α nn w n ]) ∗ (46) = [ α nn w n , α w , α w , w , . . . , α ( n − n − w n − ] , but(47) T ( (cid:126)w ∗ ) = T [ w n , w , w , . . . , w n − ] = [ α w n , α w , . . . , α nn w n − ] . Therefore, even if the trajectories are constant (or linear), the commutation does notgenerally hold.In the more general case treated in the derivation of the semi-Lagrangian exponen-tial scheme, the trajectories are nonlinear. Also, the linear operator is not necessarilydiagonal, but one could think of its diagonalized version in complex space in a similarway, for which the terms in the diagonal would be the eigenvalues of the operator.
REFERENCES[1]
J. D. Anderson and J. Wendt , Computational fluid dynamics , vol. 206, Springer, 1995.[2]
A. Arakawa and V. Lamb , Computational design of the basic dynamical processes of theUCLA general circulation model , Methods in Computational Physics, 17 (1977), pp. 173–265.[3]
R. K. Archibald, K. J. Evans, J. Drake, and J. White , Time acceleration methods foradvection on the cubed sphere , in International Conference on Computational Science,Springer, 2009, pp. 253–262.[4]
S. R. M. Barros, D. Dent, L. Isaksen, G. Robinson, G. Mozdzynski, and F. Wollenwe-ber , The IFS model: A parallel production weather code , Parallel Computing, 21 (1995),pp. 1621–1638.[5]
P. Bartello and S. J. Thomas , The cost-effectiveness of semi-Lagrangian advection , Monthlyweather review, 124 (1996), pp. 2883–2897.[6]
J. Bates, F. Semazzi, R. Higgins, and S. R. Barros , Integration of the shallow water equa-tions on the sphere using a vector semi-lagrangian scheme with a multigrid solver , MonthlyWeather Review, 118 (1990), pp. 1615–1627.[7]
A. Bermudez, M. R. Nogueiras, and C. Vazquez , Numerical analysis of convection-diffusion-reaction problems with higher order characteristics/finite elements. part i: timediscretization , SIAM Journal on Numerical Analysis, 44 (2006), pp. 1829–1853.[8]
A. Berm´udez, C. Rodr´ıguez, and M. Vilar , Solving shallow water equations by a mixedimplicit finite element method , IMA journal of numerical analysis, 11 (1991), pp. 79–97.[9]
G. Beylkin, J. M. Keiser, and L. Vozovoi , A new class of time discretization schemes for thesolution of nonlinear PDEs , Journal of Computational Physics, 147 (1998), pp. 362–387.EMI-LAGRANGIAN EXPONENTIAL INTEGRATION [10] E. Carr, I. Turner, and P. Perr , A variable-stepsize jacobian-free exponential integrator forsimulating transport in heterogeneous porous media: Application to wood drying , Journalof Computational Physics, 233 (2013), pp. 66 – 82.[11]
E. Celledoni , Eulerian and semi-Lagrangian schemes based on commutator-free exponentialintegrators , Group theory and numerical analysis, 39 (2005), pp. 77–90.[12]
E. Celledoni and B. K. Kometa , Semi-Lagrangian Runge-Kutta exponential integrators forconvection dominated problems , Journal of Scientific Computing, 41 (2009), pp. 139–164.[13]
E. Celledoni, B. K. Kometa, and O. Verdier , High order semi-Lagrangian methods forthe incompressible Navier–Stokes equations , Journal of Scientific Computing, 66 (2016),pp. 91–115.[14]
C. Clancy and P. Lynch , Laplace transform integration of the shallow-water equations. PartII: Semi-Lagrangian formulation and orographic resonance , Quarterly Journal of the RoyalMeteorological Society, 137 (2011), pp. 800–809.[15]
C. Clancy and J. A. Pudykiewicz , On the use of exponential time integration methodsin atmospheric models , Tellus A: Dynamic Meteorology and Oceanography, 65 (2013),p. 20898.[16]
S. M. Cox and P. C. Matthews , Exponential time differencing for stiff systems , Journal ofComputational Physics, 176 (2002), pp. 430–455.[17]
M. Diamantakis , The semi-Lagrangian technique in atmospheric modelling: current statusand future challenges , in ECMWF Seminar in numerical methods for atmosphere andocean modelling, 2013, pp. 183–200.[18]
T. Dubos, S. Dubey, M. Tort, R. Mittal, Y. Meurdesoif, and F. Hourdin , DYNAMICO-1.0, an icosahedral hydrostatic dynamical core designed for consistency and versatility ,Geoscientific Model Development, 8 (2015), pp. 3131–3150.[19]
D. R. Durran , Numerical methods for fluid dynamics: With applications to geophysics , vol. 32,Springer, 2010.[20]
ECMWF , PART III: DYNAMICS AND NUMERICAL PROCEDURES , IFS Documentation,ECMWF, 2017, ch. ., p. .[21]
M. Falcone and R. Ferretti , Convergence analysis for a class of high-order semi-lagrangianadvection schemes , SIAM Journal on Numerical Analysis, 35 (1998), pp. 909–940.[22]
S. N. Figueroa, J. P. Bonatti, P. Y. Kubota, G. A. Grell, H. Morrison, S. R. Barros,J. P. Fernandez, E. Ramirez, L. Siqueira, G. Luzia, et al. , The Brazilian globalatmospheric model (BAM): performance for tropical rainfall forecasting and sensitivity toconvective scheme and horizontal resolution , Weather and Forecasting, 31 (2016), pp. 1547–1572.[23]
J. Galewsky, R. K. Scott, and L. M. Polvani , An initial-value problem for testing numericalmodels of the global shallow-water equations , Tellus A, 56 (2004), pp. 429–440.[24]
F. Garcia, L. Bonaventura, M. Net, and J. S´anchez , Exponential versus IMEX high-ordertime integrators for thermal convection in rotating spherical shells , Journal of Computa-tional Physics, 264 (2014), pp. 41–54.[25]
S. Gaudreault and J. A. Pudykiewicz , An efficient exponential time integration method forthe numerical solution of the shallow water equations on the sphere , Journal of Computa-tional Physics, 322 (2016), pp. 827 – 848.[26]
A. Gelb and J. P. Gleeson , Spectral viscosity for shallow water equations in spherical geom-etry , Monthly Weather Review, 129 (2001), pp. 2346–2360.[27]
T. Haut, T. Babb, P. Martinsson, and B. Wingate , A high-order time-parallel scheme forsolving wave propagation problems via the direct construction of an approximate time-evolution operator , IMA Journal of Numerical Analysis, (2015), p. drv021.[28]
M. Hochbruck and C. Lubich , On Krylov subspace approximations to the matrix exponentialoperator , SIAM Journal on Numerical Analysis, 34 (1997), pp. 1911–1925.[29]
M. Hochbruck and A. Ostermann , Exponential integrators , Acta Numerica, 19 (2010),pp. 209–286.[30]
M. Hortal , The development and testing of a new two-time-level semi-lagrangian scheme(settls) in the ecmwf forecast model , Quarterly Journal of the Royal Meteorological Society,128 (2002), pp. 1671–1687.[31]
G. L. Kooij, M. A. Botchev, and B. J. Geurts , An Exponential Time Integrator for theIncompressible Navier–Stokes Equation , SIAM Journal on Scientific Computing, 40 (2018),pp. B684–B705.[32]
J. N. Koshyk and K. Hamilton , The horizontal kinetic energy spectrum and spectral budgetsimulated by a high-resolution troposphere–stratosphere–mesosphere GCM , Journal of theAtmospheric Sciences, 58 (2001), pp. 329–348.[33]
P. H. Lauritzen, C. Jablonowski, M. A. Taylor, and R. D. Nair , Numerical techniques P. S. PEIXOTO, M. SCHREIBER for global atmospheric models , vol. 80, Springer Science & Business Media, 2011.[34]
E. Lindborg , Can the atmospheric kinetic energy spectrum be explained by two-dimensionalturbulence? , Journal of Fluid Mechanics, 388 (1999), pp. 259–288.[35]
J. Loffeld and M. Tokman , Comparative performance of exponential, implicit, and explicitintegrators for stiff systems of ODEs , Journal of Computational and Applied Mathematics,241 (2013), pp. 45 – 67.[36]
A. Majda , Introduction to PDEs and Waves for the Atmosphere and Ocean , vol. 9, AmericanMathematical Soc., 2003.[37]
A. McDonald , Accuracy of Multiply-Upstream Semi-Lagrangian Advective Schemes II , Mon.Wea. Rev., 115 (1987), pp. 1446–1450.[38]
C. McLandress , On the importance of gravity waves in the middle atmosphere and their pa-rameterization in general circulation models , Journal of Atmospheric and Solar-TerrestrialPhysics, 60 (1998), pp. 1357–1383.[39]
G. Mengaldo, A. Wyszogrodzki, M. Diamantakis, S.-J. Lock, F. X. Giraldo, and N. P.Wedi , Current and Emerging Time-Integration Strategies in Global Numerical Weatherand Climate Prediction , Archives of Computational Methods in Engineering, (2018), pp. 1–22.[40]
P. S. Peixoto and S. R. Barros , On vector field reconstructions for semi-lagrangian transportmethods on geodesic staggered grids , J. Comput. Phys., 273 (2014), pp. 185 – 211.[41]
O. Pironneau , On the transport-diffusion algorithm and its applications to the Navier-Stokesequations , Numerische Mathematik, 38 (1982), pp. 309–332.[42]
F. Plunian, R. Stepanov, and P. Frick , Shell models of magnetohydrodynamic turbulence ,Physics Reports, 523 (2013), pp. 1–60.[43]
D. A. Randall , Geostrophic Adjustment and the Finite-Difference Shallow-Water Equations ,Mon. Wea. Rev., 122 (1994), pp. 1371–+.[44]
A. Robert , A stable numerical integration scheme for the primitive meteorological equations ,Atmosphere-Ocean, 19 (1981), pp. 35–46.[45]
A. Robert , A semi-Lagrangian and semi-implicit numerical integration scheme for the prim-itive meteorological equations , Journal of the Meteorological Society of Japan. Ser. II, 60(1982), pp. 319–325.[46]
R. Sadourny , The dynamics of finite-difference models of the shallow-water equations , Journalof the Atmospheric Sciences, 32 (1975), pp. 680–689.[47]
M. Schreiber and R. Loft , A parallel time-integrator for solving the linearized shallow waterequations on the rotating sphere , under revision in Numer. Linear Algebra Appl., (2018).[48]
M. Schreiber, P. S. Peixoto, T. Haut, and B. Wingate , Beyond spatial scalability limita-tions with a massively parallel method for linear oscillatory problems , The InternationalJournal of High Performance Computing Applications, (2017), p. 1094342016687625.[49]
M. Schreiber, N. Schaeffer, and R. Loft , Exponential Integrators with Parallel-in-TimeRational Approximations for Shallow-Water Equations on the Rotating Sphere , ParallelComputing, (2019).[50]
J. C. Schulze, P. J. Schmid, and J. L. Sesterhenn , Exponential time integration using Krylovsubspaces , International journal for numerical methods in fluids, 60 (2009), pp. 591–609.[51]
W. C. Skamarock, J. B. Klemp, M. G. Duda, L. D. Fowler, S.-H. Park, and T. D.Ringler , A Multiscale Nonhydrostatic Atmospheric Model Using Centroidal Voronoi Tes-selations and C-Grid Staggering , Mon. Wea. Rev., 140 (2012), pp. 3090–3105.[52]
A. St-Cyr and S. J. Thomas , Nonlinear operator integration factor splitting for the shallowwater equations , Applied Numerical Mathematics, 52 (2005), pp. 429–448.[53]
A. Staniforth and J. Ct , Semi-Lagrangian Integration Schemes for Atmospheric Models - AReview , Mon. Wea. Rev., 119 (1991), pp. 2206–2223.[54]
M. Tokman , Efficient integration of large stiff systems of ODEs with exponential propagationiterative (EPI) methods , Journal of Computational Physics, 213 (2006), pp. 748–776.[55]
D. L. Williamson , The evolution of dynamical cores for global atmospheric models , J. Mete-orol. Soc. Jpn., 85B (2007), pp. 241–269.[56]
D. L. Williamson, J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber , A standardtest set for numerical approximations to the shallow water equations in spherical geometry ,J. Comput. Phys., 102 (1992), pp. 211–224.[57]