Composite Grid Construction and Applications Citation Chesshire, Geoffrey Samuel (1986) Composite Grid Construction and Applications. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/K5KS-4N46. https://resolver.caltech.edu/CaltechETD:etd-03242008-145553 Abstract For many problems that involve the discrete representation of a function on a two-dimensional region, a single rectangular grid is inappropriate. However, any compact region with a smooth boundary can be discretised using a composite grid, which provides boundary-fitted coordinates and allows stretched coordinates for boundary layers and internal layers wherever these are needed. We have developed interactive software for specifying composite grids, and general-purpose subroutines for using them in discretising PDE boundary-value problems. The discretisation of an elliptic BVP with an internal layer and of a supersonic flow problem are given as examples of the uses of composite grids and the capabilities of the software. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: Applied Mathematics Degree Grantor: California Institute of Technology Division: Physics, Mathematics and Astronomy Major Option: Applied Mathematics Thesis Availability: Public (worldwide access) Research Advisor(s): Kreiss, Heinz-Otto Thesis Committee: Kreiss, Heinz-Otto (chair) Defense Date: 15 May 1986 Funders: Funding Agency Grant Number Natural Sciences and Engineering Research Council of Canada (NSERC) UNSPECIFIED Department of Energy (DOE) DE-AA03-76SF0076 Office of Naval Research (ONR) N00014-83-K-0422 NSF ATM8201207 Record Number: CaltechETD:etd-03242008-145553 Persistent URL: https://resolver.caltech.edu/CaltechETD:etd-03242008-145553 DOI: 10.7907/K5KS-4N46 Related URLs: URL URL Type Description https://doi.org/10.1137/0908074 DOI Article adapted for Appendix 2. Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 1105 Collection: CaltechTHESIS Deposited By: Imported from ETD-db Deposited On: 04 Apr 2008 Last Modified: 21 Dec 2019 02:58 Thesis Files Preview PDF (Chesshire_gs_1986.pdf) - Final Version See Usage Policy. 5MB Repository Staff Only: item control page

