# Local Source Terms

The `physics/sourceTerms`

organizational directory contains several
units that implement forcing terms. The `Burn`

, `Stir`

, `Ionize`

,
and `Diffuse`

units contain implementations in Flash-X. Two other
units, `Cool`

and `Heat`

, contain only stub level routines in their
API.

## Burn Unit

The nuclear burning implementation of the `Burn`

unit uses a
sparse-matrix semi-implicit ordinary differential equation (ODE) solver
to calculate the nuclear burning rate and to update the fluid variables
accordingly (Timmes 1999). The primary interface routines for this unit
are `physics/sourceTerms/Burn/Burn_init`

, which sets up the nuclear
isotope tables needed by the unit, and
`physics/sourceTerms/Burn/Burn`

, which calls the ODE solver and
updates the hydrodynamical variables in a single row of a single block.
The routine `physics/sourceTerms/Burn/Burn_computeDt`

may limit the
computational timestep because of burning considerations. There is also
a helper routine
`Simulation/SimulationComposition/Simulation_initSpecies`

(see
`Simulation/Simulation_initSpecies`

) which provides the properties of
ions included in the burning network.

### Algorithms

Modeling thermonuclear flashes typically requires the energy generation rate due to nuclear burning over a large range of temperatures, densities and compositions. The average energy generated or lost over a period of time is found by integrating a system of ordinary differential equations (the nuclear reaction network) for the abundances of important nuclei and the total energy release. In some contexts, such as supernova models, the abundances themselves are also of interest. In either case, the coefficients that appear in the equations are typically extremely sensitive to temperature. The resulting stiffness of the system of equations requires the use of an implicit time integration scheme.

A user can choose between two implicit integration methods and two
linear algebra packages in Flash-X. The runtime parameter
`Burn/odeStepper`

controls which integration method is used in the
simulation. The choice `odeStepper = 1`

is the default and invokes a
Bader-Deuflhard scheme. The choice `odeStepper = 2`

invokes a
Kaps-Rentrop or Rosenbrock scheme. The runtime parameter
`Burn/algebra`

controls which linear algebra package is used in the
simulation. The choice `algebra = 1`

is the default and invokes the
sparse matrix MA28 package. The choice `algebra = 2`

invokes the GIFT
linear algebra routines. While any combination of the integration
methods and linear algebra packages will produce correct answers, some
combinations may execute more efficiently than others for certain types
of simulations. No general rules have been found for best combination
for a given simulation. The most efficient combination depends on the
timestep being taken, the spatial resolution of the model, the values of
the local thermodynamic variables, and the composition. Users are
advised to experiment with the various combinations to determine the
best one for their simulation. However, an extensive analysis was
performed in the Timmes paper cited below.

Timmes (1999) reviewed several methods for solving stiff nuclear reaction networks, providing the basis for the reaction network solvers included with Flash-X. The scaling properties and behavior of three semi-implicit time integration algorithms (a traditional first-order accurate Euler method, a fourth-order accurate Kaps-Rentrop / Rosenbrock method, and a variable order Bader-Deuflhard method) and eight linear algebra packages (LAPACK, LUDCMP, LEQS, GIFT, MA28, UMFPACK, and Y12M) were investigated by running each of these 24 combinations on seven different nuclear reaction networks (hard-wired 13- and 19-isotope networks and soft-wired networks of 47, 76, 127, 200, and 489 isotopes). Timmes’ analysis suggested that the best balance of accuracy, overall efficiency, memory footprint, and ease-of-use was provided by the two integration methods (Bader-Deuflhard and Kaps-Rentrop) and the two linear algebra packages (MA28 and GIFT) that are provided with the Flash-X code.

### Reaction networks

We begin by describing the equations solved by the nuclear burning unit. We consider material that may be described by a density \(\rho\) and a single temperature \(T\) and contains a number of isotopes \(i\), each of which has \(Z_{i}\) protons and \(A_i\) nucleons (protons + neutrons). Let \(n_i\) and \(\rho_i\) denote the number and mass density, respectively, of the \(i\)th isotope, and let \(X_i\) denote its mass fraction, so that

where \(N_A\) is Avogadro’s number. Let the molar abundance of the \(i\)th isotope be

Mass conservation is then expressed by

At the end of each timestep, Flash-X checks that the stored abundances
satisfy [Eqn:mass conservation] to machine
precision in order to avoid the unphysical buildup (or decay) of the
abundances or energy generation rate. Roundoff errors in this equation
can lead to significant problems in some contexts (*e.g.*, classical
nova envelopes), where trace abundances are important.

The general continuity equation for the \(i\)th isotope is given in Lagrangian formulation by

In this equation \(\dot{R_i}\) is the total reaction rate due to all
binary reactions of the form *i(j,k)l*,

