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. For SDC coupling see the Spectral Deferred Corrections section.

The equations we integrate to do a nuclear burn are:

(1)\[\frac{dX_k}{dt} = \dot{\omega}_k(\rho,X_k,T)\]

(2)\[\frac{de}{dt} = \epsilon(\rho,X_k,T)\]

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

In this system, \(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 (see Neutrino Losses) in addition to the energy release from the changing binding energy of the fusion products.

Note

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

(3)\[\frac{dn_k}{dt} = \dot{\omega}_k(\rho,n_k,T)\]

(4)\[\frac{de}{dt} = \epsilon(\rho,n_k,T)\]

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 the system above 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 \(\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.

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.

`burner` interface

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

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().

Call the EOS directly on the input burn_t state using \(\rho\) and \(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.

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

Any reaction network must provide a righthand side and Jacobian function.

Important

Microphysics integrates the reaction system in terms of mass fractions, \(X_k\), but most astrophysical networks use molar fractions, \(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

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:

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, \(d({Y}_k)/dt\)
state.ydot(net_ienuc) : the change in the internal energy from the net, \(de/dt\)

Important

The righthand side routine is assumed to return the change in molar fractions, \(dY_k/dt\). These will be converted to the change in mass fractions, \(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, \(dX_k/dt\), then these will need to be converted to molar fraction rates for storage, e.g., \(dY_k/dt = A_k^{-1} dX_k/dt\).

Jacobian implementation

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:

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>

There are 4 different regions in the Jacobian: \(\partial \dot{\bf Y} / \partial {\bf Y}\), \(\partial \dot{\bf Y} / \partial {e}\), \(\partial \dot{e} / \partial {\bf Y}\), \(\partial \dot{e} / \partial {e}\). We discuss how these are computed and stored below:

\(\partial \dot{\bf Y} / \partial {\bf Y}\) :

This corresponds to elements \(d(\dot{Y}_m)/dY_n\)
- stored as: jac(m, n) for \(\mathrm{m}, \mathrm{n} \in [1, \mathrm{NumSpec}]\)
- computed as: each network has a function to compute these elements directly, since we need to know the stoichiometry. This is easy as it is just differentiating the righthand side with respect to species.
\(\partial \dot{\bf Y} / \partial {e}\) :

This corresponds to elements: \(d(\dot{Y}_m)/de\)
- stored as: jac(m, net_ienuc) for \(\mathrm{m} \in [1, \mathrm{NumSpec}]\)
- computed as: we directly compute the temperature derivative of the \(dY_m/dt\) expressions by computing the temperature derivative of the rates, i.e. \(d\lambda/dT\), and then evaluating the \(dY_m/dt\) using these temperature derivatives to get \(d{\dot{\bf Y}}/dT\).
  
  We then convert it to an energy derivative via the chain rule, namely:
  
  \[\frac{\partial\dot{\bf Y}}{\partial e} = \frac{1}{c_v} \frac{\partial \dot{\bf Y}}{\partial T}\]
\(\partial \dot{e} / \partial {\bf Y}\) :

This corresponds to \(d(\dot{e})/dY_n\)
- stored as: jac(net_ienuc, n) for \(\mathrm{n} \in [1, \mathrm{NumSpec}]\) :
- computed as: there are 3 different terms that make up the energy evolution:
  
  \[\frac{de}{dt} = \epsilon_\mathrm{nuc} - \epsilon_{\nu,\mathrm{weak}} - \epsilon_{\nu,\mathrm{therm}}\]
  
  The derivative of each of these (\(\partial \epsilon_* / \partial Y_n\)) are computed separately, and in different fashions:
  - \(\epsilon_\mathrm{nuc}\) : this is the energy release just from the change in mass:
    
    \[\epsilon_\mathrm{nuc} = -N_A \sum_{m=1}^{\mathrm{NumSpec}} \dot{Y}_m m_m c^2\]
    
    where \(m\) is the index of the nucleus, and \(m_m\) is the mass of that nucleus. Differentiating with respect to \(Y_n\), we have:
    
    \[\frac{\partial{\epsilon_\mathrm{nuc}}}{\partial Y_n} = -N_A \sum_{m=1}^{\mathrm{NumSpec}} \frac{\partial \dot{Y}_m}{\partial Y_n} m_m c^2\]
    
    We already have the \({\partial \dot{Y}_m}/{\partial Y_n}\), so the contribution of \(\epsilon_\mathrm{nuc}\) to each entry \(\partial (\dot{e})/\partial Y_n\) in the Jacobian is just the sum down column \(n\), weighting by \(mc^2\).
  - \(\epsilon_{\nu,\mathrm{weak}}\) : this represents the neutrino losses from weak rates. Presently this is not accounted for.
  - \(\epsilon_{\nu,\mathrm{therm}}\) : these represents the thermal neutrino losses (see Neutrino Losses). The neutrino loss functions directly provide \(\partial \epsilon_{\nu,\mathrm{therm}} / \partial \bar{A}\) and \(\partial \epsilon_{\nu,\mathrm{therm}} / \partial \bar{Z}\), so we can compute \(\partial \epsilon_{\nu,\mathrm{therm}} / \partial Y_n\) via the chain rule.
\(\partial \dot{e} / \partial {e}\) :
- stored as: jac(net_ienuc, net_ienuc)
- computed as: just like \(\partial \dot{e} / \partial Y_n\), there are 3 different terms that make up the energy evolution. The method for computing each contribution is largely the same:
  - \(\epsilon_\mathrm{nuc}\) : now we differentiate this with respect to temperature:
    
    \[\frac{\partial \epsilon_\mathrm{nuc}}{\partial T} = -N_A \sum_{m=1}^{\mathrm{NumSpec}} \frac{\partial \dot{Y}_m}{\partial T} m_m c^2\]
    
    and as before, we already have the \({\partial \dot{Y}_m}/{\partial T}\), so the contribution of \(\epsilon_\mathrm{nuc}\) to \(\partial (\dot{e})/\partial e\) is computed by summing down the last column of the Jacobian (weighting by \(mc^2\)) and adding the \(c_v\) weighting to convert from \(\partial/\partial T\) to \(\partial/\partial e\).
  - \(\epsilon_{\nu,\mathrm{weak}}\) : as above, we do not presently account for this.
  - \(\epsilon_{\nu,\mathrm{therm}}\) : as above, the neutrin loss functions directly provide \(\partial \epsilon_{\nu,\mathrm{therm}} / \partial T\).

Important

The Jacobian returned by the network is assumed to be in terms of molar fractions. However, we do convert the temperature derivative to an energy derivative already in the network by multiplying by \((1/c_v)\).

The form of the Jacobian return by the integrator looks like:

\[\begin{split}\left ( \begin{matrix} \dfrac{\partial \dot{Y}_1}{\partial Y_1} & \dfrac{\partial \dot{Y}_1}{\partial Y_2} & \cdots & \dfrac{\partial \dot{Y}_1}{\partial Y_\mathrm{NumSpec}} & \dfrac{1}{c_v} \dfrac{\partial \dot{Y}_1}{\partial T} \\ % \dfrac{\partial \dot{Y}_2}{\partial Y_1} & \dfrac{\partial \dot{Y}_2}{\partial Y_2} & \cdots & \dfrac{\partial \dot{Y}_2}{\partial Y_\mathrm{NumSpec}} & \dfrac{1}{c_v} \dfrac{\partial \dot{Y}_2}{\partial T} \\ % \vdots & \vdots & \ddots & \vdots & \vdots \\ % \dfrac{\partial \dot{Y}_\mathrm{NumSpec}}{\partial Y_1} & \dfrac{\partial \dot{Y}_\mathrm{NumSpec}}{\partial Y_2} & \cdots & \dfrac{\partial \dot{Y}_\mathrm{NumSpec}}{\partial Y_\mathrm{NumSpec}} & \dfrac{1}{c_v} \dfrac{\partial \dot{Y}_\mathrm{NumSpec}}{\partial T} \\ % \dfrac{\partial \dot{e}}{\partial Y_1} & \dfrac{\partial \dot{e}}{\partial Y_2} & \cdots & \dfrac{\partial \dot{e}}{\partial Y_\mathrm{NumSpec}} & \dfrac{1}{c_v} \dfrac{\partial \dot{e}}{\partial T} \end{matrix} \right )\end{split}\]

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.

Wrappers

To translate between the network’s righthand side and Jacobian functions and those expects by the ODE integrators, we provide a set of wrappers. These handle the conversion of variables (e.g. \(Y\) to \(X\)) and ensure that the state is thermodynamically consistent at the start.

These wrappers take an integrator state and copy back and forth to the burn_t that the networks want.

Note

In the flowcharts below, we’ll refer to the generic integrator type as int_state. The actual type will depend on the integrator used, e.g. dvode_t for VODE, rkc_t for RKC, …

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:

call clean_state on int_state
update the thermodynamics by calling update_thermodynamics. This takes both the int_state 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 \(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 wrapper

The integrator provides a wrapper that sits between the integration routines and the network’s implementation of the Jacobian. Its flow is:

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.

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.

Since \(Y_k = X_k/A_k\), we have \(\partial/\partial X_k = A_k^{-1} \partial/\partial Y_k\).

We transform by:
- multiplying all rows of the form \(\partial Y_m / \partial \star\) by \(A_m\) (where \(\star\) is either a molar fraction or \(T\)/\(e\)).
- multiplying all columns of the form \(\partial \star / \partial Y_n\) by \(1/A_n\).
add correction terms proportional to \(\partial e/\partial X_k |_{\rho, T, X_{j,j\ne k}}\) if integrator.correct_jacobian_for_const_e is 1.

The system we integrate is \((X_k, e)\), but the derivatives we took in the analytic Jacobian were in terms of \(T\) and not \(e\). So we need to correct for the fact that for some quantity \(q\),

\[\left . \frac{\partial q}{\partial X_k} \right |_e \ne \left . \frac{\partial q}{\partial X_k} \right |_T\]

If we write \(q = q(\rho, T(\rho, X_k, e), X_k)\) then we find that:

\[\left . \frac{\partial q}{\partial X_k} \right |_{\rho, e, X_{j,j\ne k}} = \left . \frac{\partial q}{\partial X_k} \right |_{\rho, T, X_{j,j\ne k}} - \frac{e_{X_k}}{c_v} \left . \frac{\partial T}{\partial X_k} \right |_{\rho, e, X_{j,j\ne k}}\]

where \(e_{X_k} = \partial e / \partial X_k |_{\rho, T, X_{j,j\ne k}}\).

This correction term is described in [43].
apply any boosting to the rates if integrator.react_boost > 0

The final form of the Jacobian is:

\[\begin{split}\left ( \begin{matrix} \dfrac{\partial \dot{X}_1}{\partial X_1} - \dfrac{e_{X_1}}{c_v} \dfrac{\partial \dot{X}_1}{\partial T} & \dfrac{\partial \dot{X}_1}{\partial X_2} - \dfrac{e_{X_2}}{c_v} \dfrac{\partial \dot{X}_1}{\partial T} & \cdots & \dfrac{1}{c_v} \dfrac{\partial \dot{X}_1}{\partial T} \\ % \dfrac{\partial \dot{X}_2}{\partial X_1} - \dfrac{e_{X_1}}{c_v} \dfrac{\partial \dot{X}_2}{\partial T} & \dfrac{\partial \dot{X}_2}{\partial X_2} - \dfrac{e_{X_2}}{c_v} \dfrac{\partial \dot{X}_2}{\partial T} & \cdots & \dfrac{1}{c_v} \dfrac{\partial \dot{X}_2}{\partial T} \\ % \vdots & \vdots & \ddots & \vdots \\ % \dfrac{\partial \dot{e}}{\partial X_1} - \dfrac{e_{X_1}}{c_v} \dfrac{\partial \dot{e}}{\partial T} & \dfrac{\partial \dot{e}}{\partial X_2} - \dfrac{e_{X_2}}{c_v} \dfrac{\partial \dot{e}}{\partial T} & \cdots & \dfrac{1}{c_v} \dfrac{\partial \dot{e}}{\partial T} \end{matrix} \right )\end{split}\]

Thermodynamics and \(e\) Evolution

The thermodynamic equation in our system is the evolution of the internal energy, \(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, \(e\) is set to the value from the EOS consistent with the initial temperature, density, and composition:

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

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

\[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 (2) requires an evaluation of the temperature at each integration step (since the RHS for the species is given in terms of \(T\), not \(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, \(c_v\), since we usually calculate derivatives with respect to temperature (as this is the form the rates are commonly provided in).

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.

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

Reaction ODE System

burner interface

Network routines

Righthand size implementation

Jacobian implementation

Wrappers

Righthand side wrapper

Jacobian wrapper

Thermodynamics and \(e\) Evolution

`burner` interface