- Email: [email protected]
Solving partial differential equations using recursive grids
__ _EB
aA
ELSEVIER
APPLIED NUMERICAL MATHEMATICS Applied
Numerical
Mathematics
14 (1994) 165-181
Solving partial differential equations using recursive grids Gregory
Kozlovsky
*
Department of Computer Science, City College of New York, New York, NY 10031, USA
Abstract In this paper we discuss numerical methods and software tools for solving partial differential equations with variable local resolution. The emphasis of our work is on finding approaches which will greatly simplify developing applications using dynamic adaptive grids. First, we investigate a class of adaptive grids that looks particularly promising for use in general-purpose solvers and propose a discretization method suitable for this class of grids. Second, we describe an interactive programming environment which can be used, among other things, as a geometrically-oriented debugger for developing applications using adaptive grids. Finally, we briefly illustrate an application of the above-mentioned methods and software to flame-flow interaction problems.
1. Introduction
When solving partial differential equations, it is often necessary to resolve small areas with high activity. ’ Blame fronts are an example of such areas. High local activity can also be found in boundary layers. It may be not computationally feasible to use grids which are uniformly spaced throughout the computational domain to resolve these small high activity areas. One increasingly common solution in this situation is the use of adaptive non-uniform grids [1,3-6,101. Adaptive non-uniform grids have different densities in different parts of the computational domain, as dictated by requirements of the accuracy of the numerical method. Because of the potential that adaptive gridding has for reducing computational costs, it is one of the forefront areas in computational physics. Many types of non-uniform grids are used in numerical computing. In general, the choice of a grid type represents a trade-off between simplicity of the grid and its ability to cover the computational domain economically with the locally required resolution. In deciding which type of grid to employ for a given problem, much depends on the geometrical configuration of the area(s) in need of high resolution. In many cases, relatively simple grids, such as tensor-product * Current address: CN Division, CERN, 1211 Geneva 23, Switzerland. E-mail: [email protected]. ’ What exactly constitutes high activity is problem-dependent. It may mean a high gradient, a large second of the solution, or something else. 016%9274/94/$07.00 0 1994 - Elsevier SSDI 0168-9274(93)E0069-E
Science
B.V. All rights reserved
derivative
166
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
grids or composite uniform grids, can do the job efficiently. We will describe, without the pretence of being comprehensive or rigorous, several typical situations where relatively uncomplicated grids can be used. For simplicity, we restrict our exposition to the two-dimensional case. l
l
l
The high activity area is centered around a singular point. A radial grid is formed by the intersection of circles centered at the point with rays emanating from it. A radial grid is actually a tensor-product grid and, as such, can be implemented using simple data structures. The density of such a grid naturally increases with proximity to the central point. The high activity area is concentrated around a line. Non-uniform tensor-product grids are well suited to this case. By suitable mapping of the computational domain the case of a simple curve can be accommodated. The high activity area can be covered by a rectangle of the size comparable to that of the high activity area. Composite uniform grids can be employed as well as non-uniform tensor-
product grids. The advantages of the former should be weighed against the increase in the complexity of computer implementation. For many important problems, however, the shape of the high activity area is more complicated than in the cases outlined above. In addition, in order to be efficient the above approaches should rely on certain geometrical shapes of the area to be refined, and a “black-box” solver must incorporate some sort of pattern recognition ability, a significant complication to be avoided. Another major consideration in grid design is the ability of the grid to change with the solution. In the simplest application of non-uniform gridding, a variably spaced grid is established initially and then held fixed through the course of the calculation. This approach to grid generation is used in many industrial finite element packages. It requires a priori knowledge of the general behaviour of the solution (usually an intuition of an experienced engineer) and may require many months of labour. Dynamic adaptive gridding methods, on the other hand, change the grid during the solution process according to the accumulated information. These methods refine the grid in the high activity areas, leave the grid coarse in the areas of little activity, and unrefine the grid in the areas where refinement is no longer needed, Automatic mesh refinement and unrefinement is absolutely necessary in case of time-dependent problems with moving areas of high activity, or when the location of the areas requiring high resolution is not known in advance. For highly nonlinear problems a large number of time steps has to be performed. A desirable property of a dynamic grid in this case is that a small change in the solution causes a small change in the grid. Tensor-product grids and, in general, composite uniform grids do not possess this property (at least when uniform subgrids are large). In this paper we investigate a class of recursive grids we call consistent recursive grids that looks particularly promising for use in general-purpose solvers. One major drawback of recursive grids is that they are more complex than tensor-product grids and, as a result, require complex data structures. In addition, numerical methods have to be reformulated to cope with the complexity of the local geometrical structure of recursive grids. Consequently, programming and, especially, debugging of applications using recursive grids become more complex. During
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
167
our research, we found out that making use of such grids possible and practical involves a whole cluster of interconnected multi-disciplinary questions including data structures, discretization methods, visualization, and interactive user interfaces. Developing PDE solvers using conventional compiled languages such as Fortran or C is time-consuming and difficult. One of the main stumbling blocks is the lack of tools for visualizing data and solutions during the development process. Engineering packages (most commonly used in structural mechanics) do have integrated interactive graphics modules which are extensively used for data debugging. Clearly, there is a similar need for interactive visualization during the process of developing new PDE solvers. As we move from relatively simple tensor-product grids to recursive grids, the need for visualization tools integrated with the program development environment becomes imperative. As a solution for this problem we propose an interactive programming environment [8] for developing adaptive grid solvers. The environment is intended to reduce drastically the amount of technical, non-problem-specific programming and can be used, among other things, as a geometrically-oriented debugger for developing applications using adaptive grids. It is currently implemented on a PC-AT and on an IBM RS-6000 workstation. In designing this environment, we aspired to provide a set of tools for the researcher rather than to create an intelligent system which “understands” PDEs and adaptive grids. A substantial effort in creating higher-level packages may thereby be saved by implementing them on top of properly designed lower-level tools. In the following sections we describe in detail consistent recursive grids, their properties, algorithms for grid handling, and data structures for grid implementation. Discretization methods on consistent recursive grids and interactive software tools are discussed only very briefly. For more details the reader is referred to [8]. In the last section, we illustrate an application of the interactive programming environment to flame-flow interaction problems.
2. Adaptive grids In selecting the type of adaptive grids suitable for a general-purpose following governing principles:
solver, we adopted the
(a) Keep the grid geometry simple and tackle the geometric complexity of the computational domain with a combination of grid refinement and special boundary finite elements. (b) Create grids which are as “tight” (economical) as possible even if it complicates programming and makes parallelization more difficult. If the library of grid management functions is reusable, the complexity of these functions is less important than in the case of ad hoc written applications. The question of parallelization can be dealt with later and, in any case, it is only one factor in the overall performance of a method. From a performance point of view, we looked for grids that could handle areas with complex geometry, have the ability to move refined area(s), be flexible enough to be used in a general-purpose environment, and allow for a high degree of local refinement. The following exposition is carried out for the two-dimensional case. The generalization for the three-dimensional case is straightforward.
168
G. Kozlovsky /Applied
Numerical Mathematics
14 (1994) 165-181
Fig. 1. Grid fragment with three grid levels. Higher grid levels are placed above the lower ones. The resulting three-dimensional picture is shown in oblique projection. The boundaries of the supercells and the lines dividing each supercell into four cells are included. Nodes are marked by black dots. Vertical lines connect co-located nodes on adjacent levels.
2.1. Definitions and terminology The main unit for building our grids is a supercell. Geometrically, a supercell is a rectangle divided by lines parallel to its sides into four equal subrectangles. Each of these subrectangles we will call a grid cell or computational cell. Grid nodes are formed by the vertices of the cells. Nine nodes are associated with each supercell. A grid starts its life as a root grid, which consists of a number of non-overlapping supercells of equal size completely covering a rectangle. ’ An elementary grid refinement is performed by creating a supercell exactly covering an existing cell (refining the cell). The refined cell is called the parent cell of the new supercell, which, in turn, is called its descendant. Refinement can later be reversed by destroying a previously created supercell (unrefining the cell). Each supercell has a grid level number associated with it and its cells. The grid level of a root supercell is zero. The grid level of a non-root supercell is equal to the level of its parental cell plus one. Nodes associated with cells on different grid levels and having identical coordinates are considered distinct and are called co-located nodes. A set of nodes co-located with a given node includes the node itself. There is only one node with the given coordinates on any grid level and, therefore, a node may be associated with more than one supercell on the same level. The grid level of a node is equal to the grid level of its associated supercells. A grid fragment with three grid levels is shown in Fig. 1. A terminal cell is a cell which does not have the descendant supercell. A terminal supercell has all four of its subcells terminal. An active supercell contains at least one terminal cell. A terminal node is a node which is not co-located with a node on a higher level.
a We are not concerned with more elaborate root grids, because a root grid is inexpensive computational domain can be economically covered using appropriate grid refinements.
and a non-rectangular
G. Kozlo~dy
/Applied
Numerical Mathematics
14 (1994) 165-181
169
It will be useful to classify the nodes by their geometrical location in an associated supercell: . l x-node: node in the center of a supercell, l c-node: node in one of the corners of a supercell, l t-node: node on a side of a supercell which is associated exclusively with this supercell (no adjacent supercell on the same grid level), l s-node: located as t-node but is associated with two adjacent supercells. Node types are illustrated in Fig. 2. Another important distinction is between internal cells and nodes, and border cells and nodes which are the cells adjacent to and the nodes located on the border of the computational domain. Internal nodes which are not t-nodes we will call r-nodes, where r stands for regular. The difference between x-nodes, c-nodes, s-nodes, and t-nodes is important from the point of view of grid data structures and basic grid algorithms. For numerical algorithms x-nodes, c-nodes, and s-nodes are all the same, except that, as we will see later, there can be more variety in the configuration of the grid around c-nodes. However, numerical methods have to be formulated differently on the r-nodes and on the internal t-nodes. Two grid cells are called neighbows or adjacent cells if they share at least one pair of co-located nodes. A pair of grid nodes are neighbours on level k if both of them have co-located nodes which are associated with the same cell of level k. Two terminal nodes are called neighbows if both nodes are co-located with nodes associated with a terminal cell. A grid is called consistent if the linear size of any adjacent terminal grid cells differs by at most a factor of two (see Fig. 31. This requirement, which is referred to as consistency, results in a certain regularity of the grid-namely, that the local structure of the grid is limited to few possible configurations. Such regularity proves advantageous from several points of view: it facilitates grid data bookkeeping in basic grid algorithms (refinement, unrefinement, grid traversal, and nearest node search), simplifies derivation and programming of numerical methods on the grid, and allows rows of stiffness matrices to be preassembled. A grid can be represented as a forest of quad-trees, with supercells as tree nodes. Unlike supercells, the grid nodes does not naturally form a tree, this being the reason why supercells are chosen as the basic building block of the grid. The tree structure of the supercells allows for certain important operations to be performed fast. Therefore, even when we use only nodes for our numerical methods, we have to keep both supercells and nodes in our data structures.
c-node
t-node
c-node
c-node
t-node
c-node
Fig. 2. Supercell
and node types.
170
G. Ko.zlousky /Applied
Numerical Mathematics
14 (1994) 165-181
Fig. 3. Consistency.
We will frequently need to refer to the relative location of neighbouring grid components (nodes, cells, and cell interfaces). We say (see Fig. 4) that nodes are located in the directions x-, x+, y-, and yf, and in the orientations x-y-, x+y-, x-y+, and xfy+ relative to the central node. It is convenient to have different generic names for directions and orientations because of their different geometrical position relative to the central node and the different role in the grid data structures. The same notation will be used for the relative positions of cells, for positions of cells relative to nodes, and vice versa. Now we can precisely define the local geometrical configuration. The geometrical configuration around a terminal node is a correspondence between the orientations and the relative levels of the terminal cells adjacent to the node, where the relative level of a cell is an integer obtained by subtracting the central node grid level from the cell grid level. A geometrical configuration is identified by its configuration number. 2.2. Properties of consistent grids Consistent definitions.
grids have the following properties,
x-y+
2-
x
most of which result directly from the
Y+
x+y+
Y-
x+y-
+ixt
Y
Fig. 4. Location
of neighbouring
nodes relative
to the central
node.
G. Kozlovshy /Applied
Numerical Mathematics
Property 2.1. There is a one-to-one correspondence x-nodes.
14 (1994) 165-181
171
between the supercells and their associated
Property 2.2. The terminal cells surrounding a terminal node are on at most two different grid levels-the level of the node and one level below it. At least one of these cells is associated with the node and is, therefore, on the same grid level. Moreover, all four terminal cells surrounding a terminal x-node or a terminal s-node are on the same level. A terminal internal t-node is surrounded by three cell-two cells associated with it and one on the level below. Only a terminal c-node can have several different grid configurations around it. Proof. Follows from the definition
of consistency.
•I
Property 2.3. A node associated with a terminal cell is either terminal or has a co-located terminal node on the next grid level. Proof. Consider the terminal node pt co-located with the given node p. According to Property
2.2 the terminal cells surrounding a terminal node can be on at most two different levels. Therefore, the highest level cell around pt can be at most one level above the level of the given terminal cell. The same is true for the levels of nodes p and pt which are associated with the aforementioned cells. 0 Property 2.4. Grid levels of two neighbouring terminal nodes differ at most by one. Proof. By definition two neighbouring terminal nodes are co-located with nodes of a terminal cell. Therefore, according to Property 2.3, the grid level of each node is either equal to the level of the cell, or is one level above it. 0 Property 2.5. An internal t-node of a consistent grid is always terminal. Proof. Given a t-node p on level 1, suppose there is a node fi on level I+ 1 co-located with node p. By the definition of an internal t-node, there is a terminal cell on level I- 1 adjacent to node p. On the other hand, there is a cell of level 1 + 1, associated with node 5, which is also adjacent to node p. Thus, the starting assumption would violate the grid consistency. •I
The following simple theorem expresses a fundamental property of consistent recursive grids, their local semi-regularity. This property distinguishes consistent recursive grids from less structured grids which can have infinitely many local configurations. Such semi-regularity proves advantageous from several points of view: it facilitates grid data bookkeeping, simplifies derivation and programming of numerical methods on the grid, and allows rows of stiffness and mass matrices to be preassembled. Theorem 2.6. The number of possible distinct geometrical configurations around an internal terminal r-node is fifteen (see Fig. 5).
G. Kozlovsky /Applied Numerical Mathematics I4 (1994) 165-181
172
Fig. 5. Nontrivial rotation.
geometrical
configurations
around
a grid node. The remaining
configurations
can be obtained
by
Proof. It follows from Property 2.2 that terminal cells of only two different sizes can surround
a terminal node and that the case where all the surrounding cells are on the level below the node level is not feasible. The theorem follows from the fact that there are four terminal cells around an internal terminal r-node. LI
2.3. Grid data structures The design of grid data structures must be a compromise between the amount of memory required and the support for the efficient execution of certain basic operations on the grids. Given the current state of computer hardware and the type of problems we deal with, speed rather than memory is the bottleneck. Therefore, in designing the data structure we gave preference to achieving higher speed over saving memory. The following is a list of basic grid operations which must be efficiently supported by grid data structures. Grid refinement/unrefinement.
Given a consistent recursive grid, a set of cells to be refined, and a set of supercells allowed to be unrefined, find a modified consistent grid in which all the necessary cell refinements are performed and as many supercells as possible are unrefined. Grid traversal. Traverse subsets of nodes in different order. Subsets of nodes can be defined by node types (such as t-nodes or r-nodes), grid levels, and node terminality. Among frequently used node orderings are the x-y, y-x orderings, and the nested dissection ordering. Nearest neighbours search. Given a node, find its neighbouring nodes and/or neighbouring cells. Given a cell, find its neighbouring cells and/or associated nodes. Point location relative to the grid. Given a point in the computational domain, find the node at the given grid level which is the closest to the point and the cell which contains the point. Support for numerical methods using hierarchical structure of the grid. 3 Given a node, find the co-located nodes on the adjacent grid levels. This operation is required, for example, in order to transfer grid vector values between different grid levels. Two arrays are used to store the information about the nodes, the supercells, and their interconnections. The array c t r e e keeps the instances of the structure 4 SUPER C E L L, and the 3 Such as the multigrid method and the nested dissection 4 The adaptive grid library is implemented in C.
method.
173
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
typedef struct supercell ( float ctr[2]; /* gn_index ctrnode; /* sc_index parcell; /* int load; /* int maxlevel; /* char pos; /* char level; /* ) SUPERCELL; typedef struct gnode c gn_index nbr[41; gn_index fine; union { gn_index coarse; sc_index scell; 1 u; unsigned char gtype; unsigned char gdir; ) GNODE;
/********** supercell structure *********/ coordinates of the super-cell center */ central node of the supercell */ index of the parent supercell */ work load associated with cell sub-tree */ maximum level of descendants *c/ position in the parent supercell */ supercell's level (root level 1s 0) */
/********* grid node structure ****f****/ /* indexes of 4 neighbouring nodes */ /* co-located next level node */ /* co-located previous level node /* associated supercell for x-node
*/ */
/* node type and inconsistency bit */ /* encoded locations of adjacent cells */
Fig. 6. Node and supercell
structures.
array g n o d e s contains the instances of the structure G N 0 D E. Each instance of these structures is identified by its index in the corresponding array. The dynamic allocation module keeps track of the free slots in c t r e e and g n o d e s. When a new supercell or node is created, a slot for it is allocated; when a supercell or a node is deleted its slot is returned to the list of the free slots. The definitions of node and supercell structures are listed in Fig. 6 with some comments. Supercell structure has references to its parent, to its children, and to its central node. Node structure has references to its immediate neighbours, to the co-located node of the parental supercell (the central node of a supercell has a reference to its associated supercell instead), and to the co-located node of a child supercell (if any). In addition, node structure contains the encoded local grid configuration and the encoded node location in its associated supercells. 2.4. Grid creation, refinement, and unrefinement Grid adaptation is done in steps which may be interspersed with numerical computations. One step of grid adaptation is a solution to the following problem. Definition 2.7 (Grid adaptation problem). Given a grid .?9’,a set of its cells 9Z’to be refined, and
a set of its supercells Z/ which are allowed to be deleted, find a minimal consistent grid in which all the cells of 9 are refined and all the supercells not belonging to ?Y are preserved. A minimal grid is a grid from which no supercells can be deleted without violating one of the above conditions. It is easy to see that because of the consistency requirement the minimal grid may contain some supercells of % marked for deletion, and some cells not belonging to 9 will be refined.
174
G. Kozlousky /Applied
Numerical Mathematics
I4 (1994) 165-181
refine(s) { For each neighbouring
cell s,, do {
If s, is one level below s, ref ine(s,) 1 create(s) 1 Fig. 7. One-cell
consistent
refinement
algorithm.
Theorem 2.8. For the given grid 22 and the given sets 2 and 2Y there exists a unique minimal grid. Proof. Suppose there are two different minimal grids ZYIand P2. Without restricting generality, we can presume that the set 2YI\PY2 is not empty and that the supercell s’ is among its supercells on the highest grid level. Then, the parent cell of S does not belong to the set 9, otherwise it would have been refined in g2. For the same reason, the supercell s’ is not among the supercells of Z? which cannot be deleted. Finally, the deletion of s’ will not cause inconsistency in ZYr,because it is absent in .Y2 and there are no supercells in .YI\3’2 which are above the grid level of s”. Therefore, contrary to the initial assumption, gI is not a minimal grid. 0
Let c r e a t e(s) be the function which creates the supercell covering the cell s. Then, the recursive function r e f i n e(s) presented in Fig. 7 refines the terminal cell s, preserving the grid consistency. Moreover, the grid remains consistent during the intermediate steps of the algorithm, because a supercell is created only after all of the neighbours on the grid level below its level have been created. Let de s t r o y(s) be the function which deletes the supercell covering the cell s. Then, the recursive function un r e f i n e(s) presented in Fig. 8 unrefines the cell s if possible without violating the grid consistency. Because of consistency, the solution to the grid adaptation problem is global in the sense that refinement or unrefinement at a particular location of the grid may depend on the situation at a different location. A naive approach to the problem can lead to repeated unrefinements and refinements of the same cell. This would not only require unnecessary
unrefine(s) boolean nounrefine = FALSE If s is non-terminal, nounrefine = TRUE For each neighbouring cell s,, do { If s, is more than one level above s, nounrefine 1 If nounref
ine is FALSE, destroy(s)
Fig. 8. One-supercell
consistent
unrefinement
algorithm.
= TRUE
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
175
adapt (1 For 1 = 0,. l,,,, do { For each supercell s E 721, refine(s) For 1 = l,,,, . .O, do { For each supercell s E L/l, unrefine(s)
Fig. 9. Optimal
grid adaptation
algorithm.
computing, but might also lead to the loss of numerical data associated with temporarily deleted nodes and cells. We propose an algorithm presented in. Fig. 9 which guarantees that this does not happen. The following notations are used: 1 is the current grid level, I,,, is the maximal grid level, &J?[and ??J1are the subsets of the sets 32 and Z, respectively, consisting of supercells on the grid level 1. Theorem 2.9. The grid adaptation algorithm is optimal in the sense that the resulting grid is minimal and that no cell refined at an intermediate stage of the algorithm is later unrefined. Proof. Let c’ be the highest level cell which was refined
and then unrefined. Because it was unrefined, c’ cannot belong to 32, and, therefore, it was refined because there was a neighbouring cell E+ on the grid level above that had to be refined. By the definition of 2, c’+ cannot be among the cells later unrefined. Hence, c’ too could not have been unrefined without violating consistency. Suppose now that the grid .?Yaresulting from the application of the grid adaptation algorithm is not minimal. Let c’ be a highest level cell which is refined in 27a but not in the minimal grid. Then, 2 does not belong to the set 33”.Therefore, c’ either was created because there was a neighbouring cell c’+ on the grid level above its level that had to be refined, or c’ belongs to the set 22 of the original grid. In the first case, by the definition of E, E+ must be refined in the minimal grid, and, therefore, c’ must also be refined in the minimal grid, which contradicts the definition of c’. In the second case, nothing prevents c’ from being deleted during the unrefinement phase of the algorithm. q 2.5. Grid traversal and subsets
Parameters of a numerical method can be node based, cell based or cell interface based. We will use the word grid component to refer to any of these three types. Most often numerical methods work with node based parameters. A grid is needed to support operations performed on data based at grid nodes, cells, or interfaces. Programs performing these operations must have the ability to traverse the grid. The order of grid traversal can be important. Grid components can be ordered by grid level, by geometrical order, by component type (such as t-nodes and r-nodes), or by any combination of the above orders.
176
G. Kozlovsky /Applied
Numerical Mathematics
14 (1994) 165-181
Two sets of tools are provided for grid traversal-the basic grid traversal functions and index lists. Basic grid traversal functions (BGTFs) traverse the grid using its internal tree structure. Index lists require some additional memory but are much more flexible than BGTFs. They allow a set of grid components to be split into any number of subsets and maintain a specific order inside every subset. Several standard splittings are provided with the system. Users can create new splittings and orderings. 2.6. Search for neighbours The grid data structure contains indexes of the neighbouring nodes on the same grid level for every direction. It also contains indexes of the co-located nodes on the adjacent levels. Therefore, all the neighbouring nodes can be found in O(1) operations. The connections between nodes and supercells are made through x-nodes. Because of this, the connections between nodes can be used to find neighbouring cells as well. 2.7. Adaptive grid library Using the data structure described, algorithms for refinement, unrefinement, grid traversal, and nearest node search for consistent recursive grids, together with necessary utilities, are implemented as a library of C functions for adaptive grid handling. The library is integrated with the programming environment and with the interactive graphics which are outlined in Section 4. At present the implementation is limited to two-dimensional grids.
3. The finite element method on recursive grids
The finite element method [12,13] is an especially attractive spatial discretization method on recursive grids because the derivation of the discrete equations on the interface between cells of different sizes is somewhat more natural than in the finite difference method, and because an h-p adaptation can potentially be employed. The treatment of the interface between computational cells of different sizes is a major problem in using recursive grids. Many authors use interpolated solution from the course grid as a boundary condition for patches of fine grid. We discretize a differential operator at all the r-nodes, including the ones on the interface between computational cells of different sizes. For such interface nodes, standard formulation of the finite element method can lead to undesirable increase in number of nodes in the discretization stencil. To solve this problem, we developed a version of the finite element method for consistent recursive grids. Our version uses discontinuous piecewise bilinear basis functions to simplify matrix generation on interfaces of computational cells of different sizes. The solution remains continuous due to additional constraints imposed on the unknowns. A particular version of the finite element method (in the weighted residual formulation) is defined by the choice of the basis functions and the weight functions.
177
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
3.1. Basis functions in space
We shall seek unknown functions in the approximate
form
(1) where +i are basis functions of space variables x and y, and a, are node based parameters to be computed. A basis function $i is a piecewise bilinear function which equals 1 at terminal node i and equals 0 at the rest of the terminal nodes. Let z be the set of elements (cells) which have node i as one of their vertices. Then 4i is non-zero only over the elements belonging to A, and is a bilinear function over each element of A. The support of the function $i, defined as the sum of the elements belonging to A, will be denoted by Si. In our case, where the elements are rectangles in two-dimensional space, the number of elements in A is four for an internal r-node and two for a t-node. Assuming that the coordinate system starts at node i, its axes are parallel to the sides of the elements, and the element j EA is located in the first quadrant, the restriction of a basis function 4i over the element j can be written as
4i=h”F’
hf-x
h;-y
I
I
(2)
where hr and hy are the x and y dimensions of the element j. The support of basis functions for several geometrical configurations of the grid is illustrated on Fig. 10. When elements of A have different sizes, or the node i is a t-node, the function 4i is apparently discontinuous. The continuity of the expansion (1) is preserved by imposing the following linear constraints on the admissible set of parameters ai ai = i(ai++ ai-)
for every internal t-node, where i+ and i- are the nodes on the ends of the cell side in the middle of which the t-node i is located. An alternative approach would be to eliminate the t-node based basis functions altogether and use continuous basis functions at the r-nodes. These functions would have larger support than basis functions c$~,which would lead to a more complicated discretization pattern and, ultimately, would also complicate programming.
Fig. 10. The support of basis functions for several geometrical marked by black disks, the remaining nodes of the discretization
configurations of the grid. The central stencil are marked by circles.
nodes
are
178
G.
Kozlovsky/Applied Numerical Mathematics 14 (1994) 165-181
3.2. Weight functions
As weight functions we take normalized linear combinations (cli=
CLi
4i
+
+
i
C
kE6
4k
of basis functions
i
where q is the set (possibly empty) of t-nodes adjacent to the node i, and pi is chosen SO that the integral of lcIiover its support Qi is unity. Note that rcIiare continuous functions of x and y. 3.3. Preassembling stiffness matrices
As was mentioned before, consistent recursive grids can have only a limited number of local geometric configurations at a node. If we use a version of the finite element method described above the number of possible discretization stencils remains equal to the number of local geometrical configurations of the grid. Consequently, rows of the stiffness matrix approximating a differential operator with constant coefficients can be preassembled for every configuration. This property leads to substantial savings of computational time.
4. The interactive
programming
environmed
The foundation of the proposed interactive programming environment [8] is an interpreter. The environment also includes an interactive graphics module and a library of functions for handling recursive adaptive grids. Integration of an interpreter with interactive graphics provides a way to connect geometrically-based human understanding with underlying data structures. This connection is important because conceptual and programming mistakes that are very difficult to find by conventional means manifest themselves through our spatial perception. We call the process of finding such mistakes geometric debugging. 4.1. The interpreter For the most part, the input language of the interpreter mimics the lexical conventions and syntax of the C programming language [7]. Most of the C fundamental data types are supported. Several specialized built-in classes and attributes have been added to the input language to facilitate work with adaptive grids. They include node- and cell-based grid vectors, grid objects, and space functions. Arrays declared as grid vectors are given special treatment by the interactive graphics module. Grid objects include data and functions to be executed upon the data every time the grid is refined or unrefined. Space functions accept a location in the computational space as their input, and output their results in a pop-up window using special output operators. Space functions are intended to be called from the interactive graphics module with a location specified by a pointing device (mouse). Let us consider grid objects in more detail. Conceptually, a grid object is an object which must be modified whenever the grid is modified. For example, when the grid is refined, the components of a grid vector corresponding to the newly created nodes have to be initialized. A
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
179
grid cell is refined either at the request of an application program or by functions of the grid management library in order to maintain the consistency of the grid. Therefore, a uniform mechanism to update grid objects is desirable. This mechanism is provided in the input language by a class of grid objects which, in addition to data, includes user-supplied functions to be called by the system when the grid is created, refined, and unrefined. These functions interface with grid management functions through a set of predefined external variables. The interpreter language is intended to call system and user functions that transform data (for example, that compute the solution for the next time step), execute an interactive graphics module to examine data, and perform simple arithmetic calculations. All the sophisticated processing and work involving a large number of operations should be carried out by compiled modules. The problem parameters and variables that are expected to be changed and/or examined interactively should be defined as interpreter variables. The data space of the interpreter is accessible both interactively and from compiled modules in C or Fortran. Program development in this environment consists of writing scripts in the interpreter input language and modules in a compiled language to be called from the interpreter. Programs are then run and geometrically debugged using the interpreter. The same interpreter can be used by application area specialists for debugging data and examining solutions. 4.2. The interactive graphics module The interactive graphics module has functions for displaying grid vectors and interactively examining spatial data with the help of a pointing device. The system maintains a list of grid vectors that is known to the interactive graphics module. Each grid vector has a set of attributes (such as the color palette and the array of isolevels) that define its visual representation by the interactive graphics module. These attributes are accessible through the interpreter. Interactive examination of the spatial data is performed with one of the grid vectors and/or the grid displayed in the graphics window. During this process the grid vector values at a location specified by the pointing device are printed in a pop-up window. If needed, these values can then be modified by overtyping them in the same pop-up window. Space functions, which are user-extensible, are also used for interactive examination of data. For example, a space function can be made to compute the numerical derivative of a grid vector at a location specified by the pointing device. Grid vector cross-sections are another visualization tool available to the user. A cross-section is the restriction of a function of two variables over the straight line interval connecting two points. The interval defining a cross-section is specified interactively with the pointing device. The scaling and the graph colour are among the grid vector attributes accessible through the interpreter. Several grid vectors can be displayed on one cross-section graph, which helps to visualize their relationship to each other.
5. Application
to flame-flow
interaction
problems
We applied the methods and programming tools described above to develop computational models that were used to investigate Bunsen burner flame tip opening and flame propagation
180
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
Fig. 11. Contours of temperature distribution for the problem of flame propagation in a periodic flow field with an adaptive grid superimposed over the computational domain. Computational domain dimensions are 0 < x < 32, 0 6 y < 10, Lewis number Le = 1.5, maximum amplitude of periodic flow A = 1.0, computed flame front velocity u = 0.73.
in a periodic flow field. The Bunsen burner model and numerical experiments with it are described in [9]. For theoretical background on the Bunsen burner flame tip opening problem see [ll]. In this section we outline the implementation of the model of premixed flame propagation in a periodic flow field [2]. This problem is intended to simulate certain basic features of turbulent combustion. A flame propagates in a two-dimensional infinite channel with a given unidirectional periodic field of velocities of a combustible mixture. A single one-step chemical reaction between fuel and oxidizer is assumed. We consider a conventional adiabatic, constant density reaction model in a nondimensionalized formulation similar to the one in [ll]. The special property of this problem is that we have to move the computational domain in order to keep the flame front inside it. The average flame front velocity is not known in advance and has to be computed numerically. The goal of the numerical experiments is to determine the dependence of the flame front velocity on the amplitude of the flow field and on the Lewis number. We discretize the problem in time by Euler’s method. The resulting equations are then discretized in space using the finite element method with an upstream correction. The refinement criterion for a terminal grid cell is computed as the maximum of the average numerical second derivatives of the temperature and the concentration at the cell’s nodes multiplied by the square of the cell linear size. The decision to refine or unrefine is made by comparing the refinement criteria with the prescribed refinement and unrefinement tolerances. A cell is refined if its refinement criterion is greater than the prescribed refinement tolerance. A supercell is unrefined if the refinement criteria for all its cells is less than the unrefinement tolerance. Typically, the unrefinement tolerance is chosen to be about one quarter of the refinement tolerance. An example of an automatically generated grid and typical results of the calculations are shown in Fig. 11.
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
181
Our experience confirmed that geometric debugging allows applications using complex grids to be developed relatively quickly compared with conventional methods of program development. Our adaptive grid method created very tight grids to resolve curved flame fronts with complex geometrical shapes. By comparison, a uniform grid with the same resolution as the maximal resolution of the adaptive grid shown in Fig. 11 would have about 20 times more nodes than the adaptive grid.
6. Acknowledgements The application of the interactive programming environment to flame-flow interaction problems was done in collaboration with Professor Gregory Sivashinsky who also provided much of the inspiration for the rest of the project. The work benefited from many valuable discussions with Professor Octavia Betancourt. This research was supported by the U.S. Department of Energy under Grant No. DE-FG02-88ER13822 and by the National Science Foundation under Grant No. CTS-9213414.
7. References Dynamic mesh adaption for unsteady nonlinear [l] F. Benkhaldoun, P. Leyland and B. Larrouturou, phenomena-application to flame propagation, in: Proceedings Second International Conference on Numerical Grid Generation in Computational Fluid Dynamics, Miami Beach, FL (1988). [2] H. Berestycki and G.I. Sivashinsky, Flame extinction by periodic flow field, SIAM J. Appl. Math. 51 (1991) 344-350. [3] M.J. Berger and J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, J. Cornput. Phys. 53 (1984) 484-512. [4] J.F. Dannenhofer and J.R. Baron, Robust grid adaptation for complex transonic flows, in: ALAA 24th Aerospace Sciences Meeting Proceedings, Paper AIAA-86-0495, Reno, NV (1986). [S] R.E. Ewing, Adaptive mesh refinements in large-scale fluid flow simulation, in: I. BabuSka, et al., eds., Accuracy Estimates and Adaptive Refinements in Finite Element Computations (Wiley, New York, 1986) 299-314. [6] J.M. Hyman and B. Larrouturou, On the use of adaptive moving grid methods in combustion problems, in: Proceedings Second Workshop on Modelling of Chemical Reaction Systems, Heidelberg (1986). 171 B.W. Kernighan and D.M. Ritchie, The C Programming Language (Prentice Hall, Englewood Cliffs, NJ, 1988). [8] G. Kozlovsky, An interactive programming environment for solving partial differential equations using adaptive grids (with application to flame-flow interaction problems), Tech. Report #CCNY-CS1291, Computer Science Department, City College of New York (1991). [9] G. Kozlovsky and G. Sivashinsky, Structure of the Bunsen flame: a numerical study, Combust. Sci. Tech. (submitted). [lo] W.C. Rheinboldt and C.K. Mesztenyi, On a data structure for adaptive finite element mesh refinements, ACM Trans. Math. Software 6 (2) (1980) 166-187. [ll] G.I. Sivashinsky, Structure of Bunsen flames, J. Chem. Phys. 62 (197.5) 638-643. [12] G. Strang and G.J. Fix, An Analysis of the Finite Element Method (Prentice-Hall, Englewood Cliffs, NJ, 1973). [13] O.C. Zienkiewicz, The Finite Element Method (McGraw-Hill, London, 1977).
aA
ELSEVIER
APPLIED NUMERICAL MATHEMATICS Applied
Numerical
Mathematics
14 (1994) 165-181
Solving partial differential equations using recursive grids Gregory
Kozlovsky
*
Department of Computer Science, City College of New York, New York, NY 10031, USA
Abstract In this paper we discuss numerical methods and software tools for solving partial differential equations with variable local resolution. The emphasis of our work is on finding approaches which will greatly simplify developing applications using dynamic adaptive grids. First, we investigate a class of adaptive grids that looks particularly promising for use in general-purpose solvers and propose a discretization method suitable for this class of grids. Second, we describe an interactive programming environment which can be used, among other things, as a geometrically-oriented debugger for developing applications using adaptive grids. Finally, we briefly illustrate an application of the above-mentioned methods and software to flame-flow interaction problems.
1. Introduction
When solving partial differential equations, it is often necessary to resolve small areas with high activity. ’ Blame fronts are an example of such areas. High local activity can also be found in boundary layers. It may be not computationally feasible to use grids which are uniformly spaced throughout the computational domain to resolve these small high activity areas. One increasingly common solution in this situation is the use of adaptive non-uniform grids [1,3-6,101. Adaptive non-uniform grids have different densities in different parts of the computational domain, as dictated by requirements of the accuracy of the numerical method. Because of the potential that adaptive gridding has for reducing computational costs, it is one of the forefront areas in computational physics. Many types of non-uniform grids are used in numerical computing. In general, the choice of a grid type represents a trade-off between simplicity of the grid and its ability to cover the computational domain economically with the locally required resolution. In deciding which type of grid to employ for a given problem, much depends on the geometrical configuration of the area(s) in need of high resolution. In many cases, relatively simple grids, such as tensor-product * Current address: CN Division, CERN, 1211 Geneva 23, Switzerland. E-mail: [email protected]. ’ What exactly constitutes high activity is problem-dependent. It may mean a high gradient, a large second of the solution, or something else. 016%9274/94/$07.00 0 1994 - Elsevier SSDI 0168-9274(93)E0069-E
Science
B.V. All rights reserved
derivative
166
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
grids or composite uniform grids, can do the job efficiently. We will describe, without the pretence of being comprehensive or rigorous, several typical situations where relatively uncomplicated grids can be used. For simplicity, we restrict our exposition to the two-dimensional case. l
l
l
The high activity area is centered around a singular point. A radial grid is formed by the intersection of circles centered at the point with rays emanating from it. A radial grid is actually a tensor-product grid and, as such, can be implemented using simple data structures. The density of such a grid naturally increases with proximity to the central point. The high activity area is concentrated around a line. Non-uniform tensor-product grids are well suited to this case. By suitable mapping of the computational domain the case of a simple curve can be accommodated. The high activity area can be covered by a rectangle of the size comparable to that of the high activity area. Composite uniform grids can be employed as well as non-uniform tensor-
product grids. The advantages of the former should be weighed against the increase in the complexity of computer implementation. For many important problems, however, the shape of the high activity area is more complicated than in the cases outlined above. In addition, in order to be efficient the above approaches should rely on certain geometrical shapes of the area to be refined, and a “black-box” solver must incorporate some sort of pattern recognition ability, a significant complication to be avoided. Another major consideration in grid design is the ability of the grid to change with the solution. In the simplest application of non-uniform gridding, a variably spaced grid is established initially and then held fixed through the course of the calculation. This approach to grid generation is used in many industrial finite element packages. It requires a priori knowledge of the general behaviour of the solution (usually an intuition of an experienced engineer) and may require many months of labour. Dynamic adaptive gridding methods, on the other hand, change the grid during the solution process according to the accumulated information. These methods refine the grid in the high activity areas, leave the grid coarse in the areas of little activity, and unrefine the grid in the areas where refinement is no longer needed, Automatic mesh refinement and unrefinement is absolutely necessary in case of time-dependent problems with moving areas of high activity, or when the location of the areas requiring high resolution is not known in advance. For highly nonlinear problems a large number of time steps has to be performed. A desirable property of a dynamic grid in this case is that a small change in the solution causes a small change in the grid. Tensor-product grids and, in general, composite uniform grids do not possess this property (at least when uniform subgrids are large). In this paper we investigate a class of recursive grids we call consistent recursive grids that looks particularly promising for use in general-purpose solvers. One major drawback of recursive grids is that they are more complex than tensor-product grids and, as a result, require complex data structures. In addition, numerical methods have to be reformulated to cope with the complexity of the local geometrical structure of recursive grids. Consequently, programming and, especially, debugging of applications using recursive grids become more complex. During
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
167
our research, we found out that making use of such grids possible and practical involves a whole cluster of interconnected multi-disciplinary questions including data structures, discretization methods, visualization, and interactive user interfaces. Developing PDE solvers using conventional compiled languages such as Fortran or C is time-consuming and difficult. One of the main stumbling blocks is the lack of tools for visualizing data and solutions during the development process. Engineering packages (most commonly used in structural mechanics) do have integrated interactive graphics modules which are extensively used for data debugging. Clearly, there is a similar need for interactive visualization during the process of developing new PDE solvers. As we move from relatively simple tensor-product grids to recursive grids, the need for visualization tools integrated with the program development environment becomes imperative. As a solution for this problem we propose an interactive programming environment [8] for developing adaptive grid solvers. The environment is intended to reduce drastically the amount of technical, non-problem-specific programming and can be used, among other things, as a geometrically-oriented debugger for developing applications using adaptive grids. It is currently implemented on a PC-AT and on an IBM RS-6000 workstation. In designing this environment, we aspired to provide a set of tools for the researcher rather than to create an intelligent system which “understands” PDEs and adaptive grids. A substantial effort in creating higher-level packages may thereby be saved by implementing them on top of properly designed lower-level tools. In the following sections we describe in detail consistent recursive grids, their properties, algorithms for grid handling, and data structures for grid implementation. Discretization methods on consistent recursive grids and interactive software tools are discussed only very briefly. For more details the reader is referred to [8]. In the last section, we illustrate an application of the interactive programming environment to flame-flow interaction problems.
2. Adaptive grids In selecting the type of adaptive grids suitable for a general-purpose following governing principles:
solver, we adopted the
(a) Keep the grid geometry simple and tackle the geometric complexity of the computational domain with a combination of grid refinement and special boundary finite elements. (b) Create grids which are as “tight” (economical) as possible even if it complicates programming and makes parallelization more difficult. If the library of grid management functions is reusable, the complexity of these functions is less important than in the case of ad hoc written applications. The question of parallelization can be dealt with later and, in any case, it is only one factor in the overall performance of a method. From a performance point of view, we looked for grids that could handle areas with complex geometry, have the ability to move refined area(s), be flexible enough to be used in a general-purpose environment, and allow for a high degree of local refinement. The following exposition is carried out for the two-dimensional case. The generalization for the three-dimensional case is straightforward.
168
G. Kozlovsky /Applied
Numerical Mathematics
14 (1994) 165-181
Fig. 1. Grid fragment with three grid levels. Higher grid levels are placed above the lower ones. The resulting three-dimensional picture is shown in oblique projection. The boundaries of the supercells and the lines dividing each supercell into four cells are included. Nodes are marked by black dots. Vertical lines connect co-located nodes on adjacent levels.
2.1. Definitions and terminology The main unit for building our grids is a supercell. Geometrically, a supercell is a rectangle divided by lines parallel to its sides into four equal subrectangles. Each of these subrectangles we will call a grid cell or computational cell. Grid nodes are formed by the vertices of the cells. Nine nodes are associated with each supercell. A grid starts its life as a root grid, which consists of a number of non-overlapping supercells of equal size completely covering a rectangle. ’ An elementary grid refinement is performed by creating a supercell exactly covering an existing cell (refining the cell). The refined cell is called the parent cell of the new supercell, which, in turn, is called its descendant. Refinement can later be reversed by destroying a previously created supercell (unrefining the cell). Each supercell has a grid level number associated with it and its cells. The grid level of a root supercell is zero. The grid level of a non-root supercell is equal to the level of its parental cell plus one. Nodes associated with cells on different grid levels and having identical coordinates are considered distinct and are called co-located nodes. A set of nodes co-located with a given node includes the node itself. There is only one node with the given coordinates on any grid level and, therefore, a node may be associated with more than one supercell on the same level. The grid level of a node is equal to the grid level of its associated supercells. A grid fragment with three grid levels is shown in Fig. 1. A terminal cell is a cell which does not have the descendant supercell. A terminal supercell has all four of its subcells terminal. An active supercell contains at least one terminal cell. A terminal node is a node which is not co-located with a node on a higher level.
a We are not concerned with more elaborate root grids, because a root grid is inexpensive computational domain can be economically covered using appropriate grid refinements.
and a non-rectangular
G. Kozlo~dy
/Applied
Numerical Mathematics
14 (1994) 165-181
169
It will be useful to classify the nodes by their geometrical location in an associated supercell: . l x-node: node in the center of a supercell, l c-node: node in one of the corners of a supercell, l t-node: node on a side of a supercell which is associated exclusively with this supercell (no adjacent supercell on the same grid level), l s-node: located as t-node but is associated with two adjacent supercells. Node types are illustrated in Fig. 2. Another important distinction is between internal cells and nodes, and border cells and nodes which are the cells adjacent to and the nodes located on the border of the computational domain. Internal nodes which are not t-nodes we will call r-nodes, where r stands for regular. The difference between x-nodes, c-nodes, s-nodes, and t-nodes is important from the point of view of grid data structures and basic grid algorithms. For numerical algorithms x-nodes, c-nodes, and s-nodes are all the same, except that, as we will see later, there can be more variety in the configuration of the grid around c-nodes. However, numerical methods have to be formulated differently on the r-nodes and on the internal t-nodes. Two grid cells are called neighbows or adjacent cells if they share at least one pair of co-located nodes. A pair of grid nodes are neighbours on level k if both of them have co-located nodes which are associated with the same cell of level k. Two terminal nodes are called neighbows if both nodes are co-located with nodes associated with a terminal cell. A grid is called consistent if the linear size of any adjacent terminal grid cells differs by at most a factor of two (see Fig. 31. This requirement, which is referred to as consistency, results in a certain regularity of the grid-namely, that the local structure of the grid is limited to few possible configurations. Such regularity proves advantageous from several points of view: it facilitates grid data bookkeeping in basic grid algorithms (refinement, unrefinement, grid traversal, and nearest node search), simplifies derivation and programming of numerical methods on the grid, and allows rows of stiffness matrices to be preassembled. A grid can be represented as a forest of quad-trees, with supercells as tree nodes. Unlike supercells, the grid nodes does not naturally form a tree, this being the reason why supercells are chosen as the basic building block of the grid. The tree structure of the supercells allows for certain important operations to be performed fast. Therefore, even when we use only nodes for our numerical methods, we have to keep both supercells and nodes in our data structures.
c-node
t-node
c-node
c-node
t-node
c-node
Fig. 2. Supercell
and node types.
170
G. Ko.zlousky /Applied
Numerical Mathematics
14 (1994) 165-181
Fig. 3. Consistency.
We will frequently need to refer to the relative location of neighbouring grid components (nodes, cells, and cell interfaces). We say (see Fig. 4) that nodes are located in the directions x-, x+, y-, and yf, and in the orientations x-y-, x+y-, x-y+, and xfy+ relative to the central node. It is convenient to have different generic names for directions and orientations because of their different geometrical position relative to the central node and the different role in the grid data structures. The same notation will be used for the relative positions of cells, for positions of cells relative to nodes, and vice versa. Now we can precisely define the local geometrical configuration. The geometrical configuration around a terminal node is a correspondence between the orientations and the relative levels of the terminal cells adjacent to the node, where the relative level of a cell is an integer obtained by subtracting the central node grid level from the cell grid level. A geometrical configuration is identified by its configuration number. 2.2. Properties of consistent grids Consistent definitions.
grids have the following properties,
x-y+
2-
x
most of which result directly from the
Y+
x+y+
Y-
x+y-
+ixt
Y
Fig. 4. Location
of neighbouring
nodes relative
to the central
node.
G. Kozlovshy /Applied
Numerical Mathematics
Property 2.1. There is a one-to-one correspondence x-nodes.
14 (1994) 165-181
171
between the supercells and their associated
Property 2.2. The terminal cells surrounding a terminal node are on at most two different grid levels-the level of the node and one level below it. At least one of these cells is associated with the node and is, therefore, on the same grid level. Moreover, all four terminal cells surrounding a terminal x-node or a terminal s-node are on the same level. A terminal internal t-node is surrounded by three cell-two cells associated with it and one on the level below. Only a terminal c-node can have several different grid configurations around it. Proof. Follows from the definition
of consistency.
•I
Property 2.3. A node associated with a terminal cell is either terminal or has a co-located terminal node on the next grid level. Proof. Consider the terminal node pt co-located with the given node p. According to Property
2.2 the terminal cells surrounding a terminal node can be on at most two different levels. Therefore, the highest level cell around pt can be at most one level above the level of the given terminal cell. The same is true for the levels of nodes p and pt which are associated with the aforementioned cells. 0 Property 2.4. Grid levels of two neighbouring terminal nodes differ at most by one. Proof. By definition two neighbouring terminal nodes are co-located with nodes of a terminal cell. Therefore, according to Property 2.3, the grid level of each node is either equal to the level of the cell, or is one level above it. 0 Property 2.5. An internal t-node of a consistent grid is always terminal. Proof. Given a t-node p on level 1, suppose there is a node fi on level I+ 1 co-located with node p. By the definition of an internal t-node, there is a terminal cell on level I- 1 adjacent to node p. On the other hand, there is a cell of level 1 + 1, associated with node 5, which is also adjacent to node p. Thus, the starting assumption would violate the grid consistency. •I
The following simple theorem expresses a fundamental property of consistent recursive grids, their local semi-regularity. This property distinguishes consistent recursive grids from less structured grids which can have infinitely many local configurations. Such semi-regularity proves advantageous from several points of view: it facilitates grid data bookkeeping, simplifies derivation and programming of numerical methods on the grid, and allows rows of stiffness and mass matrices to be preassembled. Theorem 2.6. The number of possible distinct geometrical configurations around an internal terminal r-node is fifteen (see Fig. 5).
G. Kozlovsky /Applied Numerical Mathematics I4 (1994) 165-181
172
Fig. 5. Nontrivial rotation.
geometrical
configurations
around
a grid node. The remaining
configurations
can be obtained
by
Proof. It follows from Property 2.2 that terminal cells of only two different sizes can surround
a terminal node and that the case where all the surrounding cells are on the level below the node level is not feasible. The theorem follows from the fact that there are four terminal cells around an internal terminal r-node. LI
2.3. Grid data structures The design of grid data structures must be a compromise between the amount of memory required and the support for the efficient execution of certain basic operations on the grids. Given the current state of computer hardware and the type of problems we deal with, speed rather than memory is the bottleneck. Therefore, in designing the data structure we gave preference to achieving higher speed over saving memory. The following is a list of basic grid operations which must be efficiently supported by grid data structures. Grid refinement/unrefinement.
Given a consistent recursive grid, a set of cells to be refined, and a set of supercells allowed to be unrefined, find a modified consistent grid in which all the necessary cell refinements are performed and as many supercells as possible are unrefined. Grid traversal. Traverse subsets of nodes in different order. Subsets of nodes can be defined by node types (such as t-nodes or r-nodes), grid levels, and node terminality. Among frequently used node orderings are the x-y, y-x orderings, and the nested dissection ordering. Nearest neighbours search. Given a node, find its neighbouring nodes and/or neighbouring cells. Given a cell, find its neighbouring cells and/or associated nodes. Point location relative to the grid. Given a point in the computational domain, find the node at the given grid level which is the closest to the point and the cell which contains the point. Support for numerical methods using hierarchical structure of the grid. 3 Given a node, find the co-located nodes on the adjacent grid levels. This operation is required, for example, in order to transfer grid vector values between different grid levels. Two arrays are used to store the information about the nodes, the supercells, and their interconnections. The array c t r e e keeps the instances of the structure 4 SUPER C E L L, and the 3 Such as the multigrid method and the nested dissection 4 The adaptive grid library is implemented in C.
method.
173
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
typedef struct supercell ( float ctr[2]; /* gn_index ctrnode; /* sc_index parcell; /* int load; /* int maxlevel; /* char pos; /* char level; /* ) SUPERCELL; typedef struct gnode c gn_index nbr[41; gn_index fine; union { gn_index coarse; sc_index scell; 1 u; unsigned char gtype; unsigned char gdir; ) GNODE;
/********** supercell structure *********/ coordinates of the super-cell center */ central node of the supercell */ index of the parent supercell */ work load associated with cell sub-tree */ maximum level of descendants *c/ position in the parent supercell */ supercell's level (root level 1s 0) */
/********* grid node structure ****f****/ /* indexes of 4 neighbouring nodes */ /* co-located next level node */ /* co-located previous level node /* associated supercell for x-node
*/ */
/* node type and inconsistency bit */ /* encoded locations of adjacent cells */
Fig. 6. Node and supercell
structures.
array g n o d e s contains the instances of the structure G N 0 D E. Each instance of these structures is identified by its index in the corresponding array. The dynamic allocation module keeps track of the free slots in c t r e e and g n o d e s. When a new supercell or node is created, a slot for it is allocated; when a supercell or a node is deleted its slot is returned to the list of the free slots. The definitions of node and supercell structures are listed in Fig. 6 with some comments. Supercell structure has references to its parent, to its children, and to its central node. Node structure has references to its immediate neighbours, to the co-located node of the parental supercell (the central node of a supercell has a reference to its associated supercell instead), and to the co-located node of a child supercell (if any). In addition, node structure contains the encoded local grid configuration and the encoded node location in its associated supercells. 2.4. Grid creation, refinement, and unrefinement Grid adaptation is done in steps which may be interspersed with numerical computations. One step of grid adaptation is a solution to the following problem. Definition 2.7 (Grid adaptation problem). Given a grid .?9’,a set of its cells 9Z’to be refined, and
a set of its supercells Z/ which are allowed to be deleted, find a minimal consistent grid in which all the cells of 9 are refined and all the supercells not belonging to ?Y are preserved. A minimal grid is a grid from which no supercells can be deleted without violating one of the above conditions. It is easy to see that because of the consistency requirement the minimal grid may contain some supercells of % marked for deletion, and some cells not belonging to 9 will be refined.
174
G. Kozlousky /Applied
Numerical Mathematics
I4 (1994) 165-181
refine(s) { For each neighbouring
cell s,, do {
If s, is one level below s, ref ine(s,) 1 create(s) 1 Fig. 7. One-cell
consistent
refinement
algorithm.
Theorem 2.8. For the given grid 22 and the given sets 2 and 2Y there exists a unique minimal grid. Proof. Suppose there are two different minimal grids ZYIand P2. Without restricting generality, we can presume that the set 2YI\PY2 is not empty and that the supercell s’ is among its supercells on the highest grid level. Then, the parent cell of S does not belong to the set 9, otherwise it would have been refined in g2. For the same reason, the supercell s’ is not among the supercells of Z? which cannot be deleted. Finally, the deletion of s’ will not cause inconsistency in ZYr,because it is absent in .Y2 and there are no supercells in .YI\3’2 which are above the grid level of s”. Therefore, contrary to the initial assumption, gI is not a minimal grid. 0
Let c r e a t e(s) be the function which creates the supercell covering the cell s. Then, the recursive function r e f i n e(s) presented in Fig. 7 refines the terminal cell s, preserving the grid consistency. Moreover, the grid remains consistent during the intermediate steps of the algorithm, because a supercell is created only after all of the neighbours on the grid level below its level have been created. Let de s t r o y(s) be the function which deletes the supercell covering the cell s. Then, the recursive function un r e f i n e(s) presented in Fig. 8 unrefines the cell s if possible without violating the grid consistency. Because of consistency, the solution to the grid adaptation problem is global in the sense that refinement or unrefinement at a particular location of the grid may depend on the situation at a different location. A naive approach to the problem can lead to repeated unrefinements and refinements of the same cell. This would not only require unnecessary
unrefine(s) boolean nounrefine = FALSE If s is non-terminal, nounrefine = TRUE For each neighbouring cell s,, do { If s, is more than one level above s, nounrefine 1 If nounref
ine is FALSE, destroy(s)
Fig. 8. One-supercell
consistent
unrefinement
algorithm.
= TRUE
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
175
adapt (1 For 1 = 0,. l,,,, do { For each supercell s E 721, refine(s) For 1 = l,,,, . .O, do { For each supercell s E L/l, unrefine(s)
Fig. 9. Optimal
grid adaptation
algorithm.
computing, but might also lead to the loss of numerical data associated with temporarily deleted nodes and cells. We propose an algorithm presented in. Fig. 9 which guarantees that this does not happen. The following notations are used: 1 is the current grid level, I,,, is the maximal grid level, &J?[and ??J1are the subsets of the sets 32 and Z, respectively, consisting of supercells on the grid level 1. Theorem 2.9. The grid adaptation algorithm is optimal in the sense that the resulting grid is minimal and that no cell refined at an intermediate stage of the algorithm is later unrefined. Proof. Let c’ be the highest level cell which was refined
and then unrefined. Because it was unrefined, c’ cannot belong to 32, and, therefore, it was refined because there was a neighbouring cell E+ on the grid level above that had to be refined. By the definition of 2, c’+ cannot be among the cells later unrefined. Hence, c’ too could not have been unrefined without violating consistency. Suppose now that the grid .?Yaresulting from the application of the grid adaptation algorithm is not minimal. Let c’ be a highest level cell which is refined in 27a but not in the minimal grid. Then, 2 does not belong to the set 33”.Therefore, c’ either was created because there was a neighbouring cell c’+ on the grid level above its level that had to be refined, or c’ belongs to the set 22 of the original grid. In the first case, by the definition of E, E+ must be refined in the minimal grid, and, therefore, c’ must also be refined in the minimal grid, which contradicts the definition of c’. In the second case, nothing prevents c’ from being deleted during the unrefinement phase of the algorithm. q 2.5. Grid traversal and subsets
Parameters of a numerical method can be node based, cell based or cell interface based. We will use the word grid component to refer to any of these three types. Most often numerical methods work with node based parameters. A grid is needed to support operations performed on data based at grid nodes, cells, or interfaces. Programs performing these operations must have the ability to traverse the grid. The order of grid traversal can be important. Grid components can be ordered by grid level, by geometrical order, by component type (such as t-nodes and r-nodes), or by any combination of the above orders.
176
G. Kozlovsky /Applied
Numerical Mathematics
14 (1994) 165-181
Two sets of tools are provided for grid traversal-the basic grid traversal functions and index lists. Basic grid traversal functions (BGTFs) traverse the grid using its internal tree structure. Index lists require some additional memory but are much more flexible than BGTFs. They allow a set of grid components to be split into any number of subsets and maintain a specific order inside every subset. Several standard splittings are provided with the system. Users can create new splittings and orderings. 2.6. Search for neighbours The grid data structure contains indexes of the neighbouring nodes on the same grid level for every direction. It also contains indexes of the co-located nodes on the adjacent levels. Therefore, all the neighbouring nodes can be found in O(1) operations. The connections between nodes and supercells are made through x-nodes. Because of this, the connections between nodes can be used to find neighbouring cells as well. 2.7. Adaptive grid library Using the data structure described, algorithms for refinement, unrefinement, grid traversal, and nearest node search for consistent recursive grids, together with necessary utilities, are implemented as a library of C functions for adaptive grid handling. The library is integrated with the programming environment and with the interactive graphics which are outlined in Section 4. At present the implementation is limited to two-dimensional grids.
3. The finite element method on recursive grids
The finite element method [12,13] is an especially attractive spatial discretization method on recursive grids because the derivation of the discrete equations on the interface between cells of different sizes is somewhat more natural than in the finite difference method, and because an h-p adaptation can potentially be employed. The treatment of the interface between computational cells of different sizes is a major problem in using recursive grids. Many authors use interpolated solution from the course grid as a boundary condition for patches of fine grid. We discretize a differential operator at all the r-nodes, including the ones on the interface between computational cells of different sizes. For such interface nodes, standard formulation of the finite element method can lead to undesirable increase in number of nodes in the discretization stencil. To solve this problem, we developed a version of the finite element method for consistent recursive grids. Our version uses discontinuous piecewise bilinear basis functions to simplify matrix generation on interfaces of computational cells of different sizes. The solution remains continuous due to additional constraints imposed on the unknowns. A particular version of the finite element method (in the weighted residual formulation) is defined by the choice of the basis functions and the weight functions.
177
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
3.1. Basis functions in space
We shall seek unknown functions in the approximate
form
(1) where +i are basis functions of space variables x and y, and a, are node based parameters to be computed. A basis function $i is a piecewise bilinear function which equals 1 at terminal node i and equals 0 at the rest of the terminal nodes. Let z be the set of elements (cells) which have node i as one of their vertices. Then 4i is non-zero only over the elements belonging to A, and is a bilinear function over each element of A. The support of the function $i, defined as the sum of the elements belonging to A, will be denoted by Si. In our case, where the elements are rectangles in two-dimensional space, the number of elements in A is four for an internal r-node and two for a t-node. Assuming that the coordinate system starts at node i, its axes are parallel to the sides of the elements, and the element j EA is located in the first quadrant, the restriction of a basis function 4i over the element j can be written as
4i=h”F’
hf-x
h;-y
I
I
(2)
where hr and hy are the x and y dimensions of the element j. The support of basis functions for several geometrical configurations of the grid is illustrated on Fig. 10. When elements of A have different sizes, or the node i is a t-node, the function 4i is apparently discontinuous. The continuity of the expansion (1) is preserved by imposing the following linear constraints on the admissible set of parameters ai ai = i(ai++ ai-)
for every internal t-node, where i+ and i- are the nodes on the ends of the cell side in the middle of which the t-node i is located. An alternative approach would be to eliminate the t-node based basis functions altogether and use continuous basis functions at the r-nodes. These functions would have larger support than basis functions c$~,which would lead to a more complicated discretization pattern and, ultimately, would also complicate programming.
Fig. 10. The support of basis functions for several geometrical marked by black disks, the remaining nodes of the discretization
configurations of the grid. The central stencil are marked by circles.
nodes
are
178
G.
Kozlovsky/Applied Numerical Mathematics 14 (1994) 165-181
3.2. Weight functions
As weight functions we take normalized linear combinations (cli=
CLi
4i
+
+
i
C
kE6
4k
of basis functions
i
where q is the set (possibly empty) of t-nodes adjacent to the node i, and pi is chosen SO that the integral of lcIiover its support Qi is unity. Note that rcIiare continuous functions of x and y. 3.3. Preassembling stiffness matrices
As was mentioned before, consistent recursive grids can have only a limited number of local geometric configurations at a node. If we use a version of the finite element method described above the number of possible discretization stencils remains equal to the number of local geometrical configurations of the grid. Consequently, rows of the stiffness matrix approximating a differential operator with constant coefficients can be preassembled for every configuration. This property leads to substantial savings of computational time.
4. The interactive
programming
environmed
The foundation of the proposed interactive programming environment [8] is an interpreter. The environment also includes an interactive graphics module and a library of functions for handling recursive adaptive grids. Integration of an interpreter with interactive graphics provides a way to connect geometrically-based human understanding with underlying data structures. This connection is important because conceptual and programming mistakes that are very difficult to find by conventional means manifest themselves through our spatial perception. We call the process of finding such mistakes geometric debugging. 4.1. The interpreter For the most part, the input language of the interpreter mimics the lexical conventions and syntax of the C programming language [7]. Most of the C fundamental data types are supported. Several specialized built-in classes and attributes have been added to the input language to facilitate work with adaptive grids. They include node- and cell-based grid vectors, grid objects, and space functions. Arrays declared as grid vectors are given special treatment by the interactive graphics module. Grid objects include data and functions to be executed upon the data every time the grid is refined or unrefined. Space functions accept a location in the computational space as their input, and output their results in a pop-up window using special output operators. Space functions are intended to be called from the interactive graphics module with a location specified by a pointing device (mouse). Let us consider grid objects in more detail. Conceptually, a grid object is an object which must be modified whenever the grid is modified. For example, when the grid is refined, the components of a grid vector corresponding to the newly created nodes have to be initialized. A
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
179
grid cell is refined either at the request of an application program or by functions of the grid management library in order to maintain the consistency of the grid. Therefore, a uniform mechanism to update grid objects is desirable. This mechanism is provided in the input language by a class of grid objects which, in addition to data, includes user-supplied functions to be called by the system when the grid is created, refined, and unrefined. These functions interface with grid management functions through a set of predefined external variables. The interpreter language is intended to call system and user functions that transform data (for example, that compute the solution for the next time step), execute an interactive graphics module to examine data, and perform simple arithmetic calculations. All the sophisticated processing and work involving a large number of operations should be carried out by compiled modules. The problem parameters and variables that are expected to be changed and/or examined interactively should be defined as interpreter variables. The data space of the interpreter is accessible both interactively and from compiled modules in C or Fortran. Program development in this environment consists of writing scripts in the interpreter input language and modules in a compiled language to be called from the interpreter. Programs are then run and geometrically debugged using the interpreter. The same interpreter can be used by application area specialists for debugging data and examining solutions. 4.2. The interactive graphics module The interactive graphics module has functions for displaying grid vectors and interactively examining spatial data with the help of a pointing device. The system maintains a list of grid vectors that is known to the interactive graphics module. Each grid vector has a set of attributes (such as the color palette and the array of isolevels) that define its visual representation by the interactive graphics module. These attributes are accessible through the interpreter. Interactive examination of the spatial data is performed with one of the grid vectors and/or the grid displayed in the graphics window. During this process the grid vector values at a location specified by the pointing device are printed in a pop-up window. If needed, these values can then be modified by overtyping them in the same pop-up window. Space functions, which are user-extensible, are also used for interactive examination of data. For example, a space function can be made to compute the numerical derivative of a grid vector at a location specified by the pointing device. Grid vector cross-sections are another visualization tool available to the user. A cross-section is the restriction of a function of two variables over the straight line interval connecting two points. The interval defining a cross-section is specified interactively with the pointing device. The scaling and the graph colour are among the grid vector attributes accessible through the interpreter. Several grid vectors can be displayed on one cross-section graph, which helps to visualize their relationship to each other.
5. Application
to flame-flow
interaction
problems
We applied the methods and programming tools described above to develop computational models that were used to investigate Bunsen burner flame tip opening and flame propagation
180
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
Fig. 11. Contours of temperature distribution for the problem of flame propagation in a periodic flow field with an adaptive grid superimposed over the computational domain. Computational domain dimensions are 0 < x < 32, 0 6 y < 10, Lewis number Le = 1.5, maximum amplitude of periodic flow A = 1.0, computed flame front velocity u = 0.73.
in a periodic flow field. The Bunsen burner model and numerical experiments with it are described in [9]. For theoretical background on the Bunsen burner flame tip opening problem see [ll]. In this section we outline the implementation of the model of premixed flame propagation in a periodic flow field [2]. This problem is intended to simulate certain basic features of turbulent combustion. A flame propagates in a two-dimensional infinite channel with a given unidirectional periodic field of velocities of a combustible mixture. A single one-step chemical reaction between fuel and oxidizer is assumed. We consider a conventional adiabatic, constant density reaction model in a nondimensionalized formulation similar to the one in [ll]. The special property of this problem is that we have to move the computational domain in order to keep the flame front inside it. The average flame front velocity is not known in advance and has to be computed numerically. The goal of the numerical experiments is to determine the dependence of the flame front velocity on the amplitude of the flow field and on the Lewis number. We discretize the problem in time by Euler’s method. The resulting equations are then discretized in space using the finite element method with an upstream correction. The refinement criterion for a terminal grid cell is computed as the maximum of the average numerical second derivatives of the temperature and the concentration at the cell’s nodes multiplied by the square of the cell linear size. The decision to refine or unrefine is made by comparing the refinement criteria with the prescribed refinement and unrefinement tolerances. A cell is refined if its refinement criterion is greater than the prescribed refinement tolerance. A supercell is unrefined if the refinement criteria for all its cells is less than the unrefinement tolerance. Typically, the unrefinement tolerance is chosen to be about one quarter of the refinement tolerance. An example of an automatically generated grid and typical results of the calculations are shown in Fig. 11.
G. Kozlovsky /Applied Numerical Mathematics 14 (1994) 165-181
181
Our experience confirmed that geometric debugging allows applications using complex grids to be developed relatively quickly compared with conventional methods of program development. Our adaptive grid method created very tight grids to resolve curved flame fronts with complex geometrical shapes. By comparison, a uniform grid with the same resolution as the maximal resolution of the adaptive grid shown in Fig. 11 would have about 20 times more nodes than the adaptive grid.
6. Acknowledgements The application of the interactive programming environment to flame-flow interaction problems was done in collaboration with Professor Gregory Sivashinsky who also provided much of the inspiration for the rest of the project. The work benefited from many valuable discussions with Professor Octavia Betancourt. This research was supported by the U.S. Department of Energy under Grant No. DE-FG02-88ER13822 and by the National Science Foundation under Grant No. CTS-9213414.
7. References Dynamic mesh adaption for unsteady nonlinear [l] F. Benkhaldoun, P. Leyland and B. Larrouturou, phenomena-application to flame propagation, in: Proceedings Second International Conference on Numerical Grid Generation in Computational Fluid Dynamics, Miami Beach, FL (1988). [2] H. Berestycki and G.I. Sivashinsky, Flame extinction by periodic flow field, SIAM J. Appl. Math. 51 (1991) 344-350. [3] M.J. Berger and J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, J. Cornput. Phys. 53 (1984) 484-512. [4] J.F. Dannenhofer and J.R. Baron, Robust grid adaptation for complex transonic flows, in: ALAA 24th Aerospace Sciences Meeting Proceedings, Paper AIAA-86-0495, Reno, NV (1986). [S] R.E. Ewing, Adaptive mesh refinements in large-scale fluid flow simulation, in: I. BabuSka, et al., eds., Accuracy Estimates and Adaptive Refinements in Finite Element Computations (Wiley, New York, 1986) 299-314. [6] J.M. Hyman and B. Larrouturou, On the use of adaptive moving grid methods in combustion problems, in: Proceedings Second Workshop on Modelling of Chemical Reaction Systems, Heidelberg (1986). 171 B.W. Kernighan and D.M. Ritchie, The C Programming Language (Prentice Hall, Englewood Cliffs, NJ, 1988). [8] G. Kozlovsky, An interactive programming environment for solving partial differential equations using adaptive grids (with application to flame-flow interaction problems), Tech. Report #CCNY-CS1291, Computer Science Department, City College of New York (1991). [9] G. Kozlovsky and G. Sivashinsky, Structure of the Bunsen flame: a numerical study, Combust. Sci. Tech. (submitted). [lo] W.C. Rheinboldt and C.K. Mesztenyi, On a data structure for adaptive finite element mesh refinements, ACM Trans. Math. Software 6 (2) (1980) 166-187. [ll] G.I. Sivashinsky, Structure of Bunsen flames, J. Chem. Phys. 62 (197.5) 638-643. [12] G. Strang and G.J. Fix, An Analysis of the Finite Element Method (Prentice-Hall, Englewood Cliffs, NJ, 1973). [13] O.C. Zienkiewicz, The Finite Element Method (McGraw-Hill, London, 1977).