Exascale models of stellar explosions: Quintessential multi-physics simulation

Abstract

The ExaStar project aims to deliver an efficient, versatile, and portable software ecosystem for multi-physics astrophysics simulations run on exascale machines. The code suite is a component-based multi-physics toolkit, built on the capabilities of current simulation codes (in particular Flash-X and Castro), and based on the massively parallel adaptive mesh refinement framework AMReX. It includes modules for hydrodynamics, advanced radiation transport, thermonuclear kinetics, and nuclear microphysics. The code will reach exascale efficiency by building upon current multi- and many-core packages integrated into an orchestration system that uses a combination of configuration tools, code translators, and a domain-specific asynchronous runtime to manage performance across a range of platform architectures. The target science includes multi-physics simulations of astrophysical explosions (such as supernovae and neutron star mergers) to understand the cosmic origin of the elements and the fundamental physics of matter and neutrinos under extreme conditions.

Keywords

Astrophysics adaptive mesh refinement computational software scientific application multi-physics

1. Introduction

The creation of the chemical elements of which our planet, and ourselves, are made links our existence to the lives and deaths of stars. The abundance patterns of the elements, which have accumulated gradually over nearly 14 billion years of cosmic time, reveal both the physics of stars and of atomic nuclei. Nearly a century after it was suggested that stars could be powered by the transmutation of nuclei, our understanding of this chemical evolution has matured greatly. Yet fundamental questions remain, chief among them: What is the source of the heaviest elements? These questions are the primary focus of the Exascale Computing Project (ECP) ExaStar application development project.

On the basis of fundamental nuclear physics and observed cosmic abundances (Burbidge et al., 1957; Cameron 1957), we know that roughly half of the isotopes heavier than iron were formed via a slow neutron capture process (the s-process) with most of the remaining half due to a rapid neutron capture process (the r-process), and with a scattering of rarer isotopes ascribed to what was originally thought to be a proton capture process (the p-process). We have considerable confidence that the s-process occurs within the stable burning shells of giant stars (Käppeler et al., 2011). In contrast, our incomplete understanding of the r- and p-processes represents one of the long-standing open problems in physics; addressing it requires new experimental measurements of highly radioactive isotopes coupled to mature computational models of the potential astrophysical sites.

Production of the observed pattern of r-process nuclei requires an explosive event that ejects high-entropy and neutron-rich material on a timescale of seconds. The appearance of r-process elements in the spectra of extremely old stars indicates that these elements began forming very early in Galactic history. This suggests that the r-process occurs during the death throes of short-lived, massive stars (Argast et al., 2004). These and other observations (Wasserburg et al., 1996) further suggest that the production site of lighter r-process species (mass number A < 125) is distinct from that of the main r-process (130 < A < 195).

Several astrophysical environments are suspected to produce the conditions needed for r-process nucleosynthesis. The most prominent scenarios are: (1) the collapse of oxygen/neon/magnesium or iron stellar cores, with their associated supernovae; (2) the decompression of nuclear matter ejected during the violent merger of two neutron stars (or a neutron star and a black hole); and (3) the winds from gaseous disks accreting onto black holes formed in failed supernovae or neutron star mergers. The extremes of density, velocities, and temperature realized in these environments also give rise to a rich set of other observables, including high-energy photons, neutrinos, and gravitational waves.

Simulations of each of these scenarios with quantitative physical fidelity rely on similar toolkits: magnetohydrodynamics (MHD), thermonuclear kinetics, the equation of state (EoS) of nuclear matter, radiation transport of neutrinos and photons, and a general relativistic description of gravity. Current models that rely on approximations to all of this physics already tax petascale architectures. To overcome these approximations and to perform definitive simulations at the exascale will require improved radiation transport solver technologies, elliptic and hyperbolic solvers for gravity, and techniques for load balancing multi-physics applications on a dynamically adaptive hierarchical mesh. In addition, a comprehensive suite of new physics implementations is needed, designed to take maximal advantage of hybrid multi-core and many-core architectures.

We are building a software ecosystem for astrophysical simulations that will effectively exploit exascale computational architectures to provide world-leading models of the mechanisms and observable consequences of a variety of stellar cataclysms that lead to explosive nucleosynthesis (Table 1). Our target is a set of simulations including fully self-consistent calculations of all of the main proposed r-process scenarios (core-collapse supernovae, neutron star mergers, and accreting black holes) as well as related stellar explosions: novae, X-ray bursts, and thermonuclear supernovae. Realistic simulations of these scenarios require advanced multi-group neutrino radiation transport (via variable Eddington factor methods) in true multi-dimensions, increasing the computational cost by an order of magnitude over current state-of-the-art calculations. These simulations will also have spatial resolutions via adaptive mesh refinement (AMR) fine enough to resolve the development of the magneto-rotational instability (Balbus and Hawley, 1991), that is, scales of order 100 m within simulations spanning roughly 50,000 km in radius. To fully understand the production of heavy elements in such events, the target simulations should evolve nuclear reaction networks with as many as 3000 isotopes in situ. All of these components are necessary to produce simulations of core-collapse supernovae and neutron star mergers (our primary simulation targets). Of all these components, the neutrino transport is the most computationally intensive and consumes the bulk of the execution time in each simulation. Nevertheless, effective implementations with graphics processing units (GPUs) of the remaining physics must be produced to ensure that none of these intrinsically less-expensive calculations reduce the overall performance of the code. This is most important because to allow a broad range of progenitor scenarios to be explored with the requisite level of physical fidelity, the wall clock time to complete a simulation will need to be reduced from the current level of months to weeks.

Table 1.

A comparison of the typical parameters that characterize the current state of the art for supernova simulation on petascale platforms and the ExaStar targets for our exascale challenge problem. Note that the neutrino phase-space resolution remains constant between the petascale and exascale, as the current resolution has been shown to be adequate for correct results.

Physics component or control parameter	Current production state of the art	ExaStar exascale target
Maximum AMR spatial resolution	$O$ (1 km)	$O$ (10s of meters)
Neutrino transport implementation	Multi-group flux-limited diffusion or variable Eddington factor approaches	Two-moment DG finite elements within an IMEX scheme
Neutrino phase-space resolution	One or two moments, ≈ 20 energy groups	Two moments, ≈ 20 energy groups
Nuclear kinetics	10–20 species	100–3000 species
Typical wall clock time to evolve 1 s of physical time	3–4 months	1–2 weeks
Implemented gravity solver^a	Newtonian gravity with an effective spherical GR potential	Full general relativity with dynamic spacetime coupling

^aThe implementation of a full GR solver is a stretch goal for our development effort and is not discussed in this article.

ExaStar’s approach to attack this challenge includes substantially increasing the capabilities of two widely used stellar astrophysics codes: Flash-X¹ and Castro.² Here, we discuss the development of new physics modules in the context of Flash-X, as we have primarily used Flash-X as the vehicle for development and testing. We have developed the ability to use these physics modules in the Castro code as well in most cases, and anticipate extending this coverage into Castro in the future. We also describe significant changes to the runtime we are instantiating in the Flash-X code, as well as the software engineering practices we use to ensure that these new developments are carefully tracked and tested as they are incorporated into the code base.

2. Code descriptions

2.1. Flash-X

Flash (Dubey et al., 2009; Fryxell et al., 2000) is a modular, multi-physics code with functional physics units that can be applied in a highly configurable way. Flash-X is a new version of the code derived from Flash to be used in the ECP project. Flash-X primarily relies on a finite-volume (FV) discretization, though as we describe in Section 3.2, individual physical operators are not restricted to FV methods. As shown in Figure 1, Flash-X is capable of nearly ideal weak scaling to 10,000 Message Passing Interface (MPI) ranks on the Summit computer at the Oak Ridge Leadership Computing Facility (OLCF). The scalability of Flash-X is crucial to efficiently using the available hardware resources at exascale and beyond. In addition to the MPI grid decomposition provided by the AMR framework, OpenMP threading is already implemented in many of the physics modules in Flash-X, as we describe in more detail in later sections.

2.2. AMReX

In Flash-X, the underlying AMR framework is provided by AMReX. For more details about AMReX, see the accompanying article also in this issue.

AMR is a particularly salient capability of Flash-X that lends itself to the problem of simulating astrophysical phenomena at scale. This adaptivity also provides a convenient block-structured domain decomposition wherein local physics operations can be applied on a per-block basis independent of location in the global domain. Blocks can then be individually assigned to distinct MPI ranks using a work-weighted Morton space-filling curve (see, e.g., Warren and Salmon, 1993) to achieve a similar computational load per process.

3. Physics modules

3.1. Nuclear burning

We use XNet³ (Hix and Thielemann, 1999), a coupled nuclear reaction network to calculate the rate of nuclear energy deposition to the fluid and evolve the detailed composition of the nuclear species responsible for key astronomical observables. Our implementation is generalizable to an arbitrary network and is a significant advance over the existing nuclear burn units in Flash-X and Castro that are limited to a small number of species. These larger networks improve the physical fidelity of simulations by more accurately modeling the spatial and temporal distribution of the energy release, neutronization of the fluid, and mixing of the larger number of species present. These improvements are necessary to build confidence in multi-dimensional models, which will allow more meaningful comparisons with observations. Experience has taught us that a network consisting of ≈160 species is well-suited to the targeted scientific problem of a core-collapse supernova.

