Explicit integration of extremely stiff reaction networks: asymptotic methods

In many scientific and technical contexts, one encounters phenomena that may be modeled by fluxes transferring population between sources and sinks for various species. Examples include kinetic processes that modify abundances and transfer energy in atomic, molecular and nuclear systems; geochemical, climate and other environmental systems; electrical circuits; economic models; and population dynamics. Terminology varies, but let us refer generically to these sources and sinks as boxes, and call the resulting systems of boxes connected by fluxes reaction networks. Such systems are commonly modeled by a coupled set of differential equations that describe a continuous flow of population through the boxes.

The reaction network is often classified as a stiff system, which we shall define to be a system of equations containing multiple timescales ranging over many orders of magnitude [3–6]. Most physical systems involve important processes operating on very different timescales, so realistic problems tend to be at least moderately stiff. Some, such as those encountered in many astrophysical thermonuclear burning applications, are extremely stiff, with fastest and slowest timescales in the problem differing by as much as 10–20 orders of magnitude. In stiff systems the timestep constraints are set by numerical stability requirements rather than accuracy considerations. Hence, explicit numerical integration of stiff systems is usually impractical because the maximum stable timestep is far too small for efficient solutions (see, e.g., [5, 6]). This is commonly addressed by employing implicit or semi-implicit stiff solvers that are stable, but that require time-consuming matrix solutions.

A given box in a reaction network very often is connected strongly only to a few other boxes. For example, the explosive burning conditions encountered in astrophysical novae, x-ray bursts or supernovae may require reaction networks with hundreds to thousands of nuclear isotopes. However, because Coulomb repulsion strongly favors the capture of light nuclei over heavy-ion captures, individual isotopes are typically connected directly to other isotopes through (at most) ∼ 10 reactions of consequence, and under many conditions no more than 2–3 reactions are important for a given isotope. Such restrictions on the direct box reaction coupling imply that the matrices appearing in the iterative implicit solution are sparse. Although various methods are available to deal with sparse matrices, in practice codes for solving large reaction networks have not always exploited sparseness in particularly effective ways.

For example, in astrophysical calculations with implicit solvers in large networks (hundreds of species), one often finds that 80–90% of the processor time is consumed in matrix operations [7, 8]. Efficient algorithms exist for the required matrix algebra (with incremental improvements in them over time), but the matrix nature of the core problem implies that the time required for an implicit solution grows nonlinearly with the size of the network. In typical working codes for large-scale applications, increasing the size of the network increases the time for solution, often quadratically, sometimes as much as cubically, until there are enough boxes in the network to justify the overhead of sparse-matrix methods with more favorable scaling. In thermonuclear networks, for example, it is often found that the overhead required to implement sparse-matrix iterative solutions is not justified until there are a couple of hundred boxes in the network. Thus, many of the present implicit stiff-network algorithms do not scale very gracefully to larger networks.

We are primarily interested in the most ambitious applications of large networks, where the reaction network is only a portion of a larger problem. Let us take as representative the astrophysical thermonuclear reaction networks, where a proper description of the overall problem typically requires multi-dimensional hydrodynamics or radiation hydrodynamics coupled tightly to a large thermonuclear reaction network. The hydrodynamical evolution controls the conditions in the network such as temperature and density, and the network influences the hydrodynamic evolution strongly through energy production and modification of composition variables. As a consequence of the limitations discussed in the preceding paragraphs, the solution of large networks by the usual approaches is time-consuming and a few multi-dimensional calculations have attempted to couple the element and energy production strongly to the hydrodynamics with a network of realistic complexity⁹. The most ambitious approaches use very small networks, perhaps tuned empirically to get critical quantities like energy production correct on average, coupled to the hydrodynamical simulation. In many calculations even this is not done and the network is replaced entirely by parameterization. Then a more complete network is run in a separate 'post-processing' step, where fixed thermodynamical profiles computed in the hydrodynamical simulation with the small network are used to specify the variation of variables such as temperature and density with time.

Astrophysical thermonuclear networks have been used here for illustration, but many problems of scientific and technical interest exhibit a similar complexity. Examples include astrochemical kinetics, where one must model large chemical evolution networks in contracting molecular clouds, or combustion chemistry, where chemical burning networks are strongly coupled to simulations of the dynamics of the air and fuel mixture. Physically realistic networks in such contexts would often be quite large. In the combustion of larger hydrocarbon molecules or studies of soot formation, hundreds to thousands of reacting species undergoing as many as 10 000 reactions may be encountered [6], and in supernova explosions hundreds to thousands of nuclear isotopes with tens of thousands of reaction couplings make non-zero contributions [8]. For such cases, one finds that current techniques do not allow for a coupling of realistic reaction networks to the full dynamics of the problem and often severely truncated or highly schematic reaction networks have been used in even the most realistic simulations.

To be definite, we shall assume that the coupling of reaction networks is done using operator splitting, where the hydrodynamical solver is evolved for a numerical timestep holding network parameters constant, and then the network is evolved over the time corresponding to the hydrodynamical timestep holding the new hydrodynamical variables constant. This places two basic constraints on methods:

• (i)  
  At the end of each hydrodynamical timestep the network must be advanced with new initial conditions. Thus, algorithms must be capable of rapid initialization and must not depend in a complex way on the conditions from previous time intervals.
• (ii)  
  With modern processors, existing algorithms are reasonably adequate for many post-processing calculations. In contrast, for the operator-split, parallel processing environment, which is of our interest here, it is highly desirable that solution of the network over a hydrodynamic timestep be fast enough that it does not require time substantially larger than that for the hydrodynamical solution alone.

Let us elaborate further on this second point. If a single processor were used to calculate both the hydrodynamical evolution and the network evolution in one hydrodynamical zone, the network evolution over a hydrodynamical timestep interval must be fast enough to not slow the calculation too much relative to the hydrodynamical evolution alone. If we take the point of view that we are willing to tolerate longer compute times in the interest of a much more realistic calculation, but not longer by multiple orders of magnitude, we estimate that the network must be capable of evolving over the time interval corresponding to the hydrodynamical timestep in roughly a second or less wall clock time.

