SPH-EXA: Exascale SPH Mini-App

• SPH-EXA is an exascale-oriented smoothed particle hydrodynamics mini-app that consolidates core kernels from legacy codes into a compact, performance-focused platform.
• It employs hybrid MPI+CUDA/HIP parallelism and advanced data structures to achieve scalable dynamic load balancing and fault tolerance on exascale systems.
• The software underpins high-resolution astrophysical and turbulence studies, enabling numerical convergence and reassessment of simulation claims.

Searching arXiv for recent SPH-EXA papers and foundational mini-app papers. arXiv_search(query="SPH-EXA", max_results=10, sort_by="submittedDate") SPH-EXA is an exascale-oriented smoothed particle hydrodynamics (SPH) software effort that began as an interdisciplinary co-design project to distill the core kernels and performance-critical components of representative SPH codes into a compact mini-app, and is later described as a GPU-based SPH code with hybrid MPI+CUDA or HIP parallelism, GPU-resident hydrodynamics, and scalable gravity for astrophysical applications (Cavelan et al., 2020, Cabezón et al., 13 Mar 2025). Across its mini-app and application-focused descriptions, SPH-EXA is associated with the design objectives “Bigger, Faster, Longer”: larger particle counts, significantly faster execution than existing codes, and resilience and reproducibility over long runs on exascale hardware (Guerrera et al., 2019). The project occupies a dual role in the literature: it is both a software-engineering vehicle for exascale SPH and a numerical platform for high-resolution studies of turbulence and tidal disruption events.

1. Co-design origin and stated objectives

The SPH-EXA project emerged from the observation that SPH simulations in astrophysics and CFD are computationally demanding because they combine meshless neighbor discovery, irregular data structures, multiple time stepping, and, in many applications, long-range interactions such as self-gravity (Guerrera et al., 2018). Early project papers frame SPH-EXA as a mini-app intended to consolidate the common basic SPH features of existing parent codes into a pure-SPH, Exascale-ready, optimized code path while also providing direct feedback to those parent codes (Guerrera et al., 2018).

The parent-code analysis emphasized three established implementations. In the 2020 mini-app review, these are ChaNGa, SPH-flow, and SPHYNX, whose complementary strengths include octree-based gravity and domain decomposition, high-order kernels and GPU-friendly data structures, and the Integral Approach to Derivatives (IAD) with sinc-family kernels (Cavelan et al., 2020). The co-design rationale was not merely to reimplement existing algorithms, but to isolate reusable “operands” and modules such as density evaluation, gradients, pressure forces, viscosity, time stepping, neighbor search, and optional self-gravity in a compact code base (Guerrera et al., 2019).

A common methodological feature of the early SPH-EXA literature is the use of shared benchmark problems to expose bottlenecks. The rotating square patch was implemented as a joint test for the parent codes, and an Evrard collapse case was also used to stress differences relevant to astrophysical applications (Guerrera et al., 2018). Profiling on Piz Daint for a 3D Rotating Square Patch with $10^6$ particles and 40 timesteps revealed load imbalance in MPI+OpenMP loops, divergent domain-decomposition strategies, and overhead in tree-based neighbor search and self-gravity multipole walks (Guerrera et al., 2019). These observations directly shaped the mini-app’s focus on scalability, portability, extensibility, dynamic load balancing, and fault tolerance.

2. Software architecture and execution model

In its mini-app formulation, SPH-EXA is described as a lightweight, flexible, header-only C++14 code with no external software dependencies, with parallelism selected at compile time and expressed through multiple programming models, including MPI + OpenMP, MPI + OpenMP + OpenMP 4.5 target offload, MPI + OpenMP + OpenACC, MPI + OpenMP + CUDA, and an experimental MPI + HPX configuration (Cavelan et al., 2020). A related 2019 description characterizes the mini-app as a modern C++ code base designed to encapsulate core SPH operands and classical SPH test problems under a compact implementation (Guerrera et al., 2019). By 2025, SPH-EXA is described as written in C++ with a hybrid MPI + CUDA or HIP design in which all of the core SPH kernels—density, neighbour search, IAD gradients, artificial viscosity, generalized volume elements, and force loops—have been fully ported to modern GPUs (Cabezón et al., 13 Mar 2025).

