Reactions

The nuclear network serves two purposes: it defines the fluid components used in both the equation of state and the hydrodynamics, and it evolves those components through a nuclear burning step. All of the reaction networks that Castro uses are provided by the Microphysics repository.

Note

To enable reactions in a simulation, you must compile with

USE_REACT=TRUE

This can be set in your GNUmakefile.

Microphysics comes with a general_null network. This is a bare interface for a nuclear reaction network. No reactions are enabled, and no auxiliary variables are accepted. It contains several sets of isotopes; for example, Microphysics/networks/general_null/triple_alpha_plus_o.net would describe the isotopes needed to represent the triple-\(\alpha\) reaction converting helium into carbon, as well as oxygen and iron.

Other reaction networks can be found in Microphysics/networks. To select a network, you set the make variable NETWORK_DIR to the name of the network directory.

Note

An arbitrary reaction network can be created for Castro via the pynucastro library.

The main interface for burning is in Microphysics/interfaces/burner.H:

void burner (burn_t& state, Real dt)

Here the burn_t type contains all of the information needed for the reaction network. It is similar to the equation of state eos_t.

Tip

The equation of state routines can be called directly with a burn_t in place of an eos_t.

Castro has several different modes of coupling reactions and hydrodynamics, selected through the parameter castro.time_integration_method. See Flowchart for the details.

Controlling burning

There are a number of reactions-related parameters that can be set at runtime in the inputs file. Reactions are enabled by setting:

castro.do_react = 1

(Note: turning reactions off for problems where they’re not required can help improve the efficiency).

It is possible to set the maximum and minimum temperature and density for allowing reactions to occur in a zone using the parameters:

castro.react_T_min and castro.react_T_max for temperature
castro.react_rho_min and castro.react_rho_max for density

Burning in Shocks

Burning can also be disabled inside shocks. This requires that the code be compiled with:

USE_SHOCK_VAR = TRUE

in the GNUmakefile. This will allocate storage for a shock flag in the conserved state array. This flag is computed via a multidimensional shock detection algorithm described in [17]. A zone is tagged as a shock if the following conditions are true:

\[\begin{split}\begin{align*} \nabla \cdot \ub &< 0 \\ \frac{|(\nabla p - \rho {\bf g}) \cdot \ub|}{p |\ub_\mathrm{cell}|} &> f_\mathrm{shock} \end{align*}\end{split}\]

This requires that there is compression and that the pressure jump (excluding the part of the pressure that balances gravity) is large. The runtime parameter

castro.disable_shock_burning = 1

will skip reactions in a zone where we’ve detected a shock. The runtime parameters castro.shock_detection_threshold and castro.shock_detection_include_sources will set the value of \(f_\mathrm{shock}\) and whether to subtract \(\rho {\bf g}\) from the pressure gradient.

Note

Both the compilation with USE_SHOCK_VAR = TRUE and the runtime parameter castro.disable_shock_burning = 1 are needed to turn off burning in shocks.

Reactions Flowchart

Here we describe how the burn_t is setup before the burn and how we update the castro state afterwards for both Strang and simplified-SDC.

Strang

In Castro_react.cpp, the flow is:

create burn_t burn_state
if NSE_NET is defined, initialize the chemical potentials that will be used as an initial guess for the NSE solve
- burn_state.mu_p \(= U(\mu_p)\)
- burn_state.mu_n \(= U(\mu_n)\)
- burn_state.y_e \(= 0\) (this will be filled if needed by the NSE routines)
initialize burn_state.dx – this is used for some NSE conditions.
set burn_state.success = true : we assume that the burn was successful. The integrator will set this to false is a problem occurred.
fill the thermodynamic quantities for input to the burner:
- burn_state.rho \(= U(\rho)\)
- burn_state.e \(= U(\rho e) / U(\rho)\)
- burn_state.T \(= U(T)\)
  
  Note
  
  It is assumed here that the temperature is thermodynamically consistent with the energy. For most networks, the temperature passed in will be used to set the thermodynamics in the burner.
- burn_state.xn[] \(= U(\rho X_k) / U(\rho)\)
- if NAUX_NET > 0: burn_state.aux[] \(= U(\rho \alpha_k) / U(\rho)\)
If we are doing castro.drive_initial_convection then we set burn_state.T_fixed by interpolating from the initial model.
Initialize the metadata that is used for diagnostics
Call the burner:
- We check to make sure that \(T\) and \(\rho\) are within the limits given by castro.react_T_min, castro.react_T_max, castro_react_rho_min, and castro.react_rho_max.
- The burner will set burn_state.success = false if it failed. This can happen for a number of reasons and is integrator-dependent.
  
  Note
  
  Castro will not abort by default here if the burn failed. Instead we leave it to the Retry Mechanism mechanism to attempt the step again with a smaller timestep.
Store the burning sources for plotting

We use the Reactions_Type StateData to hold the reactive sources that are output to the plotfile and the burn_weights MultiFab to hold the number of righthand side evaluations for diagnostics.

We fill these as:
- energy generation rate:
  
  \(\mathtt{reactions}(\rho e) = \dfrac{U(\rho) \, \cdot\, \mathtt{burn\_state.e}\, -\, U(\rho e)}{\Delta t}\)
- species and auxiliary creation rates (only if castro.store_omegadot = 1):
  - \(\mathtt{reactions}(\rho X_k) = U(\rho) \dfrac{\mathtt{burn\_state.xn[k]}\, -\, U(\rho X_k) / U(\rho)}{\Delta t}\)
  - \(\mathtt{reactions}(\rho \alpha_k) = U(\rho) \dfrac{\mathtt{burn\_state.aux[k]}\, -\, U(\rho \alpha_k) / U(\rho)}{\Delta t}\)
