See pages 1-last of paper.pdf


# A New Re-redistribution Scheme for Weighted State Redistribution with
Adaptive Mesh Refinement

I. Barrio Sanchez  A. S. Almgren  J. B. Bell  M. T. Henry de Frahan  W. Zhang

###### Abstract

State redistribution (SRD) is a recently developed technique for stabilizing
cut cells that result from finite-volume embedded boundary methods. SRD has
been successfully applied to a variety of compressible and incompressible flow
problems. When used in conjunction with adaptive mesh refinement (AMR),
additional steps are needed to preserve the accuracy and conservation
properties of the solution if the embedded boundary is not restricted to a
single level of the mesh hierarchy. In this work, we extend the weighted state
redistribution algorithm to cases where cut cells live at or near a coarse-
fine interface within the domain. The resulting algorithm maintains
conservation and is demonstrated on several two- and three-dimensional example
problems.

###### keywords:

State redistribution , Adaptive Mesh Refinement

††journal: Journal of Computational Physics

[inst1]organization=University of Pittsburgh, addressline=4200 Fifth Ave,
city=Pittsburgh, postcode=15260, state=PA, country=USA

[inst2]organization=Lawrence Berkeley National Laboratory, addressline=1
Cyclotron Rd, city=Berkeley, postcode=94720, state=CA, country=USA

[inst3]organization=National Renewable Energy Laboratory, addressline=15013
Denver W Pkwy, city=Golden, postcode=80401, state=CO, country=USA

##  1 Introduction

Numerical solution of partial differential equations in complex geometric
domains has a wide range of important scientific and engineering applications.
Embedded boundary methods, also known as “cut cell” methods, provide an
attractive approach for these types of problems since mesh generation is
robust and automatic. The basic idea of the embedded boundary / cut-cell
approach is to treat the domain boundary, or solid object within the domain,
as a surface intersecting a regular Cartesian grid. The process of mesh
generation for the embedded boundary approach then reduces to intersections of
the surface with regular finite volume cells.

The difficulty that arises with this process is that cells cut by the
interface may have arbitrarily small volumes. This is referred to as the
“small cell problem” and is one of the major challenges when using the
embedded boundary / cut-cell approaches. A number of different approaches have
been proposed to address the small cell problem. A novel dimensionally split
approach was proposed by Gokhale et al. [1]. $H$-box methods [2] have nice
theoretical properties but have not been extended to three dimensions. Cell
merging, in which small cells are merged with larger cells, is an intuitive
approach that has been successfully used to address the small cell problem,
and has seen significant recent progress. See, for example, [3, 4, 5, 6].

An alternative approach to cell merging that has been successfully used for
simulation of compressible flow [7, 8, 9, 10, 11, 12] and incompressible flow
[13, 14] is called flux redistribution (FRD). In flux redistribution schemes,
a conservative but potentially unstable update for the solution is first
computed using the difference of area-weighted fluxes divided by the cell
volume. The solution in cut cells is then stabilized/regularized by adjusting
the update to the solution to be stable, and redistributing the remainder of
the update to neighboring cells in a conservative manner. This approach is
attractive because it is simple to implement as a postprocessing step.
However, the downside is that there is a loss of accuracy at the cut cells and
the scheme is not linearity preserving.

Recently, the state redistribution algorithm (SRD) first proposed by Berger
and Giuliani [15] has shown great promise as an alternative to flux
redistribution. The state redistribution scheme, like flux redistribution,
starts by computing a conservative but potentially unstable update to the
solution. However, SRD updates the conserved variables and then redistributes
the updated solution, whereas FRD redistributes the update before updating the
conserved variables. State redistribution, like flux redistribution, is fully
conservative and has successfully been used to stabilize both finite volume
[15] and discontinuous Galerkin [16] solvers. Unlike the FRD approach, SRD is
linearity preserving, which leads to improved accuracy at the boundary.
Giuliani et al. [17] extend the original state redistribution algorithm to
three dimensions and introduce a number of algorithmic improvements. In
particular, the authors introduce a weighted version of the algorithm (WSRD)
that is easily generalized and demonstrate that state redistribution can
stabilize both advective and diffusive contributions to the solution update.
The authors also demonstrate that the methodology can be extended to
incompressible flows. Berger and Giuliani [18] introduce an algorithm for
selecting weights that makes WSRD monotone, total variation diminishing and
GKS stable in most situations.

The basic concept of embedded boundary approaches, namely, that they are based
on an underlying regular finite volume grid, makes them a natural candidate
for block-structured adaptive mesh refinement (AMR) of the type initially
proposed by Berger and collaborators [19, 20]. Several of the algorithms
referenced above use this type of AMR. In AMR algorithms for hyperbolic
conservation laws, the data at each level are advanced independently of the
other levels, other than using Dirichlet data provided from coarser grids as
Dirichlet “boundary conditions” where appropriate. When the data on grids at
different levels reach the same time they must be synchronized to construct a
stable, conservative composite solution. When no cut cells are present, this
synchronization has two parts. The first is re-defining the coarse solution as
the average of the fine solution where possible; this changes solution values
only on coarse cells covered by fine cells. The second part is “refluxing”,
whereby the fluxes used to update the coarse solution are effectively over-
written by the time- and space-averaged fine fluxes, thus modifying the coarse
data adjacent to, but not covered by, fine grids. Redistribution algorithms
for embedded boundary methods require additional steps in the synchronization
process when coarse/fine interfaces are close to the embedded boundary. Pember
et al. [7] developed a so-called re-redistribution algorithm that addresses
this additional synchronization for flux redistribution.

The goal of this paper is to extend the weighted state redistribution
algorithm to enable AMR without the restriction that the entire embedded
boundary be resolved at the same level. The focus here is on the AMR
synchronization process when a coarse/fine boundary is near an embedded
boundary. (If the entire embedded boundary is always at the same, typically
finest, level, then no special synchronization is required and the current
algorithm as presented in [17] is sufficient.) We will present the algorithm
for a system of hyperbolic conservation laws. The implementation of the re-
redistribution algorithm is done in the finite-volume, block-structured
adaptive mesh refinement framework of AMReX [21].

The paper is organized as follows. In Section 2, we review the embedded
boundary discretization of conservation laws and the basics of AMR. In Section
3, we review the weighted state redistribution (WSRD) algorithm and describe
how to do re-redistribution for WSRD. In Section 4, we validate our numerical
scheme using canonical two- and three-dimensional example problems and present
an application of the methodology to a more realistic flow problem. Finally,
we summarize our conclusions in Section 5. The appendices contain the details
of the flux redistribution algorithm (Appendix A), the re-redistribution
algorithm for flux redistribution (Appendix B) and the necessary modifications
to the Godunov advance in the presence of cut cells (Appendix C).

##  2 Mathematical preliminaries

In this section we review the basic concepts of embedded boundary
discretizations, discuss the redistribution approach to addressing the small
cell problem, and introduce adaptive mesh refinement (AMR). For the discussion
here, we consider a system of hyperbolic conservation laws in three
dimensions,

$\displaystyle\frac{\partial U}{\partial t}+\nabla\cdot F(U)=0$ (1)

where $U$ is the state variable and $F(U)=(F^{x}(U),F^{y}(U),F^{z}(U)\;)$ is
the flux.

###  2.1 Notation

Embedded boundary (EB) methods are used to discretize partial differential
equations (PDEs) in complex domains. The basic idea is to represent an
irregular boundary (between “fluid” and “body”, where we wish to solve the
PDEs in the fluid region only) by intersecting an explicitly specified
boundary with a uniform finite-volume Cartesian grid with mesh spacings
$\Delta x$, $\Delta y$ and $\Delta z$. This representation introduces
irregularly shaped cells only adjacent to the boundary. Following standard
notation, we define each grid cell $(i,j,k)$ to be either body, cut, or
regular. We define the geometric volume fraction, $\Lambda$, of each cell to
be the fraction of that rectangular cell volume that is inside the fluid
region: body cells have $\Lambda=0,$ regular cells have $\Lambda=1$, and for
cut cells $0<\Lambda<1.$ Area fractions $a$ are stored on each cell face,
again with values in $[0,1],$ representing the fraction of the face not
covered by the body. In the AMReX implementation used here, the EB information
needed to represent the geometry, such as face area fractions and cell volume
fractions, is precomputed and stored in a distributed database at the
beginning of the calculation. For additional details on the embedded boundary
implementation, we refer to the AMReX documentation https://amrex-
codes.github.io/amrex/docs_html/.

Going forward, we will use the notation $A$ to represent the face area of a
regular cell multiplied by the area fraction of that particular face, for
example $A^{x}_{i+\mathchoice{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},j,k}=a_{i+\mathchoice{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},j,k}\Delta y\Delta z$, and $V$ to represent
the volume of a regular cell multiplied by the volume fraction of that
particular cell, i.e. $V_{i,j,k}=\Lambda_{i,j,k}\Delta x\Delta y\Delta z$.

###  2.2 Single level algorithm

For the discussion here, we will consider a single-step integration scheme
typical of higher-order Godunov-type discretizations. Generalizations to a
method of lines approach are discussed at the end of Sec. 3. We let
${U}_{i,j,k}$ and $F^{x}_{i+\mathchoice{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},j,k}$, $F^{y}_{i,j+\mathchoice{\kern
0.0pt\raise 1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},k},$ and $F^{z}_{i,j,k+\mathchoice{\kern
0.0pt\raise 1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2}}$ represent the discretized solution $U$,
and the normal flux $F$ on $x$-faces, $y$-faces and $z$-faces, respectively.
We define a standard finite volume scheme in the form,

$\displaystyle U_{i,j,k}^{n+1}$ $\displaystyle=$ $\displaystyle
U_{i,j,k}^{n}+\Delta t\;\delta U^{c}_{i,j,k}$

where

$\displaystyle\delta U^{c}_{i,j,k}=$ $\displaystyle-$
$\displaystyle\frac{\left(A^{x}_{{i+\mathchoice{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},j,k}}F^{x}_{{i+\mathchoice{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},j,k}}-A^{x}_{{i-\mathchoice{\kern
0.0pt\raise 1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},j,k}}F^{x}_{{i-\mathchoice{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},j,k}}\right)}{V_{i,j,k}}$ $\displaystyle-$
$\displaystyle\frac{\left(A^{y}_{{i,j+\mathchoice{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},k}}F^{y}_{{i,j+\mathchoice{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},k}}-A^{y}_{{i,j-\mathchoice{\kern
0.0pt\raise 1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},k}}F^{y}_{{i,j-\mathchoice{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},k}}\right)}{V_{i,j,k}}$ $\displaystyle-$
$\displaystyle\frac{\left(A^{z}_{{i,j,k+\mathchoice{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2}}}F^{z}_{{i,j,k+\mathchoice{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2}}}-A^{z}_{{i,j,k-\mathchoice{\kern
0.0pt\raise 1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2}}}F^{z}_{{i,j,k-\mathchoice{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2}}}\right)}{V_{i,j,k}}$ $\displaystyle-$
$\displaystyle\frac{A^{f}_{i,j,k}F^{f}_{i,j,k}}{V_{i,j,k}}\;.$

