# The Cartesian Grid Active Flux Method with Adaptive Mesh Refinement 111This
work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) - 325695158. Furthermore, research was supported by a subcontract
with Boise State University, Boise, ID, under DARPA Cooperative Agreement
HR00112120003 with Embry-Riddle Aeronautical University, Daytona Beach, FL,
USA. This work is approved for public release; distribution is unlimited. The
content of the information does not necessarily reflect the position or the
policy of the Government, and no official endorsement should be inferred.

Donna Calhoun Department of Mathematics, Boise State University, Boise, ID,
USA<EMAIL_ADDRESS>Erik Chudzik Institute of Mathematics,
Heinrich-Heine-University Düsseldorf<EMAIL_ADDRESS>Christiane
Helzel33footnotemark: 3<EMAIL_ADDRESS>(corresponding author)

###### Abstract

We present the first implementation of the Active Flux method on adaptively
refined Cartesian grids.

The Active Flux method is a third order accurate finite volume method for
hyperbolic conservation laws, which is based on the use of point values as
well as cell average values of the conserved quantities. The resulting method
has a compact stencil in space and time and good stability properties.

The method is implemented as a new solver in ForestClaw, a software for
parallel adaptive mesh refinement of patch-based solvers. On each Cartesian
grid patch the single grid Active Flux method can be applied. The exchange of
data between grid patches is organised via ghost cells. The local stencil in
space and time and the availability of the point values that are used for the
reconstruction, leads to an efficient implementation. The resulting method is
third order accurate, conservative and allows the use of subcycling in time.

Keywords: Cartesian Grid Active Flux Method, Hyperbolic Conservation Laws,
Adaptive Mesh Refinement

AMS subject classification: 65M08, 65M25, 65M50

##  1 Introduction

The Active Flux method is a finite volume method for hyperbolic conservation
laws that was previously introduced by Eymann, Roe and coauthors [11, 12, 13,
18, 20]. In its original form the method is third order accurate. This is
achieved by using a continuous, piecewise quadratic reconstruction and a
sufficiently accurate quadrature rule to compute numerical fluxes. The
quadrature method, i.e. typically Simpson’s rule, requires point values of the
conserved quantities at grid cell interfaces at the current time as well as at
later time levels. These point values together with the cell average value are
also used to compute the reconstruction. While classical finite volume methods
only use cell average values of the conserved quantities as degrees of
freedom, the Active Flux method involves both point values and cell average
values degrees of freedom. This adds flexibility to the numerical method. For
linear advection and the acoustic equations in one, two or three spatial
dimensions, these point values can be updated using exact evolution formulas,
making the method truly multidimensional.

Originally, Roe and Eymann [13] used unstructured triangular grids for their
two-dimensional active flux method. In [2] and [14], two-dimensional Cartesian
grid versions of the active flux method were introduced. The Cartesian grid
method will also be used in this paper and briefly reviewed in the next
section.

While the order of convergence of a numerical scheme is a property that can be
shown in the limit when the mesh width and the time step goes to zero, for
practical computations it is desirable to obtain accurate results on
relatively coarse grids. Roe [19] argues that the accurate approximation on
coarse grids is strongly influenced by the computational stencil and that
exact evolution operators perform well in this respect. Barsukow showed that
the two-dimensional Cartesian grid Active Flux method for the acoustic
equations with an exact evolution operator as described in [1] is stationary
preserving. This means that the numerical scheme does not add dissipation to
discrete representations of all stationary states of the acoustic equations.
As a consequence such states can be computed with very high accuracy on coarse
grids, while other methods would require a much higher resolution.

In earlier related work, Lukáčová-Medvid’ová et al. [16, 17], use exact
multidimensional evolution operators as building blocks of finite volume
methods of various order. Those methods differ from the Active Flux method
mainly in the choice of the degrees of freedom.

Although the Active Flux method performs well on coarse grids, for practical
applications it might still be desirable to vary the size of the grid cells
adaptively in order to allow a higher resolution in parts of the computational
domain as needed. A possible application that could benefit from local
refinement is the propagation of a high frequency acoustic wave.

In this paper we show how the Active Flux method can be applied on adaptively
refined Cartesian grids. The local stencil of the method allows an efficient
transfer of data between the different grid patches. Our numerical results
confirm third order accuracy of the resulting method.The method is implemented
as a new solver in ForestClaw [6, 8, 9].

This paper is organised as follows. In Section 2 we briefly review the Active
Flux method for two-dimensional Cartesian grids. Section 3 describes the
extension to adaptively refined grids. In Section 4 we introduce new Active
Flux methods for advective transport and illustrate the performance of the
adaptively refined Active Flux method for a variety of test problems.

##  2 The Cartesian grid Active Flux Method

In this section we provide a brief review of the Active Flux method on a
single two-dimensional Cartesian grid. More details can be found in Barsukow
et al. [2] and Helzel et al. [14].

We consider hyperbolic conservation laws in divergence form

$\partial_{t}q+\partial_{x}f(q)+\partial_{y}g(q)=0,$

where $q:\mathbb{R}^{2}\times\mathbb{R}^{+}\rightarrow\mathbb{R}^{s}$ is a
vector of conserved quantities and
$f,g:\mathbb{R}^{s}\rightarrow\mathbb{R}^{s}$ are vector valued flux
functions. On a single patch we use a two-dimensional Cartesian grid with
equidistant mesh sizes $\Delta x$ and $\Delta y$. The grid cell $(i,j)$ is
described by
$[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}]\subset\mathbb{R}^{2},\,i,j,\in\mathbb{Z}.$
As a finite volume method, the Active Flux method computes cell averaged
values of the conserved quantities via an update of the form

$Q_{i,j}^{n+1}=Q_{i,j}^{n}-\frac{\Delta t}{\Delta
x}\left(F_{i+\frac{1}{2},j}-F_{i-\frac{1}{2},j}\right)-\frac{\Delta t}{\Delta
y}\left(G_{i,j+\frac{1}{2}}-G_{i,j-\frac{1}{2}}\right),$ (1)

where $Q_{i,j}^{n}$ is an approximation of the cell average values of the
conserved quantities in grid cell $(i,j)$ at time $t_{n}$ and
$F_{i\pm\frac{1}{2},j}$, $G_{i,j\pm\frac{1}{2}}$ are numerical fluxes at the
grid cell interfaces given by

$\begin{split}F_{i+\frac{1}{2},j}&\approx\frac{1}{\Delta t\Delta
y}\int_{t_{n}}^{t_{n+1}}\int_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}f(q(x_{i+\frac{1}{2}},y,t))dy\,dt\\\
G_{i,j+\frac{1}{2}}&\approx\frac{1}{\Delta t\Delta
x}\int_{t_{n}}^{t_{n+1}}\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}g(q(x,y_{j+\frac{1}{2}},t))dx\,dt.\end{split}$

As suggested by Eymann and Roe, we use Simpson’s rule to compute the numerical
fluxes. For fluxes $F_{i+\frac{1}{2},j}$, this leads to the formula

$\begin{split}F_{i+\frac{1}{2},j}&:=\frac{1}{36}\Big{(}f(Q_{i+\frac{1}{2},j-\frac{1}{2}}^{n})+4f(Q_{i+\frac{1}{2},j}^{n})+f(Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n})\\\
&\quad+4f(Q_{i+\frac{1}{2},j-\frac{1}{2}}^{n+\frac{1}{2}})+16f(Q_{i+\frac{1}{2},j}^{n+\frac{1}{2}})+4f(Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}})\\\
&\quad+f(Q_{i+\frac{1}{2},j-\frac{1}{2}}^{n+1})+4f(Q_{i+\frac{1}{2},j}^{n+1})+f(Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n+1})\Big{)}\end{split}$
(2)

We use an analogous formula for flux $G_{i,j+\frac{1}{2}}$. The $Q$ values in
the right hand side of (2) are approximations to point values of the conserved
quantities at the grid cell interface at times $t_{n}$, $t_{n+\frac{1}{2}}$
and $t_{n+1}$. To compute these point values, we assume that at time $t_{n}$
the average values $Q_{i,j}^{n}$ and point values at cell corners
$Q_{i\pm\frac{1}{2},j\pm\frac{1}{2}}^{n}$ and edge midpoints
$Q_{i\pm\frac{1}{2},j}^{n}$ and $Q_{i,j\pm\frac{1}{2}}^{n}$ are known. The
location of these point values are shown in Figure 1.

