HIGH-ORDER FINITE-VOLUME METHODS ON LOCALLY-STRUCTURED GRIDS

. We present an approach to designing arbitrarily high-order ﬁnite-volume spatial discretizations on locally-rectangular grids. It is based on the use of a simple class of high-order quadratures for computing the average of ﬂuxes over faces. This approach has the advantage of being a variation on widely-used second-order methods, so that the prior experience in engineering those methods carries over in the higher-order case. Among the issues discussed are the basic design principles for uniform grids, the extension to locally-reﬁned nest grid hierarchies, and the treatment of complex geometries using mapped grids, multiblock grids, and cut-cell representations.


1.
Introduction. It is often the case that one wants to compute solutions to partial differential equations containing terms of the form ∇ · F , F = (F 1 , . . . , F D ), F d = F d (x). Classical examples include Poisson's equation and time-dependent systems of conservation laws If F = F (U ), the system is hyperbolic if ∇ U (n · F ) has real eigenvalues and a complete set of eigenvectors for any constant unit vectorn. One of the oldest and most widely-used spatial discretization methods is the finite-volume method on rectangular grids. In its simplest form, a rectangular spatial domain Ω on which the PDE is to be solved is discretized into a union of control volumes V i , i ∈ Z D V i = [ih, (i + u)h] , i ∈ Z D , u = (1, 1, ..., 1).
Averages over control volumes are denoted by · i By the divergence theorem, the average of ∇ · F over a control volume is equal to the sum of differences of the averages of the normal component of F over the faces bounding the control volume where e d is the unit vector in the d th coordinate direction. A numerical approximation to ∇ · F i is obtained by replacing the integrals of the fluxes on the boundary of the control volume with quadratures. The resulting discretization of ∇ · F satisfies a discrete form of the divergence theorem: for any union of control volumes Υ, the sum of the volume-weighted averages of ∇ · F over the control volumes are equal to the sum of the approximations of the boundary integrals.
where ± = −, + if the boundary face is on the low or high side, respectively of a control volume contained in Υ. To discretize solutions to time-dependent problems, we average (2) over control volumes with the numerical fluxes on the right-hand side computed as a function of the U , and the resulting system of ODE discretized in time. Using such a method of lines approach, the fully-discrete system can be written in the form of a single-step update where U n i ≈ U i (t n ), and F d n+ 1 2 i+ 1 2 e d is the approximation to 1 ∆t implied by the time discretization of (8).
Finite-volume methods were first introduced by Lax in [25,24] for the purpose of computing discontinuous solutions to hyperbolic conservation laws. He found that the discrete divergence property (7) leads to the discrete traveling waves that satisfy the Rankine-Hugoniot jump relations. This empirical observation was later supported by the theorem in [26]: if the discrete flux function is consistent and the discrete time evolution is given by a one-step discretization of the form (8), then the computed solution converges to a weak solution of (2), provided it converges at all.
In addition to its essential role in computing discontinuous solutions to conservation laws, the finite-volume approach is useful in at least two other settings. One is the solution to marginally-resolved problems. For example, in low-Mach number combustion problems, it is possible to obtain temperatures that significantly exceed the adiabatic flame temperature if the thermodynamic energy and chemical species are not strictly conserved [35]. The second is in solving Poisson's equation with Neumann or periodic boundary conditions, for which there is a nontrivial null space consisting of constant functions. The resulting solvability condition is manifest at the discrete level in a finite-volume discretization by applying (7) to the entire computational domain.
For more than 60 years, finite volume methods have been a mainstay of applied numerical PDE, with application to complex multiphysics problems, accurate and robust treatment of discontinuities, and the use of elaborate gridding strategies, both for representing complex geometries and for solution adaptivity. However, up until about a decade ago, they were almost exclusively second-order accurate methods, at least for multidimensional problems. In the context of the rectangular grid discretization (5), this typically corresponds to the use of the midpoint rule as the quadrature method for approximating the average of the fluxes over the faces. There are a combination of software complexity and algorithmic reasons for this. One is that, for complex multiphysics problems with multiple time scales, models for different physical processes are coupled together using operator-splitting methods that are at most second-order accurate. This permits the use of different combinations of implicit or explicit time discretization methods, depending on degree of stiffness of the underlying submodel relative to the time scales of interest. Operator splitting also allows for a software design strategy that permits straightforward insertion and removal of software components to simulate different physical processes. A second issue is that the second-order spatial discretizations are robust: they lead to M-matrices for elliptic and parabolic problems; for hyperbolic problems, there are well-established strategies for suppressing oscillations at discontinuities and underresolved regions, while remaining second-order accurate in regions where the solution is smooth. More recently, this point of view has begun to change, for a number of reasons.
• In certain specialized fields such as fundamental studies of turbulent flows high-order methods (fourth-and higher-order) have been shown to be necessary in regimes in which scientists would like to simulate more general complex problems. • There are new methods for time discretization [22,9] that are higher-order accurate, while retaining the other desirable algorithmic and software properties of operator splitting. • Finite-volume methods typically suffer a loss of one order of accuracy in truncation error at places where the mesh fails to be smooth, such as refinement boundaries in structured adaptive mesh refinement, block boundaries in multiblock methods, and near the irregular control volumes arising in cut-cell approaches to complex geometries. Dropping from second-order accuracy to first order accuracy, even on a set of codimension one, can lead to a significant increase in the error; the relative impact of the loss of one order of accuracy is less as the order of accuracy of the discretization increases. • Processors used in high-performance computers have been undergoing disruptive changes in their architectures [31,40]: the ratio of total memory capacity to aggregate floating point capability is decreasing, and the cost (in time and power consumption) of data movement is increasing relative to that of performing float point operations. Higher-order accurate methods could be an effective response to these changes. They typically perform more floating point operations per unit of data traffic, and they have the potential to produce a result with a given level of accuracy with fewer computational degrees of freedom, at nearly the same throughput as a low-order method [3]. In what follows, we describe one specific approach to designing arbitrarily highorder finite-volume spatial discretizations on locally-rectangular grids. It is based on the use of a particularly simple class of high-order quadratures for computing the average of fluxes over faces. This approach has the advantage of being a variation on widely-used second-order methods, so that the prior experience in engineering those methods carries over to the higher-order case. Some of the design issues include: • High-order accurate methods on a Cartesian grid. Given averages of the solution over control volumes, how does one compute high-order accurate averages of nonlinear functions of the solution and / or its derivatives over faces ? For time-dependent hyperbolic problems, another set of issues are the development of limiters and other dissipation mechanisms to maintain stability, particularly at discontinuities and underresolved regions, while preserving high-order accuracy for resolved smooth solutions. • Adaptive mesh refinement. Some aspects of adaptive mesh refinement on nested rectangular grids are made simpler by the fact that the averages (4), (6) on a coarse grid can be computed exactly by arithmetic averaging of the corresponding quantities defined on a nested refined grid. This allows the use of uniform-grid methods to interpolate boundary conditions at refinement boundaries by averaging the fine grid solution onto the coarse grid, while still maintaining high-order accuracy. In the case of time-dependent problems in which different time steps are used at different levels of refinement, it is essential to perform interpolation in time at refinement boundaries without having to keep the solution history at multiple time levels. • Complex geometries. Mapped grids, multiblock grids, and cut-cell representations of geometry are all widely employed to handle various kinds of geometric complexity. Extending them to higher-order accuracy involves resolving a number of new conceptual and technical issues.
2.1. Cartesian grids. We will use time-dependent conservation laws of the form (2) as a motivating example. If we discretize first in space (8), we can use a highorder ODE solution algorithm to solve the resulting system of ODE, such as an explicit or semi-implicit Runge-Kutta method. In that case, it suffices to be able to compute a high-order approximation to the flux averages on faces, given the values of U . We first consider the case F (U ) = CU , C d each a constant M × M matrix. In one space dimension, we note [16,43] that the averages over cells provide values for the indefinite integral of U evaluated at cell faces Where D is a finite-difference approximation to the spatial derivative of any order of accuracy. For example, if D is the classical fourth-order centered-difference approximation to ∂ x , then This reasoning applies to constant coefficient systems in higher dimensions without modification. If U(x d ) is the indefinite integral of the average of U across rectangular patches aligned with a row of the faces {A d We can use a similar approach to compute fluxes that are a linear combination with constant coefficients of ∇U i+ 1 2 e d . For the derivative of U in the direction normal to the face, we use and approximate ∂ 2 x d by a suitably high-order difference approximation. For the case of ∂ x d U i+ 1 2 e d , d = d, we first compute U i+ 1 2 e d using (12). We then compute ∂ x d U i+ 1 2 e d using the exact relationship ∂U ∂x d i+ 1 and approximate the averages of U over (d, d ) edges of the grid on the right-hand side by applying (12) to the face averages of U . For nonlinear problems, we have to do something more complicated, due to the fact that F (U ) = F ( U ) + O(h 2 ). However, applying a nonlinear function to U evaluated at a point has zero truncation error. So we need to be able to easily compute point values from average values, and conversely. This is easily done, to any order of accuracy. We denote by q i the value of q at the center of a control volume. Then It is tedious, but straightforward, to apply the approach given above to compute the derivatives in the right hand side to the necessary order -the process can be automated using a symbolic algebra package. These formulas apply equally well to transform between face averages and point values at the centers of faces, except that the sums over the derivatives in the correction terms exclude the derivatives with respect to the direction normal to the face over which one is averaging. A special case that often arises is the calculation of the divergence where the fluxes are the product of two spatially-varying functions, e.g. the advection operator ∇ · ( uq), u = u(x). In that case, it is useful to express the average of the product as the product of the averages plus higher order correction terms. For fourth-order accuracy, such a formula is given by where D d the two-point centered-difference approximation to ∂ d at i + 1 2 e d . There are further corrections that extend this to any order of accuracy.
To see how this works in practice, we outline how to apply this to equations in which the flux F = F (W, ∇W ), where W = W (U ) is a nonlinear function of U . This is typical of nonlinear fluid equations, both in the hyperbolic case [28], where one may want to perform a change of variables to compute limiters; or in the presence of diffusive transport, in which case the diffusive contributions to the fluxes are given in terms of gradients of the derived quantities (velocity, temperature), rather than the conserved variables (momentum, total energy). Given U i , we want to compute a fourth-order accurate approximation to F d (W, ∇W ) i+ 1 2 e d . We do so in the following steps.
where D 2 d is the three-point centered-difference approximation to ∂ 2 x d . Note that we use the average of U and W applied to the average of U to compute the O(h 2 ) correction terms in this step. This keeps the size of the stencil to a minimum, while still maintaining fourth-order accuracy. It also guarantees that, if W (U ) = U , then W = U . 2. Compute W i+ 1 2 e d , ∇W i+ 1 2 e d given W i , using the methods described above. 3. We compute the fluxes by evaluating the flux at the point value of W , then apply (20) again to obtain the average over the face. We use a similar secondorder approximation to derivatives in the O(h 2 ) correction in the final step as used in (22) to control the diameter of the stencil.
Given that U i is the exact value for the average of U over the control volume, and that U is a smooth function of x and F is a smooth function of U , the error F d i+ 1 2 e d as an approximation to the face average of the fluxes is O(h 4 ), where the coefficient multiplying the leading-order term of the error is a smooth function of space. When this is substituted into the right-hand side of (7), we see that the truncation error of the divergence is O(h 4 ), due to the usual cancellation of smoothly-varying errors on successive faces.
In the discussion above, we have focused mainly on the truncation error of these methods. With respect to stability, we can use discrete Fourier analysis to prove von Neumann stability in the constant coefficient case. For the higher-order extensions described here of the classical second-order methods for second-order elliptic operators with smooth coefficients, we obtain discretizations for which geometric multigrid converges in O(1) iterations [2]. Boundary-value problems have the same well-posedness and solvability properties as for the second-order case. For hyperbolic problems, the interaction between accuracy and stability / robustness is more delicate, and is discussed in Section 2.3 below.

Mapped coordinates.
We assume that we want to compute solutions to conservation laws on a domain that is the image of a smooth mapping from the unit cube in an abstract coordinate space into physical space [14] In that case, the divergence of fluxes in physical space can be written in terms of a divergance in the mapping space.
and N is the cofactor matrix of ∇ ξ X. Conservation laws are discretized in terms of averages of ∇ x F over images of the control volumes in the coordinate space.
where h is the mesh spacing in the coordinate space, and · i , · i+ 1 2 e d denote the averages over control volumes (4) and faces (6) in the coordinate space. Using the change of variables formula we have converted the average of ∇ x · F into the average of a divergence of fluxes over cubic control volumes in the coordinate space, so that the quadratures in the previous section applies here. In particular, we can use product formulas such as (21) to compute (N T F ) d i+ 1 2 e d to fourth-or higher-order accuracy For time-dependent problems (2), we can use the change of variables formula to obtain an evolution equation for averages over the control volumes in the coordinate If F is identically equal to a constant as a function of space, then ∇ x · F ≡ 0. We want to preserve this property under discretization. This requirement is often referred to as freestream preservation. While this is trivially satisfied for the Cartesian discretization (5), it imposes a nontrivial constraint on the calculation of N i+ 1 2 e d for the right-hand side of (31) to vanish. If we use (32) to compute the face-averaged fluxes, then the higher-order terms vanish if F is a constant, and (31) becomes a constraint on the quadrature rule used to compute N i+ 1 2 e d that, if met, is sufficient to guarantee freestream preservation. By equality of mixed partials, the We can then use Stokes' theorem to compute the average of N s over faces in terms of the integrals of N over edges where E ± d,d are the low (−) and high (+) edges of A d i+ 1 2 e d corresponding to both ξ d and ξ d equal to constants. To obtain freestream preservation and p th -order accuracy, we can compute N s i+ 1 2 e d using the right-hand side of (34) with the edge integrals replaced by any quadrature rule that is p th order accurate, provided that they satisfy the antisymmetry condition corresponding to N s d,d = −N s d ,d . To complete the construction, we note that, while N is not unique, there is a distinguished choice that is a local function of X, ∇ ξ X.
Given (34), the extension of the approach used for Cartesian grids in the previous section to the mapped-grid setting is relatively straightforward. For example, to compute the right-hand side of (33) corresponding to (22)-(26), one can use a cellcentered version of the product formula (21) to make the transformation JU i → In that case, we use the chain rule ∇ x W = J −1 N ∇ ξ W combined with the product formula (21). We then use the product formula (32) to compute the average of the fluxes at faces.
2.3. Dissipation mechanisms, limiters, and positivity preservation for hyperbolic problems. Dissipative Methods. Methods such as (14) that use stencils that are centered around the face where the value is being computed correspond to even-order, centered-difference approximations. For constant-coefficient problems, as the wavelength of a Fourier mode approaches 2h, the propagation speed of the mode appoaches zero independent of the exact wave propagation velocity for the PDE, while the amplitude of the Fourier mode remains unchanged. While these methods are still Fourier stable, the near-neutral stability and large phase error of short-wavelength modes can lead to various degrees of instability when applied to variable-coefficient and nonlinear problems. For this reason, it was suggested in [23] that such methods be modified so as to be dissipative: in the long-wavelength limit, the order of accuracy of the method remains unchanged, while reducing the amplitude of all nonconstant Fourier modes. This can be done in one of two ways. The approach proposed in [23] is to add to the flux in the d direction of a 2p th -order method a term of the form Ch 2p+1 ∂ 2p+1 This corresponds to adding a dissipation term to the righthand side of (8). An alternative approach is to use an (2p + 1) st -order method with a stencil for the flux at a face centered around the cell upwind to that face relative to the propagation speed. Such a one-point upstream-centered method has a truncation error that is instrinsically dissipative, damping all nonconstant Fourier modes, including the high-wavenumber ones. For systems of equations, one computes one-point upstream-centered face values centered on both cells adjacent to the face, and then uses an approximate Riemann solver to synthesize a upwind state relative to the different wave propagation speeds. Either approach provides a suitable mechanism for damping neutrally-stable high-wavenumber modes, at least in one dimension.
There are additional complications in more than one space dimension. For solutions that depend on only one space variable, it is possible for modes with wavenumbers perpendicular to that direction of propagation to be neutrally stable, and dissipation methods based on upwinding introduce no dissipation in the transverse directions. The complete lack of dissipation in those settings can lead to nonlinear instabilities, as observed, for example, in [41] for strong shocks aligned with one of the coordinate directions. This instability is eliminated by modifying the nonlinear diffusion term introduced near shocks so that the diffusion coefficient has the same magnitude in all of the coordinate directions. This has the effect of damping the nonlinear instability transverse to the shock, while maintaining the advertised level of accuracy in smooth regions [16,28]. This suggests that any dissipation that is introduced to stabilize short-wavelength modes should introduce damping in all coordinate directions appropriate to the desired level of accuracy. Limiters. For hyperbolic conservation laws, there is a fundamental problem to representing discontinuities as discrete traveling waves on a grid. In the neighborhood of the discontinuity, it is necessary to add a regularization term in the form of a numerical viscosity with a coefficient that is O(h) in the neighborhood of the discontinuity [32,24]. Such a term being applied everywhere leads to large errors in regions where the solution is smooth or for linearly degenerate discontinuities. This was recognized as a problem early on in the development of methods for discontinuities, and led to a variety of approaches for tuning the artificial viscosity so that it was sufficiently large at nonlinear discontinuities, and small in smooth regions. Early methods of this type required the setting of various parameters in the viscous terms that were problem-dependent. A more robust class of solutions, called flux limiters, were first introduced in the 1970's [8,39,42] to address this issue. The simplest version of this approach is the Flux-Corrected Transport (FCT) limiter in which the time-averaged flux in (9) is given by 2 e d is a time-averaged flux computed using a high-order method such as the one described above, F d L i+ 1 2 e d is a flux computed using a stable first-order accurate method, such as donor-cell differencing, or the CTU method in [12,36]. The quantity η i+ 1 2 e d ∈ [0, 1] is a hybridization coefficient: if η i+ 1 2 e d = 0, we are using the high-order accurate flux; if η i+ 1 2 e d = 1, the first-order accurate flux. To understand the idea behind this approach, consider the case of a scalar equation, i.e., U (x, t) ∈ R. We note that, by modified equation analysis, the correction term The FCT approach in [42] is to compute η to be as small as possible subject to the conditions that the solution remains within prescribed bounds, e.g. the componentwise max and min of in the neighborhood of i after the low-order method has been applied for one time step, combined with the max and min of U n in the same neighborhood. Geometric limiting [39] is an approach that is specific to the interpolation-based evaluation of the fluxes. In the semi-discrete setting considered here, we apply the limiter at each evaluation of the right-hand side. For a scalar, we replace a single-valued calculation of U i+ 1 2 e d at the face with a double-valued one where U i+ 1 2 e d is computed using (11). The single-valued flux is computed as an approximate solution to the Riemann problem for the 1D equation obtained from assuming that the solution is constant in the coordinate directions d = d is constant. In the case where the sign of the wave speed is unambiguous, this simply corresponds to choosing the upwind value of (37). The prescription in [39,16] for computing the hybridization coefficient η i+ 1 2 e d is based on assuming an interpolating polynomial derived from the U i+ 1 2 e d does not introduce any new extrema not already present in the values U i . In the case where η ≡ 1, we are replacing the high-order method for computing the flux with a first-order upwind discretization in space, which introduces adequate dissipation at discontinuities. For both geometric limiting and FCT, additional conditions and dissipation mechanisms need to be applied, e.g. to enforce entropy conditions and obtain robust treatment of strong nonlinearities.
In order to obtain high-order accuracy for smooth solutions for h sufficiently small using either of these approaches, this choice of bounds assumes that the exact U i+ 1 2 e d can be bounded by local values of U . This fails at smooth extrema, where it is possible for U i+ 1 2 e d (t) to not be within the range of the local values of U . If that occurs, η > 0 may be nonzero at some faces. This leads to a low-order update to the solution, and typically a method that has O(h) errors in the solution at smooth extrema. This was recognized in [42], and a method for modifying such limiters to eliminate this loss of accuracy was proposed there. One way of summarizing this approach is given by the following [15,28,10].
• If there is no possibility of a smooth extremum, the solution is bounded by values of U in the neighborhood of the control volume, and and we use the limiters described above.
• In regions where there is a possibility of a smooth extremum, we use a different limiter, based on comparing various approximations to the second derivatives of the solution, and applying a limiter if they differ by O(1), or change sign.
This is a restatement in algebraic terms of the extremum-preserving method in [42], in the spirit of [39]. In the latter work, the hybridization coefficient in one dimension is a function of the ratios of several finite difference approximations to first derivative. The coefficient η is nonzero if these ratios are not sufficiently close to 1 or the estimates of the first derivative do not have a consistent sign. Using constraints based on the relative sizes of first derivatives are reliable indicators of whether the solution is sufficiently well-resolved on the grid so that the use of higherorder methods can be used, except at locations where the first derivative vanishes and the solution is not discontinuous. In that case, we use the same strategy, except based on several finite-difference approximations to the second derivative. In one dimension, this is the generic case, in the sense that a function with extremum at which both of the first two derivatives vanish (e.g. a quartic) is transformed by small perturbations into a solution in which there are two non-degenerate extrema.
Thus this approach to limiting using the second-derivative condition will preserve high-order accuracy for any order method. The extension to multiple dimensions is done by using the derivatives in each of the coordinate directions, with additional conditions where the solution varies in only one of the coordinate directions, or in the neighborhood of points where the solution is sufficiently close to a cubic [28]. In any case, the cost of applying these more elaborate limiters at extrema are small, since they are only computed at locations where there is a possibility of a smooth extremum, a condition which is inexpensive to compute, and occurs typically on a set of codimension one. There are two possible approaches to extending the limiters described above to the case of systems of equations. The first is to apply the limiter componentwise, possibly after first performing a nonlinear change of variables, e.g. transforming from conserved variables to primitive variables in gas dynamics. This is a straightforward generalization, but can fail for more complex systems, such as magnetohydrodynamics. The second approach is to apply the limiters in characteristic variables. To do this, one expands either δF d i+ 1 2 e d in (36) or ( U i − U i+ 1 2 e d ) in (37) as a sum of right eigenvectors, compute a limiting coefficient for each amplitude, and resum the expansion with the corrected amplitudes to obtain the limited values. In the scalar case, these coefficients are obtained as nonlinear combinations of differences. For characteristic limiting, the limiter for the amplitude of the k th eigenvector is computed as the same nonlinear combination of the k th expansion coefficient of the same derivatives as used in the scalar case. While this is a sounder approach theoretically, it also must be done with care, due to the nonlinear dependence of the eigenvectors on the solution. All of the expansions in characteristic variables used to compute a given limiting coefficient must be done with respect to eigenvectors evaluated at the same value of the solution. Positivity Preservation. Solutions to the advection equation remain nonnegative for all times if they are initially nonnegative. It is desirable that this property be preserved under discretization, for example, if other terms in a system of equation are not defined unless an advected quantity on which they depend is nonnegative. Such a property of a discretization is called positivity preservation. In the initial development of limiters in one dimension, positivity preservation was considered one of the design goals. In more than one dimension, this becomes problematic.

PHILLIP COLELLA
The only classical positivity-preserving linear low-order finite-volume method for advection by an arbitrary spatially-varying, discretely-divergence-free velocity field is donor-cell differencing, which has a CFL constraint that scales like D −1 . In high dimensions, such as for kinetic problems in 6D phase space, this is a significant time step penalty. A second issue is that positivity preservation constrains the design space for extremum-preserving limiters, even in one dimension: in order to allow for smooth extrema to be preserved, it is necessary to relax constraints on the max and min of the discrete solution so that small violations of the bounds are introduced. An alternative that has seen some success has been to use a separate post-processing step after each time step to preserve positivity (or other bounds satisfied by the PDE) [33,11,20]. At the end of every time step, one computes the increments in each control volume corresponding to the extent to which positivity has been violated. Those increments are then redistributed to nearby control volumes, with weights proportional to the mass of the quantity in that control volume, and with weights summing to one. If at the beginning of the time step, the averages were all nonnegative, there is some nearby control volumes that have a surplus of mass corresponding to the negative increment being redistributed. Thus the advected quantity is conserved, and positivity is preserved. Our experience is that, in the presence of discontinuities, redistribution cannot be used by itself, but must be applied to solutions to which limiting is also applied.
3. Adaptive mesh refinement. In the block-structured finite-volume adaptive mesh refinement algorithm in [5] the solution to a system of conservation laws is represented on a hierarchy of locally-refined grids. At each level of refinement = 0, . . . , max , a subset Ω l of the spatial domain Ω is discretized into a union of control volumes of the form (3), which are in turn organized into a disjoint union of rectangles k Ω ,k = Ω (figure 1). The grids are assumed to be nested: Ω ⊂ Ω −1 , h −1 = h n ref , where n ref ∈ N + is the refinement ratio. The control volumes defined by these grids are assumed to be nested, i.e. that any control volume in Ω that intersects Ω +1 is completely covered by Ω +1 . The additional requirement of proper nesting is usually imposed, i.e. that the distance between the boundary of Ω +1 and any control volume in Ω −1 not covered by Ω is at least some integer multiple of h . This proper nesting distance is chosen to make the calculation of ∇ · F on Ω use only those values defined on Ω { , ±1} .
For a Cartesian grid hierarchy, the calculation of ∇ · F on Ω , where F =   Computing fluxes at refinement boundaries. The flux into the coarse control volume to the right of the refined grid is given by the arithmetic average 1 2 (F L,top + F L,bot ), which is also the flux out of the adjacent coarse control volume covered by the fine region. Since arithmetic averaging of fine face averages to obtain a coarse face average is exact, the accuracy of the flux on the refinement boundary is the same as the accuracy of the component fluxes.
This algorithm satisfies the natural generalization of the discrete conservation property (7) to nested grid hierarchies. This follows easily from the replacement of the fluxes in Step 2 with averages from the next finer level, combined with the replacement of values of ∇ · F on the parts of Ω covered by Ω +1 by the arithmetic averages of ∇ · F defined on the latter. A key observation is that arithmetic averaging from fine grids to coarser grids of averages over control volumes and over faces is exact in the sense of truncation error. This implies that if the truncation error of the single-level method is O(h p ), then the truncation error in the calculation of the divergence on the hierarchy is O(h p ) away from boundaries between different levels of refinement, and O(h p−1 ) near refinement boundaries. For all faces, the error in the flux calculation is O(h p ). Near refinement boundaries, the coefficient multiplying that error is not C 1 , so the cancellation of error in the final flux divergence calculation does not hold, thus leading to a loss of one order of accuracy. For classical PDE, this loss of accuracy is often less, even in max norm. For elliptic and parabolic PDE, the solution error is smoothed by the solution operator, typically leading to an O(h p ) solution error. For hyperbolic PDE, the solution error can be O(h p ) if the refinement boundary is non-characteristic. This is suggested by a modified equation analysis similar to that used in [13], and observed in practice [2,28].
Formally, we view the solution to a system of conservation laws on a nested grid hierarchy as being defined only on the control volumes at each level that are not covered by the finer levels, i.e. on Ω − Ω +1 . Furthermore, we need to interpolate values at control volumes contained in the stencils for operators applied to level control volumes that extend outside of Ω . In the former case, we compute U in control volumes covered by the next finer level to be the arithmetic average of the values on the next finer level. Since arithmetic averaging of averages over control volumes is exact, the latter are p th -order accurate, so is the average. To interpolate in space from the next coarser level, it is convenient to use a least-squares approach to define interpolation stencils [28]. This is similar to what is done for mappedmultiblock methods, so we will defer discussion to the next section. In addition, AMR algorithms for time-dependent problems often employ refinement in time, so that values of U at intermediate corresponding to the time steps at the next finer level are required. For the Runge-Kutta methods used here, this is done by means of dense output, i.e. reconstructing a polynomial in time from the intermediate stages of a p th -order Runge-Kutta method to obtain a p th -order accurate approximation to the solution at any intermediate time. This is done on for the coarse-grid values near the refinement boundary prior to performing interpolation in space.
For mapped grids, there is a potential additional complication in defining discretizations on the AMR hierarchy, in that the discretized grid mappings and the associated quadratures on faces may not match up exactly between refinement levels. In the case where the mapping is specified analytically, and in addition that the averages J i and N s d i+ 1 2 e d can be computed analytically, arithmetic averaging between levels is exact, and there is no problem. In [6], this issue was addressed by constructing the grid at the finest possible level, then constructing the coarser levels by averaging these quantities down from the finer grids, so that the requisite consistency conditions are satisfied. This is not a practical solution in 3D, particularly when the grid hierarchy is changing as a function of time. Instead [19], we use (29) and (34) to construct values for J i and N s d i+ 1 2 e d that are consistent across refinement boundaries, and, in the case of N s d i+ 1 2 e d , satisfy freestream preservation. In the case of J , this is just the calculation of ∇ ξ · (D −1 N T X) on the AMR grid hierarchy in the Cartesian mapping space. In that case, the discrete conservation property becomes conservation of volume independent of how the grid is refined. The calculation of the face averages uses the same idea, except one dimension lower. For all edges on level that are covered by edges on the boundary of level + 1, we replace the integrals in (35) over the coarse edges by the sum of the integrals of the fine edges ( figure 3). By a similar argument to the divergence case, the resulting calculation of N s d i+ 1 2 e d satisfies the freestream condition. There is no loss of accuracy of the edge integrals, since summing integrals is exact. however, there is a loss of one order of accuracy from differencing the edge integrals. Therefore, in order to obtain p th -order order accuracy for the fluxes, the edge integrals in (35) must be computed with a (p + 1) st -order quadrature rule. 4. Complex geometries. There are two geometric challenges related to generalizing the discretization methods described above that we address here. One is representing simple domains, but where it is desirable for reasons of efficiency and accuracy to align the grid with global anisotropies, while at the same time there are topological constraints that make it impossible to do so with a single smooth mesh mapping of the type described in Section 2.2. A simple example of this are problems on thin spherical shells, such as the earth's atmosphere. Large-scale motions can be shown to have large aspect ratios [18], and therefore the it is desirable to have the grid aligned with the normal and surface directions to the sphere, i.e. have the map be the tensor product of a map to the surface of the sphere and a 1D mapping in the radial direction, with the ratio of the grid spacing in those directions comparable to the aspect ratio. It is not possible to do that with a single mapping of the form (27) without introducing singularities, such as the pole singularities in classical spherical coordinates. A widely used alternative is the cubed-sphere grid, in which the surface of the sphere is represented by a mapping from six coordinate blocks of the form (27) that are continuous at boundaries, but cannot be represented in terms of a single smooth coordinate map ( figure 4).
For highly complex geometries coming from CAD specifications of engineering devices, geophysical data, or image data, we use Cartesian grid embedded boundary (EB) / cut-cell methods [33,11,7,4,34,21,1], an irregular domain Ω is discretized by intersecting a rectangular grid with the domain to obtain a collection of control volumes. One of the principal advantages of EB methods is the ease of what corresponds to grid generation, i.e. the calculation of volume fractions, the area fractions, and the centroids. These can be computed easily for highly complex geometries, such as the exteriors of full aircraft configurations defined by surface triangulations of components [1], and microscale geometries of porous media obtained from image data [38]. This uses the cubed-sphere map, which maps the surface of a cube onto the surface of a sphere. The map is smooth, except at boundaries between blocks (represented by the heavier dark lines) corresponding to the boundaries between different faces of the sphere, where the map has discontinuous derivatives in the direction normal to the boundary. The tensor-product of a cubed-sphere map with a radial coordinate can be used to represent a thin 3D spherical shell.
For mapped multiblock methods, the basic discretization approach is that in Section 2.2, and the main issue is the definition of that discretization method near the boundaries between blocks that preserves high-order accuracy. For the embedded boundary approach, we have to define methods for computing high-order approximations to the average of ∇ · F over irregular control volumes generated where Cartesian control volumes intersect the irregular domain, as well as for faces are sufficiently close to the boundary to require an irregular stencil. In both cases, the problem is one of generalizing the constructions used in Section 2 to compute fluxes to cases for which the unknowns are expressed in terms of averages of smooth functions over control volumes that fail to be smooth rectangular discretizations of space on a set of codimension one. The approach we have taken is based on constructing polynomials whose averages over all sufficiently-nearby control volumes are known. We use more control volumes than polynomial coefficients, leading to an overdetermined linear system for the coefficients or, equivalently, an underdetermined system for the stencil coefficients that express fluxes in terms of nearby averages. To solve this overdetermined system, we use least-squares methods, with the introduction of weighting in some cases to obtain stable discretizations.

4.1.
Mapped multiblock methods. The method outlined here for computing ∇ · F on a mapped-multiblock grid is that given in [30]. In that approach, the method in Section 2.2 to computing fluxes can be used without modification for control volumes sufficiently far away from block boundaries. To compute the fluxes near to the block boundaries, we assume that there is a smooth extension of each of the mappings beyond the boundary of the block. Using this extended mapping, we can define an extended rectangular grid and associated control volumes sufficient to compute the fluxes for all of the faces surrounding the control volumes inside the block. To compute the average of the solution at i for block s, we construct a polynomial approximation to the solution U near i To compute the coefficients a p , we use the conditions that the average of the polynomial over a collection of nearby valid control volumes is equal to the average of U over those control volumes.
where the moments x p (i ,s ) are computed using quadratures. Given (i, s), we need to specify how to how to choose N (i,s) . For the case of fourth-order accuracy, we first find (x 0 , s 0 ), the valid control volume containing the center of (i, s). Then we add to that set all of the valid control volumes that share a face, edge, or corner with (x 0 , s 0 ). Finally, we add the valid control volumes that share a face opposite the shared faces from the previous set. Away from three-or-more block intersections, This stencil is a set of points of size (3 D + 2D) (figure 5). At three-or-more-block intersections, the number of points can be larger or smaller, but in any case no less than the number of coefficients {a p }, and is sufficient to construct polynomials with P = 4. In two and three dimensions, (39) is in general an overdetermined system of equations. However, it has maximal rank, is well-conditioned, and a least-squares solution can be computed. This same approach is used to compute ghost values near domain boundaries, with slightly different heuristics [30] for choosing N (i,s) .
Once the coefficents have been constructed, the ghost value can be computed as This provides us with the necessary ghost values to compute fluxes on all of the faces of each block. This leads to fluxes that are double-valued on faces at block boundaries. We can compute a single-valued flux either by averaging the two fluxes (appropriate for elliptic operators), or using some upwinding method for computing a single-valued flux (appropriate for hyperbolic operators).

4.2.
Embedded boundary methods. The discretization of space given by the intersection of Cartesian control volumes with the irregular domain leads to a finitevolume discretization corresponding to (5) (figure 6), where v is a connected component of Ω V i , a ± d = ∂v A i± 1 2 e d , and a EB = ∂v ∂Ω. The quantity κ = κ(v) is the volume fraction of v, i.e. V olume(v) = κh D . Globally, the geometry is represented as a graph: the nodes of the graph are the control volumes, and the arcs of the graph are the faces α d connecting a pair of control volumes in adjacent cells. Since the a connect exactly two adjacent control volumes, the discretization (41) satisfies a discrete conservation principle analogous to (7): the volume-weighted sum of ∇ · F over a region defined as a collection of control volumes is equal to the area-weighted sum of the flux averages ± F d d,± , F ·n EB summed over the boundary of the region. The relationship (41) is exact, and the most commonly used discretization is based on a second-order accurate approximation of the flux, combined with the analogue of the midpoint rule: where the area fractions α are defined by α d,± h D−1 = Area(a d,± ), α EB h D−1 = Area(a EB ), and x d,± , x EB are the centroids of the corresponding a's. We note that there is a loss of one order in h in the truncation error relative to what we get from the midpoint rule in control volumes that do not intersect ∂Ω. This is a loss of one order of accuracy on a set of codimension one, which results in second-order accuracy in max norm for elliptic problems, and at least in L 1 norm or possibly max norm or for hyperbolic problems, depending on whether or not the boundary is characteristic [13]. In addition, the presence of the volume fractions κ in denominators, which can take on arbitrarily small values, might be a cause for concern, both from a stability and an accuracy standpoint. This turns out not to be a problem. For elliptic problems, or parabolic problems that are discretized implicitly in time, the volume fractions in the denominator are eliminated by diagonal scaling [21,29], while still leading to stable discretizations. For time-dependent hyperbolic problems, one can still use explicit time discretizations by hybridizing with a stable, non-conservative method, and redistribute the increments corresponding to the local loss of mass in each irregular control volume to nearby control volumes in a way that is stable overall [33,11,4].
We generalize the EB approach to higher order for domains Ω that are defined in terms of an implicit function: Ω = {x : φ(x) < 0}, where φ is assumed to be smooth for the purposes of constructing high-order accurate methods. This is a representation that arises naturally from image processing based on level-set methods [27], and is increasingly an approach used in CAD systems based on constructive solid geometry. In that case, ∂Ω = φ −1 (0), and We take the Taylor expansion of of the fluxes inside the flux integrals (43) -(44).
where we use the representation (46) to define the Taylor expansion ofn and the points x d,± , x EB are points on the Cartesian face containing a d,± , and in the Cartesian control volume containing v, respectively. Thus grid generation in this extension of EB to higher order reduces to computing the moments over a d,± , a EB in (48), (49), along with moments over the volume fractions Once the moments are computed, the computation of ∇ · F v is reduced to finding suitable approximations to the ∇ q F d , or equivalently, a polynomial approximation to F d of the appropriate order. Following [37], we compute suitable approximations to these moments by applying (41) to fluxes of the form The error term in this expansion is sufficient to lead the asymptotic error term in (49) unchanged. For each R ≥ 0, the equations of the form (51) for {r : ||r|| 1 = R, d = 1, . . . , D} form a well-conditioned overdetermined system of maximal rank for the moments over v and a EB on the left-hand side, where the right-hand side consists of higher-order moments over a EB or moments of the form (50) for dimension D = D − 1. The hierarchy in R terminates at R = Q, for which the contribution to the right-hand side of integral over a EB vanishes, so that, if one can compute the moments over a d,± , we obtain a closed hierarchy of least-squares problems for the moments which can then be solved, starting at the highest moment. We can apply this argument recursively in dimensionality, until we get to the case D = 1, for which the moments can be computed analytically by finding the intersection of coordinate lines with ∂Ω.
There are two cases in which the above construction must be modified. One is in the case where a D < D calculation involves a zero set that grazes a D coordinate surface; the second is where thickness of the region {x : φ(x) > 0} is comparable to h (figure 7). In both cases, the solution is the same, i.e. to locally refine the grid in the process of generating the control volumes and moments. In the latter case, this will also lead to multiple control volumes per Cartesian cell, and a graph representation of the embedded boundary grid that is not a subgraph of the lattice Z D . Figure 7. Example of an underresolved collection of control volumes in 2D. The region above the heavy curve is contained in the domain. However, a simple representation of the geometry and topology of the embedded boundary grid based on computing a single intersection of the coordinate line with the face is inadequate: there are two disconnected components to that face that connect the single volume on top to a pair of disconnected control volumes below. We deal with this problem by locally refining the grid as indicated here to represent accurately the geometry, as indicated by the dotted lines.
The main problem for the higher-order EB algorithm formulation is finding approximations to the flux and its derivatives that are stable. There is no rigorous theory for stability, so the goal is to find a methodology and a set of heuristics for designing stable methods. One example of such a combination is contained in [17] for Poisson's equation. In that case, a least-squares algorithm similar to that in used to compute ghost cells at block boundaries in Section 4.1 is used to compute a polynomial interpolant to the flux on the faces. The essential new heuristic that leads to a stable discretization is the use of weighted least-squares, in which the control volumes nearest the face are weighted the most heavily. This has the effect of more heavily weighting the diagonal of the resulting matrix, while still providing enough degrees of freedom to obtain higher-order accuracy. 5. Conclusions. The principal focus of the discussion here has been to lay out an approach for extending classical finite-volume methods for (1), (2) to obtain discretization methods of any order of accuracy in the truncation error. It is based on the use of uniform or smoothly-varying grids to define a rectangular lattice of control volumes. Because of the grid smoothness, the truncation error in the discretization of the operator is O(h p ) provided that the face averages of fluxes are also computed to O(h p ), and that the truncation error is a smooth function of its inputs and location in space. This leads to a loss of one order of accuracy in truncation error when the latter condition is not satisfied, e.g. at domain boundaries, at boundaries between different levels of refinement in AMR, at block boundaries for mapped-multiblock methods, and at control volumes near the irregular domain boundary in embedded boundary methods. For classical PDE, this loss of accuracy in the truncation error is often smoothed out by the solution operator, leading to a solution error that can recover O(h p ) accuracy, even in max norm. The approach has been generalized to discretizations that are approriate for representing a variety problems in complex domains.
Within the context of the approach given here, there are still some possible extensions that would increase the geometric flexibility of the methods, specifically combinations of mapped and mapped-multiblock methods with embedded boundary geometries. This could be useful, for example, in representing orography in atmospheric modeling, or CFD in internal combustion engines. In these cases, a mapped-multiblock background grid would be used to represent some global feature: the thin-layer anisotropy in atmospheric modeling, or in case of the combustion problem, alignment of the grid with the cylinder wall in order to efficiently represent the large-scale viscous boundary layer. Then the embedded boundary approach can be used to represent the detailed geometry of variations in the surface of the earth, or of valves, injectors, piston geometries, and other detailed engineering features in the cylinder. Finally, for moving or free boundaries, the methods described in Section 4.2 can be applied in space-time to obtain the moment information required to derive finite-volume discretizations for the space-time divergence of (U, F ) to generalize the methods in [11,4] to any order of accuracy.