3.1.1. Methods

The fully implicit nature of XNet necessitates the choice of a suitable integration scheme. The linear system of ordinary differential equations (ODEs) defining the reaction network is quite stiff and so requires an implicit integration method. This necessitates the use of a dense linear solver that scales $O (N_{s}^{3})$ , where N_s is the number of species in network. Previous work has shown the simple first-order backward Euler method to be very efficient in advancing nuclear abundances within the constraints of the CCSN problem. With this scheme, nuclear abundances, $\vec{Y}$ , are evolved by some change, $Δ \vec{Y}$ , of the system over a timestep, Δt = t⁽ⁿ⁺¹⁾ − t⁽ⁿ⁾, according to

{\vec{Y}}^{(n + 1)} = {\vec{Y}}^{(n)} + Δ \vec{Y}

(1)

This is done using the Newton–Raphson method, based on the Taylor series expansion of ${\vec{Y}}^{(n + 1)} = {\vec{Y}}^{(n)} + f ({\vec{Y}}^{(n + 1)}) Δ t$ , where

\begin{array}{l} f (Y_{i}^{(n + 1)}) \equiv \frac{\partial Y_{i}^{(n + 1)}}{\partial t} \\ = \sum_{j} N_{j}^{i} λ_{j} Y_{j}^{(n + 1)} \\ + \sum_{j, k} N_{j, k}^{i} ρ N_{A} 〈 j, k 〉 Y_{j}^{(n + 1)} Y_{k}^{(n + 1)} \\ + \sum_{j, k, l} N_{j, k, l}^{i} ρ^{2} N_{A}^{2} 〈 j, k, l 〉 Y_{j}^{(n + 1)} Y_{k}^{(n + 1)} Y_{l}^{(n + 1)} \end{array}

(2)

N_A is Avogadro’s number, λ_j is the decay-rate, and ⟨j, k⟩ and ⟨j, k, l⟩ are cross-sections for two- and three-body reactions, respectively. $N_{j}^{i} = N_{i}, N_{j, k}^{i} = N_{i} / \prod_{m = 1}^{n_{j, k}} | N_{m} |$ , and $N_{j, k, l}^{i} = N_{i} / \prod_{m = 1}^{n_{j, k, l}} | N_{m} |$ provide for proper accounting of numbers of nuclei, where N_i specify how many particles of species i are created or destroyed in a reaction, and n_j,k, and n_j,k,l denote the number of different species destroyed in the reactions. This reduces to iteratively solving a linear system $(A \vec{x} = \vec{b})$ in the form

(\frac{I}{Δ t} - J) Δ {\vec{Y}}^{[m]} = - \frac{({\vec{Y}}^{(n + 1), [m]} - {\vec{Y}}^{(n + 1), [0]})}{Δ t} + f ({\vec{Y}}^{(n + 1), [m]})

(3)

{\vec{Y}}^{(n + 1), [m + 1]} = {\vec{Y}}^{(n + 1), [m]} + Δ {\vec{Y}}^{[m]}

(4)

for m = 0, 1, …, with a suitable initial guess

{\vec{Y}}^{(n + 1), [0]}

, where

J = \partial f ({\vec{Y}}^{(n + 1), [m]}) / \partial {\vec{Y}}^{(n + 1), [m]}

is the Jacobian of

f ({\vec{Y}}^{(n + 1), [m]})

and I is the identity matrix. Iteration continues until the solution converges according to nucleon number conservation or a predefined convergence criteria,

{‖ \vec{E} ‖}_{1} < ε

, is satisfied. We typically use a tolerance ɛ = 1 × 10⁻⁴ and define

E_{i} \equiv {\begin{array}{l} 0 & if Y_{j}^{[m]} < Y_{\min} \\ | Δ Y_{i}^{[m]} | / Y_{i}^{[m]} & otherwise \end{array}

(5)

where Y_min = 1 × 10⁻²⁵.

Each Newton–Raphson iteration requires computing the full set of abundance derivatives, calculating all reaction rates, evaluating the Jacobian, evaluating the right-hand side, and then performing one LU decomposition ( $\sim (2 / 3) N^{3}$ floating-point operations) and one back substitution (∼N² floating-point operations). Extremely large variations in timescales between strong, electromagnetic, and weak reactions make astrophysical nuclear reaction networks extraordinarily stiff (Hix and Meyer, 2006). More rigorously, a system is stiff if the ratio of the maximum and minimum real components of the eigenvalues of the Jacobian is much greater than one (values that exceed 10¹⁵ are not uncommon in astrophysics). Practically, this means that the step size is set by considerations of numerical stability rather than accuracy. Complicating the computation itself is indirect memory addressing and loop carried dependencies associated with building the Jacobian. Double-precision floating-point arithmetic is required for the calculation, as the approach to equilibrium at various stages of the burning can lead to the near-cancellation of large reaction fluxes.

Traditionally, a simple backward Euler (BE) integration has proven to be the most computationally efficient choice, so long as the time being evolved by the network is limited to less than 10⁻⁴ seconds such that the benefit of larger step-sizes from higher-order integration methods is not realized since the burning module takes many sub-timesteps (burning timesteps) to evolve through a single hydrodynamic (global) timestep.

Motivated by the results of Longland et al. (2014), we can now use an implicit ODE solver in XNet based on variable-order Backward Differentiation Formula (BDF) methods (Gear, 1971), following the implementations described in Byrne and Hindmarsh (1975) and Brown et al. (1989). The main principle behind the method is to use past behavior of the system to predict the solution for a given timestep. At first order, BDF reduces to BE, but is also capable of using trial error estimates to set the size of the next timestep and adjust the integration order if it would result in a sufficient increase in the size of the step size. Furthermore, the computational cost per timestep of BDF is similar to that of BE, regardless of order, since each step is still dominated by one LU factorization and subsequent back substitution. For a Flash-X simulation, burning would be coupled to other physics and evolved up to the global timestep, typically set by the Courant–Friedrichs–Lewy (CFL) condition and ranging from 10⁻⁷ s to 10⁻⁵ s. For this duration of time, the BDF method is able to sufficiently take advantage of higher-order integration and consequently requires between four to five times fewer timesteps. Since the computational expense per timestep is comparable between the two methods, this results in much better overall performance for the burner. A downstream effect of this substantial decrease in the number of burning timesteps taken in the most expensive cells is an improvement to the load imbalance of the burner.

3.1.2. Parallelization strategy

Within the burning module there is a set of nested do loops used to iterate over the grid zones within all blocks local to an MPI rank. Without any threading, each iteration of this loop structure sends a single zone of thermodynamic (temperature and density) and composition (mass fractions of all species) data to the burner. Inside the burner, these data are used to build a system of ODEs representing the evolution of the N_s species. The solution of this system gives the updated composition and nuclear energy release for the zone. To reduce the time spent computing the nuclear kinetics (and so the overall simulation time), we apply OpenMP directives to this set of loops. Each thread iterates over all zones within a block, “burning” each zone in turn. The number of threads determines the number of blocks iterated over concurrently. For the latter, the zones within each block are distributed among the available threads, so each thread “burns” its portion of the zones sequentially (but in parallel with other threads).

In order to take advantage of GPUs, we primarily rely on optimized batched linear algebra libraries. The batched library calls are used to solve multiple systems of ODEs (i.e., burn multiple zones) in parallel on the GPU. For each model, we initialize our simulations by allocating (for each thread) memory on the host and pinning it to a corresponding memory allocation on the GPU (pinned memory is required to allow asynchronous data transfers). This block of memory is allocated for the data structure that holds the batched thermodynamic and composition data sent to the burner each hydrodynamic timestep which includes the $N_{s}^{2}$ matrix and right-hand side in equation (3) that will be solved in batches.

When simulating stellar explosions, there is often a natural load imbalance imposed by the presence of a burning front moving through the stellar material. This process of load balancing is particularly important for nuclear burning because of there is a high degree of variability in the computational cost that is very sensitive to initial conditions. However, the mechanism we use to weight the space-filling curve is generalizable to other physics modules and provided by the AMR framework. Burning that occurs within the front takes many more sub-timesteps than burning that occurs behind it or ahead of it (where no burning occurs). Therefore, if one MPI rank is assigned several blocks where a significant amount of burning takes place and another MPI rank is assigned only blocks where no burning occurs, then a load imbalance will be present. To improve this load balancing, we weight each block based on the average number of burning timesteps taken by any zone within it. By doing so, MPI ranks assigned blocks with many burning timesteps are given fewer total blocks to compute.

For large numbers of MPI ranks, we must balance improved homogeneity in the burning cost with increased communication costs related to guardcell filling (the process of getting the values for halo cells surrounding computational domain in stencil computations) required by hydrodynamics and shock detection. By heavily weighting the block redistribution using the total sum of sub-timesteps required per block, we can reduce the load imbalance of the burning even further, but this results in a much more expensive simulation due to the additional incurred communication. On Summit, we find weighting each block by the average number of burning timesteps required per batch (i.e., the maximum number of timesteps for any cell within the batch) results in better overall load balance than other measures.

3.1.3. Example simulation

With the goal of understanding the observable impacts of the central engine and convective mixing on CCSN ejecta, we have taken a 3D model of the central engine for a 9.6 M_⊙ progenitor star evolved by the chimera code (Bruenn et al., 2020) as initial conditions to a new simulation in Flash-X that follows the progression of the explosion throughout the entire star. Using self-consistent chimera models, rather than parameterized models, provides the most physically faithful starting point currently available. Flash-X gives us the unprecedented ability to track 160 nuclear species from 1H to 64Ge to shock breakout.

In addition to hydrodynamics, self-gravity was included via Flash-X’s Newtonian multi-pole solver, which solves the Poisson equation through a multi-pole expansion. As we removed the proto-neutron star (PNS) in this simulation, Flash-X’s Helmholtz EoS was able to be used, which maintains thermodynamic consistency and includes contributions to internal energy from ions, electrons, positrons, and radiation. The inclusion of XNet in Flash-X allows us to track the composition of 160 nuclear species throughout the evolution of the explosion, which leads to a more accurate analysis of the ejecta seen at shock breakout and has traditionally not been included in these types of simulations. It is this ejecta that is injected behind the shock and propagates as “bullets” throughout the star (see Figure 2). Because these bullets exit the star first, the expanded network in XNet allows us to more accurately compare to astronomical observations of CCSNe.

3.2. Radiation transport

The toolkit for high-order neutrino radiation hydrodynamics⁴ (thornado) (Chu et al., 2019a, 2019b; Endeve et al., 2019) is being independently developed for modeling CCSNe. For ExaStar, we are interested in the capability of thornado to evolve the neutrino radiation field with spectral neutrino transport using a two-moment model, where moments of the neutrino phase-space distribution function—representing spectral energy and momentum densities—are evolved and where higher-order moments (e.g., components of the radiation pressure tensor) are approximated using algebraic expressions involving the lower-order, evolved moments to form a closed system of equations. thornado is not developed as a supernova simulation code with all the utilities needed to run at large scale, but rather as a collection of modules that can be incorporated into such simulation codes (e.g., Flash-X and Castro). For this reason, a primary focus of the thornado development is node-level performance on heterogeneous computing systems.

3.2.1. Methods

The neutrino radiation hydrodynamics equations are discretized in phase-space using a combination of FV; Leveque (2002)) and discontinuous Galerkin (DG; Hesthaven and Warburton (2008)) methods. Specifically, the convection and gravitational source terms in the hydrodynamics equations are discretized with FV methods. The semi-discrete system evolved with FV methods can then be written as a system of ODEs

