Particle simulation methods in plasma physics

Computer Physics Communications 43(1986) 89—106 North-Holland, Amsterdam

89

PARTICLE SIMULATION METHODS IN PLASMA PHYSICS James W. EASTWOOD Cuiham Laboratory, Abingdon, Oxfordshire 0X14 3DB, UK (EURA TOM/ UKAEA Ft~sionAssociation)

A brief survey of particle simulation methods is presented. Section 2 outlines the particle—mesh (PM) method for collisionless phase fluids. Section 3 describes how PM is extended to make simulation studies of dense plasmas feasible: The 3M) algorithms reduces the cost scaling to aN~from the aN~for conventional methods, particle—particle/particle—mesh (P where N~is the number of particles. Section 4 summarises particle modelling of fluids and magnetofluids. New algorithms (EPIC schemes) arise when the optimisation techniques of collisionless PM models are applied to the fluid and MHD equations. These are shown to outperform state-of-the-art finite difference schemes in both cost and performance.

I. Introduction The aspect of plasma physics where particle— mesh (PM) methods have been most frequently used over the past two decades is in the study of collisionless phase fluids. Initially, PM schemes were devised by appealing to physical argument the collisionless phase fluid equations are themselves derived from corpuscular descriptions of plasmas. Theoretical physics, rather than numerical analysis techniques, led to the detailed understanding we now have of such particle models. This understanding has led to optimisation for the (purely hyperbolic) collisionless case. In this paper, it is shown that optimal particle methods offer a cost effective solution to plasma and fluid modelling beyond the class of problem for which particle models are generally regarded as appropriate, One special feature making PM methods attractive for collisionless phase fluids is that only moments of the distribution are required for source terms in the field equations (see section 2). Consequently, no velocity space grid is required as the particle distribution carries all the required information. In addition, particle models are robust; too few particles simply enhances diffusion. Too little spatial resolution causes inverse Landau damping. Both move physical parameters to ranges accessible to the numerical model rather than causing catastrophic failure. —

Extending the PM method to describe kinetic effects (section 3) again leads to a system where only moments are required for the field source terms, but now the problem is one of insufficient spatial resolution: the P3M (particle—particle/ particle—mesh) method overcomes this by using the mesh only for the smoothly varying part of the force. The result is that Np-body Coulombic force sums can be evaluated in i22(N~)rather than operations, making simulation of large systems possible. Robustness and accuracy are the features of particle models which make particle modelling of Fokker—Planck, MHD and fluid equations attractive (section 4). Ephemeral Particle-In-Cell (EPIC) schemes automatically combine many special features of state-of-the-art finite difference schemes to give stable, accurate and cost effective numerical schemes. No attempt has been made in this paper to present an exhaustive set of references to particle methods. Two research monographs [1,2] describing particle methods and their applications provide a source of several hundred references. More recent references appear in refs. [3,4].

2. The collisionless phase fluid Particle—mesh methods for collisionless plasmas use a combination of particle and mesh repre-

OO1O-4655/86/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

90

J. W. Eastwood

/

Particle simulation methods in plasma physics

sentations. The plasma is represented by an ensemble of particles. The evolution of the system is described by the motions of these particles. The purpose of the mesh in electrostatic models is to provide a discrete representation of the charge density, the potential and the force fields. For electromagnetic models, this is extended to inelude current density, and electric and magnetic fields. As will be shown below (section 2.3), the use of fields on a mesh gives enormous computational speed gains over direct particle—particle force summation.

Eqs. (2) are represented by finite differences (fda’s). The computational box is divided into a set of cells of width H. where H L/N. At the centre of each cell is a mesh point at which values of charge density, potential and electric field are stored. Derivatives are approximated by differences of mesh values. The timestep cycle of the NGP model is as follows: i E [1, N~]}.where n a) Start: Given fx~’,v~’1~’2 denotes timelevel, time t nit. b) Assign charge to mesh.’ Charge density at mesh point p at position x~, pH is computed by summing all the charge in cell p (x~ H/2
=

=

2.1. The NGP scheme

—

We use the original and simplest PM method, nearest-grid-point (NGP) scheme [5], to illustrate

p~

=

1

Q+

~

-~

p

(3)

0.

the basic features of the PM calculation cycle. We consider a one dimensional electron gas embedded in a fixed neutralising background: This illustrates the essential features without complications of isolated boundaries and initial equilibria. PM schemes generalise trivially to two and three spatial dimensions [1—5]. The state of the 1-D NGP model at any timelevel is described by the set of particle positions and velocities (x,, v~,i E [1, N~]}.Initial condilions are prescribed by distributing the particles of mass M and charge Q throughout the computational box, x e [0, U, according to the initial distribution, f(x, v, 0), of density and velocity. Boundary conditions are prescribed for both partides and fields at the boundaries x 0 and ~ L of the computational box. For periodic conditions, particles moving through x 0 re-enter the com=

=

=

putational box at x U and vice-versa. The evolution of the PM model is computed by integrating the equations of motion of the partides: =

dx. dt

=

v;

do M—~- F(x~) QE(x,). dt =

=

(1)

particles in cell P

In practice, { ps,, p ~ [0, Ng

=

dx2

—

-p-; E ~o

=

—

~

dx

(2)

1]) are computed by

sweeping through the particles, locating the cell, p, containing each particle and incrementing p~by Q/H. c) Solve for the potential: the simplest fda to Poisson s equation gives —

P

E

2~+ [0, Ng

=

—

1~.

(4)

Eq. (4) with periodic boundary conditions may be inverted to obtain potentials {4~,,;p E [0, Ng 1]} d) Accelerate particles.’ NGP force interpolation takes the force field at any point in a cell to be equal to the value at the mesh point at the centre of that cell. Thus, for particle i in cell p the force is given by F(x,) QE(x~) Q(4~ ~~/2H (5) —

=

=

.~—

and velocities are updated using the difference approximation 2 V7~2 F” z~t/m. (6) v;’~’~ The acceleration phase, eqs. (5) and (6), are performed by sweeping through particles rather than cell-wise. e) Move particles; The new velocities { v”~1’~2 } —

The electric field E pis ofobtained from plus the charge density distribution the particles background charge by solving Poisson’s equation

—

=

J. W. Eastwood

are used to obtain new positions x’=x7+v7~1”2 ~t

/

Particle simulation methods in plasma physics

(7)

and it is at this stage that particle boundary conditions are applied.

obtained using weighted means of samples, i.e. fdvvnf_~Jdvv (N x+X/2 n~sf W(x—x’) X xf(x’, v,

2.2. The mathematical model

8f + ~f +

F Bf may d4 F=qE= -q~--, ~—

8t

ax

—

—

=

0,

t)

dx/}

=~~W(x—x,)v~

A common misconception concerning PM simulations is that of the nature of the physical process being simulated. NGP (and more sophisticated variants) provide means to simulate collisionless systems, and for a given number of simulation particles give a more accurate representation than would a set of point particles in a meshless system. The mathematical model for a 1 -D Vlasov— Poisson system plasma is [6] —

91

(13)

W is a weighting factor, A is the range of the averaging and N~is a scaling factor determined by mass (or charge) conservation. Setting A H and the weighting factor to W(x) 11(x), i.e., =

=

W(x)

=

\

/1, 0,

—H/2
(14)

reduces eq. (13) to the NGP assignment. A physical interpretation of collisionless PM

(8)

models is to view each particle as a finite sized cloud of electrons where internal degrees of free-

(9)

dom are suppressed. complete isdescription of the cloud (slab in oneThe dimension) given by its velocity, position, density profile and charge/mass

=

—

dx2 p

=

(10)

ratio. The cloud must have a finite width to soften short range interactions to the extent that they freely pass through each other, otherwise binary

(11)

collisional effects will dominate collective effects. Each cloud corresponds to N~electrons, where N~ is typically iO~,so without softening short range forces qualitatively incorrect behaviour would result. Generally, a cloud width (~H) of order the Debye length is used: This minimises collisional effects without seriously affecting collective properties. The value of charge density at a given point is given by the sum of contributions of clouds overlapping that point, so the range of neighbouring particles for purposes of finding weighted mean values may be interpreted as the cloud width (cf. eq. (13)). To avoid statistical fluctuations in such mean values, we need a large number (say 10) of clouds to overlap: This factor leads to the requirement that the number of particles per Debye slab

~o’

p 0

+

qff dv,

where f(x, v, t) is the one particle electron distnbution function. Eq. (8) is represented by particles and the integral in eq. (11) is evaluated by ‘Monte Carlo’ methods [7] in PM models. One interpretation of the particles in the PM models is simply as a mathematical artefact for transforming the awkward advection equation (8) into a set of ordinary differential equations (1). Eqs. (1) are the characteristics of eq. (8), along which f is constant: Given f at some time t the equations of motion (1) may be used to map f to some later time t’. The discrete approximation to f is a set of samples: N~

~ 8(x I

—

x.)8(v

—

vi),

(12)

1

where {x1, v.; i E [1, N~]}evolve according to eq. (1). f comprises samples of a smooth distribution function f, so estimates of moments of f are

(square or cube) is large. 2.3. PM vs. direct summation The finite sized particle interpretation of the NGP scheme begs the question, why use a mesh at

92

J. W. Eastwood

/

Particle simulation methods in plasma physics

all? The beneficial reduction in collisional effects arising from smoothing of forces by the mesh could equally well be achieved by summing interparticle forces (Coulomb’s law) between finite sized clouds. In this manner, all assignment, interpolation and mesh errors would be avoided. The answer is that the direct summation method is an untenable cost/quality compromise. The cost savings that can be achieved by employing particle—mesh schemes are readily seen by considering a simple example. The direct particleparticle (PP) scheme computes forces by pairwise sums: >f,,. Take, for example, a system of point charges where f,1 (x1 x1)/~x x, Assume + ± and X each count as one operation and ( )3/2 counts as three, then the operations count for evaluating F~for N~particles and advancing positions and velocities by one timestep is 10N~.To be compared with this is an operations count of aN~+ /3(N) for PM schemes, where the term aN~arises from charge assignment, force interpolation and integration of equations of motion. The term /3(N) is the operations count 3 for solving Poisson’s equation (~N 5N~ in 3-D). Using typical values a 20, log2N 32, N~ iO~and time, T per floating point operation of I p.s we obtain cpu time estimates of

NGP assignment (eq. (3)) as

Q p~,

=

0 +

~

—

=

=

=

PP:

+

iONp2 T

3

‘

SN log =

2N )T 4~s/timestep, 10 s/st I day/timestep. ‘

