Chebyshev–Fourier spectral methods in bipolar coordinates

Journal of Computational Physics 295 (2015) 46–64

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Chebyshev–Fourier spectral methods in bipolar coordinates Zhu Huang a , John P. Boyd b,∗ a b

School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an, PR China Department of Atmospheric, Oceanic & Space Science, University of Michigan, 2455 Hayward Avenue, Ann Arbor, MI 48109, United States

a r t i c l e

i n f o

Article history: Received 17 November 2014 Received in revised form 4 March 2015 Accepted 22 March 2015 Available online 7 April 2015 Keywords: Pseudospectral Chebyshev polynomials

a b s t r a c t Bipolar coordinates provide an efficient cartography for a variety of geometries: the exterior of two disks or cylinders, a half-plane containing a disk, an eccentric annulus with a small disk offset from the center of an outer boundary that is a large circle, and so on. A pseudospectral method that employs a tensor product basis of Fourier functions in the cyclic coordinate η and Chebyshev polynomials in the quasi-radial coordinate ξ gives easy-to-program spectral accuracy. We show, however, that as the inner disk becomes more and more offset from the center of the outer boundary circle, the grid is increasingly non-uniform, and the rate of exponential convergence increasingly slow. One-dimensional coordinate mappings significantly reduce the non-uniformity. In spite of this non-uniformity, the Chebyshev–Fourier method is quite effective in an idealized model of the wind-driven ocean circulation, resolving both internal and boundary layers. Bipolar coordinates are also a good starting point for solving problems in a domain which is not one of the “bipolar-compatible” domains listed above, but is a sufficiently small perturbation of such. This is illustrated by applying boundary collocation with bipolar harmonics to solve Laplace’s equation in a perturbed eccentric annulus in which the diskshaped island has been replaced by an island bounded by an ellipse. Similarly a perturbed bipolar domain can be mapped to an eccentric annulus by a smooth change of coordinates. © 2015 Elsevier Inc. All rights reserved.

1. Introduction Classic Chebyshev and Fourier spectral methods are most effective in multivariate applications when the geometry is a direct product of one-dimensional domains. This allows the construction of multivariate basis sets and grids as tensor products of one-dimensional basis sets and grids. In the eighteenth, nineteenth and twentieth centuries, mathematical physicists and applied mathematicians constructed many such direct product domains and the corresponding coordinate systems so as to apply the method of separation of variables. Before the computer age, this semi-analytical method reigned supreme for solving partial differential equations although there was some application of finite difference and spectral methods by hand [55,22]. Such direct product or tensor product domains are very amenable to spectral methods even when the partial differential equation at hand is not actually separable in the domain. In this article, we explore the numerical ramifications of applying Chebyshev–Fourier spectral methods in bipolar coordinates. The motive was an application to idealized ocean flow [34,12] as illustrated in Section 9. Numerous numerical issues

*

Corresponding author. E-mail address: [email protected] (J.P. Boyd).

http://dx.doi.org/10.1016/j.jcp.2015.03.062 0021-9991/© 2015 Elsevier Inc. All rights reserved.

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

47

Table 1 Notation. B (x, y )

Function whose zero isolines define the boundary

D N r R s

Total degree of a polynomial (for xm+n , D = m + n) Total number of basis functions Radial coordinate in polar coordinates Radius of the outer boundary circle of an eccentric annulus Relative displacement of the center of the smaller boundary circle from the center of the larger, outer boundary of an eccentric annulus Non-uniformity ratio: the maximum, on a circle of fixed ξ , of the Jacobian function divided by the minimum Relative radius of inner boundary circle in an eccentric annulus; the absolute radius is R  Alternative quasi-radial coordinate Quasi-radial bipolar coordinate Inner boundary is the circle ξ = ξi Outer boundary is the circle ξ = ξo Angular coordinate in polar coordinates Angular bipolar coordinate

U(ξ )

 ρ ξ ξi ξo θ

η

Fig. 1. Contours of constant quasi-radial coordinate ξ (left panel) and quasi-angular coordinate η (right panel) in the x–y plane. The dashed circle is the boundary between the interior and exterior regions; the bipolar coordinate system is defined for all (x, y ), but in this application, we employ it only for the flow exterior to the boundary at ξ = ξ0 where ξ0 is arbitrary, but chosen to be 1 for this schematic [dashed circle]. The contours of both coordinates are circles. The boundary disk is always a circle of constant ξ ≡ ξ0 . The “circle” ξ = 0 includes the antisymmetry line, x = 0, as well as infinity. ξ is everywhere negative in the left half-plane as illustrated not here but on p. 1210 of Morse and Feshbach [52].

arose in the course of that work which suggested that an article devoted primarily to computational efficiency and theory would be useful. In the process, we also shed light on the method of separation of variables in general. All notations used in this article are given in Table 1. Let the poles which give the coordinate systems its name be located at (±d, 0) in Cartesian coordinates (x, y ). Then bipolar coordinates (ξ, η) are



ξ = arctanh



2dx d2 + x2 + y 2



η = arctan

,

2dy d2 − x2 − y 2



(1)

The inverse transformation is

x=d

sinh(ξ ) cosh(ξ ) + cos(η)

,

y =d

sin(η) cosh(ξ ) + cos(η)

(2)

Other essential formulas, such as the metric factors required to transform partial differential equations from Cartesian coordinates to bipolar coordinates will be given in later sections. Note that there is no unanimity of conventions. Other authors replace η in our definition by −η or η + π , rotate the coordinates by ninety degrees, etc. Our convention is that of Morse and Feshbach [52]. Lucht [49] gives an extensive comparison of different conventions. The bipolar coordinate contours are circles as shown in Fig. 1. Each circle of constant ξ = ξ0 has a center at

xC = d cosh(ξ0 )/ sinh(ξ0 )

(3)

and radius

R = d/ sinh(ξ0 )

(4)

48

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

Each circle such that ξ(x, y ) =  for some constant  > 0 is enclosed by all circular isolines of ξ such ξ < . (Note that somewhat counterintuitively, the limit ξ → ∞ is a contour of tiny radius around the right focus, (x = d, y = 0) while ξ = 0 is the y-axis.) In mathematics/geology jargon, the ξ isolines form a layering or “foliation” of the x– y plane. The contours of constant η are also circles, but every such circle intersects every other circle at the two points they have in common, the poles (±d, 0). Bipolar coordinates have proved useful for a variety of different geometries including: 1. 2. 3. 4. 5.

Eccentric annulus. Exterior of two disks or cylinders. Single disc or cylinder and a half-plane or half-space. Infinite plane with semicircular trough. Any domain that can be smoothly mapped into one of the domains described above.

Definition 1 (Eccentric annulus). An eccentric annulus is a domain bounded by two circles, one of radius R and the other of radius R  where  < 1 and where the center of the smaller circle is displaced from the center of the larger circle by a distance Rs where by convention s > 0 regardless of the direction of the translation. If, without loss of generality, we choose the translation to be leftward along the x-axis where this axis is chosen to pass through the centers of both circles, then the eccentric annulus can be equivalently defined in bipolar coordinates as follows with specification of R,  and s:



r1 ≡

1 −  2 − s2 +

ξi = log(r1 ) log(r1 ) ξo =  + sr1

 √ ( + 1 + s) (1 +  − s) (1 −  + s) (1 −  − s)

[Coordinate foci at] x = Rd,

2 s (5) (6) d = sinh(ξo )

where ξi and ξo are the bipolar quasi-radial coordinates of the inner and outer boundaries, respectively. The center of the big circle is at

xc = R cosh(ξo )

(7)

