Smoothed Particle Hydrodynamics in pkdgrav3 for Shock Physics Simulations. I. Hydrodynamics

Christian Reinhardt; Douglas Potter; Joachim Stadel; Thomas Meier

arxiv: 2511.11425 · v5 · submitted 2025-11-14 · 🌌 astro-ph.EP · astro-ph.IM

Smoothed Particle Hydrodynamics in pkdgrav3 for Shock Physics Simulations. I. Hydrodynamics

Thomas Meier , Douglas Potter , Christian Reinhardt , Joachim Stadel This is my paper

Pith reviewed 2026-05-17 22:09 UTC · model grok-4.3

classification 🌌 astro-ph.EP astro-ph.IM

keywords Smoothed Particle Hydrodynamicstree codeself-gravityparallel computingplanetary impactsshock physicshydrodynamic simulationsGPU acceleration

0 comments

The pith

pkdgrav3 is a parallel tree-SPH code that scales large hydrodynamic simulations with self-gravity to thousands of cores and GPUs.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper introduces pkdgrav3 as a fully parallel code that pairs a hierarchical tree for gravity and neighbor searches with a modern SPH hydrodynamics solver. It targets large-scale problems such as planetary impacts where particles must track material identity and thermodynamic history through strong shocks. The implementation uses a hybrid shared and distributed memory model plus asynchronous communication to run on mixed CPU and GPU clusters. Validation against a suite of standard hydrodynamic tests shows close agreement with analytic and reference solutions. The authors note that the code has already been applied to peer-reviewed studies of self-gravitating impact events.

Core claim

pkdgrav3 combines an efficient hierarchical tree algorithm for gravity and neighbor finding with a modern implementation of Smoothed Particle Hydrodynamics optimized for massively parallel hybrid CPU/GPU architectures. Its hybrid shared/distributed memory model with asynchronous communication scales efficiently to thousands of CPU cores and GPUs. The code demonstrates excellent agreement with analytic or reference solutions across a suite of standard tests and has been used to model planetary-scale impacts where SPH's Lagrangian nature tracks material origin and thermodynamic evolution.

What carries the argument

The hierarchical tree algorithm for gravity and neighbor finding paired with an SPH solver on a hybrid shared/distributed memory model with asynchronous communication.

If this is right

The Lagrangian particle approach allows direct tracking of material origin and thermodynamic states through strong shocks in impact events.
The hybrid memory and asynchronous scheme supports efficient use of modern heterogeneous supercomputers for self-gravitating flows.
The same framework can address other astrophysical problems that require both hydrodynamics and self-gravity on large particle counts.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

If the scaling and accuracy claims hold, the code could reduce the cost of resolving fine-scale mixing and ejecta in giant impacts compared with fixed-grid methods.
The tree-based neighbor search may lend itself to adaptive resolution in regions of high density contrast without manual remeshing.
Integration with additional physics such as strength models or radiation could follow directly from the existing parallel infrastructure.

Load-bearing premise

That agreement on a standard test suite is sufficient to guarantee accuracy in highly dynamical self-gravitating planetary impact simulations at full scale without hidden implementation errors.

What would settle it

A side-by-side comparison of pkdgrav3 output against an independent high-resolution reference solution for a self-gravitating Sedov blast or shock-tube problem at planetary scales.

Figures

Figures reproduced from arXiv: 2511.11425 by Christian Reinhardt, Douglas Potter, Joachim Stadel, Thomas Meier.

**Figure 1.** Figure 1: Left: 2D illustration of the particles (gray dots) in the cell bound (black) with their kernel balls (green) and the resulting ball bounds, the box-of-balls in blue and the ball-of-balls in orange. Center: Combining two cells in the tree build results in combined bounds. The original cell bounds are gray solid lines, while the combined cell bound is shown by a black dashed line. The combination of the boxe… view at source ↗

**Figure 2.** Figure 2: Illustration of how a single time step is done (the kick phase is indicated by the yellow box). First, the density update is done for all particles (or the active particles and their neighbors in the case of FastGas, see text). Then, using the time derivatives from the last step, the velocity and internal energy are predicted (green arrows). Using these values, the force calculation, the closing kick of th… view at source ↗

