Plasma and Fusion Research: Regular Articles

Volume 9, 1403029 (2014)

Local Gyrokinetic Vlasov Simulations with Realistic
Tokamak MHD Equilibria
Motoki NAKATA, Akinobu MATSUYAMA, Nobuyuki AIBA, Shinya MAEYAMA,
Masanori NUNAMI1) and Tomo-Hiko WATANABE2)
Japan Atomic Energy Agency, Rokkasho, Aomori 039-3212, Japan
1)
National Institute for Fusion Science, Toki 509-5292, Japan2)
Department of Physics, Nagoya University, Nagoya 464-8602, Japan(Received 15 January 2014 / Accepted 21 February 2014)
A local gyrokinetic Vlasov simulation code GKV is extended to incorporate realistic tokamak equilibria
including up-down asymmetry, which are produced by a free-boundary 2D Grad-Shafranov equation solver
MEUDAS. By using a newly developed interface code IGS, two dimensional rectangular equilibrium data fromMEUDAS is converted to straight-field-line flux coordinates such as Hamada, Boozer, and axisymmetric coordinates, which are useful for gyrokinetic micro-instability and turbulent transport analyses. The developed codes
have been verified by a cross-code benchmark test using Cyclone-base-case like MHD equilibrium, where good
agreement in the dispersion relation of ion temperature gradient (ITG) driven mode has been confirmed. The
extended GKV is applied to two types of shaped plasmas expected in JT-60SA tokamak devices, i.e., ITER-likeand highly-shaped plasmas, and ITG-mode stability and residual zonal-flow level are investigated. Through thedetailed comparisons, more favorable stability properties against the ITG mode are revealed for the highly-shaped
case, where the lower ITG-mode growth rate and higher residual zonal-flow levels compared to the ITER-likecase are identified.
c 2014 The Japan Society of Plasma Science and Nuclear Fusion ResearchKeywords: gyrokinetic Vlasov simulation, realistic geometry, ITG mode
DOI: 10.1585/pfr.9.1403029

1. Introduction
Turbulent transport of heat, particle, and momentum
is a dominant transport channel in magnetically confinedplasmas [1]. Understanding physics behind the turbulent
transport and predicting quantitatively the resultant transport levels are, thus, important issues for the establishmentof self-ignited steady burning plasmas in future fusion devices such as ITER and DEMO reactors. To this end, ion
temperature gradient (ITG) driven turbulence, which is expected as one of the main causes of ion heat transport, has
extensively been studied by means of gyrokinetic simulations [e.g., see Ref. [2]]. A number of important resultson the ITG turbulent transport processes have been obtained so far, and the ion heat transport suppression by
nonlinearly generated zonal flows is one of the most important findings [3, 4]. In particular, local fluxtube gyrokinetic model [5], which is utilized for investigation of local micro-instabilities, fine-scale turbulence structures andthe related transport levels and mechanisms, is nowadays afundamental tool not only for theoretical studies, but also
for experimental analyses.
In view of turbulence control, the plasma shaping is
an eﬀective method to improve the confinement performance in tokamaks. In addition to MHD stability, microauthor’s e-mail: nakata.motoki@jaea.go.jp

instabilities, turbulent transport, and zonal-flow properties
can be strongly aﬀected by the plasma shape, such aselongation, triangularity and up-down asymmetry, throughthe change of the linear frequencies, mode-structures,and trapped/passing guiding-center orbits. There aresome earlier works looking into shaping eﬀects on microinstabilities and turbulent transport. Miller’s analytic local equilibrium model [6] is widely used in the local fluxtube model, then the fundamental dependencies with respect to the elongation and triangularity parameters havebeen investigated [7–9]. The impact of negative triangularity on electron heat transport has also been investigated by using both global and local turbulence simulations [10, 11], where the change of electron precession
drifts leads to the stabilization of trapped electron modes(TEM). Many eﬀorts have also been devoted for the crosscode benchmark [12] and validation against experimental
observations [13, 14].
The plasma shaping can also aﬀect the zonal-flow dynamics. Especially, an extension of linear zonal-flow theory [15] has revealed that the residual zonal-flow level isenhanced due to the plasma elongation [16]. The residual
zonal-flow theory including the plasma shaping has also
been applied to gyrofluid simulations [17] based on kinetic
closure model [18].
In our previous studies [19, 20] using a local fluxtube

1403029-1

c 2014 The Japan Society of Plasma
Science and Nuclear Fusion Research

Plasma and Fusion Research: Regular Articles

Volume 9, 1403029 (2014)

gyrokinetic Vlasov code GKV [21], a critical role of zonal
flows leading to successive nonlinear entropy transfer processes in the regulated ITG turbulent transport has been
elucidated, where a simple concentric circular tokamak
geometry was considered. In this study, a realistic tokamak equilibrium calculated by a free-boundary 2D GradShafranov equation solver MEUDAS [22] is newly incorporated to GKV using an interface code IGS, in order toextend the linear and nonlinear analyses to the tokamak
experimental conditions. The extension is also indispensable for the quantitative validation against experimental results and for the prediction studies towards future fusion
devices. By using the extended code, linear gyrokinetic
simulations are performed for two types of equilibria, i.e.,
ITER-like and highly-shaped cases, which are expected in
JT-60SA tokamak device, where the MHD stability has already been well investigated [23]. Then, the ITG-mode
stability and the residual zonal-flow levels are comparedin detail. These linear results provide one with fundamental insights for the later nonlinear gyrokinetic simulation
studies towards JT-60SA and ITER.
The rest of the paper is organized as follows. In Sec. 2,
the interface code IGS is presented, as well as some basic
introductions to flux coordinate systems. Then, calculation model in GKV and verification results are given inSec. 3. In Sec. 4, linear simulation results on ITG-mode
stability and residual zonal flows for JT-60SA plasmas arepresented. Finally, the summary is given in Sec. 5.

2. Realistic Tokamak Geometry
All the flux coordinates considered here are constructed by 2D data obtained from a free-boundary GradShafranov equation solver MEUDAS [22], where general
up-down asymmetric tokamak equilibria with the divertor
legs are calculated consistently with the realistic externalcoil configurations. The interface code connecting the flux
coordinates to the fluxtube one used in local gyrokinetic
simulations is briefly presented in this section.

2.1

Magnetic field representation

