This procedure is also called “padding” because one must add (N/2) zeros to the spectral coef¬-

cients before making the coef¬cient-to-grid Fast Fourier Transform.

Whatever the name, the ¬ltering or padding is always performed just before taking coef¬cient-

to-grid transforms so that the shortest un¬ltered waves on the grid have wavelength 3h.

A total of (3/2) N basis functions are used during intermediate stages of the calculation,

but only the lowest N wavenumbers are retained in the ¬nal answer at each time level.

Two alternative deliasing methods have also been tried. Patterson & Orszag (1971)

evaluated grid point sums with a phase shift. However, this trick has fallen from favor

because it is always more expensive than the Two-Thirds Rule (Canuto et al. 1988, pg. 85).

Rogallo (1977, 1981, Rogallo and Moin, 1984), who is better known for inventing the

wing used in the ¬rst modern hang-glider, showed that for a two-step time scheme, it is

possible to apply the phase shift method at negligible extra cost to reduce aliasing errors

11.6. ENERGY-CONSERVING: CONSTRAINED INTERPOLATION 213

to O([ t]2 ) of those of the standard pseudospectral method. This is very appealing, even

though a small residue of aliasing remains, but for a decade, it was never mentioned except

in Rogallo™s own work until praised by Canuto et al (1988).

Other strategies which control aliasing indirectly are described in the next section.

11.6 Energy-Conserving Schemes and

Skew-Symmetric Advection:

Interpolation with Constraints

Arakawa (1966) proposed a very different strategy: modifying the numerical algorithm so

that it exactly conserves the discretized energy. (By “energy”, we mean the total energy

of the model.) He reasoned that blow-up implies a spectacular (and wholly spurious) in-

crease in energy. If the discretized energy is constant, however, then supersonic winds are

obviously impossible.

Similarly, for Hamiltonian systems, one can use special time-marching schemes (“sym-

plectic” algorithms) which preserve the Hamiltonian structure even in discrete form. These

have been very useful in megayear-long integrations of solar system orbital dynamics

problems (Sanz-Serna and Calvo, 1994).

As noted by Canuto et al.(1988), it is often straightforward to obtain spectral schemes

which are energy-conserving. Indeed, Galerkin™s method is automatically energy-conserving

for some equations, such as for the barotropic vorticity equation with a triangular trunca-

tion of the spherical harmonic basis (Haltiner and Williams, 1980). However, for most

problems, pseudospectral schemes conserve the discrete energy only after (small) modi¬-

cations (Canuto et al., 1988, Secs. 4.5 and 4.6).

Zang (1991a) has shown that the most important such modi¬cation is to use the skew-

symmetric form of advection for the Navier-Stokes equations. A skew-symmetric matrix

or differential operator is one whose eigenvalues are all pure imaginary. Skew-symmetry

is a property of advection operators and also the operators of wave equations: the imagi-

nary eigenvalues (paired with a ¬rst order time derivative) imply that operator advects or

propagates without causing growth or decay. It is desirable to preserve skew-symmetry in

numerical approximation because pure advection and propagation is precisely what doesn™t

happen when a model is blowing up.

The usual forms of advection are

u · u, [Standard form]

ω—u+ ((1/2) u · u) , [Rotation form] (11.14)

However, Zang (1991a) shows that skew-symmetry is preserved only by the more compli-

cated form

1

1

u· [Skew-Symmetric Advection] (11.15)

u+ (u u) ,

2 2

where (u u) is a vector whose components are · (u uj ) where uj is one of the three

vector components of u. The property of skew-symmetry also depends upon boundary

conditions, but we refer the reader to his paper for details.

Blaisdell, Spyropoulos and Qin(1996), through a mixture of theoretical analysis and

experiments for Fourier pseudospectral codes for turbulence and Burgers™ equation, offer

additional support for the skew-symmetric form. They show that this form cancels much

of the differentiation error of either the conservative form or the usual non-conservative

advection form for the largest aliased components.

CHAPTER 11. ALIASING, SPECTRAL BLOCKING, & BLOW-UP

214

More recently, a wide variety of algorithms have been developed to avoid spurious

oscillations near fronts and shock waves, such as Van Leer™s MUSCL scheme, Roe™s algo-

rithm, and the Piecewise-Parablic Method (Carpenter et al., 1990). These do not attempt

to control the energy, but rather use analytical theories of shock structure to suppress the

growth of 2h-waves.

Some of these schemes have the disadvantage of being more costly than conventional

schemes. Arakawa™s ¬nite difference Jacobian requires twice as many operations as ap-

proximating the advection operation using standard differences. The Two-Thirds Rule re-

