Numerical Grid Generation
Foundations and Applications
By: Joe E. Thompson, Z.U.A. Warsi and C. Wayne Mastin
V. TRUNCATION ERROR
Difference representations on curvilinear coordinate systems are constructed by first transforming derivatives with respect to cartesian coordinates into expressions involving derivatives with respect to the curvilinear coordinates (the metric coefficients). The derivatives with respect to the curvilinear coordinates are then replaced with difference expressions on the uniform grid in the transformed region. The "order" of a difference representation refers to the exponential rate of decrease of the truncation error with the point spacing. On a uniform grid this concerns simply the behavior of the error as the point spacing decreases. With a nonuniform point distribution, there is some ambiguity in the interpretation of order, in that the spacing may be decreased locally either by increasing the number of points in the field or by changing the distribution of a fixed number of points. Both of these could, of course, be done simultaneously, or the points could even be moved randomly, but to be meaningful the order of a difference representation must relate to the error behavior as the point spacing is decreased according to some pattern. This is a moot point with uniform spacing, but two senses of order on a nonuniform grid emerge: the behavior of the error (1) as the number of points in the field is increased while maintaining the same relative point distribution over the field, and (2) as the relative point distribution is changed so as to reduce the spacing locally with a fixed number of points in the field.
On curvilinear coordinate systems the definition of order of a difference representation is integrally tied to point distribution functions. The order is determined by the error behavior as the spacing varies with the points fixed in a certain distribution, either by increasing the number of points or by changing a parameter in the distribution, not simply by consideration of the points used in the difference expression as being unrelated to each other. Actually, global order is meaningful only in the first sense, since as the spacing is reduced locally with a fixed number of points in the field, the spacing somewhere else must certainly increase. This second sense of order on a nonuniform grid then is relevant only locally in regions where the spacing does in fact decrease as the point distribution is changed.
In the following sections an illustrative error analysis is given. The general development from which this is taken appears in Ref. [17], together with references to related work.
1. Order On Nonuniform Spacing
A general onedimensional point distribution function can be written in the form

(1) 
In the following analysis, x will be considered to vary from 0 to l. (Any other range of x can be constructed simply by multiplying the distribution functions given here by an appropriate constant.) With this form for the distribution function, the effect of increasing the number of points in a discretization of the field can be seen explicitly by defining the values of at the points to be successive integers from 0 to N. In this form, N+1 is then the number of points in the discretization, so that the dependence of the error expressions on the number of points in the field will be displayed explicitly by N. This form removes the confusion that can arise in interpretation of analyses based on a fixed interval , where variation of the number of points is represented by variation of the interval . The form of the distribution function, i.e., the relative concentration of points in certain areas while the total number of points in the field is fixed, is varied by changing parameters in the function.
Considering the first derivative in one dimension:

(2) 
with a central difference for f_{ } we have the following difference expression (with =1 as noted above):

(3) 
where T_{1} is the truncation error. A Taylor series expansion then yields

(4) 
Here the metric coefficient, x_{ } is considered to be evaluated analytically, and hence has no error. (The case of numerical evaluation of the metric coefficients is considered in a later section.)
The series in (4) cannot be truncated without further consideration since the derivatives of f are dependent on the point distribution. Thus if the point distribution is changed, either through the addition of more points or through a change in the form of the distribution function, these derivatives will change. Since the terms of the series do not contain a power of some quantity less than unity, there is no indication that the successive terms become progressively smaller.
It is thus not meaningful to give the truncation error in terms of derivatives of f. Rather, it is necessary to transform these derivatives to xderivatives, which, of course, are not dependent on the point distribution. The first derivative follows from (2):

(5) 
Then

(6) 
and

(7) 
Each term in f_{ } contains three differentiations. This holds true for all higher derivatives also, so that each term in f_{ } will contain five differentiations, etc.
A. Order with fixed distribution function
From Eq. (1) we have