General expressions for the magnetic field B and thecurrent density j, as well as the ideal MHD force balance
relation, in a flux coordinate system (ρ, θ, ζ) are given by[24]
B = Ψ  ∇ζ × ∇ρ + Φ ∇ρ × ∇θ + ∇ρ × ∇α̃= β̃∇ρ + I∇θ + J∇ζ + ∇γ̃ ,4π
j = J  ∇ρ × ∇ζ + I  ∇ρ × ∇θ + ∇β̃ × ∇ρ ,cp = −I  B · ∇θ − J  B · ∇ζ + B · ∇β̃ ,B as in Eq. (1), the covariant one is shown in Eq. (2). The
toroidal magnetic flux, the poloidal one, and the toroidal
current inside the magnetic flux surface ρ are denoted by
2πΦ, 2πΨ and 2πI, respectively, while 2πJ represents
the poloidal current outside ρ. Ensuring physical singlevaluedness for B and j, the quantities α̃(ρ, θ, ζ), β̃(ρ, θ, ζ),and γ̃(ρ, θ, ζ) must be doubly periodic functions with respect to θ and ζ. It should be noted that Eqs. (1) – (4) areof invariance under the following transformation,α̃∗ = α̃ + g(ρ) ,∗β̃ = β̃ + h (ρ) ,∗γ̃ = γ̃ − h(ρ) ,where ρ is the radial coordinate representing a magnetic
flux surface, and θ and ζ are the poloidal- and toroidalangle coordinates, respectively. The prime symbol meansthe ρ-derivative: X  = dX(ρ)/dρ, and c denotes the speedof light. In addition to the contravariant representation for

(6)
(7)

0

where V  (ρ) is the radial derivative of the plasma volume
at ρ. Solvability constraints of the magnetic diﬀerential
equations, i.e., B · ∇γ̃ = 0 and B · ∇β̃ = 0 in Eqs. (1),(2) and (4), provide one with

 
Φ 
(11)
B2 = 4π2  J + q−1 I ,VΦ p = −4π2  J  + q−1 I  ,(12)V
where q(ρ) = Φ /Ψ  denotes the safety factor profile.
Although the above relations are consequence of general discussions, irrespective of the choice of (ρ, θ, ζ), one
can construct more specific coordinate systems, which areuseful for theoretical and/or numerical analyses. Following
Pustovitov’s method [25], the most important three fluxcoordinate systems, i.e., Hamada, Boozer, and axisymmetric (or natural) coordinates, are classified systematically.
First, from the definition of the flux-surface average[Eq. (9)], one finds that the Jacobian can be written by anarbitrary function f = f (ρ, θ, ζ),
√

(2)

(4)

(5)

where g and h are arbitrary flux functions depending onlyon ρ, and note that β̃ and γ̃ should be simultaneously trans√
formed by h. The Jacobian gρθζ and the flux-surfaceaverage operator · · ·  are defined as√gρθζ = (∇ρ × ∇θ · ∇ζ)−1 ,
(8)
 2π
1
√
dθdζ gρθζ X(ρ, θ, ζ) ,(9)
X = 
V
0
 2π
√
V =
dθdζ gρθζ ,
(10)

(1)

(3)



4π2 f
1
.
= 
gρθζ
V f

(13)

Note that the physical dimension of f is irrelevant in theabove definition, and this is significant feature to define
various coordinates in the same manner. The above relation is applied to Eq. (1), then we find1403029-2

Ψ f,V  f Φ f(B − ∇ρ × ∇α̃) · ∇ζ = 4π2 .
V f(B − ∇ρ × ∇α̃) · ∇θ = 4π2

(14)
(15)

Plasma and Fusion Research: Regular Articles

Volume 9, 1403029 (2014)

Equations (14) and (15) are coupled equations with respect√
to four variables, α̃, f (or gρθζ ), θ, and ζ, for given magnetic field B and the radial coordinate ρ. When any two ofthem are given, one can obtain the corresponding flux coordinate system determined by Eqs. (14) and (15). Then,
Hamada, Boozer, and the axisymmetric coordinates areclassified as follows:
Hamada (ρ, θH , ζH ) : α̃ = 0, f = Const. ,(16)

Boozer (ρ, θB , ζB ) : α̃ = 0, f = B ,

(17)

2

Axisymmetric (ρ, θA , ζA ) : α̃ = 0, f = B·∇ζA , ζA = φ , (18)where α̃ = 0 makes the magnetic field line straight, i.e.,dζ/dθ = B · ∇ζ/B · ∇θ = q(ρ), then such coordinates arecalled the straight-field-line (SFL) flux coordinates. Rigorously speaking, α̃ = s(ρ) (s is an arbitrary flux function)
also gives the SFL flux coordinates, but one can always sets(ρ) to be 0 by Eq. (5). Note also that the second condition in Eq. (18), f = B · ∇ζA , is actually equivalent to theflux-surface average of Eq. (15), and does not provide any
relevant equations for determining ζ. Thus, ζA = φ shouldbe imposed as an additional condition for the axisymmetric coordinates, where φ denotes the azimuthal angle in thecylindrical coordinates (R, Z, φ). Each coordinate system
possesses diﬀerent advantages, and their detailed properties and the components of B and j are summarized in Appendix A.
The coordinate transformation from one to another isgiven by the generating function G(ρ, θ, ζ):θ̄ = θ + Ψ G(ρ, θ, ζ) ,ζ̄ = ζ + Φ G(ρ, θ, ζ) ,the cubic-spline interpolation is applied to the 2D rectangular data of Ψ (R, Z). Then, the axisymmetric coordinates, R = R(ρ, θA ), Z = Z(ρ, θA ), ζA = φ, are first constructed from Ψ (R, Z) by solving dR/dlp = |∇Ψ |−1 ∂Ψ/∂Z,
dZ/dlp = −|∇Ψ |−1 ∂Ψ/∂R, dθA /dlp = J/qR|∇Ψ |, where dlp
is the arc-length element on the poloidal cross section. After that, Hamada and Boozer coordinates, R = R(ρ, θH,B ),Z = Z(ρ, θH,B ), ζH,B = φ + ΦG H,B (ρ, θA ) are obtained from
Eqs. (19) – (21).
The covariant metric tensor, gi j is calculated by thefollowing relation:
gi j =

∂φ ∂φ
∂R ∂R
∂Z ∂Z
+
+ R2
,
∂ui ∂u j ∂ui ∂u j∂ui ∂u j(23)

where (ui , u j ) = {ρ, θ, ζ} (the coordinates label A, B, andH are omitted). The determinant of gi j gives the Jacobian√gρθζ = det(gi j )1/2 , and the contravariant metric tensor gi j
is calculated by

1 
(24)
gil = √ 2 g jm gkn − g jn gkm ,gρθζ
where {i, j, k} and {l, m, n} are the even permutations of{ρ, θ, ζ}.
The examples of constructed flux coordinates areshown in Fig. 1. Here, the Solov’ev analytic equilibrium

(19)
(20)

where G is determined by the following magnetic diﬀerential equation,
B · ∇G(ρ, θ, ζ) = √

1
1
− √
.
gρθ̄ζ̄gρθζ

(21)

For instance, G(ρ, θB , ζB ) giving the transformation from
Hamada to Boozer coordinates is 
B2
G(ρ, θB , ζB ) =    β̃∗ (ρ, θB , ζB ) . (22)Φ p J + q−1 I

