Adaptive clipping-and-redistribution algorithms for bounded and conservative high-order interpolations applied to discontinuous and reactive flows

A new adaptive clipping-and-redistribution method is presented which provides bounds-preservation for multidimensional interpolation in the context of high-order finite-volume discretizations with adaptive mesh refinement (AMR). The underlying finite-volume method (FVM) for the computational fluid dynamics applications is fourth-order accurate for smooth solutions and utilizes AMR for computational efficiency in solving multiscale problems involving turbulence and combustion. High-order interpolation between different AMR levels is required. However, this operation often leads to numerical issues because combustion species must have physical bounds preserved. The present study overcomes two major challenges in the development of the high-order interpolation method. First, the method needs to be bound-preserving near extrema or discontinuities to prevent the emergence of unphysical oscillations while maintaining fourth-order accuracy in smooth flows. Second, the method needs to satisfy the conservation requirement in multiple dimensions, particularly in the context of curvilinear coordinate transformations. Additionally, the method is designed to be localized and computationally inexpensive. The new interpolation scheme is demonstrated by solving reacting flows, which are extremely sensitive to unphysical overshoots in conserved quantities. The test problems are shock-induced H 2 -O 2 combustion and a C 3 H 8 -air flame in a practical bluff-body combustor. Results show the method prevents new extrema near discontinuities while maintaining high-order accuracy in smooth regions. In particular, the method is extremely beneficial for combustion with stiff chemistry. With the proposed new method, even if flame fronts cross AMR interfaces or new grids are created in the vicinity of the flame, solution stability is retained.


NOTATION
The notation is generally explained as it is introduced. Bold type is reserved for real physical vectors, such as solution variables of a system of the governing equations. However, a bold symbol can also represent vectors in integer space. Matrices are named with capital letters, such as N, and the entries are named with subscripts, N i,j , where i and j are the indices for the row and the column, respectively. The vector symbol ⃗ x is used for vectors (or column matrices) depending on space dimensions. Some symbols are listed and defined here for convenience. D The spatial dimension

INTRODUCTION
Numerical strategies for finite-volume methods (FVMs) such as high-order schemes and adaptive mesh refinement (AMR) on structured grids have proven themselves valuable for improving computational efficiency. 1 However, these methods are not without challenges, especially when solving combustion physics. The purpose of AMR is to better resolve regions of interest in a solution in order to reduce the discretization error. AMR requires interpolation in both space and time at the boundaries of grid refinements, to enable local time stepping at each resolution. When the underlying FVM is high-order accurate for smooth flows, as in the present study, oscillations can appear near extrema or discontinuities due to the use of a high-order conservative interpolation scheme. These overshoots are particularly troublesome when they become unphysical, and ideally AMR would tag (mark locations for refinement) discontinuous features before they form. When solving compressible reactive flows, shock waves and flame fronts are features of interest. A high-order interpolant near or crossing these features often brings about oscillations, resulting in undesirable characteristics such as negative mass fractions, or such that the sum of all species' mass fractions are greater than unity in some localized region. Eventually, these unphysical phenomena can destroy the solution. Two types of interpolation are required for the finite-volume AMR strategy discussed herein. Interpolation is used to fill ghost cells surrounding a patch of fine grid cells, and interpolation is used to define new fine grids as the solution evolves. For both cases, standard practice is to adjust the fine grids so that discontinuities always remain well away from refinement boundaries, and the interpolations are of smooth solution profiles. Having a discontinuous feature encroach on a grid refinement boundary is avoided. However, for reasons of computational cost, this is not always practical. In many simulation strategies, regions are refined based on criteria other than following a discontinuity. For example, if reproducing an experiment, one may add fixed levels of refinement that match a physical observation window to focus computational effort, and minimize error, in that window. Strategies such as embedded/large/eddy/simulation (embedded-LES) and embedded/direct/numerical/simulation (embedded-DNS), to be discussed shortly, have similar requirements. For some problems, such as combustion described by highly sensitive reaction mechanisms, tiny overshoots and undershoots in bounded species mass fractions may be enough to corrupt the solution, in regions well away from the identified flame discontinuity. Although it is not the preferred way to apply AMR, there are strong motivations for robustly and accurately interpolating in the vicinity of discontinuities. The goal of this work is to create high-order accurate AMR methods that are more robust to tagging and refinement variability, especially when AMR interfaces are near regions with strong gradients.
To demonstrate the problem specifically, an example of high-order interpolation is shown in Figure 1 where significant erroneous features appear due to mesh refinement in patches highlighted by violet bubbles. The figure shows the distribution of CH 2 O species mass fraction from combustion of a C 3 H 8 -air mixture. Further details for this case are explained in Section 5.2, but here we mainly point out the oscillations produced by the high-order interpolation. These seemingly small noises, resulting from mesh refinement, have been shown to grow and interact with the flow dynamics, eventually leading to numerical instability. This is observed to be especially problematic for chemically reacting flows where small differences in mass factions are able to push the chemical reaction to unphysical states that lead to solution failure due to the stiff nature of the reactions. Therefore, for combustion with stiff chemistry, this necessitates a high-order interpolation approach that is not only bound-preserving but also maintains conservation while providing a high-order accurate solution in smooth regions. By contrast, for calorically perfect gases, even at high Mach numbers, it is the authors' experience that bounded interpolation is not necessary as long as discontinuities are contained in fine grids.
Another example is controlling oscillations that may appear from high-order interpolation when studying turbulent flow applications utilizing embedded-LES and embedded-DNS methods. 2 In these approaches, only a small subset of the domain has sufficient resolution for a detailed query of the physics of interest using LES or DNS modeling of turbulence. Elsewhere, the domain can be considered as approximating the influence of realistic geometries and boundary conditions using far more affordable methods. With additional control of error, for example, via data assimilation, it is a way forward when full LES or DNS is otherwise completely infeasible due to the computational expense. Since embedded-DNS and embedded-LES cases can only afford limited amounts of mesh refinement, discontinuous flow features often must span AMR interfaces. This is illustrated by Figure 2, where only a small window of a flame front is resolved using DNS. In this example, the discontinuous flame front must cross multiple different levels of grid refinement. To reduce computational time, there is also interest in allowing the flame to only develop on the coarse level for an amount of time before refining and using it as an initial condition for the embedded-DNS region. For these applications, it is important to have an algorithm which is able to mitigate the oscillations present in high-order interpolation to prevent unwanted numerical instabilities due to abrupt AMR-driven changes in grid resolution. In summary, we have devised a new high-order bound-preserving interpolation scheme, which improves the robustness of AMR for high-order FVMs when modeling combusting flows. The new high-order adaptive clipping-and-redistribution (HO-ACR) method has been implemented, verified, and validated with our in-house computational fluid dynamics (CFD) code Chord. 1,3-10 Chord employs fourth-order finite-volume discretizations in space and fourth-order Runge-Kutta methods (both fully explicit and additive implicit-explicit) for time integration. AMR in both space and time is key to enable computational efficiency for multiscale problems. Mapped structured grids (e.g., transforming a complex physical geometric domain to a Cartesian computational domain) are used to accommodate the geometric complexity and utilize AMR. Chord is capable of solving a wide range of compressible fluid dynamics with combustion and turbulence including Reynolds-averaged Navier-Stokes and large eddy simulation for turbulence. Chord is built upon the Chombo framework 11 and implemented using flat message passing interface (MPI) for parallel computing on high-performance computers with 10 5 or more cores.