The data structures and decomposition strategies reported for SPH-EXA vary by phase and application. The mini-app literature describes an H-tree, an orthogonal recursive bisection variant, that partitions the global bounding box into one subdomain per MPI rank; each leaf cell stores a range of particle indices, and each domain builds a local octree for particle lookup, yielding a tree-walk complexity of $O(N_{\rm loc}\log N_{\rm loc})$ (Guerrera et al., 2019). The 2020 review additionally reports particles stored in DFS/Morton order of the local octree for cache locality and a global octree “top” identical on every rank, updated by two MPI_AllReduce operations per timestep, with point-to-point halo exchanges for neighbor particles (Cavelan et al., 2020). In the turbulence-oriented GPU implementation, space-filling curves of Morton/Z-order type are used to partition the domain across MPI ranks, local cell-lists are built on the GPU, and communication is overlapped with computation (Cabezón et al., 13 Mar 2025). In the tidal disruption event calculations, the gravity tree is a distributed-memory Barnes–Hut variant with opening angle $\theta_0<0.5$, built on GPUs using the “Cornerstone” oct-tree construction algorithms, and tree walks for gravity and neighbour searches are fused to maximize data locality and reuse (Kubli et al., 30 Oct 2025).

Gravity support is a defining architectural feature. The mini-app papers inherited multipole gravity from ChaNGa and discussed self-gravity as an optional module (Cavelan et al., 2020). The 2025 turbulence paper reports a highly-scalable, octree-based Fast Multipole Method running on the GPU, with tree building, multipole expansions, and local essential tree exchanges all GPU-driven (Cabezón et al., 13 Mar 2025). The TDE implementation instead emphasizes a distributed Barnes–Hut tree coupled to the hydrodynamics loop, with boundary exchanges of the locally essential tree and minimal MPI traffic for tree-boundary data (Kubli et al., 30 Oct 2025). Taken together, these descriptions indicate that SPH-EXA is not restricted to a single decomposition or gravity backend; rather, it is organized around interchangeable, performance-oriented building blocks.

3. Discretization, kernels, and physical modules

The SPH discretization used in SPH-EXA is grounded in standard kernel-weighted particle sums, but the later code descriptions emphasize a modernized formulation aimed at reducing gradient error and dissipation. In the 2020 mini-app review, the density, momentum, and energy equations are written in an IAD-based conservative form with artificial viscosity, using a sinc-family kernel and the usual grad-$h$ correction $\Omega_a$ (Cavelan et al., 2020). The turbulence paper presents a more detailed Lagrangian formulation in which the momentum equation is

$$\frac{d v_{i,a}}{dt}= -\sum_b m_b \Bigl[ \frac{X_a^{2-\sigma}X_b^{\,\sigma}\,P_a}{\Omega_a\,m_a^2\,k_a}\,A_{i,ab}(h_a) +\frac{X_b^{2-\sigma}X_a^{\,\sigma}\,P_b}{\Omega_b\,m_b^2\,k_b}\,A_{i,ab}(h_b) \Bigr] +\sum_b m_b\,\Pi_{ab}^{AV}\,\widehat{x}_{ab,i},$$

with a corresponding energy equation in which $X_a=m_a/\rho^0_a$, $k_a=\sum_b X_bW_{ab}$, and $\sigma\in[0,1]$ is a local Lagrangian-compatibility parameter (Cabezón et al., 13 Mar 2025).

A central numerical component is the Integral Approach to Derivatives. For each particle $a$, SPH-EXA solves a $O(N_{\rm loc}\log N_{\rm loc})$0 linear system,

$O(N_{\rm loc}\log N_{\rm loc})$1

where

$O(N_{\rm loc}\log N_{\rm loc})$2

and uses the matrix-free form

$O(N_{\rm loc}\log N_{\rm loc})$3