and the center of the smaller circle is at R (xc − s). Flow between eccentric rotating cylinders is very useful in the theory of bearings and lubrication and has been extensively investigated [18,20,21,23,29,30,32,42–44,46]. Other applications of bipolar coordinates in the eccentric annulus include [34,51,42–44,46,66]. Chen, Tsai and Liu [17] compared conformal mapping and bipolar coordinates for solving the Laplace equation in this geometry. Convective flow and heat transfer in an eccentric circular annulus have also been studied by [33,27,59,15,64,50,19]. The axial-coordinate-independent flow between two circular cylinders which are not coaxial but have parallel axes has been the topic of many articles. Fixed cylinders are treated in [25,45]. Rotating circular cylinders have been examined in [41,61,53,57,58,60,65]. Chemical engineering plants often have mixers with two parallel rotating blades or cylinders. A bicylindrical geometry is popular for ore-crushers. Solutions to PDEs in the half-plane, x ∈ [0, ∞] ⊗ y ∈ [−∞, ∞], are discussed in many textbooks; Cartesian coordinates are usually satisfactory for problems in a half-plane. When it is necessary to impose boundary conditions on a circle as well as on the boundary of the plane, x = 0, bipolar coordinates are useful as illustrated in [52]. Oceanic flow between a long, straight stretch of coast and an offshore island is an obvious application. A modon is a pair of contra-rotating vortices, translating at a constant speed under their mutually interacting tangential velocities as a long-lived coherent structure. They are also called “Lamb–Chaplygin dipoles” after the two fluid dynamicists who discovered these exact analytic two-dimensional solutions around the turn of the twentieth century. These vortex pairs form spontaneously in turbulence and in river estuaries. Modons are of great interest in atmospheric and oceanic dynamics as well as other areas of fluid mechanics. In the nineteenth century, Helmholtz had previously discovered that a pair of equal but opposite “point vortices”, which are idealizations in which all of the vorticity is concentrated in a single point, would propagate at a steady velocity along a line perpendicular to a line segment drawn from one point vortex to its mate. (The circulation of a point vortex is finite, but the vorticity has the spatial structure of a Dirac delta function.) The Chaplygin–Lamb solutions have all the positive vorticity and all the negative vorticity distributed smoothly over two half-disks; outside the union of the two semidisks, the flow is irrotational. Although no exact solutions are known except for these two extreme cases, Eyedland and Turkington showed numerically that a complete spectrum of intermediates exist in which the vorticity of a given sign is concentrated in a roundish patch whose boundary is approximately a circle [26]. It seems likely by analogy with the ellipse-bounded modons computed by Ma and Boyd [11] that modons exist in which the flow is irrotational everywhere outside the bipolar coordinate contours, ξ = ±ξb , for all ξb ∈ [0, ∞] where the lower limit, ξb = 0, is the Chaplygin–Lamb dipole and ξb = ∞ is a pair of point vortices.

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

49

So far, no one has yet taken up the challenge of calculating such conjectured “bipolar modons”. Morse and Feshbach (pp. 1212–1213 of [52]) solve the potential flow of an inviscid, incompressible fluid, uniform far upstream, past a pair of cylinders using the bipolar harmonic series given in Section 2.2 here. It is important to note that the cylinders need not be of equal diameters. (This series also applies to a uniform electric field with boundary conditions imposed on two cylinders.) Bipolar coordinates have proved equally useful in electromagnetics. Morse and Feshbach [52] give a pithy review on pp. 1210–1215. The potential due to a point charge at (d, 0) in Cartesian coordinates and an infinite conducting plane at x = 0 can be solved by Lord Kelvin’s method of images (p. 1180 of op. cit.); if an equal but opposite charge is placed at (−d, 0), the sum of the two point charge fields will cancel identically everywhere along the plane x = 0 and thus the sum solves the problem. The isolines of the potential are in fact the contour lines of constant ξ of a bipolar coordinate system! It is straightforward, in bipolar coordinates, to compute the electric potential outside a pair of cylinders even if the cylinders are of different diameters. Morse and Feshbach give explicit, closed-form solutions in elementary functions for the case where one cylinder is at zero potential while the other is at constant potential V ; an explicit solution for the capacitance is also provided. Magnetic field lines are treated, too. Jarem [39,40] is a more recent treatment of electromagnetic scattering from cylinders using bipolar coordinates. Liemert gives a good treatment of the Green’s function for the Poisson equation, motivated by various electrostatics applications [47]. Electromagnetic cloaking is analyzed by Chen and Chan [16]. Lucht gives a careful analysis of the two-cylinder capacitor [49]. All these applications have the following property Definition 2 (Bipolar-compatible domain). A domain is “bipolar-compatible” if all parts of its boundary ∂ coincide with a coordinate surface (i.e., a curve of constant ξ or η ) of a wisely-chosen bipolar coordinate system. When a domain is bipolar-compatible, this implies that one can always apply a tensor product grid with a tensor product basis using bipolar coordinates. This is good, but it is also restrictive: all boundaries in a “bipolar-compatible” domain must be circular boundaries. (In this context, an infinite straight line which is the limit of larger and larger circles is included in the category of “boundary circles”.) However, Scharstein and Davis apply bipolar coordinates to electromagnetic scattering from a semi-circular trough in a grounded plane [56]. (The scatterer is made by cutting a strip from the infinite plane and then bridging the gap with an infinitely long trough of semicircular cross-section.) The foci of the bipolar coordinates are at the points where the semicircular trough joins the infinite plane. In their matched asymptotic expansion, the zeroth order term is an exact solution to Laplace’s equation in the form of an integral in bipolar coordinates. The semicircular trough is not a bipolar-compatible domain. It is an example of a broader category that we may dub “bipolar-useful”. The restriction to circular boundaries and bipolar-compatible domains can be relaxed through a simple, smooth changeof-coordinate, known for brevity as a “coordinate mapping”. Orszag, for example, [54], mapped an annulus with undulating boundaries into a conventional annulus by the simple radial stretching from the polar coordinates (r , θ) to (ρ , η) where η ≡ θ and

ρ≡

( ro − ri ) + (1 −  )r ro (θ) − r i (θ)

(8)

where the radii of the undulating boundaries of the original domain are r i (θ) and ro (θ) and where the new computational annulus is bounded by the unit circle and by a concentric circle of radius  . If the mapped computational domain is an eccentric annulus, then bipolar coordinates with a Fourier–Chebyshev spectral basis can be usefully employed even though the original domain may have had boundaries far from circular. For further background on bipolar coordinates, we recommend the leisurely and very readable review by Philip Lucht [49] and the highly compressed but also useful treatment of Morse and Feshbach [52]. The outline of the rest of the paper is as follows. Section 2 is a discussion of separation of variables for solving partial differential equations (PDEs). It also treats bipolar harmonics, which are explicit, individual solutions to Laplace’s equation. Their usefulness can be extended from bipolar-compatible domains to geometries that are sufficiently small perturbations of bipolar-compatible domains by coordinate mapping or by “boundary collocation”. The latter approximates the Laplace solution by a series of bipolar harmonics. In separation of variables, the coefficients of the bipolar harmonics series are determined by interpolating the boundary conditions by the series. The “boundary collocation” generalization is needed when one or more boundaries is not circular. It is then usually necessary to replace interpolation with least-squares fitting: more points on the boundary than the number of bipolar harmonics in the truncated series. Ill-conditioning rapidly becomes a serious problem as the boundary is deformed more and more from a circle. Nevertheless, Section 3 shows a case where bipolar harmonic boundary collocation triumphs. The next section defines the Chebyshev–Fourier pseudospectral method and its advantages over alternatives. Section 5 shows that tensor product spectral methods are equally useful for the “sectorial annulus” except that the basis must be Chebyshev polynomials in both coordinates. The only significant flaw is that grid squares vary wildly in size. We analyze non-uniformity in Section 6 and show in the following section that one-dimensional coordinate mappings can reduce but

50

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

not eliminate this difficulty, and also analyze rates of convergence. Lastly, we show that the Chebyshev–Fourier pseudospectral method in bipolar coordinates is very effective for the Sverdrup jump problem in Stommel’s model of the wind-driven ocean circulation in spite of narrow boundary layers, high gradient internal layers and boundary layer detachment. The three appendices provide information for extending bipolar coordinates into three dimensions, for transforming differential equations from Cartesian coordinates to bipolar and lastly for understanding the connection between bipolar coordinates in the complex-plane and the isoconvergence contours of rational Chebyshev functions for approximating functions on the real line [6]. 2. Separability: solving Laplace’s equation through bipolar coordinates 2.1. Separability The Poisson equation is



(cosh(ξ ) + cos(η))2 d2

 ∂2 ∂2 + u = f (ξ, η) ∂ξ 2 ∂ η2

(9)

(Conversion from Cartesian to bipolar coordinates is discussed in Appendix A.) This can be rewritten in a form identical to the same equation in Cartesian coordinates except for the modified inhomogeneity.



 ∂2 d2 ∂2 + u= f (ξ, η) 2 2 ∂ξ ∂η (cosh(ξ ) + cos(η))2

(10)

The Helmholtz equation is



(cosh(ξ ) + cos(η))2 d2

 ∂2 ∂2 + u + κ 2 u = f (ξ, η) ∂ξ 2 ∂ η2

(11)

Multiplying through by the same factor as for the Poisson equation gives



 ∂2 ∂2 κ 2 d2 d2 + u+ u= f (ξ, η) 2 2 2 ∂ξ ∂η (cosh(ξ ) + cos(η)) (cosh(ξ ) + cos(η))2

(12)

