<<

. 45
( 136 .)



>>

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

<<

. 45
( 136 .)



>>