Existing methods to cope with unphysical numerical solutions
Traditional CFD approaches for managing interpolation overshoots often utilize a measure to detect discontinuities and locally revert to low-order methods in their presence. 12 There is a wealth of well-developed methods on this topic in the form of flux limiters 13 and the family of essentially non-oscillatory (ENO) schemes. 14 Common methods for bounds preservation in CFD schemes exist in the form of limiters. The piecewise parabolic method (PPM) limiter 13 is well known for its application in high-order schemes, and relative simplicity. This method operates by constructing high-order approximations at cell faces, while ensuring these constructions do not generate new extrema and are monotonic. Furthermore, solution smoothness is accounted for in order to prevent excessive limiting and reduction in the order of accuracy in regions of smooth extrema. Although well-developed, in general limiters are only designed for face reconstructions and most often applied in one dimension at a time. Approaches for limiters in multiple dimensions do exist, 15 but are not easily modified to achieve strict bound-preserving cell-based interpolation.
In a similar spirit to limiters, the flux corrected transport (FCT) methods utilize low-order solutions to provide bounds for high-ordered flux reconstructions. Unlike limiters, FCT methods have been extended to operate in multiple dimensions. The FCT-like maximum preserving principle (MMP) in the work by Anderson et al. 16 enforces the flux correction by instead adjusting Discontinuous Galerkin elements. This adjustment requires satisfying flux constructions, while redistributing and preserving element mass. Despite the approach having many valuable properties, the aims are subtly different from ours and can not be directly translated into a FVM scheme.
One approach that is multidimensional, conservative, and solves the spatial interpolation problem we face is the central essentially non-oscillatory (CENO) method described by Ivan et al. 17 This approach formulates a smoothness measure to identify regions which are nonsmooth. In nonsmooth regions a low-order interpolation is used, while a high-order interpolation is employed elsewhere. Although it is an effective approach, it requires a tuned coefficient for smoothness detection. Additionally, large stencils of data can be required which necessitate expensive data communication, but further developments of this approach are able to reduce stencil sizes by more frequently moving smaller amounts of information. 18 Implementing this interpolation approach other than within the CENO scheme would be nontrivial since it requires careful consideration of limiter choice and significant changes to communication patterns.
Another approach for the multidimensional, conservative interpolation problem qis the clipping-and-redistribution algorithm described by Hilditch et al. 19 After a high-order interpolation, any solution outside physically imposed bounds is clipped off. The sum of clipped quantities is tracked and redistributed elsewhere inside the domain to maintain conservation. This approach is developed to prevent physically constrained overshoots based upon predetermined bounds. Additionally, the redistribution happens globally and thus requires solving an implicit system which is potentially expensive. Nevertheless, this work is the inspiration for the algorithm that we present. We seek to adaptively enforce bounds near discontinuities and reduce the cost of redistribution. In addition, the present algorithm operates on mapped quantities since the underlying CFD infrastructure uses mapped structured grids to accommodate complex geometries.
None of the above approaches fully match our required criteria for refining complex combusting flows. Although the interpolation method in the CENO scheme is the nearest to our goals, it is not straightforward to incorporate into our existing code. In the present study, we develop the high-order adaptive clipping-and-redistribution (HO-ACR) method and successfully demonstrate the design properties of the interpolation: conservative, solution-dependent bounds-preservation, efficiency, interpolation without compromising accuracy in smooth regions. Additionally, the HO-ACR algorithm is designed to operate as a correction to a high-order interpolation. This makes it noninvasive and relatively simple to implement independent of the scheme it is used in.
The rest of paper is organized as follows. Section 2 describes the foundational numerical methodology including that of the FVM, AMR, and unbounded interpolations used to fill ghost cells and populate newly refined grid levels. Section 3 presents the details of the newly developed HO-ACR algorithm, including the solution procedures and simple verification tests. In Section 4, the HO-ACR algorithm is implemented in Chord and verified by advection physics. Results and discussions are provided in Section 5, which both validate the methodology and demonstrate the effectiveness for complex flows. Finally, conclusions are drawn in Section 6 and future follow-up work is discussed.

BASE METHODOLOGY
Although not the focus of this paper, the high-order FVM on mapped structured grids is briefly described for convenience since the new HO-ACR method is developed, implemented, and applied in this test bed. A more complete description is available for Cartesian grids in McCorquodale and Colella 20 and for mapped grids in Guzik et al. 1,21

High-order FVMs
Consider hyperbolic PDE's of the form FVMs approximate the integral form of the conservation laws with many contiguous control volumes The underlying FVM code employs mapped grids to accommodate geometric complexities and the Cartesian grid in the computational domain also facilitates the use of AMR. The Cartesian grid, with coordinates ⃗ , is transformed into physical space, with coordinates ⃗ x by the forward mapping ⃗ x = X( ⃗ ), as illustrated in Figure 3. Using the mapping function, geometries of moderate complexity can be solved while maintaining the computational advantages of Cartesian grids including tractability of AMR.

Node
Cell grid grid

F I G U R E 3 A representation of the transformation between computational and physical space
For mapped grids, the integral equations are transformed into computational space, ⃗ , by using grid metric terms such that where the transformation matrix, N ⊺ = J ⃗ ∇ x ⃗ , describes the grid metrics and the metric Jacobian J ≡ det( ⃗ ∇ ⃗ x). Cell-averaged quantities are defined as and together with the divergence theorem yield a discretized form where the subscript d denotes the dth row of N ⊺ . The term ⟨N ⊺ d ⃗ F⟩ is devised such that free-stream preservation is maintained. 1 This formulation is exact until numerical methods are used to approximate fluxes and time derivatives.

