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 ±

2±

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