Numerical Grid Generation
Foundations and Applications
By: Joe E. Thompson, Z.U.A. Warsi and C. Wayne Mastin
IV. NUMERICAL IMPLEMENTATION
1. Transformed Eqations
In order to make use of a general boundaryconforming curvilinear coordinate system in the solution of partial differential equations, or of conservation equations in integral form, the equations must first be transformed to the curvilinear coordinates. Such a transformation is accomplished by means of the relations developed in the previous chapter and produces a problem for which the independent variables are time and the curvilinear coordinates. The resulting equations are of the same type as the original ones, but are more complicated in that they contain more terms and variable coefficients. The domain, on the other hand, is greatly simplified since it is transformed to a fixed rectangular region regardless of its shape and movement in physical space. This facilitates the imposition of boundary conditions and is the primary feature which makes grid generation such a valuable and important tool in the numerical solution of partial differential equations on arbitrary domains.
A numerical solution of the transformed problem can be obtained using standard techniques once the problem is discretized. Since the domain is stationary and rectangular, and since the increments of the curvilinear coordinates are arbitrary, the computation can always be done on a fixed uniform square grid. Spatial derivatives at nearly all field points in the transformed domain can therefore be represented by conventional finitedifference or finitevolume expressions, as discussed in the next section. In fact, the transformed problem has the appearance of a problem on a uniform cartesian grid and thus may be treated as such both in the formation of the difference equations and in the solution thereof.
The specific form of the transformed equations to be solved depends, of course, on which of the realtions in Chapter III are used, i.e., conservative or not. As an example, consider the generic convectiondiffusion equation

(1) 
Equations (III123), (III42), and Eq. (III43) may be used to transform the convective terms, the gradient, and the second divergence, respectively, and thereby yield the conservative form:

(2) 
where now the time derivative is understood to be at a fixed point in the transformed region, and the contravariant velocity components (relative to the moving grid) are given by Eq. (III122). Eq. (2) can also be written in the form

(3) 
which clearly shows the conservative form. It is the product A rather than the function A itself, which is conserved in this form. The derivative inside the j summation can be expanded and Eq. (III40) invoked to obtain the simplified form:

(4) 
which is still conservative in regard to the derivatives.
These conservative forms are in the commonly used form
where the solution vector is , the "flux" vectors ^{i} are given by the brackets in (3) and (4), and the source vector is =.
The flux vectors ^{i} contain metric derivatives and depend on time and the curvilinear coordinates through these metric elements, as well as through the solution vector and this must be taken into account in the construction of factored solution methods. A general formulation of split solution methods (encompassing both time splitting, e.g., approximate factorization, and spatial splitting, e.g., MacCormack method) in the curvilinear coordinates can, however, be formulated.
The nonconservative form of Eq. (1) follows using Eq. (III125) for the convective terms, Eq. (III65) for the gradient, and Eq. (III67) for the divergence . The resulting equation may be written

(5) 
since Eq. (III70) gives

(6) 
(The last summation in Eq. (5) is just , which vanishes for incompressible flow.) Comparison of Eq. (5) with the original equation, written in the form

(7) 
demonstrates that the equation has been complicated by the transformation only in the sense that the coefficient u_{i} has been replaced by the coefficient U^{i}+(^{2} ^{i}), and the Kroniker delta in the double summation has been replaced by g^{ij}, thus expanding that summation from three terms to nine terms, and through the insertion of variable coefficients in the last summation. This exemplifies the fact that the use of the general curvilinear coordinate system does not introduce any significant complications into the form of the partial differential equations to be solved. When it is conthat the transformed equation (5) is to be solved on a fixed rectangular field with a uniform square grid, while the original equation (7) would have to be solved on a fiels with moving curved boundaries, the advantages of using the curvilinear system are clear.
These advantages are further evidenced by consideration of boundary conditions. In general, boundary conditions for the example being treated would be of the form

(8) 
where is the unit normal to the boundary and , , and are specified. From Eq. (III79) these conditions transform to

(9) 
for a boundary on which ^{i} is constant. For comparison, the original boundary conditions (8) can be written in the form

