AN ADAPTIVE, HIGH-ORDER PHASE-SPACE REMAPPING FOR
THE TWO-DIMENSIONAL VLASOV-POISSON EQUATIONS

arXiv:1204.6512v1 [math.NA] 29 Apr 2012

BEI WANG

∗,

GREG MILLER† , AND PHIL COLELLA‡

Abstract. The numerical solution of high dimensional Vlasov equation is usually performed by
particle-in-cell (PIC) methods. However, due to the well-known numerical noise, it is challenging to
use PIC methods to get a precise description of the distribution function in phase space. To control
the numerical error, we introduce an adaptive phase-space remapping which regularizes the particle
distribution by periodically reconstructing the distribution function on a hierarchy of phase-space
grids with high-order interpolations. The positivity of the distribution function can be preserved by
using a local redistribution technique. The method has been successfully applied to a set of classical
plasma problems in one dimension. In this paper, we present the algorithm for the two dimensional
Vlasov-Poisson equations. An efficient Poisson solver with infinite domain boundary conditions is
used. The parallel scalability of the algorithm on massively parallel computers will be discussed.
Key words. Particle-in-cell (PIC) Methods, Adaptive Mesh Refinement, Phase-space Remapping, Numerical Noise, Vlasov-Poisson equation
AMS subject classifications.

1. Introduction. The Vlasov equation describes the dynamics of a species of
charged particles under electromagnetic fields. In the electrostatic case, the normalized equation reads
∂f
s
+ v · ∇x f + (−1) (E + E e ) · ∇v f = 0,
∂t

(1.1)

where f (x, v, t) is the distribution function of the species in phase space (x, v) ∈
Rd × Rd with d = 1, 2, 3. s is 0 for positive charges and is 1 otherwise. E and E e
denote the self-consistent and the external electric field, respectively. This equation
is the simplest model to study collisionless plasmas and beam propagation which is
of importance to controlled thermonuclear fusion and accelerator modeling.
The Vlasov equation is a nonlinear hyperbolic equation in phase space so methods
of solution can be guided by the well-established numerical analysis of classical partial
differential equations. Accordingly, grid methods in fluid dynamics, such as transform
methods, finite-volume methods, and semi-Lagrangian methods, can be employed.
Operator splitting was successfully applied to the solution of the Vlasov equation by
Cheng and Knorr [4] in 1970s. It reduces the solution of the multi-dimensional Vlasov
equation to a set of one-dimensional advection problems, and therefore has become
a widely used technique. Even with these well-established algorithms, performing
high-dimensional simulations using grid methods is still a challenging task. The issue
is the computational time and memory cost in dealing with the whole six dimensional
phase space. With the advances of supercomputer, grids methods have achieved
large development in the last decade. In semi-Lagrangian methods, Sonnendrucker
et al. [30] introduced the cubic spline method. Nakamura and Yabe [26] introduced
the cubic interpolated propagation method. In finite-volume methods, Fijalkow [14]
∗ Department of Applied Science, University of California, 1 Shields Avenue, Davis CA 95616.
Present Address: Princeton Institute for Computational Science and Engineering, Princeton University, Princeton, NJ 08540. E-mail: beiwang@princeton.edu
† Department of Chemical Engineering and Materials Science, 1 Shields Avenue, Davis CA 95616.
‡ Applied Numerical Algorithms Group, Lawrence Berkeley National Laboratory, MS 50A-1148,
1 Cyclotron Road, Berkeley, CA 94720.

1

2
presented the flux balance method. A similar idea is used in a high-order finite-volume
method based on mapped coordinate by Colella, Dorr and Hittinger [6]. Filbet,
Sonnendrucker, and Bertrand [16] proposed the positive and flux conservative scheme
using the idea of limiter.
A more widely used approach for the solution of the Vlasov equation is PIC methods [19, 1]. In PIC methods, the particles, a Lagrangian discretization of the distribution function, follow trajectories computed from the characteristic curves given by the
Vlasov equation, whereas the self-consistent fields are calculated on a grid. Since the
methods employ the fundamental equations without much approximation, it allows us
to observe most of the physics in a plasma system with relatively few particles. However, as with all other particle methods, PIC methods suffer from numerical noise
such that they have difficulty in simulating some problems, e.g., the problem with
large dynamic ranges in velocity space. To remedy this deficiency, there are usually
two approaches. One method is the so-called δf method [13, 28, 21], discretizing
only the perturbation δf with respect to an equilibrium state f0 based on a particle
method. The δf method has been successfully used in realistic applications, e.g., microtubulence in magnetic confined plasmas [22]. The limitation of this method is that
it can only be applied to the problems which are close to equilibrium. An alternative
approach is through periodically reconstructing the distribution function on a grid in
phase space. Such remapping technique has been used in particle methods in fluid
dynamics [8, 2], i.e., vortex methods and smoothed particle hydrodynamics (SPH),
to maintain regularity of the particle distribution and thereby improve accuracy, but
has much more limited use in PIC methods in plasma physics. It is worth mentioning
that early work of Denavit [12] and more recent work of Vadlamani [31] and Yang
[3] used the idea of remapping for PIC methods. However, They all used low order
interpolation function which results in a first order method overall.
We studied a high order remapping scheme to PIC methods for the solution of the
one-dimensional Vlasov-Poisson equations early [32]. Meanwhile, we provided a local
redistribution technique such that the positivity of the distribution function could be
preserved after high-order remapping. The initial numerical experiments on a set of
classical plasma problems in one dimension are very encouraging. Remapping significantly reduces the numerical noise and results in a more consistent second-order
convergence rate in the electric field error. We also investigated the effects of integrating mesh refinement to the uniform remapping. This is motivated by the observation
that remapping, a numerical diffusive procedure, tends to create a large number of
small-strength particles at the low density region of the distribution function. Mesh
refinement has the potential to reduce this side effect.
In this paper, we extend the algorithm to the solution of the two-dimensional
Vlasov-Poisson equations. This includes the use of an efficient Poisson solver with infinite domain boundary conditions for beam problems. High-dimensional simulations
are very expensive with respect to memory usage. We perform the simulation on a
parallel machine using domain decomposition in physical space. A scalable implementation based on domain decomposition in phase space will be discussed. We consider
two types of numerical tests: plasma problems including linear Landau damping and
the two stream instability, and a beam problem based on the paraxial model [15].
The rest of the paper is organized as follows. In §2, we first review the classical
PIC methods for the Vlasov-Poisson equations. An efficient algorithm which solves
the Poisson equation with infinite boundary conditions is described. Then we present
the high-order and positive remapping on a hierarchy of locally-refined grids in two

