arXiv CSopen access

A HIP-Compatible Accelerator Backend for Fourier-Bessel Particle-in-Cell Simulations on CPU/DCU Heterogeneous Clusters

A HIP-Compatible Accelerator Backend for Fourier-Bessel Particle-in-Cell Simulations on CPU/DCU Heterogeneous Clusters Jingliang Fana,b , Ruiqing Hec , Yang Wand , Jiandong Shanga , Hengliang Guoa , Qiang Chena,e,∗ a National Supercomputing Center in Zhengzhou, Zhengzhou University, Zhengzhou, 450001, China b School of Computer and Artificial Intelligence, Zhengzhou University, Zhengzhou, 450001, China c School of Communication and Artificial Intelligence, School of Integrated Circuits, Nanjing Institute of Technology, Nanjing, 211167, China d School of Physics and Laboratory of Zhongyuan Light, Zhengzhou University, Zhengzhou, 450001, China

arXiv:2609.06680v1 [physics.comp-ph] 6 Sep 2026

e Laboratory for Advanced Computing and Intelligence Engineering, Wuxi, 214000, China

Abstract FBPIC (Fourier-Bessel particle-in-cell) is a high-performance simulation code for relativistic plasma and accelerator physics. Its original accelerator backend relies on Numba CUDA, which limits its direct deployment on accelerators using the HIP (Heterogeneous-Compute Interface for Portability) programming environment, such as DCU (Deep Computing Unit) accelerators. In this work, we develop an accelerator backend compatible with HIP that enables FBPIC to run efficiently on DCU platforms while preserving its Python user interface and high level simulation workflow. For the evaluated LWFA (laser-wakefield acceleration) workloads, the proposed backend achieves 1.32-1.54× speedups over the original FBPIC implementation on an NVIDIA V100 GPU and enables efficient execution on the DCU platform. We also summarize the key lessons learned from porting FBPIC to the DCU platform. Multi-DCU experiments achieve a 1.88× strong-scaling speedup on four accelerators and a 2.72× increase in aggregate throughput at approximately 68% weak-scaling efficiency, with communication analysis identifying inter-node communication and synchronization as the main scalability limitations. Beyond FBPIC, the proposed approach provides a practical reference for porting and optimizing other scientific computing applications developed with Python on heterogeneous accelerator platforms. Keywords: Fourier-Bessel particle-in-cell, HIP, CuPy, DCU, Performance portability

1. Introduction The PIC (Particle-In-Cell) method has been widely used for plasma physics, accelerator physics, and laser-plasma interaction simulations[1, 2, 3]. These problems involve electromagnetic structures and charged particles evolving over multiple spatial and temporal scales, and their predictive accuracy is often limited by both numerical dispersion and computational cost. As plasma-accelerator studies move toward larger domains, longer propagation distances, and more extensive parameter scans, accelerator computing has become increasingly important for high-fidelity PIC simulations[3, 4]. The DCU (Deep Computing Unit) is a GPU-like accelerator for highperformance heterogeneous computing and supports the HIP (Heterogeneous-computing Interface for Portability) programming model[5]. Its CUDA-like programming interface provides a practical environment for porting GPU-accelerated scientific applications[6, 7]. With the emergence of such accelerator platforms, improving the performance portability of highfidelity PIC codes becomes increasingly important. FBPIC (Fourier-Bessel Particle-In-Cell) is a Particle-In-Cell code designed for relativistic plasma simulations, with particular relevance to laser-wakefield and plasma-wakefield acceleration ∗ Corresponding author

Email address: [email protected] (Qiang Chen)

problems[8]. By representing the electromagnetic fields in a quasi-cylindrical spectral basis, FBPIC can retain essential three-dimensional physical effects for nearly axisymmetric systems while requiring far fewer degrees of freedom than a full three-dimensional Cartesian PIC simulation[8]. The FourierBessel spectral formulation also avoids the numerical dispersion associated with conventional finite-difference time-domain PIC solvers, which is particularly important for simulations involving ultra-relativistic beams and laser pulses[8, 9]. The original FBPIC implementation is written in Python and relies on Numba just-in-time compilation for its CPU and CUDA GPU execution paths. Its GPU support requires both CuPy and Numba-related GPU dependencies[10, 11, 12]. This design provides a productive high-level programming model and good performance on NVIDIA GPUs, but it also tightly couples the performance-critical kernels to the CUDA-oriented Numba backend. This becomes a major limitation when targeting DCU/HIP platforms, because Numba’s ROCm target has been officially unmaintained since version 0.54.0[13]. As a result, the original Numba-based GPU path cannot be directly reused as a robust and efficient backend for DCU systems. To overcome this portability limitation, we develop a HIPcompatible accelerator backend for FBPIC that removes the dependence on Numba from its performance-critical GPU execution path. Instead of translating Python-defined kernels

through the Numba JIT infrastructure, the dominant PIC kernels are reimplemented as explicit C/C++ GPU kernels and integrated into the existing Python framework through CuPy RawKernel[14]. CuPy provides CUDA- and ROCm-backed execution environments, allowing the same Python-side dispatch mechanism to be retained across NVIDIA GPUs and HIP-compatible accelerators[12, 15]. In addition, particle and field data remain resident in CuPy device arrays throughout the time-stepping loop, so the backend can be replaced without introducing additional host-device data movement or modifying the high-level simulation workflow. The existing MPI-based longitudinal domain decomposition is preserved, so the proposed backend changes only the intra-device execution layer rather than the physical model or distributed-memory decomposition of FBPIC. The main contributions of this work are summarized as follows: A portable accelerator backend for FBPIC. We develop a HIP-compatible accelerator backend for FBPIC by replacing its performance-critical Numba CUDA kernels with explicitly implemented C/C++ GPU kernels. These kernels are integrated into the existing Python framework through CuPy RawKernel, preserving the original user interface, simulation workflow, device-resident data model, and MPI-based domain decomposition. Kernel-level optimization for DCU execution. We optimize the rewritten kernels for DCU execution through compiler- and occupancy-oriented techniques, including pointer-based memory access, selective loop unrolling, register-pressure control, and kernel-specific thread-block tuning, and investigate the performance trade-offs associated with particle locality and binning. Cross-platform numerical and performance evaluation. We verify the numerical consistency of the backend using linearwakefield and LWFA (laser-wakefield acceleration) benchmarks and demonstrate performance portability across CUDA and HIP platforms. On an NVIDIA V100 GPU, the rewritten implementation achieves 1.32-1.54× speedups over the original FBPIC implementation for the evaluated workloads. Multi-DCU scalability and communication analysis. We characterize the multi-DCU scalability of FBPIC through strong- and weak-scaling experiments and mpiP communication profiling. The implementation achieves a 1.88× strongscaling speedup on four DCUs and a 2.72× increase in aggregate throughput at approximately 68% weak-scaling efficiency, while the profiling results identify inter-node communication and synchronization as the primary limitations at larger process counts.

speed laser propagation, finite-difference field solvers may introduce numerical dispersion and numerical Cherenkov-type artifacts[16, 17]. Spectral and pseudo-spectral solvers have therefore become important alternatives. FBPIC follows this direction by adopting a Fourier-Bessel quasi-cylindrical spectral formulation, which is particularly effective for close-toaxisymmetric laser wakefield and plasma wakefield acceleration problems. Its spectral cylindrical representation provides high accuracy at substantially lower cost than full threedimensional Cartesian PIC for such geometries[8]. Modern high-performance PIC codes increasingly target many-core and GPU-based supercomputers. WarpX is an exascale-oriented electromagnetic PIC code built on the AMReX framework. Its performance-portable implementation targets NVIDIA, AMD, and Intel GPUs and has demonstrated scalability on leadership-class systems[18, 19]. PIConGPU represents another major GPU-oriented development path. It was originally developed as a CUDA-based PIC implementation for GPU clusters[20]. It was subsequently moved toward performance portability using the cupla interface over Alpaka, enabling single-source C++ kernels to target different processor architectures[21]. Smilei is a collaborative, open-source C++ PIC code co-developed by plasma physicists and HPC specialists, with an emphasis on modularity, multiphysics capabilities, and parallel performance[22]. HiPACE++ further demonstrates the effectiveness of performance-portable GPU implementations for quasi-static plasma-accelerator modeling, reporting near-optimal strong scaling from 1 to 512 GPUs[23]. These developments illustrate a broader shift in highperformance PIC software toward accelerator-aware and performance-portable implementations[18, 21, 23]. Representative implementations are predominantly developed in C++ and rely on portability layers such as AMReX or Alpaka to map computational kernels onto different processor architectures[24, 21]. FBPIC differs from these codes in that it combines a high-level Python workflow with a spectral quasicylindrical algorithm. The original implementation relies on Numba JIT compilation for performance and can run on multicore CPUs or GPUs, with GPU execution being much faster for large simulations[25]. This Python-JIT design provides productivity and flexibility, but it also limits portability. In particular, the original GPU backend is closely tied to Numba CUDA. This becomes problematic for AMD/DCU-class accelerators because Numba’s ROCm target was moved to an unmaintained status in version 0.54.0 and relocated outside the main Numba repository[13]. A recent experimental FBPIC pull request explored the use of the numba.hip interface with numba.hip.pose_as_cuda() in order to minimize changes to the existing CUDA-oriented implementation[26]. The reported tests on one GPU die of an MI250X identified several interoperability and functionality limitations and showed lower-than-expected performance[26]. These observations suggest that a minimal-change Numba-HIP compatibility layer alone may be insufficient to provide a robust and high-performance FBPIC backend for production use. CuPy provides another possible route for Python-based GPU computing. Its RawKernel interface allows user-defined kernels

