Adaptation procedure

refines flow regions with large numerical errors, and coarsens regions with small truncation errors. The cell-edge data structure is applied to the flow solver and the adaptation procedure. The source code is written in Fortran, about 1,100 lines for the adaptation procedure (greatly simplified from and 500 lines for the flow solver.

Post-processing

transforms data on the cell-edge data structure to the conventional vertex-centered or cell-centered structure, then flow visualization and data processing can be realized by using commercial softwares. The VAS2D provides a few simple but quick subroutines to investigate results, by outputting iso-contours or fringes in the PostScript format; numerical results shown in this webpage (.GIF pictures) are directly transformed from the output of the VAS2D. Post-processing will not be discussed further (although it is equally important).

The source code was originally programmed in Fortran. The  adaptation module contains about 1,100 lines for the adaptation procedure (300 lines shorter than previous versions) . A typical module in the flow solver is a few hundreds lines long. For example, the module for inviscid gas modeling using the MUSCL-Hancock scheme contains 448 lines. Main modules included in the flow solver are data reconstruction module (computing spatial gradients),  viscous module  (computing viscosity & heat conduction terms),  two inviscid modules (computing convection terms using two different schemes).

The three-dimensional adaptive solver using hexahedral cells has been developed in Fortran 90.
 

Performance on Cray C90

The flow solver used in this test is the classic Lax-Wendroff scheme. The spurious oscillations which numerically appear with discontinuities are removed by an additional conservative smoothing step. The computational efficiency of all algorithms is measured by calculating a shock wave diffracting over a 90 degree corner. The test runs 400 time steps and uses about 10,000 cells (cell number can not be fixed because of adaptation). The memory requirement for the computation is nearly 2mega words.  figure  shows the CPU time required for main subroutines of the algorithm in one time step. The flow solver uses 70% of the compute time, other time is mainly spent by grid adaptation, geometry update. The geometry update only computes cell volumes and edge lengths. According to hardware statistics, the average vector length is 117.6 (the maximum machine length is 128), therefore the algorithms are well vectorized.  The average CPU time requires about 4 microseconds  per cell per step. This speed is over one hundred times faster than that conducted in a workstation using upwind schemes with either exact or approximate Riemann solvers. The efficiency is not only due to the vectorizable data structure but also due to the Euler solver, which is a central-difference scheme with the conservative smoothing.

Note: All computations are carried out in Cray C90 in single processor mode.

Cell-edge data structure

A fundamental nature of the data structure is that every cell points to its four edges, and every edge points to its two neighboring cells. The basic information saved for a cell is its location and the indices of four edges. The four edges are uniquely ordered. In Figure, edge BC, CD, DA, AB are neighboring edges NE1, NE2, NE3, NE4 respectively. Every edge is defined as a directional segment as shown in Figure. The direction assures that the left of NE1, NE4 and the right of NE2, NE3 are the cell itself, while the right of NE1, NE4 and the left of NE2, NE3 corresponds to the  right, below, above and left neighboring cells, respectively. This strict definition reduces conditional statements in the code especially in adaptation.

The basic information saved for an edge is the locations of its two ends and the indices of its left and right cells. Note that a cell never directly points to another cell but through their common edge.  For instance, if cell i need the information of its right cell j, then it should get the right index of its NE1. It is clear that one adjacent connectivity only needs two memory reads during computation,  because of the unique definition of geometry.

Adaptation procedure

The refinement and coarsening procedures are handled separately. Both procedures have similar steps for vectorization:

1. handling the inside of cells which are flagged to be refined or coarsened,
2. handling the edges of the flagged cells,
3. arranging memory.

Step 1 is conducted based on cells, and updates all inside information, such as edge deletion or inserting new cells,  which is independent of other cells. Step 2, based on edges, renews every edge of refined or coarsened cells and its two neighboring cells, which may be done without influencing other edges.

In step 1, it is required that no two neighboring cells differ by more than one refinement level. This rule prevents pathologically large volume ratios under certain circumstances. An illustration of this rule can be found in some literature. A cell to be refined must be divided into four subcells or sons as shown in the Figure. If an edge is split into two subedges, we call the edge a mother, and the subedges daughters. There are some options in the division. One is that the node inside can be optimized to get good shapes of the sons. In present computation the node is set to the center of the father cell. Another option is the locations of new nodes on the outside edges, usually at the centers. But if an edge is a boundary, AB for instance,  we choose the new location E so that edge OE is perpendicular to the boundary if this choice does not generate extremely skewed cells. If all four sons are flagged coarse, then they are deleted and recovered to their father cell, which is coarsening.

Step 1 and step 2 change the status of some cells and edges, for example from sons to fathers. Step 3 just cleans and arranges the memory, and divides all cells and all edges into two groups, father/non-father and mother/non-mother, according to their updated statues. The non-mother edges are further divided into boundary/non-boundary edges. These classifications generate an efficient data structure for a flow solver.