Hostname: page-component-74656f9884-kv2lp Total loading time: 0 Render date: 2026-11-06T00:20:32.524Z Has data issue: false hasContentIssue false
  • English
  • Français

GPU-MPI parallelisation for QuickPIC, a three-dimensional quasi-static particle-in-cell algorithm

Published online by Cambridge University Press:  28 July 2026

Yueran Tian ,
Yueluo Wang ,
Thamine Dalichaouch ,
Viktor Decyk ,
Frank S. Tsung ,
Warren B. Mori  and
Weiming An Open the ORCID record for Weiming An [Opens in a new window]
Yueran Tian
Affiliation:
School of Physics and Astronomy, Beijing Normal University, No. 19, Xinjiekouwai St., Haidian District, Beijing 100875, PR China
Yueluo Wang
Affiliation:
School of Physics and Astronomy, Beijing Normal University, No. 19, Xinjiekouwai St., Haidian District, Beijing 100875, PR China Department of Astronomy, Peking University, No. 5, Yiheyuan Road, Haidian District, Beijing 100871, PR China
Thamine Dalichaouch
Affiliation:
Department of Physics and Astronomy, University of California Los Angeles, Los Angeles, CA 90095-1547, USA
Viktor Decyk
Affiliation:
Department of Physics and Astronomy, University of California Los Angeles, Los Angeles, CA 90095-1547, USA
Frank S. Tsung
Affiliation:
Department of Physics and Astronomy, University of California Los Angeles, Los Angeles, CA 90095-1547, USA
Warren B. Mori
Affiliation:
Department of Physics and Astronomy, University of California Los Angeles, Los Angeles, CA 90095-1547, USA Department of Electrical and Computer Engineering, University of California Los Angeles, Los Angeles, CA 90095-1547, USA
Weiming An*
Affiliation:
School of Physics and Astronomy, Beijing Normal University, No. 19, Xinjiekouwai St., Haidian District, Beijing 100875, PR China Institute for Frontiers in Astronomy and Astrophysics, Beijing Normal University, Changping District, Beijing 102206, PR China
*
Corresponding author: Weiming An, anweiming@bnu.edu.cn
Rights & Permissions [Opens in a new window]

Abstract

QuickPIC is a quasi-static PIC program for simulating plasma wakefield acceleration. It supports a hybrid parallel simulation on multiple central processing units via OpenMP and message passing interface (MPI). We now present a version of QuickPIC that runs on graphics processing units (GPU), with MPI for communication between different GPUs, based on the UCLA Parallel Partice-In-Cell framework. By rewriting relevant Fortran77 functions with CUDA C, we port the two-dimensional particle and field computations onto GPU platforms. QuickPIC-GPU supports simulation with a single GPU for up to a resolution of $1024^3$ and 25 particles per cell to avoid time spent on cross-processor communications, and simulation with multiple GPUs for larger problem sizes.

Keywords

plasma simulationintense particle beams

Information