The update is conservative by construction, and the time accuracy of the
scheme is determined by the details of the flux construction. Here,
$A^{f}_{i,j,k}$ represents the area of the boundary intersected with cell
(${i,j,k}$) and $F^{f}_{i,j,k}$ represents the flux in the normal direction at
the boundary. We note that $A^{f}_{i,j,k}$ and the unit normal to the surface
$\mathbf{n}^{f}_{i,j,k}$ (pointing from body to fluid) can be approximated
using

$\displaystyle A^{f}_{i,j,k}\mathbf{n}^{f}_{i,j,k}$
$\displaystyle=(A^{x}_{{i+\mathchoice{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},j,k}}-A^{x}_{{i-\mathchoice{\kern
0.0pt\raise 1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},j,k}})\;\mathbf{i}+\;(A^{y}_{{i,j+\mathchoice{\kern
0.0pt\raise 1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},k}}-A^{y}_{{i,j-\mathchoice{\kern
0.0pt\raise 1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2},k}})\;\mathbf{j}+\;(A^{z}_{{i,j,k+\mathchoice{\kern
0.0pt\raise 1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2}}}-A^{z}_{{i,j,k-\mathchoice{\kern
0.0pt\raise 1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\sevenrm1}\kern-1.04996pt/\kern-1.04996pt\lower
0.75346pt\hbox{\sevenrm2}}{\kern 0.0pt\raise
1.50694pt\hbox{\fiverm1}\kern-1.39998pt/\kern-1.04996pt\lower
0.75346pt\hbox{\fiverm2}}{1\\!/2}}})\;\mathbf{k}\;.$

This formula is exact for planar surfaces and gives second-order averages for
smooth surfaces.

In the absence of cut cells, where all cells have the same volume $V=\Delta
x\Delta y\Delta z$, the CFL condition for an explicit update constrains
$\Delta t$ based on the ratio of cell volume $V$ to face areas $A$ as well as
the maximum wave speed in the problem, i.e. the time step depends linearly on
the mesh spacing. The “small cell problem” discussed in the introduction
refers to the fact that when an embedded boundary representation of geometry
is used, the cut cells can have arbitrarily small volumes, making the update
in Eq. (2.2) unstable if the time step calculation does not take into account
the reduction in cell volume [13]. Reducing the time step to account for
arbitrarily small cell volumes would be too computationally expensive, thus as
discussed above there are a number of approaches to addressing the small cell
problem, all of which allow the computation of $\Delta t$ based only on the
full-cell CFL constraint. Recalling that the FRD scheme modifies the update
itself, we can write the single-level single-timestep advance using FRD as

$\displaystyle U_{i,j,k}^{n+1}$ $\displaystyle=$ $\displaystyle
U_{i,j,k}^{n}+\Delta t\;R^{FRD}(\delta U^{c}_{i,j,k})$ (3)

where $R^{FRD}$ represents the flux redistribution operator. Details of the
flux redistribution algorithm are discussed in A. Analogously, since WSRD
operates on the full solution, we can write

$\displaystyle U_{i,j,k}^{n+1}$ $\displaystyle=$ $\displaystyle
R^{WSRD}(U_{i,j,k}^{n}+\Delta t\;\delta U^{c}_{i,j,k})$ (4)

where $R^{WSRD}$ represents the state redistribution operator, which will be
described in detail in Section 3.

###  2.3 Multilevel algorithm

We consider a block-structured adaptive mesh refinement algorithm (AMR) in
which regions at a given level requiring additional resolution are covered
with grid patches of finer resolution. (For our purposes, a “level” refers to
all the grids with the same resolution, and higher levels denote finer
resolution. We call two levels “adjacent” if they are not separated by any
intermediate levels.) With this type of refinement approach, one can think of
the time stepping algorithm as a recursive procedure, starting with the
coarsest level. Data on coarse grids are advanced first, then data from the
coarse grids are interpolated to provide Dirichlet boundary conditions for the
grids at the next finer level (if there is one). We refer to boundary cells
for fine grids that are filled from coarse data as “ghost” cells. Data on the
finer grids are then advanced, possibly with smaller time steps, until they
reach the same time at the coarse grid. Advancing the fine grids with smaller
time steps than the coarse grid is referred to as subcycling. If there are
still finer levels, this process repeats recursively. Once the data on grids
at adjacent levels reach the same time, the data on the two levels are
synchronized.

There are a number of different facets to an AMR algorithm including
refinement criteria, data distribution and load balancing, that we will not
consider here. Here we will focus on how the data at adjacent levels are
synchronized to create a composite solution. We define the refinement ratio
$r_{\ell}$ to be the ratio of the mesh spacing at level $\ell$ to that at
level $(\ell+1)$, and let $n_{\ell}$ be the number of level $\ell+1$ time
steps corresponding to a single level $\ell$ time step. To advance level
$\ell,$ $0\leq\ell\leq\ell_{\mathrm{m}ax}$ the following steps are taken.

Starting with level $\ell=0,$

  * •

Fill boundary conditions for $U$ from level $\ell-1$ if level $\ell>0$, and
from the physical domain boundary conditions.

  * •

Define $U^{\ell,n+1}$ using Eq. (3) or Eq. (4) above.

  * •

Then, if $\ell<\ell_{max}$

    * –

Advance $U$ at level ($\ell+1$) for $n_{\ell}$ time steps with time step
$\Delta t^{\ell+1}=\frac{1}{n_{\ell}}\Delta t^{\ell}.$

    * –

Synchronize the data between levels $\ell$ and $\ell+1$

###  2.4 Synchronization

As mentioned in the introduction, there are several steps in the
synchronization procedure. First, we want to replace the coarse solution by
the average of the fine solution where possible. Here we define a “covered”
cell at level $\ell$ as one which is overlaid by $r_{\ell}^{d}$ level
$(\ell+1)$ grid cells, where $d$ is the spatial dimension of the problem. In
addition we define a coarse cell to be “at the coarse-fine boundary” if it
shares a face with fine cells but is not covered, and a fine cell to be “at
the coarse-fine boundary” if it shares a face with a coarse cell that is not
covered. These definitions are illustrated in Figure 1. For EB methods we
replace the solution in covered coarse grid cells with the volume-fraction-
weighted average of the fine grid solution. We refer to this step as the
average down step.

Figure 1: Schematic of the different type of cells and levels close to an EB
(black line). Blue grid lines: coarse level, red grid lines: fine level.
Orange (1) fine ghost cells; purple (2) uncovered coarse cells; burgundy (3)
valid fine cells; and magenta (4) covered coarse cells. The coarse-fine
interface lies between A and B. Note that we have not labeled all of the cells
to the left of $\overline{AB}$; all of the coarse cells left of
$\overline{AB}$ are covered and all of the fine cells left of $\overline{AB}$
are valid. Here, cells with a face that is part of $\overline{AB}$ are
considered to be at the coarse-fine boundary. Shaded blue cells indicate the
coarse cells at the coarse-fine interface; shaded red cells denote the fine
cells at the coarse-fine interface.

The other steps in the synchronization address issues that arise from
advancing the levels independently. Although the operations are defined
procedurally, we note that these operations, in combination with the level-by-
level update, implicitly define a composite discretization. (See [22] for a
more general discussion of this perspective.) One of these issues is that at
the coarse-fine boundary, the flux computed on the coarse grid does not match
the average of the fine grid fluxes. (We note that the average of the fine
grid area fractions does in fact equal the coarse grid area fraction, but the
values of the flux itself need not match.) Correcting for this mismatch is
referred to as refluxing. We compute the difference between the coarse grid
flux and the time and space average of the fine grid fluxes on every face
along the coarse-fine boundary. This difference is used to correct the
solution in the coarse grid cells at the coarse-fine boundary, i.e. coarse
cells adjacent to the faces where the fluxes are corrected but not covered by
finer grid (Figure 1).

For example, if coarse cell $(I,J,K)$ is located at the coarse-fine boundary
and cell $(I-1,J,K)$ is covered, then the refluxing step updates
$U_{I,J,K}^{\ell}$ by adding

$\displaystyle\frac{1}{V_{I,J,K}}\delta F^{x,\ell}_{I-\frac{1}{2},J,K}$

where

$\delta F^{x,\ell}_{I-\frac{1}{2},J,K}\equiv-\Delta
t^{\ell}A^{x,\ell}_{I-\frac{1}{2},J,K}F^{x,\ell}_{I-\frac{1}{2},J,K}+\sum_{n=1,n_{\ell}}\sum_{j,k}\Delta
t^{\ell+1}A^{x,\ell+1}_{i+\frac{1}{2},j,k}F^{x,\ell+1}_{i+\frac{1}{2},j,k,n}$
(5)

represents the extensive flux correction associated with edge
$({I-\frac{1}{2},J,K})$. Here $i$ is the index for fine cells at the
coarse/fine boundary and the summation in $j$ and $k$ ranges over the fine
grid cells at the coarse/fine boundary that cover cell $(I-1,J,K)$ and $n$
refers to substeps of the fine grid corresponding to a coarse grid time step.
Note that Eq. (5) incorporates the orientation of the face; if the coarse cell
at the boundary is instead on the left side of the fine grid the signs of the
flux terms are reversed.

For non-EB problems, and when cut cells are sufficiently far from coarse/fine
boundaries, refluxing is the only synchronization that is needed to make the
composite solution conservative. However, when there are cut cells near the
coarse/fine boundary, additional synchronization is required in order to
account for the quantity of each conserved variable that has effectively
crossed the coarse-fine boundary during the redistribution step which occurred
independently at the coarse and fine levels. This additional redistribution
synchronization is analogous to the refluxing procedure in that during the
redistribution step at each level we accumulate the mismatch, then once the
coarse and fine data have reached the same time, we modify the coarse data at
the coarse/fine boundary to correct the mismatch.

While this process is sufficient to make the composite solution conservative,
there is an additional step required to maintain stability. Specifically,
since the coarse cells whose data is modified by refluxing and the
redistribution synchronization step can be arbitrarily small, we need to
redistribute the corrections from these procedures. Details of these two
steps, which collectively are referred to as re-redistribution, were first
discussed in Pember et al. [7] in the context of flux redistribution; see B
for a summary of this approach.

##  3 Weighted state redistribution scheme

