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

8.4. OTHER DISCRETE SYMMETRIES 167

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,

1952.)

¬‚ow that develops from the initial condition:

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

= ’ cos(x) sin(y) cos(z) (8.15)

v

2

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

3

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

∞∞∞

(8.16)

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

CHAPTER 8. SYMMETRY & PARITY

168

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

computer.

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.

8.4. OTHER DISCRETE SYMMETRIES 169

V=1

1

2

/

=1

V=0

V=0

V

0

‚V/‚ z=0

-1

0 1

-1 V=1

x

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.

1

0.9

0.8

0.7

0.5 0.6

0.5

0.5

0.1

0 0.1 0.4 0.3 0.2

0.2 0.3 0.4

0.5 0.5

0.6

-0.5

0.7

0.8

0.9

-1

-1 -0.5 0 0.5 1

x

Figure 8.6: Solution to Richardson™s Laplace problem.

CHAPTER 8. SYMMETRY & PARITY

170

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

1

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)

u=

j=0

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.

8.5. AXISYMMETRIC & APPLE-SLICING MODELS 171

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

knowledge.”

” 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