Adaptive mesh refinement with mapped grids
Adaptive mesh refinement is effective in solving multiscale problems. When discontinuities appear in a solution, high-order methods must locally revert to low-order methods around discontinuities. The decrease in order of the scheme is necessary and adversely impacts the solution accuracy. To remedy this, locally refining the grid near the discontinuities with AMR is the best option available. Traditionally, discontinuities are well-embedded inside the fine grid. Layers of buffer cells ensure that interpolations to fill ghost cells and the fine patches are kept away from any discontinuities. However, there are situations in combusting flows and embedded-DNS where it is not practical to fully refine every discontinuity. Building upon this base methodology for AMR, we develop the HO-ACR method in Section 3 to make AMR more robust for nonideal circumstances. While the details of AMR in Chord can be found in McCorquodale and Colella 20 and Guzik et al., 1   In addition to the refinement in space, AMR can make use of subcycling for refinement in time, as detailed by McCorquodale and Colella. 20 Briefly, subcycling updates each grid level using a different time step size to reduce computational costs. Solutions on level + 1 use a time step sized as Δt +1 = Δt ∕n ref . Levels + 1 and are synchronized at the beginning of each time step on , and operate independently between these time synchronization points. During this synchronization, the solution values on level are updated to match the cell-averages from + 1, and fluxes at coarse-fine interfaces are corrected to ensure conservation. Flux correction replaces the fluxes on the coarse grid with those from the fine grid at the interface between the two levels. At intermediate steps on the fine level + 1, ghost cells are interpolated from the coarse level in space and time.

Conservative interpolation schemes for AMR
With an AMR grid, each level is intended to advance independently with data transfer between levels occurring at only refinement boundaries. For the FVM, the coarse cell-average is defined exactly by the average of fine cell-averages. From coarse to fine, data is interpolated. There are two AMR operations where interpolation is used directly. When regridding, a new level + 1 is created, and its valid region is interpolated from existing coarse data on level . This interpolation must be conservative to maintain the strict conservation property of finite-volume schemes. At AMR interfaces, a number of operations to communicate between levels and + 1 rely on having invalid ghost cells of level + 1 interpolated from level . An example grid with one layer of invalid ghost cells is show in Figure 4. This interpolation needs not be conservative because invalid ghost cells are only intermediate states used for flux reconstruction, although enforcing conservation may improve the quality of the solution. 1 These two interpolations between levels have slightly different implementations, but both follow this procedure: 1. For each coarse cell i , determine an interpolation stencil (i ). The stencil must have sufficient degrees of freedom for the desired order of accuracy.
2. Construct a scheme-order multidimensional polynomial using the data in (i ) while maintaining conservation of ⟨U⟩ i . 3. For each desired fine cell i +1 ∈  −1 (i ), evaluate the cell-averaged value using the prior reconstruction.
The interpolation procedure for both newly refined cells and invalid ghost cells as described by McCorquodale and Colella 20 is a form of k-exact reconstruction. 22 This process is prone to problematic oscillations when done at high-order, independent of whether a conservation constraint is enforced. Due to this similarity, bound-preservation can be enforced using the same method for both newly refined cells and invalid ghost cells.

Multidimensional polynomial interpolation
When interpolating an AMR solution, the idea is to construct a piecewise interpolant, U i ( ⃗ ), for a coarse cell, i. To achieve a specified order of accuracy, a multidimensional Taylor expansion about the center of cell i is used. Multidimensional notation is defined using q as an integer vector of size D. Vector powers of q of vector quantities ⃗ are expressed as A multidimensional form of the Taylor series up to order Q, based on cell centers, ⃗ i , is expressed as where U , and h is the cell spacing.
Choosing c q = 1 q! U (q) ( ⃗ i ) as the unknowns to be solved, the Taylor polynomial in Equation (7) working with cell-averaged quantities, defined in Equation (4), can be written compactly as For a unique construction, K = terms (binomial coefficients) are needed. In general, there is no obvious unique choice for how to select a stencil to satisfy exactly this many constraints.

Conservative least squares interpolation on mapped grids
Since there is no obvious choice of the stencil, it is preferred to use a stencil (i) of size N > K so long as i ∈ (i). This allows a stencil choice that appears more regular and has flexibility near boundaries. 20 In practice we use interpolation stencils as shown by McCorquodale et al. 20 For interior solutions, this stencil uses a radius of one cell from the interpolation center with an additional cell in each face-normal direction, resulting in 13 cells in two-dimensional (2D) and 33 in three-dimensional (3D). A similar stencil is used near boundaries, but shifted toward the interior to maintain the total number of cells (see Reference 20 for specific stencil construction details).
Since the stencil size is overdetermined, the least squares method must be used to create the interpolating polynomial of best fit. This still retains the desired interpolation order of accuracy and has the added benefit of allowing extra constraints. 20 For interpolation between AMR levels, preserving scalar conservation is required. To fulfill this, a constrained least squares problem of the form is formulated, where the unknowns x are a vector of the coefficients c q . The conservation constraint Bx = d is enforced using the equation while the objectives Ax − b to minimize are described by Once the interpolating function has been constructed by solving coefficients c q , fine cells are evaluated as Due to the conservation constraint, it can be shown that using the refinement operator  −1 . The values of ⟨Ĵq⟩ are computed using a fourth-order product formula and the interpolation step is shown to be freestream-preserving. 1 When filling invalid ghost cells (see Figure 4), interpolation of conservative values ⟨U⟩ is sufficient, which is solved using constant values of J in space. For regridding operations where valid cells are filled, conservation of mapped quantities ⟨JU⟩ is required. This method of least squares interpolation is said to be k-exact, 22 that is polynomial solutions up to order Q are interpolated exactly using the approach. Any least squares remainder is of the desired order for smooth solutions. Using this approach, fourth-order accurate methods with AMR have been demonstrated. 1,20

ACR ALGORITHM
The ACR method is designed as a postprocessing operation independent of the high-order interpolation. This simplifies implementation since it requires no fundamental changes to the existing high-order interpolation algorithm, and can be adaptively enabled. The method maintains conservation, avoids significant data exchange, and does not require predetermined bounds. The procedure for ACR is outlined by the following steps: 1. Interpolate the fine data from the coarse data using a high-order conservative scheme as described in Section 2.3. 2. Establish minimum and maximum bounds for interpolation values of each fine cell. These bounds are created using local coarse cell information. 3. Clip each fine cell to fall within the established bounds and track the solution quantity change. 4. Redistribute the clipped excess to local regions that can accept it while remaining inside their bounds.
A few main features deserve emphasis and clarification. The primary objective of the ACR method is to prevent new extrema while maintaining solution conservation. Notably, this does not strictly enforce monotonicity of the solution. It is important to properly establish the local clipping bounds since their selection directly impacts the resulting solution. Furthermore, the extension of the ACR to high order is desirable for maximum accuracy. At the end of the section, we provide verification of the new method.