(10) 
The transformed boundary conditions thus have the same form as the original conditions, but with the coefficient n_{j} replaced by g^{ij}/. The important simplification is the fact that the boundary to which the transformed conditions are applied is fixed and flat (coincident with a curvilinear coordinate surface). This permits a discrete representation of the derivatives A_{j} along the transformed boundary without the need for interpolation. By contrast, the derivatives A_{xj} in the original conditions cannot be discretized along the physical boundary withoutinterpolation since the boundary is curved and may be in motion.
This discussion of a generic convectiondiffusion equation and associated boundary conditions should serve to allow specific physcial equations to be transformed. References to application of these equations are gven in the surveys Ref. [1] and [5]. Several examples also appear in Ref. [2].
2. Discrete Representation of Derivatives
Approximate values of the spatial derivatives of a function which appear in the transformed equations may be found at a given point in terms of the function's value at that point and at neighboring points. As noted earlier, with the problem in the transformed space, only uniform square grids need be considered, hence the standard forms for difference representation of derivatives may be used. For example, in two dimensions the first, second, and mixed partials with respect to the curvilinear coordinates and are ordinarily represented at an interior point (i,j) by finite differences or finitevolume expressions which contain function values at no more than the nine points shown below.
This centered, ninepoint "computational molecule" is usually preferred because of the associated difference representations which are symmetrypreserving and secondorder accurate. Examples of finitedifference approximations of this type are:

(11a) 

(11b) 

(12a) 

(12b) 

(13) 
Other secondorder approximations of the mixed partial (f_{})_{ij} which use the ninepoint molecule are:

(14) 
and

(15) 
It is clear that at boundary points, where at most first partials must be represented, the computational molecule cannot be centered relative to the direction of the coordinate ^{} which is constant on the boundary (see diagram below).
There a onesided difference must be used to approximate f_{}. The secondorder formula appropriate for the boundary point indicated above is
Any standard text on the subject of finitedifference methods will provide formulas of alternate order and/or based on other computational molecules.
A finitevolume approach uses function values at gridcell centers and approximates derivatives at a cell center by line (surface in 3D) integrals about the cell boundary which are equivalent to averages over the cell. In particular, the identity

(16) 
is used, where V is the volume of D. Thus, if a function is assumed constant along a gridcell face, it is a simple matter to evaluate the line integral in (16) when D is a grid cell in transformed space. In terms of the twodimensional grid:
this approach gives

(17a) 

(17b) 
With an edge value approximated as the average of the center values of the two cells sharing that face, e.g.

(18) 
the values given by (17) are equivalent to ordinary central differences (cf. Eq. (11)) and hence are secondorder accurate. The first partials of f may also be assumed constant along each cell edge in order to derive from (16) the following approximations of second and mixed partials at a cell center:
Now, however, the averaging scheme in (18) cannot be used to approximate edge values of the derivatives without going outside the ninepoint computational molecules shown above. Instead, a secondorder accurate representation can be obtained on the ninepoint molecule using a forward (backward) assignment for the center value of a function and a backward (forward) assignment for the first partial on a given side. There are four possible schemes of this type. One uses

(20) 
to evaluate f(,) at all cell centers according to (17), and then uses

(21) 
to evaluate the second and mixed partials given in (19). This method is equivalent to a finitedifference scheme which approximates first partials by backward differences of the function, and then approximates second and mixed partials by forward differences of the first partials. Consequently, the second derivatives which result are equal to those given in Eq. (12), while the resulting representations of the two mixed partialsare unequal and only firstorder accurate. If the two mixed partials are averaged, however, the secondorder expression (15) is recovered. This is also true of the reverse scheme:

(22) 