‘fr

The enormous speed advantage of PM over pp is obtained at the cost of increased algorithm complexity but with negligible loss of accuracy in the case of more sophisticated PM schemes, 2.4. Assignment and interpolation NGP is the lowest (‘zeroth’) order member of the hierarchy of PM schemes. It gives a force between particles which varies stepwise, is zero for particles within the same cell, is not displacement invariant and which depends on direction with respect to the mesh (‘angular anisotropy’). Generalisation of NGP assignment and interpolation may be obtained as follows. Using the ‘top-hat’ function, eq. (14) and the density of particle centres n(x) ff do we may express the =

—

=

p0 +

f W( x1

—

x’)

dx’ (15)

and interpolation as ,,

F(x,)

=

N5q

~

W(x,

—

x~)E1,.

(16)

p=O

—

—

PM: (tiN

W(x~, x)

X

=

=

N

p

Higher order schemes result from replacing the NGP form of W in eqs. (15) and (16) by other appropriate functions. Criteria used in selecting functional forms for W are i) at large separations compared to H, fluctuations and angular dependence of interpartide forces become small, ii) charge assigned and force interpolated should vary smoothly as partides are displaced w.r.t. the mesh, iii) momentum should be (i) conserved computational cost, Criterion leads and to aiv)multipole expansion hierarchy of schemes [8,9], with successively higher order schemes having dipolar, quadrupole, etc., leading error and involving increasing numbers of .

mesh points. Criterion (ii) selects from those schemes satisfying (i) the subset which lead to ‘

smoothly varying forces (and thence smaller collisional effects and better energy conservation). Criterion (iii) demands that the force interpolation function be identical to the charge assignment function to avoid self-forces and possible related numerical instabilities. The cost criterion favours compact assignment schemes of product form. Applying all four criteria one is left with the cost/quality hierachy of product form charge assignment/force interpolation functions,NGP, dC, TSC, PQS [1]. The names associated with the schemes arises from a graphical description of the assignment schemes. Fig. I illustrates this for the case of TSC (Triangular Shaped Cloud) in one dimension. If one imagines that each particle carries a triangle with it, the fraction of its charge assigned to a particular mesh point is equal to the overlap of the triangle with the cell containing the

J. W. Eastwood

/

Particle simulation methods in plasma physics

93

An example of the pitfalls that may arise in choosing time integration schemes without due attention to accuracy, stability and efficiency is

____

x”~’~+v”t+~L4F”

__

1

1.1/5-5—.,

k-X-.,

=

Fig. 1. In TSC assignment and interpolation, the fraction assigned/interpolated to or from a mesh point is equal to the overlap area of the triangle with the mesh point’s cell,

~) ,

~

—

F”’].

=

—

which is identical to leapfrog! Thus, the ‘Beeman’

1(x) ~(x + t. x x ~) 5 \2 W.. ~ S i 1 where by displacement invariance, W’ — —

—~[2F”~’+ 5F” 6

(18)

.

2 =

+

F”’]

This scheme has been favoured in certain circles as a more accurate (!?) alternative to leapfrog for particle calculations. However, by eliminating v from eqs. (18) we find (19) x”~~ 2x” + x”~ F” At2/m

mesh point. Inspection of fig. 1 then reveals

W

~“

I..—1/2.X—--—--—.’

-

(17)

algorithm will give the same results as leapfrog age and four times the amount of computer time! (apart from roundoff), but require twice the stor-

