un Advection

un+1/3

Waves & Pressure

Adjustment

un+2/3 Viscosity/Diffusion

un+1

Figure 13.4: Splitting by process: Advection, adjustment of ¬elds (linear wave dynamics

and pressure), and viscosity and diffusion are split into three separate, one-after-the-other

processes.

n+2/3

’ f v n+2/3 = ’φx

n+2/3

(13.49a)

ut

n+2/3

+ f un+2/3 = ’φy

n+2/3

[Pressure Step alias (13.49b)

vt

n+2/3 n+2/3 n+2/3

“Adjustment of Fields”] (13.49c)

φt + ux + vy =0

un+1 = µ un+1 (13.50a)

t

n+1

v n+1 [Diffusion Step] (13.50b)

vt =µ

φn+1 = φn+2/3 (13.50c)

In each case, the solution at the previous fractional step provides the initial condition for

the equations that model the next physical process.

If we wish, we can apply an additional splitting to decompose the diffusion step into

a sequence of one-dimensional diffusion problems. However, it is usually a bad idea to

CHAPTER 13. SPLITTING & ITS COUSINS

264

attempt dimensional splitting for the “adjustment of ¬elds” step because the Coriolis force

and wave dynamics strongly couple the coordinates. Splitting advection into separate x,

y and z advection is unnecessary because this subproblem is normally integrated by an

explicit method.

As noted earlier, the great advantage of splitting is that we can solve each fractional

step using a different numerical method. The advective step is almost always advanced

explicitly and the adjustment of ¬elds is integrated implicitly. In engineering problems,

diffusion is usually marched in time by an implicit algorithm, but Marchuk (1974) barely

discusses diffusion because µ in his numerical weather prediction model is so small that

he can use an explicit method.

It is not even necessary to be consistent by using spectral or pseudospectral algorithms

for all sub-problems. Lambiotte et al. (1982) compute only the advection and pressure steps

pseudospectrally. Since µ is small for their high Reynolds ¬‚ow, they employ ¬nite differences

to integrate the diffusive step.

Chapter 14

Semi-Lagrangian Advection &

Other Integrating Factor Methods

“Go with the ¬‚ow”

“ New Age proverb

“... there is a mistaken belief that semi-Lagrangian schemes are only good for smooth

¬‚ows.”

ˆe

“ Staniforth & Cot´ (1991), pg. 2211

14.1 Introduction: Concept of an Integrating Factor

An “integrating factor” is a transformation of the unknowns in a differential equation. The

most familiar example yields the solution of the general linear ¬rst order ODE, which is

ux ’ q(x) u = f (x) (14.1)

Let u(x) be replaced by the new unknown v(x) where

u(x) ≡ I(x) v(x) (14.2)

where the integrating factor is

x

I(x) ≡ exp (14.3)

q(y) dy

Substitution gives the transformed problem

I vx + q(x) I v ’ q(x) I v = f (x) (14.4)

using Ix = q(x)I. Dividing by I(x) gives

x

’ (14.5)

vx = exp q(y) dy f (x)

265

CHAPTER 14. SEMI-LAGRANGIAN ADVECTION

266

Thus, the integrating factor has reduced the differential equation to a quadrature, i. e., v(x)

is the inde¬nite integral of the right-hand side of (14.5):

x z

’ (14.6)

v =A+ exp q(y) dy f (z) dz

0

where A = u(0).

Integrating factors come in all sorts of ¬‚avors, like ice cream, but we shall concentrate

on two particular classes. The ¬rst is when a time-dependent PDE can be written in the

form

(14.7)

ut = L u + M (u)

where L is a linear differential operator with spatial derivatives only and where M (u) is

an operator, possibly nonlinear, that contains everything else in the equation including the

forcing terms. The new unknown de¬ned by

u ≡ exp(Lt) v (14.8)

solves the simpli¬ed equation

(14.9)

vt = exp(’Lt) M (v)

Just as for the ¬rst order ODE, the integrating factor has purged the transformed equation

of the differential operator L except through the exponential. This is sometimes a good

idea ” and sometimes not.

The second class of integrating factors is to remove translation/advection by a change-

of-coordinates to a moving coordinate system. Formally, this change-of-coordinate can

be expressed in terms of operators, too. However, it is more convenient to perform the

change-of-coordinate directly.

The simplest case is to shift into a frame of reference which is moving at a constant

speed with respect to the original frame. For example,

(14.10)

ut + cux = f (x, t), u(0, t) = g(t), u(x, 0) = Q(x)

can be transformed by the shift

ξ ≡x’ct ” u(x, t) ≡ v(ξ, t) (14.11)

x=ξ+ct &

into

(14.12)

vt = f (ξ + ct, t), v(’ct, t) = g(t), v(ξ, 0) = Q(ξ)

For the special case f = 0, this is clearly an improvement on the original solution because

then vt = 0 which implies that v is a function of ξ only, that is, u depends on (x, t) only as

the combination x ’ ct.

Semi-Lagrangian time-marching schemes remove nonlinear advection through a sim-

ilar but more complicated transformation. Such algorithms have so many good proper-

ties that they have become standard on two important geophysical models: the NCAR

Community Climate Model (CCM) and the operational forecasting code of the European

Centre for Medium Range Forecasting, which provides ¬ve-day global predications to an

