t=15.0

-10

-12

-14

-16

0 32

16

Wavenumber k

Figure 11.4: Spectral blocking in a three-dimensional transition-to-turbulence model. As

time increases, more and more resolution is needed; by t = 27.5, the vorticity contour plot

(for this 64 — 64 — 64 run) is a sea of noise. Redrawn from Fig. 4 (“z resolution spectra”) of

Zang, Krist and Hussaini (1989).

A wave equation with constant coef¬cients will have a truncation error ET because the

exact solution is, except for special cases, an in¬nite series whereas the numerical approx-

imation includes only a ¬nite number of terms. However, the Fourier pseudospectral dis-

cretization of a linear, constant coef¬cient wave equation does not display spectral blocking

because waves of different wavenumbers are represented exactly as Fourier coef¬cients. If

the time step is suf¬ciently short, the Fourier algorithm will advect these components with-

out amplitude or phase error so that the high wavenumbers in the truncation will not grow,

but rather will remain the same in amplitude. The lack of amplitude and phase errors in

propagation and advection is one of the great virtues of spectral methods.

If the equation has spatially-varying coef¬cients, either due to nonlinearity or to spa-

tial variations in the wave environment such as a change in the depth of the water, then

even a spectral method will have aliasing error. The product of an x-dependent coef¬cient

with u(x) has a Fourier expansion with components outside the truncation limit which are

aliased to smaller wavenumbers.

Although some spectral blocking is almost inevitable in a long time integration of a

nonlinear system (unless the dissipation is large), there is also a really stupid way to gen-

erate spectral blocking: Use too large a time step. For an explicit time-marching scheme

for the diffusion equation, for example, the high wavenumbers become unstable before

the low wavenumbers as the time step is increased. It is thus possible to accurately and

stably integrate the low wavenumbers while the wavenumbers near the truncation limit

are spuriously amplifying because the time step is too long. Since these high wavenum-

ber components have exponentially small amplitudes at t = 0 (if the initial condition is

smooth), the amplitudes of the high wavenumber components will remain exponentially

small for a ¬nite time interval. However, the instability will cause the graph of the ampli-

tudes of the spectral coef¬cients (on a log-linear plot) to curl up more and more with time,

i. e., exhibit spectral blocking. This result of violating the Courant-Friedrichs-Lewy (CFL)

timestep may be dubbed “CFL Blocking” or alternatively, “Really Stupid Blocking”.

11.3. “2 H-WAVES” AND SPECTRAL BLOCKING 209

Fig. 11.5 is an illustration. There is nothing special about the partial differential equa-

tion, boundary conditions, or initial condition; spectral blocking can occur for any linear

equation when the phase speed or diffusion rate of a mode increases with wavenumber or

degree so that the CFL criterion for aj (t) becomes increasingly restrictive as j increases.

Coeffs. vs. Wavenumber Coeffs. vs. Time

0

10

k=0

k=1

Absolute value of Fourier coefficients

Absolute value of Fourier coefficients

Long

k=2

Timestep

k=3

-2

10 -2

k=4

10

15

k=5

k=

k=6

k=7

k=8

-4

10 -4

10 k=9

4

k=10 k=1

k=11

k=12

k=13

-6 -6

10 10

Short Timestep

0 5 10 15 0 0.5 1

wavenumber k time

Figure 11.5: Spectral blocking in a LINEAR differential equation due to violation of the CFL

timestep limit for LARGE wavenumbers (although the small wavenumbers are STABLE).

The low wavenumbers k = 0, 1, 2, . . . , 11 for basis functions exp(ikx) are stable (as are the

corresponding negative wavenumbers, not shown), but the three highest wavenumbers

grow exponentially with time for the long timestep used (solid curve). The dashed line

illustrates the amplitudes of the exact solution, which are independent of time. The solid

curve was generated by the adaptive fourth/¬ve order Runge-Kutta subroutine ode45 in

Matlab 5.2 with 12 subroutine calls on the time interval t ∈ [0, 1.25]. Because the high

wavenumber Fourier coef¬cients are initially very small, inaccuracies in integrating them

had only a negligible impact on the accuracy of the initial step and so the adaptive rou-

tine happily picked a time step which was unstable for the highest three wavenumbers.

