niteness” wars must be fought.

CHAPTER 14. SEMI-LAGRANGIAN ADVECTION

284

14.10.2 Spectral Off-grid

Spectrally-accurate off-grid interpolation algorithms include the following:

1. Fast Multipole Method (FMM) (Boyd, 1992c, Dutt & Rohklin, 1993, 1994)

¨

2. Partial Chebyshev expansion method (Suli & Ware, 1991, 1992, Ware, 1991,1994)

3. Fine-grid Lagrange interpolation algorithm (Boyd, 1992)

As explained at greater length in Chapter 10, all three methods are asymptotically faster

than direct summation. (See the review by Ware, 1998). However, the proportionality

factors, which vary with the desired error tolerance, are always much greater than for the

Fast Fourier Transform.

¨

The Suli-Ware algorithm is the only one that has been explicitly applied in semi-Lagrangian

calculations. It performed well in both one and two space dimensions.

Boyd™s ¬ne-grid polynomial interpolation algorithm is probably both the fastest and

the simplest to explain. The key idea is that low order polynomial interpolation can be

very accurate if performed on a suf¬ciently ¬ne grid. So, his ¬rst step is to apply the FFT

to transform the data to a ¬ner grid than used for the rest of the calculations. One can then

use ordinary Lagrange interpolation or splines on the ¬ne grid with little loss of accuracy

compared to spectral off-grid interpolation on the original coarser grid.

His algorithm has not been explicitly used in semi-Lagrangian models to date. However,

the ECMWF forecasting model and the NCAR Community Climate Model both employ

his idea implicitly. For reasons of cost, and because fast off-grid interpolations were not yet

known at the time, the programmers used various forms of low-order interpolation and

hoped for the best. (Careful testing showed that this inconsistent use of low order off-grid

interpolation in a high order model does not signi¬cantly degrade the overall forecast.)

However, these models employ dealiasing so that inner products are evaluated on a grid

which is much ¬ner than the actual number of spectral degrees of freedom. Thus, these

models implicitly use Boyd™s strategy.

The inconsistency between spectral calculation of derivatives and non-spectral off-grid

interpolation is in any event not as illogical as it seems. The reason is that to minimize

dispersion, it is important to have accurate computation of derivatives in the linear terms

which are not treated by semi-Lagrangian advection. The off-grid interpolation does not

have to approximate derivatives, but only the velocities themselves, which is easier.

There are two open questions about spectral off-grid interpolation. First, is it cost-

effective in a real climate or forecasting model, or is low order polynomial interpolation

suf¬cient? The second is: can spectral semi-Lagrangian schemes be improved, when nar-

row fronts form, by skillful use of ¬lters or sum-acceleration methods?

We discuss ¬lters brie¬‚y in Chapter 18, Sec. 20, but it is useful to at least hint at the

possibilities. In Fig. 14.9, there is a large spike of error in the immediate neighborhood

of the front in both panels; it is very dif¬cult to eliminate this because even a slight er-

ror (O(h)) in the location of the front will inevitably generate O(1) errors, although the

shape of the computed ¬‚ow may look very similar to the true ¬‚ow. Away from this small

neighborhood, however, Fig. 14.9 shows that ¬ltering can dramatically reduce errors in a

Fourier semi-Lagrangian method. There is only a little experience with ¬ltered spectral

semi-Lagrangian schemes, and ¬lters themselves are a hot research topic. The graph sug-

gests that further work on ¬ltered semi-Lagrangian algorithms would be rewarding.

14.10.3 Low-order Polynomial Interpolation

Finite difference and low order ¬nite element models commonly use polynomial interpo-

lation to approximate quantities off-grid. Cubic splines (Purnell, 1976) have been tried, but

14.10. OFF-GRID INTERPOLATION 285

. Filtered

1 .1 Unfiltered

.01 .01

.001 .001

1E-4

1E-4

π -π π

-π 0

0

Figure 14.9: Errors in the solution to the one-dimensional advection equation by Fourier

three-level semi-Lagrangian computation with N = 64 grid points and a time step „ = 1/2.

The initial condition is u(x, 0) = 1 ’ sin(x). Left: Vandeven ¬lter of order p = 6. Right:

un¬ltered numerical solution.

the most common choice for interpolating particle properties (i. e., u’ in the SL schemes

above) is cubic Lagrange interpolation. (Leslie and Purser, 1991, used quintic interpola-

tion.) There is less consistency in calculating displacements because as noted in the pre-

vious subsection, the forecast error is much less sensitive to errors in the velocity in the

¬xed-point iteration for ± than to errors in u’ , F 0 , etc. Some therefore use linear interpo-

lation of u(x ’ ±, tn ) in calculating displacements. McDonald (1986) suggested a hybrid

strategy: linear interpolation on the ¬rst ¬xed-point iteration followed by cubic interpola-

tion on the second (and usually last) iteration.

Moisture has a very patchy distribution in the atmosphere; the air above the Mediter-

