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 eﬃcient cartography for a variety of geometries: the exterior of two disks or cylinders, a halfplane 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 quasiradial coordinate ξ gives easytoprogram 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 nonuniform, and the rate of exponential convergence increasingly slow. Onedimensional coordinate mappings signiﬁcantly reduce the nonuniformity. In spite of this nonuniformity, the Chebyshev–Fourier method is quite effective in an idealized model of the winddriven 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 “bipolarcompatible” domains listed above, but is a suﬃciently 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 onedimensional domains. This allows the construction of multivariate basis sets and grids as tensor products of onedimensional 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 semianalytical method reigned supreme for solving partial differential equations although there was some application of ﬁnite 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 ramiﬁcations of applying Chebyshev–Fourier spectral methods in bipolar coordinates. The motive was an application to idealized ocean ﬂow [34,12] as illustrated in Section 9. Numerous numerical issues
*
Corresponding author. Email address:
[email protected] (J.P. Boyd).
http://dx.doi.org/10.1016/j.jcp.2015.03.062 00219991/© 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 deﬁne 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 Nonuniformity ratio: the maximum, on a circle of ﬁxed ξ , of the Jacobian function divided by the minimum Relative radius of inner boundary circle in an eccentric annulus; the absolute radius is R Alternative quasiradial coordinate Quasiradial 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 quasiradial coordinate ξ (left panel) and quasiangular 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 deﬁned for all (x, y ), but in this application, we employ it only for the ﬂow 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 inﬁnity. ξ is everywhere negative in the left halfplane 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 eﬃciency 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 deﬁnition 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 yaxis.) 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 halfplane or halfspace. Inﬁnite plane with semicircular trough. Any domain that can be smoothly mapped into one of the domains described above.
Deﬁnition 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 xaxis where this axis is chosen to pass through the centers of both circles, then the eccentric annulus can be equivalently deﬁned in bipolar coordinates as follows with speciﬁcation 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 quasiradial 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 ﬂow and heat transfer in an eccentric circular annulus have also been studied by [33,27,59,15,64,50,19]. The axialcoordinateindependent ﬂow 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 orecrushers. Solutions to PDEs in the halfplane, x ∈ [0, ∞] ⊗ y ∈ [−∞, ∞], are discussed in many textbooks; Cartesian coordinates are usually satisfactory for problems in a halfplane. 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 ﬂow between a long, straight stretch of coast and an offshore island is an obvious application. A modon is a pair of contrarotating vortices, translating at a constant speed under their mutually interacting tangential velocities as a longlived coherent structure. They are also called “Lamb–Chaplygin dipoles” after the two ﬂuid dynamicists who discovered these exact analytic twodimensional 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 ﬂuid 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 ﬁnite, 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 halfdisks; outside the union of the two semidisks, the ﬂow 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 ellipsebounded modons computed by Ma and Boyd [11] that modons exist in which the ﬂow 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 ﬂow of an inviscid, incompressible ﬂuid, 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 ﬁeld 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 inﬁnite 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 ﬁelds 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, closedform 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 ﬁeld 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 twocylinder capacitor [49]. All these applications have the following property Deﬁnition 2 (Bipolarcompatible domain). A domain is “bipolarcompatible” if all parts of its boundary ∂ coincide with a coordinate surface (i.e., a curve of constant ξ or η ) of a wiselychosen bipolar coordinate system. When a domain is bipolarcompatible, 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 “bipolarcompatible” domain must be circular boundaries. (In this context, an inﬁnite 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 semicircular trough in a grounded plane [56]. (The scatterer is made by cutting a strip from the inﬁnite plane and then bridging the gap with an inﬁnitely long trough of semicircular crosssection.) The foci of the bipolar coordinates are at the points where the semicircular trough joins the inﬁnite 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 bipolarcompatible domain. It is an example of a broader category that we may dub “bipolaruseful”. The restriction to circular boundaries and bipolarcompatible domains can be relaxed through a simple, smooth changeofcoordinate, 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 bipolarcompatible domains to geometries that are suﬃciently small perturbations of bipolarcompatible 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 coeﬃcients 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 leastsquares ﬁtting: more points on the boundary than the number of bipolar harmonics in the truncated series. Illconditioning 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 deﬁnes 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 signiﬁcant ﬂaw is that grid squares vary wildly in size. We analyze nonuniformity in Section 6 and show in the following section that onedimensional coordinate mappings can reduce but
50
Z. Huang, J.P. Boyd / Journal of Computational Physics 295 (2015) 46–64
not eliminate this diﬃculty, 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 winddriven 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 complexplane 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 modiﬁed 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 coeﬃcients 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 coeﬃcients 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 quasiradial 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 coeﬃcients of the inner series, sinh(m[ξ − ξo ])/ sinh(m[ξi − ξo ]), decrease monotonically until all coeﬃcients are identically zero at the outer boundary. Likewise, the coeﬃcients of the outer series vary monotonically from zero on the inner boundary to matching the Fourier coeﬃcients 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 nonmonotonic behavior. Nevertheless, the property of “termmonotonicity” 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 coeﬃcients by interpolation on the boundaries. Once the coeﬃcients 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 leastsquares boundary collocation. The difference between the harmonic series and the Laplacian solution at sixty points on each boundary circle was minimized in a leastsquares 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 satisﬁed 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 xaxis. 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 leastsquares problem, even with three times as many boundary points was so illconditioned 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 ﬁltering only contracted the size of the disaster without removing the unacceptability of the result. Extreme illconditioning 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 singularityfree. 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 ﬁxed at the center of the outer boundary circle at x = 1.4, y = 0. Thus inﬁnite 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 deﬁnes 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 ﬁgure except that bipolar harmonics were used instead of polar harmonics. Isolines of the solution are shown in the upper panel; inner outer coeﬃcients 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 ﬁgure 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 speciﬁcation 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 semimajor axis, parallel to the yaxis is half again as long as the semiminor 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 coeﬃcients of v (ξ = ξo , η) by trigonometric interpolation; it is suﬃcient 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
Deﬁne 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)
Deﬁne 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 coeﬃcients 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 symmetricinη 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 leastsquares sense — the Matlab backslash operator does this automatically by the QR algorithm for rectangular matrices — nearmachine 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 suﬃces 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 quasiradial 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 socalled “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 coeﬃcients f 2n (r ) and g 2n (r ) for integer n must be symmetric functions of radius r while all the odd coeﬃcients 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 nth Chebyshev coeﬃcient 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. Nonuniformity 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 nonuniform. 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 ﬁgure 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 nonuniformity. 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 inﬁnitesimal 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)
Nearlyuniform resolution requires nearlyuniform area which in turn requires a nearly uniform Jacobian. Unfortunately, the Jacobian can be very nonuniform. The ratio of the maximum of the Jacobian to its minimum on a curve of constant ξ , which we shall dub the “nonuniformity 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 base10 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 ﬂurry 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 nonuniformity. 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 nonuniformity. 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 ﬁxed basis truncation, the ﬁgure shows that the error worsens slightly with increasing island size for ﬁxed shift s of the island, but the error grows rapidly when s is increased for any island size. A very offcenter 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 nonuniform grid. 7. Onedimensional mappings 7.1. Tensorpreserving transformations A further change of coordinates, known for brevity as a “mapping”, can reduce grid nonuniformity. Bivariate mapping functions need not destroy the tensor product character of the grid; indeed, the Cartesiantobipolar transformation is an example of a tensor productpreserving mapping with bivariate functions. However, onedimensional transformations are obviously simpler and suﬃce to smooth the grid. Denoting the new computational coordinate by ζ , a onedimensional mapping of the form η = f (ζ ) always preserves the tensor product property. The standard Chebyshev–Fourier grid in ξ − η may be rather nonuniform in both directions. In the two subsections, we therefore treat ﬁrst a changeofcoordinate for the angular coordinate η , and then follow this with the introduction of a new quasiradial coordinate. 7.2. A periodicitypreserving 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 changesofcoordinates 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 quasiradial 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 nonuniformity of Chebyshev–Fourier grid is evident, both a small inner circle and a large eccentricity of the circles results in nonuniform 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 quasiradial 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 deﬁned ξ as in the usual deﬁnition 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 quasiradial 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 ﬁgure 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 quasiradial coordinates were tried, but the ﬁrst 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 ﬁxed 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 ﬁxed the angular map parameter L and increased the number of points M in ξ to compare the inﬂuence 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 nonuniform 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 ﬁxed at 100. The other ﬁxed parameters were s = 0.4, L = 5, d = 1 and = 0.2. The errors ﬂatten 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 ﬁfteen Chebyshev polynomials in ρ and ﬁfteen 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 coeﬃcients 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 ﬁnite 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 yaxis. 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 nonuniform 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 ﬂow 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 coeﬃcients computed by using the ξ − η grid while the lower right plot is the same but using the “arctan/tan map” (L = 2) and the ρ quasiradial coordinate. The upper left plot shows the numerical approximation to the Stommel ﬂow 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 ﬂow adjusts rapidly so as to satisfy the homogeneous boundary condition on the western boundary, too. The “westernintensiﬁed” 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 ﬂow. 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 ﬂow 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 timelike 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 offcenter 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 ﬂow 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 quasiradial 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 ρ quasiradial coordinate gives slightly smaller errors than the ξ coordinate. 10. Summary and open problems Bipolar coordinates made it possible to apply an eﬃcient Chebyshev–Fourier spectral method to islandbound boundary layers as detailed earlier and in our forthcoming article [34]. We also show through a mixture of review and original analysis that spectralfriendly “bipolarcompatible” domains can be extended by boundary collocation and changesofcoordinates to a much wider class of “bipolar useful” geometries. Most usefulbutnotcompatible 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 bipolarcompatible geometry. Although exotic, bipolar coordinates are actually rather easy to use. A caveat is that although changesofcoordinates (mappings) and boundary collocation and so on have greatly expanded the reach of singledomain 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 roughlycircular 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 deﬁned 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 tplane, symmetric about the real taxis, 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 xplane 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 convergencelimiting 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 inﬁnity. The generic case is therefore convergence only on the real xaxis 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 deﬁned in [6,9]. However, Alﬁmov, Usero and Vaázquez [1] and Boyd and Xu [13] independently showed that the solitary wave of the cubicallynonlinear Benjamin–Ono equation is singular only at a pair of complexconjugate singularities on the imaginary axis. [As known much earlier, the solitary wave of the usual, quadraticallynonlinear 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. Alﬁmov, 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 twodimensional elliptic problems, Numer. Algorithms 30 (2002) 199–239. L. Badea, P. Daripa, A fast algorithm for twodimensional 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 inﬁnite 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, SpringerVerlag, 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 intensiﬁcation of winddriven ocean currents: the rotatedchannel 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 ﬂow in a freedowel 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 ﬂow past two circular cylinders using a boundary element method, Comput. Fluids 24 (1995) 787–798. A. Eydeland, B. Turkington, A computational method of solving freeboundary 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 ﬂow 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 ﬂow between nonconcentric circular cylinders, J. Franklin Inst. 267 (1959) 25–34. [34] Z. Huang, J.P. Boyd, Smoothing the Sverdrup ﬂow discontinuity, J. Phys. Oceanogr. (2015), in preparation. [35] S.Z. Husain, J.M. Floryan, Spectrallyaccurate 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 ﬂow 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 ﬂow problems described by the Laplace operator, Int. J. Numer. Methods Fluids 56 (2008) 1765–1786. [38] S.Z. Husain, J.M. Floryan, Implicit spectrallyaccurate 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 ﬂuid, Proc. R. Soc. Lond. Ser. A 101 (1922) 169–174. [42] N. Joukowski, Motion of a viscous ﬂuid 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 ﬂow between eccentric rotating cylinders, Trans. ASME J. Basic Eng. 88 (1966) 717–724. [45] S. Kang, Characteristics of ﬂow over two circular cylinders in a sidebyside arrangement at low Reynolds numbers, Phys. Fluids 15 (2003) 2486–2498. [46] E. Kulinski, S. Ostrach, Journal bearing velocity proﬁles for small eccentricity and moderate modiﬁed Reynolds numbers, Trans. ASME J. Appl. Mech. 34 (1967) 16–22. [47] A. Liemert, The Green’s function of the Poisson equation on the nonconcentric 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 twocylinder 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 ﬂuid ﬂow 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 ﬂow 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 ﬁnite 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 lowfrequency scattering by a semicircular 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 ﬂuid, Mathematika 38 (1991) 63–66. [61] S. Surattana, M. Nikolay, Numerical simulation of steady viscous ﬂow 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 ﬂows 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 ﬂuid, Mathematika 42 (1995) 105–126. [66] W.W. Wood, The asymptotic expansions at large Reynolds numbers for steady motion between noncoaxial rotating cylinders, J. Fluid Mech. 3 (1957) 159–175.