State redistribution [15, 17] has recently emerged as an alternative to flux
redistribution for addressing the small cell problem associated with embedded
boundary methods. State redistribution addresses the problem by creating
logical merging neighborhoods where small cut cells are merged with their
neighbors until the neighborhood volume is greater than a specified threshold.
Unlike traditional cell-merging algorithms, a single cell can be in multiple
neighborhoods. As with FRD, the first step in advancing the solution in time
is to create a conservative but possibly unstable update for each of the state
variables using a finite volume formulation as in Eq. (2.2). The solution is
then provisionally updated using the possibly unstable update, and this new-
time provisional solution is modified by the state redistribution procedure.
This procedure itself has two steps: first, the values of the provisionally
updated state variables are used to define values associated with each
neighborhood. In the second step, the neighborhood values are used to define
new solution values on the original cells. We note that solution values in
both cut and regular cells are used to define the neighborhood solutions, and
that, as with FRD, the SRD procedure modifies the final solution in both cut
and regular cells. Unlike FRD however, the use of linear interpolation in this
second step of SRD preserves linearity of the full solution where possible.

For a wide range of applications, it would be desirable to couple state
redistribution with adaptive mesh refinement. Current (W)SRD implementations
impose a restriction that cells on or near the embedded boundary must all be
at the same level of refinement; a coarse/fine boundary cannot intersect the
embedded boundary. Here we will introduce a “re-redistribution” algorithm for
WSRD that removes this restriction. First, we review how the neighborhoods are
defined and how the solution is advanced.

###  3.1 Defining the neighborhoods

To create merging neighborhoods we first identify “small” cells, i.e. cells
with a volume fraction less than a specified target threshold, $V_{target}$.
These cells are then logically merged with neighboring cells until the
neighborhood volume, i.e. the sum of the volume of all cells in the
neighborhood reaches the target value. The neighborhoods in WSRD do not form a
tessellation of the computational domain as is done in traditional cell
merging. Instead they are logical constructs where cells can potentially be
part of several neighborhoods. Here we use

$\displaystyle V_{target}=\Delta x\Delta y\Delta z/2.$

This choice for the target threshold is informed by [15], but in practice the
appropriate threshold may depend on factors such as the details of the
numerical scheme and boundary conditions.

There are several ways to choose which neighbor a small cell merges with.
Using the simplest option, normal merging, we first merge small cells with
their (cut or regular) neighbor in the direction of the boundary normal. If
the addition of the neighbor’s volume is not sufficient for the target
threshold, additional cells are added to the merging neighborhood. In our
implementation, our second choice neighbor is based on the next largest
component of the normal. This creates an “L-shaped” neighborhood and we
automatically merge the cell in the same plane as the two other neighbors to
turn this neighborhood into a $2\times 2$ neighborhood. In 3D if the threshold
target is still not met, we add the additional cells necessary to create a
$2\times 2\times 2$ neighborhood. An alternative process, called central
merging, merges a small cell with all of its (cut or regular) neighbors in the
$3\times 3\times 3$ box centered on the small cell.

The main motivation for WSRD (as opposed to the initial formulation of SRD) is
that a large cut cell that is just below the threshold volume requires less
stabilization than a very small cell. Observing that every cell is always in
its own neighborhood, we define $M_{i,j,k}$ to be the set of cells contained
in the neighborhood associated with cell $(i,j,k)$, and
$M^{-}_{i,j,k}=M_{i,j,k}-\\{(i,j,k)\\}$, i.e. $M^{-}_{i,j,k}$ is the set of
cells contained in the neighborhood of cell $(i,j,k)$ with the exception of
cell $(i,j,k)$ itself. We define $W_{i,j,k}$ to be the set of indices
$(r,s,t)$ such that cell $(i,j,k)$ is in the neighborhood associated with
$(r,s,t)$, i.e. so that cell $(i,j,k)$ is in $M_{r,s,t}.$ We define
$W_{i,j,k}^{-}=W_{i,j,k}-\\{(i,j,k)\\}$. Finally, we define $N_{i,j,k}$ to be
the number of neighborhoods that cell $(i,j,k)$ is in, i.e. the size of
$W_{i,j,k}$.

For the weighted version of SRD, at every regular or cut cell $(i,j,k)$ we
define

${\beta}_{i,j,k}=\begin{cases}(V_{target}-V_{i,j,k})\;\;/\displaystyle\sum_{(r,s,t)\in
M^{-}_{i,j,k}}\,V_{r,s,t}\,&\text{if $V_{i,j,k}<V_{target}$}\\\ 0&\text{if
$V_{i,j,k}\geq V_{target}$}\end{cases}$ (6)

and then

${\alpha}_{i,j,k}=1\,-\,\frac{1}{N_{i,j,k}}\,\sum_{(r,s,t)\in
W^{-}_{i,j,k}}\,\beta_{r,s,t}\,.$ (7)

Equations (6) and (7) imply that $0\leq\alpha_{i,j,k},\beta_{i,j,k}\leq 1$. We
observe that $\beta_{i,j,k}$ depends on the volume fractions of the cells in
the neighborhood of $(i,j,k),$ while $\alpha_{i,j,k}$ depends on the values of
$\beta$ in cells whose neighborhoods contain $(i,j,k).$

As $V_{i,j,k}$ becomes larger, the weighted algorithm increases the dependence
of the final solution in cell $(i,j,k)$ on the provisional solution in
$(i,j,k)$, unlike the original SRD algorithm. This weighted SRD algorithm may
have somewhat different stability properties than the original SRD algorithm,
especially for one-dimensional test cases that are not representative of
Cartesian cut cell meshes in higher dimensions.

We define the weighted volume, $\widehat{V}_{i,j,k}$, of the neighborhood
associated with $({i,j,k})$ by

${\widehat{V}}_{i,j,k}=\alpha_{i,j,k}V_{i,j,k}+\beta_{i,j,k}\sum_{(r,s,t)\in
M^{-}_{i,j,k}}\,\frac{V_{r,s,t}}{N_{r,s,t}}.$ (8)

and the centroid of the neighborhood associated with $(i,j,k)$ by

$\displaystyle\begin{aligned}
({\widehat{x}}_{i,j,k},{\widehat{y}}_{i,j,k},{\widehat{z}}_{i,j,k})=\frac{1}{\widehat{V}_{i,j,k}}\biggl{(}&\alpha_{i,j,k}V_{i,j,k}(x_{i,j,k},y_{i,j,k},z_{i,j,k})\,+\\\
&\;\beta_{i,j,k}\sum_{(r,s,t)\in
M^{-}_{i,j,k}}\frac{V_{r,s,t}}{N_{r,s,t}}(x_{r,s,t},y_{r,s,t},z_{r,s,t})\biggr{)}.\end{aligned}$
(9)

where $(x_{i,j,k},y_{i,j,k},z_{i,j,k})$ is the original centroid of cell
$(i,j,k)$.

In an efficient implementation, all of the steps in this subsection can be
done when the EB information is defined at the start of a simulation or, in
the case of AMR, after any step in which changes in the grid hierarchy would
require it.

###  3.2 Advancing the solution

As noted in Sec. 2.1, we start by defining a conservative update, $\delta
U_{i,j,k}^{c}$ using Eq. (2.2). Here for simplicity we initially focus on a
single step integration scheme, in particulate a piecewise-linear variant of
Godunov’s method of Collela [23] (2D) and Saltzman [24] (3D) that has been
modified in and near cut cells. The details of how we modify the algorithm are
discussed in C. (Berger and Giuliani [15] use a similar method in two
dimensions that uses a somewhat different approach to dealing with cut cells.)

For WSRD, once we have computed $\delta U_{i,j,k}^{c}$ we define a provisional
update

$\widehat{U}_{i,j,k}=U^{n}_{{i,j,k}}+\Delta t\;\delta U_{i,j,k}^{c}$ (10)

The next step is to define the weighted solution average for each neighborhood
as

$\widehat{Q}_{i,j,k}=\frac{1}{{\widehat{V}}_{i,j,k}}\,\left(\alpha_{i,j,k}V_{i,j,k}\widehat{U}_{i,j,k}\,+\beta_{i,j,k}\,\sum_{(r,s,t)\in
M^{-}_{i,j,k}}\frac{V_{r,s,t}}{N_{r,s,t}}\widehat{U}_{r,s,t}\right).$ (11)

We observe that if $(i,j,k)$ is a regular cell, or a cut cell with
$V_{{i,j,k}}>V_{target},$ the neighborhood associated with $(i,j,k)$ contains
only itself, and $\widehat{Q}_{i,j,k}=\widehat{U}_{i,j,k}.$

To achieve second order accuracy where possible, in both the original and
weighted SRD approaches, we define a gradient,
$(\widehat{\sigma}_{x},\widehat{\sigma}_{x},\widehat{\sigma}_{x})$ in all
neighborhoods. To do so, we define a $3\times 3\times 3$ ($3\times 3$ in 2D)
box centered on $(i,j,k)$ for all $(i,j,k)$ whose neighborhoods contain at
least two cells; otherwise we simply set the gradient to zero. (Note that this
box has no relation to the size or shape of the merging neighborhoods
themselves.) We then use the least squares equation to fit the values of
$\widehat{Q}$ in all cut or regular cells in this box, treating the value of
each $\widehat{Q}_{r,s,t}$ as defined at the centroid of that neighborhood,
$(\widehat{x}_{r,s,t},\widehat{y}_{r,s,t},\widehat{z}_{r,s,t})$. If the
distance between centroids in any direction is below a given threshold (half
the mesh spacing in that direction), then the gradient stencil is increased to
have width $5$ rather than $3$ in that direction. If the increased stencils
either don’t contain enough points or don’t have sufficient separation of
coordinates, we set the gradient to 0. In order to not generate new maxima or
minima of $\widehat{Q},$ we also limit the components of the gradient with a
Barth-Jesperson-style limiter. This allows us to define a linear function for
neighborhood $(i,j,k)$ of the form

$\widehat{q}_{i,j,k}(x,y,z)=\widehat{Q}_{i,j,k}+\widehat{\sigma}_{x,i,j,k}(x-\widehat{x}_{i,j,k})+\widehat{\sigma}_{y,i,j,k}(y-\widehat{y}_{i,j,k})+\widehat{\sigma}_{z,i,j,k}(z-\widehat{z}_{i,j,k}).$
(12)

We then compute the final solution,

$U^{n+1}_{i,j,k}=\alpha_{i,j,k}\;\widehat{q}_{i,j,k}(x_{i,j,k},y_{i,j,k},z_{i,j,k})+\frac{1}{N_{i,j,k}}\sum_{(r,s,t)\in
W^{-}_{i,j,k}}\beta_{r,s,t}\;\widehat{q}_{r,s,t}(x_{i,j,k},y_{i,j,k},z_{i,j,k}).$
(13)

We observe that if a small cut cell $(i,j,k)$ is in no neighborhoods but its
own, then $\alpha_{i,j,k}=1$ and
${U}^{n+1}_{i,j,k}=\widehat{q}_{i,j,k}(x_{i,j,k},y_{i,j,k},z_{i,j,k})$
regardless of how many other cells are in its neighborhood.

###  3.3 Matrix form

We can also express the WSRD algorithm in matrix form, which will be useful in
defining the re-redistribution algorithm. We first define the matrix $A$ such
that

