Parallelization of a local gyrokinetic Vlasov
simulation code for peta-scale computing
S. Maeyama 1, T.-H. Watanabe 2, Y. Idomura 1, M. Nakata 1, M. Nunami 2, and A. Ishizawa 2
1

1

Japan Atomic Energy Agency, 2-166 Omotedate, Obuchi, Rokkasho-mura, Kamikita-gun 039-3212, Japan
2
National Institute for Fusion Science, 322-6 Oroshi-cho, Toki 509-5292, Japan

Introduction

Turbulent transport is one of the most important issues inmagnetically-confined plasma researches. Pressure gradientsof confined plasma destabilize micro-instabilities, and thedestabilized micro-fluctuations drive plasma turbulence via
nonlinear coupling between perturbations of electric potential
and plasma distribution function. The plasma turbulencecauses particle and heat transport perpendicular to magnetic
fields and degrades confinement. Typical wavelengths of theplasma turbulence are of the order of gyroradii of charged
particles, while its time scales are slower than gyrations.
To treat plasma turbulence, the gyrokinetic theory hasbeen developed. The equations describe time evolution ofgyrophase-averaged distribution functions and electricpotential with retaining finite gyroradius effects. A lot ofnumerical simulations based on the gyrokinetic equations has
been carried out and have contributed to understandings ofplasma turbulence[1]. The development of numerical
techniques for gyrokinetic simulations has continued toimprove their applicability, accuracy and efficiency.
Gyrokinetic simulations require expensive computationalresources, since they have to solve time evolution in fivedimensional phase space. We have investigated ion-scale
plasma turbulence by using our gyrokinetic simulation code
GKV[2] on tera-scale supercomputers. Plasma turbulence,
however, is originally multi-scale physics, which includes
temporal and spatial scales of electrons and ions. Problem
sizes for resolving the multi-scale turbulence are 1836
(corresponding to the ion-to-electron mass ratio) times larger
than that of ion-scale turbulence. To deal with thisnumerically challenging problem, we have to develop
numerical methods for applying the GKV code to peta-scalecomputing.

2The GKV code

The GKV code is originally developed to investigate iontemperature-gradient-driven turbulence with the adiabatic
electron approximation. Extensions of the code for treating
both of kinetic ions and electrons are in progress. Since theemployed numerical algorithms are principally the same, we
treat gyrokinetic equations with the adiabatic electron
approximation in the following manuscript.
2.1 Governing equations
The GKV code solves the so-called δf gyrokinetic
equations, where the distribution function is split into theequilibrium part FM and the perturbed part δf. Then the time
evolution of the gyrophase-averaged perturbed ion
distribution function 𝛿𝑓i̅ �𝒓, 𝑣|| , 𝜇; 𝑡� is described by thegyrokinetic Vlasov equation,𝜇∇|| 𝐵 𝜕𝜕𝑩+ �𝑣|| + 𝒗d + 𝒗E � ∙ ∇ −� 𝛿𝑓i̅ = 𝑆 + 𝐶, (1)𝑚i 𝜕𝑣||
𝜕𝑡𝐵
where v||, vd and vE are the velocity parallel to theconfinement magnetic field, the perpendicular magnetic drift
velocity and the perpendicular E×B drift velocity due toelectric potential perturbations. The term with the magnetic
moment µ and the ion mass mi represents parallel accelerationby the mirror force. The linear term associated with theequilibrium distribution S contains contributions of theparallel electric field and equilibrium pressure gradients,
which drive micro-instabilities and plasma turbulence. The
model collision operator C is friction and diffusion operatorsin velocity space (v||, µ). The perturbed electric potential φ isgiven by the gyrokinetic quasi-neutrality equation with theelementary electric charge e and the ion equilibrium
temperature Ti,
𝑒𝐹�i − M (𝜙 − 𝜙�)� 𝑑𝑣 3 = δ𝑛e ,� �𝛿𝑓(2)𝑇iwhere 𝜙� is the gyrophase-averaged potential. It should be
noted that the velocity space integral must be taken holding
particle (not gyrocenter) position fixed. The perturbed
electron density ne is assumed to be
𝛿𝑛e 𝑒(𝜙 − 〈𝜙〉)
=,
(3)
𝑛0𝑇e
where 〈⋯ 〉 denotes the flux surface average.
�

2.2 Simulation domain and boundary conditions

In the δf framework, it is assumed that a steadyequilibrium exists and satisfies the MHD equilibrium
condition. Then, we can employ magnetic coordinates asconfiguration space coordinates r = (x, y, z), where theequilibrium magnetic field is described as
𝐵0 𝜕𝒓.
(4)𝑩 = 𝐵0 ∇𝑥 × ∇𝑦 =�𝑔 𝜕𝑧

Fig. 1. An example of the flux-tube simulation domain
is plotted by black lines. Green lines represent a circulartoroidal flux surface.

