---
title: 'SPH-EXA: Exascale SPH Mini-App'
url: https://www.emergentmind.com/topics/sph-exa
type: topic
---

# SPH-EXA: Exascale SPH Mini-App

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 [2005.02656] [2503.10273]. 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 [1905.03344]. 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 [1809.08013]. 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 [1809.08013].

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 [2005.02656]. 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 [1905.03344].

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 [1809.08013]. 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 [1905.03344]. 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 [2005.02656]. 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 [1905.03344]. 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 [2503.10273].

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})\) [1905.03344]. 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 [2005.02656]. 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 [2503.10273]. 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 [2510.26663].

Gravity support is a defining architectural feature. The mini-app papers inherited multipole gravity from ChaNGa and discussed self-gravity as an optional module [2005.02656]. 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 [2503.10273]. 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 [2510.26663]. 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\) [2005.02656]. 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 [2503.10273].

A central numerical component is the Integral Approach to Derivatives. For each particle \(a\), SPH-EXA solves a \(3\times3\) linear system,
\[
\boldsymbol\nabla f_a=\mathcal{T}_a^{-1}\,\mathbf{I}_a,
\]
where
\[
\mathcal{T}_{ij,a}=\sum_b V_b\,(x_{i,b}-x_{i,a})(x_{j,b}-x_{j,a})\,W_{ab}(h_a),\qquad
I_{i,a}=\sum_b V_b\,f_b\,(x_{i,b}-x_{i,a})\,W_{ab}(h_a),
\]
and uses the matrix-free form
\[
\nabla_i f_a=\sum_b V_b\,(f_b-f_a)\,A_{i,ab}(h_a)
\]
inside the momentum and energy loops [2503.10273]. The same IAD structure appears in the TDE implementation, where it is explicitly credited with requiring only \(O(n_n)\) neighbours while yielding higher accuracy and reduced kernel noise [2510.26663].

Kernel choice is another characteristic feature. SPH-EXA uses the sinc-family
\[
W_{ab}(q)=\frac{K_n}{h_a^d}
\begin{cases}
\bigl[\mathrm{sinc}(\tfrac{\pi}{2}q)\bigr]^n & 0<q\le2,\\
0 & q>2,
\end{cases}
\qquad 3\le n\le 9,
\]
and the turbulence paper states that high-\(n\) sinc kernels suppress pairing up to very high \(k\), requiring fewer neighbours than Wendland kernels [2503.10273]. In the TDE calculations, the chosen kernel is a sixth-order sinc kernel with 100 neighbours [2510.26663].

Artificial viscosity is Monaghan-style, but with modern switching and cleaning procedures. The turbulence paper gives
\[
\Pi_{ab}^{AV}=
\begin{cases}
-\,\tfrac{\alpha_{ab}}{2}\,v_{ab}^{sig}\,w_{ab}, & x_{ab}\!\cdot v_{ab}<0,\\
0, & \text{otherwise,}
\end{cases}
\]
with \(v^{sig}_{ab}=\tfrac12(c_a+c_b)-2w_{ab}\), augmented by a linear field cleaner that replaces the naive velocity jump by a cleaned jump \(v'_a-v'_b\), and by a dynamic individual \(\alpha_a\) that relaxes between \(\alpha_{\min}\simeq0.05\) and \(\alpha_{\max}\simeq1\) through Read–Cullen-type switches [2503.10273]. Generalized volume elements are defined through
\[
X_a=\frac{m_a}{\rho^0_a},\qquad
V_a=\frac{X_a}{\sum_b X_bW_{ab}(h_a)},\qquad
\rho_a=\frac{m_a}{V_a},
\]
and are coupled to an Atwood-number-based switch \(\sigma_{ab}\) with \(At_{\min}=0.1\) and \(At_{\max}=0.2\) to suppress tensile instability while preserving conservation elsewhere [2503.10273]. For the TDE problem, the hydrodynamics is supplemented by a polytropic equation of state \(P=K\rho^\gamma\) with \(\gamma=5/3\), and relativistic apsidal precession is represented with a pseudo-Newtonian “Einstein” potential
\[
\Phi(r)=-\frac{G\,M_\bullet}{r}\Bigl(1+3\,\frac{r_g}{r}\Bigr),\qquad r_g=\frac{G\,M_\bullet}{c^2}
\]
evaluated on the GPU [2510.26663].