2. Current state of the art Considerable effort has been devoted to improving the accuracy and scalability of PIC simulations. Conventional Cartesian finite-difference PIC solvers remain widely used because of their algorithmic locality and compatibility with domain decomposition. However, for relativistic beams and high2

written in raw CUDA source to be compiled and launched from Python, while caching the compiled binary for reuse[14]. CuPy also provides experimental ROCm support. It provides a practical integration layer through which carefully written C/C++ GPU kernels can be connected to a Python scientific application. The present work focuses on this unresolved gap. While preserving the usability of FBPIC, this work implements a version that can run efficiently across multiple GPU platforms by leveraging the RawKernel interface provided by CuPy. This study complements existing research on PIC code portability and provides a practical reference for the efficient and robust migration of other Python-based scientific computing applications to emerging GPU-like accelerator platforms such as DCUs.

Charge and current deposition. Particle charge and current are projected back to the interpolation grids. Similar to gathering, each particle contributes to several radial and longitudinal grid points for every azimuthal mode. Unlike gathering, deposition performs concurrent updates to shared grid locations. Multiple particles may therefore write to the same element, requiring atomic accumulation or an equivalent conflict-resolution strategy. Spectral source correction and field advancement. The deposited sources are transformed to spectral space, corrected to satisfy the selected current-continuity treatment, and used to advance Maxwell’s equations with the pseudo-spectral analytical time-domain method. The spectral representation combines Fourier transforms along z with Hankel transforms along r. This spectral field-solver formulation avoids the numerical dispersion commonly associated with conventional finitedifference field solvers. Boundary and particle exchange. For multi-device simulations, guard-cell values are exchanged between adjacent longitudinal subdomains, deposited source quantities are accumulated across overlapping regions, and particles leaving a local subdomain are transferred to the neighboring rank. The communication volume is determined by the guard-region width and the number of migrating particles rather than the total local particle population. Local packing, unpacking, buffer initialization, and particle rearrangement remain accelerator-intensive operations and must therefore be included in the portable kernel backend.

3. Parallel implementation of FBPIC on DCU accelerators 3.1. Computational structure and porting requirements of FBPIC FBPIC is organized as a Python-controlled, object-oriented simulation framework. The top-level Simulation object owns the field, particle, diagnostic, and boundary-communication components and invokes the PIC cycle through its step method. The high-level Python layer is primarily responsible for physical configuration, data organization, and execution control, whereas the computationally intensive operations are implemented as device kernels or GPU array operations. In the original implementation, GPU execution requires both CuPy and Numba CUDA[10]. FBPIC exposes two levels of parallelism. At the inter-device level, the computational domain is decomposed along the longitudinal direction, and neighboring MPI ranks exchange field guard cells and particles that cross subdomain boundaries. At the intra-device level, particles and grid elements within each subdomain are processed concurrently by GPU threads[27, 28]. The present work preserves the existing longitudinal domain decomposition and MPI-level simulation semantics, while replacing the intra-device Numba implementation and the associated buffer-processing kernels. This separation limits platform-specific modifications to the accelerator backend and avoids changes to the physical decomposition of the original code. The five major computational stages of a typical FBPIC time step are described below, and their key characteristics are summarized in Table 1. Field gathering. The electric and magnetic fields stored on the interpolation grids are evaluated at the particle positions. Because FBPIC represents the fields as a truncated set of azimuthal modes, the local field experienced by each particle is reconstructed by accumulating contributions from all retained modes. Gathering is dominated by irregular read accesses and repeated interpolation operations. Particle push. After gathering, the particle momenta and positions are advanced according to the Lorentz force. FBPIC uses a time-centered PIC sequence in which the momentum and position updates are staggered relative to the field quantities. The update of each macroparticle is independent, which provides abundant data parallelism.

3.2. Kernel execution layer for DCU platforms DCU accelerators are programmed through a HIPcompatible heterogeneous-computing environment. HIP adopts a CUDA-like host/device model in which a kernel launch is defined by a grid of thread blocks, and the threads within a block can cooperate through synchronization and on-chip shared memory. This similarity reduces the structural effort required to translate CUDA C/C++ kernels, but it does not make Python accelerator backends automatically portable. The original Numba-based GPU path cannot be reused directly in the target environment. Numba officially classified its ROCm target as unmaintained in version 0.54.0 and moved the corresponding implementation outside the main repository. Consequently, the original FBPIC kernels cannot be used directly as a stable and performance-controllable execution path in the DCU software environment. To remove this dependency, the performance-critical kernels are rewritten as explicit C/C++ GPU kernels using the common subset of CUDA and HIP kernel syntax. The kernels retain the grid-block-thread programming model but no longer depend on Python JIT translation. More importantly, it separates kernel implementation from the Python compiler ecosystem and places the computational core closer to the native programming model of both CUDA and HIP devices. CuPy RawKernel is used as the interface between the rewritten kernels and the FBPIC Python layer. Compilation occurs on first use and is skipped when a valid cached binary is available. On CUDA builds, RawKernel provides NVRTC and NVCC 3

Table 1. computational characteristics: N p denotes the number of macroparticles in a local subdomain, Nm the number of retained azimuthal modes, Nz and Nr the longitudinal and radial grid dimensions, respectively, and S the number of grid points covered by the particle shape function. The dominant per step operations can then be grouped into five computational motifs.

Component Field gather Particle push Charge/current deposition Fourier-Hankel transforms PSATD field update

Approximate work

Parallel pattern

O(N p Nm S ) O(N p ) O(N p Nm S ) O[Nm (Nr Nz log Nz + Nz Nr2 )] O(Nm Nz Nr )

One or more threads per particle Independent particle updates Particle-parallel scattered updates Mode and grid parallel Independent spectral-grid updates

4. Lessons from the DCU port

compilation backends, whereas CuPy’s ROCm build provides a HIP-backed implementation of the same kernel compilation interface[14, 15]. This allows the Python-side kernel invocation mechanism to be retained across CUDA and ROCm builds. Thus, a kernel written within the CUDA/HIP-compatible language subset can be compiled through either the CUDA or ROCm toolchain without changing the FBPIC Python call site. Figure 1 presents the software architecture of the ported FBPIC implementation. The architecture consists of three layers: the unchanged Python application layer, a backend integration layer, and the device-kernel layer. Python application and simulation-control layer. The simulation continues to be initiated by a standard FBPIC input script. Users define the computational domain, particle species, laser parameters, boundary conditions, diagnostics, and timeadvancement settings through the original Python interface. The Simulation, Particles, Fields, and BoundaryCommunicator objects retain their roles in data ownership and PIC-cycle scheduling. Existing physical models and input scripts therefore require little or no modification. Backend dispatch and kernel-construction layer. The original Numba launch sites are redirected to a CuPy-based backend. This layer selects the appropriate kernel variant, prepares scalar parameters, determines the grid and thread-block configuration, and constructs or retrieves the corresponding RawKernel object. Platform-dependent compiler options and kernel specializations are confined to this layer. The high-level particle and field modules do not need to distinguish between a CUDA GPU and a DCU/HIP device. Device-data and kernel layer. Particle and field state remains represented by CuPy arrays throughout the simulation. The same device allocation is therefore visible to CuPy array operations, FFT routines, and RawKernel launches, avoiding redundant device copies and preserving compatibility with the existing FBPIC storage model. The architecture follows three design principles. First, the user-facing Python interface and the existing physical data model are preserved. Second, particle and field data remain device-resident throughout the PIC loop, minimizing hostdevice traffic. Third, common kernel source is shared between CUDA and HIP wherever practical, while architecture-specific variants are permitted when required for correctness or performance. Together, these principles provide a maintainable path for extending FBPIC to DCU accelerators without sacrificing its Python-based usability or preventing low-level kernel optimization.