COMPOSITE GRID CONSTRUCTION AND APPLICATIONS Thesis by Geoffrey Samuel Chesshire In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy California Institute of Technology Pasadena, California 1986 (Submitted May 15, 1986) Acknowledgements I would like to thank my advisor, Heinz Kreiss, for his many helpful suggestions and for his infinite patience. A large part of the credit belongs to Bill Henshaw, whose collaboration made this work possible. I was able to complete this thesis only with the steady encouragement and support of my family and the special friends I have met while at Caltech. 1 would also like to thank Gerald Browning, who arranged for me to work several summers at the National Center for Atmospheric Research. My work at Caltech has been supported by the Natural Science and Engineering Research Council of Canada through scholarships, the Department of Energy under contract DE-AA03-76SF0076, the Office of Naval Research under contract NO0014— 83-K-0422, and the National Science Foundation under contract ATM8201207. The computations were done on the Fluid Dynamics VAX 11/750 at Caltech. The figures were prepared using the NCAR graphics package and the manuscript was typeset using TpX. —ili- Abstract For many problems that involve the discrete representation of a function on a two-dimensional region, a single rectangular grid is inappropriate. However, any compact region with a smooth boundary can be discretised: using a composite grid, which provides boundary-fitted coordinates and allows stretched coordinates for boundary layers and internal layers wherever these are needed. We have developed interactive software for specifying composite grids, and general-purpose subroutines for using them in discretising PDE boundary-value problems. The discretisation of an elliptic BVP with an internal layer and of a supersonic flow problem are given as examples of the uses of composite grids and the capabilities of the software. Table of Contents Acknowledgments ....... 2.0.0 cece ence een ene e een e nent neta n nee eee ii ADStract oo eee cccc eeee t eee eeeen eee e ee eee eee e eee Eee lil Table of Contents 2.0... 0c ccc cette tenn e teen eee ete een e ee ne nee iv List of Figures...... 0.0 cece eee eee e eee teeter eee e nee vi List of Tables 0.00... ccc cece ee nee eee eee e eee e nee e nese eee e eens vi Introduction. 2... ee ee 1 Chapter 1 Curvilinear Grids 21... ce ee ee ees 4 Section 1.1 Notation for Curvilinear Grids ...... 2.00. ee eee ee eee 4 Section 1.2 PDEs on Curvilinear Grids 2.0... eee ee ee 7 Section 1.3 Curvilinear Grid Construction ...... 2... 0.02. ee eee 10 Section 1.3.1 Algebraic Methods of Grid Construction ......... 11 Section 1.3.2 Grid Construction Based on PDEs .............. 13 Chapter 2 Composite Grids 2.2... 0... eee tee 17 Section 2.1 Notation for Composite Grids ..... 2... ... 0.2.02 0008. 18 Section 2.2 PDEs on Composite Grids 1.2... ... 0... 0c eee eee 19 Section 2.3 Composite Grid Construction ........... 0.20.20 00008 24 Section 2.4 Composite Grids for Multigrid ...............000056. 28 Chapter 3 Composite Grid Software .. 0... . ee eee 30 Section 3.1 A Program for Composite Grid Construction ........... 30 Section 3.1.1 Construction of Curves ... 2... 00... eee eee 31 Section 3.1.2 Construction of Grids ..............--2 000004 34 Section 3.1.3 Construction of Composite Grids ............... 37 Section 3.2 Subroutines for Using Composite Grids ............... 42 Chapter 4 Applications of Composite Grids ............ 00.0000 000. 44 ~V~ Section 4.1 An Elliptic Problem With an Internal Layer ........... 44 Section 4.2 A Supersonic Flow Problem ............0-. eae 49 Appendix 1 CMPGRD: A Composite Grid Construction Program ........ 57 Section Al.1 Purpose... 0... ce ee eee 57 Section Al.2 Elements of a Composite Grid............. Lee eee 58 Section Al.2.1 First A Simple Example ................---- 59 Section Al.2.2 Regions in the Plane .... 2... .. 2. eee ees 60 Section A1l.2.3 Regions on a Cylinder or Torus ...........044. 64 Section A1l.3 PDEs on Composite Grids ........ 0.0000 eeee 65 Section Al.3.1 A Simple Example .......... 00.2 eee eee ees 65 Section A1.3.2 More General PDEs ......... 2.0.00 eee ees 67 Section Al.4 Using Program CMPGRD ........ 0... 2. ee ee ees 70 Section Al.4.1 Getting Started. 2... 0... ee ee 70 Section Al.4.2 Specifying Curves .... 0. eee ee ee 73 Section Al.4.3 Specifying Grids and Transformations........... 75 Section Al.4.4 Specifying Composite Grids ...............004 81 Section Al.4.5 Specifying Parameters ........ 0.00 cee eee eee 84 Section Al.56 Appendices ..... 0... 0c ee eens 85 Section Al.5.1 Subroutine CURVE ...........0...0...200.05. 85 Section A1.5.2 The Graphics Package... ...... cee ee 88 Section A1.5.3 The Composite Grid Data File................ 99 Appendix 2 Multigrid on Composite Meshes ........... 0000002 e eae 102 Section A2.1 Introduction .... 0... 0... ee 102 Section A2.2 Description, ©... 2... ee ee eee 104 Section A2.2.1 Multigrid. 2.0.0.0... 0... cece ees 106 Section A2.3 Numerical Results .. 20... 00.00. cee ee ee eee 111 References List of Figures Figure 2.1 2... eee ee ee Dee 21 Figure 4.1 2... ee eee 47 Figure 4.2 0. ees 50 Figure 4.3 0... ee eee eee 54 Figure 4.4 2... ee 55 Figure 4.5 2.0. ee eee ee 56 Figure ALL... ee eee eee 61 Figure Al.2 2... ee eens 61 Figure AL.3 2.0. ee eens 63 Figure A1l4 . 0... eee eee 78 Figure A2.1 2... ee ee eee 114 Figure A2.2 0... eee 115 List of Tables Table 2.1 0... eee eee 24 Table 2.2 eee ee ee ee eee 24 Table 4.1 2... ee ee eee ee 48 Table 4.2 0 eee eee eee es 52 Table A211... eee eee eee eee 112 Table A2.2 0 ee ne nee nets 113 Introduction For many problems that involve the discrete representation of a function on a two-dimensional region, a single rectangular grid is inappropriate. For example, a rectangular grid is poorly suited to the discretisation of a partial differential boundary-value problem on any region but a rectangle, as the boundary of the region will not in general correspond to a line of the grid. As a result, the accuracy of the solution may be significantly degraded, as demonstrated by Henshaw [12]. Many methods have been devised to construct a one-to-one transformation from a region bounded by piecewise smooth curves onto a rectangle. See, for example, the review article by Eiseman {9} and the conference proceedings edited by Thompson [22]. A cartesian coordinate system on the rectangle then corresponds to a curvilinear coordinate system on the region, and at least for the case of a simply-connected region, the boundary of the rectangle corresponds to the boundary of the region. A rectangular grid on the rectangle corresponds to a curvilinear grid on the region, and the boundary of the region corresponds everywhere to a gridline. One difficulty faced by all of these methods is that for a region of complex geometry they have very limited control over the distribution of coordinate lines, and hence of the gridlines. This difficulty is addressed with some success by the elliptic grid generation methods of Thompson and Mastin [15] and the adaptive grid generation methods of Brackbill and Saltzman [4]. However, even for a region of relatively simple geometry such as a simply-connected region consisting of two discs connected by a narrow strip, the density of gridlines will inevitably be greater on the narrow strip than on the discs. Another difficulty is that any transformation from a multiply-connected _2)~ region (other than a doubly-connected region) to a rectangle must have a singularity somewhere on the boundary. Both of these difficulties are consequences of using a single grid to cover the entire region. An alternative is to cover the region with a composite grid consisting of several overlapping curvilinear grids, each of which is the image, under a smooth trans- formation, of a rectangular grid, and chosen so that the boundary of the region is composed entirely of edges of the curvilinear grids. The density of gridlines on each curvilinear grid is independent of the other grids and is determined by the transformation specified for that grid. Using enough of these transformations and grids it is possible to discretise problems on complicated regions for which a sin- gle grid is hopelessly inadequate. A composite grid can be constructed for any compact region with a piecewise smooth boundary, even for a multiply-connected region. It provides boundary-fitted coordinates and allows stretched coordinates for boundary layers and internal layers wherever these are needed. The concept of a composite grid extends easily to regions with cylindrical or toroidal symmetry and to n-dimensional regions. We have developed software for the interactive construction of composite grids for compact regions in the plane and for regions with cylindrical or toroidal symme- try, and general-purpose subroutines for using them in discretising PDE boundary- value problems. These include subroutines for reading in a data file containing a description of a composite grid, for calculating the coefficients in the centered second-order discretisation of a linear combination of derivatives on a nine-point stencil, for interpolating data between the grids of a composite grid in the regions of their overlap, and for the second-order discretisation and direct solution of a second-order linear elliptic equation with mixed boundary conditions on a compos- ite grid. We have generalised the procedure of B. Kreiss [14] for construction of a pair of overlapping grids for a simply-connected region with a smooth boundary _3- so that we can construct a composite grid for any compact two-dimensional region with a piecewise-smooth boundary, even a multiply-connected region periodic in one or both directions. Further generalisations we have included allow for the use of stretched coordinates on the component grids to resolve boundary layers and in- ternal layers, the discretisation of a PDE and its boundary conditions to any order of accuracy, the interpolation of function values between the component grids to any order of accuracy, and the construction of sequences of composite grids for use with multigrid methods. The discretisations of an elliptic boundary value problem with an internal layer and of a supersonic flow problem with multiple shocks are given as examples of the uses of composite grids and the capabilities of the software. CHAPTER 1 Curvilinear Grids Numerical simulations in fluid mechanics and many other problems involving the computation of a solution field on a bounded region in the plane are often carried out on a rectangular computational grid chosen so that the boundary of the region follows lines of the grid. This is only possible for regions bounded by curves consisting of line segments intersecting at right angles. Although regions bounded by such curves are involved in many interesting fluid flow problems, numerical simulations are carried out on such regions less often for this reason than because of the ease of accurately applying boundary conditions on a grid fitted to the boundary of the region. For some applications it may not be important to use a boundary- fitted grid; a rectangular grid may suffice. However, boundary conditions for the discretisation of partial differential equations are difficult to apply when the grid is not boundary-fitted, and for a problem with a boundary layer, truncation error in the boundary condition results in an inaccurate solution. Examples of this effect are given in {12}. Boundary conditions for a partial differential equation may be accurately discretised on a region with a piecewise-smooth boundary if a boundary- fitted curvilinear grid can be found. For this reason the construction of such grids is of interest. 1.1 Notation for Curvilinear Grids We shal] consider only curvilinear grids in two dimensions. Although the con- cepts discussed here extend easily to three or more dimensions, the notation becomes _5- difficult and regions and their boundaries become difficult to describe or imagine. Since the principal use for curvilinear grids is the discretisation of PDE boundary- value problems by finite-difference methods, we shall consider only “logically rect- angular” grids, those on which the gridpoints can be labelled by two indices (, 7) so that gridlines drawn through all of the gridpoints with the same index 7 do not intersect, and likewise for gridlines of constant 7. Such curvilinear grids may be continuously deformed from a uniform rectangular grid. This leads to the following definition of a curvilinear grid. Definition 1.1.1 A curvilinear grid for a compact region D in the plane, on a cylinder or on a torus is a continuous transformation d from the unit square onto D, D = {d(r,s):0<r<1,0<s< I}, and a uniform rectangular grid G on the unit square, r—-1 g-Il ; . c= {(2 _ Jide is tenh, m—-in—-1l We shall henceforth say “grid” whenever it is clear from the context that we mean “curvilinear grid.” Just how smooth the transformation d needs to be will depend on the specific application the grid is to be used for. The boundary of a multiply-connected region, such as an annulus, consists of two or more disjoint curves. Under any continuous transformation d from the unit square onto such a region D, the boundary curves of D are the image of parts of the boundary of the unit square; other parts of the boundary of the unit square are mapped onto curves called “branch cuts” in the interior of D that connect the boundary curves. Branch cuts may also be introduced, as explained below, so that the cylinder and the torus may be represented in the plane. However, these need not add to the confusion if they are chosen to follow the above branch cuts wherever they intersect D. _6- A cylinder may be represented in the plane by a “period strip” bounded by a pair of parallel lines C, and C, which correspond to a single line, a branch cut, on the cylinder. If corresponding points on C, and Cy, are parametrised by ci(€) = a; + bé and c2{€) = a, + bé, any function on the cylinder corresponds to a periodic function in the plane, with periodicity cj. = a, — ag. It is sometimes convenient to choose a curve other than a straight line for the branch cut on the cylinder. In that case the period strip in the plane is bounded by curves C; and Cy on which corresponding points may be parametrised by functions cy and cy, which differ by the constant vector cy, = ¢i(€) — c2(€). For example, if the the image of the boundary of the unit square under the transformation d crosses the branch cut on the cylinder, the corresponding transformation from the unit square onto the period strip in the plane is discontinuous. To avoid this problem the branch cut must be chosen so that the image of the boundary of the unit square under the transformation d does not cross it; this is achieved by choosing the branch cut in the cylinder so that it follows a branch cut due to the transformation d wherever it intersects D. At the risk of some confusion, we use the same variables (x,y) to label coordinates on both the cylinder and the period strip in the plane, and use the same name for a function on the cylinder and the corresponding function on the period strip. A torus may likewise be represented in the plane by a “period box” bounded by a trapezoid. Each pair of opposite sides of the trapezoid corresponds to a single line segment, a branch cut, on the torus. More generally, if the branch cuts on the torus are chosen to be curves, then the period box in the plane is bounded by two pairs of curves (C,,C2) and (C3,C,4) on which corresponding points may be parametrised by two pairs of functions (¢,(€),c2(€}) and (c1(€),c2(€)) which differ by the constant vectors Cy. = c,(€) ~ c2(€) and c34 = ¢3(€) — c4(é). Any function on the torus corresponds to a doubly-periodic function in the plane with periods _7- Cy. and c34. The branch cuts on the torus should be chosen so that the image of the boundary of the unit square under the transformation d does not cross them. We use the same variables (x,y) to label coordinates on both the torus and the period box in the plane, and use the same name for a function on the torus and the corresponding function on the period box. 1.2 PDEs on Curvilinear Grids The principal use for curvilinear grids is in discretising PDE boundary-value problems using finite-difference methods on a region D whose boundary consists of smooth curves. To be useful for this purpose the grid must satisfy certain smooth- ness constraints that depend on the PDE. If the grid is constructed as the image under the transformation d of a rectangular grid, the constraints on the smoothness of the grid can be expressed as constraints on the transformation d. Consider a PDE boundary value problem 0 Oa where F represents the PDE in the interior of D and the boundary conditions on the boundary of D, and it is understood that higher derivatives may be involved. We can define a function u’ which represents u in cartesian coordinates on the unit square, u'(r,s) = u(d(r,s)) for (r,s) € [0,1]?. (1.2.2) If the Jacobian of the transformation d is non-singular, we can use it to express the problem F in (r,s) coordinates. Let J(r,s) = os = [ * r | , (1.2.3) Then the inverse of J(r,s) is 1 _ Jot 5¢ -——— | Ys el=(% wy. 12s (r 5) LrYs — Urts ~Ur Ly Sz By (1 2 4) _8~ With the understanding that by r,, ry, , and sy we mean the functions of (r,s) which form the inverse of the Jacobian, we can express the problem (1.2.1) in (r,s) coordinates as ] e] 6] i] F (125. + S25 ty a + S50’) -0 for (r,s) € [0,1]: (1.2.5) If the boundary-value problem involves at most first-order derivatives and rz, Ty, Sz and s, are continuous, the problem expressed in (r,s) coordinates is equivalent to the original problem. If the BVP involves nth derivatives, the transformation must satisfy the additional constraints ortid OrPds4 continuous for allp+q<n. (1.2.6) A problem may arise when parts of the boundary of the unit square do not correspond, under the transformation d, to the boundary of D. This occurs when part of the boundary of the unit square corresponds to a branch cut. If all of the required derivatives of d are continuous across the branch cut, the appropriate periodicity conditions are obvious. For example, if the branch cut is the image of the two sides r = 0 and r = 1 of the unit square and the Jacobian of d extends smoothly to a function periodic in r (with period 1) then the appropriate periodicity condition is that u’ must extend smoothly to a function periodic in r. It may be useful to note in passing that curvilinear grids may also be used for problems, such as the graphical representation of data on a region with a smooth boundary, where no differentiation is involved. For such problems the transformation may not need to be smooth at all; there may be no conditions at all on the Jacobian and its inverse. If a sufficiently fine uniform rectangular grid on the unit square and an appropri- ate difference method are used, an accurate discrete solution to the problem (1.2.1) may be obtained. The following second-order central difference approximations may be used for many problems: ot = 7,Do,u! + 8,Do.u' + O(Ar’, As’) x Ou ' ' 2 2 a, = ryDoru! + $y Dosu' + O(Ar*, As’) ¥ Oru 26-2, ' 27 772,,! 5n2 7 (r,)°(Déu') + 2r.8,(D,Du') + (s,)°( Diu’) x + [r2(Dorrz) + $2(Doerz)|(Dort'') + [r2(Dors2) + 82(Dos82)](Dosv') + O(Ar?, As?) 1.2.7 + [rz(Dorty) + 82(Dosty)|(Doru') + [rz(Dorsy) + 82(Doesy)|(Dosu') + O(Ar’, As’) O*u _ ay? = + [ry(Dorry) + Sy(Doery)|( Dor’) + [ry(Dorsy) + 8y(Dossy)|(Dosv’) (ry)°(Dyu') + 2rysy(D, Deu’) + (sy)’ (Deu) + O(Ar’, As’) where Dov (ris 83) = rien Si) — vrs 83) UIT, , $5 — Uv Ty, SG Dov (ris 83) = ( ae j 1) Dang.) x UlTivn 83) — 20(rin 83) + U(ti-1 85) Dyv(ri, 8;) — Ay? D.v(ri, $3) = Ast 1 1 Ar= — and As= ———~. mI n—-I The following second-order one-sided differences may in many instances be used instead of Do, and/or Do, for boundary conditions involving first derivatives: D_,v(r;,8;) = 2u(r;,8;) — 3u(ri-1,8;) + U(ti-25 85) Ar a vr; » 85 Q bo ; ye 2 7 ; Ds,v(ri,8;) = U(T42, 8) fens $3) — 2u(ri, 8) 2.9 2v(1;, 83) ~ 3u(r, S;-1) + v(r;, $;-2) (1 ) D_.v(ti, 8;) cs Ap D,.0(ri, $j) = —0(ri, Sj42) + Bv(Ti, 8341) — 2u(ri5 83) Ar -10- Periodicity conditions are easy to discretise on a branch cut that is the image under the transformation d of two opposite sides of the unit square and the Jacobian of the transformation extends smoothly to a function periodic in the appropriate coordinate, r or s. Periodicity is represented by the discrete condition that the index 2 or 7 is to be interpreted modulo m — 1 or n — 1: Periodicity inr: v(m, $;) = v(Ti4+m-15 55) (1.2.10) Periodicity ins: v(ri,8;) = v(t, 5;4n-1)- 1.3 Curvilinear Grid Construction Many methods have been devised to construct a one-to-one transformation from the unit square onto a bounded region D in the plane, under which a rectangular grid on the unit square corresponds to a boundary-fitted curvilinear grid on D. These methods generally fall into one of three classes: those which construct the transformation as the solution to a PDE (usually elliptic) with boundary conditions specifying that the edges of the unit square are mapped to the boundary of D, those which construct the transformation algebraically as a weighted sum of points on the boundary or of functions that describe the boundary (and in some cases other points or functions describing curves interior to D also), and those which construct the transformation from commonly known coordinate transformations by the appropri- ate choice of parameters. The latter methods readily provide transformations for such simple regions as a trapezoid or an annulus, but lack the versatility to handle regions of general shape. More sophisticated examples are the conformal Schwartz- Christoffel transformations for general polygonal regions, but these can be written in closed form only for a small class of regions. Generally, as the shape of the region becomes more complicated, these methods offer diminishing return for the effort of finding an appropriate transformation. Algebraic methods are more versatile in that they can provide transformations for regions whose boundaries need not be de- scribed in terms of simple functions. These methods can handle regions of relatively -11- complex geometry with a minimum of effort and computation. Methods based on PDEs are the most versatile and can be used for regions of quite general shape, but at the expense of the transformation being much slower to compute than one constructed by an algebraic method. Some PDE methods have the advantage of a maximum principle that ensures the resulting transformation is one-to-one. Some algebraic and PDE methods provide control, to a limited degree, over the density of the curvilinear gridlines. The following sections describe the most popular of these methods and some of their strengths and their limitations; Thompson [22] and Eiseman [9] give a more complete survey. 1.3.1 Algebraic Methods of Grid Construction Given a description of the boundary of D, and in some cases, curves in the interior of D also, the algebraic methods construct a transformation from the unit square onto D using various forms of interpolation. These methods are usually described as they apply to regions in three dimensions. However, to avoid notational complexity we shall consider only two-dimensional examples. All of these methods are quick to compute in comparison to the PDE methods. Some of them can handle regions of quite complicated shape. The simplest algebraic method is the shearing transformation d(r,s) = (1 — r)ei(s) + re2(s) from the unit square onto the region bounded by two curves C;, = {c;(s):0 <5 < 1} and the line segments joining their endpoints at s = 0 and at s = 1. This transformation may be used for any region D whose boundary consists of two line segments joining the endpoints of two curve segments, provided that functions c;, parametrising the curves C, can be found. If the curved part of the boundary of D is described only by discrete data points, the parametrisations c, may be constructed as spline functions through the data points. The main strength of this method is ~12- the speed with which the transformation may be computed, since it requires only two function evaluations. Its main weaknesses are that it can handle only regions with at least two straight sides, one opposite the other, and that the transformation resulting from a region with a convoluted boundary may not be one-to-one; it may have folds where its Jacobian derivative is singular. Most of the other algebraic methods are essentially generalisations of the shearing transformation, designed to overcome one or more of these weaknesses or to provide additional features. Hermite interpolation is a generalisation of the shearing transformation to allow control over derivatives of the transformation at the curved boundary segments of D. For example, d(r,s) = (1 — 8r? ~ 2r®)eyo(s) +. 77(3 — 2r)ego(s) + r(1 — r)?eqy(s) + 7?(1 — r)eoi(s) allows both d and @d/dr to be specified at r= O and r =1: d(0,s) = cyo(s), 2“(0.8) = €4;(s) d(1,s) = cgo(s), oc (1,3) = €91(s). This control over the derivative at the boundary is useful when coordinate transfor- mations for two adjacent regions must be found such that their Jacobian derivatives are identical along their common boundary. Hermite interpolation suffers the same weaknesses as the shearing transformation. Another generalisation of the shearing transformation, the “multisurface trans- formation” proposed by Eiseman [8], provides control over the transformation and its derivative on curves in the interior of D, so that od d(rz,8) = Cxo(s), By (Tee S) = Cz1(s) may be specified on arbitrary lines 0 = ry,7r2,...,ry = 1. The multisurface transfor- mation, like Hermite interpolation, is a linear combination of the functions ce, and C41 With polynomial coefficients; the coefficients of the polynomials are determined -13- by a linear system of equations. It maps the unit square onto the region bounded by Cy, Cy and two curves that pass through the points {czo(0):& = 1,...,N} and {cz:{1):k = 1,...,N}, respectively. It may be used for regions D of more complicated shape. Another advantage it has over the shearing transformation is additional control over the transformation (and over resolution in the r-direction) in the interior of D. A one-to-one multisurface transformation can be found in many cases where a one-to-one shearing transformation cannot. Still another generalisation of the shearing transformation is the “Coons patch,” (7| or “transfinite interpolation.” The simplest formulation of this, d(r,s) = (1 — s)ei(r) + se2(r) + (1 — r)eg(s) + rea(s) — (1 —s)(1 —r)c,(0) ~ (1 — s)re2(0) — s(1 — r)ex(1) — sree(1), maps the unit square onto the region D bounded by four curves Cy = {cxz{s):0 < s <1} whose parametrisations are chosen so that they intersect at c,(0) = c3(0), ci{1) = c4(0), c2(0) = c3(1) and c2(1) = c4(1). Like the shearing transformation, the simple Coons patch is fast to compute, requiring only four function evaluations. Its advantage is that it can handle a region whose boundary consists of four given curve segments. It has the same trouble with regions that have convoluted bound- aries; the resulting transformation may have folds where its Jacobian derivative is singular. A more general formulation of the Coons patch allows d(r,s) to be specified on more than two curves in each direction, d(ri,s) =ci(s) fort =1,...,M d(r,s;) =¢2,;{s) for 7 =1,...,N, with a specified number of derivatives continuous across these curves. This gener- alisation can handle regions with more complicated boundaries. 1.3.2 Grid Construction Based on PDEs PDE methods for construction of coordinate transformations from the unit square onto a region D are generally more versatile than algebraic methods in ~14- that they can handle regions of more complicated shape and that they offer greater control over resolution and skewness. This versatility comes at the expense of high computational cost. The most popular of these methods construct the transforma- tion as the solution to an elliptic PDE boundary-value problem. Other methods use hyperbolic PDEs (Starius [17], for example) or parabolic PDEs. The best-known elliptic method for construction of a coordinate transformation is conformal mapping. A conformal map can be found for any simply-connected region D in two dimensions, as the solution to the Cauchy-Riemann equations Or Os Os Or ax dy’ dx dy in the interior of D with boundary geometry (r,s) = W(z,y) specified, where © is a function (to be determined) that maps the boundary of D onto the boundary of the unit square. By differentiating and combining the Cauchy-Riemann equations, we get Laplace’s equation V?(r,s) = 0 in the interior of D. With the Cauchy- Riemann equations as additiona] boundary conditions these may alternatively be used in the interior of D. Conformal maps offer no control over resolution in the interior or on the boundary; gridlines nearby the boundary of D are always more dense where the boundary is concave than where it is convex. Conformal maps usually have singularities on the boundary of D where the Jacobian derivative becomes zero or infinite. For example, the conformal map from a circular disc onto the unit square is singular at the points corresponding to the corners of the unit square, and the inverse mapping is singular at the corners of the unit square. Furthermore, conformal mapping methods do not generalise to three dimensions. For these reasons conformal maps are not always appropriate. As a more useful method for generating coordinate transformations, Winslow [23| proposed using Laplace’s equations V7(r,s) = 0, as in the conformal mapping method, but relaxing the boundary conditions. The boundary condition (r,s) = W(x,y) is imposed, with the function W specified; the Cauchy-Riemann equations ~15- are not imposed at the boundary. Godunov and Prokipov [10] have shown that a maximum principle guarantees that the Jacobian of the resulting transformation is nonsingular. This approach has the advantage over conformal mapping that it allows some control over the transformation near the boundary. In particular the boundary conditions can be chosen so that corners of D correspond to corners of the unit square, so that the transformation need not have singuarities. Also, the transformation is faster to compute since it involves solving only a linear problem. This method generalises easily to three dimensions. A disadvantage is that the resulting grid is no longer orthogonal. The use of elliptic systems for generating coordinate transformations was fur- ther generalised by Thompson, Thames and Mastin [21] with the inclusion of forcing terms in the elliptic equations, and branch cuts to allow for multiply-connected re- gions. The terms f and g in the equations V?(r,s) = (f,g) allow much greater control over the density of gridlines, not only near the boundary but throughout the interior. Constraints on f and g have been found so that a maximum principle still guarantees that the Jacobian of the transformation is non-singular. The intro- duction of more than one branch cut, as neccessary for many multiply-connected regions, results in a region D with more than four corners. In such cases some of the corners of D cannot correspond to corners of the unit square, and the transfor- mation must have singularities there. Brackbill and Saltzman (4) generate coordinate systems using a variational for- mulation. They find the transformation that maximizes a functional J = 7,+A,/,+ Ati, where I, is a global measure of the smoothness of the transformation, J, is a global measure of resolution and Jj is a global measure of orthogonality. The Kuler equation of this variational problem is an elliptic system of which the system considered by Thompson is a special case. The variational formulation is used to specify the forcing functions in the elliptic problem. This method provides more -~16- systematic global control over the resolution and other properties of the grid. This systematic control makes possible the automatic adaptive generation of grids. —17- CHAPTER 2 Composite Grids Two limitations apply to the use of a single curvilinear grid on a region with a curved boundary. One is that even on a region with a fairly simple geometry the density of gridlines will inevitably be highly nonuniform. For example, a curvilinear grid for an oblong region with a narrow waist will inevitably have a much higher density of gridlines in the narrow waist than in some other parts of the region. The other limitation is that if the region D is not doubly-connected, any transformation must have a singularity somewhere on the boundary. For example, the Jacobian of any transformation from the unit square onto a disc (a simply-connected region) will be singular at some of the corners of the square; if polar coordinates are used, the Jacobian will be singular all along the side of the square corresponding to the center of the disc; if a non-periodic transformation which maps the corners of the square to points on the boundary of the disc is used, the Jacobian will be singular at all the four corners. Both of these problems can be overcome by the use of composite grids consisting of several overlapping curvilinear grids. Since the component grids of a composite grid may be much simpler than a single curvilinear grid for the same region, many of the methods of curvilinear grid construction discussed in the Chapter 1 are adequate to construct the component grids. In fact, any region may be broken up into enough sub-regions so that the shape of each is simple enough that even the least sophisticated methods may be used to construct the component grids, ~18- 2.1 Notation for Composite Grids We shall consider only composite grids in two dimensions, although all of the concepts extend to three or more dimensions. One-dimensional composite grids may be useful as models for higher-dimensional composite grids in the study of stability and convergence of PDE discretisations, but they have no practical advantage over a single stretched grid for the same interval. So we define composite grids as follows. Definition 2.1.1 A composite grid for a compact region D in the plane, on a cylin- der or on a torus is a set of N component curvilinear grids consisting of continuous transformations d, from the unit square onto sub-regions D; chosen such that the union of the D; contains D and the union of the boundaries of the D, contains the boundary of D, DC Jp, and ODC Jan, where k k D, = {di (r,s):0 <r <1,0<s<1}, and and uniform rectangular grids G, on the unit square, a-1 g-1 ; . G,; = ; ta d,...,m, 7 =1,...,nE >, m;,—-1 n,-1 for1<k<N. The component grids of a composite grid are curvilinear grids as defined in the Chapter 1. All of the remarks there about branch cuts due to the transformations and branch cuts on the cylinder or torus apply to these grids. Any region D may be broken up into enough subregions so that each may have a simple geometry. For example, a multiply-connected region D in the plane whose boundary consists of n smooth closed curves may be broken up into a rectangular grid on a rectangle that encloses D and nm annular grids each following one of the boundary curves. Those parts of the recangular grid which extend beyond the boundary of D can be ignored. The transformation used for the rectangular grid introduces no branch cuts. It can be a linear function of (r,s). The transformations used for the annular ~19- grids can be chosen to map one side of the unit square onto the boundary curve, the opposite side of the unit square onto a nearby closed curve in the interior of D and the other two sides of the unit square onto a branch cut. The Jacobian of each of these transformations extends to a function periodic in the tangential direction and the transformation can be chosen so that its Jacobian is nonsingular and continuous, even across the branch cut. 2.2 PDEs on Composite Grids Consider a PDE boundary value problem F (50) =0 for (c,y)<€D (2.2.1) where F represents the PDE in the interior of D and the boundary conditions on the boundary of D, and it is understood that higher derivatives may be involved. For each k we can define a function u;, which represents u in cartesian coordinates on the part of the unit square which is mapped by d; onto D, uz(r,s) = u(dy(r,s)) for dy(r,s) € D. (2.2.2) If the Jacobian of the transformation d; is non-singular we can use it to express the problem F in the (r,s) coordinates of each component grid. Let _ Ody (2, a, Ji (r,s) = ars) [ v. ; (2.2.3) Then the inverse of J;(r,s) is Joi(r,8) = 1 ( Ys ee | = [ % q (2.2.4) LrYs — UrLs ~Ur Ly v With the understanding that by r,, r,, 8, and s, we mean the functions of (r,s) which form the inverse of the Jacobian, we can express the boundary value problem in (r,s) coordinates as fed fe) OQ O F (5. + S29 Ty a + syst =0 for d,(r,s) € D. (2.2.5) ~20- The transformation d,; must satisfy the smoothness constraints (1.2.6). Parts of the boundary of the unit square may not correspond, under the trans- formation d;, to the boundary of D. A periodicity condition applies on those parts of the boundary of the unit square which correspond to branch cuts. To specify the problem completely we will in genera] need to impose matching conditions on the remaining parts of the boundary of the unit square. Actually since some parts of the unit square may correspond to the exterior of D and other parts of the unit square may correspond to parts of D that are represented on other grids, the PDE and boundary conditions (2.2.5) may apply only on part of the unit square, and the matching conditions need only be applied on the boundary of this part. We will in general need to apply matching conditions in the regions of overlap between the D,. In order to do this in a systematic way it helps to be able to order the transformations d,; so that wherever two or more of the D; overlap at the boundary of D, this occurs at the boundary of the D, with the highest index k of those involved. Although it is not always possible to order a given set of transformations in this way, it is always possible to choose an appropriate set of transformations that can be so ordered. To help in deciding on an appropriate sequence of transformations, imagine drawing D, in the plane, next drawing Ds, erasing that part of D, which overlaps D2, and so on, finally drawing Dy and erasing those parts of each underlying D, which it overlaps. The result will be that the entire region Dy remains, along with those parts of each underlying D; which are not covered by some D; with 7 > 1. If the sequence was chosen properly, the boundary of D is composed only of visible parts of the boundaries of the D,. Figure 2.1 illustrates how an appropriate sequence D,, Dz, Dz of subregions may be chosen for a periodic channel D. First D, is chosen to be a rectangle enclosing D. Then D, is chosen so that it follows the lower boundary of D and overlaps D,. Finally D, is chosen so that it follows the upper boundary of D and overlaps both -21- D, and Dz, with care that it does not overlap the lower boundary of D. If D2 and D, had been chosen in the reverse order, part of the upper boundary of D would have been hidden and the sequence of subregions would have been inappropriate. ™ rm) ; 4 Lo eee nee eee eee eee Ss ne The Region D Subregion D, Unit Square for G, Pp mn Pp 8 Subregions 0,0, Dy Subregion D; Unit Squore for G3 LO . Relevant Parts of Subregion D; Unit Square D,, Oz, Ds for G; B: Boundary M: Matching P: Periodicity Figure 2.1: Overlapping Subregions of a Periodic Channel Now we are ready to specify the neccessary matching conditions. First we specify conditions identifying uy with u,; on regions of overlap between Dy and the underlying regions D;. Where Dy overlaps other regions we apply the matching condition uy(r,s) = u,(d5*(dy(r, s))), where 7 is the highest index of those regions which underlie Dy at dy(r,s). The matching condition is applied only on those edges of Dy where neither boundary conditions nor periodicity conditions apply. 22- Depending on the type of PDE involved, some higher derivatives of uy and u,; may also need to be matched in the same way. In this way Dy is divided into two parts: the part (A) where the PDE and boundary conditions are applied and the part (B) of its boundary where the matching condition is applied. Matching conditions for the underlying regions are similarly specified. Where region D, overlaps other regions we apply the matching condition u,(r,s) = uj(d5"(di(r, s))), where 7 is the highest index of those regions which overlap or underlie D; at d,(r,s). The matching condition is applied on a curve (or curves) chosen so as to divide D; into three parts: the part (A) where the PDE and boundary condition are applied, the curve or curves (B) where the matching condition is applied, and a part (C) which includes any of D; that is exterior to D. Where 2 > 7 this curve will lie along edges of D; on which neither boundary conditions nor periodicity conditions apply; where z < 7 the curve lies on that part of D; where the PDE and boundary condition would otherwise be applied. By induction we can demonstrate that it is possible to cover D with overlapping parts of all the D; where the PDE is applied, each of which is separated from parts exterior to D by curves where a matching condition is applied. At this point the problem F is completely specified in (r,s) coordinates. Figure 2.1 shows how the curves may be chosen where matching conditions are applied on the three regions D,, D., Ds of the periodic channel. On D, the matching condition ui(r,s) = u2(dy’(di(r,s))) is applied on the lower halves of curves that separate parts of D, at the left and right, where the PDE applies, from a part which extends beyond the boundary or D; the matching condition ui(r,s) = us(dz'(di(r,s))) is applied on the upper halves of these curves. On Dy the matching condition u2(r, s) = uy(dj*(do(r, s))) is applied on a curve which cuts off part of the top of D2; the matching condition ug(r,s) = us(dz'(da(r,s))) is applied on the remainder of the top edge of D,. On D3; the matching condition ug(r,s) = u,(dyj*(da(r,s))) is applied on parts of the bottom edge of Dz at the -23- left and right; the matching condition ug(r,s) = u2(dz'(dg(r,s))) is applied on the remainder of the bottom edge. Each region D, is divided into parts of three types: parts (A) where the PDE and boundary condition apply, curves (B) where matching conditions apply, and parts (C) which are either redundant (covered by some other D; with jy > k) or exterior to D. The PDE and boundary conditions (2.2.5) are discretised at all of the gridpoints of G; which lie in the part of unit square corresponding to the parts (A) of D,. For example, the difference approximations (1.2.7) may be appropriate. The discretisation at some of the gridpoints near the curves (B) may involve other gridpoints not in parts (A). The value of u at these other gridpoints is given by the matching condition, approximated discretely using bivariate interpolation for- mulas. The set of gridpoints where the discrete matching condition is applied are the discrete analogue of the curves (B) of D,; gridpoints (A), where the PDE and boundary conditions are discretised, are separated by this set (B) from redundant and/or exterior points (C), where uj, is not defined. The matching condition uz (Ty8) = up (dz (di (r,))) (2.2.6) mMr-1? nme-l may be approximated at the point (r;,s;) = ( tol it.) of grid G,; using the interpolation formula uz (Ti, 83) = ‘> Ope (7, 95047 Ue (Ty, Sy.) + OC Ar’, As") (2.2.7) (49) where Ar’ = moi? As! = west (ri, sh) = (22.254) and ay, is a matrix of coefficients chosen to provide interpolation from grid Gy, onto grid Gy to the desired order of accuracy. For example, bilinear interpolation of uz(r,,s;) from grid G, takes a linear combination of u,: at the four gridpoints of G, at the corners of a rectangle with sides of length Ar! and As’ centered nearest to the point (r’, s') = d;, (dx (ri, ;)). The weights in the linear combination are the nonzero coefficients of -24- Ore (2, J,7,7') given in Table 2.1, where 6r! = (r'—r},)/Ar', 6s' = (s'—s).)/As', and (c4,,95) are chosen such that 24-1 < r'/Ar' <% and jj-1<s'/As' < 75. This gives second-order accuracy (p = g = 2) in equation (2.2.7). Biquadratic interpolation of uxz(7;,8;) from grid G, takes a linear combination of u, at the nine gridpoints of G, within a rectangle with sides of length 2Ar' and 2As' centered nearest to the point (r',s'), with the weights given in Table 2.2, where (25, 9) are chosen such that th — 1.5 <r'/Ar' < 14 —0.5 and 7j — 1.5 < s'/As' < 7, — 0.5. This gives third-order accuracy (p = q = 3) in equation (2.2.7). (a', J’) Xk ke («, J a’, J') (205 Jo) (1 — 6r')(1 — 6s’) (th + 1,36) 6r'(1 — 6s') (75,96 + 1) (1 — ér')és' (t6 + 1,96 + 1) ér'és! (7",7') are (i, 950", 9') (% — 1,36 ~ 1) 5r'(1 — 6r')6s!(1 — 6s!) /4 (7949 — 1) —(1+ ér')(1 — 6r')ds'(1 — 6s') /2 (% + 1,93 — 1) ~(1 + 6r')6r'Ss!(1 — 68')/4 (%) ~ 1,30) —6r'(1 — 6r')(1 + 68')(1 — 6s') /2 (2, Jo) (1 + 6r')(1 — 6r') (1 + 6s'}(1 ~ 6s’) (15 + 1,95) (1+ 6r')ér'(1 + 6s')(1 — 58')/2 (% — 1,95 + 1) ~6r'(1 — 6r')(1 + 68')5s'/4 (%, 96 + 1) (1 + 6r')(1 — 6r')(1 + 6s") 6s! /2 (% + 1,93 +1) (1 + br')6r'(1 + bs')bs'/4 Table 2.2: Biquadratic Interpolation Weights 2.38 Composite Grid Construction The method used here for the construction of composite grids is a generalisa- tion of the method used by B. Kreiss [14] to construct two-component composite -25- grids for simply-connected regions. Atta and Vadyak [1] also have constructed two- component composite grids, and show that composite grids are useful for three- dimensional problems. Benek et al. [2] have constructed two-component composite grids. A composite grid for a region D is a set of curvilinear grids for subregions D, whose union contains D and the union of whose boundaries contains the boundary of D. To construct a composite grid for D we (1) decide on the sequence of overlapping subregions, (2) construct curves forming the boundary of D and of the subregions, including any branch cuts, (3) use these curves to construct a curvilinear grid for each D,, and (4) decide which points on each grid are used for discretisation of the PDE and boundary conditions, which points are interpolated from other grids, and which points are not used at all. The first step is somewhat of an art, and would be difficult to automate. The second and third steps are straightforward, as they require only choosing the parametrisations c; of the curves and the transformations d, from a small class of available functions and choosing the number of gridlines in each direction. The fourth step would be tedious except that it has been automated. The details of these four steps are given below. The first step in construction of a composite grid for a region D is to decide how D may be broken up into a sequence of N overlapping sub-regions D, whose union contains D and the union of whose boundaries contains the boundary of D, so that the shape of each D, is simple enough that a grid can be constructed for each. A single grid will suffice if D is simple enough. For example, it is easy to construct a grid for a convex region with four corners (something like a rectangle). It is also easy to construct a grid for an annular region each of whose two boundary curves deviates only slightly from a circle. However, in general a composite grid consisting of two or more grids is appropriate and is much easier to construct than a single grid for the same region. or PDE boundary-value problems with boundary layers, ~26- interior layers or other regions where the solution varies rapidly, an appropriate composite grid consists of fine grids in the regions of rapid variation and coarse grids elsewhere. The choice of sub-regions will depend on the application. However, the method outlined in §2.1 applies to a rather general class of problems. In order that matching conditions may be applied, the sequence of subregions must be chosen in accordance with the principles given in §2.2. The second step is to construct curves forming the boundary of D and of the subregions. We parametrise each curve C; by a smooth function ¢; such that C; = {c;(€):0 < € < 1}. These functions may be specified analytically if a formula for the boundary curve shape is known, or by a spline function if only some points on the boundary are known. In order that each transformation d, may be smooth and its Jacobian may be nonsingular, we require that the boundary of D, and any branch cuts consist of smooth curves intersecting at four corners with included angles smaller than 7, and choose d; so that the corners of the unit square are mapped to these points. Each side of the unit square is mapped to part of a curve which falls into one of three categories: (1), those curves C? which form the boundary of D; (2), those curves C? which have been chosen as branch cuts either of D (in case D lies on a cylinder or torus; i.e., is periodic) or of a subregion (in case the subregion is annular); and (3), any other curves C/ in the interior of D where two or more subregions overlap. A boundary curve C? or branch cut C? of D may form parts of the boundary of more than one of the subregions. In such a case, since the subregions overlap to cover D, those parts of the boundaries of the subregions overlap to cover the curve. Other parts of the boundary of D, consist of curves which typically do not form any part of the boundary of another subregion. Branch cuts of D appear in the plane as pairs of identical curves, shifted with respect to each other by the corresponding period vector. It is important in the discretisation of boundary-value problems that the transformations d;, from the ~27- unit square constructed using the parametrisations of these curve should have the same periodicity and be smooth across the branch cut. To facilitate this we require that the parametrisations differ by a constant vector, the period vector. The third step is to construct each component curvilinear grid. We specify the rectangular grid G, on the unit square by choosing the numbers n, and n, of gridlines in each direction. We construct the transformation d; from the unit square onto D, using the parametrisations of the curves, segments of which form the boundary of D,. Each side of the unit square is mapped onto a segment of one of these curves, and each corner of the unit square is mapped onto an intersection of two of these curves. In order that each corner may be well-defined, we require that these intersections not be tangential. Any of the methods discussed in Chapter 1 may be used to construct the transformation. However, since the principal appli- cation for composite grids is the discretisation of PDE boundary-value problems, it is usually important to use only those methods which result in a transformation which is smooth, even across branch cuts. Since enough subregions may always be chosen so that the shape of each is very simple, even the simplest methods are usu- ally adequate to construct the transformation d,. In particular, the “Coons Patch” transformation is guaranteed to be as smooth as the boundary curve parametrisa- tions. To compute it requires only four evaluations of the curve parametrisations, so it is very fast. It is easily modified to allow stretched coordinates for boundary layers and interior layers. For these reasons it is the method of choice in most cases. The fourth step is to decide which points on each grid are used for discretisa- tion of the PDE and boundary conditions, which points are interpolated from other grids, and which points are not used at all. Points where the PDE or boundary condition may be discretised will be called discretisation points, points which are interpolated from other grids will be called interpolation points, and points which are not used at all will be called exterior points. For simplicity it is assumed that the -28- PDE and the boundary conditions are both discretised to the same order and that the order is positive and even. If the discretisation is of order p, a square of points (on the rectangular grid) with p+ 1 points on a side is used in the discretisation. The discretisation of the PDE is assumed to be centered, while the discretisation of the boundary condition is assumed to be one-sided. A point may be an discreti- sation point if all the neighbouring points required for discretisation of the PDE or boundary condition are either discretisation points or interpolation points. Inter- polation is assumed to be of order p, where p > 0. Interpolation of a point on grid k from grid k' to order p using the formula (2.2.7) requires a square of points on grid k' with p points on a side. The point on grid k may be an interpolation point if all the points in the nearest such square on grid &’ are discretisation or interpolation points. An algorithm for finding which points may be discretisation points and which points may be interpolation points, based on these criteria, is given in Chap- ter 3. This algorithm attempts to minimize the total number of discretisation and interpolation points. When the grids overlap sufficiently and the grids are ordered in accordance with the principles given in §2.2, a PDE boundary-value problem of the given order can be discretised using the discretisation and interpolation points of the resulting composite grid. 2.4 Composite Grids for Multigrid The multigrid method (see, for example, Brandt [5]) is a fast iterative method for the discrete solution of linear elliptic PDEs. It is especially useful for problems in which the large number points required for an accurate discretisation precludes the use of a direct method to solve the resulting linear system of equations. The multigrid method is usually discussed in the context of a rectangular grid. When an elliptic problem is discretised on a composite grid, it is transformed into a set of elliptic problems on rectangular grids in the (r,s) coordinates of each compo- nent grid, coupled by interpolation equations. Using this fact, Henshaw [{12], [13] (cf. Appendix 2) has implemented a multigrid algorithm for composite grids by in- corporating the interpolation equations into the multigrid iteration. This algorithm uses a finite sequence of related composite grids M’, where M! is the grid on which the elliptic problem is to be discretised, and M! for | = 2,3,... are successively coarser composite grids. Corresponding component grids at each multigrid level are constructed using the same transformation d,. The rectangular grid Gj), (on ihe unit square) of component grid k at multigrid level / contains twice the number of gridlines in each direction as the corresponding rectangular grid Gin at the next multigrid level. To be precise, the numbers of gridlines in each direction on each grid at each level are related by mi, — 1 = 2(mi*! —1) and ni, —1 = 2(nit! — 1). The multigrid algorithm requires the transfer of data from level | to levels | — 1 and 1+1 at each point where the equation or boundary condition is applied. This requirement must be taken into account when deciding which points of each grid at each level are to be used for discretisation of the PDE and boundary conditions, which points are interpolated from other grids, and which points are not used at all. Because of this, the composite grids at all levels must be constructed simulta- neously. The one-level algorithm given in Chapter 3 for finding which points may be discretisation points and which may be interpolated generalises to any number of multigrid levels. -30~ CHAPTER 3 Composite Grid Software We have developed software for the interactive construction of composite grids for compact regions in the plane and for regions with cylindrical or toroidal symme- try. It may be used to specify composite grids consisting of an arbitrary number of grids. A manual for this software is included as Appendix 1. This software has been generalised to construct sequences of composite grids for use with multigrid meth- ods. In addition we have developed general-purpose subroutines for using composite grids in discretising PDE boundary-value problems. These include subroutines for reading a data file containing a description of a composite grid, for calculating the coefficients in the centered second-order discretisation of a linear combination of derivatives on a nine-point stencil, for interpolating data between the grids of a composite grid in the regions of their overlap, and for the second-order discretisation and direct solution of a second-order linear elliptic equation with mixed boundary conditions on a composite grid. These subroutines have been generalised to act on raultigrid sets of composite grids, and Henshaw |13] (cf. Appendix 2) has generalised the elliptic solver so that it may use the multigrid algorithm. 3.1 A Program for Composite Grid Construction We have developed the interactive Fortran program CMPGRD for construction of a composite grid, consisting of as many component grids as neccessary, for any compact region in the plane bounded by piecewise-smooth curves. The region may have cylindrical or toroidal symmetry, i.e., be periodic in one or two directions. -31- The region may also be multiply-connected; the boundary of the region may consist of any number of distinct closed or periodic curves. CMPGRD may be used to construct the curves forming the boundary of the region, to construct the component grids for subregions bounded by these and other curves, and to decide which points on each grid are used for the discretisation of a PDE and its boundary conditions, which points are interpolated from other grids, and which points are not used at all. The algorithms for each of these steps are given below. 3.1.1 Construction of Curves There are two methods by which CMPGRD may be used to construct the curves C; that bound a region D and the subregions D, which overlap to cover D. They may be specified either analytically by a Fortran subroutine supplied by the user or as spline curves through data points given by the user. Segments of several curves constructed by either method may be used to specify the four parts of the boundary of any subregion D; corresponding under the transformation d, to the four sides of the unit square, and curves specified by different methods may be used for the same subregion. (CMPGRD is set up so that other methods of curve construction may be incorporated later with a minimum of difficulty.) Instructions for using CMPGRD to specify curves are given in §A1.4.2; the computations involved are explained below. To specify a curve C; analytically, the user must provide subroutine CURVE, which, given the curve index 7 and the arcwise parameter € in the range |0, 1], returns the value (x,y) and the derivative (xg, y¢) of a function ¢;(€) which parametrises C;, and an integer code p; that indicates how the curve is to be used. General specifications for subroutine CURVE are given in §A1.5.1. For example, if the boundary of the region D consists of a single closed curve Cy, this curve could be specified analytically by a periodic function c, on the interval [0,1]; CURVE -32- should return the value c;(¢) and the derivative of this function, and the integer code p, = —1 to indicate that this curve is a boundary curve. One of the subregions D, used in constructing a composite grid for D might be an annular region that follows this boundary curve. In that case we would want to find a transformation d, that maps one side of the unit square onto the curve C,, the two adjacent sides onto a branch cut, and the opposite side onto a curve in the interior of D where D, overlaps other subregions. The curve forming the branch cut might be labelled C, and the other curve C3, and these might be specified analytically by functions Cc» and cs. The branch cut should intersect C, at the point c,;{0) = c;(1) and intersect Cz at c3(0) = c3(1), so that each side of the unit square may be mapped under the transformation d, onto a part of a curve parametrised on a subinterval of [0,1]. In general, every periodic curve should be parametrised on the interval [0,1] so that it intersects a branch cut at each end of the interval. The integer code pp = 2 is used to indicate that curve Cz is a branch cut, and that whenever one side of the unit square is mapped by any transformation d, onto a segment of this curve, the opposite side of the unit square will be mapped onto the same segment of the same curve. In general, branch cuts occur as pairs of identical curves, so that if C; is a branch cut then there must be a corresponding curve C,, identical to C,; but possibly shifted with respect to it by a constant vector. This correspondence would be indicated by p; = 7 and p,; = 1, of which p2 = 2 is a trivial case. The integer code p3 = 0 is used to indicate that the closed curve Cz in the interior of D is neither a boundary curve nor a branch cut. CMPGRD may be used to specify a curve as a cubic spline through data points given by the user either in a data file or through the use of interactive graph- ics. Curves specified in this way are numbered consecutively in the order they are specified, following the highest-numbered curve that was specified analytically by subroutine CURVE. At the time a curve is to be specified, it must be declared as -33- either a boundary curve, a branch cut or an interpolation curve, and it must be declared whether the curve is to be periodic or not. Each “knot” through which the spline is drawn may be specified either directly by giving its (x,y) coordinates or indirectly by giving the coordinates of a point nearby another curve and taking the knot to be the nearest point on that curve. A nearby point on the nearest curve is determined by a global minimization over a large number (typically 100) of points on each of the previously-specified curves. The nearest point on that curve is determined using a bisection iteration. By this method it is possible to insure that the new curve intersects another particular curve, and that the intersection is at an endpoint of the new curve. In particular, if the new curve was declared to be periodic, the first knot must be a point on a branch cut, and the last knot will be taken to be the corresponding point on the opposite branch cut; if the two branch cuts coincide, the first and last knots are identical and the curve is a closed curve. At least two knots are required to specify a non-periodic curve; a cubic spline curve is constructed which passes through each of the knots in the order they were entered, with zero second derivatives at the end-points. At least four points are needed to specify a closed periodic curve, and at least two points are needed to specify a periodic curve with distinct endpoints; in either case, a cubic spline with periodic boundary conditions is constructed which passes through the knots in the order they were entered, after the last point has been corrected to be the periodic image of the first point. In all cases, the parainetrisation of the spline function is chosen so that at the knots, its parameter (independent variable) is proportional to the polygonal arc-length from the initial end-point through the knots. Since only one new curve may be specified at a time, new pairs of distinct branch cuts cannot be specified; any new branch cut is paired with itself. ~34- 3.1.2 Construction of Grids A grid for each subregion D, of D is constructed in several steps. First, the boundary of D, is specified as four segments of intersecting curves. Second, a mapping from the unit square onto D,; is chosen. Third, a stretching function is specified which maps the unit square onto itself, to allow some freedom to resolve boundary and other layers. The transformation d; is obtained by first applying the streching, if any, and then the mapping. Fourth, the number of gridlines in each direction is specified. These steps are explained below. Instructions for using CMPGRD to perform these steps are given in §A1.4.3. The boundary of the subregion D, is specified as four segments of intersecting curves; each side of the unit square is mapped under the transformation d, onto one of these segments. In order that these intersections be well-defined, only non- tangential intersections may be used to determine them. CMPGRD may be used to specify these segments in terms of data points given by the user either in a data file or through the use of interactive graphics, as follows. First, a curve C;, is chosen by giving the coordinates of a point nearby that curve. The nearest curve is determined by a global minimization over a large number (typically 100) of points on each of the curves. Then each endpoint of the segment of C; is chosen by giving the coordinates of a point nearby the intersection of C; with another curve C, or nearby an endpoint of C;. A nearby point on the nearest curve is determined by global minimization over a large number of points on each of the curves. If the given point is nearby the point c;({ ) on C;, the endpoint of the segment is taken to be the endpoint of C; nearest to c;(é ) in terms of €. Otherwise, if the given point is nearby the point c,;(£;) on C; (with 7 #2), the endpoint of the segment is taken to be the nearest intersection of C; and C,;. This intersection is determined by first finding a point c;(&) on C; nearby c;(€,) using a global minimization over points on C;, and then iterating simultaneously on € and &, by Newton’s method ~35— until the intersection c;(&;) = ¢;(€:) is accurately found. When four curve segments have been specified that together form the boundary of a subregion D, of D, we can use these to construct a mapping from the unit square onto D,. To construct the mapping it is only necessary to associate each side of the unit square with a curve segment and to choose the type of mapping to be constructed from this boundary data. First, a side of the unit square is cho- sen. Second, a curve segment is chosen by giving the coordinates of a point nearby that segment. The nearest segment is determined by a global minimization over a large number (typically 100) of points on each of the segments. Third, the orien- tation of the segment is chosen so that each endpoint of the segment corresponds to the appropriate corner of the unit square. CMPGRD checks that the segments associated with each side of the unit square are compatible with each other. In par- ticular, it checks that each corner of the unit square corresponds to the intersection of two curve segments at their endpoints, and that if one side of the unit square corresponds to a segment of a branch cut then the opposite side of the unit square corresponds to the identical segment of the opposite branch cut. When all sides of the unit square are associated with compatible curve segments, parametrisations of the curve segments may be used as boundary data to construct the the mapping from the unit square onto D,. Although many types of mappings could be con- structed, CMPGRD offers only one, a “Coons Patch” [7], also known as transfinte interpolation. . This mapping is a linear combination of points on four segments of curves which form the boundary D,. To avoid problems with notation, suppose we have numbered the curves so that segments of curves C, and C2 form two opposite sides of the region D,, and segments of Cz and Cy, form the other two sides, so that these curves intersect at 1 (£13) = €3(&31), C1 (E14) = C4(E41), (3.1.1) C2( £23) = €3(&39), €2( E24) = C4 (E42). -36- CMPGRD uses the Coons patch defined by (x,y) = (1—u)ei(Eis + (E14 ~ E13) t) + wee (E23 + (E24 — €23)t) + (1 —t)e3(€31 + (32 ~ €31)u) + tea(Ear + (E42 — 41) 4) — (1 —u)(1 — the: (13) — u(1 — t)e2(23) — (1 - u)ter (E14) — utes (Eo) (3.1.2) CMPGRD is set up so that other methods of grid construction may be incorporated later with a minimum of difficulty. The only restrictions are that the mapping may be evaluated anywhere on the unit square (not just at gridpoints), that it be smooth, and that its Jacobian be invertible wherever D; overlaps another subregion of D. If the composite grid is to be used for discretising a PDE boundary-value problem, the Jacobian must be invertible everywhere in the intersection of D, and D. CMPGRD checks the Jacobian determinant of the mapping for zeros and changes of sign by evaluating it on a 40 point by 40 point grid on the unit square; if no zeros or changes of sign are found at these points, it searches for zeros or changes of sign nearby the point where the smallest value of the determinant was found. In order to allow some freedom to construct grids that resolve boundary layers and interior layers, the transformation d; from (r,s) coordinates on the unit square to (x,y) coordinates is constructed as the composition of two functions, namely, a “stretching” function from (r,s) coordinates on the unit square to (¢,u) coordinates on the unit square (which maps the unit square onto itself), and the mapping func- tion from (t,u) coordinates to (z,y) coordinates as constructed above. Although many different stretching functions could be useful, CMPGRD provides only inde- pendent exponential stretching functions in the r and s directions: Each of these is a function of the form H(t) = Lt Shales) — 210) 14+ [a (1) — ;(0)] (3.1.3) where o;(t) = a; tanh b(t — c;) (3.1.4) -37- if the problem is not periodic in r, and o;(t) = a; tanh b(t — ¢;) + S> a;{tanh b(t —¢; ~j) + tanhb,(t—¢;+ 7)] (3.1.5) g=l if the problem is periodic in r. Using trigonometric identities, the infinite series (3.1.5) may be written as = a; sinh 2b;(t — ¢;) s(t) = a; tanh 6,( —, 3.1.6 oi(t) = ai tan e+) Soa — ¢;) + cosh 27); ( ) j=l which converges quickly and requires only three evaluations of hyperbolic functions (not including the constants cosh 276,;). Each term o;(t) has three free parameters: the weight a;, the exponent 6; and the location ¢;. Each term produces a stretching in the ¢ coordinate centered at t = c,. That is, the image in the (t,u) plane of a uniform grid in the (r,s) plane is a nonuniform rectangular grid with a maximum in density of gridlines for t near each c;. The density of gridlines in the ¢-direction at the point (t,u) is proportional to the derivative r'(t), and may be adjusted by varying the parameters a;, 0; and c; in each term. In particular, a; is approximately the ratio of the number of gridlines in the stretched region of the grid near c; to the total number of gridlines in unstretched regions, while 6; is approximately the ratio of the density of gridlines near c; to the density of gridlines in unstretched regions. If the problem is periodic in r, an infinite serics of correction terms is needed to ensure that the derivative ¢'(r) of the stretching is periodic. For example, if there is one term in the stretching (n = 1), and we choose ay = 1, 6} = 1 and c; = 0, the t coordinate is stretched at ¢ = 0; if the problem is periodic in r, the t coordinate is equally stretched at ¢ = 1. Another sum of the same form expresses s in terms of u, and together these specify (t,u) = (¢(r),u(s)), where t(r) and u(s) are the inverses of functions given by weighted sums of linear and hyperbolic tangent terms. 3.1.3 Construction of Composite Grids Once the sequences of regions D,, transformations d;, and grids G, have been -38- specified, CMPGRD may be used to decide which points on each grid are used for the discretisation of a PDE and its boundary conditions, which points are inter- polated from other grids and which points are not used at all. The algorithm for determining these points is given below. In addition, the (r,s) coordinates (in the unit square of the grid from which it is interpolated) of each interpolated point are computed, so that the coefficients a, in the interpolation formula (2.2.7) may be calculated. This information, together with the (z,y) coordinates and Jaco- bian derivative of each transformation d, at the gridpoints of G,, is sufficient to discretise a PDE boundary-value problem on the composite grid. Each point on each grid of the composite grid is used either for the discretisa- tion of a PDE or its boundary conditions, or as a point where function values there are interpolated from another grid, or it is not used at all. A gridpoint which may be used for the discretisation of a PDE or its boundary conditions will be called a discretisation point, and a gridpoint where function values are interpolated from another grid will be called an interpolation point. A point which is neither a discreti- sation point nor an interpolation point will be called an exterior point. A gridpoint may be a discretisation point if each adjacent point used in the discretisation at that point is either a discretisation point or an interpolation point. The points used in the discretisation depend on the order of the discretisation and whether the PDE or the boundary condition is discretised. It is assumed that centered differences are used for the PDE, so that, if second-order differences are used to discretise a second-order differential operator then a nine-point stencil will be required, and if fourth-order differences are used then a 25-point stencil will be required. It is assumed that one-sided differences are used for the boundary conditions, using a stencil of the same size as for the PDE. A gridpoint of G,; may be an interpolated point, interpolated from G,., if each point in Gy from which function values at that point are to be interpolated is either a discretisation point or an interpolation point. ~39— It is assumed that bilinear or biquadratic interpolation is used, so that a square of four or nine points, respectively, of G, centered nearest to the (r,s) coordinates (in G,-) of the interpolated point is required for the interpolation. The first step is to mark each gridpoint on each grid as a discretisation point, so that, unless it is later determined that the point cannot be a discretisation point, it may be used in the discretisation formulae of adjacent discretisation points on the same grid and in the interpolation formulae of interpolation points on other grids. The second step is to mark non-boundary points that lie close to a boundary side of another grid as exterior points. (Boundary points and boundary sides are those gridpoints and sides of the unit square which are mapped onto points and segments, respectively, of curves which were declared to be boundary curves.) This is done so that if a grid covers a hole in a multiply-connected region, then those of its gridpoints in the hole will be marked as exterior, as they should be. The third step is to find which grid, if any, each point on each grid can be interpolated from. Starting with the last (highest) grid and working down to the first (lowest) grid, each gridpoint on each G,; is examined to find the highest G,: with k' > k, if any, from which it can be interpolated. For each side of Gx, each point which is neither a periodic point (a point on a side of the unit square which is mapped onto a branch cut) nor a boundary point, and cannot be interpolated from a higher grid, is then examined to find the highest G, with k' < k, if any, from which it can be interpolated. The points used in these latter interpolations are marked as needed for interpolation, so that they may not be later marked exterior. The end result is that each point that can be interpolated from a higher grid is marked as interpolated from the highest such grid, each point on each non- boundary, non-branch-cut side of each grid which can only be interpolated from a lower grid is marked as interpolated from the highest such grid, and each point used in the interpolation of a point on a higher grid is marked as needed for interpolation. ~40- The fourth step is to examine each remaining discretisation point to see if it can really be a discretisation point. Each grid is swept to find any point which is currently marked as a discretisation point but some of whose neighbours, required for discretisation, are neither discretisation points nor interpolated points. Such points are then marked as exterior points. If any of these points has been marked as needed for interpolation, the composite grid will be useless and CMPGRD gives an error message to that effect. Each grid is swept repeatedly until no more points are found which must be changed to exterior points. The end result is that each point marked as a discretisation point really can be a discretisation point, but that some of the points marked as interpolated points can no longer be interpolated. However, provided there is sufficient overlap between the grids, all the interpolation points needed for discretisation at the remaining discretisation points can still be interpolated. The fifth step is to examine each interpolated point to see if it is needed for the discretisation of any discretisation point or has been marked as needed for interpola- tion. Starting with the lowest grid and proceeding to the highest, each interpolated point on each grid is examined. Those interpolated points which are not needed are marked as exterior points. Those points on higher grids required to interpolate the remaining interpolation points are marked as needed for interpolation. The end result is that all of the remaining interpolation points can really be interpolated, provided there is sufficient overlap between the grids, and that each point needed for interpolation of some point on another grid is marked as such. At this point all the points on all the grids are marked truthfully as to whether they may be discretisation points, interpolation points or neither. However, when there is a choice, it may be preferable that a point be a discretisation point rather than an interpolation point, since the system of interpolation equations may be better-conditioned if there are fewer interpolation points, and hence fewer equations. _41- The sixth step is to examine each interpolation point to see whether it could just as well be a discretisation point, and if so, then mark it as such. The end result is a composite grid with fewer interpolation points. ‘The seventh and final step is to output to a file all the information about the composite grid relevant to discretising a PDE boundary-value problem, namely: the coordinate bounds of the region D, defining a bounding rectangle to help with plotting; the number of grids; the number of gridlines in each direction on each grid; a flag array indicating for each point on each grid whether the point is a discretisation point, an interpolation point or neither; the boundary conditions on each side of each grid, indicating also which curve the side is mapped to; a list of interpolated points on each grid, indicating the grid from which it is interpolated and the (r,s) coordinates of the point on that grid; the (x,y) coordinates and the Jacobian derivative at each point on each grid. This file may be read, with the help of a utility subroutine described in the next section, by a program written to discretise a PDE boundary-value problem. The above algorithm spends most of its time in computing the (r,s) coordinates on one G, of points on another G;,. CMPGRD stores the coordinates (2,4) = di(rs,s;) = as ( t-1 j-l m,—-1'n, -1 of all the points on all grids, but must calculate the inverses of these transformations (r',s') = dy (z,y) to determine, for example, whether the gridpoint (r;,s,;) can be interpolated from G,,. This function inversion is performed in several steps. First, it is determined whether (x,y) lies within a previously-calculated rectangle which bounds D,.. If not, then the inversion does not proceed. Second, beginning with an arbitrary gridpoint (7',7') on Gy, a sequence of adjacent gridpoints of Gy is followed which finds a local minimum of |z ~ 2'| + |y — y'|. Third, if this local minimum occurs on -42- the boundary of G,,, then |z— x'|+ |y ~ y'| is minimized globally over the boundary of G, and the second step is repeated. At this point we have found the nearest point (7’,7') on Gy: in the one norm. Fourth, Newton’s method is used to solve the equation d,-(r',s') = (z,y) starting from the coordinates of the point (1’,j'), taking periodicity into account when neccessary. If the Newton interations take (r’, s') outside the unit square (beyond a small tolerance), the inversion does not proceed. 3.2 Subroutines for Using Composite Grids We have developed several general-purpose subroutines for using composite grids in discretising PDE boundary-value problems. These include subroutine CG- DATA for reading a data file containing a description of a composite grid, subroutine DISCRT for calculating the coefficients in the centered second-order discretisation of a linear combination of derivatives on a nine-point stencil, subroutine INTERP for interpolating data between the grids of a composite grid in the regions of their overlap, and subroutine SOLVER for the second-order discretisation and direct so- lution of a second-order linear elliptic equation with mixed boundary conditions on a composite grid. The use of these subroutines is outlined below. CGDATA is a Fortran subroutine that reads a file containing the description of a composite grid and stores it in a form that can be used by other composite grid subroutines. DISCRT is a Fortran subroutine that, given a point (7,7) on G, and a set of coefficients Crz, Czy, Cyy, Cx, Cy and c,, computes the weights for the nine points in the centered second-order discretisation of the operator a a? a? a ) L = Cre 5-3 + Coy} Cw arg + Ce ys + vB 3.2.1 “rr ay2 T © "Axdy “wy Oy? re Ox rey oy ¥ ( ) It| on G, using the difference approximations (1.2.7) and (1.2.8). This subroutine is meant as a prototype to help in discretising PDE problems. -43- INTERP is a Fortran subroutine that, given a function defined on the discreti- sation points of a composite grid, computes interpolated values of the function at the interpolation points using bilinear or biquadratic interpolation formulae. If the interpolation formulae for any of the interpolation points involve interpolation points on other grids, then the interpolation formulae become a coupled linear sys- tem of equations. INTERP factors this sparse system using subroutines from the Yale Sparse Matrix Package (YSMP) [11]. SOLVER is a Fortran subroutine that computes the second-order discretisation and direct solution of a second-order linear elliptic equation ; O*u 4 O7u in O*u n Ou n Ou 4 f (3.2.2 rz ~ Cy c la aT ch = ote Ox? "Oxdy Oy? “on 1 Oy Cet ) with oblique mixed boundary condition Ou Ou b, — —+bu= 2. aq + bvGy + bet g (3.2.3) using the difference approximations (1.2.7), (1.2.8) and (1.2.9), and bilinear or bi- quadratic interpolation. SOLVER factors this sparse system using subroutines from the Yale Sparse Matrix Package (YSMP) |{11]. —44- CHAPTER 4 Applications of Composite Grids Composite grids are useful in discretising a wide variety of PDE boundary-value problems by finite-difference methods. The primary advantage of using a composite grid rather than a single grid is the flexibility it allows in discretising problems on regions of complicated geometry and problems with boundary layers and internal layers. The high degree of contro! over the density of gridlines in a composite grid permits accurate discretisation of problems using in some cases only a fraction of the number of gridpoints a single grid would require. We have developed general- purpose subroutines for using composite grids in discretising PDE boundary-value problems. The discretisations of an elliptic boundary value problem with an internal layer and of a supersonic flow problem with multiple shocks are given as examples of the uses of composite grids and the capabilities of the software. 4.1 An Elliptic Problem With an Internal Layer To demonstrate how effectively composite grids may be used in PDE boundary- value problems with layers of rapid transition, we consider here a second-order elliptic equation with Dirichlet boundary condition. The coefficients in the equation are chosen so that the solution will have an interior layer, and the forcing terms and boundary data are chosen so that the exact solution is known. Numerical results show that when the internal layer is very sharp, a composite grid can be used to calculate an accurate solution, while a single rectangular grid cannot; the number of gridpoints required to resolve the solution on a single grid is prohibitive. —~45— We consider an elliptic boundary-value problem of the form a(z,y)V*w —w = f(z,y) (4.1.1) inside the square D = [—1,1]? with Dirichlet boundary condition w(z,y) = g(z,y), where a(z,y) = O(1), f(x,y) = O(1) and g{z,y) = O(1), and a(z,y) is bounded away from zero independently of € except near the circle of radius one half, where a(z,y) = O(e”). In particular, we choose 1 2 a(z,y) = (orem) , where (P(x, y)) tanh 72=4 + tanh #2 pt tan + tan Ze and P(z,y) — /_? + y?, 1+ tanh + 2+ tanh = For positive values of ¢, Q(r) is a smooth, odd function of r, so a(z,y) is a smooth Q(p) = function of x and y. For small, positive e, a(z,y) has approximately the value 9 over all of the unit square, except near the circle of radius one half where it drops to approximately 9c”. We choose the forcing function and boundary data so that the exact solution is w = cos7Q(P(z,y)). In particular we choose _ i 75 msintQ and g(z,y) =cosmQ, where Q = Q(P) and P= P(z,y). f(z,y)=-( + 2”) cos TQ — For small, non-zero values of €, w(x,y) and f(x,y) are smooth everywhere except for O(1) transitions across a layer of width O(e) near the circle of radius one half, and g({z,y) is smooth everywhere on the boundary of D. In the limit as e — 0 w(x,y) tends to the function which is one inside the circle of radius one half and zero on the rest of D. In order to resolve the solution w(z,y) and the coefficient function f(z,y) ona uniform rectangular grid on D, the numbers n, and n, of gridlines in each direction must be taken large enough so that the gridline spacing is of the order of the layer width or smaller. That is, A 2 O Ay 2 -46- and the total number of gridpoints is Nj = nzny = O(e~”), For example, if we take € Agr = Ay= 5 ande= a then the total number of gridpoints is N, = 1.6 x 10°. This number of gridpoints would be needed in order to discretise the problem using second-order centered differences and obtain an accurate solution. We can construct a composite grid, consisting of a radially stretched (nonuni- form) annular grid covering part of D that includes the layer and a uniform square grid covering the rest of D, which resolves w(z,y) and f(z,y) using a number of gridpoints which is independent of the layer width e. We construct the composite grid shown in Figure 4.1 using the two transformations di(r,s) = (27 —1,2s —1) from the unit square [0,1]? onto Dy = D and d2(r,s) = (0.2 + 0.67 (r))(cos 278, sin 27s) from the unit square onto the annulus D; = {(z,y):0.2 < P(z,y) < 0.8}, where T(r) is the inverse of the function R(t) = t+ $(tanh 4 + tanh ze) 1 + tanh 3: The radial stretching function R(t) concentrates gridlines in a layer of width ¢ around the circle of radius one half. Although a variety of functions could have been used to achieve the same result, this particular R(t) was chosen because CMPGRD readily provides stretching functions of this form. In the new coordinates (r,s) given by the transformations d, and d, the so- lution has no layer of rapid transition; the problem may be discretised and an accurate solution obtained using a number of gridpoints which is independent of é. In the scaled rectangular coordinates given by the transformation d,, equation 4.1.1 becomes 1\°/@ @ (20°) (2+ 5) Mamma = dite s)}, (A-2) —47— Tl i T oot a T tan oe MeEVIN f aac = a ae - = “ ms y Z = 4 x ort AQ /; Z S Aye ae: 7 ee ones 0 [saad 7 ¥ hk \\A ihia di Tih t re | LO i aul sho ne'u ‘ pat rT i Wy meat : a6) ‘Otge i aul tean8t BeGU care CS a inary SOY ee - a 4g e, “S NS Sy ‘< yr Ppt roe as i i Computational Grid Computational Grid L L L L -1.00 -1.00 -1.00 -1.00 L L L L -1.00 -1.00 -1..00 -1.00 L L L -1.00 -1,00 -1,00 -1¢00 L L L L -1.00 -1,00 -1.00 -1.00 Exact solution Exact solution Calculated solution Calculated solution Fig. 4.1: Comparison of a Single Rectangular Grid and a Composite Grid for an Elliptic BVP with an Internal Layer _48- where wi(r,s) = w(d,(r,s)). The coefficients of both 0?w,/Or? and 0?w,/ds? in equation 4.1.2 are O(1) and bounded away from zero independently of ¢, and the coefficients, the forcing function and the boundary values are O(1) and smooth in the part of D where the square grid is used. Therefore, the number of gridpoints needed to accurately discretise the problem on this part of D is independent of «. In the stretched polar coordinates provided by the transformation d2, equation 4.1.1 becomes R' 2 Ow» R! 1 dw, 1 2 aw, . (4.1.3) where w2(r,s) = w(de(r,s)). The coefficient of 0?w2/dr? in equation 4.1.3 is O(1) and bounded away from zero independently of ¢, and the coefficients, the forcing function and the boundary values are O(1) and smooth for all (r,s) € {0,1]?. There- fore the number of gridlines in the r-direction needed to accurately discretise the problem on this part of D is independent of ¢. Since we know the solution is inde- pendent of s, the number of gridlines in the s-direction does not affect the accuracy. Hence, the number of gridpoints needed is independent of e. Table 4.1 gives a comparison of numerical results for this problem for € in the range of + to a on a single rectangular grid and on appropriately stretched com- posite grids. It shows that a fixed number of gridpoints may be used to accurately solve the problem with a composite grid, but not with a singular rectangular grid. € Single Grid Composite Grid ee) .0062 1/10 0332 0112 1/20 4190 0196 1/40 0733 Table 4.1: Maximum Error We can discretise the problem using a number of gridpoints independent of ¢€ ~49- in the above example only because the ratio R'/Q' is O(1) even where the solution has a layer of rapid transition. That is, we were able to choose a function R(t) to specify a transformation of coordinates so that (1) the gridlines are closely spaced only in the layer where Q' is small, and (2) the density of gridlines in the layer is high enough to accurately resolve the PDE coefficients and the solution on the grid. In more general problems with layers this is possible only when the location and width scale of the layers are accurately known. However, even for a problem with a layer whose location is known only approximately, a composite grid can be used to accurately discretise the problem with significantly fewer gridpoints than a single grid would require. 4.2 A Supersonic Flow Problem A typical feature of supersonic flow is the formation of shocks. In two dimensions these are curves across which some component of the flow is discontinuous. In general these curves may intersect each other and they may appear, move and disappear as the flow evolves. Wherever the flow is smooth, it may be modelled by a nonlinear hyperbolic system. However, such a system cannot model shocks unless, for example, it is regarded as the inviscid limit of a parabolic system, or some special rules are applied wherever shocks appear. We shall take the former approach, which has the advantage of being straightforward to apply and the disadvantage of approximating shocks as layers of rapid transition instead of discontinuities. Composite grids are well suited to the discretisation of steady supersonic flow problems. By allowing the use of fine grids in the regions of shocks and relatively coarse grids where the solution is smooth they significantly reduce the total number of gridpoints needed to obtain an accurate solution. The position of the shocks cannot in general be determined a priori, so a composite grid which will resolve these shocks cannot initially be found. However, an iterative procedure may be used _50- to construct such a grid and use it to calculate an accurate steady solution. Since automatic refinement of composite grids is not currently available, this iteration requires considerable human intervention. This approach is practical because only a few iterations are needed in order to construct a composite grid on which the solution is resolved. This approach is impractical for time-dependent solutions with moving shocks since a new composite grid may need to be constructed after every few timesteps. NSASNASNANAAAAN SN SAN ANS SSSA NAAN ANA ANAS ANAS NS Outflow ——> — > Inflow DQ NN NIN NNN Figure 4.2: A Duct With a Step To demonstrate how composite grids may be used to discretise a supersonic flow problem with a steady solution, we look at the problem of supersonic flow in a duct with a step. This problem is similar to one considered by Woodward and Colella [24]. In particular we consider a duct extending from —7 to 3m in the x-direction and from 0 to 7 in the y-direction, with a curved step extending from (0,0) where it intersects the bottom of the duct at right angles, to a point near (37, 7/2) at the downstream end of the duct. The shape of the step is described by x+y = log(1 + 7¢s) for s>0. We consider the case of supersonic inflow so that we may specify all components of the flow at the upstream end of the duct. In particular we specify p = 1.4, u = 3.0, v = 0.0 and p= 0.2 at g = —7. At the downstream end the step narrows the duct 51 to approximately half its full height, so we may anticipate that the outflow velocity is higher than the inflow velocity. We assume that the flow at the downstream end is supersonic and we do not specify any component of the flow at outflow. We specify that there is no normal flow at solid boundaries. With these boundary conditions the flow evolves to a steady solution independent of the initial condition. The main feature of this steady flow is a shock intersecting the bottom wall of the duct upstream of the step and extending to the top wall of the duct, where it is reflected. The reflected shock becomes weaker toward the downstream end and reflects once more from the step before it passes out of the duct. The Euler equations for compressible flow in two dimensions are Ou OF OG _ bt Oe? dy 7° (4.2.1) where p pu, pv u= ||, Fs P pe. , G=], fn , (4.2.2) € u(p + e) v(p +e) p =density, p =static pressure, (u,v) =velocity and e =internal energy related to the pressure by the equation of state e = p/(y — 1) + o(u? + v?). We would like to use the Euler equations 4.2.1 to model supersonic flow in a duct with a step. However, several shocks develop in this flow, and the corresponding solution of the Euler equations becomes discontinuous. Once discontinuities develop in the solution, these equations can no longer be used. To circumvent this difficulty, we add a dissipation term Ou OF OG op op 2 4.2.3 31 + Da + By eV‘u ( ) and consider the limit « —> 0. For small ¢ the dissipative term has almost no effect where u is smooth, but acts to smooth out shocks, where u varies rapidly in space. To see this, consider the case of a steady solution independent of y, and suppose it has a layer of width 6 that approximates a shock. That is, |du/dz| = O(|u|/6) = _5Q- O(6~3) and |dF/dz| = O(|F\/6) = O(6"*) in the layer. This solution satisfies the equation dF /dz = ed?u/dx? so in the layer we have |du/dz| = O(|F|/e) = O(e7’). This shows that shocks are smoothed out to layers of width e. The parabolic system of equations 4.2.3 requires boundary data for all compo- nents of u on the entire boundary. The extra boundary conditions are chosen so that they do not introduce boundary layers into the solution where there were none in the solution to the Euler equations 4.2.1. In particular, we bound the normal derivative of each component for which we have no other boundary condition. To bound the normal derivative of a component of the solution at the boundary, we set the value of the component at the boundary points of each grid equal to its value on the next line of gridpoints in the interior on the same grid. The boundary conditions we use are given in table 4.2. At Inflow At Outflow On Walls p=14 Op/On bounded Op/On bounded u = 3.0 O(u,v)/An bounded (u,v)-n=0 v = 0.0 dv/dn bounded O((u,v) -t)/An bounded p = 0.2 de/On bounded Oe/On bounded Table 4.2: Boundary Conditions for Equation 4.2.3 We discretise the parabolic system 4.2.3 on a composite grid using leap-frog time differencing with lagged dissipation and centered second-order differences in (r,s) on each grid. We use this time-differencing as a relaxation method to find a steady solution. Biquadratic interpolation is used in the regions of overlap between grids. An accurate solution can be obtained if the grids are fine enough to resolve the shock layers smoothly. If the positions of the layers are unknown then fine grids must be used everywhere and the total number of gridpoints will be O(e~’). However, if the positions of the layers are known a priori then we can construct a composite grid with fine grids only where they are needed. This greatly reduces the -53- total number of gridpoints. We use an iterative procedure to calculate the steady solution to 4.2.3 in the limit « — 0. At the first step in the iteration we take a modest value for « and a reasonable initial condition such as the potential flow p=1.4 2 iy =34 —_ + U-w=3st Dae Fy) p= 0.2 chosen to satisfy the boundary conditions. Using a composite grid composed of only coarse grids we run the problem until the solution approaches a steady state. For each subsequent step of the iteration we construct a new composite grid composed of two coarse grids and one or more fine grids in the regions where the steady solution from the previous step was least accurately resolved. Then we take a smaller value for €, interpolate the final state from the previous step onto the new composite grid as an initial condition and run the problem again until the solution approaches a new steady state. These steps are repeated until a target value of € is reached. In particular, figures 4.3, 4.4 and 4.5 show the composite grids used and the steady states obtained with « = 0.5, « = 0.02 and « = 0.015, respectively. -54—- weN [ (= SS Fig. 4.3: Composite Grid and Steady State for the First Step (with ¢ = 0.5) -55— Huy i Fig. 4.4: Composite Grid and Steady State for the Second Step (with € = 0.02) —56- $0801] i sans Fig. 4.5: Composite Grid and Steady State for the Third Step (with « = 0.015) —~57- APPENDIX 1 CMPGRD: A Composite Grid Construction Program* Al.1 Purpose For many problems that involve the discrete representation of a function on a two-dimensional region, a single rectangular grid is inappropriate. For example, a rectangular grid is poorly suited to the discretisation of a PDE boundary-value problem on any region other than a rectangle, because the boundary of the region will not in general correspond to a gridline. One alternative is to cover the region with several overlapping curvilinear grids, each of which is the image, under a smooth transformation, of a rectangular grid, and chosen so that the boundary of the region is composed entirely of edges of the curvilinear grids. Using enough of these transformations and grids it is possible to discretise problems on complicated regions for which a single grid is inadequate. CMPGRD is an interactive FORTRAN program for generating composite grids in the plane, on a cylinder or on a torus. Any compact region with a smooth boundary can be covered with a composite grid, even a multiply-connected region. A composite grid for a region D in the plane is a set of smooth transformations from the unit square to regions in the plane which overlap to cover D, and uniform grids on the unit square that are mapped by these transformations onto curvilinear grids in the plane. It provides boundary-fitted coordinates and allows stretched coordinates for boundary Jayers and internal layers wherever these are needed. The * A User’s Manual for program CMPGRD -58- above definition of a composite grid extends easily to regions with cylindrical or toroidal symmetry, and to n-dimensional regions. Program-CMPGRD also could be generalised to generate n-dimensional composite grids, but at present it is restricted to two-dimensional regions that can be represented in the plane. A1.2 Elements of a Composite Grid The essential elements of a composite grid for a compact region D in the plane, on a cylinder or on a torus are a set of smooth transformations dy: [0, 1)? TR? chosen such that the union of the images D, of the unit square under d, contains the region D and the union of the boundaries of D; contains the boundary of D, and a set of uniform rectangular grids G, on the unit square. DC\JD, and aDC|JeD,, where k k D, = {dz({r,8):0<r<1,0<s<1}, and t—1 p— J ; . G, = _ r2=|,...,me, 7 = 1,..-,2e?- m,p—1 ne—1 The choice of transformations is free to the extent that there is more than one way to.choose the D, to cover D and there is more than one way to choose a smooth transformation d, under which the unit square is mapped onto D,. The choice of D,, is restricted by the requirement that the transformations d, must be smooth. Just how smooth the d; need to be will depend on the specific application the grid is to be used for. For example, if the grid is to be used to discretise a second-order partial differential equation, the jacobian derivative of d; must be continuous and nonsingular for d,; <€ D. Anticipating this use, we shall consider only composite grids composed of transformations d, with continuous, nonsingular jacobians. A consequence of this restriction is that the boundary of each D; must consist of four smooth sides and four corners with interior angles smaller than 7, but even with this restriction there is considerable freedom in the choice of the regions D,;. The _59- boundary of D consists of smooth closed curves; the boundary of D; consists of segments of these and other smooth curves. The following sections describe how these curves may be specified, how D, may be specified in terms of these curves, and how the transformations d, may be specified in terms of the boundary of D,. For the sake of clarity we will first consider only regions D in the plane. The generalisation to regions on the cylinder or torus is then straightforward. A1.2.1 First A Simple Example Most of the concepts involved in constructing a composite grid for any two- dimensional region are illustrated by the simple case of the unit disc in the plane, D = {(z,y) € R*:2? +y? < 1}. For this region we can construct, for example, a composite grid consisting of two transformations d,(r,s) = ((0.4 + 0.6r) cos 27s, (0.4 + 0.6r) sin 27's)) d,(r,s) = (—0.6 + 1.27, —0.6 + 1.2s) from the unit square [0,1]? onto the regions D, = {(z,y):0.4 <2? + y? < 1.0} D; = {(z,y):—0.6 < x < 0.6,-0.6 < y < 0.6} as shown in Figure Al.1l, and uniform grids ~—-1] 9-1 . . : Gi= ey rt=l,...,7,7 =1,...,38 6 32 7-1 7-1 . . G,= ——— , --— }:1 = 1,...,18,7 =1,...,18 12 12 on the unit square, as shown in Figure Al.2. The two regions overlap, and their union is the region D, so the two grids and transformations are a composite grid for D. This composite grid can be constructed in three steps: —-60- (1) Specify parametrisations on [0,1] for the curves forming the boundary of the images, under the proposed transformations d;,, of the unit square. c4(€) = (—0.6 + 1.2€, —0.6) ci(€) = (cos 27€,sin 27) c5(€) = (—0.6 + 1.2€, +0.6) c2(€) = (0.4cos27€,0.4 sin 27) eo(£) = (—0.6, -0.6 + 1.28) c3(€) = (0.4 + 0.6€, 0) C7(€) = (40.6, —0.6 + 1.2£) (2) Use these curves to specify the transformations dy. One option is to use a simple “Coons patch” |7] transformation to specify di(r,s) = (1 ~ s)es(r) + sea(r) + (1 — r)ea(s) + rei(s) — (1—s)(1 ~ r)es(0) — s(1 — r)es(0) + (1 — 8)reg(1) — sres(1) = [(1 — r)e2(s) + rei(s)] + [es(r) — (1 — r)e3(0) — res(1)] = (1—r)e2(s) + rey(s) = ((0.4 + 0.6r) cos 27s, (0.4 + 0.6r) sin 27:8) d.(r,s) = (1 — s)e4(r) + ses(r) + (1 — r)eo(s) + re7(s) — (1 — s)(1 — r)eg(0) — s(t — r)es(0) + (1 — s)req(1) — sres(1) = |(1 —r)es(s) + re7(s)] + (1 ~ s)ea(r) — (1 — r)e4(0) — rea(1)} + s[es(r) ~ (1 ~ r)es(0) — res(1)] = (1 —r)ec(s) + re7(s) = (~0.6 + 1.2r, -0.6 + 1.2s) (3) Specify the grids G, and G, on the unit square by giving the number of gridlines in each direction on each grid: my = 7, A1.2.2 Regions in the Plane Consider a compact region D in the plane, with a smooth boundary. In order to describe the boundaries of regions D,; whose union contains D and the union of whose boundaries contains the boundary of D, we must define three categories of curves: D, Deg Figure Al.1: The unit square is mapped by d, and d, onto D, and Dp. >” > o Figure Al.2: Grids G; and G2 are mapped onto curvilinear grids. (1) The closed curves CP? which form the boundary of D. The boundary of a simply-connected region consists of exactly one such curve; that of a doubly- connected region consists of two, and so forth. These curves will be referred to as boundary curves. (2) Curves CP of which a segment is the image of two opposite sides of the unit square under a periodic transformation d,. Every d, which maps one side of the unit square onto a segment of one of these curves must be periodic and map the opposite side of the unit square onto the same segment of that curve. These curves will be referred to as branch cuts. -62- (3) Other curves Cj. No segment of one of these curves may lie on the bound- ary of D or be the image of two opposite sides of the unit square under any transformation d,. Every curve used in constructing a composite grid falls into one of these categories. In order to specify the transformations d; in terms of segments of these curves, we must parametrise the curves. That is, we must specify continuous functions c,:{0, 1] + IR? such that For the example of the unit disc we defined parametrisations for one boundary curve CP? , one branch cut C?, and five other curves Cj, C1, C1, Cl and C}. The following three restrictions on the parametrisations c, are neccessary but not sufficient to ensure that the transformations d, may be smooth and have nonsingular jacobian derivatives: (1) The functions c, must have continuous, non-vanishing derivatives. (2) The parametrisation c;, of every closed curve C, must be periodic. Without loss of generality, we may insist that the period be one. (3) Every closed curve C, must intersect a branch cut and the intersection must not be tangential. Without loss of generality, we may insist that such an intersection occur at c,(0) = ¢;(1). In the unit disc example C? and Cj were closed curves. Their parametrisations c) and cy were periodic and they intersected curve C2 at c,(0) = c)(1) = (1,0) and c2(0) = c2(1) = (0.4,0), respectively, at right angles to curve Ch. When we have specified and parametrised enough curves so that every edge of every D, is a segment of some curve, then we are ready to specify the transforma- tions d,. To avoid problems with notation, suppose we have numbered the curves so that segments of curves C; and C, form two opposite sides of D,;, and segments of Cz and Cy, form the other two sides, as indicated in the diagram, so that these ~63- Figure A1.3: Curves Bounding a Region curves intersect at €1( £13) = €3(£31), €1( £14) = €4( £41), €2( £53) = C3(E32), C2(E24) = c4( £42). Then we must find a smooth transformation d, such that di.(r,0) = €1(€13 + (S14 — Ers)r), da (r, 1) = €2(€23 + (E24 — &23)r), d,(0,s) = e3(&31 + (E32 ~ €31)r), da (1,8) = c4(€4a1 + (€42 — €41)r). It may not be clear how best to choose such a transformation. The simplest option is a Coons patch, which is a linear combination of points on the edges of Dx, in particular dx (r,s) = (1 — s)d;(r,0) + sd (r, 1) + (1 — r)d,(0,s) + rd,(1, s) — (1 — s)(1 — r)d, (0,0) — s(1 — r)d,(0, 1) ~ (1 ~ s)rd,(1,0) — srd,(1, 1) A Coons patch will be adequate if the boundary of D, is not too convoluted and the parametrisations of the curves are appropriate; otherwise, the transformation may contain folds. This problem can always be avoided by a careful choice of enough regions D, so that the shape of each is simple enough that a Coons patch is adequate. -64- Another option is a conformal map. A conformal map from the unit square onto D, can always be found, but it is not generally possible to ensure that the corners of the unit square correspond to the corners of D;, under such a transformation. Even when this is possible, the transformation will have a singularity at any corner of D, that does not form a right angle. Related but more versatile methods exist [23], [21], [4] to determine a smooth transformation as the solution to a system of elliptic equations. However, if the boundary of D, is too convoluted, these methods will result in a highly non-uniform transformation. As a result, these methods have in practice at most limited advantage over the less sophisticated Coons patch. For the example of the unit disc we used Coons patches for both dy and dz, but conformal maps would have given exactly the same transformations in each case. A1.2.3 Regions on a Cylinder or Torus In case D lies on a cylinder or torus, we represent these as a strip and a box, respectively, in the plane. The period strip of the cylinder is specified by a pair of curves (C?,C?) parametrised by two functions c, which differ by a constant vector ci(€) _ c2(€) = C12 and have continuous, non-vanishing derivatives. The parametrisation should be chosen so that c,(0) and ¢;(1) lie on the boundary of D. The period box of the torus is similarly specified by two pairs of curves (C?,C?) and (C2, C71), each pair parametrised by two functions c, which differ by a constant vector c1(€) —e2(€) =r, 3(€) — e4(€) = Cag and have continuous, non-vanishing derivatives. In this case, the derivative of c} (£) should be periodic with period one acy. dc, —“(€+1)=—= b=)... -65- and the curves should intersect only at the four points €1(0) =¢3(0), ¢1(1) = ¢4(0), c2(0) =cs(1), ¢2(1) = e4(1). In most cases it is convenient for the curves specifying a period strip or period box to be straight lines. For example, if D is periodic in z with period a, bounded below by y = 0 and bounded above by y = 1 then we could specify the period strip by curves (C?,C?) parametrised by c1(€) = (0, é), c2(€) = (a, é). These curves do not at present fit into the category of branch cuts (those curves indicated by the superscript P) defined earlier. In order to accommodate them, that category is extended to include pairs of curves of which corresponding segments are images of opposite sides of the unit square under a transformation d, whose derivative is periodic. A1.3 PDEs on Composite Grids The primary use of composite grids is the discretisation of partial differential equations. With appropriate difference methods they can be used to discretise ellip- tic, parabolic or hyperbolic PDE boundary-value problems. While a discussion of difference methods is beyond the scope of this work, the discretisation of a partic- ular elliptic boundary-value problem is given as an example. Aside from that, the discussion is limited to how a PDE boundary-value problem may be expressed in cartesian coordinates on the unit square using the transformations of a composite grid. A1.3.1 A Simple Example Before looking at more general problems, it is instructive to consider the simple -66- example of Poisson’s equation on the unit disc with Dirichlet boundary conditions: Aus OPu ax? Oy? u=g(z,y) on OD = {(z,y):27 + y’ = 1}. = f(x,y) on D = {(z,y):2’ +y’ < 1} This problem may be discretised using the composite grid discussed-earlier, in par- ticular the transformations di (r,s) = ((0.4 + 0.6r) cos 278), (0.4 + 0.6r) sin 27s)) d2(r,s) = (—0.6 + 1.2r, -0.6 + 1.28) from the unit square onto the regions D, = {(z,y):0.4 < 2? +y’ < 1.0} Dz = {(z,y):-0.6 < x < 0.6,-0.6 < y < 0.6} a~-1ljg-1 ' . G.- pe Jt =1,...,7,7 =1,...,33 i( 6° 32 ] a~-1ly-l1\. . c= { (GF ian ian...n3} on the unit square. We represent the solution u by functions u, and wu, on the unit and grids square ui(r,s) = u(dy(r,s)), ue(r,s) = u(de{r, s)) and write the boundary value problem in terms of u; and uz. In the polar coordi- nates defined by d,, Poisson’s equation becomes 2 . 2 1 Oru, 1 Ou, 1 O*u; — Fldilr.s (a5) Or? (caarrees) ar | (503 + rat asi ~ F{dilr,s)) and in the scaled coordinates defined by d, it becomes 1\7 (@?u, Ou (5) (Se + 5) = Hailes), The boundary condition becomes ui(1,s) = g(di(1,s)). —G7— In order to complete the specification of the problem in (r,s) coordinates and to make sure that wu, and u2 agree in the region of overlap, we apply matching condi- tions us(r,8) = w2(45*(di(r,5))) on the inner edge of D,, uz(r,8) = wi (4; (da(r,5))) on all four sides of D2, and the periodicity condition Ou u3(r, 0) = ui(r, 1), a, (10) = 35 1) for 0 <r <1. A well-conditioned discretisation of the problem can be obtained using the grids G, and G, on the unit square. Using centered second-order dif- ferencing for the PDE and bilinear (second-order) interpolation for the matching conditions, a solution can be obtained which is second-order accurate in the limit that the number of gridlines in both directions on both grids tends to infinity while the two transformations remain fixed. This statement is misleading, since in prac- tice as finer grids are considered, the transformations will be changed in such a way that the amount of overlap between grids decreases inversely as the number of gridlines. In this limit biquadratic (third-order) interpolation must be used to obtain a second-order accurate solution [12, §5.1]. A1.3.2 More General PDEs Consider an arbitrary PDE boundary value problem 0 Oa where F represents the PDE in the interior of D and the boundary conditions on the boundary of D, and it is understood that higher derivatives may be involved. Suppose we have a composite grid for D consisting of n transformations d, from the unit square to regions D, = {di (r,s):0 <r <1,0<s< 1}. For each k we can define a function uz, which represents u on the part of the unit square which is mapped by d, onto D, uz(r,s) = u(d;(r,s)) for dk(r,s) € D. If the jacobian of the transformation d, is non-singular we can use it to express the problem F in (r,s) coordinates. Let Jlrs) i O(r,s) — [* Us , Then the inverse of J, (r,s) is 1 _ -1 _ Ye Le _ Ter Ty 7. (758) - ZrYs — YUr@s ( Ur Ly 7 ( Se Sy ; With the understanding that by r,, r,, 8, and s, we mean the functions of (r,s) which form the inverse of the jacobian, we can express the boundary value problem in (r,s) coordinates as 0 0 0 oO F (-$. + S25 Tua + a) =O for d,{r,s)€D. To specify the problern completely we will generally need to impose matching conditions in the regions of overlap between the D,. In order to do this in a systematic way it helps to be able to order the transformations d, so that wherever two or more of the D, overlap at the boundary of D, this occurs at the boundary of the D, with the highest index & of those involved. Although it is not always possible to order a given set of transformations in this way, it is always possible to choose an appropriate set of transformations that can be so ordered. To help in deciding on an appropriate sequence of transformations, imagine drawing D, in the plane, next _69- drawing D», erasing that part of D; which overlaps D2, and so on, finally drawing D,, and erasing those parts of each underlying D, which it overlaps. The result will be that. the entire region D,, remains, along with those parts of each underlying D; which are not covered by some D; with j > 2. If the sequence was chosen properly, the boundary of D is composed only of visible parts of the boundaries of the D,. Now we are ready to specify the neccessary matching conditions. First we specify conditions identifying u, with u, on regions of overlap between D,, and the underlying regions D,. Where D,, overlaps other regions we can apply the matching condition Un(r,s) = Uz (a; "(an (r, s))), where k is the highest index of those regions which underlie D, at d,(r,s). The matching condition is applied on those sides of D, which are specified in terms of curves C/ (other than boundary curves and branch cuts). Depending on the type of PDE involved, some higher derivatives of u,, and uz, may also need to be matched in the same way. In this way D,, is divided into two parts: the first part where the PDE and boundary conditions are applied and the second part a curve where the matching condition is applied. Matching conditions for the underlying regions are similarly specified. On region D; we apply the condition u(r, $) = u; (4;*(a.(r, s))), where 7 is the highest index of those regions which overlap or underlie D, at d;(r, s). The matching condition is applied on a curve (or curves) chosen so as to divide D; into three parts: the first part where the PDE and boundary condition are applied, the second part a curve where the matching condition is applied, and a third part which includes any of D; that is exterior to D. This curve will lie along sides of D,; which are specified in terms of of curves C’ (other than boundary curves and branch cuts) wherever ¢ > 7, while for 7 < 7 the curve lies on that part of D, where the PDE and boundary condition are applied. By induction we can demonstrate ~70— that it is possible to cover D with overlapping parts of all the D, where the PDE is applied, each of which is separated from parts exterior to D by parts (curves) where a matching condition is applied. At this point the problem F is completely specified in (r,s). Al.4 Using Program CMPGRD CMPGRD is an interactive FORTRAN program for generating composite grids in the plane, on a cylinder or on a torus. It can generate a composite grid for any compact region with a smooth boundary, or even with a boundary that is only piecewise smooth (that is, a region with corners). The program and the instructions given in this manual both rely heavily on the use of an interactive graphics terminal, but it is possible with practice, patience and a good sketch of the region, to use the program without a graphics terminal or even in a batch mode. Specification of a composite grid is accomplished in three steps. First, curves which form the boundary of the region and of the individual grids are specified. Second, grids and transformations which form the composite grid are specified in terms of these curves. Third, the composite grid is specified in terms of these grids and transformations. Although composite grids are defined in terms of smooth transformations and the program also represents them internally as such, their primary use is in discretising PDE boundary value problems. For this purpose CMPGRD produces as output data a discrete description of the composite grid. This data may then be used by other programs to provide a finite-difference approximation to a PDE boundary value problem. Al.4.1 Getting Started Program CMPGRD is designed to be used interactively. At every step where it requires input it either gives instructions or supplies a menu of possible responses. -71- When the program starts, it gives the prompt and waits for a response. If you type HELP then the program lists the menu COMP Specify a new composite grid. CURVE Specify a new curve. EXIT Exit after specifying one or more composite grids. GRID Specify a new grid. HELP Type this list. PARAMS Specify new parameters. QUIT READ Read command input from a file. SAVE Save command input in a file. and repeats the > prompt. This menu lists all of the commands accepted after the > prompt. Any other response will result in an error message suggesting that you type HELP. The most important commands (other than HELP, of course) are COMP, CURVE and GRID. The prompt will change to COMP>, CURVE>, GRID> or PARAMS> when the corresponding command is entered; instructions are given after other commands whenever input is required. Use of the commands accepted after the > prompt is explained below. COMP When one or more grids have been specified, this command is used to specify a composite grid for a region as a sequence of grids that overlap to cover the region. If the resulting composite grid is satisfactory then a discrete description of it may be saved in a data file in a form useful for discretising PDE boundary value problems. This command causes the prompt to change to COMP>, al which point a new menu of recognised commands is available. A later section of the manual is devoted to the use of this command. ~72~ CURVE Curves used as the sides of grids may be specified either by a subroutine EXIT GRID HELP PARAMS QUIT READ SAVE supplied by the user or interactively using the cursor of a graphics terminal. This command is used to specify spline curves through data points entered by the graphics cursor. This command causes the prompt ‘to change to CURVE>, at which point a new menu of recognised commands is available. A later section of the manual is devoted to the use of this command. This command is the normal way to finish a session with CMPGRD. If no composite grid data has been saved in a file yet an error message to that effect is typed and the > prompt returns. When several curves have been specified, this command is used to specify a transformation from the unit square onto a region bounded by segments of these curves and a uniform grid on the unit square that corresponds to a curvilinear grid on this region. This command causes the prompt to change to GRID>, at which point a new menu of recognised commands is available. A later section of the manual is devoted to the use of this command. This command types a menu of recognised commands with a brief explana- tion of each. This command is used to set the size of the plotting area and parameters which determine the action taken by CMPGRD when it detects errors. This command causes the prompt to change to PARAMS>, at which point a new menu of recognised commands is available. A later section of the manual is devoted to the use of this command. This command causes a session with CMPGRD to end immediately. This command is used to read command input from a file saved with the SAVE command in a previous run of CMPGRD. This command is used to save command input in a file for use in later runs of the program. The resulting file may be edited to remove erroneous input -73- or to insert other command input. The command input files saved from several runs may be combined and used as input to a subsequent run. In this way the individual curves and grids used to specify a composite grid may each be specified in separate runs of the program, followed by a run that specifies a composite grid. A1.4.2 Specifying Curves Curves are divided into three categories according to how they are used in specifying the grids and transformations of a composite grid. Some curves form parts of the boundary of the region to be covered by the composite grid. These are called boundary curves. Segments of some curves or pairs of curves are used as opposite sides of grids on which the problem is periodic. These are called branch cuts. Segments of all other curves are used only as sides of grids where they overlap other grids. Each side of a grid must be a segment of a curve that belongs to only one of these categories. Adjacent sides of a grid must be segments of different curves, so the endpoints of these segments are specified as intersections of curves. To ensure that these endpoints are well-defined and to avoid the possibility of highly-skewed grids, these intersections may not be tangential. CMPGRD provides two methods by which curves may be specified. One is to to specify a cubic spline curve passing through data points entered using the cursor of a graphics terminal or from a data file. The other is to use a subroutine (supplied by the user) which specifies parametrisations for one or more curves. Spline curves have the advantage that they are easy to use in modelling regions whose boundaries cannot be described in terms of known functions. Curves specified by subroutine have the advantage that their parametrisations may be arbitrarily smooth functions, while cubic splines have discontinuous third derivatives. Grids may be specified in terms of curves specified by either method, and segments of curves specified by _74- different methods may be used as different sides of the same grid. Curves used as opposite sides of the period strip of the cylinder or the period box of the torus must occur in pairs. Spline curves may not be used for this purpose since they are specified only one at a time. Subroutine CURVE must be supplied by the user. If the parametrisations of any curves are to be specified by subroutine, then CURVE gives these parametri- sations. CURVE also gives the number of such curves, so even if there are none subroutine CURVE must still be supplied. At the beginning of a session, CMP- GRD calls subroutine CURVE to find how many curves it specifies. Then it checks that the parametrisation of each curve has a non-vanishing derivative and that the parametrisations of pairs of branch cuts have identical derivatives. General specifi- cations for subroutine CURVE and two examples are given in an appendix. The CURVE command in the main menu (after the > prompt) is used to specify curves interactively. This command changes the prompt to CURVE>, at which point the curves currently specified are plotted for reference and the following menu of commands is available: EXIT Exit after specifying a curve. HELP Type this list. QUIT SPLINE Specify a spline curve. This menu lists all of the commands accepted after the CURVE> prompt. The use of these commands is explained below. EXIT This command is the normal way to return to the main menu and the > prompt after specifying a curve. If no new curve has been specified yet, an error message is typed and the CURVE> prompt returns. HELP This command types a menu of recognised commands with a brief explana- tion of each. QUIT This command abandons the curve currently being specified and returns to SPLINE -~75- the main menu and the > prompt. This command is used to specify a spline curve through data points entered interactively using the cursor of a graphics terminal. A spline curve is spec- ified in three steps. First, CMPGRD asks whether the curve is a boundary curve, a branch cut or a curve used only for interpolated sides of grids. Sec- ond, CMPGRD asks if the problem is periodic on the curve. In particular, the problem is periodic on any curve which joins corresponding points on a pair of branch cuts. To help avoid the mistake of trying to specify a grid two opposite sides of which are incompatible segments of a pair of branch cuts, CMPGRD will not allow a non-periodic curve to begin or end on a branch cut. To help avoid the mistake of trying to specify a grid one side of which is a segment of a closed curve which is not parametrised continuously, CMPGRD will not allow a closed curve to begin or end anywhere except on a branch cut. Third, CMPGRD asks for points through which the spline passes. Four or more points are needed to specify a closed curve; other curves need at least two. Each point chosen may be either the cursor point or the nearest point on another curve, with the restrictions noted above. Af- ter all the points on a periodic curve have been entered, CMPGRD checks that the endpoints are compatible; the first point lies on a branch cut, so CMPGRD moves the last point to the corresponding point on the opposite branch cut. A1.4.3 Specifying Grids and Transformations The curvilinear grids which make up a composite grid are the images under smooth transformations of uniform grids on the unit square. To specify a grid we must specify a uniform grid on the unit square in the (r,s) plane and a transforma- tion from the (r,s) plane to the (x,y) plane. The uniform grid on the unit square is -16- specified by the number of gridlines in each direction. The transformation from the (r,s) coordinates to (x,y) coordinates is specified in two stages. First, to allow for stretched coordinates we specify a transformation from the unit (r,s) square onto the unit (t,u) square. This transformation will be called the stretching function. Second, we specify a transformation from the unit (t;u) square onto a region D, in the (z,y) plane bounded by segments of curves. This transformation will be called the mapping function. The transformation from the unit (r,s) square onto D, is the composition of the stretching and mapping functions. Unless otherwise specified, the stretching is taken to be the identity function and the transformation is just the mapping function. Although the transformation could be specified in a single step, these two steps make stretched coordinates for boundary layers and interior layers easier to tailor to the needs of a specific problem. CMPGRD provides a choice of two types of stretching functions from the unit square in the (r,s) plane to the unit square in the (t,u) plane. One option is the identity function (t,u) = (r,s). The other option is an exponential stretching function. In this case t is a function of r only and u is a function of s only, and both of these functions have the same form. The inverses of these functions are weighted sums of a linear function and one or more terms involving hyperbolic tangent functions. In particular, ¢t is related to r by the function r(t) = t+ Dolo: (t) mA) 14+ SO le: (1) - 2: (0)] where o;(t) = a; tanh b;(t — ¢;) if the problem is not periodic in r, and o,(t) = a; tanh b;(t — ¢;) + Ss, a;{tanh b,(t — ¢; — 7) + tanhb,(t — ¢; + 3)] j=l if the problem is periodic in r. Each term o;(t) has three free parameters: the weight a;, the exponent }; and the location c;. Each term produces a stretching _77- in the t coordinate centered at t = c;. That is, the image in the (t,u) plane of a uniform grid in the (r,s) plane is a nonuniform rectangular grid with a maximum in density of gridlines for t near each c;. The density of gridlines in the t-direction at the point (t,u) is proportional to the derivative r'(t), and may be adjusted by varying the parameters a;, 6; and c; in each term. In particular, a; is approximately the ratio of the number of gridlines in the stretched region of the grid near c; to the total number of gridlines in unstretched regions, while 6; is approximately the ratio of the density of gridlines near c; to the density of gridlines in unstretched regions. If the problem is periodic in r, an infinite series of correction terms is needed to ensure that the derivative t'(r) of the stretching is periodic. For example, if there is oné term in the stretching (n = 1), and we choose a; = 1, b} =1 and c, = 0, the ¢ coordinate is stretched at t = 0, as shown in Figure A1.4; if the problem is periodic in r then the t coordinate is equally stretched at ¢ = 1. Another sum of the same form expresses s in terms of u, and together these specify (t,u) = (t(r), u(s)), where t(r) and u(s) are the inverses of functions given by weighted sums of linear and hyperbolic tangent terms. CMPGRD provides only one type of mapping function from the unit square in the (¢,u) plane to the region D; in the (x,y) plane, namely a Coons patch [7]. This mapping is a linear combination of points on four segments of curves which form the boundary D,. To avoid problems with notation, suppose we have numbered the curves so that segments of curves C; and C, form two opposite sides of the region D,, and segments of Cs and Cy, form the other two sides, so that these curves intersect at €1(£13) = eg(€31), 1 (E14) = C4 (Ear), C2( £23) = €3(€32), €2( E24) a C4 (E42). —78- t t 1 | 14 + ] | ] | 0 3 +—— +—+ > Ot —r 0 1 0 1 Non-—Periodic Periodic Exponential Stretching Functions n=1, a,=1., by, =5., ce, =0. Figure Al.4: Examples of Stretching Functions CMPGRD uses the Coons patch defined by (x,y) = (1 ~ u)er(€1s + (S14 — €1s)t) + we2(€23 + (E24 — £23)t) + (1 — thes(€ai + (&32 — €si)u) + teg( Ear + (S42 ~ &4i) u) ~ (1 —u)(1 — t)ei(fis) — wd — t)e2(&2s) ~ (1 — u)ter (E14) — utea( fos) Since the parametrisations of the curves are smooth functions this mapping is also a smooth function. However, the constraint on the curve parametrisations that their derivatives not vanish does not imply that the jacobian derivative of the mapping is nonsingular or that the mapping is locally invertible. For example, if one of the sides of D,; has zero length then the jacobian will be singular on that side. In this case the resulting grid may still be useful for discretising a PDE boundary value problem if this singularity is taken into account. However, even if none of the sides has zero length the mapping may have folds, in which case it is not locally invertible and the resulting grid cannot be used for discretising a PDE. If this problem occurs, _79- one remedy is to subdivide the region D, into two or more overlapping regions and try again. Despite this difficulty, the Coons patch mapping is still a useful tool. The GRID command in the main menu (after the > prompt) is used to specify grids and transformations. This command causes the prompt to change to GRID>, at which point the curves currently specified are plotted for reference and the following menu of commands is available: EXIT Exit after specifying a grid. HELP Type this list. - LINES Specify the number of gridlines in each direction. MAPPING Select a mapping function (t,u) --> (x,y). QUIT SEGMENT Specify a segment of a curve. SIDE Associate a side of the grid with a segment. STRETCH Specify a stretching function (r,s) --> (t,u). This menu lists all of the commands accepted after the GRID> prompt. The prompt will change to MAPPING> or STRETCH> when the corresponding command is entered; instructions are given after other commands whenever input is required. Use of the commands accepted after the GRID> prompt is explained below. EXIT This command is the normal way to return to the main menu and the > prompt after specifying a grid. If no new grid has been specified yet, an error message is typed and the GRID> prompt returns. | HELP This command types a menu of recognised commands with a brief explana- tion of each. LINES This command is used to specify the numbers m, and m, of gridlines in the r and s directions respectively on the unit square. MAPPING This command is used to select a mapping function and check that its jacobian is nowhere singular. At this time there is no choice of mapping function; a Coons Patch is always used. This command causes the prompt to change to MAPPING>, at which point the following menu of recognised commands is available: ~80- COONS Use a Coons Patch mapping function. EXIT Exit after specifying a mapping function. HELP Type this list. The COONS command is used to specify the Coons patch mapping function described above. If this mapping has a singular jacobian anywhere on the unit square then CMPGRD types a warning message to that effect. QUIT SEGMENT SIDE This command abandons the grid currently being specified and returns to the main menu and the > prompt. This command is used to specify a segment of a curve that forms one or more sides of the region in the (z,y) plane that the unit square is mapped onto. The curve and the endpoints of the segment on the curve are specified interactively using the cursor of a graphics terminal. First, the curve is chosen by entering a cursor point on or near it. The curve which passes nearest to the cursor point will be chosen. Second, the endpoints of the segment are chosen one at a time by entering cursor points at or near them. If the nearest curve to the cursor point is the curve already chosen then the endpoint will be the end of the curve nearest, in terms of its arcwise parameter, to a nearby point on the curve. Otherwise the endpoint will be a nearby intersection of the chosen curve and the curve nearest to the cursor point. This command is used to specify which segment of a curve is associated with a side of the unit square in the (r,s) plane. CMPGRD displays the unit square and the segments that have so far been specified. First, a side of the unit square is chosen by entering a cursor point on or near it. Next, the associated segment is chosen by entering a cursor point on or near it. Finally, the orientation of the segment is reversed, if necessary, so that each endpoint of the side corresponds to the proper endpoint of the segment. -81- STRETCH This command is used to select a stretching function. The only available stretching function other than the identity is an exponential stretching. This command causes the prompt to change to STRETCH>, at which point the following menu of recognised commands is available: EXIT Exit after specifying a stretching function. EXPO Specify an exponential stretching function. HELP Type this list. IDENT Use the identity stretching function. The EXPO command is used to specify an exponential stretching of the type described above. CMPGRD first prompts for the direction (r or s) in which stretching is to be specified, then for the number of terms, and finally for the three coefficients a,, b; and ¢; (discussed in §A1.4.3) for each term 7. If the problem is periodic then the exponent in each term of the stretching in the periodic direction must be at least 1.0 so that the infinite series of correction terms converges quickly. This command must be repeated if stretching in both r and s is required. A1.4.4 Specifying Composite Grids A composite grid for a region D is a set of smooth transformations from the unit square onto regions which overlap to cover D, and uniform grids on the unit square which are mapped by these transformations onto curvilinear grids. To be useful in discretising PDE boundary value problems, the grids and transformations comprising a composite grid must be given a hierarchy. That is, the composite grid is specified as a sequence of grids chosen so that wherever two or more grids overlap at the boundary of D, the boundary corresponds to an edge of the grid that occurs latest in the sequence. Think of the grids as being laid down in sequence, so that the second grid obscures part of the first grid, the third grid obscures parts of the first and second grids, and so forth. All of the last grid in the sequence remains visible. If the sequence has been chosen correctly the boundary of D will consist —~82~ only of parts of visible edges of the grids. The COMP command in the main menu (after the > prompt) is used to specify a composite grid and to save a discrete description of it in a data file in a form useful for discretising PDE boundary value problems. This command causes the prompt to change to COMP>, at which point the following menu of commands is available: COORD Find coordinates of a point on composite grid. EXIT Exit after writing out data for composite grid. GRIDS Specify the sequence of overlapping grids. HELP Type this list. QUIT SAVE Save composite grid data in a file. This menu lists all of the commands accepted after the COMP> prompt. The use of these commands is explained below. COORD EXIT GRIDS This command prompts prompts for cursor input of points (z,y) on D and returns (r,s), such that (x,y) = d,(r,s) and & is the highest index of any grid from which the point can be interpolated. The results are also put into the log file. This command is the normal way to return to the main menu and the > prompt after specifying a composite grid and saving data for it in a file. If data for a composite grid has not been saved yet, an error message is typed and the COMP> prompt returns. This command is used to specify a sequence of overlapping grids that make up a composite grid for a region D. First, CMPGRD prompts for the number of grids in the composite grid. Next, it prompts for the sequence of grids, which must conform to the principle explained at the beginning of this section. Next, it prompts for the order of the PDE discretisation and of the interpolation. The order of the discretisation must be an even positive integer, and it determines the size of the computational molecule. For example, if the order of the discretisation is two then the computational HELP QUIT SAVE -33- molecule will consist of nine points, while if the order of discretisation is four then the computational molecule will consist of 25 points. The order of the interpolation must be a positive integer, and it determines the type of interpolation and the number of points needed for interpolation. For example, if the order of the interpolation is two then bilinear interpolation can be used, which involves four points, while if the order of the interpolation is three then biquadratic interpolation will be used, which involves nine points. After the number of grids, their sequence and the order of the discretisation and of the interpolation is specified, CMPGRD determines which gridpoints will be used for the PDE and which for the matching conditions. If a valid sequence of grids was specified, all other gridpoints are either redundant or lie outside the region D. Otherwise, the regions of overlap do not separate all of the parts of the grids where a PDE could be discretised from parts that are exterior to the region D. CMPGRD plots those parts of each grid on which the PDE or the matching condition will be applied, leaving blank all unneeded parts. If the grids were specified in an invalid sequence or there was insufficient overlap between the grids then CMPGRD types an error message which reports that one or more of the grids will be null or that points needed for interpolation were deleted. In this case the composite grid cannot be used to discretise a PDE, and the sequence of grids must be respecified. This command types a menu of recognised commands with a brief explana- tion of each. This command abandons the composite grid currently being specified and returns to the main menu and the > prompt. This command is used to save data for the composite grid in a file for later use by programs that discretise PDI! boundary value problems. The format —84~ of the data in such a file is described in an appendix. A1.4.5 Specifying Parameters The PARAMS command in the main menu (after the > prompt) is used to set the size of the plotting area and parameters which determine the action taken by CMPGRD when it detects errors. This command causes the prompt to change to PARAMS>, at which point the following menu of recognised commands is available: BOUNDS Specify new bounds for the plotting region. DEBUG Specify which items to save debug data for. ERROR Specify a new error handling procedure. EXIT HELP Type this list. This menu lists all of the commands accepted after the PARAMS> prompt. The use of these commands is explained below. BOUNDS DEBUG EXIT HELP ERROR This command is used to specify new bounds for the plotting region. Ini- tially, bounds of the plotting region are determined by extrema of the curves defined by subroutine curve, supplied by the user. These bounds may be inadequate if more curves are to be specified. This command is used for debugging. Internal variables used in generating a composite grid may be saved in a readable form in the log file. In addition, the common blocks may be saved when the program exits. This command is used to return to the main menu and the > prompt. This command types a menu of recognised commands with a brief explana- tion of each. This command is used to specify a new error handling procedure. Initially the procedure for handling errors is to crash the program (to get a traceback) only on fatal errors. Other errors (such as unrecognised commands) cause an error message to be typed and contro] to return to the terminal if a -85- command input file was being read at the time. Warnings are typed, but commands continue to be taken from a command input file if any was being read at the time. Informational messages are ignored. If the program is run as a batch job it is best to set the program to crash on all errors. A1.5 Appendices The following appendices give specifications and examples of subroutine CURVE, which must be supplied by the user; of the graphics package, which is specific to the graphics terminal and operating system; and of the composite grid data file which CMPGRD produces for later use by programs that discretise PDE boundary value problems. A1.5.1 Subroutine CURVE Subroutine CURVE must be supplied by the user. The details of the subroutine will vary with the geometry and topology of the region for which a composite grid is to be generated. General specifications and two examples are given here. subroutine CURVE(k, s, x,y, xs, ys, kp) Purpose: Specify smooth curves c, for 1 < k < kmnaz, namely C* functions c,:(0.0, 1.0} IR”, their derivatives and information regarding their use. Arguments on input: k Type: Integer Range: O,1,...,kmazr If k = 0, subroutine CURVE returns as kp the highest number, kmaz, for which the subroutine defines a curve. The first call to subroutine CURVE is with k = 0, and all subsequent calls are with 1<k< kmaz. If k > 1 then k is called the curve number. —86- S Type: Real Range: (0.0, 1.0] (see below) s is the arc parameter of the curve, and should be in the range (0.0, 1.0]. However, sometimes subroutine CURVE is called with s outside this range, and although it may be a lot to ask, the definition of the curve should be extended to accommodate s at least in the range |—0.1,1.1). It is not essential that the-deriva- tive should be continuous on the extended interval. Just keep in mind that Newton’s method is used extensively and is expected to work even when it involves calling subroutine CURVE near the endpoints of the interval [0.0, 1.0]. Arguments on output: xyy Type: Real The value of the kth curve function at s. Ifk <lork > kmaz then xz and y need not be defined. xs,ys Type: Real The derivative of the kth curve function at s. If k <lork> kmazx then xs and ys need not be defined. kp Type: Integer If k = 0, set kp = kmaz, the highest number for which a curve is defined by subroutine CURVE. If 1 < k < kmaz then kp indicates how curve k is to be used: (a) If curve k is used as some part of the boundary of the region to be described, then kp = —1. (b) If the composite grid is generated for a problem which is periodic on any of the component grids and curve k is used as a branch cut on one or two corresponding sides of such a grid, kp is the number of the curve which forms the branch cut on the opposite side. In this case the derivatives of curves k and kp must be equal for each s. {c) If curve k is used only as an edge of a grid where that grid overlaps another grid, then kp = 0. Two examples of subroutine CURVE are given below. The first specifies no curves at all. The second specifies seven curves, three of which may be used to form the boundary of an annular region, including a branch cut, and four of which bound a square. 1 subroutine curve(k,s,x,y,xs,ys, kp) BWW o OMAN OO bk WON eS wow won nN NNNNNNNN EER ER Re RR Re BOoOovuoreoaoanrdaoawx one coumrnanrreronre 35 Woe Tt oO gooenan an Cy he aanaoana ~87- This subroutine specifies no curves at all. kp=0 return end subroutine curve(k,s,x,y,xs,ys,kp) The curves specified here can be used to specify a composite grid for a disc of radius 1 consisting of two grids: - (1) an annular grid on the region between circles of radii 0.4 and 1.0 centered at the origin (2) a rectangular grid on the square [-0.6,0.6]x[-0.6,0.6] parameter (twopi=2.%3.14159265359) if(k.eq.0)then Seven curves are defined. kp=7 elseif (k.eq.1)then Curve 1 is the circle of radius 1.0 centered at the origin. x=cos (twopis) y=sin(twopi*s) xs=-twopi*sin(twopis) ys= twopitcos(twopi*s) Curve 1 forms the boundary of the region. kp=-4 , elseif (k.eq.2)then Curve 2 is the circle of radius 0.4 centered at the origin. x=0.4*cos (twopi*s) y=0 .4*sin(twopi*s) xs=-twopit0.4+sin(twopi*s) ys= twopi+0.4*cos(twopi*s) Curve 2 is the edge of the annular grid which overlaps the rectangular grid. kp=0 elseif (k.eq.3) then Curve 3 joins curves i and 2 by a radial line. x=0.4+0.6%s y=0. xs=0.6 ys=0.0 Curve 3 forms two opposite sides of the annular grid. kp=3 elseif (k.eq.4)then Curve 4 is the line y=-0.6 for -0.6.le.x.le.0.6 x=1.2*8-0.6 a -0.6 =1.2 ys=0.0 Curve 4 is an edge of the rectangular grid which overlaps the annular grid. kp=0 elseif (k.eq.5)then Curve 5 is the line y=0.6 for -0.6.le.x.1e.0.6 x=1.2%s-0.6 . y=0.6 xs=1.2 ys=0.0 Curve 5 is an edge of the rectangular grid which overlaps the annular grid. : kp=0 elseif (k.eq.6)then Curve 6 is the line x=-0.6 for -0.6 le. y.le.0.6 ~88— 56 x=-0.6 57 yol .2%s-0.6 58 xs=0.0 59 ys=1.2 60 ¢ Curve 6 is an edge of the rectangular grid which overlaps the Gioc annular grid. 62 kp=0 G3 elseif (k.eq.7) then 64 ¢ Curve 7 is the line x=0.6 for -0.6.le.y.le.0.6 65 x=0.6 6G y=1.2#8-0.6 67 xs=0.0 G8 ys=1.2 69 ¢ Curve 7 is an edge of the rectangular grid which overlaps the 70°C annular grid. 71 kp=0 72 endif 73 return 74 end A1.5.2 The Graphics Package All graphic input and output is handled by the following routines, which are specific to the graphics terminal and operating system and so must be written by the user. General specifications for these subroutines are given here, along with two examples of implementations of the graphics package. subroutine SETUP (ploter) Purpose: Initialise graphics package and graphics terminal. This subrou- tine needs only to set the value of ploter. However, it may be used for any initialisation required by the graphics package and/or the graphics terminal. Calling sequence: This subroutine is called only once, before any other subroutine in the graphics package is called. ~&RQ— Arguments on output: ploter Type: Integer (a) If the terminal used does not have graphics capabilities then ploter = 0. In this case no other subroutines of the graphics package will be called. (b) If the graphics terminal used is such that parts of a plot cannot be erased by plotting them over in the background colour, then ploter = 1. This is the case for most graphics terminals. (c) If the graphics terminal used is such that parts of a plot can be erased by plotting them over in the background colour, then ploter = 2. This is case for some raster graphics termi- nals. subroutine RESET Purpose: Reset the graphics package and terminal. This subroutine may be used for any resetting required by the graphics package and/or the graphics terminal. Calling sequence: This subroutine is called only once, after which no other subrou- tine in the graphics package is called. subroutine STRTPL Purpose: Put the terminal into graphics mode. Calling sequence: This subroutine is called before each item is plotted. After the call to STRTPL and before the next call to ENDPLT, no al- phanumeric output will be sent to the terminal by subroutine REPLY and no alphanumeric input will be expected by subrou- tine RQUEST of the input/output package. subroutine ENDPLT Purpose: Take the terminal] out of graphics mode. -9Q- Calling sequence: This subroutine is called after each item is plotted. After the call to ENDPLT, alphanumeric output to the terminal from sub- routine REPLY may be sent and alphanumeric input may be expected by subroutine RQUEST of the input/output package. subroutine ERASE Purpose: Erase the graphics area of the terminal. " subroutine WINDOW (tw, za, zd, ya, yb) Purpose: Set up a linear transformation onto the graphics area. CMP- GRD often needs to draw two plots to be seen both at the same time. For this purpose the graphics area is conceptually divided up so that it contains two square regions called the primary and secondary windows. The primary window should take up most of the graphics area; the size of the secondary window is unim- portant, but the two should not overlap. Subroutine WINDOW sets up a linear transformation from [za, x6] x [ya, yb] onto one of two windows. Arguments on input: iw Type: Integer Range: 1,2 If iw = 1 then use the primary (larger) window. If cw = 2 then use the secondary (smaller) window. za,xzb Type: real Range: ra < xb xaand rb are coordinates of the left and right sides of the window under the linear transformation. ya,yb Type: real Range: ya < yb ya and yb are coordinates of the bottom and top of the window under the linear transformation. — -91- subroutine SHADE(sh) Purpose: Set colour, intensity or dash pattern for plotting. If the graphics terminal is capable of plotting in different colours, choose two different colours for sh = 0.5 and sh = 1.0, and use the back- ground colour for sh = 0.0, which is used for erasing. If the graphics terminal is not capable of plotting in different colours or intensities then set up a dash pattern distinct from a solid line for use when sh = 0.5, set the dash pattern to solid for sh = 1.0 and to null for sh = 0. Bear in mind that items plot- ted with sh = 0.5 may be plotted over with sh = 1.0, and if ploter = 2, then items plotted with sh = 1.0 may be plotted over with sh = 0.5 or sh = 0.0 and items plotted with sh = 0.5 may be plotted over with sh = 0.0. Arguments on input: sh Type: real Range: {0.0,0.5, 1.0} If sh = 0.0, set the colour to the background colour, or set the dash pattern to null. If sh = 0.5, set the colour to dim or gray, or set the dash pattern to something distinguishable from a solid line. If sh = 1.0, set the colour to bright or white, or set the dash pattern to solid. subroutine MOVE(z, y) Purpose: Move the “pen” point to the image of (z,y) under the linear transformation set up by the most recent call to subroutine WIN- DOW. The point (x,y) should be in [za, 2b] x [ya, yb], but may be outside in some unimportant cases. Arguments on input: yy Type: real Range: [za,xb) x [ya, yb) (see above) The coordinates of the point. subroutine DRAW(z, y) Purpose: Move the “pen” point to the image of (x,y) as in subroutine MOVE, and draw a line from the previous “pen” point. —-92- Arguments on input: zy Type: real Range: [ca, xb] x |ya, yb] (see above) The coordinates of the point. subroutine MARK(z, y) Purpose: Move the “pen” point to the image of (z,y) as in subroutine MOVE, and draw a small circle or other recognizable symbol there. Arguments on input: Ly Type: real Range: (xa, xb| x [ya, yb) (see above) The coordinates of the point. subroutine CURSOR(z, y, key) Purpose: Obtain a cursor report from the graphics terminal using its joy- stick or equivalent, if any, and return the result (X,Y) under the transformation set up by the most recent call to subroutine window. Also return a one-character key entered from the key- board. For example, on a Tektronix graphics terminal this can be the key that was hit in order to send the cursor report. This subroutine might be tricky to write since it may require that input from the terminal be read which is not terminated by a carriage-return. If this is too difficult or the graphics terminal has no joystick, then a simple way out is to enter (z,y) and the key directly from the keyboard using simple read statements. Arguments on output: L,Y Type: real The coordinates of the cursor under the transformation set up by the most recent call to subroutine window. key Type: character+(*) key(1:1) should be set to a one-character key entered from the the keyboard of the terminal. ‘wo examples of implementations of the graphics package are given below. The first is an implementation for use with no graphics terminal. This is not of much ~93- use, since CMPGRD relies heavily on interactive graphics capabilities. Next is an implementation for use with a Tektronix 4107 graphics terminal. The final two listings contain system-dependent subroutines used by the latter implementation. They are for use on computers running the DEC VMS and BSD 4.2 Unix operating systems, respectively. ic 2c Graphics package for use with no graphics terminal. 3.¢ 4 subroutine setup(ploter) 5 integer ploter 6 ploter=0 7 return 8 entry reset 9 entry strtpl LO entry endplt li entry erase 12 entry window 13 entry shade 14 entry move 15 entry draw 16 entry mark 17 entry cursor 18 end re} Graphics package for Tektronix 4107 terminal. subroutine setup(ploter) character buffer512,str%3 integer ploter,buflen logical alpha common/graph/alpha,xs ,x0,ys,yO, buffer, buflen data alpha,buflen/.true. ,0/ call opterm c Select TEK code, clear dialog scroll and set dialog area to 10 lines. ee Fc OU MN AOC ek WN HK 12 call encint(1i0,str,n) 13 call bufout i4 + (char (27)//'4!0'//char(27)//’LZ’//char(27)//’LL'//str(:n)) 1S call erase 16 ¢ Set ploter=2 for raster graphic display 17 ploter=2 18 return 19 end 20 21 subroutine reset 22 character buffer*512,str*3 23 integer buflen 24 logical alpha 25 common/graph/alpha,xs,x0,ys,y0, buffer, buflen 26 call erase 27 ¢ Clear dialog scroll, set dialog area to 32 lines and select ANSI code. 28 call encint(32,str,n) 29 call bufout ~94— 30 + (char (27)//°LZ’//char(27)//’LL’//str(:n)//char(27)//°%!1') 3] call wrterm(buffer, buflen) 32 call clterm 33 return 34 end 35 36 subroutine strtpl 37 return 38 end 39 40 subroutine endplt 4. ¢ Put the terminal back into alpha mode. 42 character buffer'!512 43 integer buflen 44 logical alpha 45 common/graph/alpha,xs ,x0,ys,y0, buffer, buflen 4G call bufout(char(31)) 47 call wrterm(buffer ,buflen) 48 alpha=.true. 49 return 50 end 51 52 subroutine erase 53. ¢ Erase the graphics area. an rs character buffer#512 integer buflen logical alpha common/graph/alpha,xe,x0O,ys,yO, buffer, buflen Or Cn aa aan a NI 8 call bufout(char(27)//char(i2)) 9 call wrterm(buffer, buflen) GU return 61 end G2 63 subroutine window(iw,xa,xb,ya,yb) G4 ¢ Set up a transformation from [xa,xb]x[ya,yb] onto a square. 65 character buffer*512 6G integer buflen 67 logical alpha 68 common/graph/alpha,xs,x0,ys,y0, buffer, buflen 69 if(iw.eq.i)then 70 ¢ Use the primary plotting region. 7h ia=256 72 ib=2303 73 ja=1024 74 jb=3071 75 else 76 ¢ Use the secondary plotting region. 7 ia=2816 78 ib=3839 79 ja=1536 80 jb=2559 81 endif 82 xs=(ib-ia)/(xb-xa) 83 xO=ia-xa*xs 84 ys=(jb-ja)/(yb-ya) 85 yO=ja-yatys 86 return 87 end 88 89 subroutine shade(sh) 90 ¢ Set the colour to black, light gray or white ~95- character*3 string if(sh.le.O.)then k=0 elseif (sh.1t.1.)then k=15 else k=1 endif call encint(k,string,n) call bufout(char(27)//’NL’//string(:n)) return end subroutine move(x,y) Move the pen point to (x,y) with the pen up. character buffer*512,string#5 integer buflen logical alpha common/graph/alpha,xs,xO,ys,y0, buffer, buflen call enccrd(x,y,string,n) call bufout(char(29)//string(:n)) alpha=.false. return end subroutine draw(x,y) Move the pen point to (x,y) with the pen down. character buffer#512,string*5 integer buflen logical alpha common/graph/alpha,xs ,x0,ys,y0, buffer, buflen if(alpha) then call bufout (char(29)//char(5)) alpha=.false. endif call enccrd(x,y,string,n) call bufout(string(:n)) return end subroutine mark(x,y) Move the pen point to (x,y) and draw a small box there. character buffer#512,stri#3,str2:5 integer buflen logical alpha common/graph/alpha,xs ,x0,ys,y0, buffer ,buflen call encint(8,stri,nt) call encerd(x,y,str2,n2) call bufout (char (27)//’MM’//stri(:n1)//char(27)//'LH’//str2(.n2)) return end subroutine cursor (x,y,key) Get a cursor report. character key«(+*) ,buffer*512, string: (8) integer buflen,hix,lox,sex,hiy,loy,fex logical alpha common/graph/alpha,xs,xO,ys yO, buffer ,buflen Enable GIN. call bufout(char(27)//’IEO1’) rae I] ran co i ~ an a a cn el a be oan a nx 159 160 161 162 163 164 165 166 1G8 is9 170 L7i 191 192 193 - 194 195 196 197 198 200 201 202 203 204 205 206 207 208 209 210 2it 212 -96~— call wrterm(buffer,buflen) Read GIN report. call rdterm(string,len(string)) Reset bypass mode. call bufout(char(10)) hiy=iand(ichar(string(2:2)),31) fex=ishft (iand(ichar(string(3:3)),12),-2) sex=iand(ichar(string(3:3)) ,3) loy=iand(ichar(string(4:4)) ,31) hix=iand(ichar(string(5:5)) ,31) lox=iand(ichar(string(6:6)) ,31) 1=128*hix+d*loxt+sex j=i28thiyt4+loy+fex x=(i-x0)/xs y=(j-y0)/ys key=string(1:1) return end subroutine encint (i, string ,n) Encode an integer. character string*(*) ii=iabs (i) ii=iand(ii,15) if(i.ge.O)il=ior(it,16) i2=iand(ishft (ii,-4) ,63) i3=iand(ishft(ii,-10) ,63) if(i3.ne.0)then n=3 string=char (ior (i3,64))// char (ior (i2,64))//char (ior (i1,32)) elseif (i2.ne.0)then n=2 string=char (ior(i2,64))//char (ior (ii,32)) else n=1 string=char (ior(i1,32)) endif return end subroutine enccrd(x,y,string,n) Encode an (x,y) coordinate. character string*(*), buffer '512 integer buflen,ox1,oy1 logical alpha common/graph/alpha,xs,x0,ys,y0, buffer, buflen data oxl,oyi/2#-1/ i=x*xs+x0 j=y*ys+yO nyi=iand(ishft(j,-7),31) ny2=iand(ishft(j,-2),31) ny3=iand( j , 3) nx1=iand(ishft(i,-7),31) nx2=iand(ishft(i,-2),31) nx3=iand( i , 3) n=0 if(ny1.ne.oy1)then n=n+1 string(n:n)=char(ior(ny1,32)) endif -~97- 213 n=n+1 214 string(n:n)=char(ior(ior(ishft(ny3,2) ,nx3) ,96)) 215 n=n+1 216 string(n:n)=char(ior(ny2,96)) Q17 if(nx1.ne.oxi)then 218 n=n+1 219 string(n:n)=char(ior(nx1,32)) 220 endif 221 n=n+1 222 string(n:n)=char (ior (nx2,64)) 223 oxl=nx1 224 oyt=nyl 225 return 226 end 227 228 subroutine bufout (string) 229 ¢ Buffer a string to the terminal. 230 character string«(*) ,buffer#512 231 integer buflen 23:2 logical alpha 233 common/graph/alpha,xs ,xO,ys,yO, buffer,buflen 234 do 10 i=i,len(string) 235 buflen=buflen+1 236 buffer (buflen:buflen)=string(i:i) 237 if (buflen.ge.512)call wrterm(buffer, buflen) 238 10 continue 239 return 240 end 1 subroutine opterm 2c Open the terminal for read and write. 3 integer chan,flag,lib$get_ef,sys$assign 4 common/termio/chan, flag 5 if(.not .lib$get_ef(flag))stop’Error in OPTERN’ 6 if(.not.sys$assign(‘sys$output’,chan,,))stop’Error in OPTERM' 7 return 8 end 9 10 subroutine clterm lic Close the terminal for read and write. 12 integer chan,flag,lib$free_ef ,sys¢dassgn 13 common/termio/chan, flag 14 if(.not .lib$free_ef(flag))stop’Error in CLTERH’ 15 if(.not .sys$dassgn(%val(chan)))stop’Error in CLTERM’ 16 return 17 end 18 19 subroutine rdterm(buffer, buflen) 20 ¢ Purge terminal input and read the next buflen bytes into buffer 21 ¢ with no echo and no interpretation of control characters. 22 parameter (io$_readvbik=49 , io$m_noecho=64 , io$m_purge=2048) 23 parameter (iocode=io$_readvblk.or.io$m_noecho.or.io$m_purge) 24 character buffer (*) 25 integer buflen,chan ,flag,trmblk(2) ,sys$qiow 26 common/termio/chan,flag 27 data trmblk/0,0/ 28 if(.not .sys$qiow(%val (flag) ,Aval (chan) ,4val (iocode) 29 + ,,,,4ref (buffer) , {val (buflen), ,%ref(trmblk),,)) 30 + stop’Error in RDTERM’ 3i return 32 33 34 35 36 37 38 39 40 41 42 43 44 48 46 47 48 19 20 word wv vw wo ee CO MON DA 32 ~98~ end subroutine wrterm(buffer, buflen) Write buflen bytes to the terminal from buffer with no interpretation of control characters; reset buflen to 0. parameter (io$_writevblk=48, io$m_noformat=256) parameter (iocode=io$_writevblk. or .io$m_noformat) character buffer (+) integer buflen,chan,flag,sys$qiow common/termio/chan,flag if (buflen.gt.0)then if(.not .sys$qiow(4val (flag) ,4val (chan) ,%val(iocode) ,,.,avet (buffer) ,jAval(buflen),,,,))stop’Error in WRTERM’ buf len=0 endif return end #include <sgtty.h> opterm_() /* Open the terminal for read and write. #/ {return;} clterm_() /* Close the terminal for read and write. «/ {return; } rdterm_(buffer , buflen) /* Flush terminal input and read the next buflen bytes into buffer with no echo. */ char buffer[]; long *buflen; { int n; struct sgttyb stateO,statel; if (ioct1(0,TIOCGETP,&stateO) == -1) printf("error in rdterm"); statel=state0; state! .sg_flags=stateO sg flags+CBREAK-ECHO; if (ioctl(O,TIOCSETP,&state1) == -1) printf("error in rdterm"); for ( n=0; n<*buflen; nt++ ) if (read(0,buffer+n,1) != 1) printf("error in rdterm"); if (ioct1(0,TIOCSETP,&stateO) == -1) printf("error in rdterm"); return; } wrterm_(buffer, buflen) /* Write buflen bytes to the terminal from buffer; reset buflen to 0. «/ char buffer[];long *buflen; { if (*buflen == 0 ) return; if (write(1,buffer,(int)*buflen) != (int) buflen) printf ("error in wrterm"); *buflen = 0; return; } long iand_(i,j) /* integer valued function to return logical AND of i and j */ /* i and j are assumed to be 4 byte integers <«/ long *i,*j; ¢ return(*i&*j); } long ior_(i, j) / integer valued function to return logical OR of i and j */ /* i and j are assumed to be 4 byte integers +/ ~99— 43 long +i,*j; ¢f 44 return(*i]*j); 45 } 46 long ishft_(i,nshift) 47 /* shift the 4 byte integer i by nshift bits «/ 48 /* aif nshift is negative, right shift end off zero fill +/ 49 /* if nshift is positive, left shift end around +/ = /* the routine behaves properly if magnitude of nshift > 32 +/ 5 51 long #i,*nshift; { 52 long jshift,nbits,; 53 if (nshitt<o) f 54 nbits = ('nshift < -32 ? 32 : -inshift); 55 jshift = (*i>>nbits) & (017777777777>>(nbits-1)); 5G } 57 else { 58 nbits = *nshift % 32; 59 jshift = (+i<<nbits) | ( (#i>>(32-nbits)) re) & ("(037777777777<<nbits)) ); Gi } 62 return(jshift); 63 } A1.5.3 The Composite Grid Data File CMPGRD produces an unformatted data file that contains all the information about a composite grid that is needed in order to discretise PDE boundary-value problems on it. This file may be read by a program written to discretise a PDE boundary-value problem. Its contents are explained below. One record of four words: xa, xb, ya, yb. The rectangle ra < x < 2b, ya < y < yb bounds the region D. This information is useful for plotting functions on the grid. One record of one word: ng. ng is the number of component grids. One record of 2ng words: (m,, n,z, for k = 1,...,ng). m, and n, are the number of gridlines in the r and s directions, respec- tively, on grid G,. One record of S°/2, myn, words: (((krij,, for i =1,...,m,), for j =1,...,n4), for k= 1,...,ng). kr,;,, indicates whether the point (7,7) on grid G, is a discretisation point, an interpolation point or neither. If kr,;, > 0 then the point is a discretisation point; if kr;;, <0 then the point is an interpolation point; ~100- if kr;;, = 0 then the point is not used. The discretisation and inter- polation points are numbered consecutively by |kr;,,| as first 1, then 9 and then k increase. The exception to this is that corresponding points on corresponding branch cuts are given the same value of kr,,,. In this way, when a PDE boundary-value problem is discretised on the com- posite grid, the discrete equation that applies at each discretisation and interpolation point is given a unique number, with the exception that the identical equations that apply at corresponding points on branch cuts share the same equation number. One record of 5ng words: (bc,, (tbe;,, for 7 = 1,...,4), for k = 1,...,ng). bc, and tbc;, indicate the boundary conditions on each side of each grid G,. bc, is an integer which indicates the directions, if any, in which grid G, is periodic. In particular, bc, = 0 if grid G, is not periodic; be, = 1 if grid G, is periodic in the r direction only; bc, = 2 if grid G, is periodic in the s direction only; and bc, = 3 if grid G, is periodic in both directions. If ¢bc;, #4 0 then |zbc,;| is the number of the curve that side 7 of grid G, is mapped to. ibc;, > 0 indicates that side 7 corresponds to a boundary curve, 7bc;, < 0 indicates that side 7 corresponds to a branch cut, and zbc;, = 0 indicates that side 7 is only interpolated from other grids. One record of ng words: (nb,, for k = 1,...,ng). nb, is the number of interpolated points on grid G,. One record of S772, 8nb, words: ((¢biz, 7bin, kbi,, for ¢ = 1,...,nb,), for k = 1,...,7g). (1b;.,7;.) form a list of the interpolated points on grid Gx, and kb,, indicates the grid from which each point is interpolated. Each interpo- lated point of grid G, appears in the list as (7b;,,76;,) for some 7. The point (2b;,,76:.) is interpolated from grid Gyp,,. One record of S072, 2nb, words: ((rbj,, sbi, for 1 = 1,...,nb,), for k = 1,...,ng). (rb;,, Sb;,) are the coordinates of the interpolated point (2b;,, 76;,) on the grid Gx, from which it is interpolated. The points used for interpola- tion and their interpolation weights can be computed from (rby,, sbix). One record of S°y%, 2myzn, words: (((xijz,¥cjk, for 7 = 1,...,m,), for 7 = 1,...,%%), fork = 1,...,ng). (XijksYijk) are the coordinates of the point (¢,7) on grid G, under the transformation d,. Fork =1,...,ng,forj =1,...,n,, fora = 1,...,m,, one record of four words: LT igks Yrigks USijky YSijzk- LVijky YTizks LSijks YSijk are the components of the Jacobian derivative Lg —101- ‘ ] at the point (7,7) on grid Gx. —102- APPENDIX 2 Multigrid on Composite Meshes* W.D. Henshaw and G. Chesshire Department of Applied Mathematics 217-50 California Institute of Technology Pasadena CA 91125 Abstract. The multigrid method is applied to numerical solution of elliptic equations on general composite overlapping meshes. Computational results show that good convergence rates are obtained. Key words. Composite meshes, overlapping grids, multigrid. A2.1 Introduction We describe the application of the multigrid method to the solution of elliptic partial differential equations (PDEs) on two dimensional regions which have been discretized using composite overlapping grids. A general purpose code, CGMG, has been developed (in FORTRAN) which can solve problems on composite meshes created by the grid construction program CMPGRD, Chesshire [6]. With CMPGRD the user may create a composite grid containing any number of component grids at any number of multigrid levels. CGMG can then be used to solve elliptic PDE boundary value problems. * Submitted to SIAM J. Sci. Statist. Comput. -103- A composite overlapping grid consists of a number of simpler component grids. These component grids cover a region and overlap where they meet. Functions defined on the composite mesh are matched by interpolation at the overlapping grid boundaries. The problem of generating grids for regions of complicated geometry can be difficult, especially for those grid generation algorithms which attempt to fit a single global grid. With a composite overlapping grid, however, the component grids can be generated almost independently of each other. Each component grid can be stretched and refined with little effect on the other component grids. The numerical solution of PDEs on such grids has been examined by, among others, Starius [18], [19], Reyna [16], Kreiss [14], Atta [1], Benek et al. [2], Henshaw [12] and Berger [3]. The multigrid method is a fast iterative method for the solution of elliptic problems. Multigrid utilizes a sequence of grids of varying degrees of coarseness to accelerate the convergence of the solution on the finest grid. The basic principle rests on the fact that it is possible to obtain iterative procedures (smoothers) for which the high frequency components of the solution converge rapidly. This means that after a few smoothing iterations, the part of the solution yet to converge is smooth and hence can be accurately solved for on a coarser grid. In the rest of this paper we outline the implementation of the multigrid algo- rithm on fairly general composite meshes. Multigrid on a model composite mesh is described in Sttiben and Trottenberg |20] and, for more general two component meshes, in Henshaw [12]. Numerical results are presented for some test cases. These results show that the good convergence rates expected from multigrid can be ob- tained on composite meshes. Acknowledgements. Computations were performed on the Fluid Dynamics Vax 11/750 at the California Institute of Technology. Support for this work came from the National Science Foundation under contract DMS-8312264 and the Office —104- of Naval Research under contract NO0014-83-K-0422. A2.2 Description. The multigrid code CGMG was written to solve linear, variable coefficient el- liptic PDEs of the form (14) L O*u n Oru n Oru n Cu n Ou 1 f ing a U I= Crp 575 Cry pa a na a cu >= dn? “Ardy | “Way” ar) Vay 1 oe m O (16) Bu:= 0, _ bys +bu=g onan. Periodic boundary conditions are also allowed. Currently a second order difference approximation has been implemented. Let us first describe how the above PDE boundary value problem is discretized on a composite mesh, without any references to multigrid. The grid construction program CMPGRD can be used to generate a composite overlapping mesh for the region 9. This mesh is composed of one or more component meshes. See, for example, the composite meshes of figures A2.1 and A2.2. Each component mesh is logically rectangular, although some points in overlap regions are discarded. A point on a component grid will be one of the following four types : (i) An intertor point where the PDE (1a) should be applied. (ii) A boundary point, which corresponds to 0M, where the boundary condition (1b) holds. (iii) A point where the solution is matched, by interpolation to the solution on other component grids. We sometimes call the set of all such interpolation points the interpolation boundary. (iv) A point in a region of overlap which is not used. There are a number of possible ways to generate the discrete approximations to the PDE (1a) the boundary conditions (1b) and the interpolation equations. The ap- proach we take is a mapping method and proceeds as follows. Each component grid —105- is simple enough so that it can be mapped smoothly to a unit square (coordinates (r,s)). The PDE can be written in these (r,s) coordinates. L O*u i Ou Oru 4 Ou in Ou teu=f (= ¢,,—— Cra +e, teu = 4 Ore Or? rs Oras “Os? Or “OS Ou Ou = b— +b, +b,u = Bu=6 Or + "OS + eu = 9 The coefficients in these equations depend on derivatives of the mapping. These derivatives are supplied as output from the composite grid construction program CMPGRD. The solution of a PDE on a composite mesh can then be considered as the solution of PDE’s on a sequence of unit squares. The solutions on the squares are coupled to each other through the interpolation boundaries. Let the discrete solution at point (¢,7) on component grid k be denoted by v,(i,j). Then, for each k we obtain discrete approximations to the elliptic equation and the boundary conditions of the form (2a) Livy = fre (26) Biv, = Gk. A point (7,7,k) on an interpolation boundary is interpolated from some other com- ponent grid k'. CMPGRD supplies the position (r',s') of the point (7,7,k) on grid k'. Hence, standard interpolation formulae for rectangular grids can be applied : ag! These equations (2a), (2b) and (2c) will be written as a single linear system (3) Av =f. The vector v of all unknowns will be called a composite mesh function. A is an example of a composite mesh operator, mapping one composite mesh function to another composite mesh function. ~106- The mesh equations (3) can be solved in any number of ways; using multigrid is but one possibility. If there are not too many equations, the system can be solved directly. Sparse matrix routines |11] have been used for this purpose. For large systems iterative methods become attractive. Many standard iterative methods are applicable. The matrix is not symmetric, however, so that some schemes do not apply. In general it is best to try and solve all the equations simultaneously, rather than iterating for too long on one component grid. In Stiiben and Trottenberg [20] it was shown that an iteration based on the Schwartz alternating procedure (where the equations on each component grid are solved exactly before updating interpolation boundaries) is slower and sensitive to the amount of overlap. The general principle to follow seems to be to iterate in such a way that, at any time, all equations have converged to about the same degree. A2.2.1 Multigrid. The multigrid algorithm can be applied to the solution of the mesh equations (3). Discussions of the multigrid method in general can be found, for example, in Brandt [5], or Stiiben and Trottenberg [20]. Once a composite mesh has been constructed using CMPGRD it is a simple matter to have CMPGRD generate the sequence of coarser composite meshes which are used for the multigrid algorithm. Figures A2.1 and A2.2 show some composite meshes which have been generated for multigrid. Denote the finest composite mesh, level 1, by M? and successively coarser meshes by M!,] = 2,3,.... Note that for simplicity the composite meshes at the different levels all have the same number of component meshes. The elliptic PDE boundary value problem can be discretized on each of the composite meshes M!. Let v! denote the composite mesh function for level /. Then at each level ! there will be mesh equations of the form (3) (4) Alv! = f!. -107- Now let us outline the multigrid algorithm as it applies to composite meshes. Let v'(p) be the pth iterate in the solution of the mesh equations on the finest mesh M!. Multigrid Algorithm while not converged do smooth Vv, times vi (+) — (S*)"v"(p) compute the defect and transfer to the coarser grid f? — R2(f1 — Alv}(s)) “solve” the defect equation on the coarser grid yee (A2) "£2 correct the fine grid solution from coarse grid solution vi(x*) — vi(x) + PPoty? smooth v2 times vi(p + 1) — (81) vi (+s) end while This is the basic defect correction scheme. The smoothing operator S!, the restriction operator R!~? and the prolongation operator P?~! will be described in the context of composite grids. The defect equation need only be “solved” approxi- mately. This approximate solution can be obtained by multigrid, in which case the algorithm becomes recursive. At the coarsest level the equations are usually solved directly. Smoothers. The smoothing operator S for composite grids consists of smooth- ing each component grid and updating the interpolation boundaries. A component grid may be smoothed with any of the standard smoothers that exist. The program CGMG allows the user to choose from a number of possibilities including Gauss— —108-— Seidel, Red-Black—Gauss-Seidel, Zebra line smoothers and alternating Zebra line smoothers. Each component grid can have a different smoother; the smoother can be tailored to the grid. For example, a particular component grid might be stretched to resolve a boundary layer, in which case one can use a line smoother in the appropriate direction. There is some freedom as to the order of smoothing and interpolation. One pos- sibility is to smooth all component grids before updating the interpolation bound- aries : Smooth first component grid Smooth second component grid S! = Smooth last component grid Interpolate From experience, however, it seems that a good procedure involves interpolating after each component grid is smoothed : Smooth first component grid Interpolate Ss! — Smooth second component grid . Interpolate The latter composite smoother requires more interpolations but this extra work is usually small compared to the smoothing operations. Restriction Operators (Fine to Coarse Grid Transfer). The defect com- puted on a given level is transferred to the next coarsest level by the restriction operator R'~?. (The superscript 1 — 2 indicates that this operator maps mesh functions on M' to mesh functions on M’.) A typical restriction operator de- termines the value at points on the coarse grid as soine weighted average of the surrounding points on the fine grid. We use the so called full weighting restriction. The defects in the boundary equations are averaged separately from the defects in the interior equations; boundary defects are averaged along the boundary line. Since the final stage of the smoothing operation involves an interpolation, the de- —109- fects in the interpolation equations are all zero. Hence, no defect need be transfered at these points. Prolongation Operators (Coarse to Fine Grid Transfer). The prolon- gation operator P?~} maps the coarse grid solution to the fine grid. This mapping usually takes the form of an interpolation. We use second order interpolation. In- terpolation boundary values could be corrected as well. However, once all other values have been corrected, the interpolation equations can be solved to update the interpolation boundary. It turns out that less overlap is needed on finer meshes for the latter approach. Choice of Parameters and Cycle. An important part of the multigrid algorithm is the choice of the parameters 1), 2, etc. and the choice of cycle. Cycle is the term used to denote the sequence in which the different levels of grids are traversed. The program dynamically determines the type of cycle and the values for yy and v2 in a manner similar to that described by Brandt [5]. There are two basic principles : (i) Perform smoothing iterations until the smoothing rate (the reduction in residual per iteration) becomes larger than some value 7, where n ~ .6. The smoothing typically starts out very fast but then slows down once the high frequencies have been reduced compared to the low frequencies. (ii) Return to a finer grid once the residual at this level has been reduced by a certain factor 6 with 6 = .1. | Ideally the program should automatically choose the types of component smoothers and the parameters 7 and 6. This could be done, for example, by examining com- ponent grids for stretching and by monitoring the smoothing rates. A program with this level of sophistication has not yet been developed. Currently one must make some educated guesses and then examine the results to determine any changes that might improve convergence rates. -110- Remarks on Grid Construction. For practical reasons the amount of overlap between component grids is kept fairly small. This reduces the number of compu- tational points. With sufficient overlap, points on an interpolation boundary can be explicitly interpolated from non-interpolation points on other component grids. However, as the amount of overlap decreases the interpolation points may become coupled through the interpolation equations. In this case a system of equations must be solved to obtain values on the interpolation boundaries in terms of other values. Of course, as the overlap goes to zero these equations may become singular. The behaviour of the numerical solution to a model elliptic problem, as a function of the amount of overlap, was considered in Henshaw [12]. Suppose that the amount of overlap goes to zero as the grid is refined. To maintain accuracy, the order of accuracy of the interpolation formulae must be greater than the order of accuracy of the interior formulae. In the results presented here we use second order accurate approximations to the elliptic equations and require the amount of overlap to be greater than one half a grid line. In this case the interpolation formulae should be third order accurate. To make the restriction and prologation operators simple, the composite meshes at different multigrid levels are strongly related. For example, interior points on a coarse grid coincide with interior points on finer grids. In Henshaw [12] the coarse grid was generated first and all finer grids were generated by simply doubling the number of grid lines. This gave reasonble results for composite grids with 2 component grids. However, this method does not immediately generalize to more component grids, and it tends to lead to more overlap on the finer grids than is really necessary. Hence, the grid construction program CMPGRD was designed explicitly to meet the requirements of the multigrid routine [6]. Remarks on Programming. A computer code for solving problems on com- posite meshes is somewhat more complicated than a code written for a single grid. —111- One of the major difficulties arises in the handling of all the data that descibe the composite grid. We have worked out a scheme for storing the composite grid data in an efficient and flexible manner. Physically, all the data is stored on a single array. Logically, the data are stored in a “directory-tree” fashion. Variables can be stored or accessed by name using utility routines. Arrays can be stored with a minimum of wasted space. With this storage structure, the composite grid data can be easily passed to subroutines. A2.3 Numerical Results In this section we present results from the multigrid solver CGMG. The com- posite grids used in the examples are shown in figures A2.1 and A2.2. Grid 1 (figure A2.1) has 3 levels and Grid 2 (figure A2.2) has 2 levels. The interpolation points are marked with small circles. In both examples we solve Ou d?u 4 —— + —— = (4) gat and f = -2n’ cos rzcos Ty. The boundary conditions are chosen so that the true solution is Utrue = COS TXCOS TY. Test 1 The Poisson equation (4) is solved on Grid 1. This grid might to be used to study flow around an obstacle in a varying channel. Periodic boundary conditions are applied at the left and right ends of the channel. Dirichlet boundary conditions, U = Utrue, are given at all other boundaries. For smoothers we chose Red—Black on the rectangular grid and line~Zebra smoothers for the 3 curvilinear grids. The Zebra smoothers are on lines normal to the boundaries, since this is the direction in which the grids are stretched. Smoothing was performed until the smoothing rate became greater than n = .6. The residual at level 2 had to be decreased by a factor of 6 = .01 before the program would return to level 1. These parameters resulted in —112- Grid 1 Grid 2 Iteration | r(p)/r(p-1) | ECR r(p)/r(p-1) ECR p=1 019 7 ~.091 73 p=2 065 72 033 63 p=3 1055 71 11 75 p=4 385 .89 14 17 Table A2.1: Convergence Rates a “W” cycle; the levels were traversed in the order (1,2,3,2,3,2,1). A value of 6 = .1 resulted in a “V” cycle, (1,2,3,2,1), which gave slightly inferior results. Test 2 Equation (4) is solved on Grid 2 with Dirichlet boundary conditions at all boundaries. Grid 2 is meant to represent an airfoil with a flap. Alternating Zebra smoothers were used on all component grids, since they all are stretched in both the r and s directions. A value of 7 = .36 = (.6)” was chosen since the alternating smoothers are comprised of two sweeps. Convergence results for the two tests are summarized in table A2.1. r(p) is the residual on the finest composite mesh after the pth multigrid iteration. r(p) = ||f' — Atv"(p)|| Define WU(p) to be the number of work units used for the pth iteration. A work unit is the amount of work to perform one iteration of SOR on the composite mesh. The effective convergence rate, ECR, is defined as ner(n) = (EY 9 win Theoretically the ECR for multigrid should be independent of the grid spacing h as h — 0. In contrast the (effective) convergence rate for many other standard iterative schemes deteriorates as h — 0. For optimal SOR, ECR ~ 1 — cyh, and for Gauss-Seidel, ECR = 1 — cyh?. In table A2.1 the ECR increases on the fourth iteration for Grid 1, since the solution has almost converged to single precision accuracy. -113- In tables A2.2 and A2.3 the errors in the solution to the elliptic problem are given. The problem was solved on the fine mesh (level 1) and also on each of the coarser meshes. The error on level J, e(/), is the maximum difference between the calculated and true solution. The number of points on the component grids in the r and s directions is given by n, and ng, respectively. Not all the points are used, since there is overlap between the grids. Level 1 of Grid 1 has 8431 computational points while there are 7431 computational points on the finest level of Grid 2. Contour plots of the calculated solutions are shown in figures A2.1 and A2.2. l k Ny Ns e(l) e(l) /e(1) 1 1 93 77 1.5x1073 1. 2 81 21 1.1x107° 1. 3 81 21 1.1x1073 1. 4 93 17 1.7x1073 1. 2 1 47 39 7.4x1073 5. 2 41 11 5.4x107° 4.8 3 4] 11 5.4x1078 4.8 4 47 9 7.4x107 4.4 3 1 24 20 3.7x107? 25. 2 21 6 2.6x107? 23. 3 21 6 2.6x107? 23. 4 24 5 3.3x107? 20. Table A2.2: Errors in Solution for Grid 1 l k ny Ns e(l) e(1) /e(1) 1 1 101 61 1.8x10°? 1 2 121 13 2.0x107? 1 3 111 13 5.7x107% 1. 2 1 51 31 9.8x 107? 5.4 2 61 7: 8.7x10°7 4.4 3 56 7 2.7x107? 4.7 Table A2.3: Errors in Solution for Grid 2 -114- Sones Zee A ERE AAA Py rt to , tt 4 AH Fatt = A ion H HHH H bit vf Hi i itt BUY tt tt ry ry HH} jaseea: 1 ed at . 3 ftgeeeabiey Hi <4 H elpassssess! i Hitt ioigietgyss HH : L-4 Z i ii ta bv ~ aoe On ae eeSaue Ui l qj] ee Hg ERS = ISA pd] A t { A LA L es = Le] BZ ru Fig. A2.1. Grid 1: Composite meshes for multigrid levels 1, 2 and 3 and a contour plot of the calculated solution. A ee 4G ES plot of the calculated solution. [1] [12] [13] [14] [15] —116- References im. H. Atta and J. Vadyak, A Grid Overlapping Scheme for Flowfield Com- putations about Multicomponent Configurations, AIAA J., 21 No. 9, (1983), pp. 1271-1277. J. A. Benek, J. L. Steger and F. C. Dougherty, A Flexible Grid Embedding Technique with Application to the Euler Equations, AIAA Paper No. 83-1944, Danvers, Mass., 1983. M. Berger, On Conservation at Grid Interfaces, ICASE Report No. 84-43, 1984. J. U. Brackbill and J. S. Saltzman, Adaptive Zoning for Singular Problems in Two Dimensions, J. Comp. Phys., 46 (1982), pp. 342-368. A. Brandt, Guide to Multigrid Development, in Multigrid Methods, W. Hack- busch and U. Trottenberg, eds., Springer-Verlag, Berlin (1982), pp. 220-312. G. Chesshire, CMPGRD: A Composite Grid Construction Program, To appear. S. A. Coons, Surfaces for Computer-Aided Design of Space Forms, Dept. of Mech. Eng., Mass. Inst. of Tech. , Cambridge, 1964. P. R. Eiseman, A Multi-Surface Method of Coordinate Generation, J. Comp. Phys., 33 (1979), pp. 119-151. P. R. Eiseman, Grid Generation for Fluid Mechanics Computations, Ann. Rev. Fluid Mech., 17 (1985), pp. 487-522. 5. K. Godunov and G. P. Prokipov, The Use of Moving Meshes in Gas-Dynamical Computations, USSR. Comput. Math. and Math. Phys., 12 (1972), pp. 299-319. . C, Eisenstat, M. C. Gursky., M. H. Schultz and A. H. Sherman, Yale Sparse Matrix Package, Research Reports 112 and 114, Yale University, Dept. of Com- puter Science, May 1977. W. D. Henshaw, Part I: The Numerical Solution of Hyperbolic Systems of Con- servation Laws; Part II: Composite Overlapping Grid Techniques, Ph.D. thesis, Department of Applied Mathematics, California Institute of Technology, 1985. W. D. Henshaw and G. Chesshire, Multigrid on Composite Meshes, submitted to SIAM J. Sci. Stat. Comput. B. Kreiss, Construction of a Curvilinear Grid, SIAM J. Sci. Stat. Comput., 4 No. 2 (1983), pp. 270-279. C. W. Mastin and J. F. Thompson, Elliptic Systems and Numerical Transfor- mations, J. Math. Anal. App., 62 (1978), pp. 52-62. [16} [17] 18} 19] [20] [21] ~117- L. G. M. Reyna, Ph.D. thesis, Department of Applied Mathematics, California Institute of Technology, 1982. G. Starius, Construction of Orthogonal Curvilinear Meshes by Solving Initial Value Problems, Numer. Math., 28 (1977), pp. 25-48. G. Starius, Composite Mesh Difference Methods for Elliptic Boundary Value Problems, Numer. Math., 28 (1977), pp. 243-258. G. Starius, On Composite Mesh Difference Methods for Hyperbolic Difference Equations, Numer. Math., 35 (1980), pp. 241-255. K. Stiiben and U. Trottenberg, Multigrid Methods: Fundamental Algorithms, Model Problem Analysis and Applications, in Multigrid Methods, W. Hack- busch and U. Trottenberg, eds., Springer-Verlag, Berlin (1982), pp. 1-176. J. F. Thompson, F. C.Thames and C. W. Mastin, Automatic Numerical Gener- ation of Body-Fitted Curvilinear Coordinate System for Field Containing Any Number of Arbitrary Two-Dimensional Bodies, J. Comp. Phys., 15 (1974), pp. 299-319. J. F. Thompson, ed., Numerical Grid Generation, Elsevier, New York, 1982. A. M. Winslow, Numerical Solution of the Quasilinear Poisson Equation in a Nonuniform Triangle Mesh, J. Comp. Phys., 2 (1967), pp. 149-172. P. Woodward and P. Colella, The Numerical Simulation of Two-Dimensional Fluid Flow with Srong Shocks, J. Comp. Phys., 54 (1984), pp. 115-173.