Establishing local clipping bounds
Local clipping bounds provide a feasible range that the interpolation of each fine cell must lie within to prevent generation of new extrema. For each fine cell i on AMR level + 1, an interval i +1 containing the valid interpolated solution is determined. Bounds are created from the minimum and maximum values of a neighborhood ((i +1 )) of coarse data on level about the fine cell i +1 . In this work, the neighborhood ((i +1 )) is chosen to be the same as the neighborhood from the interpolation stencils. However, it need not be identical. Clipping bounds are determined as The ACR method will ensure the interpolated fine level respects the bounds from the existing coarse level. As a result, the ACR method does not require explicit bounds to maintain physically realizable solutions, so long as the coarse solution itself is physically realizable. The only requirement on local neighborhood choice ((i +1 )) is that it must contain the coarse cell (i +1 ); however other choices do impact the behavior of the ACR method. For example, the smallest clipping neighborhood may be set to include only the coarse cell (i +1 ). This removes any flexibility in the interpolation, and reverts to a piecewise constant interpolation. In the other extreme, the clipping neighborhood is expanded to the entire domain, and the ACR method will only restrict the creation of new global extrema. Our compromise is that clipping stencils should have the same locality as interpolation stencils; smaller neighborhoods are overly restrictive, while larger neighborhoods are not restrictive enough and are computationally impractical.

High-order extension
In regions where data is smooth and in state variables with no strict physical bounds (such as velocity), high-order interpolation is expected to perform well without assistance from the ACR method. Reducing the interpolation order in these smooth regions can compromise physics that are correctly captured by the high-order scheme. 13 In such cases, interpolation of smooth data is expected to accurately generate subtle new extrema. This is illustrated in Figure 5 where high-order interpolation correctly captures a Gaussian profile, and enforcing bounds via the ACR method results in a poor quality result. For this purpose the ACR method is extended to high-order, hereafter called the HO-ACR method. The HO-ACR method is adaptive in that the operation is not engaged in smooth regions, while being engaged elsewhere. Regions are determined as smooth if all second derivative constructions share the same sign. The second derivatives are computed for each coarse cell i , and compared to each of the directly adjacent neighbors. If the entire set share the same sign, the region is determined as smooth and the HO-ACR for fine cells i +1 ∈  −1 (i ) is disabled. It is at the discretion of the researcher whether high-order detection is used for any of the state variables. HO-ACR is used for all results presented herein except for the bluff-body combustor in Section 5.2.
Computing the curvature using a first-order accurate method is sufficient. Specifically, curvature is evaluated using a second-order centered scheme as given by where h is the cell spacing. For locations near the physical boundaries, first-order one-sided formulations such as are used as required.

Local linear redistribution
Once clipping is completed, properly redistributing the clipped quantities is necessary to maintain conservation. One of our design goals for the redistribution scheme is to be computationally inexpensive by avoiding data exchange. A redistribution region is chosen such that the redistribution can be solved explicitly without additional data exchange by choosing non-overlapping redistribution regions . The natural choice for such regions is the coarse cells where data is being interpolated from, such that  = (i +1 ).
Many CFD algorithms use a mapped grid technique to accommodate moderately complex geometries, meaning that a general curvilinear coordinate transformation is used to map the physical domain to the computational domain. Operations on the mapped quantities ⟨JU⟩ require additional care to enforce bounds on the physical solution ⟨U⟩. Bounds can only be enforced on conservative quantities ⟨U⟩ because gradients in ⟨JU⟩ may appear solely as the result of grid mapping. However, regridding requires conservation of the mapped variables as in Equation (11). The mapped HO-ACR algorithm must first solve for a second-order conservative variable On Cartesian grids, the same HO-ACR algorithm can be applied directly to conservative quantities ⟨U⟩ by setting J to one, and Equation (18) then becomes exact.
The method for localized linear redistribution begins by determining the overshoot from each fine cell ⟨U⟩ i with where max i and min i are respective maximum and minimum of the previously determined clipping bounds i (e.g., an interval). The mapped quantities that are clipped off are computed by which again is a second-order operation. Using the clipping bounds, the extra mass to be redistributed in region  is given by To determine where mass may be pushed to, each fine cell has a maximum allowable difference calculated as Using this, each fine cell can be updated using the clipping-and-redistribution as In effect, this drops the order of accuracy of the interpolation to first-order in regions which are fully clipped. This mapping process only needs to be second-order, since the HO-ACR method effectively reduces the interpolation order in nonsmooth regions where clipping is engaged. Note that the smoothness check for the high-order extension must operate on ⟨U⟩.

Conservation proof
The ACR method can be shown to maintain conservation, by demonstrating Assuming the prior equation holds, and substituting it into Equation (23), conservation can be shown by which is true by the definition of Equation (21), and therefore In general ∑ i∈⟨ U⟩ i ≠ ∑ i∈ ⟨U⟩ i except for the case of Cartesian grids in physical space when J is one.

1D verification
The property of bound-preservation of the new method is verified using a range of one-dimensional (1D) profiles from smooth to nonsmooth, including six tests-(a) a step function, (b) a smooth step function, (c) a piecewise linear function, (d) a Gaussian profile, (e) a piecewise cubic function, and (f) a discrete delta function. Figure 6 shows the spatial profiles of solution variable ⟨U⟩ and compares the performance of the HO-ACR method for the six cases on a Cartesian grid (by setting J = 1 uniformly). The conservative interpolation is fourth-order accurate using a centered 5-point centered stencil. For all cases, the HO-ACR method successfully clips off the overshoots (e.g., cases a, b, c, e, f), while preserving the smooth extrema (e.g., case d). Additionally, all 1D cases maintain conservation to machine precision.
Despite the improvement in interpolation near sharp features, there is still room for improvement in the solution. In cases (c) and (e), the high-order interpolation at one coarse-cell away from the discontinuities does not correctly capture the coarse polynomial solution. This is due to the interpolation stencil spanning the discontinuities, and as a result oscillations are present in the high-order conservative interpolation. While the HO-ACR method is not able to correct the high-frequency oscillations, it still does no worse than the existing high-order interpolation. It is important to note, the objective of the HO-ACR method is to prevent new extrema near discontinuous features and the method does not enforce monotonic interpolations of non-smooth coarse data. The emphasis of the HO-ACR method is to maintain solutions that are physically bounded, but leave the quality of these solutions up to the base interpolation scheme. Finally, we demonstrate the accuracy for some cases. For the following tests, we start with exact solution cell averages on progressively finer grids, and apply the base interpolation and HO-ACR algorithms using a factor of four refinement. The resulting convergence rates for L 1 and L ∞ norms over the O(N) domain are plotted in Figure 7.
The two smooth cases, the Gaussian from Figure 6D and smooth step from Figure 6B, behave as expected in Figure 7A and indicate fourth-order accuracy as anticipated. Note that the results are nearly identical for the base and HO-ACR scheme in smooth regions with sufficient resolution. The smooth step on coarser grids is clearly under resolved, resulting in reduced convergence rates. However, the HO-ACR reduces the error slightly at these coarse resolutions. Figure 7B shows the errors for the piecewise linear and piecewise cubic cases, corresponding to the profiles in Figure 6C,E, respectively. In both of these cases interpolation errors next to the discontinuities dominate, while interpolation is exact more than two coarse cells away (note the outermost cells in Figure 6C,E). Because of the discontinuities, the base method has an O(1) overshoot, so that the L ∞ norm is constant and does not converge with refinement. Meanwhile, the HO-ACR method shows local first-order convergence. In both cases the error of the L 1 norms are one-order smaller than the L ∞ norms because a constant number of cells interpolate the discontinuity, and the L 1 norms scale these by the cell size. Thus, while the base interpolation is first-order, the HO-ACR method increases this to second order. Using the same six cases as shown in Figure 6, we demonstrate the mapped HO-ACR performance on mapped grids. The setting of Figure 8 and Figure 9 is the same as that of Figure 6-the sequence of cases (A-F), the line color and symbols, and the interpolation scheme. Cases in Figure 8 use a mapping function of x = a tan(b ) while in Figure 9 the mapping function is x = 1 2 2 . The mapped HO-ACR performs exceptionally well for both cases. The clipped profiles are bound-preserving and smooth extrema-preserving, while maintaining conservation to machine precision. As in the prior example, the cases (c) and (e) do not correctly match the coarse polynomial solution due to the high-order interpolation. Although there is room for improvement in regions away from the discontinuity, the solution is still improved as a whole over the unbounded high-order interpolation. Notably these small regions of error do not violate any existing coarse bounds of the solution.
Algorithm 1 presents the mapped HO-ACR for ⟨JU⟩ of AMR level + 1. The same algorithm can be applied directly to conservative quantities ⟨U⟩ on a Cartesian grid by setting J to be one.

