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:
Here,
In this system,
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:
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
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
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.
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 and 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 advode_t
.Save the initial thermodynamic state for diagnostics and optionally subtracting off the initial energy later.
Call the ODE integrator,
dvode()
, passing in thedvode_t
and theburn_t
— as noted above, the auxiliary information that is not part of the integration state will be obtained from theburn_t
.Convert back to a
burn_t
by callingintegrator_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
Important
Microphysics integrates the reaction system in terms of mass
fractions,
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 ifNAUX_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 temperaturestate.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 theNumSpec
species that we are evolving,state.ydot(net_ienuc)
: the change in the internal energy from the net,
The righthand side routine is assumed to return the change in molar fractions,
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 thedvode_t
update the thermodynamics by calling
update_thermodynamics
. This takes both thedvode_t
and theburn_t
and computes the temperature that matches the current state.call
actual_rhs
convert the derivatives to mass-fraction-based (since we integrate
) 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
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>
The Jacobian matrix elements are stored in jac
as:
jac(m, n)
for :jac(net_ienuc, n)
for :jac(m, net_ienuc)
for :jac(net_ienuc, net_ienuc)
:
The form looks like:
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 theburn_t
call
actual_jac()
to have the network fill the Jacobian arrayconvert the derivative to be mass-fraction-based
apply any boosting to the rates if
integrator.react_boost
> 0
Thermodynamics and Evolution
The thermodynamic equation in our system is the evolution of the internal energy,
integrator.MAX_TEMP
(defaulting to 1.e11
) by clamping the
temperature if necessary.
At initialization,
As the system is integrated,
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 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
Note also that for the Jacobian, we need the specific heat,
Note
If desired, the EOS call can be skipped and the temperature and
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
.