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

This paper introduces pkdgrav3, a high-performance, fully parallel tree-SPH code optimized for hybrid CPU/GPU architectures that accurately simulates large-scale, self-gravitating hydrodynamic systems such as planetary impacts.

The Gist

Imagine you are trying to simulate a cosmic event, like two planets smashing into each other, or a star collapsing under its own weight. To do this, you need a computer program that can track billions of tiny pieces of matter as they fly, crash, heat up, and cool down.

This paper introduces pkdgrav3, a super-powered computer program designed to do exactly that. Think of it as the "Ferrari" of astrophysics simulation software. Here is a breakdown of what it does and why it's special, using some everyday analogies.

1. The Problem: Simulating a Cosmic Dance

In the past, scientists had two main ways to simulate fluids (like gas, magma, or molten rock):

• The Grid Method: Imagine a giant chessboard covering the universe. You calculate what happens in each square. The problem? If a planet moves across the board, the squares don't move with it. It's like trying to track a runner by only looking at the tiles on the floor; you lose track of them if they run off the edge of a tile.
• The Particle Method (SPH): Instead of a grid, imagine the fluid is made of millions of tiny marbles. Each marble carries its own temperature, speed, and pressure. As they move, they carry their properties with them. This is Smoothed Particle Hydrodynamics (SPH). It's perfect for things that change shape wildly, like an exploding star or a crashing planet, because the "marbles" naturally flow wherever the material goes.

The Catch: Traditional SPH codes are slow. If you want to simulate a planet with high detail, you need billions of marbles. Calculating how billions of marbles interact with each other is like trying to introduce every person at a stadium of 100,000 people to every other person. It takes forever.

2. The Solution: pkdgrav3

The authors built pkdgrav3 to solve this speed problem. They took an existing, highly efficient gravity simulator (which calculates how things pull on each other) and added a brand-new, high-speed engine for fluid dynamics (SPH).

Here is how they made it fast, using some metaphors:

A. The "Smart Neighborhood" Search

In a standard simulation, every particle has to ask, "Who are my neighbors?" to calculate forces. If you have a billion particles, asking everyone individually is a nightmare.

• The pkdgrav3 Trick: They use a Tree Structure. Imagine organizing a city not by asking every house who lives next door, but by grouping them into blocks, then neighborhoods, then districts.
• The Analogy: If you want to know who lives near your house, you don't knock on every door in the city. You look at your block. If the block is far away, you ignore it entirely. pkdgrav3 uses a "Fast Multipole Method" (FMM) to do this instantly. It groups particles into "buckets" and only does the heavy math for buckets that are actually close. This turns a task that would take years into one that takes minutes.

B. The "Hybrid Engine" (CPU + GPU)

Modern supercomputers have two types of processors:

• CPUs: The "brains." They are good at complex, step-by-step thinking.
• GPUs: The "muscle." They are good at doing the same simple calculation millions of times at once (like painting a wall with a roller vs. a brush).
• The pkdgrav3 Trick: This code is built to use both simultaneously. It splits the work: the "brains" manage the big picture and the tree structure, while the "muscle" (GPUs) crunches the billions of particle interactions in parallel. It's like having a conductor leading an orchestra where every musician plays their part perfectly in sync.

C. The "Asynchronous" Workflow

Usually, in a computer simulation, everyone waits for the slowest person to finish before moving to the next step.

• The pkdgrav3 Trick: They use an asynchronous system. Imagine a relay race where runners don't wait for the baton to be perfectly handed off; they keep running as soon as they have the info they need. This means the computer never sits idle waiting for data. It keeps chugging along, making the simulation incredibly efficient.

3. Why Does This Matter? (The "Shock Physics")

The paper focuses on Shock Physics. This is the study of what happens when things hit each other really hard—like a meteorite hitting Earth or planets colliding in the early solar system.

• The Challenge: When things crash, they create shockwaves, heat up to thousands of degrees, and change state (solid to liquid to gas). Standard simulations often get this wrong, creating "fake" friction that stops the fluids from mixing naturally.
• The Fix: pkdgrav3 includes special math tricks (like "interface corrections") to handle the messy boundaries between different materials (e.g., rock vs. gas) without creating these fake forces. It allows the simulation to show how a moon might form from the debris of a giant collision, tracking exactly which piece of rock came from the impactor and which came from the target planet.