W(x

Lorenz force terms are usually treated implicitly: 5~’~2 v””~”~ qE” + (v”””2 + v””~’2) ~n At m 2 (20)

—

1~,( x)

=

,

—p).

—

2.5. Time integration The large number of particles involved in PM simulations favours simple low order schemes. Stability requires Q At 1 for most explicit schemes independent of their order. Limitations of storage make two timelevel second order schemes the most attractive, where greater accuracy is obtamed by reducing the timestep. Indeed, it is questionable as to whether higher order schemes do introduce any gains where fields and their derivatives are defined numerically. A further factor favouring the lower order schemes is the absence of the potentially troublesome nonphysical (parasitic) roots of higher order schemes. The most widely used integration scheme is the time symmetrical second order accurate leapfrog scheme (eqs. (6) and (7)). Leapfrog is stable for timesteps At < 2/f2, where f2 is the maximum characteristic frequency of the system. It requires simultaneous storage of particle positions and velocity at only one timelevel. Numerical errors arising from finite At appear as phase errors (oscillation frequency too large) without dissipation.

An interesting property of eq. (20) is that for large At, the component perpendicular to the magnetic field reduces to the adiabatic drift approximation (v E X B/B2). A frequency correction factor to compensate for the finite At phase errors in the cyclotron motion is commonly used in eq. (20) [1]. Explicit integration schemes, as outlined here, are limited by the fastest frequencies encountered. Usually, these are frequencies related to the electron mass. Several methods have been devised to overcome this. The first is to operate with reduced mass ratio, and then scale phenomena to the correct ratio. In a similar vein, the electron mass may be taken to zero, so electron dynamics are reduced to an Ohms law (a hybrid fluid/particle scheme). Electron subcycling (i.e. taking many electron integration steps per ion timestep), electron orbit averaging and implicit time integration (through implicit treatment of fields rather than particles) have also been proposed. We refer readers to refs. [2,21—28]for further details. =

.

=

94

J. W. /2astwood

/

Particle simulation methods in plasma phisic.c

2.6. Energy conserving schemes ‘Energy conserving’ schemes [1,2,12—14]are a class of collisionless PM schemes derived using a variational finite element method. They are optimal for the basis functions chosen, Hamiltonian in the limit At—sO, and can be applied to fully electromagnetic systems, awkward geometries and arbitrarily spaced meshes. To use illustrate the method ofplasma discretisation, we again the one-dimensional model. The action integral is a functional of potential 4) and particle coordinate, x: I[4), x]

ft’L

=

dt,

(21)

where the combined field and particle Lagrangian is [151 U

=

f [~mv2

fdft

fdx

+

—

/ ~o I ~ 2

I2

4))

—

(22)

The first integral is over phase space and the second over physical (x) space. The finite element approximation to L is obtained by replacing f by a set of samples, f (cf. eq. (12)) and replacing the potential 4) by trial functions 4)

4)

=

(23)

~4)~(t)J4’~,(x).

p Amplitudes { 4~,} correspond to mesh potential values and l’J’, correspond to charge assignment/ force interpolation functions. The discretised approximation, L to eq. (22) is N 5

I

r

I ~mv~

=

—

q~4)~~(x,)l

Eqs. (25) specify the optimal (in the sense that. it minimises I ~ I) energy conserving scheme once W is chosen. The lowest order conforming choice for W is the CIC scheme (W TI * H) where * denotes convolution. This gives =

(4)~~ —

m~1

=

—

q

~.

I)

H 4)p~l

—

H2

+

.

I 0

I [


(26)

1 ~ ft(x,)~. j (27)

,~

The right hand side of eq. (27) is CIC assignment and the left hand side is the three point fda to the Laplacian. The difference from the momentum conserving PM scheme comes in force interpolation (l.h.s. of eq. (26)), which is NGP centred at cell boundary points rather than mesh points. A consequence of the mixed CIC/NGP scheme is that momentum conservation is lost: an isolated particle placed arbitrarily on a mesh will oscillate in its own force field, but this factor is generally unimportant [14]. 2. 7. Physical properties Particle mesh models are amenable to the same analysis as the differential models from which they derive. Transform space analysis gives an exact description of the errors incurred in the force calculation cycle. Linearised analysis shows how numerical effects change wave dispersion and damping. Kinetic theory for finite sized particle ensemles describes effects of the ‘graininess’ introduced by discretisation of the phase fluid, and a simple stochastic theory accurately describes the loss of energy conservation due to non-conservative numerical errors in the force. Space precludes

i=i[

+

Jdx,[~

dW

1

2 —

p0~~I~ .

P

]

(24)

The approximate action integral, I, arising from

I leads to the Euler—Lagrange equations aI d aI 0; al =

—

~

av,

=

0.

(25)

all but a brief outline here: Details m~ybe found in refs. [1,2,15—17]. For a uniform periodic (or infinite) spatial mesh, the charge assignment step may be recognised as a convolution followed by a sampling and the steps to solve for potential, difference to find fields and interpolate forces may be seen as convolutions. Application of the convolution theorem then reveals that for a density of particle centres. n(x),

J. W. Eastwood

/

Particle simulation methods in plasma physic’s

whose transform is ñ(k), the force field (for q H 1) is

95

1,

putational costs demands as few particles as possible. Long collision times are obtained with wider

(28)

treated bywhereas particles, smallest dispersive practicalproperties particle width are better H(~AD for plasmas).

Successive terms in the r.h.s. of eq. (28) are the transforms of the operations, force interpolation, potential differencing, solving Poisson’s equation and charge assignment, respectively. Comparing eq. (28) with the expression for the correct force harmonic amplitudes

Numerically stable PM calculations generally exhibit slow secular increases of total energy [18]. The dominant contribution to this stochastic heating is the fluctuating part of the force this part corresponds to the n # 0 terms in the alias sum in eq. (28). The size of the fluctuating force is determined by the rate of decay of W, or, in x-space, the smoothness of W: this is one reason why the

=

=

i( k)

~xact

=

=

(

Wb~ i’i’(

—

~

—

k~) ~(ku)).

ik~ñ

(29)

enables the combined effect of all the errors on the forces to be evaluated. Errors in the force direction arise from potential differencing (f~ and the (non-conservative) force fluctuations are due solely from charge assignment. The sum over ii in eq. (28) is known as the alias sum. Aliasing arises because the mesh can only resolve wavelengths lying in the principal (first Brillouin) zone. Shorter wavelengths in the density distribution masquerade as wavelengths in the principal zone. Optimally matching elements of the force calculation using eqs. (28) and (29) can give orders of magnitude improvements in force accuracy without any increase in computational cost. Using eq. (28) in linearised dispersion analysis reveals that the non-physical mode coupling (aliasing) can cause numerical instability for AD ~ H. The consequence of this instability is a rapid increase in total energy, until the plasma temperature is sufficient for AD H, i.e. the mathematically convenient ‘cold’ plasma approximation is inaccessible to the PM simulation experiment, Collisions are analysed by interpreting the PM model as a finite sized particle ensemble. Collisionless behaviour is obtained if the number of particles per particle is large. For plasmas, the collision time is related to the plasma period number of particles per Debye square (in 2D) ND, and particle width w by ~,

I =

ND

w (~) ]. 21

+

(30)

In computations a balance has to be struck between collisional and dispersive properties. Corn-

—

schemes in the NGP—CIC—TSC hierarchy are generally superior to multipole schemes. Typically, energy conservation is 25 times better for CIC than for NGP [19]. 2.8. Optimal PM

The process of setting up an optimal PM simulation model comprises the following steps: a) choose At for stability and accuracy. b) choose computational box to cover longest wavelengths. c) choose W and mesh size to get best quality/cost compromise, where quality is measured by dispersion/conservation compromise. Almost invariably, the area weighting (CIC) function is the best for collisionless plasmas. d) choose N~to get acceptable collisional effects. e) given W for assignment/interpolation, choose a differencing operator D to get commensurate accuracy in the force direction, and choose an optimal G (cf. section 3.3) to offset reshaping effects of W and D for wavelengths of importance, and set G to zero for large wavenumbers to suppress the strong short wavelength alias fluctuations. The outcome is a numerical scheme which identically conserves mass, charge and momentum, and guarantees density positivity. Stability is assured if ~2At <2. Long collision times, accurate wave dispersion and good energy and angular momentum conservation are readily obtained using the Q-minimising optimisation procedure [20].