$\widehat{V}=AV,$

where $V$ is the vector of cell volumes and $\widehat{V}$ is the vector of
neighborhood volumes. This is the matrix form of Equation (8). Thus, a nonzero
entry $a_{p,q}$ of $A$ corresponds to the weighted contribution of cell-$p$ to
neighborhood-$q$. (Tying back to our previous notation, cell-$p$ might be
$(i,j,k$) and cell-$q$ might be cell $(r,s,t)$ in 3D.) We note that the matrix
$A$ can be stored in a sparse format. Specifically, for each neighborhood we
only need the values corresponding to the nonzero entries and an associated
index.

The definition of the $\alpha$’s and the $\beta$’s guarantees that the nonzero
elements of $A$ are positive. An additional property of $A$ is that the
columns sum to 1, i.e.,

$e^{T}A=e^{T},$

where $e$ is a vector of all 1’s. This property ensures that the WSRD
algorithm is conservative. Although not pursued here, this suggests that any
collection of non-negative weights that reach the target neighborhood volumes
with $e^{T}A=e^{T}$ can be used to define the SRD algorithm. This opens the
possibility of more sophisticated approaches, including ones that might depend
on the local solution.

We can now write Equation (11) as

$\widehat{Q}=\mathrm{Diag}(\widehat{V})^{-1}A\;\mathrm{Diag}(V)\widehat{U}.$

Furthermore, if the gradients are zeroed in Eq. (12) then we have simply

$U^{n+1}=A^{T}\widehat{Q}.$

Including the gradients in Eq. (12) redistributes mass to preserve linearity
but does not alter the conservation properties of the method. Continuing, we
can write (13) as

$\displaystyle
U^{n+1}=A^{T}\widehat{Q}+\left[\mathrm{Diag}(x)A^{T}-A^{T}\mathrm{Diag}(\widehat{x})\right]\widehat{\sigma}_{x}$
$\displaystyle+\left[\mathrm{Diag}(y)A^{T}-A^{T}\mathrm{Diag}(\widehat{y})\right]\widehat{\sigma}_{y}$
$\displaystyle+\left[\mathrm{Diag}(z)A^{T}-A^{T}\mathrm{Diag}(\widehat{z})\right]\widehat{\sigma}_{z},$

where $(x,y,z)$ are the centroids of the original cells,
$\widehat{x},\widehat{y},\widehat{z}$ are vectors of the components of the
weighted centroids of the neighborhoods and
$\widehat{\sigma}_{x},\widehat{\sigma}_{y},\widehat{\sigma}_{x}$ are vectors
of the components of the gradients from the neighborhood gradient
reconstruction discussed above.

###  3.4 Coupling to AMR

As discussed earlier, when cut cells are sufficiently near to the coarse/fine
boundary, the synchronization operations for AMR must include refluxing the
redistribution as well as redistributing the reflux in order to define a
conservative and stable composite solution. Borrowing some of the notation
from [7], which described re-redistribution for algorithms using flux
redistribution, here we present the re-redistribution algorithm for algorithms
using weighted state redistribution.

We first briefly introduce some additional notation that will be helpful: fine
valid cells are cells inside level $(\ell+1)$ grids, and fine ghost cells are
fine cells that are not inside level $(\ell+1)$ grids but are adjacent to fine
valid cells. Fine ghost cells overlay coarse uncovered cells at the
coarse/fine boundary.

We recall that in order to perform the standard refluxing operation at the
coarse/fine boundary we accumulate the difference of coarse and fine fluxes in
face-centered $\delta F$’s as shown in Eq. 5. In order to reflux the
redistribution, we analogously compute cell-centered $\delta R$’s that track
the extensive amount of each variable that effectively crosses the coarse/fine
interface in the redistribution step at each level. Because the data
associated with a fine ghost cell or a coarse covered cell is not part of the
composite solution, contributions to or from those cells are not correctly
accounted for in the composite solution. When redistributing at the coarse
level, $\ell$, we include in $\delta R$ the extensive contributions that went

  1. (A)

from coarse uncovered cells to coarse covered cells

  2. (B)

from coarse covered cells to coarse uncovered cells

When redistributing at the fine level, $\ell+1,$ we include in $\delta R$ the
extensive contributions that went

  1. (a)

from fine valid cells to fine ghost cells

  2. (b)

from fine ghost cells to fine valid cells

####  3.4.1 Computing $\delta R$

Here we describe how to compute $\delta R$ for the WSRD algorithm. To do so,
we first create a matrix, $\mathcal{R},$ that captures the contributions of
$\widehat{U}$ to $U^{n+1}$ in nearby cells.

If we set the gradient terms in Eqs. (12)-(13) to zero, from the matrix form
of WSRD we can write the final update (in terms of $\widehat{U}$) as

$\mathrm{Diag}(V)U^{n+1}=\mathrm{Diag}(V)A^{T}\mathrm{Diag}(\widehat{V})^{-1}A\;\mathrm{Diag}(V)\widehat{U}$

From this relation we can see that the matrix

$\mathcal{R}\equiv\mathrm{Diag}(V)A^{T}\mathrm{Diag}(\widehat{V})^{-1}A\;\mathrm{Diag}(V)\mathrm{Diag}(\widehat{U})$

specifies how values of $\widehat{U}$ affect the final update. In particular,
the row of $\mathcal{R}$ corresponding to $(I,J,K)$ specifies, in extensive
form, the contribution of each of the $\widehat{U}$’s to $U^{n+1}_{I,J,K}$. We
note that $\mathcal{R}$ is a sparse matrix with the same sparsity pattern as
$A^{T}A$, making the nonzeros easy to identify using the sparse representation
of $A$.

Let us consider the contributions from the redistribution operation at each
level. We can define $\mathrm{R}_{I,J,K}$ and $\mathrm{C}_{I,J,K}$ to be the
row and column of in $\mathcal{R}$, respectively, corresponding to cell
$(I,J,K)$. At the coarse level, we associate with uncovered cell $(I,J,K)$ the
sum of its contributions to covered cells minus the sum of contributions of
covered cells to $(I,J,K)$; i.e., for uncovered cell $(I,J,K)$ we define

$\delta
R_{I,J,K}^{\ell}=\sum_{(II,JJ,KK)\epsilon\mathcal{C}^{\ell}}\mathcal{R}^{\ell}_{\mathrm{R}_{II,JJ,KK},\mathrm{C}_{I,J,K}}-\;\mathcal{R}^{\ell}_{\mathrm{R}_{I,J,K},\mathrm{C}_{II,JJ,KK}}$
(14)

where $\mathcal{C}^{\ell}$ is the set of covered cells at level $\ell$ and
$\mathcal{R}^{\ell}$ is matrix $\mathcal{R}$ at level $\ell$. Here the first
and second terms on the right hand side correspond to contributions of type
(a) and type (b), respectively. Similarly, at the fine level, we associate
with ghost cell $(i,j,k)$ the sum of contributions to cell $(i,j,k)$ from
valid cells minus the sum of contributions of cell $(i,j,k)$ to valid cells;
i.e.,

$\delta
R_{i,j,k}^{\ell+1}=\sum_{(ii,jj,kk)\epsilon\mathcal{V}^{\ell+1}}\mathcal{R}^{\ell+1}_{\mathrm{R}_{i,j,k},\mathrm{C}_{ii,jj,kk}}-\;\mathcal{R}^{\ell+1}_{\mathrm{R}_{ii,jj,kk},\mathrm{C}_{i,j,k}}$
(15)

where $\mathcal{V}^{\ell+1}$ is the set of valid cells at level $\ell+1$ and
$\mathcal{R}^{\ell+1}$ is the matrix $\mathcal{R}$ at level $\ell+1$. In this
case, the first and second terms on the right hand side correspond to
contributions of type (c) and type (d), respectively. We note that in both Eq.
(14) and Eq. (15), $\mathcal{R}$ is sparse with a known sparsity pattern,
making the summations easy to compute.

We now need to describe how to treat the gradient terms in Eqs. (12)-(13). One
possible approach would be to capture explicitly how the gradient depends on
the $\widehat{U}$’s; however, limiting makes this a nonlinear process and the
potentially large stencil for the gradient dramatically enlarges the stencil
for the redistribution. Alternatively, we will view the gradients in Eq. (13)
as fixed and then let

$\delta\widehat{q}^{I,J,K}_{R,S,T}=\widehat{\sigma}_{x,R,S,T}(x_{I,J,K}-\widehat{x}_{R,S,T})+\widehat{\sigma}_{y,R,S,T}(y_{I,J,K}-\widehat{y}_{R,S,T})+\widehat{\sigma}_{z,R,S,T}(z_{I,J,K}-\widehat{z}_{R,S,T})$

be the contribution of neighborhood $(R,S,T)$ to the final update of cell
$(I,J,K)$ as in Eq. (13). The issue here is that we cannot explicitly link the
correction back to the $\widehat{U}$’s. One could potentially associate this
update directly with the neighborhood. This would correspond to treating the
entire $\delta q$ correction for neighborhood $(R,S,T)$ as originating from
cell $(R,S,T)$. Although such an approach would be conservative, it could
potentially generate large corrections for cells with very small volumes,
which would introduce an instability. Instead we define

$\left(\delta
q^{I,J,K}_{R,S,T}\right)_{II,JJ,KK}=\frac{1}{\widehat{V}_{R,S,T}}A_{\mathrm{R}_{R,S,T},\mathrm{C}_{II,JJ,KK}}\;V_{II,JJ,KK}\;\delta\widehat{q}^{I,J,K}_{R,S,T}\quad,$

which distributes $\delta\widehat{q}^{I,J,K}_{R,S,T}$ over the cells
$(II,JJ,KK)$ that are in neighborhodd $(R,S,T)$. If $(I,J,K)$ is uncovered and
$(II,JJ,KK)$ is covered then we add

$V_{II,JJ,KK}A_{\mathrm{R}_{R,S,T},\mathrm{C}_{II,JJ,KK}}\left(\delta
q^{II,JJ,KK}_{R,S,T}\right)_{I,J,K}-V_{I,J,K}A_{\mathrm{R}_{R,S,T},\mathrm{C}_{I,J,K}}\left(\delta
q^{I,J,K}_{R,S,T}\right)_{II,JJ,KK}$ (16)

to $\delta R_{I,J,K}^{\ell}$. Note that in Eq. (16) the indices of $A$
correspond to multiplication by $A^{T}$. Analogous corrections are also added
to $\delta R_{i,j,k}^{\ell+1}$.

Similar to FRD, we associate the WSRD redistribution correction with the
uncovered coarse cells at the coarse/fine boundary. Specifically, for
$(I,J,K)$ at the coarse/fine boundary, we define

$\delta\mathbf{R}^{\ell}_{I,J,K}=\delta R^{\ell}_{I,J,K}+\sum_{i,j,k}\delta
R^{\ell+1}_{i,j,k}$

