involves only a ¬rst derivative and demands a time step inversely proportional to kmax

itself rather than to its cube.

The exponential-of-derivatives operator looks forbidding, but it can be interpreted in

the usual way as the exponential of the matrix that discretizes the third spatial derivative,

which in turn can be interpreted as multiplication of the k-th eigenvector of the derivative

matrix by the exponential of (’t) times the k-th eigenvalue. In a Fourier basis, this is

trivial because the basis functions, exp(ikx), are the eigenvalues of the derivative operator.

It follows that for any function v(x, t) with Fourier coef¬cients vk (t),

∞

exp(i k 3 t) vk exp(ikx) (14.24)

exp(’‚xxx t) v(x, t) =

k=’∞

14.3. SEMI-LAGRANGIAN ADVECTION: INTRODUCTION 269

To put it another way, the integrating factor method for the KdV equation is equivalent to

transforming the unknowns of the ODE system which comes from the usual Galerkin or

pseudospectral discretization. In place of the Fourier coef¬cients ak (t) of u(x, t), use the

modi¬ed Fourier coef¬cients (of v(x, t)) given by

vk (t) ≡ exp(’ik 3 t) ak (t) (14.25)

Unfortunately, this is really a terrible idea when the prime target is solitary waves. The

simplest counterexample is the special case that the exact solution is a solitary wave of

phase speed c. The exact Fourier coef¬cients ” before the integrating factor ” are

ak = constant exp(’ikct) (14.26)

and thus vary rather slowly in the sense that the frequency of the highest wavenumber

in the truncation grows only linearly with kmax . In contrast, the time dependence of the

transformed coef¬cients is

vk (t) = constant exp(’ik 3 t ’ ikct) (14.27)

which varies rapidly in the sense that the highest frequency after transformation grows as

3 3

kmax . To accurately compute vkmax (t), a time step proportional to kmax is necessary ”

in other words, a time step as short as the explicit stability limit for the original, untrans-

formed KdV equation.

In discussing implicit methods, we have stressed repeatedly that their usefulness is

SITUATION-DEPENDENT. Implicit methods are good for computing KdV solitons because

they misrepresent only the high frequency transient while accurately computing the slow-

moving solitary waves. Integrating factor methods, in contrast, are bad for solitons because

the integrating factor extracts a high frequency motion that is not actually present in the

low frequency motion (i. e., the soliton) that we are computing. Thus, the usefulness

of integrating factor methods is also situation-dependent, but often inversely to implicit

methods.

Suppose our goal is to compute KdV solutions which are not solitons but rather a soup

of weakly interacting waves such that the Fourier components are approximately

ak = Ak ( t) exp(ik 3 t) (14.28)

where the Ak are amplitudes that vary slowly under weak wave-wave interactions. Now

because the integrating factor method successfully removes the high frequency oscillations,

the transformed coef¬cients are

(14.29)

vk (t) = Ak ( t)

These vary slowly with time, allowing a long time step. In contrast, implicit and Nonlinear

Galerkin methods as well as standard explicit schemes all require a very short time step

to be accurate for the high frequency motion. Weakly coupled wave systems (“resonant

triads” and their generalization) have been widely analyzed by theorists (Craik, 1985).

Thus, it is quite wrong to conclude that integrating factor methods are counter-productive

in numerical applications. Sometimes they are very valuable. The difference between the

two KdV cases is that in the ¬rst (soliton) case, the integrating factor extracts a high fre-

quency dependence that was not actually present in the slow solutions (bad!) whereas in

the second (weakly interacting sine waves) case, the integrating factor extracts high fre-

quency motion that is an essential part of the desired solution (good!).

CHAPTER 14. SEMI-LAGRANGIAN ADVECTION

270

14.3 Semi-Lagrangian Advective Schemes: Introduction

A “semi-Lagrangian” (SL) time-marching algorithm is an integrating-factor method in

which the integrating factor is an advection operator that shifts into a coordinate system

that travels with the ¬‚uid. In so doing, it borrows ideas from both the method of char-

acteristics and the Lagrangian coordinate formulation of ¬‚uid mechanics. The “total” or

“particle-following” derivative of any quantity “q”, for simplicity given in two space di-

mensions, is

Dq ‚q ‚q ‚q

≡ (14.30)

+u +v

Dt ‚t ‚x ‚y

In the moving coordinate system, the total derivative becomes the ordinary partial time