## 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 \(0.08\,\mathrm{s}\) per timestep on 1 node and approximately \(0.06\,\mathrm{s}\) on 16 nodes, while the parent codes required approximately \(0.5\)–\(0.7\,\mathrm{s}\) on 1 node and converged to \(0.2\,\mathrm{s}\) on 16 nodes, implying an approximately \(6\times\) speedup of SPH-EXA over the legacy codes at all scales [1905.03344]. The same study reports strong-scaling efficiency \(\eta_{\rm strong}=T_1/(p\,T_p)\) dropping from approximately \(0.90\) at 4 nodes to approximately \(0.30\) at 32 nodes, with good scaling up to 16 nodes and deterioration once \(N_{\rm loc}\) falls below approximately \(10^4\) particles per core; weak-scaling efficiency remains approximately \(1.0\) up to 256 nodes and approximately \(0.3\) at 512 nodes for constant \(N_{\rm loc}=3\times10^4\) per rank [1905.03344].

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 \(65{,}536\) 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 [2005.02656]. 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 [2005.02656].

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 [1905.03344]. 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 \(O(10)\,\mathrm{TFLOP/s}\) per GPU in double precision [2503.10273]. In the TDE work, SPH-EXA is reported to scale strongly and weakly to \(2.8\times10^{12}\) particles across thousands of GPUs on LUMI-G and ALPS, with peak per-GPU performance between \(1\times10^7\) and \(2.5\times10^7\) particle-timesteps per second on a GH200 GPU and a memory footprint of approximately 80–100 bytes per particle [2510.26663].

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 [1905.03344]. The 2020 review also lists fault tolerance, silent corruption detection, optimal checkpoint intervals, and mixed precision with error control among the principal future extensions [2005.02656].

## 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 [1807.11244]. With the filter
\[
\widetilde{\phi}(\mathbf{x})=\int_{\Omega}\phi(\mathbf{r})W_h(\mathbf{x}-\mathbf{r})\,d\mathbf{r}
\]
and particle motion at the filtered velocity
\[
\frac{d\mathbf{x}}{dt}=\widetilde{\mathbf{u}}_h(\mathbf{x}),
\]
the filtered momentum equation contains the subgrid stress
\[
\tau_{ij}(\mathbf{x})=
\int \rho(\mathbf{r}')
\bigl[u_i(\mathbf{r}')-\widetilde{u}_i(\mathbf{x})\bigr]
\bigl[u_j(\mathbf{r}')-\widetilde{u}_j(\mathbf{x})\bigr]
W_h(\mathbf{x}-\mathbf{r}')\,d\mathbf{r}',
\]
equivalently \(\tau_{ij}=\langle \rho u_i u_j\rangle-\langle \rho_h\rangle \widetilde{u}_i\widetilde{u}_j\) [1807.11244]. 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 [1807.11244].

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 \(\tau_{ij}\) 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 [1807.11244]. 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 \(E_0\) 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 [2503.10273]. At \(250^3\) particles it recovers a clear \(k^{-5/3}\) slope from \(k\approx5\) to \(k\approx25\); at \(500^3\) the inertial range extends to \(k\in[5,50]\); and at \(1000^3\) to \(k\in[5,100]\) [2503.10273]. 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 [2503.10273].

## 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 \(10^6\,M_\odot\) black hole from initial disruption to the moment of stream self-intersection, at resolutions up to \(10^{10}\) particles and including relativistic apsidal precession [2510.26663]. 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 \(10^8\) 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 [2510.26663]. However, as the resolution is increased through \(10^6\), \(1.6\times10^7\), \(1.28\times10^8\), \(5.12\times10^8\), and \(10^{10}\) 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 \(1.6\times10^7\) to \(10^{10}\) particles, and at \(10^{10}\) particles there is effectively no change between the incoming and outgoing stream widths [2510.26663]. 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 [2510.26663].

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 [2503.10273]. 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.

Source: https://www.emergentmind.com/topics/sph-exa