Cover
Preface
1. Introduction
2. Boundary-Conforming Coordinate Systems
3. Transformation Relations
4. Numerical Implementation
5. Truncation Error
6. Elliptic Generation Systems
7. Parabolic and Hyperbolic Generation Systems
8. Algebraic Generation Systems
9. Orthogonal Systems
10. Conformal Mapping
3. Variational Approach
4. Other Approaches
5. Correlations
Exercises
12. Appendix A
Appendix B
Appendix C
References

Numerical Grid Generation
Foundations and Applications
##### By: Joe E. Thompson, Z.U.A. Warsi and C. Wayne Mastin

In an adaptive grid, the physics of the problem at hand must ultimately direct the grid points to distribute themselves so that a functional relationship on these points can represent the physical solution with sufficient accuracy. The idea is to have the grid points move as the physical solution develops, concentrating in regions of large variation in the solution as they emerge. The mathematics controls the points by sensing the gradients in the evolving physical solution, evaluating the accuracy of the discrete representation of the solution, communicating the needs of the physics to the points, and finally by providing mutual communication among the points as they respond to the physics. The basic techniques involved then are as follows:

(1) a means of distributing points over the field in an orderly fashion, so that neighbors may be easily identified and data can be stored and handled efficiently.

(2) a means of communication between points so that a smooth distribution is maintained as points shift their position.

(3) a means of representing continuous functions by discrete values on a collection of points with sufficient accuracy, and a means for evaluation of the error in this representation.

(4) a means for communicating the need for a redistribution of points in the light of the error evaluation, and a means of controlling this redistribution.

Several considerations are involved here, some of which are conflicting. The points must concentrate, and yet no region can be allowed to become devoid of points. The distribution also must retain a sufficient degree of smoothness, and the grid must not become too skewed, else the truncation error will be increased as noted in Chapter V. This means that points must not move independently, but rather each point must somehow be coupled at least to its neighbors. Also, the grid points must not move too far or too fast, else oscillations may occur. Finally the solution error, or other driving measure, must be sensed, and there must be a mechanism for translating this into motion of the grid. The need for a mutual influence among the points calls to mind either some elliptic system, thinking continuously, of some sort of attraction (repulsion) between points, thinking discretely. Both approaches have been taken with some success, and both are discussed below. It should be noted that the use of an adaptive grid may not necessarily increase the computer time, even though more computations are necessary, since convergence properties of the solution may be improved, and certainly fewer points will be required.

With the time derivatives at fixed values of the physical coordinates transformed to time derivatives taken at fixed values of the curvilinear coordinates, no interpolation is required when the adaptive grid moves. Thus, as given by Eq. (III-116),

 (1)

where xi and i are the cartesian and curvilinear coordinates, respectively. The computation thus can be done on a fixed grid in the transformed space, without need of interpolation, even though the grid points are in motion in physical space. The influence of the motion of the grid points is registered through the grid speeds, (xi)t, appearing in the transformed time derivative. This is the appropriate approach when the grid evolves with the solution at each time step. Some methods, however, change the grid only at selected time steps, and here interpolation must be used to transfer the values from the old grid to the new since the grid movement is not continuous.

In the following discussion, the problem of grid adaption will be formulated as a variational problem, the ideas being developed first in one dimension and then extended to multiple dimensions.

A. Equidistribution

A number of studies of numerical solutions of boundary-value problems in ordinary differential equations have shown that the error can be reduced by distributing the grid points so that some positive weight function, w(x), is equally distributed over the field, i.e.,

 (2)

or, in discrete form,

 (3)

where xi is the grid interval, i.e., xi = xi+1 - xi. (The subscript here indicates position on the line in this one-dimensional case.) With this condition, the grid interval will, of course, be small where the weight function is large, and vice versa. Thus if the weight function is some measure of the error, or the solution variation, the grid points will be closely spaced in regions of large error, or solution variation, and widely spaced where the solution is smooth. (It may be more appropriate in some cases to replace the equal sign in Eq. (2) and (3) with "less than or equal", and thus to "sub-equidistribute" the weight function.)

This approach has also been applied to redistribute the grid points (or to add points) at each time step, or at certain intervals, in numerical solutions of initial/boundary-value problems in one-dimensional partial differential equations. A number of references to the use of equidistribution are cited in Ref. [45]. It can be shown that the point distribution is asymptotically optimal if some error measure is distributed evenly, and that this optimum error is rather stable under perturbations of the point distribution. Thus it is not necessary to locate the grid points with excessive accuracy.

B. Equidistribution by transformation

The nonuniform point distribution can be considered to be a transformation, x(), from a uniform grid in -space, with the coordinate serving to identify the grid points. The grid points are conveniently defined by successive integer values of , making = 1 by construction and the maximum value of , i.e., N, equal to the total number points on the line. Then x = x = x so that x represents the variation in x between grid points. Hence the equidistribution statement, Eq. (3), can be represented as

 (4)

With the weight function w taken as a function of this is just the Euler equation for the minimization of the integral

 (5)