Figure 1:  Degrees of freedom used to reconstruct the piecewise quadratic
function in grid cell $(i,j)$. Point values are marked by a dot and the cell
average is marked by a square.

Based on these known values, a two-dimensional quadratic polynomial can be
reconstructed in each grid cell. For simplicity, all Cartesian grid cells are
mapped to a reference cell $[-1,1]\times[-1,1]$. In each cell, the
reconstruction has the form

$q^{rec}(\xi,\eta)=c_{00}+c_{10}\xi+c_{01}\eta+c_{20}\xi^{2}+c_{11}\xi\eta+c_{02}\eta^{2}+c_{21}\xi^{2}\eta+c_{12}\xi\eta^{2}+c_{22}\xi^{2}\eta^{2}$
(3)

with $c_{ij}\in\mathbb{R}$, $i,j=0,1,2$ and $(\xi,\eta)\in[-1,1]\times[-1,1]$.
The precise form of the reconstructed function is obtained by interpolating
the eight known point values along the grid cell boundary and by requiring
that the average of the reconstructed function agrees in each grid cell with
the known cell average. This can conveniently be obtained by using appropriate
basis functions as introduced in [14].

We denote the (at time $t_{n}$) reconstructed continuous, piecewise quadratic
function by $q^{n}(x,y)$. The quadratic polynomial that was reconstructed in
grid cell $(i,j)$ is denoted by $q_{ij}^{n}(x,y)$.

For special linear partial differential equations, the required point values
of the conserved quantities can be calculated using an exact evolution formula
applied to the piecewise quadratic data. This is in particular the case for
linear advection and acoustics. The use of the exact evolution formula leads
to a truly multidimensional numerical method. While the point values at time
$t_{n+\frac{1}{2}}$ are only used in order to compute the numerical fluxes,
the point values at time $t_{n+1}$ are also used to compute the continuous,
piecewise quadratic reconstruction $q^{n+1}(x,y)$, which provides the initial
data for the next time step.

Since the point values are located along the grid cell boundaries, they are
used for the reconstruction in two or four grid cells. Thus, the Active Flux
method for two-dimensional Cartesian grids has four degrees of freedom per
grid cell: the cell average plus three point values along the grid cell
boundary. In our implementation, we associate point values
$Q_{i-\frac{1}{2},j}$, $Q_{i-\frac{1}{2},j-\frac{1}{2}}$,
$Q_{i,j-\frac{1}{2}}$ and the average value $Q_{i,j}$ with grid cell $(i,j)$.

###  2.1 Advection equation

For the advection equation

$\partial_{t}q+a\partial_{x}q+b\partial_{y}q=0,$ (4)

with $q:\mathbb{R}^{2}\times\mathbb{R}^{+}\rightarrow\mathbb{R}$,
$a,b\in\mathbb{R}$, the exact evolution formula has the form

$q(x,y,t+\tau)=q(x-a\tau,y-b\tau,t).$

Thus, the required point values at time $t_{n+\frac{1}{2}}$ and $t_{n+1}$ can
easily be computed by evaluating the continuous, piecewise quadratic function
$q^{n}$ at the appropriate upwind points, for example

$\begin{split}Q_{i+\frac{1}{2},j-\frac{1}{2}}^{n+\frac{1}{2}}&=q^{n}\left(x_{i+\frac{1}{2}}-a\frac{\Delta
t}{2},y_{j-\frac{1}{2}}-b\frac{\Delta t}{2}\right)\\\
Q_{i+\frac{1}{2},j-\frac{1}{2}}^{n+1}&=q^{n}\left(x_{i+\frac{1}{2}}-a\Delta
t,y_{j-\frac{1}{2}}-b\Delta t\right).\end{split}$ (5)

In the numerical method, we want to restrict the time steps, so that waves
propagate at most one grid cell per time step. This can be expressed in the
form

$\mbox{CFL}:=\mbox{max}\left(\frac{|a|\Delta t}{\Delta x},\frac{|b|\Delta
t}{\Delta y}\right)\leq 1.$ (6)

We have shown in [10] that linear stability requires a more severe time step
restriction. In practical computations we therefore use time steps which
satisfy $\mbox{CFL}\leq 0.75$.

In Section 4.1 the approach is extended to advective transport in a spatially
and temporally varying velocity field.

###  2.2 Burgers’ equation

For the two-dimensional Burgers’ equation

$q_{t}+\left(\frac{1}{2}q^{2}\right)_{x}+\left(\frac{1}{2}q^{2}\right)_{y}=0$
(7)

with $q:\mathbb{R}^{2}\times\mathbb{R}^{+}\rightarrow\mathbb{R}$, we do not
have an exact evolution formula. Instead we approximate the point values using
an approximative approach described in [10].

For smooth solutions, equation (7) can equivalently be written in the
advective form

$q_{t}+qq_{x}+qq_{y}=0,$ (8)

which suggests the implicitly defined evolution formula

$q(x,y,t)=q(x-q(x,y,t)t,y-q(x,y,t)t,0).$ (9)

Starting with an initial guess

$\left(Q_{i-\frac{1}{2},j}^{n+\frac{1}{2}}\right)^{0},\quad\left(Q_{i-\frac{1}{2},j}^{n+1}\right)^{0},$

we iteratively compute

$\begin{split}\left(Q_{i-\frac{1}{2},j}^{n+\frac{1}{2}}\right)^{\ell}&=q^{n}\left(x_{i-\frac{1}{2}}-\left(Q_{i-\frac{1}{2},j}^{n+\frac{1}{2}}\right)^{\ell-1}\frac{\Delta
t}{2},y_{j}-\left(Q_{i-\frac{1}{2},j}^{n+\frac{1}{2}}\right)^{\ell-1}\frac{\Delta
t}{2}\right)\\\
\left(Q_{i-\frac{1}{2},j}^{n+1}\right)^{\ell}&=q^{n}\left(x_{i-\frac{1}{2}}-\left(Q_{i-\frac{1}{2},j}^{n+1}\right)^{\ell-1}\Delta
t,y_{j}-\left(Q_{i-\frac{1}{2},j}^{n+1}\right)^{\ell-1}\Delta
t\right)\quad\ell=1,2,\ldots.\end{split}$ (10)

Wave speeds at other positions are computed analogously. Each iteration
improves the accuracy by one order. We start the iteration with a first order
accurate approximation. The piecewise quadratic reconstruction limits the
achievable accuracy to third order. Thus, a third order accurate approximation
can be achieved after two iterations.

An obvious initial guess would be to use the point values at the respective
location, i.e.

$\left(Q_{i-\frac{1}{2},j}^{n+\frac{1}{2}}\right)^{0}=\left(Q_{i-\frac{1}{2},j}^{n+1}\right)^{0}=Q_{i-\frac{1}{2},j}^{n}$

and analogously for all other point values along the grid cell boundaries.
These wave speeds are third order accurate in space and first order accurate
in time. However, this choice suffers from an instability, if the
characteristic speed changes sign as explained in detail in [14]. The
instability can be mitigated if data from all adjacent grid cells are used to
compute the initial guess. This leads to a stronger coupling of wave speeds
and cell average values. Here, we compute the initial guess of the wave speed
based on the neighboring cell average values, i.e. we use

$\begin{split}\left(Q_{i-\frac{1}{2},j}^{n+\frac{1}{2}}\right)^{0}&=\left(Q_{i-\frac{1}{2},j}^{n+1}\right)^{0}=\frac{1}{2}\left(Q_{i-1,j}^{n}+Q_{i,j}^{n}\right)\\\
\left(Q_{i-\frac{1}{2},j-\frac{1}{2}}^{n+\frac{1}{2}}\right)^{0}&=\left(Q_{i-\frac{1}{2},j-\frac{1}{2}}^{n+1}\right)^{0}=\frac{1}{4}\left(Q_{i-1,j}^{n}+Q_{i,j}^{n}+Q_{i-1,j-1}^{n}+Q_{i,j-1}^{n}\right)\\\
\left(Q_{i,j-\frac{1}{2}}^{n+\frac{1}{2}}\right)^{0}&=\left(Q_{i,j-\frac{1}{2}}^{n+1}\right)^{0}=\frac{1}{2}\left(Q_{i,j}^{n}+Q_{i,j-1}^{n}\right).\end{split}$

