Spatial Discretization and Dispersion: 2.1 Advection Equation
Spatial Discretization and Dispersion: 2.1 Advection Equation
Spatial Discretization and Dispersion: 2.1 Advection Equation
17
12.950 Atmospheric and Oceanic Modeling, Spring ’04 18
The phase speed �/k = c and group speed �k � = c are equal and constant
the waves are said to be non-dispersive.
Approximating the spatial derivative with the centered second order deriva
tive yields:
c i
dt ω + ∂i ω = 0
�x
or
c
d t ωi + (ωi+1 − ωi−1 ) = 0
2�x
Substituting in our solution of the form ei(kx−�t) gives:
c � ik�x �
−i� = − e − e−ik�x
2�x
ci
= − sin k�x
�x
In the limit of small �x, the frequency of long waves (k�x << 1) is:
c �x�0
�= sin k�x = ck
�x
so the approximation appears to converge. However, the numerical waves
are dispersive since the group and phase speeds are not constant with wave
number. Also, the frequency of the grid-scale waves (wavelength of 2�x �
k = η/�x) is always zero regardless of the resolution.
Now consider the fourth order accurate representation of the spatial
derivative:
�i
c 1
�
dt ω + ∂i ω − ∂ii ω = 0
�x 6
or
c
dt ω + ∂i (ωi−2 − 8ωi−1 + 8ωi+1 − ωi+2 ) = 0
12�x
Substituting in a wave solution yields the dispersion relation:
c � 2ik�x �
−i� = −e + 8eik�x − 8e−ik�x + e−2ik�x
12�x
ci
3.5
ω’ =πk’
ω’ =2s c
2 k k
2
3 ω’ =2s c (1+2/3s )
4 k k k
2 2
ω’ =2s c (1+2/3s (1+4/5s ))
6 k k k k
2.5
2
ω’=ωΔx/c
1.5
0.5
0
0 0.2 0.4 0.6 0.8 1
k’=kΔx/π
Figure 2.1: Dispersion relations for constant flow advection in the contin
uum, and using second and fourth order spatial differences.
a)
u,η u,η u,η
x
i−1 i i+1
b)
η u η u η
x
i−1 i−1/2 i i+1/2 i+1
Figure 2.2: A regular grid with the gravity wave variables u and � co-located
on an non-staggered grid (a) and on a staggered grid (b).
The pattern that emerges is that the grid-scale modes are always station
ary and that increasing the order of the scheme does not necessarily improve
the behavior of the shortest scales.
�t u = −g�x �
�t � = −H�x u
where g is gravity, H is the nominal depth of the layer, u is the flow and �
the free-surface displacement. These first order equations can be re-arranged
into a second order wave equation:
�tt � = gH�xx �
Now consider the second order spatial discretization of the gravity wave
equations on a regular grid corresponding to Fig. 2.2a:
g
�t u = − ∂i � i (2.1)
�x
H
�t � = − ∂i ui (2.2)
�x
(2.3)
We can follow the same procedure of substituting the first equation into the
second to obtain:
gH
�tt � = ∂ii � ii
�x2
The discrete operator on the right hand side in full form is:
1
∂ii � ii = (�i−2 − 2�i + �i+2 )
4
and so substituting in the solution � = �o ei(kx−�t) into the discrete wave
equation gives:
gH � −i2k�x �
−� 2 = e − 2 + e i2k�x
4�x2
gH
= (2 cos 2k�x − 2)
4�x2
4gH k�x k�x
= − 2
sin2 cos2
�x 2 2
Here, we’ve expressed the frequency as a function of k�x/2 since 0 �
sin k�x
2
� 1 and is monotonic. Note that the largest discrete wave num
ber permitted is the Nyquist mode, k = η/�x.
In order to plot these relations, we will use a non-dimensional wave num-
�k ��x
ber k � = �x and plot � gH
. The curves for the continuum and non-staggered
grid are plotted in Fig. 2.3. From these curves we note the following:
• the numerical gravity waves are dispersive while the continuum is non-
dispersive,
• the effective group speed has the wrong sign for modes higher than
�
k = 2�x (or wavelengths shorter than 4�x),
12.950 Atmospheric and Oceanic Modeling, Spring ’04
22
3.5
ω’=πk’
ω’=2sin(πk’/2)cos(πk’/2)
3 ω’=2sin(πk’/2)
2.5
ω’=ωΔx/cg
1.5
0.5
0
0 0.2 0.4 0.6 0.8 1
k’=kΔx/π
Figure 2.3: Dispersion of numerical gravity waves for the non-staggered grid
(green) and the staggered grid (red). The continuum (� = k) is plotted for
comparison (blue).
�
• the Nyquist mode (k = �x
or k � = 1) is stationary (i.e. has zero
frequency).
• when compared to the continuum we see that the numerical modes are
still dispersive,
The remarkable improvement in behavior of the staggered grid over the non-
staggered grid is due to the complete absence of operators with null-spaces
(such as the averaging operator). The formal truncation error is smaller but
the improved accuracy is not responsible for the qualitative improvements.
12.950 Atmospheric and Oceanic Modeling, Spring ’04 24
H�x 0 �t �
and taking the determinant of the operator matrix
�� ⎛�
� −i� −f gik �� ⎠
� �=0
��
f −i� 0 ⎝� = 0 �
⎜�
2
� = f 2 + gHk 2
��
Hik 0 −i� �
� �
�
A)
u,v,η u,v,η u,v,η
x
i−1 i i+1
B)
η u,v η u,v η
x
i−1 i−1/2 i i+1/2 i+1
C)
v,η u v,η u v,η
x
i−1 i−1/2 i i+1/2 i+1
D)
u,η v u,η v u,η
x
i−1 i−1/2 i i+1/2 i+1
Figure 2.4: The inertia-gravity wave variables u, v and � in the four possible
configurations (A, B, C and D) on a finite difference grid.
• C-grid model
g
�t u − f v i + ∂i � = 0
�x
�t v + f ui = 0
H
�t � + ∂i u = 0
�x
• D-grid model
g
�t u − f v i + ∂i � i = 0
�x
�t v + f ui = 0
H
�t � + ∂i ui = 0
�x
�2 4L2d 2 2
A: = 1 + s
c
f2 �x2 k k
12.950 Atmospheric and Oceanic Modeling, Spring ’04 26
r=(2Ld/Δ x)=2
12
ω=1+π2r2k2/4
2 2
ω =1+r2s c Continuum
A k k
2 2
ω =1+r s
10 B k
2 2 2
ω =c +r s
C k k
2 2 2 2
ω =c +r s c
D k k k
8
ω’=ω/f
6
B
C
4
2
A
a)
D
0
0 0.2 0.4 0.6 0.8 1
k’=kΔx/π
r=(2L /Δ x)=0.5
d
1.8
ω=1+π2r2k2/4
ω =1+r2s c
2 2 Continuum
1.6 A k k
2 2
ω =1+r s
B k
ω =c +r2s
2 2
1.4 C k
2 2 2 2
k
ωD=ck +r sk ck B
1.2
A
1
ω’=ω/f
0.8
0.6
0.4
C
0.2
b)
D
0
0 0.2 0.4 0.6 0.8 1
k’=kΔx/π
�2 4L2d 2
B: = 1 + s
f2 �x2 k
�2 2 4L2d 2
C: = c k + s
f2 �x2 k
�2 2 4L2d 2 2
D: = c k + s c
f2 �x2 k k
�
where Ld = gH/f is the radius of deformation and sk = sin k�x 2
and
k�x
ck = cos 2 and short-hand of the common factors. The coefficient in front
of the gravity-wave contribution (s2k ) is a resolution parameter:
�2
2Ld
�
r 2 =
�x
and is controls the form of the some of the dispersion relations. The dis
persion relation of the continuous and four numerical inertia-gravity waves
are plotted in Fig. 2.5. In panel (a) the resolution parameter is larger than
one meaning the deformation radius is resolved. In panel (b) the resolution
parameter is less than one meaning the deformation radius is not resolved.
In the continuum, the frequency increases monotonically with wave number,
the waves are dispersive and the group and phase speeds are pointed in the
direction of the wave vector. None of the numerical schemes behave precisely
in this manner:
• Scheme A, regardless of resolution has a false maximum at k � = 1/2
so that grid-scale waves (k � = 1) oscillate with the inertial frequency
rather than propagate as gravity waves,
• Scheme D has false maxima but is worse in that the grid-scale waves
become stationary,
• Scheme C has the wrong slope (group speed in wrong direction) if r < 1
but is qualitatively correct if r > 1,
• Scheme B is qualitatively correct for both values of r and so appears
to be the best scheme.
This analysis suggests that staggering variables in the form of the B grid is
most likely to avoid computational modes when solving these one-dimensional
shallow water equations. However, the results of an analysis in two dimen
sions are different.