derivative, and the advective terms disappear, buried in the change-of-coordinate.

The idea of tracing ¬‚ow by following the motion of individual blobs or particles of

¬‚uid is at least two centuries old. Particle-following coordinates are called “Lagrangian”

coordinates in ¬‚uid mechanics. However, exclusive use of a Lagrangian coordinate system

is usually a disaster-wrapped-in-a-catastrophe for numerical schemes because the parti-

cle trajectories become chaotic and wildly mixed in a short period of time where “chaos”

is meant literally in the sense of dynamical systems theory even for many laminar, non-

turbulent ¬‚ows.1 Semi-Lagrangian (SL) algorithms avoid this problem (and earn the mod-

i¬er “semi”) by reinitializing the Lagrangian coordinate system after each time step. The

usual choice is that the tracked particles are those which arrive at the points of a regular

Fourier or Chebyshev spatial grid at time tn+1 . Since the departure points for these trajec-

tories at t = tn are generally not on the regular grid, we follow a different set of N particles

on each time interval of length „ where „ is the time step. In contrast, in a pure Lagrangian

coordinate system, one follows the evolution of the same N particles forever.

There are two advantages to an SL method:

1 Trajectories of ¬‚uid blobs are often chaotic even when the advecting velocity ¬eld is laminar.

t5

t5

t4 t4

t3

t3

t2

t2

t1 t1

t0 t0

x x

Figure 14.1: Left panel: Fluid mechanics in Lagrangian coordinates tracks the same set

of particles throughout all time. Right panel: Semi-Lagrangian methods. A new set of

particles to track is chosen at the beginning of each time step.

14.4. ADVECTION & METHOD OF CHARACTERISTICS 271

1. If the scheme is “multiply upstream”, a term de¬ned later, SL will be both stable

and accurate for time steps far greater than usual explicit limit. To be precise, the

advective time step limit disappears, and the maximum time step is then controlled

entirely by non-advective processes like viscosity and wave propagation.

2. SL methods are very good at resolving fronts and shock waves even in the absence

of arti¬cial viscosity.

Because of these virtues, SL time-marching is now employed in the NCAR Community

Climate Model and the European Centre for Medium-Range Forecasting weather predic-

tion model. Neither is a toy model, as too often used to test algorithms; instead, both are

global, three-dimensional spherical harmonic models with lots of physics.

Simple semi-Lagrangian algorithms have at least ¬ve disadvantages:

1. Dif¬cult to parallelize (for “multiply upstream” schemes).

2. Not strictly conservative of mass, energy, etc.

3. Inherent computational dissipation (McCalpin, 1988)

4. Instability and inaccuracy for ¬‚ows over mountains.

5. More expensive per time step than most competing schemes because of the need for

• “off-grid” interpolation and

• ¬xed point iteration.

None of these drawbacks, however, has proved fatal. In particular, conservative, shape-

preserving and low dissipation schemes have been recently developed, and the mountain

¬‚ow or “orographic” problems have been solved, too.

Semi-Lagrangian method are designed to be especially good at advection. However,

the algorithms can be applied to the full Navier-Stokes equations without dif¬culty and

can be combined with other ideas. For example, in meteorology it is common to combine

semi-Lagrangian treatment of advection with a semi-implicit treatment of the viscous and

Coriolis force terms (“SL/SI” schemes).

Semi-Lagrangian time-stepping can be combined with any strategy for discretization of

the spatial coordinates.

14.4 A Brief Digression on the One-Dimensional and Two-

Dimensional Advection Equations and the Method of

Characteristics

The Navier-Stokes equations reduce to the so-called “advection equations” if we neglect

all processes except advection. Restricting ourselves to two-dimensions for simplicity, the

pure advection equations are

Du Dv

(14.31)

=0 & =0

Dt Dt

or written in full, where (u, v) are the x and y velocities,

‚u ‚u ‚u

+u +v = 0

‚t ‚x ‚y

‚v ‚v ‚v

(14.32)

+u +v = 0

‚t ‚x ‚y

CHAPTER 14. SEMI-LAGRANGIAN ADVECTION

272

Figure 14.2: Contour plot in the space-time plane of the solution u(x, t) to the one-

dimensional advection equation, ut + uux = 0. The initial condition is sinusoidal; negative

contours are dashed. Where contours intersect, the solution becomes multi-valued and

ux ’ ∞. The ¬‚ow is said to have “formed a front” or experienced “wave breaking”.