where the sum is over fine grid ghost cells $(i,j,k)$ that overlay cell
$({I,J,K})$. The quantity $\delta\mathbf{R}^{\ell}_{I,J,K}$ represents the
difference between redistribution on the coarse and fine levels in extensive
form.

####  3.4.2 Updating the solution

Once we have computed $\delta\mathbf{R}^{\ell}_{I,J,K},$ if
$\Lambda_{I,J,K}=1$ then we can simply update

$U^{\ell}_{I,J,K}:=U^{\ell}_{I,J,K}+\frac{1}{V_{I,J,K}}\,\delta\mathbf{R}^{\ell}_{I,J,K}$

to complete the synchronization so that the resulting scheme is conservative.
However, if $\Lambda_{I,J,K}<1$ adding the entire increment $\delta
R^{\ell}_{I,J,K}$ to $U^{\ell}_{I,J,K}$ could introduce an instability.
Instead the update needs to be redistributed as well. This redistribution can
be done using either the methodology used in FRD or using a WSRD approach.
Here we describe the FRD approach for this “redistribution of reflux” step.
Specifically, we update

$U^{\ell}_{I,J,K}:=U^{\ell}_{I,J,K}+\frac{\Lambda_{I,J,K}}{V_{I,J,K}}\,\delta\mathrm{R}^{\ell}_{I,J,K}$

and redistribute the remainder,
$(1-\Lambda_{I,J,K})\,\delta\mathbf{R}^{\ell}_{I,J,K}$ to neighboring cells
(neighboring in the FRD sense, see Appendix A; not a WSRD neighborhood). In
particular we define

$V_{nbh}=\sum_{II,JJ,KK\,\in\,\mathbf{N}({I,J,K})}V_{II,JJ,KK}$

where $\mathbf{N}(I,J,K)$ is the set of regular and cut cells $(II,JJ,KK)$
such that none of the indices differs from $(II,JJ,KK)$ by more than one and
where $(II,JJ,KK)$ can be reached from $(I,J,K)$ by a monotone path of regular
and cut cells in index space; i.e., a path in which none of the indices are
both incremented and decremented.

We then update cells $(II,JJ,KK)$ in $\mathbf{N}(I,J,K)$ by

$U^{\ell}_{II,JJ,KK}:=U^{\ell}_{II,JJ,KK}+\frac{(1-\Lambda_{I,J,K})\,\delta\mathbf{R}^{\ell}_{I,J,K}}{V_{nbh}}$

This redistribution process can update covered coarse cells. In that case the
correction to the covered coarse cell must be interpolated to the fine cells
that cover the coarse grid cell. In particular, if $(II,JJ,KK)$ is covered, we
use piecewise constant interpolation to update

$U^{m}_{i,j,k}:=U^{m}_{ii,jj,kk}+\frac{(1-\Lambda_{i,j,k})\,\delta\mathbf{R}^{\ell}_{i,j,k}}{V_{nbh}}\qquad\forall\,{i,j,k}\,|\,C^{m}_{i,j,k}\subset
C^{\ell}_{II,JJ,KK}\;\mathrm{for}\;m>\ell$

at all cells $(i,j,k)$ at level $(\ell+1)$ that overlay cell $(II,JJ,KK)$ at
level $\ell.$

For cut cells ($\Lambda_{i,j,k}<1$), adding the entire flux correction from
Eq. (5) can also introduce an instability. Consequently, reflux increments
also need to be redistributed. Here we have again used the FRD style
redistribution for the reflux increments. We note that, operationally, this
can be done by adding the reflux corrections to $\delta\mathrm{R}$ before
redistributing it.

Method-of-lines temporal integration

The methodology described above can be easily extended to work using a method-
of-lines discretization. Berger and Giuliani [15] show how to use state
redistribution in the context of a simple predictor / corrector temporal
discretization. Here we will derive essentially the same scheme but from a
different perspective. The alternate perspective provides useful insight into
the requirements for the synchronization.

The WSRD process can be viewed as defining an operator, $R^{WSRD}$ such that

$U^{n+1}=R^{WSRD}(U^{n}+\Delta t\;\delta U^{c})$ (17)

In the absence of slopes (Eq. (12))

$R^{WSRD}=A^{T}\mathrm{Diag}(\widehat{V})^{-1}A\mathrm{Diag}(V)$

Starting with Eq. (17) we define an effective $\delta U^{n,WSRD}$ using

$\displaystyle U^{n+1}$ $\displaystyle=R^{WSRD}(U^{n}+\Delta t\;\delta U^{c})$
$\displaystyle=U^{n}+\Delta t\;\frac{R^{WSRD}(U^{n}+\Delta t\>\delta
U^{c})-U^{n}}{\Delta t}$ $\displaystyle=U^{n}+\Delta t\;\delta U^{n,WSRD}$
(18)

We note that, when the reconstructed slopes are not considered, $R^{WSRD}$ is
linear and we can rewrite

$\delta U^{n,WSRD}=R^{WSRD}(\delta U^{c})+\frac{1}{\Delta t}(R^{WSRD}-I)U^{n}$

The properties of $A$ show that the second term on the right hand side is an
(anisotropic) numerical diffusion term with magnitude that scales with
$\frac{1}{\Delta t}$.

We can then define the predictor / corrector scheme as

$\displaystyle U^{*}$ $\displaystyle=U^{n}+\Delta t\;\delta U^{n,WSRD}$
$\displaystyle U^{n+1}$ $\displaystyle=\frac{U^{n}+U^{*}+\Delta t\;\delta
U^{*,WSRD}}{2}$ $\displaystyle=U^{n}+\Delta t\;\frac{\delta U^{n,WSRD}+\delta
U^{*,WSRD}}{2}$ (19)

This form shows that in forming $\delta R_{I,J,K}^{\ell}$ and $\delta
R_{i,j,k}^{\ell+1}$ using Eqs. (14) and (15), respectively, the contributions
need to be multiplied by $\frac{1}{2}$ to reflect their weightings in Eq.
(3.4.2).

This basic approach can be extended to higher-order temporal integration
strategies with suitable weights in defining the $\delta R$ that reflect the
underlying temporal quadrature rule. (We note that for higher-order temporal
discretizations in conjunction with AMR and subcyling, special care must be
taken in defining boundary conditions for fine grid from coarse data to avoid
order reduction. See, for example, [25]).

##  4 Numerical results

In this section, we present numerical results demonstrating the properties of
the proposed state re-redistribution method. We use two shock tube problems,
with different initial conditions and grid refinement strategies to
demonstrate the mass conservation properties (Section 4.1) and to provide a
validation case (Section 4.2). A shock interacting with a cylinder provides
additional validation in the presence of multiple AMR levels (Section 4.3).
These cases are run with the CAMR code [26]. Finally, we demonstrate the use
of state re-redistribution with a more realistic three-dimensional geometry
and present some performance metrics showing the cost savings from only
partially refining the embedded boundary (Section 4.4). This final problem is
run with the compressible reacting flow solver, PeleC [27, 28, 29], part of
the Pele suite of solvers.

###  4.1 Conservation properties

We first demonstrate that finite volume embedded boundary methods without the
re-redistribution step violate conservation when the coarse-fine interface
intersects the embedded boundary, and that the re-redistribution step restores
conservation. For this demonstration we choose very simple cases. In the
three-dimensional case, the fluid region is the interior of a cylinder whose
axis passes through $(0,0,0)$ in a [-2:2,-2:2,-2:2] domain and is 30∘ off the
x-axis in the y direction. In the two-dimensional case we take the
intersection of that cylinder with the $z=0$ plane, forming a fluid region
bounded above and below by lines slanted at 30∘ off the x-axis in the
y-direction. The radius of the cylinder is $0.172$. The initial data in both
cases is zero velocity with a discontinuous one-dimensional profile for
density and pressure, namely $\rho=0.125$ and $p=0.1$ for $x\leq 0$ and
$\rho=1.0$ and $p=1.0$ for $x>0.$ The region of the domain define by $-2\leq
x\leq 2$, $-0.125\leq y\leq 0.125$ and, in 3D, $-2\leq z\leq 2$ is refined by
a factor of two relative to the base level. Figure 2 illustrates the solution
at times $t=0$ and $t=0.4$, the bounding geometry and region of grid
refinement. Here we have intentionally presented results from a low resolution
simulation to emphasize the role of the cut cells. When adequately resolved,
at early times the initial discontinuity evolves into a rarefaction wave
expanding to the right, a contact discontinuity and a shock wave propagating
to the left. The shock reflects off the upper boundary and generates an
expansion wave at the lower boundary. Later in the evolution, the different
waves begin to interact with each other, producing a fairly complex flow even
in this relatively simple geometry.

Figure 3 shows the change in mass per time step in the entire domain. It is
hard to scale this plot in a meaningful way, so here we compare the rate of
mass loss that would occur without refluxing on the same plot as the rate of
mass loss due to not re-redistributing. In Figure 3(a) we show two-dimensional
results using the Godunov methodology and state redistribution. In Figure 3(b)
we show two-dimensional results using the method-of-lines methodology and
state redistribution in order to demonstrate the lack of sensitivity of the
results to the specific advection scheme, and that state re-redistribution
behaves correctly for a multi-step integrator. In Figure 3(c) we show two-
dimensional results using the Godunov methodology and flux redistribution to
compare the behavior of flux re-redistribution and state re-redistribution.
And finally in Figure 3(d) we show the three-dimensional results using the
Godunov methodology and state redistribution.

We note that the flux of mass per time step across the coarse/fine interface
over the time shown here ranges from roughly $10^{-7}$ to $10^{-3}$ in two
dimensions and $10^{-8}$ to $4.\times 10^{-4}$ in three dimensions. Thus we
can see that the synchronization accounts for at most a few percent of the
total mass crossing the coarse/fine interface. It should also not be
surprising that the effect of refluxing is greater than that of re-
redistribution; the refluxing corrections occur over the entire
coarse/interface while the re-redistribution only affects the solution at the
part of the coarse/fine interface that is near the embedded boundary. We also
note, by contrasting the solutions with and without the synchronization steps,
that the difference in solutions with and without the synchronization steps is
also on the order of a few percent in max norm.

We mention a few caveats here. This problem is relatively simple and run for a
relatively short time. For more complicated problems run for longer, it is
quite possible for the methodology to fail in the absence of either refluxing
or the redistribution synchronization. In addition, recall that the final step
of the synchronization methodology is the redistribution of the correction
after it is added to the coarse solution. It is straightforward to demonstrate
that if the cylinder is shifted so that there are uncovered coarse cells with
sufficiently small volume fractions (e.g., less than $10^{-4}$ in this case)
at the coarse/fine boundary, the methodology fails in the first time step with
negative energy and pressure, due to too large an update being added to that
cell in the synchronization step.

Finally, we note that while we have shown the behavior of the algorithm in
terms of mass conservation, the re-redistribution methodology acts on all
conserved variables, thus maintaining the appropriate conservation properties
of the entire solution.