4. The Results: Superpowers

The authors tested their code with a battery of standard physics puzzles (like shock tubes and exploding waves) and found it matched the theoretical answers perfectly.

But the real magic is in the Scale:

• Resolution: They can simulate 2.1 billion particles in a single run. To put that in perspective, previous simulations of similar events might have used 10 million. This is like going from a low-resolution pixelated image to a 4K ultra-HD movie. You can finally see tiny details like planetary crusts or thin atmospheres.
• Speed: They ran this massive simulation on a supercomputer in Switzerland (Piz Daint) using 256 GPU nodes, and it finished in a reasonable amount of time.

Summary

pkdgrav3 is a new, ultra-fast tool that lets scientists simulate the violent, chaotic birth of planets and stars with unprecedented detail.

• Old way: Slow, blurry, and limited to small models.
• New way (pkdgrav3): Fast, crystal clear, and capable of simulating entire planetary systems with billions of "marbles," all while running efficiently on modern supercomputers.

It's like upgrading from a flip phone to a smartphone: suddenly, you can do things that were previously impossible, opening up new ways to understand how our universe was built.

Technical Summary

Here is a detailed technical summary of the paper "Smoothed Particle Hydrodynamics in pkdgrav3 for Shock Physics Simulations I: Hydrodynamics" by Meier et al.

1. Problem Statement

Hydrodynamical simulations are critical for modeling diverse astrophysical phenomena, from galaxy formation to planetary impacts. While grid-based methods offer high efficiency in structured domains, they struggle with large density contrasts, free surfaces, and complex geometries due to advection errors and the difficulty of adaptive mesh refinement. Particle-based methods, specifically Smoothed Particle Hydrodynamics (SPH), offer a fully Lagrangian, adaptive approach ideal for these scenarios but suffer from well-known numerical issues:

• Artificial Viscosity: Required for shock capturing, it often introduces excessive dissipation in smooth flows.
• Surface Tension: Traditional SPH formulations suppress fluid mixing and hydrodynamic instabilities at contact discontinuities due to spurious surface tension.
• Scalability: Modern High-Performance Computing (HPC) architectures rely heavily on hybrid CPU/GPU systems. Many existing SPH codes do not scale efficiently to thousands of nodes or fully exploit GPU acceleration.

The authors aim to address these challenges by implementing a high-performance, fully parallel SPH module within pkdgrav3, a code already renowned for its efficient N-body gravity solver, to enable large-scale, self-gravitating shock physics simulations.

2. Methodology

The paper details the integration of a modern SPH solver into the pkdgrav3 framework, leveraging its existing hierarchical tree structure and hybrid parallelization strategies.

Code Architecture & Parallelism

• Hybrid Model: The code utilizes a shared/distributed memory model combining MPI (distributed), Boost threads (shared), and SIMD/SIMT vectorization.
• GPU Acceleration: Data is converted from an Array of Structures (AoS) to an Array of Structures of Arrays (AoSoA) format to optimize memory access for CUDA warps (32 threads).
• Tree Structure: The simulation domain is decomposed via Orthogonal Recursive Bisection (ORB) into a binary tree. This tree is reused for both gravity (Fast Multipole Method - FMM) and neighbor finding, ensuring memory balance and efficient load distribution.

SPH Implementation Details

• Formulation: The authors use the standard density-energy formulation derived from the Lagrangian (Springel & Hernquist 2002), which is compatible with common Equations of State (EOS) for planetary materials (e.g., Tillotson, ANEOS, M-ANEOS).
• Neighbor Search: A novel approach reuses the FMM gravity tree for neighbor finding. Instead of a particle-by-particle search, the code walks the tree for groups of particles (sinks) using "box-of-balls" and "ball-of-balls" bounds. This reduces the complexity and communication overhead, allowing for vectorized evaluation of interaction lists.
• Interface Correction: To mitigate spurious surface tension at material interfaces (e.g., planet cores vs. mantles), the code implements the correction method by Ruiz-Bonilla et al. (2022). This uses a kernel imbalance metric to estimate and correct temperature and pressure at interfaces before calculating density.
• Time Integration: A hierarchical kick-drift-kick leapfrog scheme is employed. It includes a predictor-corrector step to handle velocity-dependent hydrodynamic forces. The code supports individual time steps (substeps) but enforces constraints to prevent energy conservation errors when particles on vastly different time-step rungs interact.
• Viscosity: Standard artificial viscosity (Monaghan 1992) is used for shock capturing. The authors note that while viscosity switches exist, they often struggle with shock capturing; instead, they rely on high resolution to mitigate non-shock dissipation.