VERIFICATION
The new HO-ACR algorithm is verified by a set of test cases involving the advection of smooth and sharp profiles on mapped grids. The features to be verified are: (1) ensuring high-order accuracy for smooth flow problems; (2) bound-preservation for nonsmooth ones; and (3) maintaining conservation. Specifically the algorithm is examined for interpolation operations involving both the invalid ghost cells and the regridding process. The HO-ACR method is applied for invalid ghost cells and newly refined cells without fundamentally changing the base interpolation algorithm. The verification is first performed with the advection of a sharp profile in 1D, then advection of both a smooth profile and nonsmooth profile in 2D. A more complicated case is explored with the solid body rotation of smooth and nonsmooth profiles in 2D. All three features are demonstrated and assessed in each case.

Advection of a 1D sharp profile
To explore the effects of the HO-ACR method in Chord, a 1D scalar advection test using an Euler equation solver is explored. The initial condition top-hat profile is defined by end for end if end for on a 1D periodic Cartesian grid of length L = 2. This initial profile aligns the discontinuity with cell boundaries on the coarse grid with spacing Δx = 1∕16, which preserves the sharp discontinuity. The profile is advanced to t = 0.5 with speed u = 1. A time step size of Δt =5e-3 is used on the coarse level, and finer levels are subcycled in time.
To investigate the HO-ACR method, three versions of "extreme" AMR are compared using Chord with and without the HO-ACR method. Each of these chosen "extreme" cases violates standard practice by intentionally refining in the vicinity of discontinuous features, and anticipates the emergence of new oscillations in the solution without the HO-ACR method. In each case, a single level of AMR with a refinement factor of 4 is used. Results are shown in Figures 10,11, and 12. As references, solutions run entirely on the fine and coarse grids are shown. Due to the nature of the fourth-order finite volume method, small overshoots and undershoots appear around discontinuous features as they evolve in time, even on the coarse grid. Ideally the results in these test cases with "extreme" AMR should never be worse than the coarse solution, and good results should match the fine solution.
The first case (10) examines the effects of delaying regridding, by initializing on the coarse grid and first advancing for one coarse time step. Beginning at the second time step, the solution is either refined using the base interpolation algorithm of Chord or using the HO-ACR method. The solutions after the second time step are shown in Figure 10A. A significant overshoot is observed in the base algorithm; the HO-ACR method effectively reduces these overshoots. We quantify the improvement in Table 1 by measuring the fine solution's overshoot relative to the analytic solution. The base algorithm has overshoots that are an order of magnitude larger than the reference fine solution, while the HO-ACR modification has slightly smaller errors than the fine solution. Although this represents a significant improvement over the base algorithm, the HO-ACR method does not fully remove the small oscillations near the discontinuity, which is reflected in the larger L 2 norm than the fine solution.
To see the long-term impact, the solutions are advanced to 100 time steps, with results in Figure 10B. The overshoots from refining using the base algorithm are still present, while the solution using the HO-ACR solution closely matches the reference fine case, with extrema quantified in Table 2. The base algorithm has extrema that are nearly an order of magnitude larger than the fine solution, while the HO-ACR has extrema that almost match the fine solution. Comparing the difference in extrema demonstrates the improvement of HO-ACR, which is nearly two orders of magnitude smaller. To summarize, improper refinement of the initial condition by tagging one time step too late has significant and lasting consequences in the evolution of a solution. However, with the use of the HO-ACR method we are able to mitigate the effects of imperfect initial tagging. In more complex cases where discontinuous features appear spontaneously, this indicates that the HO-ACR method may be more accurate and robust with new flow features.
The second case examines the effects of ghost cell interpolation by specifying an initial profile with discontinuities on the coarse grid, shown in Figure 11. Directly ahead of the profile is a region of AMR fixed in space for x > 0.75. Figure 11A shows the solution as the leading edge enters the refinement region. The base algorithm produces a leading negative undershoot on the fine region, as well as an overshoot in the coarse region. The HO-ACR method prevents both of these new extrema, with minor oscillations from the fourth-order evolution scheme. In Figure 11B, the solution has been advanced until the profile has mostly entered the fine region. The extrema persist in the base algorithm while the HO-ACR method keeps them suppressed. Extrema at the end time are quantified in Table 3 relative to the analytic and fine solutions. The overshoots of the HO-ACR method are only marginally larger than the fine solution, while the base algorithm roughly doubles the overshoot. Surprisingly, the minimum values from the HO-ACR are even smaller than those in the fine solution. Because of the smoother solution from the HO-ACR method, the L 2 norm of its extrema is smaller than the fine solution's, and the L 1 norm nearly matches, while the base algorithm case is significantly worse. The differences in extrema between the base algorithm and the HO-ACR become more pronounced when compared relative to the fine solution. Relative to the extrema allowed by the fine solution, the HO-ACR method reduces the overshoot in L 2 and L 1 norms by over two orders of magnitude from the base algorithm. To summarize, this test demonstrates better robustness around sharp features and AMR interfaces when using HO-ACR. In realistic cases, this shows that discontinuous features can interact with AMR interfaces without causing excessive overshoots and undershoots. The third case, shown in Figure 12, defines the initial profile with immediate refinement based on the solution gradients. The mesh is regridded at the beginning of each time step over the entire run. However, the refinement is applied minimally so that only one coarse cell on each side of the discontinuity is refined. This causes insufficient spacing between the discontinuity and the interpolation stencil for filling invalid ghost cells. After the time t = 0.125, Figure 12A shows the base algorithm solution has new extrema emerge in both the coarse and fine regions, and begins to move away from the reference fine solution. The HO-ACR method prevents these new extrema and is able to closely match the reference fine solution with less than half the number of cells. This becomes more pronounced as the solutions advance to the final time t = 0.5 ( Figure 12B). At the end time Table 4 again shows the extrema of the base algorithm exceeding the analytic and fine solutions. Comparing extrema of the fine and HO-ACR solutions shows good agreement. Although we are hesitant to recommend this as a general practice, this test shows that using the HO-ACR method enables smaller regions of refinement (and reduced computational cost) with minimal loss of solution quality.

