Adaptive Mesh Refinement Using Wave-Propagation Algorithms for Hyperbolic Systems.
An adaptive mesh refinement algorithm developed for the Euler equations of gas dynamics has been extended to employ high-resolution wave-propagation algorithms in a more general framework.  This allows its use on a variety of new problems, including hyperbolic equations not in conservation form, problems with source terms or capacity functions, and logically rectangular curvilinear grids.  This framework requires a modified approach to maintaining consistency and conservation at grid interfaces, which is described in detail.  The algorithm is implemented in the AMRCLAW package, which is freely available.
Introduction
The multi-dimensional wave-propagation algorithm described in [14] is a "high-resolution" method
that is second order accurate on smooth solutions while maintaining sharp discontinuities through
the use of slope-limiters. While based on ideas developed for hyperbolic systems of conservation
laws
in the context of shock capturing, these methods apply in a more general framework that allows their
application to other hyperbolic systems which are not in conservation form. Variable-coefficient
hyperbolic systems of the form
arise, for example, in studying acoustics or elasticity in heterogeneous materials with varying material
properties (see Section 6).
For most practical problems, it is desirable to use mesh refinement to cluster grid points in
regions where they are most needed, for example around shocks or other regions where the solution
Courant Institute, New York University, 251 Mercer St., New York, NY 10012. berger@cims.nyu.edu
2 Department of Applied Mathematics and Department of Mathematics, University of Washington, Box 352420,
Seattle,
January, 1997. To appear in a special issue of SIAM J. Numer. Anal. dedicated to Ami Harten.
Revised July, 1997.
has steep gradients or complicated structure. This should be done in an adaptive manner, based on
the behavior of the solution, and for time-dependent problems the region of refinement must move
adaptively with the interesting structure. An effective adaptive mesh refinement (AMR) strategy
has been developed [1-6] that uses Cartesian grids, with refinement in both space and time, over
rectangular patches. The refinement is by an arbitrary even integer ratio (typically 4) and further
recursive refinement can be done within these patches to an arbitrary depth. This algorithm was
originally developed for the Euler equations of gas dynamics using flux-differencing methods, but
can be easily extended to other systems of conservation laws. A crucial ingredient is the manner
in which fluxes at grid-refinement interfaces are coordinated to insure that the method is globally
conservative. This is described in detail in [5] and reviewed in Section 2.1.
In this paper we show how the wave-propagation algorithm can be used in conjunction with this
AMR strategy. This algorithm is written in a more general form that does not use flux differencing
per se, though in the special case of a conservation law the method can be re-expressed in flux-
differencing form and is fully conservative. In Section 4 we show that at grid-refinement interfaces
it is possible to apply a correction procedure that maintains global conservation when applied on
a conservation law (and reduces to the flux-based approach), while retaining the more general
wave-propagation framework that allows application on more general hyperbolic systems.
We also discuss how further extensions presented in [14] can be incorporated into AMR. This
includes source terms, capacity-form differencing, and applications on mapped curvilinear grids.
These topics are discussed in Sections 5 and 6.
The domain is assumed to be rectangular, at least in computational space. Metric terms
can be included so that we are actually solving a problem on a nonrectangular physical domain
with a curvilinear grid, obtained by a smooth mapping of the Cartesian computational grid. The
grid Jacobian function can be properly incorporated using the "capacity function" described in
Section 5.2.
The wave-propagation algorithm has been implemented in a general software package (in for-
called clawpack (Conservation LAWs PACKage, a holdover from earlier versions that
applied only to conservation laws). This package is available from netlib[11]. This software
has recently been combined with Berger's implementation of AMR, incorporating the extensions
described in this paper. This is now freely available as the amrclaw package[6]. All of the numerical
results presented in this paper were obtained with amrclaw, and several such examples
are included with the package. Details on the use of this package can be found in the on-line
documentation.
2 The AMR algorithm
The adaptive mesh refinement algorithm for conservation laws is fully described in [5] and only
a brief summary will be given here. The AMR approach to adaptive mesh refinement uses a
collection of logically rectangular meshes that make up the coarse grid; refinements cover a subset
of the domain and use smaller rectangular grid patches. These fine patches can be recursively
nested until a given level of accuracy is attained. Typically, if a patch at Level L is refined in x
and y by an even integer RL , then the time step is also refined by the same factor, so that RL
time steps must be taken on the refined grid at Level L each step on the grids at Level L.
The mesh ratios \Deltat/\Deltax and \Deltat/\Deltay are then the same on all grids, ensuring stability with explicit
difference schemes.
Every K time steps on a particular grid level, all finer level grids are regenerated in order to follow
moving features of the flow. An error estimation procedure based on Richardson extrapolation
determines the regions where resolution of the solution is insufficient. This procedure compares the
solution obtained by taking 2 steps on the existing grid with one computed by taking 1 time step on
a grid that is twice as coarse in each direction. Cells where the error is greater than some tolerance
are flagged for refinement. Other criteria might be used in addition to, or instead of, this error
estimate, e.g. identifying steep gradients in some variable. A buffer zone around the flagged cells
is also flagged to insure that features of interest do not escape from the refinement region over the
next K time steps. The buffer width and K are adjustable parameters which must be coordinated.
Flagged cells are then organized into rectangular grid patches, typically containing several
hundred to several thousand grid points per patch. Note that some cells not tagged for refinement
are also included in new fine grid patches. Typically, our grid generation algorithm produces grids
with 70% of the cells within new grid patches tagged as needing refinement; the remaining 30%
are untagged but still lie within the new patch boundaries. By taking very small patches one
could avoid refining too many cells, but this must be balanced with the competing desire to create
relatively few separate patches, and to minimize computational overhead on the boundaries of fine
grids. See [1] for more details on the refinement and clustering algorithms. See Figure 6 for an
example of refinement on logically rectangular patches.
An alternative would be to use a quadtree data structure (e.g., [7], [15]) in which only the
flagged cells are refined, but the storage overhead of these data structures, typically 30 to 50 words
per cell, usually exceeds the storage overhead associated with the block structured approach, even
with the 30% additionally refined grid.
A finite volume method is used to advance the solution on the resulting grid hierarchy. Cell
averages of each variable are stored in each grid cell. When solving a conservation law, these cell
averages are updated by a flux-differencing algorithm based on fluxes through the cell edges. (In
section 3 we review wave-propagation methods and their extension to nonconservative hyperbolic
systems.) The integration proceeds by grid level. All grids on level 1 are first integrated over a
time step. Then grids at level 2 are integrated over R 1 time steps to catch up. This approach is
applied recursively on each level.
Boundary conditions on all grids are imposed using "ghost cells". The computational domain
is extended by G ghost cells in each direction, and values are assigned to the ghost cells at the start
of each time step. (G=2 for the methods described here). At a physical boundary, the user must
set the ghost cell values at each time step based on the problem specification (e.g. extrapolation
at an outflow boundary or reflection at a solid wall). However, the boundary of a fine grid may
be interior to the domain. In this case, if there is no neighboring fine grid to supply values for
the ghost cells, they are interpolated from a coarser parent grid. Space-time interpolation must in
general be used since more time steps are taken on the fine grid than on the coarse grid, and at
intermediate times there are no coarse grid values available. However, since coarse grids are always
advanced first, we have data available from both an earlier and later time on the coarse grid from
which to interpolate.
The techniques developed in this paper for adaptive refinement are most easily described in one
space dimension, where the refined grid patches become intervals. We concentrate on this case in
the development below, but all of these techniques carry over directly to multiple space dimensions,
and have been implemented in two dimensions in the amrclaw software. Two-dimensional results
are presented in Section 6.
2.1 Conservation at grid interfaces
At an interface between coarse and fine grids, we must also insure that the formulas used to update
the solution on each grid are consistent with one another. In particular, when a conservation law
is being solved we must preserve global conservation.
To illustrate the conservative flux correction required at grid interfaces, and the manner in which
this must be modifed for the wave-propagation algorithms, we will consider the one-dimensional
case with only two levels of grid and a factor of two refinement. Denote the coarse grid spacing by
h and the time step by k. The grid spacing on the fine grid is -

Figure

1(a)
shows the coarse grid over a single time step in x-t space. Figure 1(b) shows also the finer grid that
is overlaid. We will denote the values on this fine grid by - q i , and assume that there are m cells on
this grid, which ends just before cell j on the coarse grid.
On the coarse grid we compute fluxes F 0
i as indicated in Figure 1(a), and update the coarse
grid values by the flux-differencing algorithm
Note that F i represents the flux at the left edge of cell i, which is the interface between grid cells
instead of the commonly-used but more cluttered notation F i\Gamma1=2 .
On the fine grid we use the fluxes indicated in Figure 1(c). In each of the two time steps on
this grid, we will use flux-differencing of the form
2. To compute fluxes near the right edge of this grid, we will need to use
values in ghost cells which are also indicated in Figure 1(c). These values are determined using
space-time interpolation from the coarse grid values q 0
, as described above.
a)
PSfrag replacements
PSfrag replacements
c)
PSfrag replacements

