14.7. NUMERICAL ILLUSTRATIONS & SUPERCONVERGENCE 279

The crucial data is in the left panel of each graph. When the time step is large, there

is a gap in the departure points in the sense that all trajectories originate at a considerable

distance from the front. This is important because the pseudospectral interpolation, which

for these ¬gures was performed without smoothing or special tricks, becomes increasingly

inaccurate as the front is approached because of the oscillations of Gibbs™ phenomenon.

When the time step is small, some particles originate very close to the front, so the inter-

polation gives poor approximations to the velocities that advect the particles towards the

front from either side. The result is large errors.

SL schemes that use various low order polynomial interpolations suffer the same curse;

close proximity to a jump discontinuity is an ill wind for all approximation schemes. Nu-

merical experiments such as Kuo and Williams (1990) have con¬rmed that SL methods

have their smallest errors as a function of time step for rather large time steps ” a couple

of hours in meteorological applications.

Fig. 14.7 illustrates another odd property of SL schemes. If one plots only the grid point

values of a pseudospectral SL scheme, one ¬nds, for moderate „ , a good approximation to

the exact solution. However, the Fourier series can be summed to evaluate u(x) at arbitrary

points not on the canonical grid. If we do this for many intermediate points, we obtain the

very oscillatory graph shown in the right panel. The pseudospectral solution is much better

on the grid than off.

Finite element methods often exhibit a similar “superconvergence” at the collocation

2 2

1

1

0

0

1.2

1.2 0.6 0.8 1.0

0.8

0.6 1.0

x

x

Figure 14.7: Same case as the two previous ¬gures, but a “zoom” which shows only part of

the x-interval around the shock. In both panels, the dashed line is the exact solution. On

the left, the solid line-with-circles shows the semi-Lagrangian Fourier solution at the points

of the collocation grid. The solid-plus-crosses on the right is the same solution except that

the Fourier interpolant has been summed at many points between those of the collocation

grid. The semi-Lagrangian solution exhibits “overconvergence”, that is, the solution is

much more accurate AT the collocation points than BETWEEN them.

CHAPTER 14. SEMI-LAGRANGIAN ADVECTION

280

points. The existence of this phenomenon has been rigorously proved for some ¬nite ele-

ment algorithms, but not as yet for spectral.

14.8 Two-Level SL/SI Algorithms

The three-level SL/SI is simple and robust, but it has drawbacks. The ¬rst is that it re-

quires storing the solution at three different levels. The second is that the implicit terms

are treated with an effective time step of 2„ whereas a standard Crank-Nicholson scheme

would average the unknowns at time levels tn+1 and tn with a time step of „ and an error

four times smaller.

Fortunately, it is possible to replace Robert™s scheme with a two-level method that has

a much smaller proportionality constant. The key idea is to apply the three-level formula

with u’ interpreted as the solution at time level tn rather than tn’1 . The intermediate

velocity at t = (n + 1/2)„ is obtained by extrapolation from previously known values of

the velocity.

The displacements are calculated from

±j = („ /2) uext (xj ’ ±j , tn+1/2 ) (14.45)

where ± is de¬ned to be the displacement over a time interval of „ /2, subscript “ext” de-

notes a quantity which must be extrapolated and j is the index of (spatial) grid points.

The actual time-marching step for Du/Dt + G(u, t) = F (u, t) is identical in form to the

three-level scheme:

u+ + („ /2)G+ = u’ ’ („ /2) G’ + („ /2) Fext

0

(14.46)

As before, u+ = u(x, tn + „ ), but

≡ F (x ’ ±, tn+1/2 )

0

Fext

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

The most widely used extrapolation formula is

15 10 3

u(x, tn ) ’ u(x, tn’1 ) + u(x, tn’2 )

uext (x, tn+1/2 ) = (14.48)

8 8 8

which Temperton and Staniforth (1987) found preferable to the two-level extrapolation

uext (x, tn+1/2 ) = (3/2)u(x, tn ) ’ (1/2)u(x, tn’1 ) (14.49)

It is possible to extrapolate along trajectories (that is, shifting x positions by the displace-

ment ± as one moves back in time), but Temperton and Staniforth (1987) found it was more

accurate (and cheaper) to extrapolate along grid points.

Extrapolation has a bad effect on stability, as noted above. However, the two-level

