Multilevel Error Estimation and Adaptive h-Refinement for Cartesian Meshes with Embedded Boundaries

15 Pages·1.8 MB·English

by  Aftosmis, M. J.

AIAA 2002-0863
Multilevel Error Estimation and Adaptive
h-Refinement for Cartesian Meshes with
Embedded Boundaries
For permission to copy or republish, contact the American Institute of Aeronautics and Astronautics
1801 Alexander BelI Drive, Suite 500, Reston, VA 20191
Multilevel Error Estimation and Adaptive h-Refinement
for Cartesian Meshes with Embedded Boundaries
M. J.Aftosmis t M. J. Berger ¥
Mail Stop T-27B Courant Institute
NASA Ames Research Center 251 Mercer St.
Moffett Field, CA 94035 New York, NY 10012
[email protected] berger@ cims.nyu.edu
Abstract More recently, there has been interest in multilevel error indi-
cators which use the residual of a higher order interpolation
This paper presents the development of a mesh adaptation of the existing solution to estimate the local error. [131114][15]
module for a multilevel Cartesian solver. While the module
This approach is attractive since it can be combined with
allows mesh refinement to be driven by a variety of different
solution of an adjoint problem to produce a mesh adaptation
refinement parameters, a central feature in its design is the
incorporation of a multilevel error estimator based upon strategy that seeks to minimize error in a specific integral
direct estimates of the local truncation error using _:-extrapo- quantity of engineering interest. [13]'[16] The hope of these
lation. This error indicator exploits the fact that in regions of techniques is that, for the cost of the adjoint solution and
uniform Cartesian mesh, the spatial operator is exactly the Jacobian storage, one can produce a mesh that focuses on
same on the fine and coarse grids, and local truncation error reducing only those errors which are important to the prob-
estimates can be constructed by evaluating the residual on lem at hand, rather than simply equidistributing the error.
the coarse grid of the restricted solution from the fine grid. A
new strategy for adaptive h-refinement is also developed to While the research on adjoint-based mesh optimization is
prevent errors in smooth regions of the flow from being very exciting, two drawbacks make it currently unattractive
masked by shocks and other discontinuous features. For cer- of a general purpose analysis code. First, if the Jacobians are
tain classes of error histograms, this strategy is optimal for not already present in the solver, it is an expensive proposi-
achieving equidistribution of the refinement parameters on tion to form them solely to drive the adaptation. Secondly
hierarchical meshes, and therefore ensures grid converged and more fundamentally, many problems of interest have
solutions will be achieved for appropriately chosen refine-
competing requirements that cannot always be encapsulated
ment parameters. The robustness and accuracy of the adapta-
into a single output functional. For example, auser may wish
tion module is demonstrated using both simple model
to optimize a single simulation for lift, drag and pitching
problems and complex three dimensional examples using
meshes with from 10°to 107cells. moment. While this is currently an area of research interest,
results are not currently in-hand.
1 Introduction In this paper, we develop anew multilevel adaptation module
for use on adapted Cartesian meshes with embedded bound-
While error estimation using local gradient recovery tech-
aries. Figure 1shows a sample hierarchy of such meshes used
niques have long been popular in structural mechanics and
by a new, parallel multigrid solver. [17] the fact that these
other disciplines governed by elliptic systems, [1] such rigor-
coarse meshes can be automatically generated makes multi-
ous error estimates are generally not available in fluid
level error estimation a viable alternative, even for grids
mechanics where the governing equations contain non-self-
around very complex geometries. Seeking to take advantage
adjoint operators. As a consequence, error estimation and
of this fact, we explore Richardson extrapolation-like error
adaptation for fluid mechanics has evolved over a different
estimators which are modified to exploit the hierarchical
path. The simplest methods in the literature are gradient, and
undivided difference based featu re detectors.-[2114] Whil' e nature of adapted Cartesian grids. Even without the sensitiv-
ity information provided by adjoint-based approaches, the
these methods have been extremely successful for various
module will still offer huge savings over a priori mesh
classes of problems [2]-[7]their lack of formalism makes them
enrichment, and taking advantage of this is an important step.
difficult to blindly apply to problems far from established
The estimators we present are light-weight, and both fast to
experience. Indicators based upon estimates of interpolation
implement and execute and leverage much of the machinery
error, [7][8][9] lend some of the missing formalism. However,
already built into multigrid solvers.
since these essentially compute the local curvature of a repre-
sentative variable or combination of variables, they are not Since Cartesian grid refinement is restricted to cell subdivi-
clearly superior (or different from) straight feature detec- sion (h-refinement), the approach is intrinsically discrete.
tion. [5]The multilevel error estimators from the literature on There exists anecessary reliance upon arefinement threshold
Adaptive Mesh Refinement (AMR) were among the first which selects cells destined for subdivision. This work exam-
Richardson extrapolation-like error estimators [l°][q 1][121but ines the topic of threshold selection in detail. We develop a
have been only narrowly used outside the AMR community. methodology for robustly setting this threshold which also
ensures consistency. The approach both controls the level of
t Research Scientist, Senior Member AIAA
error in the domain and equidistributes the remainin error as
¥Professor, Courant Institute, Member AIAA
fast as possible within the restrictions of hierarchical (nested)
ThispaperisdeclaredaworkoftheU.S.Governmentandisnotsubjectto
copyrightprotectionintheUnitedStates.
American Institute of Aeronautics and Astronautics, AIAA Paper 2002-0863
J
/
/
Figure 1."Hierarchical sequence of partitioned meshes around tiffing body configuration used by parallel
flow solver of ref. [17].
refinement. Given a robust strategy for setting the adaptation satisfied the difference equation, the exact solution will not,
threshold, adaptive cases may be run fully automated. This and the difference will be the local truncation error,
reduces the sensitivity of the results to user experience/inter- LTEhA(X,t).
vention which has long been a drawback of adaptive meth-
ods. 1 R(u(x, t))
LTEh, k(x, t) = _c[u(x, t+ k) - u(x, t)] + h (3)
2 Multilevel Error Estimation If the problem reaches a steady state, then
u(x, t+ k) = u(x, t) and we can drop the dependence on
Multilevel Richardson-type error estimators like those in [12]
time.
are attractive because of their low cost and firm theoretical
underpinnings in smooth solutions. After developing these
for perfectly nesting Cartesian mesh hierarchies, this section LrEh(X ) = R(u(hx)) (4)
goes on to examine their application at adaptation boundaries
and cells cut by geometry. It goes on to present shock and
Equation 4 is instructive since it relates the local truncation
vortex-core detection strategies also planned as options to
error to the discrete residual operator R(°) and the cell dimen-
directly adapt to under-resolved features.
sion h. Furthermore is clear from the derivation that in multi-
ple dimensions LTE(_c) - Residual/hal where d is the number
2.1 Local Truncation Error Measurement
of spatial dimensions of the domain.
Consider the 1-D wave equation ut+A ux = 0 discretized
with forward Euler in time and arbitrary spatial differencing 2.2 Local Truncation Error on Uniform Cartesian
using atimestep k, and acell dimension h. Meshes
Equation 4 is useful since it permits direct measurement of
k!(UJ_+1- U_J)+ _I(FR, j-FL,j) = 0 (1) the local truncation error of a particular numerical method on
an actual computational grid for any problem that has an
is the discrete approximation to the continuous solution, exact solution. One such example is the supersonic vortex
u(x, t), in cell j at time n, while FL j and FR jare the numeri- model problem used in reference 18. Figure 2 outlines the
cal fluxes through the left and rig}_t bounda fie.' ofj and con-
problem and shows 4 sample telescoping meshes. The actual
tain the details of the spatial operator. The difference of the meshes used had from 130 to 7800 control volumes. After
numerical fluxes is the residual and the discrete equation can initialization with the exact solution, the local truncation
be written more compactly as
error was computed within each cell by dividing the residual
by the cell volume (following §2.1 and eq. (4)) using the
Euler solver in reference 17without limiters.
(_ +1_ j)+ - 0, (2)
The L1norm of the LTE in density was computed for each
where the discrete residual operator R(o) now contains all mesh and is displayed in figure 3as afunction of the number
details of the numerical flux balance for cell j. of cells in each mesh. Performing a regression analysis on the
Cartesian mesh data reveals an asymptotic slope of 2.11. For
The local truncation error, LTE, measures how well the dis-
smooth solutions, one can show that the global error is of the
crete equation models the actual PDE throughout the domain. same order 19.
It is defined by replacing the approximate solution _ with
the exact solution u(x, t) in equation (2). Since U7 exactly
-2-
American Institute ofAeronautics and Astronautics, AIAA Paper 2002-0863
2.3 Estimation of LTE on Nested Multilevel Meshes
Since the LTE can be computed for each cell in the domain, it
Y
provides an excellent basis for constructing an refinement
parameter within each cell. Unfortunately, since few prob-
lems of practical interest have an exact solution, we must
develop a method of estimating the LTE within each cell.
In §2.2, the fact that we could reasonably draw straight lines
through the data in figure 3 established the local and global
order of the method. We can express this more formally by
recalling that for a discrete method with local order p, there
ri
exists aconstant Clte such that
II '1
2_
LTEh(x ) ._Ct_e_ . (5)
p(r) = Pi 1 T i l-
Section 2.2 also established that for the solver in reference
Figure 2: Overview of supersonic vor-
17,p was around 2.11 when I[_1is the 1-norm. Furthermore,
tex model problem from ref. 18
since the RHS of equation 5vanishes for small h, the method
used to investigate local truncation
is obviously also consistent. Consistency implies that when h
error., at conditions: Min=3.0,
pin= 1/y, pin= 1,ri= 1,ro= 1.9. is small, the discrete solution Uj,hin cellj of a mesh will be a
good approximation to the exact solution u(}), and this
For comparison, figure 3 also contains the results of the same statement is the basis of our multilevel error estimation pro-
experiment performed with three types of body-fitted grids: cedure.
(1) regular, nearly unit aspect ratio quads, (2) right triangles
made by subdividing the quads, and (3) a "quality" triangula- Assume that the numerical method has converged on a fine
tion [2°] in which no angle was less than 29°. All methods grid with mesh spacing h so that the residual of the discrete
used the same numerical scheme. While the details of this solution will be zero for every cell j in the domain.
experiment are contained in the Appendix, figure 3 displays
Rh(Uj) = 0 (6)
the results. The asymptotic slopes of all four mesh types are
shown in the figure. We note that the slopes of all but the Since the method is consistent, the current discrete solution is
"quality" mesh are similar but that the level of error on the assumed to be close to the exact solution, U *,u(:_). An
j,h-
Cartesian grids is consistently 6 to 10 times lower than the estimate of the local truncation error on a coarse grid with
other grids over the full range of problem sizes. mesh spacing H can then be written by substituting the dis-
crete solution on the fine grid for the exact solution in equa-
tion 4.
" RH(4Uj'h) (7)
i0°_ ........ I- ' ' LTEH(X) - H
Following equation 7we obtain an estimate of the LTE on the
coarse grid by first restricting the discrete solution on the fine
lO'l'
grid Uj h to a coarser grid using the interpolation operator
E \"_ 1.02["-'--.__
1 " /_, and' then evaluating the discrete residual of the restricted
>,
solution on the coarse grid. Since Ujhsatisfied the numerical
scheme on the fine grid, the restricted solution I_U.- will
. n J, n
•_ 10_ not, in general, produce exactly zero residual on the coarse
grid. Thus, in the same way as the exact solution provided a
means of measuring LTE hin §2.1, the discrete solution on the
fine grid provides a method of estimating the error on the
m-3 _ C_e_ia_M._._
coarse grid. Since adaptively refined Cartesian meshes nest
_ "Quality"Triangulation _
• • RegularQuads exactly, the residual operator R is exactly the same on coarse
•-----oRegularRightTriangles and fine meshes 1.This permits us to apply equation 5 on a
cell-by-cell basis. The coarse grid estimates are then used to
[0 _ , , _ , , J,llO00 , , , , , ,,IO(X)O trigger cell refinement on fine grid. This method of error esti-
Numberofcontrolvolumes
mation is known asT-extrapolation within the multigrid corn-
Figure 3: L1norm of LTE in density for supersonic vortex
model problem on sequence of meshes for 4mesh. The 1. In ref. 12an estimate similar to eq.7 was used on curvilin-
asymptotic slope of each curve is labeled. Further ear meshes. However, oncurvilinear meshes, Rh_Rt-1
details of this experiment are presented inAppendix A. since the coarse grid control volumes are not the exact
union of the finegrid control volumes.
-3-
American Institute of Aeronautics and Astronautics, AIAA Paper 2002-0863
munity, [2U and the perfectly nested residual operators that
occur on Cartesian grids makes its use extremely attractive (a) Fine Mesh
here. I. I-/ .jI. i
In the present work, the restriction operator IH is a volume
R ---- ' ..I+_• . I . IR -R±
weighted average, which is the same as that used by the mul- /
tigrid smoother used for convergence acceleration. By using "J_ . r, . .
the prolongation operator from the multigrid scheme as well, *1 • I
we are able to construct refinement parameters on the fine oI • I
._..i. • , •
mesh for the cost of one restriction, one coarse mesh residual
_iII,-i
evaluation and one prolongation. Since these operations use
machinery which already exists in the multigrid solver, they Cells with _rregutar r,,n_ ,, i,_ i,. ,i.,,.
may be implemented with little effort. ._ne mesh s_enc_t stencil on coarse mesh
R ,_R± s
2A Mesh Irregularities (b) Coarse Mesh f
The preceding section examined error estimation on hierar- .[
chies of uniform Cartesian meshes. Adaptively refined Carte- / °b
sian meshes, however, introduce some real-world
i___o
complications which need discussion. As mentioned in the
final paragraphs of §2.3, the multilevel error estimation pro-
cedure relies upon use of the same spatial operator R on both
• R_R±
the fine and coarse meshes. We denote the spatial operator on
the pure Cartesian regions of the mesh as R z. While R L is
clearly in use in uniform regions of the grid, nearer to adapta-
tion boundaries and cut-cells the situation merits closer Figure 4: Illustration of difference stencils against adapta-
examination. tion boundary between cells of different levels. Cells in
the shaded region of the fine mesh restrict their solution
Figure 4 depicts asample fine grid (a) and the corresponding into the corresponding coarse mesh shown. Since the
coarse ceils (b) near asimple refinement boundary. The adap- residual operator on this coarse mesh is different from
tation boundary introduces irregularities into the connectivity that in the corresponding cells on the fine mesh, a cell-
graph which locally changes the residual operator. Thus, on by-cell interpretation of eq.5 is not valid and accurate
the fine mesh, (Fig. 4a) all the cells which are immediately
adjacent to the adaptation boundary have R _ R±. The situa-
stant extrapolation of the error estimates from neighboring
tion is the same on the coarse mesh (Fig. 4b) and since the cells.
form of Rh and that of RH must be identical for equation (5)
to hold, the multilevel error estimator from §2.3 can't be Inspection of figure 4 demonstrates that this irregularity
accurately constructed for the fine ceils shown shaded in fig- affects a larger portion of the mesh than one may initially
ure 4a. Since this difficulty is associated with mesh irregular- suspect since the cell next to an interface has its gradient cal-
culation affected. This cell in turn contaminates the stencil of
ity, it appears in the residual operator of cut-cells as well.
acell two away from the interface. Even a minor change like
While the formal basis for local truncation error estimation
this can be important if it occurs in a region of near zero mea-
makes it preferable to straight feature detection, there are
sured LTE. Additionally, since it is the coarse grid which pro-
several drawbacks that make them difficult to use - even for
vides the estimates, the pollution can extend as far as 4 cells
Cartesian multilevel meshes. First, for a reliable error esti-
from the interface on the fine grid, so that the fine grid cells
mate, the same stencil needs to be applied to both the coarse
have no reliable estimate of their own (cf. fig. 4). To produce
and fine grids. Thus, at mesh irregularities such as grid inter-
grids with a stencil of 5 regular cells in each direction, large
faces and cutcells, the procedure described above does not
regions of uniformly refined cells need to be generated in the
produce a reliable estimate. Colella et. a/[22] present a nice
mesh. Although x-estimates have been used with great suc-
treatment of this problem in which the estimate is multiplied
cess in block structured AMR methods, they are somewhat
by the number of cells of that type. In essence, this scales a
more difficult to implement on unstructured multilevel Carte-
cell's contribution to the error by the number affected cells.
sian. Asymptotic arguments contend that irregular cells are
While the local truncation error at an interface is larger than 2
confined to only O(N ) cells in a 3-D mesh with O(N 3) cells.
at a uniform cell of the finer or even the coarser level, related
In practice however, the difference between O(N 2) and O(N 3)
work on convergence theory for non-uniform grids suggests
is not as great on the coarse mesh as it is for the fine mesh,
that the estimate is still overly pessimistic. [23][241 (These
and even reasonably refined meshes can have as much as
results point out that its possible for O(h) LTE to still pre-
30% of the cells contaminated by irregularity. When the
serve O(h) global accuracy.) At present, we simply ignore
neighbors of cells directly affected are also counted, esti-
error estimates at interfaces and cut ceils, and rely on abuff-
mates on over half the mesh may have some contamination.
ering procedure to supply error values using piecewise con-
-4-
American Institute of Aeronautics and Astronautics, AIAA Paper 2002-0863
A second problem with LTE estimates is the need for a some-
what reliable solution to have already been computed. The rj-- j,j T
approximate solution should be in the asymptotic region for a
valid error estimate. This observation leads to a somewhat
which is a normalized version of the first difference parame-
different strategy than that followed by feature detection ter advocated by Warren et al.[26]
strategy. Feature detectors generally perform well even when
starting from a very coarse grid, making possible many Similarly, a second difference based parameter should remain
un-divided to give the same behavior. 2
cycles of adaptive refinement. For the LTE based refinement,
our strategy is to start with a fairly refined initial grid, and
on/y use two or three additional cycles of LTE driven adap- 62OJ- *J+'J:-2*J +*J-tJ2 = h} vzoj
tive refinement. In the examples here the initial grid uses rj- ddj ddj _j (9)
geometry based refinement only, but initial grids generated
from feature-detection based refinement could also be used. The 1-D refinement parameters in eqs.(8) and (9) can be
extended using a finite volume approach. This produces a
2.5 Prespeeified Adaptation Regions vector refinement parameter, with components for each of the
cell's dimensions. The first difference based refinement
Both LTE estimates and feature detection approaches suffer
parameter in the kthdirection of cellj is:
from a difficulty in localizing them to regions of interest. For
example, errors in the wake behind a wing persist for quite a
distance downstream. However, if a user is interested in the
rj, k h 2 (Vq_j "k')
pressure distribution on the wing, research suggests that there = j,k ¢j (it)
is no need to refine the wake region very much. _ In the
where I¢is the unit vector in the kthdirection.
present work we simply allow the user to define a prespeci-
fled (Cartesian) region in which adaptation is permitted. These vector refinement parameters can be used to drive both
While not particularly elegant, this approach is adequate in isotropic and anisotropic cell subdivision as discussed in
practice. ref.[4], and a similar approach may be used for the second
difference based parameters.
2.6 Feature Detection
Popular choices of ¢ are the density, velocity magnitude or
Feature detection attempts to use the current discrete solution
sometimes the local static pressure. In addition, combinations
to determine where mesh needs enhancement. Common
of these scalars are also possible. The investigations in sec-
schemes use first or second-order undivided differences and
tion 4use either density or velocity magnitude for detection.
gradient information of various flow quantities. They attempt
to "smooth out" computational space in hopes of producing a
uniform error distribution. 3 An Optimal Strategy for h-Refinement
Since the link between flow features and truncation-error is The adaptation strategy seeks to refine the mesh using the
not as formal as the Richardson-like LTE estimates discussed LTE estimates or other refinement parameters from §2 to
earlier, approaches in the literature vary significantly and improve the solution. Unlike unstructured or structured body-
everything from gradients and second derivatives to unscaled fitted approaches, r-refinement (redistribution) is not an
and scaled differences have been used. [2][3][4][51126] option with Cartesian grids, and cell sub-division (h-refine-
ment) is the only alternative. The adaptation strategy takes
In 1991 Warren et al. showed that since gradients and second the refinement parameters from section 2as input and returns
derivatives stay approximately constant with mesh refine- a boolean tag for each cell in the fine mesh. Subdividing each
ment, they make poor refinement parameters. A cell with a tagged cell will improve the accuracy of the discrete solution
high gradient will continue tohave ahigh gradient even after by reducing the local errors.
subdivision. From the standpoint of grid convergence, we
note that if the refinement parameter is to act as a substitute Algorithm A outlines the adaptation procedure.
for the real LTE it should have the same asymptotic behavior Algorithm A: Adaptation Strategy
as the LTE in smooth regions of the flow. For a pth order
Input: Vector of refinement parameters for each
scheme, this means that halving the mesh should reduce the
cell on fine mesh, r"
error by 2p.
For a first difference based quantity to mimic the LTE behav- Output: Vector of cells tagged for h-refinement "_.
ior of a second-order scheme, it must be scaled by the local A.1 Normalize: Normalize the refinement parameter vec-
mesh size, h. In one dimension, this leads to afirst difference tor. 39= Irjl/llrll=
based refinement parameter of the form:
Mesh optimization methods like those inreferences [15] 2. Eq. (9)differs from thatpresented inref.[26] (eq.15) since
and [16] localize the adaptation automatically when the theft parameter isre-scaled bythe local mesh dimension
user constructs afunctional thatdefines the singleoutput of and will therefore vanish faster than the LTE asthe mesh
interest. size isdecreased.
-5-
American Institute of Aeronautics and Astronautics, AIAA Paper 2002-0863
A.2 Tag: Apply the adaptation criteria a(.) to the normal- to achieve a different goal. Its easy to conceive of inverse-dis-
ized refinement parameter to produce a vector of tance weightings or error weightings that take the local char-
tagged cells. "_ = a(_') acteristics into account in order to identify those errors that
have the strongest impact the output functionals of interest.
A.3 Rules: Modify the set of tagged cells to ensure the
validity and smoothness of the output mesh. Figure 5 shows the histogram of refinement parameters in a
mesh which has achieved equidistribution, since all cells in
A.3.1 Buffer: Add buffer layers of tagged cells.
the domain have the same error, they fall in the same bin, and
-__ rb(_ )
the histogram looks like a delta function whose height is
A.3.2 Smoothness: Filter for island/void suppression. Ncells.
Like most strategies for adaptation, the paradigm ofAlg. A is
A.3.3 Validity: Ensure adaptation boundaries do not to start with some initial distribution of refinement parame-
exceed 2:1. _ ,---rl(-_) ters, and drive this distribution toward the idealized distribu-
tion shown in fig. 5.
A.4 Output final vector of tagged cells "_for subdivision.
Figure 6 shows the Gaussian-like distribution of refinement
The final set of tagged cells output in A.4 is then h-refined
and the solution vector is initialized on the new fine mesh parameters which serves as a model for the histogram prior to
h-refinement. A common approach found in the literature is
with the prolongation operator from the multigrid scheme.
to set the refinement threshold to some fraction of a standard
In the following sub-sections we present a new strategy for deviation above the mean of the distribution .[3][4][26]
reducing the refinement parameters to some predetermined
level. The method is optimal in that it accomplishes this task 3.2 A Fresh Look at Refinement Histograms
as fast as possible. In addition the method ensures that the
Figure 7shows an actual refinement histogram resulting from
permissible level of the refinement parameter chosen at the
a coarse grid simulation (3775 cells) of flow over an ONERA
outset will not vary as the mesh and solution evolve.
M6 wing at transonic conditions [25]. This plot bears little
resemblance to the idealized Gaussian-like model shown in
3.1 Equiflistribution of Refinement Parameters
figure 6.The refinement parameters lie between 0 and 1.The
Equidistribution aims at producing a mesh which contains mean value is 0.011 but the standard deviation 0.04. More-
the same level of refinement parameter in each cell. Since the over, fully 8.2% of the cells lie below the mean, and almost
refinement parameters are stand-ins for the local truncation 50% have 1_< 0.001. As a consequence of this extreme dis-
error, this goal spreads the remaining discretization error out parity in scales, setting the threshold, t, anyplace above the
evenly over the domain. As this level is reduced, the method mean only addresses the error in a handful of the most severe
is guaranteed to converge to the correct solution. This princi- cells. Only after the very worst errors are reduced by many
ple guards against over-resolving some features of the flow cell refinements will error in the bulk of the domain be
while overlooking others. addressed. In shocked flows, the refinement parameter will be
In practice, equidistribution is somewhat over conservative highest in cells with shocks or other strong non-linear fea-
most of the time. It assumes that all errors are equally impor- tures. As a consequence, this approach will inadvertently
tant to the simulation, and this is certainly not the case most result in arefinement process which over-resolves shocks and
of the time. However, without additional guidance about other severe features without ever addressing smooth regions
what is important for a particular simulation, equidistribution of the flow. Oversights such as these have been shown to pro-
simply ensures that everything in the simulation is equally duce an adaptive- p:[r2o6c]edure which can actually converge to
correct. In addition, if amethod can control the LTE distribu- the wrong solution .
tion to achieve equidistribution, then it can control the error
Ncells--
Mean/Median t E
_ Refinement
Refinement Parameter
Refinement Parhmeter
Figure 5: Idealized histogram of refinement parameters
Figure 6: Idealized distribution of refinement parameters
in a mesh which as achieved equidistribution.
prior to h-refinement. Region acontains Nacells.
-6-
American Institute of Aeronautics and Astronautics, AIAA Paper 2002-0863
40OO ' i , I l L
300O
et_ e_
-_200n
c
10(X)
I,
O_L r ..0,J2 _ 0.4I 0.6I , 0.8I ,
RcllL_lloent Pttraltl Refinement Parameter
Figure 9: h-refinement moves the cells in region ato the left
Figure 7" Histogram of adaptation parameter for a
coarse mesh simulation of flow around an ONERA according to the order of accuracy of the scheme. If a
M6 wing at transonic conditions [251. contains Nacells, then a* contains m Nacells, where m
is the number of children produced by refining acell.
In earlier work [4] we advocated a filtering approach which
aof figure 9, and each cell is subdivided into m children, then
removed cells containing shocks or other strong features in
a* wiU contain m Nanew cells.
an attempt to clean up the histogram prior to computing the
mean and standard deviation. Even this approach, however is If our goal is equidistribution, then we desire to build the
dubious given the huge disparity in scales. delta function of figure 5.An optimal method construct these
as rapidly as possible. Assuming that the histogram is mono-
A less ad-hoc approach toward compressing the scales in the
tonically decreasing to the right of the median, the new histo-
histogram is to simply take the log of the refinement parame-
gram grows most rapidly if the highest point of a* is placed
ter prior to binning up the histogram. Figure 8 shows the his-
on top of the median of the existing distribution. Since the
togram, of the same data as in figure 7 but computed using
cells in a move 2 units to the left, the threshold which builds
log2(l_'l) rather than simply I_'t. The mean value of this new
the highest peak is identically 2.
distribution is -6.4 and the standard deviation is 4.3. The
rescaled data much more closely resembles the idealized dis- Subsequent refinements will add to this same peak, and the
tribution in figure 5. target level of error will remain constant as the mesh evolves.
It is interesting to note that if the threshold is chosen above or
[40
below than this value the target error level will continue
migrate higher or lower (respectively) with subsequent cell
12(I -.'
refinements.
.._100
Just as refinement moves cells to the left on the histogram,
coarsening transfers ceils to the right. In the absence of
coarsening, these low-error cells will remain in the histogram
and appear as a "tall" to the left of the peak value. Figure 10
4O illustrates the evolving histogram. With the threshold chosen
as described above, newly refined ceils will not alter the his-
211 -15 -10 -5 0
log_2fRcfit_ment P,-u-ztm,)
Figure 8: Histogram of data from fig.7 computed using
log2(Irt) rather than the raw refinement parameter.
3.3 Optimal Threshold Selection
The motivation for choosing base 2for the logarithms used to
rescale the refinement parameters in the proceeding section
becomes clear when selecting an appropriate threshold. Near
grid convergence, each 2:1 cell refinement using apth-order
scheme will reduce the LTE by a factor of 2p. With 2:1 cell
refinement, and the present second-order scheme, the chil-
___3_em / Par\a__,r/ent!
dren of a h-refined cell will therefore get translated an abso-
lute distance of 2units to the left on these base-2 histograms. Figure 10: Evolution of a histogram for a mesh without
Figure 9illustrates this process. If there are Nacells in region coarsening. Growth in the number of ceils before the
peak will be approximately geometric.
-7-
American Institute of Aeronautics and Astronautics, AIAA Paper 2002-0863
togram to the left of the peak value, and therefore no newly
refined cells can ever end up in this tail. Since they were cells
in the original unadapted mesh, this tail contains only coarse
cells, and since coarse cells fill space very quickly, there can-
not be very many of them in the domain as compared to the
number of highly refined cells in the peak. In addition, since
the peak was built via cell subdivision, and the number of
children produced per cell is generally constant, the growth
approaching the peak from the left must be approximately
geometric.
3.4 Specifying an Error Level
Skeptics may point out that setting the adaptation threshold
to its optimal value removes the user's "control" of the adap- 4000
tive process. While this is precisely the goal of automation,
there is aclear need for auser to be able to have some control
over the level of error in the final solution.
In three dimensions, isotropic h-refinement of almost any
polyhedron produces cells very quickly. Hexahedra are
refined 8:1, and depending on the method of subdivision, tet-
rahedra multiply at about the same rate. Even allowing for
anisotropic subdivision of some cells, the 3-D growth rate is
0
such that adaptively refined cells virtually always account for
the median value in the histogram.
Since the location of the peak can be controlled by adjusting
t, the most efficient strategy is to set the desired peak location
on as coarse a mesh as possible (i.e. the first adaptation
cycle). Subsequent adaptations will then continue to build on
this peak using the optimal threshold. This allows the user to
set a desirable error level based upon the histogram of the
unadapted coarse mesh, and then drive the refinement hands-
free.
-t5 -10 -5 0
3.5 An Illustrative Example
Using the coarse mesh simulation from the base-2 histogram
in figure 8, figure 11 shows the evolution of the histogram in
this simulation over the next 5 adaptation cycles. This exam-
ple clearly shows the rapid growth of the peak in the histo-
gram confirming its approach toward equidistribution of the
refinement parameter.
Since this is a real 3-D transonic flow, several issues merit
discussion. Adaptation was driven by the undivided first dif-
ference of velocity magnitude scaled by the local mesh
dimension as presented in §2.6. The adaptive procedure in
Alg. A tagged ceils according to adaptation criteria, and these
tags were then modified to satisfy the smoothness and mesh
validity rules detailed in A.3. In addition, the adaptation crite-
ria used in this example was:
^ Irj]>I;I+2
a(rj) = 1 (11)
0 otherwise
-15 -10 -5 0
,'2q h)l_2(lclllll21ill;,lll I_l/anij
Where ]rl is the mean value of the magnitude of the refine-
Figure 11." Evolution of histogram through 5 adaptation
ment parameter, rather than the precise median as called for cycles for transonic ONERA M6 wing case using the
by the theoretical development earlier in this section. Of optimal threshold.
course for the narrow peaked histograms shown in the figure,
-8-
American Institute of Aeronautics and Astronautics, AIAA Paper 2002-0863
the mean is close to the median - it is within one unit of the This flow was computed using both the multilevel x-extrapo-
median at every cycle after the first. Nevertheless one would lation LTE estimates and scaled first-difference (eq.10) based
expect somewhat better performance if the true median was feature detection using _ = In as a refinement parameter.
used. Figure 12a displays both the mesh and solution resulting
from the LTE based adaptation, while 12b contains these
4 Numerical Examples and Discussion same plots using feature detection. In figure 12, the x-extrap-
olation analysis used an initial mesh with 8levels of geome-
The LTE estimates and feature detection approaches in §2 try based refinement and then an additional three cycles of
have been applied to both simple and complex configurations solution-based refinement following the philosophy of §2.4.
in three dimensions. This section begins by presenting results The final mesh shown contains 1.8M cells. The feature detec-
showing that, when combined with the h-refinement strategy tion based results in figure 12b began on a mesh created with
of §3, both can provide valid approaches to achieve grid-con- 5 levels of geometry-based adaptation and about 30k cells
verged solutions. We then present adapted solutions on two followed by 6cycles of adaptation. The final mesh contained
complex configurations to demonstrate the robustness of the 1.9M cells. After 5 levels of adaptive refinement, the mesh
procedure when run "hands-off' on real-world complex con- contained about 1.2M cells and integrated quantities (normal,
figurations. axial and side force) were virtually the same as those
obtained on the final (1.9M cell) mesh. The integrated quanti-
4.10NERA M6 Wing
fies for both examples changed by less than 0.1% in the last
In ref.[17] the baseline Euler solver was validated using the adaptation cycle suggesting that the results are grid con-
well-known ONERA M6 wing example. [25]This case con- verged.
siders the transonic flow over this wing at M_ =0.84 and
Lift and drag coefficients for z-extrapolation were: 0.3041
a = 3.06°. Although viscosity was obviously present in the
and 0.0117 while those the feature-detection were 0.3042 and
experiment, the case has been widely studied using inviscid
0.0116. Comparison between the two simulations reveals a
solvers, and amultitude of Euler solutions are available in the
difference of less than 0.04% in the magnitude of the total
literature for comparison.
force on the wing, and the final meshes have cell counts
within 6% of each other.
;C0n[6Urs /
\ / density contours
Figure 12a: Flow over an ONERA M6 wing at Moo---0.84 Figure 12b: Flow over an ONERA M6 wing at Moo=0.84
and cc=3.06 °, adaptation driven by x-estimates of and a = 3.06 °, adaptation driven using scaled first
LYE in density. The final mesh contains 1.8M cells differences of velocity magnitude. The final mesh
contains 1.9M cells
-9-