A Fully Explicit Integrator for Modeling Astrophysical Reactive Flows

Simulating complex astrophysical reacting flows is computationally expensive -- reactions are stiff and typically require implicit integration methods. The reaction update is often the most expensive part of a simulation, which motivates the exploration of more economical methods. In this research note, we investigate how the explicit Runge--Kutta--Chebyshev (RKC) method performs compared to an implicit method when applied to astrophysical reactive flows. These integrators are applied to simulations of X-ray bursts arising from unstable thermonuclear burning of accreted fuel on the surface of neutron stars. We show that the RKC method performs with similar accuracy to our traditional implicit integrator, but is more computationally efficient when run on CPUs.


INTRODUCTION
In an X-ray burst (XRB), a neutron star accretes He or H/He from a companion star building up a thin layer on its surface.The strong gravitational acceleration compresses this material to the point of thermonuclear runaway, producing a brief flash of X-rays (Galloway & Keek 2021).Observations suggest that the burning begins locally and spreads across the neutron star as a flame (Bhattacharyya & Strohmayer 2007).
Modeling XRB flames can be challenging due to the difference in burning and hydrodynamic timescales.Traditionally, explicit integration methods are applied to hydrodynamics while implicit integration methods are applied to the stiff reaction systems (Hujeirat & Rannacher 2001).For an ordinary differential equation (ODE) system, a first-order implicit discretization over a timestep ∆t is: For small ∆t, the change in the solution, ∆y, will be small as well, allowing the system to be Taylor-expanded around (t n+1 , y n+1 ).This requires the Jacobian matrix, J = ∂f/∂y to be computed and stored every timestep leading to computationally expensive linear system solves.In contrast, explicit methods do not need the Jacobian, instead constructing the update from the current solution.Explicit methods can use less memory, but they can be unstable when applied to the reacting system.In this paper, we study explicit methods with increased stability applied to the reacting system, with the goal of conserving memory and increasing the overall performance for large-scale simulations.Niemeyer & Sung (2014) showed that the Runge-Kutta-Chebyshev (RKC) method (Sommeijer et al. 1998) can be an efficient integrator when applied to chemical kinetics in terrestrial reactive-flow simulations.RKC is based on explicit Runge-Kutta methods, using the Chebyshev formulas to add stages to the integration to increase its stability region.At each step, an estimate of the spectral radius is needed, which is computed internally in RKC using the power method without constructing the Jacobian.
We implement RKC for astrophysical nuclear reactions in the Castro hydrodynamic code (Almgren et al. 2010(Almgren et al. , 2020)).Castro uses the AMReX adaptive mesh library (Zhang et al. 2019) and a general equation of state and nuclear reaction network, diffusion, and gravity.The advective update is fully explicit, using an unsplit piecewise parabolic method (Miller & Colella 2002).Several methods for coupling hydrodynamics and reactions are available, here we consider Strang-splitting (Strang 1968) and the simplified spectral deferred corrections (SDC) approach (Zingale et al. 2022).In all cases, the reaction system updates species and energy according to the reaction rates.
In the next section, we compare RKC and the implicit VODE integrator (Brown et al. 1989) for simulations of XRBs.We look at both CPU and GPU performance and consider both the Strang-split coupling and the simplified-SDC method.

RESULTS
We use the Castro flame wave setup (Eiden et al. 2020;Harpole et al. 2021) for our XRB flame model.The domain is 0.256 km × 1.024 km, with a base grid of 512 × 128 and two levels of refinement (jumps of 4×, then 2×), giving a 25 cm resolution overall.Reactions are modeled using the 13 isotope aprox13 network (Timmes et al. 2000).We ran eight simulations, each on eight nodes of the NERSC Perlmutter machine.CPU runs used 128 MPI tasks each with eight OpenMP threads; GPU runs used 32 NVIDIA A100 GPUs.For the RKC integrator we found it critical to scale the energy to be O(1).
Figure 1 shows the computational cost as a function of simulation time for each simulation.We see that RKC outperforms VODE when using Strang and simplified-SDC integration on CPUs.On GPUs, the performance between the RKC and VODE is almost the same, regardless of the integration method.This differs from the findings of Niemeyer & Sung (2014); Stone & Davis (2013) and may be due to the differing costs in evaluating the nuclear reaction right-hand side function or simply because our VODE implementation has been optimized more than RKC.As expected, the GPU runs overall are much faster than the corresponding CPU runs.The right panel of the figure shows the total 12 C mass from burning 4 He.All simulations show the same trend: an initial build up of carbon and then a drop as it is consumed to produce heavier elements.We also see that for this problem, the simplified-SDC method does not improve the efficiency of this simulation, consistent with (Chen et al. 2023).
Our results show that a fully explicit integration method for reactions can be more efficient than an implicit method in modeling XRBs.Finally, we note that a similar study was done looking at He detonations, but there RKC struggles to perform as well as VODE, consistent with findings for more-stiff terrestrial combustion simulations (Niemeyer & Sung 2014).This is likely because the reactions in the detonation are much more stiff and require many more stages, so the cost of evaluating the right-hand side of the reaction ODE system dominates.Future work will explore larger networks as well as higher-order stabilized explicit integrators, like the ROCK family (Abdulle 2002).
Castro is available at http://github.com/AMReX-Astro/Castro.This research used the DOE National Energy Research Scientific Computing Center and was supported by DOE Office of Nuclear Physics grant DE-FG02-87ER40317 and NSF grant PHY-2243856.

Figure 1 .
Figure 1.Number of node hours and carbon mass as a function of time for the different integrators and time-integration method on CPU/GPU.