Plasma and Fusion Research: Regular Articles

Volume 11, 2401096 (2016)

Multi-Scale Simulations of Magnetic Reconnection Using
Particle-in-Cell and Magnetohydrodynamics with Adaptive
Mesh Refinement Technique∗)Tomoya OGAWA, Shunsuke USAMI1) , Ritoku HORIUCHI1) ,Mitsue DEN2) and Kazuyuki YAMASHITA3)
Kitasato University, Sagamihara 252-0373, Japan
National Institute for Fusion Science, Toki 509-5292, Japan
2)
National Institute of Information and Communications Technology, Koganei 184-8795, Japan
3)
University of Yamanashi, Kofu 400-8510, Japan
1)(Received 9 December 2015 / Accepted 27 April 2016)
We performed magnetic reconnection simulations with our simulation code. Our code is a multi-scale threedimensional plasma simulation code in which Particle-in-Cell (PIC) method is used in a restricted region andmagnetohydrodynamics (MHD) simulation is performed in the rest of a simulation box. A hierarchical mesh isused in the MHD region, which is controlled by an adaptive mesh refinement (AMR) technique. Restriction ofthe PIC region and the adaptively controlled mesh realize treatment both of micro-scale particle dynamics andmacro-scale structure with low computational costs. In this paper, we describe methods of our simulation andresults of calculations.
c 2016 The Japan Society of Plasma Science and Nuclear Fusion Research

Keywords: multi-scale simulation, magnetic reconnection, PIC, MHD, AMR
DOI: 10.1585/pfr.11.2401096

1. Introduction
Magnetic reconnection is a multi-scale phenomenon,
in which micro-scale particle dynamics, macro-scale fluid
and magnetic structure are important. A magnetohydrodynamics (MHD) simulation is eﬃcient for treating macroscale fluid, but it needs some kind of phenomenologicalresistivity model or extension of MHD. Diﬀerent models of resistivity may lead to diﬀerent results [1]. Extended MHD simulations are used (e.g. Otto 2001 [2]),but such simulations can only treat micro-scale eﬀects partially. Particle simulations can treat particle dynamics directly, but the method needs high computational resources.
In magnetic reconnection, micro-scale dynamics is important only near reconnection points, while micro-scale dynamics is not important far from reconnection points. Usami et al. developed a multi-scale simulation code [3–5].
They use Particle-in-Cell (PIC) method in a restricted region and perform MHD simulation in other regions. The
multi-scale simulation code achieves the treatment of both
micro-scale particle dynamics and macro-scale structurewith lower computational costs than pure PIC simulations.
Ogawa et al. extended their simulation box by connecting
to a hierarchical mesh [6]. The hierarchical mesh is controlled by an adaptive mesh refinement (AMR) technique.
It can extend a simulation box with a very small increaseauthor’s e-mail: ogawa@kitasato-u.ac.jp∗) This article is based on the presentation at the 25th International Toki
Conference (ITC25).

in computational costs.
In this paper, we applied the simulation code to magnetic reconnection and show results below.

2. Method
A computational box is divided into three domains: aPIC domain, an interface domain, and an MHD domain.
The latter two domains are separated into upper and lower
regions by the centered PIC domain. A regular mesh isused in the PIC and the interface domains, and in the inner
region of the MHD domain. A hierarchical mesh is used in
the outer region of the MHD domain. Figures 1 and 2 show
the structure of domains and meshes of two calculations.
PIC simulation is calculated in the PIC and the interface domains. MHD simulation is calculated in the interface and the MHD domains. In the interface domain, both
PIC and MHD simulations are calculated simultaneously
and exchange their information as below. Two sets of fluid
quantities are derived by PIC and MHD algorithms: QPICand QMHD , respectively. Figure 3 shows time developing
procedure in the interface domain. The upper arrow represents time step in MHD and the lower arrow represents
time step in PIC. At t = ti , QPIC is derived by averaging
PIC fluid quantities over several steps and we obtain mixed
quantities as
Qmix = αQMHD + (1 − α)QPIC ,(1)where the interconnection function α is within 0 ≤ α ≤ 12401096-1

c 2016 The Japan Society of Plasma
Science and Nuclear Fusion Research

Plasma and Fusion Research: Regular Articles

Volume 11, 2401096 (2016)