The slight increase of the stencil, as introduced by this initial guess of the
wave speeds, leads to a stable approximation.

###  2.3 Acoustics

The acoustic equations are given by

$\begin{split}\partial_{t}p+c\nabla\cdot\boldsymbol{u}&=0\\\
\partial_{t}\boldsymbol{u}+c\nabla p&=0,\end{split}$ (11)

where
$\boldsymbol{u}:\mathbb{R}^{2}\times\mathbb{R}^{+}\rightarrow\mathbb{R}^{2}$
is the velocity vector,
$p:\mathbb{R}^{2}\times\mathbb{R}^{+}\rightarrow\mathbb{R}$ is the pressure
and $c\in\mathbb{R}^{+}$ is the speed of sound.

The evolution formula for the two-dimensional acoustic equations used in this
paper can be found in [13]. It is based on the observation that (11) can be
rewritten as

$\begin{split}\partial_{tt}p-c^{2}\Delta p&=0\\\
\partial_{tt}\boldsymbol{u}-c^{2}\Delta\boldsymbol{u}&=c^{2}\nabla\times\boldsymbol{w},\end{split}$
(12)

where $\boldsymbol{w}=\nabla\times\boldsymbol{u}$ is the vorticity and
$\triangle$ is the Laplacian operator. In the two-dimensional case, considered
here, $\boldsymbol{w}=(0,0,v_{x}-u_{y})^{T}$. Thus, in a flow with constant
vorticity both pressure and velocity satisfy a wave equation. Furthermore, it
is easy to verify that the vorticity is stationary, i.e. in the two-
dimensional case the relation $\partial_{t}(v_{x}-u_{y})=0$ holds.

Assuming constant vorticity, Eymann and Roe [13] derived the evolution
formulas

$\begin{split}p(t)&=M_{R}\\{p\\}+R\left(\partial_{R}M_{R}\\{p\\}-M_{R}\\{\nabla\cdot\boldsymbol{u}\\}\right)\\\
\boldsymbol{u}(t)&=M_{R}\\{\boldsymbol{u}\\}+R\left(\partial_{R}M_{R}\\{\boldsymbol{u}\\}-M_{R}\\{\nabla
p\\}\right),\end{split}$ (13)

where $R=c\cdot t$ and $M_{R}\\{f\\}$ is the spherical mean. The values of
pressure and velocity at the right hand side of the evolution equation (13)
are given initial values at time $t=0$. For a scalar function
$f:\mathbb{R}^{2}\rightarrow\mathbb{R}$, the spherical mean over a disc with
radius $R$, centred at $(x,y)$ is defined by

$M_{R}\\{f\\}(x,y):=\frac{1}{2\pi
R}\int_{0}^{2\pi}\int_{0}^{R}f(x+s\cos\phi,y+s\sin\phi)\frac{s}{\sqrt{R^{2}-s^{2}}}dsd\phi.$
(14)

In the vector valued case, the formula is applied component wise. The solution
formula can be evaluated exactly, if during each time step the previous values
of pressure and velocity are replaced by the corresponding components of the
reconstructed continuous, piecewise quadratic function $q^{n}$.

Each time step of the explicit Active Flux method is restricted so that the
circle around the edge midpoint over which the integration takes place remains
inside the two adjacent grid cells. This condition will be met if

$\max\left(\frac{c\Delta t}{\Delta x},\frac{c\Delta t}{\Delta
y}\right)\leq\frac{1}{2}.$ (15)

In [10] we showed that this necessary condition is sufficient for linear
stability on a regular Cartesian grid.

##  3 Adaptive Mesh Refinement for the Active Flux Method

We now describe the implementation of the adaptive Active Flux method as a new
solver in ForestClaw [6], a software for parallel adaptive mesh refinement
based on a quadtree approach. In ForestClaw, Cartesian grid patches occupy
quadrants in a quadtree, or multi-block forest of quadtrees. ForestClaw was
developed by Calhoun and Burstedde based on the p4est software [7].

Mesh refinement is realised by a bisection of grid patches so that a quadrant
of resolution level $\ell$ is replaced by four quadrants of resolution level
$\ell+1$. A patch of level zero would correspond to a single Cartesian grid
discretizing a single, square domain. The number of grid cells on a single
quadrant is constant for all levels, resulting in a 2:1 refinement ratio
between resolution levels. Typically $8\times 8$, $16\times 16$ or $32\times
32$ grid cells are used on a single patch, which allows a flexible change of
the resolution. Furthermore, single grid patches can efficiently be handled by
separate processors in a parallel computation. Figure 2 shows a typical
situation.

Figure 2:  Three refinement levels of quadrants of a ForestClaw mesh are
shown. Each quadrant is occupied by a logically Cartesian grid of fixed
resolution (e.g. $8\times 8$ in the figure). Each grid also contains a layer
of ghost cells (shaded region) which are used to facilitate the exchange of
data between adjacent grid patches and between local and remote processors.

A ForestClaw mesh inherits many properties from the underlying p4est mesh [7].
In particular, a ForestClaw mesh is well-balanced so that adjacent quadrants
never differ by more than one level. This is a necessary condition for the
definition of data transfer between different patches, which is organized
using ghost cells and described in more detail below. Another key feature of
the ForestClaw mesh is that it is dynamically adapted to follow solution
features of interest. At each time step, we apply coarsening and refinement
criteria to the solution on each quadrant. If the coarsening criteria is
satisfied by the solution in each quadrant in a family of four quadrants, the
four quadrants will be replaced by a single quadrant, and the solution will be
averaged from the finer grids to the new coarser solution. Otherwise, if the
solution on a quadrant satisfies the refinement criteria, the quadrant will be
subdivided into four quadrants, and the solution will be interpolated from the
coarse parent to the new finer resolution solution on each of the child
quadrants. The details of the interpolation and averaging are provided next.

###  3.1 Spatial transfer of grid cell information

Solution data in a composite ForestClaw mesh needs to be communicated between
adjacent grids sharing quadrant boundaries, and when dynamically coarsening
and refining the mesh. The Active Flux method allows a very efficient transfer
of both pointwise and cell-average information. In the following, a ”grid
patch” is both the p4est quadrant and the solution data in the quadrant. The
following three situations need to be considered:

  1. 1.

A transfer from a fine grid to a coarse grid is needed if four grid patches at
level $\ell+1$ are coarsened to a single patch at level $\ell$. The same
approach is used for the computation of ghost cell values on a patch of level
$\ell$ from a neighbouring patch of level $\ell+1$.

  2. 2.

A transfer from a coarse grid to a fine grid is needed if a patch of level
$\ell$ is marked for refinement and four patches at level $\ell+1$ need to be
reconstructed. The same approach is used for the computation of ghost cell
values for a patch at level $\ell+1$ from a neighbouring patch of level
$\ell$.

  3. 3.

For neighbouring grid patches of the same level the ghost cell information is
simply copied from the neighbouring grid patches.

We will now discuss the first two approaches in more detail.

#### Transfer from fine to coarse grids

The degrees of freedom of a coarse grid cell are computed from the degrees of
freedom of four grid cells on the finer level as illustrated in Figure 3. The
cell average of the coarse cell is the average of the four cell average values
on the fine grid. The point values are copied from the point values at the
vertices of the fine grid cells.

Figure 3: Illustration of coarsening: The cell average on the coarse grid is
the average of the four average values on the fine grid. The point values are
copied from vertices of the fine grid cells.

#### Transfer from coarse to fine grids

We use the degrees of freedom of the coarse grid cell to reconstruct a
quadratic polynomial as described in Section 2. This polynomial can be
evaluated at all the required point values along the edges of the fine grid
cells. The cell average values of the four fine grid cells are computed using
Simpson’s rule. This requires the additional computation of four point values
on the coarse grid cell at the positions of the centers of the fine grid
cells. An illustration is shown in Figure 4.

Figure 4: Illustration of refinement: Point values and cell average values on
the four fine grid cells are computed from the reconstruction of the conserved
quantities at the coarse grid cell.

Since Simpson’s rule is exact for our reconstructed function the sum of the
cell average values on the four fine grid cells agrees with the cell average
on the coarse grid cell exactly. This is a necessary property in order to
maintain conservation.

###  3.2 Subcycling for a local refinement in space and time