Figure

1: a) Coarse grid for the one-dimensional example, shown in space-time. Values of q and
fluxes F are indicated. b) Also shows the fine grid that is overlaid for and the grid interface. Values
of -
q on the fine grid are indicated. c) The flux values -
F needed on the fine grid. The ghost cells
q m+1 and -
q m+2 are also indicated. These values are needed to compute the fine grid fluxes at the
grid interface.
We now need to coordinate the values obtained on the two different grids at the final time.
First, in the coarse grid cells overlapped by fine grid cells (i we replace the value q 1
i by
the average of the fine-grid values:
This is sensible since the fine-grid values are presumably more accurate than the coarse-grid value,
and is also crucial in maintaining global conservation when regridding, if the fine grid is eliminated,
and important in maintaining stability since information must be allowed to pass from the fine to
coarse grid as well as in the other direction.
To be fully conservative, however, we also must modify the coarse-grid value q 1
j . In initially
computing q 1
j we used the coarse-grid flux F 0
j at the left edge of this cell. To be conservative, we
must instead use a left-edge flux that agrees with the flux used in determining the fine-grid values
to the left of this cell. Instead of using F 0
, we should use 1( -
j is corrected
by the difference between the two:
We then have global conservation of the total mass at multiples of the coarse grid time step in the
sense that the total mass
i-m
is conserved up to boundary effects at the farfield.
3 The wave-propagation algorithms
In one space dimension, the wave-propagation algorithm described in [14] is based on solving a
Riemann problem at each interface between grid cells, and using the resulting wave structure to
update the solution in the grid cell to each side. This is, of course, the basis for a host of methods
for conservation laws, dating back to Godunov's method[8]. The Riemann problem consists of
the original conservation law together with piecewise constant initial data defined by the two
neighboring cell values. For a wide class of conservation laws, the Riemann problem can be solved
(either exactly or approximately) and the solution is a similarity solution consisting of a set of
waves propagating at constant speeds. For a system of conservation laws, this solution can be used
to define a flux at the cell interface, yielding Godunov's method. Second order accuracy can be
achieved in various ways, (e.g. by introducing slope information), and "slope-limiters" or "flux-
limiters" are then used to give good resolution of discontinuities without spurious oscillations. See,
e.g., [12] for a general discussion of such methods.
The wave-propagation algorithms are based on using the waves directly to update cell values,
including second order corrections with "wave limiters". For conservation laws these methods can
be rewritten in conservation form by defining flux functions in terms of the waves, but they are
implemented in a way that allows their application to hyperbolic problems not in conservation
form, for which there is still a well-defined wave structure but no flux function. To illustrate this,
consider the advection equation with variable velocity,
This equation, sometimes called the "color equation", is not in conservation form. The value of q is
constant along characteristics but the integral of q is not conserved. Assume u(x) ? 0 everywhere.
At the interface x i\Gamma1=2 between cells we can define the Riemann problem as
where the function - u(x) and initial data - q are given by
where u are cell-centered values of u(x) in cells In the solution to this
Riemann problem, the wave W i simply propagates with speed u i (to the right
since u ? 0). Over time step k this wave moves distance ku i into cell i, and modifies the cell average
The first order upwind method, in wave-propagation form, is thus
More generally, if u(x) has arbitrary sign, the first order wave-propagation algorithm is
where
Alternatively, the Riemann problem might be defined by using an edge value u i\Gamma1=2 at the
interface between cells i. The update formula (7) would then become
This has advantages in two-dimensional incompressible flow (see [13]), but for illustration here we
will use the formulation (7), which is based on the cell-centered velocites (6) rather than edge
values.
High-resolution second-order corrections are easily introduced, and in fact this can be done
in a "flux-differencing" form even for nonconservative equations such as the color equation. This
simplifies the procedure for insuring conservation when applied to a conservation law. At the cell
interface x i\Gamma1=2 , we define
~
~
where ~
W i is a limited version of the wave W
~
where limiter(a; b) represents some standard limiter such as minmod or superbee. After including
these corrections, the method (7) becomes
Note that if u is constant then the advection equation (5) is a conservation law, and in this case
reduces to (assuming u ? 0, for example),
where
This is the numerical flux for a standard flux-limiter method on the advection equation (see [19],
for example). In particular, if no limiter is used and ~
reduces to the Lax-Wendroff
method.
The more general form (9) is not in flux-differencing form, but works just as effectively on the
color equation as standard flux-differencing does on the constant-coefficient advection equation.
The form (9) is easily extended to general hyperbolic systems. Consider the variable-coefficient
linear system
where now q 2 lR m and A(x) 2 lR m\Thetam is diagonalizable with real eigenvalues. Then for the
Riemann problem at x i\Gamma1=2 we decompose q
m) in such
a way that left-going waves are eigenvectors of the matrix A i\Gamma1 defined on cell traveling with
speeds - p
(eigenvalues of A i\Gamma1 ), while the right-going waves are eigenvectors of A i traveling
with speeds - p
The summation term gives the first-order upwind method, while the ~
are again the second-order
corrections defined now by
~
~
where ~
W i+1 is a limited version of the wave W p
obtained by comparing it with W p
or with W p
details, including a worked example for acoustics in a
heterogeneous medium, can be found in [14] and the documentation with [11].
The wave-propagation algorithm is extended to nonlinear systems of conservation laws using
a Roe approximate Riemann solver[16], which linearizes the problem at each cell interface in such
a way that the wave-propagation approach is guaranteed to be conservative. The general wave-propagation
algorithm is written symbolically as
represents the right-going "fluctuation" from the i'th Riemann problem, at the left
edge of Cell i, while A \Gamma \Deltaq i+1 is the left-going fluctuation from the Riemann problem at the right
edge of this cell. Each fluctuation is just the sum over all waves moving in the appropriate direction
of the wave speed multiplied by the wave strength. The notation is motivated by the fact that, for
the constant-coefficient linear system q t
and
where A being the eigendecomposition of A, so that
is the matrix of right eigenvectors. The waves W p in this case are given by W where the
vector of wave strengths is
Note that in this linear case A A. For a general conservation law, the method (11) is
conservative provided that
i.e., the fluctuations are defined by a flux-difference splitting. This is easy to see since we can then
define
so that
Using this in (11) shows that the method can be rewritten in conservation form, with flux F
4 Wave-propagation at grid interfaces
At the interface between a fine and coarse grid, the wave-propagation form can still be used to
update the values on each grid independently, using "ghost cell" values as needed near grid inter-
faces. We can still replace the new coarse grid value by an average of fine grid values in any cell
covered by a fine grid. The only tricky part in extending the AMR algorithm to wave-propagation
algorithms is the conservative correction of the coarse cells adjacent to finer grids, e.g., the value
in

Figure

1(b). Recall that with the flux-differencing algorithm this value must be modified by
the correction (4) to insure conservation, since then the flux "into" the coarse cell agrees with the
total flux "out of " the adjacent fine cells. With the wave-propagation algorithm, we must apply a
similar fix-up to insure that the waves match up in an appropriate manner to yield conservation
when conservation is expected. This turns out to be only slightly more difficult when no numerical
flux is available than for the flux-differencing form discussed in Section 2.
Note that the second order correction terms are written in flux-differencing form ~
even for nonconservative systems, and so these terms can again be corrected using (4). It is only
the first-order upwind terms written in terms of the fluctuations A \Gamma \Deltaq and A + \Deltaq that must be
handled differently, and so we concentrate on the first-order alogrithm below.
Both the difficulty and the solution can be most easily seen by examining the constant coefficient
advection on the one-dimensional grid of Figure 1, in the case u ! 0, so
that waves are moving from the coarse grid to the fine grid. Suppose we are solving the Cauchy
problem and the data has compact support while the fine and coarse grids extend off to \Gamma1 and
respectively. Then we hope to have conservation in the sense that
i-m
i-m
On the coarse grid we have
and when we sum this over the coarse grid cells we get a telescoping of the q-differences everywhere
except in the first cell, so that
Similarly, after 2 steps on the fine grid, we find that
i-m
i-m
Summing these two results we see that the method is globally conservative only if we add in a
correction to this global sum of magnitude
This can be accomplished by modifying the value q 1
by
Note that this is exactly the modification (4) obtained when the flux-differencing form is used, since
for the advection equation
Instead of relying on a flux function, we will interpret this correction in another way using wave
propagation. Figure 2(a) shows all of the waves which affect the relevant cell values at the end
of the time increment, if we only apply the wave-propagation algorithm on each grid separately.
something is wrong with this picture. If we view the initial data as defining a piecewise
constant function, then we need to solve the Riemann problem at each discontinuity and add the
total fluctuation from the Riemann problem \Deltaq) to some grid value(s) in order to
maintain conservation. In general A \Gamma \Deltaq is added to the cell on the left, and A + \Deltaq to the cell
on the right. As Figure 2(a) shows, we have failed to solve the Riemann problem between states
(the ghost-cell value on the fine grid) and q 0
j at the initial time. To restore conservation we
must solve the Riemann problem between these states and add in the resulting total fluctuation
weighted by - since the time step is - k while the cell size is h, to the cell
value
j . In terms of maintaining conservation, this sort of correction could equally well be added
to some other cell value(s) instead of to q 1
j , but this choice is clearly most reasonable from the
figure and agrees with how modifications are applied in the flux-differencing framework.
Similarly, in the second step on the fine grid we must also solve a Riemann problem between
m+1 and q 0
j and add these fluctuations into q 1
. For the advection example considered above, these
two corrections will sum to exactly the required correction (14), since this can be rewritten as
This correction is easily extended to an arbitrary hyperbolic system, since we presumably have a
Riemann solver that produces A \Gamma \Deltaq and A + \Deltaq from the two states -
m+1 and q 0
. We modify q 1
by
PSfrag replacements
PSfrag replacements

