日本原子力研究開発機構機関リポジトリ
Japan Atomic Energy Agency Institutional Repository

Title
Author(s)Citation
Text Version

Benchmarking of flux-driven full-F gyrokinetic simulations
Asahi Yuichi, Grandgirard V., Idomura Yasuhiro, Garbet X.,
Latu G., Sarazin Y., Dif-Pradalier G., Donnel P., Ehrlacher C.
Physics of Plasmas, 24(10), 102515 (2017)
Publisher

URL

https://jopss.jaea.go.jp/search/servlet/search?5060888

DOI

https://doi.org/10.1063/1.4998015ⓒ 2017 American Institute of Physics. This article may be
downloaded for personal use only. Any other use requires prior

Rightpermission of the author and AIP Publishing. The following article
appeared in Physics of Plasmas, 24(10), 102515 (2017) and may befound at https://doi.org/10.1063/1.4998015.

Japan Atomic Energy Agency

Benchmarking of flux-driven full-F gyrokinetic simulations
Yuuichi Asahi, Virginie Grandgirard, Yasuhiro Idomura, Xavier Garbet, Guillaume Latu, Yanick Sarazin, Guilhem
Dif-Pradalier, Peter Donnel, and Charles Ehrlacher

Citation: Physics of Plasmas 24, 102515 (2017);
View online: https://doi.org/10.1063/1.4998015
View Table of Contents: http://aip.scitation.org/toc/php/24/10
Published by the American Institute of Physics

Articles you may be interested in
Toroidal angular momentum balance during rotation changes induced by electron heating modulation in
tokamak plasmas
Physics of Plasmas 24, 080701 (2017); 10.1063/1.4996017
Verification of linear resistive tearing instability with gyrokinetic particle code VirtEx
Physics of Plasmas 24, 102125 (2017); 10.1063/1.4999166
A conservative scheme of drift kinetic electrons for gyrokinetic simulation of kinetic-MHD processes in toroidal
plasmasPhysics of Plasmas 24, 102516 (2017); 10.1063/1.4995455
Two-fluid tearing mode instability in cylindrical geometry
Physics of Plasmas 24, 072102 (2017); 10.1063/1.4986116
Zonal flows driven by the turbulent energy flux and the turbulent toroidal Reynolds stress in a magnetic fusion
torus
Physics of Plasmas 24, 102508 (2017); 10.1063/1.5004555
Nonlinear gyrokinetic simulation of ion temperature gradient turbulence based on a numerical Lie-transform
perturbation method
Physics of Plasmas 24, 082515 (2017); 10.1063/1.4986395

PHYSICS OF PLASMAS 24, 102515 (2017)

Benchmarking of flux-driven full-F gyrokinetic simulations
Yuuichi Asahi,1,2,a) Virginie Grandgirard,1 Yasuhiro Idomura,2 Xavier Garbet,1Guillaume Latu,1 Yanick Sarazin,1 Guilhem Dif-Pradalier,1 Peter Donnel,1and Charles Ehrlacher1
1
2

CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France
Japan Atomic Energy Agency, 178-4-4 Wakashiba, Kashiwa, Chiba, Japan(Received 28 July 2017; accepted 2 October 2017; published online 20 October 2017)
Two full-F global gyrokinetic codes are benchmarked to compute flux-driven ion temperature gradient (ITG) turbulence in tokamak plasmas. For this purpose, the Semi-Lagrangian code
GYrokinetic SEmi-LAgrangian and the Eulerian code GT5D are employed, which solve the full-F
gyrokinetic equation with a realistic fixed flux condition. The equilibrium poloidal flow profile formation processes are benchmarked and compared against the local neoclassical theory. The simulations above are carried out without turbulence, which agree well with each other and with thetheoretical estimates. Here, a lot of attention has been paid to the boundary conditions, which havehuge impacts on the global shape of radial electric field. The behaviors of micro-instabilities arebenchmarked for linear and nonlinear cases without a heat source, where we found good agreements in the linear growth rates and nonlinear critical gradient level. In the nonlinear case, initialconditions are chosen to be identical since they dominate the transient turbulence behavior. Using
the appropriate settings for the boundary and initial conditions obtained in the benchmarks above, aflux-driven ITG turbulence simulation is carried out. The avalanche-like transport is assessed witha focus on spatio-temporal properties. A statistical analysis is performed to discuss this selforganized criticality (SOC) like behaviors, where we found a 1/f spectra and a transition to 1/f3
spectra at high-frequency side in both codes. Based on these benchmarks, it is verified that the
SOC-like behavior is robust and not dependent on numerics. Published by AIP Publishing.
https://doi.org/10.1063/1.4998015
I. INTRODUCTIONFirst-principles kinetic simulations are considered as
essential tools to understand and predict plasma transportphenomena in magnetic confinement fusion devices. So far,a number of simulation studies have addressed the physicalmechanisms involved in turbulent transport, and comparedthem with experimental results (see the review paper,1 forexample). In most cases, these codes tackle local transportphenomena under fixed plasma profiles. Since they computethe perturbed part of the distribution functions only, they areoften called local df gyrokinetic simulations. Tremendous
efforts have been done for code verification, namely, benchmarking of different local df gyrokinetic codes, in order toensure better credibility.2,3Concerning the experimental situation, the local
approach may fail whenever the characteristic size of turbulence is non-negligible compared with the machine size or
the profile formation is affected by plasma turbulence itself,
possibly leading to an organization of meso-scale structures.
Robust features of non-local transport have been reported.4,5
To address this issue, a new generation of gyrokinetic codes
called full-F gyrokinetic codes have been developed, where
the fluctuations and equilibrium profiles evolve self-consistently.6–10 More importantly, full-F gyrokinetic codes cancompute turbulent heat transport in the presence of a physically relevant fixed flux condition.
a)Electronic mail: yuuichi.asahi@cea.fr

1070-664X/2017/24(10)/102515/17/$30.00

In addition, full-F gyrokinetics is an approach to treat
turbulent and collisional transport consistently. This favorable feature is rather problematic in benchmarking due to thecomplex physics involved in these simulations. In fact, there
exist few benchmark works of full-F gyrokinetic codes11 andthey are addressing fixed-gradient global simulations withadaptive heat source terms12 rather than a realistic heatsource term with a constant heating source. In order toaddress the complicated dynamics of full-F flux-driven simulations, we first decompose the problem into subsets, eachtackling a specific physics issue. For collisional (neoclassical) properties, we carried out simulations where the Ion
Temperature Gradient (ITG) mode is stable so that transport
is only due to collisions. For turbulence, we investigated linear ITG mode and nonlinear decaying ITG turbulence in theabsence of collisions (or weakly collisional). After completing these tests, we move to the benchmark of flux driven
ITG turbulence, which involves all the physical processes
listed above. The physics to be addressed mainly is theavalanche-like transport and profile corrugations and theirstatistical properties. We focus on the self-organized criticality (SOC) like behaviors, which are at the heart of the full-F
flux driven gyrokinetic simulations. It should be stressed thatthese physical processes can only be simulated in the framework of global full-F gyrokinetic models.
This paper contributes to the benchmark of long-timeflux-driven ITG turbulence simulations, which is organized
as follows: Section II describes the full-F gyrokinetic modelwith source and sink terms. In Sec. III, the differences in the24, 102515-1

Published by AIP Publishing.

102515-2

Asahi et al.

Phys. Plasmas 24, 102515 (2017)

numerical schemes are shown. The benchmarks for the neoclassical physics and zonal flow dynamics against theories
are shown in Sec. IV. The linear and nonlinear benchmarks
for gyrokinetics without source and sink terms are shown inSec. V. The benchmark for the flux driven case is summarized in Sec. VI. The results are summarized in Sec. VII.
II. MODELIn this section, we briefly describe the physical models
used in GYrokinetic SEmi-LAgrangian (GYSELA)9,13 andGT5D10 codes. Although these codes can now treat the
kinetic electrons,14 we employ the conventional versions
with adiabatic electrons in this work. Since the gyrokinetic
part is fundamentally the same, we emphasize the differences in the collision operators and sink/source models.
A. Full-F gyrokinetic models

In the modern gyrokinetic theory,15 the time evolutionof the guiding center distribution function f can be computedby the 5D gyrokinetic equation_@J f_ þ @J v k f ¼ J ðDr ð f Þ þ Kð f Þ þ Cð f Þ þ SÞ;þ r  J Rf
@t
@vk
(1)
where J ¼ m2i Bk is the Jacobian in the gyrocenter coordinate and vk is the parallel velocity. mi is the mass of ion,Bk ¼ b  B is the parallel component of B , with the unit
vector in the parallel direction b, and B ¼ B þ ðBvk =Xi Þr  b. Xi ¼ ðei BÞ=ðmi cÞ is the ion cyclotron frequency
with the charge of ion ei and the speed of light c. For theright hand side, Dr and K are a diffusion term and a Krook
operator applied on radial buffer regions, C represents a collision operator, and S is a source term. The details for the collision operator are defined in Sec. II B, and source/sink terms(diffusion or Krook) are described in Secs. II C and II D.
The nonlinear characteristics in Eq. (1) are given byB @Hc
R_  fR; H g ¼
þb  rHmi Bk @vk ei Bkc 
ei rh/ia þ mi v2k b  rb þ lrB ;¼ vk b þ

ei Bk


B
v_k  vk ; H ¼ 
 rH
mi Bk¼B  ei rh/ia þ lrB ;mi B k(2)(3)

where fF; Gg is the gyrokinetic Poisson bracket operator

Xi @F @G @F @GBþ
fF; Gg 
B @a @l @l @a
mi Bk
!
@G @F
c
 rF

rG 
b  rF  rG
@vk @vkei Bk(4)
in the gyro-center coordinates Z ¼ ðR; vk ; l; aÞ, where l is
the magnetic moment and a is the gyro-phase angle. H represents the Hamiltonian of the system defined as1
H ¼ mi v2k þ lB þ ei h/ia :
2

