methods-implementation.tex

\section{Implementation} \label{methods:code}

In this section, we discuss details of the code itself, including the
available application program interfaces (APIs), code layout, and how
the existing models and networks may be extended.  For more detailed
information, users should consult online documentation available at
\texttt{https://grackle.readthedocs.org}.  In addition to the
documentation, the \texttt{Grackle} community maintains a mailing list
where users can post questions or comments and receive help from other
users and developers.  More information can be found on the front page
of the documentation.

\subsection{Simulation Code API} \label{methods:api}

The \texttt{Grackle} library provides five main functions to the user
for use in simulation codes through C or Fortran bindings: solving the
chemistry and cooling (i.e., updating the chemical species and
internal energy) and calculating the cooling time, temperature,
pressure, and ratio of specific heats.  Before these functions can be
called, the code must be initialized with various user-specified
settings.  This initialization process is also responsible for loading
data from external files and calculating the chemistry and cooling
rate tables used by the solvers.  All \texttt{Grackle} run-time
parameters are stored within a C \texttt{struct} of type
\texttt{chemistry\_data}.  The user initializes this structure by
calling the function \texttt{set\_default\_chemistry\_parameters} and
supplying a pointer to a \texttt{chemistry\_data} structure.  The
\texttt{chemistry\_data} pointer is then attached to a globally
viewable pointer called \texttt{grackle\_data}, allowing all run-time
parameters to be accessible without having to store the struct
manually.  Once all parameters are properly set, the user must call
\texttt{initialize\_chemistry\_data} to finalize the initialization
process.  An example of this procedure is shown below:

\vspace{0.5cm}
\begin{minipage}[b]{0.5\linewidth}
\begin{verbatim}

int rval;
chemistry_data my_pars;
rval =
  set_default_chemistry_parameters(&my_pars);

grackle_data.use_grackle = 1;
grackle_data.with_radiative_cooling = 1;
grackle_data.primordial_chemistry = 3;
grackle_data.metal_cooling = 1;
grackle_data.UVbackground = 1;
grackle_data.grackle_data_file = 
  "CloudyData_UVB=HM2012.h5";

rval = initialize_chemistry_data(&my_units);

\end{verbatim}
\end{minipage}

In the above example, the variable \texttt{my\_units} is a C
\texttt{struct} that holds unit conversions from internal code units
to the CGS unit system for quantities such as density, length, and
time.  These are required in order to set up the internal unit system
for the chemistry and cooling rate tables.

Once \texttt{Grackle} has been initialized, the functionality described
above can be called by the simulation code.  The available functions
are: 1) \texttt{solve\_chemistry} to integrate chemistry and cooling
equations over a specified time step, 2)
\texttt{calculate\_cooling\_time} to calculate the cooling time
($e/(de/dt)$) for each computational element, 3)
\texttt{calculate\_temperature} to calculate the temperature from the
internal energy and chemical species densities, 4)
\texttt{calculate\_pressure} to calculate the gas pressure, and 5)
\texttt{calculate\_gamma} to calculate the ratio of specific heats.
In \texttt{Grackle}, the ratio of specific heats is only altered from that
of a monatomic gas by the presence of H$_{2}$.

For efficiency, \texttt{Grackle}'s functions are designed to operate on
multiple computational elements simultaneously.  The user provides
arrays of the required fields to \texttt{Grackle} and their values are
updated by the chemistry and cooling solvers.  Because the number of
required fields depends on the specific solver being used, \texttt{Grackle}
makes use of another C \texttt{struct} as a means of passing field
arrays to the \texttt{Grackle} functions.  The \texttt{grackle\_field\_data}
\texttt{struct} contains pointers to which can be attached the arrays
of density, internal energy, and any optional fields, such as
individual species densities or the arrays of constant heating rates
(Section \ref{section:constant-heating}).  Since arrays of the
optional fields are only accessed based on run-time parameter
settings, the user has the option of only providing the fields they
wish to use.  The field arrays can be one, two, or three-dimensional,
allowing both Lagrangian particle-based codes and Eulerian mesh-based
codes to provide fields in their native layout.  The
\texttt{grackle\_field\_data} \texttt{struct} contains entries to
specify the field dimensionality as well as to flag certain array
elements as boundary cells which are to be ignored.  An example of
calling the main chemistry solver function on a 10$^{3}$ grid with 3
boundary zones on each side is shown below:

\vspace{0.5cm}
\begin{minipage}[b]{0.5\linewidth}
\begin{verbatim}

grackle_field_data my_fields;

my_fields.grid_rank = 3;
my_fields.grid_dimension = new int[3];
my_fields.grid_start = new int[3];
my_fields.grid_end = new int[3];
for (int i = 0;i < 3;i++) {
  my_fields.grid_dimension[i] = 10;
  my_fields.grid_start[i] = 3;
  my_fields.grid_end[i] = 12;
}

my_fields.density = density_array;
my_fields.internal_energy = energy_array;
my_fields.HI_density = HI_array;
...

// 1 Myr in internal units
double dt = 3.15e7 * 1e6 /
  my_units.time_units;

rval =
  solve_chemistry(&my_units, &my_fields, dt);