Porting FBPIC to the DCU platform while achieving functional portability is insufficient to fully unleash the performance of the new accelerator architecture; therefore, this paper explores the following three aspects on the DCU platform. 4.1. Portable kernel reimplementation and compiler optimization The performance-critical Numba-CUDA kernels were systematically reimplemented as C/C++ GPU kernels using a restricted programming subset accepted by both the CUDA and HIP compilation environments. The rewritten kernel exposes information to the compiler that is difficult to express through high-level array abstractions. Device-array accesses are replaced by explicit pointer arguments, with readonly inputs declared as const T* and non-aliasing arguments marked with __restrict__ where the non-aliasing assumption is valid. Restrict-qualified pointers expose non-aliasing information to the compiler, enabling more aggressive code reordering, common-subexpression elimination, and reuse of loaded values[29, 30]. Small device helpers are inlined selectively, and thread indices, grid dimensions, and kernel arguments are made explicit. These changes can reduce redundant address calculations and enable more aggressive load reuse and commonsubexpression elimination. They also make generated instructions, register usage, memory transactions, and spill behavior more directly accessible to compiler reports and profiling tools[31]. The original FBPIC GPU backend first compiles Pythondefined kernels into PTX using Numba and subsequently loads and launches the generated PTX through CuPy[32]. Compared with launching kernels directly through the Numba runtime, this design reduces host-side launch overhead by exploiting CuPy’s lower-overhead execution interface, type-based kernel caching, and preprocessing of array metadata before kernel invocation. In the proposed implementation, the performancecritical kernels are instead expressed directly as C/C++ GPU source code and invoked through CuPy’s runtime-compilation interface. On CUDA platforms, the kernel source is compiled into device code through the CUDA runtime compilation toolchain, cached for subsequent invocations, and launched by directly passing the underlying pointers of CuPy arrays together with scalar arguments. On HIP-based platforms, the Pythonside invocation interface remains unchanged, while CuPy maps the kernel compilation and runtime operations to the HIP backend, which generates device-specific code for hip-compatible 4

Figure 1. CuPy based accelerator backend of the ported FBPIC implementation. Dashed boxes and arrows denote runtime compilation and binary caching performed upon the first invocation, whereas solid arrows indicate operations repeated at each time step. The cached binary is reused by the custom-kernel path, and both paths operate on device-resident CuPy arrays on a CUDA GPU or DCU/HIP device. Each time step concludes with synchronization and guard exchange, followed by diagnostics, before control returns to Simulation for the next time step.

or DCU accelerators. This unified execution path eliminates the intermediate Python-to-PTX translation stage and reduces the amount of framework-level processing required during kernel preparation and dispatch. Measurements on the evaluated platforms show that the revised backend further decreases firstuse compilation latency and host-side kernel-launch overhead, while preserving a consistent high-level programming interface across CUDA and HIP environments.

Registers are allocated from finite vector and scalar register files shared by the wavefronts resident on a compute unit. An increase in per-thread register usage therefore reduces the number of wavefronts that can be scheduled concurrently and may weaken the ability of the hardware to hide memory and arithmetic latency. If the compiler cannot accommodate the required live variables in the available register file, some values may be spilled to scratch memory. Because scratch memory is backed by the device memory hierarchy rather than by the register file, spilling introduces additional load and store instructions and can substantially increase kernel latency[35, 34]. Loop unrolling illustrates the resulting trade-off between instruction efficiency and resource consumption. The interpolation stencil used by the cubic gather operation has a fixed extent known at compilation time. Selectively unrolling such short loops can eliminate loop-control instructions, reduce repeated index calculations, and expose additional instruction-level parallelism to the compiler. However, aggressive unrolling may increase the number and live ranges of intermediate values, thereby increasing register pressure and potentially causing reduced occupancy or spilling[35, 34]. For this reason, the DCU implementation employs selective unrolling. Small, fixed-trip-count loops on frequently executed paths are unrolled only when the resulting register allocation remains below an occupancy-relevant threshold. Loops associated with a larger number of intermediate quantities, more complex control flow, or infrequently executed paths retain their iterative form. Variable scopes and temporary values are also restricted where possible to shorten live ranges and allow the compiler to reuse registers. The purpose of these transforma-

4.2. Occupancy optimization Most computational kernels in FBPIC exhibit substantial thread level parallelism; however, their practical performance is jointly influenced by thread block size, per-thread register usage, memory access patterns, and the overhead of atomic operations[33, 34, 35]. Accordingly, this section investigates occupancy-oriented execution configurations on the DCU platform by accounting for the distinct computational and memory access characteristics of different kernel classes. 4.2.1. Register pressure Register pressure is a major determinant of the achievable occupancy of particle kernels. In FBPIC, a particle-processing thread may simultaneously retain the particle position, momentum, charge or macroparticle weight, interpolation coefficients, multiple electromagnetic-field components, and intermediate quantities associated with coordinate transformations or particle updates. The live working set is particularly large in fieldgather and current-deposition kernels, where several interpolation weights and accumulation variables may remain active over an extended instruction sequence. 5

tions is not simply to minimize the static instruction count, but to obtain a more favorable balance between instruction-level parallelism and resident wavefront parallelism.

where a value of S k represents the best-performing configuration for a given kernel. The tuning benefit shown in Figure 2(b) is quantified as Sk =

4.2.2. Thread-block configuration optimization The thread-block configuration determines how GPU threads are organized and therefore affects wavefront utilization, memory-access coalescing, resource allocation, and latency hiding[33, 36, 37, 38]. On the target DCU architecture, threads are scheduled in wavefronts of 64 work-items. Block sizes that are multiples of 64 are therefore generally preferable because they avoid partially occupied wavefronts[5]. Nevertheless, increasing the number of threads per block does not necessarily improve performance. Large blocks may increase register and shared-memory consumption, reduce the number of concurrently resident blocks, and consequently limit the occupancy available for hiding memory and instruction latencies. Moreover, for multidimensional kernels, the block shape influences the mapping between threads and array dimensions and thus affects address continuity, memory coalescing, and atomicupdate contention. The tuning experiments were conducted using a relatively large-scale LWFA simulation. With the original FBPIC thread configuration, the simulation required approximately 150 ms per time step. This workload was selected because it provides sufficient particles and grid points to expose a large number of thread blocks, allowing the major kernels to operate in a throughput-oriented regime rather than being dominated by kernel-launch latency or insufficient parallelism. It also exercises the principal computational stages of FBPIC, including field gathering, particle pushing, charge and current deposition, data-layout transformation, and spectral field operations. The selected workload therefore provides a throughput-oriented test case for evaluating kernel-specific launch configurations under the target DCU architecture and software environment. For each kernel, a set of candidate configurations was constructed according to its indexing dimensionality and computational structure. One-dimensional particle and grid kernels were evaluated primarily using block sizes ranging from 128 to 512 threads, whereas multidimensional copy, transformation, and deposition kernels were additionally tested with different block shapes. All candidates preserved the numerical algorithm, data layout, and physical simulation parameters; only the kernel-launch geometry was modified. To exclude runtime compilation and first-launch overheads, each configuration was warmed up before measurement and repeatedly executed under identical problem sizes, numerical precision, and compiler settings. The execution time of each candidate configuration was measured using the DCU profiling tool, and the results were normalized to the best-performing configuration for each kernel. Let T k,c denote the execution time of kernel k using candidate configuration c. The normalized runtime shown in Figure 2(a) is defined as Rk,c =

T k,c . minc′ T k,c′

maxc T k,c . minc T k,c

(2)

represents the runtime ratio between the original and best candidate configurations. As shown in Figure 2, sensitivity to thread-block configuration varies considerably among kernels. Across the 11 representative kernels, the runtime ratio between the worst and best configurations ranges from 1.01 to 3.49. Memory-layout transformation and deposition kernels exhibit the highest sensitivity, whereas several particle and vector kernels maintain near-optimal performance over a relatively broad range of block sizes. The application-level benefit of the selected configurations was evaluated using the complete LWFA time-stepping loop. Relative to the launch configurations used by the original FBPIC implementation, the combined kernel-specific configuration reduced the average time per simulation step from 149.171 ms to 112.918 ms. This corresponds to a time reduction of 24.30%, an overall time-step speedup of 1.32×. 4.3. Particle locality Particle reordering is widely used in accelerator-oriented PIC implementations to restore spatial locality. By grouping particles that occupy the same or neighboring cells, subsequent particle kernels operate on spatially coherent subsets of the particle distribution. This improves the locality of field gathering. In WarpX, periodic particle sorting was shown to substantially improve cache reuse and the performance of field gathering and current deposition on GPUs[18].However, particle sorting itself introduces additional work. A conventional sorting procedure requires the spatial key of every particle to be generated, a permutation to be constructed, and multiple particle attribute arrays to be physically rearranged. Because particle positions change continuously, this organization must be periodically reconstructed. For simulations with a large number of macroparticles or species carrying many auxiliary attributes, the associated memory traffic can become non-negligible and partially offset the performance gained from improved locality. We examine the particle organization strategy used by FBPIC and develop a lightweight binning sorting alternative that avoids physically reordering all particle attributes. The following subsections first discuss the deposition algorithm, which motivates the need for spatial particle organization, and then compare the original sorting procedure with the proposed binning/counting approach. 4.3.1. Current deposition Among the major stages of a PIC timestep, particle pushing is comparatively straightforward to parallelize because the state of each particle can largely be updated independently. Field gathering is also naturally particle parallel: each execution thread reads grid values surrounding a particle and writes the interpolated fields only to that particle. Current and charge