**Figure 3.** Figure 3: L1 error norm of the linear traveling sound wave problem. In blue, the body centered cubic grid shows N −0.898 scaling with the number of grid cells in the direction of wave propagation, while in orange the primitive cubic grid shows N −0.769 scaling. As expected for traditional SPH, the scaling falls short of first-order and becomes dominated by the floating-point precision used by the code above N = 1024… view at source ↗

**Figure 5.** Figure 5: Results of the Sedov-Taylor blast wave test for different resolutions. The top frame shows the density while the bottom frame shows the radial velocity. All values are binned in 100 bins, with the marker showing the average value in the bin while the error bars show maximum and minimum values. The Sedov-Taylor blast wave (L. I. Sedov 1946; G. I. Taylor 1950a,b; L. I. Sedov 1959) is a spherical shock wave g… view at source ↗

**Figure 6.** Figure 6: Results of the Evrard collapse test for different resolutions. The top frame shows the density, the middle frame the radial velocity and the bottom frame the entropy. All values are binned in 100 bins, with the marker showing the average value in the bin while the error bars show maximum and minimum values. the Sedov-Taylor blast wave test, we use a global time step smaller than the smallest needed by any… view at source ↗

**Figure 7.** Figure 7: Top: Time evolution of the potential, kinetic, internal and total energy in Evrard’s spherical collapse test for different resolutions. The simulation with N = 1.12 billion particles is only run up to t = 0.8. Bottom: The total energy is conserved to within 0.06 %, even during shock formation, with conservation improving at higher resolution. perfectly balanced by the pressure gradient. Artificial viscosi… view at source ↗

**Figure 9.** Figure 9: Results of the hydrodynamic equilibrium test. Shown is the density, calculated by distributing the mass of each particle to each pixel according to the intersection between the pixel area and the particle’s kernel using the SPH kernel function. The left column shows the results of the equal-mass test, while the right column shows the results of the equal-spacing test. In each block, the results are for N … view at source ↗

**Figure 10.** Figure 10: Results from the Blob test. Each row shows a timeseries, from left to right: 0.25 τKH, 1.0 τKH, 1.75 τKH, 2.5 τKH and 4.0 τKH, where the top half of the frame is from the simulation using the cubic spline kernel with 32 neighbors and the bottom half is from the simulation using the Wendland C6 kernel with 400 neighbors. From top to bottom, the resolution increases, with the short sides of the domain being… view at source ↗

**Figure 11.** Figure 11: The evolution of the cloud mass fraction in the blob test. The simulations with the cubic spline kernel and 32 neighbors are shown as solid lines, while the runs with the Wendland C6 kernel and 400 neighbors are shown as dashed lines. tions for both the cubic spline kernel with 32 neighbors and the Wendland C6 kernel with 400 neighbors. At low resolutions, we get similar results as O. Agertz et al. (2007)… view at source ↗

**Figure 13.** Figure 13: Results of the Kelvin-Helmholtz instability test with a density ratio of χ = 1 after (from left to right) 0.3 τKH, 0.6 τKH, 1.0 τKH and 2.0 τKH. The top half of each row shows the results of the runs with α = 1.0 and β = 2.0, while the lower half shows the runs with α = 0.1 and β = 0.0. From top to bottom, the resolution increases, with the simulation x-axis being sampled by 512, 1024, 2048, and 4096 BCC … view at source ↗

**Figure 14.** Figure 14: Scaling of time per step on a single node. The total time per step closely follows the expected O [PITH_FULL_IMAGE:figures/full_fig_p018_14.png] view at source ↗

**Figure 15.** Figure 15: Strong scaling results with 2048 million particles. The top frame shows the time per step and the bottom frame shows the parallel efficiency. The values for 1, 2 and 4 nodes in the blue shaded region are extrapolated from the result for 8 nodes, as the simulation would not fit into memory at those node counts. tested particle counts. The last three data points are extrapolated from runs with more nodes … view at source ↗

**Figure 18.** Figure 18: Weak scaling results with 8 million particles per node. The top frame shows the time per step and the bottom frame shows the parallel efficiency. particles being run on 1024 nodes. But as shown in Figure 15, at this node count, the result has already visibly departed from perfect scaling. A production simulation would thus never be run in this configuration, but rather with a more realistic 256 nodes. Th… view at source ↗

