🔭 astrophysics

A novel algorithm for GPU-accelerated particle-mesh interactions implemented in the QUOKKA code

This paper presents a novel, GPU-optimized "particle-mesh-particle" algorithm implemented in the QUOKKA code that efficiently handles particle-mesh interactions in massively parallel hydrodynamics simulations, achieving approximately 50% weak-scaling efficiency on up to 8192 GPUs while enabling scalable studies of star formation and feedback.

Original authors: Chong-Chong He (ANU), Benjamin D. Wibking (MSU), Aditi Vijayan (ANU), Mark R. Krumholz (ANU), Pak Shing Li (SAO)

Published 2026-03-19

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Chong-Chong He (ANU), Benjamin D. Wibking (MSU), Aditi Vijayan (ANU), Mark R. Krumholz (ANU), Pak Shing Li (SAO)

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to simulate the life of a galaxy on a computer. You have billions of tiny bits of gas (the "mesh") and billions of stars (the "particles"). The stars pull on the gas with gravity, and sometimes they explode as supernovas, blowing energy back into the gas.

The problem? Doing this math on a supercomputer with thousands of graphics cards (GPUs) is like trying to organize a massive party where everyone is shouting at once.

The Old Problem: The "Neighbour Search" Bottleneck

In the past, to figure out how a star affects the gas around it, the computer had to play a game of "Who's near me?" for every single star.

The Analogy: Imagine you are at a crowded concert. To find out who is standing next to you, you have to turn your head and ask every single person in the crowd, "Are you near me?"
The Issue: If you have millions of people (stars) and thousands of computers (GPUs) working together, this "asking" process is incredibly slow. It requires constant, messy communication between the computers. It's like trying to pass a note across a stadium while everyone is screaming; the signal gets lost, and the computers spend more time talking to each other than doing the actual work.

The New Solution: The "Particle-Mesh-Particle" Algorithm

The authors of this paper invented a new way to do this math that is much faster and works perfectly on modern supercomputers. They call it the "Particle-Mesh-Particle" algorithm.

Here is how it works, using a simple analogy:

1. The "Buffer" (The Whiteboard)

Instead of stars asking neighbors who is near them, they just write their notes on a giant whiteboard (the "mesh") that everyone can see.

Step 1: Every star looks at the whiteboard cells right next to it and writes down: "I am adding mass here," or "I am blowing energy here."
The Trick: If two stars try to write on the same spot on the whiteboard at the same time, the computer uses a special "atomic addition" trick (like a magic pen that instantly adds the numbers together without the writers bumping into each other).

2. The "Ghost Zones" (The Handshake)

The whiteboard is split up among thousands of computers. Some computers have the left side, some have the right.

The Trick: The computers share a few extra rows of the whiteboard (called "ghost zones") with their neighbors. This way, if a star on the left edge writes a note that spills over to the right edge, the neighbor computer sees it immediately. This is the only time the computers need to talk to each other, and they do it in a very organized, predictable way.

3. The "Limiter" (The Safety Valve)

Sometimes, two stars might try to suck up all the gas in a single cell, or two supernovas might try to add so much energy that the math breaks (like trying to make a density negative, which is physically impossible).

The Solution: The computer has a "safety valve" (a limiter). It checks the notes on the whiteboard. If the total change is too extreme, it scales it back just enough to keep the physics realistic, ensuring the simulation doesn't crash.

4. The Update (The Final Result)

Finally, the computer takes the notes from the whiteboard and updates the actual gas and the stars.

The gas gets the new mass and energy.
The stars get the new mass (if they ate gas) or disappear (if they exploded).

Why This is a Big Deal

1. No More "Shouting": The old method required stars to constantly search for neighbors. This new method is like everyone just writing on a shared board. It's much faster.
2. Perfect for Supercomputers: Modern supercomputers (like the Frontier mentioned in the paper) have thousands of GPUs. This algorithm is designed so that all these GPUs can work in perfect sync without getting confused.
3. Real-World Results: The authors tested this by simulating:

Gas falling into black holes (Bondi accretion).
Stars exploding (Supernovae).
Whole galaxies forming.

They found that their new method is 50% efficient even when running on 8,192 GPUs simultaneously. That's like having a choir of 8,000 singers all singing in perfect harmony without anyone being out of tune.

The Bottom Line

Before this paper, simulating how stars and gas interact on the world's fastest computers was like trying to herd cats. It was slow, messy, and prone to errors.

This new algorithm is like giving every cat a whistle and a clear map. Now, the stars and gas can interact efficiently, allowing scientists to run massive, realistic simulations of galaxy formation and star birth that were previously impossible. It's a game-changer for understanding how the universe works.

1. Problem Statement

Modern astrophysical simulations of star formation and galactic feedback increasingly rely on GPU-accelerated architectures to handle massive computational loads. However, integrating particles (representing sink particles, stars, or black holes) with grid-based hydrodynamics presents significant challenges on GPUs:

Neighbour Searches: Traditional particle-based methods (like SPH) require expensive neighbour searches, which involve irregular memory access patterns and complex synchronization that are inefficient on GPUs.
Order-Dependence: Many feedback recipes (e.g., supernova energy injection) are non-commutative. On CPUs, particles can be processed in a fixed serial order to ensure reproducibility. On massively parallel GPUs, especially across distributed memory nodes, enforcing a serial order is prohibitively expensive or impossible, leading to non-deterministic and inconsistent results.
Communication Overhead: Existing particle-mesh (PM) methods often struggle with complex interaction rules (like accretion limits or feedback merging) that require global communication or complex data structures not native to GPU optimization.