96

J. W. Eastwood

/

Particle simulation methods in plasma physics

3. Plasma kinetic modelling

integration scheme, giving the timestep loop:

Computer experiments on the statistical mechanics of plasmas measuring correlation functions and thence transport coefficients would require a prohibitively fine mesh for the PM method to handle interparticle forces between neighbouring electrons or ions. Costs also militate against direct summation unless only a small number of partides (‘~500) is involved, because of the N~scaling of operation count for the long range Coulombic forces. In dense plasmas, ionic liquids (and also in stellar clusters, galaxy clusters and molecular structure), the assumption of short range correlations is often poor: It is in such dircumstances that the P3M methods [1,29—32]are appropriate, for they enable the N operations count of the PP method to be reduced to an operations count which is proportional to N~. The trick used in P3M is similar in concept to the Ewald summation method. The interparticle force, f, expressed as the sum of a short range part, f~r and a smoothly varying part, R

a) Start: (x,”, ~,n 1/2; E [1, N~]}; b) PM force calculation: p7 = p,’””2 + F mesh; i ~ [1, Np]; c) PP force calculation: p” 1/2 = ~* + ~i fsr d) Move: x,” ~ = x,” + p,” ~ e) Go to.’ (a).

—

—

f=f sr +R,

(31)

which is non-zero only for small particle separations (r < i~), is computed by direct pairwise summation. If N~is the number of neighbouring particles in the range (= re) of each particle, the operations count for this step is aN~N~. The

The momentum, p, rather than velocity is used in P3M to reduce the computational costs of step (c). Step (b) follows the PM cycle. Cost/quality compromise favours TSC for almost all P3M calculations (cf. CIC for PM). Step (c) requires a fast algorithm for locating the neighbours of each particle: Described in ref. [1] is a chaining mesh and a linked list technique which reduces the search to a negligible portion of the timestep cycle. On vector machines where indirect addressing has a large cost penalty, it becomes more efficient to sort particle coordinates rather than use linked lists. 3.2. Error analysis and optimisation

fsr’

The errors introduced into the force calculation in the PP step are generally negligible. The residual error is that due to the failure of the mesh calculated approximation F(x = x 2 x1 x1) to model exactly the desired smooth part of the force R(x) between two particles at positions x1 and x2. A measure of the error is given by 1 Q = dxif dx F(x; x1) R(x) I 2 (32) —

smoothly varying force, R, is computed using the PM method, with operation count aN~+ /3(N). Generally, a pseudo-spectral treatment of the PM part of the calculation is employed to simplify the minimization of numerical Problems 3M errors. calculations are ini) setto ting uptheeconomical P r of f small and ii) to make cutoff radius make anisotropies and fluctuations in R negligible. These have been solved using the Q-minimising optimisation method [1 15,20 30,31]. 3.1. The timestep loop The timestep loop of P3M differs from that of PM schemes only in the addition of the short range particle—particle (PP) force calculation. Again, leapfrog provides the most convenient time

—

f ~

—

where 0b and ~ are, respectively, the volumes of the computational box and potential mesh cells. Applying the power theorem [31] enables eq.(32) to be written in terms of Fourier transforms W, D and ~ (cf. eq. (28)). For a given choice of charge assignment function, W, eq. (32) reduces to Q = Q(~,b). Minimising Q w.r.t. harmonics G of the Green’s function and over a restricted class of difference operator D gives optimal combinations of elements of the PM part of the force calculation as a function of re.

J. W. Eastwood

/

Particle simulation methods in plasma physics

Subsidiary error measures can also be evaluated to identify which component of the calculation cycle is the dominant error source, thus enabling elements of the calculation to be optimally matched. Indeed, the optimisation may be cornpletely automated, reducing the choice of scheme to a choice of cycle time/particle number/accuracy compromise: This one choice alone is sufficient to determine W, re, G and mesh size. A full account of this procedure, and a 3-D implementation of the P3M algorithm may be found in refs. [1,29]. The operations count, T, for P3M is given by T= aN~+ fiN3 log N3 + yN~N~.

(33)

The first term, aNn, is the count for charge assignment, force interpolation, filling chaining mesh and linked list tables and integrating the equations of motion, the second terms is for solving the field equations (using FFT) and the third term is the count for the short range PP force summation. Constants a, /3 and y are determined by the particular choice of W, D, etc. Substituting for the number of neighbours N~in eq. (33) gives the cost per particle T

N3 log N3 N~

r

3

(34)

Minimising F w.r.t. N~gives the optimal operating point /3N3 log N3

=

y(re/N)3Nc~,

(35)

at which the potential solver time equals the short range force calculation time. In addition, eqs. (34) and (35) give the optimal number of particles per cell

1/3

/ N~\ )oPt =

~ ~‘~‘

L

]

log N3 ]i/2 (36)

and the optimal operations per particle per timestep =

a + 2[y/3re~ log N3] 1/2

(37)

Breakeven for P3M with PP schemes occurs at 300 particles. For 1000 particles, P3 M shows a

speed gain tides!

10, rising to

97

450 for 30000 par-