Unfortunately, it is not possible to apply separation of variables to the Helmholtz equation. 2.2. Solving Laplace’s equation through bipolar harmonics In bipolar coordinates, the Laplace equation has the same form as in Cartesian coordinates: u ξ ξ + u ηη = 0. (However,

η is a cyclical, angular coordinate that requires periodic boundary conditions while ordinary Dirichlet, Neumann or Robin boundary conditions are imposed on the two circles ξ = ξo and ξ = ξi .) The “bipolar harmonics” are thus the same as for Laplace’s equation on a rectangle. We shall therefore simply give the classic, textbook solution [52]. A useful preliminary is to note that any function can be split into its parts which are symmetric and antisymmetric with respect to η = 0 by

v (ξ, η) = v S + v A ,

v S = (1/2){ v (ξ, η) + v (ξ, −η)}, v A = (1/2){ v (ξ, η) − v (ξ, −η)}

(13)

It is convenient to list the symmetric and antisymmetric solutions separately. The general symmetric solution is

v S = a0 + b 0 ξ +

∞ 

[am exp(−mξ ) + bm exp(mξ )] cos(mη)

(14)

m =1

The coefficients are found by expanding the boundary conditions as Fourier series and matching each term in the Fourier series. Suppose that boundary functions for Dirichlet boundary conditions have the Fourier series in angle

v S (ξi , η) = ainner + 0

∞ 

inner am cos(mη)

(15)

outer am cos(mη)

(16)

m =1

v S (ξo , η) = aouter + 0

∞  m =1

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

51

The symmetric solution is then

v S (ξ, η) = ainner 0

 ∞   ξ − ξo ξ − ξi inner sinh(m[ξ − ξo ]) outer sinh(m[ξ − ξi ]) + aouter + a a cos(mη) + m m 0 ξi − ξo ξo − ξi sinh(m[ξi − ξo ]) sinh(m[ξo − ξi ])

(17)

m =1

The antisymmetric solution is very similar with sines replacing cosines, and also with the constant and the linear function in ξ replaced by η and ηξ :

v A = A0η + B 0ξ η +

∞ 

[ Am exp(−mξ ) + B m exp(mξ )] sin(mη)

(18)

m =1

The coefficients are found by expanding the boundary conditions as Fourier series and matching each term in the Fourier series. Suppose the antisymmetric parts of the boundary functions have the sine series

v A (ξi , η) =

∞ 

inner bm sin(mη)

(19)

outer bm sin(mη)

(20)

m =1

v A (ξo , η) =

∞  m =1

η = 0, the Laplace solution is  inner sinh(m[ξ − ξo ]) outer sinh(m[ξ − ξi ]) bm sin(mη) + bm sinh(m[ξi − ξo ]) sinh(m[ξo − ξi ])

Specializing to functions with antisymmetry about

v A (ξ, η) =

∞ 



m =1

2.3. Monotonic quasi-radial variation of individual terms in the inner and outer series It is useful to split the symmetric solution into two parts:

v S (ξ, η) = ainner 0 v inner, S =

∞ 

ξ − ξo ξ − ξi + aouter + v inner, S + v outer, S 0 ξi − ξo ξo − ξi

inner am

m =1

v outer, S =

∞ 

outer am

m =1

sinh(m[ξ − ξo ]) sinh(m[ξi − ξo ]) sinh(m[ξ − ξi ]) sinh(m[ξo − ξi ])

(21)

cos(mη)

cos(mη)

(22)

Observe that

v inner, S (ξ = ξi ) =

∞  m =1

v inner, S (ξ = ξo ) = 0

inner am

sinh(m[ξ − ξo ]) sinh(m[ξi − ξo ])

cos(mη)

(23) (24)

At the boundaries, each series is identical with the corresponding boundary Fourier series. As we pass from ξ = ξi to ξo , the coefficients of the inner series, sinh(m[ξ − ξo ])/ sinh(m[ξi − ξo ]), decrease monotonically until all coefficients are identically zero at the outer boundary. Likewise, the coefficients of the outer series vary monotonically from zero on the inner boundary to matching the Fourier coefficients of the outer boundary series. This monotonic variation from one boundary to the other is a property of an individual term, and only after the solution has been split into an inner series and an outer series. The sum of the series may exhibit very complex non-monotonic behavior. Nevertheless, the property of “term-monotonicity” is very useful in analyzing rates of convergence as we shall do later. 3. Boundary collocation An alternative strategy for solving the Laplace equation in a complicated domain is to apply the species of harmonics associated with a simpler geometry, such as using harmonics for a disk (“polar harmonics”) when the domain is actually bounded by an ellipse or is an eccentric annulus, and then determine the coefficients by interpolation on the boundaries. Once the coefficients of a harmonic series is known, one can evaluate the solution at an arbitrary point merely by summing the series [6]. This strategy is called “boundary collocation” [28]. Any complete set of Laplace solutions can be applied; the

52

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

Fig. 2. Polar harmonic boundary collocation. The solution was approximated by forty terms in a polar harmonic cosine series with least-squares boundary collocation. The difference between the harmonic series and the Laplacian solution at sixty points on each boundary circle was minimized in a least-squares sense. The radius of the outer boundary is one; the radius of the circular island is  = 1/2. The center of the island is shifted by s = 0.3 left of the center of the outer boundary circle.

Method of Fundamental Solutions (MFS) is an alternative to boundary collocation that employs a sum of Green’s functions, each with its singularity slightly outside the boundary [5], as the Laplace solutions. Often more points than harmonics are used with the interpolation conditions satisfied only in a least squares sense. One could, for example, solve problems in an eccentric cylinder domain by expanding the solution as a series of polar harmonics, which are Laplace solutions for the interior of a disk. For simplicity only, assume solutions symmetric about the x-axis. Using a polar coordinate system centered at the center of the outer boundary circle,

u polar = a N log(r ) +

N −1 

an r n−1 cos(nθ)

(25)

n =1

Unfortunately, boundary collocation with polar harmonics is a disaster for the example illustrated. Fig. 2 shows that the computed values are an order of magnitude too high with a ring of spikes around the outer boundary. The least-squares problem, even with three times as many boundary points was so ill-conditioned that it was necessary to use an SVD decomposition of the interpolation matrix, discarding all modes with singular values whose magnitude, relative to the largest singular value, was smaller than 10−13 . But this SVD filtering only contracted the size of the disaster without removing the unacceptability of the result. Extreme ill-conditioning is generic for boundary collocation, but a more fundamental problem is that the domain of convergence of the polar harmonics series is the largest disk, centered on the center of the polar coordinate system, which is singularity-free. Therefore, the disk of convergence of the polar harmonics series here is limited by the logarithmic singularity of the solution at x = 1, y = 0. The center of the polar coordinate system is here fixed at the center of the outer boundary circle at x = 1.4, y = 0. Thus infinite series converges only within the disk (x − 1.4)2 + y 2 < 0.4. The attempt to apply boundary conditions on the outer circle of unit radius, far outside the radius of convergence of the polar harmonics, is absurd. In contrast, the logarithmic solution has no pathologies on the bipolar domain, and the bipolar harmonic series therefore converges exponentially fast. Fig. 3 shows that bipolar harmonics yield geometrically converging series for both the inner and outer contours for an arbitrary but typical solution to Laplace’s equation,

u = log((x − 1)2 + y 2 )

(26)

The monotonic decrease of each individual term in the inner and outer series from one boundary to the other, demonstrated in Section 2.3, shows that the rate of convergence is worst at the boundaries, and is therefore geometric but faster on the interior. Because of this fast truncation, truncating the series to a moderate number of terms gives very high accuracy. Nevertheless, boundary collocation — that is, using a species of Laplace harmonics for the wrong geometry but subject to the right boundary conditions — is still used [35,62,4,3,2,38,36,37]. The reason is that the method is simple; it is accurate, too, for domains which are small perturbations of the domain that defines the harmonics. Thus, boundary collocation works well — with bipolar harmonics instead of polar harmonics — if offset circular islands were replaced by islands bounded by ellipses of small eccentricity. A case study is illustrated in Fig. 4 and Table 2 below. Homogeneous boundary conditions are applied on an outer circle of unit radius and on an inner boundary which is an ellipse parameterized by

x = xC +  cos(t ),

y =   sin(t )

(27)

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

53

Fig. 3. Same as previous figure except that bipolar harmonics were used instead of polar harmonics. Isolines of the solution are shown in the upper panel; inner outer coefficients am and am are shown on a logarithmic scale in the lower panel.