ranean Sea is very humid because the warm water evaporates like crazy, but the Mediter-

ranean Sea is bordered by the Sahara desert, which is not humid at all. Global spectral

models, with a resolution of only O(100 km), exhibit Gibbs™ phenomenon in the vicinity

of such moisture fronts. The resulting regions of negative water vapor are decidely un-

physical and have a bad effect on model climate. SL schemes do much better. (Rasch and

Williamson,1990b,1991 and Williamson, 1990).

The moisture transport problem has motivated studies which replace ordinary inter-

polation by a procedure in which the interpolant is adjusted by side constraints, such as

smoothness and positive-de¬niteness. Such “quasi-monotone” or “shape-preserving” SL

schemes have been investigated by Williamson and Rasch(1989), Rasch and Williamson(1990a)

and Bermejo and Staniforth (1992).

Adjusted or constrained interpolation can also be used to enforce conservation of total

mass. Priestley (1993) and Gravel and Staniforth (1994) have generalized SL to be quasi-

conservative.

CHAPTER 14. SEMI-LAGRANGIAN ADVECTION

286

14.10.4 McGregor™s Taylor Series Scheme: an Alternative for Computing

Departure Points

The only disadvantage of polynomial interpolation is that if done accurately, it is rather

expensive, especially in multiple space dimensions. McGregor (1993) suggested an alter-

native procedure, which is to calculate the trajectory by a Taylor series in time about the

departure point. The catch is that because the trajectory is the goal, the time derivatives

must be total time derivatives. Thus, in one space dimension,

„2 „3 ‚ ‚

± = ’„ u + (ut + uux ) ’ (14.53)

+u (ut + uux ) + . . .

2 6 ‚t ‚x

Calculating the local time derivatives above is very hard. However, McGregor argues

that these can be neglected compared to the advective terms so that

„2 „3 ‚

± ≈ ’„ u + uux ’ u (14.54)

(uux ) + . . .

2 6 ‚x

where the velocities are to be evaluated at the spatial location of the arrival point at the

midpoint in time (tn for the three-level scheme, tn+1/2 for a two-level method).

His method is an order of magnitude cheaper than bicubic interpolation and accurate to

second order. Adding additional terms does not increase the formal time accuracy beyond

second order, but McGregor shows that empirically, including partial third order terms

does improve accuracy considerably for ¬nite „ .

14.11 Higher Order SL Methods

Purser and Leslie (1994) have developed an ef¬cient third-order SL scheme. Earlier, Maday,

¨

Patera and Rønquist (1990) and Suli and Ware (1991, 1992) and Ware (1991, 1994) have

developed a family of algorithms which can be extended to arbitrary order. Building on

Holly and Preismann (1977) and Yang and Hsu (1991), Yang, Belleudy and Temperville

(1991) developed a ¬fth order accurate SL scheme.

The bad news is that these high order schemes are all rather complicated. The assump-

tion that a particle is advected at constant speed by the velocity at the midpoint of the

trajectory is suf¬cient to give second order accuracy in time. At higher order, the displace-

ment on t ∈ [tn’1 , tn ] is no longer the same as t ∈ [tn , tn+1 ], and a simple ¬xed point

iteration no longer gives direct access to the trajectory.

¨

The Maday-Patera-Rønquist and Suli-Ware algorithms require solving auxiliary time-

marching problems, k problems for a (k + 1)-th order method. Maday et al. show that

one can “sub-cycle”, integrating these auxiliary advection problems with a short time step

while applying the implicit time-marching formulas (for non-advective terms) only over a

long time interval.

Even so, these algorithms are uncompetitive with the second order SL schemes for me-

teorological applications. We shall therefore not describe them in detail. These are likely

to become more important as resolutions improve so that time-truncation errors begin to

catch up with space truncation errors. The latter are currently much worse for forecasting

and climate models than the temporal differencing errors.

14.12 History and Relationships to Other Methods

14.12. HISTORY AND RELATIONSHIPS TO OTHER METHODS 287

Table 14.1: A Selected Bibliography of Semi-Lagrangian Methods

SL is an abbreviation for “semi-Lagrangian”; SI for “semi-implicit”.

Reference Comments

Fjørtoft(1952, 1955) First Lagrangian algorithm (meteorology)

Wiin-Nielsen (1959) Earliest SL, at least in meteorology

Sawyer (1963) coined “semi-Lagrangian”

Holly&Preissman (1977) SL variant using Hermite interpolation (hydrology)

Haltiner&Williams(1980) Good review of (singly-upstream) SL

Robert (1981, 1982) Earliest pairing of SL with SI

Douglas&Russell(1982) Invention of an SL scheme

Pironneau (1982) Invention of SL as “transport-diffusion” algorithm

Bates&McDonald (1982) First multiply-upstream SL

Bates (1984) More ef¬cient multiply-upsteam SL

Morton (1985) Invention of SL as “Lagrange-Galerkin” scheme