(5)

Þ
Here, h/ia ¼ /da=2p stands for the gyro-averaged electrostatic potential /.
The electrostatic potential / is determined by the quasineutrality condition

q2ti
1 
r? / þ 2 /  h/if
2
kDi
kDe
ð
¼ 4pei f dð½R þ q  xÞd6 Z  nG;eq ;

r? 

(6)

where R þ q is a particle position, d6 Z ¼ m2i Bk dRdvk dlda is
the phase space volume of the gyro-center coordinates, qti is
the Larmor radius evaluated with the thermal velocity vti , kDi
and kDe are the ion and electron Debye lengths, hif is a flux
surface average, and nG;eq is the equilibrium guiding-center
density. The first term of the left hand side stands for the ion
polarization density in the long wavelength approximation.
This assumption should be valid for the core ITG turbulence,
which is characterized by k? qti  1 with the perpendicular
wavenumber k? . By coupling the gyrokinetic equation (1)and quasi-neutrality condition (6), the gyrokinetic model canbe solved self-consistently.
B. Collision operators

In this subsection, we briefly describe the collision operators used in GYSELA13 and GT5D.16 Although there aresome differences, it has been shown that the conservative
properties13,16 and some neoclassical properties17,18 are satisfied in both codes. It should be noted that the finite-Larmorradius (FLR) effect is not taken into account in the collision
operators in GYSELA and GT5D, which is equivalent to acollision operator computed in the drift-kinetic limit.
1. Collision operator in GYSELAThe version of the collision operator considered in thepresent study for GYSELA (other versions exist, such as arecently benchmarked multi-species version19,20) is given bya simplified Lenard-Bernstein collision operator,21 where
only the vk contribution is taken into account
"
!#
1 @
@f

B Dk
 V kf:(7)
Cð f Þ ¼ 
Bk @vk k@vkThe operators Dk and V k , respectively, model a diffusion
and a drag in the parallel velocity direction. The conservation properties of parallel momentum and energy are ensured
by constraining Dk to depend on l only and defining the
local fluid velocity Vk;coll and ion temperature Tcoll in an
appropriate manner.13
2. Collision operator in GT5D

GT5D is equipped with the linearized Fokker-Planck
collision operator22 as Cðf Þ ¼ CT ðf Þ þ CF . The test particle
operator CT is given by102515-3

Asahi et al.

 @  1 @2 @ 
?1 v2 f þk1 uf þ
?2 v4 f
2
@s
@u
2 @s


1 @2 
@2 2
k? v3 f ;þk2 v f þ
2
@s@u
2 @u

Phys. Plasmas 24, 102515 (2017)

CT ð f Þ ¼(8)

where s ¼ 2lB=mi and u ¼ vk  Uk are defined in a moving
frame with respect to the parallel flow velocity Uk andv2 ¼ u2 þ s. The field particle operator CF ensures the conservation properties of parallel momentum and energy. The
definitions of the collision frequencies, ?1 ; k1 ; ?2 ; k2 ;k? , and the field particle operator CF are given in the
Appendix of Ref. 16.
C. Source models

Let us describe the gyrokinetic equation neglecting allterms but the source term@f
¼ Ssrc :(9)@t
We assume that the source term can be expressed in the following way:
Ssrc ¼ S0 SE ðq; EÞSr ðqÞ;(10)

where S0 is the amplitude of the source in physical unit,
E ¼ v2k =2 þ lB is the dimensionless energy, and SE and Sr
represent source profiles in velocity space and real space,respectively. q ¼ r=a is the normalized minor radius used asa radial variable throughout the paper, with a minor radius a.
In GYSELA and GT5D, the radial and velocity space profiles (all variables are dimensionless) are given by
1q  qs  3Ls
Sr;GYS ðqÞ ¼  tanhLs2
q  3Ls  q;þtanh sLs1q  qs
Sr;GT5D ðqÞ ¼ 1  tanh
;
Ls
2

1
E 3
 exp ½E=Ts ; (11)SE;GYS ðq; EÞ ¼ pﬃﬃﬃﬃﬃﬃ 5=2
T2s
3p 2pTs
B0(12)
SE;GT5D ðq; EÞ ¼  ½fM1  fM2 ;
Bkwhere fM1 represents the Maxwellian1
exp E=2T1
fM1 ¼
3=2
ð2pT1 Þ(13)

for the volume averaged temperature T1 and fM2 is for thevolume averaged temperature T2 . If we take T2 ¼ T1  DT,
Eq. (12) can be approximated by
SE;GT5D ðq; EÞ ’

B0 @fM1
B0@lnfM1
DT ¼  fM1
DT : (14)Bk @T
Bk@T

If we neglect the difference in B0 and Bk , of the order of theinverse aspect ratio  ¼ r=R0 and use the relationshiplnfM1 ¼ 3=2lnT1  E=T1 , it leads tofMs E 3
SE;GT5D ðq; EÞ ’ C(15)

Ts Ts 2by replacing T1 by Ts. Since the difference in Bk and B is of
the order of q? and negligible, the major difference between B0
and Bk is of the order of the inverse aspect ratio , that is thedifference between B and B0. In our flux-driven test case witha relatively large aspect ratio (see Table III), the difference in
B0 and Bk is about 10% in the source region, since the source
term is strongly localized in the q < 0:3 region, wherein  issmaller than 0.1. Although this gives the minor difference inthe region with q < 0:3, it may not affect much in the region
with 0:3 < q < 1:0 since the integrated input power which
affects the dynamics in this region is the same for both codes.
On one hand, the GYSELA source model has isotropic velocity space structure, but it gives a momentum input due to thevelocity space Jacobian found in Eq. (1). On the other hand,
the GT5D source model does not impose any extra momentum
input, at the cost of anisotropic velocity space structure. The
impact of momentum input is of the order of q? , but it could
be important when discussing the parallel flow dynamics, as
found in Sec. VI. Equation (15) is the same as Eq. (11) except
for the factor C. The difference is then absorbed in the source
amplitude S0. Figure 1 shows the radial and velocity space profiles of source terms used in this study. For radial profile, we
set qs ¼ 0:01212 and Ls ¼ 0.1 in GYSELA and qs ¼ 0:31
and Ls ¼ 0.1 in GT5D. For the velocity space, we employ
Ts ¼ 1:5T0 and DT ¼ 0:1T0 in order to make the velocity
space structures the same in both codes in case we dismiss thedifference in Bk and B0. Here, T0 is the reference temperature
at the mid-minor radius (see Sec. III D for detail). These
parameters are used in the flux-driven simulations in Sec. VI.
As is clear from Eq. (11), this source term continuously
reduces the number of slow particles and increases the numberof fast particles. This process leads to plasma heating. This isdifferent from a Krook-type source term as used in Ref. 12,
which forces the temperature profile to be fixed on average.
D. Sink models

In practice, GYSELA often employs diffusion terms inthe buffer regions as a sink term (although a Scrape-OffLayer (SOL)-like outer boundary condition has been recently
developed23). However, for better comparison with GT5D,
we used a Krook type operator Kr ðf Þ ¼ ðqÞðf  fM Þ in
GYSELA as well as GT5D. The radial profile of the coefficient ðqÞ is given by"#
1q  qk ðqÞ ¼ 0 1 þ tanh
;
(16)
Lk
2
where  0 is the diffusion coefficient, qk and Lk are, respectively, the location and the steepness of the buffer regions.
For the flux-driven simulations in Sec. VI, we employ
qk ¼ 1:0, Lk ¼ 0.1, and 0 a=vti ¼ 0:1. The Krook operator
forces the distribution function to relax towards the initial
Maxwellian fM in the buffer region. Thus, turbulence isdamped in the buffer region.
III. NUMERICAL SCHEMESIn this subsection, we briefly describe the numerical
schemes used in the codes. We will focus, in particular, on102515-4

Asahi et al.

Phys. Plasmas 24, 102515 (2017)

FIG. 1. Profiles of source terms in
GYSELA and GT5D. (a) Radial structure and (b) Velocity space structure.

some differences that involve different numerical properties
which could end up with different results. For example, thenumerical schemes to solve the gyrokinetic Poisson equation, which are relevant to the boundary conditions, may
have a strong impact on the global mode structure of theelectrostatic potential. The comparisons of numerical costs
are provided in the Appendix.
A. GYSELA code

As named, GYSELA (GYrokinetic SEmi-LAgrangian)code uses a Backward Semi-Lagrangian scheme to solve
the advection term. In order to avoid a costly multidimensional interpolation, the advection operator is split into
three different operators based on Strangs’s operator
decomposition method:24 2D ðr; hÞ þ 1D ðvk Þ þ 1D ðuÞ
advection terms. Here, GYSELA employs the toroidal coordinate system ðr; h; uÞ for a spatial coordinate and ðvk ; lÞfor a velocity space coordinate. The Poisson equation issolved with Fourier expansion in toroidal and poloidal
directions and finite differences in the radial direction. The
difficulty lies in the estimate of flux surface averaged electrostatic potential h/if , since this term is nonlinear in h dueto the h dependence of the Jacobian. To tackle this, thePoisson equation is decomposed in two equations, one forthe difference between the potential and its flux surface
average and the other for the flux surface averaged component. These equations can be solved sequentially, andfinally, the electrostatic potential can be obtained. The
boundary conditions used in these equations are describedin Sec. III C.
Contrary to local flux-tube codes, the gyro-averaging
operation is not straightforward. In the flux-tube codes, thegyro-average operator can be represented by a Bessel function in the Fourier space based on the periodicity assumptions both in radial and poloidal directions (that is the local
approach). In a global code, a periodic boundary condition in
the radial direction is not allowed so that the Fourier representation of the gyro-average operator is no longer feasible.
Therefore, we have to compute the gyro-average operator
numerically in real space with some approximations. In
GYSELA, there are two possibilities to approximate thegyro-average: A Pade approximation and integration on
gyro-circles using Hermite interpolation. Here, we employ

