*******************
Reaction ODE System
*******************

.. important::

   This describes the integration done when doing Strang operator-splitting, which is the
   default mode of coupling burning to application codes.

The equations we integrate to do a nuclear burn are:

.. math::
   \frac{dX_k}{dt} = \dot{\omega}_k(\rho,X_k,T)
   :label: eq:spec_integrate

.. math::
   \frac{de}{dt} = \epsilon(\rho,X_k,T)
   :label: eq:enuc_integrate

Here, :math:`X_k` is the mass fraction of species :math:`k`, :math:`e`
is the specific nuclear energy created through reactions. Also needed
are density :math:`\rho`, temperature :math:`T`, and the specific
heat. The function :math:`\epsilon` provides the energy release from
reactions and can often be expressed in terms of the instantaneous
reaction terms, :math:`\dot{X}_k`. As noted in the previous section,
this is implemented in a network-specific manner.

In this system, :math:`e` is equal to the total specific internal
energy. This allows us to easily call the EOS during the burn to obtain the temperature.

.. note::

   The energy generation rate includes a term for neutrino losses in addition
   to the energy release from the changing binding energy of the
   fusion products.

.. index:: integrator.use_number_densities

.. note::

   By setting ``integrator.use_number_densities=1``, number densities will be
   integrated instead of mass fractions.  This makes the system:

   .. math::
      \frac{dn_k}{dt} = \dot{\omega}_k(\rho,n_k,T)
      :label: eq:spec_n_integrate

   .. math::
      \frac{de}{dt} = \epsilon(\rho,n_k,T)
      :label: eq:enuc_n_integrate

   The effect of this flag in the integrators is that we don't worry
   about converting between mass and molar fractions when calling the
   righthand side function and Jacobian, and we don't do any normalization
   requiring $\sum_k X_k = 1$.


While this is the most common way to construct the set of
burn equations, and is used in most of our production networks,
all of them are ultimately implemented by the network itself, which
can choose to disable the evolution of any of these equations by
setting the RHS to zero. The integration software provides some
helper routines that construct common RHS evaluations, like the routine
that converts a temperature update to :math:`\dot{e}`, but these calls
are always explicitly done by the individual networks rather than
being handled by the integration backend. This allows you to write a
new network that defines the RHS in whatever way you like.

.. index:: integrator.react_boost

The standard reaction rates can all be boosted by a constant factor by
setting the ``integrator.react_boost`` runtime parameter.  This will simply
multiply the righthand sides of each species evolution equation (and
appropriate Jacobian terms) by the specified constant amount.

Interfaces
==========

The interfaces to all of the networks and integrators are written in C++.

``burner``
----------

The main entry point for C++ is ``burner()`` in
``interfaces/burner.H``.  This simply calls the ``integrator()``
routine (at the moment this can be ``VODE``, ``BackwardEuler``, ``ForwardEuler``, ``QSS``, or ``RKC``).

.. code-block:: c++

    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void burner (burn_t& state, amrex::Real dt)

The input is a ``burn_t``.

.. note::

   For the thermodynamic state, only the density, temperature, and
   mass fractions are used directly--we compute the internal energy
   corresponding to this input state through the equation of state
   before integrating.

When integrating the system, we often need auxiliary information to
close the system.  This is kept in the original ``burn_t`` that was
passed into the integration routines.  For this reason, we often need
to pass both the specific integrator's type (e.g. ``dvode_t``) and
``burn_t`` objects into the lower-level network routines.

Below we outline the overall flow of the integrator (using VODE as the
example).  Most of the setup and cleanup after calling the particular
integration routine is the same for all integrators, and is handled by
the functions ``integrator_setup()`` and ``integrator_cleanup()``.

.. index:: integrator.scale_system, burn_to_integrator, integrator_to_burn
.. index:: integrator.call_eos_in_rhs, integrator.subtract_internal_energy, integrator.burner_verbose

#. Call the EOS directly on the input ``burn_t`` state using
   :math:`\rho` and :math:`T` as inputs.

#. Scale the absolute energy tolerance if we are using
   ``integrator.scale_system``

#. Fill the integrator type by calling ``burn_to_integrator()`` to create a
   ``dvode_t``.

#. Save the initial thermodynamic state for diagnostics and optionally
   subtracting off the initial energy later.

#. Call the ODE integrator, ``dvode()``, passing in the ``dvode_t`` *and* the
   ``burn_t`` --- as noted above, the auxiliary information that is
   not part of the integration state will be obtained from the
   ``burn_t``.