2.2

Numerical constructions of SFL flux coordinates

The free-boundary Grad-Shafranov equation solver
MEUDAS provides one with 2D data of Ψ (R, Z) for givenp(Ψ ) and J(Ψ ). These fundamental data are, then, input to a newly developed interface code IGS, which produces Hamada, Boozer, and the axisymmetric coordinates
numerically. The basic idea in IGS is similar to that inRef. [26]. In addition to the original data format in MEUDAS, IGS can also treat the G-EQDSK format which is
widely used in the fusion physics field. The radial coordinate is chosen as ρ = (Ψ/Ψedge )1/2 or ρ = (Φ/Φedge )1/2 ,where Ψedge and Φedge are the poloidal and toroidal magnetic fluxes at the plasma edge, respectively. In IGS,

Fig. 1 Comparison of the SFL coordinates: (a) the original
2D data Ψ (R, Z) for the Solov’ev analytic equilibriumwith the up-down asymmetry, and (b) axisymmetric, (c)Boozer, (d) Hamada coordinates constructed by IGS.

1403029-3

Plasma and Fusion Research: Regular Articles

Volume 9, 1403029 (2014)

with up-down asymmetry is considered:
1 (1 − D)R̄2 Z̄ 2 − DZ̄ 2 1 2Ψ= 2+ (R̄ − 1)2
2
Ψedge
4
E
a
νasym 2+(R̄ − 1)Z̄ ,(25)
awhere R̄ and Z̄ are the positions normalized by the major
radius on the magnetic axis Rax , and a = a/Rax , E, D, and,
νasym are shape parameters related to the inverse aspect ratio, elongation, triangularity, and up-down asymmetry, respectively. One finds that the iso-angle surfaces of θ in
Hamada coordinates are less deformed in comparison withthose in Boozer and the axisymmetric ones. Instead, in the
non-axisymmetric systems, Hamada coordinates show relatively strong deformation of iso-angle surfaces of ζ whichleads to the spectral broadening on B in the toroidal mode
numbers [24].
A numerical accuracy of the constructed flux coordinates is verified by comparisons with the analytic expression. Figure 2 shows the absolute errors indicating the deviation from the following analytic relations,√

1 ∂ρ[R, Z]
1 ∂R[ρ, θ]= −,gρθζ ∂θR ∂Zgρρ [ρ, θ] =2
∂ρ[R, Z]
∂ρ[R, Z]
∇R +
∇Z ,
∂R
∂Z

One can see that all the constructed flux coordinates
indicate the numerical accuracy up to ∼ 10−8 . Besides,
the consistency among the metric tensor and magnetic fieldcomponents, e.g., Bθ = Bθ gθθ + Bζ gζθ , have been confirmedin IGS.

2.3 Local fluxtube coordinate system
Local fluxtube coordinates (x, y, z) used in gyrokinetic
simulations are based on the SFL flux coordinates. The
definition of each coordinate is given as
x = c x (ρ − ρ0 ), y = cy q(ρ)θ − ζ , z = θ ,(28)

where (ρ, θ, ζ) are arbitrary SFL flux coordinates, and c x ,cy , and ρ0 are the normalization factors and the reference
radial position, respectively. Using the above definitions,the magnetic field B is written as−1
B = cB ∇x × ∇y, cB = Ψ  c−1
x cy .

(29)

The partial-derivative operators with respect to x, y, and z
are given by1 ∂q ŝθ ∂∂=+,∂xc x ∂ρ c x ρ ∂ζ
∂
1 ∂
= −
,
∂y
cy ∂ζ
∂
∂
∂
=
+q ,
∂z
∂θ
∂ζ

(26)
(27)

where ρ = (Ψ/Ψedge )1/2 is considered here.

(30)
(31)
(32)

where ŝ = ρq /q denotes the magnetic shear. Note that∂/∂y = 0 for axisymmetric plasmas with ∂/∂ζ = 0. One
can also obtain the contravariant metric tensor and the Jacobian as follows:g xx = ∇x·∇x = c2x gρρ ,q ŝθ ρρg xy = ∇x·∇y = c x cy
g + qgρθ − gρζ ,ρ
g xz = ∇x·∇z = c x gρθ ,⎡⎢ q ŝθ 2 ρρyy
2⎢g + q2 gθθ + gζζ
g = ∇y·∇y = cy ⎢⎢⎢⎣
ρ
 2 


2q ŝθ ρθ
2q ŝθ ρζθζ
+,g − 2qg −gρρ
q ŝθ ρθgyz = ∇y·∇z = cy
g + qgθθ − gθζ ,ρ
gzz = ∇z·∇z = gθθ ,√−1 √g xyz = (∇x × ∇y · ∇z)−1 = c−1gρθζ .
x cy

(33)
(34)
(35)

(36)
(37)
(38)
(39)

The wavenumber vector perpendicular to the magnetic
field and the diﬀerential operator along the field line aredefined by

Fig. 2 Comparison of the absolute errors in (a) Eq. (26) and(b) Eq. (27), for axisymmetric, Hamada, and Boozer coordinates, where the radial distributions for θA,H,B ={0, π/2, π, 3π/2} are plotted in each coordinate system.

k⊥ = k x ∇x + ky ∇y ,(40)k⊥2 = k2x g xx + 2k x ky g xy + ky2 gyy ,∂cB,b·∇ = √
B g xyz ∂z

(41)
(42)

where b is the unit vector parallel to the magnetic field.

1403029-4

Plasma and Fusion Research: Regular Articles

Volume 9, 1403029 (2014)

As will be shown in Sec. 3, the radial coordinate ofρ = (Φ/Φedge )1/2 and c x = a, cy = aρ0 /q(ρ0 ), cB = Bax
are basically used in GKV, where Bax denotes the magnetic
field intensity on the magnetic axis of ρ = 0, and the plasma
minor radius is defined by a = (2Φedge /Bax )1/2 . Note that,in the local fluxtube formulation shown below, all the normalized quantities and operators do not depend explicitly
on the plasma minor radius a, but on the inverse aspectratio a . Although IGS gives us the complete metric information of realistic equilibria from MEUDAS, the analytic
geometry model is still useful to subtract the shaping effects. The analytic expression of the concentric circularequilibrium model and the simplified s−α one are summarized in Appendix B.

the triad-interaction condition of k⊥ = k⊥ + k⊥ . Collisionaleﬀects are introduced in terms of a linearized model collision operator Cs .
The magnetic and diamagnetic drift frequencies, ωDs
and ω∗T s , are given by
cT sk⊥ · b × μ∇B + ms 2 b·∇b
ωDs =
es B