the Hermite interpolation approach with 8 sampling points.
The time integration is performed using a predictor-corrector
method. The detailed implementation of GYSELA is foundin Ref. 13.
B. GT5D code

In GT5D, the gyrokinetic Poisson bracket operator is discretized using the fourth-order non-dissipative and conservative finite difference scheme,10 which is a variant of theMorinishi scheme25 extended to incompressible Hamiltonian
flows, and the gyrokinetic Poisson equation Eq. (6) is computed using a toroidal mode expansion and a 2D finite element
approximation on the poloidal plane. Similar to GYSELA, thegyro-average is evaluated by a finite-point sampling technique, with 40 sampling points. GT5D solves the gyrokinetic
Vlasov equation in the cylindrical coordinate ðR; Z; uÞ andthe gyrokinetic Poisson equation in the flux coordinates
ðw; h ; uÞ, with the poloidal flux w and the straight-field-line
poloidal angle h . The velocity space coordinate is uniform inðvk ; wÞ directions, with w ¼ ð2B0 l=mi Þ1=2 ¼ ðB0 =BÞ1=2 v? .
The use of an aligned coordinate in the w direction is beneficial to capture trapped-passing boundaries.
The time integration is performed using the secondorder additive semi-implicit Runge-Kutta method and a stiff
linear term involving the parallel streaming is treated implicitly. An implicit part is solved using the generalized conjugate residual method. The detailed implementation of GT5D
is found in Refs. 10 and 16.
C. Boundary conditions

Due to the differences in the numerical schemes used tosolve the Poisson equation, GYSELA and GT5D impose different boundary conditions on electric potential /. The nonaxisymmetric components of / are forced to be zero at bothinner (if exists) and outer boundaries in both codes. The axisymmetric /00 mode is treated in a different way. In GT5D,
the Dirichlet boundary condition /00 ¼ 0 is applied at theedge, while the natural boundary condition is imposed at themagnetic axis. It should be noted that the use of cylindricalcoordinates ðR; Z; uÞ in GT5D allows one to avoid the singularity at the magnetic axis. In GYSELA, the plasma region
qmin < qð¼ r=aÞ < qmax , with some tiny value qmin is solvedin order to avoid the singularity in 1=r at the magnetic

102515-5

Asahi et al.

Phys. Plasmas 24, 102515 (2017)

FIG. 2. Mesh grids in GYSELA (a)and GT5D (b). Although emphasized
by setting qmin ¼ 0:1 in figure (a), theactual hole inside the inner boundary israther small with qmin ¼ 0:01 evenwith the Neumann condition.

axis.13 Conventionally, for the axisymmetric /00 mode, theDirichlet boundary condition is imposed at the edge
(q ¼ qmax ), while the Dirichlet or Neumann conditions canbe chosen at the inner boundary (q ¼ qmin ). The Neumannboundary condition can be applied only with qmin 102 ,which ensures that the poloidal dependence of the Jacobian
is negligible. Recently, we employed the scheme in Ref. 26
to avoid imposing any boundary condition at the innerboundary (q ¼ qmin ). In this scheme, the radial mesh widthis adjusted to eliminate the coefficients of the grids at the
inner boundary so that there is no inner boundary condition
needed. Here, the first radial grid point r0 is set as Dr=2
with the radial grid width Dr. The interpolation in the
Vlasov solver inside r0 ¼ Dr=2 is performed with the mapping from the polar to rectangular coordinates, where the
1=r singularity can be avoided in the rectangular coordinate
system.27 With this technique, the theta dependence of theJacobian is taken into account appropriately in GYSELA as
well as in GT5D, except for the linear GYSELA simulations with Neumann boundary condition as described
below.
Three different choices for the inner boundary conditions can be used in GYSELA. However, we employ the
new boundary condition (called the Null boundary condition
hereinafter) for all the simulation shown in this paper except
for linear simulations, which is closer to the natural boundary condition at the magnetic axis in GT5D. In the linear
simulations, we used the Neumann boundary condition withqmin ¼ 0:01 instead. The use of different boundary conditions can be justified, since the axisymmetric /00 ¼ 0 is filtered out in the linear simulations and is nothing to do withthe results. We solve the plasma region qmin < q < qmax ,with qmin ¼ 0 ðor 0:01 for linear simulationsÞ and qmax ¼ 1:0
for GYSELA and qmin ¼ 0 and qmax ¼ 1:0 for GT5D. Figure
2 shows examples of meshes in GYSELA and GT5D. For
simplicity, the low resolution cases with ðNr ; Nh Þ ¼ ð16; 32Þ
for GYSELA and ðNR ; NZ Þ ¼ ð16; 16Þ for GT5D are shown
here. As is easily seen from Fig. 2(a), the core region is wellresolved while the edge region may be under resolved in
GYSELA. Figure 2(a) corresponds to the Dirichlet or
Neumann boundary condition, and is no longer correct for
the Null boundary condition wherein the inner holedisappears.

As summarized in Sec. III E, we choose the resolutions
in each code to be the same in the edge region q 0:9.
D. Normalization

In order to compare the values directly, we employ acommon normalization that is different from the ones usedin both codes. The normalizations used in this paper are summarized in Table I, where the notation X^ is used for thedimensionless value of X.
The subscripts s represents the species. Here, we
employed the reference variables including a reference mass
m0 ¼ A0 mp , a reference charge e0 ¼ Z0 e, a reference magnetic field B0 on magnetic axis and a reference temperature
T0, where mp is the proton mass. The over-lined value means
the value evaluated with the volume averaged density andtemperature, and the value without it means the value at the
mid-minor radius qð¼ r=aÞ ¼ 0:5.
E. Simulation settings

Since there are many different benchmarking parameters, it appears useful to summarize all relevant physical,
geometrical, and numerical settings for the different
benchmark cases (neoclassical, zonal flow, linear microinstabilities, nonlinear decaying turbulence and nonlinear
flux-driven turbulence).
Tables II and III show the relevant numerical and physical settings used in this paper. The numerical settings are
kept the same as in conventional runs in the spirit of presenting the consistency with previous publications. On the other
hand, we use the exactly same physical parameters. In the
TABLE I. Comparisons of the normalizations in this work and the originalones in each code.
Common normalization

GYSELAGT5D

ms ¼ m0 A^ses ¼ Z0 eZ^sl ¼ qt0 l^
t ¼ ^t R0 =vt0
ns ¼ n0 n^s
Ts ¼ T0 T^s
B ¼ B0 B^

ms ¼ m0 A^ses ¼ Z0 eZ^sl ¼ qt0 l^
t ¼ ^t =Xc0
ns ¼ n0 n^s
Ts ¼ T0 T^s
B ¼ B0 B^

ms ¼ m0 A^ses ¼ Z0 eZ^sl ¼ qt0 l^
t ¼ ^t R0 =vt0
ns ¼ n0 n^s
Ts ¼ T0 T^s
B ¼ B0 B^

102515-6

Asahi et al.

Phys. Plasmas 24, 102515 (2017)

TABLE II. Numerical settings for the different benchmarking cases.
GYSELAðNr ; Nh ; Nu ; Nvk ; Nl Þ

Wedge size

Inner boundary

Time step

Collision

Source/Sink

Linear
Neoclassical
Zonal flow
Nonlinear
Flux driven

ð256; 256  512; 128  256; 128; 16Þ
(256, 256, 16, 128, 16)(256, 256, 16, 128, 32)(256, 256, 32, 128, 16)(256, 256, 32, 128, 16)

1
1/2
1
1/6
1/6

Neumann
Null
Null
Null
Null

10-40
10
5
15
15

OFF
ON
OFF
ON
ON

OFF
OFF
OFF
OFF
ON

ðNR ; Nz ; Nu ; Nvk ; Nv? Þ

Wedge size

Boundary

Time step

Collision

Source/Sink(200-400, 200-400, 16, 128, 16)(200, 200, 1, 128, 16)(200, 200, 1, 128, 32)(200, 200, 32, 128, 16)(200, 200, 32, 128, 16)

1/5–1/30
1
1
1/6
1/6

NaturalNaturalNaturalNaturalNatural

5-10
5
5
5
5

OFF
ON
OFF
ON
ON

OFF
OFF
OFF
OFF
ON

GT5D
Linear
Neoclassical
Zonal flow
Nonlinear
Flux driven

linear simulations, we increased the perpendicular grid resolutions and decreased the time step size for n ¼ 25 and 30
cases in both codes. Correspondingly, we increased the resolution in toroidal direction in GYSELA while keeping thewedge size. In contrast, we reduced the wedge size while
keeping the resolution in toroidal direction in GT5D.
IV. BENCHMARKS OF COLLISIONAL TRANSPORTAND ZONAL FLOW DYNAMICS

Let us begin with the benchmarking of collisional transport and poloidal flow dynamics, namely, neoclassical transport in the axisymmetric limit and zonal flow dynamics. The
neoclassical results can be obtained at the drift kinetic limits
as shown in Ref. 17.
A. Neoclassical benchmarks

As for the neoclassical benchmarks, we compare theneoclassical ion heat flux and the parallel flow relation at aquasi-equilibrium state. The grid points used in GYSELA
and GT5D are ðNr ; Nh ; Nu ; Nvk ; Nl Þ ¼ ð256; 256; 16; 128; 16Þ
and ðNR ; Nz ; Nu ; Nvk ; Nv? Þ ¼ ð200; 200; 1; 128; 16Þ, respectively. The q1  a=qti ¼ 150 is employed here, which is
carefully chosen to be equivalent for both codes (see the normalization defined in Sec. III D). The profiles used in thisneoclassical computation are shown in Fig. 3. As explainedin Sec. III D, the normalized variables like vti ; qti , and Ti are
the ones defined at the mid-minor radius at t ¼ 0 in thefollowing.
The temperature and density gradients are set identical,
that is, gi ¼ Lti =Ln ¼ 1, wherein the ITG mode is stable. In
addition, we consider the plasma with large aspect-ratio withTABLE III. Physical parameters at the mid-minor radius q ¼ 0:5 for the different benchmarking cases.