#. Convert back to a ``burn_t`` by calling ``integrator_to_burn``

#. Recompute the temperature if we are using ``integrator.call_eos_in_rhs``.

#. If we set ``integrator.subtract_internal_energy``, then subtract
   off the energy offset, the energy stored is now just that generated
   by reactions.

#. Normalize the abundances so they sum to 1 (except if ``integrator.use_number_density`` is set).

#. Output statistics on the integration if we set ``integrator.burner_verbose``.
   This is not recommended for big simulations, as it will output information
   for every zone's burn.

.. index:: integrator.subtract_internal_energy

.. important::

   By default, upon exit, ``burn_t burn_state.e`` is the energy *released* during
   the burn, and not the actual internal energy of the state.

   Optionally, by setting ``integrator.subtract_internal_energy=0``
   the output will be the total internal energy, including that released
   burning the burn.

Network Routines
----------------

.. important::

   Microphysics integrates the reaction system in terms of mass
   fractions, :math:`X_k`, but most astrophysical networks use molar
   fractions, :math:`Y_k`.  As a result, we expect the networks to
   return the righthand side and Jacobians in terms of molar
   fractions.  The integration wrappers will internally
   convert to mass fractions as needed for the integrators.

Righthand size implementation
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The righthand side of the network is implemented by
``actual_rhs()`` in ``actual_rhs.H``, and appears as

.. code-block:: c++

   AMREX_GPU_HOST_DEVICE AMREX_INLINE
   void actual_rhs(burn_t& state, amrex::Array1D<amrex::Real, 1, neqs>& ydot)

All of the necessary integration data comes in through state, as:

* ``state.xn[NumSpec]`` : the mass fractions.

* ``state.aux[NumAux]`` : the auxiliary data (only available if ``NAUX_NET`` > 0)

* ``state.e`` : the current internal energy. It is very rare (never?) that a RHS
  implementation would need to use this variable directly -- even though this is
  the main thermodynamic integration variable, we obtain the temperature from the
  energy through an EOS evaluation.

* ``state.T`` : the current temperature

* ``state.rho`` : the current density

Note that we come in with the mass fractions, but the molar fractions can
be computed as:

.. code-block:: c++

      amrex::Array1D<amrex::Real, 1, NumSpec> y;
      ...
      for (int i = 1; i <= NumSpec; ++i) {
          y(i) = state.xn[i-1] * aion_inv[i-1];
      }

.. warning::

   We use 1-based indexing for ``ydot`` for legacy reasons, so watch out when filling in
   this array based on 0-indexed C arrays.

The ``actual_rhs()`` routine’s job is to fill the righthand side vector
for the ODE system, ``ydot(neqs)``. Here, the important
fields to fill are:

* ``state.ydot(1:NumSpec)`` : the change in *molar
  fractions* for the ``NumSpec`` species that we are evolving,
  :math:`d({Y}_k)/dt`

* ``state.ydot(net_ienuc)`` : the change in the internal energy
  from the net, :math:`de/dt`

The righthand side routine is assumed to return the change in *molar fractions*,
:math:`dY_k/dt`. These will be converted to the change in mass fractions, :math:`dX_k/dt`
by the wrappers that call the righthand side routine for the integrator.
If the network builds the RHS in terms of mass fractions directly, :math:`dX_k/dt`, then
these will need to be converted to molar fraction rates for storage, e.g.,
:math:`dY_k/dt = A_k^{-1} dX_k/dt`.

Righthand side wrapper
^^^^^^^^^^^^^^^^^^^^^^

The integrator provides a wrapper that sits between the integration
routines and the network's implementation of the righthand side.  Its
flow is (for VODE):

#. call ``clean_state`` on the ``dvode_t``

#. update the thermodynamics by calling ``update_thermodynamics``.  This takes both
   the ``dvode_t`` and the ``burn_t`` and computes the temperature that matches the
   current state.

#. call ``actual_rhs``

#. convert the derivatives to mass-fraction-based (since we integrate :math:`X`)
   and zero out the temperature and energy derivatives if we are not integrating
   those quantities.

#. apply any boosting if ``integrator.react_boost`` > 0


Jacobian implementation
^^^^^^^^^^^^^^^^^^^^^^^

.. index:: integrator.jacobian

Either an analytic or numerical Jacobian is used for the implicit
integrators, selected via the ``integrator.jacobian`` runtime
parameter (``1`` = analytic; ``2`` = numerical).  For VODE, the
numerical Jacobian is computed internally.  For the other integrators,
a difference method is implemented in
``integration/utils/numerical_jacobian.H``.