- NSE flag (only if NSE is defined). This simply stores the value of burn_state.nse.
Update the conserved state:

Note

\(\rho\) and \(\rho \ub\) are unchanged by reactions so those variables are not updated here. They are already the “new” state.
- \(U^\mathrm{new}(\rho e) = U^\mathrm{new}(\rho) \cdot \mathtt{burn\_state.e}\)
- \(U^\mathrm{new}(\rho E) = U^\mathrm{old}(\rho E) + (U^\mathrm{new}(\rho e) - U^\mathrm{old}(\rho e))\)
- \(U^\mathrm{new}(\rho X_k) = U^\mathrm{new}(\rho) \cdot \mathtt{burn\_state.xn[k]}\)
- if NAUX_NET > 0: \(U^\mathrm{new}(\rho \alpha_k) = U^\mathrm{new}(\rho) \cdot \mathtt{burn\_state.aux[k]}\)
- if NSE_NET :
  - \(U(\mu_p) = \mathtt{burn\_state.mu\_p}\)
  - \(U(\mu_n) = \mathtt{burn\_state.mu\_n}\)

Simplified-SDC

In Castro_react.cpp, the flow is:

create burn_t burn_state
if NSE_NET is defined, initialize the chemical potentials that will be used as an initial guess for the NSE solve
- burn_state.mu_p \(= U(\mu_p)\)
- burn_state.mu_n \(= U(\mu_n)\)
- burn_state.y_e \(= 0\) (this will be filled if needed by the NSE routines)
initialize burn_state.dx – this is used for some NSE conditions.
set burn_state.success = true : we assume that the burn was successful. The integrator will set this to false is a problem occurred.
fill the conserved state – this is stored in the burn_t only when we are using simplified-SDC.
- burn_state.y[SRHO] \(= U(\rho)\)
- burn_state.y[SMX] \(= U(\rho u)\)
- burn_state.y[SMY] \(= U(\rho v)\)
- burn_state.y[SMZ] \(= U(\rho w)\)
- burn_state.y[SEDEN] \(= U(\rho E)\)
- burn_state.y[SEINT] \(= U(\rho e)\)
- burn_state.y[SFS+k] \(= U(\rho X_k)\) for \(k = 0 \ldots N_{\mathrm{spec}} - 1\)
- if NAUX_NET > 0 : burn_state.y[SFX+k] \(= U(\rho \alpha_k)\) for \(k = 0 \ldots N_{\mathrm{aux}} - 1\)
fill the thermodynamic quantities in the burn_t :
- burn_state.rho \(= U(\rho)\)
- burn_state.T \(= U(T)\) – this is mainly going to be used as an initial guess
Note

We don’t initialize burn_state.xn[] or burn_state.aux[]
- if NAUX_NET > 0: burn_state.aux[] \(= U(\rho \alpha_k) / U(\rho)\)
If we are doing castro.drive_initial_convection then we set burn_state.T_fixed by interpolating from the initial model.
Store the advective update that will be used during the SDC integration.
Compute
Initialize the metadata that is used for diagnostics
Call the burner:
- We check to make sure that \(T\) and \(\rho\) are within the limits given by castro.react_T_min, castro.react_T_max, castro_react_rho_min, and castro.react_rho_max.
- The burner will set burn_state.success = false if it failed. This can happen for a number of reasons and is integrator-dependent.
  
  Note
  
  Castro will not abort by default here if the burn failed. Instead we leave it to the Retry Mechanism mechanism to attempt the step again with a smaller timestep.
Store the burning sources for plotting

We use the Reactions_Type StateData to hold the reactive sources that are output to the plotfile and the burn_weights MultiFab to hold the number of righthand side evaluations for diagnostics.

We fill these as:
- energy generation rate:
  
  \(\mathtt{reactions}(\rho e) = \dfrac{U(\rho) \, \cdot\, \mathtt{burn\_state.e}\, -\, U(\rho e)}{\Delta t}\)
- species and auxiliary creation rates (only if castro.store_omegadot = 1):
  - \(\mathtt{reactions}(\rho X_k) = U(\rho) \dfrac{\mathtt{burn\_state.xn[k]}\, -\, U(\rho X_k) / U(\rho)}{\Delta t}\)
  - \(\mathtt{reactions}(\rho \alpha_k) = U(\rho) \dfrac{\mathtt{burn\_state.aux[k]}\, -\, U(\rho \alpha_k) / U(\rho)}{\Delta t}\)
- NSE flag (only if NSE is defined). This simply stores the value of burn_state.nse.
Update the conserved state:

Note

\(\rho\) and \(\rho \ub\) are unchanged by reactions so those variables are not updated here. They are already the “new” state.
- \(U^\mathrm{new}(\rho e) = U^\mathrm{new}(\rho) \cdot \mathtt{burn\_state.e}\)
- \(U^\mathrm{new}(\rho E) = U^\mathrm{old}(\rho E) + (U^\mathrm{new}(\rho e) - U^\mathrm{old}(\rho e))\)
- \(U^\mathrm{new}(\rho X_k) = U^\mathrm{new}(\rho) \cdot \mathtt{burn\_state.xn[k]}\)
- if NAUX_NET > 0: \(U^\mathrm{new}(\rho \alpha_k) = U^\mathrm{new}(\rho) \cdot \mathtt{burn\_state.aux[k]}\)
- if NSE_NET :
  - \(U(\mu_p) = \mathtt{burn\_state.mu\_p}\)
  - \(U(\mu_n) = \mathtt{burn\_state.mu\_n}\)