When the calculation was repeated with much more frequent calls to the subroutine ode45,

the adaptive routine was forced to use a short timestep, and all wavenumbers were accu-

rately integrated. The PDE is the linear Benjamin-Davis-Ono equation, ut = ’H(uxx )

where H is the linear operator known as the Hilbert transform; in wavenumber space,

H {exp(ikx)} = isign(k) exp(ikx). Arbitrarily, the spatial period was chosen to be 2π and

the initial condition was u(x, 0) = 2tanh(1)/(1 ’ sech(1) cos(x)).

The CFL criterion can sometimes be exceeded in the middle of an integration even

though the computation is stable at t = 0. The reason is that the advective time-stepping

limit depends on the maximum current, and this can increase during the time evolution of

a ¬‚ow. Fig. 11.6 shows contour plots of four different time levels in the evolution of a dipo-

lar vortex. At t = 10, the ¬‚ow is still perfectly smooth, but then 2h-waves develop in the

zone of intense shear between the two vortices and rapidly amplify until the calculation be-

comes nonsense. The culprit is not lack of spatial resolution or a blocked turbulent cascade,

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

210

timestep 10 timestep 12

2 2

1 1

0 0

-1 -1

-2 -2

-2 0 2 -2 0 2

timestep 14 timestep 16

2 2

1 1

0 0

-1 -1

-2 -2

-2 0 2 -2 0 2

Figure 11.6: Streamlines of the evolution of two contra-rotating vortices. (Negative con-

tours are dashed). At the tenth time step (upper left), the ¬‚ow is very smooth, but 2h-waves

appear at the twelfth step and amplify rapidly. Although spectral blocking is occurring in

the sense that the high wavenumbers are rapidly amplifying after the tenth step, the dif¬-

culty is completely cured by merely halving the time step.

however: the calculation can be extended inde¬nitely merely by halving the timestep.

Elsaesser(1966a) describes a similar example in numerical weather prediction. His Fig.

12 shows a forecast two time steps before over¬‚ow; like our Fig. 11.6 [timestep 14, lower

left panel], Elsaesser™s ¬‚ow is smooth everywhere except for one small region where 2h-

waves have appeared, and his dif¬culty was cured by shortening the time step from four

hours to three hours.

If the timestep is suf¬ciently short, however, spectral blocking can still arise, but only

due to spatially-varying coef¬cients in the differential equation. Nonlinear terms are merely

a special case of spatially-varying coef¬cients.

11.4 Aliasing Instability: History and Remedies

Three-dimensional hydrodynamic codes known as “general circulation models” (GCMs)

have become a major tool in understanding future climate, such as the effects of an increase

in greenhouse gases in the atmosphere. The ur-GCM, the ancestral Adam from whence all

later models came, was Phillips (1956). Starting from rest, and imposing no forcing except a

mean north-south temperature gradient (warm at the equator, cold at the poles), the model

spontaneously developed mean jets and travelling vortices rather closely resembling those

of the real atmosphere. Hurrah!

Unfortunately, after a few more days of model time, the winds became supersonic. The

GCM had blown up!

Phillips tried to stabilize his model by greatly decreasing both the time step and the

11.5. DEALIASING AND THE ORSZAG TWO-THIRDS RULE 211

spatial grid size, but at roughly the same time (about 20 days) the model blew up anyway.

Although his model employed ¬nite differences, he reasoned spectrally. He noticed that

the ¬rst warning of the impending disaster was the appearance of 2h-waves. In Fourier

space, this means that energy is building up near the aliasing limit,

K≡πh ” wavelength = 2h (11.12)

Now the hydrodynamic equations are quadratically nonlinear. The interaction of a wave

with |k| > K/2 with another wave of similar magnitude must generate waves with k > K.

However, these cannot be resolved on a grid with spacing h, but are instead “aliased” to

waves of lower wavenumber. Phillips conjectured that this spurious, high-to-low wavenum-

ber transfer of energy was the source of his instability.

He tested his hypothesis (Philips, 1959) by repeating his earlier experiment with a twist.

Every few steps, he calculated the Fourier transform of the grid point values of each ¬eld

and then set the upper half of the resolved spectrum to zero. That is to say, he applied the