In explicit finite volume methods for hyperbolic conservation laws, the time
step is necessarily restricted by a CFL condition, which requires that the
numerical domain of dependence must contain the true domain of dependence of
the partial differential equations [15].

The Active Flux method has a very compact stencil as explained in Section 2.
Thus, for stability it is necessary to restrict the time step in such a way
that information travels at most through one grid cell. Our results from [10]
show that the time step should be restricted by $\mbox{CFL}\leq 0.75$ for two-
dimensional advection problems and by $\mbox{CFL}\leq 0.5$ for the two-
dimensional acoustic equations.

On an adaptively refined mesh, the smallest grid cells would typically dictate
the time step restriction for the whole domain. To increase the efficiency of
the computation, local time stepping (or ”subcycling”) can be used. In a
subcycled computation, several time steps on more refined patches are taken
for one time step on the coarsest grid. Subcycling was included in the
original AMR algorithm by Berger and Oliger [3] and Berger and Colella [4] and
is a standard feature of many AMR codes, including AMRClaw [5], AMReX [21] and
many others. The local stencil of the Active Flux method allows for efficient
implementation of subcycling, which will now be described in more detail.

#### Subcycling

In order to use subcycling, we reconstruct the solution in two layers of ghost
cells surrounding each patch.

We illustrate the idea for the one-dimensional case but an extension to the
two-dimensional situation is straightforward. Assume a situation with three
different grid patches as illustrated in Figure 5.

Figure 5:  Illustration of one-dimensional grid interfaces with different
resolution

The coarsest grid, denoted as grid 1, has refinement level $\ell$. The grid
cells on this part of the domain have length $\Delta x_{\ell}$. Our grid 2 has
refinement level $\ell+1$ and the grid cell length satisfies $\Delta
x_{\ell+1}=\Delta x_{\ell}/2$. On the finest mesh, i.e. grid 3, we have
$\Delta x_{\ell+2}=\Delta x_{\ell+1}/2$. For efficiency, we will use a time
step $\Delta t_{\ell}$ on grid 1, two time steps $\Delta t_{\ell+1}=\Delta
t_{\ell}/2$ on grid 2 and four time steps $\Delta t_{\ell+2}=\Delta
t_{\ell+1}/2=\Delta t_{\ell}/4$ on grid 3.

In the ForestClaw implementation, we start the time stepping on the finest
level grids. Referring to grids 1,2 and 3 described above, the algorithm
proceeds as follows for this three-level mesh configuration.

  1. 1.

Advance the solution one step on grid 3 using time step $\Delta t_{\ell+2}$.

  2. 2.

Recursively advance the solution one step on grid 2 (using time step $\Delta
t_{\ell+1}$) and on grid 1 (using time step $\Delta t_{\ell}$).

  3. 3.

Advance the solution a second step on grid 3.

  4. 4.

Grids 2 and 3 are now time synchronized and ghost cell data is exchanged
between these levels.

  5. 5.

Advance the solution a third step on grid 3.

  6. 6.

Recursively update the solution a second step on level 2

  7. 7.

Advance the solution a fourth step on grid 3

In each grid advance, the first layer of ghost cells is updated along with all
interior cells. These ghost cell values are needed to update interior cells at
intermediate time levels (e.g. fine time levels that do not exist on coarser
levels). However, whenever two levels are time synchronized, the updated ghost
cell data is replaced by data averaged or interpolated from the neighboring
finer or coarser grids at the same time level. For those grids at the physical
boundary, physical boundary conditions are used at all time levels. The time
step on each grid uses a stable time step appropriate for that grid. These
steps are illustrated in Figure 6.

Figure 6: Subcycling algorithm for advancing three levels (with grids 1, 2 and
3) a single coarse grid time step of size $\Delta t_{\ell}$. The algorithm
starts at the lower right by advancing one step on the finest level. All
arrows of the same color are advanced recursively. The blue shaded boxes
indicate time synchronized levels where ghost exchanges between neighboring
grids takes place. Grids without time synchronized coarse grid neighbors
advance by using updated values in an extra layer of ghost cells.

The key difference between the subcycling approach described above and the
subcycling described in the original second order finite volume schemes
described by Berger, Oliger and Colella is that in the Active Flux method,
ghost cells at the intermediate time level are not filled using time
interpolation from neighboring coarser grid cells. Instead we obtain all the
required data directly from previously updated degrees of freedom of the
Active Flux method.

#### Global conservation property

To ensure the global conservation property of the method, the fluxes at grid
cell interfaces need to be defined in a unique way. At regular grid cell
interfaces the Active Flux method automatically satisfies this conservation
property. At interfaces of patches with different refinement level,
conservation needs to be enforced by some kind of “conservative fix”. We use
the classical approach described by Berger and Colella [4], and update both
the coarse and the fine grid cells using the fluxes that have been computed
for the more resolved grid.

##  4 Numerical results

In this section we show numerical results for advection, advective transport
with spatially and temporally varying velocity field, Burgers’ equation and
acoustics. We introduce new Active Flux methods for advective transport
problems and discuss the preservation of constant states.

For all computations we use subcycling and the conservative fix described
above unless otherwise noted. Each grid patch uses $16\times 16$ Cartesian
grid cells plus ghost cells. Our results confirm third order accuracy of the
Active Flux method on adaptively refined grids.

###  4.1 Convergence study for advection

We consider the advection equation (4) on the domain $[0,1]\times[0,1]$ with
initial condition

$q(x,y,0)=H(r(x,y)+r_{0})-H(r(x,y)-r_{0}),$ (16)

where $r(x,y):=\sqrt{(x-x_{0})^{2}+(y-y_{0})^{2}}$ and
$H(r):=(\tanh(r/0.02)+1)/2$ with $r_{0}=0.15$ and $x_{0}=y_{0}=0.5$. We use
the advection speeds $a=1$ and $b=0.5$ and time steps which satisfy
$\mbox{CFL}=0.6$.

Two different adaptively refined grids are considered. In the first case
refinement is allowed only along the diagonal of the domain as shown in Figure
7 (left plot). In the second case a patch is refined, if
$q_{max}-q_{min}>0.001$. In this case the refined grid will follow the
solution structure as shown in Figure 7 (right plot). By comparing the
numerical solution with the exact solution we can measure the error and
compute the experimental order of convergence (EOC). The results are shown in
Table 1 for refinement along the diagonal and for dynamic refinement which
follows the solution structure.

Table 1: Error at time $t=1$ measured in the 1-norm and convergence rates for
the advection problem with constant refinement along the diagonal (left table)
and refinement that resolves the relevant solution structure (right table).

Level | Error | EOC
---|---|---
4 | $6.461\,135\text{\times}{10}^{-05}$ | —
5 | $8.450\,336\text{\times}{10}^{-06}$ | $2.9347$
4-5 | $5.075\,896\text{\times}{10}^{-05}$ | —
5-6 | $6.607\,801\text{\times}{10}^{-06}$ | $2.9414$
Level | Error | EOC
---|---|---
5 | $8.450\,336\text{\times}{10}^{-06}$ | —
6 | $1.066\,256\text{\times}{10}^{-06}$ | $2.9865$
3-5 | $8.467\,527\text{\times}{10}^{-06}$ | —
3-6 | $1.069\,544\text{\times}{10}^{-06}$ | $2.9849$

For the academic test case with refinement along the diagonal, we see that the
changes in the grid structure did not introduce any grid-induced artifacts.
The accuracy observed on the adaptively refined grid is comparable with the
accuracy on a regular Cartesian grid on the coarser level.

If the adaptive mesh follows the solution structure, then the accuracy
obtained on the adaptively refined grid compares well with the accuracy
obtained on a regular Cartesian grid that uses the highest level of refinement
in the full domain. These results are shown in Table 1. The accuracy obtained
on the adaptive mesh with levels $3-5$ or $3-6$ compares well with the
accuracy obtained on the grids that are refined uniformly to levels $5$ or
$6$.

Figure 7: Solutions to (16) with constant refinement along the diagonal (left)
and with refinement level 3-6 (right)at time $t=0.1$.

###  4.2 Solid-body rotation

Next we consider the linear advection equation with variable coefficients

$\partial_{t}q+\partial_{x}(a(x,y)q)+\partial_{y}(b(x,y)q)=0$ (17)