Fig. 4. The striped shading shows the domain for an eccentric annulus with the inner circle replaced by an ellipse [green]. The parameters used in boundary collocation, appearing in the bipolar harmonics series, are  = 1/2, s = 3/10, xC = 1.1000, ξi = 1.42542, ξo = 0.86701 and ratio of axes in the ellipse are 3/2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 2 Maximum pointwise errors for the whole domain for a case study of boundary collocation with bipolar harmonics. N inner

N b,inner

N outer

Error ∞

30 30 30 30 60 90 60

30 60 90 120 60 90 120

30 30 30 30 30 30 30

18.2 1.4E−6 5.4E−9 5.3E−9 36.8 758 9.5E−9

which is equivalent to the implicit specification as the zero isoline of B (x, y ) = (x − xC )2 + y 2 /2 −  2 . The parameters ξo and ξi are identical to those for the eccentric annulus with given values of s and  where for the case illustrated s = 3/10 and  = 1/2. The ellipse has  = 3/2 so that the semi-major axis, parallel to the y-axis is half again as long as the semi-minor axis. For simplicity, we chose an exact solution which is symmetric about η = 0. The outer bipolar harmonic series can be computed independently of the inner series by calculating the Fourier coefficients of v (ξ = ξo , η) by trigonometric interpolation; it is sufficient for the number of collocation points N b,outer to equal the number of outer harmonics we choose, N outer . Denote the outer harmonic series by v outer :

v outer (ξ, η) = aouter 0

N outer −1 ξ − ξi outer sinh(m[ξ − ξi ]) + am cos(mη) ξo − ξi sinh(m[ξo − ξi ]) m =1

Define the inner harmonic basis by



φk (ξ, η) =

ξ −ξo ξi −ξo , sinh(m[ξ −ξo ]) , sinh(m[ξi −ξo ])

k=1 m = k − 1,

k>1

(28)

54

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

The inner collocation points are

x j = xC +  cos(t j ), y j =   sin(t j ), t j = π

j−1 N b,inner − 1

, j = 1, 2, . . . N b,inner

(29)

Define a vector whose elements are

f j = v (x j , y j ) − v outer (ξ(x j , y j ), η(x j , y j )), j = 1, 2, . . . N b,inner

(30)

The elements of the boundary collocation matrix are

H jk = φk (x j , y j )

(31)

Let a inner denote the vector of coefficients of the φ j . These can be found by solving the matrix problem

a inner = f H

(32)

The complete boundary collocation solution is then (with truncation) identical in form to the general symmetric-in-η solution given previously:

v S (ξ, η) = ainner 0

N b,inner  ξ − ξo ξ − ξi inner sinh(m[ξ − ξo ]) + aouter + am cos(mη) 0 ξi − ξo ξo − ξi sinh(m[ξi − ξo ]) m =1

N b,outer



+

outer am

m =1

sinh(m[ξ − ξi ]) sinh(m[ξo − ξi ])

cos(mη)

(33)

Note that the restriction to symmetry has been made only for simplicity. The antisymmetric solution can be computed independently of the η -symmetric solution but in the same manner. The method also generalizes to when the island is not symmetric about y = 0, but it is then necessary to lump the symmetric and antisymmetric components together and solve a single big matrix problem to impose the boundary conditions on the island. Table 2 shows the errors in the approximation of the exact Laplace equation solution. The exact solution for the case study is

v (x, y ) = log([1 − x]2 + ( y − 1/3)2 ) + log([1 − x]2 + ( y + 1/3)2 )

(34)

which has point singularities inside the island. Interpolation always fails; this is common in boundary collocation. However, when the number of points is larger than the number of harmonics, and the matrix problem is solved in a least-squares sense — the Matlab backslash operator does this automatically by the QR algorithm for rectangular matrices — near-machine precision accuracy is possible. Bipolar coordinates greatly expand the range of “small perturbations” for boundary collocation. A full treatment of boundary collocation requires a separate article, but this case suffices to illustrate the possibilities. 4. Chebyshev–Fourier pseudospectral algorithm for the eccentric annulus Because the angular bipolar coordinate η is cyclic, all functions are periodic in η with period 2π . It is natural to use a Fourier series to represent the angular dependence. Since the dependence in ξ is not periodic, it is natural to use a Chebyshev polynomial basis for the quasi-radial coordinate. The minor technical complication is that the argument of the Chebyshev polynomials must be stretched from the canonical Chebyshev interval z ∈ [−1, 1] to ξ ∈ [ξi , ξo ]. Let the argument of the Chebyshev polynomials be z. The forward and inverse mappings are

ξ≡ z≡

ξi − ξo 2 2

z+

ξi + ξo

2   ξi + ξo ξ−

ξi − ξo 2ξ − ξo − ξi = ξi − ξo

2

(35) (36) (37)

The Chebyshev polynomials are the images of Fourier cosine basis functions under a change of coordinate as expressed by the Chebyshev identity

T m (cos(t )) = cos(mt )

(38)

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

55

It is often convenient to represent the function as a fully trigonometric series:

u (x, y ) =

N M  

amn cos(mt [ξ ]) cos(nη) +

m =0 n =0

where

 t (ξ ) = arccos

bmn cos(mt [ξ ]) sin(nη)

ξi − ξo 2dx d2 + x2 + y 2

 (40)





η = arctan

,

2dy



d2 − x2 − y 2

The grid is a tensor product grid, uniform in the bipolar angular coordinate dinate t:

η j = 2π

2j −1 2N + 2

(39)

m =0 n =1

2ξ − ξi − ξo

 ξ = arctanh

N M  

, j = 1, . . . N + 1; t k = π

2k − 1 2M + 2

, k = 1, . . . M + 1

(41)

η and also in the radial trigonometric coor(42)

The Fourier–Chebyshev method for the disk is extensively discussed in Boyd and Fu [14] and (for the closely related sphere) by Livermore, Jones and Worland [48] and the papers they cite. 4.1. Annulus versus disk: logarithmic singularities and the absence of parity conditions The annulus is both simpler and more complex than a disk. The Chebyshev–Fourier series of an analytic function on a disk must satisfy so-called “parity conditions”. If we write in polar coordinates

f (r , θ) =

N  n =0

f n (r ) cos(nθ) +

N 

gn (r ) sin(nθ)

(43)

n =1

then regularity at r = 0 requires that (i) all even coefficients f 2n (r ) and g 2n (r ) for integer n must be symmetric functions of radius r while all the odd coefficients f 2n−1 (r ) and g 2n−1 (r ) are of odd parity, that is, f 2n−1 (−r ) = − f 2n−1 (r ). If the domain does not extend to r = 0 — if the geometry is an annulus — then the parity conditions disappear. The bad news is that when a disk near the center of the domain is excluded from the domain (i.e., the island), there is no reason to exclude a singularity of the solution within the island. Poisson’s equation is obviously a very simple PDE, but in an annulus (whether eccentric or centered), a small island of radius  is a singular perturbation. The effect of the island is not O ( ) when   1, but rather O (1/ log( )). For example, if the domain is the annulus between an island of radius  and a “mainland coast” which is a circle of unit radius with homogeneous boundary conditions on both circles and the inhomogeneous term in the differential equation is the constant (−1), the problem and exact solution are

ψxx + ψ y y = −1, ψ(r = 1) = 0, ψ(r =  ) = 0    1 log(1/r ) ψ= 1 − r 2 − (1 −  2 ) 4 log(1/ )

(44) (45)

The dependence on 1/ log( ) is explicit. As  → 0, the location of the singularity of log(r ) at r = 0 moves closer and closer to the endpoint r =  ; the Chebyshev series in radius √ therefore converges slower and slower as the island size  decreases. (The n-th Chebyshev coefficient is proportional to ( 8 )n for   1 [9].) Thus, the singularity slows the rate of geometric convergence for the Chebyshev–Fourier series, but does not prevent an exponential rate of convergence. 5. Sectorial annulus A “sectorial annulus” is similar to an annulus except that the domain does not completely wrap around the inner circle, but instead is restricted to

η ∈ [η1 , η2 ],

η 2 − η1 < 2 π

(46)

A Fourier expansion in η is now a bad idea because the solution in a sectorial annulus need not be a periodic function of For the sectorial eccentric annulus, the best choice of spectral basis is Chebyshev in both coordinates.

η.

6. Non-uniformity and the Jacobian of the transformation “We suspect, however, that the role of Schwarz–Christoffel mapping in mesh generation is not likely to expand in the future.” [Tobin A. Driscoll and Lloyd N. Trefethen, p. 108 of [24]]

56

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