method is stable if the extrapolation is done carefully. With complicated physics ” moun-

tains, complicated temperature dependence and so on ” it is sometimes necessary to add

extra damping, such as replacing Crank-Nicholson by the “theta” scheme (which is a lin-

ear combination of Crank-Nicholson with Backwards-Euler that weakly damps high fre-

quency). It is also a good idea to transfer as many terms as possible from the extrapolated,

0

explicit quantity Fext to the implicitly-treated term G. These empirical ¬xes are problem-

speci¬c and model-speci¬c, so we must defer the details to the papers catalogued in the

bibliographic table below.

14.9. NONINTERPOLATING SL & NUMERICAL DIFFUSION 281

The two-level scheme is in theory twice as cheap as the three-level scheme in the sense

that one can achieve the same level of time accuracy with twice as long a time step in the

two-level model as in the older algorithm. However, because the three-level scheme is

more stable and does not need time-extrapolation, it remains competitive with two-level

methods.

14.9 Noninterpolating SL Methods and Computational Dif-

fusion

The error in interpolating from the evenly spaced grid to the “feet” of trajectories is a source

of computational diffusion. McCalpin (1988) has given a thorough analysis for polynomial

interpolation schemes of various orders, but it would take us too far a¬eld to reproduce his

analysis here. Instead, we shall offer only a heuristic argument.

Consider the simplest case: the pure advection equation in one dimension. The formula

for computing the new unknown in the three-level SL scheme is then

u(x, tn+1 ) = u(x ’ 2±, tn’1 ) (14.50)

Interpolation errors will add a small random component to the displacement ±. Sometimes

the ¬‚uid blob will be advected a little too far, sometimes not far enough, because interpo-

lation errors generate inexact particle velocities and displacements. To the SL scheme, the

foot of a trajectory is not known exactly, but only to within a line segment (or circle or

sphere) which is the size of the average interpolation error. The advected ¬elds are blurred

by this randomness, and blurring is equivalent to a computational damping.

McCalpin™s analysis shows linear interpolation generates ordinary, second derivative

diffusion; cubic interpolation gives a biharmonic (fourth derivative) damping and so on.

For ocean models where the scale of the energy-containing eddies is an order of magni-

tude smaller than in the atmosphere, and the affordable number of grid points-per-eddy-

diameter correspondingly smaller, too, McCalpin argues that, arithmurgically speaking,

this computational diffusion would turn the seas into great viscous pools of honey.

Fig. 14.8 illustrates his point. The semi-Lagrangian solution is damped much more than

that of the standard pseudospectral method even though the explicit, physical viscosity is

the same.

Ritchie (1986) was thus motivated to invent a “non-interpolating” SL algorithm. The

basic idea is very simple. At large Courant number, split the advecting velocity into two

parts:

(14.51)

uj = uj,SL + uj,E

where uj,SL is chosen so that the “foot” of the trajectory is exactly on a point of the regular

grid when the particle is advected by uj,SL only. Of course, the actual foot of the trajectory

almost never coincides with a point of the regular grid, but this is not a serious limitation

since the SL algorithm is really just a transformation to a moving coordinate. There is

no law of nature or arithmurgy that says that the velocity of the coordinate change must

exactly match that of the ¬‚uid.

However, Ritchie™s modi¬ed transformation eliminates only part of the advection, so

the remainder must be treated by the standard pseudospectral method, just as if we used

an Eulerian (¬xed grid/trajectory-ignoring) coordinate system throughout. It is important

to choose uj,SL to advect to the grid point nearest the actual foot of the trajectory so that

uj,E is as small as possible. The residual Courant number is then always less than one so

CHAPTER 14. SEMI-LAGRANGIAN ADVECTION

282

Viscosity=h/2

1 1

0.9

0.5

0.8

0

0.7

-0.5

0.6

-1 0.5

-2 0 2 2 2.5 3

x x

Figure 14.8: Comparison of the semi-Lagrangian/semi-implicit Fourier method (thick

lower curve) with the standard Fourier pseudospectra method (upper curve with circles)

for Burgers™ equation, ut + uuxx = ν uxx where ν = h/2 and h = π/32 is the grid spacing.

The SL time step was „ = 1/2; that for the standard method was determined by an adap-

tive 4th/5th order Runge-Kutta subroutine. The solution from u(x, 0) = sin(x) “breaks”

