Accelerating high-order continuum kinetic plasma simulations using multiple GPUs

TL;DR

It is demonstrated that the improved compute performance can be used to explore configurations which were previously computationally infeasible such as resolving fine-scale distribution function filamentation and multi-species dynamics with realistic electron-proton mass ratios.

Abstract

Kinetic plasma simulations solve the Vlasov-Poisson or Vlasov-Maxwell equations to evolve scalar-variable distribution functions in position-velocity phase space and vector-variable electromagnetic fields in physical space. The computational cost of evolving high-dimensional variables often limits the utility of continuum kinetic simulations and presents a challenge when it comes to accurately simulating real-world physical phenomena. To address this challenge, we developed techniques that accelerate and minimize the computational work required for a scalable Vlasov-Poisson solver. We present theoretical hardware compute and communication bounds required for solving a fourth-order finite-volume Vlasov-Poisson system. These bounds are then used to inform and evaluate the design of performance portable algorithms for a multiple graphics processing unit (GPU) accelerated version of the Vlasov-Poisson solver VCK-CPU. We demonstrate that the multi-GPU Vlasov solver implementation VCK-GPU simultaneously minimizes required inter-process data transfer while also being bounded by the machine network performance limits, leaving minimal room for theoretical performance improvements. This resulted in an overall strong scaling speedup per timestep of up to 40x in three-dimensional phase space (one position, two velocity coordinates) and 54x in four dimensional phase space (two position, two velocity coordinates) and a 341x increase in simulation throughput of the GPU accelerated code over the existing CPU code. The GPU code is also able to weak scale up to 256 compute nodes and 1024 GPUs. We then demonstrate how the improved compute performance can be used to explore configurations which were previously computationally infeasible such as resolving fine-scale distribution function filamentation and multi-species dynamics with realistic electron-proton mass ratios.

Figures (18)

• Figure 1: 1D-1V fourth-order finite-volume stencil required to update the solid black cell. The black cell and a subset of the dark gray cells are used for reconstruction in Eq. \ref{['eq:fvm_poly']}. Light gray cells correspond to $C_{\boldsymbol{i}}$, defined in Eq. \ref{['eq:spatial_corrections']}. Dashed outline cells are not used to update the solid black cell. Dark gray and black cells are not needed for $C_{\boldsymbol{i}}$ since their contributions cancel out.
• Figure 2: Local moment integration kernel designs for different memory layouts. Algorithm L1 is optimized for memory contiguous in $\boldsymbol{v}$ while Algorithm L2 is optimized for memory contiguous in $\boldsymbol{x}$. Both algorithms are generally equally efficient so long as the layout of $\bar{f}$ matches the algorithm used.
• Figure 3: Local zeroth moment integration kernel performance for an $N^{d+v}$ domain. Figure \ref{['fig:moment0_layout_test']} demonstrates that in 1D-1V both Algorithms L1 and L2 are equally efficient so long they are matched with the correct memory layout. Figure \ref{['fig:moment0_perf_port']} demonstrates that Algorithm L1 can be scaled to 1D-1V and 2D-2V across multiple GPU architectures without loss of efficiency. Error bars denote the 5th, 50th, and 95th percentiles and the 95% confidence interval is shaded. Student's t-test is used to determine the number of samples required for a 95% confidence intervalstudent1908. The kernel has little noise in execution time.
• Figure 4: Walltime per Poisson solve for domains with $N$ cells along each dimension. GPU-accelerated FFT solvers outperform sparse matrix solvers for small $N$, however, the better scaling performance of AMG-based methods leads eventually to AMG methods performing better. Error bars denote the 5th, 50th, and 95th percentiles and the 95% confidence interval is shaded. Student's t-test is used to determine the number of samples required for a 95% confidence intervalstudent1908.
• Figure 5: Effective memory bandwidth of the fused stage+fast RK4 method utilizing a nested tiling strategy for the spatial parallelization. Similar portable performance is achieved in 1D-2V while in 2D-2V the reduced caching performance of the MI250x negatively impacts overall throughputsai2022. Error bars denote the 5th, 50th, and 95th percentiles and the 95% confidence interval is shaded. Student's t-test is used to determine the number of samples required for a 95% confidence intervalstudent1908. The kernel has little noise in execution time.