where \(\lambda _{kj}\) and \(\lambda _{jk}\) are the reverse (creation) and forward (destruction) nuclear reaction rates, respectively. Contributions from three-body reactions, such as the triple-\(\alpha\) reaction, are easy to append to [Eqn:binary rate]. The mass diffusion velocities \(\textbf{V}_i\) in [Eqn:isotope continuity] are obtained from the solution of a multicomponent diffusion equation (Chapman & Cowling 1970; Burgers 1969; Williams 1988) and reflect the fact that mass diffusion processes arise from pressure, temperature, and/or abundance gradients as well as from external gravitational or electrical forces.

The case \(\textbf{V}_i\equiv 0\) is important for two reasons.
First, mass diffusion is often unimportant when compared to other
transport processes, such as thermal or viscous diffusion (*i.e.*, large
Lewis numbers and/or small Prandtl numbers). Such a situation obtains,
for example, in the study of laminar flame fronts propagating through
the quiescent interior of a white dwarf. Second, this case permits the
decoupling of the reaction network solver from the hydrodynamical solver
through the use of operator splitting, greatly simplifying the
algorithm. This is the method used by the default Flash-X distribution.
Setting \(\textbf{V}_i\equiv 0\) transforms
[Eqn:isotope continuity] into

which may be written in the more compact, standard form

Stated another way, in the absence of mass diffusion or advection, any changes to the fluid composition are due to local processes.

Because of the highly nonlinear temperature dependence of the nuclear reaction rates and because the abundances themselves often range over several orders of magnitude in value, the values of the coefficients which appear in [Eqn:nucrate 1] and [Eqn:nucrate 2] can vary quite significantly. As a result, the nuclear reaction network equations are “stiff.” A system of equations is stiff when the ratio of the maximum to the minimum eigenvalue of the Jacobian matrix \({\tilde {\bf J}}\equiv\partial{{\bf f}}/\partial{{\bf y}}\) is large and imaginary. This means that at least one of the isotopic abundances changes on a much shorter timescale than another. Implicit or semi-implicit time integration methods are generally necessary to avoid following this short-timescale behavior, requiring the calculation of the Jacobian matrix.

It is instructive at this point to look at an example of how [Eqn:nucrate 1] and the associated Jacobian matrix are formed. Consider the \(^{12}\)C(\(\alpha\),\(\gamma\))\(^{16}\)O reaction, which competes with the triple-\(\alpha\) reaction during helium burning in stars. The rate \(R\) at which this reaction proceeds is critical for evolutionary models of massive stars, since it determines how much of the core is carbon and how much of the core is oxygen after the initial helium fuel is exhausted. This reaction sequence contributes to the right-hand side of [Eqn:nucrate 2] through the terms

where the ellipses indicate additional terms coming from other reaction sequences. The minus signs indicate that helium and carbon are being destroyed, while the plus sign indicates that oxygen is being created. Each of these three expressions contributes two terms to the Jacobian matrix \({\tilde {\bf J}}\)=\(\partial{{\bf f}}/\partial{{\bf y}}\)

Entries in the Jacobian matrix represent the flow, in number of nuclei per second, into (positive) or out of (negative) an isotope. All of the temperature and density dependence is included in the reaction rate \(R\). The Jacobian matrices that arise from nuclear reaction networks are neither positive-definite nor symmetric, since the forward and reverse reaction rates are generally not equal. In addition, the magnitudes of the matrix entries change as the abundances, temperature, or density change with time.

This release of Flash-X contains three reaction networks. A
seven-isotope alpha-chain (`Iso7`

) is useful for problems that do not
have enough memory to carry a larger set of isotopes. The 13-isotope
\(\alpha\)-chain plus heavy-ion reaction network (`Aprox13`

) is
suitable for most multi-dimensional simulations of stellar phenomena,
where having a reasonably accurate energy generation rate is of primary
concern. The 19-isotope reaction network (`Aprox19`

) has the same
\(\alpha\)-chain and heavy-ion reactions as the 13-isotope network,
but it includes additional isotopes to accommodate some types of
hydrogen burning (PP chains and steady-state CNO cycles), along with
some aspects of photo-disintegration into \(^{54}\)Fe. This 19
isotope reaction network is described in Weaver, Zimmerman, & Woosley
(1978).

The networks supplied with Flash-X are examples of a “hard-wired” reaction network, where each of the reaction sequences are carefully entered by hand. This approach is suitable for small networks, when minimizing the CPU time required to run the reaction network is a primary concern, although it suffers the disadvantage of inflexibility.

