- Email: [email protected]
Numerical solution of multiscale problems in atmospheric modeling
Applied Numerical Mathematics 62 (2012) 1531–1543
Contents lists available at SciVerse ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
Numerical solution of multiscale problems in atmospheric modeling Martin Schlegel a,∗ , Oswald Knoth a , Martin Arnold b , Ralf Wolke a a b
Institute for Tropospheric Research, Permoserstrasse 15, 04318 Leipzig, Germany Martin Luther University Halle-Wittenberg, Institute for Mathematics, 06099 Halle (Saale), Germany
a r t i c l e
i n f o
Article history: Available online 7 June 2012 Keywords: Time integration Runge–Kutta method Multirate scheme Flux splitting IMEX
a b s t r a c t Explicit time integration methods are characterized by a small numerical effort per time step. In the application to multiscale problems in atmospheric modeling, this benefit is often more than compensated by stability problems and stepsize restrictions resulting from stiff chemical reaction terms and from a locally varying Courant–Friedrichs–Lewy (CFL) condition for the advection terms. In the present paper, we address this problem by a rather general splitting technique that may be applied recursively. This technique allows the combination of implicit and explicit methods (IMEX splitting) as well as the local adaptation of the time stepsize to the meshwidth of non-uniform space grids in an explicit multirate discretization of the advection terms. Using a formal representation as partitioned Runge–Kutta method, convergence of order p 3 may be shown if some additional order conditions are satisfied. In a series of numerical tests, the convergence results are verified and the consequences of different splitting strategies like flux splitting and cell splitting are analysed. © 2012 Published by Elsevier B.V. on behalf of IMACS.
1. Introduction State of the art numerical models can be used to simulate a broad spectrum of atmospheric phenomena. These range from microphysical and chemical processes with characteristic times below a second and characteristic scales in the submillimeter range up to high/low pressure areas and planetary waves on scales about 1000 kilometers and characteristic times up to months. Additionally the simulation domain is often locally refined to allow for a more accurate examination of regions of interests [1,3]. These may be densely populated areas in contrast to surrounding woodland or agricultural regions, or in a global model continents in contrast to oceans. Atmospheric phenomena can be described using the advection–diffusion–reaction equation
∂ ∂ ∂ ∂ c + f (x, t , c ), ρD c = − (uc ) + ∂t ∂x ∂x ∂x ρ
(1)
with the concentration of a transported quantity c, wind speed u, density ρ , diffusion coefficient D and a reaction term f . The latter may also include sources and sinks. The advection part of this equation can be solved efficiently using explicit Runge–Kutta (ERK) time integration methods. These methods have in common that stability requirements limit the global time step to be smaller than some critical value proportional to the ratio of grid size and the magnitude of the wind speed for each cell [2,7,17]. Consequently the smallest cell, or more exactly the cell with the smallest characteristic time determines the global time step.
*
Corresponding author. E-mail addresses: [email protected] (M. Schlegel), [email protected] (O. Knoth), [email protected] (M. Arnold), [email protected] (R. Wolke). 0168-9274/$36.00 © 2012 Published by Elsevier B.V. on behalf of IMACS. http://dx.doi.org/10.1016/j.apnum.2012.06.023
1532
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
Multirate time integration methods can be employed to adapt the time step locally so that slower components take longer and fewer time steps, resulting in a moderate to substantial reduction of the computational cost, depending on the scenario to simulate [4]. Multirate approaches have been developed since the early 1980s. Osher and Sanders [11] presented a scheme allowing multirate Euler steps. More current approaches allow for a wide variety of multirate schemes by proposing generic multirate methods based on classical ERKs. In 2005 Tang and Warnecke [14] proposed different internally consistent multirate schemes. These schemes however are not mass preserving. The approach of Constantinescu and Sandu [2] yields multirate ERK schemes which are mass preserving and at most of second order accuracy. In 2007 Hundsdorfer et al. showed that multirate schemes which in absence of a slow domain reduce to multiple applications of the underlying base method cannot be mass preserving and consistent at the same time [6]. The diffusion and reaction terms in (1) on the other hand are usually stiff and consequently cannot be solved efficiently using explicit methods. Instead implicit methods must be employed. Splitting methods solving the non-stiff terms using explicit methods and the stiff terms using implicit methods are known as implicit/explicit (IMEX) methods. This class of splitting methods has proven efficient for problems which can be split into non-stiff (or mildly stiff) and stiff terms. We presented an internally consistent (in the nomenclature of [6]) generic multirate scheme for constructing third order accurate multirate methods [13]. Since for atmospheric modeling mass preservation or more generally the preservation of linear invariants is of major interest we employ a mass preserving splitting of the advection operator, thus trading consistency (i.e. truncation error O (t ) + O (x)) for mass preservation. In Section 2 we shall exemplify the construction of multirate schemes using our algorithm. The possible extension to an implicit/explicit (IMEX) multirate scheme will also be considered. Section 3 is dedicated to the formal analysis of the constructed methods. Numerical tests will be presented in Section 4. 2. A new class of hierarchical multirate IMEX schemes The multirate scheme we shall introduce is based on an implicit/explicit (IMEX) method introduced by Knoth and Wolke in [10] for the efficient solution of advection–diffusion–reaction equations in air pollution. Consider a splitting of the righthand side
w = F ( w ) + G ( w ),
w (0) = w 0 ,
w ∈ RN .
(2)
In the context of this splitting G is a non-stiff term which will be solved explicitly whereas F is a stiff term which must be solved using an implicit method. For the theoretical analysis of the method the authors assume that the implicit integration is carried out exactly. Defining the nodes of the explicit method
τi = tn + tc i , the full algorithm computing an approximate solution w n+1 for time tn+1 from an approximate solution w n at time tn is defined as follows
W 1 = wn, ri =
i −1
(3)
(ai j − ai −1, j )G ( W j ),
(4)
j =1
v i (τi −1 ) = W i −1 , ∂ vi 1 = r i + F ( v i ), ∂τ c i − c i −1
(5)
τ ∈ [τi−1 , τi ], i = 2, . . . , s + 1,
(6)
W i = v i (τi ),
(7)
w n +1 = W s +1 ,
(8)
with s denoting the number of stages of an explicit Runge–Kutta method with parameters (a, b, c ) in common notation. For simplicity we define as+1, j = b j , thus making a separate treatment of the summation stage unnecessary. We note that the cumulative implicit integration time is t. There is another class of comparable IMEX schemes where the integration for the stiff part starts in each internal stage at the old time stage tn and the initial value from the previous integration step and therefore the cumulative interval is usually larger than the explicit time step [16]. If the stiff term is identical to zero F ≡ 0 we obtain the underlying explicit method applied to w = G ( w ). To obtain a multirate scheme we again consider a splitting of the right-hand side,
w = G 1 ( w ) + G 2 ( w ),
w (0) = w 0 ,
w ∈ RN ,
(9)
but now both terms (now denoted G 1 and G 2 ) are assumed to be non-stiff with G 2 imposing a more restrictive time step limitation. Referring to the characteristic times of these methods we will also call G 1 a slow term opposed to the fast term G 2 . Now two potentially but not necessarily different explicit Runge–Kutta methods are employed. The method replacing
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
1533
the implicit integrator is called the inner method, denoted by a superscript I. The other method is called the outer method, denoted by a superscript O. The multirate algorithm reads
W 1 = wn, ri =
i −1
(10)
j =1
O aO i j − ai −1, j G 1 ( W j ),
(11)
V i ,1 = W i −1 ,
(12)
V i ,k = V i ,k−1 + t c iO − c iO−1
k −1 j =1
I akj − akI −1, j
1 c iO − c iO−1
ri + G 2 ( V i, j )
(13)
i −1 k −1 O O I O = V i ,k−1 + t ckI − ckI −1 a i j − aO akj − akI −1, j G 2 ( V i , j ), i −1, j G 1 ( W j ) + t c i − c i −1 j =1
(14)
j =1
W i = V i , s I +1 ,
(15)
w n +1 = W s O +1 = V s O +1 , s I +1 .
(16)
We call this algorithm recursive flux splitting multirate (RFSMR), [13]. The name emphasizing the splitting by fluxes comes from the original application to the advection equation in atmospheric simulations where mass preservation is essential. The algorithm (10)–(16) can also be cast as a partitioned Runge–Kutta (PRK) system:
W 1 = wn, W i = w n + t
i −1
aslow G 1 ( W j ) + afast ij i j G 2 (W j ) ,
j =1
w n+1 = w n + t
s
bslow G 1 ( W j ) + bfast j j G 2 (W j ) ,
j =1
which is the general form of an explicit PRK applied on a splitting of the right-hand side. We will make use of this notation in the following chapter. The Runge–Kutta parameters (aslow , bslow , c slow ) and (afast , bfast , c fast ) employed here depend on the parameters of the base methods as shown below. For easier association of the parameters of the partitioned method with the parameters of the underlying base methods we employ tuple indexes. The first element of a tuple is associated with a stage of the outer method, the other element corresponds to an inner method stage:
( i , k ) ↔ ( i − 1) s I + k . The following equations can be deduced by transforming algorithm (10)–(16) from incremental to common form and solving the resulting recurrence relations. The complete set of multirate PRK coefficients reads:
s = sO · sI ,
= aslow aslow (i ,k),( j ,l) = (i −1)sI +k,( j −1)sI +l
bslow ( j ,l) =
bOj
if l = 1,
0
if l > 1,
(17) aO i, j
+ (aOi +1, j
0
− aOi , j )ckI
if l = 1, if l > 1,
(18) (19)
O I O O c (slow i ,k) = c i + ck c i +1 − c i , ⎧ ⎪0
(20)
afast (i ,k),( j ,l)
(21)
=
⎨
if j > i , (c Oj +1 − c Oj )akI ,l if j = i ,
⎪ ⎩ (c O
− c Oj )blI if j < i , j +1 O O I bfast ( j ,l) = c j +1 − c j bl , slow O I O O c (fast i ,k) = c (i ,k) = c i + ck c i +1 − c i .
(22) (23)
Note that different pairs of base methods may result in equivalent multirate methods. A possible disadvantage of the RFSMR algorithm is that there is no parameter defining the time step ratio, i.e. the relative stability and corresponding to that the relative computational cost of the slow and the fast method. Instead the slow method is equivalent to the outer base method while the fast method is equivalent to the inner method executed sO times with
1534
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
Table 1 RFMSMR(RK2a) – an example for a second order multirate scheme; redundant stages omitted. 0 1
0 1/2 1/2 1
1 1/2 1/2 (RK2a)
outer base method 0 1/2 1/2 1 1
1/2 1/2 1 1 1/2
0 0 0 0 0 0 0 “slow” part
0 0
1/2 1/4 1/4 1/4 1/4 1/2 1/4 1/4 1/4 inner base method 0 1/2 1/2 1 1
1/2
1/2 1/4 1/4 1/4 1/4
1/4
1/4 1/4 1/2 1/4 1/4 1/4 1/4 “fast” part
1/4 1/4
0
Fig. 1. Schematic of recursive construction.
time step sizes of (c O − c Oi )t; i ∈ {2, . . . , sO + 1}. This disadvantage may however be easily overcome by substituting the i +1 inner method by a composition of n steps of step size t /n. For instance in Table 1, the inner method is equivalent to the outer method executed twice with a half time step, in combination of the outer method’s nodes c = (0, 1, 1) leading to a time step ratio of R = 2 for the partitioned multirate method. For practical purposes two temporal refinement levels are often insufficient. For instance many temporal refinement levels are desirable when nested spatial refinement is employed, as is common in atmospheric modeling. To obtain a multirate scheme with an arbitrary number of temporal refinement levels the recursive flux-splitting multirate approach may be applied recursively. One of the methods created in one construction stage then serves as outer method in the following construction stage. Note that the PRKs created using the RFSMR algorithm are internally consistent (i.e. c slow = c fast ) and the parameters of the fast method do not depend directly on aO or bO , comp. (21)–(23). Consequently it makes no difference if the slow or the fast method created in the previous construction stage is employed as outer method in the next stage. Like this an arbitrary number of temporal refinement levels can be realized, see Fig. 1. Arrows indicate use of variables. The ellipses marked “slow construction” denote construction of the slow method’s parameters as defined in (18)–(20); analogously the “fast construction” ellipses are correlated to (21)–(23). Complex multiphysics simulations often involve both terms which may be approximated efficiently using explicit methods and stiff terms requiring an implicit integrator. Based on an IMEX splitting, our explicit multirate algorithm may easily be extended to a multirate IMEX algorithm combining the advantages of both splitting techniques. As before consider a splitting of the right-hand side
w = G 1 ( w ) + G 2 ( w ) + F ( w ), now with a stiff term F and two non-stiff terms G 1 , G 2 with different characteristic times, G 1 being the slower one. A multirate IMEX splitting may now be obtained by hierarchical application of the IMEX splitting (3)–(8). Again we define intermediate times, now correlated to the stages of the smaller explicit time step
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
1535
Fig. 2. Illustration of intermediate points in time.
τi,k = tn + t c iO−1 + c iO − c iO−1 ckI = tn + tc ifast ,k , see also Fig. 2. First we apply the splitting on w = G 1 ( w ) + [G 2 ( w ) + F ( w )]:
W 1 = wn, ri =
i −1
(24)
j =1
O aO i j − ai −1, j G 1 ( W j ),
(25)
v i (τi ,1 ) = W i −1 ,
(26)
∂ vi 1 = O r i + G 2 ( v i ) + F ( v i ), ∂τ c i − c iO−1
τ ∈ [τi,1 , τi+1,1 ], i = 2, . . . , sO + 1,
(27)
W i = v i (τi +1,1 ),
(28)
w n +1 = W s O +1 .
(29)
The same splitting (3)–(8) is applied to (27). Note that the time dependent variable v employed in (27) now becomes a Runge–Kutta stage variable V , with V i ,1 ≈ v (τi ,1 ) and V i ,sI +1 ≈ v (τi +1,1 ).
V i ,1 = W i −1 ,
(30)
u i ,k (τi ,k−1 ) = V i ,k−1 , r i ,k =
(31)
k −1 j =1
I akj − akI −1, j
1 c iO − c iO−1
∂ u i ,k 1 = I r i ,k + F (u i ,k ), ∂t ck − ckI −1
ri + G 2 ( V i, j ) ,
(32)
τ ∈ [τi,k−1 , τi,k ], k = 2, . . . , sI + 1,
(33)
V i ,k = u i ,k (τi ,k ),
(34)
W i = V i , s I +1 .
(35)
3. Order conditions As shown in the preceding section the RFSMR algorithm may be cast in PRK form. This gives us the means to examine the order of the created multirate methods in the well established context of partitioned Runge–Kutta methods. Order conditions for PRKs are composed of the conditions for the separate methods:
Order 1:
s
b i = 1,
(36)
i =1
Order 2:
s
bi ci =
i =1
Order 3:
s
b i c i2 =
i =1 s s i =1 j =1
1 2
(37)
,
1 3
(38)
,
b i ai j c j =
1 6
,
(39)
1536
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
and additional coupling conditions (see e.g. [8]) setting the parameters of the different methods into relation:
Order 2:
s
s
bslow c ifast = i
i =1
Order 3:
s
slow bfast = i ci
i =1
bslow c ifast i
2
=
i =1
s
bfast c islow i
fast bslow afast = i ij c j
s s
i =1 j =1
(40)
,
2
1
= , 3
slow slow bfast cj = i ai j
i =1 j =1
bfast c ifast − c islow i
2
=
i =1 s
2
i =1
s s
s
1
s
bslow c ifast − c islow i
1 6
2
(41)
,
= 0,
i =1
s
fast fast bfast c i − c islow = i ci
i =1
bslow c islow c ifast − c islow = 0, i
i =1
s s
fast
bfast i ai j
slow
c fast −cj j
=
i =1 j =1 s s
s s
= 0, bslow aslow c fast − c slow i ij j j
i =1 j =1
slow fast = bfast c j − c slow i ai j j
i =1 j =1
s s
= 0. bslow afast c fast − c slow i ij j j
i =1 j =1
As the methods generated using RFSMR are internally consistent (i.e. c slow = c fast ) all coupling conditions apart from (41) are redundant. Note that if the method would additionally preserve all linear invariants (i.e. bslow = bfast ) the second order conditions (37), (40) would be satisfied automatically. We can prove order conditions for the generated partitioned methods assuming that the underlying base methods satisfy certain conditions. Theorem 3.1. If both the inner and the outer method satisfy the conditions for second order consistency (36), (37) the generated partitioned method is of second order accuracy, i.e. the fast and the slow method satisfy (36), (37). Corollary 3.2. A second order multirate scheme with an arbitrary number of temporal refinement levels can be constructed based on a set of explicit Runge–Kutta methods if all of the employed base methods satisfy the conditions for second order consistency (36), (37). Proof of Theorem 3.1. Proving this theorem is straightforward; apply definitions (17)–(23) on the order conditions and simplify the resulting expression. At some points we will employ the assumptions on the base methods. This will be denoted by under braces. Order 1 condition (36) for slow method: O I
s s
O
bslow = j
j =1
s
bOj = 1.
j =1
=1
Order 1 condition (36) for fast method: O I
s s j =1
O
bfast j
=
I
s s
c Oj +1
j =1 k =1
− c Oj
bkI
=
sO
j =1
c Oj +1
s +1
O I
s s j =1
O
bslow c slow = j j
s
bOj c Oj =
j =1
=1/2
1 2
.
=c OO
Order 2 condition (37) for slow method:
− c Oj =1
I
s
bkI .
k =1
=1
= 1.
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
1537
Order 2 condition (37) for fast method: O I
s s
O
fast bfast = j cj
j =1
I
s s j =1 k =1 O
=
I
s s O j =1 k =1
=
2
c j c Oj +1 − c Oj bkI + c Oj +1 − c Oj bkI ckI
sO O j =1
c Oj +1 − c Oj bkI · c Oj + c Oj +1 − c Oj ckI
I
c j c Oj +1
− c Oj
s
bkI +
k =1
c Oj +1
2 − c Oj
=
j =1
2
c Oj +1 c Oj − c Oj
I
s
bkI ckI
k =1
=1
sO
=1/2
+
2 1 O 2 c j +1 − 2c Oj +1 c Oj + c Oj 2
sO
1 O 2 O 2 1 O 2 1 = = c j +1 − c j c O +1 = . 2 2 s 2 j =1
2
=1
An analogous theorem can be formulated for the conditions for third order consistency (38), (39). Trying to solve for the coupling condition (41) for third order consistent partitioned methods however we find an additional condition to be satisfied by the outer method: s i −1 1 (c i +1 − c i ) (ai +1, j + ai j )c j = . i =1
j =1
3
(42)
Note that here this equation is applied to the parameters of the outer method
(a, b, c ) = aO , bO , c O . The more general notation is chosen because it will be applied to parameters of other methods as well. Condition (42) was already derived by Knoth and Wolke in [10] for the IMEX splitting introduced there. However the authors assumed that the implicit integrator (in the current context replaced by the inner explicit method) was exact. In the same paper also explicit Runge–Kutta methods satisfying this condition were presented. One of those methods and the multirate PRK resulting from it are listed in Table 2. Theorem 3.3. If both the inner and the outer method satisfy the conditions for third order consistency (36)–(39) and the outer method additionally satisfies the new order condition (42) the generated partitioned method is of third order accuracy, i.e. the fast and the slow method satisfy (36)–(41). If both of the base methods satisfy the conditions for third order consistency (36)–(39) and both methods satisfy the additional condition (42) both the fast and the slow method generated satisfy (42). The proof of this theorem follows the same idea as the proof of Theorem 3.1 but is very technical and would go beyond the scope of this paper. A last result following directly from Theorem 3.3 shall be given: Corollary 3.4. A third order multirate scheme with an arbitrary number of temporal refinement levels can be constructed based on a set of explicit Runge–Kutta methods if all of the employed base methods satisfy the conditions for third order consistency (36)–(39) and all but the last inner method satisfy (42). Consequently we are also able to construct a multirate IMEX scheme of third order accuracy employing an implicit Runge–Kutta method satisfying the classical conditions for third order consistency (36)–(39) as last inner method in a construction cascade as illustrated in Fig. 1. 4. Numerical tests In the context of partitioned methods it is of interest how exactly the splitting is defined. With a finite volume approach conservation laws
1538
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
Table 2 RFSMR(RK43) An example for a third order multirate scheme; redundant stages omitted. 0 1/2 1/2 1
1/2
−1/6
2/3 1/3 −1/3 1 1/6 1/3 1/3 1/6 base method (inner and outer) (RK43)
0 1/4 1/4 1/2 1/2 1/2 3/4 3/4 1 1
0 1/4 1/4 1/2 1/2 1/2 3/4 3/4 1 1
1/4 1/4 1/2 1/2 −1/6 1/12 1/12 1/3 1/3 1/6
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 2/3 0 1/6 0 1/6 0 −1/3 0 −1/3 0 1/3 “slow” part
1/2 1/2 1 1 1/3
0 0 0 0
0 0 0
0 0
1/6
1/4
−1/12
1/3
1/6 1/12 1/12 1/12 1/12 1/12 1/12
−1/6 1/6 1/6 1/6 1/6 1/6 1/6
1/2 1/6 1/6 1/6 1/6 1/6 1/6
1/12
1/6
1/6
1/12 1/12 1/12 1/12 1/12 1/12
0 0 0 0 0
1/4
−1/12 1/6 1/12
1/12 0 1/12 “fast” part
1/3
−1/6 1/6
1/2 1/6
1/12
1/6
1/6
1/12
0
Table 3 Multirate scheme according to Tang and Warnecke [14]. 0 1/2 1/2 1
1/2 1/4 1/4 1 0 0 1/2 0 0 slow method
0 1/2 1/2 1 1/2
1/2 1/4 1/4 1/4 fast
1/4 1/4 1/2 1/4 1/4 method
1/4
Table 4 Multirate scheme according to Constantinescu and Sandu [2]. 0 1 0 1
1 0 0 1/4
0 0 1 1/4 1/4 slow method
0 1/2 1/2 1 1/4
1/2 1/4 1/4 1/4 fast
1/4 1/4 1/2 1/4 1/4 method
1/4
∂w f (w) =− ∂t ∂x naturally lead to a so-called flux form
∂ wi 1 = F i w (t ) = − f i +1/2 w (t ) − f i −1/2 w (t ) ∂t xi with cell indexes i and the numerical fluxes f i ±1/2 at the boundaries of the cells. There are several possibilities to split the operator F = F 1 + F 2 . The most commonly used approach is the so-called splitting by cells alias splitting by components. This can be interpreted in such a way that F 1 describes the net fluxes of a subset of all cells and F 2 describes the net fluxes of the remaining cells. If not mentioned otherwise we tacitly assume a splitting by components. Alternatively a so-called splitting by fluxes can be employed. This approach can be interpreted so that F 1 describes the fluxes over a subset of all cell boundaries and F 2 describes the fluxes over the remaining cell boundaries. Opposed to the splitting by components this guarantees mass preservation independent from the specific time integration method as both F 1 and F 2 are mass preserving. However as already shown by Hundsdorfer et al. in [6] applying a splitting by fluxes will not lead to a consistent full discretization.
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
1539
Fig. 3. Decomposition of the unit interval in slow and fast blocks.
Though called recursive flux splitting multirate, the multirate schemes constructed via our algorithm may also be applied on a splitting by cells; thereby trading mass preservation for consistency. The following test setting is intended to examine the effects of these different splitting approaches when the resulting split scheme is solved using an explicit multirate method generated using RFSMR. It also is a mean to compare RFSMR to multirate methods from the literature, one given by Tang and Warnecke in [14] the other by Constantinescu and Sandu in [2]. (See Tables 3 and 4.) Both of these methods are based on RK2a (see Table 1) and converge with second order. Our test equation is one-dimensional advection with unit speed and a smooth initial profile
∂w ∂w =− , ∂t ∂x
w (0) = sin2 (2π x),
x ∈ [0, 1],
with periodic boundary conditions. The system is discretized in space employing an equidistant grid and the WENO5 approximation [9] for fluxes. In [15], Wang and Spiteri formally showed that this approximation is not linearly stable if combined with certain explicit Runge–Kutta methods, including RK2a. Numerical evidence for instability in [15], Section 4.1 was given for a Courant number of ν = 1.32. We did however not observe this behavior for the Courant numbers ν ∈ {0.2, 0.4} considered in the following tests. To examine the properties of the various multirate methods the grid is split into 5 blocks solved using the respective slow method and 5 blocks solved using the fast method, see Fig. 3. We consider the difference between the numerical solution at final time t = 1 and the exact solution. In [6], Hundsdorfer et al. presented a similar setting intended to examine the differences of consistent and mass preserving multirate schemes. As pointed out in the same paper, a true multirate scheme cannot be consistent and mass preserving at the same time. In the nomenclature of the authors a true multirate scheme is a scheme for which the respective slow method is equivalent to the base method in terms of computational cost and stability and the fast method is equivalent to the base method executed R times with time steps t / R for a given time step ratio R. First the second order schemes are examined. We choose t = 5 · 10−4 ; x = 1/800 leading to a Courant number of ν = t /x = 0.4. Figs. 4(a) and 4(b) show the difference of the numerical solution to the exact solution for each cell at t = 1. The mass preserving full discretizations lead to significant error peaks at the block interfaces indicated by the dashed vertical lines. Opposed to this the consistent discretizations show a much smoother behavior with negligible effects of the interfaces. These effects are also the reason for the mass preserving schemes showing significantly higher maximum errors compared to the latter, see Table 5. The mass defect for both consistent discretizations is in the order of magnitude 10−8 . Due to round-off errors the mass preserving discretizations also lead to mass defects about the magnitude of the machine epsilon, 10−17 . To compare the third order multirate scheme listed in Table 2 to the schemes from the literature the time step of the second order schemes is reduced by a factor of 1/2 so that the solutions are obtained with the same computational cost, measured in evaluations of the right-hand side. Again the spatial resolution is x = 1/800, Courant numbers are ν = 0.2 for the second order schemes and ν = 0.4 for the third order scheme; see Figs. 5(a) and 5(b) for results. Even for these moderate Courant numbers the third order scheme proves superior to the second order schemes, yielding errors which are 50% to 80% smaller compared to the second order schemes, comp. Table 5. With a second test setting we want to examine the properties of our multirate IMEX approach. Our test equation now is an advection–reaction equation.
∂w ∂w =− + f ( w ). ∂t ∂x The transported variable w now consists of three components which may be interpreted as three different chemical compounds. The reaction term f is the well-known Robertson example [12]:
⎞
⎞ −0.04w 1 + 104 w 2 w 3 2 4 7 f ( w ) = ⎝ f 2 ( w 1 , w 2 , w 3 ) ⎠ = ⎝ 0.04w 1 − 10 w 2 w 3 − 3 · 10 w 2 ⎠ , f 3(w1, w2, w3) 3 · 107 w 22 ⎛
f 1(w1, w2, w3)
⎛
the initial condition is given by
w (t = 0) =
max(sin(2π x), 0)2 0 0
.
1540
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
Fig. 4. Comparison of second order schemes.
Table 5 Difference of numerical to exact solution in L 1 and L ∞ norm at t = 1. Number of cells Mass preserving discretizations
Consistent discretizations
100
200
Constantinescu and Sandu (RK2a)
L1 L∞
3.33 · 10−4 6.41 · 10−4
7.68 · 10−5 16.81 · 10−5
400 1.86 · 10−5 5.85 · 10−5
800 4.58 · 10−6 24.12 · 10−6
RFSMR (RK2a) flux splitting
L1 L∞
3.00 · 10−4 6.27 · 10−4
7.31 · 10−5 21.95 · 10−5
1.86 · 10−5 10.43 · 10−5
4.79 · 10−6 59.57 · 10−6
RFSMR (RK43) flux splitting
L1 L∞
1.35 · 10−4 3.17 · 10−4
1.86 · 10−5 7.05 · 10−5
0.42 · 10−5 2.62 · 10−5
1.06 · 10−6 12.93 · 10−6
Tang and Warnecke
L1 L∞
2.71 · 10−4 5.61 · 10−4
6.39 · 10−5 10.25 · 10−5
1.61 · 10−5 2.58 · 10−5
4.06 · 10−6 6.45 · 10−6
RFSMR (RK2a) cell splitting
L1 L∞
2.65 · 10−4 5.78 · 10−4
6.32 · 10−5 10.34 · 10−5
1.60 · 10−5 2.59 · 10−5
4.08 · 10−6 6.52 · 10−6
RFSMR (RK43) cell splitting
L1 L∞
1.26 · 10−4 2.81 · 10−4
1.49 · 10−5 4.28 · 10−5
0.29 · 10−5 0.86 · 10−5
0.73 · 10−6 2.06 · 10−6
RK43 solved with t /x = 0.4, remaining methods solved with t /x = 0.2, leading to equal computational cost for each column. For better readability the exponents are uniform for each column.
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
1541
Fig. 5. Comparison of second order schemes from literature against third order multirate scheme. Solutions with the second order schemes were obtained using the half time step.
Table 6 A third order diagonally implicit Runge–Kutta scheme.
√
1/2 − 3/6 √ 1/2 + 3/6
√
1/2 − 3/6 √ 1/ 3 1/2
√
1/2 − 3/6 1/2
Again we consider the unit interval in space with periodic boundary conditions. The grid size now is 1/100, the advection operator is discretized using the third order upwind biased spatial discretization with flux limiter [5] to guarantee positivity. Explicit multirate integration is applied by solving the first 50 cells using the fast method and the remaining half using the slow method. A mass preserving splitting by fluxes is employed. For sake of comparison we also examine a “pure” (i.e. non-multirate) IMEX scheme based on the same methods. Note that formally the latter scheme should yield larger errors as the larger advection time step is employed on all cells. We employ a third order multirate IMEX scheme based on RK43 (see Table 2). Starting from this explicit multirate scheme we obtain a multirate IMEX scheme by employing the diagonally implicit Runge–Kutta method listed in Table 6 as the inner method in an additional RFSMR construction stage. To examine the error behavior of this discretization we consider the solution at t = 1. As in this case no exact solution is available we obtain our reference solution by solving the system with a time step of t = 10−5 with the pure IMEX scheme. The difference of solutions obtained with larger time
1542
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
Fig. 6. Error to reference, obtained with different time steps.
steps to this reference in L 1 norm are shown in Fig. 6. The average slope of the graphs in the double logarithmic diagram affirm the order of convergence that was shown formally in Section 3. 5. Conclusions and outlook In this paper we presented and examined a generic explicit Runge–Kutta multirate scheme, named recursive flux splitting multirate (RFSMR). The created methods can be cast in partitioned Runge–Kutta form, they are internally consistent and of second order accuracy for second order base methods. Even third order convergence can be obtained if the base methods satisfy additional order conditions. A third order multirate scheme was presented and has proven superior to second order schemes even at moderate Courant numbers. A further test was dedicated to the multirate IMEX approach. Numerical results affirm that third order convergent IMEX schemes as well as multirate IMEX schemes can be constructed using the RFSMR approach. An implementation of RFSMR in a three-dimensional air pollution model developed and employed at the Leibniz Institute for Tropospheric Research showed promising results. The multirate approach was employed only on the advection operator causing only a vanishing fraction of the total computational cost in a full scale simulation. Anyway even for simulations involving complex chemistry simulations a speedup comparable to the speedup of pure advection can be obtained. This is mainly due to the local nature of the chemistry computations. The advection operator influences chemistry in the form of a source term included in the implicit integration. Every stage of the explicit Runge–Kutta stage changes this source term and thus causes a discontinuity in the right-hand side of the equation solved implicitly. For higher order implicit solvers this makes a computationally expensive restart procedure necessary. Application of a multirate scheme leads to less frequent updates of the right-hand side which leads to fewer discontinuities and finally to an increased efficiency of the implicit solving routine. A separate paper on details of this implementation is in preparation. This will also include algorithmic details and discussions on the efficient implementation and special challenges arising in the context of parallel programming. Acknowledgement This work was funded by the German Research Foundation (DFG), grant No. WO 1204/1-2. References [1] E. Constantinescu, A. Sandu, On Adaptive Mesh Refinement for Atmospheric Pollution Models, Lecture Notes in Comput. Sci., vol. 3515, Springer, Berlin/Heidelberg, 2005, pp. 798–805. [2] E. Constantinescu, A. Sandu, Multirate timestepping methods for hyperbolic conservation laws, J. Sci. Comput. 33 (3) (2007) 239–278. [3] W. Dabberdt, M. Carroll, D. Baumgardner, G. Charmichael, R. Cohen, T. Dye, J. Ellis, G. Grell, S. Grimmond, S. Hanna, J. Irwin, B. Lamb, S. Madronich, J. McQueen, J. Meagher, T. Odman, J. Pleim, H. Schmid, D. Westphal, Meteorological research needs for improved air quality forecasting, Bull. Am. Meteorol. Soc. 85 (2004) 563–586. [4] C. Gear, D. Wells, Multirate linear multistep methods, BIT Numer. Math. 24 (4) (1984) 484–502. [5] W. Hundsdorfer, B. Koren, M. van Loon, J. Verwer, A positive finite-difference advection scheme, J. Comput. Phys. 117 (1995) 35–46. [6] W. Hundsdorfer, A. Mozartova, V. Savcenco, Analysis of explicit multirate and partitioned Runge–Kutta schemes for conservation laws, CWI Report MAS-E0715, Centrum voor Wiskunde en Informatica, Amsterdam, 2007. [7] W. Hundsdorfer, J. Verwer, Numerical Solution of Time-Dependent Advection–Diffusion Reaction Equations, Springer Ser. Comput. Math., Springer, Berlin, Heidelberg, New York, 2003.
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
[8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
1543
Z. Jackiewicz, R. Vermiglio, Order conditions for partitioned Runge–Kutta methods, Appl. Math. 45 (4) (1998) 301–316. S.-S. Jiang, C.-W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228. O. Knoth, R. Wolke, Implicit–explicit Runge–Kutta methods for computing atmospheric reactive flows, Appl. Numer. Math. 28 (2–4) (1998) 327–341. S. Osher, R. Sanders, Numerical approximations to nonlinear conservation laws with locally varying time and space grids, Math. Comp. 41 (1983) 321–336. H.H. Robertson, Numerical Analysis: An Introduction, Academic Press, 1966, pp. 178–182 (Chap. “The solution of a set of reaction rate equations”). M. Schlegel, O. Knoth, M. Arnold, R. Wolke, Multirate Runge–Kutta schemes for advection equations, J. Comput. Appl. Math. 226 (2009) 345–357. H. Tang, G. Warnecke, A class of high resolution schemes for hyperbolic conservation laws and convection–diffusion equations with varying time and space grids, SIAM J. Sci. Comput. 26 (4) (2005) 1415–1431. R. Wang, R. Spiteri, Linear instability of the fifth-order WENO method, SIAM J. Numer. Anal. 45 (5) (2007) 1871–1901. L. Wicker, W. Skamarock, A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta time differencing, Mon. Weather Rev. 126 (1998) 1992–1999. R. Zvan, P. Forsyth, K. Vetzal, Robust numerical methods for PDE models of Asian options, J. Comput. Finance 1 (1998) 39–78.
Contents lists available at SciVerse ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
Numerical solution of multiscale problems in atmospheric modeling Martin Schlegel a,∗ , Oswald Knoth a , Martin Arnold b , Ralf Wolke a a b
Institute for Tropospheric Research, Permoserstrasse 15, 04318 Leipzig, Germany Martin Luther University Halle-Wittenberg, Institute for Mathematics, 06099 Halle (Saale), Germany
a r t i c l e
i n f o
Article history: Available online 7 June 2012 Keywords: Time integration Runge–Kutta method Multirate scheme Flux splitting IMEX
a b s t r a c t Explicit time integration methods are characterized by a small numerical effort per time step. In the application to multiscale problems in atmospheric modeling, this benefit is often more than compensated by stability problems and stepsize restrictions resulting from stiff chemical reaction terms and from a locally varying Courant–Friedrichs–Lewy (CFL) condition for the advection terms. In the present paper, we address this problem by a rather general splitting technique that may be applied recursively. This technique allows the combination of implicit and explicit methods (IMEX splitting) as well as the local adaptation of the time stepsize to the meshwidth of non-uniform space grids in an explicit multirate discretization of the advection terms. Using a formal representation as partitioned Runge–Kutta method, convergence of order p 3 may be shown if some additional order conditions are satisfied. In a series of numerical tests, the convergence results are verified and the consequences of different splitting strategies like flux splitting and cell splitting are analysed. © 2012 Published by Elsevier B.V. on behalf of IMACS.
1. Introduction State of the art numerical models can be used to simulate a broad spectrum of atmospheric phenomena. These range from microphysical and chemical processes with characteristic times below a second and characteristic scales in the submillimeter range up to high/low pressure areas and planetary waves on scales about 1000 kilometers and characteristic times up to months. Additionally the simulation domain is often locally refined to allow for a more accurate examination of regions of interests [1,3]. These may be densely populated areas in contrast to surrounding woodland or agricultural regions, or in a global model continents in contrast to oceans. Atmospheric phenomena can be described using the advection–diffusion–reaction equation
∂ ∂ ∂ ∂ c + f (x, t , c ), ρD c = − (uc ) + ∂t ∂x ∂x ∂x ρ
(1)
with the concentration of a transported quantity c, wind speed u, density ρ , diffusion coefficient D and a reaction term f . The latter may also include sources and sinks. The advection part of this equation can be solved efficiently using explicit Runge–Kutta (ERK) time integration methods. These methods have in common that stability requirements limit the global time step to be smaller than some critical value proportional to the ratio of grid size and the magnitude of the wind speed for each cell [2,7,17]. Consequently the smallest cell, or more exactly the cell with the smallest characteristic time determines the global time step.
*
Corresponding author. E-mail addresses: [email protected] (M. Schlegel), [email protected] (O. Knoth), [email protected] (M. Arnold), [email protected] (R. Wolke). 0168-9274/$36.00 © 2012 Published by Elsevier B.V. on behalf of IMACS. http://dx.doi.org/10.1016/j.apnum.2012.06.023
1532
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
Multirate time integration methods can be employed to adapt the time step locally so that slower components take longer and fewer time steps, resulting in a moderate to substantial reduction of the computational cost, depending on the scenario to simulate [4]. Multirate approaches have been developed since the early 1980s. Osher and Sanders [11] presented a scheme allowing multirate Euler steps. More current approaches allow for a wide variety of multirate schemes by proposing generic multirate methods based on classical ERKs. In 2005 Tang and Warnecke [14] proposed different internally consistent multirate schemes. These schemes however are not mass preserving. The approach of Constantinescu and Sandu [2] yields multirate ERK schemes which are mass preserving and at most of second order accuracy. In 2007 Hundsdorfer et al. showed that multirate schemes which in absence of a slow domain reduce to multiple applications of the underlying base method cannot be mass preserving and consistent at the same time [6]. The diffusion and reaction terms in (1) on the other hand are usually stiff and consequently cannot be solved efficiently using explicit methods. Instead implicit methods must be employed. Splitting methods solving the non-stiff terms using explicit methods and the stiff terms using implicit methods are known as implicit/explicit (IMEX) methods. This class of splitting methods has proven efficient for problems which can be split into non-stiff (or mildly stiff) and stiff terms. We presented an internally consistent (in the nomenclature of [6]) generic multirate scheme for constructing third order accurate multirate methods [13]. Since for atmospheric modeling mass preservation or more generally the preservation of linear invariants is of major interest we employ a mass preserving splitting of the advection operator, thus trading consistency (i.e. truncation error O (t ) + O (x)) for mass preservation. In Section 2 we shall exemplify the construction of multirate schemes using our algorithm. The possible extension to an implicit/explicit (IMEX) multirate scheme will also be considered. Section 3 is dedicated to the formal analysis of the constructed methods. Numerical tests will be presented in Section 4. 2. A new class of hierarchical multirate IMEX schemes The multirate scheme we shall introduce is based on an implicit/explicit (IMEX) method introduced by Knoth and Wolke in [10] for the efficient solution of advection–diffusion–reaction equations in air pollution. Consider a splitting of the righthand side
w = F ( w ) + G ( w ),
w (0) = w 0 ,
w ∈ RN .
(2)
In the context of this splitting G is a non-stiff term which will be solved explicitly whereas F is a stiff term which must be solved using an implicit method. For the theoretical analysis of the method the authors assume that the implicit integration is carried out exactly. Defining the nodes of the explicit method
τi = tn + tc i , the full algorithm computing an approximate solution w n+1 for time tn+1 from an approximate solution w n at time tn is defined as follows
W 1 = wn, ri =
i −1
(3)
(ai j − ai −1, j )G ( W j ),
(4)
j =1
v i (τi −1 ) = W i −1 , ∂ vi 1 = r i + F ( v i ), ∂τ c i − c i −1
(5)
τ ∈ [τi−1 , τi ], i = 2, . . . , s + 1,
(6)
W i = v i (τi ),
(7)
w n +1 = W s +1 ,
(8)
with s denoting the number of stages of an explicit Runge–Kutta method with parameters (a, b, c ) in common notation. For simplicity we define as+1, j = b j , thus making a separate treatment of the summation stage unnecessary. We note that the cumulative implicit integration time is t. There is another class of comparable IMEX schemes where the integration for the stiff part starts in each internal stage at the old time stage tn and the initial value from the previous integration step and therefore the cumulative interval is usually larger than the explicit time step [16]. If the stiff term is identical to zero F ≡ 0 we obtain the underlying explicit method applied to w = G ( w ). To obtain a multirate scheme we again consider a splitting of the right-hand side,
w = G 1 ( w ) + G 2 ( w ),
w (0) = w 0 ,
w ∈ RN ,
(9)
but now both terms (now denoted G 1 and G 2 ) are assumed to be non-stiff with G 2 imposing a more restrictive time step limitation. Referring to the characteristic times of these methods we will also call G 1 a slow term opposed to the fast term G 2 . Now two potentially but not necessarily different explicit Runge–Kutta methods are employed. The method replacing
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
1533
the implicit integrator is called the inner method, denoted by a superscript I. The other method is called the outer method, denoted by a superscript O. The multirate algorithm reads
W 1 = wn, ri =
i −1
(10)
j =1
O aO i j − ai −1, j G 1 ( W j ),
(11)
V i ,1 = W i −1 ,
(12)
V i ,k = V i ,k−1 + t c iO − c iO−1
k −1 j =1
I akj − akI −1, j
1 c iO − c iO−1
ri + G 2 ( V i, j )
(13)
i −1 k −1 O O I O = V i ,k−1 + t ckI − ckI −1 a i j − aO akj − akI −1, j G 2 ( V i , j ), i −1, j G 1 ( W j ) + t c i − c i −1 j =1
(14)
j =1
W i = V i , s I +1 ,
(15)
w n +1 = W s O +1 = V s O +1 , s I +1 .
(16)
We call this algorithm recursive flux splitting multirate (RFSMR), [13]. The name emphasizing the splitting by fluxes comes from the original application to the advection equation in atmospheric simulations where mass preservation is essential. The algorithm (10)–(16) can also be cast as a partitioned Runge–Kutta (PRK) system:
W 1 = wn, W i = w n + t
i −1
aslow G 1 ( W j ) + afast ij i j G 2 (W j ) ,
j =1
w n+1 = w n + t
s
bslow G 1 ( W j ) + bfast j j G 2 (W j ) ,
j =1
which is the general form of an explicit PRK applied on a splitting of the right-hand side. We will make use of this notation in the following chapter. The Runge–Kutta parameters (aslow , bslow , c slow ) and (afast , bfast , c fast ) employed here depend on the parameters of the base methods as shown below. For easier association of the parameters of the partitioned method with the parameters of the underlying base methods we employ tuple indexes. The first element of a tuple is associated with a stage of the outer method, the other element corresponds to an inner method stage:
( i , k ) ↔ ( i − 1) s I + k . The following equations can be deduced by transforming algorithm (10)–(16) from incremental to common form and solving the resulting recurrence relations. The complete set of multirate PRK coefficients reads:
s = sO · sI ,
= aslow aslow (i ,k),( j ,l) = (i −1)sI +k,( j −1)sI +l
bslow ( j ,l) =
bOj
if l = 1,
0
if l > 1,
(17) aO i, j
+ (aOi +1, j
0
− aOi , j )ckI
if l = 1, if l > 1,
(18) (19)
O I O O c (slow i ,k) = c i + ck c i +1 − c i , ⎧ ⎪0
(20)
afast (i ,k),( j ,l)
(21)
=
⎨
if j > i , (c Oj +1 − c Oj )akI ,l if j = i ,
⎪ ⎩ (c O
− c Oj )blI if j < i , j +1 O O I bfast ( j ,l) = c j +1 − c j bl , slow O I O O c (fast i ,k) = c (i ,k) = c i + ck c i +1 − c i .
(22) (23)
Note that different pairs of base methods may result in equivalent multirate methods. A possible disadvantage of the RFSMR algorithm is that there is no parameter defining the time step ratio, i.e. the relative stability and corresponding to that the relative computational cost of the slow and the fast method. Instead the slow method is equivalent to the outer base method while the fast method is equivalent to the inner method executed sO times with
1534
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
Table 1 RFMSMR(RK2a) – an example for a second order multirate scheme; redundant stages omitted. 0 1
0 1/2 1/2 1
1 1/2 1/2 (RK2a)
outer base method 0 1/2 1/2 1 1
1/2 1/2 1 1 1/2
0 0 0 0 0 0 0 “slow” part
0 0
1/2 1/4 1/4 1/4 1/4 1/2 1/4 1/4 1/4 inner base method 0 1/2 1/2 1 1
1/2
1/2 1/4 1/4 1/4 1/4
1/4
1/4 1/4 1/2 1/4 1/4 1/4 1/4 “fast” part
1/4 1/4
0
Fig. 1. Schematic of recursive construction.
time step sizes of (c O − c Oi )t; i ∈ {2, . . . , sO + 1}. This disadvantage may however be easily overcome by substituting the i +1 inner method by a composition of n steps of step size t /n. For instance in Table 1, the inner method is equivalent to the outer method executed twice with a half time step, in combination of the outer method’s nodes c = (0, 1, 1) leading to a time step ratio of R = 2 for the partitioned multirate method. For practical purposes two temporal refinement levels are often insufficient. For instance many temporal refinement levels are desirable when nested spatial refinement is employed, as is common in atmospheric modeling. To obtain a multirate scheme with an arbitrary number of temporal refinement levels the recursive flux-splitting multirate approach may be applied recursively. One of the methods created in one construction stage then serves as outer method in the following construction stage. Note that the PRKs created using the RFSMR algorithm are internally consistent (i.e. c slow = c fast ) and the parameters of the fast method do not depend directly on aO or bO , comp. (21)–(23). Consequently it makes no difference if the slow or the fast method created in the previous construction stage is employed as outer method in the next stage. Like this an arbitrary number of temporal refinement levels can be realized, see Fig. 1. Arrows indicate use of variables. The ellipses marked “slow construction” denote construction of the slow method’s parameters as defined in (18)–(20); analogously the “fast construction” ellipses are correlated to (21)–(23). Complex multiphysics simulations often involve both terms which may be approximated efficiently using explicit methods and stiff terms requiring an implicit integrator. Based on an IMEX splitting, our explicit multirate algorithm may easily be extended to a multirate IMEX algorithm combining the advantages of both splitting techniques. As before consider a splitting of the right-hand side
w = G 1 ( w ) + G 2 ( w ) + F ( w ), now with a stiff term F and two non-stiff terms G 1 , G 2 with different characteristic times, G 1 being the slower one. A multirate IMEX splitting may now be obtained by hierarchical application of the IMEX splitting (3)–(8). Again we define intermediate times, now correlated to the stages of the smaller explicit time step
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
1535
Fig. 2. Illustration of intermediate points in time.
τi,k = tn + t c iO−1 + c iO − c iO−1 ckI = tn + tc ifast ,k , see also Fig. 2. First we apply the splitting on w = G 1 ( w ) + [G 2 ( w ) + F ( w )]:
W 1 = wn, ri =
i −1
(24)
j =1
O aO i j − ai −1, j G 1 ( W j ),
(25)
v i (τi ,1 ) = W i −1 ,
(26)
∂ vi 1 = O r i + G 2 ( v i ) + F ( v i ), ∂τ c i − c iO−1
τ ∈ [τi,1 , τi+1,1 ], i = 2, . . . , sO + 1,
(27)
W i = v i (τi +1,1 ),
(28)
w n +1 = W s O +1 .
(29)
The same splitting (3)–(8) is applied to (27). Note that the time dependent variable v employed in (27) now becomes a Runge–Kutta stage variable V , with V i ,1 ≈ v (τi ,1 ) and V i ,sI +1 ≈ v (τi +1,1 ).
V i ,1 = W i −1 ,
(30)
u i ,k (τi ,k−1 ) = V i ,k−1 , r i ,k =
(31)
k −1 j =1
I akj − akI −1, j
1 c iO − c iO−1
∂ u i ,k 1 = I r i ,k + F (u i ,k ), ∂t ck − ckI −1
ri + G 2 ( V i, j ) ,
(32)
τ ∈ [τi,k−1 , τi,k ], k = 2, . . . , sI + 1,
(33)
V i ,k = u i ,k (τi ,k ),
(34)
W i = V i , s I +1 .
(35)
3. Order conditions As shown in the preceding section the RFSMR algorithm may be cast in PRK form. This gives us the means to examine the order of the created multirate methods in the well established context of partitioned Runge–Kutta methods. Order conditions for PRKs are composed of the conditions for the separate methods:
Order 1:
s
b i = 1,
(36)
i =1
Order 2:
s
bi ci =
i =1
Order 3:
s
b i c i2 =
i =1 s s i =1 j =1
1 2
(37)
,
1 3
(38)
,
b i ai j c j =
1 6
,
(39)
1536
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
and additional coupling conditions (see e.g. [8]) setting the parameters of the different methods into relation:
Order 2:
s
s
bslow c ifast = i
i =1
Order 3:
s
slow bfast = i ci
i =1
bslow c ifast i
2
=
i =1
s
bfast c islow i
fast bslow afast = i ij c j
s s
i =1 j =1
(40)
,
2
1
= , 3
slow slow bfast cj = i ai j
i =1 j =1
bfast c ifast − c islow i
2
=
i =1 s
2
i =1
s s
s
1
s
bslow c ifast − c islow i
1 6
2
(41)
,
= 0,
i =1
s
fast fast bfast c i − c islow = i ci
i =1
bslow c islow c ifast − c islow = 0, i
i =1
s s
fast
bfast i ai j
slow
c fast −cj j
=
i =1 j =1 s s
s s
= 0, bslow aslow c fast − c slow i ij j j
i =1 j =1
slow fast = bfast c j − c slow i ai j j
i =1 j =1
s s
= 0. bslow afast c fast − c slow i ij j j
i =1 j =1
As the methods generated using RFSMR are internally consistent (i.e. c slow = c fast ) all coupling conditions apart from (41) are redundant. Note that if the method would additionally preserve all linear invariants (i.e. bslow = bfast ) the second order conditions (37), (40) would be satisfied automatically. We can prove order conditions for the generated partitioned methods assuming that the underlying base methods satisfy certain conditions. Theorem 3.1. If both the inner and the outer method satisfy the conditions for second order consistency (36), (37) the generated partitioned method is of second order accuracy, i.e. the fast and the slow method satisfy (36), (37). Corollary 3.2. A second order multirate scheme with an arbitrary number of temporal refinement levels can be constructed based on a set of explicit Runge–Kutta methods if all of the employed base methods satisfy the conditions for second order consistency (36), (37). Proof of Theorem 3.1. Proving this theorem is straightforward; apply definitions (17)–(23) on the order conditions and simplify the resulting expression. At some points we will employ the assumptions on the base methods. This will be denoted by under braces. Order 1 condition (36) for slow method: O I
s s
O
bslow = j
j =1
s
bOj = 1.
j =1
=1
Order 1 condition (36) for fast method: O I
s s j =1
O
bfast j
=
I
s s
c Oj +1
j =1 k =1
− c Oj
bkI
=
sO
j =1
c Oj +1
s +1
O I
s s j =1
O
bslow c slow = j j
s
bOj c Oj =
j =1
=1/2
1 2
.
=c OO
Order 2 condition (37) for slow method:
− c Oj =1
I
s
bkI .
k =1
=1
= 1.
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
1537
Order 2 condition (37) for fast method: O I
s s
O
fast bfast = j cj
j =1
I
s s j =1 k =1 O
=
I
s s O j =1 k =1
=
2
c j c Oj +1 − c Oj bkI + c Oj +1 − c Oj bkI ckI
sO O j =1
c Oj +1 − c Oj bkI · c Oj + c Oj +1 − c Oj ckI
I
c j c Oj +1
− c Oj
s
bkI +
k =1
c Oj +1
2 − c Oj
=
j =1
2
c Oj +1 c Oj − c Oj
I
s
bkI ckI
k =1
=1
sO
=1/2
+
2 1 O 2 c j +1 − 2c Oj +1 c Oj + c Oj 2
sO
1 O 2 O 2 1 O 2 1 = = c j +1 − c j c O +1 = . 2 2 s 2 j =1
2
=1
An analogous theorem can be formulated for the conditions for third order consistency (38), (39). Trying to solve for the coupling condition (41) for third order consistent partitioned methods however we find an additional condition to be satisfied by the outer method: s i −1 1 (c i +1 − c i ) (ai +1, j + ai j )c j = . i =1
j =1
3
(42)
Note that here this equation is applied to the parameters of the outer method
(a, b, c ) = aO , bO , c O . The more general notation is chosen because it will be applied to parameters of other methods as well. Condition (42) was already derived by Knoth and Wolke in [10] for the IMEX splitting introduced there. However the authors assumed that the implicit integrator (in the current context replaced by the inner explicit method) was exact. In the same paper also explicit Runge–Kutta methods satisfying this condition were presented. One of those methods and the multirate PRK resulting from it are listed in Table 2. Theorem 3.3. If both the inner and the outer method satisfy the conditions for third order consistency (36)–(39) and the outer method additionally satisfies the new order condition (42) the generated partitioned method is of third order accuracy, i.e. the fast and the slow method satisfy (36)–(41). If both of the base methods satisfy the conditions for third order consistency (36)–(39) and both methods satisfy the additional condition (42) both the fast and the slow method generated satisfy (42). The proof of this theorem follows the same idea as the proof of Theorem 3.1 but is very technical and would go beyond the scope of this paper. A last result following directly from Theorem 3.3 shall be given: Corollary 3.4. A third order multirate scheme with an arbitrary number of temporal refinement levels can be constructed based on a set of explicit Runge–Kutta methods if all of the employed base methods satisfy the conditions for third order consistency (36)–(39) and all but the last inner method satisfy (42). Consequently we are also able to construct a multirate IMEX scheme of third order accuracy employing an implicit Runge–Kutta method satisfying the classical conditions for third order consistency (36)–(39) as last inner method in a construction cascade as illustrated in Fig. 1. 4. Numerical tests In the context of partitioned methods it is of interest how exactly the splitting is defined. With a finite volume approach conservation laws
1538
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
Table 2 RFSMR(RK43) An example for a third order multirate scheme; redundant stages omitted. 0 1/2 1/2 1
1/2
−1/6
2/3 1/3 −1/3 1 1/6 1/3 1/3 1/6 base method (inner and outer) (RK43)
0 1/4 1/4 1/2 1/2 1/2 3/4 3/4 1 1
0 1/4 1/4 1/2 1/2 1/2 3/4 3/4 1 1
1/4 1/4 1/2 1/2 −1/6 1/12 1/12 1/3 1/3 1/6
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 2/3 0 1/6 0 1/6 0 −1/3 0 −1/3 0 1/3 “slow” part
1/2 1/2 1 1 1/3
0 0 0 0
0 0 0
0 0
1/6
1/4
−1/12
1/3
1/6 1/12 1/12 1/12 1/12 1/12 1/12
−1/6 1/6 1/6 1/6 1/6 1/6 1/6
1/2 1/6 1/6 1/6 1/6 1/6 1/6
1/12
1/6
1/6
1/12 1/12 1/12 1/12 1/12 1/12
0 0 0 0 0
1/4
−1/12 1/6 1/12
1/12 0 1/12 “fast” part
1/3
−1/6 1/6
1/2 1/6
1/12
1/6
1/6
1/12
0
Table 3 Multirate scheme according to Tang and Warnecke [14]. 0 1/2 1/2 1
1/2 1/4 1/4 1 0 0 1/2 0 0 slow method
0 1/2 1/2 1 1/2
1/2 1/4 1/4 1/4 fast
1/4 1/4 1/2 1/4 1/4 method
1/4
Table 4 Multirate scheme according to Constantinescu and Sandu [2]. 0 1 0 1
1 0 0 1/4
0 0 1 1/4 1/4 slow method
0 1/2 1/2 1 1/4
1/2 1/4 1/4 1/4 fast
1/4 1/4 1/2 1/4 1/4 method
1/4
∂w f (w) =− ∂t ∂x naturally lead to a so-called flux form
∂ wi 1 = F i w (t ) = − f i +1/2 w (t ) − f i −1/2 w (t ) ∂t xi with cell indexes i and the numerical fluxes f i ±1/2 at the boundaries of the cells. There are several possibilities to split the operator F = F 1 + F 2 . The most commonly used approach is the so-called splitting by cells alias splitting by components. This can be interpreted in such a way that F 1 describes the net fluxes of a subset of all cells and F 2 describes the net fluxes of the remaining cells. If not mentioned otherwise we tacitly assume a splitting by components. Alternatively a so-called splitting by fluxes can be employed. This approach can be interpreted so that F 1 describes the fluxes over a subset of all cell boundaries and F 2 describes the fluxes over the remaining cell boundaries. Opposed to the splitting by components this guarantees mass preservation independent from the specific time integration method as both F 1 and F 2 are mass preserving. However as already shown by Hundsdorfer et al. in [6] applying a splitting by fluxes will not lead to a consistent full discretization.
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
1539
Fig. 3. Decomposition of the unit interval in slow and fast blocks.
Though called recursive flux splitting multirate, the multirate schemes constructed via our algorithm may also be applied on a splitting by cells; thereby trading mass preservation for consistency. The following test setting is intended to examine the effects of these different splitting approaches when the resulting split scheme is solved using an explicit multirate method generated using RFSMR. It also is a mean to compare RFSMR to multirate methods from the literature, one given by Tang and Warnecke in [14] the other by Constantinescu and Sandu in [2]. (See Tables 3 and 4.) Both of these methods are based on RK2a (see Table 1) and converge with second order. Our test equation is one-dimensional advection with unit speed and a smooth initial profile
∂w ∂w =− , ∂t ∂x
w (0) = sin2 (2π x),
x ∈ [0, 1],
with periodic boundary conditions. The system is discretized in space employing an equidistant grid and the WENO5 approximation [9] for fluxes. In [15], Wang and Spiteri formally showed that this approximation is not linearly stable if combined with certain explicit Runge–Kutta methods, including RK2a. Numerical evidence for instability in [15], Section 4.1 was given for a Courant number of ν = 1.32. We did however not observe this behavior for the Courant numbers ν ∈ {0.2, 0.4} considered in the following tests. To examine the properties of the various multirate methods the grid is split into 5 blocks solved using the respective slow method and 5 blocks solved using the fast method, see Fig. 3. We consider the difference between the numerical solution at final time t = 1 and the exact solution. In [6], Hundsdorfer et al. presented a similar setting intended to examine the differences of consistent and mass preserving multirate schemes. As pointed out in the same paper, a true multirate scheme cannot be consistent and mass preserving at the same time. In the nomenclature of the authors a true multirate scheme is a scheme for which the respective slow method is equivalent to the base method in terms of computational cost and stability and the fast method is equivalent to the base method executed R times with time steps t / R for a given time step ratio R. First the second order schemes are examined. We choose t = 5 · 10−4 ; x = 1/800 leading to a Courant number of ν = t /x = 0.4. Figs. 4(a) and 4(b) show the difference of the numerical solution to the exact solution for each cell at t = 1. The mass preserving full discretizations lead to significant error peaks at the block interfaces indicated by the dashed vertical lines. Opposed to this the consistent discretizations show a much smoother behavior with negligible effects of the interfaces. These effects are also the reason for the mass preserving schemes showing significantly higher maximum errors compared to the latter, see Table 5. The mass defect for both consistent discretizations is in the order of magnitude 10−8 . Due to round-off errors the mass preserving discretizations also lead to mass defects about the magnitude of the machine epsilon, 10−17 . To compare the third order multirate scheme listed in Table 2 to the schemes from the literature the time step of the second order schemes is reduced by a factor of 1/2 so that the solutions are obtained with the same computational cost, measured in evaluations of the right-hand side. Again the spatial resolution is x = 1/800, Courant numbers are ν = 0.2 for the second order schemes and ν = 0.4 for the third order scheme; see Figs. 5(a) and 5(b) for results. Even for these moderate Courant numbers the third order scheme proves superior to the second order schemes, yielding errors which are 50% to 80% smaller compared to the second order schemes, comp. Table 5. With a second test setting we want to examine the properties of our multirate IMEX approach. Our test equation now is an advection–reaction equation.
∂w ∂w =− + f ( w ). ∂t ∂x The transported variable w now consists of three components which may be interpreted as three different chemical compounds. The reaction term f is the well-known Robertson example [12]:
⎞
⎞ −0.04w 1 + 104 w 2 w 3 2 4 7 f ( w ) = ⎝ f 2 ( w 1 , w 2 , w 3 ) ⎠ = ⎝ 0.04w 1 − 10 w 2 w 3 − 3 · 10 w 2 ⎠ , f 3(w1, w2, w3) 3 · 107 w 22 ⎛
f 1(w1, w2, w3)
⎛
the initial condition is given by
w (t = 0) =
max(sin(2π x), 0)2 0 0
.
1540
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
Fig. 4. Comparison of second order schemes.
Table 5 Difference of numerical to exact solution in L 1 and L ∞ norm at t = 1. Number of cells Mass preserving discretizations
Consistent discretizations
100
200
Constantinescu and Sandu (RK2a)
L1 L∞
3.33 · 10−4 6.41 · 10−4
7.68 · 10−5 16.81 · 10−5
400 1.86 · 10−5 5.85 · 10−5
800 4.58 · 10−6 24.12 · 10−6
RFSMR (RK2a) flux splitting
L1 L∞
3.00 · 10−4 6.27 · 10−4
7.31 · 10−5 21.95 · 10−5
1.86 · 10−5 10.43 · 10−5
4.79 · 10−6 59.57 · 10−6
RFSMR (RK43) flux splitting
L1 L∞
1.35 · 10−4 3.17 · 10−4
1.86 · 10−5 7.05 · 10−5
0.42 · 10−5 2.62 · 10−5
1.06 · 10−6 12.93 · 10−6
Tang and Warnecke
L1 L∞
2.71 · 10−4 5.61 · 10−4
6.39 · 10−5 10.25 · 10−5
1.61 · 10−5 2.58 · 10−5
4.06 · 10−6 6.45 · 10−6
RFSMR (RK2a) cell splitting
L1 L∞
2.65 · 10−4 5.78 · 10−4
6.32 · 10−5 10.34 · 10−5
1.60 · 10−5 2.59 · 10−5
4.08 · 10−6 6.52 · 10−6
RFSMR (RK43) cell splitting
L1 L∞
1.26 · 10−4 2.81 · 10−4
1.49 · 10−5 4.28 · 10−5
0.29 · 10−5 0.86 · 10−5
0.73 · 10−6 2.06 · 10−6
RK43 solved with t /x = 0.4, remaining methods solved with t /x = 0.2, leading to equal computational cost for each column. For better readability the exponents are uniform for each column.
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
1541
Fig. 5. Comparison of second order schemes from literature against third order multirate scheme. Solutions with the second order schemes were obtained using the half time step.
Table 6 A third order diagonally implicit Runge–Kutta scheme.
√
1/2 − 3/6 √ 1/2 + 3/6
√
1/2 − 3/6 √ 1/ 3 1/2
√
1/2 − 3/6 1/2
Again we consider the unit interval in space with periodic boundary conditions. The grid size now is 1/100, the advection operator is discretized using the third order upwind biased spatial discretization with flux limiter [5] to guarantee positivity. Explicit multirate integration is applied by solving the first 50 cells using the fast method and the remaining half using the slow method. A mass preserving splitting by fluxes is employed. For sake of comparison we also examine a “pure” (i.e. non-multirate) IMEX scheme based on the same methods. Note that formally the latter scheme should yield larger errors as the larger advection time step is employed on all cells. We employ a third order multirate IMEX scheme based on RK43 (see Table 2). Starting from this explicit multirate scheme we obtain a multirate IMEX scheme by employing the diagonally implicit Runge–Kutta method listed in Table 6 as the inner method in an additional RFSMR construction stage. To examine the error behavior of this discretization we consider the solution at t = 1. As in this case no exact solution is available we obtain our reference solution by solving the system with a time step of t = 10−5 with the pure IMEX scheme. The difference of solutions obtained with larger time
1542
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
Fig. 6. Error to reference, obtained with different time steps.
steps to this reference in L 1 norm are shown in Fig. 6. The average slope of the graphs in the double logarithmic diagram affirm the order of convergence that was shown formally in Section 3. 5. Conclusions and outlook In this paper we presented and examined a generic explicit Runge–Kutta multirate scheme, named recursive flux splitting multirate (RFSMR). The created methods can be cast in partitioned Runge–Kutta form, they are internally consistent and of second order accuracy for second order base methods. Even third order convergence can be obtained if the base methods satisfy additional order conditions. A third order multirate scheme was presented and has proven superior to second order schemes even at moderate Courant numbers. A further test was dedicated to the multirate IMEX approach. Numerical results affirm that third order convergent IMEX schemes as well as multirate IMEX schemes can be constructed using the RFSMR approach. An implementation of RFSMR in a three-dimensional air pollution model developed and employed at the Leibniz Institute for Tropospheric Research showed promising results. The multirate approach was employed only on the advection operator causing only a vanishing fraction of the total computational cost in a full scale simulation. Anyway even for simulations involving complex chemistry simulations a speedup comparable to the speedup of pure advection can be obtained. This is mainly due to the local nature of the chemistry computations. The advection operator influences chemistry in the form of a source term included in the implicit integration. Every stage of the explicit Runge–Kutta stage changes this source term and thus causes a discontinuity in the right-hand side of the equation solved implicitly. For higher order implicit solvers this makes a computationally expensive restart procedure necessary. Application of a multirate scheme leads to less frequent updates of the right-hand side which leads to fewer discontinuities and finally to an increased efficiency of the implicit solving routine. A separate paper on details of this implementation is in preparation. This will also include algorithmic details and discussions on the efficient implementation and special challenges arising in the context of parallel programming. Acknowledgement This work was funded by the German Research Foundation (DFG), grant No. WO 1204/1-2. References [1] E. Constantinescu, A. Sandu, On Adaptive Mesh Refinement for Atmospheric Pollution Models, Lecture Notes in Comput. Sci., vol. 3515, Springer, Berlin/Heidelberg, 2005, pp. 798–805. [2] E. Constantinescu, A. Sandu, Multirate timestepping methods for hyperbolic conservation laws, J. Sci. Comput. 33 (3) (2007) 239–278. [3] W. Dabberdt, M. Carroll, D. Baumgardner, G. Charmichael, R. Cohen, T. Dye, J. Ellis, G. Grell, S. Grimmond, S. Hanna, J. Irwin, B. Lamb, S. Madronich, J. McQueen, J. Meagher, T. Odman, J. Pleim, H. Schmid, D. Westphal, Meteorological research needs for improved air quality forecasting, Bull. Am. Meteorol. Soc. 85 (2004) 563–586. [4] C. Gear, D. Wells, Multirate linear multistep methods, BIT Numer. Math. 24 (4) (1984) 484–502. [5] W. Hundsdorfer, B. Koren, M. van Loon, J. Verwer, A positive finite-difference advection scheme, J. Comput. Phys. 117 (1995) 35–46. [6] W. Hundsdorfer, A. Mozartova, V. Savcenco, Analysis of explicit multirate and partitioned Runge–Kutta schemes for conservation laws, CWI Report MAS-E0715, Centrum voor Wiskunde en Informatica, Amsterdam, 2007. [7] W. Hundsdorfer, J. Verwer, Numerical Solution of Time-Dependent Advection–Diffusion Reaction Equations, Springer Ser. Comput. Math., Springer, Berlin, Heidelberg, New York, 2003.
M. Schlegel et al. / Applied Numerical Mathematics 62 (2012) 1531–1543
[8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
1543
Z. Jackiewicz, R. Vermiglio, Order conditions for partitioned Runge–Kutta methods, Appl. Math. 45 (4) (1998) 301–316. S.-S. Jiang, C.-W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228. O. Knoth, R. Wolke, Implicit–explicit Runge–Kutta methods for computing atmospheric reactive flows, Appl. Numer. Math. 28 (2–4) (1998) 327–341. S. Osher, R. Sanders, Numerical approximations to nonlinear conservation laws with locally varying time and space grids, Math. Comp. 41 (1983) 321–336. H.H. Robertson, Numerical Analysis: An Introduction, Academic Press, 1966, pp. 178–182 (Chap. “The solution of a set of reaction rate equations”). M. Schlegel, O. Knoth, M. Arnold, R. Wolke, Multirate Runge–Kutta schemes for advection equations, J. Comput. Appl. Math. 226 (2009) 345–357. H. Tang, G. Warnecke, A class of high resolution schemes for hyperbolic conservation laws and convection–diffusion equations with varying time and space grids, SIAM J. Sci. Comput. 26 (4) (2005) 1415–1431. R. Wang, R. Spiteri, Linear instability of the fifth-order WENO method, SIAM J. Numer. Anal. 45 (5) (2007) 1871–1901. L. Wicker, W. Skamarock, A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta time differencing, Mon. Weather Rev. 126 (1998) 1992–1999. R. Zvan, P. Forsyth, K. Vetzal, Robust numerical methods for PDE models of Asian options, J. Comput. Finance 1 (1998) 39–78.