A Set of Benchmark Tests for Validation of 3D Particle In Cell Methods
S. O'Connor, Z. D. Crawford, J. Verboncoeur, J. Lugisland, B. Shanker
11 A Set of Benchmark Tests for Validation of 3DParticle In Cell Methods
S. O’Connor,
Student Member, IEEE,
Z. D. Crawford,
Student Member, IEEE,
J. Verboncoeur,
Fellow, IEEE,
J.Lugisland,
Fellow, IEEE,
B. Shanker,
Fellow, IEEE
Abstract —While the particle-in-cell (PIC) method is quite ma-ture, verification and validation of both newly developed methodsand individual codes has largely focused on an idiosyncraticchoice of a few test cases. Many of these test cases involveeither one- or two-dimensional simulations. This is either due toavailability of (quasi) analytic solutions or historical reasons. Ad-ditionally, tests often focus on investigation of particular physicsproblems, such as particle emission or collisions, and do notnecessarily study the combined impact of the suite of algorithmsnecessary for a full featured PIC code. As three dimensional (3D)codes become the norm, there is a lack of benchmarks test thatcan establish the validity of these codes; existing papers eitherdo not delve into the details of the numerical experiment orprovide other measurable numeric metrics (such as noise) thatare outcomes of the simulation.This paper seeks to provide several test cases that can be usedfor validation and bench-marking of particle in cell codes in 3D.We focus on examples that are collisionless, and can be run witha reasonable amount of computational power.Four test cases are presented in significant detail; these include,basic particle motion, beam expansion, adiabatic expansion ofplasma, and two stream instability. All presented cases arecompared either against existing analytical data or other codes.We anticipate that these cases should help fill the void of bench-marking and validation problems and help the development ofnew particle in cell codes.
Index Terms —Particle in Cell, PIC, Finite Element Method,Maxwell Solvers, Vlasov Equations, Plasma, Space Charge
I. I
NTRODUCTION S IMULATION of space charge and plasma is an integralpart of many scientific and engineering processes, withspecific applications of pulsed power, particle accelerators,directed energy, integrated chip manufacture, satellites, andmedicine, to name a few examples. A number of methodshave been used to simulate the underlying physics; these rangefrom molecular dynamics, many-body, global-model, 2-fluid,magneto-hydrodynamics (MHD), etc. [1]. The key with all ofthese methods is the self-consistent evolution of the fields andplasma. In the area of kinetic computational plasma physics,this means evolving the particle distribution function in amanner consistent with dynamic electric and magnetic fields.One such method, the Particle-in-Cell (PIC), does this by (a)mapping all particle-field interactions through a mesh, and(b) using macro-particles as statistically significant markersof the distribution function’s phasespace [2], [3]. PIC has
S. O’Connor, Z. D. Crawford, J. Verboncour, J. Lugisland, B. Shanker arewith the Department of Electrical and Computer Engineering, Michigan StateUniversity, East Lansing, MI, 48824.E-mail: [email protected] been around for several decades, originating in the 1960’s.Much of the early work focused on using finite-differencetime-domain (FDTD) for electromagnetic fields solutions ona regular grid, more specifically, Yee Cells [4]. This methodhas become popular due to its simplicity and scalabilty [5].Additionally, these methods are local, which makes parallelimplementation relatively straight forward. While efficient,curved surfaces, or any geometry that is not aligned with aregular structured grid, is poorly represented. Significant workhas been done to adapt Cartesian grids with better modelboundaries. The most popular is the Dey-Mittra adaption ofFDTD to PIC [6]. It should be noted that there have been anumber of other advances in conformal FDTD [7], but theyhave not made their way into PIC codes. An alternative isto use an axisymmetric methods when possible by treating athree dimensional problem in a two-dimensional setting. Anexcellent example is XOOPIC, which is open source and hasbeen used extensively [8].Since the introduction of FDTD, the fields community hasinvested time and resources to developing high fidelity fieldsolvers. These include finite element method (FEM) and dis-continuous Galerkin methods (DG). Both can provide higherorder representation of fields and geometry. Additionally, onehas better understanding of efficient time marching schemes[9]. Here, our focus is finite element based Maxwell solver,that have been a) adapted to PIC and b) shown to be devoidof several problems associated with charge conservation.The challenge with using unstructured grids, using FEM orDG, is developing a consistent scheme to solve the coupledproblem of evolving fields and particle population. Giventhe relative lack of maturity of FEM-PIC, it is important tounderstand how the nuances of these algorithms affect theaccuracy of the physics that is predicted. Small variationsin method can have significantly different prediction; for in-stance, time stepping method, solving the first order Maxwell’sequations vs. solving the second order wave equation, orderof spatial basis functions, order of time basis function, etc.,can yield different predictions. Perhaps, the most critical ofthese nuances is charge conservation in that it explicitlyties together the evolution of the fields and the particles.Additionally, while there has been extensive efforts for FEMin the context of source-free electromagnetics, FEM-PIC ismuch less developed. The most recent of these is the work in[10] where a charge conserving scheme for Whitney basis waspresented. This has been developed further by Teixeira [11].This begs the question as to how PIC codes are validated.Correctness of a simulation has typically been done through a r X i v : . [ phy s i c s . c o m p - ph ] J a n a number of methods such as simple property test such asenergy conservation, to more rigorous tests such as rates ofconvergence. These test are described in more rigor in [12].The most rigorous tests require a comparison with analyticalsolutions. But, this is much harder to come by in PIC settings,as analytic solution for most problems do not exist. Most testsrely either demonstrating proper behavior, or comparison toanother code. As a result, there is a need for a set of testexamples that can be used to validate/benchmark and see howthese variations affect the solution of the solver. Validationand verification methods play import roles in quantifying error,numerical noise, and comparison of method performance. Oflate, in the PIC community, there has been a push for morerobust bench marking and validation of (PIC) codes [12]–[14]. Much of this current literature focuses on dischargeand particle emission and collisions. Less attention is paidto test examples where the grid-particle interaction dominatethe numerical effects. This is the void we seek to fill.In collision-less PIC, there are a number of metrics one maywant to measure such as field noise, charge conservation andenergy fluctuation. A set of benchmark tests are needed tomeasure each of these metrics. In this paper, we have deviseda set of test cases to measure various metrics of particle incell codes to both verify the code and allow for comparisonsbetween codes. The four cases presented here have (quasi)-analytic solutions, and can be verified against other methodssuch as an electrostatic PIC codes. We provide details of eachcase such that comparison and validation is easy. Each ofthese cases are a 3-D example that can be set up with simpleboundary conditions – perfect electric conductor (PEC) or firstorder absorbing boundary condition (ABC). The domain sizeto capture the physic of the problem is limited such that largeamount of computational resources are not needed, but can bescaled in a parallel context, to test algorithm scaling.II. P ROBLEM D ESCRIPTION
Next, we present a succinct outline of the problem descrip-tion and implementation details for the test cases presented.Consider a domain Ω enclosed by a boundary ∂ Ω . The domainconsists of a free space permittivity and permeability definedby (cid:15) and µ , and can contain multiple charge species. Both thefields and charge distribution evolve over time self consistentlywith accordance to both the Vlasov Eq. (1) and Maxwell’s Eq.(2), ∂ t f ( t, r , v ) + v · ∇ f ( t, r , v )+ (1) qm [ E ( t, r ) + v × B ( t, r )] · ∇ v f ( t, r , v ) = 0 . Here, f ( t, r , v ) is the phase space distribution function.Typically, one does not solve Eq. (1) directly. The usualapproach, that we follow as well, is to define charge andcurrent via the first and second moment of this distributionfunction; i.e., ρ ( t, r ) = q (cid:82) Ω f ( t, r , v ) d v and J ( t, r ) = q (cid:82) Ω v ( t ) f ( t, r , v ) d v . These are then solved self-consistentlywith Eq. (2) and the particle location are updated using theLorentz force F = q ( E + v × B ) . − ∂ B ( t, r ) ∂t = ∇ × E ( t, r ) (2a) ∂ D ( t, r ) ∂t = ∇ × H ( t, r ) − J ( t, r ) (2b) ∇ · B ( t, r ) = 0 (2c) ∇ · D ( t, r ) = ρ ( t, r ) (2d)We assume that the boundary can be subdivided into eitherDirichlet Γ D , Nuemann Γ N or impedance Γ I , boundary con-ditions as defined by, ˆ n × E ( r , t ) = Ψ D ( r , t ) on Γ D , (3a) ˆ n × B ( r , t ) µ = Ψ N ( r , t ) on Γ N , (3b) ˆ n × B ( r , t ) µ − Y ˆ n × ˆ n × E ( r , t ) = Ψ I ( r , t ) on Γ I . (3c)where ˆ n is an outward pointing normal to Γ D , Γ N or Γ I .We define Y = (cid:112) (cid:15) /µ as the free space surface admittance,and Ψ D ( r , t ) , Ψ N ( r , t ) , Ψ I ( r , t ) are boundary condition func-tions. Initial conditions, must satisfy Gauss’s electric and mag-netic law, such that ∇· D (0 , r ) = ρ (0 , r ) and ∇· B (0 , r ) = 0 .The algorithm’s preservation of these conditions can then betested by our problems. A. The Discrete Domain
To solve the system, we use a 3D finite element Maxwellsolver. Consider a domain Ω broken into an irregular non-overlapping tetrahedral elements with boundaries ∂ Ω . In eachtetrahedral element the field quantities E and B can be repre-sented using Whitney spaces. The moments of the distributionfunction can be represented using using delta functions, shownin Eq. (4a),(4b), ρ α ( t, r ) = q α N p (cid:88) p δ ( r − r p ( t )) (4a) J α ( t, r ) = q α N p (cid:88) p v p ( t ) δ ( r − r p ( t )) (4b)for each charge species α , where q α is the charge of the macro-particle and N p is the number of macro-particles. The systemproceeds as follows. The Whitney basis functions are usedto convert Maxwell’s continuous equation to discrete ones.Current and charge are mapped to the grid via Whitney-1 (cid:0) W (1) i (cid:1) and Whitney-0 (cid:0) W (0) j (cid:1) basis sets, and the solutionto Faraday’s and Ampere’s law uses a leap frog scheme.Details to pertinent to each stage of the scheme can befound in [10], [11]. Using E ( t, r ) = (cid:80) N e i =1 e i ( t ) W (1) i ( r ) and B ( t, r ) = (cid:80) N f i =1 b i ( t ) W (2) i ( r ) and Galerkin testing results in, B n +1 / = B n − / − ∆ t [ D curl ] · E n (5) E n +1 = E n + ∆ t [ (cid:63) (cid:15) ] − · (cid:18) [ D curl ] · [ (cid:63) µ − ] · B n +1 / − J n +1 / (cid:19) (6) where N e and N f are the number of edges and faces, ∆ t isthe time step size, E n = [ e n , e n , ..., e nN e ] , B n +1 / = [ b n +1 / , e n +1 / , ..., b n +1 / N f ] , J n +1 / = [ j n +1 / , j n +1 / , ..., j n +1 / N e ] with E n , B n +1 / , and J n +1 / being the field coefficients attime step n and n +1 / . The Hodge and discrete curl matricesin Eq. (5) and Eq. (6) are, [ (cid:63) (cid:15) ] i,j = (cid:104) W (1) i ( r ) , ε · W (1) j ( r ) (cid:105) (7) [ (cid:63) µ − ] i,j = (cid:104) W (2) i ( r ) , µ − · W (2) j ( r ) (cid:105) (8) [ D curl ] i,j = (cid:104) ˆn i , ∇ × W (1) j ( r ) (cid:105) (9)where, [ D curl ] is the discrete curl matrix, [ (cid:63) µ − ] is the Hodgestar matrix, [ (cid:63) (cid:15) ] is the Hodge epsilon matrix. The current ismapped using a charge conserving method presented in [10],[11], j n +1 / i = q p ∆ t (cid:90) r n +1 p r np W (1) i ( r ) · d l (10) = q p ∆ t [ W (0) i ( r np ) W (0) j ( r n +1 p ) − W (0) i ( r n +1 p ) W (0) j ( r np )] , (11)where W (0) i ( r ) is the barycentric coordinate for node i withina cell. III. B ENCHMARK T ESTS
In this section, we present four cases that could be used forvalidation. In each of these cases, we will present a numberof metrics and discuss what would happen if some of thesemetrics were not to be satisfied. Each test is designed tohighlight an aspect of the method. The first test case, basicparticle motion, highlights the particle push’s performance.The expanding particle beam is setup to verify current map-ping, discrete Gauss’s magnetic and electric law, and thediscrete continuity equation. The adiabatic expanding plasmatest case allows for multi-species evolving distributions andcan be scaled up to test the scalability of methods. Two streaminstability displays more dynamic features of the method andhas analytic solution based on the rise time of the simulation.
A. Basic Particle Motion
The first tests of any PIC method should start with validatingthe particle motion. To this end, we set up two tests, 1) aparticle moves under the influence of an external field dueto Newton’s laws and 2) the motion of a particle is trackedunder the influence of both a background external field andgrid interaction solved by Maxwell’s equation and Newton’slaws. We use Newton’s law to obtain position and velocityat every time point. Such a test verified that the trajectory ofthe particle — as mapped on a grid — is computed correctly.Given that this trajectory can be obtained analytically, one canobtain error bounds and relate them to time step size. Next this test can be rerun assuming radiation, that is solv-ing Newton and Maxwell’s equation self consistently. If thesolutions are consistent, the two path should be approximatelyidentical as the radiated fields are significantly smaller than theDC fields. All results that bear this exertion are shown in Fig3. The parameters chosen are such that force due to the particlemotion is significantly less than the force due to the DCbackground fields. This choice causes the main source of errorin the particle trajectories to be dominated by the particle pushand not the field solver. This setup also presents an opportunityto verify current continuity and the discrete Gauss’s law in asetup where the particle path is for all intents and purposes isknown. Although the field values are small, they are non-zero.This case is simple enough that it can be performed on a smallenough meshes that any errors are easily identified. a) External E-Field:
The first term of the Lorentz forcecan be checked by placing a particle at rest in an uniformbackground electric field of V/m and letting the particleaccelerate as shown in Fig. 1. The macro-particle has themass and charge of a single electron. We calculate the analyticmotion and the error we get from the particle push method. b) External B-Field:
To check the v × B term we, applya uniform external magnetic flux density B = 6 . · − ˆ zW b/m with no electric field. A single particle with the massand charge of an single electron is given an initial velocity v = 3 · ˆ y m/s and position r = [0 . , . , which resultsin a gyro-radius of . m, as shown in Fig. 2. We analyze theerror in Fig. 3 by comparing four test cases; the motion due toa fully coupled system labeled mesh , the background fieldswith the current time step size, background fields with thetime step size 10 time large, and background fields with thetime step 10 time smaller. The error is calculated by findingthe distance between the analytic and simulated locations at agiven step for each case. The two dips in the error stem fromthe particle starting at the same location — zero error — anddrifting array from the analytic circle and once the rotation iscomplete ending up at the start location. In Fig. 4, we comparethe average error over one cycle of cyclotron motion vs. thenumber of step per cycle. Fig. 1. Particle acceleration due to electric field.Fig. 2. Cyclotron particle motion of an electron in a magnetic field.Fig. 3. Cyclotron motion error from just Boris push with several differenttime step sizes, as well as influence from the mesh. Fig. 4. Error as time step size decreases for Boris push.
B. Expanding Particle Beam
In the second test, a beam of electrons are injected into acylindrical drift tube and is evolved over time. The reasons forthis test are 1) it has an approximate solution found in [15];2) charge conservation is easy to measure and any error showsup in the solution as higher noise levels or spatial striationsin the beam [16].We create a cylindrical drift tube geometry as shown in Fig.5, with a radius of . m and a length of . m. We set allexterior surfaces to perfect conductors (PEC) with the interiorvolume set to vacuum. The drift tube is discretized such thatthe average edge length is . mm. The beam has a voltageand current of 7.107 kV and 0.25 A. This setup amounts to abeam radius of mm centered in the cylinder with a startingvelocity of v = [0 . , . , · ] m/s. The beam is rampedup by changing the macroparticle weight from 0 to 52012electrons per macroparticle linearly over the first 2 ns of thesimulation. The parameters are listed in Table I. The current TABLE IE
XPANDING P ARTICLE B EAM P ARAMETERS
Parameter Value
Cavity Radius 20 mmCavity Length 100 mmBoundary Conditions PEC v p · m/s v p /c r b in this example was chosen such that a converged solutionsof the expansion could be represented with a coarse mesh toreduced the simulation time needed to solve the problem. Wealso compare the results with the 2-D FDTD code XOOPICin r - θ mode. As is evident, the agreement between the twocodes — despite obvious advantages afforded by asymmetryin XOOPIC — agree well with each other. The trajectories ofthe particle is shown in Fig. 5, and there is excellent agreement Fig. 5. XOOPIC vs FEMPIC for use in particle in cell between FEMPIC and XOOPIC. While FEMPIC and XOOPICagree with each other, they differ slightly from the analyticdata, although there is good agreement on the overall trend.The analytic solution derived in [15] is also shown in Fig. 5,it assumes an idealized beam, with uniform density radiallyand ignores the perfect electric conductor (PEC) boundaryconditions on the entrance and exit surfaces of the beam. Wepostulate that this leads to the slight disagreement. Note, itis difficult, if not impossible to create the exactly identicalconditions for either the numerical simulation or analyticalanalysis. Nerveless, we include this data to show that the sametrends are captured.Next, the kinetic and potential energy in both the electricand magnetic field obtained by FEMPIC and XOOPIC agreeas shown in Fig. 6. We also examined field noise. Note, theboundaries are perfect electrically conducting, and there areno losses in this problem. This results in cavity modes thatare excited. The radial component of the electric field, E r ,at halfway through the tube at a radial distance of cmfrom the center of the tube is shown for three cases in Fig.7. Each case uses a different random number generator seedfor the location of where the particle in the beam is injected.This allows us to see if the noise is random or does it fallinto certain modes. Various modes of the cylindrical cavityare shown in Fig. 8. These show up upon taking the Fouriertransform of the noise after the system reaches steady state at . ns and the DC component had been removed. Non-randominjection does artificially reduces the noise of the system. Wealso measured the discrete continuity equation as a functionof time in Fig. 9, it is evident charge is conserved to machineprecision. C. Adiabatic Expanding Plasmas
In the third test, we simulate an adiabatically expandingplasma ball with a Gaussian distribution in the radial direction.This experiment has an approximate analytic solutions [17]that can used for validation. These analytic solutions havebeen compared with experimental results in [18]. In earlier
Fig. 6. Particle Beam Energy with a 2ns turn on timeFig. 7. Noise measured from three beam test cases with different randomnumber generatorFig. 8. Fourier transform of the electric field for three beam tests that haddifferent random number generator seeds.
Fig. 9. Charge Conservation for particle beam test case measured by | ( ρ n +1 − ρ n ) / ∆ t − [ ∇· ] J n +1 / | / | ( ρ n +1 − ρ n ) / ∆ t | /N part Fig. 10. Distribution of electrons radially at two instants in time with analyticsolutions. examples, the analytical examples are obtained in 1-D and2-D, making it difficult to show a proper expansion with 3-D numerical code [11]. As analytical data is available in 3-D, we design the test to reflect the experiment done in [18]but change the density such that the Debye length can fullyresolved with modest computational resources. The expandingplasma has two species, Sr+ ions with an initial temperatureof 1 K and electrons at 100 K. The initial number density ofboth species is e particle per m − , and the initial σ = 0 . msuch that the system is initial charge neutral everywhere. Otherdetails of the simulation parameters can be found in TableII. When the simulation starts, the electrons expand outwardfaster than ions, due to the higher temperature, creating a radialelectric field. After the initial expansion outward, the growingelectric field slows the expansion of the electrons. At the sametime the ions expand outward pulled by the electrons. Thisexample can be scaled – by changing the density of charges– such that various levels of computational power are needed TABLE IIA
DIABATIC E XPANDING P LASMAS
Parameter Value
Mesh Radius 6mmBoundary Conditions First order ABC T ion T electron Sr + Macro-Particle Size 52012.58Min Edge Length 1.529mmMax Edge Length 6.872mm thus making it a good test of the scalability of a code. We usean first order absorbing boundary condition (ABC) to truncatethe domain and place the boundary conditions sufficiently farfrom the initial distribution to capture the initial expansion ofthe electrons but not exit the mesh. Early electron expansioncan be validated in small geometries, which we do here. Theanalytic distribution of electrons at initial conditions and attime 168 ns, are show in Fig. 10. This was derived in [18],and is succinctly presented here. The density function f α ofspecies α can be expanded as, f α ( t, r , v ) = n (cid:20) m α πk b T α (0) (cid:21) / exp (cid:26) − r σ α ( t ) − m α v k b T α ( t ) (cid:27) . (12)The radial distribution is a symmetric Gaussian distributionfunction with a deviation σ α ( t ) is a function of time, σ α ( t ) = σ α (0)(1 + t /τ exp ) . (13)The velocity distribution depends on the temperature T α ( t ) ofspecies α , T α ( t ) = T α (0)1 + t τ exp , (14)Both σ α ( t ) and T α ( t ) depend on the characteristic expansiontime τ exp τ exp = (cid:115) m α σ α (0) k b [ T e (0) + T i (0)] . (15)Integrating velocity space, we can write the number density η ( t, r ) as a function of time and radius from center of thesphere, η ( t, r ) = η (cid:20) T α ( t ) T α (0) (cid:21) / e − r σ α . (16)This derivation was used in creating Fig. 10. D. Two Stream Instability
Next, we present the two stream instability test case fromBirdsall and Langdon [2]. We shall follow its prescription andparaphrase it here. Two streams of particles are setup, eachwith a density of n and the same sign. This test case waschosen due to analytic solutions associated with the growthrate of the instability. While it is inherently a 1-D problemand has been done many times before by others, it is a classicproblem that can adapted to 3-D and provides insight into codeperformance. Typically, this test is performed by using periodicboundary conditions. Alternatively, when that functionality is TABLE IIIT WO S TREAM I NSTABILITY P ARAMETERS
Parameter Value
Boundary Conditions PeriodicSolver Type electrostaticspecies type electronsdensity . · cm − v · m/sMesh Edge Length 5 cmNumber particles 1000macro-particle size 9000.0length of domain 2 mFig. 11. Two stream phase diagram for electrostatic finite element solver at184 ns into the simulation with 400 particles. not available, one can use perfectly electric conductor (PEC),pack the volume with charge with the two velocity setup,and emit particles from both PEC surfaces with the samevelocity, such that you get a similar system. PEC providesimage charge and thus retains the same physical features of theproblem in a physically realizable configuration. The specificexample parameter details can be found in Table III. Thetwo-stream example shown here produces the standard swirlson the phases diagram shown in Fig. 11. The energy of thesystem, shown in Fig. 12, is conserved. The small bumps inenergy in the field solve is caused from electrons moving tooclose to one another causing spikes in electric field values.IV. S UMMARY
In this paper, we have provided four tests that can be usedto validate 3D PIC codes. These tests have been developed as(a) there is significant recent interest in developing PIC codesto fully resolve complex boundaries, as properties of theseboundaries, and (b) as there is a void in examples with detailsthat can be used to validate these codes. In all examples that wehave developed, there exist quasi-analytic solutions, as well asopen source codes that can be used for benchmarking results.We have also provided sufficient detail so as to facilitate theuser to examine all possible metrics that may subtly affect thephysics.
Fig. 12. Energy plot for two stream instability with 3D periodic boundaryconditions. A CKNOWLEDGMENT
This work was supported by SMART Scholarship program.We thanks the MSU Foundation for support through theStrategic Partnership Grant. The authors would also like tothank the HPCC Facility, Michigan State University, EastLansing, MI, USA. R