A Comparison of Massively Parallel Performance Portable Particle-in-Cell schemes for electrostatic kinetic plasma simulations

This paper evaluates the performance and portability of various Poisson solvers, including FFT, PCG, FEM, and the novel Particle-in-Fourier (PIF) schemes, within the IPPL library for electrostatic PIC simulations on diverse GPU architectures, finding that while FFT is fastest, the PIF scheme offers excellent scalability as a high-fidelity alternative.

The Gist

Imagine you are trying to simulate a massive crowd of people (particles) moving through a room, where everyone is pushing and pulling on each other based on their position. This is essentially what scientists do when they simulate plasma (a super-hot, electrically charged gas) to understand how it behaves.

This paper is a "race report" comparing different ways to calculate the forces between these particles to see which method is the fastest and most reliable on the world's most powerful supercomputers.

Here is the breakdown of the race using simple analogies:

The Setting: The "Particle-in-Cell" Loop

Think of the simulation as a game played in rounds. In every round, the computer does four things:

1. Scatter: It takes the positions of the particles and paints their "charge" onto a grid (like a chessboard).
2. Solve: It calculates the electric field (the push/pull force) everywhere on that grid based on the painted charges. This is the main event of the race.
3. Gather: It reads the force from the grid and tells each particle how to move.
4. Push: The particles move to their new spots.

The authors tested four different "Solvers" (methods to calculate step 2) to see which one wins.

The Four Racers

1. The FFT Solver (The Fast Sprinter)

• How it works: This method uses a mathematical trick called the "Fast Fourier Transform." Imagine trying to solve a puzzle by instantly seeing the whole picture in a mirror rather than looking at one piece at a time. It's incredibly fast.
• The Catch: It only works if the room has "periodic" boundaries. Think of it like a video game world where if you walk off the right edge, you instantly appear on the left. It can't handle walls or open doors.
• The Result: It was the absolute fastest in terms of raw time. However, on one specific supercomputer (Alps), it stumbled because the "particle movement" part of the loop got stuck, slowing down the whole race.

2. The PCG Solver (The Reliable Workhorse)

• How it works: This method breaks the grid down into tiny squares and solves the math step-by-step, like a detective checking every clue one by one. It uses a "Preconditioned Conjugate Gradient" approach.
• The Catch: It is much slower than the FFT (about 10 times slower in raw time), but it is very flexible. It can handle walls (Dirichlet) or open spaces (Neumann), not just the "wrap-around" video game world.
• The Result: It scales well (gets faster as you add more computers), but it takes longer to finish the job.

3. The FEM Solver (The High-Accuracy Architect)

• How it works: This is the "Finite Element Method." Instead of a rigid grid, it treats the space like a flexible mesh that can bend and fit complex shapes. It's like using a tailor-made suit instead of a boxy, off-the-rack shirt.
• The Catch: Like the PCG, it is slower than the FFT. It also struggles a bit with communication between computers because it has to constantly check the edges of its flexible mesh.
• The Result: It's great if you need high precision or complex shapes, but it's not the speed champion.

4. The PIF Solver (The New Contender)

• How it works: This is the "Particle-in-Fourier" scheme. Instead of painting particles onto a grid first, it projects them directly into "frequency space" (a mathematical representation of waves). It's like skipping the map entirely and navigating by the rhythm of the waves.
• The Catch: It requires special math (Non-Uniform FFTs) to handle particles that aren't perfectly aligned.
• The Result: It is more expensive (slower) than the FFT, but it is incredibly stable and accurate. It doesn't suffer from the "ghosting" or "aliasing" errors that happen when you try to fit a round particle onto a square grid. It scales beautifully on all machines, meaning it gets faster very efficiently as you add more power.

The Race Track (The Supercomputers)

The authors ran these tests on three different "tracks" (supercomputers) with different engines:

• Alps (Switzerland): Uses Nvidia's newest chips.
• LUMI (Finland): Uses AMD chips.
• JUWELS Booster (Germany): Uses older Nvidia chips.

The Winner's Podium

• Raw Speed: The FFT Solver wins hands down, but only if your problem fits its strict rules (periodic boundaries) and you aren't using the specific Alps machine where a technical glitch slowed it down.
• Flexibility: The PCG and FEM solvers are the best choice if your simulation has walls or complex shapes. They are slower but get the job done where the FFT can't go.
• High Fidelity: The PIF solver is the new star. While it takes a bit longer than the FFT, it offers the best balance of speed, stability, and accuracy. It's like a sports car that is slightly slower than a Formula 1 car but handles the turns much better and is safer to drive.

The Bottom Line

The paper concludes that there is no single "best" solver.

• If you need speed and have simple boundaries, use FFT.
• If you need flexibility (walls, complex shapes), use PCG or FEM.
• If you need high accuracy and stability without the errors of standard methods, PIF is an excellent, scalable alternative.

The authors also noted that they are currently working on fixing the "particle update" glitch on the Alps supercomputer and improving the "preconditioning" (a way to speed up the math) for the FEM solver to make them even faster in the future.

Technical Summary

Technical Summary: A Comparison of Massively Parallel Performance Portable Particle-in-Cell schemes for electrostatic kinetic plasma simulations

