Development of numerical techniques toward extreme scale fusion
plasma turbulence simulations
Yasuhiro Idomura1 , Susumu Yamada1 , Toshiyuki Imamura2 , Motoki Nakata1 , and Masahiko Machida1121

Japan Atomic Energy Agency, JapanThe University of Electro-communications, Japan

Introduction

The magnetic confinement fusion is one of the mostpromising approaches to a fusion reactor, and in the next
generation experiment, the International Thermonuclear
Experimental Reactor (ITER), the scientific and technological feasibility of a magnetic fusion reactor will bedemonstrated. Fusion plasmas are extremely complex
physical system consisting of multiple fluids (electronsand multiple species ions) coupled through electromagnetic fields and weak Coulomb collisions. Because of thiscomplexity, computer simulations have been established as
essential tools in wide spectrum of fusion science. Among
several critical issues, the turbulent transport is the mostdemanding issue, because five dimensional (5D) gyrokinetic simulations [1, 2] are needed for studying turbulent fusion plasmas. So far, we have been developing a
global Gyrokinetic Toroidal 5D Eulerian code (GT5D) [3].
GT5D is based on the so-called global full-f approach, inwhich the equilibrium part f0 and turbulent perturbation
δf of particle distribution are solved simultaneously using
the same first principles in a full torus computational domain. In contrast to conventional δf simulations, whichsolve only δf with fixed f0 by assuming complete spatiotemporal scale separation between f0 and δf , full-f simulations can simulate relevant multi-scale phenomena suchas interactions between macro-scale mean flows and microscale turbulence and evolutions of temperature profiles affected by turbulent energy transport in a self-consistentmanner. The full-f approach requires robust and accurate
numerical treatments, which can resolve small amplitude
perturbations (δf /f0 ∼ 1%) for time scales much longer
than a turbulent correlation time. Such numerical requirements were satisfied by developing a Non-Dissipative Conservative Finite Difference (NDCFD) scheme [4]. The accuracy of GT5D was demonstrated through quantitative
verification studies [3]. In addition, validation studies [5, 6]
captured qualitative transport properties in the experiment, such as the stiffness of temperature profiles, intermittency of avalanche-like heat transport, spontaneous
plasma rotation, and the plasma size scaling of the heatdiffusivity. Although the rapid progress of computing
power enhanced capabilities of GT5D, we need to establish the path towards next generation Peta-scale and Exascale computing for simulating the ITER, which is several
times larger than existing fusion devices. In such extreme
scale computing, severe requirements on the memory usage and the parallel efficiency are expected. In this work,we present numerical techniques required for extreme scalefusion plasma turbulence simulations.

2Calculation models

We consider the electrostatic ion turbulence describedby gyrokinetic ions and adiabatic electrons in a tokamak
configuration, which consists of nested magnetic surfaces
with an axisymmetric torus topology. GT5D is based on
the modern gyrokinetic theory [1], which describes dynamics of the plasma particle distribution f in 5D phase space
Z = (t; R, vk , µ). Here, R denotes the configuration space,v is the velocity of particles, vk = b · v and v⊥ = |b × v|
are the velocities in the parallel and perpendicular direc2
tion to the magnetic field, µ = mi v⊥/2B is the magnetic
moment, mi is the mass of ions, and B is the magnetic
field. A conservative form of the gyrokinetic equation isgiven as∂J f∂+ ∇ · (J Ṙf ) +(J v̇k f ) = J C(f ) + J Ssrc ,∂t∂vk(1)

where J is the Jacobian of coordinates, and Ṙ and v̇kare Hamiltonian flows, which are calculated using themagnetic field B and the turbulent electrostatic potential φ, and satisfy an incompressible condition, ∇ · (J Ṙ) +∂vk (J v̇k ) = 0, by their definition. The rhs shows a source
term Ssrc and a Coulomb collision term C(f ), which isgiven by a linear convection-diffusion operator based on
the linear Fokker-Planck collision model [7]. The turbulent
electrostatic potential φ is solved by using the gyrokinetic
Poisson equation,−∇⊥ · P1 ∇⊥ φ + P2 (φ − hφif ) = g − n0 ,(2)where P1 , P2 are coefficients given by plasma parameters,
g is the ion density computed from f , n0 is the equilibrium
electron density, and h·if is an averaging operator over themagnetic surface.
GT5D consists mainly of three solvers, the Gyrokineticsolver (4D convection operator in the lhs of Eq.(1)), theCollision operator (2D convection-diffusion operator in therhs of Eq.(1)), and the Field solver (Eq.(2)). The conservative gyrokinetic equation (1) is discretized using the 4th order accurate NDCFD, while the collision operator uses the6th order accurate centered finite difference. In the timeintegration, the Courant-Friedrichs-Lewy (CFL) condition
is determined by the linear part of the 4D convection operator, which involves the thermal motion of particles in thedirection of magnetic field. Therefore, we adopt the 2ndorder accurate additive semi-implicit Runge-Kutta method(ASIRK) [8], and treat the stiff linear term in an implicit
manner. The gyrokinetic Poisson equation (2) is solved
using a Fourier mode expansion in the azimuthal angle oftorus and a 2D finite element approximation on the crosssection of torus.

3

Numerical techniques