Type
Research Article
Information
Journal of Plasma Physics , Volume 92 , Issue 4 , August 2026 , E103
Check for updates
Creative Commons
Creative Common License - CCCreative Common License - BYCreative Common License - NC
This is an Open Access article, distributed under the terms of the Creative Commons Attribution-NonCommercial licence (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original article is properly cited. The written permission of Cambridge University Press or the rights holder(s) must be obtained prior to any commercial use.
Copyright
© The Author(s), 2026. Published by Cambridge University Press

1. Introduction

Plasma/laser wakefield acceleration uses high-energy particle/laser beams to generate oscillations in plasma, which in turn creates longitudinal electromagnetic fields that accelerate charged particles (Tajima & Dawson Reference Tajima and Dawson1979; Chen et al. Reference Chen, Dawson, Huff and Katsouleas1985). If the magnitude of the group velocity of the high-energy beam is close to the speed of light, then the phase velocity of the plasma waves would also be close to the speed of light, therefore maintaining acceleration over a relatively long distance (Joshi & Katsouleas Reference Joshi and Katsouleas2003). It has been proved to be a highly feasible and efficient alternative from conventional accelerators, as shown in various experiments (Blumenfield et al. Reference Blumenfield2007; Litos et al. Reference Litos2014; Corde et al. Reference Corde2015; Adli et al. Reference Adli2018).

To predict the effects of plasma acceleration experiments and improve the design of plasma wakefield accelerators, computer simulations are essential. The particle-in-cell (PIC) algorithm (Dawson Reference Dawson1983; Birdsall & Langdon Reference Birdsall and Langdon2018) is a mainstream method for simulating plasma-based acceleration. It discretises electromagnetic fields on grid points and self-consistently tracks the motion of particles between grids, and is known for its robustness and efficiency (Joshi, Corde & Mori Reference Joshi, Corde and Mori2020). In experiments, the particle density in plasmas far exceeds the number of particles that a computer can track (Dawson Reference Dawson1983). Thus, each ‘particle’ in a PIC simulation represents many real particles in the plasma.

Generally, the PIC algorithm involves four steps:

  1. (i) Initialise the position and momenta of the particles within the computational domain together with the initial electromagnetic fields.

  2. (ii) Calculate the charge and current densities at grid points (or midpoints) according to the particle positions and momenta.

  3. (iii) Solve Maxwell’s equations to update electromagnetic fields (usually using spectral (Huang et al. Reference Huang, Decyk, Ren, Zhou, Lu, Mori, Cooley, Antonsen and Katsouleas2006; Decyk Reference Decyk2007), finite-element (Bai et al. Reference Bai, Cao, Chu and Zhang2018) or finite-difference (Fonseca et al. Reference Fonseca2002) methods).

  4. (iv) Interpolate electromagnetic fields at particle locations and update particle positions and momenta over a time step.

Steps 2–4 are repeated until the program reaches the desired simulation duration.

Popular PIC codes include OSIRIS (Fonseca et al. Reference Fonseca2002), WarpX (Myers et al. Reference Myers2021), Smilei (Derouillat et al. Reference Derouillat2018), EPOCH (Arber et al. Reference Arber2015), etc. However, these ‘fully explicit’ (traditional) (Joshi et al. Reference Joshi, Corde and Mori2020) PIC codes are considered to have too high a computational load due to the full electromagnetic solver (Huang et al. Reference Huang, Decyk, Ren, Zhou, Lu, Mori, Cooley, Antonsen and Katsouleas2006). Therefore, codes such as WAKE (Mora & Antonsen Reference Mora and Antonsen1997), QuickPIC (An et al. Reference An, Decyk, Mori and Antonsen2013), QPAD (Li et al. Reference Li, An, Decyk, Xu, Hogan and Mori2021), HiPACE++ (Diederichs et al. Reference Diederichs, Benedetti, Huebl, Lehe, Myers, Sinn, Vay, Zhang and Thévenet2022), WAND-PIC (Wang et al. Reference Wang, Khudik, Kim and Shvets2022) and LCODE (Sosedkin & Lotov Reference Sosedkin and Lotov2016) use a method called the quasi-static approximation: since particle beams evolve at a much longer time scale than the plasma wavelength, these codes separate the beam and plasma evolution time scales, and ‘pause’ the beam evolution while calculating multiple steps of plasma response (Huang et al. Reference Huang, Decyk, Ren, Zhou, Lu, Mori, Cooley, Antonsen and Katsouleas2006). This approximation brings a considerable speed up to PIC algorithms. QuickPIC, for example, can be 10 squared $10^2$ to 10 Superscript 4 $10^4$ times faster than a full PIC code (e.g. OSIRIS) (An et al. Reference An, Decyk, Mori and Antonsen2013), while maintaining accuracy for many problems.

Also, advances in computer hardware can bring further speedup. To adapt to evolving computer architectures, PIC codes have been parallelised on both central processing unit (CPU) (Liewer & Decyk Reference Liewer and Decyk1989) and graphics processing unit (GPU) (Decyk & Singh Reference Decyk and Singh2011) platforms. Parallelised PIC programs typically divide the simulation domain into subregions assigned to different processors. Each processor handles particles and fields within its subregion, and particle and field data are transferred across boundaries when needed. These programs show good scaling when run on multiple processors, and can be further accelerated by moving from CPU to GPU platforms.

Currently, many GPU-based fully explicit PIC codes (PIConGPU (Burau et al. Reference Burau, Widera, Hönig, Juckeland, Debus, Kluge, Schramm, Cowan, Sauerbrey and Bussmann2010); Picador (Bastrakov et al. Reference Bastrakov, Donchenko, Gonoskov, Efimenko, Malyshev, Meyerov and Surmin2012); WarpX (Myers et al. Reference Myers2021), etc.) have been developed, and these codes all support the use of multiple GPUs through message passing interface (MPI) communication. Also, the GPU versions of WarpX and Picador both show considerable speedup compared with their CPU versions. However, to our knowledge, the only quasi-static PIC code with GPU acceleration, before our work, is HiPACE++ (Diederichs et al. Reference Diederichs, Benedetti, Huebl, Lehe, Myers, Sinn, Vay, Zhang and Thévenet2022). It achieves parallelisation via the AMReX framework, in which low-level functions are compiled into either CPU vectorised loops or GPU kernels, depending on the target platform. Also, it only supports one GPU for a two-dimensional (2-D) slice, though multiple GPUs can be used for different longitudinal domains via pipelining. Different from HiPACE++, we used CUDA C to directly write GPU kernel/host functions, for finer control of low-level details; also, we maintained the MPI parallelisation of the CPU program (with each slice assigned to one or more MPI rank(s), and assigned a GPU to each MPI rank, so that multiple GPUs can be used for one 2-D slice (QuickPIC also supports longitudinal parallelisation via pipelining).

QuickPIC is based on the UPIC (UCLA parallel PIC) framework (Decyk Reference Decyk2007). UPIC is an open-source framework developed by V.K. Decyk et al. at University of California, Los Angeles (UCLA) to provide optimised components for large-scale parallel PIC codes. It offers common algorithms (e.g. particle momentum/position updates, charge density calculations, field solvers and MPI-based data transfers) and templates for building customised PIC codes. The framework supports multiple programming paradigms, with Fortran 77 for low-level functions, Fortran 90 for encapsulation and C/Fortran 2003 for high-level logic. UPIC also supports MPI-based distributed memory parallelisation, dividing the simulation domain into subregions with an equal number of particles (or volume) assigned to different ranks; shared-memory parallelisation inside a rank is achieved via OpenMP for CPU and CUDA for GPU platforms.

To adapt QuickPIC for GPU-MPI parallelisation, we referenced the PIC-skeleton code within the UPIC framework (Decyk Reference Decyk2015; Decyk Reference Decyk2023). The PIC-skeleton provides a simplified, modular structure for implementing GPU-accelerated PIC algorithms, focusing on core computational components while omitting complex input/output and analysis features. This framework is designed for educational and debugging purposes, allowing developers to isolate and optimise critical functions such as particle pushing, charge deposition and field solving.

In this article, we first review the basic equations in QuickPIC and the structure of the program. Then, we present the revisions of the 2-D loop in QuickPIC. These revisions realise the complete GPU-MPI hybrid parallelism of the 2-D loop which processes the plasma wakefield, and the storage of all wakefield data on GPUs. After revision methods, we report the correctness and speed of the new program under different conditions. Finally, we summarise the performance of the new code and discuss potential directions for future development.

2. QuickPIC algorithm review

In solving for particle movement and field evolution, QuickPIC adopts a coordinate system left parenthesis x comma y comma s equals z comma xi equals c t minus z right parenthesis $(x,y,s=z,\xi =ct-z)$ , where s $s$ represents the penetration depth of the high-energy particle beam, and xi $\xi$ denotes the distance relative to the beam head; x $x$ and y $y$ are 2-D orthogonal coordinates in directions perpendicular to beam propagation z $z$ ; c $c$ denotes the speed of light and t denotes time. Particle movements in plasma are described by x $x$ , y $y$ and xi $\xi$ , while movements of the beam particles are described by x $x$ , y $y$ , xi $\xi$ and s $s$ . Under the quasi-static assumption, we consider beam evolution along s $s$ to be far slower than the plasma response along xi $\xi$ . Therefore, as shown in figure 1, QuickPIC separates beam and plasma computations into 3-D (outer) and 2-D (inner) loops. In each 3-D step, the simulation box is sliced along xi $\xi$ , and a slice of unperturbed plasma is initialised in front of the beam. Then, the slice travels ‘backwards’ in xi $\xi$ , while the evolution in transverse directions is calculated in 2-D steps. In each 2-D step, particles interact via the electromagnetic fields they produce. First, charge and current densities of the particles are deposited onto nearby grid points; then, the fields at grid points are evaluated using a simplified version of Maxwell’s equations, neglecting terms involving partial differential slash partial differential s $\partial / \partial s$ and partial differential squared slash partial differential s squared $\partial ^{2}/ \partial s^{2}$

(2.1) bold nabla Subscript up tack Baseline times bold italic upper E Subscript up tack Baseline equals minus StartFraction partial differential Over partial differential xi EndFraction upper B Subscript z Baseline ModifyingAbove z With caret comma \begin{equation} \boldsymbol{\nabla} _{\perp } \times \boldsymbol{E}_{\perp } = -\frac {\partial }{\partial \xi } B_z \hat {z} ,\end{equation}
(2.2) bold nabla Subscript up tack Baseline times left parenthesis upper E Subscript z Baseline ModifyingAbove z With caret right parenthesis equals minus StartFraction partial differential Over partial differential xi EndFraction left parenthesis bold italic upper B Subscript up tack Baseline minus ModifyingAbove z With caret times bold italic upper E Subscript up tack Baseline right parenthesis comma \begin{equation} \boldsymbol{\nabla} _{\perp } \times ( E_z \hat {z} ) = -\frac {\partial }{\partial \xi } ( \boldsymbol{B}_{\perp } - \hat {z} \times \boldsymbol{E}_{\perp } ) ,\end{equation}
(2.3) bold nabla Subscript up tack Baseline times bold italic upper B Subscript up tack Baseline minus upper J Subscript z Baseline ModifyingAbove z With caret equals StartFraction partial differential Over partial differential xi EndFraction upper E Subscript z Baseline ModifyingAbove z With caret comma \begin{equation} \boldsymbol{\nabla} _{\perp } \times \boldsymbol{B}_{\perp } - J_z \hat {z} = \frac {\partial }{\partial \xi } E_z \hat {z} ,\end{equation}
(2.4) bold nabla Subscript up tack Baseline times left parenthesis upper B Subscript z Baseline ModifyingAbove z With caret right parenthesis minus bold italic upper J Subscript up tack Baseline equals StartFraction partial differential Over partial differential xi EndFraction left parenthesis bold italic upper E Subscript up tack Baseline plus ModifyingAbove z With caret times bold italic upper B Subscript up tack Baseline right parenthesis comma \begin{equation} \boldsymbol{\nabla} _{\perp } \times ( B_z \hat {z} ) - \boldsymbol{J}_{\perp } = \frac {\partial }{\partial \xi } ( \boldsymbol{E}_{\perp } + \hat {z} \times \boldsymbol{B}_{\perp } ) ,\end{equation}
(2.5) bold nabla Subscript up tack Baseline bold dot bold italic upper E Subscript up tack Baseline minus rho equals StartFraction partial differential Over partial differential xi EndFraction upper E Subscript z Baseline comma \begin{equation} \boldsymbol{\nabla} _{\perp } \boldsymbol{\cdot }\boldsymbol{E}_{\perp } - \rho = \frac {\partial }{\partial \xi } E_z ,\end{equation}
(2.6) bold nabla Subscript up tack Baseline bold dot bold italic upper B Subscript up tack Baseline equals StartFraction partial differential Over partial differential xi EndFraction upper B Subscript z Baseline comma \begin{equation} \boldsymbol{\nabla} _{\perp } \boldsymbol{\cdot }\boldsymbol{B}_{\perp } = \frac {\partial }{\partial \xi } B_z ,\end{equation}
(2.7) bold nabla Subscript up tack Superscript 2 Baseline left parenthesis phi minus upper A Subscript z Baseline right parenthesis equals minus left parenthesis rho minus upper J Subscript z Baseline right parenthesis comma \begin{equation} \boldsymbol{\nabla} _{\perp }^2 \left ( \phi - A_z \right ) = - \left ( \rho - J_z \right ) ,\end{equation}

in which bold nabla Subscript up tack Baseline equals left parenthesis partial differential slash partial differential x right parenthesis ModifyingAbove x With caret plus left parenthesis partial differential slash partial differential y right parenthesis ModifyingAbove y With caret $\boldsymbol{\nabla} _{\perp }=({\partial }/{\partial x})\hat {x}+({\partial }/{\partial y})\hat {y}$ , rho $\rho$ represents charge density, phi $\phi$ represents electric potential, upper A $A$ represents magnetic vector potential, upper J $J$ represents electric current density, and bold italic upper E Subscript up tack $\boldsymbol{E}_{\perp }$ and bold italic upper B Subscript up tack $\boldsymbol{B}_{\perp }$ represent the electric field intensity and magnetic induction intensity in the x upper O y $xOy$ plane, respectively. These functions are solved using the spectral method; fields go through fast Fourier transform (FFT) (with sine and cosine as the orthogonal bases) into functions of wavenumbers left parenthesis k Subscript x Baseline comma k Subscript y Baseline right parenthesis $(k_x,k_y)$ , and are solved as linear equations, then transformed back to functions of left parenthesis x comma y right parenthesis $(x,y)$ . At last, the solved fields are used to update the positions and momenta of the particles. After each 2-D step, the state of the slice represents the state of the plasma particles in the corresponding xi $\xi$ position after the whole 3-D step. The complete 3-D–2-D loop is shown in figure 1. For parallel computation, each 2-D slice is assigned to one or more MPI ranks. If multiple MPI ranks are used, the slice is separated equally (in area) among them. Parallelisation along the xi $\xi$ direction is also supported via pipelining: the simulation area can be divided along xi $\xi$ into several pipeline stages, each containing an equal number of 2-D slices. One or more MPI ranks can be assigned to a pipeline stage, in which 2-D slices are processed sequentially, and each 2-D slice is dealt with in the aforementioned manner.

Figure 1.

Flow chart of 3-D–2-D loops in QuickPIC.

A flowchart illustrating the 3-D2-D loops in QuickPIC for plasma/laser wakefield acceleration.
Figure 1. Long description

The flowchart illustrates the 3-D2-D loops in QuickPIC for plasma/laser wakefield acceleration. The process begins with the initialization of the beam and plasma. Beam particles are then deposited, followed by a call to the 2D routine. The 2D routine involves depositing plasma particles, solving for fields, and pushing plasma particles. After the 2D routine ends, the process continues with pushing beam particles and renewing the plasma. This sequence repeats until the 3D routine ends. The moving window concept is depicted with a shift by a distance delta s, and the 2D routine is applied within a slice delta xi.

3. The GPU-MPI parallelisation approach for QuickPIC

3.1. Target analysis

QuickPIC employs a three-layer structure:

  1. (i) Bottom layer: Fortran77 functions for basic computations. These functions are based on legacy PIC codes at UCLA, and Fortran77 is chosen for its high performance; it is observed that Fortran77 performs better than any other language (at the time when UPIC and QuickPIC were developed) for it is rather primitive, and compilers can transform and optimise it well. Also, Fortran77 can be called by many other languages. However, it is inflexible and unsafe, and thus these functions need wrappers for more complex operations.

  2. (ii) Middle layer: Fortran 2003 classes encapsulating low-level details. These classes include arrays that store particle/field information, and define wrapper functions that call one or more Fortran77 functions, but offer simpler arguments. Also, these wrapper functions can check array sizes (which are hidden at higher layers), check if arrays have been allocated and make use of pointer arrays (which, if used at the Fortran77 layer, would degrade performance).

  3. (iii) Top layer: high-level classes combining middle-layer components for convenient calling in the main loop. Both the middle layer and the top layer are object oriented. The procedures in these two layers do not perform any substantial calculations; they merely make decisions and delegate the work to underlying Fortran77 functions (Decyk Reference Decyk2007).

Among all high-level classes, the field2d (electromagnetic fields) and species2d (particles) classes, which are called in the 2-D loop, dominate computational load. The field2d class relies on three mid-layer classes: fft2d (FFTs), fpois2d (Poisson equation solver) and ufield2d (field storage and arithmetic operations). The species2d class depends on the part2d class for particle-related computations, such as charge density calculation, momentum updates and particle migration across MPI boundaries. In addition, the part2d class relies on the fdist2d class for particle array initialisation. The mid-layer classes fft2d, fpois2d and part2d also use the ufield2d class, for they all include read/write operations on 2-D electromagnetic field data. Figure 2 shows the various main classes in QuickPIC, their relationships, and the main GPU revision target classes-fft2d, fpois2d, ufield2d, part2d and fdist2d, in which underlying Fortran 77 functions were replaced with CUDA C functions.

Figure 2.

Class diagram of QuickPIC. Classes outlined in red are the main GPU revision targets.

Flowchart illustrating the structure of QuickPIC with main GPU revision targets highlighted.

QuickPIC discretises the simulation domain into grids to compute and store field-related information. Grid cells are distributed evenly across non-shared-memory MPI ranks. Particles are inherently discrete, but are grouped by larger subregions called tiles – each containing multiple grid cells. Particles in the same tile are stored next to each other in a particle array, for faster r/w in particle-related operations (Decyk & Singh Reference Decyk and Singh2011).

For a domain with n Subscript x Baseline times n Subscript y $n_x\times n_y$ fine grids, and a $a$ MPI ranks on one or more physical computing nodes (each rank assigned with b $b$ CPU cores) assigned to the task, n Subscript x Baseline times n Subscript y Baseline divided by a $n_x\times n_y/a$ grids are assigned to each rank. The ufield2d class manages electromagnetic field data (e.g. bold italic upper E $\boldsymbol{E}$ and bold italic upper B $\boldsymbol{B}$ ) at these grid points. Above this, every m Subscript x Baseline times m Subscript y $m_x\times m_y$ grid cells form a tile. The total number of tiles in a rank is m Subscript x Baseline 1 Baseline times m Subscript y Baseline 1 $m_{x1}\times m_{y1}$ , where m Subscript x Baseline 1 Baseline equals left floor left parenthesis n Subscript x Baseline minus 1 right parenthesis divided by m Subscript x Baseline right floor plus 1 $m_{x1}=\lfloor (n_x -1)/m_x\rfloor +1$ and m Subscript y Baseline 1 Baseline equals left floor left parenthesis n Subscript y Baseline divided by a minus 1 right parenthesis divided by m Subscript y Baseline right floor plus 1 $m_{y1}=\lfloor (n_y/a -1)/m_y\rfloor +1$ ( left floor x right floor $\lfloor x \rfloor$ denotes the largest integer that is smaller than or equal to x $x$ ). A particle located near the grid cell left parenthesis n comma m right parenthesis $(n,m)$ belongs to the tile at row left floor n divided by m Subscript x Baseline right floor $\lfloor n/m_x \rfloor$ and column left floor m divided by m Subscript y Baseline right floor $\lfloor m/m_y\rfloor$ . The part2d class managess particle data (e.g. positions and momenta). Both field and particle data are stored in the shared memory of the rank. Each core within a rank processes m Subscript x Baseline 1 Baseline times m Subscript y Baseline 1 Baseline divided by left parenthesis a bold dot b right parenthesis $m_{x1}\times m_{y1}/(a\boldsymbol{\cdot }b)$ tiles; tasks are assigned to these cores via OpenMP. While the number of grid cells per tile is fixed, the particle count per tile varies dynamically during simulation. Thus, the part2d class includes:

  1. (i) part: a 2-D array (particle count times $\times$ attributes per particle) for storing particle positions and momenta.

  2. (ii) ppart: a 3-D array (number of tiles times $\times$ maximum particles per tile times $\times$ attributes per particle) for tracking particles within each tile.

  3. (iii) kpic: an array recording the number of particles in each tile.

These arrays (and other particle-related arrays) are allocated assuming a maximum number (larger than the number of particles initialised) of particles for each MPI rank, to allow for particles moving between MPI ranks, and a rank dealing with more particles than the number at initialisation.

The loops in the underlying Fortran 77 functions (e.g. particle updates, charge density calculations) were replaced with GPU-parallel kernel functions.

3.2. Mapping simulation objects to CUDA architecture

To achieve efficient GPU parallelisation, we established a hierarchical correspondence between QuickPIC’s simulation objects and CUDA’s computational units (threads, blocks, grids). Each MPI rank drives one GPU, which is the common practice of PIC codes with GPU support (Burau et al. Reference Burau, Widera, Hönig, Juckeland, Debus, Kluge, Schramm, Cowan, Sauerbrey and Bussmann2010; Bastrakov et al. Reference Bastrakov, Donchenko, Gonoskov, Efimenko, Malyshev, Meyerov and Surmin2012; Myers et al. Reference Myers2021; Xiong et al. Reference Xiong, Huang, Yuan, Kuang, Sharma, Jiang, Wu and Yu2025). Tiles are assigned to CUDA blocks and particles/grid points are distributed across threads. Take the function gpu_deposit, which calculates charge density at grid points, as an example

none

In this way, a block with n threads deals with a tile, and the i $i$ th thread in the block deals with the i $i$ th, left parenthesis n plus i right parenthesis $(n+i)$ th, left parenthesis 2 n plus i right parenthesis $(2n+i)$ th… particles in the tile (Decyk & Singh Reference Decyk and Singh2011), as shown in figure 3. This method has also been applied to other PIC codes (Xiong et al. Reference Xiong, Huang, Yuan, Sharma, Kuang, Jiang and Yu2024). For current deposition (calculation of current density at grid points), the function gpu_amjdeposit has a similar structure, only that the speed of the particle is also read from ppart[].

Figure 3.

Mapping of CUDA architecture to particle simulation.

A diagram illustrating the mapping of CUDA architecture to particle simulation, showing the structure of a simulation box divided into blocks and tiles, and the assignment of particles to threads within a block.
Figure 3. Long description

A diagram of the mapping of CUDA architecture to particle simulation. Panel A: The simulation box is divided into multiple blocks, each containing several tiles. Each block is labeled sequentially starting from block 0. Panel B: A detailed view of a single block corresponding to a tile. The block contains multiple particles, each assigned to a specific thread. The threads are labeled sequentially starting from thread 0. The particles are listed in order, with each particle assigned to a thread in a cyclic manner.

Field interpolation in QuickPIC-GPU is done in the function for pushing particles, gpu_ppush2l. This function first loops over grid points in a tile, to load global fields onto block-level shared memory; then, after a block-level synchronisation, loops over particles in a tile, to calculate field values at particle locations, and update particle positions and momenta

none

Also, since particles in the same tile are stored next to each other, when they move from one tile to another, they need to be reordered. The pmove() subroutine, which handles this, will be discussed later in the article.

For field computations without tile-level structure, the block number is set to be the maximum number of tasks that can be processed in parallel (for FFT processes, this is usually the grid size in the x or y direction; for other field-related processes, this is usually the total 2-D grid size) divided by the total number of threads in a block (currently set to 64), to best utilise GPU computational power, and minimise the task size of a single thread.

In both particle and field computations, CUDA block- and grid- level synchronisations are needed. For example, block-level synchronisation is needed before block-level shared-memory results in particle deposition being aggregated to global memory; grid-level synchronisation is needed when all blocks need to finish aggregating block-level particle deposition results before moving on to other tasks that would require these results. For block-level synchronisation we use the CUDA call __syncthreads(); grid-level synchronisation is implicitly done when the next task is another kernel call (it would wait in a queue of kernel calls before the GPU is finished with all the calls before it) and a memory copy between the CPU and GPU on the same rank (it would wait until all kernel calls before it are finished on the GPU, and block any kernel calls after it until it is done), and explicitly done by calling cudaDeviceSynchronize() before and after MPI transfers.

3.3. Atomic operations

When porting QuickPIC to GPUs, managing memory access order across threads became critical. Parallelising serial loops into kernel functions risks concurrent thread access to shared-memory addresses, leading to race conditions. For example, in the gpu_qdeposit() kernel, each thread calculates the charge density contribution of a particle to grid points. If multiple particles contribute to the same grid point, threads may simultaneously read, modify and write to the same memory location. Without synchronisation, threads could read unmodified values, resulting in incomplete updates (e.g. only one particle’s contribution recorded instead of the cumulative sum). Problems alike exist in current deposition, field interpolation, etc. To resolve this, CUDA-specific atomic operations (e.g. atomicAdd(), atomicMax()) were employed. These operations serialise access to shared/global memory locations, ensuring that only one thread modifies a target address at a time.

3.4. Host functions

Host functions in CUDA bridge the gap between high-level Fortran 2003 classes and low-level GPU kernels. These functions manage kernel launches, parameter passing and error handling. For example, the fgpu_deposit() host function encapsulates the gpu_deposit kernel

none

in which dimGrid denotes the number of blocks in a CUDA device (GPU), and dimBlock denotes the number of threads in each block.

Besides, memory management host functions handle GPU memory allocation, deallocation and data transfer.

3.5. Adaptations involving multiple kernel functions and MPI

3.5.1. Rewriting of pmove()

The pmove subroutine handles particle sorting: it sorts the particle array ppart[] once every 2-D step, to ensure that particles in the same tile are stored together, and can be quickly loaded for other particle operations, particle migration across tiles and MPI ranks. The ppart[] array is arranged as follows: the simulation area is separated into mxyp1 tiles, for each tile a maximum of nppmx particles are allowed, and the number of attributes of each particle is dim; therefore, the total size of ppart is dim*nppmx*mxyp1; the attribute l of particle m in tile n can be accessed via ppart[l+m*ndim+n*nppmx*dim], and the actual number of particles in tile n is stored in the nth place in the array kpic[]. Pmove mainly consists of three phases: firstly, the Fortran 77 function ppporder2la loops through all the particles in each tile. Each particle is determined whether it is leaving its current tile, and, if it does, in which direction it is going; the leaving particles and their destinations are stored in the array ppbuff[]; and the vacancies they left (positions in the particle array) is recorded in the ihole[] array. Furthermore, if a particle is determined to be crossing the MPI rank boundary, it is copied from ppbuff[] to another special buffer array. Secondly, particles crossing MPI rank boundaries are transferred via MPI send/receive operations. Thirdly, the ppporder2lb function loops through ihole[], and fills vacancies in each tile with incoming particles read from ppbuff; if incoming particles for a tile are fewer than the holes that departed particles left, the extra holes are filled with particles stored at the end of this tile; if there are more incoming particles than holes in a tile, the extra particles are attached to the end of the storage space allocated for this tile. At the same time, kpic[] is updated (Decyk & Singh Reference Decyk and Singh2011).

In the original code, ppporder2la and ppporder2lb contain multiple OpenMP loops with inter-dependent results. To adapt this for GPUs, each OpenMP loop was converted into a separate CUDA kernel. When these kernels are called by the host, they queue up to be executed sequentially (the next kernel is only executed after all the threads in the GPU have finished executing the last kernel), to ensure that inter-dependent results (e.g. number and positions of leaving particles in a tile) are calculated before they are read by a different thread in a different block. Buffer data are transferred directly between GPUs via CUDA-aware MPI.

3.5.2. Rewriting of FFT transpose

In QuickPIC, for multiple MPI ranks, the 2-D fields are partitioned in the y $y$ direction, and continuously stored in x $x$ , by default. For example, for 2 MPI ranks and a 2-D slice with boundaries left parenthesis minus x Subscript b Baseline comma x Subscript b Baseline right parenthesis $(-x_b,x_b)$ and left parenthesis minus y Subscript b Baseline comma y Subscript b Baseline right parenthesis $(-y_b,y_b)$ , 2-D data in left brace left parenthesis x comma y right parenthesis vertical bar minus x Subscript b Baseline less than x less than x Subscript b Baseline comma minus y Subscript b Baseline less than y less than 0 right brace $\{(x,y)|-x_b\lt x\lt x_b,-y_b\lt y\lt 0\}$ are stored in one rank, and 2-D data in left brace left parenthesis x comma y right parenthesis vertical bar minus x Subscript b Baseline less than x less than x Subscript b Baseline comma 0 less than y less than y Subscript b Baseline right brace $\{(x,y)|-x_b\lt x\lt x_b,0\lt y\lt y_b\}$ are stored in the other rank. Therefore, 2-D FFT in QuickPIC involves three steps: firstly, 1-D transform is performed in x $x$ : for the bit-reversal operation, each data point is assigned to a thread; for the iterative combination part, each row in the direction to be transformed is assigned to a block and each combination in an iteration step is assigned to a thread in the block; secondly, the whole field matrix is transposed across the simulation box, so that the field is now continuously stored in y $y$ ; finally, 1-D transform is performed in y $y$ ; the parallelisation method is the same as in the first step (Decyk & Singh Reference Decyk and Singh2011). In order to perform a matrix transpose across multiple ranks, 2-D data needs to be first transposed locally. Some data would be in a different rank when flipped across the main diagonal of the matrix; these are then sent to the corresponding rank(s). Other data would remain in the same rank, and need not be transferred. Similar to pmove, we rewrote the MPI communication functions into the separate C file, and linked them to the CUDA host functions for FFT transpose.

4. Performance evaluation

To evaluate accuracy and speed of the revised program, we tested the program before/after revision using the same input parameters (except for particle density and resolution). All simulations in this section ran on CPU/GPU clusters provided by the Beijing Super Cloud Computing Centre. For CPU clusters, each computing node is equipped with 2 AMD EPYC 7452 processors with 32 cores each; for GPU clusters, at most 8 NVIDIA A100 GPUs (PCIe, 40 GB, 2 computing nodes with 4 GPUs each) can be allocated for a task. The simulation object is a nonlinear plasma wakefield created by a particle beam. The beam is initially placed at the origin, with a distribution of n Subscript b Baseline e Superscript minus r squared divided by sigma Super Subscript r Super Superscript 2 Baseline e Superscript minus z squared divided by sigma Super Subscript z Super Superscript 2 $n_be^{-r^{2}/\sigma _{r}^{2}}e^{-z^{2}/\sigma _{z}^{2}}$ ; the plasma is initially evenly distributed with a density of n Subscript p $n_p$ . The normalised parameters are: n Subscript b Baseline divided by n Subscript p Baseline equals 10.0 $n_b/n_p=10.0$ , k Subscript p Baseline sigma Subscript r Baseline equals 0.2 $k_p\sigma _r=0.2$ and k Subscript p Baseline sigma Subscript z Baseline equals 1.0 $k_p\sigma _z=1.0$ ; k Subscript p Baseline equals omega Subscript p Baseline divided by c $k_p=\omega _p/c$ is the plasma wavenumber, and omega Subscript p Baseline equals StartRoot n Subscript p Baseline e squared divided by m Subscript e Baseline epsilon 0 EndRoot $\omega _p=\sqrt {n_pe^2/m_e\epsilon _0}$ is the plasma frequency. The program ran a single 3-D step, and the number of 2-D steps is equal to the resolution in the z $z$ direction.

Figure 4.

Simulated with the program before(left)/after(right) revision, electric field in direction z ( xi $\xi$ ), distributed across the x– xi $\xi$ plane.

Two heat maps showing electric field distribution across a plane.

Figure 5.

Comparison of simulation result before/after revision, electric field in direction z ( xi $\xi$ ), distributed along the z axis. Simulations were conducted with 64 CPU cores (for the old program) and 1, 2 and 4GPUs (for the new program).

A line graph comparing electric field in direction z along the z axis using different computational methods.

Figure 4 shows distribution of upper E Subscript z $E_z$ , simulated by the program before and after revision; figure 5 shows a comparison between on-axis lineouts. We see that there is excellent agreement between these results, and differences are negligible.

To test the speed performance of the new program, we measured the time taken to run a 2-D loop, with 4/9/25/100 particles per cell, and resolutions of 512 cubed $512^3$ (9*9*9 in the input file) and 1024 cubed $1024^3$ (10*10*10 in the input file). Each combination of ppc and resolution was measured with 8/16/32 CPU cores and 1/2/4 GPUs, though all set-ups except the largest (100 particles per cell and 1024 cubed $1024^3$ resolution) fit on a single GPU. The results are shown in figure 6. For CPU cases, the more threads we use, the less time the simulation requires; however, for more than half of the GPU cases, one GPU takes less time than 2 or 4 GPUs; also, except for the 100 particles per cell and 1024 cubed $1024^3$ resolution case, which cannot fit on a single GPU, one GPU takes less time than all CPU set-ups.

Figure 6.

Measurement of simulation 2-D time in seconds. Subplots (a), (b), (c) and (d) show set-ups with 4/9/25/100 particles per cell, respectively. Execution time is measured with 8 (blue, plain)/16 (orange, plain)/32 (green, plain) CPU threads and 1 (blue, slashed)/2 (orange, slashed)/4 (green, slashed) GPUs. In each subplot, results with a simulation resolution of 512 cubed $512^3$ are shown on the left side, and results with a simulation resolution of 1024 cubed $1024^3$ are shown on the right side.

Four bar graphs compare execution times for different grid sizes and particle per cell setups.

Also, since most problems fit on a single GPU, we also calculated and compared the node-hours spent by 32 CPU threads and 1 GPU. For CPU tasks, the number of node hours is calculated as time taken times $\times$ number of CPU cores/64; for GPU tasks, the number of node hours is calculated as time taken times $\times$ number of GPUs/4; as each CPU computation node has 64 cores, and each GPU computational node has 4 GPUs. The results are shown in figure 7.

Figure 7.

Number of node hours spent in D simulation. Subplots (a), (b), (c) and (d) show set-ups with 4/9/25/100 particles per cell, respectively. Execution time is measured with 32 (green, plain) CPU threads and 1 (red, slashed) GPU. In each subplot, results with a simulation resolution of 512 cubed $512^3$ are shown on the left side, and results with a simulation resolution of 1024 cubed $1024^3$ are shown on the right side.

Four bar graphs compare node hours spent in D simulation using CPU and GPU for different particles per cell and grid sizes.
Figure 7. Long description

Four bar graphs depict the number of node hours spent in D simulation for different setups with varying particles per cell. Panel A shows setups with 4 particles per cell, Panel B with 9 particles per cell, Panel C with 25 particles per cell, and Panel D with 100 particles per cell. Each panel compares execution time using 32 CPU threads (green, plain bars) and 1 GPU (red, slashed bars). The left side of each panel shows results with a simulation resolution of 512 cubed, and the right side shows results with a simulation resolution of 1024 cubed. The y-axes represent node hours, with separate scales for each resolution. The graphs indicate that GPU execution consistently requires fewer node hours compared to CPU execution across all setups and resolutions.

Among all cases that can be run on 1 GPU, the computational time of 32 CPU cores, measured in node hours, is 2.27–3.37 times longer than the computational time of 1 GPU. At the higher resolution ( 1024 cubed $1024^3$ ), using 1 GPU instead of multiple CPU cores also lead to a higher speedup, measured in node hours.

It would also be worthy to compare the speed of QuickPIC-GPU with that of a fully explicit PIC code with GPU support, e.g. WarpX. We measured the average number of particle advances per nanosecond, using a single V100 GPU and a resolution of 256 cubed $256^3$ (other simulation parameters were the same as the ones given above), to compare with the average number of particle advances per nanosecond when using a single V100 GPU, given in the article (Myers et al. Reference Myers2021). QuickPIC-GPU gives approximately 0.5 particle advances per nanosecond, while WarpX gives approximately 1 particle advance per nanosecond. This difference in update time could have several origins: QuickPIC-GPU uses a 2-D push, not a 3-D one, but the field solve is more complicated due to the quasi-static approximation; QuickPIC-GPU uses a spectral solver, while WarpX uses a finite difference time domain solver. We also note that this does not indicate a comparison of the total simulation time needed for the same physical process, since quasi-static approximation allows for much longer time steps while maintaining accuracy.

5. Discussion

As mentioned in § 3, when using multiple GPUs, the FFT transpose algorithm involves a large amount of data transfer across MPI ranks, in multiple steps. Therefore, we evaluated the computational cost of this transfer by measuring the time spent by these MPI calls, in 2-GPU and 4-GPU runs. The results are shown in figure 8. For comparison, the total time of 1-GPU runs (cross-rank transfers do not exist in these cases) is also shown. We see that MPI transfers take up more than half of the total time, for 2-GPU and 4-GPU runs; after subtracting the time cost of MPI in FFT transpose, the more GPUs we use, the faster other tasks are processed.

Figure 8.

Time spent on FFT cross-rank data transfer/other processes. Subplots (a), (b), (c) and (d) show set-ups with 4/9/25/100 particles per cell, respectively. Execution time is measured with 1/2/4GPUs. The FFT MPI transfer time is shown in orange, while time spent on other processes is shown in blue. In each subplot, results with a simulation resolution of 512 cubed $512^3$ are shown on the left side, and results with a simulation resolution of 1024 cubed $1024^3$ are shown on the right side.

Four bar graphs showing time spent on processes with different particles per cell and GPU configurations.

Thus, we propose that the field solver in QuickPIC needs to be updated for better performance on multiple GPUs. Suppose a $a$ MPI ranks are used; in the current algorithm, to transpose across the whole simulation box, each partition is divided into a $a$ squares in the direction different from the direction in which the simulation box is partitioned. Then, in a loop that runs a $a$ times, these squares are transposed locally, and then exchanged with squares that are in the symmetric position across the main diagonal of the matrix. Thus, the more MPI ranks used, the more steps FFT transpose would take. For the CPU version, this was not an eminent problem; MPI transfers across CPU threads (which are mostly on the same device) do not take long, compared with other processes (which were not so highly parallel as in the GPU version). However, in the case of multiple GPUs, as seen in figure 8, FFT transpose takes up more than half of the total time.

In the future, there are two optimisation plans that we shall attempt; the first is to use cuFFT, the library for GPU-accelerated FFT algorithms developed by NVIDIA, which currently supports distributing calculation across multiple GPUs, instead of the old FFT algorithm originally designed for CPUs. This would require a change in the parallelisation method of QuickPIC-GPU, and a switch from sine/cosine to e Superscript i x $e^{ix}$ / e Superscript minus i x $e^{-ix}$ as the orthogonal basis. Instead of having each GPU connects to a different CPU, and the CPUs are connected via MPI, cuFFT would require using its own application programming interface to connect multiple GPUs to a single CPU. We shall test if this is applicable to QuickPIC. The second plan is to use finite-difference field solvers that do not include a transpose across the whole simulation box, and thus require less data transfer. Also, we plan to try using the NVIDIA Collective Communication Library (NCCL) for cross-rank transfers, to see if this transfer speed can be further improved.

Besides FFT transpose, other processes also have the potential to be optimised. Currently, the functions in the 3-D loop still run on CPUs, and before and after the call to the 2-D loop in each 3-D step, required data (charge density distribution, electromagnetic fields, etc.) are transferred between the CPU and the GPU; in the future, we plan to also move 3-D particle and field solvers to GPUs, to reduce data transfer and optimise the 3-D loop.

6. Conclusions

We have developed QuickPIC-GPU, a new version of QuickPIC, that runs parallel on one or more GPUs. Referring to the UPIC framework, we use CUDA C functions to replace the Fortran 77 subroutines at the bottom layer, and MPI to communicate between multiple GPUs, each connected to a different CPU. With the same simulation parameters, the results of the new program matches excellently with the old program, using up to 4 GPUs; we see significant speedup (up to 3.37 times $\times$ ), when the simulation time (measured in node hours) of 1-GPU runs is compared against that of multiple (up to 32) CPU threads.

When using multiple GPUs, we note that FFT transpose takes up more than half of the total time; in future revisions, we shall attempt to reduce this cost by using cuFFT or finite-difference field solvers. The GPU-MPI parallelisation of 3-D loop functions, and additions of other functions (for example, simulation of laser wakefields; a CPU version is already under development), are also considered for future work.

Acknowledgements

Work supported by Strategic Priority Research Program of the Chinese Academy of Sciences, Grant No. XDB0530000, NSFC Grant 12575156 and 12075030, Beijing Normal University Scientific Research Initiation Fund for Introducing Talents 310432104 and U.S. Department of Energy (DOE), Offices of Advanced Scientific Computing Research and of High Energy Physics (HEP) through the Scientific Discovery through Advanced Computing (SciDAC) project, ComPASS 4, under Contract Nos. DE-AC02-06CH11357 and DE-AC02-05CH11231, and a FNAL subcontract 644405, and by a DOE HEP Grant DESC0010064 and by the U.S. National Science Foundation through award No. 2108970.

Editor Luís O. Silva thanks the referees for their advice in evaluating this article.

Declaration of interests

The authors report no conflict of interest.

References

Adli, E. et al. 2018 Acceleration of electrons in the plasma wakefield of a proton bunch. Nature 561, 363–367.10.1038/s41586-018-0485-4CrossRefGoogle ScholarPubMed
An, W., Decyk, V.K., Mori, W.B. & Antonsen, T.M. 2013 An improved iteration loop for the three dimensional quasi-static particle-in-cell algorithm: QuickPIC. J. Comput. Phys. 250, 165–177.10.1016/j.jcp.2013.05.020CrossRefGoogle Scholar
Arber, T.D. et al. 2015 Contemporary particle-in-cell approach to laser-plasma modelling. Plasma Phys. Control. Fusion 57, 113001.10.1088/0741-3335/57/11/113001CrossRefGoogle Scholar
Bai, J., Cao, Y., Chu, Y. & Zhang, X. 2018 An improved immersed finite element particle-in-cell method for plasma simulation. Comput. Math. Appl. 75, 1887–1899.10.1016/j.camwa.2017.08.001CrossRefGoogle Scholar
Bastrakov, S., Donchenko, R., Gonoskov, A., Efimenko, E., Malyshev, A., Meyerov, I. & Surmin, I. 2012 Particle-in-cell plasma simulation on heterogeneous cluster systems. J. Comput. Sci-NETH. 3, 474–479.10.1016/j.jocs.2012.08.012CrossRefGoogle Scholar
Birdsall, C.K. & Langdon, A.B. 2018 Plasma Physics Via Computer Simulation. CRC Press.10.1201/9781315275048CrossRefGoogle Scholar
Blumenfield, I. et al. 2007 Energy doubling of 42 gev electrons in a metre-scale plasma wakefield accelerator. Nature 445, 741–744.10.1038/nature05538CrossRefGoogle Scholar
Burau, H., Widera, R.ée, Hönig, W., Juckeland, G., Debus, A., Kluge, T., Schramm, U., Cowan, T.E., Sauerbrey, R. & Bussmann, M. 2010 PIConGPU: a fully relativistic particle-in-cell code for a GPU cluster. IEEE Trans. Plasma Sci. 38, 2831–2839.10.1109/TPS.2010.2064310CrossRefGoogle Scholar
Chen, P., Dawson, J.M., Huff, R.W. & Katsouleas, T. 1985 Acceleration of electrons by the interaction of a bunched electron beam with a plasma. Phys. Rev. Lett. 54, 693–696.10.1103/PhysRevLett.54.693CrossRefGoogle ScholarPubMed
Corde, S. et al. 2015 Multi-gigaelectronvolt acceleration of positrons in a self-loaded plasma wakefield. Nature 524, 442–445.10.1038/nature14890CrossRefGoogle Scholar
Dawson, J.M. 1983 Particle simulation of plasmas. Rev. Mod. Phys. 55, 403–447.10.1103/RevModPhys.55.403CrossRefGoogle Scholar
Decyk, V.K. 2007 Upic: A framework for massively parallel particle-in-cell codes. Comput. Phys. Commun. 177, 95–97.10.1016/j.cpc.2007.02.092CrossRefGoogle Scholar
Decyk, V.K. 2015 Skeleton particle-in-cell codes on emerging computer architectures. Comput. Sci. Eng. 17, 47–52.10.1109/MCSE.2014.131CrossRefGoogle Scholar
Decyk, V.K. 2023 Gpu skeleton code.Google Scholar
Decyk, V.K. & Singh, T.V. 2011 Adaptable particle-in-cell algorithms for graphical processing units. Comput. Phys. Commun. 182, 641–648.10.1016/j.cpc.2010.11.009CrossRefGoogle Scholar
Derouillat, J. et al. 2018 Smilei: a collaborative, open-source, multi-purpose particle-in-cell code for plasma simulation. Comput. Phys. Commun. 222, 351–373.10.1016/j.cpc.2017.09.024CrossRefGoogle Scholar
Diederichs, S., Benedetti, C., Huebl, A., Lehe, R., Myers, A., Sinn, A., Vay, J.-L., Zhang, W. & Thévenet, M. 2022 HiPACE++: a portable, 3D quasi-static particle-in-cell code. Comput. Phys. Commun. 278, 108421.10.1016/j.cpc.2022.108421CrossRefGoogle Scholar
Fonseca, R.A. et al. 2002 Osiris: A three-dimensional, fully relativistic particle in cell code for modeling plasma based accelerators. In Computational Science — ICCS (ed. Peter M. A. Sloot, Alfons G. Hoekstra, C. J. Kenneth Tan & Jack J. Dongarra), pp. 342–351. Springer.Google Scholar
Huang, C., Decyk, V.K., Ren, C., Zhou, M., Lu, W., Mori, W.B., Cooley, J.H., Antonsen, T.M. & Katsouleas, T. 2006 QUICKPIC: a highly efficient particle-in-cell code for modeling wakefield acceleration in plasmas. J. Comput. Phys. 217, 658–679.10.1016/j.jcp.2006.01.039CrossRefGoogle Scholar
Joshi, C., Corde, S. & Mori, W.B. 2020 Perspectives on the generation of electron beams from plasma-based accelerators and their near and long term applications. Phys. Plasmas 27, 070602.10.1063/5.0004039CrossRefGoogle Scholar
Joshi, C. & Katsouleas, T. 2003 Plasma accelerators at the energy frontier and on tabletops. Phys. Today 56, 47–53.10.1063/1.1595054CrossRefGoogle Scholar
Li, F., An, W., Decyk, V.K., Xu, X., Hogan, M.J. & Mori, W.B. 2021 A quasi-static particle-in-cell algorithm based on an azimuthal Fourier decomposition for highly efficient simulations of plasma-based acceleration: QPAD. Comput. Phys. Commun. 261, 107784.10.1016/j.cpc.2020.107784CrossRefGoogle Scholar
Liewer, P.C. & Decyk, V.K. 1989 A general concurrent algorithm for plasma particle-in-cell simulation codes. J. Comput. Phys. 85, 302–322.10.1016/0021-9991(89)90153-8CrossRefGoogle Scholar
Litos, M. et al. 2014 High-efficiency acceleration of an electron beam in a plasma wakefield accelerator. Nature 515, 92–95.10.1038/nature13882CrossRefGoogle Scholar
Mora, P. & Antonsen, T.M. Jr 1997 Kinetic modeling of intense, short laser pulses propagating in tenuous plasmas. Phys. Plasmas 4, 217–229.10.1063/1.872134CrossRefGoogle Scholar
Myers, A. et al. 2021 Porting WarpX to GPU-accelerated platforms. Parallel Comput. 108, 102833.10.1016/j.parco.2021.102833CrossRefGoogle Scholar
Sosedkin, A.P. & Lotov, K.V. 2016 LCODE: A parallel quasistatic code for computationally heavy problems of plasma wakefield acceleration. Nucl. Instr. Meth. Phys. Res. Sect. A. 829, 350–352.10.1016/j.nima.2015.12.032CrossRefGoogle Scholar
Tajima, T. & Dawson, J.M. 1979 Laser electron accelerator. Phys. Rev. Lett. 43, 267–270.10.1103/PhysRevLett.43.267CrossRefGoogle Scholar
Wang, T., Khudik, V., Kim, J. & Shvets, G. 2022 Accelerated three-dimensional quasistatic particle-in-cell code. Phys. Rev. Accel. Beams 25, 104603.10.1103/PhysRevAccelBeams.25.104603CrossRefGoogle Scholar
Xiong, Q., Huang, S., Yuan, Z., Sharma, B., Kuang, L., Jiang, K. & Yu, L. 2024 GPIC: a set of high-efficiency CUDA Fortran code using GPU for particle-in-cell simulation in space physics. Comput. Phys. Commun. 295, 108994.10.1016/j.cpc.2023.108994CrossRefGoogle Scholar
Xiong, Q.Y., Huang, S.Y., Yuan, Z., Kuang, L., Sharma, B., Jiang, K., Wu, H.H. & Yu, L. 2025 Applications of multi-GPU computing in large-scale fully kinetic particle-in-cell simulations of space plasma. AIP Adv. 15, 025201.10.1063/5.0251687CrossRefGoogle Scholar
Figure 0

Figure 1. Figure 1 long description.Flow chart of 3-D–2-D loops in QuickPIC.

Figure 1

Figure 2. Class diagram of QuickPIC. Classes outlined in red are the main GPU revision targets.

Figure 2

Figure 3. Figure 3 long description.Mapping of CUDA architecture to particle simulation.

Figure 3

Figure 4. Simulated with the program before(left)/after(right) revision, electric field in direction z (ξ$\xi$), distributed across the x–ξ$\xi$ plane.

Figure 4

Figure 5. Comparison of simulation result before/after revision, electric field in direction z (ξ$\xi$), distributed along the z axis. Simulations were conducted with 64 CPU cores (for the old program) and 1, 2 and 4GPUs (for the new program).

Figure 5

Figure 6. Measurement of simulation 2-D time in seconds. Subplots (a), (b), (c) and (d) show set-ups with 4/9/25/100 particles per cell, respectively. Execution time is measured with 8 (blue, plain)/16 (orange, plain)/32 (green, plain) CPU threads and 1 (blue, slashed)/2 (orange, slashed)/4 (green, slashed) GPUs. In each subplot, results with a simulation resolution of 5123$512^3$ are shown on the left side, and results with a simulation resolution of 10243$1024^3$ are shown on the right side.

Figure 6

Figure 7. Figure 7 long description.Number of node hours spent in D simulation. Subplots (a), (b), (c) and (d) show set-ups with 4/9/25/100 particles per cell, respectively. Execution time is measured with 32 (green, plain) CPU threads and 1 (red, slashed) GPU. In each subplot, results with a simulation resolution of 5123$512^3$ are shown on the left side, and results with a simulation resolution of 10243$1024^3$ are shown on the right side.

Figure 7

Figure 8. Time spent on FFT cross-rank data transfer/other processes. Subplots (a), (b), (c) and (d) show set-ups with 4/9/25/100 particles per cell, respectively. Execution time is measured with 1/2/4GPUs. The FFT MPI transfer time is shown in orange, while time spent on other processes is shown in blue. In each subplot, results with a simulation resolution of 5123$512^3$ are shown on the left side, and results with a simulation resolution of 10243$1024^3$ are shown on the right side.