The MA28 sparse matrix package used by Flash-X is described by Duff,
Erisman, & Reid (1986) and is used as an external library. This package,
which has been described as the “Coke classic” of sparse linear algebra
packages, uses a direct – as opposed to an iterative – method for
solving linear systems. Direct methods typically divide the solution of
\(\tilde{{\bf A}} \cdot {\bf
x} = {\bf b}\) into a symbolic LU decomposition, a numerical LU
decomposition, and a backsubstitution phase. In the symbolic LU
decomposition phase, the pivot order of a matrix is determined, and a
sequence of decomposition operations that minimizes the amount of
fill-in is recorded. Fill-in refers to zero matrix elements which become
nonzero (*e.g.*, a sparse matrix times a sparse matrix is generally a
denser matrix). The matrix is not decomposed; only the steps to do so
are stored. Since the nonzero pattern of a chosen nuclear reaction
network does not change, the symbolic LU decomposition is a one-time
initialization cost for reaction networks. In the numerical LU
decomposition phase, a matrix with the same pivot order and nonzero
pattern as a previously factorized matrix is numerically decomposed into
its lower-upper form. This phase must be done only once for each set of
linear equations. In the backsubstitution phase, a set of linear
equations is solved with the factors calculated from a previous
numerical decomposition. The backsubstitution phase may be performed
with as many right-hand sides as needed, and not all of the right-hand
sides need to be known in advance.

MA28 uses a combination of nested dissection and frontal envelope decomposition to minimize fill-in during the factorization stage. An approximate degree update algorithm that is much faster (asymptotically and in practice) than computing the exact degrees is employed. One continuous real parameter sets the amount of searching done to locate the pivot element. When this parameter is set to zero, no searching is done and the diagonal element is the pivot, while when set to unity, partial pivoting is done. Since the matrices generated by reaction networks are usually diagonally dominant, the routine is set in Flash-X to use the diagonal as the pivot element. Several test cases showed that using partial pivoting did not make a significant accuracy difference but was less efficient, since a search for an appropriate pivot element had to be performed. MA28 accepts the nonzero entries of the matrix in the \((i, j, a_{i,j}\)) coordinate system and typically uses 70\(-\)90% less storage than storing the full dense matrix.

#### Two time integration methods

One of the time integration methods used by Flash-X for evolving the reaction networks is a 4th-order accurate Kaps-Rentrop, or Rosenbrock method. In essence, this method is an implicit Runge-Kutta algorithm. The reaction network is advanced over a timestep \(h\) according to

where the four vectors \(\Delta^i\) are found from successively solving the four matrix equations

\(b_i\), \(\gamma\), \(a_{ij}\), and \(c_{ij}\) are fixed constants of the method. An estimate of the accuracy of the integration step is made by comparing a third-order solution with a fourth-order solution, which is a significant improvement over the basic Euler method. The minimum cost of this method \(-\) which applies for a single timestep that meets or exceeds a specified integration accuracy \(-\) is one Jacobian evaluation, three evaluations of the right-hand side, one matrix decomposition, and four backsubstitutions. Note that the four matrix equations represent a staged set of linear equations (\(\Delta_4\) depends on \(\Delta_3 \ldots\) depends on \(\Delta_1\)). Not all of the right-hand sides are known in advance. This general feature of higher-order integration methods impacts the optimal choice of a linear algebra package. The fourth-order Kaps-Rentrop routine in Flash-X makes use of the routine GRK4T given by Kaps & Rentrop (1979).

Another time integration method used by Flash-X for evolving the
reaction networks is the variable order Bader-Deuflhard method (*e.g.*,
Bader & Deuflhard 1983). The reaction network is advanced over a large
timestep \(H\) from \({\bf y}^n\) to \({\bf y}^{n+1}\) by
the following sequence of matrix equations. First,

Then from \(k=1,2,\ldots,m-1\)

and closure is obtained by the last stage

This staged sequence of matrix equations is executed at least twice with \(m=2\) and \(m=6\), yielding a fifth-order method. The sequence may be executed a maximum of seven times, which yields a fifteenth-order method. The exact number of times the staged sequence is executed depends on the accuracy requirements (set to one part in 10\(^6\) in Flash-X) and the smoothness of the solution. Estimates of the accuracy of an integration step are made by comparing the solutions derived from different orders. The minimum cost of this method — which applies for a single timestep that met or exceeded the specified integration accuracy — is one Jacobian evaluation, eight evaluations of the right-hand side, two matrix decompositions, and ten backsubstitutions. This minimum cost can be increased at a rate of one decomposition (the expensive part) and \(m\) backsubstitutions (the inexpensive part) for every increase in the order \(2k+1\). The cost of increasing the order is compensated for, hopefully, by being able to take correspondingly larger (but accurate) timestep. The controls for order versus step size are a built-in part of the Bader-Deuflhard method. The cost per step of this integration method is at least twice as large as the cost per step of either a traditional first-order accurate Euler method or the fourth-order accurate Kaps-Rentrop discussed above. However, if the Bader-Deuflhard method can take accurate timesteps that are at least twice as large, then this method will be more efficient globally. Timmes (1999) shows that this is typically (but not always!) the case. Note that in , not all of the right-hand sides are known in advance, since the sequence of linear equations is staged. This staging feature of the integration method may make some matrix packages, such as MA28, a more efficient choice.