(1)

6

Figure 2. Kernel-specific thread-block tuning results for FBPIC on the DCU platform using a large-scale LWFA workload. (a) Normalized execution times of candidate block configurations for representative kernels. The color scale denotes the runtime relative to the minimum runtime of each kernel, and orange boxes indicate the selected configurations. Candidates are ordered according to the total number of threads per block and then by block shape. (b) worst-to-best speedup and the selected block configuration for each kernel.

ing all particle contributions as independent global-memory updates. Contributions from the particles associated with a cell are first accumulated into thread-local variables, after which the accumulated quantities are written to the surrounding grid points. Atomic operations remain necessary because the deposition stencils of neighboring cells overlap, but a potentially large number of particle contributions can be locally aggregated before the corresponding global updates are issued. Particles located within the same cell reuse closely related grid data and share the same deposition neighborhood, allowing their contributions to be processed in a structured manner.

deposition, in contrast, perform the reverse particle-to-grid operation and are therefore more challenging on massively parallel architectures. Contributions from multiple particles may overlap on the same grid locations, creating concurrent updates that must be handled without introducing data races[39, 40, 18]. A straightforward particle-parallel implementation can resolve these conflicts using global-memory atomic operations. Although atomics preserve correctness, their performance depends strongly on the spatial distribution of particles. When many concurrently processed particles contribute to overlapping grid points, the corresponding memory updates contend for the same locations, limiting the amount of effective parallelism. Alternatively, private or shared-memory accumulation buffers can reduce the frequency of global atomic updates, but their storage requirements increase rapidly with the grid stencil, number of field components, and number of particles processed concurrently. Excessive use of such buffers can increase shared-memory or register pressure and consequently reduce kernel occupancy. FBPIC adopts a different organization that exploits particle locality at the cell level. Particles are first classified according to the cell involved in their deposition stencil. The resulting cell indices are ordered, and a prefix-sum array records the range of particles belonging to each cell. The GPU deposition kernel can therefore assign work at the cell level rather than treat-

4.3.2. Particle sorting and binning Binning and cell-based particle organization are commonly used in GPU PIC implementations to trade sorting overhead for improved locality[40, 18]. We investigated a binning/counting strategy for the DCU backend. Instead of constructing a complete ordering and physically rearranging every particle attribute, particles are first mapped to spatial bins. A bin may correspond to an individual grid cell or, more generally, to a super-cell composed of several neighboring cells. A counting kernel determines the number of particles assigned to each bin, and an exclusive prefix scan converts these counts into offsets in an auxiliary index space. Particle indices are then placed into their corresponding bin ranges. Subsequent kernels can traverse 7

normalized amplitude a0 = 0.01, a waist of 20 µm, and a pulse length cτ = 6 µm. Each simulation was advanced for 1500 time steps using a moving window that propagated at the speed of light. The physical parameters, grid, particle loading, time step, and diagnostics were identical across the three implementations. Following the acceptance criterion of the official benchmark, we normalized the maximum absolute error by the maximum analytical field amplitude. For Ez , the resulting errors were 7.408%, 7.409%, and 7.410% for Original FBPIC, Ported CUDA, and Ported DCU, respectively, below the prescribed 8% limit. The corresponding Er errors were 6.508%, 6.509%, and 6.507%, below the 11% limit. Thus, all three implementations satisfied the analytical verification criteria. We separately quantified implementation-level differences using the relative L2 metric

particles according to these ranges while leaving the original particle attribute arrays unchanged. The principal advantage of this approach is that the volume of data physically moved during particle organization is substantially reduced. The method reorganizes a compact array of particle indices rather than repeatedly permuting all position, momentum, weight, and auxiliary particle arrays. Its complexity is also closer to linear in the number of particles and bins, consisting primarily of particle classification, histogram construction, prefix scan, and index generation. The counting stage introduces concurrent updates when multiple particles are assigned to the same bin and therefore generally requires atomic increments or an equivalent privatized histogram scheme. The prefix scan and construction of the binned-index array introduce additional global-memory passes. Moreover, because particle attributes remain in their original locations, processing particles through an index array introduces a level of indirection that can weaken the memory-access regularity of particle attributes compared with a fully reordered structure. Our measurements show that, despite reducing the amount of full-array particle movement, the binning/counting implementation results in a total timestep cost comparable to that of the original physical-reordering scheme. This observation indicates that the cost removed from particle-array permutation is largely replaced by histogram construction, prefix processing, index generation, and indirect memory accesses in subsequent kernels. More importantly, it demonstrates that minimizing sorting complexity in isolation does not necessarily minimize end-to-end PIC runtime.

ϵ2 (u) =

∥u − uorig ∥2 , ∥uorig ∥2

(3)

where uorig denotes the field produced by Original FBPIC. For Ported CUDA, ϵ2 was 5.85 × 10−6 for Ez and 1.04 × 10−4 for Er . For Ported DCU, the corresponding values were 7.98×10−6 and 1.25 × 10−4 . These values are substantially smaller than the differences between the numerical and analytical solutions. Figure 3 compares the two-dimensional distributions of Ez and Er . The original CUDA implementation, the ported CUDA implementation, and the DCU implementation reproduce the same longitudinal oscillation, radial field distribution, and spatial localization as the analytical solution. In particular, the positions of the accelerating and decelerating phases in Ez , together with the alternating radial structure of Er , remain unchanged after the kernel reimplementation. No backend dependent displacement, phase shift, or distortion of the wakefield structure is observed. The field distributions produced by the ported backend on CUDA and DCU are also visually indistinguishable from those obtained with the original FBPIC implementation. The off-axis lineouts at r = 5.25 µm in Figure 4 provide a more sensitive comparison. The three numerical curves overlap throughout the displayed longitudinal interval, including the extrema and zero crossings. Their similar deviations from the analytical solution indicate that the backend replacement introduced no additional systematic error pattern at this resolution.

5. Verification of numerical consistency Before evaluating performance, we assessed whether replacing the Numba-based GPU kernels with explicit C/C++ kernels changed the numerical results of FBPIC. We used two complementary tests. The linear-wakefield benchmark compares each implementation with an analytical solution, whereas the nonlinear LWFA benchmark evaluates cross-backend consistency during coupled particle-field evolution. Hereafter, the upstream Numba-CUDA implementation is denoted as Original FBPIC, the ported backend running on NVIDIA GPUs as Ported CUDA, and the ported backend running on the DCU/HIP platform as Ported DCU. 5.1. Linear wakefield benchmark We first evaluated the numerical accuracy and cross-backend consistency of the three implementations using the FBPIC linear-wakefield verification test[41]. This test exercises the complete PIC cycle by simulating a linear laser-driven plasma wakefield and comparing the longitudinal and radial electric fields, Ez and Er , with analytical reference solutions derived in the linear wakefield regime [42]. The benchmark used Nm = 2 azimuthal modes and a linearly polarized Gaussian laser pulse. The domain contained Nz = 800 longitudinal and Nr = 120 radial grid points, with longitudinal and radial extents of 40 µm and 60 µm, respectively. The plasma density was 8 × 1024 m−3 . The laser had a

5.2. LWFA benchmark We further evaluated the numerical consistency of the ported backend using a representative LWFA simulation. Compared with the linear wakefield benchmark, this case involves a substantially higher laser intensity and a stronger plasma response, thereby providing a more representative application level test of the coupled particle and field evolution in FBPIC. The simulation is based on the representative LWFA input example distributed with FBPIC[43]. The computational domain contains Nz = 800 longitudinal and Nr = 50 radial grid points, with Nm = 2 azimuthal modes. The longitudinal simulation window extends from −10 µm to 30 µm, while the radial extent is 8

Figure 3. Two-dimensional r-z distributions of (a-d) the longitudinal electric field Ez and (e-h) the radial electric field Er for the Nm = 2 linear-wakefield benchmark after 1500 time steps. From left to right, the columns show the analytical solution, Original FBPIC, Ported CUDA, and Ported DCU. A common symmetric color scale is used within each row.