Fig. 5.√Error in interpolation of four typical functions by polynomial interpolation at the√roots of the Chebyshev polynomial of degree 20. f 1 (x) = (1 + √ 2 2 x/2)/ 1 − x2 ; f 2 (x) = cos(10(x + π /30)); f 3 (x) = (1 + x/2)/ 2 √ approximation √− x ; f 4 (x) = (1 + x/2)/ 5 − x . All four are analytic everywhere on the interval x ∈ [−1, 1] except that f 1 is singular proportional to 1/ 1 − x at the right endpoint and a weaker branch point proportional to 1 + x at the left endpoint.

Fig. 6. Expansion of f = cos(4π x) cos(4π y ) [left graph] in the region bounded by the bipolar coordinate surfaces ξi = 2.2924 and ξo = 0.8670 as a Chebyshev series in ξ combined with a Fourier cosine series in η . Right: errors in the spectral expansion using 31 Chebyshev polynomials in ξ and 102 cosine functions in η . The error is very small everywhere that is green; other colors show that the error is very non-uniform. The Jacobian of the transformation from Cartesian coordinates to bipolar coordinates varied from 0. 027 to 6.0, a factor of 225. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

One of the most desirable features of spectral methods is their uniformity: the error in a Fourier or Chebyshev series oscillates on a short spatial scale, but the amplitude of the oscillations (“envelope”) does not greatly vary except when the function f (x) represented by the Chebyshev series is a singular function. Fig. 5 shows the Chebyshev error for four representative f (x). The function f 1 is strongly singular at x = 1 with a weaker branch point at x = −1, so it is not surprising that the error has a large spike near the right endpoint and a small spike near the left endpoint. The errors for the functions smooth on the interval show no strong concentration. A comprehensive theory of Chebyshev and Fourier interpolation error envelopes is given in [7]. Unfortunately, bipolar coordinates induce a considerable degree of non-uniformity. Fig. 6 shows the errors in the Chebyshev–Fourier expansion of f = cos(4π x) cos(4π y ). Although f (x, y ) oscillates uniformly in all directions, the error is concentrated in only part of the domain. An infinitesimal cell of area d A is changed by a ratio equal to the Jacobian J of the transformation, which is the square of the metric factor h = d/(cosh(ξ ) + cos(η)):

d A = dxdy

(47)

2

= h dξ dη =

(48) d

2

(cosh(ξ ) + cos(η))2

dξ dη

(49)

Nearly-uniform resolution requires nearly-uniform area which in turn requires a nearly uniform Jacobian. Unfortunately, the Jacobian can be very non-uniform. The ratio of the maximum of the Jacobian to its minimum on a curve of constant ξ , which we shall dub the “non-uniformity ratio”, is

U(ξ ) ≡

maxη∈[−π ,π ] ( J )(ξ ) minη∈[−π ,π ] ( J )(ξ )

(50)

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

57

Fig. 7. Isolines of the base-10 logarithm of the error norm (maximum pointwise error) for Chebyshev–Fourier interpolation of f = cos(4π x) cos(4π y ) in bipolar coordinates, plotted versus the radius of the island  and the leftward shift s of the center of the island from the center of the larger circular boundary circle, which is of unit radius. The white region is (s +  ) ≥ 1 where the boundary circles overlap and domain is not an eccentric annulus.

=

(cosh(ξ ) + 1)2 (cosh(ξ ) − 1)2

≈ (16 ξ −4 +

(51)

  16 −2 26 ξ + + O ξ2 ) 3 45

ξ 1

(52)

Bipolar coordinates are an example of a conformal mapping of a rectangular domain to an eccentric annulus. There was a brief flurry of interest is using conformal mapping as a universal tool to generate orthogonal grids for complicated domains. However, this agenda spawned many disappointments, prompting the gloomy words of Driscoll and Trefethen in their book on conformal mapping quoted above. The villain is non-uniformity. Conformality is so restrictive that there is only limited freedom to control grid density. Wildly varying grid density is not a vice peculiar to bipolar coordinates. Most of the separable, tensor product coordinate systems are also conformal mappings and share the curse of non-uniformity. Fig. 7 shows how the interpolation error varies with the radius of the inner boundary  and the translation of the center of the “island” leftward from the center of the outer, larger boundary circle. Computations were restricted to (s +  ) ≤ 1 because otherwise the boundary circles overlap and a continuous ring no longer exists. For a fixed basis truncation, the figure shows that the error worsens slightly with increasing island size for fixed shift s of the island, but the error grows rapidly when s is increased for any island size. A very off-center annulus — a very eccentric annulus — is a stiff numerical challenge. Fig. 8 shows the Chebyshev–Fourier tensor grid in different circles. We can see that a small radius of the inner circle and large eccentricity of the circles results in a non-uniform grid. 7. One-dimensional mappings 7.1. Tensor-preserving transformations A further change of coordinates, known for brevity as a “mapping”, can reduce grid non-uniformity. Bivariate mapping functions need not destroy the tensor product character of the grid; indeed, the Cartesian-to-bipolar transformation is an example of a tensor product-preserving mapping with bivariate functions. However, one-dimensional transformations are obviously simpler and suffice to smooth the grid. Denoting the new computational coordinate by ζ , a one-dimensional mapping of the form η = f (ζ ) always preserves the tensor product property. The standard Chebyshev–Fourier grid in ξ − η may be rather non-uniform in both directions. In the two subsections, we therefore treat first a change-of-coordinate for the angular coordinate η , and then follow this with the introduction of a new quasi-radial coordinate. 7.2. A periodicity-preserving mapping of the angular coordinate Since η is cyclic, it is obviously desirable to use a mapping in the angular coordinate that does not wreck the convergence of the η -Fourier series. Boyd developed a comprehensive theory of changes-of-coordinates that (i) respect the periodicity of an angular coordinate and (ii) concentrate grid points where needed. One useful map is the “arctan/tan map” [8]



η j = 2 arctan L tan



π [1 − (2 j − 1)/(2N )] 2



,

j = 1, . . . N

(53)

Fig. 9 shows the grid distributions for different choices of the arctan/tan map parameter L (in η ) and both ξ and exp(ξo − ξ ) as quasi-radial coordinate. Changes in both coordinates have visible effects on the uniformity of the grid.

58

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

Fig. 8. Twenty grid points in angle η were combined eight different points in radius were used to generate each of these tensor product grids in bipolar coordinates. The non-uniformity of Chebyshev–Fourier grid is evident, both a small inner circle and a large eccentricity of the circles results in non-uniform grid.

The effects of the arctan/tan map parameter L on the interpolation error are shown in Table 3. To isolate the effects of the mapping of the angular coordinate, we intentionally employed a large number of interpolation points in the quasi-radial direction. Error is minimized when L = 2. Fig. 9 also shows that the grid distribution for L = 2 is more uniform. 7.3. A logarithmic radial coordinate We have defined ξ as in the usual definition of bipolar coordinates. However, for large ξ , the radii of the circles of constant ξ are R / sinh(ξ ) ≈ ( R /2) exp(−ξ ). This exponential shrinking of the coordinate circles degrades the rate of convergence for many functions. We therefore also used the quasi-radial coordinate ρ where

ξ = ξo − log(ρ )

⇔

ρ = exp(ξo − ξ )

(54)

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

59

Fig. 9. Grid distributions for different arctan/tan map parameters L. Both radial coordinates are shown: ξ as the black stars and exp(ξo − ξ ) as the red squares. The island offset is s = 0.4, d = 1 and the relative island radius is  = 0.2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 3 The column shows the variation of the interpolation error for a representative function, f (x, y ) = cos(4π x) cos(4π y ), as the map parameter L of the arctan/tan mapping of the angular η is varied. Two different quasi-radial coordinates were tried, but the first two digits of the errors were identical for both radial coordinates, so only a single column of errors is presented. The number of points on the tensor product grid was fixed at 80 points in the ξ direction, 100 points in the η direction, s = 0.4, d = 1 and  = 0.2. L

Either ξ coordinate

1 2 3 4 5 6 10

7.8e−11 2.7e−14 2.9e−14 2.9e−14 9.3e−13 1.2e−9 2.3e−4

Next we fixed the angular map parameter L and increased the number of points M in ξ to compare the influence of different ξ coordinates. Table 4 lists the interpolation error for different M. We can see that the coordinate exp(ξo − ξ ), instead of ξ , offers faster convergence. 8. Convergence rate of bipolar Fourier series One of the great virtues of Chebyshev and Fourier expansions on an interval or a rectangle is spatial uniformity [7,63]. Unfortunately, tensor product spectral expansions in bipolar coordinates can be rather non-uniform as we now demonstrate with case studies. For example, the bipolar harmonic expansion of x, a function of (ξ, η) equal at all points to the usual Cartesian coordinate x at that point, is given explicitly in [52]:

60

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

Table 4 Variation of the interpolation error of f (x, y ) = cos(4π x) cos(4π y ) as the number of points M in ξ was varied; the number of points in η was fixed at 100. The other fixed parameters were s = 0.4, L = 5, d = 1 and  = 0.2. The errors flatten out for large M due to hitting the “round off plateau” explained in [9]; the corrupted values are in italics and would be smaller if higher precision arithmetic were used. All calculations used the default sixteen decimal digit precision in Matlab. M

ξ coordinate

exp(ξo − ξ ) coordinate

20 25 30 35 40 45 50 60

7.9e−3 2.8e−4 7.3e−6 1.4e−7 2.1e−9 2.1e−11 9.3e−13 9.3e−13

1.3e−4 4.2e−7 1.1e−9 1.4e−12 9.3e−12 9.3e−13 9.3e−13 9.3e−13

Table 5 Maximum pointwise errors over the entire domain as s varies. The approximations were computed by interpolation with M = N = 15, that is, with fifteen Chebyshev polynomials in ρ and fifteen cosine functions in the angular coordinate. The upper half of the table used the unmapped angle η while the lower half shows the L ∞ errors using the arctan/tan map with L = 2 in the coordinate mapping. f (x)

s = 0.1

s = 0.2

s = 0.3

s = 0.4

s = 0.5

Unmapped η coordinate x2 y2 exp(−x2 ) exp(− y 2 ) exp(−x2 − y 2 ) sech( y ) x2 sech(2 y ) cos(4π x) cos(4π y )

2.4114e−13 2.1294e−13 1.0175e−13 1.0863e−11 2.5535e−15 6.2481e−12 3.7159e−09 7.2885e−06

1.1022e−08 9.5496e−09 3.8307e−09 4.3707e−07 4.5707e−11 2.336e−07 9.0927e−05 0.17665

7.252e−06 6.1748e−06 1.5537e−06 0.00021076 1.6911e−08 9.934e−05 0.020625 1.2972

0.00093566 0.00077243 8.0473e−05 0.014619 2.7984e−07 0.006042 0.65841 2.9452

0.059835 0.046622 0.0013982 0.27035 1.8521e−05 0.11292 7.5507 2.7706

3.5905e−09 7.1133e−09 4.5232e−09 1.5819e−08 1.8211e−08 7.2886e−09 3.3275e−08 0.0055695

1.0784e−09 1.4402e−09 2.5212e−09 1.1897e−08 1.0087e−08 5.5882e−09 2.605e−05 0.25299

3.5973e−10 3.3957e−10 3.9934e−08 1.9288e−05 2.1578e−09 1.724e−05 0.02241 1.5141

1.911e−05 1.7152e−05 8.3983e−06 0.016108 9.6897e−09 0.0065766 1.6891 2.7973

L = 2 in arctan/tan mapping of x2 y2 exp(−x2 ) exp(− y 2 ) exp(−x2 − y 2 ) sech( y ) x2 sech(2 y ) cos(4π x) cos(4π y )

η coordinate

x = dsign(ξ )

1+2

3.6389e−09 1.0448e−08 5.0261e−09 1.5846e−08 1.914e−08 7.3464e−09 8.1499e−09 0.00099233

∞ 

 m

(−1) [exp(−m|ξ |)sign(ξ )] cos(mη)

(55)

m =1

A linear univariate polynomial is the most trivial bivariate function one could imagine; it is also a solution to the Laplace equation. When ξ > 0, the series converges exponentially fast. On the symmetry axis x = 0, however, the rate of convergence is terrible. The coefficients are not decreasing at all. The alternating signs indicate the singularity producing this slow rate of convergence must be at η = ±π . (For a detailed analysis of the connection between the signs of a spectral series and singularity location, see [10].) Observe that as ξ tends to zero, the corresponding coordinate isoline is a circle of increasing radius, and thus variation over a finite range in the angle η is mapped to a longer and longer arc; in the limit ξ = 0, the η Fourier series must approximate the sum of the spectral series on the entire y-axis. This slow convergence as ξ → 0 is one negative characteristic of bipolar series. Even when the geometry is an annulus so that the interesting range of ξ is bounded away from zero by ξo , the value of ξ on the outer boundary, the rate of convergence is still sensitive to geometry. Table 5 shows that for a wide range of functions, error grows dramatically as the island is shifted by a larger distance s from the center of the outer circle. 9. Sverdrup boundary layers attached to a circular island in a circular ocean The Stommel problem is a sort of counterpoint to previous sections because the whole point is that the Stommel solution is very non-uniform over the domain.

δ(ψxx + ψ y y ) + ψx = −1

(x, y ) ∈ ,

ψ = 0 on ∂

(56)

where δ  1 is nondimensional bottom friction and the constant on the right is the curl of the wind stress, idealized as uniform over the ocean basin. The domain and its boundary ∂ vary from ocean to ocean and with the degree of

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

61

Fig. 10. Numerical solution of the Stommel flow for a circular ocean with an offset circular island, driven by a uniform curl of the wind stress, on two different Chebyshev–Fourier grids. The upper right plot shows the spectral coefficients computed by using the ξ − η grid while the lower right plot is the same but using the “arctan/tan map” (L = 2) and the ρ quasi-radial coordinate. The upper left plot shows the numerical approximation to the Stommel flow on the 60 × 60 grid. The lower left panel shows the maximum pointwise error [error in the L ∞ norm] on an m × m grid as a function of m. The errors are computed in the usual way, when an analytical exact solution is unavailable, by treating a higher resolution numerical solution as the “exact” solution.

idealization; here, the domain is restricted to a circular island in a circular ocean, the eccentric annulus geometry. (This is idealized but oceanographically illuminating.) A “Sverdrup” balance applies over most of the ocean,

ψx = −1

(57)

with the homogeneous boundary condition applied only on the eastern boundary of a latitudinal strip of ocean. Stommel showed that the western shores of an ocean basin had a boundary layer (“Stommel layer”) of thickness O (δ) where the flow adjusts rapidly so as to satisfy the homogeneous boundary condition on the western boundary, too. The “westernintensified” boundary layer jets are the simplest models for the Gulf Stream, Kuroshio and Alghulas currents. Only later it was realized that islands and peninsulas also disrupt Sverdrup flow. Just north of the northernmost tip of an island or peninsula, the boundary condition is applied far from the island on the eastern shore of the enclosing ocean. Slightly to the south, the boundary condition must be applied on the island or peninsula itself. The Sverdrup flow is discontinuous. The singularity is smeared out in a boundary layer, but one of complex structure. East of the tip, the boundary layer is a Stommel layer, but one that becomes thicker and thicker as the tip is approached. West of the tip, the boundary layer detaches from the coast and erupts westward into the interior of the ocean. The physics is diffusive with longitude as the time-like variable, and the layer width steadily increases westward from the island. Sverdrup jump layers and Stommel layers are a severe numerical challenge when the friction δ  1. A circular island in a circular ocean is obviously an idealization, but not atypical of oceanographic theory and related numerical experiments. Bipolar coordinates allow inexpensive but very accurate Chebyshev–Fourier numerical solutions in [34] and here in Fig. 10. In the case illustrated, the (nondimensional) radius of the outer circular ocean basin is one while the radius of the inner circular island radius is  = 1/5 and the off-center shift of the island is s = 1/2. We imposed ψ = 0 on both the ocean basin and island and applied the Chebyshev–Fourier spectral method in bipolar coordinates to solve the Stommel flow with δ = 1/100. (The coordinate system allows the island to be offset from the center of the surrounding annular ocean.)

62

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

