We now consider the amplitude and phase properties of the split-explicit and linearly implicit methods. The analytical solution of
˙
y=M y with M from (21) is
y(tm) = exp(∆tM)y(tm−1).
When we apply the methods to the test equation the amplification matrix A advances the solution to the next time level, i.e.
Ym =AYm−1.
If there is no diffusion M has one eigenvalue which results from the advection part and two eigenvalues incorporating advection and acoustics. Let λ be the eigenvalue of exp(∆tM) which originates from the pure advection part. Then |λ| is the amplification factor and
arctanImλReλ Imλ
is the relative phase speed. Previously we only considered the maximum of the moduli of the eigenvalues of the amplification matrices A in order to obtain the regions where the corresponding methods are stable. Each of these amplification matrices has nine eigenvalues because we consider three-stage methods which are applied to 3×3-systems.
One of these nine eigenvalues corresponds to the advection part which is advanced one time step, i.e. it corresponds to λ of exp(∆tM). For this eigenvalue amplitude and rela-tive phase speed are analogously defined as for λ. These values are compared with the amplitude and relative phase speed of the analytical solution.
0 0.5 1 1.5 2 0
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Amplitude for 4∆x Wave
Cadv
Amplitude
RK3 explPeer ROS3Pw implPeer2 implPeer3
0 0.5 1 1.5 2
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Amplitude for 4∆x Wave
Cadv
Amplitude
ROS3Pw implPeer2 implPeer3
0 0.5 1 1.5 2
0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6
Amplitude for 7∆x Wave
Cadv
Amplitude
RK3 explPeer ROS3Pw implPeer2 implPeer3
0 0.5 1 1.5 2
0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6
Amplitude for 7∆x Wave
Cadv
Amplitude
ROS3Pw implPeer2 implPeer3
0 0.5 1 1.5 2
0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6
Amplitude for 10∆x Wave
Cadv
Amplitude
RK3 explPeer ROS3Pw implPeer2 implPeer3
0 0.5 1 1.5 2
0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6
Amplitude for 10∆x Wave
Cadv
Amplitude
ROS3Pw implPeer2 implPeer3
Figure 6: Amplitude for the 4∆x wave (top), the 7∆x wave (middle) and the 10∆x wave (bottom) for split-explicit methods and linearly implicit methods with partial Jacobian (left) and for linearly implicit methods with simplified Jacobian (right). The analytic amplitude in black.
0 0.5 1 1.5 2 0.5
1 1.5 2 2.5
Relative phase speed for 4∆x Wave
Cadv
Relative phase speed
RK3 explPeer ROS3Pw implPeer2 implPeer3
0 0.5 1 1.5 2
0.5 1 1.5 2 2.5
Relative phase speed for 4∆x Wave
Cadv
Relative phase speed
ROS3Pw implPeer2 implPeer3
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5
Relative phase speed for 7∆x Wave
Cadv
Relative phase speed
RK3 explPeer ROS3Pw implPeer2 implPeer3
0 0.5 1 1.5 2
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5
Relative phase speed for 7∆x Wave
Cadv
Relative phase speed
ROS3Pw implPeer2 implPeer3
0 0.5 1 1.5 2
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2
Relative phase speed for 10∆x Wave
Cadv
Relative phase speed
RK3 explPeer ROS3Pw implPeer2 implPeer3
0 0.5 1 1.5 2
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2
Relative phase speed for 10∆x Wave
Cadv
Relative phase speed
ROS3Pw implPeer2 implPeer3
Figure 7:Relative phase speed for the 4∆xwave (top), the 7∆x wave (middle) and the 10∆x wave (bottom) for split-explicit methods and linearly implicit methods with partial Jacobian (left) and for linearly implicit methods with simplified Jacobian (right). The analytic relative phase speed in black.
Figure 6 shows the amplitudes for the considered methods and three wave num-bers, Figure 7 shows the relative phase speeds. Because the split-explicit peer method explPeer is the underlying method for implPeer2 and we consider the amplitude and phase properties of the advection part only, implPeer2 with the partial Jacobian has the same amplitude and phase errors as explPeer, i.e. it adopts the good amplitude and phase properties of explPeer. We can see that the split-explicit methods and the linearly implicit peer methods with partial Jacobian have quite similar amplitude and phase properties while ROS3Pw with partial Jacobian has worse properties and becomes un-stable for smaller CFL numbers in comparison to the other methods. When the simplified Jacobian is used for the linearly implicit methods implPeer2 has the best properties. im-plPeer3 has worse amplitude properties than implPeer2 while ROS3Pw has larger phase errors. We will see the effects of these properties in Section 6.4. We can also see that the differences between the methods decrease if waves with larger wavelengths are considered because for them the eigenvalues are closer to 0.
5 Implementation
5.1 The compressible Euler equations
We consider the two-dimensional dry compressible Euler equations in conservative form with diffusion:
∂ρ
∂t =−∂ρu
∂x − ∂ρw
∂z , (24)
∂ρu
∂t =−∂ρuu
∂x − ∂ρwu
∂z − ∂p
∂x+νρ (∂2u
∂x2 +∂2u
∂z2 )
, (25)
∂ρw
∂t =−∂ρuw
∂x − ∂ρww
∂z − ∂p
∂z −ρg+νρ (∂2w
∂x2 +∂2w
∂z2 )
, (26)
∂ρθ
∂t =−∂ρuθ
∂x −∂ρwθ
∂z . (27)
Here u and w are the horizontal and vertical winds, ρ is the density, θ the potential temperature, g the acceleration of gravity, ν the diffusion coefficient andp the pressure which is given diagnostically by the equation of state
p=
(Rdρθ pκ0
)1−1κ
(28) where Rd is the gas constant for dry air, κ = Rd/cp, cp the heat capacity of dry air at constant pressure and p0 is the pressure at ground level. The red terms belong to the advection part, the black terms are the acoustics and the diffusion is marked in blue. 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 in the cell centers as illustrated by Figure 8.
Because the atmosphere nearly is in hydrostatic equilibrium and the change of pres-sure is small in comparison to its absolute value we compute with its deviation from a background state which is in hydrostatic equilibrium, i.e. this background state satisfies
∂p
∂z =−ρg.
So instead of computing with the black terms in (26) we implement
−∂(p−p)
∂z −(ρ−ρ)g.
In order to obtain linear systems of equations when using implicit methods like the trapezoidal rule for the fast part of split-explicit methods we need a linearized pressure as shown for the one-dimensional case in Section 4.1. Therefore we use the Exner pressure
π=
(Rdρθ p0
)1−κκ
(29) and obtain
− Rd
1−κπ∂ρθ
∂x and − Rd
1−κπ∂ρθ
∂z
Figure 8: The positions of the variables on the Arakawa C grid.
instead of
−∂p
∂x and −∂p
∂z,
i.e. we use the equations (24) to (27) as the right-hand side while the Jacobian originates from
∂ρ
∂t =−∂ρu
∂x − ∂ρw
∂z , (30)
∂ρu
∂t =−∂ρuu
∂x − ∂ρwu
∂z − Rd
1−κπ∂ρθ
∂x +νρ (∂2u
∂x2 + ∂2u
∂z2 )
, (31)
∂ρw
∂t =−∂ρuw
∂x − ∂ρww
∂z − Rd
1−κπ∂ρθ
∂z −ρg+νρ (∂2w
∂x2 +∂2w
∂z2 )
, (32)
∂ρθ
∂t =−∂ρuθ
∂x −∂ρwθ
∂z . (33)
As for the pressurep in the vertical momentum equation we compute without the back-ground state, i.e. with
− Rd
1−κπ∂(ρθ−ρθ)
∂z − Rd
1−κ(π−π)∂ρθ
∂z −(ρ−ρ)g.
For simplicity of notation we will use the formulation where the background state is not subtracted in the remainder of this thesis.
In Section 6.8 we will also present the application of the methods to the dam-break problem which is a test case for the shallow water equations. These can easily be gener-ated from the compressible Euler equations by settingθ = 1, removing the gravity term from the vertical momentum equation (26) and changing the constants in the equation of state for the pressure. In this formulation the equations for ρ and ρθ are redundant and the shallow water equations read
∂ρ
∂t =−∂ρu
∂x − ∂ρw
∂z ,
∂ρu
∂t =−∂ρuu
∂x − ∂ρwu
∂z − ∂p
∂x,
∂ρw
∂t =−∂ρuw
∂x − ∂ρww
∂z − ∂p
∂z
with
p= g 2ρ1−1κ.
Now ρ is the depth of the water and u and w are the velocities.
The Euler equations are implemented in MATLAB on an Intel Core 2 Quad Q9550
@ 2833 Mhz with 3 GB RAM.