Figure 4. Longitudinal lineouts at r = 5.25 µm for the Nm = 2 linear-wakefield benchmark: (a) Ez and (b) Er . The curves show Original FBPIC, Ported CUDA, Ported DCU, and the analytical solution.

20 µm. The plasma electron density is 4 × 1018 cm−3 , with two macroparticles per cell in both the longitudinal and radial directions and four particles in the azimuthal direction. A Gaussian laser pulse with normalized amplitude a0 = 4, waist w0 = 5 µm, and duration τ = 16 fs is used to drive the wakefield. A 40 µm linear density up-ramp is applied at the entrance of the plasma, and the simulation employs a moving window propagating at the speed of light. The interaction length is 50 µm. The same physical parameters, discretization, and diagnostic configuration are used for the original FBPIC implementation, the ported implementation on CUDA, and the ported implementation on the DCU/HIP platform. Figure 5 compares the longitudinal electric field Ez and electron density ne obtained from the three implementations at iteration 1750, corresponding to t = 291.87 fs. As shown in Figure 5(a)-5(c), the three implementations produce nearly identical longitudinal electric-field distributions. Both the position and spatial extent of the wakefield are preserved, including the rapidly oscillating longitudinal field behind the laser pulse and the strong field variation associated with the plasma response. No visible shift in the wakefield phase or change in its spatial morphology is observed between the original FBPIC, ported CUDA, and ported DCU results.

The electron-density distributions in Figure 5(d)-5(f) provide a complementary comparison of the particle dynamics. The laser pulse produces a pronounced electron-depleted region surrounded by a compressed density structure. The location, shape, and extent of this density modulation are consistently reproduced by all three implementations. In particular, the lowdensity cavity and the high-density electron accumulation near its boundary occur at essentially the same longitudinal and radial positions. Since the electron-density distribution is determined by repeated field gathering, particle pushing, and charge and current deposition over many time steps, its close agreement provides a sensitive application-level indication that the rewritten particle kernels preserve the behavior of the original implementation. Figure 6 compares the signed transverse Cartesian electric field E x in the reconstructed x-z plane. The spatial region occupied by the laser pulse and its internal oscillation pattern are consistent across the three implementations. The corresponding E x values of ϵ2 were 5.19 × 10−4 for Ported CUDA and 4.86 × 10−4 for Ported DCU. No backend-dependent phase displacement is apparent at the plotted resolution.

9

Figure 5. Reconstructed x-z cross-sections at y = 0 for the LWFA benchmark at iteration 1750 (t = 291.87 fs). Panels (a-c) show the longitudinal electric field Ez , and panels (d-f) show the electron-density estimator ne = −ρ/e. The columns correspond to Original FBPIC, Ported CUDA, and Ported DCU. A common color scale is used within each row. The density panels use an asinh normalization to display the background plasma and compressed-density structures on the same scale.

Figure 6. Signed transverse Cartesian electric field E x in the reconstructed x-z plane at y = 0 and iteration 1750 (t = 291.87 fs): (a) Original FBPIC, (b) Ported CUDA, and (c) Ported DCU. A common symmetric color scale is used in all panels.

ducted on an NVIDIA Tesla V100 GPU, whereas the DCU experiments were performed on the Sugon 8000 (Dengfeng) supercomputing system equipped with Hygon BW1000 heterogeneous high-performance computing accelerators. The BW1000 is a general-purpose DCU accelerator designed for high-performance computing and heterogeneous workloads and provides 64 GB of HBM2e memory. At the system level, the Sugon 8000 provides a large-scale heterogeneous computing environment interconnected through a scaleFabric highspeed network with native RDMA capability and supported by the ParaStor distributed storage system. This infrastructure provides the inter-node communication and storage environment used in the multi-DCU scalability experiments. The NVIDIA platform used CUDA Toolkit 12.6, including CUDA compiler version 12.6.20. The DCU platform used DTK 25.04.

Table 2. Hardware and software configurations used in the performance evaluation. Item

NVIDIA CUDA platform

Hygon DCU platform

Device

NVIDIA Tesla V100 PCIe 32 GB HBM2 7.0 TFLOPS 900 GB/s CUDA Toolkit 12.6

Hygon BW1000 DCU

Memory Peak FP64 performance Peak memory bandwidth Software stack

64 GB HBM2e 30 TFLOPS 1.8 TB/s DTK 25.04

6. Performance results 6.1. Experimental platform The hardware platforms used in the performance evaluation are summarized in Table 2. The CUDA experiments were con10

DTK (DCU Toolkit) is the software stack provided for Hygon DCU accelerators and includes the HIP-compatible programming environment, compiler toolchain, runtime libraries, and accelerator libraries required for heterogeneous computing. On both platforms, the ported FBPIC backend retains the same Python-side execution model based on CuPy, while the rewritten C/C++ kernels are compiled and executed through the corresponding CUDA or DTK/HIP backend. The NVIDIA V100 platform was used to compare the original FBPIC implementation with the ported CUDA backend, thereby evaluating whether the backend reimplementation affects performance on the original CUDA platform. The Sugon 8000 platform was used for the DCU performance evaluation and multi-accelerator scalability studies described in the following sections. 6.2. Single accelerator performance To evaluate the single-accelerator performance of the proposed backend, we employ a LWFA benchmark, which is also used in the subsequent performance and scalability evaluations unless otherwise specified. LWFA represents a typical plasmaacceleration problem involving the interaction between an intense laser pulse and plasma. As the laser propagates through the plasma, it drives plasma-electron oscillations and generates a wakefield that can be exploited to accelerate charged particles. The simulation involves the major computational components of FBPIC, including electromagnetic field evolution, particle pushing, and particle-mesh interpolation and current deposition, and therefore provides a representative workload for assessing both computational and communication performance in realistic plasma-acceleration applications. For the single-accelerator evaluation,we compare three execution configurations: the original FBPIC implementation on an NVIDIA V100 GPU (Original FBPIC), the ported backend running on the same V100 GPU (Ported CUDA), and the ported backend running on a single DCU accelerator (Ported DCU). The average execution time per simulation step is used as the performance metric. Three problem sizes are considered by varying the longitudinal and radial grid resolutions while keeping the remaining physical and numerical parameters unchanged. The largest case uses Nz = 2452 and Nr = 512. The second case reduces the longitudinal resolution to Nz = 1226 while retaining Nr = 512, and the third case further reduces the radial resolution to Nr = 256. These configurations provide workloads with different computational granularities and therefore allow the influence of the rewritten kernel execution path to be examined over a range of problem sizes. As shown in Figure 7, the ported CUDA backend consistently outperforms the original FBPIC implementation on the NVIDIA V100. For the 2452 × 512 case, the average time per step decreases from 216 ms to 164 ms, corresponding to a speedup of 1.32×. For the 1226 × 512 case, the runtime is reduced from 157 ms to 109 ms, yielding a 1.44× speedup. The largest relative improvement is observed for the 1226 × 256 case, for which the runtime decreases from 80 ms to 52 ms and

Figure 7. Single-accelerator performance of the LWFA benchmark for three grid configurations. The bars show the average time per simulation step for Original FBPIC and Ported CUDA on an NVIDIA V100 GPU and for Ported DCU on a single DCU accelerator. The line shows the V100 speedup of Ported CUDA relative to Original FBPIC, defined as T Original /T Ported CUDA . The ported CUDA backend achieves speedups of 1.32×, 1.44×, and 1.54× for the three problem sizes.

the speedup reaches 1.54×. These reductions correspond to approximately 24.1%, 30.6%, and 35.0% of the original V100 execution time, respectively. The increasing speedup as the problem size decreases suggests that the revised backend improves not only kernel execution but also the fixed overhead associated with the original Numba-based kernel path. Such overhead represents a larger fraction of the time step for smaller workloads, making the benefit of the direct C/C++ kernel implementation and CuPy-based invocation path more visible. The Ported DCU configuration requires 112, 90, and 50 ms per step for the three problem sizes, respectively. These results demonstrate that the same ported backend can execute the LWFA workload efficiently on the target HIP-compatible accelerator without changing the high-level FBPIC simulation workflow. Since the V100 and DCU measurements are obtained on different accelerator architectures, their absolute runtimes are reported as a platform-level comparison rather than as an architecture-normalized speedup. Overall, the results show that extending FBPIC to the DCU platform does not compromise its CUDA performance; instead, the rewritten backend improves execution efficiency on the original NVIDIA platform while providing efficient execution on the target DCU system. 6.3. Strong scaling performance on DCUs To characterize the multi-accelerator scalability of the ported backend, we first perform a strong-scaling experiment using the LWFA benchmark. The global problem size is fixed at Nz = 2453, Nr = 512, and Nm = 3, while the number of MPI ranks is increased from 1 to 8. Each MPI rank is mapped to one DCU. A finite-order stencil with norder = 32 is used to enable longitudinal domain decomposition and guard-cell exchange, while all other physical and numerical parameters are kept unchanged. 11