c(ms 2 +μB) =K x k x + Ky ky ,
(44)
es cB
⎡
⎞⎤
⎛
⎜⎜⎜ ms 2 +2μB 3 ⎟⎟⎟⎥⎥⎥
cT s ⎢⎢⎢⎢⎜
ω∗T s =− ⎟⎟⎠⎥⎥⎦ k⊥ · b × ∇ ln ns⎢1 + ηs ⎜⎝es B ⎣
2T s
2
⎞⎤
⎡
⎛
2
cT s ⎢⎢⎢⎢ 1
1 ⎜⎜⎜⎜ ms  +2μB 3 ⎟⎟⎟⎟⎥⎥⎥⎥=−+− ⎟⎠⎥⎦ ky ,(45)⎢⎜
es cB ⎣ Lns LTs ⎝2T s2

3. Local Gyrokinetic Vlasov Simulation Code GKV

where ηs = Lns /LTs with Lns = −(d ln ns /dx)−1 and LTs =−(d ln T s /dx)−1 . The geometric coeﬃcients K x and Ky aredefined as follows:

In this section, we present the simulation model implemented in GKV, as well as the verification and benchmark test with realistic geometry from the MHD equilibrium solver.

3.1

g xz g xy − g xx gyz ∂ ln B ∂ ln B
−
,
∂z
∂y
B2 /c2Bg xz gyy − g xy gyz ∂ ln B ∂ ln B
+
.
Ky =
∂z
∂x
B2 /c2B

Kx =Calculation modelCalculation model in GKV is briefly summarized
below including electromagnetic eﬀects. One can alsofind the implementation of VMEC equilibria for 3dimensional helical/stellalator plasmas in Ref. [27]. The
wavenumber-space representation of electromagnetic gyrokinetic Vlasov equation for the particle species of s isgiven byc  ∂+  b·∇ + iωDs δhsk⊥ −b· k⊥ × k⊥ δψsk⊥ δhsk⊥
∂t
B Δ




es FMs ∂+ iω∗T s δψsk⊥ + Cs δhsk⊥ ,
=
(43)
Ts
∂t
where δhsk⊥ = δhsk⊥ (z,  , μ) stands for the non-adiabatic(g)part of the gyrocenter perturbed distribution function δ fsk⊥ ,(g)i.e., δhsk⊥ = δ fsk⊥ + es J0 (k⊥ ⊥ /Ωs )δφk⊥FMs /T s , and ms ,es , T s , and Ωs = es B/ms c are the particle mass, the electric charge, the equilibrium temperature, and the gyrofrequency, respectively. The parallel velocity  and the magnetic moment μ are used as the velocity-space coordinates,
where μ is defined by μ ≡ ms 2⊥ /2B with the perpendicular velocity ⊥ . The potential fluctuations averaged overthe gyrophase are denoted by δψsk⊥ = J0 (k⊥ ⊥ /Ωs )[δφk⊥ −( /c)δA k⊥ ], where J0 is the zeroth-order Bessel function,and the former and latter terms mean the electrostatic andelectromagnetic parts. Since we often focus on the finite
but low-β plasmas, the parallel magnetic field fluctuation
δB k⊥ is ignored here. The equilibrium part of the distribution function is given by the local Maxwellian distribution, i.e., FMs = ns (ms /2πT s )3/2 exp[−(ms 2 + 2μB)/2T s ],where ns represents the equilibrium density. The symbol
Δ appearing in the nonlinear term of Eq. (43) means double summations with respect to k⊥ and k⊥ , which satisfy

(46)
(47)

The electromagnetic potential fluctuations are determined by the Poisson-Ampère equations:
 


2
−2
k⊥ + λD δφk⊥ = 4πes dJ0s (k⊥ ⊥ /Ωs ) δhsk⊥ , (48)s4πk⊥2 δA k⊥ =c es d  J0s (k⊥ ⊥ /Ωs ) δhsk⊥ ,

(49)

s


where λD = ( s 4πns e2s /T s )−1/2 is the Debye length. Note
1, thethat, in the adiabatic electron limit with k⊥ ρte
gyrocenter perturbed distribution for electrons is approxi
(g)mated by d δ fek⊥  −ene (δφk⊥ − δφk⊥ z δky ,0 )/T e , whereρts = ms cts /|es |Bax means the gyroradius evaluated withaverthe thermal speed ts = (T s/ms )1/2 , and the field-line√√
age is defined by Xk⊥ z = dz g xyz Xk⊥ / dz g xyz .
By using the above gyrokinetic Vlasov and Poisson Ampère equations, one can derive another important equation describing the balance and transfer of the entropy(g)variable defined as δS sk⊥ ≡  d |δ fsk⊥ |2 /2FMs z . The entropy balance equation [e.g., see Ref. [28]] provides us
with many physical insights closely associated with the
turbulence saturation, and is formulated with the relation
among the development of fine-scale fluctuations on the
distribution function, the turbulent transport, and the collisional dissipation, as well as the linear and nonlinear entropy transfer dynamics. The detailed numerical analysesof the entropy balance and the transfer processes in ITG
turbulence are shown in, e.g., Ref. [19–21, 29].

3.2 Code verification and benchmark test
The implementation of realistic MHD equilibria on
GKV is verified by the benchmark test with other gyrokinetic codes, GENE [30], GS2 [31], and GKW [32]. In

1403029-5

Plasma and Fusion Research: Regular Articles

Volume 9, 1403029 (2014)

Fig. 3 Cross-code benchmark result on the ITG-mode dispersion relations with the analytic s − α equilibrium andCyclone-base-case like MHD equilibrium produced by
CHEASE.

this benchmark test, the same geometry data for Cyclonebase-case like MHD equilibrium from an MHD solver
CHEASE [33] , which is used in Ref. [34], is also applied
to GKV. The fluxtube geometry on the axisymmetric coordinates with ρ = (Ψ/Ψedge )1/2 is considered here. The
growth rate and the real frequency of electrostatic ITG
modes with the adiabatic electrons are compared in Fig. 3,
where the detailed physical parameters are shown in Ref.
[34]. In addition to the realistic MHD equilibrium, the results with analytic s−α models are also compared. The implementation of the realistic MHD geometry to GKV and
IGS has successfully been verified with good agreement in
the dispersion relation of the ITG-mode.

4. Linear Gyrokinetic Analyses with
JT-60SA Tokamak Equilibria
The extended GKV with IGS is applied to linear stability and residual zonal-flow analyses for the typical plasmas expected in JT-60SA tokamak device.

4.1

ITER-like and highly-shaped plasmas

