. 37
( 136 .)


Figure 8.3: Grids for Fourier series.

Although this invariance is not an automatic property of basis functions, it is straight-
forward to take linear combinations of Chebyshev polynomials (or whatever) which are
invariant under C4 .
A real world example of such a function is the third derivative of the similarity solu-
tion for a laminar boundary layer at separation (Boyd, 1997a). Functions which have C3
symmetry in the complex x-plane or are the products of such functions with a power of x
include the Airy functions Ai(x) and Bi(x) and the general, unseparated similarity solution
for boundary layer ¬‚ow.
Rotational symmetry is sometimes found in connection with parity and/or re¬‚ection
symmetry. Fig. 8.4 shows a familiar example: a snow¬‚ake, which is not only invariant
under the “rotation” or “cyclic” group C6 but also possesses the property that each of its
six points is symmetric about its midpoint. In the language of group theory, snow¬‚akes
belong to the “dihedral” group D6 . The signi¬cance of “dihedral” as opposed to “cyclic”
symmetry is the we can reduce the basis set still further as shown below:

basis set is {1, cos(nx), sin(nx), cos(2nx), sin(2nx), . . . }
Cn :
basis set is {1, cos(nx), cos(2nx), . . . }
Dn :

where x is an angular coordinate. Unfortunately, dihedral symmetry is rare: Travelling
waves will disrupt this symmetry unless we shift to a coordinate system that is travelling
with the waves.
Two- and three-dimensional symmetries also exist, and some are not simple combina-
tions of one-dimensional symmetries. A famous example is the Taylor-Green vortex, which
has been widely studied as a model for a three-dimensional turbulent cascade. This is the