all-or-nothing ¬lter

±

|k| < K/2

ak

ak ’ [f iltered] (11.13)

|k| > K/2

0

It worked! With this ¬ltering, it was possible to integrate his model for arbitrarily long

periods of time.

In a note less than a full page long, Orszag (1971a) pointed out that Phillips™ ¬ltering

was wasteful: if only the upper third of the spectrum were ¬ltered, aliasing would still be

eliminated. (To put it more precisely, aliasing would still occur, but only into wavenumbers

purged by the ¬ltering. Thus, all the surviving wavenumbers are alias-free.)

Akio Arakawa (1966) proposed another solution. Aliasing caused the discretized en-

ergy to grow without bound when it should be conserved. Arakawa therefore developed

a modi¬ed ¬nite difference formula for the Jacobian operator (advection) which exactly

conserved the discretized energy. This also prevented aliasing instability.

Yet another remedy is to use special ¬nite difference or spectral algorithms that have

been modi¬ed to better cope with fronts or shock waves. (The “synoptic scale” vortices in

Phillips™ GCM should form narrow fronts as they mature; it was during this “frontogene-

sis” stage that Phillips™ model went bad.)

All three strategies ” ¬ltering, energy-conservation, and shock-resolution ” have their

uses. However, if a ¬‚ow is being well-resolved by an exponentially-convergent spectral

model, shouldn™t the energy be conserved to a high approximation anyway? Are these

alias-¬xers really necessary? Moreover, it turns out that these remedies offer seductive

possibilities for self-deception. We will return to these dif¬cult questions after ¬rst describ-

ing a couple of these procedures in more detail in the next two sections.

11.5 Dealiasing and the Orszag Two-Thirds Rule

Phillips suggested a dealiasing procedure: apply a spatial ¬lter so as to suppress all waves

with wavelengths between 2h and 4h. For spectral methods, we ¬lter by simply deleting

all the offending wavenumbers just before each transform to grid point values.

Orszag (1971a) showed that purging half the spectrum is wasteful. If we ¬lter all waves

such that |k| > (2/3)K, then the quadratic interaction of two wavenumbers p, q that survive

the ¬lter will alias only to wavenumbers that are purged. The reason is that the aliasing

shift in k must always be by a multiple of 2K (Fig. 11.7). The conclusion is that it is only

necessary to ¬lter waves with wavelengths between 2 h and 3 h to suppress aliasing.

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

212

k ∈ [-(4/3) K,-K]

(2/3) K k ∈ [K,(4/3) K]

k

K

FILTERED UNFILTERED 0

-K

- (2/3) K

Figure 11.7: The “Aliasing Wheel”. The wavenumber k for Fourier basis functions

exp(i k x) is depicted as an angle in a polar coordinate system, scaled so that the basis

set k ∈ [’K, K] spans the full circle. Orszag™s “Two-Thirds Rule” ¬lter has no effect on

waves in the range k ∈ [’(2/3)K, (2/3)K], but all Fourier coef¬cients with |k| > (2/3)K

are set equal to zero at the end of each time step.

Quadratically nonlinear terms have Fourier components that span the range k ∈

[’(4/3)K, (4/3)K], which is represented by the double-headed arrow that completes one-

and-one-third complete turns of the circle. As shown by the diagram, wavenumbers

on k ∈ [K, (4/3)K] are aliased to k ∈ [’K, ’(2/3)K] (lower shaded band). Similarly,

wavenumbers on k ∈ [’(4/3)K, ’K] are aliased to k ∈ [(2/3)K, K] (upper shaded area).

The ¬ltering removes all effects of aliasing because the misrepresented wavenumbers ”

the shaded bands ” lie within the “death zone” of the ¬lter. The dotted horizontal ray is

the branch cut: all k > K or < ’K will be aliased.

Rule-of-Thumb 10 (TWO-THIRDS RULE) To obtain an alias-free computation on a grid of N

points for a quadratically nonlinear equation, ¬lter the high wavenumbers so as to retain only

(2/3)N un¬ltered wavenumbers.

This precept is alternatively (and confusingly) known as the “Three-Halves Rule” because to

obtain N un¬ltered wavenumbers, one must compute nonlinear products in physical space on a