In this study, two typical plasmas expected in JT60SA are considered. One is an ITER-like plasma withsingle-null separatrix [Figs. 4 (a) and 4(b)], and the otheris the highly-shaped plasma with quasi-double-null separatrix [Figs. 4 (c) and 4(d)], where the basic equilibrium
parameters are summarized in Table 1. The highly-shaped
plasma shows higher plasma current Ip , normalized-β βN ,elongation κ, and triangularity δ. Besides, the global shaping factor S = q95 Ip /Bax a is 1.45 times larger than that in
the ITER-like plasma, where q95 = q(ρ = 0.95). Note alsothat the highly-shaped plasma is almost up-down symmetupupdown, δ95  δdown
ric (κ95  κ9595 ), while the ITER-like plasma
up
up
down
, δ95 < δdown
shows the up-down asymmetry (κ95 < κ95
95 ).
The q- and ŝ-profiles for both cases are plotted in Fig. 5.
The ideal MHD stabilities have been well examined inFig. 4 Comparisons of the original 2D data Ψ (R, Z) for (a)ITER-like and (b) highly-shaped equilibria produced by
MEUDAS, and of [(b), (d)] constructed Boozer coordinates by IGS, where Zax means the vertical elevation of the magnetic axis. The flux surfaces for ρ ={0.30, 0.50, 0.75} are shown by red thick lines in (b) and(d).

Ref. [23], where the ballooning stability diagrams for boththe above ITER-like and highly-shaped plasmas have been
identified.
In the following, the electrostatic ITG-mode stability and the residual zonal-flow levels are examined, where
gyrokinetic ions and adiabatic electrons are considered.
The fluxtube geometry on Boozer coordinates with ρ =(Φ/Φedge )1/2 is used, and we look into three diﬀerent radial
positions, i.e., ρ = {0.30, 0.50, 0.75}, which correspondto diﬀerent flux-surface shapes. In order to focus on theimpact of local shaping eﬀects on the ITG-mode stability,
the normalized density and temperature gradient parame-

1403029-6

Plasma and Fusion Research: Regular Articles

Volume 9, 1403029 (2014)
Table 1 Equilibrium parameters.

ITER-likehighly-shaped

Rax [m]
3.10
3.18

a[m]
1.32
1.61

Bax [T]
2.51
2.61

Ip [MA]
2.49
5.00

q95
4.32
3.97

βN
2.02
3.11

up

κ95
1.63
1.82

down
κ95
1.83
1.81

up

δ95
0.29
0.46

δdown
95
0.44
0.43

S
3.25
4.72

Fig. 6 Comparison of structures along the magnetic field lines at diﬀerent ρ for [(a),(b)]the magnetic field intensity B, [(c),(d)] the drift
frequency ωDi (k x ρti = 0, ky ρti = 0.5,  /ti = 2.5, μBax /T i = 0), and [(e),(f)]the squared perpendicular wavenumber k⊥2 (k x ρti = 0, ky ρti =0.5), where the circular cases with same local values of q, ŝ and are also plotted by dashed lines. IT: ITER-like, HS: highlyshaped.

4.2 Structures along confinement magnetic
field

Fig. 5 Radial profiles of the safety factor q and the magnetic
shear ŝ for the ITER-like (IT) and highly-shaped (HS)cases.

ters are assumed to be constant in the radial region of interest, where Rax /Lni = 2.0 and Rax /LTi = 6.0 with T i = T e .
These are not the exact values obtained from equilibrium
calculations, but are still feasible. Note also that these parameters appear only in ω∗Ti [see Eq. (45)] which does notdepend on geometrical quantities.

Comparisons of the confinement field and linear advection operators, which are closely associated with the
ITG-mode stability, are first presented in this section. In
Fig. 6, the magnetic field intensity B, the drift frequency
ωDs , and the squared perpendicular wavenumber k⊥2 along
the magnetic field line are compared between the ITERlike and highly-shaped plasmas, where the circular caseswith same q, ŝ, and = ρ a shown in Appendix B are
also plotted. One finds that the deviation of B from thecircular case becomes more significant for the outer sideof the plasma [Figs. 6 (a) and 6(b)], where the magnetic
field strength increases (decreases) in the low (high) fieldside of z ∼ 0 (z ∼ π). These features are also seen in
the drift frequency ωDs for ρ = 0.30 and ρ = 0.50, while
the qualitatively diﬀerent structures with the multiple setsof nodes appear for ρ = 0.75, and the deviation from thecircular case is more pronounced for the highly-shapedcase [Figs. 6 (c) and 6(d)]. Similarly, in the squared perpendicular wavenumber k⊥2 , the deviation becomes drastic
for both the ITER-like and highly-shaped cases [Figs. 6 (d)and 6(e)], and then the slightly stronger up-down asymmetry is seen in the ITER-like case, e.g., the asymmetry at
ρ = 0.5 and ρ = 0.75 is more apparent.
The comparison of the trapped-passing boundaries at

1403029-7

Plasma and Fusion Research: Regular Articles

Volume 9, 1403029 (2014)

Fig. 7 Trapped-passing boundaries for diﬀerent μ-coordinatesin (a) ITER-like and (b) highly-shaped equilibria, where
the circular results are also plotted by dashed line. The
highly-shaped case shows slightly wider region in comparison with the ITER-like case.

ρ = 0.75 is shown in Fig. 7, where the circular cases arealso plotted for the ITER-like plasma. As compared withthe circular cases with the sinusoidal boundary shapes, thenarrower trapped region with sharper boundaries is found.
Also, the trapped region becomes slightly wider for thehighly-shaped case, and these features indicate an impact
on the trapped-electron mode (TEM) stability.

4.3

Fig. 8 Comparisons of linear ITG-mode dispersion relations
among (a) ITER-like (IT) and (b) highly-shaped (HS)cases, where the circular results are shown by dashedlines. The highly-shaped case shows lower growth rate
in comparison with the ITER-like case.

Comparison of ITG-mode stability

In this section, we make detailed comparisons on the
ITG-mode stability between the ITER-like and highlyshaped equilibria. The linear growth rate γITG and the real
frequency ωr are shown in Fig. 8, where the corresponding circular cases are also plotted for comparisons. In the
ITER-like case [Fig. 8 (a)], we see that, for ρ = 0.30 andρ = 0.50, the maximum growth rate is similar to that in thecircular case, while the profiles are shifted to the higher-ky
region. The up-shift of γITG -profile is more significant forρ = 0.75, where the maximum growth rate also increases ascompared to the circular case. The real frequency ωr also
shows the up-shift feature with increasing ρ. The similar
qualitative features on the spectral up-shift of γITG and ωr
are also identified for the highly-shaped case [Fig. 8 (b)],whereas the magnitude of the growth rate is lower thanthat in the ITER-like case for all three radial positions.
Since the low-wavenumber modes around ky ρti ∼ 0.5 are,in general, considered to make dominant contributions to
the ITG-driven turbulent transport, the stabilization in thelower-ky region are expected to reduce the transport levels.
A qualitative explanation on the physical mechanism
leading to relatively lower growth rate in the highly-shapedcase is given as follows. By considering a fluid limit
ignoring the wave-particle resonance and finite gyroradius eﬀects, the ITG-mode growth rate is evaluated byγITG  [ω∗Pi ωDi /(1 + τe k⊥2 ρ2ti )]1/2 for ω∗Pi ωDi > 0, whereω∗Pi = −ky (1 + ηi )cT i /ei cB Lni  0, τe = T e /T i , and ωDi (z)and k⊥2 (z) depend on the metric components [see. Eqs. (41)and (44)]. The reduction of |ωDi | and the increase of k⊥2
due to the shaping eﬀect are, thus, responsible for the decrease of γITG . Actually, for ωDi < 0, the magnitude of|ωDi |/(1 + τe k⊥2 ρ2ti ) in the highly-shaped case is lower thanthat in the ITER-like one for all three radial positions.
In Fig. 9, linear eigenmode structures along the fieldlines are compared, where the eigenmodes giving the maximum growth rate are plotted. Reflecting the diﬀerent
structures of B, ωDs , and k⊥2 shown in Fig. 6, the deviationof the eigenmode profile from the circular case becomeslarger as ρ increases. One also finds that the highly-shapedcase indicates more strongly localized peak near θ ∼ 0, as
compared with the ITER-like (and also circular) case.