(From the calculus of variations, the function x() for which the integral F(x,x )d is an extremum is given by the solution of the differential equation . This equation is called the Euler's variational equation.) The integral (5) can be taken to represent the energy of a system of springs, with spring constants w(), spanning each grid interval, considering all the points to have been expanded from a common point so that x is the extension of the spring at . The grid point distribution resulting from the equidistribution thus represents the equilibrium state of such a spring system, i.e., the state of minimum energy. Since x represents the distance between grid points, this variational problem can also be interpreted as the minimization of the cumulative spacing between the grid points in the least-squares sense, subject to the weight function w().

If the weight function is taken to be a function of x, instead of , then the integral for which Eq. (4) is the Euler equation is

 (6)

The variational problem in this case is the least-squares minimization over the grid of the cummulative grid point spacing weighted by the weight function.

Integration over , as in both these cases, constitutes a summation over the grid points, with x representing the spacing between grid points. In the first case above, i.e., Eq. (5), the weight function w(), being a function of is associated with the grid points themselves, not with their locations. In the second case, Eq. (6), however, the weight function w(x) is associated with the locations of the grid points, rather than directly with the points. Since there is a relation x() representing the locations of the grid points, any weight function can obviously be transformed from one argument to the other. However, in deriving the Euler equations for a variational problem it is only the direct dependence that is considered in the partial derivatives or , i.e., whether the weight function is determined by the identity of the grid point or by the location of the grid point, w(x), although implicit differentiation is used in the total derivatives and .

The constant in Eq. (4) can be evaluated by normalizing x to the interval (O,L). If is normalized to (1,N) we have from Eq. (4),

and hence
 (7)

so that

 (8)

Since x=1/x , the transformation is then determined by

 (9)

Thus

 (10)

so that Eq. (2) is realized by taking equal increments in , i.e., varying by equal increments between grid points as was stated initially. From Eq. (8) the grid point spacing is given by

 (11)

An alternative viewpoint results from integrating over x, instead of over , i.e., summing over the grid intervals rather than over the grid points. Since identifies the grid points, x represents the change in ,i.e., the number of grid points per unit distance, and hence is the grid point density. Eq. (4) is now the Euler equation for minimization of the integral

 (12)

Here the integral in question is F(x, )dx, so that the Euler equation is given by .

Since x can be considered to represent the point density, this variational problem represents a minimization over the field of the density of grid points in the least-squares sense, subject to the weight function, and thus produces the smoothest point distribution attainable. Here the weight function w(x) is associated with the grid point locations, not directly with the points. If the weight function is associated with the points themselves, rather than the locations, then w = w() and the integral for which Eq. (4) is the Euler equation is

 (13)

This variational problem is the least-squares minimization over the field of the cumulative point density weighted by the weight function.

The constant in Eq. (4) is evaluated in this form by writing,

so that with the normalization as defined above,

 (14)

The transformation then is given by

 (15)

Thus

 (16)

so that again Eq. (2) is realized by taking equal increments in . The point spacing is now given by

 (17)

The grid and solution may be determined separately, perhaps even in an iterative fashion. However, the transformation allows the grid and solution to be dynamically coupled so that both evolve together. With the spring analogy approach, Eq. (8) supplies the following differential equation for the grid:

 (18)

which supplies an additional differential equation to be solved simultaneously with the differential equation system of the physical problem at hand, with the grid point location x as the independent variable. Similarly, with the smoothness approach, the differential equation for the grid is

 (19)

Eq. (18) and (19) really differ only by the way the constant is evaluated, i.e., whether by integration over or over x. This is a real difference in implementation, though, since integration over is dependent on the grid, but integration over x is not. Thus, with the spring analogy approach, the weight function is associated with the grid points, i.e., with , and the grid adjusts to achieve a uniform value of wx. The uniform value reached, however, is dependent on the grid since the right-hand side of Eq. (18) is dependent on the point distribution. In contrast, in the smoothness approach, where the weight function is associated with the the grid points, i.e., with x, the grid adjusts to achieve a specified uniform value of wx, since the right-hand side of Eq. (19) is an integral in physical space, independent of the grid. In the first approach, the points move to change the spacing x between points, while in the second the points move to change the point density x. (Note that Eq. (4) can also be written as x/w=constant.) Either approach is viable, unless it is intended that the uniform value of wx be fixed beforehand, as would be the case if the weight function is taken to be representative of truncation error and a certain bound is to be imposed on this error. The smoothness approach, i.e., integration over x, has been the most widely used because it is natural in most physical problems to associate the weight function with some physical property which varies in space.

Implementation of the two forms proceeds as follows: The form based on the grid point density is implemented using Eq. (15). With the solution u() known on the current grid points at a given time step, the weight function is evaluated at each point and then the integral in the denominator of Eq. (15) is evaluated by numerical quadrature, i.e., by summing the product wx over the grid points using coefficients in the summation appropriate to whatever type of numerical quadrature is intended. The integral in the numerator is similarly evaluated out to values of the upper limit x that produce the successive integral values of which define the grid points. Thus we have xi defined by

 (20)

These values of xi then are the new grid point locations, and the solution proceeds to the next time step.

The spring analogy form, however, requires iteration. Here we have, from Eq. (8), the point locations xi defined by

 (21)

With the solution known at a given time step, the weight function is evaluated at each grid point, and the integral in the denominator is evaluated numerically as before. Then the integral in the numerator is evaluated with the upper limit set at the successive integral values of as indicated, and this defines a changed point distribution, xi. The complication here is that the integral in the denominator, i.e., the constant in Eq. (4), depends on the point distribution, amounting to a sum of 1/w over the points since =1 by construction regardless of the distribution.

(By contrast, the corresponding integral in Eq. (20), i.e., the constant in Eq. (4), does not depend on the point distribution, being simply an integral of a function in physical space.) Therefore, this integral must be re-evaluated using the changed point distribution.

The integral in the numerator is then also re-evaluated for each point, thus changing the point distribution again. This process must be continued until convergence before the final new point distribution is obtained. The solution then proceeds to the next time step. The necessity for iteration with the spring analogy form clearly makes this form more difficult to implement than the grid point density form. Since no particular advantages of the former have been noted, preference naturally falls to the latter.

A number of examples of both the point density form and the spring analogy form, as well as other applications of the use of one-dimensional equidistribution are cited in the survey of adaptive grids given as Ref. [45].

C. Weight functions

As noted above, the effect of the weight function w is to reduce the point spacing x where w is large, and therefore the weight function should be set as some measure of the solution error, or as some measure of the solution variation. The simplest choice is just the solution gradient, i.e.,

 (22)

In this case, Eq. (4) becomes

which then reduces to

With the solution gradient as the weight function the point distribution adjusts so that the same change in the solution occurs over each grid interval, as illustrated below:

This choice for the weight function has the disadvantage of making the spacing infinitely large where the solution is flat, however.

A closely-related choice, also based on the solution gradient, is the form

 (23)

An increment of arc length, ds, on the solution curve u(x) is given by

so that this form of the weight function may be written

and then Eq. (4) becomes

which reduces to

Thus, with the weight function defined by Eq. (23), the grid point distribution is such that the same increment in arc length on the solution curve occurs over each grid interval. For the curve shown above this gives the following point distribution:

Unlike the previous choice, this weight function gives uniform spacing when the solution is flat. The concentration of points in the high-gradient region, however, is not as great. This concentration can be increased, while still maintaining uniform spacing where the solution is flat, by altering the weight function to

 (24)

where is a parameter to be specified. Considering u to be plotted against x/, we have for an increment of arc length on this solution curve

so that this weight function is equivalent to

and Eq. (4) becomes

which reduces to

Thus we have equal increments of arc length on the solution curve with u plotted against x/ in this case. Now division of the abscissa by a for a flat curve would simply reduce the spacing by the same factor. However, since the slope steepens as the curve is compressed to the left by this change of scale, the effect on the spacing where the curve is not flat will be a greater reduction in spacing.

In fact, since the 1 in the weight function given by Eq. (24) tends to produce equal spacing, while the 2ux2 tends to produce concentration in the high-gradient regions, with infinite spacing in flat regions, this weight function involves a weighted average between the tendency toward equal spacing and that toward concentration entirely in the high-gradient regions. The larger the value of , the stronger will be the concentration in the high-gradient regions and the wider the spacing in the flat regions.

Now a disadvantage of all the above forms of the weight function is that regions near solution extrema, i.e., where ux=0 locally, are treated similar to flat regions, as is illustrated below for the form given by Eq. (22):

Although the distributions produced by the solution arc length forms, Eq. (23) and (24), would have closer spacings near the extrema, the effect is still the same, i.e., to concentrate points only near gradients, not extrema.

Concentration near solution extrema can be achieved by incorporating some effect of the second derivative uxx into the weight function. A logical approach is to include this effect through consideration of the curvature of the solution curve:

If the weight function is taken as

 (25)

then points will be concentrated in regions of high curvature of the solution curve, e.g., near extrema, with a tendency toward equal spacing in regions of zero curvature, i.e., where the solution curve is straight (not necessarily flat). This weight function, however, has the serious disadvantage of treating high-gradient regions with little curvature essentially the same as regions where the curve is flat. Thus in the curve shown above, nearly all the points would be concentrated near the maximum in the curve, with very wide spacing in the high-gradient regions on both sides.

A combination of the weight functions given by Eq. (24) and (25) provides the desired tendency toward concentration both in regions of high gradient and near extrema. The effect of the inclusion of the curvature is illustrated below (cf. Ref. [37]) with the function following):

 (26)