18-nation consortium.

However, integrating factor schemes can also be a trap. The crucial distinction is that

integrating factors are useful only when the integrating factor removes an important part

of the dynamics. We will explain this through a couple of bad examples in the next section.

(This section can be skipped without loss of continuity.)

In the remainder of the chapter, we describe semi-Lagrangian schemes. We begin with

second order schemes, which are the easiest to describe, and then explain a number of

elaborations.

14.2. MISUSE OF INTEGRATING FACTOR METHODS 267

14.2 Misuse of Integrating Factor Methods

Our ¬rst counterexample is the ODE which was used earlier to explain the concept of

the slow manifold:

(14.13)

ut + »u = f ( t)

Assume »/ >> 1, i. e., that the ratio of the “fast” time scale to the slow time scale is

large, and also that the real part of » is positive so that the homogeneous solution decays

exponentially with time. We shall further assume that the goal is to compute the slow

manifold of this equation, that is, the part which varies only on the same slow O(1/ ) time

scale as the forcing f .

Unfortunately, the time step for an explicit scheme is restricted by the fast time scale of

the homogeneous solution exp(’»t). Implicit methods and Nonlinear Galerkin methods

both remove this restriction, but have their own disadvantages. Is an integrating factor a

viable alternative?

De¬ne a new unknown by

(14.14)

u = exp(’» t)v(t)

The ODE becomes

(14.15)

vt = exp(» t)f ( t)

Hurrah! The term »u has been removed and with it the usual CFL time-step limitation. To

show this, let us analyze the Forward Euler method as applied to (14.15).

This scheme with time step „ and with superscript n denoting the time level is

v n+1 ’ v n

” v n+1 = v n + „ exp(» n „ )f ( n „ ) (14.16)

= exp(» n „ )f ( n „ )

„

From an arbitrary starting value v 0 and for the special case f ( t) = F , a constant,

= v 0 + „ F {1 + exp(» „ ) + exp(2 » „ ) + . . . + exp([n ’ 1] » „ )}

vn

1 ’ exp(n » „ )

= v0 + „ F (14.17)

1 ’ exp(» „ )

Because of the integrating factor transformation, v n is stable for arbitrarily large time steps

„ . To be sure, the solution is blowing up exponentially, but since the transformed forcing

is growing exponentially, a similar growth in v is okay.

The same solution in terms of the original unknown u is, using v 0 ≡ u0 ,

exp(’n » „ ) ’ 1

un = exp(’n »„ )u0 + „ F (14.18)

1 ’ exp(» „ )

In the limit that „ ’ 0 and n ’ ∞ such that t = n „ grows large, the Forward Eu-

ler/integrating factor solution simpli¬es to

’1

≈

un „F

1 ’ (1 + » „ + O(»2 „ 2 ))

F

≈ {1 + O(»„ )} (14.19)

»

CHAPTER 14. SEMI-LAGRANGIAN ADVECTION

268

where F/» is the slow manifold for this problem. Thus, for suf¬ciently small time step, the

scheme does successfully recover the slow part of the solution as t ’ ∞.

However, the whole point of the integrating factor is to use a long time step without the

bother of implicit or Nonlinear Galerkin methods. Since the stability limit of the Forward

Euler scheme without integrating factor is „max = 1/», it is interesting to see how well the

integrating factor variant does for this „ . Unfortunately, substituting „ = „max into (14.18)

gives

F 1 F

un ∼ t ’ ∞, „ = 1/» (14.20)

= 0.58 ,

» ’1 + exp(1) »

The good news is that the large time solution is proportional to the true slow manifold,

F/»; the bad news is that the amplitude of the computational solution is barely half that of

the true solution even when „ is no larger than the largest time step that would be stable

in the absence of the integrating factor. When „ >> 1/», which is of course what we really

want if the integrating factor scheme is to compete with implicit schemes,

F

un ∼ t ’ ∞, »„ >> 1 (14.21)

» „ exp(’» „ ),

»

which is exponentially small compared to the correct answer.

Thus, the integrating factor method is a dismal failure. If the time step „ = r „max

where „max = 1/» is the stability limit for the Forward Euler method (without integrating

factor), then the Forward Euler/integrating factor algorithm gives a slow manifold which

is too small by a factor of r exp(’r). In contrast, implicit methods like Crank-Nicholson

and Backwards-Euler and also Nonlinear Galerkin schemes are accurate for the slow part

of the solution even for large r.

A second bad example is the Fourier pseudospectral solution to the KdV equation

(14.22)

ut + uux + uxxx = 0

when the goal is to model the formation and dynamics of solitary waves. For a soliton of

phase speed c, the coef¬cient of exp(ikx) varies with time as exp(’ikct) which grows only

linearly with k. However, the equation also admits small amplitude waves whose phase

speeds grows as k 3 , so the stability limit for an explicit time-marching scheme is inversely

3

proportional to kmax . The KdV equation is very “stiff” in the sense that one must use a

very short time step for stability reasons compared to the slow time scale of the dynamics

(i. e., the solitons) that one is actually interested in.

The remedy ” well, sort of a remedy ” is to use the integrating factor exp(’ ‚xxx t) to

remove the third derivative from the transformed equation:

(14.23)

vt + v exp(’‚xxx t) vx = 0