Figure

2: a) Waves arising from solving all Riemann problems on the fine and coarse grids sepa-
rately. b) The wave arising from solving the Riemann problem between the ghost cell value - q 0
and the coarse grid value q 0
must also be included in order to maintain conservation.
A similar modification must be made in each of the R time steps on the refined grid within the
single coarse grid step, where R is again the refinement ratio. The fix-up algorithm thus takes the
for
solve the Riemann problem with data -
m+1 and q 0
j to compute A \Gamma \Deltaq and A + \Deltaq,
update q 1
\Deltaq).
For the case of a conservation law, this will restore conservation and in fact agrees with the flux
function modification (4) if numerical fluxes are defined by (13). Note that in view of (12), in
the conservation law case the coarse grid value q 1
j is simply updated by the entire flux difference
need to actually solve the Riemann problem. However, implementing
it as presented above gives a uniform and general formulation.
The above idea extends directly to two space dimensions (and also to three dimensions). The
multi-dimensional wave-propagation algorithm consists of solving one-dimensional Riemann problems
normal to each cell interface. This defines waves and fluctuations exactly as in one dimension.
These are used as in the one-dimensional algorithm, and are also used to define "transverse correc-
tions" by essentially solving a Riemann problem in the transverse direction using the fluctuations
as data. These corrections are needed to give second order accuracy in multi-dimensions as well
as to increase the stability limit to allow Courant numbers close to 1. The corrections are fully
described in [14]. For our present purposes we need only note that these corrections modify the
correction fluxes ~
F (and corresponding y-fluxes ~
G) and are in flux-differencing form, so that they
are automatically corrected at grid interfaces in the step (4) where corrections are made due to
the fluxes. Numerical experiments demonstrate full second order accuracy for both conservation
laws and nonconservative equations, when the above approach is used (along with the second-order
correction terms discussed in Section 3).
5 Further extensions
The amrclaw software contains some further extensions of the wave-propagation algorithms discussed
in [14] and implemented in clawpack. These will be only briefly described here, with
emphasis on new issues that arise in connection with mesh refinement.
5.1 Source terms
Consider the hyperbolic equation
The source terms / can be handled using a standard fractional step method, also called a "splitting
method". In this approach, we alternate between solving the homogeneous hyperbolic equation,
ignoring the source term, and solving the ODE q In the context of AMR, we need to apply
this sequence within each time step on each grid. We must implement this carefully in order to
avoid generating excessive noise at the interface which can contaminate the solution. Whenever
we solve a Riemann problem we must ensure that the data on each side contains the same total
contribution from source terms in order to avoid the generation of spurious waves. This must be
carefully observed in computing ghost cell values for the fine grid at intermediate times, and also
in solving the additional Riemann problem need for the conservative fix described in Section 4.
Recall that this Riemann problem is between the ghost cell values at each intermediate time and
the coarse grid value at the original time, but when source terms are included this coarse grid value
must be modified to incorporate the correct source contributions.
These interface details will be presented below. First we discuss the basic fractional step method
that we use on each grid. In each time step we first solve the homogeneous hyperbolic equation over
time \Deltat and then use the resulting solution as initial data for the source-term equation over time
\Deltat. This is a so-called "first-order" splitting (or "Godunov splitting"), as opposed to the "Strang
splitting" in which one advances first by a half time step on one equation, then by a full time step
on the other equation, and ends with a half time step on the first equation again[17], [18]. Formally
the Strang splitting can give second order accuracy in situations where only first order accuracy is
achieved with our choice, but in practical problems where methods of this nature are useful, the
differences in resolution actually seen are generally negligible. This is because applying the Strang
splitting over N time steps is equivalent to starting with a half time step with one operator, then
alternating with \Deltat steps with each until the N'th step, where we finish with a half step of the first
operator again. The change made by this minor modification of shifting one half a time step from
the beginning to the end of the computation can be formally O(\Deltat), and hence reduce the global
accuracy to first order, but clearly will not degrade the overall resolution of the solution to any
degree. The solution may simply be shifted by O(\Deltat) relative to its correct location, for example.
This is very different from the sort of errors introduced by a genuine first order method, e.g., by
using a first order upwind method in place of the high-resolution flux-limiter method. See [14] for
more discussion of this point.
Using the simpler splitting has a number of advantages in the context of coupling source terms
with the hyperbolic solver, particularly with AMR. For one thing it is less expensive, since the
source terms are advanced only once per time step instead of twice. (Of course one could combine
the half time steps together in the Strang splitting as alluded to above, but this is impractical
when variable time steps are used, particularly with adaptive refinement.) The specification of
boundary conditions is also simplified. The user-supplied routine that extends values from the
computational grid to the ghost cells is called at the beginning of the time step, producing values
in the ghost cells that guarantee the physical boundary conditions will be satisfied. For example,
the normal momentum must be negated at a solid wall boundary for the Euler equations. Since
the hyperbolic equation is advanced first, these conditions are used immediately. Then the source
terms are advanced, typically a local ODE solve in each grid cell which does not require using ghost
cell values. We do not need to worry about solving the ODE's in ghost cells since these values are
replaced immediately at the start of the next time step.
Applying the hyperbolic solver first also simplifies the modifications needed at the grid interfaces
to minimize noise generation, and the algorithm is a fairly simple extension of what has been
presented already. Let H(k) represent the solution operator for the hyperbolic equation and S(k)
the solution operator for the source terms. Then the fractional step method over one time step on
a single grid takes the form
Now consider a refined grid, with refinement ratio R. As before, let -
q represent the solution on the
fine grid, where the time step is - k, and let q be the coarse grid solution with time step
k.
The algorithm is then:
Coarse grid update:
# Initialize coarse grid value needed for conservation fix-up:
Advance fine grids:
for
Space-time interpolate the ghost cell values -
m+2 using q 0 and q 1
Conservation fix-up:
Solve the Riemann problem with data -
m+1 and q 0;N
j to compute A \Gamma \Deltaq and A + \Deltaq
Apply source terms to coarse grid value:
This is basically a direct extension of the algorithm described previously except for the final
step, which ensures that the Riemann problem solved in the correction phase has consistent data
at each fine time step.
Numerical tests have shown that this approach performs quite well. Attempts with other styles
of splitting, e.g., the Strang splitting, or the Godunov splitting with the order of H and S reversed,
were much less successful. To gain a better appreciation of how mismatches at the interface can
generate noise, the reader is encouraged to explore various strategies on the simple scalar equation
with constant advection speed u and data q that is initially constant in space. If a numerical ODE
method is used to solve q than using the trivial exact solution operator), then any
approach will lead to the generation of noise at the interface simply because this ODE will be solved
more accurately on the fine grid than on the coarse grid, leading to jumps in q at the interface and
hence to waves propagating at speed u. The approach outlined above will only generate O(k 2 ) noise,
as long as the ODE solver is at least first order accurate, whereas any other approach considered
would generate O(k) noise.
5.2 Capacity-form differencing
The wave-propagation algorithm in [14] is also described in a more general form that applies to a
quasilinear equation of the form
in one dimension, with obvious generalization to more dimensions. Here the function -(x) is called
the "capacity function" since it represents in some way the capacity of the medium to hold the
quantity q. Working in this form is particularly useful for equations in the conservation form
where it is -(x)q(x; t) that is the conserved quantity rather than q alone, while the flux is defined
in terms of q. As a one-dimensional example consider flow through a variable-area duct where
-(x) is the cross-sectional area and q represents concentration per unit volume. In this case - i h
is the volume of the i'th grid cell. Flow in porous media is another example, where - represents
the porosity in one or more dimensions. This form also arises in using mapped curvilinear grids, in
which case - is the Jacobian of the grid transformation. More details are presented in [14].
The idea of capacity-form differencing is to replace the update formula (11) by
In the case of a conservation law, assuming (12) holds, this guarantees conservation of h
In
the definition of ~
must also incorporate - i into the second order corrections, replacing (10)
by
~
h- i
~
Using capacity-form differencing with AMR is direct, provided we insure that conservation is
maintained in transferring information between grids. Returning to Figure 1(b), suppose we have
capacities - i defined on the coarse grid and -
defined on the fine grid. In a region where the two
grids overlap, we assume that these values are consistent. For example, in Figure 1(b) we should
have
so that the capacity of this coarse grid cell agrees with the total capacity of the two fine cells:
h-
When updating the coarse grid value q 1
by the average of the fine grid values - q 2
we must weight by the capacity functions and replace (3) by
-m
Finally, in the grid-interface correction (15), we replace h in the denominator by - j h.
6 Numerical Results
We present AMR computations on some examples from [14] to illustrate that source terms, nonconservative
hyperbolic systems, capacity form differencing, and curvilinear grids can all be successfully
handled. Further examples can be found on the amrclaw webpage[6], including some animations.
Example 6.1. We repeat Example 3.8.3 from [14], which consists of the linear equations of
acoustics (a hyperbolic system of 3 equations) with a discontinuity in the sound speed across a line
oblique to the grid. A plane wave strikes the interface at some angle, leading to transmitted and
reflected waves. The time-evolution is best seen in Figure 4.
The acoustics equations for the pressure perturbation p and velocities u and v can be written
where
The coefficients are the density ae(x; y) and bulk modulus of elasticity K(x; y). In the example ae
has a discontinuity across the interface while K is constant. The Riemann solvers for this system
in the wave-propagation form are given in [14].