4. Fluids and magnetofluids Particle methods for fluids and magnetofluids are in many respects similar to those for collisionless plasmas. The incompressible inviscid vortex model [32,33] differs from the PM model of section 2 in that velocity, rather than forces are the relevant vector field. A mathematically identical system arises in the guiding centre model of a plasma [36]. The earliest successful schemes are the PlC models [37—40].Again, this success may be attributed to the physical nature of the errors: mass positivity and conservation of mass, momentum and energy are assured, with numerical viscous heating being properly accounted for in theinenergy equation. Although PlC is not diffusive mass transport, it is in other variables. This is of little consequence in short timescale shock calculations but is a serious drawback in smooth flow studies. Several authors have devised particle fluid codes which retain velocity as a particle, rather than a field, attribute [41—43].These models develop distributions of velocities throughout space, whereas the fluid velocity is a single valued function of position. To prevent multistreaming, numerical viscosity is added to the equations of motion to redistribute momentum locally. It may be argued that the PlC fluid codes which retain velocity as a particle attribute are not cost effective, as they carry extra information which has to be suppressed by artificial viscosity. Even more difficult to justify on a cost basis are the meshless fluid particle methods such as smooth particle hydrodynamics [3,44]and Larson’s method [45]. Like their meshed counterparts, they use the number of particles/particle to estimate the volume per particle, and so require a relatively large numberIf of particles to overcome statistical fluctuations. long range forces are present one is faced with N scaling (cf. section 3). It is for these reasons that I believe these meshless methods are effective only for the special cases (e.g. accretion with voids) for which they were originally devised.

98

J. W. Eastwood

/

Particle simulation methods in plasma physics

Another class of ‘meshless’ particle fluid methods are the free-Lagrangian methods [46]. These are not meshless, but are meshed Lagrangian methods with dynamical reasoning. The ‘particles’ are the nodes of the dynamical mesh. They have more in common with methods such as those described in refs. [47,48] and the above comments on meshless systems do not apply to these. PlC methods have been largely discounted by many researchers as too expensive and inaccurate for fluid and MHD calculations, except in a few special circumstances. However, applying the optimisation techniques developed for Vlasov and Liouville equation solvers (sections 2 and 3) to fluid and MHD equations shows this view to be ill founded. We give the acronymn EPIC to the optimised methods [49]: if the mesh is taken to be principal carrier of information, EPIC are Ephemeral Particle-In-Cell methods, whereas if the particles are the principal information carriers, they are Extended Particle-In-Cell methods. Most of the remainder of section 4 is devoted to the former case. 4.1. Ephemeral Particle-In-Cell (EPIC)

conceptual devices and EPIC schemes for hyperbolic equations appear similar to the characteristic finite element methods [52—54]. 4.1.1. Coordinate projection For low dissipation systems, it is natural to think in terms of the motions of parcels of fluid. This, the Lagrangian view, simplifies the transformation of the governing (almost) hyperbolic equations to particle equations. We assume here that particles move with the fluid velocity, but this is not necessary (cf. section 4.3.3). We recall that if a parcel of fluid initially at x~ moves to x at later time t, so that x=x(x0,

then the infinitesimal volume d’r() and mapped to dT

=

D

-

X0

iS

(39)

dT0.

where =

D1

~

.

EPIC extends the ideas of optimised PM meti. ods to more general hyperbolic and mixed hyperbolic/parabolic equations for fluid flows: prognostic equations are transformed to coordinate projection in space and time, a mesh (or finite element net) is used to describe the displacement field, and minimising errors in principal variables leads to prescriptions for ‘sharpening operators’ [1]. The finite element formulation employed in EPIC is advantageous in the MHD context: it decouples parallel and perpendicular transport and assures non-linear stability in the semidiscrete limit [50 51] Finite element and particle—mesh methods were initially combined in an elegant Haniiltonian formulation by Lewis (section 2.6 and refs. [1,2,12— 14]). In (magneto-)fluid dynamics, the need to use particles as the principal carriers of information becomes less compelling, in which case the partides can become ephemeral. In 1-D, or with fractional timestepping, integrals become easy to do exactly, in which case particles becomes merely

(38)

t),

t0

.

(40)

=

.

D is the displacement gradient matrix with elements D and D J is its determinant. Similarly, the transformation of a vector element of surface area dS() dS is described by ,

—*

dS

=

D J.

dS~•D

(41)

It follows that the density p satisfies mass conservation ~ D (42) =

in the absence of source terms, and the magnetic field B obeys flux conservation D —

—

•

~X0,

—

( )

in the absence of resistivity. Sources s add to e.g. eq. (42) the term D

D Is dt’,

(44)

1)

which we neglect in order to keep the exposition simple, but which will be illustrated in examples below (section 4.4). Eqs. (38), (39) and (41) provide the means for,

J. W. Eastwood

/

Particle simulation methods in plasma physics

respectively, transforming point, volume and surface conservation laws from Eulerian to point Lagrangian form. They can be equally well applied to incompressible flows, reduced MHD, etc. 4.1.2. Discretisation We project eq. (42) onto the set of trial functions { Wk ( x, t)) to give the equations for the amplitudes (p }

tive definite quadratic quantities can be shown to be non-increasing [50,51]. We specialise to case (ii) since arbitrarily moving elements introduce complications but no fundamental differences. The lowest order elements which conform for both advection and diffusion terms are piecewise linear. In 1-D elements of unit width Wk are the CIC assignment functions: Wk(x)

Jwk[P—P

0IDI”] dT= _fWk(dT_0~ (45) where ~ is the finite element approximation to p, i.e. p = ~ + , = p14)1, (4)~} are basis functions and the summation convention is used. Depending on the choice of time-dependence of { Wk), eq. (45) describes (i) a Lagrangian finite element method [48], (ii) an Eulerian EPIC scheme or (iii) a moving mesh EPIC scheme. For (i) (Wk) satisfies Wk(x) = Wk(xo) and (45) becomes (Jwk(x)4)l(x)

dT)P/=

f Wk(xo)P(xo) d~, (46)

while for (ii) awk/at = 0 so

(Jwk4)/1

dT)Pl=

1

mk

W(x 1—

W(x)={0

xk),

lxi,

IxI~1, otherwise.

The right-hand-side of eq. (47) becomes = JWk(x + E)W~(x)dxp0~,

(49)

(50)

where ~ is the displacement of point x0 in the time interval (t, t + At). Representing ~ by piecewise linear functions reduces the integral in eq. (50) to a piecewise cubic given by the overlap integral of “triangle” functions (fig. 2). The mass matrix on the left-hand-side of eq. (47) is tridiagonal with stencil (1/6,2/3,1/6), so it is easily inverted to give p at time t + At. This gives the timestep loop:

(48)

v) go to (i).

to be solved for unknowns { p1 }. mk may be interpreted as the mass associated with node k. For (i) node masses are fixed and the mass matrix Akl changes at each timestep (because the geometry and connectivity of the elements vary), viceversa for (ii), while both Ak! and mk change in (iii): (iii) corresponds to adaptive elements which e.g. align themselves according to the direction of the magnetic field locally [55]. In each case, summing over k shows the scheme is conservative, For the remainder of the paper we shall assume the Ritz—Galerkin approximation i.e. { Wk } and (4)k) are the same functions. Eq. (45) then minimises the mean-square error in p, and in addition non-linear instability is suppressed because posi-

—

mk

(47) Eqs. (46) and (47) may each be written as Akipi =

=

i) ii) iii) iv)

W’~(x)p(x0) d’r0.

99

t:=t+At, find displacements ~, compute node masses {mk} (eq. (50)), solve matrix equation p = A”m to get (p1(t + At)),