Linear
Neoclassical
Zonal flow
Nonlinear
Flux driven

R0 =Lti

R0 =Lns ¼ Ti =Te





q

s^

6.92
6.0
107
6.92
10.0

2.22
6.0
107
2.22
2.22

1.0
1.0
1.0
1.0
1.0

0
0.01-10
0.18
0.02
0.02

0.18
0.1
1.4
0.18
0.18

1.4
1.4
0
1.4
1.4

0.78
0.78
0
0.78
0.78

1
a ¼ R0 =a ¼ 1:738=0:3477 ¼ 5, in order to mimic the large
aspect ratio assumption used in the neoclassical theory.
Here, a represents the inverse aspect ratio. In this configuration, we can focus on the neoclassical dynamics (without turbulence) and direct comparisons against neoclassical theory
are possible. The low resolution in the toroidal direction
reflects the axisymmetrical nature of this test.
For the collisionality dependency of the ion thermal diffusivity, we use a well-known formula derived by ChangHinton28
pﬃﬃﬃ
vi =vi; GB ¼ 2 2 K2 ;! 
0:66 þ 1:881=2  1:54
B20
K2 ¼
1=2
B2 f
1 þ 1:03 þ 0:31
! " 
 2 1 #
0:79
0:74 3=2
B20
B
;
þ 1=2

B2 f
1 þ 0:74 3=2
B20 f

(17)
with the gyro-Bohmpﬃﬃcoefficientvi; GB ¼ q2ti vti =R0 andﬃ
 ¼ r=R0 . The factor 2p2ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
results from our definition of thethermalvelocityv¼
T=m, whereas Chang-Hinton usesti
i
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
vTi ¼ 2Ti =m. The effective collision frequency is given by ¼ qR0 =ð3=2 sii vti Þ. Figure 4 shows the collisionality
dependence and its comparison with Chang-Hinton’s theory.
The neoclassical ion heat diffusivity vNeo
¼ QNeoii =ðni rTi Þ
is defined using the heat flux driven by the curvature drift
ð

Neo
3
d vf vD  rr :
(18)
Qi 
f

Concerning experimentally relevant plasmas where thecollisionality  varies a lot in the radial direction, it is quite
important to verify the applicability of the codes for a wide
range of  . We found good agreement for a wide range of , namely, banana and plateau regimes as shown in Fig. 4.
From the first order gyrokinetic theory, the radial force
balance equation can be derived
dlnhni if
 dlnhTi if
ei Er
mi Xi dw
þ 1  kð ; Þ;
hUk if þ¼
Ti
Ti RBt dr
dr
dr
(19)

102515-7

Asahi et al.

Phys. Plasmas 24, 102515 (2017)

FIG. 3. (a) Initial radial profiles ofdensity, temperature, safety factor, and(b) temperature and density gradientsused for neoclassical benchmarks.

where Er is the radial electric field, Uk the parallel mean
flow velocity, R ¼ R0 þ r cos h is the plasma major radius,
Bt the toroidal magnetic field, and w the poloidal flux. The
coefficient k is a coefficient of the neoclassical poloidal flow.
Based on Hinton’s theory,29 the coefficient (the H-H formula) is
!
1=2


1:17  0:35
1 þ 2 3 : (20)
 2:12 3kð ; Þ ¼
1=2
1 þ 0:7
Similar to Ref. 16, we compare the radial profile of theelectric field against the theoretical estimate. Since neoclassical theory gives an estimate of the poloidal flow, we canestimate the radial electric field based on the radial force balance equation Eq. (19). Figure 5 shows the radial electricfield computed from the electrostatic potential and from Eq.
(19) for  ¼ 0:1 case. We drew two theoretical estimatesusing profiles from GYSELA and GT5D, respectively. The
time evolution of profiles is very close with each other but
not exactly the same, especially when looking at the derivatives like temperature gradients. We found good agreementswith the H-H formula in both codes, confirming that the

FIG. 4. The collisionality dependence of the ion heat diffusivity vi with theradial average for r=a ¼ 0:4  0:6. The dashed curve shows the analyticalestimation by Chang-Hinton’s formula.

collision operators in both codes produce the appropriate
neoclassical poloidal flows.
B. Zonal flow dynamics

In the collisionless regime, the dynamics of the zonal
component /00 can be calculated analytically.30,31 Namely,
the time oscillatory components or geodesic acoustic modes(GAMs) are decaying towards the residual poloidal flow. In
this section, we computed the time evolution of zonal flow
component /00 and compared its damping rate and residuallevel against theory. The initial distribution function is given
by f ¼ feq ð1 þ e sin ð2pðr=2aÞÞÞ with a perturbation amplitude e ¼ 105 . As shown in the work by Sugama and
Watanabe, the time evolution of the zonal-flow potential isgiven by
/00 ðtÞ ¼ /00 ð1Þ þ /00 ð0Þ  /00 ð1Þ cos ðxG tÞexp ðctÞ;(21)
where /00 ð1Þ ¼ /00 ð0Þ=ð1 þ 1:6q2 =1=2 Þ is the residual
zonal flow derived by Rosenbluth-Hinton.30 xG and c are thefrequency and damping rate of GAM. Since the analyticaltheories30,31 are derived in the local limit, we employ radially uniform profiles to eliminate global effects. The safety

FIG. 5. Radial profile of electric field Er against the theoretical estimate
at tR0 =vti ¼ 130 from GYSELA and GT5D. Solid lines correspond to
EHHr
GYSELA results and dashed lines correspond to GT5D results. Theoreticalare estimated using the profiles of GYSELA and GT5D,
values EHHr
respectively.

102515-8

Asahi et al.

factor profile is given by qðqÞ ¼ 1:5 ð8q 2 ½0; 1Þ, and the
aspect ratio is set as a ¼ a=R0 ¼ 0:18. The density and temperature profiles are almost constant with R0 =LTi ¼ R0 =Ln ¼ 107 . The grid points used in GYSELA and GT5D are
ðNr ; Nh ; Nu ; Nvk ; Nl Þ ¼ ð256; 256; 16; 128; 32Þ and ðNR ; Nz ;Nu ; Nvk ; Nv? Þ ¼ ð200; 200; 1; 128; 32Þ, respectively. Theq1
 ¼ 150 is used here.
Figure 6 shows the time evolution of /00 mode at
q ¼ 0:5 and corresponding analytical results. They show
good agreements with the theoretical estimates of frequencies, damping rates, and residual flow levels. Table IV showsthe estimates of frequencies and damping rates from simulations and theories. From simulations, the GAM frequency isestimated by FFT with the time window of 0 tR0 =vti 50,
and the damping rate is estimated for the time window of
10 tR0 =vti 20. The shorter time window for the damping
rate estimate is chosen to avoid the time-evolution of theradial wavenumber kr of /00 mode.
For the wavenumber used to estimate the damping rate
analytically, we employed the same initial radialpﬃﬃﬃ wavenumpﬃﬃﬃ
ber of /00 mode, that is, kr qTi ¼ ð2p=2aÞ  2qti ¼ 2pq
¼ 0:0296.
Difficulty of this benchmark lies in the time-evolution of kr,
which is largely affected by boundary conditions. As describedin Eq. (2.10) in Ref. 31, the damping rate of GAM depends on
kr. Thus, it is important to give the appropriate boundary conditions for two codes. In addition, it would be helpful to investigate the point near the mid-minor radius q ¼ 0:5, in order tominimize the influence from the boundaries.
V. BENCHMARKS OF MICROINSTABILITY

To benchmark microinstability, we considered the
Cyclone-Base case (CBC),2 which is a standard benchmark
case for local df gyrokinetic simulations. Contrary to the
local df gyrokinetic simulations, we have to take into accountPhys. Plasmas 24, 102515 (2017)
TABLE IV. Estimates of frequencies and damping rates from simulations
and the theories.30,31

GYSELA
GT5D
Sugama-Watanabe

cR0 =vti

xG R0 =vti

0.031
0.036
0.039

2.644
2.643
2.690

the plasma profiles. The initial plasma distribution function
is given by a local Maxwellian. The density ni, temperature
Ti, safety factor q, and the magnetic field B are given by thefollowing form:Dnq  0:5;(22)
ni ðqÞ ¼ n0 exp  tanh
Ln
Dn =a

DTi
q  0:5
Ti ðqÞ ¼ Ti0 exp 
;
(23)
tanh
LTi
DTi =aqðqÞ ¼ 0:854 þ 2:184q2 ;(24)

B0 R0 reh þ eu ;
R qR0

(25)

B¼

where q ¼ r=a is the normalized minor radius, B0 is themagnetic field at the magnetic axis, and eh and eu are poloidal and toroidal unit vectors. We choose Dn =a ¼ 0:3 andDTi =a ¼ 0:3 for density and temperature profiles. Figure 7
shows the initial radial profile of ion temperature Ti, density
ni, and safety factor q.
The electron and ion temperature profiles are the sameso that s ¼ Ti =Te ¼ 1 everywhere. At the mid-minor radius
q ¼ 0:5, the profile gives the CBC parameters as R0 =LTi
¼ 6:92; R0 =Ln ¼ 2:22; s ¼ Ti =Te ¼ 1:0, q ¼ 1.4, and s^¼ 0:78. Although the normalized Larmor radius for original
DIII-D shot 8149932 is given by q1
 ¼ 180 with R ¼ 1:71 m