3
dimensions. §3 discusses the parallel implementation of the algorithm. We show the
numerical results in §4. Conclusion and future research will be given at the end.
2. Algorithms.
2.1. PIC methods. PIC methods are based on the Lagrangian description of
the Vlasov equation
df (X, V , t)
= 0,
dt

(2.1)

where the characteristics (X(t), V (t)) are the solution of the equation of motion:
dX
dV
s
= V (t),
= (−1) (E(X, t) + E e (X, t))
dt
dt

(2.2)

with initial conditions X(t = 0) = x and V (t = 0) = v.
In the beginning, the distribution function is approximated by a collection of
point particles,
X
f (x, v, t = 0) ≈
qk δ(x − xk )δ(v − v k ),
(2.3)
k

where (xk , v k ) is a initial particle location at the cell center of a grid in phase space
(quite start). qk = f (xk , v k , t = 0)hx hy hvx hvy is the weight of a particle. Then each
particle follows a trajectory described by the equation of motion,
dqk
= 0,
dt

dX̃ k
(t) = Ṽ k (t),
dt

dṼ k
s
(t) = (−) (Ẽ k (t) + E ek (t)),
dt

(2.4)

where X̃ k (t = 0) = xk and Ṽ k (t = 0) = v k .
At any time that a smooth representation of the distribution function is required,
we approximate the function with a collection of finite size particles, where the exact
delta function is replaced by a smoothed delta function. That is,
X
f (x, v, t) ≈
qk δεx (x − X̃ k (t))δεv (v − Ṽ k (t)), t > 0.
(2.5)
k

The smoothed delta function satisfies
Z
δε (y)dy = 1

(2.6)

R2

and
 
1
Y
1
yd
δε (y) =
u
,
εd
εd

(2.7)

d=0