Advection of a smooth profile
In order to ensure the HO-ACR correctly allows for smooth interpolation on 2D mapped grids, a A Gaussian profile in density is specified by where  1 and (1.0, 0.5), respectively. An initial profile is specified on a coarse grid and advected in time. The grid mapping is defined as To test invalid ghost cell interpolation, a fixed region of AMR (also called fixed AMR throughout the paper) is added as shown in Figure 13A. The initial profile does not initially fall inside the refined region, but is advected into it over time.  Testing regridding is done by initially starting with the coarse resolution only, then refining the solution profile based on the gradient of with a level of AMR after a specified time of 0.25 time units, as shown in Figure 13B. For a smooth profile, the HO-ACR scheme is designed to allow for new interpolated extrema to maintain the fourth-order accuracy of the underlying finite-volume scheme. The profile is advected for a total of 0.5 time units, after which solution errors are calculated by comparing the numerical and exact solutions. Tables 5 and 6 show the measured solution errors in L 1 , L 2 , L ∞ for grids sizes 64 2 , 128 2 , 256 2 , and 512 2 , and the convergence rates in L 1 , L 2 , L ∞ between each two consecutive grids, respectively. Clearly, the measures verify the high-order accuracy for both the fixed AMR and regridded AMR cases using a smooth profile.
Conservation is measured by computing a relative difference in the volume weighted sum of the solution. Ideally, these differences should be within machine precision of zero between the beginning and the end of a run. Using the fixed AMR case, Chord with the HO-ACR method is measured to have a relative difference in conservative values of 2.104e-14, whereas base Chord has value of 2.326e-14 on the 256 2 case. This small difference shows that the HO-ACR method is conservative. For the regridded AMR case, the relative conservation difference for the 256 2 mesh is measured as 3.416e-13 and 2.777e-13 using Chord with and without the HO-ACR method, respectively.

Advection of a sharp profile
In this test, using a similar setup to the Gaussian profile, a sharp profile is advected in a periodic domain allowing for verification that the HO-ACR method is bound-preserving. The domain is of size [0, 1] D where D is 2. The sharp profile is the "balls and jack's" inspired by Anderson et al. 16 This profile is composed of three features prescribed in density 0 + Δ where 0 = 1.4 and Δ = 0.14, as shown in Figure 14. Relative to a profile center of (0.25, 0.375), a smaller circular shell is located at distance (0.1, −0.1) from the center with an outer radius of 0.07 and an inner radius of 0.04. A larger circular shell is located at distance (0.1, 0.1) from the profile center with an outer radius of 0.10 and an inner radius of 0.07. Finally, the jack base is placed at the distance (−0.11, −0.11) from the profile center, with one arm specified by Beginning with the initial conditions as shown in Figure 14 and a fixed region of AMR, the profile is advected forward in time for 0.492843 time units. Figure 15 shows the fixed region of mesh refinement and compares the solution profiles obtained by Chord on the left and the inclusion of HO-ACR algorithm on the right. To emphasize, the HO-ACR method establishes the bounds from the solution on the coarse grid without AMR, serving as the basis to determine overshoots. Despite only bounding invalid ghost cells, a significant reduction in regions that experience overshoots beyond the coarse regions is seen when using the HO-ACR method. Although the improvement may not be immediately present based on the number of points producing new extrema, the magnitude of these overshoots is greatly reduced as shows by the tabulated data in Figure 15. Some new extrema are to be expected, as the high-order finite-volume algorithm includes other unbounded operations.
To test the HO-ACR method for the regridding process, the initial sharp profile in Figure 14 is advected forward in time for 0.492843 time units. A new region of refinement which tracks the profile is then added as shown in Figure 16, where the refinement spans the periodic boundary. The figure compares the solution profiles obtained by Chord without (a) and with (b) the HO-ACR algorithm. Clearly, the HO-ACR algorithm prevents many of the new extrema when regridding sharp solutions as in Figure 16. Notice from the HO-ACR solution profile that a small number of values appear to violate the bounds. However, this is due to the fact that some fourth-order operations have happened between the final application of the HO-ACR and the output of the solution. These operations will naturally produce new extrema beyond prior solution bounds. One example is the conversion from one set of cell-averaged quantities to another, 23 such as done for computing primitives from conservative values or conservative values from the mapped values. Table 7 shows the solution growth of extrema, their L 1 , and L 2 -norms. For the fixed AMR test, the HO-ACR method is shown to have reduced the magnitude of the new maximum by a factor of 16, minimum by a factor 1.5, and only modest reduction in the L 1 and L 2 -norms. More impressively, the adaptive test reduces the magnitude of the new maximum by a factor of 8, and the minimum by a factor of 25. Decreases of over an order of magnitude are seen by the L 1 and L 2 -norms. As in the prior section, conservation is measured and the HO-ACR method is found to have negligible impact on the total solution conservation. The relative conservation difference of the 256 2 fixed AMR case is measured as 3.318e-15 and 5.214e-15 using Chord with and without the HO-ACR method, respectively. Similarly, the regridded case measures relative conservation of 6.099e-13 and 1.054e-13 using Chord with and without the HO-ACR method.