4.1.3. Quadrature in space In the timestep loop above, particles enter only as conceptual devices. For two and three dimensional calculations, where a uniform orthogonal mesh is used, this one dimensional scheme can be employed using a fractional timesteping scheme where the displacement is approximated by a Sequence of one dimensional displacements. If fractional timesteping is inappropriate, then the overlap integrals may be evaluated using quadrature. Setting 4~= W, in eq. (47) and using the trapezium rule gives ~ (x,) W,,(x,) p1 = ~ Wk (x, ) p (x0, ) ~.

(~

~)

(51)

~oo

J. W. Lastwood

/

Particle simulation method.,

There are M = 1/3 uniformly spaced points in

ii,

plasma physki,

W

1 (X0)

/‘Nko

each unit width element, where x, = i~(i.e. one quadrature point coincides with each element node for all M). This choice of quadrature assures stability. To evaluate the r.h.s. of eq. (51) (‘assignment’), we set ~(x0) = x(x) to give x0, = x(x,) and = I a~0/a~ I &~,. Quadrature on the l.h.s. changes the tridiagonal stencil for the one dimensional uniform net to ((1 ~2 )/6. (2 + ~2)/3 (1 82)/6). Inspection of the l.h.s. of eq. (51) shows it to be identical to the move and assign steps of PM models. The difference in the Ephemeral PlC schemes is that new particles are generated at each step, where for reasons of stability and economy, they are chosen to arrive at uniformly spaced positions at the end of each timestep.

N _____________

1

NN

____

~

N

k-I~I~

0

k__i_tk_

—

~,

—

Fig. 2. Contributions to A51(~) arc given by the overlap integral of ~F~(x,,) with the displaced and deformed triangle function W5(x~

—

[36—40]follows because, in effect, it lumps the mass matrix. This large diffusion is absent in EPIC. Time splitting in PlC cn be interpreted as a consequence of applying one term (explicit) quadrature in time to terms such as eq. (44). Generalisation to higher order then follows trivially.

4.2. ‘EPIC and other schemes 4.3. A comparison of schemes

4.2.1. FCT Steps (iii) and (iv) in the above timestep loop (section 4.1.2) are like flux-corrected transport [56] in that mk may be interpreted as a low order approximation to new densities pk(t + At) and inverting the mass matrix as antidiffusion. Explicitly, eq. (48) with linear elements becomes (52) 2Pk + Ph_i). PA

mA

=

~~(Pk~i

which has the form of a diffusion equation with coefficient 1/6. —

4.2.2. Godunov The “Lagrange then remap” procedure has been shown to be effective in finite difference schemes right-hand side of eq. (50) amounts to a mapping for compressible flows [57,58]. We remark that the of a density distribution obtained from a Lagrangian procedure back onto a fixed mesh in one dimension (cf. fig. 2). PlC The right-hand side of eq. (51) is the ‘move and assign’ operation in the jargon of PM. It may be interpreted as: take particle I at x 0, = x1 move it to x1 and assign mass 8m1 = p(x0~)~0, to node k according to the assignment function Wk. The diffusive nature of conventional PlC 4.2.3.

431 Linear analysis

Linear analysis gives the variation with wavenumber of numerical dispersion and dissipation. The best schemes generally have small damping for wavenumbers with small phase errors, and large dissipation errorstoare large.for a Applying the where usual phase procedure EPIC problem with uniform advection (~ = constant) and node spacing, we derive first the linearised equation, a convolution, I( p

—

q)pq(t + At)

=

J( p

—

q

—

flpq(t)’

(53)

where for CIC functions WA ,,

=

{(4 (2 0.

— —

2z2[2 I z ])/6. Z )~/6, —

:1 ~ I, I ~ Z ~ 2. otherwise. (54)

Fourier transforming eq. (53) gives the dispersion relation A=e~’~=J(k,fl/i(k,O),

(55)

—

which is to be compared with the exact result w At

=

k~.

(56)

J. W. Eastwood / Particle simulation methods in plasma physics

Linear stability is assured since A lies in the unit circle for all ~ The dispersion of EPIC with linear basis functions (CIC) is compared with that of Lax— Wendroff and upwinding in fig. 3a and with FCT [56] and the high order Godunov scheme PPM [57] in fig. 3b. All curves are drawn for a Courant number C = 0.25: larger timesteps are more flattering to EPIC, but we have to choose a value for which the other schemes are stable. Lax—Wendroff has large phase errors and small damping so we expect (and find cf. below) spunous wiggles and incorrect propagation speeds. Upwind has smaller phase errors but much larger dissipation. FCT has about the same phase errors as EPIC, but greater damping at shorter wavelengths: PPM has slightly larger phase errors and dissipation intermediate between EPIC and FCT. For wavelengths greater than 4 node spacings (wavenumbers <1/2), fig. 3b shows there is little to choose between EPIC, FCT and PPM on the basis of linear analysis. Fig. 4 shows phase errors and damping for EPIC with quadrature for 1—4 particles per ele-

Exact

—

EPIC

—-—

Lax -Wendroft Upwind

a Phase

——

‘

—

02~

—

——

.—

~

I

0

0.5

0.2

_.~

~—

~

~

05

1 0

Wavenumber

.

-0.2

1,0

Damping

,

-

—-

—

Exact

FCT

‘

—EPIC

——‘—PPM

b

ment. The additional dissipation introduced by the “graininess” of the particle discretisation uniformly decreases as the number of particles in-

0.2

integration caseerrors exceptdiffer for Mlittle = 1. from the exact crease. Phase

0

4.3.2. Kinematic test problem . A simple test problem to compare schemes is given by transporting a Gaussian hump through a unit length periodic box using a velocity field u(x) that provides both compression and expansion. Take u=2—sin2~x, N=100 nodes in a spatial period, C = 0.25, and a Gaussian density hump placed initially with its peak at x = —0.2 (p(t = 0) = exp(—[x + 1/5]2/a2), a = 0.1). The compression reduces the effective a of the Gaussian to approximately 3 node spacings at smallest (proportionately increasing its height to conserve total mass). The compressed p Fourier transforms to a Gaussian, the amplitude of which is <2 x iO’”°of peak value for wavelengths less than two mesh-spacings, i.e. p can always be resolved to within round-off error (or 4-byte arithmetic) by the chosen mesh.

—— —

101

Phase

——— —~

-

— —-

i

—— —

-

——

c

-

0

0.2

0.5 -

1.0

Damping

7 I

~

C

________________________________________ I

0.5

1.0

Wavenumber 02

-

Fig. 3. (a) and (b) real (phase) and imaginary (damping) parts of the frequency versus wavenumber for uniform advection. The curves are identified in the legends.

Fig. 5 shows the outcome for the five schemes of section 4.3.1. The broken curves are the initial conditions and the solid ones the result after

102

J. W. Easiwood Phase

/

Particle simulation methods in plasma physics

-+v-=R~’-~-—~() at ax ax2 on the unit periodic interval with initial conditions 1.’ = sin 2’rx. We consider three cases i) conservative ii) non-conservative and iii) non-conser-

.-

0.2—

0.1

‘~

i

-

—

I

0.5 02

Damping

1 1

~

C

0.5 Wavenumber

vative with quadrature andofineq. each treat the diffusion term in onspace, the r.h.s. (57)case by the new timelevel implicit approximation. The projection equations to be solved to obtain node values v~at the new timelevel from v 0(x0) at the old timelevel are for case i) =