\end{verbatim}
\end{minipage}

An added benefit of this approach is that adding new
features which use additional fields will not require a change in the
function signature.  This will, in theory, allow \texttt{Grackle} to maintain
backward compatibility indefinitely.  We note that versions of \texttt{Grackle}
prior to 3.0 did not make use of the \texttt{grackle\_field\_data}
\texttt{struct} and instead required all fields to be provided as
individual arguments.  We acknowledge that the release of
\texttt{Grackle} 3.0 constitutes a significant change to the API, but
one that will ultimately provide more stability moving forward.

\subsection{Pygrackle - Grackle's Python API}

As described above, \texttt{Grackle}'s native API is provided through C and Fortran
bindings.  This is particularly useful for simulation codes, as they are
typically written in either C or Fortran compatible languages, but it provides
a barrier to entry for experimental work and testing of \texttt{Grackle} functionality.
Python is a high-level, interpreted language, increasingly used in scientific
computing both as ``glue" code as well as a mechanism for authoring production
scientific codes.  For instance, in 2016, one of the Gordon Bell Prize nominees
utilized the Python code PyFR to demonstrate extreme scaling of finite element
calculations \citep{pyfr}.

To facilitate access to \texttt{Grackle} functionality, we provide Python bindings to
it, designed to make interacting with chemical rates and rate equations
available to researchers at all stages.  This enables rapid iteration over
different rate coefficients, different initial state vectors, and over
arbitrary time periods.  The bindings are written in Cython
\citep{behnel2010cython} with minimal overhead from Python operations.  Below,
we describe two particular aspects of \texttt{Pygrackle}.

\subsubsection{Fluid Container} \label{sec:pyfluid}

\texttt{Pygrackle} provides an object called a ``fluid container''.  When data is passed
to \texttt{Grackle} during the course of a simulation, it may be sent as a single zone,
as multiple zones that are organized in a 3D array, or as a 1D ``pencil beam'' of
data.  The fluid container object is designed to mimic this, enabling
individuals to ``create" a collection of fluids either by reading them from
disk or constructing them in-memory from NumPy arrays \citep{numpy}.  This
enables calculations of cooling time, chemical evolution, etc., from
analytically-defined gas parcels and distributions.  While this will not take
into account hydrodynamics (unless a package spiritually similar to \texttt{Grackle},
but for hydrodynamics, is released) it can provide useful and scientifically
relevant information.

The fluid container provides several high-level functions, such as computing
the cooling time, the pressure, and so forth, most of which are usually used
only internally in \texttt{Grackle}.  Additionally, this object provides compatibility
with \texttt{yt} \citep{2011ApJS..192....9T}, enabling individuals to load data
with \texttt{yt} and then evolve it through \texttt{Grackle}.  One possible use case for
this is to load in a dataset and use \texttt{Pygrackle} to compute different cooling
times for a collection of gas identified in \texttt{yt} based on different
metallicity assumptions, different chemical rate coefficients, and different
radiation backgrounds.

\subsubsection{Evolution Models} \label{sec:pyevolve}

In addition to providing access to \texttt{Grackle}'s primary
functionality, \texttt{Pygrackle} also features a set of convenience functions
to evolve a fluid container forward in time following simple models.
These functions return a dictionary of arrays of all fluid quantities
(i.e., species densities, internal energy, temperature, etc.) for time
values of the evolution.  These functions can be used in semi-analytic
models that require knowledge of the thermal evolution of gas under
different conditions.  Examples of use are discussed in
Section \ref{sec:testing}.

The simplest of these evolves a fluid container assuming a constant
density model.  The \texttt{evolve\_constant\_density} function takes
an initialized fluid container and repeatedly calls \texttt{Grackle}'s
\texttt{solve\_chemistry} function until a specified stopping time or
temperature has been reached.  For each iteration, the timestep is
taken to be a fraction of the local cooling time, specifiable by the
user and defaulting to 0.01.

The second of these functions models the evolution of a parcel of gas
collapsing due to self-gravity.  The \texttt{evolve\_freefall} function
closely follows the one-zone free-fall collapse model introduced by
\citet{2000ApJ...534..809O} and modified by
\citet{2005ApJ...626..627O} to include the effects of thermal pressure
support.  The gas density evolves following the modified collapse
model of \citet{2005ApJ...626..627O}, given by
\begin{equation}
\frac{d \rho}{dt} = \frac{\rho}{t_{\rm col}},
\end{equation}
where $t_{\rm col}$ is the collapse time-scale expressed as
\begin{equation}
t_{\rm col} = \frac{t_{\rm ff}}{\sqrt{1-f}},
\end{equation}
and the free-fall time is given by
\begin{equation}
t_{\rm ff} = \sqrt{\frac{3 \pi}{32 G \rho}}.
\end{equation}
Thermal pressure forces, which act to slow the collapse of the cloud,
are modeled by the factor $f$, which is expressed as
\begin{equation}
f = \left\{
\begin{array}{lr}
0, & \gamma < 0.83,\\
0.6 + 2.5(\gamma - 1) - 6.0(\gamma - 1)^2, & 0.83 < \gamma < 1,\\
1.0 + 0.2(\gamma - 4/3) - 2.9(\gamma - 4/3)^2, & \gamma > 1,
\end{array}
 \right.
\end{equation}
where the effective adiabatic index, $\gamma$, is
\begin{equation}
\gamma \equiv \frac{\partial \ln p}{\partial \ln \rho}.
\end{equation}
As the density increases, the internal energy is altered by a
combination of adiabatic compression and radiative cooling, computed
by \texttt{Grackle}, and is given by
\begin{equation}
\frac{de}{dt}= -p \frac{d}{dt} \frac{1}{\rho} - {\Lambda},
\label{eq:energy}
\end{equation}
where the pressure is given by
\begin{equation}
p = \frac{\rho k T}{\mu m_{\rm H}},
\end{equation}
the specific internal energy is given by Equation \ref{eqn:e-T},
and $\Lambda$ is the
radiative cooling rate.  For each iteration of this function, the
timestep is taken to be a fraction of the local free-fall time,
specifiable by the user and defaulting to 0.01.  An example use of
this function is shown in Section \ref{sec:free-fall-test}.

%\subsection{Code Structure} \label{Code_Structure}

\subsection{Adding New Models/Rates}

Adding new rates to \texttt{Grackle} involves modifying the values stored in
the rate coefficients and additionally defining a new network for both
the chemistry and if desired also the cooling. As it stands, some
direct modification of the code structure is required (although see \S
\ref{Future_Directions} for future improvements in this area), and we
give more details of that in this section.

The structure of the current code base is set up to realize either an
equilibrium chemistry model using cooling tables derived from \texttt{Cloudy},
or a non-equilibrium cooling model based on a six, nine or twelve
species chemistry model.  To implement a different equilibrium cooling
model, the process is relatively straightforward: \texttt{Grackle} reads
the \texttt{Cloudy} cooling and interpolates the data along one, two or three
dimensions as appropriate. By substituting in a different cooling
table, with the correct format, the behavior of the cooling for a
given species can be easily changed.

To modify the non-equilibrium cooling behavior, more work is
required. \texttt{Grackle} assumes that non-equilibrium chemistry models are
hierarchical, with the nine species model composed of the six species
model plus three additional species and likewise for the twelve species model.
To expand the network to include, say, a fourth model containing the twelve species 
model as a subset then the procedure is straightforward. However, if a different
chemical network is envisaged containing some species already in the
twelve species model and some not, then this will require modifying
the existing hierarchical structure.  The next sections provide some
details of this process.

\subsubsection{Updating the rate coefficients}

The rates are allocated and declared in the
\texttt{initialize\_chemistry} function. To update the rate
coefficients (e.g., based on more recent experimental data), simply
make use of the already existing rate coefficient arrays (\texttt{k}
array) as declared in \texttt{initialize\_chemistry}. The rates are
populated in \texttt{calc\_rates\_g}. Within this function the rates
are fully described. For example, \texttt{k1} is the rate coefficient
for the collisional ionization of neutral hydrogen by electrons
(\texttt{k1}: H + e$^{-}$ $\rightarrow$ H$^{+}$ + 2e$^{-}$). If a
significant number of new rate coefficients are needed then the most
expedient approach would be to insert a preprocessor directive into
\texttt{calc\_rates\_g} in which the appropriate function call can be
inserted.  For any additional rate coefficients that are required, the
corresponding \texttt{k} array needs to be declared and allocated in
\texttt{initialize\_chemistry}.  Furthermore, the interpolation of the
new rates will need to be added in the function
\texttt{lookup\_cool\_rates1d}.


\subsubsection{Updating the chemistry network}

To update the chemical network to either include or exclude reactions,
a new rate network will be required. As a template, the function
\texttt{step\_rate} can be used.  If the new network is simply an
addition to the existing network (e.g., a 15 species model) then the
easiest option is to simply augment this network with the three extra
species using the appropriate interactions. The species will then be
evolved until they converge. More complex additions would require
creating a new network update routine using \texttt{step\_rate}
as a template.

\subsubsection{Updating the cooling model}

Apart from the chemical network, the cooling model may also be
modified. As discussed previously, if the intention is to implement a
new cooling table then the changes are straightforward.  For changes
to the equilibrium network (say, to modify the cooling rate due to, for
example, emission line cooling from a given species), this is handled
in \texttt{cool1d\_multi\_g}. As discussed in \S \ref{Cooling}, line
emission cooling is determined using collisional excitation,
collisional ionization and recombination rates. If the intention is to
update/modify existing rates then the cooling rates are also set in
\texttt{calc\_rates\_g} and can be modified there. If new cooling
rates are required from another species whose cooling properties are
important for the network then the rates can be added there also. Any
new arrays required in this case will also need to be declared and
initialized in \texttt{initialize\_chemistry} in a way similar to the
rate coefficients. The rates can then be interpolated to the required
temperature values in \texttt{calc\_rates\_g}, following the examples
there. Finally, once the new cooling rate has been determined its
values needs to be added to the \texttt{edot} array, which tracks
the total cooling/heating rate.