inside the momentum and energy loops (Cabezón et al., 13 Mar 2025). The same IAD structure appears in the TDE implementation, where it is explicitly credited with requiring only $O(N_{\rm loc}\log N_{\rm loc})$4 neighbours while yielding higher accuracy and reduced kernel noise (Kubli et al., 30 Oct 2025).

Kernel choice is another characteristic feature. SPH-EXA uses the sinc-family

$O(N_{\rm loc}\log N_{\rm loc})$5

and the turbulence paper states that high-$O(N_{\rm loc}\log N_{\rm loc})$6 sinc kernels suppress pairing up to very high $O(N_{\rm loc}\log N_{\rm loc})$7, requiring fewer neighbours than Wendland kernels (Cabezón et al., 13 Mar 2025). In the TDE calculations, the chosen kernel is a sixth-order sinc kernel with 100 neighbours (Kubli et al., 30 Oct 2025).

Artificial viscosity is Monaghan-style, but with modern switching and cleaning procedures. The turbulence paper gives

$O(N_{\rm loc}\log N_{\rm loc})$8

with $O(N_{\rm loc}\log N_{\rm loc})$9, augmented by a linear field cleaner that replaces the naive velocity jump by a cleaned jump $\theta_0<0.5$0, and by a dynamic individual $\theta_0<0.5$1 that relaxes between $\theta_0<0.5$2 and $\theta_0<0.5$3 through Read–Cullen-type switches (Cabezón et al., 13 Mar 2025). Generalized volume elements are defined through

$\theta_0<0.5$4

and are coupled to an Atwood-number-based switch $\theta_0<0.5$5 with $\theta_0<0.5$6 and $\theta_0<0.5$7 to suppress tensile instability while preserving conservation elsewhere (Cabezón et al., 13 Mar 2025). For the TDE problem, the hydrodynamics is supplemented by a polytropic equation of state $\theta_0<0.5$8 with $\theta_0<0.5$9, and relativistic apsidal precession is represented with a pseudo-Newtonian “Einstein” potential

$h$0

evaluated on the GPU (Kubli et al., 30 Oct 2025).

4. Scalability, load balancing, and fault tolerance

The performance literature on SPH-EXA is unusually detailed because the code originated from performance analysis as much as from numerical-method development. In the 2019 rotating square patch study on Piz Daint, the SPH-EXA mini-app achieved approximately $h$1 per timestep on 1 node and approximately $h$2 on 16 nodes, while the parent codes required approximately $h$3–$h$4 on 1 node and converged to $h$5 on 16 nodes, implying an approximately $h$6 speedup of SPH-EXA over the legacy codes at all scales (Guerrera et al., 2019). The same study reports strong-scaling efficiency $h$7 dropping from approximately $h$8 at 4 nodes to approximately $h$9 at 32 nodes, with good scaling up to 16 nodes and deterioration once $\Omega_a$0 falls below approximately $\Omega_a$1 particles per core; weak-scaling efficiency remains approximately $\Omega_a$2 up to 256 nodes and approximately $\Omega_a$3 at 512 nodes for constant $\Omega_a$4 per rank (Guerrera et al., 2019).

The 2020 mini-app review extends these results to a production supercomputer setting. Weak scaling on Piz Daint is reported for 32 million particles per node from 1 to 2,048 nodes, corresponding to a total of $\Omega_a$5 million particles, with a parallel efficiency of 67% at 2,048 nodes; strong scaling is reported for 267 million particles from 64 to 2,048 nodes, with average time per iteration decreasing nearly linearly up to 1,024 nodes and still improving at 2,048 nodes (Cavelan et al., 2020). The same paper identifies the dominant losses at scale as build-tree and domain-decomposition overheads, and notes that the average halo-to-node ratio rises from approximately 30% at 64 nodes to more than 200% at 2,048 nodes (Cavelan et al., 2020).