4.4 Comparison of linear zonal-flow dynamics
Finally, we present the comparisons of linear zonalflow evolutions in this section. Figure 10 shows time evolutions in the linear zonal-flow damping and the radialwavenumber (k x ρti ) dependence of the residual zonal-flow
levels for ρ = 0.30, ρ = 0.50, and ρ = 0.75, where the circular cases are also plotted for the ITER-like cases. We seethat the geodesic acoustic mode (GAM) frequency does

1403029-8

Plasma and Fusion Research: Regular Articles

Volume 9, 1403029 (2014)

Fig. 9 Comparison of eigenmode structures giving the maximum growth rate for [(a), (c), (e)] ITER-like (IT) and [(b), (d), (f)] highlyshaped (HS) cases, where “Re” and “Im” represent the real and imaginary parts, respectively, and the circular cases are also shownby dashed lines.

Fig. 10 [(a), (b), (c)] time evolutions in the linear zonal-flow damping for k x ρti = 0.1 at ρ = 0.30, ρ = 0.50, and ρ = 0.75, and [(d), (e), (f)]the radial wavenumber dependencies of the residual zonal-flow level evaluated at t = 70, where “R-H” shows the analytic circular
result of (1 + 1.6q2 / 1/2 )−1 [15]. The higher residual zonal-flow level is confirmed in the highly-shaped case for all three radial
positions.

not change much in all the cases, while the residual zonalflow levels are enhanced as the shaping eﬀects become
more significant [Figs. 10 (a) – 10(c)], i.e., “circular case”< “ITER-like case” < “highly-shaped case”, and “innerside (ρ = 0.30)” < “mid-radius (ρ = 0.50)”< “outer side(ρ = 0.75)”. Particularly, for k x ρti = 0.1 which is a typical wavenumber for nonlinearly generated zonal flows in
the ITG turbulence, approximately 2 times higher residuallevel is identified in the highly-shaped case.
Another important finding is diﬀerence in the radial
wavenumber dependence of the residual level [Figs. 10 (d)– 10(f)]. As derived in Ref. [15], the k x -dependence of theresidual zonal-flow level in circular tokamaks appears onlyin the higher order in . Indeed, such weak k x -dependence
is confirmed in the present numerical simulations with thecircular equilibrium. In contrast, the strong k x -dependenceof the residual levels is revealed for both the ITER-likeand highly-shaped cases, where almost linear dependencies with −0.08k x for ρ = 0.30 and ρ = 0.50, and with
−0.16k x for ρ = 0.75 are identified for k x ρti > 0.5.

5. Summary
Realistic tokamak geometries have been successfully implemented to a local fluxtube gyrokinetic Vlasovcode GKV, by using a newly developed interface code1403029-9

Plasma and Fusion Research: Regular Articles

Volume 9, 1403029 (2014)

IGS which can construct straight-field-line flux coordinates systems from MHD equilibria produced by a freeboundary 2D Grad-Shafranov equation solver MEUDAS.
The accuracy of the flux coordinates, i.e., the axisymmetric, Boozer, and Hamada, constructed by IGS has been verified with the analytic Solov’ev equilibrium model. As another verification of GKV with IGS, cross-code benchmark
test has been carried out by using Cyclone-base-case like
MHD equilibrium, and good agreement in the dispersion
relation of ITG-mode are confirmed.
By using the extended GKV, the ITG-mode stability
and linear zonal-flow dynamics are investigated for two
types of shaped plasmas expected in JT-60SA tokamak device, i.e., ITER-like and highly-shaped plasmas in which
the MHD stability has been well examined. The detailed
comparisons including the conventional circular equilibrium reveal that both the ITER-like and highly-shaped
cases show the up-shift of γITG - and ωr -profiles to thehigher-ky region in comparison with the circular cases. The
maximum growth rate is similar to or less than the circularcase for the inner sides of plasma (ρ = 0.30 and ρ = 0.50),
while the outer side with ρ = 0.75 show the higher maximum growth rate at higher ky . It is also found that, whenthe ITER-like and highly-shaped cases are compared, thelatter case shows relatively lower growth rate for all the radial positions examined here, where the maximum growing
eigenmode structure around θ ∼ 0 becomes sharper. Linear zonal-flow damping simulations clarify that the residual level is significantly enhanced by the shaping eﬀects,
and the radial wavenumber dependence becomes stronger
in both the ITER-like and the highly-shaped cases.
From these analyses, it is expected that the highlyshaped equilibrium shows less ITG-driven turbulent transport than that in the ITER-like one, due to the reduction of
the ITG-mode growth rate and the enhancement of residual zonal-flow levels. It should be noted that since theactual zonal-flow generation in ITG turbulence also depends on the turbulence intensity, which is well correlated
to γITG /k⊥2 [35, 36], the reduction of the ITG-mode growthrate can also decrease the source of zonal flows. However,
the enhancement ratio of the residual zonal-flow level fromthe ITER-like to highly-shaped cases exceeds the reduction
ratio of the ITG-mode growth rate. The detailed comparisons on the nonlinear zonal-flow generations and turbulent
transport levels are demanded in a future work, and will bereported elsewhere.

Acknowledgments

The Authors would like to thank Dr. S. Brunner, Dr.
J. Dominski, and Dr. X. Lapillonne for providing the
CHEASE equilibrium data and benchmark results, andalso thank Dr. Hideo Sugama, Dr. M. Honda, and Dr. M.
Yoshida for fruitful discussions.
Numerical simulations were performed by HELIOSat Computational Simulation Centre of International Fu-