In order to improve the parallel efficiency, the following
key techniques are developed. The first one is a multilayer hybrid parallelization model [9], in which multidimensional domain decomposition is designed following
the physical symmetry properties of each solver, and all
the communications within and among three solvers are
implemented on a single layer of the hierarchical network
consisting of multiple MPI layers and a SMP layer. This
implementation avoids collective communications among
all the MPI processes, and suppresses the size of each
MPI communicator below ∼ 100 up to ∼ 106 MPI processes. In addition to this global network design, we apply
novel communication overlap techniques to each communications. The communication overlap techniques use either non-blocking communications with additional control
sequences or communication threads in a MPI/OpenMP
hybrid parallel model. Our communication overlap techniques are based on conventional MPI libraries, and workon most of existing platforms based on dedicated and commodity networks.
On the other hand, the following low memory usage
techniques are developed. Firstly, in a flat-MPI implementation, the sizes of buffers used in domain decomposition
and in MPI processes show an explosive growth beyond∼ 104 processes. To avoid this issue, a MPI/OpenMP hybrid parallel model is adopted, and the number of MPI
processes is reduced by an order of magnitude. Secondly,
the size of the 2D finite element matrix in the Field solver
exceeds ∼ 100GB for ITER-size parameters. A filtering
technique [10] is developed based on the physical property of fusion plasma turbulence, which is aligned to themagnetic field reflecting a resonance condition. By applying the filtering technique, the size of finite element
matrix is reduced by three orders of magnitude, and theField solver is implemented using a direct solver in LAPACK. Finally, the implicit part of ASIRK is given bya huge asymmetric matrix with 17 stencils (5 points ineach direction of 4D grids). Since the problem size is typically ∼ 1010 , it is difficult to store such a huge matrix
even with sparse matrix storage formats, and therefore,
general iterative matrix solver libraries are not available.
To overcome this difficulty, massively parallel Krylov subspace solvers with optimized domain decomposition aredeveloped by implementing the matrix as finite difference
equations. The original solver was based on a generalized
conjugate residual (GCR) method [11]. In recent work [12],it is extended by adopting preconditioning techniques andadvanced iterative solver algorithms such as GMRES [11]and BiCGstab [11]. By using extended solvers, both thenumber of iterations and the computational cost are reduced, and the processing efficiency on scalar processors
with relatively low memory bandwidth is improved.

Table 1. The memory usage, the sustained performance,
and the ratio of communication cost observed in strong
scaling of GT5D on K with 8 SMP threads. The problem
size is 240 × 64 × 240 × 128 × 32 for 5D grids and 192 ×384 for 2D finite elements. ( ) shows the data without
communication overlap techniques.

Cores
2,048
4,096
8,192
16,38432,768Mem.(TB)
1.8(2.2)
2.0(2.4)
2.5(2.9)
3.2(3.6)
5.0(5.4)

Perf.(TF)
4.3 (3.6)
8.4 (6.3)
16.4(13.5)
29.7(24.8)
54.5(47.6)

Comm.(%)
2.9(14.0)
4.1(13.6)
5.6(18.3)
5.4(17.8)
8.9(23.4)

indicates that the memory usage is significantly reducedby the MPI/OpenMP hybrid parallel model, where thenumber of MPI processes is reduced by 1/8 on K. The
peak ratio of sustained performance with (without) thecommunication overlap techniques is varied from ∼ 13.1%(∼ 10.8%) at 2,048 cores to ∼ 10.4% (∼ 9.1%) at 32,768
cores. Thanks to the multi-layer hybrid parallelization
model, GT5D shows reasonably good strong scaling even
without the communication overlap techniques. However,by applying the communication overlap techniques, the ratio of communication cost is dramatically reduced and issuppressed below ∼ 10% up to 32,768 cores.
5Conclusion

Key requirements in next-generation Peta-scale or Exascale fusion plasma turbulence simulations are lower memory usage and extreme parallelism. To achieve these requirements, novel numerical techniques are developed on
GT5D in various levels ranging from physics models tocomputational models. The performances of these techniques show promising features on K.
Acknowledgement
This work is supported by the MEXT grant for HPCI
Strategic Program Field No.4: Next-Generation IndustrialInnovations. A part of the results is obtained by early
access to the K computer at the RIKEN-AICS. One ofthe authors (Y.I.) is supported also by the MEXT (GrantNo. 22866086) and by the G8 Research Councils Initiatives(the Nuclear Fusion Simulations @ Exa-scale project).
References[1] X. Garbet, et al., Nucl. Fusion 50, 043002(2010).
[2] Y. Idomura, et al., C. R. Physique 7, 650 (2006).
[3] Y. Idomura, et al., Comput. Phys. Commun. 179, 391(2008).
[4] Y. Idomura, et al., J. Comput. Phys. 226, 244 (2007).

4Strong scaling of GT5D on K

Table 1 shows a summary of performances observed instrong scaling of GT5D on K. In this problem size, thesize of finite element matrix becomes ∼ 6.5 GB, while itis reduced to ∼ 74 MB with the filtering technique. At
2,048 cores (256 nodes), the total memory usage is ∼ 1.8TB. The memory usage is significantly increased with thenumber of MPI processes, and at 32,768 cores, it is increased by ∼ 2.77 times. Since the memory usage is determined almost by the number of MPI processes, this data

[5] Y. Idomura, et al., Nucl. Fusion 49, 065029 (2009).
[6] S. Jolliet, et al., Nucl. Fusion 52, 023026 (2012).
[7] X. Q. Xu, et al., Phys. Fluids B 3, 627 (1991).
[8] X. Zhong, J. Comput. Phys. 128, 19 (1996).
[9] Y. Idomura and S. Jolliet, in Proceedings of SC11 (2011).
[10] S. Jolliet, et al., J. Comput. Phys. 231, 745 (2012).
[11] Y. Saad, Iterative Methods for Sparse Linear Systems,
Second Edition, SIAM (2003).
[12] S. Yamada, et al., in this conference.