Figure

4(a) shows a contour plot of the initial pressure, a cosine hump as in [14] moving towards
the upper right. The dashed line shows the location of the discontinuity in sound speed. The
heterogeneous material is described by a density and bulk modulus of elasticity, and here the bulk
replacements

Figure

3: Pressure contours for a plane wave hitting a discontinuity in sound speed in the acoustics
equations, at time 0:6. The finer grids are shaded, with darker shading indicating finer levels.
Three levels are used with a 76 \Theta 50 grid at the coarsest level and refinement by a factor of 2 in
each level.
modulus is taken to be constant while the density is discontinuous, leading to the discontinuity in
sound speed. In each grid cell the density is defined as the cell average of the true density over that
cell. In [14] it is shown that the wave-propagation algorithm handles this problem well even when
the discontinuity in density is not aligned with the grid. Figure 3 shows an amrclaw calculation
where the coarsest grid is 76 \Theta 50 and two levels of refinement are used with
that the finest grid compares with the resolution seen in the lower plots of Figure 8 of [14]. The
tolerance used here was chosen in such a way that the weaker reflected wave has only been refined
to Level 2. A smaller tolerance would cause Level 3 refinement of this wave as well. Here we have
used a ratio of 2 for refinement in each level for demonstration purposes.
Example 6.2. The previous example does not fully test the new interface conditions between
the fine and coarse grids in the nonconservative case. These acoustics equations fail to be in
conservation form only along the interface where the density is discontinuous, and the wave stays
embedded in Level 3 grids as it moves along this interface. As a more severe test we repeated this
computation with a simple change in the error estimation procedure so that points are flagged for
refinement only if x ! 0:6. For x ? 0:6 there is only the coarsest grid, so the wave moves from the
initial fine grids onto the coarse grid as time advances. Figure 4 shows a sequence of times ending
with the time shown in Figure 3(b). Some smearing of the wave is seen on the coarser grid, which
is inevitable, but no difficulties are observed along the discontinuity in density.
Example 6.3. We repeat Example 3.10.1 from [14]. This is the advection equation
for a tracer q(x; t) in a density-stratified flow over a hump. We take u(x; y); v(x; y)) to be a fixed
(a) 1.2 1.40.10.30.50.70.9pressure at
(c) 1.2 1.40.10.30.50.70.9pressure at
1.2 1.40.10.30.50.70.9pressure at