(a) $t=0.0$.

(b) $t=0.4$.

Figure 2: Density profiles for the two-dimensional simulation of flow in a
rotated channel with fixed mesh refinement.

(a) Two-dimensional simulation with Godunov methodology and state
redistribution.

(b) Two-dimensional simulation with method-of-lines methodology and state
redistribution.

(c) Two-dimensional simulation with Godunov methodology and flux
redistribution.

(d) Three-dimensional simulation with Godunov methodology and state
redistribution.

Figure 3: Change in mass per timestep as function of time

###  4.2 Sod shock tube validation

Using the same two-dimensional geometry described in Section 4.1, we perform a
validation test case using the Sod shock tube problem. For this problem the
initial discontinuity is normal to the centerline of the channel, i.e., the
left and right states are defined with respect to $x^{\prime}>0$ where
$x^{\prime}=x\cos(\theta)+y\sin(\theta)$ and $\theta$ is the angle of the
channel. This case is similar to validation case considered by Gulizzi et al.
[6]. The AMR refinement is based on a dynamic tagging criteria: additional
levels are added where the local change in density between adjacent cells is
greater than 0.05. Simulations were performed with one level of AMR and three
base grid resolutions: $\Delta x_{0}=0.0125,0.00625,$ and $0.003125$. Coarse-
fine interfaces intersect the EB in this problem, and these intersections are
dynamic, i.e., they move depending on the flow features. These simulations
used the method-of-lines methodology and the proposed state re-redistribution.

The simulation results are presented in Figure 4. The data is taken through a
line in the center of the cylinder. There is good agreement across all flow
variables and mesh resolutions with the exact solution. As the base grid is
refined, the numerical solution converges to the exact solution. As expected
from previous results using the state re-redistribution, mass is conserved
throughout the simulation.

(a) Density

(b) Pressure

(c) Velocity

(d) Temperature

Figure 4: Profiles for the Sod shock tube problem in a rotated channel at
$t=0.1$.

###  4.3 Flow past a cylinder

Here we consider two-dimensional flow past a cylinder and contrast the results
with cases in which the cylinder is fully refined and with experimental data.
This case has been extensively used in other works to validate numerical
approaches [6, 30, 31, 32, 33, 34, 35] by comparing to experimental results
[36]. For this case, the domain is $[0,1]\times[0,1]$ with the center of the
cylinder located at $(0.5,0.5)$. A planar shock wave is initialized at $x=0.2$
and travels towards the cylinder at a Mach number $M=2.81$. The initial right
state, upstream of the shock, is $\rho=1.0$, $u=v=0$, and
$p=\nicefrac{{1}}{{\gamma-1}}$, where $\gamma=1.4$. The left state, downstream
of the shock, is given by the shock jump relations [37]. Outflow conditions
are specified at the $x$-boundaries and walls are specified at the
$y$-boundaries. The base grid for the simulation is $256\times 256$ cells. Two
levels of static refinement are used. We intentionally refine only the top
half of the cylinder at the first level, and slightly less at the finest,
i.e., second, level. The locations of the refinement regions are shown in
Figure 5(a). This problem differs from the earlier cases in that there are two
levels of refinement, and the flow is around an object rather than within a
constrained geometry. Comparisons are made to a case where the refinement
region encapsulates the cylinder entirely, instead of bisecting it, as shown
in Figure 5(d).

Simulation results are shown in Figures 5 and 6. For the case where the top
half of the cylinder if refined, the solution in the refined region is similar
to the case with full refinement around the cylinder. The simulation results
also retain all the features, e.g., the reflected shock (R.S.), the primary
Mach stem (M.S.1), the primary contact discontinuity (C.D.1), and the primary
triple point (T.P.1), observed in the experimental results, Figure 7.

(a) $t=0.0$, showing full domain and refinement levels.

(b) $t=0.068$, close-up view.

(c) $t=0.088$, close-up view.

(d) $t=0.0$, showing full domain and refinement levels.

(e) $t=0.068$, close-up view.

(f) $t=0.088$, close-up view.

Figure 5: Profiles of density for the case of a shock interacting with a
cylinder. Top row: partial refinement around the cylinder; bottom row: full
refinement around the cylinder.

(a) $t=0.068$.

(b) $t=0.088$.

Figure 6: Profiles of density gradient magnitude for the case of a shock
interacting with a cylinder. Top row: partial refinement around the cylinder;
bottom row: full refinement around the cylinder. Figure 7: Comparison between
the experimental Schlieren photograph [36] and the numerical solution at
$t=0.088$ (case where only half the domain is refined, the section of the
domain where the refinement is performed is shown here on the bottom of this
figure). Note: the figure from the reference is mirrored compared to the case
description.

###  4.4 Flow in a compression ignition engine

For this case, we use the new state re-redistribution scheme for three-
dimensional simulation of a production-level case, a compression ignition
engine. These simulations were performed with PeleC.

####  4.4.1 Case description

The case chosen for this demonstration is a large scale simulation of a
simplified piston-bowl geometry, at realistic engine conditions, with a dual
pulse fuel injection. In contrast with previous cases described in this work,
this case involves solving the compressible Navier-Stokes equations, including
viscous stresses, thermal conduction, and species diffusion [27]. We note
that, similarly to what was shown for state redistribution in [17], the
inclusion of the diffusive fluxes does not change the re-redistribution
algorithm.

This case is ideally suited to demonstrating the advantages of adaptive
refinement crossing the the EB as it contains four different jets interacting
with a geometry, where, for most of the simulation run time, it is not
necessary to have full EB refinement. This case forms the basis of the
performance benchmark used by PeleC on high performance computing systems
[27]. The geometry contains a 24mm diameter cylinder and a piston head,
leading to a 2.29mm top volume height and a 3.92mm cavity extending into the
piston head. The geometry and jet positions are shown in Figure 8. This small
volume, representative of a piston near top-dead-center, forms the
computational domain. For these demonstration simulations, a relatively coarse
Cartesian base grid is chosen: 128 $\times$ 128 $\times$ 32 cells.

Figure 8: Computational domain, EB geometry, and jet locations for the
compression ignition engine case.

The initial conditions in the chamber are a mixture of methane (CH4) and
vitiated air at an equivalence ratio $\phi=0.5$. The initial pressure and
temperature in the domain are 6079500.0 Pa and 900K, respectively. Turbulent
velocity fluctuations are imposed on the initial velocity field to represent a
low-swirl flow.

A diesel surrogate, $n$-dodecane (NC12H26), is directly injected through four
jets, angled at a $45^{\circ}$ angle, at the $xy$-plane at the top of the
domain, 10 mm downstream of the actual diesel nozzle. A 35 transported species
mechanism is used for the finite rate chemistry. The diameter of the jets is
680 $\mu$m. They are assumed to be in a fully gaseous phase (fully evaporated
liquid jets). The mixture composition of the gaseous jets are 45% of
$n$-dodecane and 55% of the initial chamber mixture. They are injected at
$470$K and a normal velocity of 28 m/s. A precursor turbulent pipe flow is
used to superimpose turbulent fluctuations on the jet velocities. The
injection duration is 0.5 ms and the jet injection velocity is ramped down
with a hyperbolic tangent function.

In previous version of this test case, as reported elsewhere [27], we had the
following options that avoided having coarse-fine interfaces cross the EB: (1)
refine the EB entirely to the finest level present in the domain, or (2) force
the EB to be refined only up to a specified level, typically less than the
finest level present in the domain and do not allow for refined sections at
the EB. Option 1 leads to a computationally unfeasible number of cells, even
on the largest supercomputers. For example, if one is resolving the ignition
kernels and therefore requiring six levels of AMR on a base grid of 512
$\times$ 512 $\times$ 128 cells, even without including the required cells on
the EB, there would 68 billion cells in the interior of the domain. Including
the cells on the EB would lead to $\mathcal{O}(10)$ times more cells. Option 2
is the one currently pursued for production simulations of this case but leads
to underresolved physics at the EB, particularly when the jets start
interacting with the geometry. A third option is enabled through this work:
refining the EB using different resolutions to respond to the physical
phenomena at the EB.

####  4.4.2 Results

We assess the three refinement options using this case by comparing flow
features and performance metrics. For the cases presented here, in order to
maintain a reasonable computational cost while still demonstrating the
capabilities of the proposed scheme, the simulations use one level of adaptive
mesh refinement and chemical reactions are not evaluated (though all the
relevant species are transported). Simulations were conducted to $t=0.9$ms.
They were performed on 4 nodes (36 Intel Xeon-Gold Skylake-6154 processor
cores per node) on NREL’s Eagle supercomputer [38].

Flow field comparisons for the three options for the magnitude of velocity and
temperature are shown in Figure 9 and for the mass fractions of methane and
$n$-dodecane in Figure 10. From these flow field visualizations, it is clear
that the option that forces the EB to be refined at a level coarser than the
finest available level (middle row in the figures) shows some notable
differences compared to the the other two cases (full refinement at the EB and
adaptive refinement at the EB). For this option, as the jet interacts with the
wall, starting around $t=0.45$ms, the jet features appear more diffuse. The
jet spreads further along the EB (over the inflection point and into the
bottom of the piston bowl). The case using adaptive refinement at the EB (top
row) presents flow feature identical to the fully refined case (bottom row).
Table 1 shows the performance metrics of these cases at two different times
($t=0.45$ms and $t=0.9$ms): the number of cells on the finest level, the
percentage of the domain that is covered by the finest levels, and the time
per time step. At $t=0.45$, just as the jets reach the walls, the time per
time step for the case with adaptive refinement at the EB and the case with a
coarse EB refinement are similar and the percent of the domain that is covered
by the finest level is approximately the same. The case with full refinement
of the EB is approximately four times more computationally expensive, while
achieving the same flow field results. At $t=0.9$ms, after the jets have
significantly interacted with the walls, the case with adaptive refinement at
the EB has increased in computational cost by $50\%$ but still only has about
half the cells on the finest level compared to the case with full EB
refinement and is a 3.2 times faster computationally.

(a) $t=0.2$ms

(b) $t=0.45$ms

(c) $t=0.9$ms

Figure 9: Pseudocolor vizualization of the velocity magnitude and temperature
fields. Top row: using state re-redistribution, middle row: refining the EB to
one level coarser than the finest level, bottom row: refining the EB to the
maximum number of AMR levels.

(a) $t=0.2$ms

(b) $t=0.45$ms

(c) $t=0.9$ms