We take as representative the multidimensional, adaptive-mesh, explicit hydrodynamical FLASH code [9] applied to Type Ia supernova simulations on large parallel systems. The explicit hydrodynamical timestep will be limited overall by the Courant time (roughly, because stability requires a hydrodynamical timestep not larger than the sound-crossing time for the zone), and more stringently in zones of rapid burning where temperature and density may be changing rapidly. In current Type Ia supernova sub-grid model simulations, the Courant time would typically be 10⁻⁴ s or smaller over most of the grid for the timescale relevant for the main part of the explosion, with rapid nuclear burning and associated temperature changes limiting the hydrodynamical timestep to 10⁻⁸–10⁻¹⁰ s for some ranges of times. For qualitative estimates, let us take as representative that a typical network integration for a single hydrodynamical timestep will be over an interval ∼ 10⁻⁸ s during the time of strong burning and ∼ 10⁻⁶ s over much of the approach to equilibration after strong burning.

In FLASH, many spatial zones will be assigned to a single Message Passing Interface rank on a parallel system. Therefore, in the absence of node-level parallelism (for example, with OpenMP), the network must be capable of calculating a number of hydrodynamic time intervals in a second or less if we wish to calculate an independent network integration for each zone. Let us take for estimation purposes that we wish to be able to reliably integrate the network in 1000 independent zones over a time interval of say 10⁻⁶ s on a single processor in 1 s wall clock time, with each network containing several hundreds of isotopes. This places extremely strong startup and speed constraints on the required network. The explicit algorithms discussed here are capable of perhaps 10⁴ network timesteps per second on a single processor with the present technology for a network with ∼ 150 isotopes, so our goals require a network algorithm that can integrate a time interval of the order of 10⁻⁶ s in not more than ∼ 10–100 timesteps, implying average stable and accurate timesteps at least as large as 10⁻²–10⁻³ times the elapsed integration time for the corresponding hydrodynamical evolution. See figure 1.

Zoom In Zoom Out Reset image size

Such large timesteps are often possible with implicit and semi-implicit algorithms, but those methods are inefficient at computing each timestep; explicit methods can compute a timestep efficiently, but timesteps this large are unthinkable with a normal explicit algorithm because they would be unstable in most realistic situations. In this and the other two papers [1, 2] of this series we shall demonstrate stabilization methods for explicit integration that realize such competitive integration timesteps in a variety of examples. Thus we shall reopen the discussion of whether explicit methods, with their faster computation of timesteps and more favorable scaling with network size, are practical for large, stiff networks.

The general task is to solve efficiently N coupled ordinary differential equations

$$\begin{eqnarray} \frac{\mathrm{d}y_i}{\mathrm{d}t} &=& F_i(y,t) = \sum_j F_{ij}(t) \nonumber\\ &\equiv& {F^+_{i}} (t)- {F^-_{i}}(t) = {F^+_{i}} (t) - k_i(t) y_i(t) \end{eqnarray} \noindent \tag{ 1 }$$

subject to appropriate boundary conditions. In this expression, the y_i(i = 1···N) describe the dependent variables (typically measures proportional to number densities for species), t is the independent variable (the time in our examples), the fluxes between species i and j are denoted by F_ij, and k_i(t) is the effective rate for all processes depleting the species i. The sum for each variable i is over all species j coupled to i by a non-zero flux F_ij, and for later convenience we have decomposed the flux into a component F⁺_i that increases the abundance of y_i and a component F⁻_i = k_iy_i that depletes it. For an N-species network there will be N such equations in the population variables y_i, generally coupled to each other because of the dependence of the fluxes on the different y_j. (For notational simplicity, we will not always display the i index explicitly on the right side in our equations.) To keep the discussion general, the variable y_i(t) will be used in most of our equations, but for the specific astrophysical examples that follow we shall replace the generic population variables y_i with the mass fraction X_i, which satisfies

$$\begin{equation} X_i = \frac{n_iA_i}{\rho N{_{\mbox{\scriptsize A}}} },\quad \sum_i X_i =1, \end{equation} \tag{ 2 }$$

where $N{_{\mbox {\scriptsize A}}}$ is Avogadro's number, ρ is the total mass density, A_i is the atomic mass number and n_i is the number density for the species i.

3.1. Example: stiffness and stability in the carbon–nitrogen–oxygen cycle

The carbon–nitrogen–oxygen (CNO) cycle that powers main-sequence stars more massive than the Sun provides a graphic illustration of stiffness in a relatively simple system of large physical significance. The CNO cycle, displayed in figure 2, uses isotopes of carbon, nitrogen and oxygen as catalysts to convert hydrogen into helium. We shall illustrate by considering the primary part of the CNO cycle illustrated on the right side of the figure. If the thermonuclear network corresponding to the main part of the CNO cycle is integrated under typical CNO cycle temperature and density conditions by explicit forward Euler methods using standard rates [10] for the reactions and constant timesteps, the integration is stable for timesteps less than or equal to 285.7 s, but becomes catastrophically unstable for a timestep of 285.8 s or more. This instability threshold is precisely two over the fastest rate for the transitions in the network, which corresponds to the β-decay of ¹⁵O to produce ¹⁵N.

Zoom In Zoom Out Reset image size

We now show that this instability for forward differencing of the CNO cycles arises because rapidly decreasing small populations can become negative in an explicit integration if the timestep is too large. These negative populations export unphysical negative population that can destabilize the system because they can lead to exponentially growing solutions in small components that ultimately couple to the larger components. Let us elaborate through the use of a simple model illustrated in figures 3–5, which will generalize a discussion that may be found in [6].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. A simple model for one form of stiffness

In figure 3, we assume a coupled system in which there are many components varying on a relatively long timescale (illustrated schematically by the curves at the top) and a single component that is small and exponentially decreasing on a much faster timescale, with a typical behavior $y(t) = \exp (-kt)$ so that dy(t)/dt = −ky(t). Let us consider timesteps that remain small compared with the larger timescales in the system (so that for a single step, we are in an approximately adiabatic situation and the populations varying on long timescales may be considered frozen), but comparable with or larger than the timescale set by 1/k.