and a ¼ 0:625 m, we consider q1
 ¼ 150 with R ¼ 1:01 m
and a ¼ 0:365 m to reduce the computational costs.
A. Linear benchmarks

FIG. 6. The time evolution of the normalized E  B flows in RosenbluthHinton test. The residual value is compared with Rosenbluth-Hinton theory.
The GAM frequencies and damping rates are compared with theoretical predictions by Sugama and Watanabe.

As for linear benchmarks, we compared the linear growthrate and real frequency of ITG modes. As discussed in Ref. 13,
the separation between linear and non-linear terms is not possible in a full-F code. In order to avoid the coupling of different
toroidal modes, we filtered out all the toroidal modes except theinitial one. The filtering is performed just after solving thequasi-neutrality equation. For example, Fig. 8 contains 6 different linear simulations initialized with n ¼ 5; 10; 15; 20; 25; 30.
The corresponding poloidal mode number is computed by using
the relationship kh qti ¼ nq0 q a=r0 with r0 ¼ a=2 andq0 ¼ qðr0 Þ. Since the high toroidal mode number is equivalent
to a tiny structure both in toroidal and poloidal directions, wedoubled the numbers of grid points both in toroidal and poloidal
directions for the simulations with n 25.
Figure 8 shows the linear growth rate and real frequencyof ITG mode computed at the mid-minor radius q ¼ 0:5
from 6 different runs. Both the linear growth rate and real
frequency show quantitative agreements between GYSELA
and GT5D.

102515-9

Asahi et al.

Phys. Plasmas 24, 102515 (2017)

FIG. 7. (a) Initial radial profiles ofdensity, temperature, safety factor, and(b) temperature and density gradients
based on Cyclone-base case.

the vk direction. In the collisionless limit, the fine structuresin the vk direction generated by parallel phase mixing willreach the grid scale at a recurrence time. Therefore, the solution will be unphysical after the recurrence time in the collisionless limit.
The initial condition is given by
0Fs ¼ Fs;eq @1 þ e

nmax mX
max ðnÞ
X

1
cos ðmh þ nu þ dmn ÞA;

n¼1 m¼mmin ðnÞ(26)

FIG. 8. Linear growth rate and real frequency obtained by GYSELA (blue)
and GT5D (red).

The global structures of the linear eigenmodes have also
been compared. Figure 9 shows the global mode structure of
the linear eigenmode with n ¼ 15. A similar shape of global
eigenmode is found, including the radial peak position andthe ballooning angle.
B. Nonlinear decaying turbulence

In this section, we present the benchmarks for nonlinear
decaying turbulence, i.e., without sink and source terms.
Weak collision with  0:01 is added in order to avoid
unphysical solution coming from the absence of diffusion inwhere mmin ðnÞ ¼ max½1; nq  dm and mmax ðnÞ ¼ min½mmax ;nq þ dm, with nmax ¼ 18; mmax ¼ 70, and dm ¼ 8. The perturbed amplitude is set as e ¼ 106 . Since the decaying turbulence simulations are sensitive to initial conditions, we use
exactly the same initial conditions for the two codes. The
grid points used in GYSELA and GT5D are ðNr ; Nh ; Nu ;Nvk ; Nl Þ ¼ ð256; 256; 32; 128; 16Þ and ðNR ; Nz ; Nu ; Nvk ; Nv? Þ¼ ð200; 200; 32; 128; 16Þ, respectively. In order to save thecomputational cost, we employed one sixth of a torus
based on convergence tests.33 The time step size is set as
DtXi ¼ 15 in GYSELA and 5 in GT5D.
The time evolution of the turbulent heat diffusivity vi
and the normalized temperature gradient R0 =Lti is plotted inFig. 10. The turbulent ion heat diffusivity vi ¼ Qi =ðni rTi Þ
is defined using the heat flux driven by E  B drift velocity

QTurb
i

ð

3



d vf vEB  rr

:

(27)

f

FIG. 9. Linear eigen functions of themost unstable mode with n ¼ 15
obtained by (a) GYSELA and (b)GT5D.

102515-10

Asahi et al.

Phys. Plasmas 24, 102515 (2017)

FIG. 10. Time evolution of the turbulent heat diffusivity vi (a) and normalized ion temperature gradient R0 =Lti(b) averaged for 0:45 q 0:55. The
heat flux level converges to the same
nonlinear critical gradient R0 =Lti
5:5, where the ion heat transport isalmost quenched.

As shown in Fig. 10(a), the time scale for the initial transientand final steady phase shows good agreements in GYSELA
and GT5D. Although the initial saturation level is not exactly
the same, the asymptotic state transport levels are close tozero in both codes, implying that turbulence is almost
quenched. Because of the chaotic nature of plasma turbulence, a perfect agreement in time evolution cannot be
expected after the initial transient phase, that is, the turbulent
phase. Thus, we compare well-defined properties like thecritical temperature gradient and the relaxation process
toward the asymptotic state. The steady state temperature
gradients are R0 =Lti 5:5 in Fig. 10(b), which is consistent
with the effective critical gradients found in Ref. 2. The nonlinear threshold value is confirmed from Fig. 11, whichshows the temperature gradient and turbulent transport values at successive time points during profile relaxation. Bothcodes show a relaxation toward the temperature gradient
R0 =Lti 5:5, where the turbulent transport is almost
quenched. Figure 12 shows the time evolution of the root
mean square of the electrostatic potential ei /=Ti , where

FIG. 11. Time evolution of the turbulent heat diffusivity vi and normalized
ion temperature gradient R0 =Lti averaged for 0:45 q 0:55. The heat flux
level converges to the same nonlinear critical gradient R0 =Lti 5:5, where
the ion heat transport is almost quenched.

zonal components, GAM components and other contributions are separately plotted. Here, we found that the time
evolutions of zonal components are roughly the same anddominate turbulence in the asymptotic state. We can also
confirm that the behavior of the turbulent components arequite similar, both in the linear phase with 50 tvti =R0 70
and the asymptotic state phase with tvti =R0 160.
Figure 13 shows contours of the n 6¼ 0 components ofthe electrostatic potential ei /=Ti in the decaying phase. The
amplitude of the normalized potential is about ei /=Ti 0:03
for both codes. Since the formation of the mean electric field
has a large impact on transport properties, we analyzed thespatio-temporal evolutions of the radial electric field Er andturbulent ion heat flux vturb
as shown in Fig. 14. As expected,i
the heat transport is quenched almost everywhere in bothcodes after tvti =R0 160. The quench of turbulent heattransport is presumably the consequence of the turbulence
regularization by the radial electric field shear found in Figs.
14(c) and 14(f). The temperature gradient is kept almost constant with R0 =Lti 5:5 after the initial transient phase [seeFigs. 14(b) and 14(e)].

FIG. 12. Time evolutions of the root mean square of the electrostatic potential ei /=Ti at the mid-minor radius q ¼ 0:5. The zonal component /0;0 isrepresented by the line with squares, the GAM component P/1;0 is represented by the dotted line and the summation of other modes m;n /m;n areplotted with lines.

102515-11

Asahi et al.

Phys. Plasmas 24, 102515 (2017)

FIG. 13. Contours for n 6¼ 0 components of the electrostatic potential at tvti =R0 ¼ 100 from (a) GYSELA and (b) GT5D.

FIG. 14. Spatio-temporal evolutions of the ion turbulent heat diffusivity vi =ðq2ti vti =R0 Þ [(a) and (d)], the normalized temperature gradient R0 =Lti [(b) and (e)]and the radial electric field shear dEr =dr [(c) and (f)] in GYSELA (upper row) and GT5D (lower row).

In order to compare the global mode structure quantitatively, we computed the radial and time correlation functions~ ¼ /  / . Weof electrostatic potential fluctuations /n¼0
employed the following definition for the correlation function:ðð~ ðt; r; u; h ¼ 0Þi
~ t þ tlag ; r þ rlag ; u; h ¼ 0/
dt drh/uðð;C/ ðtlag ; rlag Þ ¼~ 2 ðt; r; u; h ¼ 0Þi
dt drh/u(28)

102515-12

Asahi et al.

Phys. Plasmas 24, 102515 (2017)

FIG. 15. (a) Temporal and (b) radial
self-correlation functions computedfrom n 6¼ 0 components of electrostatic
potential ei /=Ti for GYSELA and
GT5D, during the time period of
tvti =R0 ¼ 160  180.

