appendix.tex

\appendix
\section{Interpolation methods}
\label{app:interpolation}

In this appendix, we provide details for the various interpolation
methods available in the code.  We assume throughout that we are
dealing with Cartesian coordinates and cells of equal sizes, although
we allow for an arbitrary refinement factor $r$ between cells at
different levels.

\vspace{0.5cm}\noindent
{\bf SecondOrderA} 

This interpolation algorithm is generally second-order, but has
monotonicity constraints as described below.
In one dimension, we define the parent values as $Q_{-1}$, $Q_0$, and
$Q_{+1}$, where the central parent cell has a left edge at $x_0$ and
width $\Delta x$. We first linearly interpolate the parent values to
the cell edges: $Q_{-1/2} = (Q_0 + Q_{-1})/2$, and similarly for
$Q_{1/2} = (Q_0 + Q_1)/2$, and then compute a monotonic slope: $\Delta
Q_0 = \minmod{Q_{1/2} - Q_0}{Q_0 - Q_{-1/2}}$ where

\begin{equation}
\minmod{a}{b} = \left\{ \begin{array}{ll}
0 & {\rm if} \quad ab < 0 \\
\min{ \left( | a |, | b | \right) } {\rm sign} (a) & {\rm otherwise}
\end{array}\right.
\end{equation}

This slope is then used to compute the interpolated subgrid values.
Defining $q_i$ to be the subgrid values at cell centers for the $r$
cells corresponding to parent cell $Q_0$ (for refinement factor $r$),
we can write:

\begin{equation}
q_i = Q_0 + \frac{i+(1-r)/2}{r} f_x
\end{equation}

where $f_x = 2 \Delta Q_0$ (this notation is used for consistency with
the 2 and 3 dimensional cases).

In two dimensions, the procedure is very similar in that we linearly
interpolate parent values to the four cell corners: $Q_{-1/2, -1/2}$,
$Q_{-1/2, 1/2}$, $Q_{1/2, -1/2}$, $Q_{1/2, 1/2}$ (by averaging parent
cell-centered parent values).  We then compute monotonic slopes across
the two diagonals, since for a linear function extrema occur at
corners:

\begin{eqnarray}
\Delta Q_0 & = & \minmod{Q_{1/2, 1/2} - Q_0}{Q_0 - Q_{-1/2, -1/2}} \\
\Delta Q_1 & = & \minmod{Q_{-1/2, 1/2} - Q_0}{Q_0 - Q_{1/2, -1/2}}
\end{eqnarray}

We then translate these slopes into the grid axes with $f_x = \Delta
Q_0 -  \Delta Q_1$ and $f_y = \Delta Q_0 + \Delta Q_1$ so that the
interpolation itself can be written simply as

\begin{equation}
q_{i,j} = Q_0 + \frac{i+(1-r)/2}{r} f_x + \frac{j+(1-r)/2}{r} f_y
\end{equation}

Finally, we write down the three-dimensional version -- unfortunately,
here the number of monotonicity constraints is four (the 4 diagonals
across the 8 opposing corners of the cube), while the number of slopes
is three, so the problem is over-constrained.  Somewhat arbitrarily,
we adopt the following procedure.  As before, we compute $\Delta Q_0$,
$\Delta Q_1$, $\Delta Q_2$, and $\Delta Q_3$ with the minmod limiter
across the 4 diagonals.  We then define

\begin{equation}
s =\frac{\Delta Q_1 + \Delta Q_2 + \Delta Q_3}{\Delta Q_0}
\end{equation}

which is the value of $\Delta Q_0$ that the other slopes ($\Delta
Q_1$, $\Delta Q_2$, and $\Delta Q_3$) imply, normalized by the desired
value of $\Delta Q_0$ itself.  If $0 < s < 1$, then no adjustment
needs to be made, as the monotonicity constraints are met.  If these
equalities are not met, then, if $s<0$ we define $\chi_n = 1$ if
$\Delta Q_n/\Delta Q_0 < 0$, and $\chi_n = 0$ otherwise (if $s > 1$, this is
reversed, so $\chi_n = 1$ if $\Delta Q_n/\Delta Q_0 > 0$, and 0
otherwise).  These weights are used to determine which slopes to
modify to match the $\Delta Q_0$ constraint.  We then compute the
amount of adjustment required:

\begin{equation}
f = - \frac{\Delta Q_0 s^\prime + \sum_n (1-\chi_n) \Delta Q_n}{\sum \chi_n \Delta Q_n + \epsilon}
\end{equation}

where $s^\prime = \min(\max(s,0), 1)$ and $\epsilon$ is a small number
to prevent numerical errors. We then use this adjustment fraction to
compute the new, adjusted $\Delta Q_n$ (for $n = 1, 2, 3$) that match
the $\Delta Q_0$ constraint as closely as possible with

\begin{equation}
\Delta Q_n = f \chi_n \Delta Q_n + (1-\chi_n) \Delta Q_n
\end{equation}

Finally, these are converted to the grid axes with $f_x = \Delta Q_2 +
\Delta Q_3$, $f_y = \Delta Q_1 + \Delta Q_3$, $f_z = \Delta Q_1 +
\Delta Q_2$.  We then do the interpolation itself with

\begin{equation}
q_{i,j,k} = Q_0 + \frac{i+(1-r)/2}{r} f_x + \frac{j+(1-r)/2}{r} f_y + \frac{k+(1-r)/2}{r} f_z
\end{equation}

% --------------

\vspace{0.3cm}\noindent
{\bf SecondOrderB} 

We also provide a variant on the above procedure, with two changes.
The first is that the slopes across the diagonals ($\Delta Q_0$, etc.)
are computed directly rather than with the minmod limiter
(e.g. $\Delta Q_0 = (Q_{1/2,1/2,1/2} - Q_{-1/2,-1/2,-1/2})/2$).  To
ensure positivity in the resulting interpolation when applied to
positive conserved quantities, the slopes are limited so that the
smallest corner value is 0.2 of the cell center value.  The procedure
described above is then applied to turn these four slopes into a
linear interpolation.

% --------------

\vspace{0.3cm}\noindent
{\bf SecondOrderC} 

This is a completely different second-order interpolation scheme,
which is based on Cloud-In-Cell (CIC) interpolation \citep{Hockney88}.
In one dimension, we define the parent values as $Q_0$, and $Q_{+1}$,
where the left parent cell has a cell center at $x_0$ and width
$\Delta x$.  Then, the interpolated value for a subgrid cell $q_i$
with a cell left edges at $x_i = x_0 + i \Delta x^p/r$, where $i$ runs
from 0 to $r-1$ for refinement factor $r$, is given by:

\begin{equation}
q_i =  \frac{2r - 1 - 2i}{2r} Q_0 + \frac{1+2i}{2r} Q_{+1}
\end{equation}

The extension to two and three dimensions is straightforward, with
weights in the other dimensions computed in the same way and then
multiplied to get a total of 4 and 8 weights for the 2 and 3
dimensional cases, respectively.  This scheme preserves monotonicity,
but is not conservative.

% --------------

\vspace{0.3cm}\noindent
{\bf ThirdOrderA} 

This interpolation method provides third-order accuracy based on the
Triangular Shaped Cloud (TSC) methodology \citep{Hockney88}.  As
usual, in one dimension, we define the parent values as $Q_{-1}$,
$Q_0$, and $Q_{+1}$, where the central parent cell has a left edge at
$x_0$ and width $\Delta x$.  Then, the interpolated value for a
subgrid cell $q_i$ with a cell left edges at $x_i = x_0 + i \Delta
x^p/r$ where $i$ runs from 0 to $r-1$ for refinement factor $r$, is
given by:

\begin{equation}
q_i = a_i  Q_{-1} + b_i Q_0 + c_i Q_{+1}
\end{equation}

and the weights are given by:

\begin{equation}
a(i) =  \frac{(r-i)^3 - (r-i-1)^3}{6r^3}; \qquad c(i) = \frac{3i^2 + 3i + 1}{6r^3}
\end{equation}

and $b(i) = 1/r - a(i) - c(i)$.  The extension to two and three
dimensions is straightforward with the weights computed in the same
way as for the one-dimensional case but then multiplied to determine
the 9 and 27 factors necessary for the 2 and 3 dimensional cases,
respectively.

One problem with this interpolation technique is that it is not
conservative.  In particular, the sum of the interpolated subgrid
values:

\begin{equation}
\tilde{Q_0} = \sum_{i=0}^{r-1} q_i
\end{equation}
is not, in general, equal to $Q_0$.  We can retain conservation by
adding the factor $(Q_0 - \tilde{Q_0})/r$ to the interpolated values
(or by multiplying the interpolated values by the ratio
$Q_0/\tilde{Q_0}$).  Unfortunately, the result of this procedure does
not preserve monotonicity -- it can introduce local minima and maxima
at parent cell boundaries.  We can then attempt to correct that by
taking a weighted average between the interpolated values on the
subgrid and the parent value, with weights computed such that the
interpolated values at the edge of the parent cell are not local
maxima compared to the interpolated values in the neighboring parent
cell.


% --------------

\vspace{0.3cm}\noindent
{\bf FirstOrderA}

Finally, for completeness, we include a first-order accurate piecewise
constant interpolator for which, using the same definitions as in the
previous case, we take $q_i = Q_0$, for $i = 0$ to $r-1$. 

\vspace{1cm}