To advance the solution for the fast component from t_n to t_n+1 by the explicit Euler method, we take a timestep h and extrapolate the solution using the derivative dy/d t = −ky evaluated at (t_n,y_n). The triangle of height Δy in figure 3 summarizes, where (t_n+1,y_n+1) represents the numerical solution and $(t_{n+1},y_{n+1} )_{\mbox {\scriptsize exact}}$ is the exact solution. From this construction, we observe that [6]:

• (i)  
  For h < 1/k, the forward Euler approximant for step n + 1 will yield a value of y_n+1 that lies between 0 and y_n, so that |y_n+1| < |y_n|.
• (ii)  
  For 1/k < h < 2/k, the sign of y_n+1 will be negative, but again |y_n+1| < |y_n|.
• (iii)  
  For h ⩾ 2/k the sign of y_n+1 will be negative and |y_n+1| ⩾ |y_n|.

For the forward Euler approximation to $y(t) = \exp (-kt)$ we have y_n+1 = (1 − hk)y_n, so that by iterating for m successive steps of fixed size h (with mh still considerably less than the longer timescales in the system, so that the adiabatic approximation remains valid), we obtain

$$\begin{equation} y_{n+m} = (1-hk)^m y_n. \end{equation} \tag{ 3 }$$

This converges toward the correct value of zero at larger times only if the product hk lies between 0 and 2, implying that the maximum value of h that yields a convergent solution is bounded by h < 2/k. Thus, in figure 3 the maximum stable value of h is less than 2/k and any point (t_n+1,y_n+1) extrapolated by forward Euler integration that lies in the shaded zone (that is, y_n+1 ⩽ − y_n) will be unstable under the iteration (3), diverging to infinity rather than converging to the correct value of zero. Because the small, fast components are coupled to the other components of the network, this divergence will quickly destabilize the entire network. Let us illustrate with a concrete example. Figure 4 corresponds to the forward Euler solution of y_n+1 = (1 − hk)y_n with k = 100 and timesteps ranging between h = 0.008 and 0.03. In this example, we observe that

• (i)  
  For h < 2/k, the solution converges to zero.
• (ii)  
  For h = 2/k, the solution oscillates between positive and negative values of the same absolute value and neither converges nor diverges.
• (iii)  
  For h > 2/k, the solution diverges to ±∞ under successive iterations, with the divergence exhibiting exponential behavior for larger times.

Figure 5 applies the preceding model of the origin of the stiffness instability to the CNO cycle. As noted earlier, (for the parameters assumed) ordinary explicit Euler integration is highly unstable for a timestep larger than 285.7 s. This critical timestep is exactly two divided by the fastest rate parameter in the system, which is that for the β-decay of ¹⁵O to ¹⁵N. Expressing the network in matrix form and diagonalizing indicates that the origin of this instability lies in eigenvalues that exceed unity in absolute value if $h>2/k_{\mbox {\scriptsize max}}$. Generalizing the discussion associated with equation (3), for a coupled set of equations a finite-difference iteration for m steps entails a matrix raised to the power m applied to the original y vector. This is guaranteed to converge only if no eigenvalue of the matrix exceeds unity in magnitude.

Less abstractly, figure 5 indicates that the origin of this instability for standard explicit integration is the tendency of the ¹⁵O population to become negative for large explicit timesteps. This is a more complex system than the previous simple decaying-exponential example because the ¹⁵O population is depleted by the β-decay but also replenished by the (p,γ) proton capture reaction on ¹⁴N (see figure 2 and the inset diagram on the left side of figure 5). Nevertheless, we see that the origin of the stiffness instability is very similar to that illustrated in the preceding simple example: a small box population becomes unphysically negative because of taking too large a timestep, so in the next timestep the offending box exports an unphysical flux of negative population, triggering a divergence that rapidly compromises the entire network.

Motivated by the properties of stiffness illustrated in the previous examples, we first introduce a rather crude approximation for stabilizing explicit integration in stiff networks. The results of the previous section indicate that one form of stiffness instability is generated by box populations that become very slightly negative and destabilize the network when they are exported as fluxes to other boxes. Thus, we invoke a simple flux-limiting prescription that if the population of a box becomes negative we do not change the population itself (which would quickly violate conservation of probability) but we suppress all export of that negative population to other boxes. Formally, we require for all computed outgoing fluxes that $F_{ij} \rightarrow \max (F_{ij},0) \equiv \tilde F_{ij}$. In shorthand, we refer to this as suppression of negative flux and refer to the resulting algorithm as the flux-limited forward difference (FLFD) algorithm. We shall illustrate using the (explicit) forward Euler method supplemented by the flux-limiting prescription. Because we do not alter populations but only restrict their flow by this algorithm, it conserves probability. We may expect that some populations in the network will now be in error because of the flux-suppression criterion. However, such effects tend to involve the smallest populations (because they are the ones most easily made negative by the numerical error), so we may expect that this approximation could be a good one for the larger populations, with the error concentrated in the smallest populations. In the next section, we test this idea on a realistic thermonuclear network.

We illustrate the FLFD algorithm by application to astrophysical thermonuclear networks under nova conditions.

5.1. Hot carbon–nitrogen–oxygen burning

Hydrogen burning in novae occurs through the hot CNO cycle where, unlike for the normal CNO cycle depicted in figure 2, the rates for proton capture exceed those for β-decay and the reactions in figure 2 break out into a much broader range of reactions. Some representative isotopic abundances under nova conditions calculated using the FLFD algorithm are displayed in figure 6(a), with the results of a standard implicit calculation shown as dashed lines and FLFD calculations shown as symbols. An adaptive timestep was used in the FLFD integration, with the timestep adjusted to keep the populations transferred between boxes in a timestep for some key populations within a prescribed range. Note the very good agreement between implicit and explicit flux-limited methods over six orders of magnitude in the mass fractions in figure 6(a).

Zoom In Zoom Out Reset image size

5.2. Stiffness and stability