Fig. 3 Time developing procedure in the interface domain. The
upper and lower rightward arrows are time in MHD andPIC schemes, respectively. Scales on these arrows represent time steps for each scheme.

Fig. 1 Illustration of domains and mesh of calculation 1. Side
lengths of a simulation box are 64 × 128 × 1. Note thatcell sizes drawn in the figure are not real. Actual cellsizes are 0.25 and 0.5 as described in left side of the figure. The maximum resolution corresponds to a regular
256 × 512 × 4 mesh.

Fig. 2 Illustration of domains and mesh of calculation 2. Side
lengths of a simulation box are 128 × 128 × 1. Cell sizes
are the same as in Fig. 1. The maximum resolution corresponds to a regular 512 × 512 × 4 mesh.

and takes a large value near the MHD domain and a smallvalue near the PIC domain. MHD quantities are replacedby the Qmix and the MHD time is advanced from ti to ti +ΔtMHD , where ΔtMHD is a time step for the MHD scheme.
The PIC macroscopic quantities such as the electric andmagnetic fields, the density, and the pressure are replacedby the Qmix . PIC particles are re-distributed at every step

based on the Qmix assuming a Maxwell distribution. At ti <t < ti + ΔtMHD , an interpolation of QMHD (ti ) and QMHD (ti +ΔtMHD ) is used as QMHD in equation (1). See Usami et al.
for details [3–5].
The hierarchical mesh is controlled by an AMR technique. The constitution unit of the mesh is “block” which
consists of 23 cells with physical quantities, pointers to sixneighbors and eight children, and some additional quantities that are needed for control of the mesh structure ormodels. The blocks are connected to each other and construct a mesh structure as a Fully Threaded Tree [7]. The
hierarchical mesh has cell sizes of the box size divided by2l , where l is a natural number. The initial mesh should beset manually. After starting calculation, mesh structure isdecided by criteria for choosing refine or unrefine. We use
expression (2) to detect shock wave and expression (3) todetect contact discontinuity:|Pi − P j |> εsmin(Pi , P j )andui > u j ,(2)|Pi − P j |< εsmin(Pi , P j )and|ρi − ρ j |> εc ,min(ρi , ρ j )(3)

where P, u, and ρ are the pressure, the velocity component
in a direction, and the density, respectively, εs and εc are
the thresholds for shock wave and contact discontinuity,
respectively, and subscriptions i and j denote the cells having smaller and larger coordinate value, respectively [7].
The values of εs and εc should be decided empirically. We
set εs = 0.5 and εc = 0.1 in present calculations. If expression (2) or (3) is satisfied for a pair of cells in a block,then the block is divided into eight children having onehalf the length of their parent. In addition, we also divide
surrounding blocks of the block that has been divided bythe criteria.

3. Models

2401096-2

The MHD basic equations arePlasma and Fusion Research: Regular Articles⎧∂ρ⎪⎪⎪= −∇ · (ρu)⎪⎪⎪∂t⎪⎪⎪⎪⎪⎪
∂(ρu)
B2
⎪⎪
⎪
= −∇ · (ρuu − BB) − ∇ P +⎪⎪
⎪2
⎨ ∂t
⎪
⎪⎪
∂B⎪
⎪
⎪
= ∇ × (u × B)
⎪⎪
⎪∂t
⎪
⎪


⎪
⎪
⎪
∂eB2⎪⎪⎪=−∇·e+P+u − B(u · B) ,⎪⎩ ∂t2

Volume 11, 2401096 (2016)(4)

where ρ, u, B, P, and e are the mass density, the fluidvelocity, the magnetic field, the pressure, and the energydensity, respectively. The governing equations in PIC algorithm are⎧∂B⎪⎪⎪= −∇ × E⎪
⎪⎪∂t⎪
⎪⎪⎪
⎪
∂E⎪⎪
⎪
= ∇ × B − CJ J⎪⎪
⎪
∂t⎪
⎪
⎪⎪⎪
⎪
⎨ ∇·B=0(5)⎪⎪⎪∇ · E = C J ρq⎪
⎪
⎪
⎪
⎪
⎪
d(γk vk )qk⎪⎪⎪(E + vk × B)=⎪⎪
⎪dtm⎪k⎪⎪⎪⎪⎪dxk⎪⎪= vk ,⎩dt
where E, J, and ρq are the electric field, the current density,and the charge density, respectively, γk , vk , qk , mk , and xk
are the Lorenz factor, the velocity, the charge, the mass,and the position of the k-th particle, respectively, J and ρqare given by⎧⎪⎪qk vk S (x − xk )
J=
⎪⎪
⎪
⎪⎨k⎪
⎪⎪
⎪
ρ =qk S (x − xk ),⎪⎪⎩ q(6)kwhere S is a form function of super-particles, and the constant C J is defined as
CJ =