(8) 
Therefore if the number of points in the grid is increased while keeping the same relative point distribution, it is clear that each term in f_{ } will be proportional to 1/N^{3} and each term in f_{ } will be proportional to 1/N^{5}, etc.
It then follows that the series in Eq. (4) can be truncated in this case, so that the truncation error is given by the first term, which is, using Eq. (6),

(9) 
The first two terms arise from the nonuniform spacing, while the last term is the familiar term that occurs with uniform spacing as well.
From (9) it is clear that the difference representation (3) is secondorder regardless of the form of the point distribution function, in the sense that the truncation error goes to zero as 1/N^{2} as the number of points increases. This means that the error will be quartered when the number of points is doubled in the same distribution function. Thus all difference representations maintain their order on a nonuniform grid with any distribution of points in the formal sense of the truncation error decreasing as the number of points is increased while maintaining the same relative point distribution over the field.
The critical point here is that the same relative point distribution, i.e., the same distribution function, is used as the number of points in the field is increased. If this is the case, then the error will be decreased by a factor that is a power of the inverse of the number of points in the field as this number is increased. Random addition of points will, however, not maintain order. In a practical vein this means that with twice as many points the solution will exhibit onefourth of the error (for secondorder representations in the transformed plane) when the same point distribution function is used. However, if the number of points is doubled without maintaining the same relative distribution, the error reduction may not be as great as onefourth.
From the standpoint of formal order in this sense there is no need for concern over the form of the point distribution. However, formal order in this sense relates only to the behavior of the truncation error as the number of points is increased, and the coefficients in the series may become large as the parameters in the distribution are altered to reduce the local spacing with a given number of points in the field. Thus, although the error will be reduced by the same order for all point distributions as the number of points is increased, certain distributions will have smaller error than others with a given number of points in the field, since the coefficients in the series, while independent of the number of points, are dependenton the distribution function.
B. Order with fixed number of points
An alternate sense of order for point distributions is based on expansion of the truncation error in a series in ascending powers of the spacing, x_{ }, with the number of points in the grid kept fixed and the point distribution changed to decrease the local spacing. From Eq. (9) secondorder requires that

(10) 
This is a severe restriction that is unlikely to be satisfied. This is understandable, however, since with a fixed number of points the spacing must necessarily increase somewhere when the local spacing is decreased.
The difference between these two approaches to order should be kept clear. The first approach concerns the behavior of the truncation error as the number of points in the field increases with a fixed relative distribution of points. The series there is a power series in the inverse of the number of points in the field, and formal order is maintained for all point distributions. The coefficients in the series may, however, become large for some distribution functions as the local spacing decreases for any given number of points. The other approach concerns the behavior of the error as the local spacing decreases with a fixed number of points in the field. This second sense of order is thus more stringent, but the conditions seem to be unattainable.
2. Effect of Numerical Metric Coefficients
The above analysis has assumed the use of exact values of x_{ } the metric coefficient. If the metric coefficient is evaluated numerically, we have, in place of Eq. (3), the difference expression

(11) 
The Taylor expansion yields
or

(12) 
The coefficient of f_{xx} here is the difference representation of x_{} while that of f_{xxx} reduces to a difference expression of . We thus have T_{2} given by the last two terms of T_{1} and the first term of T_{1} has been eliminated from the truncation error by evalutating the metric coefficient numerically rather than analytically.
Thus the use of numerical evaluation of the coordinate derivative, rather than exact analytical evaluation, eliminates the f_{x} term from the truncation error. Since this term is the most troublesome part of the error, being dependent on the derivative being represented, it is clear that numerical evaluation of the metric coefficients by the same difference representation used for the function whose derivative is being represented is preferable over exact analytical evaluation. It should be understood that there is no incentive, per se, for accuracy in the metric coefficients, since the object is simply to represent a discrete solution accurately, not to represent the solution on some particular coordinate system. The only reason for using any function at all to define the point distribution is to ensure a smooth distribution. There is no reason that the representations of the coordinate derivatives have to be accurate representations of the analytical derivatives of that particular distribution function.
We are thus left with truncation error of the form

