Abstract

A new time-splitting method for the integration of the compressible Euler equations is presented. It is based on a two-step peer method, which is a general linear method with second-order accuracy in every stage. The scheme uses a computationally very efficient forward–backward scheme for the integration of the high-frequency acoustic modes. With this splitting approach it is possible to stably integrate the compressible Euler equations without any artificial damping. The peer method is tested with the dry Euler equations and a comparison with the common split-explicit second-order three-stage Runge–Kutta method by Wicker and Skamarock shows the potential of the class of peer methods with respect to computational efficiency, stability, and accuracy.

1. Introduction

In atmospheric models the highest-frequency modes are often not the physical modes of interest. On the other hand, severe stability constrains for the numerical integrator arise from those meteorologically irrelevant modes. A common strategy to avoid this problem is a splitting approach: the differential equation is split into two parts. The slow part is integrated with one numerical method and a time step size restricted by the Courant–Friedrichs–Lewy (CFL) number of the low-frequency modes. For the integration of the high-frequency modes a simpler method is used together with smaller time steps so that these small time step sizes satisfy the CFL condition dictated by the high-frequency modes.

A widespread method was introduced by Klemp and Wilhelmson (1978) that uses a leapfrog time discretization for the slow part and the computationally very efficient forward–backward scheme (FB) by Fischer (1965) for the high-frequency modes. Skamarock and Klemp (1992) introduced the divergence damping which stabilizes split-explicit schemes by damping acoustic modes. Then Wicker and Skamarock (1998) improved the time-splitting idea by using a second-order two-stage Runge–Kutta method (RK2) for the low-frequency modes instead of the leapfrog method together with the same forward–backward scheme for the high-frequency modes and divergence damping for stabilization. They showed that RK2 is stable and as accurate as the leapfrog method but computationally more efficient. Wicker and Skamarock (2002) replaced the two-stage Runge–Kutta method by a three-stage Runge–Kutta method (RK3) in order to use this method together with higher-order spatial discretizations (RK2 does not run stably with even-order spatial discretizations) and to improve stability and accuracy. Purser (2007) showed that RK3 is second order in time for nonlinear problems. These two split-explicit Runge–Kutta methods by Wicker and Skamarock are widely used [e.g., in the Weather Research and Forecasting Model (WRF) of the National Center for Atmospheric Research (NCAR; Skamarock et al. 2005] and in the Consortium for Small-Scale Modeling (COSMO) operational model provided by the German Weather Service [Deutscher Wetterdienst (DWD)].

The disadvantage of these Runge–Kutta methods is the fact that the high-frequency modes still constrain the maximal time step size if no additional damping term such as divergence damping is used. In a linear test equation even with the use of the analytical solution for the high-frequency modes (i.e., with an infinite number of small time steps) there is the constraint csΔtx < π, where cs is the speed of sound. In contrast the CFL number of advection for RK3 is 3, which means that in case of maximal wind speeds below of 190 m s−1 the acoustic modes constrain the maximal time step size independently of the number of small time steps per large time step. One approach to avoid this problem if one does not want to use an additional damping term is to use another fast integrator. We found out that using a method introduced by Shchepetkin and McWilliams (2005) removes the stability constraints caused by the acoustic modes if the number of small time steps per large time step is sufficient large. On the other hand this fast integrator is at least twice as expensive as the forward–backward integrator and the scheme is only stable for large time steps (there is a region of medium time steps for which the scheme becomes unstable). Another approach was done by Kar (2006) who combined RK3 together with the trapezoidal scheme for the fast modes instead of the forward–backward scheme and used the same time step sizes for both parts what is possible because the trapezoidal scheme is implicit. Tests with shallow-water equations by Kar (2006) showed that this approach was computationally more efficient than RK3 with FB and only little less accurate. But the use of an implicit integrator requires the solution of three elliptic equations per time step, which makes the method much harder to implement.

Our aim is the construction of a split-explicit method, which is stable without the use of any artificial damping. Furthermore, we want to use the computationally very efficient forward-backward scheme FB as integrator for the high-frequency modes like Klemp, Wilhelmson, Wicker, and Skamarock did before. For this purpose we consider a subclass of general linear methods. General linear methods were introduced by Burrage and Butcher (1980), a good overview can be found in Hairer et al. (1993). Peer methods are a very comprehensive class of general linear methods; they include the common Runge–Kutta methods and linear multistep methods. They can be interpreted as cyclic multistage multistep methods, which means that each of the stages of a peer method is a linear multistep method. The new feature of peer methods is that they possess several stages like Runge–Kutta methods, but all these stages have the same properties, no extraordinary solution variable is used. These methods combine positive features of both Runge–Kutta and linear multistep methods. Because of the generality of this class of methods they showed to be applicable in many different kinds of problems. Linear-implicit peer methods were introduced by Weiner et al. (2004). Because of the same order in every stage they have no order reduction even for very stiff systems and are implemented in the KARDOS code (Gerisch et al. 2008). We consider explicit peer methods that were successfully used by Weiner et al. (2008) and Schmitt et al. (2008) among others. With explicit peer methods we want to construct a method that is as computationally efficient and accurate as RK3 and the maximal time step size should not be restricted by the acoustic modes. In contrast to RK3 all these aims should be properties of the pure split-explicit peer method without the addition of any artificial damping term.

This paper is organized as follows: we present a new view on split-explicit Runge–Kutta methods together with a von Neumann stability analysis for a linearized test equation in section 2. In section 3 we define the class of peer methods, apply them to split differential equations, and present stability results and order conditions together with the parameters of a peer method. Tests with the compressible Euler equations are the content of section 4. Finally, we give some conclusions and an outlook in section 5.

2. The time-splitting idea for Runge–Kutta methods

We consider the numerical solution of autonomous split differential equations

 
formula

in [t0, te] where f represents the slow part and g represents the fast part.

As a preliminary step to the introduction of split-explicit peer methods we consider split-explicit Runge–Kutta methods in a more technical way than Wicker and Skamarock did because we do not restrict ourselves to one integrator for the fast part of (1) but our formulation of the split-explicit scheme contains the fast part as an initial value problem. With this point of view we can focus on the properties of the main solver for the slow part (i.e., on the underlying Runge–Kutta and peer methods). Furthermore, this approach will simplify the derivation of order conditions for split-explicit peer methods.

To solve split differential equation (1) we propose a scheme where an explicit Runge–Kutta method is used for the slow part f of the right-hand side while the solution of the ordinary differential equation (ODE) with a constant slow part, that is, = c + g(y), is defined implicitly by a differential equation so that in the absence of g the explicit Runge–Kutta method is recovered. For simplicity of notation we will denote the weights bi as the s + 1st row of 𝗔 from the Butcher tableau; for example, for RK3 the Butcher tableau is

 
formula

and therefore we have aij = 0 except a21 = ⅓, a32 = ½, and a43 = 1. The nodes c1 = 0, c2 = ⅓, c3 = ½, and c4 = 1 are the points on the time axis where the numerical solution approximates the analytical solution and Δt is the length of the integration interval. The ith stage (for 2 ≤ is + 1) of such a method is the following:

 
formula
 
formula
 
formula

The considered split-explicit Runge–Kutta methods of Wicker and Skamarock (1998) and Wicker and Skamarock (2002) fit into this formulation if the FB scheme is used for the integration of (2). Another approach was performed by Knoth and Wolke (1998) where the integration of (2) does not start at but at for the ith stage. This approach has the advantage that the sum of the integration intervals is not Σci but only 1 (in case of nonnegative ascending nodes). Therefore, it is computationally more efficient, but on the other hand it has worse stability properties. However, we will focus on the generalization of the Wicker and Skamarock approach in the remainder of this paper.

With this splitting approach we now want to discuss stability properties. Therefore, we consider a linearized version of the one-dimensional Euler equations. Because we will later consider Runge–Kutta methods together with divergence damping we include such a damping term in our linear test equation:

 
formula

This linear test equation describes the essential difficulties that arise when integrating the one-dimensional Euler equations. The variable u can be interpreted as a wind speed and p as a pressure. Here U is the constant wind speed to obtain a linear equation and cs is the speed of sound and therefore csU. The term νuxx is the divergence damping as presented by Skamarock and Klemp (1992), which becomes the second derivative of the wind in the one-dimensional case. With this linear equation we perform our stability analysis: we discretize them in space and time and use the common von Neumann stability idea to discuss the resulting Fourier modes separately.

A detailed description of this procedure can be found in Wicker and Skamarock (1998) for a staggered grid with a third-order upwind scheme. We use the same schemes as Wicker and Skamarock (2002) for spatial discretization, including the third- and fifth-order upwind schemes. After spatial discretization we obtain a split ODE

 
formula

where 𝗟 is an antidiagonal matrix (for ν = 0), which comes from the fast part of the linear test equation and 𝗡 is a diagonal matrix that represents the advection part. While the eigenvalues of the fast part 𝗟 are purely imaginary in the absence of divergence damping and undergo a little shift into the left half-plane if divergence damping is used, the slow part 𝗡 has eigenvalues that lie in the left half-plane with comparable imaginary and real parts. The left panel of Fig. 1 shows these eigenvalues for all wavelengths.

Fig. 1.

(left) Eigenvalues of the slow part of the spatial discretized linear test equation (for third-order upwind in black and for fifth-order upwind in gray) and (right) maximum amplification factor for RK3 only for this slow part and third-order upwind. Unstable regions with an eigenvalue > 1 are in gray. The contour interval is 0.05 (in the stable regions).

Fig. 1.

(left) Eigenvalues of the slow part of the spatial discretized linear test equation (for third-order upwind in black and for fifth-order upwind in gray) and (right) maximum amplification factor for RK3 only for this slow part and third-order upwind. Unstable regions with an eigenvalue > 1 are in gray. The contour interval is 0.05 (in the stable regions).

At first we consider stability results for the pure advection equation (i.e., 𝗟 = 0). We apply the second-order RK3 method introduced by Wicker and Skamarock (2002) to the linear advection equation. The right panel of Fig. 1 shows the maximum amplification factor for different wavenumbers 2π/k and Courant numbers UΔtx. Obviously instabilities occur for the high-frequency modes.

We now apply RK3 to the full test equation (5) with a time step size Δt. The implicitly defined part in (2) is integrated analytically and requires the solution of a linear ODE. For given CFL numbers UΔtx and csΔtx we apply RK3 to (5) and compute the stability matrix and its maximal absolute eigenvalues for wavelengths k between 2Δx and 20Δx. Figure 2 shows the contours of the maximal absolute eigenvalues in dependency of Cadv = UΔtx (vertical axis) and Csound = csΔtx (horizontal axis). The thick line has a slope of 1/12 and only the stability region below of this line is relevant because we suggest U < cs/12 in practical applications as Wicker and Skamarock (1998) did. Figure 2 has to be interpreted in that way that, for a given point on the horizontal axis, the stability between this point and the corresponding point on the thick line above is relevant, because for these coordinates of CFL numbers the wind speed U varies between 0 and cs/12 and all these wind speeds can occur in practical applications. Another time step size Δt can be considered by starting from another point of the horizontal axis. So the diagrams are valid for a fixed value of cs.

Fig. 2.

Maximum amplification factor for the linear test equation with third-order upwind spatial discretization for RK3 (left) without and (right) with divergence damping. Unstable regions with an eigenvalue > 1 are in gray. The contour interval is 0.1 (in the unstable regions).

Fig. 2.

Maximum amplification factor for the linear test equation with third-order upwind spatial discretization for RK3 (left) without and (right) with divergence damping. Unstable regions with an eigenvalue > 1 are in gray. The contour interval is 0.1 (in the unstable regions).

As we can observe RK3 is stable in the whole region of interest only if it uses divergence damping. We used a fixed value of 0.025 = νΔτx2, the same as Wicker and Skamarock (1998). This means that the original differential equation is artificially altered by a term νuxx, which depends on the actual spatial and temporal discretization and increases with decreasing time step size. Such a method can never converge in time if the spatial step size is kept constant. Without this damping term, RK3 is unstable in the whole region where Csound > π and therefore the acoustic modes restrict the maximal time step size because csU. These results are the main motivation for the construction of time-splitting methods, which do not suffer from stability restrictions due to acoustic modes, but on the other hand, do not need any artificial damping term. This aim will be reached by split-explicit peer methods.

3. Order and stability investigations for peer methods

Explicit two-step peer methods for first-order differential equations = f (y) were considered by Biermann (2005) and Weiner et al. (2008), and in parallel form also by Weiner et al. (2006) and Schmitt et al. (2008). In this paper we also imply values of the numerical solution from the current time step. For constant step sizes these methods are given by

 
formula

with the following notations:

 
formula

and 𝗔, 𝗕, 𝗥, 𝗦 ∈ ℝs×s, where 𝗥 and 𝗦 are strictly lower triangular matrices to obtain explicit methods.

For the integration of = f (y) the use of the numerical solution from the actual time level (i.e., 𝗦 ≠ 0) does not lead to a more general class of peer methods than the case 𝗦 = 0 because every peer method with 𝗦 ≠ 0 can be generated by a peer method with 𝗦 = 0. This can be seen by subtracting 𝗦𝗬m from (6) and then multiplying both sides with (𝗜 − 𝗦)−1. Therefore, in former applications of peer methods 𝗦 = 0 was used. However, because of the following splitting approach the use of 𝗦 ≠ 0 leads to a more general class of split-explicit peer methods and therefore we use these additional degrees of freedom.

To solve the split differential equation in (1) we propose a scheme where the explicit peer method is used for the slow part f, while the solution of = c + g(y) is defined implicitly by a differential equation so that in the absence of g the explicit peer method is recovered. A further generalization is done by allowing the lengths of the integration intervals αi to be independent of the nodes ci in contrast to split-explicit Runge–Kutta methods where αi = ci. The ith stage (for 1 ≤ is) of such a method reads as the following:

 
formula

with initial value

 
formula

and define

 
formula

In case of g = 0 such a method is obviously an explicit peer method with 𝗔 = (aij), 𝗕 = (bij), 𝗥 = (rij), and 𝗦 = (sij).

Every split-explicit Runge–Kutta method can be written as a split-explicit peer method. For example RK3 reads

 
formula

as a peer method.

The split-explicit scheme with the leapfrog method for the slow part, used by Skamarock and Klemp (1992) and Skamarock and Klemp (1994), can also be written as a peer method with

 
formula

We now derive order conditions for split-explicit peer methods. Consistency of split-explicit peer methods is discussed by considering the local residuals obtained by substituting the analytical solution into the method. We are interested only in the error of the slow part of the split differential equation in (1) and therefore we assume that we can solve (7) exactly. Then the residuals are defined by

 
formula

We define the order of consistency p of a method in (9) by the order of the local residuals, that is, the method has order of consistency p if

 
formula

holds.

Because a peer method has the same order in every stage the derivation of order conditions is simple and similar to the derivation of order conditions for linear multistep methods (i.e., in practice this is done by Taylor expansions as shown in the appendix). We have the following result.

A method (9) has order of consistency p ≤ 2 if

 
formula

with

 
formula

for i = 1, … , s.

As the proof in the appendix shows it is not possible to construct stable split-explicit peer methods with the approach in (8) with order of consistency of 3 because of the order condition α = 0. This comes from the absence of function evaluations in (8), the inclusion of function evaluations allows to solve the order conditions for the order of 3 with α ≠ 0 and therefore the construction of peer methods with the order of 3 and higher. However, we are interested in methods with the order of 2 and therefore we consider the simpler case where no function evaluations are used as initial values for the fast ODE in (7) in the remainder.

The order conditions 𝗔𝗕(k) = 0 are the traditional order condition for peer methods while the conditions are new order conditions that result from the splitting approach. We now want to solve the order conditions to obtain split-explicit peer methods that are of order p = s − 1. For this aim we use the coefficient matrix 𝗔 to fulfill the traditional order conditions as was done by Biermann (2005), Weiner et al. (2006), Weiner et al. (2008), and Schmitt et al. (2008). Because we only want methods of order s − 1, 𝗔 has s degrees of freedom remaining, which are expressed with a vector β for simplicity of notation. The matrix 𝗕 is used to fulfill the consistency condition , where and the s − 1 additional split order conditions. Therefore, 𝗕 is fully determined by the order conditions in contrast to former applications of peer methods where the coefficients of 𝗕 were available for stability optimizations or superconvergence conditions. With the notations

 
formula

one obtains methods with order of consistency p = s − 1 ≤ 2 by

 
formula

This choice of parameters is possible if the nodes are distinct from each other because then the Vandermonde matrix 𝗩1 is regular.

After considering order conditions we now want to discuss stability properties. Therefore we apply the split-explicit peer methods to the linear test equation in (5) considered in section 2. Figure 3 shows the results for a peer method (with its coefficients given at the end of this section) with third-order and fifth-order upwind spatial discretization.

Fig. 3.

Maximum amplification factor for (top) the linear advection equation and (bottom) the full linear test equation for the peer method with (left) third-order and (right) fifth-order upwind spatial discretization. Unstable regions with an eigenvalue > 1 are in gray. The contour interval at the top is 0.05 and bottom 0.1.

Fig. 3.

Maximum amplification factor for (top) the linear advection equation and (bottom) the full linear test equation for the peer method with (left) third-order and (right) fifth-order upwind spatial discretization. Unstable regions with an eigenvalue > 1 are in gray. The contour interval at the top is 0.05 and bottom 0.1.

For the pure advection equation (Fig. 3), the peer method has stability properties comparable to RK3. The resulting CFL numbers for the peer method are approximately 1.7 for third-order and 1.5 for fifth-order upwind discretization. For the full test equation the peer method is stable in the whole region of interest (the instabilities above the thick line correspond to wind speeds that do not occur in practical applications) and therefore it fulfills the aim that the acoustic modes should not restrict the maximal time step size. Remember that this is possible without the introduction of divergence damping.

We now describe how we got the coefficients of the considered peer method. Our aim was to construct a split-explicit peer method in (9) with s = 3 stages, order p = s − 1 = 2, and stability properties comparable to RK3. For these aims we have the coefficient matrices 𝗔, 𝗕, 𝗥, and 𝗦 and the vectors c and α. As described above we use s(s − 1) coefficients of 𝗔 to satisfy the classical peer order conditions and the remaining s degrees of freedom are expressed with a vector β. The additional split order conditions are fulfilled by 𝗕. Furthermore, we define cs = 1 so that the last stage is an approximation to the analytical solution at time points tm. Therefore we have s(s + 2) − 1 = 14 remaining parameters to find a good method.

For this optimization we use a Monte Carlo search strategy in this 14-dimensional subspace of peer methods. Our main goal is the construction of a method that has at least the stability properties of RK3. So during the first phase of the optimization process the optimization function consists of the Courant numbers for third-order and fifth-order upwind discretization and of a measure for the stability region below of the line with the slope ratio 1/12. We compute the stability regions by using the full linear test equation where the integration of the fast part is not done analytically but numerically with the FB scheme and 10 steps per large time step. If a peer method with Courant numbers of at least 1.7 (1.4) is found and the method is stable up to a value of 9 on the horizontal axis, which means that the maximal amplification factor is smaller than 1 in 90% of the region of interest, the optimization process enters its second phase. Now the optimization function contains the error constants and ‖α1 while the achieved stability properties remain as side conditions.

Finally we derived the method presented in this article. Its coefficients are

 
formula

It has the following stability properties where the indices of the Courant numbers denote the spatial order (the corresponding values for RK3 are in parentheses for comparison):

 
formula

The expense for the integration of the fast part of the split differential equation depends on the sum of the integration intervals αi. For the peer method ‖α1 = 1.44, while for RK3 c1 + c2 + c3 = 11/6, which is 27% more. The moduli of the eigenvalues of (𝗜 − 𝗦)−1𝗕 must not exceed 1 for zero stability of the peer method, they are

 
formula

4. Numerical tests

So far we presented a peer method with its properties and stability diagrams and compared it with RK3. Now we want to compare both methods in a practical application. We consider the two-dimensional dry Euler equations in conservative flux form (with divergence damping):

 
formula
 
formula
 
formula
 
formula

Here u and w are the horizontal and vertical winds, ρ is the density, p is the pressure, θ is the potential temperature, and g is the acceleration of gravity. The prognostic variables are ρ, ρu, ρw, and ρθ. The pressure p is given diagnostically by the equation of state:

 
formula

where R is the gas constant for dry air, κ = R/cp, cp is the heat capacity of dry air at constant pressure, and p0 is the pressure at ground.

We use a finite-volume spatial discretization on an Arakawa C grid, so the winds are defined on the cell edges while all scalar variables are defined at the cell centers as illustrated by Fig. 4. The time splitting is done in the manner that the large time step Δt is split into ns equal small time steps Δτ = Δt/ns. Then the nonbold terms are evaluated in every small time step while the bold terms are kept constant. So the nonbold derivatives can simply be evaluated with one difference per derivative. This corresponds to the integration of (7) and is done with the FB scheme. As first step ρu and ρw are updated from (11) and (12) simultaneously with values ρ and p from the last small time step and frozen bold derivatives. Then ρ and ρθ are computed with these updated winds from (10) and (13). Finally p is evaluated from ρθ by the equation of state in (14). For the computation of the bold derivatives we need to compute the fluxes with the third- or fifth-order upwind scheme. With one of these schemes we compute the winds at the centers and vertices and the potential temperature on the edges. The fluxes in (13) are computed in the way that θ = ρθ/ρ is computed on the cell edges with the upwind scheme and updated by the underlying method (e.g., Σaijθj for a split-explicit Runge–Kutta method), and then kept frozen for the ns small time steps. Obviously the evaluation of the bold fluxes is much more expensive because of the expensive high-order upwind-biased schemes.

Fig. 4.

The positions of the variables on the Arakawa C grid.

Fig. 4.

The positions of the variables on the Arakawa C grid.

We present the results for two test cases that were used by Wicker and Skamarock (1998) and Wicker and Skamarock (2002). The first test is the dry rising bubble. As an initial condition we have an adiabatic atmosphere with the exception of an initial thermal of 2°C with radius and height of 2 km; that is, the potential temperature is perturbed by

 
formula

where

 
formula

xc = 0 km and zc = xr = zr = 2 km as described by Bryan and Fritsch (2002). We use a 20 km × 10 km domain with a spatial resolution of 125 m and the third-order upwind scheme. A uniform horizontal flow of 20 m s−1 (from the left) and periodic boundary conditions cause a lateral transport of the bubble. After the integration period of 1000 s the bubble should be located in the center of the domain and remain symmetric. Figure 5 shows the contours of the vertical wind and the potential temperature after 1000 s. The integration was performed with a time step size of 2 s and ns = 10. The maximal (horizontal) wind speed that occurred was approximately 29 m s−1, so that this time step size corresponds to a Courant number of 0.46. Both numerical solutions are in agreement on the main lines with the exact solution (exact with respect to time). They show no numerical oscillations and are nearly perfect symmetric, only the vertical winds show tiny asymmetries, which come from the horizontal flow and disappear if no background wind is used. So both methods, RK3 and the peer method, are able to integrate the compressible Euler equations stably and approve the theoretical results from the linear stability theory in practice. But RK3 has to use divergence damping (with the damping coefficient 0.025 as mentioned at the end of section 2) while this is not necessary for the peer method.

Fig. 5.

(top) Vertical wind and (bottom) potential temperature after 1000 s for the rising bubble problem computed with (left) RK3 and divergence damping and (right) the peer method without divergence damping, both tests with a time step size of 2 s, ns = 10, and with third-order upwind spatial discretization. The contour interval at the top is 1.5 m s−1 and bottom 0.25 K.

Fig. 5.

(top) Vertical wind and (bottom) potential temperature after 1000 s for the rising bubble problem computed with (left) RK3 and divergence damping and (right) the peer method without divergence damping, both tests with a time step size of 2 s, ns = 10, and with third-order upwind spatial discretization. The contour interval at the top is 1.5 m s−1 and bottom 0.25 K.

Furthermore, the peer method is also appropriate to low Mach number problems.: It however runs stably even for very high numbers of ns (e.g., ns = 30). In practice such a method is useful if the speed of sound dominates the wind speed and therefore the CFL condition for the fast waves is much more restrictive than the pure advective CFL condition. For example, in this test case the maximal horizontal wind speed is 29 m s−1; this is more than one magnitude smaller than the speed of sound. The resulting CFL condition of the peer method with the third-order upwind scheme for pure advection is

 
formula

and therefore the maximal time step size (for a spatial resolution of 125 m) is

 
formula

But computations with a time step size of 7.4 s result in a Courant number for the sound waves of

 
formula

where we assume 340 m s−1 for the speed of sound and the factor 2 comes from two-dimensionality of the considered test (i.e., from diagonally propagating sound waves). Therefore, the number of small time steps must be ns = 29 at least when using FB as fast integrator. We tested the peer method with ns = 30 and a time step size of 7 s what corresponds to a Courant number of 1.6. Despite this large time step the solution (Fig. 6) looks quite reasonable. It shows no oscillations and is nearly symmetric.

Fig. 6.

(left) Vertical wind and (right) potential temperature after 1000 s for the rising bubble problem computed with the peer method without divergence damping and with a time step size of 7 s, ns = 30, and third-order upwind spatial discretization. The contour interval on the left is 1.5 m s−1 and on the right 0.25 K.

Fig. 6.

(left) Vertical wind and (right) potential temperature after 1000 s for the rising bubble problem computed with the peer method without divergence damping and with a time step size of 7 s, ns = 30, and third-order upwind spatial discretization. The contour interval on the left is 1.5 m s−1 and on the right 0.25 K.

The second test problem is related to the first one. Again there is a balanced initial atmosphere, but for the density current test the bubble is colder than the surrounding air. The temperature is perturbed by

 
formula

with L from (15), but with zc = 3 km and xr = 4 km. Because the bubble is colder, it sinks and after crashing the bottom (with free-slip vertical boundary conditions) several large eddies form. The domain has a width of 36 km and the uniform horizontal flow of 20 m s−1 (from the right) should provide a symmetric solution after 900 s (half orbit). The horizontal resolution is 200 m and we use the fifth-order upwind scheme for the spatial discretization. This test was performed with a time step size of 2 s which corresponds to a Courant number of 0.5 because of the strong winds up to 50 m s−1. The fast part is integrated with 1/3 s (i.e., with ns = 6). Again only the Runge–Kutta method uses divergence damping. Both solutions show three eddies at each side where the backward directed (left) eddies are better pronounced (see Fig. 7). This phenomenon was also documented by Wicker and Skamarock (2002). Nevertheless, the solutions differ slightly from that by Wicker and Skamarock (2002). This comes from two reasons: on one hand the spatial order drops at the boundaries (up to order 1 for the pure upwind scheme at the outermost grid cell) in order to use smaller stencils there. A more sophisticated approach might produce a better solution. On the other hand, we did not add diffusion in contrast to Wicker and Skamarock (2002) because we want to see the results for the pure Euler equations. In fact, there is no solution for the pure Euler equations for this problem because discontinuities form in a finite time. The dissipation in the advection scheme produces a viscous solution for the inviscid Euler equations. So there is damping in this formulation for both methods and in contrast to the first test the damping from the advection scheme has an important influence on the solution. In any case this test shows again that the differences between the results of the two methods are marginal, the peer method without divergence damping is as stable as RK3 with divergence damping even when high wind speeds up to 50 m s−1 occur.

Fig. 7.

Perturbation after 900 s for the density current problem computed with (top) RK3 and divergence damping and (bottom) the peer method without divergence damping, both with a time step size of 2 s, ns = 6, and fifth-order upwind spatial discretization. The contour interval is 0.5 K.

Fig. 7.

Perturbation after 900 s for the density current problem computed with (top) RK3 and divergence damping and (bottom) the peer method without divergence damping, both with a time step size of 2 s, ns = 6, and fifth-order upwind spatial discretization. The contour interval is 0.5 K.

5. Conclusions

We presented a new methodology to describe time-splitting methods. The basic principle of the presented approach is the assumption that one part, the fast part, of the split differential equation can be solved analytically so that stability and order investigations can be made for the underlying method, which solves the slow part. With this methodology we considered split-explicit Runge–Kutta methods and reproduced the known stability results.

Then we introduced the class of peer methods and used them as underlying method for the solution of split differential equations. We presented order conditions and stability results and finally derived a split-explicit peer method, which is as accurate and efficient as the common split-explicit Runge–Kutta method, RK3. Particularly linear stability analysis showed that the peer method even without any artificial damping is as stable as RK3 with divergence damping.

The computational effort of the peer method is comparable to RK3: both have s = 3 stages and therefore need three evaluations of the expensive, slow part of the right-hand side per large time step. Because the sum of the fast integration intervals for the peer method ‖α1 = 1.44 is about 20% smaller than the sum of the nodes of RK3 (1/3 + 1/2 + 1 = 11/6) the expense for the evaluation of the fast part is a little bit lower for the peer method. In general peer methods have a higher overhead than Runge–Kutta methods because they use more linear combinations of the numerical solution and its function evaluations but for PDEs this fact is a negligible disadvantage because of the very expensive right-hand side.

One disadvantage of the peer method might be the memory requirement. Because only the secondary diagonal of the Butcher tableau of RK3 has nonzero entries RK3 needs the storage capacity for three numerical solutions and function evaluations: one for the initial value, another for the function evaluation, and the third for the updated fast part. In contrast, the peer method needs twice of this storage: s = 3 for the linear combinations eiT𝗕Ym−1 + eiT𝗦Ym, which are successive updated to Ymi, that is, in the ith stage the integration of the fast part starts with the initial value eiT𝗕Ym−1 + eiT𝗦Ym and every update is stored at the same place so that in the end eiT𝗕Ym−1 + eiT𝗦Ym is overwritten by Ymi, the same holds for the function evaluations so that the storage for six solutions is needed at all. Because in practical applications problems are solved in a massively parallel environment that usually uses domain decomposition techniques every CPU has to know only the values on his small part of the whole domain and therefore the requirement for twice of the memory capacity is no disadvantage, there are enough reserves.

The application of the split-explicit methods to the compressible Euler equations in conservative flux form with a finite-volume spatial discretization confirmed the linear stability results: both methods integrate the two considered test problems stably. So we found a split-explicit method for the elastic equations that is as accurate, stable, and efficient as RK3 but without the need of damping. The peer method is stable even for a very high number of small time steps and therefore it is appropriate for solving low Mach number problems. The use of a time step size that nearly corresponds to the Courant number from linear stability theory still produces an acceptable solution.

Important considerations for further research are the linear-implicit peer methods. They were presented by Podhaisky et al. (2005) among others. For the order conditions no exact Jacobian is needed and because of the full order in every stage no order reduction occurs. In the considered tests they showed their potential compared to established Rosenbrock-type solvers. It would be interesting to know whether they can be an attractive alternative to existing implicit solvers of the compressible Euler equations with respect to stability, accuracy, and efficiency.

Acknowledgments

This work is a contribution to the project MetStröm (Skalenübergreifende Modellierung in der Strömungsmechanik und Meteorologie) funded by the German Research Foundation (DFG). The authors thank the referees for their valuable comments, which helped to significantly improve the presentation of the paper.

REFERENCES

REFERENCES
Biermann
,
K.
,
2005
:
Explizite Zweischrittmethoden mit Peer-Variablen. M.S. thesis, Institute of Mathematics, University of Halle, 55 pp. [Available from Martin-Luther-Universität Halle-Wittenberg, Naturwissenschaftliche Fakultät III, Institut für Mathematik, 06099 Halle, Germany.]
.
Bryan
,
G. H.
, and
J. M.
Fritsch
,
2002
:
A benchmark simulation for moist nonhydrostatic numerical models.
Mon. Wea. Rev.
,
130
,
2917
2928
.
Burrage
,
K.
, and
J. C.
Butcher
,
1980
:
Non-linear stability of a general class of differential equation methods.
BIT
,
20
,
185
203
.
Fischer
,
G.
,
1965
:
Comments on “Some problems involved in the numerical solutions of tidal hydraulics equations.”.
Mon. Wea. Rev.
,
93
,
110
111
.
Gerisch
,
A.
,
J.
Lang
,
H.
Podhaisky
, and
R.
Weiner
,
2008
:
High-order linearly implicit two-step peer-finite element methods for time-dependent PDEs.
Appl. Numer. Math.
,
59
,
624
638
.
Hairer
,
E.
,
S. P.
Nørsett
, and
G.
Wanner
,
1993
:
Solving Ordinary Differential Equations I.
2nd ed. Springer-Verlag, 528 pp
.
Kar
,
S. K.
,
2006
:
A semi-implicit Runge–Kutta time-difference scheme for the two-dimensional shallow-water equations.
Mon. Wea. Rev.
,
134
,
2916
2926
.
Klemp
,
J. B.
, and
R. B.
Wilhelmson
,
1978
:
The simulation of three-dimensional convective storm dynamics.
J. Atmos. Sci.
,
35
,
1070
1096
.
Knoth
,
O.
, and
R.
Wolke
,
1998
:
Implicit-explicit Runge–Kutta methods for computing atmospheric reactive flows.
Appl. Numer. Math.
,
28
,
327
341
.
Podhaisky
,
H.
,
R.
Weiner
, and
B. A.
Schmitt
,
2005
:
Rosenbrock-type ‘Peer’ two-step methods.
Appl. Numer. Math.
,
53
,
409
420
.
Purser
,
R. J.
,
2007
:
Accuracy considerations of time-splitting methods for models using two-time-level schemes.
Mon. Wea. Rev.
,
135
,
1158
1164
.
Schmitt
,
B. A.
,
R.
Weiner
, and
S.
Jebens
,
2008
:
Parameter optimization for explicit parallel peer two-step methods.
Appl. Numer. Math.
,
59
,
769
782
.
Shchepetkin
,
A. F.
, and
J. C.
McWilliams
,
2005
:
The regional oceanic modeling system (ROMS): A split-explicit, free-surface, topography-following-coordinate oceanic model.
Ocean Modell.
,
9
,
347
404
.
Skamarock
,
W. C.
, and
J. B.
Klemp
,
1992
:
The stability of time-split numerical methods for the hydrostatic and the nonhydrostatic elastic equations.
Mon. Wea. Rev.
,
120
,
2109
2127
.
Skamarock
,
W. C.
, and
J. B.
Klemp
,
1994
:
Efficiency and accuracy of the Klemp–Wilhelmson time-splitting technique.
Mon. Wea. Rev.
,
122
,
2623
2630
.
Skamarock
,
W. C.
,
J.
Dudhia
,
D.
Gill
,
D.
Barker
,
W.
Wei
, and
J.
Powers
,
2005
:
A description of the advanced research WRF version 2.
NCAR Tech. Note NCAR/TN-468+STR, National Center for Atmospheric Research, Boulder, CO, 88 pp
.
Weiner
,
R.
,
B. A.
Schmitt
, and
H.
Podhaisky
,
2004
:
Parallel ‘peer’ two-step W-methods and their application to MOL-systems.
Appl. Numer. Math.
,
48
,
425
439
.
Weiner
,
R.
,
S.
Jebens
,
B. A.
Schmitt
, and
H.
Podhaisky
,
2006
:
Explicit parallel two-step peer methods.
Reports of the Institute of Numerical Mathematics, Rep. 10, 18 pp. [Available online at Martin-Luther-Universität Halle-Wittenberg, Naturwissenschaftliche Fakultät III, Institut für Mathematik, 06099 Halle, Germany.]
.
Weiner
,
R.
,
B. A.
Schmitt
,
H.
Podhaisky
, and
S.
Jebens
,
2008
:
Superconvergent explicit two-step peer methods.
J. Comput. Appl. Math.
,
223
,
753
764
.
doi:10.1016/j.cam.2008.02.014
.
Wicker
,
L. J.
, and
W. C.
Skamarock
,
1998
:
A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta time differencing.
Mon. Wea. Rev.
,
126
,
1992
1999
.
Wicker
,
L. J.
, and
W. C.
Skamarock
,
2002
:
Time-splitting methods for elastic models using forward time schemes.
Mon. Wea. Rev.
,
130
,
2088
2097
.

APPENDIX

Proof of Order of Consistency

For simplicity of notation we consider only values from the actual time step (i.e., rij and sij). The inclusion of values from other time levels is straightforward. For Zi(t) we obtain the following expansion:

 
formula

Therefore the residuals Δi = y(tmi) − Zi(αiΔt) satisfy the following:

 
formula

By computation one can see that the terms up to order p ≤ 2 vanish if the predetermined order conditions are fulfilled. But the condition at (αiΔt3/6)g2[y(tm)](tm) together with the order conditions for the order of 2 lead to the conditions αi = 0.

Footnotes

Corresponding author address: Stefan Jebens, Leibniz Institute for Tropospheric Research, Permoserstr. 15, D-04318 Leipzig, Germany. Email: jebens@tropos.de