One of such particle models is the particle-in-cell model. This model derives its name from a special technique used to calculate the electric field in a self-consistent manner by solving the Poisson equation on a finite grid. This tech-nique, reducing the computational complexity from O(n2) to O(nlogn) will be explained in subsection 4.2.2. As self-consistent modeling started to emerge in the second half of 20th century [86–88], it became a routine scientific technique
1As long as the particle mass does not change, i. e. no relativistic effects or inelastic collisions are accounted for.
4.2. Particle-in-cell model with the increase of computational power, enabling us to simulate large scale phenomena [89, 90].
In this part of the thesis, the basic principles of the model will be explained as well as some approaches that are common for not only this model type. The magnetic field in this description will be considered to have only the time-invariant external component Bext(r).
4.2.1 Plasma characterization
Basic quantities needed to assess the physical quantities simulated by a particle model are particle positions and velocities. The particle positions are tracked in either two (ri = (yi, zi)) or three dimensions (ri = (xi, yi, zi)), while particle velocities are recorded usually in all three dimensions to describe the Larmor rotation properly, vi = (vx, i, vy, i, vz, i). Particles can be of arbitrary number of species, each carrying the electric chargeqi and having mass ofmi.2 Additionally, supporting quantities are useful for effective model run, such as local charge densityρ(r)and potential ϕ(r) needed for the calculation the electric fieldE(r). The important feature of the model is the simulation grid. Since charge density and several other quantities will not be treated as continuous, a discretization is to be applied. The simplest type of the discretization of the simulation space is to divide it by a rectangular (or cuboid) grid. An example of this approach in 2D is in the fig. 4.1 (left panel). This brings a transition of the continuous ρ(r) to discrete values of ρ(ix, iy, iz), where i[] stand for grid edge boundaries with inter-edge distance d[]. Since this discretization exists principally on account of field calculation, the spacing is defined by the Debye length of simulated plasma as λD is the stability limit for PIC algorithm [91]. As indicated by the fig. 2.2, δx/λD = δy/λD = δz/λD ≡ δ/λD = 1 is usually sufficient to obtain smooth potential resolution. To properly describe potential structures with L ≈ λD, lower values are necessary, makingd = 1 the upper bound of grid spacing.
Suitable boundary conditions have to be chosen both with respect to the par-ticles and with respect to fields, however, this depends notably on the simulated phenomenon and/or geometry.
The grid presented in this example is rectangular. However, this is not an inherent feature of PIC models. Cylindrical grids are advantageous for simulation of axisymmetric problems, triangular grids provide more accurate description of objects surfaces inside the plasma volume, while making the field calculation more complex.
At this moment, it is useful to elaborate more on the representation of par-ticles. The number of particles in a cubic centimetre of a plasma with n0 = 1×1018m−3 is 1012. Such a number would mean that the computer performing the simulation would use roughly 48 TB of operational memory just to store the information about positions and velocities (in 3D), while the plasma density is still quite low for typical conditions in the SOL of tokamaks. This is obviously
2It is to be noted that indices in italics, such asi in mi refer to i-th particle species. The particular species such as ions and electrons are labelled by upright indices (e. g. mi,me).
Figure 4.1: Left: The grid in the two-dimensional PIC simulation. Particles have arbitrary position in the space, however, the charge density ρ is divided into a finite grid. The grid also contains projections of object immersed into the plasma.
Right: The contribution of single particle electric chargeqito the total density is divided into neighbouring nodes A–D via weighting specified by gray areas WA–WD (note that the weights are drawn in opposite corners). 3D situation is similar, with the difference that the rectangular grid is replaced with the cuboid one.
out of capabilities of contemporary technology, but fortunately, it is not necessary to simulate all the particles in the plasma. The key is to merge the particles into a reasonably small number of macroparticles in a grid cell represented by a Debye cube (a cube with side equal toλD), while keeping the key plasma characteristics.
In the aforementioned plasma with Te= 10 eV, the number of particles in a such cube is ND ≈ 1.3×104. There are several conditions implied on the minimal value of ND sufficient for stability of the simulations [92, p. 5], which is around ND = 10. Additional constraint is the possibility of artificial diagnostics of the system, since the quantities are evaluated with uncertainty ∼√
N /N, which falls below 10 % for N = 100. This aspect can be however treated by time averaging.
The mass mmp and chargeqmp of this substitute is a sum of masses and charges of individual particles from which is the macroparticle formed.
The first part of the model’s namesake is to assign the charge density to grid points. To do so, the technique named cloud-in-cell [91] is applied. The macroparticle itself is not seen a charged point mass, but rather as a cloud with the dimension of one cell. In this way there is a straightforward division of the charge between corner points of the cell in which is the particle contained at specific time as can be seen in the fig. 4.1. The particle position divides the grid into 4 rectangles (in 2D; 6 cuboids in 3D) each of them carrying a part of the charge qW defined by its area W; qW = qmpW/(δyδz). These weighted contributions are calculated for all macroparticles in the system and assigned to respective corner points forming the charge density grid ρ(ix, iy, iz).
4.2. Particle-in-cell model
4.2.2 Electric field
The technique devised to calculate the electric field that was coined by Evans and Harlow in 1957 [86]. In order to do so, the computationally demanding procedure ofE calculation via the Coulomb law (4.2) is replaced by a much faster calculation of the electric potential ϕ via Poisson equation (2.11). The field problem is now reduced to calculation of the charge density ρ(r), because E(r) can be easily calculated as
E(r) =−∇ϕ(r) . (4.3) Having introduced the transition betweenρ(r)andρ(ix, iy, iz)in the previous subsection, the information needed for calculation of particle movement is now complete. Last remaining step is to solve the Poisson equation.
While to do so is faster than looping through all available pairs of particles, it is still a demanding task. Various methods are employed, usually based on finite differences, elements or volumes method, ranging from slow Gauss relaxation through more optimized iterative solvers to direct solvers or spectral methods.
The selection of suitable solver depends on both grid size and situation geometry as well as the dimensionality of the simulation – large 3D grids for example forbid usage of direct solution due to extreme memory demands.
4.2.3 Particle motion
Particle motion itself is provided by the integration of the equation of motion.
For PIC, the leap-frog [93] method is sufficient, while keeping minimal memory equations and having the numerical stability for a reasonably large time step
∆t. The transition from continuous to discrete values is provided by the finite differences method:
midvi
dt =Fi −→ mivi,new−vi,old
∆t =Fi,old, dri
dt =vi −→ ri,new−ri,old
∆t =vi,new.
(4.4)
In the leap-frog method, the velocity and location of i-th particle are tracked with a time shift of ∆t/2, i. e. that when the position rnew is being calculated, the velocity vnew from previous half-step is used.
To perform calculations using the finite differences method, some value of the time step has to be chosen. There are two particular time intervals that limit the time step length. They are
ttrans ≈ dλD
⟨ve⟩ =q
√︃ε0me
3e2n0, the time particle needs to cross the cell, and (4.5) tgyro ≈ 2
ωc, e = 2me|e|
B , the limit of iteration per one gyration orbit, (4.6) where the gyration limit time tgyro is the stability criterion time [92, p. 56] for leap-frog algorithm. The cell size is usually taken as a fraction of Debye length
(dλD), since its magnitude depends on the size of potential structures expected in the system.
In simple PIC models with no wave representation, these limits are sufficient for a smooth calculation as well as stable numeric integration. Each particle species has its own critical time, but for the model operation, the time step defined by the lightest particle species defines the global limit. It is clear that these critical times are given mostly by electron mass, since other parameters are either fundamental constants or essential plasma characteristics. It is common for PIC models to use higher value of me than it is in reality, since for most studies, it is not necessary to simulate the minuscule finite Larmor effects of electrons.
The ratio mi/me ≈ 200 [40, 44, 94], instead of the real value mi/me = 3670, is usually a fair compromise.
4.2.4 Boundary conditions and particle injection
Due to a finite volume simulated by the model, it is necessary to have a strat-egy describing the process of handling the situation when a particle crosses the boundary (located at r =rb ∈∂V) of the system as well as the boundary condi-tions applied on the potential in the simulation volume. Multiple boundaries are usually present inside the simulated volume – apart from the boundaries of the system itself, the object immersed into the plasma have boundaries on their own.
Boundary conditions on the potential are derived from conditions usually applied on solution of partial differential equations and are usually implemented into Poisson equation solver. Most common conditions are:
Dirichlet ϕ(rb) = f(rb) – the value of potential is prescribed on the boundary by an arbitrary function (usually constant).
Neumann dϕdr|r=rb =g(rb)– the value of electric field is prescribed on the bound-ary by an arbitrbound-ary function(usually whenE(rb) = 0). It is possible to apply this condition only to one of E(rb) components.
Cauchy ϕ(rb) = f(rb) and dϕdr|r=rb = g(rb) – in the case of both prescribed potential and electric field on the boundary.
These conditions can be different for the boundaries present in the model, since e. g. behaviour of boundary of conducting objects is different from the be-haviour of neighbouring plasma. Similarly, particles should be treated differently according to the physical origin of the boundary itself. Possible treatments are:
Periodic boundary Particles crossing the boundary on one side of the simu-lated volume should be re-injected on the opposite side.
Particle sink After crossing this type of boundary, the particle is discarded and recorded according to the nature of the boundary origin (e. g. the current is recorded, a surface charge is added, etc.).
Treating particles this way does not bring any new ones into the system and does not replenish the supply of particles in case any of them is lost. In models
4.2. Particle-in-cell model of the sheath, one or more boundaries are usually representing the bulk plasma, which is the source of particles inside the simulation volume in reality. Then the third particle boundary condition, theparticle source, is to be applied.
In order to inject the particles accurately, several conditions have to be met.
Most important one is that the energy distribution function is physical, i. e. the velocities of injected particles have the appropriate physical meaning including the orientation of the velocity vector. Second condition is that the injection leads to an appropriate increase the density of particles specified by the plasma model, and the final one is the correct spatial distribution during the injection. In case of simulation of plasma in magnetic field, an injection box is usually added to the boundary from which the particles are being injected. This is due to the fact that during Larmor gyration, the particles would cross the boundary multiple times and thus be lost, since this boundary usually acts a particle sink as well.
4.2.5 Model run
Since the kinetic simulation is performed in time steps, there is a model control loop that is being evaluated in each iteration. Usually, it proceeds as follows:
1. Particle pushingSolution of the Newton equation of motion (4.1). Eval-uation of the boundary conditions on particles.
2. Charge weighting Calculation of the charge density.
3. Field calculationSolution of the Poisson equation.
4. Force calculation Calculation of forces acting on particles.
Before the loop gets off, model has to be initialized. This process includes a setup of the grid as well as objects inside the simulated volume, the assembly of the Poisson equation matrix and initial injection of particles. Afterwards, the model loop begins. Since there are no particles in the system at the beginning, a certain time is needed for their numbers to saturate. In PIC models with mag-netic field, this time is given by velocity perpendicular to the most prominent particle sink (usually representing the wall). If the saturation point is reached, artificial diagnostics and averaging can be performed. Particle model is available to track almost every information about particles at the cost of operational mem-ory. Most common tracked variables include potential in the plasma, particle velocities and particle and energy fluxes on surfaces of immersed objects.
5. SPICE model
SPICE (sheath particle-in-cell) model family is a set of PIC models developed for investigation of finite Larmor effects on deposition of particle deposition into gaps between PFC tiles by J. P. Gunn and R. Dejarnac [50, 95, 96] and later on by M. Komm [52, 59, 60, 97]. It is based upon previous PIC models developed by J. P. Gunn [68, 98, 99]. The author of this thesis joined the development of the team with the aim on implementation of the parallel solver of the Poisson equation.
In this chapter, key features of SPICE2 and SPICE3 will be explained and the work on the Poisson solver equation will be summarized.
5.1 Model highlights
There are several key features of SPICE model family shared between both 2D and 3D alternatives of the code. First, the plasma in the simulation is considered to be collisionless. This regime is reached by meeting certain conditions related to the collision frequency (5.1)
νcoll = πe4n0
(4πε0)2(kBTe)3/2 ln Λ, (5.1) where the Coulomb logarithm ln Λ [24, p. 176–181] describes the mean impact parameter. Using the mean free path defined as Lmfp = vth/νcoll and the total system size L, we can write these conditions as
Lmfp >√︁
L2+ (2πL)2, the mean free path should be longer than the
traversed distance, and (5.2)
νcoll <min(2πωc, νpl, . . .), the collision frequency should be smaller
than any other characteristic frequency. (5.3) The traversed distance has two components since the simulation box of sizeL is traversed by a gyrating particle with mean velocity cs. Linear approximation is used, but due to the strong electric field, the actual value can be much lower.
The second criterion sets the limit for fast particles that cross the box within one gyration period and relates the collisional frequency to other oscillations present in the system (such as those given by the plasma frequency νpl).
When simulating the plasma withTe≥10 eV, both conditions (5.2) and (5.3) are effortlessly met as can be seen in the fig. 5.1. Any simulation falling within the colored area for a certainTecan be considered collisionless. The curves shown in the figure were calculated for fixed B = 1 T and selected values of Te. The operational space grows with B magnitude.
5.1.1 Input and normalization
Both 2D3V (SPICE2) and 3D3V (SPICE3) (three/two spatial dimensions and three velocity dimensions) variants of the code exist, sharing most of the key
5.1. Model highlights
Figure 5.1: Collisionless regime vs. SPICE operational parameters. Filled areas indicate validity of the collisionless approach for certain simulation size Land input densityn0. Curved part of the boundary is given by the condition (5.2), while the straight part is given by (5.3) and extends the operational space with increasingB.
characteristic with appropriate implementation. The model itself is written in Fortran 90 and can be compiled with various compilers (Intel, PGI, Pathscale [97]), although currently, the Intel Parallel Studio with ifort/mpiifort is the most commonly used option. The code works with dimensionless quantities on Cartesian grid (rectangular and cuboid). The magnetic field is homogeneous and it is specified by its magnitude and direction [97]. Direct Poisson solver is implemented in SPICE2 and SPICE3 in both single-threaded [100] and multi-threaded (this thesis) variants, while 2D code uses a direct solver, and the 3D variant is based on multigrid approach.
To run the simulation, a proper input should be specified. SPICE family derives all its operational parameter from a few quantities that are specified in the input file. The main parameters describing the plasma are:
Magnetization parameter: ξ = rL
λD (5.4)
Temperature ratio: τ = Ti
Te (5.5)
Mass ratio: µ= mi
me (5.6)
B direction: αB = (αxz, αyz) (5.7)
Both the Larmor radius and the Debye length are calculated in the cold ion approximation, since it is simple to calculate the hot ion values when the τ is known.
Geometry of the situation has to be specified as well. The grid is Cartesian with base unit of d/λD = 1 in all directions. However, this can be overridden to simulate finer grids. Upon this grid with dimensions lx, ly and lz, several types of objects can be built.
Many simulations of different plasma conditions can be simulated while having similar plasma response. To address multiple issues at once and to minimize the rounding error when treating numbers with different orders of magnitude, the dimensionless normalization is a useful tool for code optimization and for smooth review of results, since the normalization coefficients are typical values of normalized quantities (marked as ′). In SPICE models, the normalization is performed as follows [97]:
Time step: ∆t′ = ∆t·ωc, i (5.8)
Particle position: r′ =r/λD (5.9)
Input velocity: v′in=vin/cs, cold (5.10)
Particle velocity: v′p =vp/(ωc, iλD) (5.11)
Potential: ϕ′ =ϕ·e/kBTe (5.12)
Particle mass: m′p =mp/mi, main (5.13)
Particle charge: qp′ =qp/qi, main (5.14)
Species density: n′i =ni/n0 (5.15)
In the particle mass normalisation, a particular ion species is selected as the main ion, with respective massmi, main. This facilitates the introduction of mul-tiple ion species.
To adhere to time step conditions, the criterion intervals are modified as follows: The cell crossing time relates to the cell diagonal, thus d term in (4.5) becomes λD√︂
[dx]2+d2y+d2z. The stability criterion time (4.6) is shortened to 1/10th of the gyration period for better resolution of electron trajectories. It is to be noted that these critical time step intervals depend on ion to electron mass ratio µ. For most of simulations performed by SPICE, the value of µ = 200 is chosen, since the duration of the simulations shrinks substantially while keeping the physical relevance. However, some situations require the real value ofµ, which in case of deuterium plasma is µ= 3670.
The duration of the simulation runtsim is defined in multiples of ion perpen-dicular travel time
tz = lz
vi,∥sinαB, (5.16)
and to reach the saturation of particle count within the simulation, several ion passes are required, with tsim= 5tz being sufficient enough.
5.1.2 Particle injection
Since a new SPICE simulation starts with an empty simulation space, proper care should be taken when injecting particles into it. The injection takes place in additional space called the injection box. The importance of the injection box lies in the fact that ions injected exactly at z =lz would be unintentionally discarded when they cross the boundary for after a few time steps after their injection (see fig. 5.2). Addition of the injection box prevents this situation to happen. The
5.1. Model highlights length of the injection box is ∼rLcosαB and the box transforms to an injection plane for B ∥z.
Two scenarios of particle injection are supported. First one is intended for the study of plasma in direct contact with the wall represented by an infinite plane. In this case, the injection box is located at z = lz and the simulation box is periodic in x and y. Objects and the plane z = 0 serve as particle sinks.
Second injection scenario is intended for simulation of plasmas where the B field is parallel to the surface. To prevent formation of shadows, particles are injected from both planes z = 0 and z =lz. Zero drift velocity is considered. Periodicity is the same as in the one injection plane scenario. The overview of these scenarios is shown in the fig. 5.2.
Figure 5.2: Both SPICE models support two scenarios of particle injection. Either the particles are injected from an injection box located at z = lz (var. A), or they are injected from two planes located at z = 0 and z =lz (var. B). The injection scenario B is supported only for B inz direction. There are three types of boundary conditions applied on particles, periodic, source and sink. Three dimensional situation is similar, boundary conditions in thex–zplane are the same as in y–z plane.
Injected particles have to have a proper energy distribution. Multiple variants of energy distribution function (EDF) are supported via input files. Currently, two distribution functions are supplied with SPICE. The first is the straightfor-ward Maxwellian EDF. While being approximately fulfilled in the bulk plasma, both the presheath and the magnetic field disrupt this distribution.
First-choice EDF for the simulation of the sheath is the Chung-Hutchinson distribution function [101]. This function is numerically calculated by and pre-set functions are available for several values ofτ in range from0.1to5. In the initial development by J. P. Gunn, the injection schema divides the particle velocity into components related to the magnetic field. The parallel velocity component v∥ is injected with the modified velocity distribution function function calcu-lated by the QPIC, which is a 1D presheath simulation model [99]. Components perpendicular to B (v⊥1 and v⊥2) are Maxwellian and the phase of the Larmor gyration motion is uniformly distributed. Then, the velocities are transformed
into Cartesian coordinate system, where B = (Bx, By, Bz) is inclined with an angle αB =(︁
αBx, αBy, αBz)︁
as defined in [97]:
vx =v∥cosαBx − v⊥1cosαBy +v⊥2cosαBxcosαBz
√︁cos2αBx + cos2αBy , vy =v∥cosαBy +v⊥1cosαBx−v⊥2cosαBycosαBz
√︁cos2αBx+ cos2αBy , vz =v∥cosαBz +v⊥2
√︂
cos2αBx+ cos2αBy,
(5.17)
which assures the uniform distribution of the gyration phase. The particle is positioned on a randomly generated place around z = zinj which lies within the injection box. The position in x and y coordinates is uniformly distributed between 0 and lx (ly respectively). To avoid artificial modes in the simulation, z coordinate of i-th particle is randomized [97]:
zi =zinj+v∥cosαBzrand(0, ∆t) +qiv⊥1mi
√︂
cos2αBx + cos2αBy, (5.18) by a small deviation in the parallel direction and with the assumption that the particle with mass mi and qi performs the Larmor gyration around a particular field line. After the particle is generated, its motion is guided by a standard leap-frog algorithm described in subsection 4.2.3.
5.1.3 Parallelization scheme
To speed up the simulation, it is beneficial to distribute the load between several processors or program threads. The SPICE code family uses the former approach with distributed memory parallelization via Message Passing Interface (MPI).
To compile the sources, compilers provided by various implementations of MPI (BullxMPI, OpenMPI and Intel MPI) were successfully utilized, while the latter two are now the preferred choice.
Optimal scheme has to be chosen to perform the calculation efficiently. Previ-ously, the parallelization by grouping of particles was used, however, a substantial disadvantage of this approach is that each thread requires the complete informa-tion about working field (the potential distribuinforma-tion and objects locainforma-tion), which is rather large (especially in 3D). This scheme was later replaced by the domain decomposition. The domains are even slices (rectangular or cuboid) of the sim-ulation volume specified in the input file by a number of cuts in each direction.
Limits of this scheme lie within the size of the slice, which should be around 128×128 cells in 2D. SPICE scales efficiently up to64processors, going further, the inter-process communication slows down the simulation beyond ideal scaling.