at t = 1, that is, would develop an in¬nite slope at that time were it not for the viscosity.

The solutions are shown at t = 2. In spite of the long time step, which tends to slow the

diffusion of the νuxx term (which is treated implicitly), the SL/SI algorithm is much more

strongly damped than the standard pseudospectral method, which treats the viscosity hon-

estly by using a short time step. The SL/SI damping comes primarily from interpolation

errors.

that the integration is stable with very long time steps where the residual Courant number

is

„

CoE ≡ max(|uE |) (14.52)

hmin

where hmin is the smallest distance between any two points on the regular pseudospectral

grid.

Ritchie™s later papers (1987, 1988, 1991) show that the “non-interpolating” SL method

works just ¬ne with both spectral and ¬nite difference spatial discretizations, even in three-

dimensional models. However, his idea has only partially replaced older SL algorithms.

There are two reasons. First, treating some of the advection by a non-SL scheme brings

back some of the disadvantages and poor front resolution of conventional methods. Sec-

ond, for atmospheric forecasting models, the resolution is suf¬ciently high and the length

of the integration is suf¬ciently short that damping is not a problem. However, climate

models, which are run at lower spatial resolution for long time integrations, may bene¬t

from non-interpolating SL.

For example, his comparisons with non-SL models showed that at a resolution of “T106”,

then the standard for forecasting, both interpolating and non-interpolating SL models lost

only about 1% of the initial energy. This was actually better than the standard non-SL

14.10. OFF-GRID INTERPOLATION 283

spherical harmonics model of the European Centre for Medium-Range Weather Forecast-

ing, which lost about 3% of its energy due to an (unphysical) biharmonic damping which

is included to suppress small scale noise and give better front resolution. The SL models,

in spite of their inherent implicit computational dissipation, were actually less damping

than the non-SL algorithm because an explicit arti¬cial viscosity was unnecessary for the

SL schemes.

At lower resolution, as in a climate model or an ocean model, Robert™s SL algorithm

is more disspative because the interpolation errors worsen as the grid coarsens. At T63,

Ritchie found that his non-interpolating model lost only 1% energy, but the interpolating

SL lost 5%.

Smolarkiewicz and Rasch (1990) have generalized non-interpolating SL schemes. They

show that their generalization can be applied to many existing, non-SL algorithms to relax

time step limitations while retaining good properties of the methods such as that of being

“positive-de¬nite monotone”.

The conclusion is that non-interpolating SL may be very useful for ocean and climate

modelling, or whenever the resolution is suf¬ciently poor that interpolation errors create

too much damping. However, the errors are a function of the order of the interpolation as

well as the grid spacing. If one uses high order interpolation, such as fast spectral off-grid

interpolation, or even quintic polynomial interpolation, instead of the cubic interpolation

which is now the standard in atmospheric models, the computational diffusion is greatly

reduced. Thus, it is only in the context of speci¬c problems with speci¬c grid size and a

particular interpolation order that one can intelligently answer the question: Are interpo-

lating or non-interpolating schemes better?

14.10 Off-Grid Interpolation: Shape-Preserving, Conserva-

tion and All That

14.10.1 Off-Grid Interpolation: Generalities

Off-grid interpolation is rather expensive, if done in the obvious way. For spectral methods,

off-grid interpolation cannot be done by the Fast Fourier Transform because the departure

points are spaced irregularly between the points of the canonical grid. Summing an N -term

spectral series at each of N points costs O(N 2 ) operations.

There are two subproblems which require off-grid interpolation: (i) calculating the dis-

placements and (ii) calculating properties at departure points. The interpolation of mid-

trajectory velocities can be done by a lower order method than interpolating velocities,

temperatures, ozone concentration and so on at the departure points. The reason is that

the displacement is the velocity multiplied by the time step, and so errors in the displace-

ment are smaller than those in u itself.

It is much more important to interpolate accurately when calculating the properties of

a ¬‚uid blob at its departure point. The reason is that errors in this step are not multiplied

by the time step, but rather are errors in the answer itself. An error of a tenth of a degree

in interpolating the temperature at the departure point will be carried through the entire

remainder of the computation.

The third comment is that SL schemes do not automatically conserve either mass or

energy, nor preserve the positive de¬niteness of densities like mass and water vapor which

are always non-negative. It is possible to modify SL procedures to recover many of these