**Figure 17.** Figure 17: Weak scaling results with 2 million particles per node. The top frame shows the time per step and the bottom frame shows the parallel efficiency. We also present weak scaling results. For this, we use the same models from Section 4.1 again, but we choose the global step size such that it satisfies the CFL condition for all resolutions directly, removing the need for the code to do substeps and thus also … view at source ↗

read the original abstract

We present pkdgrav3, a high-performance, fully parallel tree-SPH code designed for large-scale hydrodynamic simulations including self-gravity. Building upon the long development history of pkdgrav, the code combines an efficient hierarchical tree algorithm for gravity and neighbor finding with a modern implementation of Smoothed Particle Hydrodynamics (SPH) optimized for massively parallel hybrid CPU/GPU architectures. Its hybrid shared/distributed memory model, combined with an asynchronous communication scheme, allows pkdgrav3 to scale efficiently to thousands of CPU cores and GPUs. We validate the numerical accuracy of pkdgrav3 using a suite of standard tests, demonstrating excellent agreement with analytic or reference solutions. The code was already used in several peer-reviewed publications to model planetary-scale impacts, where SPH's Lagrangian nature allows accurate tracking of material origin and thermodynamic evolution. These examples highlight pkdgrav3's robustness and efficiency in simulating highly dynamical, self-gravitating systems. pkdgrav3 thus provides a powerful, flexible, and scalable platform for astrophysical and planetary applications, capable of exploiting the full potential of modern heterogeneous high-performance computing systems.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Desk Editor's Note private letter to a colleague

pkdgrav3 adds a working SPH module to the established pkdgrav tree code for hybrid CPU/GPU runs, with standard validation and some prior use in impact papers.

read the letter

The main takeaway is that this group has folded a modern SPH implementation into pkdgrav3, keeping the tree for gravity and neighbor finding while adding hybrid memory handling and async communication to run on large CPU/GPU clusters. It builds directly on the long pkdgrav lineage rather than starting from scratch, which is the practical move here. The code already shows up in a few published planetary impact studies, and the abstract says the usual hydro tests line up with analytic or reference solutions. That combination gives it some real-world grounding beyond just the test suite. The scaling claims to thousands of cores and GPUs sound plausible given the architecture choices, and the Lagrangian tracking of material is a natural fit for the impact work they mention. The soft spot is exactly the one the stress-test note flags: standard shock-tube and blast-wave tests do not strongly exercise the coupled tree-SPH behavior in the messy, self-gravitating regimes that matter for planetary collisions. Neighbor finding, artificial viscosity, and asynchronous hand-offs could hide small conservation or mixing errors that only appear at scale or in highly dynamical cases. The prior publications help, but they do not replace a clear set of coupled gravity-hydro benchmarks in the paper itself. This is a tool paper aimed at computational planetary scientists and astrophysicists who already run or want to run large SPH-plus-gravity simulations on heterogeneous hardware. Readers who need an alternative to existing codes and care about scaling on modern machines will get the most out of it. It is solid enough on implementation and prior use that a serious referee should see it; the details on the SPH kernel choices and any extra coupled tests would be the main things to check. I would send it to peer review.

Referee Report

1 major / 2 minor

Summary. The manuscript presents pkdgrav3, a high-performance fully parallel tree-SPH code for large-scale hydrodynamic simulations including self-gravity. It builds on the pkdgrav framework by combining an efficient hierarchical tree algorithm for gravity and neighbor finding with a modern SPH implementation optimized for hybrid CPU/GPU architectures and asynchronous communication to enable scaling to thousands of cores. Numerical accuracy is validated via a suite of standard tests that demonstrate excellent agreement with analytic or reference solutions. The code has already been applied in peer-reviewed publications to model planetary-scale impacts, where the Lagrangian nature of SPH enables accurate material tracking and thermodynamic evolution.

Significance. If the implementation and validation hold, this work supplies a scalable, flexible platform for astrophysical and planetary simulations on modern heterogeneous HPC systems. The reuse of established tree infrastructure and the demonstrated prior application to self-gravitating impacts are concrete strengths that could benefit the community working on shock physics and large-scale dynamical problems.

major comments (1)

