J. Plasma Fusion Res. SERIES, Vol. 9 (2010)

Development of full particle-in-cell simulation code withadaptive mesh refinement technique
Toseo MORITAKA 1,3), Masanori NUNAMI 2,3) and Hideyuki USUl 1,3)
1) Graduate school of Engineering,

Kobe University, Kobe 657-8501, Japan
Toki 509-5292, Japan
3) Japan Science and Technology Agency(JST), CREST, Kawaguchi, 332-0012, Japan
2) National Institute of Fusion Science,(Received: 20 November 2009 / Accepted: 21 March 2010)

A numerical scheme for interface between different hierarchy levels with different grid sizes andtime step intervals is designed for a full particle-in-cell simulation code with adaptive mesh refinement
technique. Physical information in the interface region is exchanged recursively over all hierarchy levels.
The exchange processes are performed at the synchronizing time steps with the target hierarchy level byusing integer and half integer time steps in leap flog method. In addition, we employ overlap grids at theinterface region in which electromagnetic field and particles are advanced under both grid conditions
independently. Test simulations of light and Langmuir wave propagations indicate that the physicalinformation can exchange smoothly through the interface region by means of the combination of these
two schemes.
Keywords: multi-scale phenomena, particle-in-cell simulation, adaptive mesh refinement

1. Introduction

dynamically applied according to the local conditions.
Fine grids suitable to the local high density region are
applied only there and other region is simulated by using moderate size grids. Therefore, increment of thenumerical cost due to the localized region is not serious by adopting the AMR technique. The AMR technique has been rarely employed in PIC codes except
for a few examples[4] while that is actively employedin fluid codes.
We have been developed a full PIC simulation
code with the AMR technique in order to consider
a next generation interplanetary flight system called"Magneto plasma sail(MPS)". The propulsive force
is obtained from electromagnetic interaction between
solar wind and artificial magnetosphere around thespacecraft[5] . The magnetosphere is created by super
conducting coils inside the spacecraft and inflated byplasma injection largely enough to interact the solar
wind sufficiently. Spatial scale of the original magnetosphere is attempted to be roughly comparable tothe electron gyration scale (a few kilo meters or less)and the injection plasma density could be 109 - 10 times
larger than the solar wind density. Therefore the solar
wind interaction and the magnetic inflation for MPS
are typical examples of the multi-scale phenomena to
be demonstrated by means of the AMR PIC simulation.
A difficulty of the AMR technique is exchangeof physical information among different hierarchy domains with different grid sizes and time step intervals.
One can not reference the neighbor grids in the otherhierarchies directly for finite difference scheme. Therefore an additional numerical scheme for the interface
is needed for the appropriate information exchange.

Full particle-in-cell (PIC) simulation[1] is a efficient simulation method to investigate microscopic
plasma phenomena. Electron as well as ion kinetic
properties are included in full PIC simulation andwidely employed to demonstrate various phenomenain space and laboratory plasmas.
A difficulty of the PIC simulation is their considerable numerical costs. Spatial grid size and time step
interval are determined by Debye length and Courant
condition. These scales are usually much smaller thantypical scales of plasma phenomena to be investigated.
This difficulty becomes critical in multi-scale phenomena controlled by multiple plasma processes on different typical scales.
The difference of these typical scales comes fromtheir physical origins and physical conditions. The
physical origins include electron gyration motion, ion
inertia, global external forces and so on. Dynamical interaction between geomagnetospheric convection driven by the global interaction with solar windand magnetotail reconnection triggered by electron kinetic effects is a typical example of such multi-scale
phenomena[2]. In addition, these typical scales are
depend on the local conditions. When a simulation
domain includes localized high density regions, for example, the grid size defined there by Debye length
is too small to treat the surrounding environment inreasonable numerical costs. Laser induced plasma expansion is an example of such cases[3] .
Adaptive mesh refinement(AMR) technique is a
useful method to investigate such multi-scale phenomena. In the AMR simulation, various size grids areauthor's e-mail: moritaka.toseo@cs.kobe-u.ac.jp586©2010 by The Japan Society of Plasma
Science and Nuclear Fusion Research