on $[-1,1]\times[-1,1]$. We define the advection field for solid body rotation
using the stream-function $\Psi(x,y)=\frac{\pi}{2}(x^{2}+y^{2})$ and obtain

$\begin{split}a(x,y)&=-\partial_{y}\Psi(x,y)=-\pi y,\\\
b(x,y)&=\phantom{-}\partial_{x}\Psi(x,y)=\phantom{-}\pi x.\end{split}$ (18)

The velocity field is divergence free, thus Equation (17) is equivalent to

$\partial_{t}q+a(x,y)\partial_{x}q+b(x,y)\partial_{y}q=0.$ (19)

Furthermore, the solution at time $t=2n$, $n\in\mathbb{N}$ agrees with the
initial values since the flow has simply made $n$ complete rotations. The
characteristics satisfy the ordinary differential equation

$\displaystyle x^{\prime}(t)$ $\displaystyle=a(x(t),y(t))$ (20) $\displaystyle
y^{\prime}(t)$ $\displaystyle=b(x(t),y(t))$

with initial values $x(0)=x_{0}$, $y(0)=y_{0}$. We are interested in
$(x(-\tau),y(-\tau))$ with $\tau=\Delta t/2,\Delta t$. Those values can be
computed using the analytical solution

$\begin{split}x(t)&=-y_{0}\sin(\pi t)+x_{0}\cos(\pi t)\\\
y(t)&=\phantom{-}y_{0}\cos(\pi t)+x_{0}\sin(\pi t).\end{split}$ (21)

This allows us to evaluate the conserved quantities at all required nodes of
Simpson’s quadrature formula by tracing back the characteristics. For the flux
$F_{i+\frac{1}{2},j}$, we use Equation (2) with, for example,

$\displaystyle f(Q_{i+\frac{1}{2},j}^{n+\frac{1}{2}})$
$\displaystyle=a(x_{i+\frac{1}{2}},y_{j})Q_{i+\frac{1}{2},j}^{n+\frac{1}{2}}=-\pi
y_{j}Q_{i+\frac{1}{2},j}^{n+\frac{1}{2}}$ (22) $\displaystyle
g(Q_{i,j+\frac{1}{2}}^{n+\frac{1}{2}})$
$\displaystyle=b(x_{i},y_{j+\frac{1}{2}})Q_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}=\pi
x_{i}Q_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}$

and

$\displaystyle Q_{i+\frac{1}{2},j}^{n+\frac{1}{2}}$
$\displaystyle=q^{n}\left(-y_{j}\sin\left(-\pi\frac{\Delta
t}{2}\right)+x_{i+\frac{1}{2}}\cos\left(-\pi\frac{\Delta
t}{2}\right),y_{j}\cos\left(-\pi\frac{\Delta
t}{2}\right)+x_{i+\frac{1}{2}}\sin\left(-\pi\frac{\Delta t}{2}\right)\right)$
(23) $\displaystyle Q_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}$
$\displaystyle=q^{n}\left(-y_{j+\frac{1}{2}}\sin\left(-\pi\frac{\Delta
t}{2}\right)+x_{i}\cos\left(-\pi\frac{\Delta
t}{2}\right),y_{j+\frac{1}{2}}\cos\left(-\pi\frac{\Delta
t}{2}\right)+x_{i}\sin\left(-\pi\frac{\Delta t}{2}\right)\right)$

and analogously for all the other nodes.

We compare numerical solutions of (19) with initial condition (16) using
refinement levels 3-6 and 3-7 after a half rotation with the exact solution
and compute the error as well as the EOC. Results are shown in Table 2 (left).

Table 2: Error at time $t=2$ measured in the 1-norm and EOC for solid-body
rotation (left) using the approach of Section 4.2 and (right) using the
approach of Section 4.3.

Level | Error | EOC
---|---|---
6 | $4.599\,185\text{\times}{10}^{-06}$ | —
7 | $5.826\,258\text{\times}{10}^{-07}$ | $2.9807$
3-6 | $4.674\,628\text{\times}{10}^{-06}$ | —
3-7 | $5.924\,630\text{\times}{10}^{-07}$ | $2.9801$
Level | Error | EOC
---|---|---
6 | $5.791\,702\text{\times}{10}^{-06}$ | —
7 | $7.347\,671\text{\times}{10}^{-07}$ | $2.9786$
3-6 | $5.792\,193\text{\times}{10}^{-06}$ | —
3-7 | $7.351\,174\text{\times}{10}^{-07}$ | $2.9781$

Again, the accuracy of the computations on the adaptively refined grids
compares well with the accuracy obtained on the equidistant grids with highest
resolution. Figure 8 shows the numerical solution with refinement level 3-6
after a half and a full rotation.

Figure 8: Solution to (19) with initial condition (16) after a half (left) and
a full rotation (right).

###### Remark 4.1

The Active Flux method for solid body rotation described in Section 4.2
preserves constant states on regular Cartesian grids.

Proof: We consider constant data at time $t_{n}$, i.e.
$q(x,y,t_{n})=C\in\mathbb{R}$. Then all the point values at time $t_{n}$,
$t_{n+\frac{1}{2}}$ and $t_{n+1}$ are equal to $C$ and the numerical fluxes
have the form

$\begin{split}F_{i+\frac{1}{2},j}&=-\frac{C\pi}{6}\left(y_{j-\frac{1}{2}}+4y_{j}+y_{j+\frac{1}{2}}\right)\\\
G_{i,j+\frac{1}{2}}&=\frac{C\pi}{6}\left(x_{i-\frac{1}{2}}+4x_{i}+x_{i+\frac{1}{2}}\right).\end{split}$

Thus, the fluxes in the $x$-direction only depend on $y$ and the fluxes in the
$y$-direction only depend on $x$. Consequently, all the flux differences
cancel and the cell average values remain constant.

$\square$

###  4.3 The swirl flow problem

Now we consider advective transport in a velocity field that depends on space
and time. For this case, we use

$\partial_{t}q+\partial_{x}\left(a(x,y,t)q\right)+\partial_{y}\left(b(x,y,t)q\right)=0.$
(24)

The velocity field is again defined via a stream function, which now has the
form

$\Psi(x,y,t)=\frac{1}{\pi}\sin^{2}(\pi x)\cdot\sin^{2}(\pi y)\cdot\cos(\pi
t).$

We use finite difference approximations to compute

$\begin{split}a(x,y,t)&=-\frac{\Psi(x,y+\Delta y/2,t)-\Psi(x,y-\Delta
y/2,t)}{\Delta y},\\\ b(x,y,t)&=\phantom{-}\frac{\Psi(x+\Delta
x/2,y,t)-\Psi(x-\Delta x/2,y,t)}{\Delta x}.\end{split}$ (25)

We trace the characteristics back in space and time using the classical fourth
order accurate Runge-Kutta method backwards in time.

For the flux computation of $F_{i+\frac{1}{2},j}$ we use the approximation

$\begin{split}&\frac{1}{\Delta t\Delta
y}\int_{t_{n}}^{t_{n+1}}\int_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}a(x_{i+\frac{1}{2}},y,t)q(x_{i+\frac{1}{2}},y,t)dydt\\\
&\approx\frac{1}{\Delta t\Delta
y}\int_{t_{n}}^{t_{n+1}}\left(a(x_{i+\frac{1}{2}},y_{j},t)\int_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}q(x_{i+\frac{1}{2}},y,t)dy\right)dt\end{split}$
(26)

and analogously for the flux $G_{i,j+\frac{1}{2}}$. Using again Simpson’s rule
we obtain

$\begin{split}F_{i+\frac{1}{2},j}&=\frac{1}{6}\Bigg{(}\frac{1}{6}a(x_{i+\frac{1}{2}},y_{j},t_{n})\cdot\left(Q_{i+\frac{1}{2},j-\frac{1}{2}}^{n}+4Q_{i+\frac{1}{2},j}^{n}+Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n}\right)\\\
&\hskip
28.45274pt+\frac{4}{6}a(x_{i+\frac{1}{2}},y_{j},t_{n+\frac{1}{2}})\cdot\left(Q_{i+\frac{1}{2},j-\frac{1}{2}}^{n+\frac{1}{2}}+4Q_{i+\frac{1}{2},j}^{n+\frac{1}{2}}+Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}\right)\\\
&\hskip
28.45274pt+\frac{1}{6}a(x_{i+\frac{1}{2}},y_{j},t_{n+1})\cdot\left(Q_{i+\frac{1}{2},j-\frac{1}{2}}^{n+1}+4Q_{i+\frac{1}{2},j}^{n+1}+Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n+1}\right)\Bigg{)}\end{split}$
(27)