Figure 10: Pseudocolor vizualization of the CH4 and NC12H26 ($n$-dodecane) mass fraction fields. Top row: using state re-redistribution, middle row: refining the EB to one level coarser than the finest level, bottom row: refining the EB to the maximum number of AMR levels. Max level at EB |  $\\#$ cells on finest level |  $\%$ finest level | time per step (s)
---|---|---|---
at $t=0.45$ms: |  |  |
1 (adaptive) | 491520 | 11.7 | 10.1
0 | 327680 | 7.8 | 10.4
1 (full) | 1982464 | 47.3 | 44.0
at $t=0.9$ms: |  |  |
1 (adaptive) | 1150976 | 27.4 | 14.4
0 | 319488 | 7.6 | 10.0
1 (full) | 2068480 | 49.0 | 45.3
Table 1: Performance comparison at $t=0.45$ms and $t=0.9$ms. Top row: using
state re-redistribution, middle row: refining the EB to one level coarser than
the finest level, bottom row: refining the EB to the maximum number of AMR
levels.

##  5 Conclusions

Weighted state redistribution is a recently developed approach for addressing
the “small cell” problem that arises in embedded boundary discretizations.
When coupled with adaptive mesh refinement, additional synchronization steps
are needed when the interfaces between different levels of refinement are
close to the embedded boundary to ensure that the overall method is both
stable and conservative. Consequently, previous implementations of (weighted)
state redistribution with AMR have required the entire embedded boundary to be
resolved at the same level, typically the finest level in the problem.

In this paper, we have developed a synchronization procedure for WSRD that
removes this restriction. The resulting methodology was tested in two and
three space dimensions for systems of hyperbolic conservation laws. We
demonstrated that the algorithm is conservative and stable in fairly simple
geometries using both a single-step piecewise-linear Godunov method and a
predictor/corrector method of lines temporal integration scheme. We also
presented validation cases for a Sod shock tube in slanted channel and shock
reflection from a sphere. Finally, we demonstrated the methodology on a more
complex model for a compression ignition engine that includes multiple species
and diffusive transport. In this more complex example we illustrate the
potential computational savings resulting from not requiring the entire
embedded boundary to be resolved at the finest level.

The methodology here can be extended to more general compressible flows. In
particular, following the overall approach outlined in [17], the approach
could be extended to include reactions. This approach can also be extended to
incompressible flows and more general low Mach number models. We also note
that there are a number of potential variations to the basic WSRD approach
that could be considered.

##  6 Funding

The work of JBB was supported by the U.S. Department of Energy, Office of
Science, Office of Advanced Scientific Computing Research, Applied Mathematics
Program under contract No. DE-AC02-05CH11231. The work of ASA was supported in
part by the U.S. Department of Energy (DOE) Office of Advanced Scientific
Computing Research (ASCR) via the Scientific Discovery through Advanced
Computing (SciDAC) program FASTMath Institute. The work of IBS and WZ was
supported by the Exascale Computing Project (17-SC-20-SC), a collaborative
effort of the U.S. Department of Energy Office of Science and the National
Nuclear Security Administration. This work was authored in part by the
National Renewable Energy Laboratory, operated by Alliance for Sustainable
Energy, LLC, for the U.S. Department of Energy (DOE) under Contract No. DE-
AC36-08GO28308. Funding was provided by U.S. Department of Energy Office of
Science and National Nuclear Security Administration. This research was
supported by the Exascale Computing Project (17-SC-20-SC), a collaborative
effort of the U.S. Department of Energy Office of Science and the National
Nuclear Security Administration. A portion of the research was performed using
computational resources sponsored by the Department of Energy’s Office of
Energy Efficiency and Renewable Energy and located at the National Renewable
Energy Laboratory. The views expressed in the article do not necessarily
represent the views of the DOE or the U.S. Government. The U.S. Government
retains and the publisher, by accepting the article for publication,
acknowledges that the U.S. Government retains a nonexclusive, paid-up,
irrevocable, worldwide license to publish or reproduce the published form of
this work, or allow others to do so, for U.S. Government purposes.

##  Appendix A Flux redistribution scheme

Flux redistribution is one approach to dealing with the small cell problem.
The basic idea is to use the conservative but unstable update to the solution,
$\delta U^{\mathrm{c}}_{i,j,k}$, from Eq. (2.2) to compute a stable but non-
conservative update, $\delta U^{\mathrm{nc}}_{i,j,k}$. A portion of these two
updates is used to update the cell and the remainder is “redistributed” to
neighboring cells so that the resulting algorithm is conservative and stable.
We first define $\delta U^{\mathrm{c}}_{i,j,k}$ using Eq. (2.2) for each valid
cell in the domain.

Then, as done in Trebotich and Graves [14], for each cell cut by the EB
geometry, we compute the non-conservative update,

$\delta
U^{\mathrm{nc}}_{i,j,k}=\dfrac{\sum\limits_{ii,jj,kk\in\mathbf{N}(i,j,k)}\Lambda_{ii,jj,kk}\;\delta
U^{\mathrm{c}}_{ii,jj,kk}}{\sum\limits_{ii,jj,kk\in\mathbf{N}(i,j,k)}\Lambda_{ii,jj,kk}}$

where $\mathbf{N}(i,j,k)$ is the set of cells $(ii,jj,kk)$ such that none of
the indices differs by more than one from $(i,j,k)$ and where $(ii,jj,kk)$ can
be reached from $(i,j,k)$ by a monotone path of valid cells in index space;
i.e., a path in which none of the indices are both incremented and
decremented. (The approach taken here is not unique. See Pember et al. [7] and
Colella et al. [8] for an alternative approach to computing the
nonconservative updates.)

For each cell cut by the EB geometry, we compute the initial convective update
$\delta U^{EB,*}_{i,j,k}$ using

$\delta U^{EB,*}_{i,j,k}=\Lambda_{i,j,k}\delta
U^{\mathrm{c}}_{i,j,k}+(1-\Lambda_{i,j,k})\delta U^{\mathrm{nc}}_{i,j,k}$ (20)

Eq. (20) provides a stable but non-conservative update for cell-${i,j,k}$. A
direct computation shows that for cut cells, the extensive difference between
$\delta U^{EB}$ and the conservative update is given by

$\displaystyle\delta M_{i,j,k}$
$\displaystyle=V_{i,j,k}(1-\Lambda_{i,j,k})[\delta
U^{\mathrm{c}}_{i,j,k}-\delta U^{\mathrm{nc}}_{i,j,k}]$ (21)
$\displaystyle=\Delta x\Delta y\Delta
z\;\Lambda_{i,j,k}(1-\Lambda_{i,j,k})[\delta U^{\mathrm{c}}_{i,j,k}-\delta
U^{\mathrm{nc}}_{i,j,k}]$

For the method to be conservative, we need to redistribute $\delta M$ to the
neighboring cells. This is done by modifying Eq. (20) using

$\delta U^{EB}_{ii,jj,kk}=\delta
U^{EB,*}_{ii,jj,kk}+\sum_{i,j,k|ii,jj,kk\in\mathbf{N}(i,j,k)}w^{i,j,k}_{ii,jj,kk}\delta
M_{i,j,k}\ $ (22)

where the weights, $w^{i,j,k}_{ii,jj,kk}$, are given

$w^{i,j,k}_{ii,jj,kk}=\dfrac{1}{\sum\limits_{ii,jj,kk\in\mathbf{N}({i,j,k})}V_{ii,jj,kk}}=\frac{1}{\Delta
x\,\Delta y\,\Delta
z}\dfrac{1}{\sum\limits_{ii,jj,kk\in\mathbf{N}({i,j,k})}\Lambda_{ii,jj,kk}}$
(23)

We note that a number of alternative weightings for the redistribution are
possible. For example, in gas dynamics a density weighted redistribution given
by

$w^{i,j,k}_{ii,jj,kk}=\dfrac{\rho_{ii,jj,kk}}{\sum\limits_{ii,jj,kk\in\mathbf{N}({i,j,k})}\rho_{ii,jj,kk}V_{ii,jj,kk}}$

is often used.

##  Appendix B Flux re-redistribution

This appendix assumes that the reader has first read the introduction to
Section 3.4. Here we give the details for how to compute $\delta R$ for flux
redistribution algorithms. The approach described here is taken from Pember et
al. [7] and the reader is referred to that paper for a more detailed
description of the algorithm.

On the coarse grid, for all coarse uncovered cells $({I,J,K})$ at the
coarse/fine boundary, we define

$\delta R^{\ell,u\rightarrow
c}_{I,J,K}=\sum_{II,JJ,KK\,\in\,\mathbf{N}({I,J,K})\,\cap\,\mathbf{C}^{\ell}}V_{II,JJ,KK}\,w^{\ell,{I,J,K}}_{II,JJ,KK}\,\delta
M^{\ell}_{I,J,K}$

which accounts for all contributions of type (a), and

$\delta R^{\ell,c\rightarrow
u}_{I,J,K}=\sum_{(II,JJ,KK)\,|\,{I,J,K}\,\in\,\mathbf{N}(II,JJ,KK)\;\mathrm{and}\;(II,JJ,KK)\in\mathbf{C}^{\ell}}V_{I,J,K}\,w^{\ell,II,JJ,KK}_{I,J,K}\,\delta
M^{\ell}_{II,JJ,KK}$

which accounts for all contributions of type (b). We note here that any given
coarse cell can be have contributions of both type (a) and type (b). Here
$\mathbf{C}^{\ell}$ is the set of covered coarse cells at level $\ell.$

We also need to capture the effect of redistribution on the fine grid. At all
fine ghost cells, we define

$\delta R^{\ell+1,g\rightarrow
v}_{i,j,k}=\sum_{ii,jj,kk\,\in\,\mathbf{N}({i,j,k})\,\cap\,\mathbf{V}^{\ell+1}}V_{ii,jj,kk}w^{\ell+1,{i,j,k}}_{ii,jj,kk}\delta
M^{\ell+1}_{i,j,k}$

which accounts for all contributions of type (c), and

$\delta R^{\ell+1,v\rightarrow
g}_{i,j,k}=\sum_{ii,jj,kk\,|\,{i,j,k}\,\in\,\mathbf{N}(ii,jj,kk)\;\mathrm{and}\;ii,jj,kk\,\in\,\mathbf{V}^{\ell+1}}V_{{i,j,k}}w^{\ell+1,ii,jj,kk}_{{i,j,k}}\delta
M^{\ell+1}_{ii,jj,kk}$

which accounts for all contributions of type (d). Here $\mathbf{V}^{\ell+1}$
is the set of valid cells at level $\ell+1$

We then define

$\delta\mathbf{R}_{I,J,K}^{\ell}\ =\delta R^{\ell,u\rightarrow
c}_{I,J,K}-\delta R^{\ell,c\rightarrow
u}_{I,J,K}+\sum_{i,j,k\,|\,C^{\ell+1}_{i,j,k}\,\subset\,\mathbf{C}^{\ell}_{II,JJ,KK}}(\delta
R^{\ell+1,v\rightarrow g}_{i,j,k}-\delta R^{\ell+1,g\rightarrow v}_{i,j,k})$

where the sum is over all fine grid ghost cells that overlay cell $({I,J,K})$.
The quantity $\delta R^{\ell}_{I,J,K}$ represents the difference between
redistribution on the coarse and fine grid in extensive form.

