cos(mθ)

φm (r, θ) ≡ Wn (r)

m

(18.15)

n sin(mθ)

where in Verkley™s notation

0,m

Wn (r) ≡ rm P(n’m)/2 (2r2 ’ 1),

m

m = 0, 1, 2, . . . ; n = m, m + 2, m + 4, . . .

(18.16)

0,m

where Pk (s) is the Jacobi polynomial of degree k of order (0, m) in its argument s. (Note

that because s ≡ 2r2 ’ 1, the Jacobi part of the basis function is of degree (n ’ m) in r.)

This basis function explicitly enforces all the constraints of the parity theorem; each basis

function is symmetric in r for even m and odd in r for odd m. The factor of rm ensures that

each basis function has an m-th order zero at the origin. Verkley(1997a, Table 1) provides

the explicit form of these basis functions for small n and m.

The Jacobi polynomials are orthogonal so that

1

±,β ±,β

Pk (s)Pk (s)(1 ’ s)± (1 + s)β ds = 0, (18.17)

k=k

’1

Changing variables from s to r and setting ± = 0, β = m gives

1

0,m 0,m

Pk (2r2 ’ 1)Pk (2r2 ’ 1)r2m rdr = 0, (18.18)

k=k

0

or in other words,

1

m m

(18.19)

Wn (r)Wn (r) rdr = 0, n=n

0

These Jacobi polynomials of order (0, m) are “one-sided” in the sense that the weight

function in Eq.(18.18) diminishes rapidly (as r2m ) as the origin is approached, but varies

only slightly near the outer edge of the disk. If the Jacobi polynomials oscillated uniformly

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

388

in r like the Robert functions, they would be nearly-dependent and ill-conditioned like the

Robert functions. The orthogonality constraint, however, forces the polynomials of differ-

ent degree to be as independent as possible. This requires that the polynomials oscillate

mostly near r = 1: the roots of the basis functions move closer and closer to the outer

boundary for ¬xed degree as the order m increases.

The one-sideness is the magic. The basis functions of high wavenumber m are the

villians that describe the fast tangential advection or diffusion that limits the timestep.

However, the radial parts of these basis functions, if one-sided Jacobi polynomials, have

only negligible amplitude close to the origin where the grid points are close. The rapid

advection or diffusion is suppressed, and a much longer timestep is possible.

As with spherical harmonics, discussed later, there are two options for truncating the

basis: “rectangular” and “triangular”. A rectangular truncation employs the same num-

ber of radial basis functions for each angular wavenumber m. It is the easiest to program,

but has fallen from grace in spherical coordinates. The triangular truncation decreases the

number of radial basis functions by m with each increase in m until the highest wavenum-

ber has but a single radial basis function. It has been shown for spherical harmonics that

this gives the most uniform resolution (and allows the longest timestep) of any truncation

that includes a given maximum zonal wavenumber M . This property of “equiareal resolu-

tion”, de¬ned formally below in Sec. 18.13, has not been proved for the One-Sided Jacobi

basis, but it seems plausible that a triangular truncation is preferable here, too.

The usual strategy of choosing the roots of the polynomials as the interpolation points

(0,m)

are different for each polar wavenumber

does not work here because the zeros of the Pk

m. To take transforms in θ, we need the unknowns de¬ned on a single set of radial grid

points. Therefore, one is forced to abandon a pseudospectral strategy and instead apply

a Galerkin method. The integration points are the Legendre-Radau grid which includes

r = 1 but not r = 0. To avoid quadrature errors, it is necessary to use more grid points than

unknowns, but the transform from radial grid point values to coef¬cients can be still be

expressed as the multiplication of a vector of grid point values or coef¬cients by a matrix,

albeit now a rectangular matrix.

These basis functions, unlike the spherical harmonics they mimic, are not eigenfunc-

tions of Laplace operator. However, both Matsushima and Marcus(1995) and Verkley(1997a)

derive recurrence relations which allow the Laplacian to be inverted by solving a penta-

diagonal matrix problem. Similarly, derivatives can be evaluated at an operation count

directly proportional to the number of spectral coef¬cients in the truncation. Overall, the

ef¬ciency of this basis is roughly the same as for spherical harmonics. Indeed, the algorith-

mic connection with spherical harmonics is so close that Verkley™s shallow-water solver

for the disk is merely a rewrite of a previous code that solved the same hydrodynamic

equations on the surface of a sphere via spherical harmonics.

Unfortunately, the Jacobi basis shares the vices of spherical harmonics. First, there is no