and

$\begin{split}G_{i,j+\frac{1}{2}}&=\frac{1}{6}\Bigg{(}\frac{1}{6}b(x_{i},y_{j+\frac{1}{2}},t_{n})\cdot\left(Q_{i-\frac{1}{2},j+\frac{1}{2}}^{n}+4Q_{i,j+\frac{1}{2}}^{n}+Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n}\right)\\\
&\hskip
28.45274pt+\frac{4}{6}b(x_{i},y_{j+\frac{1}{2}},t_{n+\frac{1}{2}})\cdot\left(Q_{i-\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}+4Q_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}+Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}\right)\\\
&\hskip
28.45274pt+\frac{1}{6}b(x_{i},y_{j+\frac{1}{2}},t_{n+1})\cdot\left(Q_{i-\frac{1}{2},j+\frac{1}{2}}^{n+1}+4Q_{i,j+\frac{1}{2}}^{n+1}+Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n+1}\right)\Bigg{)}.\end{split}$
(28)

###### Theorem 4.2

The Active Flux method with fluxes of the form (27), (28) and $a,b$ as defined
in (25), provides a method for (24) that preserves constant states on regular
Cartesian grids.

Proof: We consider constant data at time $t_{n}$, i.e.
$q(x,y,t_{n})=C\in\mathbb{R}$. Then all the point values at time $t_{n}$,
$t_{n+\frac{1}{2}}$ and $t_{n+1}$, which are obtained by tracing back the
characteristics, are also equal to $C$. The finite volume update now reduces
to

$\begin{split}Q_{i,j}^{n+1}=Q_{i,j}^{n}&-\frac{\Delta t}{\Delta
x}\frac{C}{6}\Big{(}a(x_{i+\frac{1}{2}},y_{j},t_{n})+4a(x_{i+\frac{1}{2}},y_{j},t_{n+\frac{1}{2}})+a(x_{i+\frac{1}{2}},y_{j},t_{n+1})\\\
&\hskip
28.45274pt-a(x_{i-\frac{1}{2}},y_{j},t_{n})-4a(x_{i-\frac{1}{2}},y_{j},t_{n+\frac{1}{2}})-a(x_{i-\frac{1}{2}},y_{j},t_{n+1})\Big{)}\\\
&-\frac{\Delta t}{\Delta
y}\frac{C}{6}\Big{(}b(x_{i},y_{j+\frac{1}{2}},t_{n})+4b(x_{i},y_{j+\frac{1}{2}},t_{n+\frac{1}{2}})+b(x_{i},y_{j+\frac{1}{2}}t_{n+1})\\\
&\hskip
28.45274pt-b(x_{i},y_{j-\frac{1}{2}},t_{n})-4b(x_{i},y_{j-\frac{1}{2}},t_{n+\frac{1}{2}})-b(x_{i},y_{j-\frac{1}{2}},t_{n+1})\Big{)}.\end{split}$

For the terms at time level $t_{n}$ we obtain

$\begin{split}&-\frac{C\Delta t}{6\Delta
x}\left(a(x_{i+\frac{1}{2}},y_{j},t_{n})-a(x_{i-\frac{1}{2}},y_{j},t_{n})\right)-\frac{C\Delta
t}{6\Delta
y}\left(b(x_{i},y_{j+\frac{1}{2}},t_{n})-b(x_{i},y_{j-\frac{1}{2}},t_{n})\right)\\\
&=-\frac{C\Delta t}{6\Delta x\Delta
y}\Big{(}-\Psi(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}},t_{n})+\Psi(x_{i+\frac{1}{2}},y_{j-\frac{1}{2}},t_{n})+\Psi(x_{i-\frac{1}{2}},y_{j+\frac{1}{2}},t_{n})-\Psi(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}},t_{n})\Big{)}\\\
&\quad-\frac{C\Delta t}{6\Delta y\Delta
x}\Big{(}\Psi(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}},t_{n})-\Psi(x_{i-\frac{1}{2}},y_{j+\frac{1}{2}},t_{n})-\Psi(x_{i+\frac{1}{2}},y_{j-\frac{1}{2}},t_{n})+\Psi(x_{i-\frac{1}{2}},y_{j+\frac{1}{2}},t_{n})\Big{)}\\\
&=0\end{split}$

In the same way the terms at time $t_{n+\frac{1}{2}}$ and $t_{n+1}$ cancel and
we obtain $Q_{i,j}^{n+1}=Q_{i,j}^{n}=C$ for all $i,j$.

$\square$

Now we consider Cartesian grids with adaptive mesh refinement. Without loss of
generality we restrict considerations to the situation shown in Figure 9.

Figure 9:  Coarse grid cell with neighboring fine grid cells.

We again consider constant data at time $t_{n}$, i.e. $q(x,y,t_{n})=C$. The
update of the small grid cells preserves by Theorem 4.2 constant states. For
the large cell update, the left flux $F_{i-\frac{1}{2},j}$ is the sum of two
fluxes used to update the two small cells, i.e.

$\begin{split}F_{i-\frac{1}{2},j}&=\frac{1}{2}\left(F_{i-\frac{1}{2},j-\frac{1}{4}}+F_{i-\frac{1}{2},j+\frac{1}{4}}\right)\\\
&=\frac{1}{2}C\Big{[}\frac{1}{6}\left(a(x_{i-\frac{1}{2}},y_{j-\frac{1}{4}},t_{n})+4a(x_{i-\frac{1}{2}},y_{j-\frac{1}{4}},t_{n+\frac{1}{2}})+a(x_{i-\frac{1}{2}},y_{j-\frac{1}{4}},t_{n+1})\right)\\\
&\hskip
28.45274pt+\frac{1}{6}\left(a(x_{i-\frac{1}{2}},y_{j+\frac{1}{4}},t_{n})+4a(x_{i-\frac{1}{2}},y_{j+\frac{1}{4}},t_{n+\frac{1}{2}})+a(x_{i-\frac{1}{2}},y_{j+\frac{1}{4}},t_{n+1})\right)\Big{]}\\\
&=\frac{1}{2}C\frac{1}{6}\Big{[}-\frac{\Psi(x_{i-\frac{1}{2}},y_{j},t_{n})-\Psi(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}},t_{n})}{\Delta
y/2}-4\frac{\Psi(x_{i-\frac{1}{2}},y_{j},t_{n+\frac{1}{2}})-\Psi(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}},t_{n+\frac{1}{2}})}{\Delta
y/2}\\\
&-\frac{\Psi(x_{i-\frac{1}{2}},y_{j},t_{n+1})-\Psi(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}},t_{n+1})}{\Delta
y/2}-\frac{\Psi(x_{i-\frac{1}{2}},y_{j+\frac{1}{2}},t_{n})-\Psi(x_{i-\frac{1}{2}},y_{j},t_{n})}{\Delta
y/2}\\\
&-4\frac{\Psi(x_{i-\frac{1}{2}},y_{j+\frac{1}{2}},t_{n+\frac{1}{2}})-\Psi(x_{i-\frac{1}{2}},y_{j},t_{n+\frac{1}{2}})}{\Delta
y/2}-\frac{\Psi(x_{i-\frac{1}{2}},y_{j+\frac{1}{2}},t_{n+1})-\Psi(x_{i-\frac{1}{2}},y_{j},t_{n+1})}{\Delta
y/2}\Big{]}\\\
&=\frac{1}{6}C\Big{[}-\frac{\Psi(x_{i-\frac{1}{2}},y_{j+\frac{1}{2}},t_{n})-\Psi(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}},t_{n})}{\Delta
y}\\\
&-4\frac{\Psi(x_{i-\frac{1}{2}},y_{j+\frac{1}{2}},t_{n+\frac{1}{2}})-\Psi(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}},t_{n+\frac{1}{2}})}{\Delta
y}-\frac{\Psi(x_{i-\frac{1}{2}},y_{j+\frac{1}{2}},t_{n+1})-\Psi(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}},t_{n+1})}{\Delta
y}\Big{]}\\\
&=C\frac{1}{6}\left(a(x_{i-\frac{1}{2}},y_{j},t_{n})+4a(x_{i-\frac{1}{2}},y_{j},t_{n+\frac{1}{2}})+a(x_{i-\frac{1}{2}},y_{j},t_{n+1})\right).\end{split}$