Validation section: the manuscript relies on a suite of standard hydrodynamics tests (shock tubes, blast waves, hydrostatic equilibria) that do not strongly couple self-gravity with SPH in highly dynamical, multi-scale regimes. Because the central claim includes robustness for planetary impact simulations, where tree-SPH neighbor finding, artificial viscosity, and asynchronous communication could affect conservation or material tracking, explicit tests or quantitative metrics for the coupled gravity-hydro system are needed to substantiate the accuracy assertions.

minor comments (2)

Implementation section: specify the exact SPH formulation (e.g., artificial viscosity switch parameters, kernel choice, and time-stepping criteria) to allow direct reproduction of the reported test results.
Figures: ensure all test-result panels include quantitative error measures or direct overlays against analytic solutions rather than qualitative visual agreement alone.

Simulated Author's Rebuttal

1 responses · 0 unresolved

We thank the referee for their constructive review and positive assessment of the manuscript's significance. We address the single major comment below and agree that strengthening the presentation of coupled gravity-hydrodynamics validation will improve the paper.

read point-by-point responses

Referee: Validation section: the manuscript relies on a suite of standard hydrodynamics tests (shock tubes, blast waves, hydrostatic equilibria) that do not strongly couple self-gravity with SPH in highly dynamical, multi-scale regimes. Because the central claim includes robustness for planetary impact simulations, where tree-SPH neighbor finding, artificial viscosity, and asynchronous communication could affect conservation or material tracking, explicit tests or quantitative metrics for the coupled gravity-hydro system are needed to substantiate the accuracy assertions.

Authors: We agree that the primary validation suite emphasizes isolated hydrodynamics tests to verify the SPH implementation. Self-gravity is computed with the established hierarchical tree algorithm inherited from the pkdgrav framework, whose gravitational accuracy has been documented in prior literature. The full coupled system is exercised in the planetary impact applications already cited in the manuscript, where conservation of energy, momentum, and material tracking were implicitly confirmed by the physical outcomes. To make this explicit, we will revise the validation section to include a new subsection presenting quantitative results from at least one self-gravitating dynamical test (e.g., a Jeans instability or collapsing cloud problem). We will report conservation errors, comparison against reference solutions, and brief discussion of how tree-based neighbor finding and asynchronous communication affect accuracy in multi-scale regimes. These additions will directly address the referee's concern without altering the core hydrodynamics focus of Part I. revision: yes

Circularity Check

0 steps flagged

No circularity: validation against external benchmarks

full rationale

This is a software implementation and validation paper for the pkdgrav3 tree-SPH code. The central claims rest on describing the code architecture and showing agreement with independent analytic or reference solutions from a standard test suite (shock tubes, blast waves, hydrostatic equilibria). No mathematical derivations, fitted parameters renamed as predictions, or self-referential equations appear. Self-citations, if present, support prior code history but are not load-bearing for any result; the validation is externally falsifiable against published benchmarks outside the paper's own fitted values or implementations. The derivation chain is therefore self-contained with no reductions to inputs by construction.

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 0 invented entities

The paper is a code presentation rather than a derivation of new physics, so the central claim rests primarily on the correctness of the software implementation and standard SPH assumptions rather than new free parameters or invented entities.

axioms (1)

domain assumption Standard SPH hydrodynamics equations and artificial viscosity formulations hold for the intended shock physics regime.
Invoked implicitly when claiming validation on standard tests for hydrodynamic accuracy.

pith-pipeline@v0.9.0 · 5504 in / 1219 out tokens · 35794 ms · 2026-05-17T22:09:56.424174+00:00 · methodology

discussion (0)

Lean theorems connected to this paper

Citations machine-checked in the Pith Canon. Every link opens the source theorem in the public Lean library.

IndisputableMonolith/Foundation/DimensionForcing.lean; IndisputableMonolith/Cost/FunctionalEquation.lean reality_from_one_distinction; washburn_uniqueness_aczel; alexander_duality_circle_linking unclear

?

unclear
Relation between the paper passage and the cited Recognition theorem.

We validate the numerical accuracy of pkdgrav3 using a suite of standard tests, demonstrating excellent agreement with analytic or reference solutions... tree-SPH code... FMM gravity solver... neighbor search... artificial viscosity

What do these tags mean?

