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

Technical Deep Dive: Smoothed Particle Hydrodynamics in pkdgrav3 for Shock Physics Simulations

Problem and Context

The paper presents a high-performance implementation of Smoothed Particle Hydrodynamics (SPH) within the pkdgrav3 N-body code. SPH is a Lagrangian particle-based method well-suited for simulating highly dynamic, strongly deforming flows, making it ideal for modeling planetary-scale collisions where tracking the motion and mixing of individual material parcels is essential.

Key challenges addressed include:

• Efficiently finding neighboring particles for SPH calculations
• Coupling SPH with the Fast Multipole Method (FMM) gravity solver in pkdgrav3
• Enabling vectorized particle-particle interaction processing on both CPUs and GPUs
• Achieving excellent performance and scaling to take advantage of modern HPC systems

Methodology

SPH Implementation in pkdgrav3

• Adopts a standard SPH formulation derived from the Lagrangian
• Uses the cubic spline and Wendland C2/C6 kernel functions
• Implements an efficient neighbor search method that reuses the pkdgrav3 FMM tree structure
• Includes an interface correction to address density discontinuities at material interfaces
• Employs a hierarchical kick-drift-kick leapfrog scheme for time integration
• Supports both evolved and smoothed density formulations

Equation of State

• Provides built-in support for the ideal gas EOS
• Integrates the generalized EOSlib interface for more advanced EOS models (e.g. Tillotson, ANEOS, REOS.3, SCvH)

Performance Optimizations

• Vectorizes particle-particle interactions on both CPUs and GPUs
• Dynamically adjusts bucket/group sizes to handle steep density gradients
• Leverages pkdgrav3's highly scalable parallel architecture with MPI, OpenMP, and SIMD

Results and Validation

The authors demonstrate the accuracy and performance of their SPH implementation through a comprehensive suite of standard hydrodynamic test cases:

1. Linear Traveling Sound Wave: Shows SPH convergence scaling close to first-order, limited by floating-point precision at high resolutions.
2. Sod Shock Tube: High-resolution simulations closely match the analytic solution, with the pressure blip at the contact discontinuity greatly reduced.
3. Sedov-Taylor Blast Wave: Captures the shock front and post-shock features accurately at high resolutions, with good energy conservation.
4. Evrard Collapse: Matches 1D reference solutions at high resolutions, demonstrating the coupling between SPH and gravity.
5. Gresho-Chan Vortex: Highlights the importance of artificial viscosity settings, with lower values producing better results.
6. Hydrostatic Equilibrium: Shows that the Wendland C6 kernel performs better than the cubic spline in preserving sharp interfaces.
7. Kelvin-Helmholtz Instability: Demonstrates the impact of artificial viscosity and surface tension on the development of fluid instabilities.

Performance and Scaling

The authors evaluate the performance and scaling of their SPH implementation on the Piz Daint supercomputer:

• Single-node scaling follows the expected O(N^4/3) runtime complexity, with the tree-based neighbor search and gravity computations dominating the runtime.
• Strong scaling results show excellent parallel efficiency, above 80% up to 512 nodes for a 2 billion particle simulation.
• Weak scaling experiments maintain over 75% parallel efficiency with 8 million particles per node.

Conclusion

The authors have developed a high-performance SPH implementation within the pkdgrav3 code that can efficiently leverage modern HPC resources, including hybrid CPU-GPU architectures. The validation tests demonstrate the accuracy and robustness of their approach, making pkdgrav3 a powerful platform for astrophysical and planetary applications involving highly dynamic, self-gravitating systems.