the inter-rank overhead progressively decreases. An additional change occurs at eight ranks in the evaluated system: the 2- and 4-rank configurations execute within a single compute node, whereas the 8-rank configuration spans two nodes. The sharp degradation at this point therefore suggests an additional contribution from inter-node communication and synchronization. This interpretation is examined quantitatively using the mpiP profiles in Section 6.5. Overall, for the fixed LWFA workload considered here, four DCUs provide the minimum execution time, while further decomposition produces an unfavorable computation-to-communication ratio. 6.4. Weak scaling performance on DCUs The weak-scaling behavior of the ported FBPIC implementation is evaluated using the same LWFA workload. The number of MPI ranks is increased from 1 to 8, with one DCU assigned to each rank. To keep the local longitudinal workload approximately constant, the global grid size Nz is increased proportionally with the number of ranks, taking values of 1200, 2400, 4800, and 9600 for 1, 2, 4, and 8 ranks, respectively. The longitudinal domain extent is increased accordingly, while the grid spacing and the remaining physical and numerical parameters are kept unchanged. Thus, each MPI rank contains approximately 1200 longitudinal grid points throughout the experiment. As in the strong-scaling measurements, 20 warm-up steps are performed before timing, followed by 500 measured steps. The maximum average step time among all MPI ranks is used as the parallel runtime. The weak-scaling efficiency is defined as

Figure 8. Strong-scaling performance of the ported FBPIC implementation for the LWFA benchmark on the DCU platform. The average time per simulation step and the corresponding strong-scaling efficiency are shown as functions of the number of MPI ranks, with one DCU assigned to each rank. The global problem size is fixed at Nz = 2453, Nr = 512, and Nm = 3. The dashed horizontal line denotes the ideal strong-scaling efficiency of 100%. The minimum measured runtime is obtained with four DCUs, whereas the eight-DCU configuration exhibits a pronounced loss of efficiency.

Each configuration is first executed for 20 warm-up time steps to exclude initialization and runtime-compilation overheads, followed by 500 measured time steps. Device-side CuPy events are used for timing. Because the progress of the distributed simulation is determined by the slowest rank, the maximum of the average step times over all MPI ranks is reported as the parallel execution time. The strong-scaling efficiency for N DCUs is defined as T1 , Es (N) = NT N

Eweak (p) =

(4)

T1 , Tp

(5)

where T 1 and T p denote the average time per time step obtained using one and p MPI processes, respectively. The corresponding normalized aggregate throughput is defined as

where T 1 is the average time per step on one DCU and T N is the corresponding time using N DCUs. Figure 8 shows that increasing the number of DCUs initially reduces the time per step, but the improvement is substantially below ideal strong scaling. The runtime decreases from approximately 177 ms on one DCU to 122 ms on two DCUs, corresponding to a speedup of approximately 1.45× and a parallel efficiency of about 72%. With four DCUs, the runtime is further reduced to approximately 94 ms, giving the best measured timeto-solution and a speedup of approximately 1.88×. The corresponding efficiency decreases to approximately 47%, indicating that parallel overhead is already significant even though additional computational resources still reduce the overall execution time. Increasing the configuration from four to eight DCUs reverses this trend. The average time per step rises to 187.7 ms, and the strong-scaling efficiency falls to 11.8%. Thus, the eightDCU configuration is slower than the four-DCU case and provides no time-to-solution advantage over the single-DCU execution for this problem size. Under strong scaling, the longitudinal extent assigned to each rank decreases as the rank count increases, whereas guard-cell exchange and synchronization are still required at every subdomain boundary. Consequently, the amount of useful computation available to amortize

S weak (p) = p

T1 = pEweak (p). Tp

(6)

Figure 9 shows moderate degradation when scaling within a single node. The average time per step increases from approximately 45 ms with one rank to about 60 ms with two ranks, corresponding to a weak-scaling efficiency of approximately 74%. At four ranks, the runtime increases only moderately further, to approximately 65 ms, while the efficiency remains approximately 68%. The corresponding normalized aggregate throughputs are approximately 1.48× and 2.72× for two and four ranks, respectively. These results indicate that the implementation retains useful weak-scaling behavior up to four DCUs despite the additional boundary communication and synchronization introduced by domain decomposition. A substantially larger degradation occurs when the execution is extended to eight ranks. The average time per step increases to 125.3 ms and the weak-scaling efficiency decreases to 35.5%. Although the number of DCUs and the global longitudinal problem size both double from four to eight ranks, the 12

GB. Likewise, the number of MPI_Isend calls is approximately 6300 per interface, and the mean message size remains close to 8.4 MB. This behavior is consistent with FBPIC’s one-dimensional longitudinal decomposition: adding ranks increases the number of subdomain interfaces, while the communication cost associated with each interface remains approximately unchanged. Under strong scaling, the local particle and grid workload decreases with increasing rank count, whereas the per-interface communication requirement does not decrease at the same rate. The resulting reduction in the computationto-communication ratio explains the efficiency loss already observed within a single node. The transition from four to eight ranks introduces an additional penalty because execution changes from intra-node to inter-node communication. The mean MPI time increases from 45.2 s to 94.2 s in the strong-scaling experiment, coinciding with the increase in the average step time from approximately 94 ms to 187.7 ms. This behavior indicates that the internode execution regime introduces substantial additional MPIassociated overhead. Since mpiP measures time spent inside MPI routines, these data do not by themselves separate network transfer latency from synchronization and waiting effects, but they clearly show that the inter-node transition is associated with the observed performance reversal. Figure 10(c) and 10(d) further show that the 8-rank configurations exhibit pronounced rank-level differences in the fraction of time spent in MPI. In the strong-scaling run, the MPI-time fraction ranges from 25.2% to 68.2%, while in the weak-scaling run it ranges from 14.0% to 39.6%. Moreover, the two groups of ranks located on different nodes exhibit systematically different MPI-time fractions. This nonuniformity indicates unequal MPI-associated waiting or progress across ranks and can further amplify synchronization overhead. The available mpiP data, however, are insufficient to distinguish whether the imbalance originates primarily from communication topology, process placement, or differences in computational progress. The weak-scaling results exhibit a related but less direct trend. Because the local longitudinal workload is kept approximately constant, the computation available to amortize communication does not decrease with rank count. Consequently, weak scaling remains relatively efficient within a single node: the efficiencies at two and four ranks are approximately 74% and 68%, with normalized aggregate throughputs of 1.48× and 2.72×, respectively. Nevertheless, the mean MPI time per rank increases from 30.6 s at two ranks to 54.4 s at four ranks and 130.2 s at eight ranks. When execution extends to two nodes, the weak-scaling efficiency falls to 35.5%, while the normalized aggregate throughput increases only from 2.72× to 2.84×. Although the MPI-time fraction decreases slightly from 33.1% at four ranks to 29.3% at eight ranks, the absolute MPI time increases by approximately 2.4×. The lower percentage therefore does not indicate reduced MPI overhead; rather, the total application time grows even more rapidly. These observations are also consistent with the execution characteristics of FBPIC’s domain-decomposed parallelization. Inter-device execution introduces boundary communication and synchronization that are absent in the single-device

Figure 9. Weak-scaling performance of the ported FBPIC implementation for the LWFA benchmark on the DCU platform. The average simulation time per step and weak-scaling efficiency are shown for 1, 2, 4, and 8 MPI ranks, with one DCU per rank. The global longitudinal grid size Nz is increased proportionally from 1200 to 9600 to maintain an approximately constant local workload per rank.

normalized aggregate throughput increases only from approximately 2.72× to 2.84×, an improvement of only about 4.4%. Hence, the additional four accelerators contribute little additional aggregate simulation throughput for this configuration. 6.5. Communication analysis of Multi-DCU scalability To further identify the factors limiting multi-DCU scalability, we profiled the strong- and weak-scaling configurations using mpiP, a lightweight statistical MPI profiling tool[44]. mpiP defines AppTime as the wall-clock time from the end of MPI_Init to the beginning of MPI_Finalize, and MPI time as the wallclock time accumulated inside MPI calls. Profiles were collected at 2, 4, and 8 ranks for both scaling modes, with one rank mapped to each DCU. The 2- and 4-rank jobs ran within one node, whereas the 8-rank jobs used two nodes with four ranks per node. Because the single-rank case involves no interrank communication, the communication analysis begins at two ranks. For strong scaling, the increase in MPI-associated overhead closely correlates with the loss of parallel efficiency. As the rank count increases from 2 to 4 and 8, the mean cumulative MPI time per rank rises from 30.8 s to 45.2 s and 94.2 s (Figure 10(a)), while the fraction of application time spent in MPI increases from 20.2% to 33.7% and 53.4% (Figure 10(b)). Over the same range, the strong-scaling efficiency decreases from approximately 72% to 47% and 11.8%, as shown in Figure 8. Thus, the computational work removed by domain decomposition is progressively offset by communication and synchronization overhead. Table 3 further shows that this degradation is not caused by increasing message size. The aggregate sent volume grows from 53.2 GB at two ranks to 159.1 GB at four ranks and approximately 371 GB at eight ranks, whereas the volume per neighboring interface remains nearly constant at about 53 13

