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:

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

(2)\[\frac{de}{dt} = f(\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 \(f\) 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 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} = f(\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 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 \(\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.

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

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void burner (burn_t& state, 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.

The overall flow of the integrator is (using VODE as the example):

Call the EOS directly on the input burn_t state using \(\rho\) and \(T\) as inputs.
Fill the integrator type by calling burn_to_integrator() to create a dvode_t.
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.
subtract off the energy offset—we now store just the energy released in the dvode_t integration state.
convert back to a burn_t by calling integrator_to_burn
normalize the abundances so they sum to 1.

Note

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, \(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

void actual_rhs(burn_t& state, Array1D<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:

Array1D<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\)

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

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 \(X\)) and zero out the temperature and energy derivatives if we are not integrating those quantities.
apply any boosting if react_boost > 0

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:

void actual_jac(burn_t& state, MathArray2D<1, neqs, 1, neqs>& jac)

The Jacobian matrix elements are stored in jac as:

jac(m, n) for \(\mathrm{m}, \mathrm{n} \in [1, \mathrm{NumSpec}]\) : \(d(\dot{Y}_m)/dY_n\)
jac(net_ienuc, n) for \(\mathrm{n} \in [1, \mathrm{NumSpec}]\) : \(d(\dot{e})/dY_n\)
jac(m, net_ienuc) for \(\mathrm{m} \in [1, \mathrm{NumSpec}]\) : \(d(\dot{Y}_m)/de\)
jac(net_ienuc, net_ienuc) : \(d(\dot{e})/de\)

The form looks like:

\[\begin{split}\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 )\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.

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.

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 \(e\) Evolution

The thermodynamic equation in our system is the evolution of the internal energy, \(e\).

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.

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)\]

In the integration routines, this is termed the energy offset.

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

\[e(t) = e_0 + \int_{t_0}^t f(\dot{Y}_k) dt\]

As noted above, upon exit, we subtract off this initial offset, 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.

Reaction ODE System

Interfaces

burner

Network Routines

Righthand size implementation

Righthand side wrapper

Jacobian implementation

Jacobian wrapper

Thermodynamics and \(e\) Evolution

`burner`