We employed the Chebyshev–Fourier pseudospectral method with two different coordinate systems, the same pair of options discussed earlier: (i) ξ − η and (ii) the “arctan/tan map” in η with ρ = exp(ξo − ξ ) as the quasi-radial coordinate. The three pseudocolor plots in Fig. 10 were generated by employing a 60 × 60 pseudospectral grid. The line graph shows how the error in the L ∞ norm varies with grid size. The ρ quasi-radial coordinate gives slightly smaller errors than the ξ coordinate. 10. Summary and open problems Bipolar coordinates made it possible to apply an efficient Chebyshev–Fourier spectral method to island-bound boundary layers as detailed earlier and in our forthcoming article [34]. We also show through a mixture of review and original analysis that spectral-friendly “bipolar-compatible” domains can be extended by boundary collocation and changes-of-coordinates to a much wider class of “bipolar useful” geometries. Most useful-but-not-compatible domains are small perturbations of compatible domains. We also provide guidelines for applying Chebyshev and Fourier spectral methods to the eccentric annulus, sectorial annulus and other bipolar-compatible geometry. Although exotic, bipolar coordinates are actually rather easy to use. A caveat is that although changes-of-coordinates (mappings) and boundary collocation and so on have greatly expanded the reach of single-domain spectral methods, spectral algorithms fail when the domain is too distorted. It is far easier to map to an elliptical island offset from the center of a roughly-circular ocean from bipolar coordinates than from polar coordinates because the physical domain is a large perturbation of an annulus, but only a small perturbation of an eccentric annulus. Much future research using bipolar coordinates will therefore focus on nearly bipolar domains. Acknowledgements This work was supported by the National Science Foundation through grants OCE 0451951 and OCE 1059703. Zhu Huang was supported by NSFC grant key program (No. 51236006) of China. Appendix A. Derivative transformations Repeated application of the chain rule gives the formulas necessary to rewrite in bipolar coordinates equations given in Cartesian coordinates.

ux = uy =

u ξ cosh (ξ ) cos (η) + u ξ + sinh (ξ ) sin (η) u η

(A.1)

d − sinh (ξ ) sin (η) u ξ + u η + u η cosh (ξ ) cos (η)

(A.2)

d

1 u xx = − 2 − (cosh (ξ ))2 uηη − 2 (cosh (ξ ))2 sin (η) uη cos (η) d

− uξ ξ (cosh (ξ ))2 (cos (η))2 + (cosh (ξ ))2 (cos (η))2 uηη − 2 sinh (ξ ) sin (η) uξ η cosh (ξ ) cos (η) + sinh (ξ ) uξ cosh (ξ ) − sin (η) uη cosh (ξ ) − 2 uξ ξ cosh (ξ ) cos (η) − 2 cosh (ξ ) sinh (ξ ) uξ (cos (η))2 − 2 sinh (ξ ) sin (η) uξ η

− uξ ξ + sin (η) uη cos (η) + uηη − (cos (η))2 uηη − sinh (ξ ) uξ cos (η) 1 u y y = 2 (cosh (ξ ))2 (cos (η))2 uηη − 2 (cosh (ξ ))2 sin (η) uη cos (η) d

− uξ ξ (cosh (ξ ))2 (cos (η))2 + (cosh (ξ ))2 uξ ξ − 2 sinh (ξ ) sin (η) uξ η cosh (ξ ) cos (η) + sinh (ξ ) uξ cosh (ξ ) + 2 uηη cosh (ξ ) cos (η) − 2 cosh (ξ ) sinh (ξ ) uξ (cos (η))2 − sin (η) uη cosh (ξ ) + uξ ξ (cos (η))2 − 2 sinh (ξ ) sin (η) uξ η + sin (η) uη cos (η) − sinh (ξ ) uξ cos (η) + uηη − uξ ξ



Appendix B. Extensions to three dimensions Bipolar coordinates can be extended to three dimensions in three ways: Bicylindrical, toroidal and bispherical coordinates [52]. All are compatible with tensor product spectral algorithms even for nonseparable partial differential equations. However, within this group of three coordinate sets, the Laplace equation is separable only in toroidal coordinates and the Helmholtz equation is separable in none.

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

63

Appendix C. Bipolar coordinates and rational Chebyshev functions TB(x; L ), a basis for an unbounded interval The rational Chebyshev functions TB(x; L ), introduced by Boyd [6] as an elaboration of a strategy invented by Grosch and Orszag [31], are a good spectral basis for problems on x ∈ [−∞, ∞]. These functions are defined by

TB(x; L ) ≡ cos(n t [x]),

x = L cot(t ),

n = 0, 1 , 2 , . . .

(C.1)

In the same way that the domain of convergence of a power series is a disk in the complex plane, the Fourier series of a function f (t ) converges within the largest strip in the complex t-plane, symmetric about the real t-axis, which is free of singularities [i.e., poles, branch points, and essential singularities of f (x)]. The domain of convergence of the rational Chebyshev series in the x-plane is the image of the strip of convergence under the mapping that transforms a Fourier series in t to a rational Chebyshev series in x. Theorem 1 of [6] asserts the domain of convergence is the exterior of a bipolar coordinate contour line, ξ = ±ξc for some ξc . The required coordinate system is rotated by ninety degrees from that employed here. For the TB convergence isolines, the foci are at (x) = ± L where L is the “map parameter” appearing in the cotangent mapping. In general, more than one sing ξ contour has a pole or branch point of f (x); let us denote these contours ξ j . The convergence-limiting contour is sing

ξc = min |ξ j j

|

(C.2)

The convergence domain for the TB series of f (x) is thus the exterior of the largest bipolar coordinate isoline which is free of singularities of the function f (x). Engineering problems on an unbounded interval usually have singularities at infinity. The generic case is therefore convergence only on the real x-axis which is the exterior of a constant ξ circle in the limit ξ → 0. This is the isoline ξ = 0 in the rotated bipolar coordinate system. The rate of convergence is “subgeometric” in the sense defined in [6,9]. However, Alfimov, Usero and Vaázquez [1] and Boyd and Xu [13] independently showed that the solitary wave of the cubically-nonlinear Benjamin–Ono equation is singular only at a pair of complex-conjugate singularities on the imaginary axis. [As known much earlier, the solitary wave of the usual, quadratically-nonlinear Benjamin–Ono equations has the same property: it is the reciprocal of a quadratic polynomial with roots that are pure imaginary.] The rational Chebyshev series converge everywhere in the complex plane excluding two disks bounded by bipolar coordinate surfaces. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]

[11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26]

G.L. Alfimov, D. Usero, L. Vazquez, On complex singularities of solutions of the equation Hu x − u + u p = 0, J. Phys. A, Math. Gen. 33 (2000) 6707–6720. L. Badea, P. Daripa, A fast algorithm for two-dimensional elliptic problems, Numer. Algorithms 30 (2002) 199–239. L. Badea, P. Daripa, A fast algorithm for two-dimensional elliptic problems, Numer. Algorithms 30 (2003) 199–239. L. Badea, P. Daripa, On a Fourier method of embedding domains using an optimal distributed control, Numer. Algorithms 32 (2003) 261–273. T. Betcke, L.N. Trefethen, Reviving the method of particular solutions, SIAM Rev. 47 (2005) 469–491. J.P. Boyd, Spectral methods using rational basis functions on an infinite interval, J. Comput. Phys. 69 (1987) 112–142. J.P. Boyd, The envelope of the error for Chebyshev and Fourier interpolation, J. Sci. Comput. 5 (1990) 311–363. J.P. Boyd, The arctan/tan and Kepler–Burger mappings for periodic solutions with a shock, front, or internal boundary layer, J. Comput. Phys. 98 (1992) 181–193. J.P. Boyd, Chebyshev and Fourier Spectral Methods, Dover, New York, 2001, 680 pp. J.P. Boyd, A proof, based on the Euler sum acceleration, of the recovery of an exponential (geometric) rate of convergence for the Fourier series of a function with Gibbs Phenomenon, in: J.S. Hesthaven, E.M. Ronquist (Eds.), Spectral and Higher Order Methods for Partial Differential Equations: Proceedings of ICOSAHOM ’09 Conference, Trondheim, Norway, June 22–26, in: Lecture Notes in Computational Science and Engineering, vol. 76, Springer-Verlag, New York, 2011, pp. 131–140. J.P. Boyd, H. Ma, Numerical study of elliptical modons by a spectral method, J. Fluid Mech. 221 (1990) 597–611. J.P. Boyd, E. Sanjaya, Geometrical effects on western intensification of wind-driven ocean currents: the rotated-channel Stommel model, coastal orientation, and curvature, Dyn. Atmos. Ocean. 65 (2014) 17–38. J.P. Boyd, Z. Xu, Numerical and perturbative computations of solitary waves of the Benjamin–Ono equation with higher order nonlinearity using Christov rational basis functions, J. Comput. Phys. 231 (2012) 1216–1229. J.P. Boyd, F. Yu, Comparing six spectral methods for interpolation and the Poisson equation in a disk: radial basis functions, Logan–Shepp ridge polynomials, Fourier–Bessel, Fourier–Chebyshev, Zernike polynomials, and double Chebyshev series, J. Comput. Phys. 230 (2011) 1408–1438. Z. Cao, K. Luo, H. Yi, T.H. P, Energy conservative dissipative particle dynamics simulation of mixed convection in eccentric annulus, Int. J. Heat Mass Transf. 74 (2014) 60–76. H. Chen, C.T. Chan, “Cloaking at a distance” from folded geometrics in bipolar coordinates, Opt. Lett. 34 (2009) 2649–2651. J. Chen, M. Tsai, C. Liu, Conformal mapping and bipolar coordinate for eccentric Laplace problems, Comput. Appl. Eng. Educ. 17 (2009) 314–322. V. Chernyavsky, Exact solution for creeping cylindrical flow in a free-dowel bearing, Dokl. Phys. 53 (2008) 19–22. C. Cho, K. Chang, K. Park, Numerical simulation of natural convection in concentric and eccentric horizontal cylinder annuli, ASME J. Heat Transf. 104 (1982) 624–630. R. Churchill, J. Brown, R. Verhey, Complex Variables and Applications, McGraw–Hill Book Company, New York, 1976. R. Dai, Q. Dong, A. Szeri, Flow between eccentric rotating cylinders: bifurcation and stability, Int. J. Eng. Sci. 30 (1992) 1323–1340. D.N. de Garrs Allen, Relaxation Methods in Engineering and Science, McGraw–Hill, 1954. R. Diprima, J. Stuart, Flow between eccentric rotating cylinders, Trans. ASME J. Lubr. Technol. 94 (1972) 266–274. T.A. Driscoll, L.N. Trefethen, Schwarz–Christoffel Mapping, Cambridge Monographs on Applied and Computational Mathematics, vol. 8, Cambridge University Press, 2002, 150 pp. L. Elliot, D. Ingham, T. Bashir, Stokes flow past two circular cylinders using a boundary element method, Comput. Fluids 24 (1995) 787–798. A. Eydeland, B. Turkington, A computational method of solving free-boundary problems in vortex dynamics, J. Comput. Phys. 78 (1988) 194–214.