T. Moritaka et al., Development of Full Particle-in-Cell Simulation Code with Adaptive Mesh Refinement Techniqueis advanced by referring to the neighbor grid data.
However information exchange for FTDT calculation
between different hierarchy levels is asynchronous dueto the distinction of the time step intervals. The asynchronous information exchange could cause numerical noise at the interface because the interface region
acts as an inappropriate boundary (for instance, fixedboundary for the domain with shorter time step interval).
Figure 1 shows a schematic picture of the calculation order and corresponding time steps in the hierarchy domains of levels Land L+ 1. Here the first
integer time step and the next integer time step inlevel L are denoted as nand n+l, respectively. Field
data can exchange at the integer time step in the lower
level L (steps nand n+l in Fig.l). Field data of lower
level L at the half integer step (step n+l/2) also send
to the higher level L+ 1 as the field data at integer
time step by separating the calculation for electricfield (Eq.2) into two steps i.e., E(n) --+ E(n + 1/2)
and E(n + 1/2) --+ E(n + 1). In this way, electromagnetic field is advanced appropriately referring the fielddata in the other hierarchy domain except at the half
integer time steps in the higher level (steps n+ 1/4 andn+3/4).
In addition, overlap grids are introduced at theinterface region as shown in Fig.2. The electromagnetic field in the overlap grids are advanced independently in both hierarchy domains by using their own
grid sizes and time step intervals. The numerical error caused at the edge of the overlap grids without thesynchronized information exchange (steps n+l/4 andn+3/4) does not act on the external region until thenext synchronizing time step (steps n+l/2 and n+l,
respectively). The incorrect information in the overlap region is wasted and updated by the informationof the other hierarchy domains in the next informationexchange.
The AMR simulation with many hierarchy levels
includes various time step intervals. The calculations
should be made in appropriate order over all hierarchy
levels. The calculations in each level are separated
into two parts. The separated processes for level L
are given as follows,

The outline of this manuscript is as follows; The
fundamental equations to be solved in each hierarchy domain are given in Sec.2. Numerical schemes
about electromagnetic field and super particles forthe interface between the different hierarchy domains
are explained in Sec.3 and Sec.4, respectively. These
schemes are evaluated by some test simulations in eachsection. Conclusion and future work are described inSec.5.

2. Fundamental
structureequations

and

data