The Flash-X runtime parameter `Burn/odeStepper`

controls which
integration method is used in the simulation. The choice
`odeStepper = 1`

is the default and invokes the variable order
Bader-Deuflhard scheme. The choice `odeStepper = 2`

invokes the fourth
order Kaps-Rentrop / Rosenbrock scheme.

### Detecting shocks

For most astrophysical detonations, the shock structure is so thin that
there is insufficient time for burning to take place within the shock.
However, since numerical shock structures tend to be much wider than
their physical counterparts, it is possible for a significant amount of
burning to occur within the shock. Allowing this to happen can lead to
unphysical results. The burner unit includes a multidimensional shock
detection algorithm that can be used to prevent burning in shocks. If
the `Burn/useShockBurn`

parameter is set to `.false.`

, this
algorithm is used to detect shocks in the Burn unit and to switch off
the burning in shocked cells.

Currently, the shock detection algorithm supports Cartesian and 2-dimensional cylindrical coordinates. The basic algorithm is to compare the jump in pressure in the direction of compression (determined by looking at the velocity field) with a shock parameter (typically 1/3). If the total velocity divergence is negative and the relative pressure jump across the compression front is larger than the shock parameter, then a cell is considered to be within a shock.

This computation is done on a block by block basis. It is important that
the velocity and pressure variables have up-to-date guard cells, so a
guard cell call is done for the burners only if we are detecting shocks
(*i.e.* `useShockBurning = .false.`

).

### Energy generation rates and reaction rates

The instantaneous energy generation rate is given by the sum

Note that a nuclear reaction network does not need to be evolved in order to obtain the instantaneous energy generation rate, since only the right hand sides of the ordinary differential equations need to be evaluated. It is more appropriate in the Flash-X program to use the average nuclear energy generated over a timestep

In this case, the nuclear reaction network does need to be evolved. The energy generation rate, after subtraction of any neutrino losses, is returned to the Flash-X program for use with the operator splitting technique.

The tabulation of Caughlan & Fowler (1988) is used in Flash-X for most
of the key nuclear reaction rates. Modern values for some of the
reaction rates were taken from the reaction rate library of Hoffman
(2001, priv. comm.). A user can choose between two reaction rate
evaluations in Flash-X. The runtime parameter `Burn/useBurnTable`

controls which reaction rate evaluation method is used in the
simulation. The choice `useBurnTable = 0`

is the default and evaluates
the reaction rates from analytical expressions. The choice
`useBurnTable = 1`

evaluates the reactions rates from table
interpolation. The reaction rate tables are formed on-the-fly from the
analytical expressions. Tests on one-dimensional detonations and
hydrostatic burnings suggest that there are no major differences in the
abundance levels if tables are used instead of the analytic expressions;
we find less than 1% differences at the end of long timescale runs.
Table interpolation is about 10 times faster than evaluating the
analytic expressions, but the speedup to Flash-X is more modest, a few
percent at best, since reaction rate evaluation never dominates in a
real production run.

Finally, nuclear reaction rate screening effects as formulated by
Wallace *et al.* (1982) and decreases in the energy generation rate
\(\dot {\epsilon}_{\rm nuc}\) due to neutrino losses as given by
Itoh *et al.* (1996) are included in Flash-X.

### Temperature-based timestep limiting

When using explicit hydrodynamics methods, a timestep limiter must be used to ensure the stability of the numerical solution. The standard CFL limiter is always used when an explicit hydrodynamics unit is included in Flash-X. This constraint does not allow any information to travel more than one computational cell per timestep. When coupling burning with the hydrodynamics, the CFL timestep may be so large compared to the burning timescales that the nuclear energy release in a cell may exceed the existing internal energy in that cell. When this happens, the two operations (hydrodynamics and nuclear burning) become decoupled.

To limit the timestep when burning is performed, an additional
constraint is imposed. The limiter tries to force the energy generation
from burning to be smaller than the internal energy in a cell. The
runtime parameter `Burn/enucDtFactor`

controls this ratio. The
timestep limiter is calculated as

where \(E_{nuc}\) is the nuclear energy, expressed as energy per
volume per time, and \(E_{int}\) is the internal energy per volume.
For good coupling between the hydrodynamics and burning,
`enucDtFactor`

should be \(< 1\). The default value is kept
artificially high so that in most simulations the time limiting due to
burning is turned off. Care must be exercised in the use of this
routine.