ωpe 2 −1 c −3· n0,ωce
ωce(7)

where ωpe is the electron plasma frequency, ωce is the electron gyrofrequency, n0 is the reference value of particle
number density, and c is the light speed. Normalization
constants of the MHD and PIC algorithms are diﬀerent,
thus we perform unit transformation. See Usami et al.
(2013) for details [4].
Figures 1 and 2 show the structure of domains andmeshes in two calculations: calculation 1 with the box sizeof 64 × 128 × 1 and calculation 2 with that of 128 × 128 × 1(length is normalized by c/ωce ). The two calculations differ only in box sizes. Note that cell sizes drawn in thesefigures are not real. The cell sizes are 0.5 or 0.25, thus
the maximum resolution of calculations 1 and 2 corresponds to regular 256 × 512 × 4 and 512 × 512 × 4 meshes,
respectively. In both calculations, the PIC domain is located on 46 < y < 82, the interface domain is located onFig. 4 Plots of density (top), y-component of velocity (middle),and x-component of magnetic field (bottom) for initialcondition on a line with x = constant and z = constant.
Color represents cell size in initial mesh structure; blackand red indicate cell sizes of 0.5 and 0.25, respectively.

44 < y < 46 and 82 < y < 84, and the MHD domain
is located on outer regions. Cell sizes of the PIC and theinterface domains and the inner part of the MHD domain
are always 0.25. Cell sizes of the outer region of the MHD
domain are 0.25 or 0.5. In initial mesh, cell size is 0.25
in 36.75 < y < 91.25 and is 0.5 in other regions. In outerregion, there are no shock waves and no contact discontinuities in the present model, thus the mesh shows fewchange in its structure.
Initial condition is given at t = 0 as
⎧
2 y − yc⎪⎪a+(1−a)1−tanhρ=ρ⎪0⎪⎪W⎪⎪⎪⎪⎪
⎪
u=0⎪
⎨

(8)
⎪
2 y − yc⎪⎪a+(1−a)1−tanh
P=P⎪
0
⎪
⎪W
⎪⎪⎪
⎪
⎪
⎪ B = B tanh y − yc ,⎩0Wwhere the density and the pressure on y = yc are ρ0 = 1
and P0 = 2/3, respectively, the constant for magnetic
field is B0 = (1, 0, 0) (magnetic field is normalized bycωce me /qe , where me is the electron mass and qe is theelectron charge). The number ratio of background to foreground (Harris-equilibrated) particles is a = 0.25, and thecenter position and the thickness (the length in y-direction)of a current sheet are yc = 64 and W = 4.77, respectively.
Figure 4 shows plots of the initial condition together with
the initial cell sizes.
We give electric field on y-boundary cells to drive
plasma inflow as⎧⎪πt E0⎪⎪⎪1−cos(t < D)
⎨ 2D(9)Ez = ⎪⎪⎪⎪⎩(t ≥ D),
E02401096-3

Plasma and Fusion Research: Regular Articles

Volume 11, 2401096 (2016)

where the terminal strength of electric field is E0 = −0.6
and the increasing duration is D = 12.5 (time is normalized
by c/(ωce vA ), where vA is the Alfvén speed). x- and zboundaries are periodic.

4. Results
Figures 5 and 6 show the results of calculation 1. Figure 5 displays color contours of density, x-component ofvelocity, x- and y-components of magnetic field from left
to right at t = 80 on the upper panels and at t = 119 onthe lower panels. Figure 6 is plots of density, y-component

Fig. 5 Results of calculation 1. Color contours of density, xcomponent of velocity, x-component of magnetic field,
and y-component of magnetic field (from left to right) on
the plane of z = 0.625 at t = 80 (upper) and t = 119(lower).

Fig. 6 Results of calculation 1. Plots of density (top), ycomponent of velocity (middle), and x-component ofmagnetic field (bottom) on the line of x = 32, z = 0.625
at t = 119. Color has the same meaning in Fig. 4.

