Parallel Particle-In-Cell Plasma Simulations
on GPU, Vector, and Scalar System
Tooru Sugiyama1
11

Earth Simulator Center, Japan Agency for Marine-Earth Science and Technology (JAMSTEC), Japan

Introduction

Particle simulation model is quite useful tool to investigate particle kinetics in plasma physics. The motions ofthe charged particles are self-consistently calculated in theelectromagnetic field [1]. The simulation step in an explicit Particle-in-cell (PIC) model consists of three parts;
particle acceleration and motion where the Lorentz force
is evolved on the particle’s position by interpolating theelectromagnetic field defined on grids; field solver whereMaxwell’s equations are solved on grids; and velocity moment calculation where the particle flux and density arecalculated from the particle velocity and position. The
number of the super-particles is huge, for example, overhundred particles are allocated in each simulation cell.
Resultantly, the calculation cost of the PIC simulation isquite higher than that of the Magneto Hydro Dynamics(MHD) model. The highest cost in the particle simulationis the velocity moment calculation. This is because theparticles are randomly distributed in the simulation system and they are related with tidily located grids. Therefore, this procedure also requires random access to thephysical memory on the computational system. Then, thecache on the CPU chip does not work well in this procedure. However, the slow data transfer band width between
CPU and DDR3 memory ( 36 GB/sec. ) doesn’t become
the bottleneck of the calculation speed. Because the number of the moment calculation procedure is quite larger
than that of the data transport. This coarse-grained modelis appropriate to the calculation on GPU and Vector-typesystem.
The GPU chip on nVIDIA C2070 board has 448 CUDA
cores, and then we need to use a parallelized algorithm tomake a good performance on the chip. The vector typesystem also has a large vector register on the chip, where
we need also to use a parallelized algorithm to make agood performance. Here we use the particle and domain
decomposition methods for parallel calculation. Since theparticles are randomly distributed in the simulation system as noted above, the moment calculation is essentiallythe atomic calculation. To avoid this memory conflict,some algorithms are introduced [3] [4] [5]. Here we haveapplied these algorithms and implement them on GPU,
Vector-type, and Scalar-type systems. We have measured
the calculation times.
2Simulation Model

We treat an ion beam instability using a threedimensional hybrid simulation model where the ions aretreated as macro-particles while the electrons are treatedas a charge-neutralizing massless fluid [2].

Fig. 1. Poster of JSST 2012.

In the moment calculation, we count how many particles
are located in a cell with the particle shape function (Fig.
1), which is a similar procedure as the counting histogram.
Note that the figure shows the first order accuracy model in
two-dimensional system for convenience. When we countthe number of particles using multi-processors, cores, or on
the vector registers, (hereafter thread as GPU language)
with the particle decomposition method, the memory conflict takes place, because each thread treats each particle.
Suppose that some particles located in a certain cell andthese particles are counted by different threads. Thesethreads try to modify the memory for counting in the cell
if these particles locate in the same cell. For the safety
counting, the read-modify-write procedures must be operated without interruption by other threads. If not, thecounting results are not correct. That is, while one processor counts a particle moment, other processors should wait
until the counting finishes. This waiting reduces the parallel performance on multi-thread system, especially on thehighly parallelized system. To avoid the conflict, we apply some algorithms, those are, SORT, RETRY and LIST
algorithm.
3

Performance Measurements

The detailed results for the performance measurements
are reported both on the followed full paper and the conference presentation.
4Discussion

We have introduced SORT, RETRY and LIST algorithms in order to perform plasma PIC simulation on
multi-core and vector-type systems. The most complicated procedure in PIC simulation on highly parallelized
platform is the moment calculation. Here, we apply thesemethods for the moment calculation.

The most useful point using the SORT method is thatthe cache works well on both calculation of the moment
and velocity advances. This is because the same field
and moment values are used for the sorted particles. On
the other hands, the most useful point using RETRY andLIST methods is that the particle sorting procedures areexcluded. One of the reasons to perform the PIC simulation is that particle trajectories are followed in the selfconsistent electromagnetic-fields in order to investigate theparticle origin of the attracted particles. By sorting theparticle ID number with the particle position in the simulation system, it makes complicated to follow the particle
trajectories.
Here, we have not arranged the simulation programs todeal with the general number of system size. Indeed, thesystem size is strongly related with the number of threads
in one thread-block on GPU as 1024, and also the numberof the vector-register as 256. Therefore, there would exist
the optimized number for the system size (e.g. the number of grids, the number of particles per cell) in order tomaximize the performance.
In addition, we have realized one challenging problem on
multi-core systems. The sorted particle ID-list created bythe multi-core system is different in each run, even though
the same program are used. It certainly leads the roundoff error in calculating the particle moments in a different
order of the particle ID number. If the results are sensitive
to the round-off error, we never have the same results. Now
we don’t have any solutions about the error, but the ion
beam instability tested here is acceptable for the error.
References[1]

C. K. Birdsall, & A. B. Langdon, 1991, Institute of Physics
Publishing, Bristol and Philadelphia.

[2]

Fujimoto, M., 1992, PhD thesis, Institute of Space andAstronautical Science, Japan.

[3]Imamura, T., 2005, IPSJ journal, 46, SIG 7 (ACS 10),52-62.

[4]

Sugiyama, T., et al. 2004, IPSJ journal, 45, SIG 6 (ACS6), 171-175.

[5]Sugiyama, T., 2011, JAMSTEC Rep. Res. Dev., 12, 13-25.