where hiu stands for the average in toroidal angle u. As
shown in Fig. 15, both codes show the quite similar correlation length around Drlag =qti 6 and correlation time around
Dtlag vti =R0 2.
VI. BENCHMARKS FOR FLUX DRIVEN SIMULATIONSIn the flux driven simulations, we employ basically thesame physical parameters as for decaying turbulence simulation except for a larger temperature gradient, R0 =Lti ¼ 10 at
the mid minor radius q ¼ 0:5. The input power is set to
2MW, whose radial and velocity space profiles are given inSec. II C. The Krook-type sink operator is applied in thebuffer region near the edge as described in Sec. II D. The
grid points used in GYSELA and GT5D are ðNr ; Nh ; Nu ;Nvk ; Nl Þ ¼ ð256; 256; 32; 128; 16Þ and ðNR ; Nz ; Nu ; Nvk ; Nv? Þ¼ ð200; 200; 32; 128; 16Þ, respectively. As well as the nonlinear decaying turbulence case, we employed the one sixth
wedge torus model. The time step size is set as DtXi ¼ 10 in
GYSELA and 5 in GT5D. The gradient of initial ion temperature profile R0 =Lti ¼ 10 is above the linear and nonlinear
thresholds with R0 =Lti 4:5 and R0 =Lti 6:0, respectively.
The nonlinear threshold value can be identified from Fig. 11,
where both GYSELA and GT5D show similar values,
namely, R0 =Lti 5:5. In this section, we plot the turbulent
ion heat flux Qi rather than the heat diffusivity vi in order todifferentiate the time evolution of turbulent transport andtemperature gradient. As found in Eq. (27), a heat diffusivity
vi is affected by both heat flux Qi and temperature gradient
R0 =Lti . In flux driven simulations, these quantities evolve
self-consistently and exhibit avalanche-like properties. In
addition, the time-averaged temperature profile spontaneously corrugates leading to zonal flows and the interplay
between zonal flows and avalanches gives rise to the socalled E  B staircase.4,5,34,35
Figure 16 shows the spatio-temporal evolutions of theion heat flux Qi, the normalized ion temperature gradient
R0 =Lti and the radial electric field shear dEr =dr. In a source
free region (q ¼ 0:5–0:9), the turbulent heat transport shows
the avalanche-like feature in both codes [see Figs. 16(a) and16(d)]. Correspondingly, the ion temperature gradient profile
in Figs. 16(b) and 16(e) exhibits the near threshold value
around R0 =Lti 6. As found in the time evolution of heatdiffusivity [Figs. 16(a) and 16(d)] and radial electric field
shear [Figs. 16(c) and 16(f)], there are inward and outward
avalanches depending on the sign of the electric field shear.
In the positive shear region with q > 0:6 avalanches propagate outward, while they propagate inward in the negative
shear region with q < 0:6. Both codes agree well on the
directions of avalanche propagation and their relationshipwith the radial electric field shear. Several discussions on
this issue are found in Refs. 16, 36, and 37.
Figure 17 shows (a) the steady state profile of the turbulent ion heat flux Qi =½ni Ti vti q2ti =a2 , (b) the normalized ion
temperature gradient R0 =Lti , (c) the radial electric fieldei R0 Er =Ti , and (d) the parallel flow velocity Uk =vti . The time
average is taken over tvti =R0 ¼ 300–750. First of all, we findquite similar radial profiles of the turbulent heat flux.
Although there are some differences in the inner region with
0 q 0:5, the temperature gradient in this region is close
to the marginal state [Fig. 17(b)] where the transport levelshould be very sensitive to the temperature gradient. In themarginal state, it has been often reported5,34 that the temperature gradient profile, poloidal flow profile, and heat transport interplay via meso-scale structure organization. In thisregime, a difference in the temperature gradient profile corrugations easily propagates to a transport level. Anotherimportant aspect is that the temperature gradient is kept
almost constant in this region, which is close to the criticalgradient level with R0 =Lti 6. The constrained profiles by
the nonlinear threshold as in Fig. 17(b) and the staircase
meso-scale organization are important characteristics of fullF simulations. As shown in Fig. 17(c), the radial electric field
determined by radial force balance shows a global V-shapewith a negative peak around the mid minor radius which isconsistent with the discussions in Ref. 16. The radial electricfield then determines the directions for avalanche propagation. The parallel flow velocity Uk in Fig. 17(d) represents
the intrinsic rotation. Although the levels of rotation speeds
are not exactly the same, both codes show co-rotation in thecore region with q ¼ r=a 0:2 and counter-rotation in theouter region with 0:5 q 1:0. The larger difference inthe core region could be explained by the difference in the
inner boundary condition, the Null (resp. natural) boundary
condition in GYSELA (resp. GT5D). Another important difference is the source term, where GYSELA imposes the102515-13

Asahi et al.

Phys. Plasmas 24, 102515 (2017)

FIG. 16. Spatio-temporal evolutions of the ion turbulent heat flux Qi =½ni Ti vti q2ti =a2  [(a) and (d)], the normalized temperature gradient R0 =Lti [(b) and (e)] andthe radial electric field shear dEr =dr=ðvti B0 =qti Þ [(c) and (f)] in GYSELA (upper row) and GT5D (lower row).

FIG. 17. The quasi-steady state radial
profiles (time average over tvti =R0¼ 300–750) of (a) the ion energy flux
Qi, (b) the normalized ion temperature
gradient R0 =Lti , (c) the radial electricfield Er, and (d) the parallel flow velocity Uk .

102515-14

Asahi et al.

Phys. Plasmas 24, 102515 (2017)

FIG. 18. The statistical analysis of the
turbulent heat flux Qi and R0 =Lti which
are evaluated for 300 tvti =R0 750
with radial average in source free
region (0:5 q 0:8). The frequency
spectra of (a) Qi and (b) R0 =Lti , andthe probability density function (PDF)
of (c) Qi and (d) R0 =Lti are shown.

momentum input in addition to the energy input (see discussions in Sec. II C for detail). It should also be noted that thehigher resolution of GYSELA near core could have some
impacts (the differences in the boundaries are detailed inSec.
Ð III C). Ð In the consequence, the net rotation, Uk; net¼ drrUk = drr, is slightly different, where GYSELA givesa large negative net rotation Uk; net =vti 0:0069 and GT5D
gives a small negative net rotation Uk; net =vti 0:0003.
The dynamics of rotation will be detailed with respect to themomentum transport processes in a separate publication.
Finally, statistical properties are assessed for 300 tvti =R0 750 with a radial average in the source free region(0:5 q 0:8). Figures 18(a) and 18(b) show the frequency
Fourier spectrum of ion turbulent energy flux Qi =QGB andnormalized temperature gradient R0 =Lti , with the gyro-Bohm
coefficient QGB ¼ ni Ti vti q2ti =a2 . The signal is Fourier transformed first at each radial point and then radially averaged.
As discussed in Ref. 38, the frequency Fourier spectrum
exhibits a 1/f decay at intermediate frequencies, which arefound in both codes [See Fig. 18(a)]. In addition, there
appears to be a transition from 1/f to 1/f3 decay at high frequencies. The transition frequency which may be interpreted
as the transition from avalanche-like dynamics to vortex
dynamics is the same in the two codes.
A 1/f decay is often considered as a typical feature of aself-organized criticality (SOC) like behavior, similar tothe sand pile automaton model by Bak et al.39 This agreement in the statistical properties indicates that both codes
exhibit a SOC-like behavior, which stems from the nonlinear interactions between the background profile and theplasma turbulence, which is the heart of full-F gyrokinetic
simulations. These findings are non-trivial since they

confirm that this SOC-like behavior is robust and notdependent on numerics.
The burst characteristics can also be extracted from theprobability density function (PDF) shown in Figs. 18(c) and18(d). The PDF is obtained from the histogram with anappropriate normalization so that the integral over the range
is equal to one. Although the PDF of normalized temperature
gradient R0 =Lti in Fig. 18(d) shows the Gaussian-like characteristic, the heat flux Qi in Fig. 18(c) shows the long-tails inboth codes indicating the non-Gaussian transport phenomena. It should still be noted that the heat flux stems fromcross-correlated fluctuating quantities (temperature andE  B velocity) so that it can show non-Gaussian PDF as
long as both temperature and E  B velocity exhibit
Gaussian PDF.40
VII. SUMMARY

In this work, we have performed benchmark tests ofglobal ITG turbulence with a fixed heat source. For this purpose, we employed two different full-F gyrokinetic codes
GYSELA and GT5D. The successful benchmark of full-F
global gyrokinetic codes gives confidence in the transportphenomena described by full-F models in the presence of a
fixed heat source. Since the flux-driven ITG turbulence
involves a wide range of physics, i.e., turbulence, neoclassical and their interplay, it is not easy to perform benchmarking straight away for a flux driven case. Thus, we firstdecomposed the problem into small subsets that correspondto idealized physical situations. Basically, the equilibrium
poloidal flow is described by the neoclassical physics andthe time evolution of turbulence is described by the gyrokinetic model. For the benchmarking of neoclassical physics,

102515-15

Asahi et al.

we choose the plasma parameters such that ITG modes arestable. For the benchmarking of turbulence, we investigatedthe linear ITG modes and decaying ITG turbulence in the
weak collisional regime. After getting a reasonable agreement, we carried out a flux-driven benchmark with theappropriate numerical and physical settings.
First, we carried out the benchmark of neoclassical theory and zonal flow damping tests. Neoclassical theory predicts the ion heat transport level and the mean poloidal
velocity. In order to compare the simulation results against
theory, we carried out a collisionality scan for a wide rangeof  ¼ 0:01  10:0 and confirmed that both codes give theion heat transport predicted by theory. In combination withthe radial force balance relation, we confirmed that the poloidal velocity agrees well with the theoretical estimate. The
short and long time behaviors of the zonal flows are also
compared against theory. For short time behavior, it is confirmed that the GAM frequency and the damping rate agreewith the theoretical estimate. The long time behavior called
the residual zonal flow is also investigated and comparedwith the Rosenbluth-Hinton theory. Both codes show theconvergence of the zonal flow velocity to the level predictedby theory. Thanks to the benchmarking above, we got readyto carry out the long-time flux-driven simulations with collisions. The long time zonal flow behavior is particularly
important since the flux-driven simulations require robust
long-time behaviors.
Second, we carried out benchmarks for linear ITG
modes and nonlinear decaying ITG turbulence to test the
gyrokinetic model. In the linear benchmark, we found good
agreements for the growth rates and real frequencies. Also,
we found a similar shape of the global eigenmodes. The
analysis for the nonlinear decaying ITG turbulence shows
good agreements in the critical temperature gradients.
Although pure collisionless situations would be preferable, aweak collision has been added in order to avoid unphysicalsolutions due to the absence of diffusion in vk direction,
whose impact on the turbulence is confirmed to be negligible
in the decaying turbulence case. It was also found that theradial structures of the zonal flow and the temperature gradient have quite similar features. In the asymptotic state, turbulence heat transport is almost quenched and the temperature
gradient is below the critical level.
Finally, a benchmarking of long-time flux-driven simulations has been performed. As described above, this kind of
simulation includes every single piece of physics benchmarked above. The spatio-temporal evolutions of ion turbulent heat flux showed an avalanche-like transport in bothcodes. Correspondingly, we found avalanche-like features inthe temperature gradient and electric field shear. We also
found an ion temperature profile constrained by the nonlinear
threshold of Lti value. In addition, we compared the radial
electric field to find a global V-shape with a negative peakaround the mid-minor radius. Although the amplitude is not
the same, we found co-current rotation near the core andcounter-current rotation near the edge in both codes. We
found similar statistical properties, like 1/f decay and its transition to 1/f 3 decay in both codes. These features can be considered as indicative of a self-organized criticality (SOC)