Basic model is a three dimensional electromagnetic explicit particle-in-cell simulation. Electromagnetic field and super particles are advanced by using
leap flog method[6] for Maxwell equations
1~t
B(n + "2) = B(n) - 2(\7 x E(n))

(1)

c2

1
E(n + 1) = E(n) + 47r ~t(\7 x B(n + "2))1+J(n +"2)

(2)

1
~t
B(n + 1) = B(n + "2) - 2(\7 x E(n + 1))(3)

and Boris-Bunemann method[l] for Newton-Lorentz
equation, respectively, where n, ~t and c are time
step, its interval and light speed. The charge conservation scheme[7] is employed for the local electrostatic
correction in the complicated grid alignment in the
AMR simulation.
The AMR simulation is structured hierarchically.
Grid size and time step interval are defined according to the hierarchy levels, where high and low levels
correspond to the fine and coarse grid systems, respectively. As the system evolves, some grids in ahierarchy domain (Level L) are removed according to
the local physical conditions and corresponding grids
in the higher hierarchy domain (Level L+ 1) are newly
produced as the divided grids.
The dynamical grid system is organized by using fully threaded tree(FTT) structure[8]. Instead ofthree dimensional data matrix, the grid data is given
by one dimensional array of structures which contains
physical value at the grid and some pointers. The
pointers represent the relationship to the other grids,
such as 'neighbor' for the surrounding grids, 'parent'
for the original coarse grids in the lower hierarchy and'child' for the refined grids in the higher hierarchy.
These relationships are varied dynamically by reconnecting the pointers.

First half; B(n) --+ B(n + 1/2), E(n) --+ E(n +1/2), data exchange of B(n + 1/2), E(n + 1/2)
with level L+l and E(n + 1/2) --+ E(n + 1)
Second half; B(n + 1/2) --+ B(n + 1) and data exchange of B(n + 1), E(n + 1) with levels L-l and
L+1.
Then, the calculation order is defined recursively as
follows,

3. Interface for electromagnetic field
3.1 Numerical scheme+-- from the first half calculation of level L-l1 calculate first half of level L

Maxwell equations are solved by finite differencetime domain (FTDT) method. Electromagnetic field

587

T. Moritaka et al., Development of Full Particle-in-Cell Simulation Code with Adaptive Mesh Refinement Technique

Chart oftime e\ olution

n

n+1/4

n+ll2

I Level L I
n+3/4

actual calculation arean+l

Ioverlap region I

time step _1-1- - - - - - - - - - - - - + 1 - - .
interpolation

Level L

1····.····.····I Level L+1 I

I~

I

I

··.····.····,.......;---,.-+-~-+-~T""""""1

L...L..............--...---I........................."""'"--...........

?verlap region--+------+----..-------+AU2~ ....i....L...i.

I

I

actual calculation area

Fig. 2 Schematic of overlap region in the interface of two
different hierarchy domains of level L(blue boxes)
and level L+1(orange boxes).

Level L+1

local magnitude of the magnetic field. Coarse (level0), middle(level 1) and fine(1evel 2) grid regions arecorrespond to the regions with strong(B z ~ 0.70Bo),
moderate(0.31Bo ~ B z ~ 0.70Bo) and weak(B z :::;
0.31Bo) magnetic field respectively, where B o is theinitial amplitude of the wave. In other words, thegrids are refined as the magnetic field gradient becomes steep. Theoretically, the initial wave propagates at light speed in x direction without any modification of the wave structure. The refined areas should
also move with the propagating waves.
Figure 3 shows the magnetic field profiles at t =2.66Po, 3.33Po and 4.84Po observed in the AMR simulation using above-described interface schemes, where
Po is the initial wave period. Blue, purple and red
points denote observed magnetic field in the hierarchy
domains of level 0, levelland level 2, respectively.
The wave profile is remarkably smooth at the interface region. The wave propagates in x direction without any numerical noise.
Numerical error in the FTDT calculation aboutwave propagation appears as a phase error. The phaseerror is estimated by comparing with the theoreticalwave structure. If the wave propagates at light speed
without any modifications, the equation

Synchronization

Fig. 1 Schematic of time step calculation using the leap
flog method in two different hierarchy domains oflevel L(top) and level L+1(bottom). The calculations flow rightward with the time evolution in thisfigure. Green, orange and purple arrows correspondto the calculations of B(n) ---+ B(n + 1/2),E(n) ---+
E(n + 1) and B(n + 1/2) ---+ B(n + 1), respectively. Red boxes indicate synchronized data exchange processes between level L and level L+1.

2 --+ goto first half calculation of level L+ 1+-- from the first half calculation of level L+ 1
3 calculate second half of level L--+ goto second half calculation of level L-1.
To advance the electromagnetic field from step n tostep n+1 in the AMR simulation with hierarchy levels
0,1 .. ·M, we start with the first half calculation forlevel M (from step n to step n+ 1/ 2M +1 ) and continue
calculation in the above-described order until after thesecond half calculation for level O. For example, thecalculation order for the case of M=2 is obtained as,(4)step n : (2-1)--+(2-2)--+(1-1)--+(2-1)--+(2-2)--+(1-2)--+(0-1)--+(2-1)--+(2-2)--+(1-1)--+(2-1)--+(2-2)--+(1-2)--+(0-2) : step n+l

is fulfilled. Figure 4 shows time evolutions of theBz(x = ct, t) normalized by B o until t = 2.5Po in
three different numerical schemes for the interface region. Those are only synchronization at half integer
time step(blue line), only introduction of overlap grids(green line) and both of these schemes (red line). It
is shown that both of these schemes are needed to
avoid the growth of the phase error modifying thewave structure.
Figure 5 is a comparison with the uniform gridcase until t = 5.0Po. Red line represents the profile ofelectric field resulting from the AMR case using bothof the numerical schemes. Green and blue lines represent that resulting from the coarse and fine uniform

where (m-1) and (m-2) denote the first and second
half processes for hierarchy level m, respectively.

3.2

Test simulations

Test simulations of light wave propagation are
performed to evaluate the numerical schemes at theinterface. At initial, cosine waves of electric fieldE y(x) and magnetic field B z (x) in vacuum are assumed. Periodic boundary condition is employed inall directions. Spatial grids are refined according to588

T. Moritaka et al., Development of Full Particle-in-Cell Simulation Code with Adaptive Mesh Refinement Technique
wave profile

phase error

1.5 ,..---...,...------r"--,.....--...,...------r"--,.....------,

(I)

(IT )

1.03

(Ill )
1.02

Lv.O
05

LV.l

1.01
0

co

o

~
co

Lv.2

0

~

];

N

co

Lv.l

-0.5
-1

0.98

-1.5 '--_.....1.-_--'-_ _'--_.....1.-_--'-_ _' - - - - '
o
0.1
0.2
0.3
0.4
0.5
0.6
0.7

0.5

x (normalized)

1.5

2.5

normalized time

Fig. 3 Magnetic field profile in AMR case at 2.66(1),3.33(II) and 4.84(II1) periods. Grid are refined according to the local magnitude of the magnetic field
Bz. Coarse grids (level 0) are refined when B z :s 0.7
and the fine grids (level 1) are refined again (level2) when B z :s 0.31. Red line and Green lines represent initial profile and theoretical wave profiles ineach time step, respectively

Fig. 5 Time evolutions of the magnetic field at the sampling position of x = ct. Red line represents theresult of AMR case. Green and blue lines represent
the result of uniform coarse and fine grid cases, respectively.

4. Interface region for super particles
4.1 Numerical scheme

phase error
1.12Super particles in PIC code can move throughthe boundary of each hierarchy domains. In addition,neighbor grids are referred to calculate the electromagnetic force which acts to the super particles and current density due to the particle motion because the super particles have their own sizes of a few grids. Thus
consistent exchange of physical information throughthe interface region is important in these calculation
processes as in FTDT calculation of electromagnetic
field.
Basic concept is same as the overlap region forelectromagnetic field. That is, the physical information to be referred until the next synchronization stepis ensured in the overlap grids. First, super particles
near the interface are copied to the overlap region of
the other hierarchy domains as shown in Fig.6. By using the copied particles, particle motion and resulting
current density are calculated in each hierarchy domain independently according to their own grid sizes
and time step intervals. In leap flog method, particle
velocity is defined at half integer time steps. Thus velocity of the copied particles is estimated by temporalinterpolation in order to synchronize the half integer
time step in the destination hierarchy level.
One calculation cycle for particle motion and current density in the higher level is made without thesynchronization with the lower hierarchy level as in
the calculations for electromagnetic field. The overlap
region should contain all grids which would be referring in these calculations. Displacement of particles
during one time step is at most one grid because of
the Courant condition, and surrounding one or two
grids are referred to estimate the electromagnetic force
acting to the particle and current density due to the

1.08
1.04
0

.;e
u
];
N

co

0.99

Lv.O

0.96
0.92

normalized time

Fig. 4 Time evolutions of the magnetic field at the sampling position of x = ct. Blue, green and red lines
represent the results of the case with the synchronization at the half integer step, with the overlap grids and with both numerical schemes, respectively.

grid cases, respectively. The grid size in the formercase is the same as that of the hierarchy level 0 in
the AMR case. In the uniform grid cases, the phaseerror increases according to the grid size but the fluctuating range is stable. The phase error in the AMRcase evolves in a similar manner to that in the coarse
uniform grid case at initial because sampling points
x = ct belong to the coarse grid area of level O. The
phase error decrease in time under the influence of
the accurate wave propagation in the fine grid area.
It is concluded that accuracy for the calculation ofelectromagnetic field with leap flog method could beenhanced without additional numerical noise by using
AMR technique with the combined interface schemeof the synchronization and the overlap region.

589

T. Moritaka et al., Development of Full Particle-in-Cell Simulation Code with Adaptive Mesh Refinement Techniquecalculation a~ea of Level X

AMR vs uni[orm(Lv.O)

overlap region

1.5

cn)

(I)

D

D

cm)

CN)

0.5

Super particles are copied
to the overlap region

0

~
CiI

overlap region

calcu'lation area of Leve!1 X+1-0.5

Particle motion and current
density are calculated withthe copied particles

-1

Lv.O

Lv.l

-1.. 5

10

Remaining particles at the

20

30

40

50

60

x/Rd

overlap region are removed.

Fig. 7 Electric field profiles at tW pe =(I)O.9, (II)2.7,(III)4.5 and (IV)6.3 in the case of partially refined grid(red) and uniform coarse grid(blue). Red
vertical lines represent the boundary between thecoarse(left) and refined(right) grid regions in the
nonuniform grid case. Green line denotes initial
profile.

Fig. 6 Schematic picture of copy process of super particles
in the overlap region. Top two figures represent the
copy process between level X and level X +1. Red
and blue circles indicate super particle in the hierarchy domains of level X and level X +1, respectively.
Translucid circles are copied particles. Bottom two
figures indicate particle motions including that of
the copied particles.

lines in Fig. 7). The unrefined grid size is the same as
the grid size in the uniform grid case. The electric field
are oscillating with a period of 27rw;/ as predicted theoretically. The profiles are almost identical in coarse
grid area. Numerical error at the interface region andtheir influences on the other region are hardly foundin the simulation results. Figure 8 is for a comparison
with the uniform fine grid case. The grid size is thesame as the refined grid in the nonuniform grid case.
The observed electric field profiles are almost same in
the fine grid area. Again numerical error is not foundat the interface region.
These results indicate that physical information
about plasmas and electromagnetic field can exchange
consistently through the interface region by using thedeveloped numerical scheme. However, there is slightdifference between the wave amplitudes in the fine andcoarse grid cases. Damping of the electric wave isfound to be obvious in the coarse grids. This difference
might come from the number of super particles in eachgrid rather than the interface scheme.
In the present scheme, total number of super particles is constant thus number of super particles per
gird decreases in the fine grid. The partcles number
per grid can affect wave-particle interaction in phasespace and resulting Landau damping. One could solve
this problem by using particle revision[9] in copy process of super particles at the interface region.

displacement. Therefore width of the overlap region
for the copied particles is needed to be more than two
grids. These calculation is combined to the numerical scheme for electromagnetic field in the overlap region. Further overlap grids are required to calculate
the time evolution of the electromagnetic field aroundto the particle.

4.2

Test simulations

Test simulation of Langmuir wave propagation is
performed to examine the numerical scheme for electromagnetic field and super particles at the interface.
Initial density, velocity and electric field profiles correspond to the Langmuir wave with a wave number k
as,neno + nsin(kx)

27rn

(5)

Ex

Tcos(kx)(6)

Vex

wpeEo . (k x,)---SIn(7)

27rnoewhere E, V e and n e denote electric field, electron net
velocity and electron density. Temperature and ion
density profiles are uniform at initial. The wave number and the fluctuation of electron density are set as
k = O.IIX;1 and n = O.lno, respectively. The simulation results in the uniform grid and partially refined
grid cases are compared.
Figure 7 shows the observed electric field from
uniform grid case(blue lines) and partially refinedcase(red lines) at tw pe = 0.90(1), 2.7(11), 4.5(111) and6.3(IV). In the latter case, grids in the range ofxl AD=33.0-57.6 are refined (between two red vertical

5. Summary
We have developed a numerical scheme for theinterface between different hierarchy levels in a fullparticle-in-cell simulation code with adaptive mesh refinement technique. Each hierarchy level has theirown grid size and time interval. Grids in the hier-590

T. Moritaka et al., Development of Full Particle-in-Cell Simulation Code with Adaptive Mesh Refinement Technique

AMR vs unijol'm(Lv.l)

to poor number of super particles per refined grid. Inorder to examine wave-particle interaction accurately,
particle division and corresponding modification of the
shape factor are needed in the next step.

1. 5 r - - - - - - - . - - - - , - - - - , . - - , - - - - - , - - - - - - , - - - - - , - - - ,

(I)

(II)

(ill)

(N)

0.5
0

Acknowledgement~
CiI

This study is supported by the JST jCREST. The
computation in the present study was performed with
the KDK system of Research Institute for Sustainable
Humanosphere (RISH) at Kyoto university.

-0.5

-1

-1. 5

10

20

30

40

50

60

[1] C.K. Birsall, and A.B. Langdon, Plasma Physics via
Computer Simulation, (Adam Hilger, New York,1991).
[2] A.T.Y. Lui, IEEE Trans. Plasma Sci 28, 1854 (2000).
[3] M. Nunami and K. Nishihara, Journal of Plasma andFusion Research SERIES 8, 815 (2009).
[4] K. Fujimoto, Physics of Plasmas 13, 072904 (2006).
[5] R.M. Winglee, et aI, Journal of Geophysical Research
105, 67 (2000).
[6] K.S. Yee, IEEE Trans. Antennas Prop 14, 302 (1966).
[7] T.Zh.Esirkepov, Computer Physics Communications
135, 144 (2001).
[8] A.M. Khokhlov, Journal of Computational Physics
143, 519 (1998).
[9] G. Lapenta, Journal of Computational Physics 181,317 (2002).

x/Rd

Fig. 8 Electric field profiles at tW pe =(1)0.9, (11)2.7,(111)4.5 and (IV)6.3 in the case of partially refined grid(red) and uniform fine grid(blue). Red
vertical lines represent the boundary between thecoarse(left) and refined(right) grid regions in the
nonuniform grid case. Green line denotes initial
profile.

archy domains are removed or produced dynamicallyaccording to local physical conditions. Consistent information exchange among the different hierarchy domains is needed to avoid numerical error by adopting
adaptive mesh refinement technique.
The information at the interface region is exchanged with the other hierarchy domains by using
overlap grids. The electromagnetic field and superparticles are copied to the overlap grids in the otherhierarchy levels and advanced independently by using
their own grid systems. The information exchange
is performed at the synchronizing time step between
different time step intervals by using integer and half
integer time steps in leap flog method. Incorrect information obtained at another time step remains inside
the overlap region and renewed at the next synchronization time step. The calculation order is definedrecursively over all hierarchy levels.
Test simulations are performed to evaluate the validity of these numerical schemes. In the simulationof light wave propagation, phase error due to finite
difference time domain calculation is reduced by using AMR technique and the interface scheme without any additional numerical noise. Langmuir waves
are observed to propagate smoothly through the interface between unrefined and refined grid regions.
These results suggest physical information can exchange with the other hierarchy domains consistently
through the interface region by using the developed
numerical scheme.
Super particles are distributed to the refined grids
at the overlap region in the present copy process. Current density and number density are correctly estimated by the distributed particles. However the effective resolution of the velocity space is reduced due

591