numerical-analysis.tex

\section{Analysis}
\label{sec.num.analysis}

\subsection{Inline analysis with yt}

Detailed analysis of simulation results requires both the tools to ask
sophisticated questions of the data and the ability to process vast
quantities of data at high time-cadence.  As simulations grow in size
and complexity, storing data for post-processing simply becomes
intractable.  To cope with this, we have instrumented \enzo\ with the
ability to conduct analysis during the course of a simulation.  This
enables analysis with extremely high time cadence (as often as every
subcycle of the finest refinement level), without attempting to write
an entire checkpoint output to disk.  The current mechanism for
conducting analysis in \enzo\ during the course of the simulation
utilizes the same computional resources as are used by the simulation
itself by transferring their usage from \enzo\ to the analysis
routines; this is often referred to as \textit{in situ} analysis or
visualization.  Utilizing a dynamically-scheduled second set of
computing resources, often referred to as \textit{co-scheduled}
analysis or visualization, provides greater flexibility and overall
throughput at the expense of simplicity.

We expose \enzo's mesh geometry and fluid quantities to the analysis
platform \texttt{yt} \citep{2011ApJS..192....9T, 2011arXiv1112.4482T}.
At compile time, \enzo\ is (dynamically or statically) linked against
the Python and NumPy libraries necessary to create proxy objects
exposing the mesh geometry, fluid quantities and particle arrays as
NumPy arrays.  This information is then passed to a special handler
inside \texttt{yt}.  \texttt{yt} interprets the mesh and fluid
information and, without saving data to disk, constructs a native
representation of the in-memory state of the simulation that appears
identical to an on-disk simulation output.  \enzo\ then executes a
user-provided analysis script, which is able to access the in-memory
simulation object.  Once the analysis script has returned control to
\enzo, the simulation proceeds.  This process can occur at either the
top of the main ``EvolveHierarchy'' loop or at the end of a timestep
at the finest level, and the frequency with which it is called is
adjustable by a run-time parameter.  During the course of conducting
analysis, the simulation is halted until the conclusion of the
analysis.

Most analysis operations that can be performed on data sets that
reside on disk can be performed on in-memory data sets in \texttt{yt}.
This includes projections (i.e., line integrals, both on- and
off-axis), slices, 1-, 2-, and 3-D fluid phase distributions,
calculation of derived quantities and arbitrary data selection.  As of
version 2.5 of \texttt{yt}, the Rockstar phase-space halo finder
\citep{2013ApJ...762..109B} can be executed through \texttt{yt} on
in-memory \enzo\ data, and so can the Parallel HOP halo finder
\citep{1998ApJ...498..137E, 2010ApJS..191...43S}.  Operations that
currently cannot be conducted on in-memory datasets are those that
require spatial decomposition of data.  For instance, calculating
marching cubes on a data object with \texttt{yt} is a fully local
operation and can be conducted \textit{in situ}.  However, calculating
topologically-connected sets requires a spatial decomposition of data
and thus cannot be conducted \textit{in situ}.  This prohibition
extends to halo finding operations other than Rockstar, most
multi-level parallelism operations, and volume rendering.

Where microphysical solvers or other operator-split physics
calculations can be done in Python, \texttt{yt} can serve as a driver
for these calculations.  A major feature set that is currently being
developed is to pass structured (i.e., non-fluid) information back
from \texttt{yt} into \enzo.  For instance, this could be the result
of semi-analytic models of the growth and evolution of star clusters,
galaxy particle feedback parameters that have been influenced by
merger-tree analysis, or even spectral energy distributions that are
calculated within \texttt{yt} and provided to \enzo\ as input into
radiative transfer calculations.  Future versions will include this,
as well as the ability to dynamically allocate computional resources
to \texttt{yt} such that the simulation may proceed asynchronously
with analysis (\textit{co-scheduled} analysis).  With this
functionality will also come the ability to dynamically partition
data, such that spatially-decomposed operations such as volume
rendering become feasible during the course of a simulation.

\subsection{Tracer particles}

One of the inherent drawbacks of a grid-based fluid method is the
inability to follow the evolution of a single parcel of fluid as it
travels through the simulation volume.  To address this, \enzo\ has
the capability to introduce Lagrangian ``tracer particles'' into a
calculation either at the beginning of the simulation or when
restarting the calculation.  These tracer particles are put into the
simulation in a rectangular solid volume with uniform, user-specified
spacing.  Each particle's position and velocity is updated over the
course of a single timestep $\Delta t$ as follows:

\begin{eqnarray}
\label{eqn.drifttrace}
x^{n+1/2} & = & x^n + (\Delta t/2) v^{\rm interp,n} \nonumber \\
v^{n+1} & = & v^{\rm interp,n+1} \\
x^{n+1} & = & x^{n+1/2} + (\Delta t/2) v^{n+1} \nonumber
\end{eqnarray}

This is essentially a drift-kick-drift particle update from time $n$
to time $n+1$ -- however, instead of computing an acceleration at the
half-timestep $t^n + \Delta t/2$ (as is done for massive particles --
see Equation~\ref{eqn.driftkick}), the particle's velocity is updated
both at the beginning of each timestep and at the half-timestep by
linearly interpolating the cell-centered baryon velocity to the
position of the particle ($v^{\rm interp}$), and assigning it to that
value.

\enzo\ saves tracer particle data at user-specified intervals,
independent of the intervals at which regular data sets are written.
The data written out typically includes the tracer particle's unique
ID and position, as well as the velocity, density, and temperature of
the gas at its location, but is easily extensible to output any
grid-based quantity that a user requires.  This capability has been
used quite effectively in several papers, including
\citet{2010ApJ...715.1575S} and \citet{2012ApJ...748...12S}. It is
crucial to keep in mind that tracer particles model a fixed number of
Lagrangian fluid trajectories. The fluid on the grid, however, models
the motion of all the mass and represents the average quantity in a
grid cell's volume. Consequently, after a period of evolution, tracer
particles -- even if they initially had the same density distribution
as the gas -- will not have the same density distribution as the
fluid. For example, they tend to accumulate at stagnation points of
the flow, and care has to be taken in using these particles
appropriately. Tracer particles are very useful in studies such as the
variety of histories of the hydrodynamic quantities in Lagrangian
fluid elements and when evaluating complex chemical and cooling models
in regard to the simpler ones used in the actual numerical evolution.

\subsection{Shock finding}

Identification of shocks and their pre- and post-shock conditions can
be accomplished through a combination of either temperature or
velocity jumps with dimensionally split or unsplit search methods.
The primary method used in \enzo~is the dimensionally unsplit
temperature jump method described in detail in
\citet{2008ApJ...689.1063S}.  We briefly outline the method here.

For every cell, we first determine whether it satisfies the following
conditions necessary to be flagged as a shock:

\begin{eqnarray}
\div \myvec{v} < 0,\\
\grad T \cdot \grad S > 0,\\
T_2 > T_1,\\
\rho_2 > \rho_1,
\end{eqnarray}

where $\myvec{v}$ is the velocity field, $T$ is the temperature,
$\rho$ is the density, and $S=T/\rho^{\gamma-1}$ is the entropy.
$T_2$ and $T_1$ are the post-shock~(downstream) and
pre-shock~(upstream) temperatures, respectively. An optional
temperature floor may be chosen, which is useful for situations such
as cosmological simulations without a radiation background where
underdense gas in the intergalactic medium (IGM) can cool
adiabatically to unphysically low temperatures.

Once a cell is flagged, the local temperature gradient is calculated,
which is then used to traverse cells parallel to the gradient to
search for the first pre- and post-shock cell that do not satisfy the
above conditions.  If during the search a cell satisfying the
conditions is found to have a more convergent flow, that cell is
marked as the center, and the search is started again. Using the
temperature values from each of these cells, the Mach number is then
solved using the Rankine-Hugoniot temperature jump conditions:

\begin{equation}
\frac{T_2}{T_1} = \frac{(5\Mach^2 - 1)(\Mach^2 + 3)}{16\Mach^2},
\end{equation}

where $\Mach$ is the upstream Mach number.

Shock finding in the context of AMR is applied grid-by-grid.  If a
search for pre/post-shock cells goes outside the bounds of the grid
ghost zones, the search is stopped for that particular shocked
cell. In most situations this is adequate since the hydrodynamic shock
is captured in fewer than the number of ghost zones.  Shock finding
can be run either upon data output or at each step in the evolution of
an AMR level if the Mach and pre/post-shock quantities are needed for
additional physics modules.