Fast Fourier Transform for the Jacobi polynomials, so all radial transforms from grid point

values to spectral coef¬cients must be done by Matrix Multiplication Transform (MMT).

This costs O(N 2 ) for each wavenumber m versus only O(N log2 (N )) for an FFT. (Trans-

forms in θ are still done by FFT, thank goodness.)

Second, the messy recurrence relations and basis functions require some investment of

both learning time and programming time. (We may hope that someday the necessary

subroutines will become widely available in software libraries, and then this objection will

disappear).

There is one profound difference from spherical harmonics, however. The One-sided

Jacobi basis has the usual high density of grid points near r = 1 where grid points are sepa-

rated in r by O(1/N 2 ). With a triangular truncation so that the maximum polar wavenum-

ber M is equal to N , the tangential distance between grid points on the outermost ring is

18.5. RADIAL BASIS SETS AND RADIAL GRIDS 389

Chebyshev Expansions of J0

0

10

-2

10

-4

10

T (2 r - 1)

j

-6

10

2

-8

Tj(2 r - 1),

10

T (r)

-10 2j

10

-12

10

0 2 4 6 8 10

degree j

Figure 18.2: Chebyshev expansions of the Bessel function J0 (r). The slowly-convergent

series (dashed) is the expansion in Shifted-Chebyshev polynomials: Tj— (r) ≡ Tj (2r ’ 1).

The two good techniques are (i) an expansion in Shifted-Chebyshev polynomials with a

QUADRATIC argument Tj (r 2 ) = Tj (2r 2 ’ 1) and (ii) an expansion in the EVEN Chebyshev

polynomials without a shift: T2j (r). These two options generate only a single curve be-

cause (i) and (ii) are in fact IDENTICAL because of the identity Tj (2r2 ’ 1) = T2j (r).

also O(1/N 2 ). This implies that the One-sided Jacobi basis has not completely solved the

dif¬culty of a restrictive timestep; indeed, if the time scales for tangential advection or dif-

fusion near the origin are similar to the time scales for radial advection or diffusion near

the boundary, then the switch from a Chebyshev basis to the Jacobi basis will not increase

the time step at all.

Nevertheless, this basis has worked well in Matsushima and Marcus(1995) and Verkley

(1997a,b).

18.5.2 Boundary Value & Eigenvalue Problems on a Disk

When there is no timestep, the pressure to use One-Sided Jacobi basis functions in radius is

greatly reduced. Orszag(1974) and Boyd(1978c) have shown that satisfying all the pole con-

ditions is not terribly important to numerical ef¬ciency; a regular Chebyshev basis seems

to work just ¬ne. (Their papers actually solve problems on the sphere, but because the

polar caps are well-approximated by small disks, their conclusions should apply to polar

coordinates, too.) Matsushima and Marcus(1995) show by solving the eigenproblem whose

exact solution is the Jm Bessel function that the Jacobi basis can reach a given error toler-

ance with only half as many coef¬cients as the Parity-Restricted Chebyshev polynomials

when m = 50. However, in a Fourier-Chebyshev series in two dimensions, the Fourier

series in m usually converges rapidly, so the Chebyshev basis is converging slowly only

for wavenumbers whose amplitude is negligible anyway.

Our recommendation for boundary and eigenvalue problems is to try a Chebyshev

basis ¬rst. If the speed of the resulting (rather easy-to-program) code is unsatisfactory,

the One-Sided Jacobi basis can be used as a fallback with the Chebyshev program used to

check the newer and much more intricate Jacobi code.

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

390

Table 18.2: Cylindrical or Polar Coordinates: Unbounded Domain or External to a Cylinder

References Comments

Zebib(1987a) Stability of ¬‚ow past a cylinder

Don&Gottlieb(1990) unsteady ¬‚ow past a cylinder

Deane & Kevrekidis spectral elements; low order basis of

&Karniadakis & Orszag(1991) empirical orthogonal eigenfunctions

Mayer&Powell(1992) Eigenvalue calculation of stability of trailing

vortex through domain truncation

Mittal & Balachandar(1996) ¬‚ow past elliptic cylinder in elliptic coordinates

good discussion of boundary conditions,

blending in¬‚ow & out¬‚ow

Blackburn&Henderson(1996) ¬‚ow past a vibrating cyclinder

Matsushima&Marcus(1997) Unbounded domain including the origin;

2D and 3D vortex ¬‚ows

18.5.3 Unbounded Domains Including the Origin in Cylindrical Coor-

dinates