d_{t} \bar{u} = T_{u} (\bar{u}, Φ)

(6)

where

\bar{u}

represents the cell averages of the solution vector u across the computational domain, the right-hand side represents the combined FV representation of convection and gravitational terms, and Φ is the Newtonian gravitational potential. The moment equations—including the coupling between the fluid and the radiation field—are discretized with DG methods (see Chu et al. (2019a, 2019b); Laiu et al. (2020, 2021) for details). The semi-discrete system evolved with DG methods can be written as the following systems of ODEs

d_{t} u = C_{u} (M, u)

(7)

\begin{array}{l} d_{t} M = T_{M} (M) + C_{M} (M, u) \end{array}

(8)

where

u

and

M

, respectively, represent all the fluid and radiation degrees of freedom evolved with the DG method. The right-hand side of equation (7) corresponds to the DG discretization of the sources and sinks due to neutrino-matter interactions, and the transport operator

T_{M}

and the collision operator

C_{M}

on the right-hand side of equation (8) correspond to streaming terms and neutrino-matter interaction terms appearing in the moment equations, respectively. The use of DG methods to solve the moment equations for neutrino transport equations is, in part, motivated by the fact that DG methods recover the correct asymptotic behavior in the so-called diffusion limit, which is characterized by frequent scatterings (Larsen and Morel, 1989; Adams, 2001).

We pause here to point out the relative size of the solution vectors $u$ and $M$ in the DG method in equations (7) and (8). With the approximation space for the DG method constructed by tensorization of one-dimensional polynomials of degree k and the computational domain covered by N_x d_x-dimensional elements, we have $u \in ℝ^{N_{u}}$ , where $N_{u} = 6 \times {(k + 1)}^{d_{x}} \times N_{x}$ . Similarly, using N_ɛ elements to cover the energy dimension of momentum space and N_Sp neutrino species, we have $M \in ℝ^{N_{M}}$ , where $N_{M} = 4 \times N_{Sp} \times {(k + 1)}^{d_{x} + 1} \times N_{x} \times N_{ε}$ . In our target application we have N_Sp = 6 and N_ɛ ∈ [10, 20], so that $N_{M} ≫ N_{u}$ .

Since neutrino-matter interactions are frequent in high-density regions of the computational domain, and thus introduce stiffness, we use implicit–explicit (IMEX; Ascher et al. (1997); Pareschi and Russo (2005)) time integration to stably integrate with timesteps that can greatly exceed the shortest timescales associated with neutrino-matter coupling using implicit methods, while we integrate the convection and transport operators ( $T_{u}$ and $T_{M}$ ) using explicit methods. The first-order accurate (in time) instantiation of this IMEX scheme, evolving equations (6)–(8) from tⁿ to tⁿ⁺¹ = tⁿ + Δt, can be written as

\begin{array}{c} {\bar{u}}^{*} = {\bar{u}}^{n} + Δ t T_{u} ({\bar{u}}^{n}, Φ^{n}) \end{array}

(9)

\begin{array}{l} M^{*} = M^{n} + Δ t T_{M} (M^{n}) \end{array}

(10)

\begin{array}{l} u^{n + 1} = u^{*} + Δ t C_{u} (M^{n + 1}, u^{n + 1}) \end{array}

(11)

\begin{array}{l} M^{n + 1} = M^{*} + Δ t C_{M} (M^{n + 1}, u^{n + 1}) \end{array}

(12)

The extension to higher-order temporal accuracy with multi-stage, diagonally implicit Runge–Kutta methods (Pareschi and Russo, 2005) is relatively straightforward, but is omitted here for brevity. Since the radiation transport timestep is smaller and independent of the global timestep (i.e., the timestep used by the main Flash-X evolution loop), we subcycle in time the radiation update in Flash-X until reaching the desired timestep, updating boundary cell information between each subcycle iteration.

3.2.2. Implementation

Notice that the implicit solve modeling neutrino-matter interactions in the coupled system in equations (11) and (12) requires a mapping (reconstruction) of the fluid variables from the FV representation ${\bar{u}}^{*}$ to the DG representation $u^{*}$ . A corresponding inverse map (projection) from $u^{n + 1}$ to ${\bar{u}}^{n + 1}$ is also needed to initiate the next timestep. To this end, we borrow ideas from Dumbser et al. (2014), and consider DG elements as composed of FV subcells. Specifically, consider an element $K \subset ℝ^{d_{x}}$ . On K , the DG representation is denoted $u_{h} (x, t) \in V_{h}^{k} (K)$ , where the approximation space $V_{h}^{k} (K)$ is constructed from the tensor product of one-dimensional polynomials of maximal degree k

u_{h} (x, t) = \sum_{j = 1}^{N_{K}} u_{j} (t) ψ_{j} (x)

(13)

where

ψ_{j} \in V_{h}^{k} (K)

and

N_{K} = {(k + 1)}^{d_{x}}

is the number of degrees of freedom per spatial element. Note that the expansion coefficients

u_{j}

from all elements comprise the solution vector

u

in equations (11) and (12). Next, the element K is divided into a subgrid S of N_K nonoverlapping FV cells S_j,

S = \cup_{j = 1}^{N_{K}} S_{j}

, so that K \ S =∅; see Figure 3 for an example in two spatial dimensions (d = 2) and k = 1. On the subgrid, the representation is given by piecewise constants (cell averages)

u_{h} (x, t) = \sum_{j = 1}^{N_{K}} χ (S_{j}) {\bar{u}}_{j} (t)

(14)

where

{\bar{u}}_{j} (t)

is the cell average in subgrid cell S_j, and χ(S_j) is the indicator function on S_j. Thus, the expansion coefficients

{\bar{u}}_{j}

in equation (14), collected across all elements K , comprise the solution vector

\bar{u}

in equation (9).

Then, knowing the cell averages ${\bar{u}}_{j}$ after the explicit step with the FV method in equation (9) (we suppress superscripts indicating time from here on), we reconstruct the DG representation $u_{h}$ by requiring that

\frac{1}{| S_{i} |} \int_{S_{i}} u_{h} (x) d x = \frac{1}{| S_{i} |} \int_{S_{i}} u_{h} (x) d x \equiv {\bar{u}}_{i}

(15)

∀ S_i ∈ S . Here |S_i| is the volume of subcell S_i. With the definition in equation (13), this gives rise to the linear system

R^{- 1} u_{K} = {\bar{u}}_{K}

(16)

where

u_{K} = {u_{j}}_{j = 1}^{N_{K}}

{\bar{u}}_{K} = {{\bar{u}}_{j}}_{j = 1}^{N_{K}}

, and the components of the N_K × N_K inverse of the reconstruction matrix R are

R_{i j}^{- 1} = \frac{1}{| S_{i} |} \int_{S_{i}} ψ_{j} (x) d x

(17)

The reconstruction step then amounts to inverting to find

u_{K} = R {\bar{u}}_{K}

(18)