Dynamic load balancing is treated as a multi-level problem. A two-level approach described for the mini-app combines process-level redistribution of H-tree leaves through the DLS4LB library and thread-level weighted dynamic scheduling based on past-step timing, with the goal of reducing imbalance without increasing global communication volume because subdomains are reassigned only along tree edges (Guerrera et al., 2019). The 2025 turbulence implementation reports re-decomposition when the particle count per rank deviates by more than approximately 10%, GPU-resident neighbor search and force loops, and weak scaling above 90% in 1k–8k GPU runs on LUMI up to 50 billion particles; strong scaling remains in the 70–80% range down to 100k particles per GPU, and neighbor search, density, gradient, and force loops achieve $\Omega_a$6 per GPU in double precision (Cabezón et al., 13 Mar 2025). In the TDE work, SPH-EXA is reported to scale strongly and weakly to $\Omega_a$7 particles across thousands of GPUs on LUMI-G and ALPS, with peak per-GPU performance between $\Omega_a$8 and $\Omega_a$9 particle-timesteps per second on a GH200 GPU and a memory footprint of approximately 80–100 bytes per particle (Kubli et al., 30 Oct 2025).

Fault tolerance and reproducibility are explicit components of the exascale design. Planned mechanisms include multilevel checkpoint/restart with in-memory snapshots, node-local SSD or burst-buffer dumps, and periodic full checkpoints to the parallel file system, as well as silent data corruption detection via checksums and algorithmic checks such as global conservation of mass and momentum within tolerance (Guerrera et al., 2019). The 2020 review also lists fault tolerance, silent corruption detection, optimal checkpoint intervals, and mixed precision with error control among the principal future extensions (Cavelan et al., 2020).

5. Turbulence, filtering, and LES-consistent formulations

A major numerical theme connected to SPH-EXA is the treatment of turbulence. The theoretical basis is clarified by the filtered-SPH analysis of “SPH Consistent With Explicit LES,” which derives the locally averaged Navier–Stokes equations and shows that, with a proper choice of velocity smoothing, SPH and LES solve the same filtered equations, differing only in representation: integro-differential for SPH and differential for LES (Chola, 2018). With the filter

$$\frac{d v_{i,a}}{dt}= -\sum_b m_b \Bigl[ \frac{X_a^{2-\sigma}X_b^{\,\sigma}\,P_a}{\Omega_a\,m_a^2\,k_a}\,A_{i,ab}(h_a) +\frac{X_b^{2-\sigma}X_a^{\,\sigma}\,P_b}{\Omega_b\,m_b^2\,k_b}\,A_{i,ab}(h_b) \Bigr] +\sum_b m_b\,\Pi_{ab}^{AV}\,\widehat{x}_{ab,i},$$0

and particle motion at the filtered velocity

$$\frac{d v_{i,a}}{dt}= -\sum_b m_b \Bigl[ \frac{X_a^{2-\sigma}X_b^{\,\sigma}\,P_a}{\Omega_a\,m_a^2\,k_a}\,A_{i,ab}(h_a) +\frac{X_b^{2-\sigma}X_a^{\,\sigma}\,P_b}{\Omega_b\,m_b^2\,k_b}\,A_{i,ab}(h_b) \Bigr] +\sum_b m_b\,\Pi_{ab}^{AV}\,\widehat{x}_{ab,i},$$1

the filtered momentum equation contains the subgrid stress

$$\frac{d v_{i,a}}{dt}= -\sum_b m_b \Bigl[ \frac{X_a^{2-\sigma}X_b^{\,\sigma}\,P_a}{\Omega_a\,m_a^2\,k_a}\,A_{i,ab}(h_a) +\frac{X_b^{2-\sigma}X_a^{\,\sigma}\,P_b}{\Omega_b\,m_b^2\,k_b}\,A_{i,ab}(h_b) \Bigr] +\sum_b m_b\,\Pi_{ab}^{AV}\,\widehat{x}_{ab,i},$$2