v0 At/2,

-0.2 Fig. 4. As fig. 3, for EPIC with M —1, 2, 3. 4 particles per element. The curves labelled ~ are for exact spatial integration.

I WW+R” j

A

/

aw.aw At___~_~~L dxv ax ax f

/

(58)

=fWk(x)vodxo

and for case ii) transport through a distance of 3 periods (about 2000 steps): broken and solid should be identical in the absence of error. Fig. 5a and b show that Lax—Wendroff and (conservative) upwinding perform badly. Fig. Sc illustrates the profile steepening that is characteristic of FCT. From fig. Sd PPM has more accurate slopes than FCT, but fig. Se reveals that EPIC clearly gives the best approximation. The cpu times taken were 0.21 for (a), 0.26 for (b), 0.44 for (c) and 1.39 for (d), all measured relative to (e) EPIC. Thus EPIC is of comparable speed to the schemes it outperforms. Fig. Sf shows that EPIC with quadrature retains the good phase error properties of the exact spatial integration, but shows increased dissipation. The good stability properties of EPIC allow large timesteps to be taken: One consequence of this is that the number of steps (and cost and total numerical dissipation) to reach a given physical time can be made smaller than comparable explicit finite difference schemes. Fig. 6 illustrates this for m = 4 particles (quadrature points) per element. Cost are inversely proportional to Courant number, 4.3.3. Nonlinear test problem We apply the linear basis function EPIC scheme to Burgers’ equation

~ = v,~At,

awk aw WkW+R~’ =

At~

f Wk(x)vO dx.

dxjv1 (59)

-‘

Case iii) is as case ii), but with ~ replaced by x and integrals replaced by trapezium rule (“partide”) approximations. Cases i) and ii) give almost identical results. For large Reynolds number R, the sinusoid rapidly becomes a “sawtooth”, then decays on the resistive timescale. The significant spatial harmonic of highest frequency in the problem has wavenumber <1.SR. Fig. 7a shows that using case ii) when R = 100, solutions v are indeed accurately resolved and smooth. At higher R wiggles develop, but note that unlike many finite difference schemes which are non-linearly unstable at such R, due to aliasing error [50,51], EPIC remains stable. The appearance of oscillations is the Gibbs phenomenon, caused by neglect of significant sub-grid-scale harmonic content. Appropriate sub-grid-scale physics models can be used with EPIC (as with any other scheme) to remove these wiggles. The extra. dissipation in case iii) suppresses some of the wiggles. More interestingly, the robustness of the particle advection allows much

J. W. Eastwood / Particle simulation methods in plasma physics

103

I—.’ I

I

a t

II

I

/

I

/

b

I -

-

05 .

II II I

-

III %

I

j

0.5

-

I

/ ___________

I

—,

___________

O.~

___________

-0.5

r\

0

05

1’~

Ii

i

i

I

I

I

I

‘

—

II -0.5

0

0.5

-0.5

I

0

-

0.5

/ I i

I

f

/

_______________________ _____________

-0.5

0

_______________________ _____________

0.5

—— _____________

~0.5

0

_______________________ I

0.5

Fig. 5. Initial (dashed curve) and final (solid curve) profiles of a Gaussian convected by a sinusoidally varying velocity field as described in section 4.3.2, for (a) Lax—Wendroff, (b) upwind, (c) FCT, (d) PPM and linear basis function EPIC schemes (e) with exact spatial integration and (f) with M = 8 particles per element.

104

J. W. Ea.,twood

/

Particle simulation methods in plasma physics C=0.5

II

C=1

2~-~-~-

/A 0.1

0 5

__

0.5

_________________ I/ II

C

I

0

05

=

2

_______________

C

____________________=

4

______________

Fig. 6. The Gaussian convected by a sinusoidally varying field for linear basis function EPIC with M = 4 and Courant numbers as shown on the curves.

larger timesteps to be used. Fig. 8 shows the solution remaining accurate as well as stable for Courant numbers up to four! R

=

100

R

=

1

000

Fig. 8. Sequences showing evolution of solutions to Burgers’ equation using case (iii) scheme of section 4.3.3. with = 4 4, V M = 100, particles per element, predictor x = i’ ~t, R i0 and Courant numbers as shown. Curves are drawn at intervals of 40/C timesteps.

5. Final remarks

R

=

10 000

B

100

000

Fig. 7. Sequences showing overall reduction in amplitude as time increases of the solutions to Burgers’ equation described in section 4.3.3. Linear basis function EPIC is used with predictor f r’ ~t; C = 0.5, N = 100, and the curves are drawn at intervals of 100 ~t. Reynolds numbers. R. are as indicated,

integral part of the arsenal of the computational Particle methods will undoutedly remain an physicist. The particle—mesh method for collisionless systems has now a well developed ‘finite sized particle’ theory. Wave dispersion and growth, drag and diffusion coefficients and conservation properties show agreement between theory and experiment. Well defined methods for analysis and optimisation of PM algorithms now exist; It is no longer justifiable to use NGP and simple finite differences for Poisson’s equation except for gross velocity space instabilities. The P3M algorithms are a further result of the analysis techniques developed for PM models. P3M opens a whole new area of previously inaccessible study; Reducing the N~body force sum from an e2(N) to an ~2(N~) operation count process has made possible the numerical study of

J. W. L’astwood

/ Particle simulation methods in plasma physics

large scale correlation in systems with long range forces. The fluid PlC scheme, despite being only ‘zeroth’ order and lacking the theoretical backing that PM models have, has remained in use for over two decades: Its greatest strength is in handung particularly awkward supersonic flows. The relative ease with which PlC handles multimaterial boundaries, voids and open surfces far outweighs the difficulties of numerical diffusion and relatively high computational cost. The meshless fluid particle models seem less likely to overcome computational cost barriers in the case where long range forces are present. EPIC, which extends the PM optinMsation techniques to fluids and magnetofluids with dissipation, outperforms state-of-the-art finite difference schemes. EPIC schemes are consistent, given unconditionally stable kinetics and the Galerkin variants can be shown to be free from non-linear instabilities. Optimal accuracy arises because EPIC causes errors in principal variables to be minimised. Indeed, the key ideas behind FCT, PPM and particle methods are embraced by the EPIC formalism. A great attraction of particle models is the natural way in which they incorporate the conservation laws: The usually difficult advection terms appear as movement of discrete ‘parcels’ of conserved quantities around the computational box. The close parallels between particles and the underlying corpuscular nature of the continuum system being modelled mean that numerical errors are physical in character. Robustness, simplicity and accuracy have made particle methods a natural choice for collisionless phase fluids. The more recent developments outlined here are likely to lead to particles becoming a natural choice for all fluid and magnetofluid flow modelling where advection is a dominant process.