where and are parameters to be specified. Clearly, concentration near high gradients is emphasized by large values of , while concentration near extrema (or other regions of large curvature) is emphasized by large .

Another approach to the inclusion of the second derivative is simply to take the weight function as

 (27)

where and are non-negative parameters to be specified.

With this form, (cf. Ref. [46] we have by Eq. (15), with and ,

 (28)

so that

 (29)

Then with R1 defined as

 (30)

we have

 (31)

Since = 1/N, where I+1 is the number of points on the coordinate line, the maximum percentage change in the solution over a grid interval,

 (32)

is related to the ratio R1, which measures the relative emphasis put on concentration of points according to the solution gradient by

 (33)

A guide for the choice of to limit the maximum percentage solution change over an interval to a value r can then be obtained using an equality in Eq. (33) with R1 from Eq. (30) and neglecting the effect of the term:

 (34)

The smallest possible value of r is 1/N.

With the second derivative term included, the value of can be continually updated to keep the same relative emphasis on concentration according to this term, as measured by the ratio R2.

 (35)

The transformation can then be written as

 (36)

where R2 is considered to be constant. In this form, the transformation appears as the weighted average of one based on the solution gradient and one related to the second derivative.

The replacement of Eq. (24) with the form given by Eq. (27), with = 0, still leaves a reasonable form for the weight function, but the clear association with the geometric properties of the solution curve are lost. In this case the weight function corresponding to Eq. (23) would, after substitution in Eq. (4), leads to the condition

which corresponds to an equal distribution of the distance between points on the solution curve along a right-angle path formed by x and u from one point to the next. While this distance has some indirect relation to arc length on the solution curve (the chord length being the hypotenuse of the right triangle formed by this x and u), the direct association with arc length would seem to be preferable. Following the same reasoning, the use of solution curve curvature, rather than simply the second derivative, is also preferable. Therefore, the form given by Eq. (26) is probably more appropriate than that of Eq. (27). A number of other variations have been used, of course, as is noted in Ref. [45].

Since the numerical evaluation of higher derivatives can be subject to considerable computational noise, the use of formal truncation error expressions as the weight function is usually not practical, hence the emphasis above on solution gradients and curvature. Some problems may arise even with solution curvature, i.e., with second derivatives, in rough transits. It is common in any case to limit the grid point movement at each time step and/or to smooth the new point distribution.

For systems of equations involving more than one physical variable, one approach is to use the most rapidly-varying or dominant physical variable in the definition of the weight function. Another is to use some average of the variations of the several variables. It is also possible to use entirely different grids for different physical variables, with values transfered among the grids by interpolation. Examples of each of these approaches are cited in Ref. [45] and [5].

A. Adaption along fixed lines

In multiple dimensions, adaption should in general occur in all directions in a mutually dependent manner. However, when the solution varies predominately in a single direction, one-dimensional adaption of the forms discussed above can be applied with the grid points constrained to move along one family of fixed curvilinear coordinate lines, and applications of this approach are noted in Ref. [45].

The fixed family of lines is established by first generating a full multi-dimensional grid by any of the grid generation techniques discussed in the earlier chapters, with the curvilinear coordinate lines of one family therein then being taken as the fixed lines. The points generated for this initial grid, together with some interpolation procedure, e.g., cubic splines, serve to define the fixed lines along which the points will move during the adaption. The one-dimensional adaption discussed above is then applied with x replaced by arc length along these lines.

Examples (cf. Ref. [46]) of application of the point density form discussed above in this manner are shown in the following figures. The first figure shows an adaptive grid for a combustion problem, where the adaption is along fixed radial lines. The flame front is clearly visable here because of the strong concentration of points therein:

The oscillations evident with the fixed grid are removed by the grid adaption. An extension of this problem appears next with a flowing gas. This gives an example of the use of separate adaptive grids for different physical variables of the problem, one for the combustion and one for the fluid mechanics, with values transferred between the two grids by interpolation.

Adaption through the spring analogy is illustrated next with adaption along fixed lines between the body and outer boundary in a hypersonic flow problem (cf. Gnoffo in Ref. [45]). Here the concentration of points makes the shock location evident in the grid:

Another obvious application of adaption along fixed lines is adaption of boundary points along a fixed boundary in two dimensions (cf. Nakamura in Ref. [45]). An example of such adaption along a boundary as a shock forms appears below:

One step beyond this one-dimensional adaption along fixed lines is the application of successive one-dimensional adaptions separately in each of the curvilinear coordinate directions. This proceeds in the same manner as for the adaption on the fixed lines, simply using the latest grid to re-define the coordinate lines to serve as the "fixed" lines in the next direction of adaption, cf. Ref. [56] and [57]. In the latter a torsion spring analogy is used, as well as the tension springs discussed above, incorporating resistance to movement away from orthogonality. This is done in effect by adding the term v()(x-xo)2 to the integral of Eq. (5), where v() is a second weight function and xo is the arc length location of the intersection of the normal from the adjacent grid line with the line on which the adaption is occurring.

The final grid in the one-dimensional adaption discussed above will, of course, be the result of the grid point movement along the one family of fixed lines, and therefore the smoothness of the original grid may not be preserved as the grid adapts. Some restrictions on the point movement have generally been necessary in order to prevent excessive grid distortion.

In multiple dimensions, in general it is desirable to couple the adaption in the different directions in order to maintain sufficient smoothness in the grid. One approach to such coupling is to generate the entire grid anew at each stage of the adaption from some basic grid generation system, be it algebraic or based on partial differential equations. The structure of the grid generation system serves to maintain smoothness in the grid as the adaption proceeds. In this approach, which is analogous to the one-dimensional equidistribution discussed above, the new point locations are determined directly from the grid generation system, and then the grid point speeds, , for use in the transformed time derivatives, Eq. (1), are calculated from the change in the point locations by difference expressions. Another approach is to determine the grid point speeds directly through some process and then to calculate the new point locations by integrating these point speeds.

D. Weight functions

The one-dimensional weight function, Eq. (23), based on arc length on the solution curve can be generalized to higher dimensions as follows: Consider a hyperspace of dimensionality one greater than that of the physical space, with the solution, u, being the extra coordinate. Let the unit vector in the solution direction be , this being orthogonal to the physical space. Then the position vector in this hyperspace is given by

 (37)

where is the position vector in physical space. Now, following Eq. (III-5), the covariant metric element, denoted Gij, in the hyperspace will be

 (38)

where gij is the metric element in physical space. Now

 (39)

so that

and then

 (40)

It can be shown that

 (41)

(This has been verified for one and two dimensions.)

In one dimension this reduces to the expression for arc length on the solution curve, i.e.,

In two dimensions Eq. (41) gives an expression for area on the solution surface:

 (42)

Thus the extension of the one-dimensional weight function based on arc length on the solution curve to two dimensions is that based on area on the solution surface:

 (43)

The extension of this form to three dimensions would also seem logical, but has not been verified.

3. Variational Approach

Considering the grid from a continuous viewpoint, it occurs that something should be minimized by the grid rearrangement, and thus a variational approach is logical. This is the natural extension of the equidistribution concept discussed above to multiple dimensions. The development in this section is a generalization of that in Ref. [47]. (cf. Ref. [1] for earlier related work.)

A. Variational formulation

The variational formulation for multiple dimensions can be constructed in analogy with the one-dimensional equidistribution discussed in Section 1. Thus in general a weighted integral measure of the accumulation of some grid property Q, either over the grid points, i.e.,

 (44)

or over the physical field, i.e.,

 (45)

where w is the weight function, will be minimized. The resulting Euler equations then will constitute the grid generation system. In formulating the variational problem there are basically three decision points.

First, if the integration is taken over then the integral represents a summation over the grid points, while integration over x represents a summation over cell volumes in physical space. With integration over it is thus the accumulation of some property over the grid points that is minimized, while with integration over x the accumulation over the physical cell volumes is minimized.

The second question concerns the weight function. If the weight function is directly dependent on , then the weight is associated with the grid points, while with weight functions dependent directly on x the weight is associated with location in physical space. As noted in Section 1 it is this direct dependence of the weight function that figures in the partial derivatives and in the Euler equations, the fact that a change of variable could be effected by the transformation x() notwithstanding. In most applications the weight function will be based on some solution gradient and hence will be naturally taken as a function of position in physical space, x.

Finally, there is the choice of what property is to be accumulated to be minimized. This choice depends, of course, on what is expected from the grid. Among the grid properties that might be considered are the following in computational space (integration over grid points, i.e., d):

(1). square of cell volume:

(2). inverse cell volume:

(3). sum squares of cell edge lengths (average of squares of diagonal lengths):

(4). cell area squared/volume ratio:

(5). cell skewness based on edge tangents:

(6). cell skewness based on face normals:

In two dimensions the two orthogonaly properties, (5) and (6), are equivalent.

These six properties correspond in order to the use of the following properties in physical space, where the integration is over the physical field (dx):

(1). inverse point density:

(2). square of point density:

(3).

(4).

(5).

(6).

Similar representations of other grid properties can also be considered, of course. The one-dimensional forms of properties (1) and (3) in the computational space reduce to , while those of properties (2) and (4) become 1/x . Therefore, in analogy with the one-dimensional equidistribution in Section 1, a weight function with properties (1) and (3) that is a function of x should actually be squared in the integral (cf. Eq. (6)), i.e.,

 (46a)

while w(x) with properties (2) and (4) appears as (cf. Eq. (12))

 (46b)

Similarly, weight functions that are functions of should appear as (cf. Eq. (5) and (13))

 (47a)
 (47b)

The construction for integration in the physical space is analogous, but noting that (1) and (3) correspond to 1/x, while (2) and (4) correspond to , in one dimension (cf. (5), (6), (13) and (12), respectively):

 (48a)
 (48b)

The grid for which the weighted accumulation of the property Q is minimized is obtained, by the calculus of variations, as the solution of the Euler variational equations for the integral I. If the integration is over these equations are

 (49)

where F is the integrand of the integral I. With integration over x the variational equations are

 (50)

These partial differential equations then constitute the generation system for the grid. Note that the equations resulting from Eq. (50) must be transformed using the relations in Chapter III so that the curvilinear coordinates become the independent variables. The equations given by Eq. (49), however, will already be in this form.

A grid generation system which involves competitive emphasis on various grid properties can be constructed by casting the integral to be minimized as a weighted average of several of the above integrals, each of which represents an accumulation of a different grid property. Since the various grid properties do not all have the same dimensions, it is necessary to scale the various integrals involved, as is done below for the Brackbill-Saltzman construction.

There clearly is no unique construction of the variational formulation for adaptive grids, and this is an area that is not yet fully developed. The constructions given later in this chapter are logical and illustrative of the procedure, but should not be considered definitive.

B. Euler equations

The derivation of the Euler equations, hence the grid generation system, is straightforward but may be algebraically involved. The following developments simplify the derivation somewhat. Consider first the integral over the grid points

 (51)

where g is the covariant metric tensor, with elements gij defined by Eq. (III-5), and w(x) is a weight function dependent on x. The Euler equations then are given by Eq. (49). As shown in Appendix B, the Euler equations produce the following generation system (with written as F'):

 (52)

where

 (53)

Here the gradient of the weight function in the last term is expressed using Eq. (III-42), with i given by Eq. (III-33). It should be noted that if the weight function in the integral (51) had been defined as a function of instead of x, a result different from Eq.(52) would have been obtained for the generation system (cf. Eq. (9) of Appendix B). The two-dimensional form of Eq. (52) is given as Eq. (10) of Appendix B.

With the variational problem formulated in the physical space, and the weight functions dependent on x, we have the integral

 (54)

where G is the contravariant metric tensor, i.e., with elements gij from Eq.(III-37). Then from the Euler equations given by Eq. (50), cf. Appendix B, the generation system is (with written as F'),

 (55)

with

C. Brackbill-Saltzman construction

As noted in Chapter V there is a need for smoothness in the grid in order to reduce certain terms in the truncation error of a solution done on the grid. The quantity the extension to multiple dimensions of the used above with the smoothness form in Eq. (12). Therefore to maximize the smoothness of the grid it is logical to minimize the integral of this quantity over the physical field:

 (57)

This amounts to a minimization of the linear point density in the least-squares sense. The property used here is that given as (4) on p.396, which corresponds to the ratio of the squares of the cell face areas to the cell volume when the accumulation is over the grid points, as given by property (4) on p. 395. The corresponding integral over the grid points is

 (58)

Substitution of F from Eq. (57) into Eq. (19) of Appendix B then yields the elliptic grid generation system

 (59)

Thus the smoothest grid is that for which the curvilinear coordinates satisfy Laplace's equation.

Emphasis on orthogonality and/or on concentration of grid lines can also be incorporated into the grid generation system by basing the system on the Euler equations for additional variational principles. Orthogonality can be emphasized by minimizing the integral Io defined with property (6) on p. 396 as

 (60)

since each of these dot products vanishes for an orthogonal grid. (Recall that i is normal to the coordinate surface on which i is constant, cf. Chapter III.) The inclusion of the g3/2, the cube of the Jacobian of the transformation, as a weight function in Io is somewhat arbitrary, and causes orthogonality to be emphasized more strongly in the larger cells. With the accumulation over the grid points, this corresponds to the use of the square of the dot product of the cell face normals in the variational statement (property (6) on p.395). The corresponding integral over the grid points is

 (61)

Finally, concentration can be emphasized by minimizing the integral Iw defined by

 (62)

where w(x) is a specified weight function. This causes the cells to be small where the weight function is large, and uses property (1) on p. 396, i.e., the inverse point density. With the accumulation over the grid points this corresponds to the use of the square of the cell volume (property (1) on p. 395), and the integral over the grid points is

 (63)

The grid generation system is obtained by minimizing a weighted sum I of these three integrals:

 (64)

where N is a characteristic number of points, L is a characteristic length, and W is the average weight function over the field:

 (65)

with V being the volume of the field. This sealing in the weighted sum is obtained as follows: From the above expressions for Is, Io, and Iw we have

Therefore, the three terms in Eq. (64) should stand in the ratios given. In two dimensions the factors on Io and Iw both become (N/L)4, since the Jacobian is then proportional to (L/N)3, rather than to (L/N)3. The characteristic length and number of points might logically be taken as the cube roots of the volume and the total number of points in the field, respectively, in three dimensions, the square root being used in two dimensions.

Emphasis is varied among the competing features of smoothness, orthogonality, and adaptivity by the choice of the coefficients o and w. For example, a large o will result in a grid that is nearly orthogonal, at the cost of smoothness and concentration, with an analogous effect of w. The Euler equations for this variational problem, which will be the weighted sums of those for the individual integrals, form the system of partial differential equations from which the coordinate system is generated. These equations will be quasilinear, second-order partial differential equations, with coefficients which are quadratic functions of the first derivatives, and are derived in general as described in the preceeding section and Appendix B.

Clearly the integral I3, Eq. (57), is the multi-dimensional generalization of the one-dimensional smoothness integral I3, Eq. (12), without the weight function, and the integral Iw, in Eq. (63), is the extension of the one-dimensional spring analogy integral I1, Eq. (5), to multiple dimensions, with the spring extension x generalizing to the volume, i.e., (the Jacobian in three dimensions, area in two). This variational approach thus is a generalization of the one-dimensional equidistribution discussed above to multiple dimensions. All of the discussion of weight functions given above in regard to equidistribution therefore has relevance here to the weight function of the integral Iw, Eq. (63). (The role of the constant in the equidistribution weight function, e.g., the 1 in Eq. (23), etc., which tends to produce a linear transformation, is taken by the smoothness integral Is of Eq. (57), which tends to produce an equally-spaced grid in multiple dimensions.)

For the three integrals given by Eq. (58), (61), and (63) we have, respectively, with (i,j,k) cyclic,

 (66) (67) (68)

Here, of course, from Eq. (III-14),

In two dimensions, g13 = g23 = 0 and g33 = 1, so that these functionals reduce to

 (69) (70) (71)

(Here an additive constant in Fs has been dropped since only derivatives of F contribute to the Euler equations.) Then using Eq. (1) of (Appendix B) the two-dimensional generation system based on concentration alone is

 (72)

and the generation system based only on orthgonality is

 (73)

The generation system based on smoothness (from Eq. (66)) is more complicated, but may be constructed from the relations given in Appendix B. The complete generation system then is obtained as the linear combination of the concentration system, Eq. (72), the orthgonality system, Eq. (73), and the smoothness system.

In Ref. [47] this combination is written in the form

 (74a) (74b)

where

with

Here the coefficients subscripted s, o, and w, arise from the smoothness, orthogonality, and concentration integrals, respectively. The coefficients 'w and 'o are, taking account of the scaling discussed above in connection with Eq. (64),

 (75a) (75b)

In one dimension, with y = x = 0, we have

Also, for the smoothness integral:

For the concentration integral:

and for the orthogonality integral:

Then

Now also

and, taking the -direction to be the one of interest, we also have y = 0.

The generation system in one dimension then reduces to

The can be made a part of so that the one-dimensional generation system finally is

 (76)

with, for the scaling,

 (77)

This then is the differential equation that can be applied on a boundary curve,interpreting x as arc length and as the curvilinear coordinate that varies along the particular boundary.

In three dimensions the calculation of the required partial derivatives of F in Eq. (52) for the concentration integral, i.e., Fw given by Eq. (68), may be expedited by noting that since ,

where cij is the signed cofactor of gij. These second derivatives vanish if k=i or l=j, and are equal to gmn otherwise, where (i,k,m) and (j,l,n) are cyclic, the sign being negative when the progression from i to k is opposite to that from j to l.

D. Applications

The dynamically-adaptive grid is applied by constructing the partial differential equations which constitute the grid generation system from the Euler equations as discussed above. These equations are solved numerically by replacing all derivatives with difference expressions (typically second-order, central differences) in the same manner as discussed in Chapter IV. As noted in Chapter III, the time derivatives in the equations of the physical problem to be solved on the grid are transformed according to Eq. (III-116), with the result that, the grid point speeds appear in the difference equations of the physical problem. The grid is re-generated at each time step, and these grid point speeds are determined from difference representations between time steps. Although the difference equations for the grid and those for the physical solution could be iterated together at each time step, the more common procedure is to solve each separately at each time step.

Grid points on boundaries may, of course, be held fixed, but it is more appropriate in most cases to allow the points to move along the boundary to adapt as in the field. This can be accomplished either by using Neumann boundary conditions in the grid generation systems, i.e., making the system orthogonal at the boundary (cf. Chapter VI), or by applying the one-dimensional form of the grid generation equations, Eq. (76), in terms of arc length, along the boundary.

Some rather spectacular two-dimensional results of the grid adapting to a reflected shock are shown below for supersonic internal flow over a step. The formation and multiple reflections of the shock are made evident by the grid adaption into the shock as it develops. Here the magnitude of the pressure gradient was used in the weight function, and both smoothing and bounding was applied to the weight function to control grid distortion.

E. Extensions

In two dimensions, the departure of the grid from conformality can be controlled by basing F on the Cauchy-Riemann conditions:

 (79)

and some applications are noted in Ref. [45].

Finally, another very useful addition to this composite variational system is a control on the grid point movement, which can be incorporated by taking F as

 (80)

where ui is the fluid velocity and is the grid speed. With represented by a difference form, and with the fluid velocity evaluated at the previous time step, this F can be considered to be a function of xi. Again some applications are noted in Ref. [45]. The following example shows the effectivenss of such control of the grid point movement:

4. Other Approaches

Several other approaches are discussed in Ref. [45], three of which follow here.

A. Attraction-Repulsion

Another approach to adaptive grids is to let the grid points all move as if under the mutual influence of forces between all points. Here instead of generating new grid point locations through the solution of partial differential equations, the grid points move directly under the influence of mutual attraction or repulsion between points. This is accomplished by assigning to each point an attraction proportional to the difference between the magnitude of some measure of error (or solution variation) and the average magnitude of this measure over all the points. This causes points with values of this measure that exceed the average to attract other points, and thus to reduce the local spacing, while points with a measure less than the average will repel other points and hence increase the spacing.

This attraction is attenuated by an inverse power of the point separation distance in the transformed field. The collective attraction of all other points is then made to induce a velocity for each grid point. Since each point is influenced by all other points, this is effectively a type of elliptic generation system. Details of implementation are given in Ref. [48] and other references cited therein.

Smoothing through the addition of diffusion -- like terms in the calculation of the grid evolution from the grid speeds has also been used. Reflections in boundaries in the transformed field are used to provide smooth grid motion near and on the boundaries. Since the transformed field is rectangular, this reflection is not complicated by the shape of the physical boundaries. A means of including terms that will induce rotational motion into the grid has been devised to cause the grid lines to align with lines of high gradients such as shocks.

This procedure does not exercise any control over either the smoothness or orthogonality of the grid, so that distortion is possible. Collapse of points into each other is, however, impeded because attraction will become repulsion as the points approach each other, since the measure which drives the motion will drop below the average as the spacing decreases. Collapse is further impeded by the fact that the grid velocity decreases with the spacing. It has been found necessary to apply some limits and some damping of the grid speeds to prevent grid oscillation and distortion. In practice, the computed grid speeds are scaled so that the maximum over the field is a set value, but with the maximum sealing also limited. Provision is also made for exponential damping of the grid speeds according to the ratio of the maximum Jacobian to a specified value.

Since this procedure has all grid points moving to cause some measure to approach uniformity over the field, it can be considered an iterative approach to the equidiatribution of this measure over the field. This occurs because the grid ceases to move when the measure is uniform, i.e., when the local value is equal to the average value everywhere. Therefore, the grid can be considered to move so as to minimize the variation in the measure over the field.

B. Reaction analogy

A different, but somewhat related, approach was noted in Ref. [45] and [5] based on a chemical reaction analogy. Here each grid interval is taken to represent a species concentration, and the reaction rate constants are made dependent on the difference between a local error measure for one grid interval compared with another. Each grid interval then is coupled with every other grid interval through reaction rate equations, so that each interval grows at the expense of others, and vice versa. A system of ordinary differential equations is solved for the intervals. This approach, as given, is somewhat inefficient, since there is no provision for limiting the effect to the nearer points. With each point affected equally by all other points, the number of ordinary differential equations to be solved is equal to the square of the total number of points.

The rate constants also contain factors designed to limit the range of variation of the grid intervals. The two-dimensional form given involves essentially applying the one-dimensional form separately along each family of curvilinear coordinate lines, with spacing in one cartesian coordinate being adjusted along one family of curvilinear lines, and the other cartesian coordinate being adjusted along the other family.

C. Moving finite elements

The moving finite element method of Miller (Ref. [49] -- [50]) is a dynamically-adaptive finite element grid method in which the grid point locations are made additional dependent variables in a Galerkin formulation. The solution is expanded in piecewise linear functions, in terms of its values at the grid points and those of the grid point locations on each element. The residual is then required to be orthogonal to all the basis functions for both the solution and the grid. The grid point locations are thus obtained as part of the finite element solution. An internodal viscosity is introduced to penalize the relative motion between the grid points. This does not penalize the absolute motion of the points. An internodal repulsive force was also introduced to maintain a minimum point separation. Both of these effects are strong but of short range. A small long range attractive force is also introduced to keep the nodes more equally spaced in the absence of solution gradients. Small time steps are used in the initial development of the solution. The results show that the oscillations typically associated with shocks with fixed grids are removed with the adaptive grid, and that dispersion and dissipation are essentially eliminated. An order-of-magnitude increase in stability was also realized over conventional methods.

5. Correlations

The ultimate answer to numerical solution of partial differential equations may well be dynamically-adaptive grids, rather than more elaborate difference representations and solution methods. It has been noted by several authors that when the grid is right, most numerical solution methods work well. Oscillations associated with cell Reynolds number and with shocks in fluid mechanics computations have been shown to be eliminated with adaptive grids. Even the numerical viscosity introduced by upwind differencing is reduced as the grid adapts to regions of large solution variation. The results have clearly indicated that accurate numerical solutions can be obtained when the grid points are properly located.

It is also clear that there is considerable commonality among the various approaches to adaptive grids. All are essentially variational methods for the extremization of some solution property. The explicit use of varational principles allows effective control to be exercised over the conflicting requirements of smoothness, orthogonality, and concentration, and this is probably the most promising approach in multiple dimensions.

The adaptive grid is most effective when it is dynamically coupled with the physical solution, so that the solution and the grid are solved for together in a single continuous problem. The most fruitful directions for future effort thus are probably in the development and direct application of variational principles and in intimate coupling of the grid with the physical solution.

Exercises

1. Show that Eq. (4) is the Euler equation for the minimization of the integrals (5), (6), (12), and (13). Hint: For (6) note that in the term , w must be differentiated with respect to implicity, i.e., w =wxx . A similar situation occurs with (13). Note, however, that implicit differentiation is not to be used in the term for (5) or in for (12).

2. Show that Eq. (4) is also the Euler equation for the integrals

3. With the weight function given by w(x)=sin(x/L), find the grid point locations from Eq. (20). Note the concentration near x=L/2 where the weight function has its maximum value.

4. For u(x)=(L/)sin(x/L), obtain the point distribution from Eq. (20) using the weight functions from Eq. (22), (23), (25), (26) and (27). Use ==1. Plot and compare.

5. Show that the average of the squares of the diagonal lengths is .

6. Verify the correspondence between the six grid properties listed on p. 395 with the six listed on p. 396. Hint: Recall that dx=d.

7. Verify that the one-dimensional forms of the first four properties on pp. 395-396 are as stated on p. 396-397. Hint: In one dimension take

8. Show that Eq. (59) is the Euler equation resulting from the integral given by Eq. (57).

9. Verify Eq. (72) and (73).

10. Show that with

the generation system is

Hint: Use Eq. (19) of Appendix B.

11. Show that with the generation system consists of Laplace equations in the computational space.

12. Show that with

the generation system is

13. Show that with , (i,j,k) cyclic, the generation system is

Hint: Note that and .