# HERACLES++: A multidimensional Eulerian code for exascale computing

Lou Roussel-Hard, Edouard Audit, Luc Dessart, Thomas Padioleau, Yushan Wang

► To cite this version:

Lou Roussel-Hard, Edouard Audit, Luc Dessart, Thomas Padioleau, Yushan Wang. HERACLES++: A multidimensional Eulerian code for exascale computing. Astronomy & Astrophysics - A&A, 2025, 696, pp.A201. 10.1051/0004-6361/202453566. hal-05005204

HAL Id: hal-05005204

https://hal.science/hal-05005204v1

Submitted on 26 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# HERACLES++: A multidimensional Eulerian code for exascale computing

L. Roussel-Hard<fnref n="1" />,<fnref n="2" />, E. Audit<fnref n="1" />, L. Dessart<fnref n="2" />, T. Padioleau<fnref n="1" /> and Y. Wang<fnref n="1" />

Received XXXX; accepted YYYY

## ABSTRACT

Numerical simulations of multidimensional astrophysical fluids present considerable challenges. However, the development of ex-ascale computing has significantly enhanced computational capabilities, motivating the development of new codes that can take full advantage of these resources. In this article, we introduce HERACLES++, a new hydrodynamics code with high portability, optimized for exascale machines with different architectures and running efficiently both on CPUs and GPUs. The code is Eulerian and employs a Godunov finite-volume method to solve the hydrodynamics equations, which ensures accuracy in capturing shocks and discontinuities. It includes different Riemann solvers, equations of state, and gravity solvers. It works in Cartesian and spherical coordinates, either in 1D, 2D, or 3D, and uses passive scalars to handle gases with several species. The code accepts a user-supplied heating or cooling term to treat a variety of astrophysical contexts. In addition to the usual series of benchmarking tests, we use HERACLES++ to simulate the propagation of a supernova shock in a red-supergiant star envelope, from minutes after core collapse until shock emergence. In 1D, the results from HERACLES++ are in agreement with those of V1D for the same configuration. In 3D, the Rayleigh-Taylor instability develops and modifies the 1D picture by introducing density and composition fluctuations as well as turbulence. The focus on a wedge, rather than the full solid angle, and the ability to run HERACLES++ with a large number of GPUs allow for long-term simulations of 3D supernova ejecta with a sub-degree resolution. Future developments will extend HERACLES++ to a radiation-hydrodynamics code.

**Key words.** computational fluid dynamics – shock waves – instabilities – supernova: general

## 1. Introduction

Massive stars end their lives in the gravitational collapse of their degenerate core and the formation of a compact object, possibly followed by an energetic explosion. This phenomenon is dubbed a core-collapse supernova (SN; Baade & Zwicky 1934). Most of the binding energy released by the collapsing envelope is emitted as neutrinos (Arnett et al. 1989). If revived, the SN shock born at core bounce propagates through the progenitor envelope and its highly stratified composition structure. During its expansion through the progenitor envelope, hydrodynamical instabilities grow, leading to chemical mixing on small and large scales. Spatially resolved SN remnants provide direct evidence for this chemical mixing, as well as the heterogeneous density structure of SN ejecta, with well-documented examples such as SN 1987A or Cas A (Fesen et al. 2006; Abellán et al. 2017; Milisavljevic et al. 2024). In type II SNe, which arise from the explosion of generally red-supergiant (RSG) stars, the main instability is Rayleigh-Taylor (Falk & Arnett 1973). Its development is influenced by the post-shock neutrino-driven convection that develops as part of the explosion (Kifonidis et al. 2000, 2003), as well as the structured convective envelopes of these extended progenitor stars (Goldberg et al. 2022a,b). The nonlinear growth of the Rayleigh-Taylor instability in SN ejecta has been extensively studied (e.g., see Chevalier 1976) or more recently Duffell 2016. Most of these studies were initially focused on SN 1987A (Ebisuzaki et al. 1989; Mueller et al. 1991; Fryxell et al. 1991;

Hachisu et al. 1994) and have been extended in recent years to type II SNe in general (e.g., Wongwathanarat et al. 2015).

When studying the development of instabilities in SN ejecta, the structure and chemical composition of the progenitor is critical. The Rayleigh-Taylor instability tends to develop at the He- and the O-He interfaces (Herant & Benz 1991; Mueller et al. 1991; Herant & Woosley 1993), and the resulting structures take the form of fingers and mushrooms, altering the density and chemical profiles obtained in spherical symmetry (e.g., Wongwathanarat et al. 2015). The modeling of such fluids in 3D is computationally intensive and it is therefore desirable to develop efficient and robust numerical tools to tackle these scientific problems.

More and more computational capacities are available with the creation of new supercomputers and the development of new hardware. Exascale supercomputers operate with a massively parallel architecture and have a computing power exceeding $10^{18}$ floating-point operations per second (FLOPS). The exascale era offers the ability to solve large and complex problems, opening up new possibilities for research by combining unprecedented computing power and data analysis. Indeed, the computational resources rely less and less on CPUs and benefit more and more from GPUs. Practically, among the world's fastest supercomputers, more than 75% of the available FLOPS are delivered by GPUs<fnref n="1" />. GPUs have become a crucial driver of computational performance in high-performance computing (HPC) and differ-

<footnote n="1">¹ Université Paris-Saclay, UVSQ, CNRS, CEA, Maison de la Simulation, 91191, Gif-sur-Yvette, France
e-mail: lou.roussel-hard@cea.fr</footnote>

<footnote n="2">² Institut d'Astrophysique de Paris, CNRS-Sorbonne université, 98 bis boulevard Arago, 75014 Paris, France</footnote>

<footnote n="1">¹ Computed from www.top500.org</footnote>
ent vendors (e.g., AMD, NVIDIA, Intel) are developing their own types of hardware, which can be used with their own programming languages (e.g., HIP/ROCM, CUDA, OneApi). To fully exploit these possibilities, scientists must develop software and numerical approaches that use the benefits of exascale hardware and ensure portability across different architectures.

To meet these challenges, and especially performance portability, we have developed the hydrodynamics code HERACLES++ using the Kokkos library (Trott et al. 2022). The unified C++ code is portable to different architectures, offering a level of abstraction that enhances user-friendliness. The code uses the Message Passing Interface (MPI) to deal with distributed memory parallelism and Kokkos for shared memory computations. With this approach, the code can efficiently use most of the CPUs and GPUs available on the market, including NVIDIAs and AMDs.

The usage of Kokkos is slowly spreading and different codes have been recently developed using the Kokkos-C++ framework. For example, IDEFIX (Lesur et al. 2023) and ARK (Padioleau et al. 2019; Bourgeois et al. 2024) are fixed-grid magnetohydrodynamics (MHD) codes designed to study protoplanetary discs or convection, respectively. Other codes, such as Athena++ (Stone et al. 2020) and the new revision, K-ATHENA++ (Grete et al. 2020), or DYABLO (Delorme et al. 2022), implement an adaptive mesh refinement (AMR) framework for problems ranging from solar physics to cosmology.