Solid body rotation
Another standard test problem for the HO-ACR method uses 2D solid body rotation to advect a passive scalar quantity with different localized initial profiles. The velocity field is defined by the profile u = (y, −x), where the rotational speed is = 20 in the clockwise direction. The initial profiles (Anderson et al. 24 ) are specified by = 0 + Δ ( cone (r cone ) + bump (r bump ) + cyl (r cyl )), where 0 = 0, Δ = 1, and r cone , r bump , r cyl specify the distances from the initial profile center as (cos( ∕6), sin( ∕6)), (cos(3 ∕2), sin(3 ∕2)), and (cos(5 ∕6), sin(5 ∕6)), respectively. The cone, bump, and notched cylinder profiles are given as and shown as evaluated on the initial mesh in Figure 17A. The structured grid is defined by an annulus centered at the origin with outer radius This test is motivated by tagging regions other than where the solution discontinuities are, and forcing the discontinuities to pass through the refined regions. For this imperfect application of AMR, the primary objective for the scalar solution is that it is stable and produces no new extrema as it is regridded and advects across AMR boundaries. This case is examined using three patches of counter rotating AMR at a refinement factor of 2. These patches are initially centered about ∕2, 7 ∕6, and 11 ∕6 with widths of ∕8, respectively marked a, b, and c in Figure 17A. The AMR patches are rotated in the opposite direction of the scalar advection by an angular velocity of −5 ∕6, and regridded every 10 time steps. Initial profiles are defined only in coarse regions with this choice of AMR, whereas the final profiles primarily fall inside refined regions, indicated with black boxes in Figure 17. Results from using the base algorithm of Chord and with the addition of the HO-ACR method are examined and shown in Figure 17C and D, respectively. A reference using a fine grid for the entire domain is shown in Figure 17B. In each of these cases values outside the analytic solution bounds of [0, 1] are highlighted, with overshoots in the yellow-orange color, and undershoots in magenta.
The reference fine solution, in Figure 17B, shows some small new extrema near the discontinuous features of the notched cylinder due to the fourth-order flux calculations in Chord. Even with imperfect AMR, the cone and bump profiles do not show any overshoot. This is expected since the scheme's numerical dissipation smooths continuous profiles over time, and AMR interpolation is designed for smooth solutions. With this in mind, the continuous cone and bump regions of the solutions show no distinguishable differences with or without the HO-ACR method in Figure 17C, which indicate it is not changing the algorithm substantially in smooth regions. However, the more extreme case of the notched cylinder using the counter rotating AMR with Chord in Figure 17C, shows significantly larger overshoots and undershoots. These errors are the most pronounced at the edges of the notched cylinder, especially near the inner and outer boundaries, parallel to the velocity. These larger extrema are quantified in Table 8, labeled as the "base" columns, where the peak errors are seen to be an order of magnitude larger than what is observed in the fine solution. In contrast, the HO-ACR method shown in Figure 17D greatly reduces the extrema in the solution and closely resembles the fine solution. For better visual comparison between the two cases, Figure 18 shows solutions in a window containing only the notched cylinder. The extrema are measured in Table 8, (columns labeled as "HO-ACR") where we see that the peak value errors are an order of magnitude smaller than the base case, and comparable to the fine case. The same trend is apparent in the

RESULTS AND DISCUSSION
The verified HO-ACR algorithm is applied to solve two types of cases. In the first case, a H 2 bubble is advected through a strong shock in 2D. This presents numerical challenges, including properly tracking of the thin H 2 -O 2 flame fronts and the strong discontinuities of shock waves. The second case is the C 3 H 8 -air combustion in a bluff body combustor in which the chemistry is extremely stiff. In both these cases, the HO-ACR method is applied to invalid ghost cells and newly refined cells with no fundamental changes to the base algorithm. The investigation has demonstrated that the HO-ACR method has played a critical role in ensuring the stability and accuracy in the process of obtaining solutions. Immediately what follows is the description of each problem. Then, results are presented and discussed.

H 2 -O 2 shock-induced combustion
The configuration is shown in Figure 19 and follows the setup used in a previous study. 7 A Mach 2 steady planar shock is located 0.75 cm to the right of the origin and parallel to the y-axis. A hydrogen bubble is located just upstream of the shock. The upper and lower boundaries of the y-direction are periodic, while x-min boundary condition is Dirichlet and the x-max is extrapolated. The initial velocities U u = 1.24 × 10 5 cm∕s and U s = 4.34 cm∕s are the velocities upstream and downstream of the shock. The hydrogen mass fraction is initialized to with r c the radius of the bubble and the coordinate (x 0 , y 0 ) the center of the bubble. The coefficient C 2 determines the sharpness of the interface between the H 2 bubble and the surrounding air. The case parameters are C 2 = 3 × 10 −3 cm −1 , r c = 0.28 cm, (x 0 , y 0 ) = (0.4, 0.75) cm, and = 1.5 cm. The mass fractions for the surrounding air are simplified to c N 2 = 0.767 and c O 2 = 0.233. The H 2 -O 2 combustion is a 9-species, 19-reaction mechanism described by Billet et al. 25 A traditional approach for this case would use AMR to completely refine around the discontinuities. Instead, an imperfect AMR tagging scenario that is representative of what might occur in an embedded-DNS simulation is considered. Assume only a small window of refinement can be afforded in the solution process, and due to this limitation significant features will cross AMR interfaces. Additionally, refinement is only added once the flow is developed. Ideally, the HO-ACR algorithm will handle this imperfect tagging; that is improve accuracy over a solution that uses no AMR whatsoever while maintaining solution stability.
In total, three simulations are performed. Two cases are run for comparison: one with base Chord (without the HO-ACR method), and the other with the HO-ACR method. The base mesh has 512 × 256 cells. Both cases begin with only the base mesh until t = 5.3 s. Then two AMR levels with a refinement ratio of 2 are added based on the gradients of density and pressure for the lower half of the domain, with a small buffer wrapping around the periodic boundary. For comparison, the third case is run by using AMR on the entire domain from the initial solution as traditional applied, which is considered as the "ideal" case. The goal here is to demonstrate that the HO-ACR method enables us to obtain a quality solution even if imposing AMR in arbitrary space-time and crossing strong discontinuities. This strategy may be performed for embedded-DNS or embedded-LES where only a portion of the space-time mesh is appropriately refined due to computational expense. Figure 20 shows the solution of H 2 O at t = 5.3 s. In the figure, there are four subfigures: (a) with the "ideal" case where AMR has been used for the entire run; (b) showing an initially coarse solution from Chord with refinement added; (c) and (d) comparing the closeups of solutions between Chord without and with the HO-ACR algorithm for the domain regridded on the bottom half, as indicated by the black rectangles in the prior subfigures. After the initial refinement as seen in the subfigure (c), negative species mass fractions, as indicated by magenta, are large, and positive overshoots beyond the coarse solution are shown in yellow-orange. While mass fractions of species must be always positive to be physically meaningful, Chord does have a small tolerance for negative mass fractions to cope with the numerical instability in the nonlinear solvers employed during the solution process. 10 In the subfigure (D), the magnitude of the negative mass fractions is reduced by an order of magnitude, which is sufficient to have kept the overall algorithm stable and allow the solution to develop properly. Note that for intermediate chemistry-dependent operations, Chord first renormalizes the mass fractions to avoid nonphysical quantities. The correction first bounds each species mass fraction between 0 and 1 and then renormalizes so that the sum of all species mass fractions equals one. 26 At the end of each time step, a similar procedure is applied to the conserved solution state except that the bounds are expanded by a small tolerance. We find the correction process handles small nonphysical perturbations from sources such as time integration well, but large errors from grid refinement are too drastic to correct, motivating the need for the HO-ACR method. Overall, reducing nonphysical overshoots at the source (in this case, from interpolation), vastly mitigates the impact of renormalization on the behavior of the dynamical system. Although not the focus of this work, balancing chemistry is an important discussion, and we refer readers to works by Christopher et al. 10 and Owen et al. 27 for details of the approach Chord uses. Figure 21 shows the H 2 O after continuing to evolve the solution to t = 10 s. The layout of this figure is the same as Figure 21. Although not matching the "ideal" solution, the partially added refinement is seen to begin capturing fine scale features not visible in the coarse solution while using significantly less computational resources. This provides an ability to assess the validity of embedded-DNS or -LES, although we refrain from making judgment here. Instead, we focus on the quality of the interpolations. The magnitude of the solution difference with and without the HO-ACR method grows smaller in comparison to that of the earlier time. There is a significant overshoot along the AMR interface that is effectively prevented by the HO-ACR algorithm. Despite what may appear to be minor differences, this is of significant impact for chemically reacting flows with stiff reactions where minor errors can push the solution to unphysical states that are unsolvable. Running further in time, this overshoot would eventually destabilize the solution if the HO-ACR method is not applied.

