18.1 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 Burn_init, which sets up the nuclear isotope tables needed by the unit, and Burn, which calls the ODE solver and updates the hydrodynamical variables in a single row of a single block. The routine Burn_computeDt may limit the computational timestep because of burning considerations. There is also a helper routine Simulation/SimulationComposition/Simulation_initSpecies (see Simulation_initSpecies) which provides the properties of ions included in the burning network.

18.1.1 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. The runtime parameter 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 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. 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 code.

18.1.2 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

$$X_i = \rho_i/\rho = n_i A_i/(\rho N_A) ,$$ (18.1)

where $N_A$ is Avogadro's number. Let the molar abundance of the $i$th isotope be

$$Y_i = X_i/A_i = n_i/(\rho N_A) .$$ (18.2)

Mass conservation is then expressed by

$$\sum_{i=1}^N X_i = 1 .$$ (18.3)

At the end of each timestep, FLASH checks that the stored abundances satisfy (18.3) 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

$${dY_i\over dt} + \nabla \cdot \left ( Y_i \textbf{V}_i \right ) = \dot{R_i} .$$ (18.4)

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

$$\dot{R_{i}} = \sum_{j,k} Y_{l} Y_{k} \lambda _{kj}(l) - Y_{i} Y_{j} \lambda _{jk}(i) ,$$ (18.5)

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 (18.5). The mass diffusion velocities $\textbf{V}_i$ in (18.4) 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 distribution. Setting $\textbf{V}_i\equiv 0$ transforms (18.4) into

$${dY_i\over dt} = \dot{R_i} ,$$ (18.6)

which may be written in the more compact, standard form

$$\dot {{\bf y}} = {\bf f} ({\bf y}) .$$ (18.7)

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 (18.6) and (18.7) 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 (18.6) 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 (18.7) through the terms

$$\dot {Y} (^4He) = - Y(^4He) Y(^{12}C) R + \ldots$$

$$\dot {Y} (^{12}C) = - Y(^4He) Y(^{12}C) R + \ldots$$ (18.8)

$$\dot {Y} (^{16}O) = + Y(^4He) Y(^{12}C) R + \ldots ,$$

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}}$

$$J(^4He,^4He) = - Y(^{12}C) R + \ldots \hskip 0.5in J(^4He,^{12}C) = - Y(^4He) R + \ldots$$

$$J(^{12}C,^4He) = - Y(^{12}C) R + \ldots \hskip 0.5in J(^{12}C,^{12}C) = - Y(^4He) R + \ldots$$ (18.9)

$$J(^{16}O,^4He) = + Y(^{12}C) R + \ldots \hskip 0.5in J(^{16}O,^{12}C) = + Y(^4He) R + \ldots .$$

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

18.1.2.1 Two linear algebra packages: MA28 and GIFT

As mentioned in the previous section, the Jacobian matrices of nuclear reaction networks tend to be sparse, and they become more sparse as the number of isotopes increases. Since implicit or semi-implicit time integration schemes generally require solving systems of linear equations involving the Jacobian matrix, taking advantage of the sparsity can significantly reduce the CPU time required to solve the systems of linear equations.

The MA28 sparse matrix package used by FLASH is described by Duff, Erisman, & Reid (1986). 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 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.

GIFT is a program which generates Fortran subroutines for solving a system of linear equations by Gaussian elimination (Gustafson, Liniger, & Willoughby 1970; Müller 1997). The full matrix $\tilde{{\bf A}}$ is reduced to upper triangular form, and backsubstitution with the right-hand side b yields the solution to $\tilde{{\bf A}} \cdot {\bf x} = {\bf b}$. GIFT generated routines skip all calculations with matrix elements that are zero; in this restricted sense, GIFT generated routines are sparse, but the storage of a full matrix is still required. It is assumed that the pivot element is located on the diagonal and no row or column interchanges are performed, so GIFT generated routines may become unstable if the matrices are not diagonally dominant. These routines must decompose the matrix for each right-hand side in a set of linear equations. GIFT writes out (in Fortran code) the sequence of Gaussian elimination and backsubstitution steps without any do loop constructions on the matrix $A(i,j)$. As a result, the routines generated by GIFT can be quite large. For the 489 isotope network discussed by Timmes (1999), GIFT generated $\sim$ 5.0$\times$10$^7$ lines of code! Fortunately, for small reaction networks (less than about 30 isotopes), GIFT generated routines are much smaller and generally faster than other linear algebra packages.