of velocity, and x-component of magnetic field from top
to bottom at t = 119. Figures 7 and 8 show the results ofcalculation 2. Figure 7 displays color contours of density
on the left panels and x-component of velocity on the rightpanels at t = 96 on the upper panels and at t = 119 on
the lower panels. Figure 8 is the same plots as Fig. 6 att = 119, however, this figure displays the results from calculation 2. In both calculations, inflows run smoothly fromupper and lower boundary to the centered PIC domain viathe MHD domain with the hierarchical mesh and the interface domain, and then approach the center. Reconnection
with single x-point occurs in calculation 1. On the otherhand, two x-points appear in the case of calculation 2. Thisis because the current sheet in calculation 2 is longer thanthat in calculation 1. This result is consistent with Ji &Fig. 7 Results of calculation 2. Color contours of density (left)
and x-component of velocity (right) on the plane of z =0.625 at t = 96 (upper) and t = 119 (lower).

Fig. 8 The same plots as Fig. 6 but for calculation 2 on the lineof x = 64, z = 0.625.

2401096-4

Plasma and Fusion Research: Regular Articles

Volume 11, 2401096 (2016)

Daughton (2011) qualitatively [8]. Fluctuations in the PIC
domain are well suppressed in the interface domain.
We estimate computational gain by restriction of PIC
region. The calculation 1 spent memory size of 105.2 GB
and CPU time of 23.3 hours. For comparison, let us assume that the whole region is calculated by a pure PIC
simulation. Such a calculation will need large memory anda long period of time. Furthermore, memory size and CPU
time in PIC simulations are not proportional to the y-lengthof computational box because the number density of particles depends on the upstream direction. In the upstream
far away from the current sheet, the plasma density is approximately uniform. Thus, we perform two smaller calculations, in which box sizes are 64 × 40 × 1 and 64 × 48 × 1(grid numbers are 256 × 160 × 4 and 256 × 192 × 4), andthen estimate costs in the case of a 64 × 128 × 1 box(256 × 512 × 4 grids) by assuming that the cost increment
is proportional to the elongation of y-length. Spent memory and CPU time in the former calculation are 101.6 GBand 36565 sec. On the other hand, spent memory and CPU
time in the latter calculation are 121.2 GB and 49522 sec.
The box elongation of 8 in y-direction brings enhancements of memory size of 19.6 GB and CPU time of12957 sec. Therefore, we can estimate computational costsin the case of 64 × 128 × 1 as memory size of 317.2 GB andCPU time of 49.7 hours. This indicates that our method reduced memory usage to about one-third and CPU time toabout one-half in comparison with a pure PIC calculation.

5. Conclusion

simulation code to magnetic reconnection. PIC region isrestricted near the current sheet. Inflow from y-boundary
has propagated smoothly to the PIC domain through the hierarchical mesh region and the interface domain. Magnetic
reconnection with one or two x-points occurred on a current sheet. PIC noises have been well suppressed in the interface domain. By restricting the PIC region, we reduced
memory usage to about one-third and CPU time to aboutone-half comparison with a pure PIC calculation. The
outer hierarchical mesh will be able to extend a computational box with low cost enhancement. Thus, we can expect that future simulations with larger boxes will achieve
higher computational gain. By this feature, we plan to connect our code with a global simulation involving a wholesystem, e.g., an Earth’s magnetosphere simulation.

Acknowledgements
This work was supported by JSPS Grant-in-Aid for
Scientific Research (B) 23340182 and the General Coordinated Research at the National Institute for Fusion Science(NIFS15KNTS039).
[1] M. Ugai, Phys. Fluids B 4, 2953 (1992).
[2] A. Otto, J. Geophys. Res. 106, 3751 (2001).
[3] S. Usami et al., Commun. In Comput. Phys. 11, 1006 (2012).
[4] S. Usami et al., Phys. Plasmas 20, 061208 (2013).
[5] S. Usami et al., J. Phys. Conf. Ser. 561, 012021 (2014).
[6] T. Ogawa et al., JPS Conf. Proc. 1, 016013 (2014).
[7] A.M. Khokhlov, J. Comput. Phys. 143, 519 (1998).
[8] H. Ji and W. Daughton, Phys. Plasmas 18, 111207 (2011).

We applied our multi-scale three-dimensional plasma

2401096-5