Figure 10. mpiP-based analysis of communication and synchronization overhead in the multi-DCU scaling experiments. (a) Mean cumulative MPI time per rank for the strong- and weak-scaling configurations with 2, 4, and 8 MPI ranks. (b) Fraction of application time spent in MPI routines for the corresponding strongand weak-scaling configurations. (c) Rank-resolved MPI-time fraction for the 8-rank strong-scaling run. (d) Rank-resolved MPI-time fraction for the 8-rank weakscaling run. The 2- and 4-rank configurations execute within a single compute node, whereas the 8-rank configuration spans two nodes with four ranks per node. MPI time denotes the wall-clock time accumulated inside MPI routines and may include communication, synchronization, and waiting overhead. Table 3. Point-to-point communication characteristics from mpiP profiles for the strong- and weak-scaling experiments. The per-interface communication volume is calculated as the aggregate sent volume divided by p − 1, where p − 1 denotes the number of internal interfaces in the longitudinal domain decomposition.

Scaling

p

Sent (GB)

Per interface (GB)

MPI_Isend calls

Mean size (MB)

Strong Strong Strong Weak Weak Weak

2 4 8 2 4 8

53.2 159.1 371.0 53.2 159.1 370.7

53.2 53.0 53.0 53.2 53.0 53.0

6,300 18,900 44,100 6,300 18,900 44,100

8.38 8.40 8.40 8.38 8.37 8.37

7. Conclusion

case. It is therefore most effective when the computational workload assigned to each accelerator is sufficiently large to amortize these costs, or when distributed execution is required because the problem cannot fit within the memory of a single device. In the strong-scaling experiment, increasing the number of DCUs progressively reduces the local workload and therefore moves the simulation toward a regime in which the approximately fixed communication cost per neighboring interface can no longer be effectively amortized. In the weak-scaling experiment, preserving the local problem size prevents this reduction in computational granularity, but communication and synchronization overhead remain exposed and increase substantially once execution extends across nodes.

This paper presented a HIP-compatible accelerator backend for FBPIC that decouples its performance-critical execution path from Numba CUDA while preserving the existing Python interface, device-resident data flow, and MPI-based domain decomposition. The dominant particle and field operations were reimplemented as explicit C/C++ GPU kernels and integrated through CuPy RawKernel, providing a common execution path for CUDA and HIP-compatible platforms. In the linear-wakefield and nonlinear LWFA benchmarks, the ported implementations reproduced the numerical behavior of the original FBPIC code. The relative L2 differences did not exceed 5.19 × 10−4 for the reported field comparisons, and no backenddependent phase shift or structural distortion was observed. The performance results show that the increased portability does not compromise execution efficiency on the evaluated CUDA platform. On an NVIDIA V100 GPU, the ported back14

end achieved speedups of 1.32×, 1.44×, and 1.54× for the three LWFA problem sizes. On the target DCU platform, kernelspecific thread-block tuning reduced the average time per simulation step from 149.171 ms to 112.918 ms, corresponding to a 24.30% reduction and a 1.32× speedup. The multi-DCU experiments achieved their best time-to-solution on four accelerators, with a strong-scaling speedup of 1.88×. Under weak scaling, four DCUs delivered 2.72× the single-DCU aggregate throughput at approximately 68% efficiency. Scaling to eight DCUs produced limited or negative gains. The mpiP profiles attribute this degradation primarily to inter-node communication, synchronization, and the increasingly unfavorable ratio between local computation and the approximately fixed communication cost per subdomain interface. Although the implementation and tuning experiments in this study were conducted on NVIDIA V100 and DCU accelerators, the underlying approach is not restricted to these particular devices. Replacing a vendor-dependent Python JIT backend with explicitly implemented kernels, preserving deviceresident data, and applying compiler- and occupancy-oriented optimization are broadly applicable to other CUDA- and HIPcompatible accelerators. The proposed design therefore provides both a practical accelerator backend for FBPIC and a reference strategy for porting other Python-based scientific applications whose performance-critical components depend on platform-specific JIT frameworks. Future work will focus on reducing multi-node communication overhead, overlapping communication with computation, improving deviceaware MPI data exchange, and introducing architecture- and workload-aware kernel configuration mechanisms to extend efficient FBPIC execution to larger DCU systems.

[2] R. A. Fonseca, L. O. Silva, F. S. Tsung, V. K. Decyk, W. Lu, C. Ren, W. B. Mori, S. Deng, S. Lee, T. Katsouleas, J. C. Adam, OSIRIS: A Three-Dimensional, Fully Relativistic Particle in Cell Code for Modeling Plasma Based Accelerators, in: Computational Science— ICCS 2002, Springer Berlin Heidelberg, 2002, pp. 342– 351. [3] V. K. Decyk, T. V. Singh, Particle-in-Cell algorithms for emerging computer architectures, Computer Physics Communications 185 (2014) 708–719. [4] J.-L. Vay, A. Almgren, J. Bell, L. Ge, D. Grote, M. Hogan, O. Kononenko, R. Lehe, A. Myers, C. Ng, J. Park, R. Ryne, O. Shapoval, M. Thévenet, W. Zhang, WarpX: A new exascale computing platform for beam–plasma simulations, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 909 (2018) 476–479. [5] Z. Liu, M. Hao, W. Zhang, G. Lu, X. Tian, S. Yang, M. Xie, J. Dai, C. Yuan, D. Wang, H. Yang, Optimizing depthwise separable convolution on DCU, CCF Transactions on High Performance Computing 6 (2024) 646–664. [6] Advanced Micro Devices, Inc., HIP documentation, 2026. URL: https://rocm.docs.amd.com/projects/HIP/ en/latest/, accessed 2026-09-03. [7] Y. M. Tsai, T. Cojean, T. Ribizel, H. Anzt, Preparing Ginkgo for AMD GPUs—A Testimonial on Porting CUDA Code to HIP, in: Euro-Par 2020: Parallel Processing Workshops, Springer International Publishing, 2021, pp. 109–121.

Acknowledgements

[8] R. Lehe, M. Kirchen, I. A. Andriyash, B. B. Godfrey, J.L. Vay, A spectral, quasi-cylindrical and dispersion-free Particle-In-Cell algorithm, Computer Physics Communications 203 (2016) 66–82.

This work is supported by the National Key Research and Development Program of China (2024YFB4504103), the Fund of Laboratory for Advanced Computing and Intelligence Engineering (2025-ZZKY-016), and the National Natural Science Foundation of China (12574380). This work is also supported by Jiangsu Province Engineering Research Center of IntelliSense Technology and System.

[9] B. B. Godfrey, J.-L. Vay, I. Haber, Numerical stability analysis of the pseudo-spectral analytical time-domain PIC algorithm, Journal of Computational Physics 258 (2014) 689–704. [10] FBPIC contributors, FBPIC 0.27.0 documentation: Installation on a local computer, 2026. URL: https://fbpic. github.io/install/install_local.html, accessed 2026-09-03.

Data and Code Availability The source code developed in this study is available in the Mendeley Data repository at https://doi.org/10.17632/ r66r2vzcjc.1. The repository is currently under embargo and will become publicly available after the embargo period.

[11] S. K. Lam, A. Pitrou, S. Seibert, Numba: a LLVMbased Python JIT compiler, in: Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, ACM, 2015, pp. 1–6.

References

[12] R. Okuta, Y. Unno, D. Nishino, S. Hido, C. Loomis, CuPy: A NumPy-Compatible Library for NVIDIA GPU Calculations, in: Proceedings of Workshop on Machine Learning Systems (LearningSys) in the Thirty-first Annual Conference on Neural Information Processing Systems (NIPS), 2017.

[1] C. K. Birdsall, A. B. Langdon, Plasma Physics via Computer Simulation, Institute of Physics Publishing, Bristol and Philadelphia, 1991.

15

[13] Numba Development Team, Deprecation Notices: ROCm target, Numba 0.54.1 documentation, 2021. URL: https://numba.readthedocs.io/en/0.54.1/ reference/deprecation.html, accessed 2026-07-20.