The FLASH runtime parameter algebra controls which linear algebra package is used in the simulation. algebra = 1 is the default choice and invokes the sparse matrix MA28 package. algebra = 2 invokes the GIFT linear algebra routines.

18.1.2.2 Two time integration methods

One of the time integration methods used by FLASH 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

$${\bf y}^{n+1} = {\bf y}^n + \sum_{i=1}^4 b_i \Delta_i ,$$ (18.10)

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

$$(\tilde{{\bf 1}}/\gamma h - \tilde{{\bf J}}) \cdot \Delta_1 = {\bf f} ({\bf y}^n)$$ (18.11)

$$(\tilde{{\bf 1}}/\gamma h - \tilde{{\bf J}}) \cdot \Delta_2 = {\bf f} ({\bf y}^n + a_{21}\Delta_1) + c_{21}\Delta_1/h$$ (18.12)

$$(\tilde{{\bf 1}}/\gamma h - \tilde{{\bf J}}) \cdot \Delta_3 = {\bf f} ({\bf y}^n + a_{31}\Delta_1 + a_{32}\Delta_2) + (c_{31}\Delta_1 + c_{32}\Delta_2)/h$$ (18.13)

$$(\tilde{{\bf 1}}^\gamma h - \tilde{{\bf J}}) \cdot \Delta_4 = {\bf f} ({\bf y}^n + a_{31}\Delta_1 + a_{32}\Delta_2) + (c_{41}\Delta_1 + c_{42}\Delta_2 + c_{43}\Delta_3)/h .$$ (18.14)

$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 makes use of the routine GRK4T given by Kaps & Rentrop (1979).

Another time integration method used by FLASH 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,

$$h = H/m$$

$$(\tilde{{\bf 1}} - \tilde{{\bf J}}) \cdot \Delta_0 = h {\bf f} ({\bf y}^n)$$ (18.15)

$${\bf y}_1 = {\bf y}^n + \Delta_0 .$$

Then from $k=1,2,\ldots,m-1$

$$(\tilde{{\bf 1}} - \tilde{{\bf J}}) \cdot {\bf x} = h {\bf f}({\bf y}_{k}) - \Delta_{k-1}$$

$$\Delta_k = \Delta_{k-1} + 2 {\bf x}$$ (18.16)

$${\bf y}_{k+1} = {\bf y}_k + \Delta_k ,$$

and closure is obtained by the last stage

$$(\tilde{{\bf 1}} - \tilde{{\bf J}}) \cdot \Delta_m = h [ {\bf f} ({\bf y}_m) - \Delta_{m-1} ]$$

$${\bf y}^{n+1} = {\bf y}_m + \Delta_m .$$ (18.17)

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) 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 Equations Eqn:BD 1 - Eqn:BD 3, 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 runtime parameter 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.

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

18.1.4 Energy generation rates and reaction rates

The instantaneous energy generation rate is given by the sum

$$\dot {\epsilon}_{\rm nuc} = N_A \sum_i {dY_{i}\over dt} .$$ (18.18)

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 program to use the average nuclear energy generated over a timestep

$$\dot {\epsilon}_{\rm nuc} = N_A \sum_i {\Delta Y_i \over \Delta t} .$$ (18.19)

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 program for use with the operator splitting technique.

The tabulation of Caughlan & Fowler (1988) is used in FLASH 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. The runtime parameter 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 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.

18.1.5 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. 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 enucDtFactor controls this ratio. The timestep limiter is calculated as

$$\Delta t_{burn} = {\tt enucDtFactor} \cdot \frac{E_{int}}{E_{nuc}}$$ (18.20)

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.