Robert&Yee&Ritchie(1985) SL with SI

Pudykiewicz et al.(1985) Convergence of ¬xed point iteration

Temperton & Staniforth(1987) Independent invention of two-level SL

McDonald & Bates (1987) 2d order accuracy through extrapolated velocities

Ritchie (1986) Invention of “non-interpolating” SL

Ritchie (1987,1988,1991) Non-interpolating SL applied to a sequence of

increasingly sophisticated spectral models

McCalpin (1988) Analysis of dissipation inherent

in SL interpolation

McDonald&Bates(1989) Showed SL with time-splitting was inaccurate

Bates&Semazzi& Finite difference SL shallow water model

Higgins&Barros(1990) Accuracy improved by avoiding time-splitting

ˆe

Cot´ &Gravel&Staniforth(1990) Finite-element SL with “pseudostaggering”

Williamson&Rasch(1989) SL with shape-preserving (positive de¬nite)

Rasch&Williamson(1990a) advection

Rasch&Williamson(1990b,1991) moisture transport in climate & weather models

Williamson(1990)

Purser&Leslie(1988,1991,1994) Not spectral but improvements in SL

Leslie&Purser(1991) latest article is 3rd order in time

Yang&Belleudy Fifth order SL

&Temperville(1991)

Yang&Hsu(1991) Improved Holly-Preismmann & Yanget al. schemes

Smolarkiewicz & Rasch(1990) Generalization of non-interpolating SL; monotone

Maday&Patera&Rønquist(1990) SL of arbitrarily high time-order

¨

Suli&Ware(1991,1992) Fourier and Chebyshev SL/SI, proofs and examples

Ware (1991,1994) independent invention of arbitrary time-order SL

¨

Baker&Suli&Ware(1992) & fast off-grid interpolation scheme

ˆe

Staniforth& Cot´ (1991) Highly readable review

ˆe

Gravel&Staniforth&Cot´ (1993) Stability analysis

CHAPTER 14. SEMI-LAGRANGIAN ADVECTION

288

Table 14.1: Bibliography of Semi-Lagrangian Methods (continued)

Reference Comments

Bermejo&Staniforth(1992) Quasi-monotone SL

Priestley(1993) Quasi-conservative SL

Gravel&Staniforth(1994) Quasi-conservative SL for shallow-water system

Inexpensive O(„ 2 ) formula for departure pts.

McGregor(1993)

Pellerin&Laprise Experiments with simple advection/

&Zawadzki(1995) condensation system

ˆe

Cot´ &Gravel Extra time level scheme to defeat

&Staniforth(1995) orographic resonance

Ritchie&Tanguay(1996) Comparison of spatially averaged Eulerian

& Lagrangian treatments of mountains

¨

Bottcher (1996) Modi¬ed exponential splines for interpolation

Makar&Karpik(1996) Periodic B-spline interpolation on sphere for SL;

comparisons with spectral methods

The Lagrangian (particle-tracking) formalism is a couple of centuries old. Fjørtoft (1952,

1955) showed that particle-tracking could be used for weather forecasting. However, We-

lander (1955) showed through a vivid and widely- reproduced graph that blobs of ¬‚uid

were rapidly stretched and deformed into a very convoluted shapes. Today, we would

say that ¬‚uids can be chaotically mixed, in the sense of “chaos” theory, and a blob of ¬‚uid

initially bounded by a small square is bounded, for later times, by a “fractal” curve of

ever-increasing length.

Wiin-Nielsen (1959) therefore introduced a “semi-Lagrangian” scheme: just as in mod-

ern schemes, each set of particles is tracked only over a single small time interval of length

„ . A new set of particles is chosen for later time steps: those which will arrive at the

points of an evenly spaced grid. Later meteorological papers by Krishnamurthi (1962) and

Sawyer (1963) and others established that SL algorithms are useful, but did not identify

overwhelming advantages.

The method become much more popular when Robert (1981, 1982) showed how to

combine semi-Lagrangian advection with semi-implicit treatment of the other terms. This

launched an explosion of activity in meteorology which has continued to the present.

Bates and McDonald (1982) made another important advance by introducing “multiply-

upstream” CFL methods. These allowed the use of a time step which far exceeded the

advective limits of conventional explicit schemes.

Temperton and Staniforth (1987) and McDonald and Bates (1987) independently de-

vised two-level schemes of second order time accuracy. (Robert™s method required three

time levels while the earlier two-level algorithms of Bates and McDonald were only ¬rst or-

der.) These allow a doubling of the time step with no loss of accuracy. However, two-level

schemes are somewhat less stable than Robert™s method, so several later articles have de-

scribed slight modi¬cations, such as a “theta” treatment of the implicit schemes, a damping

of divergence, and so on to render such algorithms more stable and robust.

Ritchie (1986) developed the “non-interpolating” version of SL, and then applied it to

a hierarchy of increasingly complicated models in later papers. Smolarkiewicz and Rasch