equivalently $$\frac{d v_{i,a}}{dt}= -\sum_b m_b \Bigl[ \frac{X_a^{2-\sigma}X_b^{\,\sigma}\,P_a}{\Omega_a\,m_a^2\,k_a}\,A_{i,ab}(h_a) +\frac{X_b^{2-\sigma}X_a^{\,\sigma}\,P_b}{\Omega_b\,m_b^2\,k_b}\,A_{i,ab}(h_b) \Bigr] +\sum_b m_b\,\Pi_{ab}^{AV}\,\widehat{x}_{ab,i},$$3 (Chola, 2018). That paper further states that, although the Smagorinsky model is common in LES, approximate deconvolution methods are the most natural choice for SPH, and it emphasizes a fundamental consequence of the filtered formulation: the equations are nonconservative in nature and a correct Lagrangian cannot be easily constructed (Chola, 2018).

The explicit implications for SPH-EXA are stated in implementation terms. For a large-scale SPH-EXA implementation, one must store both filtered fields and closures, compute the pairwise forces arising from $$\frac{d v_{i,a}}{dt}= -\sum_b m_b \Bigl[ \frac{X_a^{2-\sigma}X_b^{\,\sigma}\,P_a}{\Omega_a\,m_a^2\,k_a}\,A_{i,ab}(h_a) +\frac{X_b^{2-\sigma}X_a^{\,\sigma}\,P_b}{\Omega_b\,m_b^2\,k_b}\,A_{i,ab}(h_b) \Bigr] +\sum_b m_b\,\Pi_{ab}^{AV}\,\widehat{x}_{ab,i},$$4 explicitly and symmetrize them to guarantee linear and angular momentum conservation, close the non-conservative energy budget with turbulence-dissipation source terms, and provide data structures that accommodate both particle averages and reconstruction stencils (Chola, 2018). This suggests a route by which SPH-EXA could incorporate explicit LES closures without abandoning its particle framework, but only by accepting that filtered SPH is not identical to the baseline conservative Lagrangian formulation.

The practical turbulence results reported for SPH-EXA are correspondingly significant. The 2025 study on subsonic turbulence states that SPH has traditionally been hindered by $$\frac{d v_{i,a}}{dt}= -\sum_b m_b \Bigl[ \frac{X_a^{2-\sigma}X_b^{\,\sigma}\,P_a}{\Omega_a\,m_a^2\,k_a}\,A_{i,ab}(h_a) +\frac{X_b^{2-\sigma}X_a^{\,\sigma}\,P_b}{\Omega_b\,m_b^2\,k_b}\,A_{i,ab}(h_b) \Bigr] +\sum_b m_b\,\Pi_{ab}^{AV}\,\widehat{x}_{ab,i},$$5 errors, inaccurate gradient evaluations, and excessive dissipation, and reports that SPH-EXA reproduces the inertial-range behavior of AREPO and GIZMO in the subsonic regime (Cabezón et al., 13 Mar 2025). At $$\frac{d v_{i,a}}{dt}= -\sum_b m_b \Bigl[ \frac{X_a^{2-\sigma}X_b^{\,\sigma}\,P_a}{\Omega_a\,m_a^2\,k_a}\,A_{i,ab}(h_a) +\frac{X_b^{2-\sigma}X_a^{\,\sigma}\,P_b}{\Omega_b\,m_b^2\,k_b}\,A_{i,ab}(h_b) \Bigr] +\sum_b m_b\,\Pi_{ab}^{AV}\,\widehat{x}_{ab,i},$$6 particles it recovers a clear $$\frac{d v_{i,a}}{dt}= -\sum_b m_b \Bigl[ \frac{X_a^{2-\sigma}X_b^{\,\sigma}\,P_a}{\Omega_a\,m_a^2\,k_a}\,A_{i,ab}(h_a) +\frac{X_b^{2-\sigma}X_a^{\,\sigma}\,P_b}{\Omega_b\,m_b^2\,k_b}\,A_{i,ab}(h_b) \Bigr] +\sum_b m_b\,\Pi_{ab}^{AV}\,\widehat{x}_{ab,i},$$7 slope from $$\frac{d v_{i,a}}{dt}= -\sum_b m_b \Bigl[ \frac{X_a^{2-\sigma}X_b^{\,\sigma}\,P_a}{\Omega_a\,m_a^2\,k_a}\,A_{i,ab}(h_a) +\frac{X_b^{2-\sigma}X_a^{\,\sigma}\,P_b}{\Omega_b\,m_b^2\,k_b}\,A_{i,ab}(h_b) \Bigr] +\sum_b m_b\,\Pi_{ab}^{AV}\,\widehat{x}_{ab,i},$$8 to $$\frac{d v_{i,a}}{dt}= -\sum_b m_b \Bigl[ \frac{X_a^{2-\sigma}X_b^{\,\sigma}\,P_a}{\Omega_a\,m_a^2\,k_a}\,A_{i,ab}(h_a) +\frac{X_b^{2-\sigma}X_a^{\,\sigma}\,P_b}{\Omega_b\,m_b^2\,k_b}\,A_{i,ab}(h_b) \Bigr] +\sum_b m_b\,\Pi_{ab}^{AV}\,\widehat{x}_{ab,i},$$9; at $X_a=m_a/\rho^0_a$0 the inertial range extends to $X_a=m_a/\rho^0_a$1; and at $X_a=m_a/\rho^0_a$2 to $X_a=m_a/\rho^0_a$3 (Cabezón et al., 13 Mar 2025). The paper attributes this result to a combination of integral-derived gradients, dynamic artificial-viscosity switching with linear-field cleaning, pairing-resistant kernels, generalized volume elements, and local Lagrangian compatibility, and explicitly states that no single improvement suffices alone (Cabezón et al., 13 Mar 2025).