Figure 8.4: Illustration of the discrete rotation group Cn [no mid-sector symmetry] and the
dihedral group Dn [symmetric about some axis in each sector] for various n. (After Weyl,

¬‚ow that develops from the initial condition:

u = sin(x) cos(y) cos(z)

= ’ cos(x) sin(y) cos(z) (8.15)
w= cos(x) cos(y) sin(z)
with boundary conditions of periodicity of 2π in all three coordinates and constant viscos-
ity. One can show that the symmetries inherent in the initial condition ” u is antisymmet-
ric about x = 0 and symmetric about y = 0 and z = 0 ” are preserved for all time, even
though the equations of motion are nonlinear. Consequently, the solution for all times can
be expanded as
∞ ∞ ∞
u = umnp (t) sin(mx) cos(ny) cos(pz)
m=0 n=0 p=0
v = vmnp (t) cos(mx) sin(ny) cos(pz)
m=0 n=0 p=0
w= wmnp (t) cos(mx) cos(ny) sin(pz)
m=0 n=0 p=0

The symmetries inherent (8.16) allow a reduction of a factor of eight in the total number
of basis functions since we need only half of the possible one-dimensional basis functions

in each of x, y, and z as factors for our three-dimensional functions. This alone is a huge
savings, but Appendix I of Brachet et al. (1983) shows that the ¬‚ow is also symmetric with
respect to rotations through 180 degrees about any of the three axes: (i) x = y = π/2 (ii)
x = z = π/2 and (iii) y = z = π/2. In mathematical terms, this implies that 7 out of every 8
coef¬cients in (8.16) are zero: m, n, and p must differ by an even integer (including 0). The
¬‚ow is also invariant under a rotation by 90 degrees about the vertical axis x = y = π/2,
but the computer code did not exploit this ¬nal symmetry.
The total number of degrees of freedom is reduced by a factor of 64. The result is a three-
dimensional code which, at 643 unknowns per ¬eld, is equivalent in accuracy to a general
spectral code with 256 grid points in each of x, y, and z. Even 643 taxed the memory of
the best supercomputer of 1983, the Cray-1, so that the coef¬cients and the corresponding
grid point representations would not both ¬t in core memory at the same time, forcing an
intricate data shuttle to disk storage.
Brachet et al. (1983) note: “The present work has been possible because of the develop-
ment of new algorithms that take advantage of the symmetries of the Taylor-Green ¬‚ow”.
The “new algorithm”, given in their Appendix III, is an ef¬cient FFT for real Fourier series
whose non-zero coef¬cients are all odd degree such as {1, cos(x), cos(3x), cos(5x), . . . }.1
Boyd (1986c) exploits a two-dimensional symmetry that has no counterpart in one di-
mension. In addition to parity in both x and y, the solution of the nonlinear PDE known
as “Bratu™s equation” is invariant under the C4 rotation group, which is equivalent to the
replacement of x by y and y by x. His article shows how to exploit the combined dou-
ble parity/C4 symmetries to reduce the basis set by a factor of eight, which allowed him to
solve the nonlinear two-dimensional BVP to ¬ve decimal places on a 1985-vintage personal
Similarly, L. F. Richardson, who made great contributions to numerical analysis, the
theory of war and con¬‚ict, meteorology, and fractals, solved a two-dimensional partial
differential in 1908 without the aid of a digital computer. He attacked Laplace™s equation
in the unit square subject to the boundary condition V (x, y = ±1) = 1 on the top and
bottom boundaries while V (x = ±1, y) = 0 on the sides. As shown in Fig. 8.5, the square
can be divided into eight triangles by horizontal, vertical, and diagonal lines. The vertical
and horizontal lines are lines of symmetry: the solution is symmetric with respect to both
x = 0 and y = 0. The diagonals are lines of antisymmetry, but there is a subtlety: it is
V (x, y) ’ 1/2, not V (x, y) itself, which is antisymmetric with respect to the diagonals:

V (x + δ, x) ’ 1/2 = ’ {V (x, x + δ) ’ 1/2} (8.17)

where δ is arbitrary and similarly for the other diagonal (Fig. 8.6).
Marcus(1984a,b, 1990) employs a “re¬‚ect-and-shift” symmetry for the velocities that
halves the number of Fourier modes (and the cost!).
Weyl(1952) and Hargittai and Hargittai(1994) are good popular treatments of symmetry
in nature and art.

1 Kida (1985) set a new “record” for symmetry by identifying a ¬‚ow in which storage requirements are only
1/192 of what they would be for a general Fourier series with the same resolution.




‚V/‚ z=0

0 1
-1 V=1

Figure 8.5: Richardson™s Laplace equation problem. By exploiting symmetry; it is suf¬cient
to solve the problem only in the shaded triangle with the indicated boundary conditions.
Symmetry then gives the solution in the other seven triangles.

0.5 0.6


0 0.1 0.4 0.3 0.2
0.2 0.3 0.4

0.5 0.5
-1 -0.5 0 0.5 1

Figure 8.6: Solution to Richardson™s Laplace problem.

8.5 Axisymmetric, Hemispheric & Apple-Slicing Models
If a physical problem does not have exact symmetries like the Taylor-Green problem, it
may still be possible to obtain big savings by making the approximation that various sym-
metries exist. Although unacceptable in operational forecasting, approximate symmetries
are widely used in meteorological research.
The crudest approximation is axisymmetry: The ¬‚ow is independent of longitude. In
spectral terms, the three-dimensional atmospheric ¬‚ow can always be represented without
error as a Fourier series in longitude with Fourier coef¬cients that are functions of latitude,
height, and time. The assumption of axisymmetry is equivalent to truncating the Fourier
series at N = 0, i. e. keeping only the constant term. The longitude-independent approxi-
mation is terrible for day-to-day forecasting, but gives useful insights into climate.
Another simpli¬cation is the hemispheric model: A grid that spans only the northern
hemisphere is justi¬ed by the ¬ction that the ¬‚ow has equatorial parity. “Symmetry about
the equator” has a special meaning to a meteorologist: From the Navier-Stokes equations,
one can show that it is not possible for all three ¬‚ow variables to be simultaneously sym-
metric about the equator. For example, the (linearized) x-momentum equation is

ut ’ f v = ’ px (8.18)

Since the Coriolis parameter is antisymmetric about the equator, it follows that v must have
parity opposite to that of u and p. Hemispheric models using spherical harmonics therefore
use symmetric basis functions for u, vertical velocity w, and p, but employ antisymmetric
harmonics to expand the north-south velocity v and the potential vorticity q.
A third approach is “apple-slicing”: Assuming that the spectrum is dominated by one
longitudinal wavenumber, and keeping only that wavenumber and its harmonics in the
basis set. The theory of “baroclinic instability” shows that when the global weather pat-
terns are expanded in a Fourier series, the coef¬cients do not monotonically decrease with
zonal wavenumber k but rather show a large peak in the range of k ≈ 6 ’ 8. A crude-but-
quick spectral approach is to keep just zonal wavenumbers from the set {0, 8, 16, 24, . . .
The reason for the name “apple-slicing” is that this sort of truncation is equivalent to
assuming that the earth can be cut into eight wedges, bounded by meridians, with identical
¬‚ow on each wedge. Consequently, solving for the weather on a single wedge gives the
¬‚ow on all the other wedges. The “apple-slicing” approximation is equivalent to assuming
that the weather is invariant under rotations through 360/m degrees for some integer m.
The spectral series in longitude is then

umj (φ, z, t) eimj» (8.19)

Fig. 8.7 illustrates this idea.
Simmons and Hoskins (1975) show that their “apple-sliced” model with m = 8, which
they used for preliminary experiments, needed only 0.084 seconds on the CDC7600 per
time step versus the 0.67 seconds/time step of their full model, which kept all spherical
harmonics up to a zonal wavenumber of 21.
Apple-slicing is a very good way of testing new paramerization schemes, vertical reso-
lution and so on. Half a dozen “apple-sliced” experiments have the same cost as a single
full-basis simulation. The disadvantages are obvious: it only applies to special problems
where the spectrum is dominated by intermediate wavenumber, and it is usually a rather
crude approximation.

Figure 8.7: Schematic of “apple-slicing” on the sphere.
(a) Side view; the dynamics in the shaded sector [which is the computational domain] is
assumed to repeat periodically in all the other sectors of the sphere.
(b) A polar projection with the north or south pole at the center of the diagram, illustrating
C5 symmetry in a contour plot of a variable such as pressure or zonal velocity.
Chapter 9

Explicit Time-Integration Methods

“ . . . In my experience, many people who do computing are reluctant to look at numbers.
At Stanford, the general level of our students has been pretty high, but . . . their main
weakness is in their inability to look at outputs and extract the meaningful information . . .
In fact, they are generally less ef¬cient than the assistants I used to have . . . in the ACE days
[1948], in spite of having far superior mathematical quali¬cations. Most of those assistants
had experience with desk computers and had learned to “look at numbers”.”
“I certainly do not want to suggest that the way to acquire this habit is to serve an
apprenticeship on desk computers, but we have yet to learn how to instill the relevant
” James Wilkinson, interviewed in BYTE, Feb. ™85

9.1 Introduction
An “implicit” time-marching scheme is one which requires solving a boundary value prob-
lem at every time step. “Explicit” schemes give the solution at the next time level in terms
of an explicit formula: un+1 = stuff where “stuff” is terms that can be evaluated by using


. 37
( 136 .)