2. Methodology: The "Particle-Mesh-Particle" Algorithm

The authors introduce a novel Particle-Mesh-Particle (PMP) algorithm designed specifically for GPU-native execution. This method avoids direct particle-particle searches and instead uses a buffer-based deposition scheme.

Core Steps of the Algorithm:

Particle to Buffer Mesh (Deposition):
- A buffer mesh (with ghost zones) is allocated, matching the hydrodynamic grid dimensions.
- Particles deposit their interaction quantities (mass, momentum, energy changes) into this buffer.
- Multiple particles contributing to the same cell are accumulated using atomic additions. This avoids race conditions without requiring complex locking mechanisms.
Inter-GPU Summation (Communication):
- The buffer mesh is summed across distributed memory ranks. Ghost zones on one rank are added to the real zones of neighboring ranks.
- Key Optimization: This communication pattern is identical to the standard ghost-zone filling used in hydrodynamic updates, requiring only one communication step per cycle. It avoids expensive all-to-all particle exchanges.
Application of Limiters:
- A cell-wise limiter function is applied to the summed buffer to enforce physical constraints (e.g., ensuring density remains positive after accretion, or capping energy injection).
- This step resolves the "order-dependence" issue. By aggregating all contributions first and then applying a deterministic limiter, the result is independent of the order in which particles were processed.
Update Particles:
- Particle properties (mass, momentum) are updated based on the limited changes. This ensures conservation laws (e.g., total mass) are maintained even if accretion is capped.
Update Mesh:
- The limited changes are applied to the real hydrodynamic mesh state.

Implementation Details:

Sink Particles: Used for star formation. The algorithm implements Bondi-Hoyle accretion with a kernel radius ( $r_K$ ) and enforces the Truelove condition to prevent artificial fragmentation.
Supernova (SN) Feedback: Used for galaxy-scale simulations. The algorithm adapts the TIGRESS prescription but modifies it to be GPU-friendly by avoiding homogenization of the deposition zone and applying limiting to the collective action of multiple SNe rather than sequentially.
Stochastic Star Formation: A stochastic model generates stellar populations based on a Chabrier IMF, spawning high-mass stars (which will explode as SNe) and low-mass stellar particles.

3. Key Contributions

Novel Algorithm: The PMP scheme is the first to enable efficient, scalable, and deterministic particle-mesh coupling on massively parallel GPU architectures without particle neighbour searches.
Deterministic Results: By aggregating contributions before applying limits, the algorithm guarantees reproducible results regardless of particle processing order or MPI rank distribution.
GPU-Native Implementation: The method is implemented in Quokka, a radiation-magnetohydrodynamics code. It utilizes atomic operations and standard ghost-zone communication, minimizing memory transfers between CPU and GPU.
Extensibility: The framework is general and can be extended to other feedback mechanisms (stellar winds, protostellar outflows) and N-body calculations.

4. Results and Validation

The authors validated the algorithm through a comprehensive suite of tests:

Isothermal Sphere Collapse: Demonstrated excellent agreement with the analytic Shu solution for self-similar collapse, with relative errors in density and velocity < 2%.
Bondi & Bondi-Hoyle Accretion:
- Reproduced the exact analytic accretion rate in the unresolved regime ( $r_B \ll \Delta x$ ).
- Showed convergence in the resolved regime ( $r_B \gg \Delta x$ ) and acceptable accuracy in marginally resolved cases.
- Tested various kernel sizes ( $2\Delta x$ to $4\Delta x$ ) with minimal impact on accuracy.
Single Supernova Feedback:
- Compared against the TIGRESS prescription. The new method correctly recovers the empirical radial momentum vs. ambient density relation across different resolution regimes (Ejecta, Sedov-Taylor, and Momentum-Conserving).
- Avoided the "overcooling" issues seen in pure thermal injection models, thanks to an improved ODE solver for temperature evolution.
Multiple Supernova Feedback:
- Tested scenarios with two overlapping supernovae. The limiter successfully handled the superposition of feedback, maintaining convergence in terminal radial momentum within 20% of the high-resolution reference solution, even when deposition kernels overlapped.
Weak Scaling Performance:
- Benchmarked on the Frontier supercomputer (up to 8192 GPUs).
- Achieved ~50% weak-scaling efficiency from 2 to 1024 nodes.
- The particle update time remained less than 10% of the total hydrodynamic update time, demonstrating that the particle module does not bottleneck the simulation.

5. Significance

This work represents a paradigm shift in how particle-mesh interactions are handled in astrophysical simulations.

Scalability: It enables the first efficient, large-scale simulations of star formation and feedback on GPU-based supercomputers, overcoming the communication and synchronization bottlenecks that previously limited such studies.
Robustness: The algorithm solves the critical issue of non-deterministic results in parallel feedback injection, a major hurdle in previous GPU implementations.
Future Impact: By providing a robust, scalable framework, this method paves the way for high-resolution, full-physics simulations of galaxy evolution and the interstellar medium on next-generation exascale computing systems.

In summary, the paper presents a highly optimized, mathematically rigorous, and practically validated algorithm that successfully bridges the gap between complex particle physics and the architectural constraints of modern GPU supercomputing.