J. Dargent, C. Riconda, M. Grech, Smilei: A collaborative, open-source, multi-purpose particle-in-cell code for plasma simulation, Computer Physics Communications 222 (2018) 351–373.

[14] CuPy Development Team, CuPy 14.2.0 documentation: cupy.rawkernel, 2026. URL: https: //docs.cupy.dev/en/v14.2.0/reference/ generated/cupy.RawKernel.html, accessed 202609-03.

[23] S. Diederichs, C. Benedetti, A. Huebl, R. Lehe, A. Myers, A. Sinn, J.-L. Vay, W. Zhang, M. Thévenet, HiPACE++: A portable, 3d quasi-static particle-in-cell code, Computer Physics Communications 278 (2022) 108421.

[15] CuPy Development Team, CuPy 14.2.0 documentation: Using CuPy on AMD GPU (ROCm), 2026. URL: https: //docs.cupy.dev/en/stable/install.html, accessed 2026-09-03.

[24] W. Zhang, A. Myers, K. Gott, A. Almgren, J. Bell, AMReX: Block-structured adaptive mesh refinement for multiphysics applications, The International Journal of High Performance Computing Applications 35 (2021) 508– 526.

[16] B. B. Godfrey, Numerical Cherenkov instabilities in electromagnetic particle codes, Journal of Computational Physics 15 (1974) 504–521.

[25] FBPIC contributors, FBPIC 0.27.0 documentation, 2026. URL: https://fbpic.github.io/, accessed 2026-0903.

[17] J.-L. Vay, C. Geddes, E. Cormier-Michel, D. Grote, Numerical methods for instability mitigation in the modeling of laser wakefield accelerators in a Lorentz-boosted frame, Journal of Computational Physics 230 (2011) 5908–5929.

[26] A. Sinn, [DO NOT MERGE] test on AMD GPUs, fbpic/fbpic Pull Request #740, 2025. URL: https://github.com/fbpic/fbpic/pull/740, open, unmerged experimental pull request; head commit e3e299bc3d5de76c3bb8781b8eb70665206b149a; accessed 2026-09-03.

[18] A. Myers, A. Almgren, L. D. Amorim, J. Bell, L. Fedeli, L. Ge, K. Gott, D. P. Grote, M. Hogan, A. Huebl, R. Jambunathan, R. Lehe, C. Ng, M. Rowan, O. Shapoval, M. Thévenet, J.-L. Vay, H. Vincenti, E. Yang, N. Zaïm, W. Zhang, Y. Zhao, E. Zoni, Porting WarpX to GPUaccelerated platforms, Parallel Computing 108 (2021) 102833.

[27] FBPIC contributors, Parallelization of FBPIC, 2026. URL: https://fbpic.github.io/overview/ parallelisation.html, FBPIC 0.27.0 documentation; accessed 2026-09-03. [28] S. Jalas, I. Dornmair, R. Lehe, H. Vincenti, J.-L. Vay, M. Kirchen, A. R. Maier, Accurate modeling of plasma acceleration with arbitrary order pseudo-spectral particlein-cell methods, Physics of Plasmas 24 (2017) 033115.

[19] L. Fedeli, A. Huebl, F. Boillod-Cerneux, T. Clark, K. Gott, C. Hillairet, S. Jaure, A. Leblanc, R. Lehe, A. Myers, C. Piechurski, M. Sato, N. Zaïm, W. Zhang, J.-L. Vay, H. Vincenti, Pushing the frontier in the design of laserbased electron accelerators with groundbreaking meshrefined particle-in-cell simulations on exascale-class supercomputers, in: SC22: International Conference for High Performance Computing, Networking, Storage and Analysis, IEEE, 2022, pp. 25–36.

[29] NVIDIA Corporation, CUDA C++ Programming Guide, 2026. URL: https://docs.nvidia.com/cuda/ cuda-programming-guide/, accessed 2026-09-03. [30] Advanced Micro Devices, Inc., HIP C++ Language Extensions, 2026. URL: https://rocm.docs.amd. com/projects/HIP/en/develop/how-to/hip_cpp_ language_extensions.html, accessed 2026-09-03.

[20] H. Burau, R. Widera, W. Hönig, G. Juckeland, A. Debus, T. Kluge, U. Schramm, T. E. Cowan, R. Sauerbrey, M. Bussmann, PIConGPU: A fully relativistic particle-incell code for a GPU cluster, IEEE Transactions on Plasma Science 38 (2010) 2831–2839.

[31] S. Ryoo, C. I. Rodrigues, S. S. Stone, J. A. Stratton, S.-Z. Ueng, S. S. Baghsorkhi, W.-m. W. Hwu, Program optimization carving for GPU computing, Journal of Parallel and Distributed Computing 68 (2008) 1389–1401. [32] FBPIC contributors, FBPIC GPU utility source: fbpic/utils/cuda.py, GitHub source file at commit 2c4b9ca4c20847672fb7a86d61df3aeb766a9e8d, 2026. URL: https://github.com/fbpic/fbpic/blob/ 2c4b9ca4c20847672fb7a86d61df3aeb766a9e8d/ fbpic/utils/cuda.py, accessed 2026-07-20.

[21] E. Zenker, R. Widera, A. Huebl, G. Juckeland, A. Knüpfer, W. E. Nagel, M. Bussmann, Performanceportable many-core plasma simulations: Porting PIConGPU to OpenPower and beyond, in: High Performance Computing, volume 9945 of Lecture Notes in Computer Science, Springer International Publishing, 2016, pp. 293–301.

[33] S. Hong, H. Kim, An analytical model for a GPU architecture with memory-level and thread-level parallelism awareness, ACM SIGARCH Computer Architecture News 37 (2009) 152–163.

[22] J. Derouillat, A. Beck, F. Pérez, T. Vinci, M. Chiaramello, A. Grassi, M. Flé, G. Bouchard, I. Plotnikov, N. Aunai, 16

[34] Advanced Micro Devices, Inc., HIP Performance Guidelines, 2026. URL: https://rocm.docs.amd.com/ projects/HIP/en/latest/how-to/performance_ guidelines.html, accessed 2026-09-03.

[44] J. S. Vetter, C. Chambreau, mpiP: A lightweight mpi profiler, 2020. URL: https://software.llnl.gov/ mpiP/, mpiP 3.5 User Guide; accessed 2026-09-04.

[35] A. Fanfarillo, N. Curtis, Register pressure in AMD CDNA2 gpus, 2023. URL: https: //gpuopen.com/learn/amd-lab-notes/ amd-lab-notes-register-pressure-readme/, originally published 2023-05-17; updated 2024-06-26; accessed 2026-09-03. [36] W. Hu, L. Han, P. Han, J. Shang, Automatic Thread Block Size Selection Strategy in GPU Parallel Code Generation, in: Parallel Architectures, Algorithms and Programming, Springer Singapore, 2021, pp. 390–404. [37] C. M. Pereira, A. L. Pinheiro, R. Schirru, Automatic block dimensioning on GPU-accelerated programs through particle swarm optimization, Information and Software Technology 123 (2020) 106299. [38] M. Lurati, S. Heldens, A. Sclocco, B. van Werkhoven, Bringing Auto-Tuning to HIP: Analysis of Tuning Impact and Difficulty on AMD and Nvidia GPUs, in: Euro-Par 2024: Parallel Processing, Springer Nature Switzerland, 2024, pp. 91–106. [39] G. Stantchev, W. Dorland, N. Gumerov, Fast parallel Particle-To-Grid interpolation for plasma PIC simulations on the GPU, Journal of Parallel and Distributed Computing 68 (2008) 1339–1349. [40] X. Kong, M. C. Huang, C. Ren, V. K. Decyk, Particlein-cell simulations with charge-conserving current deposition on graphic processing units, Journal of Computational Physics 230 (2011) 1676–1685. [41] FBPIC contributors, FBPIC linear-wakefield verification test: tests/test_linear_wakefield.py, 2026. URL: https://github.com/fbpic/fbpic/blob/ 2c4b9ca4c20847672fb7a86d61df3aeb766a9e8d/ tests/test_linear_wakefield. py, gitHub source file at commit 2c4b9ca4c20847672fb7a86d61df3aeb766a9e8d; accessed 2026-09-04. [42] E. Esarey, C. B. Schroeder, W. P. Leemans, Physics of laser-driven plasma-based electron accelerators, Reviews of Modern Physics 81 (2009) 1229–1285. [43] FBPIC contributors, FBPIC laser-wakefield acceleration example: lwfa_script.py, 2026. URL: https://github.com/fbpic/fbpic/blob/ 2c4b9ca4c20847672fb7a86d61df3aeb766a9e8d/ tests/test_linear_wakefield. py, gitHub source file at commit 2c4b9ca4c20847672fb7a86d61df3aeb766a9e8d; accessed 2026-09-04. 17