Some codes are more specifically dedicated to the study of core-collapse SNe (for a review, see Müller 2020). These codes were developed initially in 2D and in the last decade extended to 3D with different strategies and assets. Their aim is to explore the neutrino-driven explosion mechanism of massive stars with a variety of physical modules (e.g., neutrino transport schemes, neutrino opacities, nuclear networks, general relativity) and spatial grids (e.g. spherical polar, dendritic, Yin-Yang, Cartesian with AMR). These codes include PROMETHEUS and variants (Fryxell et al. 1991; Mueller et al. 1991; Rampp & Janka 2002; Wongwathanarat et al. 2010, 2013, 2015, 2017; Bollig et al. 2021; Gabler et al. 2021), CHIMERA (Bruenn et al. 2020; Mezzacappa et al. 2020), FORNAX (Skinner et al. 2019; Burrows et al. 2019; Vartanyan et al. 2019; Burrows et al. 2020), COCONOT-FMT (Müller et al. 2019), or FLASH (Couch et al. 2015; O'Connor & Couch 2018). TESS and variants (Duffell & MacFadyen 2011, 2013; Duffell 2016; Mandal & Duffell 2023) use a variety of expanding grids (moving Voronoi or Cartesian meshes) to study the ejecta as they turn into SN remnants at late times.

In our group, we have developed and used the radiation-hydrodynamics code HERACLES (González et al. 2007) for the study of SN ejecta. The focus has not been the explosion mechanism but instead the evolution of the SN ejecta on a longer timescale. We have studied in 3D the impact of ejecta structure (i.e., clumping) on type II SN light curves (Dessart & Audit 2019). We have also studied the interaction of SN ejecta with circumstellar material, both in spherical symmetry (Dessart et al. 2015) and in 2D (Vlasis et al. 2016). One asset of our approach is the post-processing of the radiation-hydrodynamics simulations with detailed nonlocal thermodynamic equilibrium (NLTE) radiative transfer for the production of not just light curves but also spectra with CMFGEN (Hillier & Miller 1998; Dessart & Hillier 2005; Hillier & Dessart 2012; Dessart et al. 2015, 2017). Here, we present our current effort at switching from our previous FORTRAN-based CPU-based version of HERACLES to the new code HERACLES++, written in C++, making use of the Kokkos library, and adapted for exascale architectures.

In this paper, we first present the hydrodynamic equations solved for in HERACLES++ as well as the physics treated (Sect.

2). We then describe their implementation and the numerical scheme used for the resolution of the equations (Sect. 3). Some benchmarking tests are presented in Sect. 4, followed by a more astrophysics-oriented test in Sect. 5 with the modeling of a RSG explosion from 500 s after the onset of the explosion until shock breakout. We present our conclusions in Sect. 6. We also include several appendix sections to present the results from resolution tests (Appendix A), performance tests (Appendix B), as well as additional benchmarking tests (Appendix C).

## 2. Formulation of hydrodynamical equations and physics

In this initial version of HERACLES++, we solve the equations of hydrodynamics and ignore radiative transfer. We include various equations of state (EoSs), gravity, as well as the possibility of a local heating and/or cooling term. The code is designed to work in Cartesian or spherical coordinates, and in 1D, 2D, or 3D.

### 2.1. Hydrodynamics

The equations of hydrodynamics in their conservative form with gravity are given by

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0 \qquad (1a) $$

$$ \frac{\partial(\rho\mathbf{u})}{\partial t} + \nabla \cdot (\rho\mathbf{u} \otimes \mathbf{u} + \mathbb{P}) = -\rho g \qquad (1b) $$

$$ \frac{\partial E}{\partial t} + \nabla \cdot [\mathbf{u}(E + P)] = -\rho g \cdot \mathbf{u} + \dot{e}_{hc} \qquad (1c) $$

$$ \frac{\partial(\rho f_x)}{\partial t} + \nabla \cdot (\rho u f_x) = 0, \qquad (1d) $$

where $\rho$ is the density, $\mathbf{u}$ the velocity, $\mathbb{P}$ the pressure tensor, $E$ the energy per unit volume, $E$, is the sum of the kinetic energy, $e_k = \rho u^2/2$, and the internal energy, $e$. The user-defined, power source term $\dot{e}_{hc}$ is used to handle a heating or a cooling process. The last equation (Eq. 1d) can be added to track the distribution of elements in the course of the simulation. $f_x$ is the mass fraction of species $x$ and defined such that $\sum_x f_x = 1$.

### 2.2. Gravity terms

Different types of gravity are implemented in HERACLES++. The first option is uniform gravity, which is defined by

$$ g = g_0, \qquad (2) $$

where $g_0$ is the constant and uniform gravitational acceleration. It can be used in all geometries.

A second option treats the gravity from a central mass, $M_c$, in spherical geometry. The gravitational acceleration is then given by

$$ g = g(r) = -\frac{GM_c}{r^2} e_r, \qquad (3) $$

where $G$ is the gravitational constant and $r$ is the radius.

Self-gravity of the outer layers can also be accounted for in the approximation of a spherically symmetric matter distribution. In this case, we have:

$$ g = g(r, t) = -\frac{GM(r, t)}{r^2} e_r, \qquad (4) $$

where $M(r, t)$ is the total mass within the sphere of radius $r$ at time $t$.
## 2.3. Equation of state

HERACLES++ can work with a generic EoS. In the present work, two EoSs are used. The first one treats the ideal gas, for which the pressure of the gas and the sound speed are given by

$$P = P_g = (\gamma - 1)e \quad \text{and} \quad c_s = \sqrt{\frac{\gamma P_g}{\rho}}, \qquad (5)$$

where $\gamma$, the adiabatic index, can be chosen by the user.

The second EoS treats together an ideal gas in equilibrium with radiation. In this case, the total internal energy and pressure are the sum of the gas and the radiative components:

$$e = \frac{\rho k_b T}{\mu(\gamma - 1)} + a_r T^4 \qquad (6a)$$

$$P = P_g + P_r \quad \text{with} \quad P_r = \frac{1}{3} a_r T^4, \qquad (6b)$$

where $k_b$ is the Boltzman constant, $\mu$ the mean molecular weight, and $a_r$ is the radiative constant. The sound speed is then given by

$$c_s = \sqrt{\frac{\gamma/(\gamma - 1) + 20\alpha + 16\alpha^2 P_g}{1/(\gamma - 1) + 12\alpha \rho}}, \qquad (7)$$

where $\alpha = P_r/P_g$.

For this EoS, the pressure cannot be computed directly from the internal energy, so we first need to determine the temperature. This is done using the internal energy equation 6a and a Newton-Raphson approach, which converges in approximately ten iterations. This EoS is used for simulations of stellar explosions, which are characterized by radiation-dominated conditions. This also enables the continuation with radiation hydrodynamics of an originally hydrodynamics-only simulation (the switch should be done when somewhere on the grid the conditions are ripe for photon decoupling from the gas).

## 2.4. Energy source terms: Heating and cooling

A term, $\dot{e}_{hc}$, may be added to the energy equation (Eq. 1c) in order to study specific astrophysical conditions. It may be positive (heating) or negative (cooling). It allows one for example to treat radioactive decay or radiative losses. The user can define a custom function fit to the physical setup being studied.

When studying SN ejecta, unstable isotopes decay and release neutrinos, $\gamma$ rays, and positrons in the process. At present, we treat the $^{56}$Ni decay chain, with daughter isotopes $^{56}$Co and $^{56}$Fe (Nadyozhin 1994). The associated term in the energy equation is given by

$$\dot{e}_{hc}(t) = \lambda \left\{ \left[ Q_{Ni} \left( \frac{\tau_{Co}}{\tau_{Ni}} - 1 \right) - Q_{Co} \right] e^{-t/\tau_{Ni}} + Q_{Co} e^{-t/\tau_{Co}} \right\} \rho f_{Ni}, \qquad (8)$$

where $Q_{Ni} = 1.75$ MeV and $Q_{Co} = 3.73$ MeV are the total energy ($\gamma$-ray and positrons) emitted per decay of $^{56}$Ni and $^{56}$Co, and $\lambda$ is given by

$$\lambda = \frac{1}{56m_u (\tau_{Co} - \tau_{Ni})}. \qquad (9)$$

The lifetimes of $^{56}$Ni and $^{56}$Co are given by $\tau_{Ni} = 8.8$ d and $\tau_{Co} = 111.3$ d. At the present time, we assume that this decay power is fully absorbed locally.

For radioactive decay heating, this emissivity term depends only on the total yield of $^{56}$Ni and its spatial distribution (i.e., the total number of $^{56}$Ni nuclei at all radii). However, $\dot{e}_{hc}$ may also depend on the temperature, as in the cases of nuclear burning or radiative cooling, and in this case we must iterate to find the proper $\dot{e}_{hc}$ and temperature for the corresponding EoS.

## 2.5. Monitoring of internal energy in highly supersonic regions

When the internal energy becomes too small compared to the total energy, it becomes difficult to track it precisely using the total energy conservation (Eq. 1c). This can be an important issue when, for example, the temperature is needed. Instead, we introduce the following equation to correct its evolution:

$$\frac{\partial e}{\partial t} + \nabla \cdot (\boldsymbol{u}e) + P\nabla \cdot \boldsymbol{u} = 0. \qquad (10)$$

The evolution of the internal energy is then given by

$$e(t + dt) = e(t) + \alpha \delta^1 e + (1 - \alpha) \delta^2 e, \qquad (11)$$

where $\delta^1 e$ is the evolution of the internal energy given by the integration of Eq. 1c and $\delta^2 e$ is determined by the integration of Eq. 10, using the velocities from the standard hydrodynamical step.

The condition for $\alpha$ is

$$\alpha = \begin{cases} 1 & \text{if } e \ge \varepsilon e_k \\ \left(\frac{e}{\varepsilon e_k}\right)^2 & \text{if } \varepsilon e_k > e > \frac{\varepsilon e_k}{10}, \\ 0 & \text{if } e \le \frac{\varepsilon e_k}{10} \end{cases}, \qquad (12)$$

where $\varepsilon$ is a criterion fixed by the user (we typically choose $\varepsilon = 10^{-6}$).

# 3. Numerical methods and implementation

HERACLES++ is designed to offer flexibility in selecting the physics modules to be used, to accommodate various numerical schemes as well as switching easily between different geometries and dimensions. We present below the implementation of a Godunov-type finite-volume method.

## 3.1. Numerical discretization

Our notation for space-time discretization is such that the position at the center of the cell is denoted by $r_i$ and the mesh nodes are given by $r_{i\pm1/2}$. The cell size is $\delta r$. The time interval between $t^n$ and $t^{n+1}$ is $\delta t$. The notation $U_i^n$ represents the average quantity at time $t^n$ for the cell centered at $r_i$. This notation is trivially extended to all geometries and in several dimensions. When using spherical coordinates, $i$ refers to the radius coordinate, $j$ to the polar-angle coordinate, and $k$ to the azimuthal-angle coordinate.

### 3.1.1. Mesh

In HERACLES++, the computational domain is defined as a logically Cartesian mesh that can be either in Cartesian or spherical coordinates and in 1D, 2D, or 3D. The surface elements (dS) and volumes (dV) of cells are computed based on the chosen geometry. The grid spacing in each direction can be either regular or logarithmic. Users also have the option to implement a custom
grid spacing. It is also possible to initialize the simulation from an HDF5 file containing both the physical quantities and the grid spacing.

In addition, the code includes a moving grid feature that allows one to shift the grid along a specified coordinate axis during the simulation. This functionality allows one to track regions of interest defined by the user. In the context of SN ejecta, the grid can be allowed to move in the radial direction in order to follow the expanding ejecta. With this option, the grid size can be reduced or resolution can be enhanced in specific regions. Furthermore, the elimination of cells in the inner ejecta reduces the constraints on the CFL conditions.

### 3.1.2. Momentum equation discretization in spherical coordinates

Equations 1 are written in a conservative form. However, in spherical coordinates, the tensor operator (Eq. 1b) presents geometrical source terms. We use a classical discretization method that conserves angular momentum. The azimuthal component, $\phi$, of the momentum equation is then written as:

$$ \frac{\partial \rho u_{\phi}}{\partial t} + \frac{1}{r^2} \frac{\partial^2 \rho u_{\phi} u_r}{\partial r} + \frac{1}{r \sin \theta} \frac{\partial \sin \theta \rho u_{\phi} u_{\theta}}{\partial \theta} + \frac{1}{r \sin \theta} \frac{\partial (\rho u_{\phi}^2 + P_{\phi\phi})}{\partial \phi} = - \left( \frac{\rho u_{\phi} u_r}{r} + \frac{\cot \theta \rho u_{\phi} u_{\theta}}{r} \right). \quad (13) $$

The source terms are discretized as follows:

$$ \frac{\rho u_{\phi} u_r}{r} = \frac{r_{i+1/2} - r_{i-1/2}}{(r_{i+1/2} + r_{i-1/2}) dV} \left( [\rho u_{\phi} u_r]_{i+1/2} dS_{i+1/2}^r + [\rho u_{\phi} u_r]_{i-1/2} dS_{i-1/2}^r \right) \frac{\cot \theta \rho u_{\phi} u_{\theta}}{r} = \frac{\sin \theta_{j+1/2} - \sin \theta_{j-1/2}}{(\sin \theta_{j+1/2} + \sin \theta_{j-1/2}) dV} \times \left( [\rho u_{\phi} u_{\theta}]_{j+1/2} dS_{j+1/2}^{\theta} + [\rho u_{\phi} u_{\theta}]_{j-1/2} dS_{j-1/2}^{\theta} \right) \quad (14) $$

where $dS^r$ is the surface elements in the radial direction, $dS^\theta$ in the polar direction, and $dV$ the volume element. The values at the cell interfaces are computed using outputs from the Riemann solver used for the cell update.

The corresponding equation for the polar component, $\theta$, is

$$ \frac{\partial \rho u_{\theta}}{\partial t} + \frac{1}{r^2} \frac{\partial^2 \rho u_{\theta} u_r}{\partial r} + \frac{1}{r \sin \theta} \frac{\partial \sin \theta (\rho u_{\theta}^2 + P_{\theta\theta})}{\partial \theta} + \frac{1}{r \sin \theta} \frac{\partial \rho u_{\theta} u_{\phi}}{\partial \phi} = - \left( \frac{\rho u_{\theta} u_r}{r} - \cot \theta \frac{\rho u_{\phi}^2}{r} \right), \quad (15) $$

with the first source term discretized as Eq. 14 and the second:

$$ \cot \theta \frac{\rho u_{\phi}^2}{r} = \frac{([\rho u_{\phi}^2]_{k+1/2} + [\rho u_{\phi}^2]_{k-1/2}) \cos[(\theta_{j+1/2} + \theta_{j-1/2})/2]}{2 \sin[(\theta_{j+1/2} + \theta_{j-1/2})/2]} \times \frac{(dS_{j+1/2}^{\theta} - dS_{j-1/2}^{\theta})}{2 dV}. \quad (16) $$

Other source terms are discretized in a similar manner.

### 3.2. Algorithmic steps

HERACLES++ was designed to easily incorporate a variety of numerical schemes. In this paper, we concentrate on a Godunov-type approach (Godunov & Bohachevsky 1959) and more specifically on the Monotone Upstream-centered Schemes for Conservation Laws (MUSCL)-Hancock method (van Leer 1977a,b, 1979, 1984; Toro 2009) to solve the hydrodynamical equations. This method is a simple but practical choice for achieving second-order accuracy in numerical simulations. Its ability to accurately represent gradients and discontinuities, its robustness, and its conservation properties make it a reliable tool to treat complex flows involving strong shocks.

The first step of the MUSCL scheme is to reconstruct time and space extrapolated variables at the cell interfaces. The space reconstruction is done with a classical slope-limiter approach:

$$ \tilde{U}_{L}^{n} = \tilde{U}_{i}^{n} - \frac{\delta r_{i}}{2} \sigma \quad \text{and} \quad \tilde{U}_{R}^{n} = \tilde{U}_{i}^{n} + \frac{\delta r_{i}}{2} \sigma, \quad (17) $$

where $\tilde{U}_L$ and $\tilde{U}_R$ are the reconstructed primitive variables $(\rho, u, P)$ at the left and right sides of the cell, and $\sigma$ the slope defined using a slope limiter. Three different slope limiters are implemented in HERACLES++: Minmod (Sandham & Yee 1989), Van Leer, and Van Albada (Berger & Aftosmis 2005; Toro 2009). The user can choose the slope limiter at runtime and it is also possible to remain at first order. From these primitive reconstructed variables, we can obtain the reconstructed conservative variables, $U_{L,R}$.

The time reconstructed values at the half time step are obtained in the following way:

$$ U_{L,R}^{n+1/2} = U_{L,R}^{n} + \frac{\delta t}{2\delta r_{i}} [F(U_{L}^{n}) - F(U_{R}^{n})], \quad (18) $$

where $F$ is the flux. The fluxes are computed from these space-time extrapolated values using the analytical flux function. In 3D, the reconstruction is done in a similar manner (e.g., see Sect. 16.5 of Toro 2009).

The evolution to the next time step is made using the Godunov scheme with the Riemann solver:

$$ U_{i}^{n+1} = U_{i}^{n} + \frac{\delta t}{\delta r} (\mathcal{F}_{i-1/2}^{n+1/2} - \mathcal{F}_{i+1/2}^{n+1/2}), \quad (19) $$

where $\mathcal{F}_{i-1/2}^{n+1/2}$ is the flux computed with the chosen Riemann solver using the space-time extrapolated values on each side of the interface. HERACLES++ currently implements three choices of Riemann solvers: the classical HLL and HLLC solvers (Toro 2009) and an all-mach regime solver that is much more accurate for low Mach regions (Bourgeois et al. 2024). The user can choose the solver at runtime and the structure of the code makes it easy to implement a new one.

Three boundary conditions are implemented: periodic, null gradient, and reflective. A specific module also allows for the implementation of user-defined boundary conditions.

### 3.3. Implementation

HERACLES++ is designed with an emphasis on both performance and portability. This enables the use of a single source code that can be compiled across various architectures (GPUs, CPUs), while ensuring optimal performance on each. The code is modular, making it easy to implement and switch between different slope limiters and Riemann solvers for the flux calculation.
### 3.3.1. Parallelism

The code can leverage MPI for distributed memory parallelism and the Kokkos library (Trott et al. 2022) for shared memory parallelism. The global computational domain is partitioned into equal subdomains, with each subdomain assigned to a separate MPI process. Within each subdomain, the Kokkos library is used to take advantage of shared memory parallelism.

HERACLES++ has a modular structure, with each module corresponding to a specific step in the hydrodynamic scheme. A module consists of a Kokkos kernel, which is a C++ functor, and is executed using the Kokkos parallel patterns, such as Kokkos::parallel_for or Kokkos::parallel_reduce. For example, Eq. 19 is implemented as:

```
Kokkos::parallel_for(
    "Parallel_for_example",
    Kokkos::MDRangePolicy<Kokkos::Rank<3>>
    ({0, 0, 0}, {nx, ny, nz}),
    KOKKOS_LAMBDA(int i, int j, int k) {
        F_L = compute FL(i,j,k);
        F_R = compute_FR(i,j,k);
        dv = compute_dV(i,j,k);
        U_new(i,j,k) = U(i,j,k) + dt/dv * (F_L-F_R));
```

where U_new, U, and S are all of the Kokkos::View type, which is the fundamental data type used in the Kokkos library. The MDRangePolicy specifies how to iterate over a multidimensional container.

### 3.3.2. Inputs and outputs

The input-output (I/O) operations in HERACLES++ are handled by PDI,<fnref n="2" /> a library designed to decouple high-performance simulation codes from I/O management. The simulation generates a significant amount of data, which is output by PDI into HDF5 files. We can initialize a simulation by directly assigning values to the physical fields or by reading an HDF5 file containing the values for these fields. This approach allows us to restart a simulation from a specific checkpoint by loading the corresponding data file.

Another advantage of using PDI is that the I/O configuration is fully decoupled from the simulation itself. This allows us to adjust the I/O settings in a separate YAML file without the need to recompile HERACLES++.

## 4. Benchmarking tests

We have conducted a large number of tests to validate and benchmark HERACLES++. These include the standard suite of hydrodynamical tests for shocks, advection, implosion, or explosion. In Appendix C, we present results for a subset of such tests, in particular those that were performed in 2D and 3D. In this section, we present results for the tests that have a more direct relevance to the context of multidimensional simulations of stellar explosions. Hence, we present results for a 2D Kelvin-Helmholtz instability test in Sect. 4.1, for a 3D Rayleigh-Taylor instability test in Sect. 4.2, and for a 3D Sedov blast wave test in Sect. 4.3. We then move on to discuss a more astrophysics-oriented test in Sect. 5 with the simulations in 1D and 3D of the propagation of a SN shock in a RSG star envelope.

### 4.1. Two-dimensional Kelvin-Helmholtz instability test

The Kelvin-Helmholtz instability test is useful to investigate the development of turbulence arising from a velocity shear within a fluid. Following Lecoanet et al. (2016), the initial conditions in the (x, z) plane for the density, $\rho$, the x-velocity, $u_x$, and the z-velocity, $u_z$, are given by

$$ 
\begin{aligned}
\rho &= 1 + \frac{1}{2} \left[ \tanh\left(\frac{z-z_1}{a}\right) - \tanh\left(\frac{z-z_2}{a}\right) \right], \\
u_x &= \tanh\left(\frac{z-z_1}{a}\right) - \tanh\left(\frac{z-z_2}{a}\right), \\
u_z &= 0.01 \sin(\pi x) \left[ \exp\left(-\frac{(z-z_1)^2}{\sigma^2}\right) + \exp\left(-\frac{(z-z_2)^2}{\sigma^2}\right) \right],
\end{aligned} 
\quad (20) $$

with $a = 0.05$, $\sigma = 0.2$, $z_1 = 0.5$, and $z_2 = 1.5$. The initial pressure, $P_0$, is 10 and we assume an ideal gas with $\gamma = 5/3$. The boundary conditions are periodic in every direction.

Figure 1 shows the evolution of the density at three times ($t = 0, 3, 6$). The instability results in the creation of small, alternate vortices, which cascade into turbulence. In Fig. 2, we quantify the growth rate of the instability-generated structures by evaluating the evolution of the quantity $\partial \ln u_z^{\text{max}}/\partial t$, with $u_z^{\text{max}}$ the maximum velocity at each time, and for different resolutions: 512×256 (blue), 1024×512 (red), and 2048×1024 (green). The linear phase of the instability ends at $t \sim 2.6$ and the dashed black line represents the best exponential fitting for this phase, corresponding to $0.01e^{1.9t}$. The growth rate flattens after that time and remains roughly constant until the end of the simulation, largely independent of the resolution.

We can compare our results to those published by Skinner et al. (2019). The linear evolution that they obtain ends around a time of 2.5 with a value of 1.5, which are similar to our results. Their growth rate is between 1.5 and 2, whereas ours is a little lower between 1.2 and 1.7. When comparing Fig. 1 with the counterpart in their work, we observe fewer structures and a more diffusive scheme, which can explain the lower values in our growth rate. In both simulations, the resolution modestly affects the development of turbulence. Indeed, we only obtain a slight increase in the growth rate when employing a higher resolution.

### 4.2. Three-dimensional Rayleigh-Taylor instability test

The first three-dimensional test that we present is the development of the nonlinear R-T instability. The initial conditions are taken from Skinner et al. (2019). The gravitational acceleration is constant with $g_z = -2$ and the initial density of the fluid is stratified such that:

$$ \rho(z) = \rho_0 \left(1 - \frac{\gamma - 1}{\gamma} \frac{\rho_0 g_z}{P_0}\right)^{\frac{1}{\gamma - 1}}, \quad \rho_0 = \begin{cases} \rho_h = 3 & \text{if } z \ge 0 \\ \rho_l = 1 & \text{otherwise} \end{cases}, \quad (21) $$

where $\gamma = 5/3$. The initial pressure is chosen to ensure that the initial conditions are in hydrostatic equilibrium:

$$ P = P_0 \left(\frac{\rho}{\rho_0}\right)^\gamma \quad \text{with} \quad P_0 = 2\pi(\rho_h + \rho_l)gL. \quad (22) $$

At $t = 0$, the interface between the two fluids, which is originally at $z = 0$, is displaced by an amount, $h(x,y)$, given by

$$ h(x, y) = \frac{1}{H} \sum_{8 \le k_x^2 + k_y^2 \le 16} \left[ a \cos(k_x X) \cos(k_y Y) + b \cos(k_x X) \sin(k_y Y) + c \sin(k_x X) \cos(k_y Y) + d \sin(k_x X) \sin(k_y Y) \right], $$

<footnote n="2">² PDI stands for Portable Data Interface. Information can be found at https://pdi.dev/master/</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1: Color map of the density for the 2D Kelvin-Helmholtz instability test and shown at times 0 (left), 3 (middle) and 6 (right). The grid resolution corresponds to $[n_x, n_z] = [2048, 1024]$. (See Sect. 4.1 for discussion.)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2: Instability growth rate for the 2D Kelvin-Helmholtz instability test shown in Fig. 1. This growth rate is estimated from the evolution of the quantity $\partial \ln u_z^{\max}/\partial t$ and for three different grid choices, corresponding to “Low” (blue), “Med” (red), and “High” (green) resolution (i.e., $[n_x, n_z] = [512, 256]$, [1024, 512], and [2048, 1024]). The dashed line corresponds to the best-fitting exponential to the growth rate in the linear regime, which breaks down at $t \sim 2.6$.</figcaption></figure>

(23)

where $X = 2\pi x/L$, $Y = 2\pi y/L$ and $k_x, k_y$ are the wave numbers with values in the range zero to four. The coefficients, $a, b, c$, and $d$, are sampled from a uniform distribution randomly sampled from -1 to 1, excluding the boundaries. The normalized coefficient is $H = \sqrt{(1/4)(a^2 + b^2 + c^2 + d^2)/(h_{rms})}$ and $h_{rms} = 3 \times 10^{-4}L$.

We tracked the evolution of the fluid with a passive scalar field initialized to unity for positive $z$ and zero otherwise. We assumed periodic boundary conditions in the $x$ and $y$ directions and reflective conditions in the $z$ direction. The size of the domain is $x, y \in [-L/2, L/2]$ and $z \in [-L, L]$, with $L = 10$. The number of grid zones in all three directions is $[n_x, n_y, n_z] = [512, 512, 1024]$.

Figure 3 shows the evolution of the passive scalar at $t = 5$ (left) and $t = 10$ (right). The initial perturbation of the interface between the two fluids generates small buoyant bubbles that grow, occupying a greater fraction of the volume as time passes.

We can quantify the mixing between the two fluids by evaluating the horizontally averaged concentration of the passive scalar versus $z$ defined as:

$$ \langle f_x \rangle = \frac{1}{L^2} \iint f_x dx dy. \qquad (24) $$

Article number, page 6 of 19

The vertical profile of the quantity, $\langle f_x \rangle$, is shown at multiple epochs in Fig. 4. Initially, the interface between the two fluids at $z \sim 0$ is sharp. However, as time passes, the mixing layer thickens and eventually encompasses the full height of the simulated domain. This mixing occurs on large scale and corresponds to advection of material with distinct $f_x$, whose average value covers from 0 to 1 over the full height.

### 4.3. Three-dimensional off-centered Sedov blast wave test

The test of the Sedov blast wave (Sedov 1946) is useful to describe the energy-conserving phase of a SN remnant. We performed this test using a grid in spherical geometry in order to test the precision of our numerical scheme for a non-Cartesian grid. The problem consists of a stationary background medium with a uniform density, $\rho_0 = 1$, and a low pressure, $P_0 = 0.01$. For our test, the domain was a spherical wedge with $r \in [0.1, 1.3]$ and $(\theta, \phi) \in [\pi/4, 3\pi/4]$. The blast was created with a total energy peak of $E_1 = 1$ at the position $r_0 = (r_0, \theta_0, \phi_0) = (0.7, \pi/2, \pi/2)$. The boundary conditions are transmissive in every direction, although the simulation is stopped before the blast wave reaches any boundary. The spherical polar grid counts 256 zones in all three directions.

The distance between the shock and $r_0$ is given analytically by

$$ d_{sh}(t) = \xi_0 \left( \frac{E_1 t^2}{\rho_0} \right)^{1/5}, \qquad (25) $$

with $\xi_0 \approx 1.15$ for an adiabatic EoS with $\gamma = 5/3$ (Rosswog 2022). It therefore predicts that the initial energy injection develops into an expanding spherical blast that migrates away uniformly from $r_0$.

Figure 5 shows a meridional slice of the density at $t = 0.08$ together with the analytical solution given by Eq. 25. Other slices through the 3D domain yield a similar pattern. Hence, the simulation is in good agreement with the analytical prediction, apart from a small offset in the direction where the shock is aligned with the grid. In such locations, the numerical diffusion is lower, causing a notorious problem for numerical schemes (Elling 2009).
<figure><img src="image_4.png" /><figcaption>Fig. 3: Passive scalar concentration, $f_x$, for the 3D Rayleigh-Taylor instability test at a time of 5 (left) and 10 (right). The Cartesian grid is characterized by $[n_x, n_y, n_z] = [512, 512, 1024]$. (See Sect. 4.2 for discussion.)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4: Multiepoch vertical profile of the horizontally averaged passive scalar ($f_x$) obtained in the 3D Rayleigh-Taylor instability test (see also Fig. 3, and Sect. 4.2 for discussion).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5: Meridional slice of the density for the 3D off-centered Sedov blast wave test. The time is 0.08 and the cut is done along $\phi = \pi/2$. The analytical prediction given by Eq. 25 is shown as a black line. (See Sect. 4.3 for discussion.)</figcaption></figure>

# 5. Astrophysical test: The case of a red-supergiant star explosion and the development of Rayleigh-Taylor instability structures

## 5.1. Preamble

We now turn to the more astrophysics-oriented context of a spherical explosion occurring in a RSG star. This initial work is used here for the benchmarking of the code but it fits within a larger project of simulating the various instabilities taking place in the extended envelopes of RSG stars before, during, and after their explosion. These so-called type II SNe are the most

frequently encountered type of SN from massive star explosions (see, e.g., Perley et al. 2020). Furthermore, detailed NLTE radiative-transfer calculations have demonstrated that both small and large-scale structures of SN ejecta have a significant impact
<figure><img src="image_7.png" /><figcaption>Fig. 6: Initial conditions for our HERACLES++ simulations. We show the density (top left), temperature (top right), velocity (bottom left), and the mass fraction of H, He, O, Si, and $^{56}$Ni, as a function of Lagrangian mass. The time is 500 s after the explosion trigger.</figcaption></figure>

on SN light curves and spectra (Dessart et al. 2018; Dessart & Audit 2019; Ergon & Fransson 2022) so progress in that sector is much needed in order to improve the physical consistency of radiative-transfer and radiation-hydrodynamics simulations and the robustness of their inferences.

For this work, we used the initial ejecta conditions of Dessart et al. (2024), which are based on a subset of evolutionary calculations for binary massive stars performed by Ercolino et al. (2024). Binarity, which occurs with a high frequency in massive stars (Sana et al. 2012), can lead to a wide range of H-rich envelope masses, *M* (H-env), in the pre-SN star because of mass transfer through Roche-lobe overflow. However, as long as a residual H-rich envelope remains, the He-core mass as well as the surface radii are largely unaffected by this mass loss, so that these objects should explode with similar properties (explosion energy and nucleosynthesis), merely modulated by the differing *M* (H-env). Dessart et al. (2024) indeed showed that through variations in *M* (H-env) alone, the full diversity of type II SN light curves may be reproduced. Chemical mixing and clumping may, however, differ widely between these ejecta due to the differences in their progenitor-envelope properties, with strong implications for their observables. Investigating this diversity will be the subject of a forthcoming study.

We limited ourselves to the model logP3p45 (this stands for the log of the initial system orbital period of 2818 d) of Dessart et al. (2024), whose initial orbit was so wide that it entirely avoided binary mass transfer and essentially evolved as a single star. The initial, primary mass is 12.8 $M_{\odot}$ and the system's original mass ratio is 0.95. The primary dies with a mass of 10.84 $M_{\odot}$, a He-core mass of 3.98 $M_{\odot}$, a small CO core of only about 2.2 $M_{\odot}$, and a surface radius of 898 $R_{\odot}$. The model also exhibits a stratification in composition typical of massive, RSG stars at death. The envelope shows a succession of layers of distinct composition with a massive H-rich envelope of about ~6.9 $M_{\odot}$ and a degenerate core of ~1.5 $M_{\odot}$.

<figure><img src="image_8.png" /><figcaption>Fig. 7: Variation of the radial grid spacing, $dr$, as the function of the radius, $r$. The grid used in the simulation with HERACLES++ is regular until $10 R_{\odot}$ and logarithmic beyond.</figcaption></figure>

In practice, we started our simulations from the output of the explosion simulation carried out with the 1D, Lagrangian radiation-hydrodynamics code V1D (Livne 1993; Dessart et al. 2010a,b), at a time of 500 s relative to the onset of the explosion trigger. The total energy injected in the star for the explosion was designed to yield a kinetic energy at infinity of $1.1 \times 10^{51}$ erg so by this time of 500 s, the shock has already reached the outer layers of the He core$^3$. Using HERACLES++, we first evolved this structure in 1D (Sect. 5.3) and then in 3D (Sect. 5.5) in order to study the development of Rayleigh-Taylor instabilities and their impact on chemical mixing and clumping.

## 5.2. Initial conditions and numerical setup

The ejecta properties used as initial conditions for the HERACLES++ simulations are shown in Fig. 6. The shock is lo-

<footnote>$^3$For details about the preSN evolution and explosion simulation, we refer the reader to Dessart et al. (2024) and Ercolino et al. (2024).</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 8: Ejecta properties for our RSG star explosion model as computed by HERACLES++ (solid) and V1D (dashed). For representative epochs between 2.8 hr and 1.6 d, we show the density (top row), the temperature (middle row), and the velocity (bottom row) vs. the Lagrangian mass (left column) and radius (right column). (See Sect. 5.3 for discussion.)</figcaption></figure>

cated at ~4.1 $M_\odot$ with a peak velocity around 12 000 km s⁻¹. Beyond the shock, the velocity is essentially zero and $\rho$ and $T$ slowly decrease until the progenitor surface at 10.84 $M_\odot$. Figure 6 also shows the mass fraction of the main elements (i.e., H, He, O, Si, and $^{56}$Ni). The position of each layer in the star reflects previous episodes of nuclear burning or mixing, as well as the ashes of explosive nucleosynthesis in the innermost ejecta layers. The extended H-rich envelope is essentially unchanged from its original composition (i.e., corresponding here to solar metallicity), whereas the He-core layers exhibit a stratification of ashes from H, He, C, O, Ne, and Si burning. The explosion energy at 500 s is primarily stored as kinetic energy, with around $9 \times 10^{50}$ erg.

Since our interest is to investigate the impact of fluid instabilities on the distribution of elements, we introduced five passive scalars to track the material originally present in the H-rich envelope and in the He-, O-, Si-, or $^{56}$Ni-rich shells. In the following sections, when referring to the distribution of O, we really mean

the distribution of the material in the O-rich shell (in practice, this includes the O/Si, O/Ne/Mg, and the O/C shells).

At the start of the HERACLES++ simulation, we remapped the V1D input (which is on a Lagrangian grid) onto the HERACLES++ grid, which we chose to be a hybrid linear-log grid. It has a relatively high resolution out to the shock location (in order to preserve the initial sharpness of the shock and the total initial energy of the explosion) and a decreasing resolution beyond. In Sect. 5.3, we constructed a regular radial grid with 10000 points to produce a grid that presents no major resolution offset with the V1D grid, which is Lagrangian and thus in mass space (i.e., the goal in this test is to confront the numerical solution of the hydrodynamical equations between the two codes, rather than assess the influence of resolution or numerical diffusion etc). In Sect. 5.5, we created a radial grid with a total of 1024 points, with 200 points for the regular part out to $10 R_\odot$ (corresponding to a Lagrangian mass of ~4.2 $M_\odot$). Figure 7 illustrates the variation in the radial increment, $dr$, with radius. In Appendix A, we
<figure><img src="image_10.png" /><figcaption>Fig. 9: Evolution of the total energy, $E_{tot}$, the kinetic energy, $E_{kin}$, the radiative energy, $E_{rad}$, and the internal (gas) energy, $E_{gas}$, during our HERACLES++ simulation from 500 s until shock breakout. At shock breakout, $E_{kin} \sim E_{rad}$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10: Estimated growth rate of the R-T instability in our RSG explosion simulation, shown at multiple epochs and vs. Lagrangian mass. (See Sect. 5.4 for discussion.)</figcaption></figure>

consider additional radial resolutions (i.e., $n_r$ of 256, 512, and 768) as part of a resolution study.

We performed hydrodynamical simulations with HERACLES++, using spherical coordinates, first in 1D (Sect. 5.3) and then in 3D (Sect. 5.5). We used the perfect gas+radiation EoS (Eq. 6b; we used $\gamma = 5/3$ and $\mu = 0.62$ for the ideal-gas part) and introduced $^{56}$Ni-decay heating (Eq. 8; this was done for consistency with V1D but it plays little role here since we stopped the simulation at about 1.5 d, which is only a quarter of the $^{56}$Ni lifetime). A point mass gravity (Eq. 3) was used to account for the presence of the compact remnant. We used a routine to track precisely the internal energy (Sect. 2.5). The boundary conditions were reflexive in the inner radius and transmissive at the outer radius. The slope limiter used was Minmod and the Riemann solver was HLLC. For consistency, this general setup was the same for all similar simulations.

## 5.3. Results in spherical symmetry and comparison with V1D

In this section, we first present a comparison of the results between V1D and HERACLES++, and thus in spherical symmetry.

Figure 8 shows the evolution of $\rho$, $T$, and $u$ over time as the shock propagates through the progenitor envelope and until its emergence at $R_\star$ at $\sim 1.6$ d.<fnref n="4" /> We show the variations in these quantities as a function of both the Lagrangian mass and the radius. There is a good agreement between the results of the two codes.

A pair of reverse and forward shocks forms when the SN shock reaches the outer edge of the He core very soon after the start of the simulation. The reverse shock then propagates slowly inward (in mass space) and will eventually reach the innermost ejecta layers after about 10 d (as is indicated in the V1D simulation; not shown here). The reverse shock plows through the extended envelope more swiftly but it undergoes a strong deceleration because of the relatively flat density profile of the massive extended H-rich envelope. Over the course of the simulation, the ejecta density drops from $10^{-4}$ to $10^{-7}$ g cm$^{-3}$, primarily through the expansion of the dense inner regions over a timescale of days. In contrast, the outer regions, located at large distances from the explosion site, are shocked last and undergo limited expansion. A high-density plateau forms around the peak density so that most of the shocked-envelope mass is located between 500 and 900 $R_\odot$ at shock emergence. A small and decreasing mass is present below the reverse shock owing to the large progenitor envelope mass relative to the He-core mass.

The velocity profile mirrors that of the density. From the time of the explosion, the peak velocity decreases from 7500 to 3000 km s$^{-1}$. Over time, the shock front remains steep while propagating through the envelope. The main difference between the two codes is the prediction of a temperature drop interior to the reverse shock. This feature is also present in MESA but absent in V1D because of the use of a temperature (or internal energy) floor in the EoS (for a similar pattern, see Fig. 28 of Paxton et al. 2015). There is a good agreement between the two codes for the temperature in the shocked regions, which are radiation-dominated (i.e., $P_r \gg P_g$).

The deceleration of the forward shock is also reflected by the decrease in the total kinetic energy over time. Figure 9 shows the evolution of the various energy components versus time from the start of the simulation at 500 s until shock breakout (for each energy component, the summation is over the full grid). The total energy is primarily kinetic initially, but is transformed into internal energy (i.e., primarily radiation) as the shock propagates. At shock breakout, the total energy is split essentially equally between radiative energy and kinetic energy. This radiation-dominated shock is at the origin of the burst of radiation taking place at shock emergence, which is not modeled in this hydrodynamics-only simulation.

## 5.4. Stability analysis

Before discussing the 3D simulations of our explosion model, we conduct an approximate local stability analysis using the method of Benz & Thielemann (1990) (see also Mueller et al. (1991)). It estimates the growth rate of the Rayleigh-Taylor instability in the compressible case as:

$$ \sigma = \frac{c_s}{\gamma} \sqrt{\mathcal{P}^2 - \gamma \mathcal{P} \mathcal{R}}, \qquad (26) $$

<footnote n="4">⁴ For this figure and to ensure the Eulerian grid is as fine in the inner regions as the Lagrangian grid in V1D, the HERACLES++ simulation employs 10 000 radial zones.</footnote>
<figure><img src="image_12.png" /><figcaption>Fig. 11: Multiepoch meridional slices for the 3D RSG-star explosion simulation with HERACLES++. The panels show the radial velocity (a), the normalized density, $\rho/\langle\rho\rangle$ (b), the polar velocity, $u_\theta$ (c), and the He mass fraction, $X_{\text{He}}$ (d). Within each panel, each quadrant corresponds to a cut at $\phi = \pi/2$ and progresses clockwise from 1, to 4, 8, and $12 \times 10^4$ s. (See Sect. 5.5.2 for discussion.)</figcaption></figure>

with the following definitions:

$$ \mathcal{R} = \frac{\partial \ln \rho}{\partial r} \quad \text{and} \quad \mathcal{P} = \frac{\partial \ln P}{\partial r}, \qquad (27) $$

where $c_s$ is the adiabatic sound speed accounting for both the gas and the radiation (Eq. 7). Regions with the largest growth rate with be the most unstable. This estimate is only indicative since it applies to the linear regime of the instability.

Similarly, the total time-integrated growth at a given mass is given by

$$ \frac{\xi}{\xi_0} = \exp\left(\int_{t_0}^{t} \sigma dt\right), \qquad (28) $$

where $\xi$ ($\xi_0$) is the amplitude of the instability at time $t$ ($t_0$).

Figure 10 shows the integrated growth rate as a function of the Lagrangian mass at selected epochs (in the incompressible case, $\xi/\xi_0$ is qualitatively similar but peaks at smaller values). As is expected, the most unstable region is located at the interface between the highly bound He-core and the loosely bound H-rich envelope, at $\sim 4 M_\odot$, and thus where there is a density cliff.

The small peak at $\sim 1.5 M_\odot$ is an artifact of the initial ejecta density structure reflecting the regions of the original thermal bomb. This analysis indicates that the Rayleigh-Taylor instability should develop most strongly at the H-He interface in 3D simulations, as was previously found (see, for example, Mueller et al. 1991).

## 5.5. Results from the 3D simulation

### 5.5.1. Setup

The 3D simulation with HERACLES++ was initialized with the same 1D input described in Sect. 5.2. This 1D profile was remapped onto the same radial grid and copied across all angular directions ($\theta, \phi$) both in the range $[\pi/4, 3\pi/4]$. We used the same setup as in Sect. 5.2 and added reflexive boundary conditions for the angular directions. We chose a grid with $[n_r, n_\theta, n_\phi] = [1024, 512, 512]$. Simulations using other resolutions in both radius and angle are presented in Appendix A. In the preliminary tests for the 3D simulations, we employed a va-
<figure><img src="image_13.png" /><figcaption>Fig. 12: Evolution of the angle-averaged density (top), velocity dispersion (middle), and composition (bottom; the solid [dashed] line corresponds to H [He]) from the initial time until shock emergence and versus Lagrangian mass. We show the same epochs as in Fig. 11. In the top panel, the final density obtained in the 1D counterpart is shown as a dashed line. (See Sect. 5.5.2 for discussion.)</figcaption></figure>

variety of initial perturbations of the density and velocity to act as seeds for the Rayleigh-Taylor instability. We find, however, that the instability grows from numerical noise and is also aided by the small-scale high-frequency oscillations present in the V1D input (i.e., in the post-shock region at the highest velocities). So, no initial perturbation of the fluid variables was introduced for the tests presented in this work.

### 5.5.2. Results

In this section, we present a 3D simulation of the perturbed 1D explosion model presented in Sect. 5.2. To facilitate the discussion of the results, we first present the evolution of various fluid quantities within meridional slices through the 3D domain. Then, we also discuss angle-averaged quantities for the same epochs but shown versus the Lagrangian mass rather than the radius. In the following, “<>” is meant as an average over the polar and azimuthal angles at a given radius or Lagrangian mass.

Figure 11 shows the radial velocity, the logarithm of the normalized density, the polar velocity, and the He mass fraction in the $\phi = \pi/2$ plane at 1, 4, 8, and $12 \times 10^4$ s – each quadrant corresponds to a time ordered chronologically from twelve noon. These slices are representative of the overall evolution of the 3D ejecta and other choices for these slices (i.e., a different azimuth, or cuts instead in the $(r, \phi)$ plane exhibit similar properties).

As in 1D, we see that a pair of reverse and forward shocks form soon after the start of the simulation. However, in contrast to the 1D simulation, fluctuations in the velocity, the density, and the composition quickly appear throughout the region in between the shocks. These shocks are clearly visible in Fig. 11a, bracketing a turbulent flow that is located between about 450 and 900 $R_{\odot}$ at $12 \times 10^4$ s. The reverse shock remains roughly spherically symmetric but the forward shock is corrugated, showing variations in extent at the 5–10% level. The structured, unstable flow extends all the way to the forward shock. This probably arises from the strong deceleration of the forward shock as it propagates through the extended RSG star envelope (compare the first and last quadrants in Fig. 11a; see also Fig. 8).

Figure 11b indicates that the density variations are maximum at the forward shock, reaching a factor of a few the angular average at that radius, and to a lesser extent at the reverse shock. This is partially an artifact from showing the density normalized to the average and the fact that the forward shock is corrugated. Within the shocked region, the density enhancements appear as filaments surrounded by more extended, lower-density regions that are close to the average density.

The meridional slice of the polar velocity (Fig. 11c) shows that the horizontal velocity is a sizable fraction of the radial velocity, with maximum values reached at the forward shock. In the turbulent region below the forward shock, the angular velocity is subsonic and lower than about 1000 km s⁻¹ at all times ($c_s \sim 2000$ km s⁻¹ in this turbulent region) – the angle-averaged polar and azimuthal velocities are essentially zero.

The morphology of the instability-generated structures is more vividly illustrated by the He distribution, which is the dominant element in the He-rich shell (at the 95% level) and the second-most abundant in the H-rich shell (at the 28% level). In Fig. 11d, we see the development of He-rich fingers, plumes, and mushrooms into the H-rich envelope, stretching out as time passes, and reaching out to the forward shock at the last time shown. These structures are well resolved in our calculations, which use 512 polar and azimuthal zones for full coverage in each direction of ninety degrees (see Appendix A for a comparison of the structure with different radial and angular resolutions). The He mass fraction in the most radially elongated structures is, however, much smaller than the original He mass fraction in the He-rich shell. These structures have probably expanded laterally and mixed with the H-rich material, which may in part arise from numerical diffusion.

Figure 12 shows the evolution (at the same epochs as in Fig. 11, apart from the initial time) of the angle-averaged den-
<figure><img src="image_14.png" /><figcaption>Fig. 13: Evolution of the angular properties of the quantity log(ρ/⟨ρ⟩) at a fixed Lagrangian mass (4.5 M<sub>☉</sub>). In the top row and from left to right, we show color maps of radial sections through the 3D domain at 127, 323, 500, and 635 R<sub>☉</sub> and at times of 1, 4, 8, and 12 × 10<sup>4</sup> s, respectively. In the bottom row, we show the fractional shell volume per normalized density bin for each of the shells shown in the top row.</figcaption></figure>

sity, velocity dispersion,<fnref n="5" /> and composition from the initial time until shock emergence and versus Lagrangian mass. Although the density contrast is strongly reduced because of expansion, we see that in 3D this contrast is further reduced so that the average density is essentially constant throughout the shocked envelope. It no longer shows a jump at the Lagrangian mass corresponding to the outer edge of the He core (see also Fig. 8). High-density extrema are still present at the forward and reverse shocks but they no longer have the 4π lateral coherence scale enforced in 1D, which tends to produce a softer density variation at the shocks. Both shocks stretch further in and out than what is obtained in the 1D counterpart. The velocity dispersion (middle panel of Fig. 12) shows a similar behavior to the polar velocity in the meridional slice shown in Fig. 11, with values that asymptote at shock emergence to 500-1000 km s⁻¹. Perhaps the most striking impact of the instabilities taking place until shock emergence is the strong chemical mixing (bottom panel of Fig. 12). In 1D, the original composition profile would remain unchanged when viewed relative to the Lagrangian mass coordinate. In 3D, the composition is instead progressively altered with the mean mass fraction of H increasing inward and the mean mass fraction of He increasing outward. This pollution stretches all the way to the reverse shock for H (into ejecta regions that are H

deficient) and to the forward shock for He (into regions that are already He rich). At the end of the simulation at 12 × 10⁴ s, the reverse shock barely reaches the O-rich shell where we expect a secondary Rayleigh-Taylor instability to develop (Mueller et al. 1991).

To characterize more quantitatively the lateral structure in the 3D simulation, Fig. 13 shows the properties of a slice of the quantity log(ρ/⟨ρ⟩) at different radii and times, and starting at 10⁴ s. We selected a slice at 120 R☉ initially (around a Lagrangian mass of 4 M☉), which is located just below the forward shock at 137 R☉ (the choice of the radius at the different times is made such that we track the same slice of the ejecta as it advects outward). The top row shows color maps of the density variations relative to the polar and azimuthal angles, whereas the bottom row shows histograms of the fractional volume (within that slice or shell) per bin in log(ρ/⟨ρ⟩). At 10⁴ s, structures are within 50% of the average, with naturally a larger volume occupied by the underdense material (the histogram is slightly skewed relative to the mean density). At the next epoch, the density contrast between underdense and overdense material has grown, with the density maximum reaching four times the average density. The histogram is now highly skewed with about 65% of the volume occupied by material about 50% below the average density. The structures are randomly distributed in angle (and remain so at all times). As time progresses and the ejecta expand further, the his-

<footnote n="5">⁵ $\langle u_{\text{disp}} \rangle \equiv \sqrt{\langle u_r^2 \rangle - \langle u_r \rangle^2 + \langle u_\theta^2 \rangle + \langle u_\phi^2 \rangle}$.</footnote>
<figure><img src="image_15.png" /><figcaption>Fig. 14: 3D rendering of the helium mass fraction in the SN test for two different viewing angles. The time is $12 \times 10^4$ s.</figcaption></figure>

togram becomes more symmetric and the entire material is found within 50% of the mean density. There is thus only a modest residual density contrast in the shocked envelope by the time the shock emerges from the progenitor star.

Finally, we show a 3D rendering of the He mass fraction in Fig. 14 for two lines of sight looking radially inward from outside as well as from the side. This gives an alternative view of the random distribution of the advected He-rich structures into the H-rich envelope of the progenitor star as well as the filamentary nature of the most He-rich material advecting out. The low He-mass fraction of the structures located at the largest radii suggests some level of numerical diffusion.

## 6. Conclusion

In this paper, we have presented the new multidimensional hydrodynamics code HERACLES++, which is designed for computing at exascale, with a parallelization that works both for CPUs and GPUs, and equipped with a good portability across a variety of architectures. The code is written in C++ and makes use of the Kokkos library. HERACLES++ is a Eulerian code that works in Cartesian and spherical coordinates and in 1D, 2D, or 3D. Passive scalars are employed to handle the evolution of a nonuniform initial composition. The mesh may be irregular, adopting for example a combination of a linear grid and a logarithmic grid. The code currently handles constant, point-mass, and interior-mass gravity. Two EoSs are possible with either a pure ideal gas or a mixture of ideal gas and radiation. The code accepts a user-supplied heating or cooling term in order to treat astrophysical conditions such as radioactive decay in SN ejecta or radiative losses from the photosphere of a star. The code was benchmarked with a variety of tests. We show here the results for some usual 1D tests (e.g., a 1D Sod shock tube or a 2D implosion) as well as a 2D Cartesian Kelvin-Helmholtz test, a 3D Cartesian Rayleigh-Taylor test, and a 3D Spherical off-center blast wave test. Performance tests (see Appendix B) were also carried out and reveal a satisfactory scalability of the code.

We also conducted a more astrophysics-oriented test with the study of the propagation of a SN shock through the envelope of

a RSG star, using as initial conditions a model of an exploding massive star (Ercolino et al. 2024; Dessart et al. 2024). This is a useful test, which has been studied in detail in the context of SN 1987A by Fryxell et al. (1991) and Mueller et al. (1991). In 1D, HERACLES++ produces results for the evolution of the density, temperature, or velocity that are in excellent agreement with those obtained with the 1D Lagrangian (radiation) hydrodynamics code V1D (Livne 1993). Because this configuration is Rayleigh-Taylor unstable, we repeated this calculation in 3D and evolved the model until the SN shock emerged at the surface of the progenitor star. We find that the instability develops from the earliest times, eventually extending throughout the shocked region, and leading to filamentary structures or fingers, corresponding to variations in density of the order of a few at several $10^4$ s and dropping to a few tens of percent by the time of shock emergence. The instability also erases the density jump that persists in 1D counterparts at the base of the H-rich layers (i.e., the interface between the H-rich envelope and the He core in the pre-SN star). There is dispersion in velocity both in the radial and angular directions, which may survive as turbulence in the ejecta after shock emergence. The most vivid illustration of this instability can be seen in the impact it has on the composition, with advection of H-rich gas down to the H-free innermost regions and advection of He-rich gas from the He-rich shell throughout the H-rich envelope. This behavior qualitatively reproduces previous studies on similar shocked envelopes.

This benchmarking test will serve as a basis of a forthcoming study on the shock propagation in supergiant stars differing in H-rich envelope mass as well as the long-term influence of $^{56}$Ni-decay heating in such shocked envelopes. Future work will also extend HERACLES++ to a full radiation-hydrodynamics code.

The code is available upon request on GitLab.<fnref n="6" />

*Acknowledgements.* This work was supported by the Paris Region. This project was granted access to the HPC resources of CINES under the allocations 2023-cin4698 and 2024-c1615137 made by GENCI.

<footnote n="6">⁶ https://gitlab.maisondelasimulation.fr</footnote>
## References

Abellán, F. J., Indebetouw, R., Marcaide, J. M., et al. 2017, ApJ, 842, L24

Arnett, W. D., Bahcall, J. N., Kirshner, R. P., & Woosley, S. E. 1989, ARA&A, 27, 629

Baade, W. & Zwicky, F. 1934, Proceedings of the National Academy of Science, 20, 254

Benz, W. & Thielemann, F.-K. 1990, ApJ, 348, L17

Berger, M. & Aftosmis, M. 2005, 43rd AIAA Aerospace Sciences Meeting and Exhibit - Meeting Papers

Bollig, R., Yadav, N., Kresse, D., et al. 2021, The Astrophysical Journal, 915, 28

Bourgeois, R., Tremblin, P., Kokh, S., & Padioleau, T. 2024, Journal of Computational Physics, 496, 112594

Bruenn, S. W., Blondin, J. M., Hix, W. R., et al. 2020, ApJS, 248, 11

Burrows, A., Radice, D., Vartanyan, D., et al. 2019, Monthly Notices of the Royal Astronomical Society, 491, 2715-2735

Burrows, A., Radice, D., Vartanyan, D., et al. 2020, MNRAS, 491, 2715

Chevalier, R. A. 1976, ApJ, 207, 872

Couch, S. M., Chatzopoulos, E., Arnett, W. D., & Timmes, F. X. 2015, ApJ, 808, L21

Delorme, M., Durocher, A., Brun, A. S., & Strugarek, A. 2022, in SF2A-2022: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. J. Richard, A. Siebert, E. Lagadec, N. Lagarde, O. Venot, J. Malzac, J. B. Marquette, M. N'Diaye, & B. Briot, 209-213

Dessart, L. & Audit, E. 2019, A&A, 629, A17

Dessart, L., Audit, E., & Hillier, D. J. 2015, MNRAS, 449, 4304

Dessart, L., Gutiérrez, C. P., Ercolino, A., Jin, H., & Langer, N. 2024, A&A, 685, A169

Dessart, L. & Hillier, D. J. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 332, The Fate of the Most Massive Stars, ed. R. Humphreys & K. Stanek, 427

Dessart, L., Hillier, D. J., & Audit, E. 2017, A&A, 605, A83

Dessart, L., Hillier, D. J., & Wilk, K. D. 2018, A&A, 619, A30

Dessart, L., Livne, E., & Waldman, R. 2010a, MNRAS, 408, 827

Dessart, L., Livne, E., & Waldman, R. 2010b, MNRAS, 405, 2113

Duffell, P. C. 2016, ApJ, 821, 76

Duffell, P. C. & MacFadyen, A. I. 2011, ApJS, 197, 15

Duffell, P. C. & MacFadyen, A. I. 2013, ApJ, 775, 87

Ebisuzaki, T., Shigeyama, T., & Nomoto, K. 1989, ApJ, 344, L65

Elling, V. 2009, Acta Mathematica Scientia, 29, 1647

Ercolino, A., Jin, H., Langer, N., & Dessart, L. 2024, A&A, 685, A58

Ergon, M. & Fransson, C. 2022, A&A, 666, A104

Falk, S. W. & Arnett, W. D. 1973, ApJ, 180, L65

Fesen, R. A., Hammell, M. C., Morse, J., et al. 2006, ApJ, 645, 283

Fryxell, B., Mueller, E., & Arnett, D. 1991, ApJ, 367, 619

Gabler, M., Wongwathanarat, A., & Janka, H.-T. 2021, Monthly Notices of the Royal Astronomical Society, 502, 3264-3293

Godunov, S. K. & Bohachevsky, I. 1959, Matematičeskij sbornik, 47(89), 271

Goldberg, J. A., Jiang, Y.-F., & Bildsten, L. 2022a, ApJ, 929, 156

Goldberg, J. A., Jiang, Y.-F., & Bildsten, L. 2022b, ApJ, 933, 164

González, M., Audit, E., & Huynh, P. 2007, A&A, 464, 429

Grete, P., Glines, F. W., & O'Shea, B. W. 2020, K-Athena: a performance portable structured grid finite volume magnetohydrodynamics code

Hachisu, I., Matsuda, T., Nomoto, K., & Shigeyama, T. 1994, A&AS, 104, 341

Herant, M. & Benz, W. 1991, ApJ, 370, L81

Herant, M. & Woosley, S. E. 1993, in American Astronomical Society Meeting Abstracts, Vol. 182, American Astronomical Society Meeting Abstracts #182, 19.04

Hillier, D. J. & Dessart, L. 2012, MNRAS, 424, 252

Hillier, D. J. & Miller, D. L. 1998, ApJ, 496, 407

Kifonidis, K., Plewa, T., Janka, H. T., & Müller, E. 2000, in Nuclear Astrophysics, 2000, ed. W. Hillebrandt & E. Müller, 60

Kifonidis, K., Plewa, T., Janka, H. T., & Müller, E. 2003, A&A, 408, 621

Lecoanet, D., McCourt, M., Quataert, E., et al. 2016, MNRAS, 455, 4274

Lesur, G. R. J., Baghdadi, S., Wafflard-Fernandez, G., et al. 2023, A&A, 677, A9

Liska, R. & Wendroff, B. 2003, SIAM J. Sci. Comput., 25, 995

Livne, E. 1993, ApJ, 412, 634

Mandal, S. & Duffell, P. 2023, SPROUT: A moving mesh hydro code using a uniformly expanding Cartesian grid

Mezzacappa, A., Endeve, E., Messer, O. E. B., & Bruenn, S. W. 2020, Living Reviews in Computational Astrophysics, 6, 4

Milisavljevic, D., Temim, T., De Looze, I., et al. 2024, ApJ, 965, L27

Mueller, E., Fryxell, B., & Arnett, D. 1991, A&A, 251, 505

Müller, B. 2020, Living Reviews in Computational Astrophysics, 6

Müller, B., Tauris, T. M., Heger, A., et al. 2019, Monthly Notices of the Royal Astronomical Society, 484, 3307

Nadyozhin, D. 1994, The Astrophysical Journal Supplement Series, 92, 527

O'Connor, E. P. & Couch, S. M. 2018, The Astrophysical Journal, 854, 63

Padioleau, T., Tremblin, P., Audit, E., Kestener, P., & Kokh, S. 2019, ApJ, 875, 128

Paxton, B., Marchant, P., Schwab, J., et al. 2015, The Astrophysical Journal Supplement Series, 220, 15

Perley, D. A., Fremling, C., Sollerman, J., et al. 2020, ApJ, 904, 35

Rampp, M. & Janka, H. T. 2002, A&A, 396, 361

Rosswog, S. 2022, Modelling astrophysical fluids with particles

Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444-446

Sandham, N. & Yee, H. 1989

Sedov, L. I. 1946, Journal of Applied Mathematics and Mechanics, 10, 241

Skinner, M. A., Dolence, J. C., Burrows, A., Radice, D., & Vartanyan, D. 2019, ApJS, 241, 7

Stone, J. M., Tomida, K., White, C. J., & Felker, K. G. 2020, ApJS, 249, 4

Toro, E. 2009, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction (Springer Berlin Heidelberg)

Trott, C. R. et al. 2022, IEEE Transactions on Parallel and Distributed Systems, 33, 805, conference Name: IEEE Transactions on Parallel and Distributed Systems

van Leer, B. 1977a, Journal of Computational Physics, 23, 263

van Leer, B. 1977b, Journal of Computational Physics, 23, 276

van Leer, B. 1979, Journal of Computational Physics, 32, 101

van Leer, B. 1984, SIAM Journal on Scientific and Statistical Computing, 5, 1

Vartanyan, D., Burrows, A., & Radice, D. 2019, MNRAS, 489, 2227

Vlasis, A., Dessart, L., & Audit, E. 2016, MNRAS, 458, 1253

Wongwathanarat, A., Hammer, N. J., & Müller, E. 2010, A&A, 514, A48

Wongwathanarat, A., Janka, H. T., & Müller, E. 2013, A&A, 552, A126

Wongwathanarat, A., Janka, H.-T., Müller, E., Pllumbi, E., & Wanajo, S. 2017, ApJ, 842, 13

Wongwathanarat, A., Müller, E., & Janka, H. T. 2015, A&A, 577, A48
<figure><img src="image_16.png" /><figcaption>Fig. A.1: Impact of the radial resolution on the results of the 3D SN shock simulation presented in Sect. 5.5.2. We show color maps of the quantity $log(\rho/\langle\rho\rangle)$ in a meridional slice at $\phi = \pi/2$ and corresponding to a time of $10^5$ s. Clockwise from twelve noon, we show the results for increasing radial resolution by varying $n_r$ from 256 to 512, 768, and 1024 – the angular resolution is kept fixed with $n_\theta = n_\phi = 256$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. A.2: Angular average of the radial velocity versus radius for the simulations that differ in radial resolution and shown in Fig. A.1.</figcaption></figure>

## Appendix A: Resolution study

We now present the results from resolution tests. We use the same setup as in Sect. 5.5 but now vary the number of radial zones or the number of angular zones to assess the impact on the R-T instability-generated structures at the time near shock emergence.

Figure A.1 shows meridional slices of the quantity $log(\rho/\langle\rho\rangle)$ at $10^5$ s with four simulations in which the number of radial zones was increased from 256, to 512, 768, and 1024, while the number of angular zones is kept fixed to 256 in both polar and azimuthal directions. Several effects are visible. First, depending on the resolution the initial structure is not equally well resolved so the initial model is resampled into the HERACLES++ grid with a resolution-dependent accuracy. At low resolution, both the reserve and forward shocks propagate to smaller and larger radii (see also Fig. A.2), respectively (for the forward shock, this feature is probably aided by the poorer resolution at larger radii). The post-shock structures are also more radially elongated at

<figure><img src="image_18.png" /><figcaption>Fig. A.3: Same as Fig. A.1 but for different angular resolutions. Clockwise from twelve noon, the number of polar and azimuthal zones increases from 128 to 256, 512, and 768 (with $n_\theta = n_\phi$). The number of radial zones is kept fixed at 768.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. A.4: Same as Fig. A.2 but for the simulations differing in angular resolution and shown in Fig. A.3.</figcaption></figure>

lower resolution. However, for $n_r$ of 768 and 1024, the final structures look similar and the overall dynamics (e.g., the position of the forward and reverse shocks) is similar.

Figure A.3 is a counterpart of Fig. A.1 but now showing the results for both $n_\theta$ and $n_\phi$ increased from 128 to 256, 512, and 768, while the number of radial zones is kept fixed at 768. In contrast to the radial-resolution dependence, the models with higher angular resolution produce a pair of reverse and forward shocks that reach smaller and bigger radii, respectively (see also Fig. A.4. The Rayleigh-Taylor instability-generated structures have a smaller angular scale at higher angular resolution although they are clearly present and to some extent resolved even with the lower resolution employed here. This low resolution corresponds to 0.7 deg and is already greater than typically employed in 3D simulations of core-collapse SNe (for example, an angular resolution of 1 deg is used in Wongwathanarat et al. 2015).

This study helps us define the size of the simulation for the astrophysical test of the code (Sect. 5). We needed to balance good resolution of the physics while performing simulations of reasonable size, given the available computational time, resources, and our ability to analyze them. We chose to work with
a grid having $n_r = 1024$ and $n_\theta = n_\phi = 512$, corresponding to a total of about $2.7 \times 10^8$ cells.

## Appendix B: Performance evaluation

We conduct a weak scaling test to assess the ability of the code to improve its performance when the workload and computing resources are increased in the same proportion. The test is performed on the French supercomputer Adastra<fnref n="7" /> hosted at CINES in Montpellier, France. This cluster is equipped with AMD MI250X GPUs, with four GPUs per node. Each GPU has a 128 GB memory and a theoretical peak memory bandwidth of 3276 GB/s. The MI250X GPU is composed of two Graphic Complex Die (GCD) and we assign one MPI process to each GCD. We also use the 5.5.1 version of rocm with the HIPCC compiler driver. The Kokkos version is 4.4.01 and we enable the Kokkos option Kokkos_ARCHIVE90A. The performance is evaluated for the temporal loop without the initialization and I/O. It shows the number of cells updated per second given by $perf = (n_x n_y n_z n_{itr})/t_{elaps}$, with $n_{itr}$ the number of iterations and $t_{elaps}$ the elapsed time for the calculation.

For this analysis, we used the 3D Rayleigh-Taylor test in Cartesian coordinates discussed in Sect. 4.2. Initially, we do not employ any passive scalar. We maximize the size of the problem with one GCD, which is $5 \times 10^7$ cells. After determining the maximum size on a single GPU, we repeated the calculation, each time doubling both the size of the problem and the number of GPUs. The results are shown in Fig. B.1. The ideal weak scaling is calculated based on the performance of a single MPI process and is represented as a line from this assumption. Our results slightly differ from this line.

We then compared with a test based on the configuration for the propagation of a SN shock in a stellar envelope and discussed in Sect. 5. The performance is now significantly lower with $1.3 \times 10^{10}$ cells s⁻¹ compared to $2.3 \times 10^{10}$ cells s⁻¹ in the Cartesian case and 512 MPI processes. This decrease is due to the additional geometrical terms appearing in the equations when spherical coordinates are used as well as the use of an EoS with a mixture of ideal gas and radiation. The use of passive scalars also comes with a sizable computational cost, further reducing the performance in HERACLES++.

We have also tested HERACLES++ on a variety of other architectures including Nvidia GPU's from the IDRIS computing center and various kinds of CPUs (from standard notebook to supercomputers). Thanks to the usage of Kokkos, HERACLES++ can run efficiently on all these systems just by changing the compiling options.

## Appendix C: Additional hydrodynamic tests

### Appendix C.1: Sod shock tube test

The classical shock tube test is designed to demonstrate the code's shock-capturing capabilities. We adopted 1D Cartesian coordinates, an ideal gas with $\gamma = 1.4$, and transmissive boundary conditions. The spatial domain is $x \in [0, 1]$, sampled with 100 grid points. The initial conditions were taken from Toro (2009) and such that

$$ (\rho, u, P)(x) = \begin{cases} (1.0, 0, 1.0) & \text{if } x < 0.5 \\ (0.125, 0, 0.1) & \text{if } x > 0.5. \end{cases} \qquad (C.1) $$

<figure><img src="image_20.png" /><figcaption>Fig. B.1: Weak scaling test showing the performance of the code when the size of the problem and the number of MPI processes are increased in the same proportion. We show this performance for different configurations (Cartesian vs spherical coordinates; with and without passive scalars).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. C.1: Results for the 1D Sod shock tube test. We show the profiles at $t = 0.2$ for the density (top), the velocity (middle), and the pressure (bottom) obtained with HERACLES++ (solid) and predicted analytically (dashed).</figcaption></figure>

This test is not particularly challenging for modern hydrodynamics codes but it is interesting because it exhibits all fundamental hydrodynamical waves. It is also well documented. The exact solution is known analytically and consists of a propagating rarefaction wave on the left side, a rightward-propagating contact discontinuity, and a rightward-propagating shock.

<footnote n="7">⁷ Adastra Top500: https://www.top500.org/system/180051/</footnote>
<figure><img src="image_22.png" /><figcaption>Fig. C.2: Evolution of the passive scalar $f_x$ for the 2D Rayleigh-Taylor instability test. From left to right, we show the results at $t = 0, 1.5$, and 3. (see Appendix C.2 for details.)</figcaption></figure>

Figure C.1 shows the analytic solution and the results from HERACLES++. The code correctly captures the different waves. The contact wave is the most diffuse, captured within approximately four grid points, whereas the shock wave spans around two grid points. The rarefaction wave is well described by our numerical scheme.

## Appendix C.2: 2D Rayleigh-Taylor test

In Sect. 4.2 we presented results for the Rayleigh-Taylor instability test in 3D. Here, we consider a similar problem but in 2D and with different initial perturbations, which we took from Skinner et al. (2019).

The fluid is contained in a box of size defined by $x \in [-0.5, 0.5]$, $z \in [-1, 1]$, and with $n_x = 1024$ and $n_z = 2048$. The gravitational acceleration is constant with $g_z = -0.5$ and the initial density of the stationary fluid is given by

$$ \rho(x, z) = \begin{cases} 2 & \text{if } z \ge 0, \\ 1 & \text{otherwise.} \end{cases} \qquad (C.2) $$

The initial pressure is given by the conditions for hydrostatic equilibrium:

$$ P(z) = P_0 + \int_{-1}^{z} \rho g_z dz, \qquad (C.3) $$

with $P_0 = 10/7 + 1/4$. For this test we adopted an ideal gas with $\gamma = 1.4$. At time $t = 0$, the interface between the two fluids was perturbed according to

$$ h(x) = h_0 \cos(\kappa x), \qquad (C.4) $$

where $h_0 = 0.01$ and $\kappa = 4\pi$. The sharp interface is then smoothed into an hyperbolic tangent profile with characteristic length 0.005. The evolution of the fluid is tracked using a passive scalar $f_x$, which is set to zero where $z < 0$ and 1 elsewhere.

Figure C.2 shows the increasing mixing between the two fluids as time progresses, as well as the formation of the characteristic mushroom structures typical of the Rayleigh-Taylor instability. With the heavier fluid above, the slightest disturbance of the interface amplifies exponentially. Figure C.3 compares the evolution of the size of the mixing layer according to the analytical theory: $h(t) = h_0 \cosh(\sqrt{Ag\kappa t})$, with $A = 1/3$, the Atwood number. The result shows a great agreement with the linear theory up to $t \sim 1.5$, when the Rayleigh-Taylor structures start to develop.

<figure><img src="image_23.png" /><figcaption>Fig. C.3: Amplitude of the interface perturbation for the single mode 2D Rayleigh-Taylor test and comparison with the analytical solution.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. C.4: Results for the 2D implosion test. We show a color map of the density at $t = 2.5$, with the $x = z$ symmetry axis in black dashed line. (see Appendix C.3 for details.)</figcaption></figure>

## Appendix C.3: Liska-Wendroff implosion test

The 2D implosion test of Liska & Wendroff (2003) consists of Sod-like initial conditions, but with a 45 deg rotation. The
boundary conditions are reflective in every direction and the size of the domain is $(x, z) \in [0, 0.3]$, which we cover with 1024 zones in each direction. We again adopted an ideal gas with $\gamma = 1.4$. Following Liska & Wendroff (2003), the initial conditions of the stationary fluid are

$$ (\rho, P)(x, z) = \begin{cases} (1, 1) & \text{if } x+z \ge 0.015, \\ (0.14, 0.125) & \text{otherwise.} \end{cases} \qquad (C.5) $$

This setup results in a shock moving away from the initial discontinuity, eventually reflecting off the boundaries and interacting with itself. The solution includes a low-density jet that travels along the symmetry axis $x = z$. This problem is very sensitive to the preservation of the reflection symmetry with the jet failing to form if symmetry is violated.

We show the results for this test in Fig. C.4 and confirm the propagation of the jet along the symmetry axis. The distance traveled by the jet is a useful probe of the numerical diffusion. Indeed, at much higher resolutions, the jet is found to propagate further (not shown).

## Appendix C.4: Multidimensional Sedov blast wave

In Sect. 4.3 we present a 3D version of the off-center Sedov blast wave in the spherical coordinate. We also test the Sedov explosion in spherical 1D and 2D, with a centered explosion. We compare with the analytical solution.

For the 1D simulation, the domain is $r \in [0, 1]$ (which we sample with 400 grid points) and the explosion is initiated at $r = 0$. The initial conditions are $\rho = 1$, $u = 0$ and $\gamma = 1.4$. $E = E_0 = 10^{-12}$ everywhere except at the point of explosion where $E = E_1 = 1$ at the point of the explosion. The boundary conditions are reflexive at the center and transmissive at the maximum radius. The 2D counterpart is started from the results of the 1D simulation at $t = 0.3$, with the fluid properties copied across all polar directions. The spherical grid is such that $r \in [0.3, 1.3]$, $\theta \in [\pi/4, 3\pi/4]$, with 400 grid points in each direction.

Figure C.5 compares the density and total energy at $t = 0.7$ for the 1D and 2D simulations (we also added the analytical solution for each quantity). The position of the shock predicted by HERACLES++ is similar in 1D and 2D and it also agrees with the analytical solution. Figure C.6 shows the density at $t = 0.7$. The shock remains spherical and its outer radius follows the analytical prediction (Eq. 25, with $\xi_0 \approx 1.15$).

<figure><img src="image_25.png" /><figcaption>Fig. C.5: Comparison of the density (top) and total energy (bottom) at $t = 0.7$ for the Sedov blast wave test performed with HERACLES++ in 1D (blue) and 2D (pink; we show the result along $\theta = \pi/2$). We also show the analytic solution (black dashed line).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. C.6: Color map of the density at $t = 0.7$ for the 2D Sedov blast wave test (the analytic solution is shown as a black dashed line).</figcaption></figure>