Low energy ion and electron beams, produced by ion sources and electron guns, find their use in surface modification, nuclear medicine and injection into high-energy accelerators. Simulation of particle dynamics is a necessary step for optimization of beam parameters. Since such simulations require significant computational resources, parallelization is highly desirable to be able to accomplish them in a reasonable amount of time. From the implementation standpoint, dynamically typed interpreted languages, such as Python 3, allow high development speed that comes at cost of performance. It is tempting to transfer all computationally heavy tasks on a GPU to alleviate this drawback. Using the example of a charged particles dynamics simulation problem, various GPUparallelization technologies available in Python 3 are compared in terms of ease of use and computational speed. The reported study was funded by RFBR according to the research project № 18-32-00239\18. Computations were in part held on the basis of the heterogeneous computing cluster HybriLIT (LIT, JINR)
The problem of charged particle dynamics simulation has been a challenge for scientists and
engineers in many fields, such as plasma physics, accelerator engineering, beam transport, ion sources,
high current facilities etc. Simulation allows to predict conditions and results of experiments,
facilitating experiment facilities design. In the case of an ion source, it is necessary to simulate an
order of 1012 particles [
A program utilizing a GPU for computing has components running on the CPU and GPU. The
GPU component performs the compute-intensive and parallelizable operations, and the CPU
component performs all the other tasks, such as managing the computation, configuring the GPU and
data I/O. It is necessary to explicitly manage data transfer between the CPU and GPU components of
the program. To manage the communication between the GPU and the program there are proprietary
and open APIs: Nvidia CUDA[
We tested 3 popular libraries for GPU computing available for Python 3: Numba[
We chose a simple problem as our example: simulating the classical dynamics of a set of N identical charged particles in a 2-dimensional volume of space. Namely, each particle is represented as its 2 cartesian coordinates, 2 velocity components, charge and mass. Dynamics equation are dx = v ∙ dt, mdv= F ∙ dt The integration was done using the basic Euler method with a constant time step.
xt+1 = xt + vt ∙ h, vt+1 = vt + at ∙ h The force acting on each particle simply as sum of Coulomb forces from the other particles.
qiqj⃗r⃗⃗i⃗j ⃗F⃗⃗⃗i = ∑
rij3 i≠j
Coulomb’s constant is 1 here, as we use arbitrary units for simulation. The particles are contained in a two-dimensional square volume. On reaching the boundary, a particle gets perfectly elastically reflected. Initially the particle positions are evenly non-randomly distributed on a spatial grid, while the velocity components are uniformly distributed in the interval (-Vmax, Vmax). Observing the coordinate and velocity distribution of the particles over time allows to debug the simulation implementations. Proceedings of the VIII International Conference "Distributed Computing and Grid-technologies in Science and
Education" (GRID 2018), Dubna, Moscow region, Russia, September 10 - 14, 2018
To compare different GPU acceleration libraries for simulation a one-step function was implemented with each library using analogous code. Then the performance of simulation for each library was measured depending on the scale of simulation. At each scale the number of particles and volume dimensions were set, then the initial particle positions and velocities were generated. The data was copied to GPU memory, and one batch of simulation steps was performed and timed, each step corresponding to one step of numerical integration. The batch size was varied with scale, to keep the batch processing time at a consistent level. The batch needed to be long enough to allow measuring performance accurately, but short enough to measure and plot at many scales.
For each batch, the total simulation time was measured using the standard python time.perf_counter() function. This was then converted into steps per second, particle updates per second and particle interactions per second, which are more comparable across scales.
The development and performance testing were carried out using Google Colaboratory on a Nvidia Tesla K80 GPU, Driver Version 396.44. Library versions are given in Table 1. In this test both Numba and OpenCL used CUDA as the back-end interface to the GPU. Numpy runs on the CPU and provides the baseline performance to compare GPU speed-up.
The plots of performance results are shown in Figure 2. At small scales (n < 102) GPU steps per second are constant, while CPU stays ahead. This suggests that start-up overhead dominates perparticle computation. At scales of 102 < n < 104 it can be observed that particle updates per second is near-constant, that is to say, computation time grows as O(n). As each particle update consists of O(n) force computations to derive acceleration, we can conclude that at these scales all update streams run perfectly in parallel. As a check, one can compare the update thread count (n) to the number of threads the GPU can execute in parallel. This can be derived from the GPU multiprocessor count and max threads per multiprocessor. For the Tesla K80, this is 13 × 2048 = 26624 which agrees with the critical point that can be seen at ~2 ∙ 104. At large scales of > 104 the GPU is fully utilized and particle interactions per second tend to a constant value, signifying O(n2) performance.
We compared several GPU parallelization technologies available for Python 3. All provide considerable (~103) speed-up compared to CPU computing using Numpy. OpenCL seems to outperform Numpy and CUDA. However, this test did not check all factors that may affect GPU performance. Among those outside of scope of this paper are: data type and shape of arrays being processed, thread block size and shape, using scalar or constant parameters for mass and charge, library version, driver version, GPU model, optimized block caching, optimizing the algorithm.
The reported study was funded by RFBR according to the research project № 18-32-00239\18. Computations were held partially on the basis of the heterogeneous computing cluster HybriLIT (LIT, JINR).