Numerical simulation of charged particle dynamics is required for the development of electron beam ion sources (EBIS). Special attention should be paid to the formation and evolution of different instabilities that can manifest in such sources. Understanding these processes will allow the efficiency of ion sources to be improved. The simulation of two-stream instability emerging between slow ions and fast electrons is particularly important. To perform the required simulation, “ef” and “ef_python” applications are under development. In these applications simulation is performed using the particlein-cell (PIC) method, which contains four main parts: a particle mover, a particle-to-grid scatter, a field solver, field-to-particle interpolation. Each part's performance depends on the number of particles, as well as on the granularity of the simulation spatial mesh. The field solver usually represents a major portion of computational difficulty. This report describes the efforts to improve the performance of the ef_python application using CuPy, AMGX and PyAMG libraries. With minimal changes in the source code, CuPy allowed us to port all simulation computations to CUDA. In addition, to accelerate the field solver, algebraic multigrid methods, implemented by the PyAMG library on CPU and the AMGX library on GPU, were utilized. Major speed-up was achieved, especially when running simulations with very high-resolution spatial grids on high-performance GPUs.
Numerical simulation of charged particle dynamics is required for the development of electron
beam ion sources (EBIS) [
To perform the required simulation, “ef” and “ef_python” applications are under development
[
A previous study using a 2-dimensional dynamic particle model as a sample showed that both
OpenCL and CUDA could dramatically improve particle simulation performance compared to
NumPy [
The idea of the particle-in-cell method is to simulate individual positions and velocities of particles over small time steps, while using a grid to approximate collective parameters, such as charge density and current density created by particles. The simulation domain is defined by a cuboid regular mesh and a time step. As part of the problem definition, Ef and Ef_pyhton allow defining conducting volumes with a fixed potential within the simulation domain. Only simple volume shapes, such as a cylinder or a sphere, are supported right now. The same shapes define particle sources that generate new particles with a normal velocity distribution within their volumes at each time step. It is possible to add externally defined electric and magnetic fields.
Each simulation step within Ef_python consists of the following operations: advance particle positions and momenta, generate new particles, calculate charge density, compute electric potential, calculate electric field.
To compute the electric field acting on each particle, linear interpolation of the electric field from the 8 surrounding grid nodes is used. This is not the only possible interpolation kernel. Grid-toparticle and particle-to-grid kernels are an important area of possible algorithmic optimization to improve simulation accuracy and performance in this project.
A “particle mover” denotes the simulation component responsible for updating particle
positions and velocities at each time step. The choice of the particle mover is important for assuring
conservation of important invariants, such as energy. As the particle mover in ef_python, the
wellestablished leapfrog second-order explicit method, known as the Boris scheme, is used [
Particles that leave the simulation volume or collide with inner conducting volumes are absorbed and removed from simulation. New particles generated by particle sources are added to the system and their momentum is simulated half a time step back as it is required by the “leapfrog” particle mover.
The charge of each particle is linearly divided into the 8 surrounding grid nodes. This firstorder scheme is sometimes called cloud-in-cell. The collective current is not used in ef or ef_python right now, as the magnetic field created by particles is considered negligible.
To simulate the electromagnetic field created by particles, Maxwell’s equations on the grid can be approximately solved by many algorithms based on 3 main approaches: finite difference methods (FDM), finite element methods (FEM) and spectral methods. In ef and ef_python, the FDM method is used to solve the Poisson’s equation and compute the electric potential created by the simulated particles, while considering the effect of conducting volumes within the simulation domain. The electric field is then simply a numerically computed gradient of the potential on the grid.
The contribution of each step of simulation to the total execution time depends on different simulation parameters and implementation details. Profiling tools allowed us to find bottlenecks and improve performance using different simulation setups. The field solver is usually the major component of simulation in terms of resource usage.
As a result of profiling, several unreasonably slow function implementations were fixed, including rewriting the FDM equation sparse matrix generator using diagonal SciPy matrix classes, caching which mesh nodes are inside conducting volumes, delegating field-to-particle interpolation to SciPy, and rewriting particle-to-field weighting for better NumPy parallelization.
In order to safely conduct the major code restructuring and reimplementation required for this project, over 200 unit tests were created, providing total code coverage of 91%.
The point of major computation complexity of ef_python is the field solver. As such, it was the first component picked for acceleration. The SciPy conjugate gradient solver for sparse matrices was used in ef_python before these improvements.
Multigrid methods are methods of numerically
solving differential equations on a hierarchy of grids
with increasing coarseness, allowing one to faster reduce
the large-wavelength error [
CuPy is a library aiming to help easy acceleration of the NumPy-based code on GPU by providing in-place replacements for most NumPy functions and operators. By using CuPy we accelerated most of the simulation operations in ef_python with minor code changes.
Grid interpolation was not available in CuPy and had to be rewritten as a custom CUDA kernel. The gradient function from Numpy that was used to compute the electric field from the potential was not implemented Figure 2 Comparison of the simulation time for the in CuPy. It had to be reimplemented for the case of regular grid steps using CuPy two-stream instability problem with higher precision operators. of the field solver and a smaller number of particles
Simulation classes were rewritten to use class injection for the convenient runtime NumPy/CuPy and PyAMG/AMGX selection. This allowed using the selected simulation backend with minimal code duplication. Only the field solver and spatial mesh classes required explicit implementation-dependent subclasses.
are shown in Figure 2 and Figure 3.
In addition to unit tests, sample simulation was regularly performed to detect errors during development. It was the simulation of a radially symmetric charged particle beam starting from a cylindrical emitter. The shape of the beam caused by electrostatic particle repulsion can be predicted analytically. The comparison between the theoretical shape and the numerical solution is shown in Figure 1.
The simulation of two-stream instability was performed with ef using OpenMP with 1 and 12 threads to provide a baseline of comparison. The performance was compared to ef_python using CuPy and PyAMGX running on Nvidia Tesla K40 and Nvidia Tesla K80 graphics cards. The results of this comparison
With minimal changes to the program, CuPy allowed accelerating most of the simulation operations on GPU through CUDA. In addition, to improve the field solver, algebraic multigrid methods, implemented by the PyAMG library on CPU and the AMGX library on GPU, were utilized. Major speed-up was achieved, especially when running simulations with very high-resolution spatial grids on high-performance GPUs.
Using multiple GPUs is an important future prospect for the ef_python application. This will require using MPI for process coordination and data exchange. AMGX already has MPI-based multiGPU support, but the PyAMGX interface library will have to be extended to add it in. Another promising avenue of future research is using OpenCL as an acceleration backend, as it can utilize nonNvidia GPUs.
Profiling and testing were partly performed on the HybriLIT heterogeneous computing
platform [
The reported study was funded by RFBR according to the research project № 18-32-00239.