Figure

4: Pressure contours for a plane wave hitting a discontinuity in sound speed in the acoustics
equations, at four different times. In this test refinement was allowed only for x ! 0:6 so that the
wave moves out of the refined region.
velocity field, chosen so that
while aeq is the conserved quantity.
In the test problem we consider flow over a hump with the bottom topography given by
in the domain 1). The velocity field is chosen by using the
"stream function"
to define
so that condition (21) is satisfied. Dividing by ae gives the velocity field. Note that (u; v) is not
divergence free and / is not a stream function for this velocity, though it is true that contours of
constant / give streamlines of the flow. We use the density profile
as in [14], so that the velocities increase exponentially with y.
In the test below we use 2:5. Figure 5 shows the initial data and
exact solution at time data consisting of a circular blob of tracer:
The problem is solved on a curvilinear grid. The irregular region of the x-y plane that can be
mapped smoothly to a rectangle. Then equation (20) can be transformed to an advection equation
on the rectangle and solved on a uniform Cartesian grid in this computational - \Gamma j space. Here
we use "Grid 2" of [14], with the mapping
Capacity-form differencing is used where being the Jacobian of the grid mapping,
as explained in [14]. Figure 6 shows computed results with 3 levels of refinement and refinement
ratio 2 in each case. The finest grid has the same resolution as the calculation shown in Figure 10

