. 46
( 136 .)


grid of (3/2)N points.
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

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
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.

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-
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
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.

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

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.

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-
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!
h=100 m


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.


. 46
( 136 .)