The results of the preceding section indicate that even a relatively crude limit on propagation of negative population can remove large amounts of stiffness instability and leads to an algorithm useful for realistic applications, even for the extremely stiff networks common in astrophysics. But we will now demonstrate that one can do much better than that by exploiting the algebraic structure of the differential equations to make a better approximation than just constraining the sign of propagating populations.

The first step is that we must look more deeply at the stiffness instability. As we now discuss, for the sort of large and very stiff networks that we are addressing here, there are several fundamentally different sources of stiffness instability that are often not clearly distinguished in the literature. The examples discussed to this point (and in many textbooks) emphasize the type of instability associated with small quantities that should strictly be non-negative becoming negative because of an overly ambitious numerical integration step. The discussion of the FLFD algorithm in section 4 illustrates that this type of instability can be removed by approximations as simple as not permitting unphysical negative quantities to influence the rest of the network. However, there are other stiffness instabilities that may be initiated even when no population variables become negative in an integration step. In this kind of instability, we end up having to take the difference of large numbers to obtain a result very near zero. The numerical errors that ensue in a standard explicit approach can then accumulate rapidly and destabilize the network, even before any abundances become negative. This may still be viewed as a stiffness instability because it results from a numerical integration trying to deal with very different timescales, but the origin of these timescales is different from that discussed above. In this case, the disparate timescales are the very rapid reactions driving the system to equilibrium contrasted with the very slow timescale associated with equilibrium itself (which tends to infinity).

As we now consider, this distinction is essential to what follows because these stiffness instabilities have essentially different solutions. Furthermore, we shall find that the second kind of instability can be divided into two subclasses requiring different stabilizing approximations, and that the approximations that we shall introduce in these cases will also naturally take care of the first class of stiffness instabilities because they will as a matter of course prevent the occurrence of negative probabilities in the network.

6.1. The approach to equilibrium

We shall use 'equilibrium' in a broad sense to mean a condition where the populations in a network are being strongly influenced by competition between terms of opposite sign on the right sides of the differential equations (1) governing their evolution. How do we measure the degree of equilibration in a large, stiff network? In terms of the coupled set of differential equations describing the network, we may distinguish two qualitatively different conditions:

• (i)  
  a macroscopic equilibration that acts at the level of individual differential equations;
• (ii)  
  a microscopic equilibration that acts at the level of individual terms within a given differential equation.

Let us consider each of these cases in turn. The differential equations that we must solve take the general form of equation (1), dy_i/dt = F⁺_i − F⁻_i, where the total flux has been decomposed into a component F⁺_i increasing the population y_i and a component F⁻_i depleting the population y_i in a given timestep.

6.1.1. Macroscopic equilibration.

One class of approximations that we will investigate depends upon assuming that F⁺_i − F⁻_i → 0 (asymptotic approximations) or F⁺_i − F⁻_i →  constant (steady-state approximations). We shall refer to these conditions as a macroscopic equilibration, since they involve the entire right side of a differential equation in equation (1) tending to zero or a finite constant¹⁰. We shall introduce approximations exploiting this whereby whole differential equations are removed from the numerical integration for a network timestep in favor of algebraic approximate solutions for that timestep. Such approximations do not reduce the number of equations to integrate, but they reduce the number of equations integrated numerically by forward difference. They reduce stiffness for any remaining equations that are integrated numerically in the timestep because removing the equations satisfying these conditions tends to reduce the disparity in timescales for the remaining equations.

6.1.2. Microscopic equilibration.

In equation (1), F⁺_i and F⁻_i for a given species i generally each consist of a number of terms depending on the other populations in the network,

$$\begin{eqnarray} \frac{\mathrm{d}y_i}{\mathrm{d}t} &=&F^+_{i} - F^-_{i} \nonumber\\ &=& (f_1^+ + f_2^+ + \cdots)_i - (f_1^- + f_2^- + \cdots)_i\nonumber\\ &=& (f_1^+- f_1^-)_i + (f_2^+ - f_2^-)_i + \cdots = \sum_j (f^+_j - f^-_j)_i. \end{eqnarray} \tag{ 4 }$$

At the more microscopic level, groups of individual terms on the right side of equation (4) may come approximately into equilibrium (the sum of their fluxes tends to zero), even if macroscopic equilibration conditions are not satisfied. The simplest possibility for this microscopic equilibration is that individual forward–reverse reaction pairs such as $A+B \rightleftharpoons C$, which will contribute flux terms with opposing signs on the right sides of differential equations in which they participate, come approximately into equilibrium¹¹. Then we may consider an algebraic approximation that removes groups of such terms from the numerical integration, replacing their sum of fluxes identically with zero. This will not generally reduce the number of equations to be integrated numerically in a network timestep, but it can reduce dramatically the stiffness of those equations by removing terms with fast rates from the equations, thereby reducing the disparity between the fastest and slowest timescales in the system.

Such considerations will be the basis of the partial equilibrium methods that will be discussed in depth in the third paper in this series [2]. There we shall also demonstrate two important general conclusions: (i) Approximations based on microscopic equilibration are much more efficient at removing stiffness than those based on macroscopic equilibration because they target more precisely the sources of stiffness in the network. (ii) The most powerful approach will be to use macroscopic and microscopic approximations simultaneously in the same set of equations, because they can complement each other in removing stiffness from the equations to be integrated numerically.

6.2. A quantitative measure of microscopic equilibration

The FLFD approximation described in section 5 illustrates the basic principle that explicit integration in non-equilibrium situations can be stabilized by forbidding the propagation of negative flux. However, the FLFD algorithm represents only a zero-order solution to this problem (a yes/no decision on whether flux components are allowed to propagate) that can be improved substantially by exploiting the structure of the coupled equations to replace the numerical solution with an algebraic approximation that is strictly non-negative. Although the approximations that we now discuss have been implemented in some form in the earlier literature [6, 13–15], we shall demonstrate that our implementation appears to be much more successful than previous applications to large, extremely stiff networks, and we shall reach different conclusions about these methods than those reached in earlier publications.

7.1. Some explicit asymptotic approximations

The differential equations that we must solve take the form given by equation (1). Generally, F⁺_i and F⁻_i for a given species i each consist of a number of terms depending on the other populations in the network. The depletion flux for the population of species i will be proportional to y_i,