Figure

5: Density-stratified flow over a hump. Left: initial data is 1 inside the circular region and
elsewhere. The region where the solution is 1 at time Computed
results on a 200 \Theta 100 Cartesian grid. (From [14])
Conclusions
An adaptive mesh refinement algorithm developed for the Euler equations of gas dynamics has been
extended to employ high-resolution wave-propagation algorithms in a more general framework. In
particular, we have discussed the modifications needed to allow the application of this method to
hyperbolic problems which are not in conservation form, following the approach of [14]. This has
been done in a way that still maintains conservation when applied to a conservation law, in spite
of the fact that a wave-propagation approach is used to update cell averages rather than standard
flux-differencing. A generalization of this has also been presented for fractional step methods on
hyperbolic equations with source terms.
We have also discussed extensions to capacity-form differencing. This formulation is useful
in applying the methods on curvilinear grids, where refinement is done on logically rectangular
patches. This is illustrated in Section 6 for one example with the advection equations. In the
future we intend to explore the use of curvilinear grids more fully for nonlinear conservation laws
such as the Euler equations, since body-fitted curvilinear grids are crucial for some applications.
Adaptive mesh refinement is even more crucial in solving three-dimensional problems than in two
space dimensions. Recently clawpack has been extended to three dimensions [9], [10] and work is
underway to extend the amrclaw version. The underlying AMR code for the Euler equations has
previously been successfully extended to three dimensions (e.g., [1], [20]) and the ideas presented
in this paper carry over easily as well.
title
(b)
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.10.50.70.9PSfrag replacements

