<<

. 58
( 136 .)



>>

rather with u for that same ¬‚uid particle at earlier time levels as shown schematically in
Fig. 14.3. Advection is invisible in (14.37) because it is implicit in the coordinate transform,
that is, in the spatial shifts represented by ± in (14.39).
This advection-hiding has two major advantages. First, since there is no advection in
the time extrapolation formula (14.37), there is no advective stability limit either. (The ex-
plicitly treated term F (u, t) may impose its own maximum time step, but in meteorological
applications, this is huge compared to the advective limit.) Second, the SL scheme is much
better at resolving steep fronts and shock waves than conventional schemes. The reason is
that the frontogenetic process, advection, is computed exactly by the method of character-
istics. In contrast, the treatment of an advective term like u ux in a conventional, non-SL
method depends on how accurately ux is approximated by spatial ¬nite differences or a
pseudospectral formula ” poorly in the neighborhood of a front.
SL schemes, alas, require two additional steps which are unnecessary in other algo-
rithms. First, one must evaluate u and other variables at x = xj ’ ±j , which generally is
not on the canonical spectral grid. This problem of “off-grid” interpolation is suf¬ciently
important to warrant an entire section below. However, ef¬cient spectral off-grid interpo-
lation methods are now available.
Second, one must solve the algebraic equation (14.39) at each grid point to compute
the corresponding displacement. The good news is that because a ¬‚uid blob only moves
a short distance over the short time interval of a single time step, a simple ¬xed point
iteration is usually suf¬cient. That is, denoting the iteration level by superscript m:

= „ u(xj ’ ±j , tn )
m+1 m
(14.40)
±j




u+
t n+1
u0
tn
u-
t n-1 ±

x

Figure 14.3: Schematic of a particle trajectory (solid curve) which arrives at one of the
points of the regular, evenly spaced computational grid at t = tn+1 . The semi-Lagrangian
scheme approximates this trajectory by the dashed line, using the velocity at the midpoint
of the trajectory to estimate the displacement ± during one time step. The SL algorithm
connects the values of u at three points along a trajectory (circles) rather than at a ¬xed
value of x, as in a conventional scheme.
14.6. MULTIPLY-UPSTREAM SL 275

0
At the beginning of the time integration, one can take ±j = 0; at later time steps, one can
initialize the iteration with ±j (tn’1 ), that is, the displacement at the previous time level.
Experience and theory both suggest that one never needs to take more than two ¬xed point
iterations except possibly at the start of the time integration where one might need three.
Each iteration increases the accuracy of the displacement by one more power of „ . The
iteration is guaranteed to converge so long as
„ | ux | < 1 (14.41)
This condition seems to be almost always met in practice though one might worry a bit
in the neighborhood of fronts where the slope ux is very large (Pudykiewicz, Benoit and
Staniforth, 1985).
The three-level SL/SI scheme generalizes easily to higher dimensions. In two dimen-
sions, let ±(x, y) and β(x, y) denote the displacements in the x and y directions, respec-
tively. The three-level SL/SI method is
≡ u(xi , yj , tn + „ )
u+
ij

≡ u(xi ’ ±ij , yj ’ βij , tn )
u0
ij
u’ ≡ u(xi ’ 2 ±ij , yj ’ 2βij , tn ’ „ )
ij
(14.42)
and similarly for the y-momentum equation. The displacements are the solutions to
±ij = „ u(xij ’ ±ij , yij ’ βij , yij , tn ) βij = „ v(xij ’ ±ij , yij ’ βij , yij , tn )
&
(14.43)
where v is the y velocity.


14.6 Multiply-Upstream: Stability with a Very Long Timestep
We begin with three useful de¬nitions.
De¬nition 31 (FOOT of a TRAJECTORY/DEPARTURE POINT)
The spatial location of the beginning of the trajectory of a ¬‚uid particle at some time t, usually tn’1
for a three-level SL scheme or tn for a two-level method, is called the “FOOT of the TRAJECTORY”
or equivalently the “DEPARTURE POINT” of the trajectory. The goal of the ¬xed point iteration
of the SL method is to compute the displacement and through it the “foot”.
De¬nition 32 (MULTIPLY-UPSTREAM)
A semi-Lagrangian method is “MULTIPLY-UPSTREAM” if the set of grid points used to interpo-
late the velocities in the computation of the “foot” of a trajectory is allowed to vary, relative to the
arrival point of the trajectory, so as to center on the foot.
De¬nition 33 (COURANT NUMBER)
The Courant number is a nondimensional parameter which is the ratio of the time step to the
smallest distance between two points on the spatial grid, multiplied by the maximum velocity:

Co ≡ Umax (14.44)
min(h)
where Umax is the maximum ¬‚uid velocity. Strictly speaking, this is an advective Courant number;
one can measure the maximum time step for wave propagation problems by de¬ning a “wave”
Courant number with the maximum advecting velocity replaced by the maximum phase speed c.
In either case, explicit methods are stable only when the Courant number is O(1) or less. In the
rest of this chapter, we shall use only the advective Courant number.
CHAPTER 14. SEMI-LAGRANGIAN ADVECTION
276



SINGLE GRID
t


x

MULTIPLY UPSTREAM
t


x

Figure 14.4: Schematic comparing the off-grid interpolation schemes of multiply upstream
and single grid semi-Lagrangian methods for a large Courant number, i. e., when the
particle travels through several grid points during the time interval t ∈ [tn , tn+1 ]. The
circles show the points whose grid point values at time level tn’1 are used to interpolate
velocities and other quantities at the foot of the trajectory. The trajectory is the heavy curve.