(23) 
Expressions (12) and (14) are similarly recovered from the other two possibilities (Eq. 20a, 21a, 22b, and 23b or Eq. 20b, 21b, 22a, and 23a). The symmetrypreserving form (13) can be recovered by averaging the averaged mixed partial obtained in one of the first two schemes mentioned and that obtained in one of the remaining two.
The manner in which boundary conditions are treated in a finite volume approach depends on the type of conditions imposed. When Dirichlet conditions are prescribed, it is advantageous to treat the boundary as the center line (plane in threedimensions) of a row of cells straddling the boundary. The centersof these cells then fall on the physical boundary where the function values are known. When Neumann or mixed conditions are given, however, the boundary is best treated as coincident with cell faces.
Suppose, for example, that boundary condition (9) is to be imposed at the cell edge =j1/2 indicated below.
The edge value of f_{i,j1/2} cannot be approximated by the usual averaging scheme (illustrated by Eq. (18)) since there is no cell center at =j1. It can, however, be found in terms of neighboring cellcentered function values by using boundary condition (9) in connection with the forward/backward scheme used to approximate second derivatives at the cell centers.
Considering the scheme represented by Eq. (20) and (21), the values of f along the cell edges shown above are:
It follows from Eq. (17) that the first partials of f at the cell center are
Eq. (21a,b) then give f_{} and f_{} along the cell edges enclosing (i,j) in terms of f_{i1,j}, f_{i1,j+1}, f_{i,j}, f_{i,j+1}, f_{i+1,j}, x_{i} and x_{i+1}. In particular,
Substitution of these expressions into boundary condition (9) then determines the edge value x_{i} as
In this way, f, and hence f_{} and f_{}, are found on all boundarycell edges in terms of cellcentered values of f.
The finitedifference and finitevolume techniques described thus far are appropriate for representing all derivatives with respect to the curvilinear coordinates, even those appearing in the metric quantities. In fact, as it is shown later in this chapter and in chapter V, the metric quantities should be represented numerically even when analytical expressions are available. One might have, for example,