where u is any interpolation function and ε is the stencil size. Usually, the stencil
size for the smoothed delta function in physical space εx is chosen as the same as
the mesh spacing of the Poisson solver. The typical interpolation function for PIC
methods is the first-order interpolation function
(
1 − |z| 0 ≤ |z| ≤ 1
u1 (z) =
(2.8)
0
otherwise.
The flow of a PIC scheme is

4
• Assign particle charges on a grid in physical space,
X
ρ̃(xj , t) =
qk δεx (xj − X̃ k (t)),

(2.9)

k

where j ∈ Z2 are the node index of the Cartesian grid in physical space.
The grid size is chosen as the same as the stencil size of the smoothed delta
function εx . In the case of the first-order interpolation function, for each
node j, the sum is restricted to the particles with |xj − X̃ k (t)| ≤ εx .
• Solve the Poisson equation on the grid with a second-order finite-difference
method:
− (4H φ)j = −

1
X
φj+ed − 2φj + φj−ed
s
= (−1) ρ̃j + ρbackground (2.10)
ε2xd
d=0

and
d

Ẽ j =

φj−ed − φj+ed
.
2εxd

(2.11)

ρbackground is the background charge density if applicable. With given
boundary conditions, the discrete Poisson equation is usually solved by a
fast Poisson solver, such as FFTs or multigrid methods.
• Interpolate the calculated field back to the particle locations with the same interpolation function in equation (2.7). It is worth mentioning that a different
interpolation function will introduce self-force errors [7].
• Integrate the equation of motion numerically, for example, using the secondorder Runge-Kutta method.
2.2. Solving the Poisson equation with infinite domain boundary conditions. To model beam problems, the Poisson equation with infinite domain boundary
conditions needs to be solved. We compute the solution using a new version of the
James-Lackner method [20] by McCorquodale et al. [23, 24]. This method solves two
Dirichlet boundary problems plus a boundary to boundary convolution.
We briefly describe the algorithm below. Assume D0 is the support domain of the
right-hand side ρ, we can solve the Poisson equation with infinite domain boundary
conditions on a slightly larger domain D1 > D with inhomogeneous Dirichlet boundary conditions. The boundary value can be calculated by Green’s function convolution
from the source ρ to the domain boundary ∂D1 . The volume source ρ to boundary
∂D1 convolution is relatively expensive, in particular for 3D problems. Instead of
using a volume to boundary convolution, we can compute the boundary value by
performing a boundary ∂D1 to boundary ∂D2 convolution and solving another Poisson equation on a domain D2 > D1 > D with Dirichlet boundary condition. The
procedure of James’ algorithm is (Figure (2.1)) :
• Step 1: Solve the Poisson equation on domain D1 with homogeneous Dirichlet
boundary conditions
− ∆φ1 = ρ

on

D1 ,

φ1 = 0 on

∂D1 .

(2.12)

• Step 2: Calculate the surface charge on ∂D1
∂ρ1 =

∂φ1
∂n

on

∂D1 .

(2.13)

5
• Step 3: Perform a boundary to boundary convolution from ∂ρ1 to ∂D2
Z
∂φ2 =
G(x − y)∂ρ1 (y)dAy
(2.14)
∂D1

using fast multiple methods [17].
• Step 4: Solve another Poisson equation on domain D2 with inhomogeneous
Dirichlet boundary conditions
− ∆φ2 = ρ

on

D2 ,

φ2 = ∂φ2

∂u2 = G · ∂ρ
u1 = 0

∂ρ =

on

∂D2 .

(2.15)

u2 = ∂u2

∂u1
∂n

u1
∆u1 = ρ

Fig. 2.1. James Algorithm Left: Solve for u1 on D1
perform convolution Right: Solve u2 on D2

u2
∆u2 = ρ

Middle: Calculate surface charge and

2.3. Particle Remapping. The convergence of particle methods for the onedimensional Vlasov-Poisson equations has been investigated by Cottet and Raviart
[9]. Their result shows that particle overlapping and regularization are important for
the convergence of the methods. Specifically, the truncation error of a particle method
is amplified by a time dependent exponential term. Based on Cottet and Raviart’s
work, we extend the error analysis to PIC methods [32]. Our result is one order higher
in the truncation error. However, as in Cottet and Raviart’s analysis, the truncation
error is amplified by a time dependent exponential term. The analysis motivates
the use of remapping technique, a widely used strategy in particle methods in fluid
dynamics, to control the exponential error. The basic idea of remapping is simple.
Since particles will gradually move away from the exact trajectories due to numerical
error, we can reduce the displacement by periodically reproducing the distribution
function f (x, v, t) on a grid by interpolation. A new set of particles, which are created
from the grid representation, then replace the distorted particle distribution. The
later step is identical to the initial step of PIC methods that we initialize the particle
positions and weights in equation (2.3). The error due to remapping will depend on
the order of the interpolation function.
In the previous work [32], we successfully applied the remapped PIC method to
the one-dimensional Vlasov-Poisson system. The remapping scheme was extended in
three aspects compared with the standard scheme. First, we used high-order interpolation functions which improve accuracy but do not preserve positivity. Second,
we preserved the positivity of a high-order interpolation by redistributing the excess
charge into its local neighborhood. The local redistribution algorithm is based on
the mass redistribution idea of Chern and Colella [5], which is first applied to enforce
positivity preservation by Hilditch and Colella [18]. Third, instead of reinitializing

6
on a uniform grid, we reproduced the distribution function on a hierarchy of locallyrefined grids. Remapping on a hierarchy of locally-refined grids significantly reduces
the number of small-strength particles located at the tail of the distribution function.
The high-order, positive, and adaptive remapping scheme in high dimensional phase
space is described below.
2.3.1. High-order Remapping. The overall accuracy introduced by remapping will be one order lower than the order of the interpolation function since we lose
one order of accuracy in the evolution step. For example, the interpolation function
with second-order accuracy only results in a first-order method overall. In this paper,
we consider an interpolation function with third-order accuracy derived by Monaghan
[25]. The function in one-dimensional can be expressed as

3s3
5s2

0 ≤ s = |x|
1 − 2 + 2
h ≤1
2
(2.16)
W4 (x, h) = 12 (2 − s) (1 − s) 1 ≤ s = |x|
h ≤2


0
otherwise.
The one-dimensional expression can be generalized to four dimensions by tensor product,
W 4 (xi − xk ) =

3
Y

W4 (xdi − xdk , hd ),

(2.17)

d=0

where hd is the remapping mesh spacing in phase space. i and k denote the index for
the cell-centered grid and the particles, respectively.
This function, called a modified B-spline, conserves the total charge and represent
a quadratic polynomial exactly. In addition, the first- and the second-order derivative
of W4 (x, h) are continuous. The smoothness property of this modified B-spline is particularly good for scattered data interpolation. However, as with all other high-order
interpolation functions, W4 is not positivity preserving. An interpolation function
without positivity might create nonphysical negative charge. This should be avoided
in simulations.
2.3.2. Positivity. The positivity preserving algorithm is based on the mass redistribution idea of Chern and Colella [5], first applied to enforce positivity preservation by Hilditch and Colella [18]. In the algorithm, we redistribute the undershoot of
cell i
δfi = min(0, fin )

(2.18)

to its neighboring cells i + ` in proportion to their capacity ξ
n
ξi+` = max(0, fi+`
).

(2.19)

The distribution function is conserved, which fixes the constant of proportionality
n+1
n
fi+`
= fi+`
+

ξi+`
neighbors
P

δfi

(2.20)

ξi+k

k6=0

for ` 6= 0 such that cell i + ` is a neighbor of cell i. Superscript n and n + 1 denote
the interpolated value before and after redistribution, respectively.
The drawback of this approach is that positivity is not guaranteed in a single
pass. One might have to apply the method iteratively. In practice, however, we find
a few iterations are sufficient.

7
2.3.3. Mesh Refinement. Mesh refinement is an attractive option in improving the efficiency of phase-space remapping. The distribution function in phase space
is inhomogeneous, for example, a Maxwellian distribution in velocity space. When
we represent the system by particles, it is best that we can have each particle carries a similar amount of weights. Remapping on a hierarchy of locally-refined grids
is a good strategy for creating a such set of particles. From another point of view,
remapping through interpolation is a numerically diffusive procedure. This results
in a large number of small-strength particles near the tail of the distribution function. The situation becomes worse as we apply remapping frequently. Remapping
on a hierarchy of locally-refined grids, with a coarser grid covering the tail of the
distribution function, can reduce the number of those small-strength particles. In the
following, we present the algorithm of remapping with mesh refinement. In designing
the algorithm, we have two guiding principles: the total charge should be conserved
and the overall accuracy on the field needs to be maintained.
Before explaining the algorithm, we introduce the definition of a composite grid.
We define a hierarchy of cell-centered grids Ω` , where 0 ≤ ` ≤ `max . Ω`=0 is the
coarsest grid that covers the whole problem domain. The finer grids Ω`>0 are constructed as a union of cell-centered rectangles (see Figure 2.2). The mesh spacing of
each level is h` = h`−1 /r `−1 , where r `−1 is the refinement ratio of level ` − 1. In four
dimensions, h` ∈ R4 and r ` ∈ Z4 . The composite grid consists of valid grids at all
levels, where a valid grid is defined as a region not overlain by a finer grid. That is,
Ωc =

`X
max


`
Ω` \ P`+1
(Ω`+1 ) ,

(2.21)

`=0
`
where P`+1
is the operator projecting from level ` + 1 to level `.
At the beginning, a set of particles are created from the cell center of the composite
grid. In the remapping step, each particle first finds the valid cell in the composite
grid it belongs to. One particle can only belong to a single valid cell. If the cell is
far enough away from a coarse fine-interface, the charge can be interpolated on the
grid as in equation (2.16). If the cell is near a coarse-fine interface such that the
interpolation stencil intersects the coarse-fine interface, special care must be taken.
First, we interpolate the charge on the surrounding cells as usual. After deposition,
we know that not all deposited cells are valid cell. We need a further step to transfer
the charge from invalid cell to valid cell. There are two cases depending on where the
invalid cell is located. If the invalid cell is in a coarser level and it is covered by the
valid cells of a finer level, we transfer the deposited charge from the coarser level to
the finer level through interpolation. On another hand, if the invalid cell is outside the
grid of the current level, the charge is transfered by projection. Figure 2.2 shows the
algorithm in two dimensions. The four-dimensional case can be generalized easily. It is
worth mentioning that we lose one order of accuracy in interpolating the coarser level
charge into the finer level. However, since the coarse-fine interface is in co-dimension
one, the expected accuracy in the field, e.g., second-order, will be preserved in L∞
norm error. Our current implementation doesn’t have time-dependent adaptivity.
This feature can be incorporated by selecting some refinement criterion, for example
each cell in phase space has similar number of particles.

3. Parallel Implementation and Issues. The parallel implementation of the
algorithm is straightforward based on domain decomposition in physical space. The
physical space is decomposed into M patches. Given a parallel machine with N

8
Ω`c

Ω`c

Ω`+1
c

Ω`+1
c

Fig. 2.2. Cross signs denote the particle locations. The valid and the invalid deposited cells
are denoted by filled circles and open circles, respectively. The refinement ratio is r 0 = (2, 2) in the
plots. Left: Particle is at the coarser level side. Right: Particle is at the finer level side. Cross
signs denote the particle locations.

processors, each patch is assigned to a processor cyclically. Particles are assigned to
patches according to their physical space positions. Using MPI, patches communicate
with each other through ghost-cells and particles move between patches. This is the
default implementation in Chombo software [27].
The current implementation is not a scalable algorithm (weak scaling) because of
decomposition in physical space only. The potential issue is that when the problem size
increases, since the number of processors is scaled in proportion to the problem size in
physical space, the computational time and memory usage will increase in proportion
to the problem size in velocity space. In the worst case, the processor will be out of
memory. An alternative implementation is based on domain decomposition in phase
space. In this implementation, it might be the case that particles belong to the same
cell in physical space (equation (2.9)) are distributed on different processors. Since the
Laplacian operator is linear, we can choose to solve the Poisson equation separately
on different processors. The total fields are then obtained by MPI Allreduce.
4. Numerical Tests. We demonstrate PIC methods with adaptive phase-space
remapping on a set of classical plasma and beam problems in two dimensions, including linear Landau damping, the two stream instability, and beam propagation in the paraxial model. For satisfying the overlapping condition, we choose
εx /hx = εy /hy = 2, where εx = (εx , εy ).
h

We use Richardson extrapolation for error estimate. If Ẽ is the electric field
computed with the initial phase space discretization (hx , hy , hvx , hvy ) and integration
2h

step size ∆t, and Ẽ computed with (2hx , 2hv , 2hvx , 2hvy ) and 2∆t, the relative
solution error in direction d is defined as
ehd = |Ẽdh − Ẽd2h |.

(4.1)

q is the order of the method and is calculated by

q = min log2
d

||e2h
d ||
||ehd ||


.

(4.2)

9
4.1. Linear Landau Damping. The initial distribution for linear Landau
damping is
f0 (x, y, vx , vy ) =

1
exp(−(vx2 + vy2 )/2)(1 + α cos(kx x) cos(ky y)),
2π

(4.3)

where α = 0.05, kx = ky = 0.5 and vmax = 6.0. The physical domain is (x, y) ∈
[0, L = 2π/kx ] × [0, L = 2π/ky ] with periodic boundary conditions. Particle charges
with strength less than 1.0 × 10−9 are ignored. In the simulation, we apply remapping
every 5 PIC time steps.
In the first test, we are interested in the evolution of the amplitude of the electric
field. According to Landau’s theory, the electric field is expected to decrease exponentially with damping rate γ = −0.394. The behavior of exponential decay has been
observed by many other authors, mostly calculated by grid methods [26, 16, 10, 11].
We initialize the problem on two levels of grids with base level at hx = hy =
L/32, hvx = hvy = vmax /16. The velocity space is refined on sub-domain v ∈ [−3, 3] ×
[−3, 3] with a refinement ratio 2. The PIC step size is dt = 1/8. We compare the
simulation with and without remapping in Figure (4.1). In the case with remapping,
the computed damping rate is very close to the theoretical value. The simulation
without remapping fails to track the exponential decay.
1

1

γ=-0.394
hx=L/32, hv=hmax/16

0.1

0.1

0.01

0.01

0.001

0.001

0.0001

0.0001

1e-05

γ=-0.394
hx=L/32, hv=hmax/16

1e-05
0

2

4

6

8

10

t
(a) with remapping

12

14

16

0

2

4

6

8

10

12

14

16

t
(b) without remapping

Fig. 4.1. The amplitude of the electric field for 2D linear Landau damping problem. Scales
(hx , hvx ) above denote the particle grid mesh spacing at the base level, where hx = hy and hvx = hvy .
With remapping, the computed damping rate is very close to the theoretical value γ = −0.394.
Without remapping, the simulation fails to track the exponential decay after a few damping circles.

In the second test, we compare the electric field errors and corresponding convergence rates with and without remapping in Figure (4.2) and (4.3). We see that
remapping significantly reduces the electric field errors and improves their corresponding convergence rates.
4.2. The Two Stream Instability. The initial distribution for the two stream
instability is
f0 (x, y, vx , vy ) =

1
exp(−(vx2 + vy2 )/2)(1 + α cos(kx x))(1 + 5vx2 ),
12π

(4.4)

where α = 0.05, kx = 0.5 and vmax = 9.0. The physical domain is (x, y) ∈ [0, L =
2π/kx ] × [0, L = 2π/ky ] with periodic boundary conditions. Particle charges with
strength less than 1.0 × 10−9 are ignored. As in linear Landau damping problem, we

10

hx=L/32, hv =hmax/16
hx=L/64, hvx=hmax/32
hx=L/128, hvx=hmax/64

error

0.001

x

0.0001
1e-05

7

convergence rate

0.01

1e-06

hx=L/64, hv =hmax/32
hx=L/128, hvx=hmax/64

6

x

5
4
3
2
1
0
-1

0

2

4

6

8

10 12 14 16

0

2

4

6

t

8

10

12

14

16

t

(a) L∞ errors of Ẽx

(b) Convergence rates

Fig. 4.2.
Error and convergence rate plots for 2D linear Landau damping problem with
remapping. Scales (hx , hvx ) above denote the particle grid mesh spacing at the base level, where
hx = hy and hvx = hvy . Second-order convergence rates are obtained. (a) the L∞ norm of the
electric field errors on three different resolutions. (b) the convergence rates for the errors on plot
(a).

7

error

0.001
0.0001
1e-05

hx=L/32, hv =hmax/16
hx=L/64, hvx=hmax/32
hx=L/128, hvx=hmax/64
x

1e-06
0

2

4

6

8

10 12 14 16

t
(a) L∞ errors of Ẽx

convergence rate

0.01

hx=L/64, hv =hmax/32
hx=L/128, hvx=hmax/64

6

x

5
4
3
2
1
0
-1
0

2

4

6

8

10

12

14

16

t
(b) Convergence rates

Fig. 4.3. Error and convergence rate plots for 2D linear Landau damping problem without
remapping. Scales (hx , hv ) above denote the particle grid mesh spacing at the base level, where
hx = hy and hvx = hvy . The errors without remapping are much larger than the case with remapping
(see Figure (4.2)). (a) the L∞ norm of the electric field errors on three different resolutions. (b)
the convergence rate for the errors on plot (a).

apply remapping every 5 PIC time steps. The velocity space is refined on sub-domain
v ∈ [−4.5, 4.5] × [−4.5, 4.5] with a refinement ratio 2.
We compare the electric field errors and their convergence rates with and without
remapping as before. Figure (4.5) shows the L∞ norm of the errors at the case without
remapping in three different resolutions. The corresponding convergence rates are
shown on the right of the error plots. Second-order convergence rates are lost at the
early time of the simulation. Comparing with the results with remapping in Figure
(4.4), we see that remapping extends the second-order convergence rates to longer
times.
We also compare the projected distribution function on plane (x, vx ) at the same
instant time t = 20 by both methods in Figure (4.6). For visualization purpose, in the
case without remapping, we interpolate the particle-based distribution function to a

11
grid in phase space. We see that the classical PIC method results in a noisy solution
(see Figure (4.6(b))). Figure (4.6(a)) shows the distribution function computed by
the PIC method with remapping.
1

convergence rate

0.1

error

7

hx=L/32, hv=hmax/16
hx=L/64, hv=hmax/32
hx=L/128, hv=hmax/64

0.01
0.001
0.0001
1e-05

hx=L/64, hv=hmax/32
hx=L/128, hv=hmax/64

6
5
4
3
2
1
0
-1

0

5

10

15

20

25

0

5

10

t

15

20

25

t

(a) L∞ errors of Ẽx

(b) Convergence rates of the errors on the left

Fig. 4.4.
Error and convergence rate plots for the two stream instability with remapping.
Scales (hx , hv ) above denote the particle grid mesh spacing at the base level, where hx = hy and
hvx = hvy . The PIC step size is dt = 1/8 at the lowest resolution. (a) the L∞ norm of the electric
field errors on three different resolutions. (b) the convergence rate for the errors on plot (a).

7

convergence rate

1

error

0.1
0.01
0.001
hx=L/32, hv=hmax/16
hx=L/64, hv=hmax/32
hx=L/128, hv=hmax/64

0.0001
1e-05
0

5

10

15

t
(a) L∞ errors of Ẽx

hx=L/64, hv=hmax/32
hx=L/128, hv=hmax/64

6
5
4
3
2
1
0
-1

20

25

0

5

10

15

20

25

t
(b) Convergence rates of the errors on the left

Fig. 4.5. Error and convergence rate plots for the two stream instability without remapping.
Scales (hx , hv ) above denote the particle grid mesh spacing at the base level, where hx = hy and
hvx = hvy . The PIC step size is dt = 1/8 at the lowest resolution. (a) the L∞ norm of the electric
field errors on three different resolutions. (b) the convergence rate for the errors on plot (a).

4.3. Semi-Gaussian Beam. The paraxial model is an approximation to the
steady-state Vlasov-Maxwell equation in three dimensions. The K-V distribution is
a measure solution of the paraxial model. Given an arbitrary initial distribution, we
can focus a beam with the same matching forces for the K-V beam using the concept
of equivalent beam. Here, we consider an initial semi-Gaussian beam focused by an
uniform electric field using the concept of equivalent beam. The model has been
considered by many authors [29, 10, 11].
In the test, the beam is composed of ionized potassium. The physical parameters
are the following: current I = 0.2A, beam velocity vb = 0.63 × 106 m/s, and the radius

12

(a) with remapping at time t = 20

(b) without remapping at time t = 20

Fig. 4.6.
Comparison of F (x, vx ), the distribution function projected on space (x, vx ), at
the same instant of time t = 20 with (Left) and without remapping (Right) for the two stream
R R∞
instability problem. The projected value is F (x, vx ) = 0L −∞
f (x, y, vx , vy )dydvy . The grid-based
distribution function is obtained by reproducing the particle-based distribution function through a
second-order interpolation. We initialize the distribution function on two levels of grids, with base
level at hx = hy = L/64, hvx = hvy = vmax /32. The grid is refined by factor of 2 in velocity space
on sub-domain v ∈ [−4.5, 4.5]×[−4.5, 4.5]. The classical PIC method results in a noisy solution with
large errors in maximum. Both numerical noise and errors in maximum are significantly reduced by
using remapping. The negative minimum in the case without remapping is a superficial effect due
to project 4D data to 2D using a high-order interpolation for visualization purpose.

of the beam a = 0.02m. We choose the tune depression η = 1/2. For the normalization
of the paraxial model, we refer to the work of Filbet and Sonnendrucker [15]. We use
normalization parameters (x0 , v0 ) = (a, x2avb ). This results in the normalized Poisson
system
Z
− 4φ =
f dv, −∇φ = E,
(4.5)
R2

with initial semi-Gaussian distribution
(
2
f0 (x, y, vx , vy ) =

4(1−η )
πη 2

exp(−(vx2 + vy2 )/2), x2 + y 2 ≤ 1
0,
otherwise

(4.6)

and the external matching field
E e (x, y, t) = −

4
(xex + yey ).
η2

(4.7)

For the numerical parameters, we choose (hx = |Lx |/128, hy = |Ly |/128) and (hvx =
vmax /128, hvy = vmax /128), where (Lx , Ly ) = (−10, 10) and vmax = 10. The PIC
time step is dt = 0.00052925.
Figure (4.7) shows the projection of the distribution function on planes (x, vx )
with and without remapping, respectively. The simulation with remapping gives a
well-resolved result which preserves the positivity of the distribution function. Meanwhile, we show the root mean square (RMS) quantities of the semi-Gaussian beam

13
in Figure (4.8). Although the RMS quantities of the semi-Gaussian beam are oscillatory, they remain close to the quantities of the associated K-V beam and they
converge as the resolution increases. As mentioned by the other authors [10], the
oscillatory behavior is due to the fact that the semi-Gaussian is not exactly a steady
state distribution.
5. Conclusion. In this paper, we have presented the adaptive remapped PIC
method to the high-dimensional Vlasov equation and demonstrated in linear Landau
damping, the two stream instability, and the beam propagation problems. The new
method reduces the numerical noise significantly. There are two extensions of the
current research. The first will be the development of a scalable algorithm based
on domain decomposition in phase space. The second will be the introduction of
time adaptivity to the current algorithm that the hierarchy of locally-refined grids are
dynamically created from the particle distribution at every remapping step.
Acknowledgments. This work was supported by the U.S. Department of Energy Office of Advanced Scientific Computing Research under contract number DEAC02-05CH11231 at the Lawrence Berkeley National Laboratory. G. H. Miller was
supported by DOE contract number de-sc0001981.
REFERENCES
[1] C. K. Birdsall and A. B. Langdon, Plasma physics via computer simulation, Institute of
Physics Publishing, Bristol and Philadelphia, 1991.
[2] A. K. Chaniotis, D. Poulikakos, and P. Koumoutsakos, Remeshed smoothed particle hydrodynamics for the simution of viscous and heat conducting flows, J. Comput. Phys., 182
(2002), pp. 67–90.
[3] Yang Chen and Scott E. Parker, Coarse-graining phase space in ∆f particle-in-cell simulations, Phys. Plasmas, 14 (2007), p. 082301.
[4] C. Z. Cheng and G. Knorr, The integration of the Vlasov equation in configuration space, J.
Comput. Phys., 22 (1976), pp. 2330–2351.
[5] I.-L. Chern and P. Colella, A conservative front tracking method for hyperbolic conservation
laws, tech. report, Lawrence Livermore National Laboratory, 1987. UCRL-97200.
[6] P. Colella, M. R. Dorr, J. A. F. Hittinger, and D. F. Martin, High-order finite-volumn
methods in mapped coordinates, Submitted to J. Comput. Phys., (2010).
[7] P. Colella and P. C. Norgaard, Controlling self-force errors at refinement boundaries for
AMR-PIC, J. Comput. Phys., 229 (2010), pp. 947–957.
[8] G.-H. Cottet and P. D. Koumoutsakos, Vortex Methods: Theory and Practice, Cambridge
University Press, Cambridge CB2 2RU, UK, 2000.
[9] G.-H. Cottet and P. A. Raviart, Particle methods for one-dimensional Vlasov-Poisson
equations, SIAM J. Numer. Anal., 21 (1984), pp. 52–76.
[10] N. Crouseilles, M. Gutnic, G. Latu, and E. Sonnendrucker, Comparision of two eulerian
solvers for four-dimensional vlasov equation: Part ii, Communications in nonlinear science
and numerical simulation, 13 (2008), pp. 94–99.
[11] N. Crouseilles, G. Latu, and E. Sonnendrucker, A parallel vlasov solver based on local
cubic spline interpolation on patches, J. Comput. Phys., 228 (2009), pp. 1429–1446.
[12] J. Denavit, Numerical simulation of plasma with periodic smoothing in phase space, J. Comput. Phys, 9 (1972), pp. 75–98.
[13] A. M Dimits and W. W. Lee, Partially linearized algorithms in gyrokinetic particle simulation,
J. Comput. Phys., 107 (1993), pp. 309–323.
[14] E. Fijalkow, A numerical solution to the Vlasov equation, Comput. Phys. Comm., 116 (319–
328), p. 1999.
[15] F. Filbet and E. Sonnendrucker, Modeling and numerical simulation of space charged
dominated beams in the paraxial approximation, Math Mod. Meth. Appl. Sci., 16 (2006),
pp. 763–791.
[16] F. Filbet, E. Sonnendrucker, and P. Bertrand, Conservative numerical schemes for the
Vlasov equation, J. Comput. Phys., 172 (2001), pp. 166–187.

14

Fig. 4.7. Time evolution of the projection of the distribution function on (x, vx ) at time
t = 0, 0.3176, 0.5239, 0.7410, 0.9527 (from top to bottom). The projected value is F (x, vx ) =
RLR∞
0 −∞ f (x, y, vx , vy )dydvy . The grid-based distribution function is obtained by reproducing the
particle-based distribution function through a second-order interpolation. The columns on the left
and on the right are simulations running with and without remapping, respectively. The classical
PIC method results in a noisy solution with large errors in maximum.

15
0.51

0.51
K-V beam
SG beam with remapping
SG beam w/o remapping

0.5

K-V beam
SG beam with remapping
SG beam w/o remapping

0.5

0.49

0.49

0.48

0.48

0.47

0.47

0.46

0.46
0

0.2

0.4

0.6

0.8

1

1.2

0

0.2

0.4

t

0.6

0.8

1

1.2

t

(a) RMS of x (64)

(b) RMS of x (128)

1.1

1.1
K-V beam
SG beam with remapping
SG beam w/o remapping

1.08

K-V beam
SG beam with remapping
SG beam w/o remapping

1.08
1.06

1.06

1.04
1.04
1.02
1.02
1
1

0.98

0.98

0.96
0

0.2

0.4

0.6

0.8

t
(c) RMS of vx (64)

1

1.2

0

0.2

0.4

0.6

0.8

1

1.2

t
(d) RMS of vx (128)

Fig. 4.8.
Comparison of RMS quantities of semi-Gaussian beam with remapping, without
remapping, and its equivalent K-V beam. We show the quantities in two different resolutions (64
vs 128 grid size in each dimension). Although the values are oscillatory, they remains close to the
equivalent K-V beam and they converge to the equivalent K-V beam as the resolution increases.

[17] Leslie Greengard, The rapid evaluation of potential fields in particle systems, PhD thesis,
Yale University, 1987.
[18] J. Hilditch and P. Colella, A projection method for low Mach number fast chemistry reacting flow, in Proc. AIAA Aerospace Sciences Meeeting, Reno, NV., 1997.
[19] R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles, McGraw-Hill,
New York, 1981.
[20] R. A. James, The solution of Poisson’s equation for isolated source distributions, J. Comput.
Phys., 25 (1977), pp. 71–93.
[21] W. W. Lee, T. G. Jenkins, and S. Ethier, A generalized weight-based particle-in-cell simulation scheme, Comput. Phys. Comm, 182, pp. 564–569.
[22] Z. Lin, T. S. Hahm, W. W. Lee, M. Tang W, and R. B. White, Turbulent transport reduction
by zonal flows: massively parallel simulations, Science, 281, pp. 1835–1837.
[23] P. McCorquodale, P. Colella, G. T. Balls, and S. B. Baden, A scalable parallel Poisson
solver in three dimensions with infinite-domain boundary conditions, In 7th International
Workshop on High Performance Scientific and Engineering Computing, (2005), pp. 814–
822.
[24]
, A local corrections algorithm for solving Poisson’s equation in three dimensions, Comm.
App. Math. and Comp. Sci., 2 (2007), pp. 57–81.
[25] J. J. Monaghan, Particle methods for hydrodynamics, Comput. Phys. Rep., 3 (1985), pp. 71–
124.
[26] T. Nakamura and T. Yabe, Cubic interpolated propagation scheme for solving the hyperdimensional Vlasov-Poisson equation in phase space, Comput. Phys. Comm., 120 (1999),
pp. 122–154.
[27] Applied Numerical Algorithms Group of Lawrence Berkeley National Laboratory,

16
Chombo software. https://commons.lbl.gov/display/chombo/.
[28] S. Parker, W. Lee, and R. Santoro, Gyrokinetic simulation of itg driven turbulence in 3d
toroidal geometry, Phys. Rev. Lett., 71 (1993), pp. 2042–2045.
[29] E. Sonnendrucker, F. Filbet, A. Friedman, E. Oudet, and J.-L. Vay, Vlasov simulations
of beams with a moving grid, Comput. Phys. Comm., 164 (2004), pp. 390–395.
[30] E. Sonnendrucker, J. Roche, P. Bertrand, and A. Ghizzo, The semi-Lagrangian method
for the numerical resolution of Vlasov equations, J. Comput. Phys., 149 (1998), pp. 201–
220.
[31] Srinath Vadlamani, Scott E. Parker, Yang Cheng, and Charlson Kim, The particlecontinuum method: an algorithmic unification of particle-in-cell and continuum methods,
Comput. Phys. Comm., 164 (2004), pp. 209–213.
[32] B. Wang, G. Miller, and P. Colella, A particle-in-cell method with adaptive phase-space
remapping for kinetic plasmas, SIAM in Sci. Comput., 33 (2011), pp. 3509–3537.