Fig. 14.4 illustrates this de¬nition. One major reason for the slow acceptance of semi-
Lagrangian schemes, which were used as early as Wiin-Nielsen (1959), Krishnamurti(1962)
and Sawyer (1963), is that the early algorithms interpolated the departure points of the
trajectory through xj by using (uj’1 , uj , uj+1 ). That is to say, the “single grid” SL methods
used the same set of points centered on the arrival point xj irregardless of whether the foot
of the trajectory was near or far, left or right, of the arrival point.
The price for this in¬‚exibility is that when the Courant number is greater than one,
single grid schemes use extrapolation rather than interpolation to estimate velocities, etc., at
the foot of the trajectory ” and are unstable. Consequently, the old-style semi-Lagrangian
algorithms had time step restrictions which are just as strict as for conventional explicit
methods.
In contrast, multiply-upstream schemes are stable and accurate even for very large
Courant numbers. Indeed, as we shall see show in the next section, something very wierd
and wonderful happens: multiply-upstream methods become more accurate as the time
step increases!
Spectral and pseudospectral schemes are always multiply-upstream.2 Thus, spectral SL
algorithms can be used at high Courant number.


2 This
assertion is true if the interpolation to compute the foot is done spectrally; some weather forecasting
models inconsistently use low order polynomial interpolation to compute the foot, even though the rest of the
computation is spectral; such mixed schemes are not automatically multiply-stream.
14.7. NUMERICAL ILLUSTRATIONS & SUPERCONVERGENCE 277

14.7 Numerical Illustrations: Higher Accuracy for Longer
Time Step and Superconvergence
Fig. 14.5 compares the standard Fourier pseudospectral method with the semi-Lagrangian
method for the one-dimensional advection equation. This solution “breaks” at t = 1 and
has a jump discontinuity for all larger times. The time step „ = 1/2 for the semi-Lagrangian
scheme ” much longer than is stable for the standard, explicit Fourier method ” is such
that the numerical integration jumps from a very smooth ¬‚ow (a sine wave!) to a front in
just two time steps! Nevertheless, the semi-Lagrangian method gives a much smoother ap-
proximation to the front than the standard Eulerian coordinates algorithm, which blindly
approximates all derivatives by the usual Fourier pseudospectral formulas in complete
ignorance of both the method of characteristics and the hyperbolic character of the differ-
ential equation.
Fig. 14.6 shows the departure points and errors for the semi-Lagrangian calculation,
same case as the previous ¬gure but at a later time when the errors have presumably wors-
ened. Even so, the error is quite small everywhere except for a narrow neighborhood of the
front. However, large relative error in a thin zone around the front is a weakness of almost
all front-resolving algorithms.
The bottom panels in Fig. 14.6 shows the same solution using ten times as many time
steps. This should reduce the error, since the method is second order, by a factor of 100.
Instead, the error is considerably worse. Why does a larger time step give smaller error?


1
1
0.8

0.6
0.95
0.4

0.2
0.9
0
-0.2
0.85
-0.4

-0.6 0.8
-0.8
-1 0.75
-2 0 2 2.2 2.6 3
x
x

Figure 14.5: Comparison of semi-Lagrangian/semi-implicit Fourier method (thick solid
curve) with standard Fourier pseudospectra method (thin line with circles) for the one-
dimensional advection equation [inviscid Burgers™ equation] , ut + uux = 0 for 64 grid
points, evolving from u(x, 0) = sin(x) and illustrated at t = 1. The right graph is just a
“zoom” or blow-up of the left plot. The semi-Lagrangian time step was „ = 1/2; the time
step for the standard method was determined by an adaptive 4th/5th order Runge-Kutta
subroutine and is not known except that it was much shorter than for SL scheme. The
solution “breaks” at t = 1, that is, would develop an in¬nite slope at that time were it not
for the viscosity.
CHAPTER 14. SEMI-LAGRANGIAN ADVECTION
278



timestep= 1/2
Particle Origins: u(x,t=2),
2 .1
ξ(x,t)




.01
1




.001

0

.0002 π
π/2
π/2 π
x
0 0 x


timestep= 1/20
ξ(x,t) u(x,t=2):
2
„=1/20
.1


u

.01
1




.001



.0002 π
π/2
0
π 0
π/2 x
0


Figure 14.6: Semi-Lagrangian Fourier calculation of a moving shock wave in Burgers™ equa-
tion (in the limit of vanishing viscosity) with 64 points on x ∈ [0, 2π]. u(x, 0) = sin(x). The
solution breaks, i. e., ¬rst develops a near-discontinuity, at t = 1; the solution is shown at
t = 2. Left panels: the origins of the points which lie on the pseudospectral grid at t = 2,
i. e., the graph is a plot of ξ(x, t) ≡ xj ’ ±j versus xj . Right panels: error at t = 2. (A
spike of O(1) error at one point in both graphs has been truncated.) In the upper pair of
plots, the time step is „ = 1/2. Since the grid spacing is about 1/10, the Courant number
is roughly 5. Bottom panels: same as top except that the time step has been reduced by a
factor of 10. The main changes are (i) the gap in particle origins is much smaller and (ii) the
error has increased. Note that in the upper panels (where Co ≈ 5), the error is smaller than
1/100 almost everywhere; with the short time step, the error is larger than 1/100 almost

<<

. 58
( 136 .)



>>