Phys. Plasmas 24, 102515 (2017)

behavior, which are the heart of the full-F flux driven gyrokinetic simulations. These simulations confirm that this SOClike behavior is robust and independent of numerics.
In addition to the presented comparisons, it may be useful to report the difficulties we encountered in this benchmarking work. Indeed, there are mainly three important
points: boundary conditions for the electrostatic potential,
initial conditions, and sink terms.
1. Boundary conditions are critical in the Rosenbluth-Hinton
test (see Sec. IV B). The radial wavenumber kr of zonal
mode can evolve in time which results in a time evolutionof GAM damping rate [see Eq. (2.10) in Ref. 31 for the
dependence of GAM damping rate on kr]. Since boundary
conditions have large impacts on the time evolution of kr,
it is important to give the appropriate boundary conditions
for the two codes.
2. Nonlinear simulations are sensitive to initial conditions.
Especially, the amplitude and time behavior of avalanchelike transport in flux-driven simulations are chaotic. This
complex dynamics is sensitive to initial conditions.33 In
addition, GYSELA and GT5D employ different numerical schemes so that we cannot expect an exact agreementin the spatio-temporal turbulence behavior. Nevertheless,one may expect the system to be insensitive to initial conditions after several confinement time, but it is generally
too costly and cannot be tackled in the present study.
Rather, we compare the statistical properties like PDF andfrequency spectra, for which we can expect an agreement
as discussed in Ref. 33. This is the reason why it is important to analyse statistical properties in flux-driven simulation benchmarking. The decaying turbulence is also
largely affected by initial conditions, where we compare
well-defined properties like the critical temperature gradient and the relaxation process toward the asymptotic
state.
3. The sink model implementation is also a matter.
Comparing two different flux-driven simulations with diffusion or Krook operators by GYSELA, it turned out thatthese operators give different toroidal flow and radial
electric field profiles. In fact, a Krook-type sink termworks not only as a heat sink but also work as a momentum source (sink) term when a negative (positive)
momentum flux is absorbed. Therefore, using the same
sink operator is also quite important to get the similar
toroidal flow profile. As a matter of fact, much effort isnow devoted to the development of physically more relevant boundaries, including the modeling of the open field
line region (Scrape-Off Layer).23
As is clear, all the points presented above should beappropriately treated for the benchmark of flux-driven
simulations.
ACKNOWLEDGMENTS

All the simulations were performed on Helios at IFERC
and Marconi at CINECA. The authors thank Ch. Passeron
for technical support on performing large-scale simulations
on these supercomputers. This work was partly supported by102515-16

Asahi et al.

Phys. Plasmas 24, 102515 (2017)

the MEXT, Grant for HPCI Strategic Program Field No. 4:
Next-Generation Industrial Innovations, and Grant for PostK priority issue No. 6: Development of Innovative Clean
Energy. This work has also received funding from theEuropean Union’s Horizon 2020 research and innovation
programme under Grant Agreement No. 633053. The views
and opinions expressed herein do not necessarily reflect
those of the European Commission.
APPENDIX: COMPARISON OF NUMERICAL COSTS

It is good to know the difference in numerical costs ofthese two codes to understand more about the algorithm efficiencies, i.e., Semi-Lagrangian scheme in GYSELA and
Finite-Difference scheme in GT5D. In order to consider theefficiency from the view point of time-to-solution, we
employ Lattice Update Per Second (LUPS) as a metric rather
than Floating-point Operations Per Second (FLOPS). LUPS
is computed by the number of grids divided by an elapsed
time for updating each grid point.
Both GYSELA and GT5D codes are highly parallelized
and based on a hybrid MPI/OpenMP programming
model.41,42 The parallelization strategies are to some extent
similar where both codes define a MPI communicator foreach value of the magnetic moment l. This strategy reflects
the underlying physics that the magnetic moment l is aparameter in Vlasov equation if we dismissed the collisional
effects [see Eqs. (1)–(3)]. Thus, the number of sub-domains
in l direction pl is equivalent to Nl the mesh size of l direction. Inside each l communicator, the two dimensional
domain decomposition is applied to the poloidal planeðNx ; Ny Þ. If we call the number of sub-domains in (x, y)
directions ðpx ; py Þ, the original domain ðNx ; Ny Þ is decomposed into a sub-domain sized ðNx =px ; Ny =py Þ. Here, (x, y)
represents (r, h) in GYSELA and (R, Z) in GT5D as
described in Sec. III C. The total number of MPI processes
can be written as NMPI ¼ px  py  Nl .
To investigate the strong scaling, we scan with the number of sub-domains ranged ðpr ; ph Þ ¼ ð1–4; 1–4Þ in
GYSELA and ðpR ; pZ Þ ¼ ð2–8; 2–8Þ in GT5D. The subdomain in l direction and the number of threads ðNl ; Nthread Þ
are kept constant with ð16; 18Þ in GYSELA and ð16; 4Þ in
GT5D. As a testbed, we employ the Intel Xeon Processor
E5–2697 V443 (Broadwell family) equipped with 18 cores.
Hyperthreading is switched off in this test. The different configuration of GYSELA and GT5D in the number of cores per
CPU is coming from the measured performance of GT5D
(GYSELA uses 18 cores per CPU, while GT5D uses 16
cores per CPU). It is also possible to use the same parallelization in GT5D as GYSELA, but we found a slightly higher
performance with 16 cores per CPU (more MPI processes
and less OpenMP threads). This can be partially explainedby a higher cache locality with the MPI-domain decomposition inside a CPU. It should be noted that our codes havebeen optimized for Sandy Bridge architecture, not for
Broadwell architecture which is used in this study.
Therefore, the results shown here cannot be considered as
fully optimized, but moderately optimized. Figure 19 shows
the strong and weak scaling of both codes obtained from the

FIG. 19. Mega Lattice Update Per Seconds (MLUPS) of both codes. For
problem sizes, we employed the small one for q1
 ¼ 100 with
ðNx ; Ny ; Nu ; Nvk ; Nl Þ ¼ ð128  128  32  128  16Þ and the middle one
for q1
 ¼ 150 with ðNx ; Ny ; Nu ; Nvk ; Nl Þ ¼ ð256  256  32  128  16Þ.
For a parallelization configuration, we employ ðpr ; ph ; Nl ; Nthread Þ¼ ð1–4; 1–4; 16; 18Þ for GYSELA and ðpR ; pZ ; Nl ; Nthread Þ ¼ ð2–8; 2–8;
16; 4Þ for GT5D. One node includes two Broadwell CPUs.
1flux-driven simulations with q1
 ¼ 100 and q ¼ 150 on
44
the Marconi supercomputer. The input power is set as11.33 MW for q1
 ¼ 100 and 2 MW for q ¼ 150, which
also scales with plasma size in order to keep the plasma
parameter as close as possible. The q1
 ¼ 100 case isresolved with a mesh size ðNx ; Ny ; Nu ; Nvk ; Nl Þ ¼ ð128
128  32  128  16Þ and the q1
 ¼ 150 case is resolved
with a mesh ðNx ; Ny ; Nu ; Nvk ; Nv? Þ ¼ ð256  256  32  128
16Þ. The time step size is set as DtXi ¼ 10 in GYSELA
and 5 in GT5D for each case. The elapsed time for Lattice
Update is measured by the total elapsed time excluding diagnostics and initialization divided by the number of iterations.
In GYSELA, the 4D advection solver occupies about
65% of the overall execution cost, which consists of the 2D(r, h) þ 1D ðvk Þ þ 1D ðuÞ interpolation (50%) and the transpose communication (15%). The most time-consuming partof the rest is the collision operator which corresponds toabout 20% of the total cost. For the interpolation, GYSELA
employs the cubic spline interpolation method, which has agood conservation property at the cost of non-locality. Dueto the non-local feature, a single MPI process should have anaccess to all the values of a distribution function f in adimension where the interpolation takes place. Therefore, forthe 2D interpolation in (r, h) directions, a transpose communication is required to construct the original (r, h) domainfrom each sub-domain in (r, h) directions kept by a single
MPI process.
In GT5D, the iterative solver for linear advection termcosts about 90% of the total elapsed time, which is solved
implicitly with the 4th-order 4D Finite-Difference scheme.
The collision operator costs about 5% of the total cost. The
asymmetric block diagonal matrix is solved with the generalized conjugate residual (GCR) method. Since the matrix
solver requires several iterations, this is the most time
consuming part, which consists of the 4th-order 4D

102515-17

Asahi et al.

Finite-Difference computation (35%) and the vector computations (55%) such as DAXPY and DDOT. The FiniteDifference part includes the halo-communication (5%) andthe vector computation part includes the inter-node allreduce communication (5%).
As shown in Fig. 19, it is found that GYSELA is fasterthan GT5D about 30%–40% for small number of nodes and
15% for 128 nodes. This represents the better overall performance of GYSELA and better scalability of GT5D. The better
performance of GYSELA can be interpreted as higher computational intensity of the advection solver, which could be an
important advantage for Semi-Lagrangian method. Higher
computational intensity could be especially important on the
recent architectures like Broadwell which have a high computational performance and a relatively low memory bandwidth.
1
The GT5D weak scaling from q1
 ¼ 100 to q ¼ 150 showsa super scaling, which stems from the difference in the numberof iterations required by GCR method. It should be noted thatthe number of iterations cannot be fixed and different plasma