3. Key Contributions

1. High-Performance SPH Module: The first implementation of a fully coupled SPH module within pkdgrav3, optimized for hybrid CPU/GPU architectures.
2. Optimized Neighbor Search: A tree-walking algorithm that processes groups of particles simultaneously, significantly reducing the computational cost of neighbor finding compared to traditional particle-by-particle searches.
3. Robust Interface Handling: Integration of the interface correction method, allowing for more accurate modeling of multi-material planetary collisions without the severe surface tension artifacts typical of standard SPH.
4. Scalability: Demonstration of the code's ability to scale to thousands of CPU cores and GPUs, handling simulations with billions of particles.
5. Open Source: The code is released under GPLv3, with source code and documentation available, promoting reproducibility in planetary science.

4. Results

The authors validated the code through a comprehensive suite of hydrodynamic tests and performance benchmarks on the Piz Daint supercomputer (CSCS, Switzerland).

Hydrodynamic Tests

• Linear Sound Wave: Demonstrated convergence rates close to first-order ($N^{-0.898}$) for Body Centered Cubic (BCC) grids, consistent with theoretical expectations for SPH.
• Sod Shock Tube: High-resolution simulations ($N=2048$) reproduced the analytical solution nearly perfectly, with the "pressure blip" at the contact discontinuity vanishing at high resolution.
• Sedov-Taylor Blast Wave: The code successfully captured the spherical shock wave. Energy conservation was maintained below 0.2% when using a global time step, though individual time steps introduced larger errors (a known trade-off).
• Evrard Collapse: A self-gravitating collapse test showed that energy conservation improves with resolution. At $N=1.12$ billion particles, the solution closely matched reference 1D codes.
• Gresho-Chan Vortex: Tests showed that while high artificial viscosity suppresses mixing, increasing resolution can compensate for this, allowing the vortex to remain stable.
• Blob Test & Kelvin-Helmholtz Instability: These tests highlighted the limitations of SPH in fluid mixing. While high resolution allowed the cloud to fragment, the material did not fully "mix" with the wind due to SPH's entropy conservation properties (lack of artificial diffusion). However, the results were significantly better than traditional SPH implementations, particularly with the Wendland C6 kernel and 400 neighbors.

Performance and Scaling

• Weak Scaling: The code achieved >75% parallel efficiency when scaling up to 1024 nodes with 8 million particles per node.
• Strong Scaling: With 2 billion particles, the code achieved a speedup of 678x on 1024 nodes (66% efficiency). Even with smaller particle counts (1 million particles on 1024 nodes), a speedup of 7x was observed.
• Complexity: The time per step scales as $O(N^{4/3})$ due to the interplay between $O(N)$ operations and the $N^{1/3}$ scaling of the time step size required by the CFL condition.

5. Significance

This work represents a major advancement in computational astrophysics and planetary science:

• Unprecedented Resolution: The ability to run simulations with billions of particles (e.g., a 2.1 billion particle simulation of Jupiter's core disruption) allows researchers to resolve fine-scale features like planetary crusts, oceans, and tenuous atmospheres that were previously unresolvable.
• Planetary Impact Modeling: The code is specifically tailored for shock physics in planetary collisions. Its Lagrangian nature allows for the precise tracking of material origin and thermodynamic evolution, crucial for understanding moon formation, debris disks, and giant impacts.
• HPC Utilization: By fully leveraging modern heterogeneous HPC systems (CPU + GPU), pkdgrav3 enables the exploration of vast parameter spaces for planetary formation and evolution, moving beyond low-resolution models.
• Future Applications: The framework sets the stage for future extensions, including strength/failure models, radiative cooling, and double-precision arithmetic, further bridging the gap between numerical simulation and physical reality in shock physics.

In summary, pkdgrav3 provides a powerful, flexible, and scalable platform that combines the robustness of SPH for free-surface flows with the efficiency of modern HPC, enabling a new generation of high-resolution astrophysical simulations.