The analytic Jacobian is specific to each network and is provided by
``actual_jac(state, jac)``.  It takes the form:

.. code-block:: c++

   template<class MatrixType>
   AMREX_GPU_HOST_DEVICE AMREX_INLINE
   void actual_jac(const burn_t& state, MatrixType& jac)

where the ``MatrixType`` is most commonly ``MathArray2D<1, neqs, 1, neqs>``

The Jacobian matrix elements are stored in ``jac`` as:

* ``jac(m, n)`` for :math:`\mathrm{m}, \mathrm{n} \in [1, \mathrm{NumSpec}]` :
  :math:`d(\dot{Y}_m)/dY_n`

* ``jac(net_ienuc, n)`` for :math:`\mathrm{n} \in [1, \mathrm{NumSpec}]` :
  :math:`d(\dot{e})/dY_n`

* ``jac(m, net_ienuc)`` for :math:`\mathrm{m} \in [1, \mathrm{NumSpec}]` :
  :math:`d(\dot{Y}_m)/de`

* ``jac(net_ienuc, net_ienuc)`` :
  :math:`d(\dot{e})/de`

The form looks like:

.. math::
   \left (
   \begin{matrix}
      \ddots  & \vdots                          &          & \vdots \\
      \cdots  & \partial \dot{Y}_m/\partial Y_n & \cdots   & \partial \dot{Y}_m/\partial e    \\
              & \vdots                          & \ddots   & \vdots  \\
      \cdots  & \partial \dot{e}/\partial Y_n   & \cdots   & \partial \dot{e}/\partial e   \\
   \end{matrix}
   \right )

.. note::

   A network is not required to provide a Jacobian if a numerical
   Jacobian is used.

.. important::

   The integrator does not zero the Jacobian elements.  It is the responsibility
   of the Jacobian implementation to zero the Jacobian array if necessary.


Jacobian wrapper
^^^^^^^^^^^^^^^^

The integrator provides a wrapper that sits between the integration
routines and the network's implementation of the Jacobian.  Its
flow is (for VODE):

.. note::

   It is assumed that the thermodynamics are already correct when
   calling the Jacobian wrapper, likely because we just called the RHS
   wrapper above which did the ``clean_state`` and
   ``update_thermodynamics`` calls.

.. index:: integrator.react_boost

#. call ``integrator_to_burn()`` to update the ``burn_t``

#. call ``actual_jac()`` to have the network fill the Jacobian array

#. convert the derivative to be mass-fraction-based

#. apply any boosting to the rates if ``integrator.react_boost`` > 0





Thermodynamics and :math:`e` Evolution
======================================

The thermodynamic equation in our system is the evolution of the internal energy,
:math:`e`.  During the course of the integration, we ensure that the temperature stay
below the value ``integrator.MAX_TEMP`` (defaulting to ``1.e11``) by clamping the
temperature if necessary.

At initialization, :math:`e` is set to the value from the EOS consistent
with the initial temperature, density, and composition:

.. math::

   e_0 = e(\rho_0, T_0, {X_k}_0)

As the system is integrated, :math:`e` is updated to account for the
nuclear energy release (and thermal neutrino losses),

.. math:: e(t) = e_0 + \int_{t_0}^t \epsilon(\dot{Y}_k) dt

.. note::

   When the system is integrated in an operator-split approach, the
   energy equation accounts for only the nuclear energy release and
   not pdV work.

If ``integrator.subtract_internal_energy`` is set, then, on exit, we
subtract off this initial $e_0$, so ``state.e`` in the returned
``burn_t`` type from the ``actual_integrator`` call represents the
energy *release* during the burn.

Integration of Equation :eq:`eq:enuc_integrate` requires an evaluation
of the temperature at each integration step (since the RHS for the
species is given in terms of :math:`T`, not :math:`e`).  This involves
an EOS call and is the default behavior of the integration.

Note also that for the Jacobian, we need the specific heat,
:math:`c_v`, since we usually calculate derivatives with respect to
temperature (as this is the form the rates are commonly provided in).

.. index:: integrator.call_eos_in_rhs

.. note::

   If desired, the EOS call can be skipped and the temperature and
   $c_v$ kept frozen over the entire time interval of the integration
   by setting ``integrator.call_eos_in_rhs=0``.

.. index:: integrator.integrate_energy

We also provide the option to completely remove the energy equation from
the system by setting ``integrator.integrate_energy=0``.