After obtaining the DG representation $u$ from the coupled implicit update in equations (11) and (12), we can easily compute the cell averages for the FV representation on the subgrid in each element

{\bar{u}}_{i} = \frac{1}{| S_{i} |} \int_{S_{i}} u_{h} (x) d x = \frac{1}{| S_{i} |} \sum_{j = 1}^{N_{K}} \int_{S_{i}} ψ_{j} (x) d x u_{j}

(19)

∀ S_i ∈ S , which can be written as

{\bar{u}}_{K} = P u_{K}

(20)

where the components of the N_K × N_K projection matrix P are

P_{i j} = \frac{1}{| S_{i} |} \int_{S_{i}} ψ_{j} (x) d x

(21)

From the definitions above, it is obvious that RP = I, that is, we do not lose any accuracy when switching between DG and FV representations.

Due to the nature of neutrino-matter interactions, the coupled implicit system given by equations (11) and (12) is completely local in position space. However, inelastic scattering and other neutrino interactions result in global coupling in momentum space as well as across neutrino species. Moreover, since neutrino opacities and associated interactions kernels depend nonlinearly on local thermodynamic conditions, and because of the presence of blocking factors inside integral collision operators to account for Pauli’s exclusion principle, which applies to neutrinos (fermions), an iterative procedure is needed for the implicit step. On the one hand, nonlinearity and global coupling in momentum space make the implicit solve one of the most compute intensive aspects of neutrino radiation hydrodynamics simulations. On the other hand, spatial locality makes this part of the solve more easily parallelizable.

We have developed and tested various strategies for solving this nonlinear system. With a slight abuse of notation, we let $u$ represent the fluid fields in a single spatial point. Similarly, we let $U$ denote all the radiation moments associated with a single spatial point, that is, $U \in ℝ^{4 \times N_{Sp} \times (k + 1) \times N_{ε}}$ . We let $U = {(u, U)}^{T}$ denote the total solution vector. Then, the solution to equations (11) and (12) can be formulated as a fixed-point problem (Laiu et al., 2020)

U = G (U) = (\begin{matrix} g (U) \\ U (u, U) \end{matrix})

(22)

where

g

and

U

are fixed-point maps for

u

and

U

, respectively. We first consider solving equation (22) using the coupled fixed-point method, which solves for the total solution vector U by taking

U^{[k + 1]} = G (U^{[k]}), k = 1, \dots

(23)

where the initial guess U^[1] can be taken as the output from the explicit step in equations (9) and (10). In practice, we employ Anderson acceleration (Anderson, 1965; Walker and Ni, 2011) to improve the convergence of fixed-point iteration. Benefits of using fixed-point iteration over Newton’s method on the coupled problem in equation (22) are (i) it does not require the construction of the Jacobian matrix, which is only available in an approximate sense due to tabulated opacities, and (ii) it does not require inversion of a dense linear system in each iteration. While the fixed-point method requires more iterations than Newton’s method to converge, we generally observe that the fixed-point algorithm is faster than Newton’s method for solving equation (22) (Laiu et al., 2021).

The coupled fixed-point method requires evaluation of the neutrino opacities in each iteration (i.e., each time $u$ is updated). Neutrino opacities are obtained from table interpolations, which we have found to constitute a computational bottleneck. Tabulated neutrino opacities and supporting infrastructure routines (e.g., interpolation and I/O) are provided by the WeakLib⁵ library, which we discuss in more detail in Section 3.4 in the context of the EoS. To reduce opacity evaluations, we have developed the second method, a nested fixed-point iteration scheme with two layers, which can be formulated as (Laiu et al., 2020)

\begin{array}{l} u = g (\hat{U} (u)) \end{array}

(24)

\hat{U} (u) = U (u, \hat{U} (u))

(25)

where the outer layer, equation (24), is a fixed-point problem on the matter states

u

, and the inner layer is a fixed-point problem on the radiation moments

\hat{U}

for fixed

u

. The iterative scheme for this nested approach then becomes

u^{[k + 1]} = g (\hat{U} (u^{[k]})), k = 1, \dots,

(26)

where

\hat{U} (u^{[k]}) = U^{[k, *]}

is the limit point of the inner iteration sequence

U^{[k, ℓ + 1]} = U (u^{[k]}, U^{[k, ℓ]}), ℓ = 1, \dots

(27)

We use Anderson acceleration also for this nested scheme to accelerate both the inner and outer iteration sequences separately.

As a third option, we have developed an alternative to the nested fixed-point algorithm, where we have replaced the inner fixed-point iteration sequence with a Newton iteration sequence. This algorithm can be formulated as (Laiu et al., 2021)

u = g (\hat{U} (u))

(28)

ℱ (u, \hat{U} (u)) = 0

(29)

where the outer layer is equivalent to the nested fixed-point method in equation (24), while the inner layer has been reformulated as a root-finding problem on the radiation moments

\hat{U}

for fixed fluid state

u

. Here, the formulation of the outer iteration is identical to equation (26), but

\hat{U} (u^{[k]}) = U^{[k, *]}

is instead the limit point of the Newton iteration sequence

U^{[k, ℓ + 1]} = U^{[k, ℓ]} + δ U^{[k, ℓ]}, ℓ = 1, \dots,

(30)

where

δ U^{[k, ℓ]} = - {[\frac{\partial ℱ}{\partial U} (u^{[k]}, U^{[k, ℓ]})]}^{- 1} ℱ (u^{[k]}, U^{[k, ℓ]})

(31)

Since the matter state is kept fixed in equation (30), the Jacobian matrix $\partial ℱ / \partial U$ can be computed analytically, without opacity interpolations. The motivation for replacing the inner layer is to take advantage of the faster convergence rate of Newton’s method, at the expense of more dense linear algebra.

3.2.3. Example core-collapse supernova simulation

To demonstrate the coupling between the FV hydrodynamics in Flash-X and the DG neutrino transport in thornado, we have performed a core-collapse supernova simulation in spherical symmetry, initiated from a 15 M_⊙ progenitor from Woosley and Weaver (1995). This simulation evolved electron neutrinos and antineutrinos (N_Sp = 2), using the SFHo EoS (Steiner et al., 2013) and the neutrino opacity set from Bruenn (1985), and we used the nested solver in equations (24) and (25) to couple the neutrinos with matter.

Figure 4 shows the evolution of the electron fraction versus radius during a critical epoch in the simulation; namely from 15 ms before to 100 ms after core bounce. (The time of bounce, t_b, is defined as the time when the central density reaches its maximum value.) The simulation exhibits the characteristic features of spherically symmetric supernova simulations (e.g., Liebendörfer et al., 2001). During the collapse phase (t < t_b; blueish lines), the core electron fraction is reduced due to electron capture processes. Once the central density reaches sufficiently high values (>10¹² g cm⁻³), neutrinos become trapped in the fluid, and the electron fraction remains relatively unchanged. After bounce and shock formation (t > t_b; reddish lines), once the shock reaches low enough densities, neutrinos can again escape from the core—producing a “burst” of electron neutrinos. This results in a reduction in the electron fraction below the shock. At t = 100 ms, the initial bounce shock has stalled at about 100 km from the center of the star. We note that this model was computed without including relativistic observer corrections in the neutrino transport equations. This results in a lower core electron fraction than is seen in relativistic models (Lentz et al., 2012). We are currently working to extend our neutrino transport scheme to include these relativistic effects.

3.2.4. Parallelization strategy

We rely on a combination of preprocessor macros, GPU compiler directives (OpenMP or OpenACC), and generic interfaces (see, e.g., Listing 1) to optimized linear algebra libraries (e.g., MAGMA; Tomov et al., 2010) to fully accelerate all components of the DG-IMEX scheme in thornado. For example, consider the code in Listing 2 that performs an array permutation during the explicit update or the opacity interpolation kernel in Listing 3, needed in the implicit part of the scheme described in Section 3.2.1. With this approach, depending on the available hardware and programming environment, we are able to choose at compile time whether to use OpenMP with GPU offload, OpenACC, or more traditional OpenMP for CPU multi-core to execute the nested loop in parallel.

Listing 1.

Example of interface for matrix–matrix multiplication to be performed by an optimized linear algebra library. The TARGET attribute is added to the dummy arguments so the C address can be ascertained and passed to C routines with the ISO_C_BINDING interface.

Listing 2.

Example of array permutation using tightly nested loops that are common throughout thornado.

Listing 3.

Example of opacity interpolation kernel.

Due to the large number of degrees of freedom evolved, radiation transport is usually the computational bottleneck in astrophysical simulations that include it. For this reason, it is important to have a highly optimized implementation. The stand-alone nature of thornado makes these fine-grained optimizations much simpler to test and later integrate into the larger ExaStar ecosystem. For example, in Figure 5, we show some preliminary profiling results for a CCSN-like benchmark problem run with thornado on a 3D grid, where we have adapted the benchmark to be more representative of workloads that one might encounter in a large-scale simulation with Flash-X. Specifically, we solve the problem in a 3D box on a Cartesian mesh with 12 elements in each dimension; this is a fairly good representation of a single AMR block in a Flash-X simulation that would be local to one MPI rank. This benchmark was run on a single node of the Summit computer at the OLCF. Each node is composed of two IBM POWER9 CPUs and six NVIDIA Volta GPUs, but we limit our tests to a single CPU core and one GPU, since this would be representative of the resources available to a single MPI rank in a large-scale simulation.