sion Energy Research Centre, Oakleaf-FX at the University of Tokyo, and BX900 at Japan Atomic Energy Agency.
This work is supported in part by the HPCI Strategic Program for Innovative Research, and in part by the NIFS Collaborative Research Program Nos. NIFS12KNTT014 and
NIFS12KNTT015, and in part by the MEXT, Grant Nos.
22760660 and 23246163.

Appendix A. Hamada, Boozer, andAxisymmetric Coordinates
The fundamental properties of each coordinate system
are briefly summarized here.

A.1 Hamada coordinates (α̃ = 0, f = Const.)V√gH (ρ) =
,
4π2

Ψ
,
V
Φ
B ζH (ρ) = B · ∇ ζH = 4π2  .
V

(A.1)

B θH (ρ) = B · ∇θH = 4π2(A.2)

From Eqs. (1), (2), (4), (11), and (12), one finds  B2B · ∇γ̃H = B2  2 − 1 ,BB · ∇β̃H = 0 ⇒ β̃H = β̃H (ρ) .

(A.3)
(A.4)

The equation (A.4) is the necessary and suﬃcient condition for Hamada coordinates, but one can always set h (ρ)ρ
in Eq. (6) such that β̃∗H = 0, i.e., h(ρ) = − 0 dρ̄ β̃H (ρ̄). Besides, no h(ρ) exists to vanish γ̃∗H = γ̃H (ρ, θH , ζH ) − h(ρ)
at any points in (θH , ζH ) for a given ρ. Then, the covariant
components of B are given as
∂γ̃∗H
,
∂ρ
∂γ̃∗
BθH (ρ, θH , ζH ) = I + H ,∂θH∂γ̃∗BζH (ρ, θH , ζH ) = J + H ,∂ζHBρ (ρ, θH , ζH ) =

(A.5)
(A.6)
(A.7)

where γ̃∗H is determined by Eq. (A.3) (Note that B · ∇γ̃H =B·∇γ̃∗H ). Also, the contravariant components of j are givenasj θH (ρ) = J  , j ζH (ρ) = I  .

(A.8)

In Hamada coordinates, both the magnetic field andcurrent density lines become straight on the flux surface.
Although the Jacobian depending only on ρ is significantly
simplified (especially, for the flux-surface average), the covariant components of B become rather complicated.

A.2 Boozer coordinates (α̃ = 0, f = B2 )1403029-10

 
2
V B
√
gB (ρ, θB , ζB ) = 2 2 ,
4π B

(A.9)

Plasma and Fusion Research: Regular Articles
Ψ  B2,
V  B2
Φ B2
B ζB (ρ, θB , ζB ) = B · ∇ ζB = 4π2   2 .
V BVolume 9, 1403029 (2014)

B θB (ρ, θB , ζB ) = B · ∇θB = 4π2

BθA (ρ, θA , φ) = I +(A.10)

From Eqs. (1), (2), (4), (11), and (12), one finds
B2
B · ∇β̃B = p 1 −  2 ,
BB · ∇γ̃B = 0 ⇒ γ̃B = γ̃B (ρ) .

(A.12)Bρ (ρ, θB , ζB ) = β̃∗B ,(A.13)BθB (ρ) = I ,(A.14)BζB (ρ) = J ,(A.15)

where β̃∗H is determined by Eq. (A.11) (Note that B · ∇β̃H =B · ∇β̃∗H ). The contravariant components of j are given as
∗

∗

∂β̃
∂β̃
, j ζB (ρ, θB , ζB ) = I  +
.
∂ζB
∂θB
(A.16)In Boozer coordinates, only the magnetic field line becomes straight on the flux surface, but the covariant components of B are simplified. This feature is advantageous
to solve numerically the guiding-center equations of motion.

A.3

Axisymmetric coordinates (α̃ = 0, f =B·∇ζA , ζA = φ)V   −2  2 Φ R2
,
R R =
J
4π2
J
B θA (ρ, θA , φ) = B · ∇θA =,qR2
JB φ (ρ, θA , φ) = B · ∇ φ = 2 .
R√

gA (ρ, θA , φ) =(A.17)(A.18)

From Eqs. (1), (2), (4), (11), and (12), one findsR−2
B · ∇β̃A = p 1 −  −2 ,R  B2R−2
B · ∇γ̃A = B2  2 −  −2
BR

(A.19).

(A.20)

In contrast to Hamada and Boozer ones, the axisymmetric coordinates show the finite B · ∇β̃A and B · ∇γ̃A . Itis, however, noted that when the poloidal magnetic fieldstrength Bpol = Ψ  |∇ρ|R−1 is negligibly small, the axisymmetric coordinates coincide with Boozer one, i.e.,B · ∇γ̃A → 0 (B2 → B2tor = J 2 R−2 ). The covariant components of B are given as
Bρ (ρ, θA , φ) = β̃A +

qΨ 2
∂γ̃A
=−
∇ρ·∇θA ,(A.21)∂ρJ(A.22)(A.23)

where, ∇ρ·∇φ = 0, ∇θA ·∇φ = 0, and ∂/∂φ = 0 are used.
Finally, the contravariant components of j are given as
j θA (ρ) = J  , j φ (ρ, θA , φ) = I  +(A.11)The Eq. (A.12) is the necessary and suﬃcient condition for
Boozer coordinates, but one can always set h(ρ) in Eq. (7)such that γ̃∗B = 0. The covariant components of B are, then,given asj θB (ρ, θB , ζB ) = J  +Bφ (ρ) = J ,∂γ̃A qΨ 2∇ρ·∇ρ ,
=
∂θA
J

∂β̃
.
∂θA

(A.24)In the axisymmetric coordinates, only Bζ and jθ are simplified, i.e., Bζ = Bζ (ρ) and jθ = jθ (ρ), and this coordinates
are widely used for theoretical and numerical analyses oftokamak plasmas.

Appendix B. Metric Components for
Circular and s-α Models
Analytic expressions of the concentric circular equilibrium and the simplified s−α one are summarized here.
The detailed descriptions are shown in Ref. [34].

B.1 Concentric circular equilibrium
For toroidal plasmas with concentric circular flux surfaces, the axisymmetric coordinates (ρ, θA , φ) are useful,where ρ = (Φ/Φedge )1/2 = r/a with a = (2Φedge /Bax )1/2
is considered here. By using geometric polar-angle fromthe magnetic axis, θ0 , in (ρ, θ0 , φ), the SFL poloidal-angle
coordinate θA is written by⎡⎤! θ "⎥⎥⎢⎢ 1 −  12
0 ⎥
−1 ⎢
⎢
⎥⎥⎥ .
θA = 2 tan ⎢⎢⎣
tan
(B.1)
1+
2 ⎦
The magnetic flux surfaces for given q are determined bydΨ/dr = Bax r/q̄, where q̄ = (1 − 2 )1/2 q. Using thesedefinitions, one can calculate the magnetic field intensity,
its derivatives, and the contravariant metric components asfollows: (Note that q0 = q(ρ0 ) and q̄0 = q̄(ρ0 ).)
⎞1
⎛
2 ⎟2
⎜⎜⎜
Bax
⎟
,
(B.2)
B = ⎜⎝1 + 2 ⎟⎟⎠q̄0 1 + cos θ0#∂BB=
sin2 θ0cos θ0 −∂x
Bax Rax (1 + cos θ0 )1− 2⎤
(1 − ŝ + 2 q20 /q̄20 ) ⎥⎥⎥
⎥⎦ ,
−
(B.3)
q̄20 + 2 q̄0
B sin θ0
∂B
=
,
(B.4)
2
∂z
q0 1 − cos θ0 cos θA√g xyz = q0 Rax (1 + cos θ0 )2 ,
(B.5)