6. Astrophysical applications and resolution-driven revisions

SPH-EXA is used in the literature not only as a performance vehicle but as the enabling numerical platform for problem classes that are strongly resolution-limited. The most prominent example is the 2025 tidal disruption event study, which simulated the disruption of a solar-like star by a $X_a=m_a/\rho^0_a$4 black hole from initial disruption to the moment of stream self-intersection, at resolutions up to $X_a=m_a/\rho^0_a$5 particles and including relativistic apsidal precession (Kubli et al., 30 Oct 2025). The simulations were specifically designed to revisit a long-standing debate over whether shock-induced dissipation during the stream’s return to pericentre, before self-intersection, is efficient enough to circularize the debris.

The resolution study reported by Kubli et al. is central to the scientific conclusion. At $X_a=m_a/\rho^0_a$6 particles, corresponding to the highest-resolution SPH simulations in the pre-existing literature, the returning debris stream shows significant in-plane spreading through pericentre, consistent with earlier interpretations of strong dissipation (Kubli et al., 30 Oct 2025). However, as the resolution is increased through $X_a=m_a/\rho^0_a$7, $X_a=m_a/\rho^0_a$8, $X_a=m_a/\rho^0_a$9, $k_a=\sum_b X_bW_{ab}$0, and $k_a=\sum_b X_bW_{ab}$1 particles, the stream width and dissipated-energy fractions converge monotonically, artificial-viscosity tests show that shock dissipation falls by more than two orders of magnitude from $k_a=\sum_b X_bW_{ab}$2 to $k_a=\sum_b X_bW_{ab}$3 particles, and at $k_a=\sum_b X_bW_{ab}$4 particles there is effectively no change between the incoming and outgoing stream widths (Kubli et al., 30 Oct 2025). The paper therefore concludes that the previously inferred strong pericentre-passage shocks and stream widening were numerical artifacts, and that debris circularisation is instead likely mediated by the originally proposed stream–stream collision scenario (Kubli et al., 30 Oct 2025).

A comparable revision of prior expectations appears in the turbulence work. The 2025 turbulence paper states that, for the first time, a purely Lagrangian SPH code can match state-of-the-art mesh and moving-mesh methods in subsonic turbulence, both in slope and inertial-range extent (Cabezón et al., 13 Mar 2025). In combination with the TDE result, this establishes a recurrent pattern in the SPH-EXA literature: extreme-scale execution is not treated as an end in itself, but as the prerequisite for numerical convergence in regimes where earlier SPH calculations were limited either by dissipation, gradient noise, or insufficient particle number. A plausible implication is that SPH-EXA’s scientific significance lies as much in enabling convergence-driven reassessment of established claims as in demonstrating exascale software techniques.