Figure

(a) Computed density at time stratified flow over a hump. Contour levels
are at 0:05; 0:95. Compare to Figure 10 in [14]. (b) A blow-up of the grids near the top
of the hump.

Acknowledgments

.
We are indebted to Smadar Karni for pointing out the difficulty of maintaining conservation at
grid interfaces with the wave-propagation framework. This led to the development of some of the
ideas presented here.
This work was supported in part by NSF Grants DMS-9505021, DMS-96226645, and DOE
Grants DE-FG03-96ER25292, DE-FG02-88ER25053, DE-FG02-92ER25139, and AFOSR grant
F49620-94-0132.



--R

Three dimensional adaptive mesh refinement for hyperbolic conservation laws
Adaptive mesh refinement for hyperbolic partial differential equations
Automatic adaptive grid refinement for the Euler equations
Adaptive mesh refinement for hyperbolic partial differential equations
Local adaptive mesh refinement for shock hydrodynamics
A test version is available at http://www.
An adaptively-refined Cartesian mesh solver for the Euler equa- tions



available from netlib.
Numerical Methods for Conservation Laws

Wave propagation algorithms for multi-dimensional hyperbolic systems
Solution of the Euler and Magnetohydrodynamic Equations on Solution-Adaptive Cartesian Grids

On the construction and comparison of