The flux-surface label x, the field-line label y and the fieldaligned coordinate z correspond to the toroidal coordinates (r,θ, ζ) as x = r - r0, y = r0(qθ - ζ)/q0, z = θ in a large-aspect-ratio
tokamak with concentric circular magnetic flux surfaces,
where q is the so-called safety factor and the quantities with
subscript 0 denotes the value at the center of the simulation
domain. While plasma turbulence has short perpendicular
wavelength, its structure elongates in the direction parallel tothe magnetic field. Therefore, a long and thin simulation
domain along magnetic field lines is suitable for capturing thenature of plasma turbulence with reducing computationalcosts. An example of this flux-tube simulation domain isshown in Fig. 1, which is written by the projection of a box
with short lengths in x and y and a long length in z. The fluxtube model is widely used to analyse turbulent transport in thelocal approximation limit, where the equilibrium quantities
are given by local values.
We assume statistically periodic boundary conditions in x
and y and apply the Fourier decomposition asFig. 2. Strong scaling of the wave-number-spacedecomposition (where Nx = 512, Ny = 512, Nz = 32, Nv||
= 32, Nµ = 16). (k, z, v||, µ) coordinates are decomposed
into 4 (and 8, 16, 32) × 2 × 2 × 2 subdomains and 8
symmetric multi-processing threads are employed persubdomain, which corresponds to 256 - 2048 cores.

𝜙(𝑥, 𝑦, 𝑧) = � � 𝜙𝑘𝑥 ,𝑘𝑦 (𝑧)exp�𝑖(𝑘𝑥 𝑥 + 𝑘𝑦 𝑦)� . (5)𝑘𝑥 𝑘𝑦Additionally, there is the physical periodicity in the poloidal
angle θ as φ(r, θ, ζ) = φ(r, θ+2π, ζ). This leads the modified
periodic boundary condition along the field-aligned
coordinate z,(6)𝜙𝑘𝑥,𝑘𝑦 (𝑧) = Θ𝜙𝑘𝑥+Δ,𝑘𝑦 (𝑧 + 2𝜋),
where the connection phase Θ and connection wave number ∆depends on the poloidal wave number ky[3].
Eqs. (1)-(3) are numerically solved in (kx, ky, z, v||, µ)space except the E×B advection. Since direct calculations of
nonlinear convolutions in wave number space arecomputationally too expensive, the nonlinear E×B advection
term is evaluated in the real space and transformed back to
the wave number space by means of the fast Fouriertransform (FFT) algorithms and the 3/2 de-aliasing rule.

3Parallelization for peta-scale computing

To attain good performance on a distributed-memory
system, domain decomposition by using message passing
interface (MPI) library is necessary, as well as threadparallelization. The original GKV code decomposes fivedimensional phase space in three directions (z, v||, µ), which is
not enough for peta-scale supercomputers. In order toadvance to peta-scale computing, we additionally decompose
wave-number space k = (kx, ky). By using the onedimensional FFT, the E×B advection term is evaluated in thefollowing way:
1D inverse FFT for ky with the decomposition in kx.
Data distribution as the data in kx becomes local.
1D inverse FFT for kx with decomposition in y.
Calculation of the E×B advection term in real space (x, y).

1D FFT for x with the decomposition in y.
Data distribution as the data in y becomes local.
1D FFT for y with the decomposition in kx.
Efficiency of this parallelization method is examined via astrong scaling test. Computations were carried out on theFX10 super computer in the University of Tokyo. In Fig. 2,
relative speed of the computation, which is proportional to
inverse of the elapsed time, is plotted as a function of thenumber of the wave-number-space decomposition. It
demonstrates that the wave-number-space parallelization has
good strong scaling in this range. We note that perpendiculargrid number required for multi-scale turbulence simulations islarger than the case shown above, e.g., Nx = Ny = 4096. The

result shown in Fig. 2 suggests applicability of the wavenumber-space parallelization on peta-scale supercomputers.
When we employ hundreds of thousands of cores, however,
MPI communication time will become crucial, since data
distributions in wave-number space are collective communications. We have to check computation-to-communication
ratio carefully, and may have to employ masking of MPI
communication time by overlapping between computations
and communications.

4Conclusion

We have presented a local gyrokinetic Vlasov simulation
code to treat turbulent transport in magnetically-confinedplasmas. Calculations of time evolution of distributionfunctionsinfive-dimensional
phase
space
are
computationally expensive, even when a flux-tube simulation
domain is employed to reduce numerical costs. The
simulation code is massively parallelized by domain
decomposition in five-dimensional phase space and threadparallelization. First scaling test of the wave-number-spacedecomposition demonstrates good strong scaling.

Acknowledgements
This work is performed with supports of HPCI Strategic
Programs for Innovative Research and JAEA-NIFS
collaboration program.

References
[1]

[2]

[3]

X. Garbet, Y. Idomura, L. Villard, and T.-H. Watanabe,
Gyrokinetic simulations of turbulent transport, Nucl. Fusion, 50,043002, 2010.
T.-H. Watanabe, and H. Sugama, Velocity–space structures ofdistribution function in toroidal ion temperature gradientturbulence, Nucl. Fusion, 46, 24, 2006.
M. A. Beer, S. C. Cowley, and G. W. Hammett, Field-aligned
coordinates for nonlinear simulations of tokamak turbulence,Phys. Plasmas, 2, 2687, 1995.