parameters result in different numbers of iterations. The better
scalability of GT5D can be explained by communicationcomputation overlapping based on OpenMP42 (allocating themaster thread to communication and others to computation),
which is applied to hide the halo-communication cost in the
Finite-Difference part. In order to reduce the inter-node allreduce communication cost, communication avoiding algorithms45 could be applied in the future. Similar strategy forcommunication-computation overlapping could be applied to
GYSELA, by pipelining the advection part to overlap the small
chunks of transpose communications and advection
computations.
1

X. Garbet, Y. Idomura, L. Villard, and T. Watanabe, Nucl. Fusion 50,
043002 (2010).
2
A. M. Dimits, G. Bateman, M. A. Beer, B. I. Cohen, W. Dorland, G. W.
Hammett, C. Kim, J. E. Kinsey, M. Kotschenreuther, A. H. Kritz, L. L.
Lao, J. Mandrekas, W. M. Nevins, S. E. Parker, A. J. Redd, D. E.
Shumaker, R. Sydora, and J. Weiland, Phys. Plasmas 7, 969–983 (2000).
3
W. M. Nevins, J. Candy, S. Cowley, T. Dannert, A. Dimits, W. Dorland,
C. Estrada-Mila, G. W. Hammett, F. Jenko, M. J. Pueschel, and D. E.
Shumaker, Phys. Plasmas 13, 122306 (2006).
4
G. Dif-Pradalier, P. H. Diamond, V. Grandgirard, Y. Sarazin, J. Abiteboul,
X. Garbet, P. Ghendrih, A. Strugarek, S. Ku, and C. S. Chang, Phys. Rev.
E 82, 025401 (2010).
5
G. Dif-Pradalier, G. Hornung, X. Garbet, P. Ghendrih, V. Grandgirard, G.
Latu, and Y. Sarazin, Nucl. Fusion 57, 066026 (2017).
6
S. Ku, C. Chang, and P. Diamond, Nucl. Fusion 49, 115021 (2009).
7
F. Jenko, Comput. Phys. Commun. 125, 196–209 (2000).
8
S. Jolliet, A. Bottino, P. Angelino, R. Hatzky, T. Tran, B. Mcmillan, O.
Sauter, K. Appert, Y. Idomura, and L. Villard, Comput. Phys. Commun.
177, 409–425 (2007).
9
V. Grandgirard, Y. Sarazin, P. Angelino, A. Bottino, N. Crouseilles, G.
Darmet, G. Dif-Pradalier, X. Garbet, P. Ghendrih, S. Jolliet, G. Latu, E.
Sonnendr€ucker, and L. Villard, Plasma Phys. Controlled Fusion 49, B173
(2007).
10
Y. Idomura, M. Ida, T. Kano, N. Aiba, and S. Tokuda, Comput. Phys.
Commun. 179, 391–403 (2008).
11
G. L. Falchetto, B. D. Scott, P. Angelino, A. Bottino, T. Dannert, V.
Grandgirard, S. Janhunen, F. Jenko, S. Jolliet, A. Kendl, B. F. McMillan,
V. Naulin, A. H. Nielsen, M. Ottaviani, A. G. Peeters, M. J. Pueschel, D.
Reiser, T. T. Ribeiro, and M. Romanelli, Plasma Phys. Controlled Fusion
50, 124015 (2008).
12
X. Lapillonne, B. F. McMillan, T. G€
orler, S. Brunner, T. Dannert, F.
Jenko, F. Merz, and L. Villard, Phys. Plasmas 17, 112321 (2010).

Phys. Plasmas 24, 102515 (2017)13

V. Grandgirard, J. Abiteboul, J. Bigot, T. Cartier-Michaud, N. Crouseilles,
G. Dif-Pradalier, C. Ehrlacher, D. Estève, X. Garbet, P. Ghendrih, G. Latu,
M. Mehrenberger, C. Norscini, C. Passeron, F. Rozar, Y. Sarazin, E.
Sonnendr€
ucker, A. Strugarek, and D. Zarzoso, Comput. Phys. Commun.
207, 35–68 (2016).
14
Y. Idomura, J. Comput. Phys. 313, 511–531 (2016).
15
A. J. Brizard and T. S. Hahm, Rev. Mod. Phys. 79, 421–468 (2007).
16
Y. Idomura, H. Urano, N. Aiba, and S. Tokuda, Nucl. Fusion 49, 065029(2009).
17
S. Satake, Y. Idomura, H. Sugama, and T.-H. Watanabe, Comput. Phys.
Commun. 181, 1069–1076 (2010).
18
G. Dif-Pradalier, P. H. Diamond, V. Grandgirard, Y. Sarazin, J. Abiteboul,
X. Garbet, P. Ghendrih, G. Latu, A. Strugarek, S. Ku, and C. S. Chang,Phys. Plasmas 18, 062309 (2011).
19
D. Estève, X. Garbet, Y. Sarazin, V. Grandgirard, T. Cartier-Michaud, G.
Dif-Pradalier, P. Ghendrih, G. Latu, and C. Norscini, Phys. Plasmas 22,122506 (2015).
20
D. Estève, Y. Sarazin, X. Garbet, V. Grandgirard, S. Breton, P. Donnel, Y.
Asahi, C. Bourdelle, G. Dif-Pradalier, C. Ehrlacher, C. Emeriau, P.
Ghendrih, C. Gillot, G. Latu, and C. Passeron, “Self-consistent gyrokinetic
modeling of neoclassical and turbulent impurity transport,” Nucl. Fusion(published online).
21
A. Lenard and I. B. Bernstein, Phys. Rev. 112, 1456–1459 (1958).
22
X. Q. Xu and M. N. Rosenbluth, Phys. Fluids B: Plasma Phys. 3, 627–643(1991).
23
G. Dif-Pradalier, E. Caschera, P. Ghendrih, Y. Asahi, P. Donnel, X.
Garbet, V. Grandgirard, G. Latu, C. Norscini, and Y. Sarazin, Plasma
Fusion Res. 12, 1203012 (2017).
24
G. Strang, SIAM J. Numer. Anal. 5, 506–517 (1968).
25
Y. Morinishi, T. Lund, O. Vasilyev, and P. Moin, J. Comput. Phys. 143,
90–124 (1998).
26
L. Ming-Chih, Numer. Methods Partial Diff. Eqs. 17, 199–203 (2001).
27
N. Bouzat, C. Bressan, V. Grandgirard, G. Latu, and M. Mehrenberger,“Targeting realistic geometry in Tokamak code Gysela,” Proceedings ofthe ESAIM (to be published).
28
C. S. Chang and F. L. Hinton, Phys. Fluids 25, 1493–1494 (1982).
29
F. L. Hinton and R. D. Hazeltine, Rev. Mod. Phys. 48, 239–308 (1976).
30
M. N. Rosenbluth and F. L. Hinton, Phys. Rev. Lett. 80, 724–727(1998).
31
H. Sugama and T.-H. Watanabe, J. Plasma Phys. 72, 825–828 (2006).
32
C. Greenfield, J. DeBoo, T. Osborne, F. Perkins, M. Rosenbluth, and D.
Boucher, Nucl. Fusion 37, 1215 (1997).
33
S. Jolliet and Y. Idomura, Nucl. Fusion 52, 023026 (2012).
34
F. Rath, A. G. Peeters, R. Buchholz, S. R. Grosshauser, P. Migliano, A.
Weikl, and D. Strintzi, Phys. Plasmas 23, 052309 (2016).
35
A. G. Peeters, F. Rath, R. Buchholz, Y. Camenen, J. Candy, F. J. Casson,
S. R. Grosshauser, W. A. Hornsby, D. Strintzi, and A. Weikl, Phys.
Plasmas 23, 082517 (2016).
36
M. Kikuchi and M. Azumi, Rev. Mod. Phys. 84, 1807–1854 (2012).
37
B. F. McMillan, S. Jolliet, T. M. Tran, L. Villard, A. Bottino, and P.
Angelino, Phys. Plasmas 16, 022310 (2009).
38
Y. Sarazin, V. Grandgirard, J. Abiteboul, S. Allfrey, X. Garbet, P.
Ghendrih, G. Latu, A. Strugarek, G. Dif-Pradalier, P. Diamond, S. Ku, C.
Chang, B. McMillan, T. Tran, L. Villard, S. Jolliet, A. Bottino, and P.
Angelino, Nucl. Fusion 51, 103023 (2011).
39
P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59, 381–384 (1987).
40
B. A. Carreras, C. Hidalgo, E. Sanchez, M. A. Pedrosa, R. Balbın, I.
Garcıa-Cortes, B. van Milligen, D. E. Newman, and V. E. Lynch, Phys.
Plasmas 3, 2664–2672 (1996).
41
J. Bigot, V. Grangirard, G. Latu, C. Passeron, F. Rozar, and O. Thomine,in Proceedings of the ESAIM (2013), Vol. 43, pp. 117–135.
42
Y. Idomura, M. Nakata, S. Yamada, M. Machida, T. Imamura, T.
Watanabe, M. Nunami, H. Inoue, S. Tsutsumi, I. Miyoshi, and N. Shida,Int. J. High Performance Comput. Appl. 28, 73–86 (2014).
43
See http://ark.intel.com/products/91755/Intel-Xeon-Processor-E5-2697v4-45M-Cache-2_30-GHz for IntelV XeonV Processor E5-2697 v4 (45M
Cache, 2.30 GHz).
44
See
http://www.hpc.cineca.it/news/marconi-new-cineca-supercomputerfor MARCONI: the new Cineca supercomputer.
45
M. Hoemmen, “Communication-avoiding Krylov subspace methods,”
Ph.D. thesis (EECS Department, University of California, Berkeley,2010).
R

R