$$\begin{equation} {F^-_{i}} = (k_1^i + k_2^i + \cdots + k_m^i)y_i \equiv k^iy_i, \end{equation} \tag{ 5 }$$

where the kⁱ_j are rate parameters (in units of time⁻¹) for each of the m processes that can deplete y_i, which may depend on the populations y_j and on local thermodynamic variables such as temperature and density. The characteristic timescales τⁱ_j = 1/kⁱ_j will differ by many orders of magnitude in the systems of interest, implying that the equations are stiff. From equation (5), we may define the effective total depletion rate kⁱ for y_i at a given time, and a corresponding timescale τⁱ as

$$\begin{equation} k^i \equiv \frac{{F^-_{i}}}{y_i}\quad\tau^i = \frac{1}{k^i}, \end{equation} \tag{ 6 }$$

permitting equation (1) to be written as

$$\begin{equation} y_i = \frac{1}{k^i} \left( {F^+_{i}} - \frac{\mathrm{d}y_i}{\mathrm{d}t} \right). \end{equation} \tag{ 7 }$$

Thus, in a finite-difference approximation at timestep t_n we have

$$\begin{equation} y_i(t_n) = \frac{{F^+_{i}} (t_n)}{k^i(t_n)} - \frac{1}{k^i(t_n)} \left. \frac{\mathrm{d}y_i}{\mathrm{d}t} \right|_{t=t_n}. \end{equation} \tag{ 8 }$$

We now define the asymptotic limit for the species i to be F⁺_i ≃ F⁻_i, implying from equation (1) that dy_i/dt ≃ 0. In this limit, equation (8) gives a first approximation y⁽¹⁾_i(t_n) and local error E⁽¹⁾_n, respectively, for y_i(t_n) as

$$\begin{equation} y^{(1)}_i (t_n) = \frac{{F^+_{i}}(t_n)}{k^i(t_n)}\quad E_n^{(1)} \equiv y(t_n) - y^{(1)}(t_n) = -\frac{1}{k(t_n)} {\frac{\mathrm{d}y}{\mathrm{d}t}}(t_n). \end{equation} \tag{ 9 }$$

For small dy_i/dt, we may then get a correction term by writing the derivative term in equation (8) as

$$\begin{equation} \frac{\mathrm{d}y}{\mathrm{d}t}(t_n) = \frac{1}{\Delta t} \left( y_i(t_n) - y_{i}(t_{n-1}) \right) +\frac{1}{\Delta t} \left( E_n^{(1)} - E_{n-1}^{(1)}\right) + O(\Delta t), \end{equation} \tag{ 10 }$$

where O(x) denotes the order of x. If we retain only the first term and approximate y(t) by y⁽¹⁾(t), we obtain that

$$\begin{eqnarray} \left. \frac{\mathrm{d}y}{\mathrm{d}t}\right|_{t=t_n} & \simeq & \frac{1}{\Delta t} \left( y_i^{(1)}(t_n) - y_{i}^{(1)}(t_{n-1}) \right)= \frac{1}{\Delta t} \left( \frac{{F^+_{i}}(t_n)}{k^i(t_n)} -\frac{{F^+_{i}} (t_{n-1})}{k^i(t_{n-1})} \right). \end{eqnarray} \tag{ 11 }$$

Therefore, the estimate for y is improved to

$$\begin{equation} y_n^{(2)} \simeq \frac{F^+_n}{k_n} - \frac{1}{k_n\Delta t} \left(\frac{F^+_n}{k_n} - \frac{F^+_{n-1}}{k_{n-1}} \right), \end{equation} \tag{ 12 }$$

where we now employ compact index notation y_n ≡ y_i(t_n), and so on, and have dropped the species index i to avoid notational clutter. Because we are approximating the derivative term, we expect that equation (12) is valid only if the second term is small, implying that our approximation becomes more valid if kΔt is large.

The preceding discussion implements an asymptotic approximation in a very simple way. Other, more sophisticated, asymptotic approximations may be derived. For example, if we retain the full expression for the derivative in equation (10) and substitute in equation (8), we obtain