64

Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64

[27] E.E. Feldman, R.W. Hornbeck, J.F. Osterle, A numerical solution of laminar developing flow in eccentric annular ducts, Int. J. Heat Mass Transf. 78 (1982) 194–214. [28] B.A. Finlayson, The Method of Weighted Residuals and Variational Principles, 2nd ed., SIAM, New York, 2013. [29] J. Frene, M. Godet, Flow transition criteria in a journal bearing, Trans. ASME J. Lubr. Technol. 96 (1974) 135–140. [30] W. Graebel, Advanced Fluid Mechanics, Elsevier, Amsterdam, 2007. [31] C.E. Grosch, S.A. Orszag, Numerical solution of problems in unbounded regions: coordinate transforms, J. Comput. Phys. 25 (1977) 273–296. [32] B. Hamrock, Fundamentals of Fluid Film Lubrication, McGraw–Hill Book Company, New York, 1994. [33] J.F. Heyda, A Green’s function solution for the case of laminar incompressible flow between non-concentric circular cylinders, J. Franklin Inst. 267 (1959) 25–34. [34] Z. Huang, J.P. Boyd, Smoothing the Sverdrup flow discontinuity, J. Phys. Oceanogr. (2015), in preparation. [35] S.Z. Husain, J.M. Floryan, Spectrally-accurate algorithm for moving boundary problems for the Navier–Stokes equations, J. Comput. Phys. 229 (2001) 2287–2313. [36] S.Z. Husain, J.M. Floryan, Gridless spectral algorithm for Stokes flow with moving boundaries, Comput. Methods Appl. Mech. Eng. 198 (2008) 245–259. [37] S.Z. Husain, J.M. Floryan, Immersed boundary conditions method for unsteady flow problems described by the Laplace operator, Int. J. Numer. Methods Fluids 56 (2008) 1765–1786. [38] S.Z. Husain, J.M. Floryan, Implicit spectrally-accurate method for moving boundary problems using immersed boundary conditions concept, J. Comput. Phys. 227 (2008) 4459–4477. [39] J.M. Jarem, Rigorous coupled wave analysis of bipolar cylindrical systems: scattering from inhomogeneous dielectric material, eccentric, composite circular cylinders, Prog. Electromagn. Res. 43 (2003) 181–237. [40] J.M. Jarem, Rigorous coupled wave analysis of bipolar cylindrical systems: scattering from inhomogeneous dielectric material, eccentric, composite circular cylinders, J. Electromagn. Waves Appl. 18 (2004) 73–74. Extended abstract. [41] G.B. Jeffery, The rotation of two circular cylinders in a viscous fluid, Proc. R. Soc. Lond. Ser. A 101 (1922) 169–174. [42] N. Joukowski, Motion of a viscous fluid contained between rotating eccentric cylindrical surfaces, in: Proceedings of the Khar’kov Mathematical Society, 1887, pp. 34–37. [43] N. Joukowski, S. Chaplygin, Friction of a lubricated layer between a shaft and its bearing, Tr. Otd. Fiz. Nauk Obshch. Lyub. Estest. 13 (1904) 24–36. [44] M. Kamal, Separation in the flow between eccentric rotating cylinders, Trans. ASME J. Basic Eng. 88 (1966) 717–724. [45] S. Kang, Characteristics of flow over two circular cylinders in a side-by-side arrangement at low Reynolds numbers, Phys. Fluids 15 (2003) 2486–2498. [46] E. Kulinski, S. Ostrach, Journal bearing velocity profiles for small eccentricity and moderate modified Reynolds numbers, Trans. ASME J. Appl. Mech. 34 (1967) 16–22. [47] A. Liemert, The Green’s function of the Poisson equation on the non-concentric annular region, J. Electrost. 72 (2014) 347–351. [48] P.W. Livermore, C.A. Jones, S.J. Worland, Spectral radial basis functions for full sphere computations, J. Comput. Phys. 227 (2007) 1209–1224. [49] P. Lucht, Bipolar coordinates and the two-cylinder capacitator, available at http://user.xmission.com/~rimrock/Documents/Bipolare, 71 pp. [50] F. Mahfouz, Natural convection within an eccentric annulus at different orientations, J. Thermophys. Heat Transf. 26 (2012) 665–672. [51] A.J. Moghadam, A.B. Rahimi, A singular perturbation solution of viscous incompressible fluid flow between two eccentric rotating cylinders, Arab. J. Math. 3 (2014) 63–78. [52] P.M. Morse, H. Feshbach, Methods of Theoretical Physics, 2nd ed., Feshbach Publishing, Wellesley, Massachusetts, 2005, 2000 pp. (in two volumes), reissue of book published in 1953. [53] M. Nakanishi, T. Kida, Unsteady low Reynolds number flow past two rotating circular cylinders by a vortex method, in: Proceedings of the 3rd ASME/JSME Joint Fluid Engineering Conference, San Frantsisco, California, USA, July 1999, 2000. [54] S.A. Orszag, Spectral methods for problems in complex geometries, J. Comput. Phys. 37 (1980) 70–92. [55] L.F. Richardson, The approximate arithmetical solution by finite differences of physical problems involving differential equations with an application to stresses in a masonry dam, Philos. Trans. R. Soc. Lond. A 210 (1910) 307–357. [56] R.W. Scharstein, A.M.J. Davis, Matched asymptotic expansion for the low-frequency scattering by a semi-circular trough in a ground plane, IEEE Trans. Antennas Propag. 48 (2000) 201–211. [57] V. Sennitskii, On the propelling of a pair of rotating circular cylinders in a liquid, Din. Sploš. Sredy 47 (1980) 145–153. [58] V. Sennitskii, On the drag force acting on a pair of circular cylinders streamed by water, Din. Sploš. Sredy 52 (1981) 178–182. [59] M. Shaarawi, H. Abualhamayel, E. Mokherimer, Developing laminar forced convection in eccentric annuli, Heat Mass Transf. 33 (1999) 353–362. [60] S. Smith, The rotation of two circular cylinders in a viscous fluid, Mathematika 38 (1991) 63–66. [61] S. Surattana, M. Nikolay, Numerical simulation of steady viscous flow past two rotating circular cylinders, Suranaree J. Sci. Technol. 13 (2006) 219–233. [62] J. Szumbarski, J.M. Floryan, A direct spectral method for determination of flows over corrugated boundaries, J. Comput. Phys. 153 (1999) 378–402. [63] L.N. Trefethen, Approximation Theory and Approximation Practice, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, 2013. [64] S. Wang, An experimental and numerical study of natural convection heat transfer in horizontal annuli between eccentric cylinders, J. Therm. Sci. 4 (1995) 38–43. [65] E. Watson, The rotation of two circular cylinders in a viscous fluid, Mathematika 42 (1995) 105–126. [66] W.W. Wood, The asymptotic expansions at large Reynolds numbers for steady motion between non-coaxial rotating cylinders, J. Fluid Mech. 3 (1957) 159–175.