In addition to showing a significant overall speedup of OpenACC relative to using only the CPU, we were able to identify that the optimal choice of nonlinear solver depends on the details of the hardware architecture. Another example of how we can use the stand-alone thornado results to guide our optimization efforts is demonstrated in Figure 6. Here, we were able to discover that the implementation of OpenMP offloading by the IBM XL compiler was about three times slower than the OpenACC implementation by PGI. Preliminary analysis suggests this may be accounted for by different choices for the GPU kernel launch configurations (i.e., fewer thread blocks and more registers), but this warrants more study before we can draw any definitive conclusions.

3.3. Magnetohydrodynamics

The nuclear kinetics and neutrino transport modules described above are coupled to an accurate solver for the dynamics of the plasma. In the ExaStar application, magnetohydrodynamics (MHD) is handled by the Spark code being developed independently in Flash-X framework to address a broad range of astrophysical MHD problems. Below we describe the general method implemented in Spark for solving the MHD equations and our approach for parallelization. For a detailed description of Spark, we refer the reader to Couch (2020).

3.3.1. Solution approach

For the ExaStar application, Spark solves the equations of ideal MHD

\begin{array}{l} \frac{\partial ρ}{\partial t} + \nabla \cdot (ρ v) = 0 \\ \frac{\partial (ρ v)}{\partial t} + \nabla \cdot [ρ v v^{T} - B B^{T} + I (P + \frac{B^{2}}{2})] = - ρ \nabla Φ \\ \frac{\partial E}{\partial t} + \nabla \cdot [(E + P + \frac{B^{2}}{2}) v - (v \cdot B) B] = - ρ v \cdot \nabla Φ \end{array}

(32)

where ρ is the mass density, v is the velocity vector, B is the magnetic field, P is the pressure, Φ is the gravitational potential, and E is the total specific energy. The total specific energy includes contributions from kinetic motion, thermal motion, and magnetic fields:

E = ρ e + (1 / 2) (ρ v^{2} + B^{2})

The Spark MHD solver (Couch, 2017) implements the cell-centered method of generalized Lagrangian multipliers (GLM; Dedner et al., 2002; Mignone et al., 2010) to control the growth of the divergence of the magnetic field. The GLM–MHD scheme employs hyperbolic advection and parabolic damping of divergence errors in order to avoid expensive elliptic divergence cleaning (e.g., Jiang and Wu, 1999) or complicated staggered-mesh constrained transport (e.g., Gardiner and Stone, 2005; Lee et al., 2009; Lee, 2013). Spark implements the GLM–MHD scheme via a FV approach using high-order primitive reconstruction, multiple Riemann solvers for flux calculation, and method-of-lines time integration using multi-stage strong stability preserving (SSP) low-storage Runge–Kutta (RK) integrators of second- and third-order (e.g., Gottlieb and Shu, 1998). For realistic, non-Gamma-law equations of state, such as are required for the ExaStar application, Spark avoids the approach of Colella and Glaz (1985) used in other Flash-X solvers and adopts the volumetric internal energy, ρe, as an auxiliary thermodynamic primitive variable instead, as in Almgren et al. (2010). We have found this approach to be generally more robust but either method should be mathematically equivalent (e.g., Zingale and Katz, 2015). By default, Spark uses fifth-order WENO reconstruction (Borges et al., 2008), second-order time integration, and the Harten-Lax-van-Leer-Discontinuities (HLLD) approximate Riemann solver which includes MHD waves (Miyoshi and Kusano, 2005).

3.3.2. Parallelization strategy

Parallelization of Spark is achieved by two levels of domain decomposition realized via a hybrid MPI+OpenMP approach. AMR grid patches are distributed among MPI ranks, and communication of “halo” or “guardcell” data necessary to complete the solution update described above is handled by the AMR package (i.e., AMReX). The default two-stage Runge–Kutta time integration in Spark nominally requires an update of the guardcell data prior to each stage, but we avoid this in Spark by doubling the size of the halo region around each AMR patch and communicating all the required data for a full timestep update at once. This approach requires duplicated computation to update some of the guardcell data for each patch, but avoids a global synchronization step for the intermediate Runge–Kutta stage and the need for a global storage space for intermediate solution data. The net benefit of this approach is to a large extent dependent on the size of the AMR patch, with larger patches better able to amortize the cost of the additional communication and calculation.

Within and AMR patch, OpenMP is used to achieve thread-based, shared-memory parallelism. While the overall hyperbolic update is completed in a directionally unsplit fashion for each Runge–Kutta stage, the calculation of the fluxes required for the update is carried out dimension-by-dimension. In Spark, each thread stores a one-dimensional auxiliary array to hold the data necessary to compute the fluxes for a spatially contiguous line of zones in each direction. These arrays are sized so as to fit in the outermost level of cache so that cache reuse is maximized during the calculation of the fluxes. The final hyperbolic update during each stage is also threaded using OpenMP. We have also modified Spark to use OpenMP offloading features to accelerators. This is done within the AMR patch for both the flux calculation and the hyperbolic update steps.

3.4. Equation of state

The MHD equations are closed with the addition of an EoS, which relates, for instance, the pressure of the fluid to its density and temperature. A polytropic EoS can be used to describe simple gases

P = (γ - 1) ρ ϵ

(33)

where γ is the ratio of specific heats. However, in the problems for which the ExaStar code base is being developed, the use of a simple EOS is insufficient. In these dense astrophysical environments, the pressure is dominated by the contribution of degenerate fermions (electrons and neutrons), so more complicated expressions must be used. Because the EoS is called many times during stellar explosion problems, the time spent in the EoS routines can become large.

These non-trivial EOSes typically require integrals of transcendental functions in evaluating the thermodynamic quantities of interest, and routines for these calculations typically are designed for accuracy and thermodynamic consistency at the expense of speed. A common way to reduce this time is by creating a table of desired quantities using the full, detailed EoS and then using a scheme that interpolates thermodynamic quantities from the table. The Helmholtz EoS FLASH module (suitable for electron-degenerate stellar environments) uses a biquintic Hermite interpolant over a table of the Helmholtz free energy

F = ϵ - T S

(34)

d F = - S d T + \frac{P}{ρ^{2}} d ρ

(35)

All partial derivatives of the pressure, entropy, and internal energy are formed analytically. Searches through the free energy table are avoided by computing hash indices from the values of any given $(T, ρ, \bar{Z} / \bar{A})$ triple (where $\bar{Z}$ and $\bar{A}$ are the average atomic number and the average atomic mass, respectively). The high order of the interpolation is generally necessary to ensure accuracy and thermodynamic consistency and is significantly more expensive than a lower-order method. However, by using GPUs to accelerate this scheme, we have been able to achieve considerable performance gains (Messer and Papatheodore, 2018).

The neutron-degenerate environment found in the PNS in core-collapse supernovae is an even more complex situation, not the least because matter undergoes a first-order phase transition at very high densities (several times 10¹⁴ g/cm³) from a mix of nuclei and nucleons to uniform nuclear matter. Deploying a more realistic EoS also has a non-trivial computational cost. For this reason, we use the tabulated equations of state provided by the WeakLib library. In addition to the tabulations themselves, WeakLib includes useful auxiliary routines for manipulation of the EoS (e.g., interpolation, numerical differentiation, and I/O). Each EoS tabulation has an accompanying set of neutrino opacities with similar supporting auxiliary functions. The interpolation and differentiation routines have all been ported to be callable from inside OpenMP target offload and OpenACC parallel regions, or they can also be invoked to launch independent GPU kernels to operate on a vector of thermodynamic states. This method of batching states together has proven most effective in thornado benchmarks (see Section 3.2.4), but requires some code restructuring to avoid point-wise EoS calls.

3.5. Gravity

For self-gravitating systems, the need arises to solve in a distributed manner the underlying elliptic equations governing the determination and evolution of the gravitational potential or its relativistic counterpart. In the Newtonian case, a solution of the three-dimensional, elliptic Poisson equation gives the gravitational potential, and extensions to more general relativistic treatments of gravity can be made in stages. For centrally concentrated problems (like CCSNe), FLASH uses the multi-pole expansion of the potential which is suitable for mass distributions dominated by high degrees of spherical symmetry. Broadly speaking, the calculation can be divided into three parts: determining the center of mass, computing the multi-pole moments, and finding the final gravitational potential. The center of mass and multi-pole moment calculations both involve local computation followed by an MPI collective. The multi-pole moment collective is more expensive by orders of magnitude due to the larger amount of data, so it is the primary target of our optimization efforts to this point. Our approach involves splitting the multi-pole calculation into two components: the first performs the local calculation to obtain the moments and makes a non-blocking collective call. Other point-wise physics can then be performed (e.g., local transport or nuclear burning). Once these local calculations have finished, the collective will have made progress. Once the communication is finished, each rank computes the gravitational potential at each point. We continue to fine-tune this approach, and the orchestration system (OS) described in Section 4 will significantly enhance its efficacy.

4. Orchestration