(24) 
3. Special Points
Many of the expressions given in the previous section break down at socalled "special points" in the field where special attention is required in the approximation of derivatives. These points commonly arise when geometrically complicated physical domains are involved. As indicated in Chapter II, special points can occur on the domain boundary and on interfaces between subregions of a composite curvilinear coordinate system. They may be recognized in physical space as those interior points having a nonstandard number of immediate neighbors or, equivalently, those points which are vertices, or the center, of a cell with either a nonstandard number of faces or a vertex shared by a nonstandard number of other cells. (In two dimensional domains, ordinary interior points have eight immediate neighbors [refer to figure on p.141]; standard twodimensional interior grid cells have four sides and share each vertex with three other cells [see diagram on p. 143].) Boundary points are not special unless they are vertexcentered and have a nonstandard number of immediate neighbors (other than five in two dimensions see diagram on p. 142 for an ordinary boundary point) and then are special only when their assocciated boundaryconditions contain spatial derivatives. Some examples of special cellcentered points and special vertexcentered points are shown below.
When a finitedifference formulation is used, the usual approach, as described in Section 2, can be followed at a special point P if the transformed equations and difference approximations at that point are rephrased in terms of suitable local coordinates. The local system is chosen so as to orient and label only the surrounding points to be used in the needed difference expressions. Choices appropriate to various special points are listed in Tables 1, 2, and 3.
The difficulties encountered at special points in a finitevolume approach are clearly seen by considering the image in the transformed plane. The first pair of diagrams below, for example, shows that at centers of cells having the usual number of faces but sharing a vertex with a nonstandard number of cells, such difficulties amount to mere bookkeeping complications when only first partials must be approximated. Equations (17) and (18) still apply, but the indices must be defined to correctly relate the cell centers on the two sides of an interface. The following diagrams
also illustrate the breakdown at all special cellcentered points of the previouslydescribed finitevolume schemes for approximating second and mixed partial derivatives. This is because the forward/backward orientation of the coordinate system in one segment cannot be consistently followed across the interface adjacent to, or intersecting, the special points. The second pair of diagrams displays the additional complication associated with grid cells having a nonstandard number of edges. Such a cell can occur on an interface between segments of a composite grid which are joined between grid lines. When the segments are transformed to their respective images, the separate pieces of the special grid cell cannot be joined without distorting them. It is thus unclear how to evaluate the volume and the outward normals of that transformed cell in order to use identity (16) in the transformed plane. Consequently, at special points of this type and at all special points where second derivatives must be approximated, the governing equations are best represented locally in the physical plane where such ambiguities do not exist.
Treatment in physical space involves approximation of the original equations by means of identity (16). Thus, for a twodimensional Nsided cell of area A with cartesian centroid P = (p_{1},p_{2}), vertices i=1,2,...,N, and edges s^{i} joining V^{i} and V^{i+1} (V^{N+1}=V^{1}) along which a function f and its first partial derivatives are constant, this approach gives
where the superscripts on f and its derivatives indicate the point or face of evaluation. As in the previous section, an obvious way to approximate f^{si} is to average the center values of the two cells sharing edge s^{i}. This same averaging scheme cannot be repeated to approximate , and , however, without rejecting the recommended strategy of avoiding use of values at points which are not immediate neighbors of the point at which a quantity is being evaluated. Instead, we propose the averaging technique:
where the vertex values are obtained by applying identity (16) to auxiliary cells formed by joining the midpoints of the edges of each cell to the cell center. To make this more precise, let V be a vertex common to Q cells and label the cell faces emanating from V as k_{i} with midpoints
Then if is the center of the cell having edges k^{i} and k^{i+1}, and if
the first partial derivatives of f at V may be approximated by
where A is the area of the 2Qfaced auxiliary cell M^{1}P^{1}M^{2}P^{2}...M^{Q}P^{Q}M^{1} indicated in the following diagram.
This technique is applicable to all grid cell centers; however, it is recommended for use only at points where the methods developed in section 2 break down, since the difference representations associated with those methods are simpler.
4. Metric Identities
When the transformed equations are in conservative form, it is possible for the metric coefficients to introduce spurious source terms into the equations, as has been noted in several works cited in Ref. [1] and as discussed also in Ref. [11] and [12]. This is because the metric coefficients have been included in the operand of the differential operators and if the differencing of these coefficients does not numerically satisfy identities (III40) and (III120), the numerical representations of derivatives of uniform physical quantities are nonvanishing.
For example, if the quantity A is constant, the conservative form for the gradient, Eq. (III42) gives
which is precisely Eq. (III40). Relations (III43)  (III45) similarly reduce to (III40) when is uniform. Therefore, Eq. (III40), or equivalently Eq. (III21), is a metric identity which must be satisfied numerically in order that the conservative expressions for the gradient, divergence, curl, and Laplacian, etc., vanish when the physical variable is uniform. This consideration does not arise with the nonconservative forms since the quantity A is differentiated directly in those expressions.
Another metric identity which must be satisfied numerically arises when the grid is timedependent. This may be seen by considering a generic conservation equation of the form
The conservative relation (III121) transforms this to

(25) 
where now the time derivative is understood to be at a fixed point in the transformed space. If A and are both constants, then Eq. (25) gives
which vanishes according to Eq. (III40). Expansion of the lefthand summation subject to Eq. (III40) then reveals the additional identity to be satisfied:

(26) 
which is just Eq. (III120). This equation, therefore, is that which should be used to numerically determine updated values of the Jacobian, . For if is instead updated directly from the new values of the cartesian coordinates, spurious source terms will appear.
The following example provides a simple illustration of differencing schemes which do, and do not, satisfy the metric identities. The conservative expression for a first derivative in two dimensions is given in Eq. (III96) as

(27) 
which for uniform f reduces to

(28) 
Suppose that f_{x} is to be represented at the center of the cell shown below.
The differencing scheme should satisfy
One possible candidate is the sequence of central differences represented by

(29a) 