Problem Statement
The paper addresses the challenge of selecting optimal numerical schemes for electrostatic kinetic plasma simulations, specifically solving the Vlasov-Poisson system, within the context of modern heterogeneous High-Performance Computing (HPC) architectures. As supercomputing shifts toward heterogeneous systems combining CPUs and GPUs from multiple vendors (e.g., AMD, Intel, Nvidia), simulation codes must be both massively parallelized and performance-portable. The authors investigate which Particle-in-Cell (PIC) field solvers and alternative schemes offer the best balance of accuracy, scalability, and portability across different hardware platforms. The study focuses on the electrostatic Vlasov-Poisson system, a coupled nonlinear high-dimensional PDE system critical for beam physics and plasma simulations.

Methodology
The authors utilize the IPPL (Independent Parallel Particle Layer) library, a C++ framework designed for performance-portable, dimension-independent particle-mesh methods. IPPL leverages Kokkos for hardware abstraction, heFFTe for FFTs, and MPI for inter-node communication.

The study benchmarks four distinct schemes implemented within IPPL using the Landau damping mini-application as a test case. The schemes compared are:

1. FFT Solver: A pseudo-spectral solver using Fast Fourier Transforms, suitable for periodic and free-space boundaries.
2. PCG Solver: A matrix-free, second-order finite difference Preconditioned Conjugate Gradient solver, capable of handling Dirichlet, Neumann, and periodic boundaries.
3. FEM Solver: A matrix-free Finite Element Method solver using first-order Lagrange basis functions (providing second-order accuracy), capable of handling unstructured meshes and various boundary conditions.
4. PIF (Particle-in-Fourier) Scheme: A novel approach that avoids real-space grid interpolation. It projects particle charges directly to Fourier space using Non-Uniform FFTs (NUFFTs), solves the field equations spectrally, and interpolates forces back to particles. This method aims to eliminate aliasing and improve conservation properties.

Experimental Setup
The evaluation was conducted on three distinct supercomputing architectures:

• Alps (CSCS, Switzerland): 2688 nodes with GH200 (Nvidia Grace-Hopper) chips.
• LUMI-G (CSC, Finland): 2978 nodes with AMD MI250X GPUs.
• JUWELS Booster (Jülich, Germany): 936 nodes with Nvidia A100 GPUs.

Two problem sizes were tested for strong scaling:

• Case A: $512^3$ grid with 8 particles per cell (~1 billion particles).
• Case B: $1024^3$ grid with 8 particles per cell (~8.5 billion particles).

Weak scaling studies were also performed, increasing problem size proportionally with GPU count to maintain constant work per GPU.

Key Results

• Absolute Performance: The FFT solver consistently demonstrated the lowest absolute runtime across all architectures, benefiting from the optimizations in the heFFTe library. However, its applicability is restricted to periodic or free-space boundaries.
• Scalability and Bottlenecks:
  • On Alps (GH200), the FFT solver's strong scaling was limited by the particle-update kernel. Enabling CUDA IPC on these nodes increased the runtime of particle updates, causing the FFT phase (which was very fast) to no longer be the bottleneck; instead, particle communication became the dominant factor, degrading scaling efficiency below 50% after 32 GPUs for Case A.
  • The PCG and FEM solvers, where the linear solve or scatter/gather operations are the primary bottlenecks, masked the particle-update costs and exhibited better scaling on Alps.
  • The PIF scheme showed excellent scalability, maintaining efficiency above 50% up to 512 GPUs on LUMI and JUWELS Booster, and up to 1024 GPUs on Alps for the larger problem size.
• Solver Comparison:
  • PCG and FEM: These solvers were approximately an order of magnitude slower in absolute time than the FFT solver but scaled similarly in the electrostatic PIC context. The FEM solver incurred slightly higher costs due to conditional logic required for boundary conditions in matrix-vector evaluations and halo communication.
  • PIF: While costlier than the FFT-based PIC scheme, PIF offered high-fidelity results with excellent conservation and stability. On larger node counts, PIF sometimes outperformed the FFT solver in absolute time on specific architectures due to the FFT's particle-update bottleneck.
• Preconditioning: The study noted that matrix-free preconditioning for FEM remains a challenge. A naive implementation using conditional logic reduced CG iteration counts but increased runtime per iteration due to GPU sub-optimality with conditionals.

Significance and Claims
The authors claim that this work provides a comprehensive characterization of different PIC and PIF schemes on massively parallel, heterogeneous supercomputers. The primary contributions are:

1. Portability Demonstration: The study validates that the IPPL library can successfully deliver performance-portable particle-mesh simulations across diverse hardware (AMD and Nvidia GPUs) without requiring architecture-specific code rewrites.
2. Solver Selection Guidance: The paper clarifies the trade-offs between solvers. While FFT is superior in speed for periodic problems, PCG and FEM offer necessary flexibility for complex boundary conditions. PIF is identified as a high-fidelity alternative that avoids aliasing and offers competitive scalability, making it a viable option for simulations where conservation and stability are paramount.
3. Identification of Bottlenecks: The work highlights that in highly optimized FFT implementations on specific architectures (like GH200), the particle-update kernel can become the limiting factor, a nuance often overlooked when focusing solely on field solvers.

The paper concludes that the choice of solver should be tailored to the specific physical constraints (e.g., boundary conditions) and the target hardware, noting that higher-order schemes like PIF can fully exploit modern massively parallel architectures despite their computational cost. Future work is identified as improving matrix-free preconditioning for FEM to further enhance its performance.