In high performance computing (HPC), two schools of thoughts have existed for almost the entire history of the field regarding how to deal with the orthogonal axes of complexity: the target system and the domain science. One school has been represented by the science code developers who have taken on both axes of complexity in order to avoid hard dependencies on third-party tools. The other school is represented by the abstraction and tool developers that seek to hide the platform complexity from the scientific components of the software. Historically, the first approach has been dominant in the absence of a strong incentive to consider high-level abstractions in the code. This has changed with increases in platform heterogeneity, and several abstraction tools have come into existence with varying degrees of success (Bianco and Benedicic, 2020; Edwards and Sunderland, 2012; Hornung and Keasler, 2020). Opportunities also exist to design between these two extremes a domain-specific performance portability system that is suitable for tuning codes to platforms without having to program at low levels of the code for every specific platform. The ExaStar team is designing one such system that will allow Flash-X to make good, efficient use of resources on heterogeneous platforms, including those where host-to-device data movement is necessary (O’Neal et al., 2020). This OS is composed of three sets of tools: the configuration tools (CFT), code translators (CT), and the orchestration runtime subsystem (OR). The CFT are used to configure an instance of an application to be run on a specific platform and to communicate to the OR where data movements will be needed between devices. The maintained source code is augmented with performance hints in the form of Flash-X native directives. Additionally, we have developed a limited DSL for expressing recipes for overall control flow of the simulation. In Flash, this functionality was provided by the driver which invoked all physics operators in a predetermined order. If one wanted a different order one would provide an alternative implementation of the driver. Several alternative implementations have been included in Flash distribution. We have introduced the recipes because control flow now may involve offloading work and data to different devices, therefore, simple drivers are no longer adequate. The CThave two tasks: (1) parse the Flash-X directives for optimization within the physics operators; and (2) to translate the recipe into a data flow description of the execution suitable for the OR. The OR then uses this generated information to orchestrate computation and data movement between devices as needed. The overarching design principle of the OS is to enable both domain-specific knowledge and platform knowledge to be utilized in application configuration without placing undue burden on either the tool developers or the domain experts (Figure 7).

The OS takes the separation of concerns to its logical conclusion. It exposes the hierarchy of data motion and allocations deep within the components of the applications and separates them from the arithmetic and logic implementing the numerical method. On the platform hardware side, they expose the hierarchical parallelism of a given platform’s hardware through appropriately designed interfaces. The objective of this exercise is to enable the software to take maximal advantage of the exposed parallelism in hardware without having to alter the arithmetic. The directives serve the purpose of letting the CTs know which kind of transformations are safe for specific devices. Some program synthesis is involved in the process that can take several passes of code generation. Similar to the code expression, the offline tools themselves subscribe to separation of concerns in the sense that each tool addresses itself to a very focused program synthesis step.

The configuration toolchain takes its inspiration from the setup tool of Flash (Dubey et al., 2009) which has been used to implement composability at the level of major functionalities in the code. For example, an application could be configured to use plain hydrodynamics with a uniformly discretized mesh and constant gravity, while another application may use magnetohydrodynamics with an adaptive mesh, self-gravity, and nuclear burning. The setup tool uses a configuration domain-specific language (DSL) to encode metadata about each component and subcomponent of the code in a completely distributed fashion in config files. The metadata about one component is completely oblivious of the metadata about other components, though it has some knowledge of the components that are available in the code. The setup tool recursively parses the config file of each mentioned component, assimilates the information, and assembles a consistent, fully configured instance of an application along with its build system and static runtime environment. The setup tool is a modestly sized Python script that also has the ability to dictate the order of indices in array references. The meta-information encoded in config files includes state variables needed by the components; of these, the state variables that need resolution at fine-coarse boundaries; other components that are necessary for this component to work correctly; other components that should not be included in the same application configuration; runtime parameters; and many more.

We begin by casting every non-trivial function in the physics operators of the code as a collection of code blocks. Some code blocks may be declarations, some may implement the control logic, and some will implement the numerics of the function. Sometimes arithmetic and logic blocks cannot be separated out; some code blocks have both. These code blocks become components in a hierarchical composability through the use of key-value dictionary where values are code snippets of arbitrary length and complexity. Keys are user defined with a provision for multiple alternative definitions which lets them mimic the template meta-programming in C++ where a single expression of an algorithm can have specializations through alternative definitions of the keys. For more detailed description of key-value feature of the CFT, see Couch et al. (2021).

The information encoded in the directives is utilized by the CTto determine which code blocks can be treated as kernels on devices such as accelerators, which code blocks can have overlapped computation with other code blocks, which code blocks can be queued sequentially so that data movement is minimized, etc. Note that the absence of such metadata has no implication for correct execution. Instead, what this decomposition aims for is better optimization with richer metadata. While the present implementation primarily requires hand coding of this information, we expect to be able to utilize third-party code analysis tools to infer this information automatically in the future. This division of labor between multiple tools is made not only to achieve good encapsulation and modularity but also to ensure that no single tool is too complex or requires too much intelligence. Our design approach enables extraction of all useful information for orchestrating computations during runtime at the time of application configuration. This is in contrast to general solutions that delay the orchestration decision due to the conservative assumption that information about the tasks and their dependencies is only fully known at runtime. While this approach may not squeeze every last bit of performance from the hardware, it has the virtue of being simple and maintainable with no runtime overheads. Also, by design, every step in program synthesis is human-readable to facilitate both correctness-debugging and performance-debugging.

The design goal of the OR is to develop a base set of runtime elements, such as thread teams, that can be composed at runtime into thread team configurations in ways that maximize the efficient use of the node. For a given simulation, the OR operates with N thread teams that are created when the simulation starts and persist until the simulation terminates, where each team is allowed to simultaneously use at most M threads at any given time. The thread teams are run in cycles such that for each cycle a team is assigned an action routine and a data type. Examples of data types are tiles, blocks, and data packets of blocks. Upon starting a cycle, the team is given data items of the appropriate type, and the threads of the team work in a coordinated fashion so that the given action is applied to each data item once. The cycle ends once the team has been informed that no more data items will be given and the action has been applied to all given data items. In our design, the pairing of the action with data items is viewed as a task so that thread teams implement task-based parallelism via thread-level parallelism. While the threads in the team execute code in the host, they can also be used to launch kernels in accelerators and therefore make use of finer-grained parallelism. Ideally, the data type assigned to a thread team for a given action will be chosen based on the hardware that the action routine has been written to use. For instance, a tile would be a sensible data type for executing an action routine that uses the CPU for heavy computation or a data packet of blocks for a routine that launches kernels on a remote device with its own memory system. While it might have been possible to adapt existing tools such as Legion for some of the functionality we are aiming for, we have opted for an end-to-end domain-specific solution that has simple enough tools that even a small team can customize and maintain with their own code. Additionally, our OS is absolutely language agnostic, even though we are exercising it with Fortran.

Note that no part of our design limits itself to using manually built task graphs nor requires extreme simplicity in them. The interfaces are designed to be robust enough that the OR can also execute complex task graphs. By disassociating the executor of the graph from the generator, we have ensured that in the future, we will be able to use automation in task graph generation without having to alter the mechanics of the OR. All the heavy lifting of analysis and scheduling can be done separately in the generator which never needs to interface with the executor. We believe that this model of combination of program synthesis with orthogonal composition concepts will have longevity because of its ability to adapt to new systems incrementally.

5. Software engineering

This section highlights the design, development, testing, verification, and validation practices of the Flash project. For more details about the evolution of the code and its practices see Dubey et al. (2013a, 2013b, 2016).

5.1. Design

The highest level software architecture of the code recognizes that there are two types of capabilities in a multi-physics scientific code that have fundamentally different functionalities and lifecycles. The infrastructure and the book-keeping parts of the code need to be more stable and robustly designed to provide a strong backbone for the code. These functionalities cover discretization methods, coarse-grain data movement between nodes, I/O, runtime environment, and orchestration of physics solvers. These components are also at the heart of any extensibility in the physics capabilities of the code. The second type of capabilities is the physics and numerical solvers. These can be made to act as clients to the infrastructure without any need to know the details of infrastructure implementation. This approach is necessary because these components are likely to be subject of research and more likely to undergo non-trivial changes during the lifetime of the code. Therefore, they should be able to localize their modifications without impacting other parts of the code. These design principles are firmly encoded in Flash’s architecture (Dubey et al., 2009; 2019) and have been behind its extensibility and growth into multiple communities.

5.2. Development

Flash-X is hosted on GitHub, and Git is used for code management. The landing page: https://github.com/ECP-Astro/Flash-X has brief documentation of code and guidelines for development. Specific GitHub features used for collaboration and development are Issues and Pull Requests. Issues are used to keep track of tasks, enhancements, and bugs. There are tags to assign the nature of the issue and also various developers can discuss and plan out action or resolution items for individual issues. Pull request is a method of submitting contribution to the main code. Robust, stable, and quality code requires that each pull request be peer-reviewed and follow development guidelines to be able to be merged to the main code. Flash-X uses three major branches namely “master,” “staged,” and “development” with the following rules established:

All feature branches are based of master.

When development on a feature branch is finished, the feature branch is first merged into the development branch. The merge can be done directly by the author or by a pull request assigned to a different developer.

When the author has successfully executed on the development branch all manual tests related to changes made in the feature branch and believes that the feature branch is ready for inclusion in Flash-X, a pull request is issued for merging the feature branch into the staged branch.

One can trigger a code review at this point if so desired.