This is the same flux that would be used on a regular Cartesian grid. Thus,
the update of the small grid cells as well as the update of the coarse grid
cells preserve by Theorem 4.2 constant states. We summarize our result.

###### Corollary 4.3

The Active Flux method described above preserves constant states on Cartesian
grids with adaptive mesh refinement when used with global time stepping (e.g.
no subcycling).

With subcycling, the sum of the fluxes from the fine grid cells and that is
used to update the coarse grid cell would contain components at the
intermediate time that are in general not balanced by the remaining coarse
grid fluxes.

To perform a numerical convergence study, we consider smooth initial values of
the form

$q(x,y,0)=\exp(-100((x-0.5)^{2}+(y-0.25)^{2}))+\exp(-100((x-0.5)^{2}+(y-0.75)^{2}))$
(29)

on the domain $[0,1]\times[0,1]$ with periodic boundary conditions.

Figure 10: Solution of the swirl problem with initial data of the form (29) at
times $t=0$ (left), $t=0.5$ (middle) and $t=1$ (right) for the swirl problem
with smooth initial condition.

Numerical results are shown in Figure 10 and the results of a convergence
study with and without subcycling is shown in Table 3. The results again
confirm third order convergence. The use of subcycling did not significantly
influence the accuracy.

Table 3: Error at time $t=1$ measured in the 1-norm and EOC for the smooth
swirl problem (left) with subcycling, (right) without subcycling.

Level | Error | EOC
---|---|---
4 | $5.722\,927\text{\times}{10}^{-06}$ | —
5 | $7.268\,088\text{\times}{10}^{-07}$ | $2.9771$
3-4 | $5.758\,274\text{\times}{10}^{-06}$ | —
3-5 | $7.346\,642\text{\times}{10}^{-07}$ | $2.9705$
Level | Error | EOC
---|---|---
4 | $5.722\,927\text{\times}{10}^{-06}$ | —
5 | $7.268\,088\text{\times}{10}^{-07}$ | $2.9771$
3-4 | $5.724\,930\text{\times}{10}^{-06}$ | —
3-5 | $7.305\,236\text{\times}{10}^{-07}$ | $2.9703$

To measure the error in the approximation of constant states by using
subcycling, we compute the solution at time $t=1$ using the same velocity
field but with constant initial values that are equal to one in the whole
domain. Refinement with levels $3-4$ is used along the diagonal as shown in
Figure 7 (left). With subcycling we observe an error of size $10^{-13}$.
Without subcycling the error is about $10^{-16}$, i.e. agrees with the
expected truncation error.

We now also consider the swirl problem with piecewise constant initial values
of the form

$q(x,y,0)=\begin{cases}1&0\leq x\leq 0.5\\\ 0&0.5<x\leq 1\end{cases}$

on the domain $[0,1]\times[0,1]$ with periodic boundaries on top and bottom
and with zero-order extrapolation on the left and right boundary using
refinement level 3-6. Results are shown in Figure 11.

Figure 11:  Solution of the swirl problem with piecewise constant initial
values at times $t=0$ (left), $t=0.5$ (middle) and $t=1$ (right).

Subcycling was used although this slightly perturbs the exact preservation of
constant states.

While the method from Section 4.2 preserves constant states for solid body
rotation due to the special form of the velocity field, we can also compute
fluxes for solid body rotation using the more general constant state
preserving approach presented in this section. Numerical results are shown in
Table 2 (right). The error is slightly larger, due to the use of the
approximation described in Equation (26), but the method clearly preserves
third order accuracy.

###  4.4 Burgers’ equation

It is straight forward to apply the adaptive Active Flux method to scalar
nonlinear hyperbolic problems. We consider the Burgers’ equation (7) on the
domain $[0,1]\times[0,1]$ with initial values

$q(x,y,0)=\sin(2\pi x)\sin(2\pi y)+0.1.$ (30)

To check the accuracy we compute numerical solutions at an early time, here we
use $t=0.05$, at which the solution structure is still smooth. The time steps
satisfy $\mbox{CFL}\leq 0.5$. Furthermore, we use subcycling and the
conservative fix. To test the accuracy of the adaptive method we enforce
refinement along the diagonal from the upper left patch to the lower right
patch. This test (not shown here) confirms third order accuracy also for this
nonlinear problem.

At later times shocks arise and we therefore also used the bound preserving
limiter introduced in [10]. Note that the characteristic speed changes sign
which has been observed to lead to some numerical difficulties as explained in
[14, 10]. Our unlimited method shows some unphysical oscillations along the
curve where the characteristic speeds changes sign as can be seen in Figure 12
(top). The use of the bound preserving limiter avoids these inaccuracies. In
this simulation grids on levels $3-5$ are used.

As a refinement criteria we search for steep gradients, so a patch is refined,
if

$\frac{q_{i+1,j}-q_{i-1,j}}{2\Delta x}\geq 15\text{ \ or \
}\frac{q_{i,j+1}-q_{i,j+1}}{2\Delta y}\geq 15\ ,$

for any $i,j$. We also use the refinement criteria to determine whether we
need to apply limiting.

Figure 12: Solution to (30) at times $t=0.2$ (left), $t=0.6$ (middle) and
$t=1.0$ (right) without limiter (top) and with limiter (bottom) for refinement
level $3-5$.

###  4.5 Acoustics

####  4.5.1 Convergence study for acoustics

In order to investigate the order of convergence of the adaptive Active Flux
method for acoustics, we consider a test problem from Lukácǒvá et al. [16],
for which the exact solution is explicitly known. In this test problem the
acoustic equations (11) with initial values of the form

$\displaystyle p(x,y,0)$ $\displaystyle=-\frac{1}{c}\left((\sin(2\pi
x)+\sin(2\pi y)\right)$ (31) $\displaystyle u(x,y,0)$ $\displaystyle=0$
$\displaystyle v(x,y,0)$ $\displaystyle=0$

are considered on the domain $[-1,1]\times[-1,1]$. The speed of sound is set
to $c=1$ and periodic boundary conditions are imposed. The exact solution has
the form

$\begin{split}p(x,y,t)&=-\frac{1}{c}\cos(2\pi ct)\left(\sin(2\pi x)+\sin(2\pi
y)\right)\\\ u(x,y,t)&=\frac{1}{c}\sin(2\pi ct)\cos(2\pi x)\\\
v(x,y,t)&=\frac{1}{c}\sin(2\pi ct)\cos(2\pi y).\end{split}$

We compute numerical solutions at time $t=1$ using time steps, which satisfy
the condition $\mbox{CFL}\leq 0.4$. For this smooth solution, adaptive mesh
refinement does not really make sense. However, in order to test the accuracy
of the adaptive method, we require refinement along the diagonal of the domain
independently of the solution structure. Each patch is a Cartesian grid with
$16\times 16$ grid cells.

Figure 13: Pressure and grid patches for the acoustics problem from Lukácǒvá
et al. [16] with refinement levels 2-3 (left) and 3-4 (right) at $t=1.0$.

We compare grids with refinement level $2-3$ and $3-4$ with results obtained
on a uniform grid with level $2$ and $3$. The results of our convergence study
for pressure and velocity are shown in Table 4.

Table 4: Error at time $t=1$ measured in $\|\cdot\|_{1}$-norm and EOC for the Lukácǒvá test problem. Level | Error | EOC
---|---|---
| $p$ | $u,v$ | $p$ | $u,v$
2 | $2.231\,166\text{\times}{10}^{-04}$ | $2.520\,880\text{\times}{10}^{-05}$ | — | —
3 | $2.793\,825\text{\times}{10}^{-05}$ | $2.648\,027\text{\times}{10}^{-06}$ | $2.9975$ | $3.2509$
2-3 | $1.641\,979\text{\times}{10}^{-04}$ | $3.328\,607\text{\times}{10}^{-05}$ | — | —
3-4 | $2.057\,597\text{\times}{10}^{-05}$ | $3.869\,926\text{\times}{10}^{-06}$ | $3.1045$ | $2.8564$

The accuracy obtained on the adaptively refined mesh is comparable with the
accuracy obtained on a regular grid which uses the coarser grid in the whole
domain. This test shows that the the accuracy is maintained at the interface
between the coarse and the fine grids.

####  4.5.2 High frequency acoustics

Now we consider a test problem where adaptive mesh refinement allows a more
efficient computation of the solution structure. We consider the acoustic
equations (11) with initial condition

$p(x,y,0)=2+\exp(-100(r(x,y)-0.5)^{2}))\sin(100r(x,y)),\ \ u(x,y,0)=0,\ \
v(x,y,0)=0$ (32)