High resolution schemes using flux limiters for hyperbolic conservation laws
Some Aspects of the Computational Dynamics of Colliding Flows in Astrophysical Nebulae
--TR

--CTR
Lars Ferm , Per Ltstedt, Space---Time Adaptive Solution of First Order PDES, Journal of Scientific Computing, v.26 n.1, p.83-110, January   2006
K. W. Liang , P. Lin , R. C. E. Tan, Numerical solution of quenching problems using mesh-dependent variable temporal steps, Applied Numerical Mathematics, v.57 n.5-7, p.791-800, May, 2007
Zhijun Tan , Zhengru Zhang , Yunqing Huang , Tao Tang, Moving mesh methods with locally varying time steps, Journal of Computational Physics, v.200 n.1, p.347-367, 10 October 2004
Ding , C. Shu, A stencil adaptive algorithm for finite difference solution of incompressible viscous flows, Journal of Computational Physics, v.214 n.1, p.397-420, 1 May 2006
James A. Rossmanith, A wave propagation method for hyperbolic systems on the sphere, Journal of Computational Physics, v.213
R. Blikberg , T. Srevik, Load balancing and OpenMP implementation of nested parallelism, Parallel Computing, v.31 n.10-12, p.984-998, October - December 2005
Ragnhild Blikberg , Tor Srevik, Nested parallelism: Allocation of threads to tasks and OpenMP implementation, Scientific Programming, v.9 n.2,3, p.185-194, August 2001
Yuan , Tao Tang, Resolving the shock-induced combustion by an adaptive mesh redistribution method, Journal of Computational Physics, v.224 n.2, p.587-600, June, 2007
Siegfried Mller , Youssef Stiriba, Fully Adaptive Multiscale Schemes for Conservation Laws Employing Locally Varying Time Stepping, Journal of Scientific Computing, v.30 n.3, p.493-531, March     2007
Nail K. Yamaleev , Mark H. Carpenter, On accuracy of adaptive grid methods for captured shocks, Journal of Computational Physics, v.181 n.1, p.280-316, September 1 2002
Alexei Y. Poludnenko , Alexei M. Khokhlov, Computation of fluid flows in non-inertial contracting, expanding, and rotating reference frames, Journal of Computational Physics, v.220 n.2, p.678-711, January, 2007
software;finite-volume methods;hyperbolic conservation laws;gas dynamics;high resolution;godunov;adaptive mesh refinement;acoustics