matches: The paper's claim is directly supported by a theorem in the formal canon.
supports: The theorem supports part of the paper's argument, but the paper may add assumptions or extra steps.
extends: The paper goes beyond the formal theorem; the theorem is a base layer rather than the whole result.
uses: The paper appears to rely on the theorem as machinery.
contradicts: The paper's claim conflicts with a theorem or certificate in the canon.
unclear: Pith found a possible connection, but the passage is too broad, indirect, or ambiguous to say the theorem truly supports the claim.

Reference graph

Works this paper leans on

4 extracted references · 4 canonical work pages

[1]

Mignone and Jonathan C

Agertz, O., Moore, B., Stadel, J., et al. 2007, Monthly Notices of the Royal Astronomical Society, 380, 963, doi: 10.1111/j.1365-2966.2007.12183.x Alonso Asensio, I., Dalla Vecchia, C., Potter, D., & Stadel, J. 2023, Monthly Notices of the Royal Astronomical Society, 519, 300, doi: 10.1093/mnras/stac3447 Asphaug, E., Collins, G., & Jutzi, M. 2015, in Aste...

work page doi:10.1111/j.1365-2966.2007.12183.x 2007
[2]

https: //ui.adsabs.harvard.edu/abs/1946JApMM..10..241S Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics (Academic Press). https://ui.adsabs.harvard.edu/abs/1959sdmm.book.....S Sod, G. A. 1978, Journal of Computational Physics, 27, 1, doi: 10.1016/0021-9991(78)90023-2 Springel, V. 2010, Monthly Notices of the Royal Astronomical Society, 4...

work page doi:10.1016/0021-9991(78)90023-2 1959
[3]

https://ui.adsabs.harvard.edu/abs/1993A&A...268..391S Stewart, S. T. 2020, Zenodo, doi: 10.5281/zenodo.3866507 Stewart, S. T., Davies, E. J., Duncan, M. S., et al. 2019, Zenodo, doi: 10.5281/zenodo.3478631 Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, The Astrophysical Journal Supplement Series, 178, 137, doi: 10.1086/5887...

work page doi:10.5281/zenodo.3866507 2020
[4]

http://www.gnu.org/s/parallel Taylor, G. I. 1950a, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 201, 159, doi: 10.1098/rspa.1950.0049 Taylor, G. I. 1950b, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 201, 175, doi: 10.1098/rspa.1950.0050 Teyssier, R. 2002, Astronomy ...

work page doi:10.1098/rspa.1950.0049 1950

[1] [1]

Mignone and Jonathan C

Agertz, O., Moore, B., Stadel, J., et al. 2007, Monthly Notices of the Royal Astronomical Society, 380, 963, doi: 10.1111/j.1365-2966.2007.12183.x Alonso Asensio, I., Dalla Vecchia, C., Potter, D., & Stadel, J. 2023, Monthly Notices of the Royal Astronomical Society, 519, 300, doi: 10.1093/mnras/stac3447 Asphaug, E., Collins, G., & Jutzi, M. 2015, in Aste...

work page doi:10.1111/j.1365-2966.2007.12183.x 2007

[2] [2]

https: //ui.adsabs.harvard.edu/abs/1946JApMM..10..241S Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics (Academic Press). https://ui.adsabs.harvard.edu/abs/1959sdmm.book.....S Sod, G. A. 1978, Journal of Computational Physics, 27, 1, doi: 10.1016/0021-9991(78)90023-2 Springel, V. 2010, Monthly Notices of the Royal Astronomical Society, 4...

work page doi:10.1016/0021-9991(78)90023-2 1959

[3] [3]

https://ui.adsabs.harvard.edu/abs/1993A&A...268..391S Stewart, S. T. 2020, Zenodo, doi: 10.5281/zenodo.3866507 Stewart, S. T., Davies, E. J., Duncan, M. S., et al. 2019, Zenodo, doi: 10.5281/zenodo.3478631 Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, The Astrophysical Journal Supplement Series, 178, 137, doi: 10.1086/5887...

work page doi:10.5281/zenodo.3866507 2020

[4] [4]

http://www.gnu.org/s/parallel Taylor, G. I. 1950a, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 201, 159, doi: 10.1098/rspa.1950.0049 Taylor, G. I. 1950b, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 201, 175, doi: 10.1098/rspa.1950.0050 Teyssier, R. 2002, Astronomy ...

work page doi:10.1098/rspa.1950.0049 1950