(13) 
when the metric coefficient is evaluated numerically. As noted above, the last term occurs even with uniform spacing. The first term is proportional to the second derivative of the solution and hence represents a numerical diffusion, which is dependent on the rateofchange of the grid point spacing. This numerical diffusion may even be negative and hence destabilizing. Attention must therefore be paid to the variation of the spacing, and large changes in spacing from point to point cannot be tolerated, else significant truncation error will be introduced.
3. Evaluation of Distribution Functions
In Ref. [17] and Ref. [18] several distribution functions are evaluated on the basis of the size of the coefficients in the error expression. Some of this evaluation procedure is illustrated in the exercises. It appears that the following conclusions can be reached on basis of these comparisons:
(1) The exponential is not as good as the hyperbolic tangent or the hyperbolic sine. (Implementation procedures for all three of these are given in Chapter VIII.)
(2) The hyperbolic sine is the best function in the lower part of the boundary layer. Otherwise this function is not as good as the hyperbolic tangent.
(3) The error function and the hyperbolic tangent are the best functions outside the boundary layer. Between these two, the hyperbolic tangent is the better inside, while the error function is the better outside. The error function is, however, more difficult to use.
(4) The logarithm, sine, tangent, arctangent, inverse hyperbolic tangent, quadratic, and the inverse hyperbolic sine are not suitable.
Although, as has been shown, all distribution functions maintain order in the formal sense with nonuniform spacing as the number of points in the field is increased, these comparisons of particular distribution functions show that considerable error can arise with nonuniform spacing in actual applications. If the spacing doubles from one point to the next we have, approximately, x_{ } = 2x_{ }  x_{ } = x_{ } so that the ratio of the first term in Eq. (13) to the second is inversely proportional to the spacing x_{ }. Thus for small spacing, such a rateofchange of spacing would clearly be much too large. Obviously, all of the error terms are of less concern where the solution does not vary greatly. The important point is that the spacing not be allowed to change too rapidly in high gradient regions such as boundary layers or shocks.
4. TwoDimensions Forms
The twodimensional transformation of the first derivative is given by

(14) 
where the Jacobian of the transformation is

(15) 
With twopoint central difference representations for all derivatives the leading term of the truncation error is

(16) 
where the coordinate derivatives are to be understood here to represent central difference expressions,e.g.,
These contributions to the truncation error arise from the nonuniform spacing. The familiar terms proportional to a power of the spacing occur in addition to these terms as has been noted.
Sufficient conditions can now be stated for maintaining the order of the difference representations, with a fixed number of points in each distribution. First, as in the onedimensional case, the ratios
must be bounded as x_{ }, x_{ }, y_{ }, y_{ } approach zero. A second condition must be imposed which limits the rate at which the Jacobian approaches zero. This condition can be met by simply requiring that cot remain bounded, where is the angle between the and coordinate lines. The fact that this bound on the nonorthogonality imposes the correct lower bound on the Jacobian follows from the fact that implies

(17) 
With these conditions on the ratios of second to first derivatives, and the limit on the nonorthogonality satisfied, the order of the first derivative approximations is maintained in the sense that the contributions to the truncation error arising for the nonuniform spacing will be secondorder terms in the grid spacing.
The truncation error terms for second derivatives that are introduced when using a curvilinear coordinate system are very lengthy and involve both second and third derivatives of the function f. However, it can be shown that the same sufficient conditions, together with the condition that
remain bounded, will insure that the order of the difference representations is maintained.
It was noted above that a limit on the nonorthogonality, imposed by (17), is required for maintaining the order of difference representations. The degree to which nonorthogonality affects truncation error can be stated more precisely, as follows. The truncation error for a first derivative f_{ x} can be written