quires dumping (1/3) of the numerical degrees of freedom in each dimension, doing the

work of a grid which is say, 963 to obtain only a 643 set of nonzero spectral coef¬cients.

Is the additional cost worth it? The answer, to borrow the title of a science ¬ction novel

by the Strugasky brothers, is De¬nitely Maybe.

Optimist™s Viewpoint: Energy-conserving and front-resolving schemes are interpolation-

with-physics-constraints. Instead of merely using the grid point values of a function, as if

it were completely arbitrary, these schemes exploit the fact that the function is not arbi-

trary, but rather is constrained by physics. Building constraints into interpolation schemes

is perfectly legitimate and should lead to an improved approximation.

Pessimist™s Viewpoint: The Arakawa scheme guarantees that the numerical solution

will stay bounded, but it does not guarantee that it is accurate. In Phillips (1956), for example,

the calculation blew up only at the stage where the baroclinic cells were forming fronts.

Since his coarse grid could not have resolved the fronts anyway, it can be argued that the

blow-up of his model really did him a favor by aborting the computation as soon as it

became inaccurate.

Similarly, dealiasing is a dumb idea: if the nonlinear interactions are so strong that a lot

of energy is being aliased, one probably needs to double the grid spacing.

Likewise, the PPM method and its cousins guarantee smoothness but not accuracy.

Who is right?

11.7 Energy-Conserving and Shock-Resolving Schemes: Dis-

cussion

The workshop published by Straka et al.(1993) computed a two-dimensional density cur-

rent using a wide variety of schemes. The comparison supports both the optimistic and

pessimistic viewpoints.

Fig. 11.8 compares a standard, second order ¬nite difference method at various resolu-

tions. The most important point is that one needs a grid spacing h = 100m or less to get a

decent solution; the plot for h = 200m is ¬lled with noise while the numerical solution for

h = 400m is completely useless.

Fig. 11.9 compares the high accuracy solution (the h = 25m solution from the previous

¬gure, which Straka et al. use as a benchmark) with standard ¬nite differences and the PPM

scheme at h = 400. The optimist is encouraged because the PPM solution is still stable,

smooth, and resembles the accurate solution even though the grid spacing is very coarse.

The PPM™s incorporation of wave physics has made it far more successful at this resolution

than the ¬nite difference scheme, which is junk. The spectral model was unstable at this

resolution.

However, the pessimist can ¬nd solace, too, because the smoothness and stability of the

PPM solution is deceiving: it is not a very good numerical solution. Two smaller vortices

and several other features in the right half of the disturbance have been simply combined

into a single smooth vortex by PPM. The truth is that h = 400m is simply too coarse to

resolve the smaller features.

11.7. ENERGY-CONSERVING SCHEMES: DISCUSSION 215

h=25m h=133m

h=33m h=200m

h=50m h=266m

h=66m h=400m

h=100m h=533m

Figure 11.8: Comparison of second order ¬nite difference solutions at different grid spac-

ings h. The algorithm has no particular front-resolving or monotonicity properties. From

Straka et al.(1993) with permission of the authors and the American Meteorological Society.

At h = 200m, the spectral solution is noisy, but much more accurate than any of the

others in computing the total energy and enstrophy, even though it is not exactly energy-

conserving. PPM is smoother, but only because this algorithm has strong computational

dissipation, losing almost half the enstrophy at the time shown in the ¬gures versus only

about 13% for the spectral code. At h = 100m, the spectral and PPM schemes are both

accurate, but the spectral method is considerably better than PPM at resolving the smaller

scale features and preserving energy and enstropy (Fig. 11.10).

The moral of this and similar comparisons is the following

Rule-of-Thumb 11 (DEALIASING/ENERGY-CONSERVING)

1. For well-resolved ¬‚ows, dealiasing and energy-conserving and front-resolving algorithms are

unnecessary. A spectral method will conserve energy to a very high accuracy and will faith-

fully resolve a shock zone or frontal zone with only two or three grid points within the zone.

2. For marginally-resolved ¬‚ows, dealiasing and other remedies may prevent blow-up and yield

acceptable results when blind application of spectral or other methods yields only a sea of

noise or worse. However, some ¬ne structure in the solution may be missed, and the special

algorithms can do nothing about this.

3. For poorly-resolved ¬‚ows, nothing helps, really, except using more grid points. Energy-

conserving and front-tracking schemes may generate smooth solutions, but these will only

vaguely resemble the true solution.

Special algorithms and tricks are “Plan B”, fallbacks to implement when the straightforward

algorithm isn™t working well and the computation already ¬lls all the available computer memory.