(29b) 
The resulting expressions for the mixed partials are
which are indeed equal and thus satisfy identity (28). An alternate choice might be to use central differences for the second differentiation as in Eq. (29a), while approximating the required edge values of the first partials by the average of the values at the adjacent nodes, e.g.
The nodal values are reasonably represented by central differences such as
This scheme cannot possibly satisfy (28), however, since the points used to represent are:
while those needed to evaluate are:
It should be noted that the representations in both of these schemes are consistent and of the same formal order of accuracy. Also, if the metric coefficients at the grid points were evaluated and stored, it would perhaps be natural to follow the second approach, using averages of the metric coefficients at the intermediate points. This, however, is not acceptable since it fails to satisfy the metric identity involved and thus would introduce spurious nonzero gradients in a uniform field.
This example suggests one basic rule that should always be followed: Never average the metric coefficients. Rather, average the coordinate values themselves, if necessary, and then calculate the metric derivatives directly. Alternatively, a coordinate system can be generated with mesh points at all of the halfinteger points, as well as at the integer points used in the physical solution. The metric coefficients can then be evaluated directly by differencing between neighboring points, even at the halfinteger points. For example,
This approach was used in Ref. [13] and problems with the metric identities were thereby eliminated.
It is also possible to construct difference representations which do not involve any averaging and yet still do not satisfy the metric identities; schemes which use unsym metric differences are an example. Fortunately, most reasonable symmetric expressions without averaging do satisfy the identities.
In the representation of the Laplacian using Eq. (III71), ^{2}^{i} should be calculated using Eq. (III74), rather than using derivatives of the metric tensor elements.
Caution is required even when the coordinate transformation is known explicitly. In that case, the metric coefficients can be evaluated analytically, but the metric identities will not in general be satisfied numerically when these coefficients are differenced. This is true even in the simple case of cylindrical coordinates as the following example shows. With
the partials of y(,) are
If the first partial derivative f_{x} is represented as in Eq. (27) and the difference in Eq. (29a) is used, but the first partials y_{} and y_{} are represented exactly, e.g.
the bracket in Eq. (27) evaluated at for uniform f becomes
which does not vanish identically. Thus the metric identity (28) is not satisfied when the metric derivatives, y_{} and y_{}, are evaluated analytically. But it was shown above that the difference form used here, Eq. (29a), does in fact satisfy the metric identity (28) when the metric derivatives are evaluated numerically without averaging.
The use of exact analytical expressions for the metric coefficients therefore does not necessarily increase the accuracy of the difference representations, and may actually degrade the accuracy. In Chapter V it is shown further that a detrimental contribution to the truncation error can be removed by evaluating the metric coefficients numerically rather than analytically. Accuracy in the representation of the metric coefficients thus has no inherent value in and of itself. Rather, it is the accuracy of the overall difference representation that is important.
To summarize, the metric identities can often be numerically satisfied through careful attention to the evaluation of the metric coefficients. These coefficients should be expressed as differences, not by analytical expressions. They should be evaluated directly from coordinate values wherever they are needed and should never be averaged, since the use of averaged values will almost certainly result in failure to satisfy the metric identities. Intermediate coordinate values needed to construct differences which are compatible with the metric identities can be obtained by averaging the coordinates at neighboring grid points or by using a grid with twice as many points in each direction as to be used in the actual solution.
The metric identities may become more difficult to satisfy numerically in three dimensions and in schemes involving higherorder operators or unsymmetric difference expressions, as may be needed at boundaries or near the special points discussed previously. When exact satisfaction is not achieved, the effects of the spurious source terms can be partially corrected, as discussed in Ref. [12], by subtracting off the product of the metric identities with either a uniform solution or the local solution. The former amounts to using a kind of perturbation form, while the latter is, in effect, expansion of the product derivatives involving the metric coefficients and retention of the supposedly vanishing terms, thus putting the equations into a weak conservation law form. Thus the gradient could be written, using Eq. (III42), as
where A_{O} is either the local value of A or a uniform value. In view of Eq. (III40), this modification does not change the analytical expression for the gradient. Difference representations based on this form of the gradient will clearly vanish for uniform A equal to A_{O}. Analogous weaklyconservative expressions for all the other derivative operations of Chapter III can be inferred immediately. Subtraction of the product of the identity and the uniform free stream solution was used in Ref. [14], because of the difficulty in satisfying the metric identities exactly with fluxvectorsplitting involving directional differences.
All codes should be checked for the presence of spurious source terms arising from the metric identities by running with uniform nonzero values for all of the dependent variables. If such a test run produces any changes at all, some failure to satisfy the metric identities has escaped detection (assuming the code is free from errors), and the difference representations should be modified or a change should be made to the weaklyconservative form described above.
5. Implementation Procedure
When a coordinate system has been generated, the values of the cartesian coordinates, x_{i} will be available as functions of the curvilinear coordinates, ^{i}, with i = 1,2,3. Although these relations might be in the form of analytical equations in the event that the coordinate system was generated by some analytical means, a more common result is a set of values generated by a numerical solution. By definition the curvilinear coordinates take on integer values at the grid points (^{i}=0,1,2,...N_{i} where N_{i}+1 is the i total number of points in the ^{i} direction). Thus the values of x_{1},x_{2},x_{3} will be available at each grid point _{1},_{2},_{3}.
Difference expressions, such as Eq. (24), are then used at each grid point to evaluate the components of the three covariant base vectors, from Eq: (III3):
As discussed previously, the metric derivatives should not be averaged, but rather should always be evaluated directly from differences between grid points. Therefore, it may be necessary in some difference formulations to have coordinate values available at points between the grid points on which the solution is to be represented. In that case, the coordinate values at such points should be generated either by averaging the coordinate values between adjacent main points or by generating the coordinate system with twice as many grid points in each direction as will be used in the solution representation.
The nine elements of the covariant metric tensor can then be evaluated at each point from Eq. (III5):
(Only six of these elements are distinct, of course, since the tensor is symmetric, so only six dot products actually need to be evaluated.) The Jacobian is then evaluated at each point using Eq. (III16):
Next the three components of each of the three contravariant base vectors are evaluated at each point from Eq. (III33):
and the nine elements of the contravariant metric tensor are evaluated at each point from Eq. (III37):
Again only six elements are distinct.
All quantities involved in the transformed derivative operations are now available at each point. Recall that if conservative forms are to be used, the product ^{i} may be stored at each point, being evaluated from
to avoid the need for multiplication of ^{i} by in all the operations.
In transforming the physical partial differential equations, the gradient, divergence, curl, and Laplacian operations will have been replaced by either the conservative expressions, Eq. (III42)(III45), or by the nonconservative expressions, Eq. (III65)(III71). Derivatives occurring individually will have been replaced by the expressions given by Eq. (III50)(III52). Finally, derivatives occurring in boundary conditions will have been replaced by the expressions in Eq. (III78), (III79) or (III82). Integrals will have been replaced by the relations given by Eq. (III83)(IIIS7). Thus, with the metric quantities evaluated at each point, as discussed above, all quantities involved in the difference representations of the transformed partial differential equations are available.
As was noted in Chapter III, the use of the conservative forms of the gradient, divergence, curl, Laplacian, etc. in the partial differential equations is equivalent to formulation of difference equations from the integral form of these equations. Hence finitevolume formulations may be set up directly from the partial differential equations by using the conservative forms for the derivative operators involved.
It should be pointed out again that the transformed partial differential equations are of the same form and type as the original equations, and are more complicated only in the sense of having variable coefficients, crossderivatives, and more terms. The field on which these equations are solved is rectangular and the grid is fixed, uniform and square. Therefore all numerical solution algorithms that have been developed for partial differential equations on cartesian coordinate systems are applicable to these transformed equations, and all the simplifications that result from the use of uniform square grids are in order, as well.
Exercises
1. Verify Eq. (2).
2. Apply Eq. (16) to the (i,j) cell in the diagram below this equation to obtain Eq. (17). In (16) interpret the gradient as
where ^{ } and ^{ } are unit vectors in the and directions, respectively, in the transformed space. The normal is ^{ } or ^{ }, as appropriate, on each of the faces of the cell. Recall ==1.
3. Following the procedure given with Eq. (20) and (21), obtain Eq. (12) from Eq. (16) applied to the (i,j) cell. Show that the two mixed partials obtained in this manner are not equal, but that their average gives Eq. (15).
4. Verify the boundary value, f_{i,j1/2} given on p.147
5. In cylindrical coordinates show that the conservative expression for does not vanish for uniform f when the metric coefficients are evaluated analytically.