Rigorous charge conservation for local electromagnetic field solvers

manyhuntingUrban and Civil

Nov 16, 2013 (4 years and 8 months ago)


Computer Physics Communications 69(1992) 306316 Computer Physics
North-Holland Communications
Rigorous charge conservation for local electromagnetic
field solvers
John Villasenor
Jet Propulsion Laboratory,4800 Oak (Jim e Dr.,Pasadena,(A 9/109,USA
Oscar Buneman
Department of Electrical Lngmeering,5tanjord Unoersity,,Stanjord,(A 94305-4055.(ISA
Received 23 January 1991:in final lorni 12 June 1991
in this paper we present a rigorous method for finding the electric field in two-andone-half and three dimensional
electromagnetic plasma simulations without resorting to coniputationally expensive transforms.A finite grid interpretation
of the divergence equation ~J =  ~ip/~)i is offered which allows the current density and thus ness local electric and
magnetic field strengths to he determined directly from knowledge of charge motion.
1.Introduction ence form,already express such local updating
Field equations are most commonly solved by
transform methods in simulations.The ready  = V xE,(1)
availability of efficient transform codes encour-
ages one to use such spectral methods.So does
our custom to describe,interpret,and understand = V X B  J.(2)
field phenomena in terms of waves.
However,transform methods are global:all (using units such that  = 1,,u = I.and c= 1).By
of the field information,near and far,contributes staggering B-component and E-component data
to each single field harmonic.Computers of the suitably in both space and time (leapfrogging).
immediate future demand minimization of data these equations provide new data from old,using
paths and therefore call for local methods.Even present values only immediately to the east,west,
with hypercube topology,ideally suited to FFT north,and south (see fig.5),and above and below
implementation,the physical length of single data in three dimensions [1.71.
links will ultimitely become the bottleneck.Meth- By contrast,the two divergence equations,
ods which allow one to update fields purely from
local data should therefore receive renewed at- V E = p.(3)
tention and it is gratifying to note that two of
Maxwells equations,when cast into finite differ- V B = 0,(4)
00l0-4655/92/$05.00 © 1992  Elsevier Science Publishers By.All rights reserved
J.Villasenor,0.Buneman/Rigorous charge conservation for local EMfield solvers 307
require distant information (spatial boundary on gently and (2) that enough steps are taken to
conditions) for their solution,particularly in the let the initial radiation escape across the compu-
case of static conditions.Coupled with V x E= 0 tational boundary.(Non-reflecting boundary con-
and V X B = J,they then lead to elliptic equations ditions,such as Lindmans [51become important
for the components (or for the associated poten- here.)
tials).Transform methods,or other global,direct In a two-dimensional Nx N grid this number
non-iterative poisson solvers would therefore of steps is of the order N with a computational
seem inevitable,effort of order N2 per step.The transform effort
Fortunately,the dynamic problem of field evo- is of the order N2 log2 N.However,since the
lution can be solved by means of the evolution- local Maxwell algorithm is so much simpler than
ary,local equations (1) and (2) alone,after impos- even the most elegant transform algorithm,the
ing the divergence equations (3) and (4) only as transient,local method wins over the global
initial conditions.One readily checks that V.B transform method even on present-day computers
remains zero if it was so initially and that this is on which the length of data paths is not yet an
rigorously true also for the simple space-and- issue.
time-centered finite difference versions of a/at,In the first proposal for solving the field prob-
curl and div that go with the data arrangement lem locally (Buneman [1]) strict charge conserva-
shown in fig.5.One also checks that V E re- tion was to be achieved with a technique called
mains p by virtue of the conservation condition zero-order current weighting.In this method,a
charge crossing a cell boundary furnishes an im-
 (5\pulse of current proportional to the value of the
   charge.While charge is conserved,the current
discontinuities introduced generate large amounts
on the charge density.of noise.The nature of this noise is discussed in
What we shall be concerned with in this paper Birdsall and Langdon [2].Morse and Nielson [3]
is to make sure that the finite-difference imple- discuss first-order current weighting,where the
mentation of the charge conservation law (5) is charge is taken to be evenly distributed across a
also consistent with the finite difference version square of unit size and the particle mover calcu-
of a/at,div and curl as applied to E and B.In lates currents by breaking a general translation
other words,if we generate E time step after into two orthogonal moves of magnitude z~xand
time step from the finite-difference versions of ~y.Most recently,Marder [4] described a simula-
(1) and (2) and calculate V.E,the resulting p tion method which allows eq.(5) to be approxi-
should satisfy (5) exactly in this finite-difference mately satisfied,and the error,or pseudo-cur-
version.rent is kept small.Both our particle mover and
Enforcing the initial conditions on the diver- those discussed by Morse and Nielson are able to
gences of Eand B at time t = 0 now remains as a rigorously satsify eq.(5);the crucial difference is
once-only task.Given the initial p (and some that we do not break up the charge motion into
boundary conditions),one can either use a trans- orthogonal moves.
form or other non-local method for this first step,
or one can let p evolve from zero to its initial
state transiently in a certain number of prelimi- 2.Calculating the fluxes and currents
nary steps,starting with source-free E and B that
satisfy the boundary conditions.Such E and B Instead of updating the longitudinal and
are in many cases known analytically.transverse parts of the E-field separately,as is
Experiments performed on a personal corn- often done in field solvers which employ spatial
puter by one of us have shown that the static field Fourier transforms,we update E entirely from
is approached satisfactorily in transience pro- Maxwells aE/at = V x B  J.In the process the
vided (1) that the charge distribution p is switched curl B term will,of course,only cause changes in
308 J.Villasenor,0.Buneman/Rigorous charge conservation for local EM field solvers
the transverse part of E.However,when J has a
longitudinal part  which by virtue of V J =
ap/at means when the charge density changes
 we get changes in the longitudinal part of E,______________
i e its flux by Just the change in p as it should
be In other words there is no need to calculate ______________ __________
the longitudinal part of E separately except at
the very begtnntng (at time t = 0) This supposes
that the ftnite difference form of current depost
tion is consistent with the continuity equation,i.e.__________________
that the flux of J represents the change in the
charge distribution p rigorously.
The transverse part of E is affected also by the
magnetic field B.The record of B has to be
updated at every step by the transverse part of E,
in accordance with Maxwells aB/at =  V x E.
In doing this,one accounts fully for both TE
modes as well as TM modes which are excited by
Fig.1.Particle-in-cell.A charge with square cross-section
the motion of the charges.. .
moves in a simulation space divided into cells also having
We start the simulation by adopting the part!- square cross-section.
cle-in-cell method  creating a grid consisting of
squares of unit size,and placing on it unit square
charges.The charges may be located at arbitrary the moves for purposes of illustration;in the
places,and can be assigned either positive or actual simulation the time step must satisfy the
negative values.The total amount of charge is Courant condition,~t < ~x/~/i for a square
assumed to be uniformly distributed over its sur- mesh in two dimensions.This limits the charge
face,and as fig.1 shows,each charge will in displacements to less than ~x/~ per step.
general contribute to the charge density in four The simplest and most common type of charge
cells.The fractions contributed to each of the motion involves four boundaries.The action of
four cells are proportional to the rectangular the charge sweeping across the cell boundaries
areas of the charge square intercepted by the gives rise to currents J~,~ in v2 as mdi-
cells.This is area weighting (see ref.[2]).cated in fig.2.
With each step of the simulation,each charge Defining the local origin for any charge to
is instructed to move a certain distance in a he the intersection of cell boundaries nearest to
certain direction.Since current density is in the the charge center at the start of its motion and
discrete case represented by motion of charge the coordinates x and y to be the location of the
into or out of a cell,it can be seen by the cell center relative to this origin,we can then
divergence equation (5) that each boundary will write the currents (assuming unit charge) as:
be swept over by an area of the square charge I
that exactly corresponds to the current in a given j51 = ~X(1 y  (6)
direction into or out of that cell.Although most = ~x(~ +y + by),(7)
moves will affect the current on only four bound-
aries,motion of a single charge can in general 1vt = x  (8)
cause currents across four,seven,or ten bound- j.,= +x + ~©~x).(9)
aries and it is thus necessary to account for all of  - -
these possibilities in the particle mover.In the These equations give the amount of charge cross-
following figures showing how these cases are ing a cell boundary.For instance,for the ~nt
treated,we have exagerrated the magnitude of equation,the depth of charge moved is ~y.the
J.Villasenor,0.Buneman/Rigorous charge conservation for local EMfield solvers 309
_______ ~/___ IJy2
_______ _______ _______
__________ /____
I ~
4 J Jx2 J~a
______ j ~
______________ _____________ ______________ _________ 1J~ 2
J 91 1
I~ 12
-3 -~
~Y2 ~
Boundaries experiencing
current flow K
Loc& Origins for sections 1 ~nd 2 of
4 Loc& Origin
Fig.2.Four-boundary case.In the simplest,most common
type of move,motion of the charge will only create a current Boundaries Experiencing
across four cell boundaries.A move as shown will create the current flow
four fluxes J~1,J~2,J~and J52 as given in eqs.(6)(9).The Fig.3.Seven-boundary case.A charge can also move in a way
coordinates (x,y) describing the location of the charge center which affects the flux across seven boundaries.The total move
at the start of the move are measured relative to the local of Ax,Ay can be treated as a four-boundary move of Ax1,
origin. Ay1 followed by another four-boundary move of Ax2,Ay2.
See eqs.(10)(15) in text.
width is ~  x at the start of the move and complex case is the seven-boundary move,shown
decreases linearly to ~ x  z~xat the end;the below in fig.3.
average width,which is relevant for linear rno- The seven boundary case can be treated as two
tion,is ~  x  ~x.Note that the main compu- four-boundary moves  the first with the charge
tational effort is four multiplications,no more center starting at x,y and moving a distance
than in the area weighting procedure which is ~,~y1 and the second with the charge moving
necessary when the longitudinal electric field is from x~,Yi a distance.~x2,~y2.There are
calculated from charge densities.For the four- actually four cases of seven boundary moves,
boundary case (but not for the seven- and ten- depending on the direction of charge motion.
boundary cases discussed below) our currents,Equations for the case where the right edge of
given in eqs.(6)(9),are the same as those calcu- the charge comes to rest with x > 1 are given
lated by Morse and Nielson.In general charges below:
may move in a way that they will affect the
current on more than four boundaries.A more ~~x1= 0.5 X,(10)
310 J.Villasenor,0.Buneman/Rigorous charge conservation for local EM field sohi ers
~Yi=(~Y/~X) ix,,(II)
xi = 0.5,(12)
4 J>~2
__________ Jy2
~x2=~x~xi,(14) __________.I
______ I
~y2 =  ay,.(15) 
9 Jx2
The currents contributed by the first part of the  _________
move,from x,y to x1,Yi are indicated in fig.3,/7 4J~2
by ~ r2 ~ and J~.2.The second portion of j52,_.V
the move starts at x 1,y 1 and contributes the I ~iu~1~
currents J~~j2,J~ and J~2.Note that the  i I
boundary connecting the local origins referenced
in the first and second portions of the move, ___________
.9 __j
respectively,receives two current contributions.~xi xl
In unusual cases,a charge can influence the
current on ten boundaries.Just as a seven bound -_____________ _____________ _____________
ary move could be analyzed in terms of two
four-boundary moves,a ten-boundary move can
be decomposed into three distinct four-boundary ~2 Y
moves.There are eight possible cases of a ten- 2
boundary move,depending as in the seven ~ ~i~
boundary case on the exact path which the charge I ~ ~ Ag
follows.For a charge moving as shown in fig.4,
we have the following moves:Ax
Move Startingposition Ax.Ay Fig.4.Ten-boundary case.The most complex move will affect
x.v Ax.A I current across ten boundaries.This move is treated as three
2 x1,y1 Ax,,Ay successive four-boundary moves,as described by eqs.(18)
x2,~ Ax~,Av~ (27).
move denoted by J and the contributions from
~xi =0.5 x,(18) the third move byJ.
~yi=(~y/~x) ~xi,(19)
x~= 0.5,(20)
3.Field updating
~ Yi = 0.5  y  ~ yt,(22) In a two-and-one-half dimensional simulation
(one in which currents and fields in the z direc-
= (~x/~y)~y2,(23) tion are calculated but variations a/az are as-
= ~  0.5,(24) sumed to be zero) the curl equations (1) and (2),
can be expressed as:
hB~ aE.hE.
(26) =--~--~,(28)
cit dy hx
~y3=~y~yi~y2.(27) bE~ hB~
Again,the currents./,~.,,~ J~,,and J~2 are =  J~. (29)
shown,with the contributions from the second
J.Villasenor,0.Buneman/Rigorous charge conservation for local EMfield solvers 311
hE These finite-difference equations can also be in-
=  -.~  (30) terpreted as statements of the integral form of
Maxwells equations.For instance,eq.(32) states
Using the staggered grid and compass directions that the change of E-flux Out of the east face of a
as given in fig.5 and refs.[2,7] the above equa- cube erected over the square in fig.5 equals the
tions are implemented in the following finite-dif- circulation of B around the face,minus the cur-
ference form,where the ratio ôt/~x relates the rent through the face.Once the updated E and
grid size and time scale,and must be chosen to B fields have been calculated,the locations of the
satisfy the Courant condition:particles must be incremented.The choice of
method for determining the force due to fields
B~~ B~td located on cell boundaries on particles which are
at arbitrary locations is extremely important.If
= ~_((E~otth  E~0u1th) (Eeast  E~~t)),the self-force (the electric force a particle ex-
x 31~ erts on itself) is not zero,the grid will contain a
- - potential well at each cell center.Several meth-
ods of interpolating field values to charge loca-
E~W E~= ( Bfb~~~ Bs0~~th) ~J (32) tions have been tried,including twisted plane and
linear interpolations as well as area weighting.It
was found that area weighting (just like that used
rnew uoid i east west\~ J f33\
 =  ~  )  s ~ ) to determine charge density in fig.1) sucessfully
5(1-Id) 6~
(11,j) (iI,)) (i,~) (i,j)
x.x ~(1,J),,p(l,J) ~.~ x
\s/f ~ 
vvesL East
5(1~~~-i) E~1~°J~~°
Fig.5.Standard cell at location (i,J) showing field locations and indexing conventions.Currents J and electric fields E lie in the
plane of the paper and are calculated at cell edges with positive direction indicated by the arrows;the magnetic field B is
perpendicular to the plane of the paper and is given at cell corners.Charge density p and electric Dotential 4 are associated with
the cell center.
312 J,Villasenor,0.Buneman/Rigorous charge conservation for local EM field solvers
eliminates self-force [2].The weights must be lnitielize fields,
applied to field components which have been perticle positions
averaged to the field centers,end velocities
4.Implementation Move perticles,obt6in
boundery current
contributions (see Fig.7)
A simulation program implementing this rigor-
ous charge conservation was written for a grid
size of 128 x 128.An equal number of positive I
I Updete B end E using
and negative charge are placed at random loca- I equetions (3 l)-(33)
tions throughout the grid,and the initial potential ______________________
is calculated using two-dimensional Hartley trans-
forms [6] to solve Poissons equation.When mov -__________________________
ing each particle,the program first uses its initial Teet Ax end ~y for eech cherge
and final positions to determine if the move will by eree weighting field velues
effect only four boundaries.If,as is usually the
case it is a four-boundary move,eqs.(6)(9) are
applied and the program proceeds to the next.,...
Fig.6.Simulation flow chart.General flow of a simulation
charge.If the move is more complicated,a see- program showing initialization,following by repeated passes
ond routine checks for a seven-boundary case.If through particle mover and field updating.
necessary,a ten-boundary routine is called.Al-
though the number of operations involved in ad-
vancing a particle increases with the complexity
of the move,it should be noted that the majority one step because of the initial curl-free E.After
of moves will call the simple four-boundary case,eqs.(32) and (33) had been used to advance the
requiring only two shifts,four multiplies,and six electric field to its new state,the divergence of
additions.The relative frequency of occurrence the new E-field was then calculated as
of four-,seven-,and ten-boundary cases obviously E~ + E~0nh  E~0h1O1)/~x and compared with
depends on the velocity distribution of the simu- the p deduced using the actual new positions of
lation,but even in a system where many particles the particles.It was found that the two methods
approach the velocity of light ten-boundary moves agreed to within roundoff,confirming that local
will remain unusual.field updating can be accomplished at great time
After all the particles have been moved and savings and without any loss of accuracy.
their contributions to current density calculated,
the program updates the fields with finite-dif-
ference equations (31 )(33) and then finds the 5.Extension to three dimensions
force on each particle through the area-weighting
technique described earlier.It should be empha- In three dimensions,area weighting be-
sized that the E field generated at each stage of comes volume weighting.This can be inter-
the simulation using this technique is identical to preted as making each particle into a uniformily
that produced by area weighting;no sacrifice in charged cube of the same size as a cell,with the
accuracy is made.Figure 6 shows the simulation center at the nominal particle position.In gen-
flow chart,and the actual particle mover is pre- eral,the particle cube will occupy parts of eight
sented in fig.7.neighboring cells.The cell mesh will cut the cube
In our testing,the initial E-field was created up into eight rectangular blocks.Figure 8 shows a
from a charge distribution and the initial B-field particle cube penetrating into the cell nearest the
was zero,and as expected,B remained zero after viewer.
J.Villasenor,0.Buneman/Rigorous charge conservation for local EMfield solvers 313
The cells are indexed by the integral position Table I
coordinates x = i,y =j,z = k of their centers.Fraction of the particle in each of the eight cells
Given a particle (i.e.the center of the particle Weightingfraction Cell
cube) at location x = i + ~,y =j + ~,z = k + z,(1 f)(1  s~)(1 C) i,j,k
where ~,ij,~ lie between 0 and 1,table 1 shows ~ (i ~)(l C) i +i,j,k
what fractions of the particle lie in each of the (1f) ai (1 C) i,1 + 1,k
eight cells.(tf)(1~)C ~,j,k +1
(1f)aiC i,j+1,k+1
These are the eight weights which must be ~(1 C + 1,j,k + 1
applied to data recorded at the cell centers (such ~~(1 C) i + 1,j + 1,k
as charge density or electric potential) when de- ~~ ~ +1,j+1,k + 1
positing quantities into the arrays or when inter-
polating from array data.Volume weighting in
this manner amounts to trilinear interpolation:particle on the four x-facing boundaries are (1 
linear in x times linear in y times linear in z.~ (1  ~),(1  ~ ~,i~(1  h),~ ~,and similarly
Each particle cube straddles twelve cell faces,for the y-facing and z-facing boundaries.
four for each of the three orientations:x-facing,Consider a particle which moves straight from
y-facing and z-facing.The areas covered by the (i + ~,I + ~,k + ~) to (i + ~2 ~ + 112~k + ~2)
Four Boundary
eZ~~FourBoundary~> Routine  particle
Yes starts at (x,y)
FourBoundary FourBoundary
c~ Seven-Boundary T2~- Routine  particle Routine  particle
Yes starts at (x,y) starts at (x 1,y 1)
FourBoundary FourBoundary Four Boundary
Routine  perticle Routine  particle Routine  particle
sterts at (x,y) starts at (x 1,y 1 )
starts at (~,y 2)
Fig.7.Particle mover.Flow graph of particle mover involving either one,two,or three calls to the four.boundary routine
depending on the complexity of the move.
314 J.Villasenor,0.Buneman/Rigorous charge conseri ation for local EMfield soli ers
linearly with time,covering a displacement ~Xx=,,,,>x ~ ______________
~2~I ~Y=~12~Th ~z=~2~i in the time
between t =  ~- and t = + ~.Its midway posi-
tion,reached at t = 0.is given by
The total flux transported into the cell indexed ~
i + 1 j + I k + I across its x facing boundary z~/,J~
shared with the cell indexed i,j + I,k + 1 is ~~///Y,~~
= f~2(~+ t ~y)(~+ t ~z) ~x dt
 1/2
= ~x~+ ~x ~y ~z/l2.(35)
This is the contribution to J~at location j +  Fig.8.Three-dimensional version of the particle-in-cell
j + I,k + 1.The other three contributions are:
~x(1   ~x ~y ~z/l2 i.e.the difference between the particle fractions
at + I k + 1,protruding into the cell before and after the
- move.This,then,confirms rigorous charge con-
 ~)  ~x ~y ~z/l2 servation.
at i + ~,I + 1,k,In 3D field solvers the data on E1,E~.E~,B~,
-  B.,B are staggered both in space and in time
~x( 1  ~)(1  ~) + ~x ~ ~z/l2 [11.E1 is on record midway between cell centers
at i + ~- I k.(36) in the x-direction (as it would be if it were
2 obtained by simple differencing of potentials
The four contributions to J~and J~are obtained recorded at the cell centers),i.e.at locations i +
from eq.(36) by the cyclic rotation I,k.This E~is (but for the constant e,,) also the
   electric flux through the face centered at that
i,~x,~ ~y,~ ~k,~z,~ ~x,~ location.(See fig.5 for the 2D version of this
(37) method.) Here we make the approximation that
variations within a cell dimension can be treated
While the first term in each contribution is a as linear,so that central values are avarages.This
plausible generalisation to three dimensions,from flux,according to Maxwell,is incremented in a
what was obtained in eqs.(6)(9) for two dimen- unit time step by the circulation of B/~x11around
sions the ~x ~ ~z/12 terms are new..the cell face and decremented by the current
Addition of the three fluxes into the cell in- through the cell face,exactly the current which
dexed i + 1,1 + 1,k + I yields we have calculated.
-- - The circulation of B is obtained from B~,and
~x~+ ~ + ~ + ~x ~y ~z/4 B values recorded in the middles of the edges of
= (~+rz~x)(~+ ~y)(~+ ~z) the face.B data are kept on the complementary
- - - mesh to that for E data,the mesh obtained by
(i ~x)(~  ~y)(~ ~z),(38) interchanging the roles of centers and corners of
J.Villasenor,0.Buneman/Rigorous charge conservation for local EMfield solvers 315
our original mesh.The face centers of the corn- occupied by any particle one needs to know in
plementary mesh are the edge centers of the which cell of the complementary mesh the parti-
original mesh and vice versa.Thus the data for cle (i.e.the center of the particle cube) is located.
the circulation of E,required to update B fluxes,The displacements ~,sj,~ are measured from a
are also just available where needed.corner of the complementary mesh which is a
One important advantage of our method of center of the original mesh.If,during a move,a
field solving is that no separate arrays are re- particle remains in the same complementary cell,
quired for recording charge densities or current its ~,s~,and ~ will remain between 0 and 1,and
densities.The currents created by the particles the same eight cells of the original mesh will be
can be added directly to the E field values in the occupied by the particle.Currents will flow only
E update step.If,at certain stages of a run,a through the twelve faces that separate those eight
charge density record is desired for diagnostic cells.
purposes,it can be obtained by a single sweep If,however,a particle leaves its complemen-
over the E data,adding the fluxes which enter tary cell during its move,new cells of the original
each cell through its six faces.The scheme was mesh will be affected and currents will flow
tested in three dimensions in just this manner,through more than 12 faces.In two dimensions
checking that the change in charge density due to we identified the corresponding situation as
the movement of some particles was as it should seven-boundary and ten-boundary cases.In
be.Agreement was exact to round-off.three dimensions we have found it expedient not
The complementary mesh has another signifi- to calculate exactly how many more cell faces will
cance:In order to determine which eight cells are receive currents.Instead,we have automated the
Subroutine Subroutine Subroutine
xsplit.Y~Pit zSplit
5:::?n0 ~:::~
2nd particle 2nd particle 2nd particle
Subroutine to
deposit currents
on cell faces
Fig.9.Flowchart showing how crossings of the complementary mesh boundaries are handled.
3 I 6 J.Villasenor,0.Buneman/Rigorous charge eonseri ation for local EM field so/lees
particle splitting procedure illustrated in the flow aries,and the majority of them will be simple
chart of fig.7 by passing each particle through four-boundary cases.Local current densities due
four nested subroutines,as shown in fig.9.to each move can be determined by application of
The creation of pairs is done in each case by eqs.(6)(9) and the electric field then can be
the procedure shown in eqs.(10)(15),with a found by direct application of Maxwells equa-
third co-ordinate to be included.tions instead of through the Poisson solutions
used hitherto.
Using the techniques outlined in this paper,it
is possible to run an entire simulation without This work was supported in part by NASA
resorting to Poissons equation while still rigor- under contract number NAS8-35350 and by a
ously satisfying the divergence equation.Initial- National Science Foundation Fellowship.
ization can be accomplished either by letting the
potential evolve using the method described in
the introduction or by using transforms,and all References
steps thereafter rely on current information cal-
culated directly from charge motion.This rigor- Ill C).Buneman and W.B.Pardo.Relativistic Plasmas (Be-
ous use of local-only information greatly simpli- njamin,New York,1968) p.205.
12] C.K.Birdsall and A.B.Landon,Plasma Physics via Com-
fies updating of fields and is ideally suited to..
puter Simulation (Mc-Graw-l-lill.New York.1985).p.362.
parallel processors.As a consequence of the 13] R.L.Morse and C.W.Nielson.Phys.Fluids 14 (1971) 830.
Courant condition,it will always be possible to 14] B.Marder,J.Comput.Phys.68 (1987) 48.
describe the motion of an arbitrary number of 151 E.Lindman.J.Comput.Phys.IS (1975) 66.
particles over any amount of time as the superpo- 16] RN.Bracewell.The Hartley Iransform (Oxtord Univ.
Press,New York,1986).
sition of the appropriate number of moves.In a......
~7i 0.Buneman.Numerical simulation ot Eields,in:Wave
non-relativistic setting in two dimensions,only a Phenomena.L.Lam and H.C.Morris,eds (Springer,
small percentage of moves will affect ten bound- Berlin,1988),fig.11th.