Filtering and skew-symmetric advection are often useful; conserving energy is helpful; dealiasing is

a last resort because of its high cost in two or three dimensions. Front-resolving schemes are quite

useful for shocks.

CHAPTER 11. ALIASING, SPECTRAL BLOCKING, & BLOW-UP

216

PPM: h=400 m

Accurate Soln.

2d order FD: h=400 m

Figure 11.9: Comparison of the piecewise-parabolic method (PPM) at h = 400 [top] with

the standard second order ¬nite difference model [bottom] at the same resolution and with

the highly accurate reference solution [middle]. Although the PPM solution is smooth and

has captured some of the features of the true solution, it has only two vortices instead of

three. Redrawn from Straka et al.(1993).

A turbulence calculation almost by de¬nition is poorly resolved, so dealiasing is not

particularly rare for modelling such ¬‚ows. Even in turbulence, however, dealiasing is often

unnecessary; even in turbulence, dealiasing is not a cure for a ¬‚ow that is badly underre-

solved.

Aliasing error is likely to be important only when, as in Phillips™ pioneering model, the

calculation is becoming highly inaccurate. Aliasing blow-up has saved many modellers

from believing lousy simulations. There have been many energy-conserving calculations

published ” it would be tactless to name names! ” which have smooth solutions that are

complete nonsense.

11.8 Aliasing Instability: Theory

A number of papers have tried to identify speci¬c mechanics for blow-up such as Briggs,

Newell and Sarie (1981). Aoyagi (1995) is a recent article with many references. These stud-

ies have successfully identi¬ed particular modes of instability for particular algorithms.

However, the studies have not been successful in developing any sort of general theory.

It may be that model blow-up, like any other train wreck, creates casualties of many sorts

and not any speci¬c, aliasing-related injury.

Phillips (1959), who explained the blow-up of his model in terms of aliasing and the

resulting spurious, wrong-way transfer of energy, believed that he had diagnosed a speci¬c

illness with a speci¬c cure. He did not think that aliasing instability was merely a lack of

resolution because he repeated his GCM calculation with a much shorter time step and

spatial grid spacing and the model still blew up at about the same time!

11.8. ALIASING INSTABILITY: THEORY 217

SPECTRAL:

h=100 m

ACCURATE

SOLN.

PPM:

h=100 m

Figure 11.10: Comparison of spectral and piecewise-parabolic method (PPM) at h = 100m

with each other and the reference solution. From Strata et al.(1993).

However, his experiments do not completely resolve the issue. When a shock wave

forms, the solution is smooth for a ¬nite time interval and then develops a jump discon-

tinuity which ” in the absence of suf¬cient viscosity ” the numerical method cannot re-

solve. The transition from exponential accuracy to blow-up may be rapid, occurring on

a time scale short compared to the total time of integration. Even with an abrupt onset,

blow-up may be simply the result of underresolution of small scale features. 1 As we saw

earlier, spectral blocking can also be the result of too long a time step (although Phillips™

tests excluded this in his model.)

Spectral blocking would seem to support Phillips™ belief that aliasing instability is not

merely underresoution. If the ¬‚ow is evolving narrower and narrower fronts, one would

expect that the Fourier spectrum would ¬‚atten out, rather than develop a curl at high

wavenumber. Actually, in simple one-dimensional PDEs like the Regularized Long Wave

equation and Burgers™ equation, this “¬‚atten-to-blow-up” scenario is observed. On a log-

linear plot, the coef¬cients asymptote to a straight line for all t, implying that the conver-

gence rate is always geometric, but the slope of the asymptotic line becomes ¬‚atter and

¬‚atter. When the plotted coef¬cients asymptote to a horizontal line, i. e., the spectrum

resembles “white noise”, disaster is quick and spectacular.

In Phillips™ simulation, and also in the Zang-Krist-Hussaini transition-to-turbulence ex-

periment illustrated earlier, however, spectral blocking develops while the large scale fea-

tures are still well-resolved and the amplitude of the “blocked” wavenumbers is small

compared to that of low wavenumbers. Doesn™t this imply a numerical instability rather

than underresolution as the cause of blocking? Not necessarily.

Theories of two-dimensional and quasi-geostrophic turbulence predict that (i) the ¬‚ow

1 Cloot, Herbst and Weideman, 1990, give a good discussion of a one-dimensional wave equation that develops

unbounded solutions in ¬nite time. Their adaptive Fourier pseudospectral method increases N close to the time

of the singularity. This is a good strategy for ¬ghting blow-up when the model physics is taking a rapid turn for

the nasty.

CHAPTER 11. ALIASING, SPECTRAL BLOCKING, & BLOW-UP

218