The reason that this drastically simpli¬ed system is interesting is that it can be solved ex-

actly by the “method of characteristics”. Semi-Lagrangian time-marching schemes, even

with all the complications of pressure and viscosity, etc., merely mimic the method of char-

acteristics.

The key observation is that the advection equations state that the velocities are con-

served following the motion. In other words, if we track the fate of a particular blob of

¬‚uid, perhaps marked with a dab of red dye, the velocity of that blob will remain constant

forever. Consequently, the trajectory of that blob and of all blobs is a straight line whose

slope in each spatial dimension is proportional to its velocity with respect to that coordi-

nate. (Fig. 14.2 is a one-dimensional illustration.) If the blob originates at x = ξ and y = ·

at t = 0, then its trajectory is

(14.33)

x = ξ + Q(ξ, ·) t, & y = · + S(ξ, ·) t

where the initial conditions are

(14.34)

u(x, y, t = 0) = Q(x, y), & v(x, y, t = 0) = S(x, y)

The exact solution for the velocities is

(14.35)

u(x, y, t) = Q(ξ[x, y, t], ·[x, y, t]), & v(x, y, t) = S(ξ[x, y, t], ·[x, y, t])

In other words, the initial conditions furnish the exact solution for all time if the arguments

of the functions (Q, S) are replaced by the “characteristic coordinates” or simply “charac-

teristics”, (ξ, ·).

14.5. THREE-LEVEL, 2D ORDER SEMI-IMPLICIT 273

Unfortunately, Eq.(14.33) gives explicit solutions only for (x, y) as a functions of the

characteristics. To evaluate the solutions at a particular point (x, y, t), we need it the other

way around. Thus, it is necessary to solve (14.33) to determine the departure point of a

given blob where the “departure point“ or “foot of the trajectory” is the location (ξ, ·) of

the blob at t = 0. However, the equations that must be solved for the characteristics are

algebraic equations, which are much simpler than differential equations.

Semi-Lagrangian time-marching schemes also require computing the departure point

of a blob by solving a system of algebraic equations. However, this is fairly easy because

the blob moves only a short distance during a time interval of length „ where „ is the time

step. The change of coordinates (14.33) does not trivialize the dynamics of the full Navier-

Stokes equations, as it does for the pure advection equations, but it does remove advection

just as thoroughly.

14.5 The Simplest SL Scheme: Three-Level, Second Order

Semi-Implicit Method

A semi-Lagrangian scheme by itself is not very useful in meteorology because it only re-

moves the advective time step limit, which is much less severe than the restriction imposed

by the propagation of gravity waves. When combined with a semi-implicit treatment of the

non-advective terms, however, SL becomes very powerful as ¬rst demonstrated by Robert

(1981, 1982).

Let the problem be

Du

(14.36)

+ G(u, t) = F (u, t)

Dt

where D/Dt denotes the total, particle-following derivative including the advective terms,

G is the sum of all terms which are to be treated implicitly while F represents all the terms

which can be treated explicitly. The basic three-level, second order SL/SI algorithm is

u + ’ u’ 1

+ (G+ + G’ ) = F 0 (14.37)

2„ 2

where the plus and minus signs denote the quantities evaluated at t = tn±1 and super-

script 0 denotes a quantity evaluated at t = tn . In form, this looks like a standard Crank-

Nicholson/explicit method except that (i) the advective terms are missing (ii) the CN-like

step is over a time interval of 2„ instead of „ where „ is the time step.

The magic is in the coordinate transform. For simplicity, assume one space dimension.

Let ±(x) denote the displacement, over the time interval t ∈ [tn , tn+1 ] of a particle that

arrives at point x at t = tn+1 . Then for any vector appearing in (14.37),

≡ u(x, tn + „ )

u+

≡ u(x ’ ±, tn )

u0

u’ ≡ u(x ’ 2 ±, tn ’ „ )

(14.38)

When the spatial dependence is discretized, u, x, and ±(x) all turn into vectors. The dis-

placement is de¬ned to be the solution of the algebraic equation

±j = „ u(xj ’ ±j , tn ) (14.39)

where j is the integer index of the grid points.

CHAPTER 14. SEMI-LAGRANGIAN ADVECTION

274

The novelty of Eq.(14.39) is that when we move back one time level, we must simulta-

neously shift our spatial position by the displacement ±j . The predicted value u(xj , tn+1 )

is combined not with u at the same spatial location at the two previous time levels, but