Matsushima&Marcus(1997) have extended their earlier work by applying an algebraic

change-of-coordinate in r to map an unbounded domain into a disk. They then apply the

One-Sided Jacobi basis described above. Like the rational basis T Bj for the in¬nite interval

described in Chapter 17, the Matsushima-Marcus basis functions are rational functions of

r ∈ [0, ∞].

The alternative is the even T B functions for even polar wavenumber m and the T B2j’1 (r)

for odd m: a Parity-Restricted Rational Chebyshev basis. The simpler functions are prefer-

able for boundary and eigenvalue problems, but the Matsushima-Marcus basis may give

faster convergence and allow a longer timestep.

Flows exterior to a cylinder in an unbounded domain have been the subject of much

study, both analytical and numerical. Exterior ¬‚ows are free of the pole problem, but im-

posing boundary conditions on a ¬‚ow at large radius can be tricky. We have therefore

collected some illustrations in Table 18.2.

18.6 Annular Domains

Flow in an annulus that does not include the origin is free from the pole problem. We

recommend a Fourier basis in the angle θ and a standard Chebyshev basis in the radial

coordinate r. In three dimensions, one should use a Chebyshev basis in the axial coordinate

z unless periodic boundary conditions in z are imposed, in which case a Fourier basis is

greatly preferable.

There is little difference between annular ¬‚ow and channel ¬‚ow except that cylindri-

cal coordinates introduce an r-dependent “metric” factor into the Laplace operator, etc.,

whereas the Laplace operator has constant coef¬cients in Cartesian coordinates. The usual

methods for separable partial differential equations, discussed in Chapter 15, Sec. 11, still

apply.

For most practical purposes, however, an annular channel is simply a straight channel

which is periodic in the downstream direction. A standard Fourier/Chebyshev basis is

best. Some representative works are listed in Table 18.3.

18.7. SPHERICAL COORDINATES: AN OVERVIEW 391

Table 18.3: Annular Flows

References Comments

Marcus(1984a,1990) Taylor-Couette ¬‚ow

Le Qu´ r´ &Pecheux(1989,1990)

ee axisymmetric convection; bifurcations

Chaouche(1990a) axisymmetric ¬‚ows;

Chaouche et al.(1990) in¬‚uence matrix method

Randriamampianina(1994) 3D ¬‚ow; vorticity-vector-potential

18.7 Spherical Coordinates: An Overview

The pattern of latitude and longitude lines on the surface of the sphere is, near the poles,

the same as that of a polar coordinate system in a disk. So, it is hardly surprising that there

are many similarities between cylindrical coordinates and spherical coordinates. However,

there are important differences.

One is that a sphere has two poles rather than one. A second, more important difference

is that the surface of a sphere has no boundaries. This has the profound consequence that

spherical basis sets automatically and individually satisfy the behavioral boundary condi-

tions on the sphere. In contrast, it is necessary to impose numerical boundary conditions at

the boundary circle of the disk.

A third difference is that for the sphere, there is an “obvious” basis set, the spherical har-

monics, which are nearly ideal. Spherical harmonics give equiareal resolution, exponential

convergence, and trivial inversion of the Laplace operator, which is the eigenoperator for

these functions. In contrast, the obvious basis for the disk, the cylindrical harmonics, are

Bessel functions of radius with only a poor, algebraic rate of convergence. One can match

most of the good properties of spherical harmonics by using the One-Sided Jacobi basis for

the disk, but these are not eigenfunctions of the Laplace operator.

The fourth difference is that the number of simulations in a cylinder is fairly modest

compared to the large body of work that has been done in spherical coordinates for weather

forecasting, climate modelling, mantle convection and stellar ¬‚ows. For this reason, we

shall discuss spectral methods on the sphere in great detail.

18.8 The Parity Factor for Scalars: Sphere versus Torus

Topologically, the sphere is a two-dimensional manifold of genus zero while the torus

(“doughnut”) is of genus one. Ironically, however, it is easier to compute on the surface

of a torus than the surface of a sphere.

Science ¬ction writers have described toroidal planets (Boyd, 1981c, 1984c) and plasma

physicists solve ¬‚ows in tori to model the ring-shaped fusion generators known as toka-

maks (Schnack et al., 1984). In this section, we will discuss a two-dimensional basis for the

surface of a torus. However, one may de¬ne a three-dimensional orthogonal coordinate

system (“toroidal coordinates”, Morse and Feshbach, 1953) in which the third coordinate

is radial distance from the surface of a torus to its centerline.