on the domain $[-1.5,1.5]\times[-1.5,1.5]$. Initially a circular shaped
acoustic wave with high frequency pressure oscillations is given. This leads
to acoustic waves moving outwards and inwards. The inward moving acoustic wave
gets reflected in the center of the domain (around the time shown in the third
plot) and afterwards propagates outwards. Adaptive mesh refinement is used to
resolve these high frequency waves. Figure 14 shows solutions at times
$t=0,0.3,0.6,0.9$ with refinement level $3-6$. A patch is refined if
$|p_{max}-p_{min}|\geq 0.001$.

Figure 14: Plots of pressure of the high frequency acoustic wave at times
$t=0.0$ (top left), $t=0.3$ (top right), $t=0.6$ (bottom left) and $t=0.9$
(bottom right). In the bottom left plot the pressure at the center is much
larger and therefore a different color map is used to visualize the solution
structure.

In Figure 15 shows scatter plots of the solution at two different times. There
are no visable spurious grid effects apart from some smearing of the solution
structure caused by the boundary of the computational domain. For this
simulation periodic boundary conditions were used.

Figure 15: Scatter plots of pressure for the high frequency acoustic wave at
time $t=0.3$ and $t=0.9$.

###  4.6 Approximation of steady states

Finally, we consider the approximation of a stationary vortex as described by
Barsukow et al. [1, 2]. The initial values have the form

$\begin{split}p(r,0)&=0\\\
\boldsymbol{u}(r)&=\boldsymbol{n}\left\\{\begin{array}[]{ccc}5r&:&0\leq r\leq
0.2\\\ 2-5r&:&0.2<r\leq 0.4\\\ 0&:&r>0.4,\end{array}\right.\end{split}$

with $r=\sqrt{x^{2}+y^{2}}$, $\boldsymbol{n}=(-\sin\phi,\cos\phi)^{T}$,
$\phi\in[0,2\pi)$ and $\boldsymbol{u}=(u,v)^{T}$. In order to test how well
the method preserves the steady state, we compute the numerical solution at
time $t=100$. The AMR computation uses grids of level 3-5 with $16\times 16$
grid cells on each patch.

In [1, 2], the authors showed that the Cartesian grid Active Flux method, with
the evolution operator presented in [2], is stationary preserving. The method
used here does not have this property. Nevertheless, the third order accuracy
of the Active Flux method together with adaptive mesh refinement leads to
accurate approximations as shown in Figure 16.

Figure 16:  Computation of the stationary vortex at time $t=100$ using an
adaptively refined grid of level 3-5. The left plot shows $|\boldsymbol{u}|$,
the right plot shows a scatter plot of $|\boldsymbol{u}|$.

We also tested the evolution operator from [2]. While we could confirm the
preservation of steady states on a regular Cartesian grid, using the solver
with AMR leads to an instability.

## Conclusions

We showed that the Active flux method can be used on Cartesian grids with
adaptive mesh refinement and subcycling. The transfer of grid information
between different Cartesian grid patches can be implemented without loss of
third order accuracy by making use of the degrees of freedom of the Active
Flux method. Our approach benefits from the local stencil of the Active Flux
method.

For advective transport, new Active Flux methods have been presented which
preserve constant states on regular Cartesian grids and grids with adaptive
mesh refinement without subcycling. In practical computations with subcycling
we also observed good accuracy although constant states are not exactly
preserved.

The AMR concept of the Active Flux method can also be used for two-dimensional
linear hyperbolic systems as illustrated for the acoustic equations. However,
the preservation of steady states, a property that was recently shown for the
Active Flux method on regular Cartesian grids [1], does not carry over to
adaptively refined meshes with the approaches presented in this paper.

## Data availability statement

For our implementations we used ForestClaw, which is publicly available on
GitHub. We added the Active Flux method as a new solver in ForestClaw.

## References

  * [1] W. Barsukow.  Stationarity preservation properties of the active flux scheme on cartesian grids.  Commun. Appl. Comput., 2020.
  * [2] W. Barsukow, J. Hohm, C. Klingenberg, and P.L. Roe.  The active flux scheme on Cartesian grids and its low Mach number limit.  J. Sci. Comput., 81(1):594–622, 2019.
  * [3] Marsha J. Berger and Joe Oliger.  Adaptive mesh refinement for hyperbolic partial differential equations.  J. Comput. Phys., 53:484–512, 1984.
  * [4] M.J. Berger and P. Colella.  Local adaptive mesh refinement for shock hydrodynamics.  J. Comput. Phys., 82:484–512, 1984.
  * [5] M.J. Berger and R.J. LeVeque.  Adaptive mesh refinement using wave-propagation algorithms for hyperbolic systems.  SIAM J. Numer. Anal., 35:2298–2316, 1998.
  * [6] C. Burstedde, D. Calhoun, K. Mandli, and A.R. Terrel.  Forestclaw: Hybrid forest-of-octrees amr for hyperbolic conservation laws.  Parallel Computing: Accelerating Computational Science and Engineering (CSE), 25:253–262, 2014.
  * [7] C. Burstedde, L.C. Wilcox, and O. Ghattas.  p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees.  SIAM J. Sci. Comput., 33:1103–1133, 2011.
  * [8] Donna Calhoun and Carsten Burstedde.  Forestclaw software.  http://www.forestclaw.org/ForestClaw/index.html.
  * [9] Donna Calhoun and Carsten Burstedde.  ForestClaw : A parallel algorithm for patch-based adaptive mesh refinement on a forest of quadtrees.  arXiv:1703.03116, 2017.
  * [10] E. Chudzik, C. Helzel, and D. Kerkmann.  The Cartesian grid active flux method: linear stability and bound preserving limiting.  Appl. Math. Comput., 393:125501, 19, 2021.
  * [11] T.A. Eymann and P.L. Roe.  Active Flux Schemes.  AIAA 2011-382.
  * [12] T.A. Eymann and P.L. Roe.  Active Flux Schemes for Systems.  AIAA 2011-3840.
  * [13] T.A. Eymann and P.L. Roe.  Multidimensional Active Flux Schemes.  AIAA Conference Paper, June 2013.
  * [14] C. Helzel, D. Kerkmann, and L.Scandurra.  A new ADER method inspired by the active flux method.  J. Sci. Comput., 80(3):1463–1497, 2019.
  * [15] R.J. LeVeque.  Finite Volume Methods for Hyperbolic Problems.  Cambridge University Press, 2002.
  * [16] M. Lukáčová-Medvid’ová, K.W. Morton, and G. Warnecke.  Evolution galerkin methods for hyperbolic systems in two space dimensions.  Mathematics of Computation, 69:1355–1384, 2000.
  * [17] M. Lukáčová-Medvid’ová, J. Saibertová, and G. Warnecke.  Finite volume evolution Galerkin methods for nonlinear hyperbolic systems.  J. Comput. Phys., 183:533–562, 2002.
  * [18] P. Roe.  Is discontinuous reconstruction really a good idea?  J. Sci. Comput., 73(2-3):1094–1114, 2017.
  * [19] P. Roe.  Designing CFD methods for bandwidth - a physical approach.  Computers & Fluids, 214:104774, 2021.
  * [20] P.L. Roe, J. Maeng, and D. Fan.  Comparing Active Flux and Discontinuous Galerkin Methods for Compressible Flow.  2018 AIAA Aerospace Science Meeting.
  * [21] Weiqun Zhang, Ann Almgren, Vince Beckner, John Bell, Johannes Blaschke, Cy Chan, Marcus Day, Brian Friesen, Kevin Gott, Daniel Graves, Max P. Katz, Andrew Myers, Tan Nguyen, Andrew Nonaka, Michele Rosso, Samuel Williams, and Michael Zingale.  AMReX: a framework for block-structured adaptive mesh refinement.  J. Open Source Software, 4(37), 2019.