See Section 3.4.2 for how $\delta R^{\ell}_{I,J,K}$ is used to update the
coarse and fine solutions.

##  Appendix C Piecewise-linear Godunov with embedded boundaries

We use a piecewise linear Godunov method based on the unsplit discretization
of Colella [23] and Saltzman [24]. The method uses a fourth-order
reconstruction of slopes within a cell. The solution is predicted to cell
faces and half-time levels using a characteristic tracing procedure.
Derivatives transverse to the face are treated in conservation form and act as
sources in the characteristic tracing.

In the version we use with cut cells, we view $U$ as being defined at cell
centers (as opposed to centroids) even for cut cells. When we use these values
to reconstruct slopes within each (regular or cut) cell, we use the original
fourth-order slopes if the fourth-order stencil does not cross any faces with
zero area. If the fourth-order computation is disallowed, we reduce to second-
order slopes in that direction if the second-order stencil does not cross any
faces with zero area. If the second-order computation is disallowed the slopes
in that direction are zeroed.

We then compute fluxes only at faces with nonzero area. If, for a given cell,
in trying to compute the flux on a face, the area fractions of all of the
faces on which we need values to compute the transverse derivatives that
contribute to that flux are not nonzero, then we do not include any transverse
terms. For example, if we are predicting from a cell center to an $x$-face,
unless both of the $y$ faces and both of $z$ faces have nonzero area
fractions, we do not include transverse derivatives in the computation of the
flux on the $x$-face.

Once we have predicted values on either side of faces with nonzero area
fractions, we then compute the flux at each face using a two-shock approximate
Riemann solver. We view the flux as being defined at the center of the face,
but when we compute the divergence of fluxes needed in Eq. (2.2) to update the
solution, we use bilinear interpolation to interpolate the fluxes from face
centers to face centroids.

## References

  * [1] N. Gokhale, N. Nikiforakis, R. Klein, A dimensionally split Cartesian cut cell method for hyperbolic conservation laws, J. Comput. Phys. 364 (2018) 186–208.
  * [2] M. J. Berger, C. Helzel, R. J. LeVeque, H-box methods for the approximation of one-dimensional conservation laws on irregular grids, SIAM J. Numer. Anal. 41 (2003) 893–918.
  * [3] B. Muralidharan, S. Menon, A high-order adaptive cartesian cut-cell method for simulation of compressible viscous flow over immersed bodies, Journal of Computational Physics 321 (2016) 342–368.  doi:https://doi.org/10.1016/j.jcp.2016.05.050.
URL https://www.sciencedirect.com/science/article/pii/S0021999116301954

  * [4] R. Saye, Implicit mesh discontinuous galerkin methods and interfacial gauge methods for high-order accurate interface dynamics, with applications to surface tension dynamics, rigid body fluid–structure interaction, and free surface flow: Part i, Journal of Computational Physics 344 (2017) 647–682.
  * [5] R. Saye, Implicit mesh discontinuous galerkin methods and interfacial gauge methods for high-order accurate interface dynamics, with applications to surface tension dynamics, rigid body fluid–structure interaction, and free surface flow: Part ii, Journal of Computational Physics 344 (2017) 683–723.  doi:https://doi.org/10.1016/j.jcp.2017.05.003.
URL https://www.sciencedirect.com/science/article/pii/S0021999117303728

  * [6] V. Gulizzi, A. S. Almgren, J. B. Bell, A coupled discontinuous galerkin-finite volume framework for solving gas dynamics over embedded geometries, Journal of Computational Physics 450 (2022) 110861.
  * [7] R. B. Pember, J. B. Bell, P. Colella, W. Y. Curtchfield, M. L. Welcome, An adaptive cartesian grid method for unsteady compressible flow in irregular regions, Journal of computational Physics 120 (2) (1995) 278–304.
  * [8] P. Colella, D. T. Graves, B. J. Keen, D. Modiano, A cartesian grid embedded boundary method for hyperbolic conservation laws, Journal of Computational Physics 211 (1) (2006) 347–366.
  * [9] X. Y. Hu, B. Khoo, N. A. Adams, F. Huang, A conservative interface method for compressible flows, Journal of Computational Physics 219 (2) (2006) 553–578.
  * [10] R. Klein, K. Bates, N. Nikiforakis, Well-balanced compressible cut-cell simulation of atmospheric flow, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 367 (1907) (2009) 4559–4575.
  * [11] D. Graves, P. Colella, D. Modiano, J. Johnson, B. Sjogreen, X. Gao, A cartesian grid embedded boundary method for the compressible navier–stokes equations, Communications in Applied Mathematics and Computational Science 8 (1) (2013) 99–122.
  * [12] L. Schneiders, D. Hartmann, M. Meinke, W. Schröder, An accurate moving boundary formulation in cut-cell methods, Journal of Computational Physics 235 (2013) 786–809.  doi:https://doi.org/10.1016/j.jcp.2012.09.038.
URL https://www.sciencedirect.com/science/article/pii/S0021999112005839

  * [13] A. S. Almgren, J. B. Bell, P. Colella, T. Marthaler, A cartesian grid projection method for the incompressible euler equations in complex geometries, SIAM Journal on Scientific Computing 18 (5) (1997) 1289–1309.
  * [14] D. Trebotich, D. Graves, An adpative finite volume method for the incompressible Navier-Stokes equations in complex geometries, Communications in Applied Mathematics and Computational Science (2015) 43–82.
  * [15] M. Berger, A. Giuliani, A state redistribution algorithm for finite volume schemes on cut cell meshes, Journal of Computational Physics 428 (2021) 109820\.
  * [16] A. Giuliani, A two-dimensional stabilized discontinuous galerkin method on curvilinear embedded boundary grids, SIAM Journal on Scientific Computing 44 (1) (2022) A389–A415.
  * [17] A. Giuliani, A. S. Almgren, J. B. Bell, M. J. Berger, M. H. de Frahan, D. Rangarajan, A weighted state redistribution algorithm for embedded boundary grids, Journal of Computational Physics (2022) 111305.
  * [18] M. Berger, A. Giuliani, A new provably stable weighted state redistribution algorithm (2023).  arXiv:2308.16332.
  * [19] M. J. Berger, P. Colella, Local adaptive mesh refinement for shock hydrodynamics, Journal of computational Physics 82 (1) (1989) 64–84.
  * [20] J. Bell, M. Berger, J. Saltzman, M. Welcome, Three-dimensional adaptive mesh refinement for hyperbolic conservation laws, SIAM Journal on Scientific Computing 15 (1) (1994) 127–138.
  * [21] W. Zhang, A. Myers, K. Gott, A. Almgren, J. Bell, AMReX: Block-structured adaptive mesh refinement for multiphysics applications, The International Journal of High Performance Computing Applications 35 (6) (2021) 508–526.  doi:10.1177/10943420211022811.
  * [22] J. Bell, M. Day, Adaptive methods for simulation of turbulent combustion, in: Turbulent Combustion Modeling: Advances, New Trends and perspectives, Springer, 2011, pp. 301–329.
  * [23] P. Colella, Multidimensional upwind methods for hyperbolic conservation laws, Journal of Computational Physics 87 (1) (1990) 171–200.
  * [24] J. Saltzman, An unsplit 3d upwind method for hyperbolic conservation laws, Journal of Computational Physics 115 (1) (1994) 153–168.  doi:https://doi.org/10.1006/jcph.1994.1184.
URL https://www.sciencedirect.com/science/article/pii/S0021999184711843

  * [25] M. Emmett, E. Motheau, W. Zhang, M. Minion, J. B. Bell, A fourth-order adaptive mesh refinement algorithm for the multicomponent, reacting compressible navier–stokes equations, Combustion Theory and Modelling 23 (4) (2019) 592–625.
  * [26] CAMR code.
URL https://github.com/AMReX-Codes/CAMR

  * [27] M. T. Henry de Frahan, J. S. Rood, M. S. Day, H. Sitaraman, S. Yellapantula, B. A. Perry, R. W. Grout, A. Almgren, W. Zhang, J. B. Bell, J. H. Chen, PeleC: An adaptive mesh refinement solver for compressible reacting flows, The International Journal of High Performance Computing Applications 37 (2) (2022) 115–131.  doi:10.1177/10943420221121151.
URL https://doi.org/10.1177/10943420221121151

  * [28] H. Sitaraman, N. Brunhart-Lupo, M. H. de Frahan, S. Yellapantula, B. Perry, J. Rood, R. Grout, M. Day, R. Binyahib, K. Gruchalla, Visualizations of direct fuel injection effects in a supersonic cavity flameholder, Phys. Rev. Fluids 6 (2021) 110504.  doi:10.1103/PhysRevFluids.6.110504.
URL https://link.aps.org/doi/10.1103/PhysRevFluids.6.110504

  * [29] H. Sitaraman, S. Yellapantula, M. T. Henry de Frahan, B. Perry, J. Rood, R. Grout, M. Day, Adaptive mesh based combustion simulations of direct fuel injection effects in a supersonic cavity flame-holder, Combustion and Flame 232 (2021) 111531.  doi:https://doi.org/10.1016/j.combustflame.2021.111531.
URL https://www.sciencedirect.com/science/article/pii/S0010218021002741

  * [30] R. Ni, J. Li, X. Zhang, X. Zhou, X. Cui, An immersed boundary-material point method for shock-structure interaction and dynamic fracture, Journal of Computational Physics 470 (2022) 111558.
  * [31] C. Jiang, J. Pan, Y. Zhu, J. Li, E. K. Quaye, Development and verification of a high-speed compressible reactive flow solver in openfoam, Journal of Computational Science 63 (2022) 101780.
  * [32] J. Yang, Y. Liu, H. Lomax, Computation of shock wave reflection by circular cylinders, AIAA journal 25 (5) (1987) 683–689.
  * [33] A. Chaudhuri, A. Hadjadj, A. Chinnayya, On the use of immersed boundary methods for shock/obstacle interactions, Journal of Computational Physics 230 (5) (2011) 1731–1748.
  * [34] J. Zółtak, D. Drikakis, Hybrid upwind methods for the simulation of unsteady shock-wave diffraction over a cylinder, Computer Methods in Applied Mechanics and Engineering 162 (1-4) (1998) 165–185.
  * [35] G. Whitham, A new approach to problems of shock dynamics part i two-dimensional problems, Journal of Fluid Mechanics 2 (2) (1957) 145–171.
  * [36] A. Bryson, R. Gross, Diffraction of strong shocks by cones, cylinders, and spheres, Journal of Fluid Mechanics 10 (1) (1961) 1–16.
  * [37] E. F. Toro, Riemann solvers and numerical methods for fluid dynamics: a practical introduction, Springer Science & Business Media, 2013.
  * [38] National Renewable Energy Laboratory, Eagle, National Renewable Energy Laboratory (2022).
URL https://www.nrel.gov/hpc/eagle-system.html