Acknowledgement I would like to thank Dr. Wayne Arter for his help in preparing figures of section 4 and for his comments on the manuscript.

105

References [1] R.W. Hockney and J.W. Eastwood, Computer Simulation using Particles (McGraw-Hill, New York, 1981). [2] C.K. Birdsall and A.B. Langdon, Plasma Physics via Cornputer Simulation (McGraw-Hill, New York, 1985). [3] J.J. Monaghan, Comput. Phys. Rep. 3(1985) 73. [4] V.K. Decyk, Comput. Phys. Rep. 4 (1986) 245. [5] R.W. Hockney, Meth. Comput. Phys. 9 (1970) 135. [61P.C. Clemmow and J.P. Dougherty, Electrodynamics of Particles and Plasmas (Addison-Wesley, Reading, 1969). [7] J M Hammersley and D.C. Handscombe, Monte Carlo Methods (Methuen, London, 1964). [81 W.L. Kruer, J.M. Dawson and B. Rosen, J. Comput. Phys. 13(1973)342. [9] J.W. Eastwood and R.W. Hockney, J. Comput. Phys. 16 (1979) 342. [10] D Beeman, 1. Comput. Phys. 20 (1976) 130. [11] J.P. Boris, Proc. 4th Conf. Nurn. Sim. Plasma (NRL, Washington, 1970) p. 3167. [12] [13] [14] [15] [16]

HR. Lewis, J. Comput. Phys. 6 (1970) 136. HR. Lewis, Meth. Comput. Phys. 9 (1970) 309. A.B. Langdon, 3. Comput. Phys. 12 (1973) 247. J.W. Eastwood, J. Comput. Phys. 18 (1975) 1. A.B. Langdon, J. Comput. Phys. 6 (1970) 247.

[17] C.K. Birdsall, A.B. Langdon and H. Okuda, Meth. Cornput. Phys. 9 (1970) 241. [18] R.W. Hockney, J. Comput. Phys. 8(1971)19. [19] A. Pieravi C.K. Birdsall, UCB/ERL College of and Engineering, Univ.Report of Calif. BerkeleyM78/32, (1978); Proc. 8th Conf. Nurn. Sim. Plasma, paper P0-9, Monterey, CA (1978). [20] J.W. Eastwood, Comput. Meth. Class. Quant. Phys. (Advance, London, 1976) p. 196. [211 R.J. Mason, J. Comput. Phys. 41(1981) 233. [22] B.I. Cohen, R P. Freis and V Thomas, J. Comput. Phys 45 (1982) 345. [23] B.I. Cohen, A.B. Langdon and A. Friedman, 3. Comput. Phys. 46 (1982) 15. [24] A.B. Langdon, B.I. Cohen and A. Friedman. J. Comput. Phys. 51 (1983) 107. [25] A.B. Langdon and D.C. Barnes, in: Multiple Time Scales. eds. lU. Brackhill and B.I. Cohen (Academic Press, New York, 1985) p. 336. [26] B.I. Cohen. in: Multiple Time Scales, eds. lU. Brackhill and B.C. Cohen (Academic Press, New York, 1985) p. 331. [27] R.J. Mason, in: Multiple Time Scales, eds. J.U. Brackhill and B.I. Cohen (Academic Press, New York, 1985) p. 233. [28] D.C. Barnes, T. Kamirnura, TN. Le Boeuf and T. Tajima, J. Comput. Phys. 52 (1983) 480. [29] R.W. Hockney, S.P. Goel and 3W. Eastwood, Chern. Phys. Lett. 21 (1973) 589. [30] lW. Eastwood, Comput. Meth. Class. Quant. Phys. (Advance, London, 1976) p. 206. [311 lW. Eastwood, R.W. Hockney and D.N. Lawrence, Cornput. Phys. Cornmun. 19 (1980) 215. [32] G. Efstathiou and J.W. Eastwood, Mon. Not. Roy. Astron. Soc. 194 (1981) 503.

106

J. W. Eastwood

/

Particle simulation methods in plasma physic’s

[33] R. Bracewell, The Fourier Transform and its Application (McGraw-Hill, New York, 1965). [34] i.P. Christiansen, I. Comput. Phys. 13(1973)363. [35] J.P. Christiansen and NB. Zabusky, J. Fluid Mech. 61 (1973) 219. [36] W.W. Lee and H. Okuda, I. Cornput. Phys. 3 (1978) 139. [37] F.H. Harlow, Meth. Comput. Phys. 3(1964)319. [38] A.A. Amsden, Los Alarnos Scientific Report LA-3466 (1966). [39] F.H. Harlow, A.A. Amsden and JR. Nix, I. Comput. Phys. 20 (1976) 119. [40] T.L. Cook, RB. Dernuth and F.H. Harlow, J. Cornput. Phys. 41(1981) 51. [41] R.L. McCrory, R.L. Morse and K.A. Taggart, Nuci. Sci. Eng. 64 (1977) 163. [42] V.M. Marder, Math. Cornput. 29 (1975) 434. [43] iN. Leboeuf, T. Tajirna and J.M. Dawson, I. Comput. Phys. 31(1979) 379. [44] R.A. Gingold and ii. Monaghan, Mon. Not. Roy. Astron. Soc. 181 (1977) 375. [45] RB. Larson, J. Comput. Phys. 27 (1981) 39~7, [46] R.A. Clark, in: Num. Meth. Fluid Dyn. II, eds. K.W. Morton and M.J. Baines (Oxford Univ. Press, Oxford, 1986) p. 255.

[47] Mi. Fritts and J.P. Boris, I. Comput. Phys. 31(1979)173. [48] lW. Eastwood, J. Cornput. Phys. (1986) in press. [49] lW. Eastwood and W. Arter, in: Num. Meth. Fluid Dyn. II, eds. K.W. Morton and Ml. Baines (Oxford Univ. Press, Oxford, 1986) p. 581. [50] lW. Eastwood and W. Arter, Culham Preprint CLM-P763 (1985); Phys. Rev. Lett. (submitted). [51] W. Arter and J.W. Eastwood, Transp. Theory Stat. Phys. (1986) in press. [52] G. Labadie, J.P. Benque and B. Latteux, in: Proc. 2nd Intern. Conf. Nurn. Meth. Laminar Turb. Flow, eds. C. Taylor and BA. Schrefler (Pineridge Press, Swansea, 1981) p. 681. [53] K.W. Morton, in: Num. Meth. Fluid Dyn., eds. K.W. Morton and Mi. Baines (Academic Press, London, 1982) p. 1 [54] 0. Pironneau, Numer. Math. 38 (1982) 309. [55] lW. Eastwood and Lee-Hsaio, Culham Laboratory Report CLM-R253 (HMSO, 1985). [56] J.P. Boris and DL. Book, I. Comput. Phys. 20 (1976) 397. [57] P. Collela and P.R. Woodward, J. Comput. Phys. 54 (1984) 174. [58] B. van Leer, J. Comput. Phys. 32 (1979) 101.