(18) 
where T_{ } and T_{ } are the truncation errors for the difference expressions of f_{ } and f_{ }. Now all coordinate derivatives ncan be expressed using direction cosines of the angles of inclination, _{ } and _{ }, of the and coordinate lines. After some simplification, the truncation error has the form

(19) 
Therefore the truncation error, in general, varies inversely with the sine of the angle between the coordinate lines. Note that there is also a dependence on the direction of the coordinate lines. To further clarify the effect of nonorthogonality, the truncation error terms arising from nonuniform spacing are considered.
The contribution from nonorthogonality can be isolated by considering the case of skewed parallel lines with x_{ } = x_{ } = x_{ } = y_{ } = y_{ } = 0 as diagrammed below:
Here (16) reduces to
Since , this may be written

(20) 
This first term occurs even on an orthogonal system and corresponds to the first term in (13). The last two terms arise from the departure from orthogonality. For <= 45° these terms are no greater than those from the nonuniform spacing. Reasonable departure from orthogonality is therefore of little concern when the rateofchange of grid spacing is reasonable. Large departure from orthogonality may be more of a problem at boundaries where onesided difference expressions are needed. Therefore grids should probably be made as nearly orthogonal at the boundaries as is practical. Note that the contribution from nonorthogonality vanishes on a skewed uniform grid.
Exercises
1. Verify Eq. (4).
2. Derive Eq. (6) and (7) by repeated differentiation of Eq. (5).
3. Verify Eq. (12).
4. Show that the coefficient of f_{xxx} in Eq. (12) can be reduced to a difference representation of .
5. (a) Show that with an exponential distribution function,
the ratio of the second term in Eq. (9) to the third term for very small spacing, s, at = 0 is approximately equal to 1/Ns at = 0 and to 1 at = N. Hint: Note that s = (x_{ })_{O} approaches zero as a approaches infinity, and that for large , /(e^{}1) approaches 1/e^{}.
(b) Show also that the average value of this ratio over the field is [Ns ln (1/Ns)]^{1}. Hint: Note that
(c) Finally, show that the first term in Eq. (9) causes a fractional error of approximately 1/6N^{2}ln^{2} (1/Ns) in f_{x} that does not vary over the field. (Recall that this term can be eliminated by using numerical metrics, however.)
6. Show that with a hyperbolic sine distribution function,
the ratio of the second term in Eq. (9) to the third term for very small spacing, s, at = 0 vanishes at = 0 and is approximately equal to 1 at = N. Show also, however, that the maximum value of this ratio occurs near /N = 0.9/ln (2/Ns) and is approximately equal to 1/2Ns. Finally, show that the average value of the ratio over the field is equal to[Ns ln(2/Ns)]1. Hint: See the preceding exercise. (Note that this distribution gives a smaller error due to the rateofchange in the spacing than does the exponential distribution of the preceding exercise and is particularly advantageous near = 0 where the spacing is the smallest.)
7. Show that with a hyperbolic tangent distribution function,
the ratio of the second term in Eq. (9) to the third term for very small spacing, s, at = 0 is approximately equal to 1/2Ns at = 0 and vanishes at = N. Show also that the average of this ratio over the field is the same as for the hyperbolic sine distribution of the preceding exercise. This distribution is thus also superior to the exponential distribution.
8. With the distribution function of the form of Eq. (1), show that the truncation error in Eq. (3) is a power series in inverse powers of N. (Hint: see Ref. [17]).
9. Verify Eq. (17).
10. Expand the differences f_{ } and f_{ } of Eq. (14) in Taylor series about the grid point x_{ij}. Substitute these expansions back in Eq. (14) thereby verifying Eq. (16). Certain identities will be useful, such as
11. Use the identity cos to verify the inequality in (17):
12. Use the following relations to write the truncation error in Eq. (18) in the form of Eq. (19).