2
2 
q
q
q̄0 sin θA,(B.6)
g xx = 02 , g xy = 02 ŝθA −q0 1 − 2q̄0q̄0q0 sin θag xz = −,(B.7)
q̄0 1 − 2q2
sin θA2q̄0gyy = 02 ŝ2 θ2A −
ŝθA
q
1
− 2
q̄0
0
⎡ 2
⎤
2
⎢⎢ q̄ (1 + 2 /q̄20 )sin2 θA ⎥⎥⎥⎥⎦ ,+ ⎢⎢⎣ 02+(B.8)
q0 (1 + cos θ0 )2 (1 − 2 )21403029-11

Plasma and Fusion Research: Regular Articles

Volume 9, 1403029 (2014)⎡sin θA q̄20
1
1 ⎢⎢⎢ q0⎢⎣− ŝθAg =+ 2
2
Rax q̄01−
q0 (1 + cos θ0 )2
yz+

sin2 θA,(1 − 2 )2⎡
1 ⎢⎢⎢ q̄20zz
g = 2 ⎣⎢ 2
Rax q0B.2

2⎤
1
sin2 θA ⎥⎥⎥+⎦⎥ .
2 (1 + cos θ )2
(1 − 2 )20(B.9)(B.10)

s-α equilibrium model

The so-called s−α equilibrium model used in gyrokinetic simulations is a simplified model of the above concentric circular one. By expanding Eqs. (B.1) – (B.10) withrespect to , but keeping the lowest or the first order terms,
one can obtain the expressions in s − α model as follows:(Note that q0 = q̄0 and θA = θ0 .)Bax cos θ0∂B=−, (B.11)∂x
Rax∂B
q0 Rax Bax√= Bax sin θ0 ,,(B.12)
g xyz =∂zBg xx = 1, g xy = ŝθ0 , g xz = 0 ,(B.13)
1
1
yy
2 2
yz
zz
, g = 2 2 . (B.14)
g = 1 + ŝ θ0 , g =Rax
Rax
B = Bax (1− cos θ0 ),It should be noted that, in contrast to the concentric circular model, the s−α model no longer holds the consistency
between the metric components and the magnetic field intensity given as B2 = c2B [g xx gyy − (g xy )2 ].
[1] W. Horton, Rev. Mod. Phys. 71, 735 (1999).
[2] X. Garbet, Y. Idomura et al., Nucl. Fusion 50, 043002(2010).
[3] A.M. Dimits, G. Bateman et al., Phys. Plasmas 7, 969(2000).
[4] P.H. Diamond, S.-I. Itoh et al., Plasma Phys. Control. Fusion 47, R35 (2005).
[5] M.A. Beer, S.C. Cowley et al., Phys. Plasmas 2, 2687
(1995).
[6] R.L. Miller, M.S. Chu et al., Phys. Plasmas 5, 973 (1998).
[7] R.E. Waltz and R.L. Miller, Phys. Plasmas 6, 4265 (1999).
[8] E.A. Belli, G.W. Hammett et al., Phys. Plasmas 15, 092303(2008).
[9] E.A. Belli and J. Candy, Phys. Plasmas 17, 112314 (2010).
[10] Y. Camenen, A. Pochelon et al., Nucl. Fusion 47, 510(2007).
[11] A. Marinoni, S. Brunner et al., Plasma Phys. Control. Fusion 51, 055016 (2009).

[12] A. Burckel, O. Sauter et al., J. Phys. Conf. Ser. 260, 012006(2010).
[13] T.L. Rhodes, C. Holland et al., Nucl. Fusion 51, 063022(2011).
[14] R.V. Bravenec, J. Candy et al., Phys. Plasmas 18, 122505(2011).
[15] M.N. Rosenbluth and F.L. Hinton, Phys. Rev. Lett. 80, 724(1998).
[16] Y. Xiao and P.J. Catto, Phys. Plasmas 13, 082307 (2006).
[17] O. Yamagishi and H. Sugama, Phys. Plasmas 19, 092504(2012).
[18] H. Sugama, T.-H. Watanabe et al., Phys. Plasmas 14,022502 (2007).
[19] M. Nakata, T.-H. Watanabe et al., Phys. Plasmas 19,
022303 (2012).
[20] M. Nakata and Y. Idomura, Nucl. Fusion 53, 113039(2013).
[21] T.-H. Watanabe and H. Sugama, Nucl. Fusion 46, 24(2006).
[22] M. Azumi et al., in Proc. 4th Int. Symp. Comput. Methods
Applied Sci. Engineering, Paris (North-Holland, Amsterdam, 1980) p.335.
[23] N. Aiba, S. Tokuda et al., Plasmas Fusion Res. 2, 010(2007).
[24] N. Nakajima, J. Todoroki et al., Kakuyugo Kenkyu 68, 395(1992).
[25] V.D. Pustovitov, J. Plasma Fusion Res. SERIES 1, 456(1998).
[26] T. Takeda and S. Tokuda, J. Comput. Phys. 93, 1 (1991).
[27] M. Nunami, T.-H. Watanabe et al., Plasma Fusion Res. 5,016 (2010).
[28] H. Sugama, T.-H. Watanabe et al., Phys. Plasmas 16,112502 (2009).
[29] S. Maeyama, A. Ishizawa et al., Comput. Phys. Comm.
184, 2462 (2013).
[30] F. Jenko and W. Dorland, Plasma Phys. Control. Fusion 43,
A141 (2001).
[31] W. Dorland, F. Jenko et al., Phys. Rev. Lett. 85, 5579
(2000).
[32] A.G. Peeters, Y. Camenen et al., Comput. Phys. Comm.
180, 2650 (2009).
[33] H. Lütjens, A. Bondeson et al., Comput. Phys. Comm. 97,
219 (1996).
[34] X. Lapillonne, S. Brunner et al., Phys. Plasmas 16, 032308(2009).
[35] M. Nunami, T.-H. Watanabe et al., Phys. Plasmas 19,
042504 (2012).
[36] M. Nunami, T.-H. Watanabe et al., Phys. Plasmas 20,
092307 (2013).

1403029-12