C 3 H 8 -air combustion in a bluff body combustor
The premixed combustion of C 3 H 8 -air is simulated in a bluff body combustor. 28,29 As shown in Figure 22, the triangular bluff body lies in a rectangular channel. At the inlet, a gaseous mixture of 4.01% C 3 H 8 , 22.36% O 2 , and 73.64% N 2 by mass fraction flows into the domain at a velocity of 15.7 m/s. The mixture has a temperature of 310 K and a pressure of 101,325 P. The flow has a bulk Mach number of 0.053 and bulk Reynolds number of 50,000. The reaction mechanism for C 3 H 8 -air is a stiff system including 25-species and 66-reactions. 30 To ensure stability, only the ACR method is explored for this application and the high-order extension is disabled. We revert to the low order (strictly bound-preserving) ACR method for this case, because the reacting species are highly sensitive to any overshoots allowed by the HO-ACR method.
The initial conditions in the domain are set to the same values as the inlet, with the exception of a small region surrounding the bluff body. Near the bluff body, the gas mixture is initialized to 12% CO 2 , 6.54% H 2 O, 5.14% O 2 , and 73.62% N 2 , with a temperature of 1300 K. This hot spot provides the ignition of the combustion. A base mesh of 11,008 cells are used, with two levels of refinement at a ratio of 2. The second level is fixed behind the bluff body and spans the experimental observation window. The third level dynamically tracks the flame and boundary layers. The bluff body case uses the conforming mapped multiblock technique to fit a structured mesh to the geometry. The initial grid boxes are shown in Figure 23. Species OH and CH 2 O are good flame markers. Therefore, during the run time, cells defining the third level are tagged based on c OH × c CH 2 O > 2 × 10 −9 . Figure 24 shows the CH 2 O contour in the region immediately downstream of the trailing edge of the bluff body. On the left, the solution obtained from base Chord after regridding displays small artifacts arising from interpolation overshoots as highlighted by the violet bubbles. On the right, the presence of these noises in the solution is diminished greatly by instead applying the ACR method. Similarly, Figure 25 compares the OH mass distributions between the base and the ACR algorithms after regridding for the same physical location and solution time. While the differences in OH are barely visible in Figure 25, after 80 time steps, these grow into quite large perturbations as shown in Figure 26. The effectiveness of using the ACR to suppress the interpolation overshoots is clearly seen in Figure 26, as highlighted by gray bubbles. Through the numerical experiments, the ACR algorithm is found to be necessary to ensure the overall simulation stability. With the ACR method, a full simulation of the 3D version of the bluff body case was performed. The z-direction is defined with width of 76.2 mm, and enforced with periodic boundaries. The base mesh is approximately 330,000 cells. One additional level with a refinement factor of 2 is added. Figure 27 shows the instantaneous isosurfaces of the vorticity, temperature, mass fractions of H 2 O and OH, respectively, for the 3D C 3 H 8 -air flame at a solution time of t = 246 ms. Without using AMR, it is much more expensive to capture the flow and flame dynamics near the bluff body as shown in Figure 27. Furthermore, without using the ACR method, this 3D combustion simulation has numerical difficulty where the flame crosses AMR interfaces. One could refine the flame everywhere in the domain to avoid numerical issues from AMR, but this would require significantly more computing resources for little benefit-high fidelity is only required in the experimental observation window.

CONCLUSIONS AND FUTURE WORK
In the present study, the HO-ACR algorithm has been implemented, verified, and validated in a fourth-order finite-volume CFD software solving complex fluid dynamics problems with the presence of chemical reactions, shock waves, and realistic geometries. For such problems, AMR is essential for capturing multiscale flow features with computational efficiency. However, the interpolation operation required by the AMR process often induces unwanted overshoots near discontinuities which destabilizes the overall simulation, particularly for flows with combustion. This is one of the main motivations for the development of the HO-ACR algorithm. The HO-ACR method has been verified to ensure high-order accuracy for smooth flow problems, preserve bounds for nonsmooth ones, and maintain conservation on mapped grids. On challenging problems such as the shock-induced H 2 -O 2 combustion, the HO-ACR method ensures better bounds-preservation in solutions even if a fine level is arbitrarily imposed across or near discontinuities. Furthermore, the application of the method to a C 3 H 8 -air premixed flame in the bluff body combustor has demonstrated that the engagement of HO-ACR is imperative to ensure the solution stability. Note that the chemical kinetics used in this combustion problem, involving 24-species and 66-reactions, is very stiff and the bluff body geometry increases the stiffness through the nonlinear interactions between the fluid and flame dynamics. Future work will be focused on extending applications of the HO-ACR method to the computational geometries represented by general mapped multi-block domains. This would require having some form of bound-preserving property in the interpolation operation across blocks, where the grid metrics may also be discontinuous. Interpolation between blocks requires a somewhat different approach for bounds preservation. Interpolation of ghost cells between different meshes can not reasonably have a conservation constraint, but with careful consideration may use clipping and redistribution for maintaining bounds across different mappings.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.