$$\begin{eqnarray*} &&y^{(2)}_n = \frac{F^+_n}{k_n} -\frac{1}{k(t_n) \Delta t} \left( y(t_n) -y(t_{n-1}\right) - \frac{1}{k(t_n) \Delta t} \left( -E_n^{(1)} + E_{n-1}^{(1)}\right). \end{eqnarray*}$$

Setting y(t_n) = y⁽²⁾_n and solving for y⁽²⁾_n gives

$$\begin{equation} y_n^{(2)} = \frac{1}{1+k_n\Delta t} \left(y_{n-1} + F^+_n \Delta t \right), \end{equation} \tag{ 13 }$$

where a term of the order of (Δt)²/(1 + k(t_n)Δt) has been discarded. Another approach is to use a predictor–corrector scheme [6, 13, 15]. We may solve equation (1) by finite difference, replacing quantities on the right side of equation (1) with averaged quantities:

$$\begin{equation} \frac{y_n - y_{n-1}}{\Delta t} = \frac{1}{2} (F^+_n + F^+_{n-1}) -\frac{y_n+ y_{n-1}}{\tau_{n} + \tau_{n-1}}, \end{equation} \tag{ 14 }$$

where for notational convenience we have traded the rate parameters k_i for the corresponding timescales τ_i ≡ 1/k_i. Solving this equation for y_n gives

$$\begin{equation} y_{n} = \frac{(\tau_n + \tau_{n-1} - \Delta t)y_{n-1} + \frac{1}{2} (\tau_n + \tau_{n-1})(F^+_n + F^+_{n-1})\Delta t}{\tau_n + \tau_{n-1} + \Delta t} , \end{equation} \tag{ 15 }$$

which is not easy to use because the τ_n and F⁺_n are implicit functions of the y_n that we seek. We obtain a usable explicit algorithm by solving an approximate form of equation (15) to get a first guess $y{^{\mbox {\scriptsize p}}}$ for y_n (the predictor step), using that approximate result to estimate values of τ_n and F⁺_n, and then using these values obtained from the predictor step to solve equation (15) (the corrector step). Specifically,

• The predictor estimate $y{^{\mbox {\scriptsize p}}}$ is obtained by setting τ_n = τ_n−1 ≡ τ₀ and F⁺_n = F⁺_n−1 ≡  F⁺₀ in equation (15).
• The corrector value $y{^{\mbox {\scriptsize c}}}$ is obtained by substituting $F_{{\mbox{\scriptsize p}}}^{+}$ for F⁺_n, and $\tau_{{\mbox{\scriptsize p}}}$ for τ_n in equation (15).

This yields a predictor and a corrector [6]

$$\begin{eqnarray} y{^{\mbox{\scriptsize p}}} &=& \frac{y^{0}(2\tau_{0}-\Delta t)+2 F_{0}^{+}\tau_{0} \Delta t}{2\tau_{0} + \Delta t} \qquad \mbox{(Predictor)}\nonumber\\ y{^{\mbox{\scriptsize c}}} &=& \frac{y^{0}(\tau{_{\mbox{\scriptsize p}}} + \tau_{0} - \Delta t) + \frac{1}{2}\Delta t (F{_{\mbox{\scriptsize p}}^{+}} + F_{0}^{+})(\tau{_{\mbox{\scriptsize p}}} + \tau_{0})}{\tau{_{\mbox{\scriptsize p}}} + \tau_{0} + \Delta t} \qquad \mbox{(Corrector)} \end{eqnarray} \noindent \tag{ 16 }$$

where '0' denotes initial quantities and 'p' denotes quantities computed using the results of the predictor step. We shall test the asymptotic approximations (12), (13) and (16) in the examples below.

7.2. Mathematical properties of asymptotic approximations

7.3. An asymptotic flux-limiting algorithm

7.4. A simple adaptive timestepper

We shall be comparing explicit and implicit methods using codes that are at very different stages of development and optimization. We assume that for codes at similar levels of optimization the primary difference between explicit and implicit methods would be in the extra time spent in implicit-method matrix operations. Hence, if the fraction of time spent on linear algebra is f for an implicit code, we assume that an explicit code at a similar level of optimization could compute a timestep a factor of F = 1/(1 − f) faster.

We have used the implicit backward Euler code XNet [11] to estimate the factors F for various networks [17, 18]. For small networks, where dense solvers are most appropriate, the matrix build and solve corresponds typically to 50–75% of the total computational cost, implying that F ∼ 2–4. For networks with hundreds of species, where sparse-matrix solvers are more appropriate, the matrix build and solve typically consumes 80–90% of the computational cost, implying that F ∼ 5–10. Based on these results, we adopt the factors F displayed in table 1 for the networks investigated here. Then we may compare roughly the speed of explicit versus implicit codes (possibly at different levels of optimization) by multiplying F by the ratio of implicit to explicit integration steps required for a given problem. This procedure has obvious uncertainties, and likely underestimates the speed of an optimized explicit versus optimized implicit code [2], but will give a useful lower limit on how fast the explicit calculation can be.

First we shall establish that the asymptotic algorithm is capable of correct integration of coupled sets of extremely stiff equations. We have tested this for a variety of calculations in two ways:

• (i)  
  comparisons with the results from standard implicit codes;
• (ii)  
  for some smaller networks where the corresponding integration time is not too prohibitive, comparison with explicit forward-Euler calculations made with timesteps short enough to be stable.

Our general finding is that the asymptotic algorithm outlined above gives the same results as standard implicit and explicit codes, even for the stiffest networks found in astrophysics applications, provided that the numerical timesteps are limited sufficiently to ensure conservation of overall probability in the network at the desired level of precision.

9.1. Tests against fully explicit calculations

Figure 7 illustrates an asymptotic approximation calculation that gives results essentially identical to the results from exact numerical integration. In this example, the sum of the mass fractions was constrained to deviate from unity by not more than 1% over the entire range of the calculation by requiring that if it deviates by more than 10⁻⁸ from unity in any one timestep, the timestep is reduced in size. We shall generally use this global 1% criterion for the examples presented in this paper. Higher precision may be obtained by tightening this constraint, but this is typically already more conservative than justified since input parameter uncertainties in realistic large thermonuclear networks and uncertainties in the coupled hydrodynamics may each be considerably larger than 1%.

Zoom In Zoom Out Reset image size

Having established that explicit asymptotic approximations can give correct results even for extremely stiff networks, we now turn to the question of their efficiency (and further tests of accuracy) by examining calculations of various astrophysical thermonuclear networks using this approximation. In this section, we shall use the simplest asymptotic approximation corresponding to equation (12). In section 10, we shall test the alternative asymptotic approximations, of equations (13) and (16).

9.2. Explicit asymptotic integration of the pp-chains

The solar pp-chains, which convert hydrogen into helium, provide a striking example of stiffness in a simple network¹². In figure 8, we illustrate integration of the pp-chains at a constant temperature and density characteristic of the core in the present-day Sun, using the asymptotic method and the implicit backward-Euler code XNet [11]. We see that the explicit asymptotic integration gives results for the mass fractions in rather good agreement with the implicit code over 20 orders of magnitude, and is generally taking timesteps dt ∼ 0.1t that are comparable with those for the implicit code over the entire range of integration. (The asymptotic method required 333 total integration steps versus 176 steps for the implicit code.)

Zoom In Zoom Out Reset image size

As shown in section 3, the maximum stable timestep for a standard explicit integration method may be estimated as the inverse of the fastest rate contributing to the network. This is illustrated by the dashed blue curve in figure 8(b). At late times the explicit integration timesteps are ∼ 10²¹ times larger than the maximum stable timestep for a normal explicit integration. The calculation illustrated in figure 8 takes a fraction of a second on a 3 GHz processor with the explicit asymptotic method (as does the implicit solver). In contrast, from figure 8(b) we estimate that a standard explicit method taking the largest stable fully explicit timestep would require ∼ 10²¹ s of processor time to compute this simple seven-element network to hydrogen depletion.

9.3. Simulations under nova conditions

The preceding example entailed an exceedingly stiff but rather small network. Let us now turn to an example that is highly stiff and involves hundreds of isotopes. In figure 9(a), we illustrate a calculation using the explicit asymptotic algorithm and a hydrodynamical profile shown in figure 9(b) that is characteristic of hot-CNO burning in a nova outburst. The explicit asymptotic timesteps are displayed in figure 10(a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In these simulations, we see that the asymptotic network takes stable and accurate timesteps corresponding to dt ∼ 0.1t over most of the time integration, except in the region of sharp temperature rise and strong burning, where dt ∼ (0.01–0.001)t. Over most of the range of integration after burning commences, the asymptotic solver timesteps (solid red curve in figure 10(a)) are a million or more times larger than the maximum stable timestep for a purely explicit method (blue dashed curve in figure 10(a)); at late times this disparity increases and by the end of the calculation the asymptotic timesteps are approximately 10¹⁰ times larger than would be stable for a normal explicit integration.

The generally large explicit asymptotic timesteps over the entire integration range illustrated in figure 10(a) are greater than or equal to those for a typical implicit code, as may be seen by comparing with the implicit (backward Euler) calculation timestepping curve shown in dotted green. In this calculation the implicit method required 1335 integration steps while the explicit asymptotic calculation required only 935 steps, and furthermore we expect that the explicit timesteps can be computed more quickly than the implicit timesteps. For a network with 134 isotopes, an optimized explicit code should calculate a timestep at least five times faster than typical implicit codes (table 1). These large explicit timesteps are possible because during the simulation many isotopes become asymptotic but very few reactions reach partial equilibrium, as illustrated in figure 10(b). Similar results for nova simulations were obtained with the asymptotic method in [17, 18] using a different nova hydrodynamical profile and a different reaction library. Allowing for uncertainties in the present comparisons that use differently optimized codes and somewhat different timestepping criteria, we conclude that the explicit asymptotic method may be intrinsically at least several times faster than a state-of-the art implicit code for simulations under nova conditions.

9.4. Scaling with network size

For explicit methods the computing time should scale linearly with the size of the network if the network is relatively sparse. This paper and the authors of [1, 2] use explicit methods modified by additional bookkeeping to remove stiffness, raising the issue of whether additional overhead destroys the linear scaling. In figure 11, we illustrate scaling of integration time with network size for the nova simulation of figure 9. The behavior is seen to be approximately linear, as expected for an explicit algorithm since no matrix inversions are required. The quasi-steady-state (QSS) method discussed in [1] requires an amount of overhead comparable to the asymptotic method and so should also be linear. The partial equilibrium method discussed in [2] adds some additional computation and logic to the asymptotic or QSS methods, but again no matrix operations are used, so it should not alter appreciably the linear scaling. Our tests so far (on networks as large as 20 species) suggest that the partial equilibrium bookkeeping overhead indeed is incremental and does not alter the network-size scaling substantially.

Zoom In Zoom Out Reset image size

All the results presented to this point have used the simplest asymptotic formula defined in equation (12). This section compares the results from equation (12) to those obtained with the more sophisticated formulae (13) and (16), for the case of an alpha network with constant T₉ = 5 and ρ = 1 × 10⁸ g cm⁻³. In figure 12, we display a composite of the computed mass fractions and the timestepping for the asymptotic approximations corresponding to these three formulae. We find that the more sophisticated asymptotic approximations in equations (13) and (16) give results for abundances that are very similar to those from equation (12), but can in some cases give more favorable timestepping by factors of several. These results are representative of tests on a variety of networks and we conclude that for typical thermonuclear networks equations (12), (13) and (16) yield similar results, except for possible differences by factors of up to 2–3 in computational speed.

Zoom In Zoom Out Reset image size

In previous sections evidence has been presented that, well-removed from equilibrium, asymptotic methods can provide stable and accurate integration of the stiffest large networks with timesteps that are comparable with those employed in standard implicit and semi-implicit stiff solvers. In practice, for astrophysical thermonuclear networks this means that timesteps are typically from 0.1 to 0.001 of the current time over most of the integration range, except for brief time periods where very strong fluxes are being produced and timesteps may need to be shorter to maintain accuracy. Since explicit methods can generally compute each timestep substantially faster than for implicit methods, this suggests that such methods offer a viable alternative to implicit solvers under those conditions.

However, the preceding statements are no longer true when substantial numbers of reaction pairs in the network begin to satisfy microscopic equilibrium conditions (according to the criteria given in section 6.2). Then the typical behavior for asymptotic approximations is for the timestep to become constant or only slowly increasing with integration time. (Indeed, we have already seen an indication of this behavior at late times in figure 12(b).) Figure 13 illustrates for two different networks coupled to a single-zone hydrodynamical simulation. We see that the 19-isotope network is only marginally able to keep up with the hydrodynamical timesteps, whereas the 150-isotope network lags orders of magnitude behind for hydrodynamical integration times later than about 10⁻⁴–10⁻⁵ s.

Zoom In Zoom Out Reset image size

The reason for this loss of timestepping efficiency for the asymptotic method is the approach to equilibrium at late times; this is documented in the inset plots for figure 13, which show the fraction of reactions in the respective networks that satisfy equilibrium conditions. The reason why the 19-isotope network is able to keep up much better than the 150-isotope network is also clear from the inset plots: the 19-isotope network slowly approaches 30–40% equilibration in this region, but the 150-isotope network, with much faster reactions because it includes fast proton and neutron reactions not found in the 19-isotope network, quickly reaches 70% equilibration. Since in most applications for extremely stiff astrophysical networks the physical phenomena require integration over many decades of time, this lag of the asymptotic timestepping as equilibrium is approached is disastrous for such approximations and they quickly lose out to implicit methods, which can continue to take large integration steps even nearing equilibrium (although they are less efficient than explicit methods at computing each timestep).

Another example of the failure of the asymptotic approximation to generate competitive timesteps is illustrated in figure 14. The asymptotic mass fractions calculated in figure 14(a) are quite accurate in comparison with the results from standard implicit codes, but the timestepping illustrated in figure 14(b) is not competitive with implicit methods at late times. For logt ∼ − 2, the asymptotic integration timesteps are only of the order of 10⁻⁸ s, whereas the implicit code is taking timesteps larger than 10⁻³ s at this point.

Zoom In Zoom Out Reset image size

The reason for this failure of the asymptotic approximation can be seen clearly in figures 14(b) and (d). In the region of rapid temperature rise illustrated in figure 14(c), the rates in the network—especially the photodisintegration reactions that are the inverses of the ⁴He capture reactions—increase dramatically. This causes the maximum stable explicit timestep to become much smaller, as illustrated in figure 14(b) by the dashed blue curve. The steadily increasing explicit timestep (solid green curve) intersects this maximum stable explicit timestep curve near logt = −8. But, from figure 14(d) we see that essentially no isotopes in the network satisfy the asymptotic condition at this point. Thus the integration is forced to use a standard explicit method and the timestep must (to maintain stability) follow the decreasing dashed blue curve until around logt = −4, when significant numbers of isotopes finally begin to satisfy the asymptotic condition and the explicit asymptotic algorithm is able begin to take timesteps larger than the explicit limit. However, at this point the asymptotic timesteps are already orders of magnitude smaller than those an implicit method would use, and the asymptotic method is able to increase the timestep by only about two orders of magnitude before the system reaches equilibrium. Thus, for the entire time range from logt ∼ − 8 until logt ∼ − 2 the explicit asymptotic method computes the network accurately but its timesteps lag many orders of magnitude behind those from standard implicit methods.

As we shall explain in considerable detail in the third paper of this series [2], the reason for the loss of efficiency for asymptotic methods as equilibrium is approached is that the asymptotic approximation removes a large amount of stiffness associated with macroscopic equilibration, but near (microscopic) equilibrium a fundamentally new source of stiffness becomes important and it is not generally removed by the asymptotic approximation. Indeed, we see clearly from figure 14(d) that the fraction of reactions in the network that satisfy equilibrium conditions increases rapidly beginning at logt ∼ − 7 and reaches unity by logt ∼ − 3. In this third paper, we shall describe a new implementation of partial equilibrium methods that can be used in conjunction with asymptotic methods to turn equilibrium from a liability into an asset and increase the explicit timestepping by orders of magnitude in the approach to equilibrium. In that paper, we will give examples suggesting that these methods are capable of giving timestepping competitive with that of implicit methods across the entire range of interesting physical integration times for a variety of extremely stiff reaction networks.

Explicit numerical integration can compute a timestep faster than implicit methods, and the time to compute a network explicitly scales linearly and therefore more favorably with network size than for implicit codes. Nevertheless, previous discussions of numerical integration for very stiff systems have concluded rather uniformly that standard explicit methods are not competitive with implicit methods for stiff networks because they are unable to take large enough stable timesteps. Improvements in explicit methods based on using asymptotic and steady-state limiting solutions to remove stiffness from the network have had some success for systems of moderate stiffness such as various chemical kinetics problems. However, for extremely stiff networks such as those encountered commonly in astrophysical thermonuclear calculations, it has been concluded in the previous literature that such methods are not competitive, failing even to give correct results, with timesteps that are far too short to be useful even if they gave correct results [6, 13]. This paper has presented evidence strongly challenging all these conclusions.

We have clearly identified three fundamentally different sources of stiffness in large networks, only the first of which is commonly emphasized in the literature:

• (i)  
  Situations where small populations can become negative if the explicit timestep is too large, with the propagation of this anomalous negative population leading to exponentially growing terms that destabilize the network.
• (ii)  
  Situations where the right sides of the differential equations expressed as d Y =F = F⁺ − F⁻ approach a constant derived from the difference of two large numbers (the total flux in F⁺ and the total flux out F⁻), and numerical errors in taking this (often small) difference destabilize the network if the timestep is too large.
• (iii)  
  Situations where on the right sides of the differential equations expressed in the form of equation (4) the net flux in specific forward–reverse reaction pairs (f⁺_i − f⁻_i) tends to zero as the system approaches equilibrium, leading to large errors if the timestep is too large because the net flux is derived from the difference of two large numbers and the timescale equilibrating the populations is short compared with the desired numerical timestep.

We have shown that these distinctions are important because different sources of stiffness require different approximations for their removal in an algebraically stabilized explicit integration.

Using the extremely stiff systems characteristic of astrophysical thermonuclear networks as a stringent test, we have shown that asymptotic methods are very successful at removing the first two types of stiffness, and give correct results, even for the stiffest of thermonuclear networks, provided that adequate attention is paid to conservation of probability (equivalently, nucleon number) in the network. Furthermore, we have shown various examples of stable and accurate timestepping with these methods in extremely stiff systems that are competitive with that of standard implicit codes, corresponding in some simple but physically important networks to integration timesteps that are as much as 20 orders of magnitude larger than the maximum timestep that would be stable in a standard explicit method. However, we have also shown that such methods give correct results but fail to exhibit competitive timestepping when the system approaches microscopic equilibrium and the third type of stiffness instability begins to dominate. In a following paper [2], we shall provide evidence for competitive timestepping, even in the approach to equilibrium, if the explicit asymptotic method is supplemented by partial equilibrium approximations designed specifically to deal with the third type of stiffness instability.

Taken together, this paper and the following ones on QSS methods [1] and partial equilibrium methods [2] present compelling evidence that algebraically stabilized explicit integration methods are capable of timesteps competitive with implicit integration methods for a variety of highly stiff reaction networks [19]. Since explicit methods can execute a timestep faster than an implicit method in a large network, our results suggest that algebraically stabilized explicit algorithms may be capable of performing as well as, or even substantially outperforming, implicit integration in a variety of moderate to extremely stiff applications. Because of the highly favorable linear scaling for explicit methods, this fundamentally new view of the efficacy of explicit integration for stiff equations may be particularly important for applications in any field where it is imperative that more realistic—and therefore much larger—networks be used in complex physical simulations.

We thank Tony Mezzacappa for useful discussions, Austin Harris for help with some of the calculations, Eric Lingerfelt for programming assistance and Christian Cardall for a careful reading of the manuscript. This research was sponsored by the Office of Nuclear Physics, US Department of Energy.