If all changes made in the feature branch are made to a portion of the code that is not yet under automated testing and if it is believed that the changes cannot break any automated tests, then the pull request can be managed fully by the creator and new tests may be added to the test suite.

If the pull request is accepted, then the merge can be made and the developer that executes the merge must launch the Jenkins staged test suite (Flash-Staged).

If all tests in Flash-Staged have passed, then the developer who performed the merge into the staged branch issues a pull request from the feature branch to the master branch. This same individual accepts the request and subsequently launches the Jenkins master test suite (Flash-Master).

The developer who triggered the Flash-Master test run follows up to confirm that the master branch was not broken by the feature branch.

If at any point during the previous steps errors are discovered that need to be fixed, the best way to go about it is this: Fix code in the feature branch; test there; then go back to point 2 above. Repeat until no test errors.

The development and staged branches should never be merged into any other branch.

5.3. Testing, verification, and validation

We view testing as critical to our software development process, and use continuous integration (CI) where possible for early detection of bugs. Testing is performed on the Flash-X code’s master, staged, and development branches mentioned earlier. The test suite includes tests at various granularities. Flash-X’s classification of tests differs from the standard classification. The unit tests can either exercise a single functionality or use multiple functionalities to pinpoint error in a single functionality using a scaffolding method described in Dubey and Wan (2018). For example, the Eos and Grid modules have single function unit tests. They are needed to exercise hydrodynamics which uses an analytically known solution for a spherical moving out shock as unit test. Even though this test cannot be run without including the Grid and Eos modules, it does pinpoint any failure to hydrodynamics if the other two units tests passed. For integration tests that exercise interoperability among components, results are compared against the baseline (simulation problem developer verified results at the time of introduction of new tests). An error is reported if the results do not match the baseline results. For more details on the methodology for building Flash-X’s test suite, see Dubey et al. (2015).

Flash-X does CI using Jenkins and the Flash-X team-developed repos called FlashTestSuite and flashTest. Because the complete test suite is quite large and not very suitable for CI, a quick running subset is identified and used for accepting pull requests. The full suite runs daily on virtual machines with multiple different software and compiler configurations.

6. Summary and future work

We have developed a strategy for incorporating fully functional applications into the ExaStar software ecosystem as modules that can be used within a larger multi-physics framework (e.g., Flash-X). A major advantage to this approach is that we can easily incorporate within the encompassing framework any future improvements to the stand-alone physics codes in physical fidelity or computational performance. For nuclear burning, we are working to offload additional kernels in the ODE integrator to the GPU in XNet. There is also ongoing work to implement in thornado $O (v / c)$ corrections to the radiation transport and port these additional kernels to GPUs. We will also continue to benefit from improvements to the AMReX framework, including, for example, those that could take advantage of the octree assumption in Flash-X to reduce overhead in grid functions. Taken in aggregate, this strategy ultimately results in a maintainable software ecosystem for exascale and beyond that enables improved scientific simulations across a variety of architectures.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the US Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan ().FLASH was developed, in part, by the DOE NNSA ASC- and DOE Office of Science ASCR-supported Flash Center for Computational Science at the University of Chicago. Funding for Flash-X and ExaStar was supported by the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the US Department of Energy under Contract No. DE-AC05-00OR22725.

ORCID iDs

J. Austin Harris

Anshu Dubey

Notes

Author biographies

J. Austin Harris is a staff scientist in the National Center for Computational Sciences at Oak Ridge National Laboratory. He received his PhD in physics from the University of Tennessee in 2015. His research interests include stellar astrophysics, reaction kinetics, and high performance computing.

Ran Chu is a PhD student in computational astrophysics at the University of Tennessee. She works as a research assistant at Oak Ridge National Laboratory since 2014. She received her master’s degree in nuclear physics in 2016. Her research interests include supernova simulation, neutrino transport, and developing next-generation supernova simulation code.

Sean M. Couch is an assistant professor at Michigan State University in the Department of Physics and Astronomy and the Department of Computational Mathematics, Science, and Engineering. He specializes in large-scale numerical simulations of astrophysical processes, especially aspects of core-collapse supernovae and massive stars. He received his PhD from the University of Texas at Austin in 2010 and was a Hubble Fellow at the Flash Center of the University of Chicago Department of Astronomy and Astrophysics from August 2011 through September 2014. Couch was then a Senior Postdoctoral Scholar at the TAPIR Walter Burke Institute for Theoretical Physics at Caltech before joining MSU in June 2015.

Anshu Dubey is a computer scientist in the mathematics and computer science division at Argonne National Laboratory and a Senior Scientist in the Department of Computer Science at the University of Chicago. She is the chief software architect for FLASH, multi-physics multiscale HPC software that is used by multiple science and engineering domains as their community code. She is interested in all aspects of HPC scientific software, with special emphasis on design, productivity, and sustainability issues.

Eirik Endeve is a staff scientist in the Multiscale Methods group at the Oak Ridge National Laboratory. He obtained his PhD in 2003 from the Institute of Theoretical Astrophysics at the University of Oslo, Norway. Eirik’s research interests include the development of numerical algorithms for hyperbolic conservation laws, including kinetic equations, and their application to problems in computational astrophysics.

Antigoni Georgiadou is an OLCF Postdoctoral Research Associate in the Advanced Computing for Nuclear, Particles, and Astrophysics Group at Oak Ridge National Laboratory. She received her PhD degree in Applied and Computational Mathematics from Florida State University in 2019. Her research interests include optimization, stellar evolution, cosmology, and high performance computing.

Rajeev Jain received his first masters in Structural Engineering from Arizona State University in 2009 and his second masters in Computer Science from University of Chicago in 2020. He has been working at Argonne National Lab. (ANL), since 2009, developing simulation software tools tailored for high performance systems. Over the past 11 years at ANL, he has been a part of development of tools for mesh generation, simulation of nuclear reactors, urban infrastructure, and astrophysics. He also develops testing infrastructure and optimization routines for exascale deep-learning-enabled precision medicine for cancer treatment.

Daniel Kasen is a faculty scientist at Lawrence Berkeley National Laboratory and a professor in the departments of physics and astronomy at UC Berkeley. He obtained his PhD from UC Berkeley. His research focuses on modeling supernovae, neutron star mergers, and other astrophysical transients.

M. Paul Laiu is a staff mathematician in the Multiscale Methods Group at Oak Ridge National Laboratory. He received his Ph.D. degree in Electrical and Computer Engineering from University of Maryland, College Park in 2016. His research interests include numerical optimization, surrogate modeling, and numerical schemes for multiscale kinetic PDEs.

O.E. Bronson Messer is Distinguished Scientist and Director of Science at the Oak Ridge Leadership Computing Facility. He received his Ph.D. in physics from the University of Tennessee in 2000. His research interests include supernovae, stellar evolution, radiative transfer, astrophysical combustion, high performance computing, and numerical relativity.

Jared O’Neal is a scientific software developer in Argonne National Laboratory’s Mathematics and Computer Science division. He currently works on computational mathematics and science applications with an eye on improving processes to ensure productivity while maintaining scientific rigor. Jared received an M.A. in physics from Boston University in 2001 and an M.M.S. in mathematics from The Ohio State University in 2017.

Michael A. Sandoval is a PhD student in computational astrophysics at the University of Tennessee. He received his bachelor’s degree in physics from San José State University in 2015. His research interests include stellar astrophysics, high performance computing, and data visualization.

Klaus Weide received a doctoral degree in physics from the University of Göttingen in Germany in 1992. He has worked as a staff member at the Flash Center at the University of Chicago since 2006, where he shared responsibility for developing and maintaining the FLASH code, adding new features, investigating and improving performance, and supporting internal and external users. As a Sr. Research Analyst, a joint appointee at Argonne National Laboratory, and a member of the Exascale Computing Project, he now focuses on portable refactoring of simulation code for new computer architectures.

References

Adams

(2001) Discontinuous finite element transport solutions in thick diffusive problems. Nuclear Science and Engineering 137(3): 298–333.

Almgren

Beckner

Bell

, et al. (2010) CASTRO: a new compressible astrophysical solver. I. Hydrodynamics and self-gravity. The Astrophysical Journal 715: 1221–1238.

Anderson

(1965) Iterative procedures for nonlinear integral equations. Journal of the ACM 12(4): 547–560. DOI: 10.1145/321296.321305.

Argast

Samland

Thielemann

F-K

, et al. (2004) Neutron star mergers versus core-collapse supernovae as dominant r-process sites in the early galaxy. Astronomy & Astrophysics 416: 997–1011.

Ascher

Ruuth

Spiteri

(1997) Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Applied Numerical Mathematics 25: 151–167.

Balbus

Hawley

(1991) A powerful local shear instability in weakly magnetized disks. I - Linear analysis. II - Nonlinear evolution. The Astrophysical Journal 376: 214–233.

Bianco

Benedicic

(2020) Gridtools. Available at: https://github.com/GridTools/gridtools (accessed 14 October 2020).

Borges

Carmona

Costa

, et al. (2008) An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws. Journal of Computational Physics 227: 3191–3211.

Brown

Byrne

Hindmarsh

(1989) Vode: a variable-coefficient ode solver. SIAM Journal on Scientific and Statistical Computing 10(5): 1038–1051. DOI: 10.1137/0910062.

10.

Bruenn

(1985) Stellar core collapse - Numerical model and infall epoch. The Astrophysical Journal Supplement Series 58: 771–841.

, et al. (2020) CHIMERA: a massively parallel code for core-collapse supernova simulations. The Astrophysical Journal Supplement Series 248(1): 11. DOI: 10.3847/1538-4365/ab7aff.

12.

Burbidge

Fowler

, et al. (1957) Synthesis of the elements in stars. Reviews of Modern Physics 29: 547–650.

13.

Byrne

Hindmarsh

(1975) A polyalgorithm for the numerical solution of ordinary differential equations. ACM Transactions on Mathematical Software 1(1): 71–96. DOI: 10.1145/355626.355636.

14.

Cameron

AGW

(1957) Nuclear reactions in stars and nucleogenesis. PASP 69: 201–222

, et al. (2019a) Realizability-preserving DG-IMEX method for the two-moment model of fermion transport. Journal of Computational Physics 389: 62–93. DOI: 10.1016/j.jcp.2019.03.037.

, et al. (2019b) Thornado-transport: IMEX schemes for two-moment neutrino transport respecting Fermi-Dirac statistics. Journal of Physics: Conference Series 1225: 012013. DOI: 10.1088/1742-6596/1225/1/012013.

17.

Colella

Glaz

(1985) Efficient solution algorithms for the Riemann problem for real gases. Journal of Computational Physics 59: 264–289.

18.

Couch

(2017) The mechanism(s) of core-collapse supernovae. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 375: 20160271.

19.

Couch

(2020) Spark: a new solver for self-gravitating compressible magnetohydrodynamics implemented in the FLASH simulation framework. in prep.

, et al. (2021) Towards performance portability in the spark astrophysical magnetohydrodynamics solver in the {Flash-X} simulation framework. in prep.

, et al. (2002) Hyperbolic divergence cleaning for the MHD equations. Journal of Computational Physics 175: 645–673.

, et al (2013a) Evolution of FLASH, a multi-physics scientific simulation code for high-performance computing. The International Journal of High Performance Computing Applications 28(2): 225–237. DOI: 10.1177/1094342013505656.

, et al. (2013b) The software development process of FLASH, a multiphysics simulation code. In: SE-CSE, San Francisco, USA, 18 May 2013.

, et al. (2016) Software process for multicomponent multiphysics codes. In: Carver

Hong

Thiruvathukal

(eds) Software Engineering for Science. London: Taylor & Francis, Chap 1.

, et al. (2009) Extensible component-based architecture for FLASH, a massively parallel, multiphysics simulation code. Parallel Computing 35(10–11): 512–522.

(2019) The dividends of investing in computational software design: a case study. The International Journal of High Performance Computing Applications 33: 322–331. DOI: 10.1177/1094342017747692.

27.

Dubey

Wan

(2018) Methodology for building granular testing in multicomponent scientific software. In: SE4Science, held in conjuction with ICSE 2018. Gothenburg, Sweden, 2 June 2018.

, et al. (2015) Ongoing verification of a multiphysics community code: FLASH. Software: Practice and Experience 45(2). DOI: 10.1177/1094342013505656.

, et al. (2014) A posteriori subcell limiting of the discontinuous Galerkin finite element method for hyperbolic conservation laws. Journal of Computational Physics 278: 47–75.

30.

Edwards

Sunderland

(2012) Kokkos array performance-portable manycore programming model. In: Proceedings of the 2012 international workshop on programming models and applications for multicores and manycores, pp. 1–10. ACM.

, et al. (2019) Thornado-hydro: towards discontinuous Galerkin methods for supernova hydrodynamics. Journal of Physics: Conference Series 1225: 012014. DOI: 10.1088/1742-6596/1225/1/012014.

, et al. (2000) Flash: an adaptive mesh hydrodynamics code for modeling astrophysical thermonuclear flashes. The Astrophysical Journal Supplement Series 131: 273–334. DOI: 10.1086/317361 .

33.

Gardiner

Stone

(2005) An unsplit Godunov method for ideal MHD via constrained transport. Journal of Computational Physics 205: 509–539.

34.

Gear

(1971) The automatic integration of ordinary differential equations. Communications of the ACM 14(3): 176–179. DOI: 10.1145/362566.362571.

(1998) Total variation diminishing Runge-Kutta schemes. Mathematics of Computation of the American Mathematical Society 67: 73–85.

36.

Hesthaven

Warburton

(2008) Nodal Discontinuous Galerkin Methods: Algorithms, Analysis and Applications. New York, NY: Springer.

37.

Hix

Meyer

(2006) Thermonuclear kinetics in astrophysics. Nuclear Physics A 777: 188–207. DOI: 10.1016/j.nuclphysa.2004.10.009.

(1999) Computational methods for nucleosynthesis and nuclear energy generation. Journal of Computational and Applied Mathematics 109: 321–351.

39.

Hornung

Keasler

(2020) Raja: managing application portability for next-generation platforms. https://computation.llnl.gov/projects/raja-managing-application-portability-next-generation-platforms.

40.

Jiang

G-S

C-c

(1999) A high-order WENO finite difference scheme for the equations of ideal magnetohydrodynamics. Journal of Computational Physics 150: 561–594.

41.

Käppeler

Gallino

Bisterzo

, et al. (2011) Thesprocess: nuclear physics, stellar models, and observations. Reviews of Modern Physics 83: 157–193.

42.

Laiu

Endeve

Chu

, et al. (2021) A DG-IMEX method for two-moment neutrino transport: nonlinear solvers for neutrino-matter coupling. Astrophysical Journal Supplement Series 253: 52.

43.

Laiu

Harris

Chu

, et al. (2020) Thornado-transport: Anderson- and GPU-accelerated nonlinear solvers for neutrino-matter coupling. Journal of Physics: Conference Series 1623: 012013. DOI: 10.1088/1742-6596/1623/1/012013.

44.

Larsen

Morel

(1989) Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes II. Journal of Computational Physics 83: 212–236.

45.

Lee

(2013) A solution accurate, efficient and stable unsplit staggered mesh scheme for three dimensional magnetohydrodynamics. Journal of Computational Physics 243: 269–292.

46.

Lee

Ramirez-Ruiz

López-Cámara

(2009) Phase transitions and he-synthesis-driven winds in neutrino cooled accretion disks: prospects for late flares in short gamma-ray bursts. The Astrophysical Journal 699: L93–L96.

47.

Lentz

Mezzacappa

Messer

OEB

, et al. (2012) On the requirements for realistic modeling of neutrino transport in simulations of core-collapse supernovae. The Astrophysical Journal 747: 73.

48.

Leveque

(2002) Finite Volume Methods for Hyperbolic Problems. Cambridge: Cambridge University Press.

49.

Liebendörfer

Mezzacappa

Thielemann

F-K

, et al. (2001) Probing the gravitational well: no supernova explosion in spherical symmetry with general relativistic Boltzmann neutrino transport. Physical Review D 63: 103004.

50.

Longland

Martin

José

(2014) Performance improvements for nuclear reaction network integration. Astronomy & Astrophysics 563: A67. DOI: 10.1051/0004-6361/201321958.

51.

Messer

OEB

Papatheodore

(2018) Exploiting hierarchical parallelism in an astrophysical equation of state using openacc and openmp. Advances in Parallel Computing 32: 517–526.

52.

Mignone

Tzeferacos

Bodo

(2010) High-order conservative finite difference GLM-MHD schemes for cell-centered MHD. Journal of Computational Physics 229: 5896–5920.

53.

Miyoshi

Kusano

(2005) A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics. Journal of Computational Physics 208: 315–344.

54.

O’Neal

Wahib

Dubey

, et al. (2020) Domaine-specific Runtime to Orchestrate Computation on Heterogeneous Platforms. Argonne National Laboratory.Technical Report ANL-20/77.

55.

Pareschi

Russo

(2005) Implicit-explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation. Journal of Scientific Computing 25: 129–155.

56.

Steiner

Hempel

Fischer

(2013) Core-collapse supernova rquations of state based on neutron star observations. The Astrophysical Journal 774: 17.

57.

Tomov

Dongarra

Baboulin

(2010) Towards dense linear algebra for hybrid GPU accelerated manycore systems. Parallel Computing 36(5–6): 232–240. DOI: 10.1016/j.parco.2009.12.005.

58.

Walker

(2011) Anderson acceleration for fixed-point iterations. SIAM Journal on Numerical Analysis 49(4): 1715–1735.

59.

Warren

Salmon

(1993) A parallel hashed oct-tree n-body algorithm. In; Proceedings of the 1993 ACM/IEEE conference on supercomputing, supercomputing ’93. NY, USA, 31 March 1993: ACM, pp. 12–21. DOI: 10.1145/169627.169640.

60.

Wasserburg

Busso

Gallino

(1996) Abundances of actinides and short-lived nonactinides in the interstellar medium: diverse supernova sources for the [ITAL]r[/ITAL]-processes. The Astrophysical Journal 466: L109–L113.

61.

Woosley

Weaver

(1995) The evolution and explosion of massive stars. II. Explosive hydrodynamics and nucleosynthesis. The Astrophysical Journal Supplement Series 101: 181.

62.

Zingale

Katz

(2015) On the piecewise parabolic method for compressible flow with stellar equations of state. The Astrophysical Journal Supplement Series 216: 31.