A new hybrid kinetic electron model for full-f gyrokinetic simulations

Journal of Computational Physics 313 (2016) 511–531

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

A new hybrid kinetic electron model for full- f gyrokinetic simulations Y. Idomura Japan Atomic Energy Agency, Wakashiba 178-4, Kashiwa, Chiba 277-0871, Japan

a r t i c l e

i n f o

Article history: Received 28 September 2015 Received in revised form 8 January 2016 Accepted 23 February 2016 Available online 27 February 2016 Keywords: Gyrokinetics Full- f model Kinetic electrons Ion temperature gradient mode Trapped electron mode Neoclassical transport

a b s t r a c t A new hybrid kinetic electron model is developed for electrostatic full- f gyrokinetic simulations of the ion temperature gradient driven trapped electron mode (ITG-TEM) turbulence at the ion scale. In the model, a full kinetic electron model is applied to the full- f gyrokinetic equation, the multi-species linear Fokker–Planck collision operator, and an axisymmetric part of the gyrokinetic Poisson equation, while in a non-axisymmetric part of the gyrokinetic Poisson equation, turbulent fluctuations are determined only by kinetic trapped electrons responses. By using this approach, the so-called ω H mode is avoided with keeping important physics such as the ITG-TEM, the neoclassical transport, the ambipolar condition, and particle trapping and detrapping processes. The model enables full- f gyrokinetic simulations of ITG-TEM turbulence with a reasonable computational cost. Comparisons between flux driven ITG turbulence simulations with kinetic and adiabatic electrons are presented. Although the similar ion temperature gradients with nonlinear upshift from linear critical gradients are sustained in quasi-steady states, parallel flows and radial electric fields are qualitatively different with kinetic electrons. © 2016 Elsevier Inc. All rights reserved.

1. Introduction Turbulent fusion plasmas contain multi-scale phenomena covering wide spatio-temporal scales from the fast Cyclotron motion to slow profile evolutions. The issue of the scale gap between the fast Cyclotron motion and micro-turbulence was successfully resolved by the development of the gyrokinetic theory [1] and gyrokinetic simulations [2–4], which are standard approaches for analyzing turbulent transport in fusion plasmas. On the other hand, another scale gap between micro-turbulence and slow profile evolutions has been treated by conventional δ f approaches, where the total particle distribution f is separated into a macroscopic collisionless or collisional equilibrium distribution f 0 and a microscopic turbulent perturbation δ f , and only the latter is evolved by assuming a complete scale separation. However, the validity of this scale separation is not so clear, when an interaction between turbulence and plasma profiles plays an essential role as in avalanche like non-local transport phenomena. In order to resolve this fundamental issue, in recent years, the so-called full- f or total- f approaches have been developed. This approach evolves the total particle distribution f following turbulent and collisional transport processes described by the same first principle. Because of this multi-scale and multi-physics properties, full- f simulations [5–7] disclosed rich physics such as self-organized critical phenomena [8] in avalanche-like non-local transport, turbulence regulation by mean flows or radial electric fields given by the neoclassical radial force balance, and intrinsic rotation induced by non-diffusive momentum transport. With increasing computing power, longer time

E-mail address: [email protected]. http://dx.doi.org/10.1016/j.jcp.2016.02.057 0021-9991/© 2016 Elsevier Inc. All rights reserved.

512

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

scale simulations and systematic parameter scans have been enabled with full- f approaches. In Ref. [9], a confinement time scale full- f gyrokinetic simulation was presented, and formations of temperature profiles and intrinsic rotation profiles were discussed based on the power balance and the toroidal angular momentum conservation. Validation studies on transport scalings with respect to the plasma size [10–13], the collisionality [12], and the heating power [13] were performed, and qualitative features of experimental transport scalings were recovered. In addition to core plasma turbulence, full- f gyrokinetic simulations have been applied to edge plasma turbulence [14] and to energetic particles driven geodesic acoustic modes (EGAMs) [15,16]. Although applications of full- f gyrokinetic simulations have been greatly expanded, their electron models have been limited to adiabatic electrons, and electron transport in full- f gyrokinetic simulations has been an open issue for a long time. In order to resolve this critical issue, in this work, we develop a new kinetic electron model for full- f gyrokinetic simulations. In δ f flux-tube simulations, an electromagnetic model with full kinetic electrons [17] has been established as a standard model. However, an electromagnetic full- f model is not available yet because of the following difficulties. One is the so-called cancellation problem, where a skin term in the Ampere’s law has to be exactly cancelled with an adiabatic part of the guiding-center current. Although this issue was resolved for δ f approaches, which is consistent with a linear approximation of the skin term, one may need to compute more accurate pull-back transform in full- f approaches. The other issue is a gyrokinetic MHD equilibrium. Since the total particle distribution f involves equilibrium quantities, one may need to consider a gyrokinetic MHD equilibrium, so that it is consistent with various kinetic effects such as trapped particles, finite orbit widths, finite Larmor radii (FLR), and a bootstrap current. Before starting electromagnetic full- f gyrokinetic simulations, we need further investigations on these issues, and therefore, in this work, we limit ourselves to electrostatic models. It is well known√that if one uses a full kinetic electron model, the so-called ω H mode [2] or a high frequency mode with ω H = (k /k⊥ ) mi /me i appears in the electrostatic limit of the kinetic Alfven wave. Since ω H becomes comparable to the ion Cyclotron frequency i at long wavelength, the ω H mode is unphysical from the viewpoint of the gyrokinetic ordering ω/i  1 [1]. This unphysical mode is unfavorable also from the numerical viewpoint, because it limits the time step width extremely short. In order to avoid the ω H mode, a kinetic trapped electron model [18–20], a bounce-averaged trapped electron model [19], and a fluid-kinetic hybrid electron model [21] are often used. Although these models eliminate the ω H mode from electrostatic ion temperature gradient driven trapped electron mode (ITG-TEM) turbulence simulations by assuming adiabatic passing electrons, the following drawbacks arise. Trapping and detrapping processes due to electron transport, turbulent pitch angle scattering, and collisional pitch angle scattering are not described, and one may need artificial boundary conditions at trapped-passing boundaries. The neoclassical transport and the ambipolar condition are not correctly described without passing electrons, and thus, one cannot determine a relevant kinetic equilibrium in full- f approaches. The ambipolar condition is important also for the toroidal angular momentum conservation. The above contradicting issues are resolved by a new hybrid kinetic electron model. In the model, we apply a full kinetic electron model to the full- f gyrokinetic equation, the multi-species linear Fokker–Planck collision operator, and an axisymmetric part of the gyrokinetic Poisson equation, while in a non-axisymmetric part of the gyrokinetic Poisson equation, turbulent fluctuations are computed using kinetic trapped electrons responses. By taking this approach, one can avoid the ω H mode with keeping important physics such as the ITG-TEM, the neoclassical transport, the ambipolar condition, and particle trapping and detrapping processes. The model is verified through test calculations on collisional temperature relaxation, neoclassical transport, zonal flow damping, and linear stability of ITG-TEMs. After these verification tests, the model is applied to a full- f gyrokinetic simulation, and influences of kinetic electrons are discussed from comparisons of flux driven ITG turbulence simulations with kinetic electrons and with adiabatic electrons. The reminder of the paper is organized as follows. In Sec. 2, the multi-species full- f gyrokinetic equations with the new hybrid kinetic electron model is presented. In Sec. 3, a series of verification tests on the hybrid kinetic electron model are shown. In Sec. 4, flux driven ITG turbulence simulations with kinetic electrons and with adiabatic electrons are compared, and influences of kinetic electrons are discussed. In this section, the numerical convergence is also shown. Finally, a summary is given in Sec. 5. 2. Calculation model 2.1. Gyrokinetic equation We consider the electrostatic ITG-TEM turbulence in an axisymmetric tokamak configuration. GT5D is based on the modern gyrokinetic theory [1], in which the gyrokinetic equation is simply given as the Liouville equation,

D fs Dt

∂ fs + { f s, Hs} ∂t ∂ fs ˙ ∂ fs ∂ fs = +R· + v˙ ∂t ∂R ∂ v = 0, ≡

using the gyro-center Hamiltonian,

(1)

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

Hs =

1 2

513

ms v 2 + μ B + q s φα ,

(2)

and the gyrokinetic Poisson bracket operator,

{F , G} ≡

s B



∂ F ∂G ∂ F ∂G − ∂α ∂μ ∂μ ∂α



+

  ∂F ∂G c − · ∇ F − ∇ G b · ∇ F × ∇G, ∗ ms B  ∂ v ∂ v q s B ∗ B∗

(3)

in the gyro-center coordinates Z = (t ; R, v  , μ, α ). Here, f s denotes the guiding-center distribution function, R is the position of the guiding center, v is the velocity of the guiding center, v  = b · v and v ⊥ = |b × v| are the velocities in the parallel and perpendicular direction to the magnetic field, μ = ms v 2⊥ /2B is the magnetic moment, α is the gyro-phase angle, B = Bb is the magnetic field, b is the unit vector in the parallel direction, ms and q s are the mass and charge of the particle species s, respectively, c is the velocity of light, s = q s B /ms c is the cyclotron frequency, B ∗ = b·B∗ is a parallel component of B∗ = B + ( B v  /s )∇ × b, φ is the electrostatic potential of turbulent fields, and the gyro-averaging operator is defined as ·α ≡ ·dα /2π . The nonlinear characteristics Z˙ = {Z, H s } are given as

B∗

∂H c + b×∇ H ms B ∗ ∂ v  e s B ∗   c = v b + b× q s ∇φα + ms v 2 b·∇ b + μ∇ B , ∗ qs B 

˙ ≡ {R, H s } = R

v˙ ≡ { v  , H s } = −

=−

B∗

ms B ∗ B∗

ms B ∗

(4)

· ∇H · (q s ∇φα + μ∇ B ) .

(5)

By using the phase space volume conservation,

∇ · (Js R˙ ) + ∂ v  (Js v˙  ) = 0,

(6)

the gyrokinetic equation (1) yields its conservative form,

∂ ∂ Js f s + ∇ · (Js R˙ f s ) + (Js v˙  f s ) = 0, ∂t ∂ v

(7)

where Js = m2s B ∗ is the Jacobian of the gyro-center coordinates. By adding a multi-species Coulomb collision operator C ( f s , f s ) and a source term S src ,s , a conservative gyrokinetic equation used in GT5D is written as,

 ∂ ∂ Js f s + ∇ · (Js R˙ f s ) + (Js v˙  f s ) = Js C ( f s , f s ) + Js S src ,s . ∂t ∂ v 

(8)

s

The multi-species collision operator is given by the zeroth order equipartition operator and the linear Fokker–Planck collision operator [22]. The former is essential in full- f approaches, in which an equilibrium distribution or a local Maxwellian distribution of each species may have different temperature at each time step, and the resulting interaction among multiple transport channels through the equipartition process affects relative importance of each transport channel and the total power balance. Details of the multi-species collision operator are summarized in Appendix A. 2.2. Gyrokinetic Poisson equation



The self-consistency is imposed by the gyrokinetic Poisson equation or the quasineutrality condition, s q s δn s = 0, where  the perturbed density δns is defined with respect to the initial density ns0 satisfying q n = 0. In this work, we define s s0 s the ion density using a linear polarization density with a long wavelength approximation, k2⊥ ρti2  1,



δni =

δ f i δ([R + ρ ] − x)d6 Z +

1 4π q i

∇⊥ ·

ρti2 λ2Di

∇⊥ φ,

(9)

where δ f s is the perturbed distribution defined against the initial distribution f s0 , R + ρ is the particle position, ρ = b × v/s is the Larmor radius, d6 Z = m2s B ∗ dRdv  dμdα is the phase space volume of the gyro-center coordinates,

ρts = v ts /s is the Larmor radius evaluated with the thermal velocity v ts = ( T s /ms )1/2 , λ Ds = ( T s /4π ns q2s )1/2 is the Debye length. Although linear growth rate spectra of the TEM expand up to a short wavelength region, k⊥ ρti > 1, and are often connected with the electron temperature gradient driven mode (ETG), multi-scale ITG-TEM-ETG turbulence simulations using the Cyclone parameters [23] have shown that the wavenumber spectra of turbulent transport fluxes peak at long wavelength, k⊥ ρti < 0.3 [24]. Therefore, the long wavelength approximation may be valid at least for ITG dominant parameters such as the electrostatic Cyclone case.

514

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

In calculating the electrostatic ITG-TEM turbulence, it is important to eliminate the ω H mode. For this purpose, we compute turbulent fluctuations using kinetic trapped electrons and an adiabatic passing electrons response. The lowest cost kinetic electron model for ITG-TEM turbulence simulations may be given by a bounce-averaged drift-kinetic model [19]. However, this model is formulated by assuming deeply trapped particles in the collisionless limit, and it is difficult to describe the neoclassical transport and de-trapping (trapping) processes of barely trapped (passing) electrons in nonlinear global gyrokinetic simulations. Another kinetic electron model is a kinetic trapped electron model, in which the drift-kinetic (or gyrokinetic) equation is solved only for trapped electrons in a straightforward manner by assuming adiabatic passing electrons. Although this model has been applied to linear [18,25,26] and nonlinear [20] global gyrokinetic simulations, the above issues have not been resolved yet. In this work, we propose a new hybrid kinetic electron model by taking account of the following effects. First, we solve the full- f gyrokinetic equation (8) using a full kinetic electron model including both trapped and passing electrons. This treatment enables us to compute collisional processes and thus the neoclassical transport without any artificial boundary condition at trapped-passing boundaries. In nonlinear simulations, this treatment can also describe trapping and detrapping processes due to the parallel nonlinearity, which changes the pitch angle, and the radial transport of electrons between magnetic surfaces with different local aspect ratios and thus trapped-passing boundaries. Second, we use only the perturbed density of trapped electrons for calculating the electrostatic potential of finite n modes, while the radial electric field E r is computed using the perturbed density of both trapped and passing electrons. This treatment is essential for eliminating the ω H mode with keeping the ambipolar condition. If one uses only the perturbed density of trapped electrons, the ambipolar condition is violated as shown in Sec. 3.3. Since the ambipolar condition is connected with the toroidal angular momentum conservation, this feature is of critical importance for computing toroidal rotation. Third, in the n = 0 mode or E r , we neglect m = 0 modes, which have finite k and are subject to the ω H mode. Here, m and n denote the poloidal and toroidal mode numbers  defined using a mode expansion with respect to the (geometry) poloidal angle θ and the toroidal angle ϕ , A (θ, ϕ ) = m,n A m,n e imθ+inϕ . It is noted that in principle, the toroidal mode coupling between m = 0 and m = 0 modes cannot be neglected in a toroidal geometry. However, if one keeps m = 0 modes, the energy conservation is quickly violated due to the ω H mode as shown in Sec. 3.4. On the other hand, if one computes the n = 0 mode using only kinetic trapped electrons responses as in the n = 0 modes, the ω H mode is avoided, but the ambipolar condition and the toroidal angular momentum conservation are not guaranteed as described above. Therefore, we believe that this approximation is somewhat artificial but most physically sound. A drawback of this approximation is that the eigenfrequency of geodesic acoustic modes (GAMs) is changed, while the residual zonal flows are not affected as shown in Sec. 3.4. The influences of this approximation on gyrokinetic simulations will be discussed in Sec. 4.1. The resulting kinetic electron density is defined as



δne,n =0 =



δne,n=0 =

δ f e,t ,n =0 δ([R + ρ ] − x)d6 Z − α p ne

qe φn =0 Te

(10)

,

δ f e,n=0 δ([R + ρ ] − x)d6 Z ,

(11)

where the second term of Eq. (10) denotes an adiabatic passing electrons response, and α p is the flux-surface averaged fraction of passing electrons. The perturbed electron distribution is decomposed into trapped and passing components, δ f e = δ f e,t + δ f e, p , which are defined by contributions from velocity space grids satisfying 12 me v 2 + μ B > μ B max and 1 m v2 2 e 

+ μ B < μ B max , respectively. Here, B max is the maximum magnetic field of the magnetic surface, on which the velocity space grid is defined. By substituting Eqs. (9), (10), and (11), the quasineutrality conditions for n = 0 and n = 0 respectively yield

−∇⊥ ·

ρti2 λ2Di

∇⊥ φn =0 +

αp λ2De

φn =0

 = 4π qi δ f i ,n =0 δ([R + ρ ] − x)d6 Z + qe δ f e,t ,n =0 δ([R + ρ ] − x)d6 Z , −∇⊥ ·

ρti2 λ2Di

∇⊥ φn=0 = 4π



(12)

qs

δ f s,n=0 δ([R + ρ ] − x)d6 Z .

s=i ,e

(13)





Here, if one takes account of the quasi-neutrality in the initial condition, f s0 δ([R + ρ ] − x)d6 Z = 0, full- f like s qs expressions of Eqs. (12) and (13) are obtained by simply replacing δ f s by f s . Eq. (13) keeps all m modes, because finite m components of δ f i and δ f e contribute to m = 0 component of φn=0 through the toroidal mode coupling. After solving Eq. (13), we decompose φn=0 into m = 0 and m = 0 components, and apply only the former to Eq. (8). In the lhs of Eqs. (12) and (13), the gyrokinetic Poisson operators including the ion polarization density and the adiabatic passing electron density are defined using the initial density and temperature. By taking the time derivative and the flux-surface average of Eq. (13) and substituting Eq. (8), the ambipolar condition is derived as

−

ρti2 ∂ E r λ2Di ∂ t

= 4π

 s=i ,e

˙ · ∇ r ) − S src ,s q s δ f s (R

 gf

,

(14)

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

515

where the gyro/flux-surface average operator is defined as,





6

 Ag f =

A (Z)δ(R + ρ − x)d Z

(15)

, f

and · f is the flux-surface average operator. It is noted that in the above kinetic electron model, the turbulent or n = 0 component of the electrostatic potential is determined by using an adiabatic passing electron response. However, in the gyrokinetic equation (8), passing electrons are not forced to be adiabatic, and non-adiabatic passing electrons responses and thus finite passing electrons transport may arise from various kinetic effects such as the resonant wave particle interaction near the mode rational surfaces satisfying ω/k ∼ v te and the interaction between passing and trapped particles through collisional and turbulent pitch angle scattering. The resulting passing electrons fluxes are subject to the ambipolar condition (14), and contribute to evolutions of the total electron distribution f e . Therefore, in this model, we take account not only of the neoclassical transport but also of the turbulent transport of passing electrons induced by passive responses to turbulent fluctuations. 2.3. Conservation laws The gyrokinetic equation (8) yields the following energy conservation law,



Hs

s=e ,i

dE f ld,s ∂ fs 6 dE kin,s dE col dE src + = + d Z= ∂t dt dt dt dt    d 1 dE kin 2 = m s v  + μ B f s d6 Z , dt

dE f ld dt dE col dt dE src dt

s=e ,i

=

dt



2

q s φα

s=e ,i

=

 

∂ fs 6 d Z, ∂t

q s φα C ( f s , f s )d6 Z ,

(16)

(17)

(18)

(19)

s=e ,i s =e ,i

=



H s S src ,s d6 Z .

(20)

s=e ,i

A kinetic energy part of the collision term vanishes due to the energy conservation of the collision operator. However, a field energy part, Eq. (19), remains because of a mismatch of FLR effects between the gyrokinetic equation and the drift-kinetic collision operator. This contribution to the energy balance was shown to be negligible [5]. The toroidal angular momentum conservation law [27] is derived as

     ∂ ms v  bϕ f s   1 ∂  ˙ s v  bϕ f s + · J Rm ∂t J ∂R df df s=e ,i      

∂ Hs qs ˙ + fs − − ms v  bϕ S src ,s df = 0, f s R · ∇ψ ∂ ϕ df c df

(21)

where ψ is the poloidal flux, bϕ is the covariant toroidal component of b, and the drift-kinetic/flux-surface average operator is defined as,





 A df =

A (Z)δ(R − x)d6 Z

.

(22)

f

Here, the collision term is cancelled by its particle and momentum conservation properties, which are defined in the driftkinetic limit. 2.4. Numerical methods The gyrokinetic Poisson bracket operator in Eq. (3) is discretized in the cylindrical coordinates ( R , ζ, Z , v  , μ) using the 4th order non-dissipative conservative finite difference scheme [28,29], where ζ = −ϕ . In the scheme, the Morinishi scheme [30] was extended for the gyrokinetic equation by developing incompressible discretization of the 5D Hamiltonian flows (4) and (5), which satisfies the phase space conservation (6). Since the scheme is a variant of the centered finite difference, it keeps the toroidal symmetry and exactly satisfies the local conservation of f s . The gyrokinetic Poisson equation for n = 0 (12) is computed in the straight-field line coordinates (r , θ ∗ , ϕ ) using a field-aligned Poisson solver [31], while

516

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

Fig. 1. Strong scaling of GT5D on the K-computer at the Riken (Sparc64VIIIfx and Tofu interconnect). The dotted line shows ideal scaling extrapolated from the data point at 24,576 cores. The benchmark test is performed for ITER-size tokamak N R × N ζ × N Z × N v  × N μ = 768 × 64 × 768 × 128 × 32 ∼ 1.5 × 1011 . The maximum performance of ∼682 TFlops (∼6.3 s/step) is achieved using 589,824 cores, and the parallel efficiency is estimated as ∼99.99989% based on the Amdahl’s law.

Eq. (13) is computed in the flux coordinates (r , θ, ϕ ). The cylindrical coordinates and the flux coordinates are connected via the FLR operator [29]. The linear Fokker–Planck collision operator is discretized in the velocity coordinates ( v  , μ) using the 6th order centered finite difference scheme, and a field particle operator is implemented based on Ref. [32], which guarantees an exact conservation of the particle, momentum, and energy (see Appendix A). The time integration is performed using the 2nd order additive semi-implicit Runge–Kutta method (ASIRK) [33] and a stiff linear term involving the parallel streaming is treated implicitly using a generalized conjugate residual (GCR) method [34]. Here, the explicit and implicit operators are separated by perturbed and unperturbed parts of the Hamiltonian (2), and both are given by the Poisson bracket operator (3), which keeps relevant conservation properties at each time integration stage. The above conservation and symmetry properties are essential for robust and accurate computation of the present full- f model. The code is highly parallelized using a multi-dimensional domain decomposition, which is designed based on physical symmetry properties of the gyrokinetic operator (μ symmetry), the collision operator (R symmetry), and the Poisson operator (toroidal symmetry) [35]. We apply the same domain decomposition model both for ions and electrons, and do not use a particle species decomposition, because the cost for solving the implicit operator is an order of magnitude different between ions and electrons, while the memory sizes are comparable. The domain decomposition model is implemented using hybrid parallelization consisting of multi-layer MPI communicators and multi-core OpenMP parallelization. In addition, a novel computation and communication overlap technique [36] is developed using communication threads, which are implemented with a heterogeneous OpenMP programing model. The strong scaling of GT5D is dramatically improved by this latency hiding technique, and on the K-computer, an excellent strong scaling is achieved up to ∼0.6 million cores (see Fig. 1). 3. Verification tests 3.1. Benchmark parameters In this work, we consider deuterium plasmas in a circular concentric tokamak configuration with R 0 /a = 2.79, a/ρti = 150, and q(r ) = 0.85 + 2.18(r /a)2 , which has the following Cyclone like parameters [23] at r s = 0.5a:  = r s / R ∼ 0.18, q(r s ) ∼ 1.4, sˆ (r s ) = [(r /q)dq/dr ]r =r s ∼ 0.78, ne ∼ 4.6 × 1019 m−3 , R 0 / L n = 2.22, T e ∼ T i ∼ 2 keV, R 0 / L te = R 0 / L ti = 6.92. Here, q is the safety factor,  is the local inverse aspect ratio, ne is the electron density, T e and T i are the electron and ion temperature, L n , L te , and L ti are the corresponding scale lengths, and the initial distribution f s0 is given by a local Maxwellian distribution. We use the hybrid electron model with a mass ratio of mi /me = 100, which could lead to a few percent errors normalized electron and ion collisionality parameters at r = r s are √ in the linear growth rate of TEM [37].√The √ ν∗e = (4/3 π )qR 0 νˆ ei / v te  3/2 ∼ 0.034 and ν∗i = (2 2/3 π )qR 0 νˆ ii / v ti  3/2 ∼ 0.024, respectively. It is noted that ν∗e and thus electron neoclassical transport properties are not affected by the heavy electron model. Although various fixed [5] and adaptive [11,12] source models are implemented in GT5D, in this benchmark, we use only a fixed heat source model, S src ,s = νsrc,s A src,s (r )( f M1 − f M2 ). On the other hand, in all linear and nonlinear calculations, a Krook type sink operator, S snk,s = νsnk,s A snk,s (r )( f s − f s0 ), is used to impose a L-mode like fixed boundary condition. Here, νsrc,s and νsnk,s are the heating rates, A src ,s and A snk,s are the deposition profiles, and local Maxwellian distributions with different temperatures, f M1 , f M2 , are chosen to impose fixed power input P in without particle and momentum input. By using the sink profile A snk,s , which is localized at r /a = 0.9 ∼ 1.0, the density, the parallel flow, and the temperature at the edge are imposed as ne ∼ 4 × 1019 m−3 , U e = U i = 0, and T e = T i ∼ 1 keV with νsnk,e = νsnk,i = 0.1v ti /a. The particle, energy, and heat fluxes are defined as s ≡  f s R˙ · ∇ r  g f , Q s ≡ [ms v 2 /2 + μ B ] f s R˙ · ∇ r  g f , and qˆ s ≡ Q s − 5/2T s s , where neoclassical and turbulent components are

˙ respectively. given by fluxes induced by axisymmetric and non-axisymmetric components of the characteristics R, From convergence tests in Sec. 4.1, we have chosen the standard numerical parameters as followings. The gyrokinetic solver uses L R = L Z = 320ρti ( R − R 0 = −1.07a ∼ 1.07a, Z = −1.07a ∼ 1.07a), N R = N Z = 160 ( R =  Z = 2ρti ), where

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

517

Fig. 2. Radial profiles of (a) parallel flows and (b) ion temperature observed in relaxation processes of a two species plasma with 50% of hydrogen and 50% of deuterium, which have different initial profiles.

N p , L p , and  p = L p / N p are the grid number, the system size, and the grid size in the p direction, respectively. The system size and grid width in ζ is changed depending on the problem. Temperature relaxation tests, neoclassical tests, and zonal flow damping tests are performed in the 4D axisymmetric limit with N ζ = 1, linear ITG calculations for a single toroidal mode n use wedge torus configuration with L ζ = 2π /n and N ζ = 16, and nonlinear ITG simulations use a 1/6 wedge torus configuration with L ζ = 2π /6 and N ζ = 32 (n = 0, 6, · · · , 48 or kθ (r s )ρti = 0 ∼ 0.9), which is chosen based on convergence tests in Ref. [11]. The field solver uses N r × N θ = 128 × 256 bi-quadratic spline finite elements, where the field-aligned Poisson solver keeps m = m0 ± 20 poloidal modes around the resonant condition m0 ∼ nq(r ). The velocity space resolution is changed depending on the problem. Since linear ITG calculations show linear eigenfunctions in the velocity space, we use relatively low resolution parameters with ( N v  , N v ⊥ , v ,max / v ts , v ⊥,max / v ts ) = (64, 16, 5, 5). On the other hand, we use high resolution parameters with ( N v  , N v ⊥ , v ,max / v ts , v ⊥,max / v ts ) = (128, 32, 5, 5) in collisionless zonal flow damping tests. This resolution is relaxed by introducing the multi-species collision operator and neoclassical tests are performed using ( N v  , N v ⊥ , v ,max / v ts , v ⊥,max / v ts ) = (80, 20, 5, 5). Based on the convergence of neoclassical tests, numerical parameters of turbulence simulations are set with the same resolution, but with extended v  space, ( N v  , N v ⊥ , v ,max / v ts , v ⊥,max / v ts ) = (96, 20, 6, 5). The FLR operator performs gyro-average using sampling points with N α = 40. The CFL condition is determined by kinetic electrons. From a convergence property of the GCR solver in the implicit part of the ASIRK, the time step width 1 is chosen as t = 2− , where the residual errors of implicit matrices (Eqs. (38) and (39) in Ref. [29]) for ion and electron i calculations are converged to ∼10−15 with ∼10 and ∼150 iterations, respectively. 3.2. Temperature relaxation in multi-species plasmas In order to verify the multi-species collision operator, we compute temperature and momentum relaxation of multispecies plasmas, in which we assume adiabatic electrons and two species ions with uniform density profiles and different parallel flow and temperature profiles at t = 0 (see Fig. 2). Here, we consider axisymmetric plasmas given by the Cyclone like parameters except for the initial density, parallel flow, and temperature profiles. In the simulation, the initial profiles are relaxed towards the same parallel flow and temperature profiles via collisional momentum and energy exchange processes. The collisional energy exchange occurs mainly through the equipartition operator (A.3). On the other hand, the collisional momentum exchange is described by the multi-species linear Fokker–Planck collision operator, because the operator is defined using a non-shifted Maxwellian distribution by assuming the low flow ordering, U i / v ti  1. Fig. 3 shows conservation properties of the multi-species collision operator in temperature relaxation tests for two species plasmas with hydrogen–hydrogen (H-H), hydrogen–deutrium (H-D), and hydrogen–tritium (H-T) ions (50% each). The particle number is exactly conserved for each species, while the momentum and energy exchanges are balanced between two species to satisfy the conservation laws (A.38), (A.39), and (A.40). The momentum and energy exchange processes are saturated in a collision time, which is estimated as τHH ∼ 1.3 ms, τDH ∼ 1.9 ms, and τTH ∼ 2.3 ms, respectively. It should be noted that in GT5D, the collision operator is defined by using a local Maxwellian distribution f Ms , which is given by the moments of f s at each time step, and the perturbed distribution ˜f s defined with respect to f Ms is small in the above tests. Therefore, without the equipartition operator, the temperature relaxation is very weak. Although the equipartition process is normally neglected in δ f approaches, it is essential in full- f simulations. For example, even if one performs numerical experiments with ion heat sources, the equipartition process generates electron heat sources, and determine a total balance of each transport channel and thus temperature difference in a steady state. 3.3. Neoclassical transport In this section, we examine the neoclassical transport in axisymmetric ion–electron plasmas given by the Cyclone like parameters. In this test, we verify the test particle operator (A.4) and the field particle operator (A.34), while the temperature relaxation simulation in the previous section examines mainly the equipartition operator (A.3). In the neoclassical

518

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

Fig. 3. Conservation properties of the multi-species collision operator in temperature relaxation tests for two species plasmas with hydrogen– (H-H), hydrogen–deutrium (H-D), and hydrogen–tritium (H-T) ions (50% each). (a), (b), and (c) show variations of particle number, δ N (t ) = hydrogen    dt d6 Z s C ( f s , f s ), parallel momentum, δ M (t ) = dt d6 Zms v  s C ( f s , f s ), and energy, δ E (t ) = dt d6 Zms v 2 /2 s C ( f s , f s ), due to collisional processes, respectively. The collision time is estimated as τHH ∼ 1.3 ms, τDH ∼ 1.9 ms, and τTH ∼ 2.3 ms, respectively.

Fig. 4. Time evolutions of the ion and electron particle fluxes at r ∼ 0.5a in the neoclassical simulation. In computing the radial electric field, (a) uses a full kinetic electrons response, while (b) uses a kinetic trapped electron response and an adiabatic passing electron response. The ambipolar condition between ion and electron particle fluxes, i = e , is satisfied in (a), while those fluxes converge at different levels in (b).

simulation, both ions and electrons are initialized by local Maxwellian distributions with the same density and temperature profiles, which generate initial radial electric fields due to differences in finite orbit width effects between ions and electrons. The initial radial electric fields show GAM oscillations, and after the damping of GAMs, the system gradually approach a quasi-steady state, where the ion and electron particle fluxes are balanced following the ambipolar condition (20). In order to examine the ambipolar condition, we compare two neoclassical simulations with the new hybrid kinetic electron model and with the conventional trapped kinetic electron model, in which the radial electric field is computed only by kinetic trapped electrons responses. In the former (see Fig. 4(a)), the ambipolar condition, i = e , is satisfied after a collision time τii ∼ 2.1 ms, while in the latter (see Fig. 4(b)), i and e converge at different levels. In both cases, Eq. (8) is solved for a full kinetic electron distribution, and e is computed including passing electrons. From this result, we choose to use a full kinetic elections response to determine the radial electric field E r following Eq. (13). In the quasi-steady state, transport properties are well described by the neoclassical theory. In Fig. 5, the radial profiles of the particle flux, the ion and electron heat fluxes, the radial electric field, and the bootstrap current at t = τii are com-

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

519

Fig. 5. Comparisons of (a) the particle flux, (b) the ion and electron heat fluxes, (c) the radial electric field, and (d) the bootstrap current between the Hirshman–Sigmar’s (H-S) moment approach [38] and the neoclassical simulation with a full kinetic electrons response. Radial profiles at τii ∼ 2.1 ms are plotted.

pared against the Hirshman–Sigmar’s moment approach [38], and reasonable agreements are observed. Here, a theoretical estimation of the radial electric field is given by the neoclassical force balance

Er = −



1

|∇ψ|

A 21 + A 22

− A 21

T e ∂ ln p e qe

∂r

F

 BU i  f

− A 22

T i ∂ ln p i qi

∂r

+ A 23

T e ∂ ln T e qe

∂r

+ A 24

T i ∂ ln T i qi

∂r

 ,

(23)

where p s = ns T s , A i j is given by elements of the matrix A = (L − M)−1 · M, and the friction and viscosity matrices, L, M, are defined for the present ion–electron plasma following Ref. [38]. 3.4. Zonal flow damping In this section, we discuss zonal flow damping tests. We consider axisymmetric plasmas with the Cyclone like parameters except for the density and temperature profiles. The equilibrium profiles of ions and electrons are given by the same uniform density and temperature profiles, and the initial perturbation is given as δ f i = f i0 × 10−5 [1 + cos(π r /a)]/2 following Ref. [29]. Fig. 6 shows comparisons of two kinetic electron models with m = 0 and all m modes. In the former, we apply only a m = 0 component of φn=0 to Eq. (8), while the latter use all m components. In the all m case, m = 0 modes have finite k , and passing electrons produce the ω H mode, which is observed as high frequency oscillations in Fig. 6(a). In order to resolve such a high frequency perturbation, one needs to use an extremely short time step width. However, even with t = 0.5i , the energy conservation breaks down quickly (see Fig. 6(b)). Form the viewpoints of both the computational cost and the numerical accuracy, it is difficult to perform confinement time scale full- f simulations in the presence of the ω H mode. On the other hand, in the m = 0 case, a good energy conservation is satisfied by using much larger time step width, t = 2i (see Fig. 6(c)). A drawback of this approximation is a change in the frequency and damping rate of GAMs. If one takes the m = 0 limit of the Sugama–Watanabe (S-W) theory, electrons responses to GAMs disappear and the frequency (the damping rate) of GAMs become lower (higher). In Fig. 7, this is confirmed by the zonal flow damping tests with adiabatic electrons, and both the m = 0 and all m cases show good agreements with the S-W theory. In the m = 0 case, the zonal flow damping tests with adiabatic and kinetic electrons agree with each other, which is consistent with the lack of passing electrons responses to GAMs. In all three cases, the residual E × B flows are converged to the Rosenbluth–Hinton (R-H) theory [40]. Therefore, influences of GAMs on turbulent transport may be limited. This point will be examined in nonlinear ITG turbulence simulations in Sec. 4.

520

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

Fig. 6. Comparisons of zonal flow damping tests with m = 0 and all m cases. (a) shows time evolutions of E × B flows. The energy conservation is plotted 1 1 for (b) all m and (c) m = 0 cases. The time step widths used are (b) t = 2− and (c) t = 0.5− , respectively. i i

Fig. 7. Comparisons of zonal flow damping tests against the Sugama–Watanabe (S-W) theory [39] and the Rosenbluth–Hinton (R-H) theory [40]. The S-W theory is plotted for m = 0 and all m cases with kr = π /a. Time evolutions of E × B flows are plotted for the zonal flow damping tests with adiabatic electrons (ade) and kinetic electrons (kie) cases. The adiabatic electrons cases are shown for m = 0 and all m cases, while the kinetic electrons case is plotted only for m = 0 case. The frequencies and damping rates show good agreements with the S-W theory for both m = 0 and all m cases, respectively. In all cases, the residual E × B flows are converged to the R-H theory.

3.5. Linear ITG-TEM calculations In this section, we calculate the eigenfrequencies and linear growth rates of ITG-TEMs, and compare them against the Gyrokinetic Toroidal 3D particle code GT3D [41]. GT3D has the capability of calculating ITG-TEMs using a kinetic trapped electron model, and its accuracy was verified against the local ballooning eigenvalue codes FULL, HD7, the Gyrokinetic Toroidal particle Code GTC, and the X-point included Gyrokinetic Code XGC1 in Refs. [25,42]. In the present benchmark calculations, main differences between calculation models of GT5D and GT3D are that GT3D uses a real mass ratio and a Pade approximation for the treatment of the ion polarization density in a short wavelength region. In Fig. 8, GT3D shows a transition from the ITG mode to the TEM at kθ (r s )ρti = nq(r s )ρti /r s ∼ 0.5 (r s = 0.5a), and the unstable region of the TEM spreads over the short wavelength region. Although GT5D also shows the similar transition in the eigenfrequency, the TEM in the short wavelength unstable region is stabilized by a Taylor expansion of the ion polarization density, and the transition point is slightly lower than GT3D. In Fig. 9, the L ti dependence of the most unstable mode at kθ (r s )ρti ∼ 0.28 is plotted. The linear critical gradient of the ITG mode with adiabatic electrons is observed at R 0 / L ti ∼ 4. On the other hand, such a

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

521

Fig. 8. The wavenumber dependency of (a) the eigenfrequency and (b) the linear growth rate is plotted for the Cyclone like parameters (R 0 / L n = 2.22, R 0 / L te = R 0 / L ti = 6.92). The collisionless results obtained from GT5D with adiabatic electrons (ade) and kinetic electrons (kie) are compared against GT3D. The kinetic electrons result with the multi-species collision operator (kie, Collision) is also plotted.

Fig. 9. The L ti dependence of (a) the eigenfrequency and (b) the linear growth rate is plotted for the Cyclone like parameters at kθ (r s )ρti = nq(r s )ρti /r s ∼ 0.28 (r s = 0.5a). The collisionless results obtained from GT5D with adiabatic electrons (ade) and kinetic electrons (kie) are compared against GT3D. Also shown are the collisional case (kie, Collision) and the case where both ion and electron temperature profiles are varied with L te = L ti .

linear critical gradient disappears with kinetic trapped electrons, and the transition from the ITG mode to the TEM occurs for R 0 / L ti < 5.5, when the electron temperature gradient is fixed at R 0 / L te = 6.92. These collisionless results agree well between GT5D and GT3D. On the other hand, when both the ion and electron temperature gradients are decreased with L ti = L te , the TEM is stabilized, and the linear critical gradient of the ITG mode is found at R 0 / L ti ∼ 3.3, which is lower than that with adiabatic electrons. In the collisional case, only the TEM suffers from collisional stabilization [43], and the transition point becomes lower at R 0 / L ti ∼ 4.4. 4. Nonlinear turbulence simulations 4.1. Nonlinear convergence tests In this section, we present an ITG turbulence simulation using the new hybrid kinetic electron model (kinetic electron case), and compare its physics properties against an ITG turbulence simulation with the conventional adiabatic electron model (adiabatic electron case). Simulation parameters are chosen based on the Cyclone like parameters. In both cases, on-axis ion heating with P in = 4 MW and without particle and momentum sources is imposed for r /a < 0.4. The edge temperatures in kinetic and adiabatic electron cases are respectively chosen as T e = T i ∼ 1 keV and ∼0.8 keV, so that the core ion temperature becomes similar values. In fact, the adiabatic case is quoted from a series of power scan in Ref. [13], which use more unstable initial ion temperature profile with R 0 / L ti = 10 and higher collisionality with ν∗i ∼ 0.072. The higher collisionality was needed for keeping the positivity of f i in the higher P in regime, while the collisionality dependence of transport levels was shown to be weak in the collisionless regime [12]. The total computational costs of the kinetic and adiabatic electron cases for ∼10 ms (∼5 × 105 steps with t = 2i for the kinetic electron case and ∼2 × 105 steps with t = 5i for the adiabatic electron case) are respectively ∼4 M CPU hours and ∼0.3 M CPU hours on the Helios supercomputer (SandyBridge-EP processors connected with a fat-tree InfinibandQDR network). Detailed cost comparisons are given in Table 1. The computational cost of the gyrokinetic equation is increased by ∼3.7 times mainly because of the GCR solver, in which the adiabatic electron case computes ∼35 iterations for ions and the kinetic electron case computes ∼10 iterations for ions and ∼150 iterations for electrons, respectively. The computational cost of the collision operator is increased by ∼4.5 times, because the multi-species collision operator computes collision

522

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

Table 1 Typical computational costs of the gyrokinetic equation (Gyrokinetic), the collision operator (Collision), and the gyrokinetic Poisson equation (Field) are compared between the conventional adiabatic electron 1 1 model with t = 5− and the new hybrid kinetic electron model with t = 2− i i . The computational costs are measured for simulations shown in Fig. 10 using 500 nodes of the Helios supercomputer (1000 SandyBridge-EP processors, InfinibandQDR interconnect). Adiabatic electrons

Kinetic electrons

Gyrokinetic (ms/step) Collision (ms/step) Field (ms/step)

791.9 121.2 42.0

2925.4 540.0 70.9

Total (ms/step) 1 Total (ms/− i )

955.1 191.0

3536.2 1768.1

Fig. 10. The plasma profiles observed at t ∼ 10 ms. (a) and (b) show the parallel flow and the ion temperature in the adiabatic election cases. (c), (d), and (e) show the electron density, the parallel flow, and the ion and electron temperature in the kinetic electron cases. Black dotted lines show the initial profiles. The adiabatic electron cases show comparisons between simulations with and without n = 0 and m = 0 modes. The kinetic electron cases show grid convergence between ( N R , N Z , N ζ , N v  , N v ⊥ ) = (160, 160, 32, 96, 20) and (240, 240, 32, 96, 20).

terms for i–i, i–e, e–i, and e–e. The computational cost of the gyrokinetic Poisson equation is increased by ∼1.7 times, since the gyro-averaged electrostatic potential is computed twice for ions and electrons. The resulting total computational cost per time step is increased by ∼3.7 times, which corresponds to the cost increase by ∼9.3 times per unit time.

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

523

Fig. 11. The radial profiles of particle fluxes averaged over t = 0.15 ∼ 10 ms. The passing electrons flux shows outward transport and resonant structures at mode rational surfaces, which are shown by black arrows. The trapped electrons flux shows inward pinch in the outer region.

The first convergence study is to examine influences of n = 0 and m = 0 modes, which are neglected in the new hybrid kinetic electron model. Since kinetic electron cases with n = 0 and m = 0 modes including the ω H mode are prohibitive from the viewpoints of a computational cost and numerical accuracy, we compare adiabatic electron cases, where the radial electric field is defined with all-m modes and only with m = 0 mode. Figs. 10(a) and (b) show comparisons of parallel flows and ion temperature profiles observed in the quasi-steady state. While zonal flow damping rates are different between these cases, turbulent transport and the resulting plasma profiles are not affected by this approximation. In Ref. [12], it was shown that changes in zonal flow damping rates due to collisions affect the ITG turbulence only near critical temperature gradient, and the collisionality scaling of the ITG turbulence is weak at far above critical temperature gradients, because in such a strong turbulence regime, driving effects of zonal flows due to Reynolds stress is almost balanced with nonlinear damping effects, and linear zonal flow damping effects play minor roles. This may also explain that the change in linear zonal flow damping rates in the case with m = 0 mode has minor influences on turbulent transport. Another important convergence study is to check grid convergence. However, the large computational costs of kinetic electron cases prohibit us to repeat the long time numerical experiment with higher grid resolution. In order to resolve this issue, we develop a refinement restart technique, in which (distributed) restart data is converted to new restart data with an arbitrary grid resolution and domain decomposition. This technique is useful also for accelerating slow profile relaxation by using lower resolution computation, when such transient phenomena are unimportant. In the present convergence study, the grid resolution of the restart data at t ∼ 8 ms is converted from ( N R , N Z , N ζ , N v  , N v ⊥ ) = (160, 160, 32, 96, 20) to ( N R , N Z , N ζ , N v  , N v ⊥ ) = (240, 240, 32, 96, 20), and plasma profiles after an ion–ion collision time τii ∼ 2.1 ms are compared. Figs. 10(c), (d) and (e) show comparisons of the electron density, the parallel flow, and the ion temperature observed at t ∼ 10 ms. The results indicate that the kinetic electron case is converged with the standard grid resolution. In the following, we discuss comparisons of physics properties between the kinetic and adiabatic electrons cases. 4.2. Particle transport One of the most significant differences between the kinetic and adiabatic electron cases is the existence of particle transport. Fig. 11 shows particle fluxes averaged over the whole nonlinear phase t = 0.15 ∼ 10 ms. The trapped electrons flux shows inward pinch in the outer region, where a fraction of trapped electrons increases. The inward pinch of trapped electrons may be consistent with the turbulent equipartition (TEP) theory [44,45] and the former quasilinear calculation [46]. On the other hand, the passing electrons flux shows outward transport with an oscillatory profile, which is characterized by large peaks in the core region, finer corrugations around mid minor radius, and a continuous structure in the outer region. These structures are consistent with the distribution of mode rational surfaces satisfying m/n = q(r ), which becomes dense towards the outer region in normal shear tokamaks. It is also found that the positions of several large peaks coincide with relatively low order mode rational surfaces. From these features, we conclude that the passing electrons flux is induced by resonant wave particle interaction at mode rational surfaces satisfying ω/k ∼ v te . This kind of passing electrons transport induced by an interaction between electrons parallel motion and the ITG turbulence was reported also in Ref. [47]. The resonant structure of passing electrons transport produces radially corrugated electron density and temperature profiles in Figs. 10(c) and (e). In Ref. [48], it was reported that this kind of corrugated electron profiles at mode rational surfaces are pronounced with kinetic electrons. 4.3. Momentum transport Another remarkable feature of the kinetic electron case is that the direction of parallel flows are reversed from the adiabatic electron case. Since the present simulations do not have momentum input in the plasma core, this difference

524

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

Fig. 12. The ion component of the toroidal angular momentum conservation in (a) the adiabatic electron case and (b) the kinetic electron case. The time average data for the initial flow generation phase t = 0.15 ∼ 2 ms is plotted. “torque”, “stress”, “field”, “current”, and “source” respectively correspond to the first, second, third, fourth, and fifth terms of Eq. (21), in which “collision” is exactly cancelled between ions and electrons.

Fig. 13. The radial electric fields ((a), (b)) and the corresponding E × B shearing rates ((c), (d)) in the adiabatic electron case ((a), (c)) and the kinetic electron case ((b), (d)). The time average data for t = 8 ∼ 10 ms is plotted. The E r profile observed in GT5D (GT5D) and that computed by the neoclassical force balance (23) (NC) are compared.

is determined by turbulent momentum transport. In order to understand the flow generation process, we show the ion component of the toroidal angular momentum conservation (21) averaged over the initial flow generation phase (t = 0.15 ∼ 2 ms) in Fig. 12. In both cases, the toroidal angular momentum conservation is well satisfied and errors are negligibly small. Compared with the adiabatic electron case, the kinetic electron case is characterized by a large contribution from the field term  f s ∂ϕ H s , which is interpreted as the off-diagonal component of the Maxwell stress tensor in the no-FLR limit [49]. The result suggests that the phase relation between the perturbed ion distribution δ f i and the toroidal electric field ∂ϕ φα is changed by kinetic electrons. The other qualitative differences are the radial ion current induced by the particle transport, and the collisional momentum transfer, which is balanced with that of electrons by the momentum conservation of the collision operator. 4.4. Radial electric fields The above qualitative differences in plasma profiles affect the radial electric field, which is one of the most important suppression mechanisms of turbulent transport. Fig. 13 shows the time averaged profiles of the radial electric field E r and the corresponding shearing rate ω E × B ≡ (cr /qB 0 )∂r (qE r /r ). In both cases, the E r profiles in the initial phase (the linear

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

525

Fig. 14. The spatio-temporal evolutions (a) the ion energy flux in the adiabatic electron case, and (b) the ion energy flux and (c) the electron energy flux in the kinetic electron case. The flux is normalized by the ion gyro-Bohm diffusivity χGB = v ti ρti2 / L n .

phase after the initial GAM damping) is characterized by the similar global negative profiles, while their amplitudes are different because of the different initial ion temperature profiles in Figs. 10(b) and (e). However, their nonlinear evolutions are quite different. In Fig. 13(a), the adiabatic electron case keeps the global negative profile, and a fine zonal structure is generated in the core region. The global negative structure is consistent with the neoclassical force balance (23), while in Fig. 13(c), the correlation between the fine zonal structure and the neoclassical force balance is weak in the core region. Therefore, it is considered that the fine zonal structure is generated by turbulence driven zonal flows. The magnitudes of global and zonal E r shear are comparable to the linear growth rate, ω E × B ∼ γ . On the other hand, in Fig. 13(b), the E r profile in the kinetic electron case is reversed from the initial one, and a fine zonal structure is observed over significant radii. Since the global positive profile agrees with the neoclassical force balance, the reverse of the global structure is determined by the sign of parallel flows. In Fig. 13(d), it is shown that unlike the adiabatic electron case, the fine zonal structure is explained by the neoclassical force balance, where the fine E r shear is almost balanced with the gradient of corrugated electron density profile. The resulting E r shear is dominated by this fine structure, and greatly exceeds the linear growth rate, ω E × B  γ . 4.5. Heat transport In both cases, the heat transport is characterized by active bursts of avalanche like non-local transport. In Fig. 14, spatiotemporal evolutions of the energy flux are showing such avalanche like non-local transport. Although they are characterized

526

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

Fig. 15. The radial profiles of (a) the ion heat diffusivity χi in the adiabatic electron case, and (b) the ion heat diffusivity χi and the electron heat diffusivity χe in the kinetic electron case. In (b), the passing and trapped electron components of χe are also plotted. The time average data for t = 8 ∼ 10 ms is plotted.

Fig. 16. The radial profiles of (a) R 0 / L ti in the adiabatic electron case, and (b) R 0 / L ti and R 0 / L te in the kinetic electron case. The time average data for t = 8 ∼ 10 ms is plotted.

by ballistic propagation with the similar propagation speed and width, their directions are different between the kinetic and adiabatic electron cases. As reported in Ref. [5], the radial electric field structure determines the direction of avalanches propagation. In the adiabatic electron case, the avalanches propagate outward (inward) in the outer (core) region, where the shear of global E r structure is positive (negative). On the other hand, in the kinetic electron case, the avalanches propagate inward (outward) in the outer (core) region. In addition, the electron energy flux in the core region shows smaller propagation widths, reflecting resonant structures of the electron particle flux. The time averaged heat diffusivity is shown in Fig. 15, where the definition of the heat diffusivity is modified as χs ≡ qˆ s /(ne | T s /r |) with  T s = T s (r − 0.1a) − T s (r + 0.1a) and r = 0.2a, so that the divergence of χs with corrugated electron temperature profiles is avoided. The time average is taken for t = 8 ∼ 10 ms, where the adiabatic electron case is fully saturated to satisfy the power balance, while the kinetic electrons case is not fully saturated (the total energy flux is ∼3.5 MW at r /a ∼ 0.5 against P in = 4 MW). In the adiabatic electron case, χi increases towards the edge, while the kinetic electron case shows a flat χi profile in the source free region. This difference is produced by the reduction of the ion energy flux due to slow evolution of the ion temperature profile in the outer region. In the kinetic electron case, ∼20% of the energy flux is carried by electrons, which is consistent with collisional energy transfer from ions to electrons at the different temperatures in Fig. 10(e). Both trapped and passing electrons contributions are comparable, while the former becomes dominant in the outer region, where the trapped electrons fraction increases. The resulting temperature gradient profiles are shown in Fig. 16. In spite of the above physics differences, in the source free region for r /a = 0.5 ∼ 0.8, both cases show the similar ion temperature gradients with nonlinear upshift from the linear critical gradients, which are limited by a nonlinear critical gradient R 0 / L ti ∼ 6 [13]. 5. Summary In this work, we have developed a full- f gyrokinetic model with kinetic electrons for simulating the ITG-TEM turbulence at the ion scale. The model consists of the multi-species full- f gyrokinetic equations with the new hybrid kinetic electron model and the multi-species linear Fokker–Planck collision operator. In the new hybrid kinetic electron model, turbulent fluctuations are determined by the non-axisymmetric part of the gyrokinetic Poisson equation with kinetic trapped electrons responses, while all the remaining part is computed using a full kinetic electron model including both trapped and passing electrons. By using this approach, the ω H mode is eliminated from electrostatic full- f gyrokinetic simulations with keeping important physics such as the ITG-TEM turbulence, the neoclassical transport, the ambipolar condition, particle trapping

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

527

and de-trapping processes, and passing electrons transport due to passive responses to turbulent fluctuations. Although the model is designed by paying special attention to the neoclassical theory, which is essential in full- f approaches, it can be applied also to δ f simulations. The model is verified through a series of benchmark calculations on collisional temperature relaxation in multi-species plasmas, neoclassical electron transport, zonal flow damping, and linear ITG-TEM stability, and it has been demonstrated that the model successfully captures essential physics ingredients for full- f ITG-TEM simulations. Another important point is that the model enables full- f ITG-TEM simulations with reasonable cost increase, about ∼10 times from the conventional adiabatic electron model. Flux driven ITG turbulence simulations with the Cyclone like parameters are performed using the new hybrid kinetic electron model and the conventional adiabatic electron model, and from their comparisons, influences of kinetic electrons are clarified. It is found that even in the ITG turbulence, passing electrons transport is important, and its resonant structure produces corrugated electron density and temperature profiles, which lead to fine zonal E r structures through the neoclassical force balance. It should be stressed that this E × B flow generation mechanism is qualitatively different from turbulence driven zonal flows in the adiabatic electron case. Another remarkable feature is that turbulent momentum transport leads to parallel flows and the resulting E r profile in the opposite direction. This is mainly attributed to the field term in the toroidal angular momentum conservation. In spite of these qualitative differences in the electron density, the electron temperature, the parallel flow, and the radial electric field, the resulting ion temperature gradients are sustained below the similar nonlinear critical gradient at R / L ti ∼ 6. This may suggest the validity of former scaling studies on the ion heat transport and the energy confinement time based on ITG turbulence simulations with adiabatic electrons. However, former results on momentum transport and intrinsic rotation may be changed, when kinetic electrons are taken into account. More detailed physics studies on the ITG-TEM turbulence is out of scope of this work and will be reported in elsewhere. The above simulation results disclosed that in addition to trapped electrons, passing electrons are important in the ITG turbulence. This is qualitatively consistent with the flux-tube calculations in Ref. [47]. However, in the hybrid kinetic electron model, passing electrons transport is computed using passive responses to turbulent fluctuations, and a quantitative impact of this approximation is still need to be estimated. In the future work, this issue will be addressed by comparing passing electrons transport between δ f flux-tube simulations with the hybrid kinetic electron model and with the full kinetic electron model. Another issue is that the hybrid kinetic electron model cannot treat a contribution from the ETG mode, which is driven by passing electrons. This limitation should be carefully taken into account when experimental data is analyzed via the hybrid kinetic electron model. Acknowledgements The computation in this work was performed on the BX900 at the JAEA, the Helios at the IFERC, the FX10 at the Univ. Tokyo, and the K-computer (hp150027) at the Riken. The author would like to thank Dr. S. Jolliet for implementing the multi-species linear Fokker–Planck collision operator on GT5D. He is also grateful to Drs. Y. Asahi, S. Maeyama, and M. Nakata for providing useful information on kinetic electrons transport in the GKV code and to Dr. N. Hayashi for his support to this work. This work is supported by the MEXT, Grant for HPCI Strategic Program Field No. 4: Next-Generation Industrial Innovations, Grant for Post-K priority issue No. 6: Development of Innovative Clean Energy, and Grant No. 22686086. Appendix A. Multi-species collision operator We consider a collision operator in the drift-kinetic limit by neglecting FLR corrections, and therefore, the distribution function in the followings can be replaced by the guiding-center distribution function. In this approximation, local conservation of density, momentum, and energy, which is critical for full- f approaches, can be strictly imposed. The general collision operator between species a and b may be written as

C ab ( f a , f b ) = C ab ( f Ma , f Mb ) + C ab ( ˜f a , f Mb ) + C ab ( f Ma , ˜f b ) + C ab ( ˜f a , ˜f b ),

(A.1)

where both distribution functions f a and f b are decomposed into f a = f Ma + ˜f a and f b = f Mb + ˜f b using the Maxwellian distribution functions f Ma and f Mb . If ˜f a and ˜f b are small such that we can neglect the nonlinear term C ab ( ˜f a , ˜f b ), the last

E term in Eq. (A.1) can be neglected, and the linearized collision operator consists of the equipartition operator C ab , the test T F particle operator C ab , and the field particle operator C ab :

L E T ˜ F ˜ C ab ( f a , f b ) = C ab ( f Ma , f Mb ) + C ab ( f a ) + C ab ( f b ),

(A.2)

E where C ab does not vanish, when species a and b have different temperature,



E C ab ( f Ma ,

f Mb ) = −2 f Ma 1 −

Ta Tb





−1

νˆ ab αab

2 )− φˆ  (xb )(1 + αab

ˆ xb ) φ( xb

 .

(A.3)

Here, the operator is defined using non-shifted Maxwellian distribution functions by assuming the low flow ordering. It should be noted that this operator is normally neglected in δ f approaches. However, in GT5D, the Maxwellian distribution functions f Ma and f Mb are defined by density and temperature at each time step, and this operator dictates the lowest order energy exchange among multiple species.

528

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531 T C ab is given by a test particle operator proposed in Ref. [22],

T ˜ T0 C ab ( f a ) = Qab C ab Qab ˜f a .

(A.4)

T0 v Here, C ab consists of the pitch angle scattering operator L and the energy scattering operator C ab ,

T0 v C ab ( g ) = ν Dab ( v )L( g ) + C ab ( g ),

1

L( g ) =



2

1

v2



sin θ v

νab ( v )

∂v

2

xa3

νab ( v ) = 2νˆ ab

G (xb ) xa3

∂g ∂θ v

v 4 f Ma

ˆ xb ) − G (xb ) φ(

ν Dab ( v ) = νˆ ab

νˆ ab ≡

∂

sin θ v ∂θ v 1 ∂

v (g) = C ab



 + 

∂ ∂v

2

1

∂ g

sin2 θ v ∂ ϕ v 2 g f Ma

(A.5)

 ,

(A.6)



,

(A.7)

(A.8)

,

(A.9)

,

4π nb qb2 qa2 ln ab ma2 v 3T a

(A.10)

,

where g denotes an arbitrary function of the velocity coordinates ( v , θ v , ϕ v ), ab is the Coulomb logarithm, xb ( v ) = v / v T b , √ αab = v T a / v T b , v T b = 2T b /mb , and

x

2

ˆ x) = √ φ(

π

d ye − y , 2

(A.11)

0

G (x) =

ˆ x) − xφˆ  (x) φ(

(A.12)

.

2x2

The operator Qab is defined as

Qab g = g + (θab − 1)Pa g , where



⎡

1 ma

Ta

θab = ⎣ 

Ta ma

+

1 mb

Tb mb

+

(A.13)

 ⎤1/2

 ⎦

(A.14)

,

and the projection operator Pa = P1a + P2a is given as

P1a g = f Ma

ma

P2a g = f Ma ua [ g ] =

δ Ta[g] Ta

Ta

δ Ta[g]



1



(A.15)

ma v 2

Ta

2T a

−

3



2

(A.16)

,

d3 v v  g ,

na

=

ua [ g ] v  ,

1 na

d3 v

2



(A.17) ma v 2

3

2T a

−

3 2

 g.

(A.18)

T0 In the velocity coordinates ( v  , μ, α ), which is used in GT5D, C ab is written as

T0 ˜ T0 C ab ( f a ) = C ab (1) ˜f a + νab

2˜ 2˜ ˜ ∂ ˜f a ∂ 2 ˜f a ab ∂ f a ab ∂ f a ab ∂ f a + ν⊥ + νab2 + ν⊥ + ν⊥ , 2 2 ∂ v ∂s ∂ s∂ v  ∂ v ∂ s2

(A.19)

where s = 2μ B /ms , v 2 = v 2 + s, and T0 C ab (1) = 2νˆ ab αab φˆ  (xb ), ab

ν = v  F ab (xb ),

(A.20) (A.21)

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

2sH ab (xb )

ν⊥ab = 2sF ab (xb ) +

+ 4G ab (xb ),

v2

529

(A.22)

2

νab2 =

v

H ab (xb ) + G ab (xb ), v2 4sv  ab ν⊥ = 2 H ab (xb ), v 2 4s ν⊥ab2 = 2 H ab (xb ) + 4sG ab (xb ), v    ma  ˆ xb ) − xb φˆ  (xb ) , F ab (xb ) = ν0ab ( v ) 1 + φ( mb

H ab (xb ) =

2 1 2

(A.24) (A.25) (A.26)





1

G ab (xb ) =

(A.23)

ˆ xb ) − G (xb ) , ν0ab ( v ) v 2 φ(

(A.27)





ˆ xb ) − 3G (xb ) . ν0ab ( v ) v 2 φ(

(A.28)

T ˜ T0 ˜ In the single species limit, one has C aa ( f a ) = C aa ( f a ) and F aa = 0, and Eq. (A.19) recovers the single species operator in Ref. [5]. F The original definition of the field particle operator C ab is given by

F ˜ T T C ab ( f b ) = − V ab [ ˜f b ]C ab ( f Ma ma v  / T a ) − W ab [ ˜f b ]C ab ( f Ma (xa2 − 3/2)),

where

Tb V ab [ ˜f b ] ≡



γab Tb

W ab [ ˜f b ] ≡

ηab



d3 v

γab = T a

d v

f Mb



Ta

T C ba

˜f b

d3 v ma v 



˜f b

3

f Mb

f Mb

mb v  Tb

 (A.30)

,

T C ba ( f Mb (xb2 − 3/2)),



T C ab

f Ma

ma v 

(A.31)

 (A.32)

,

Ta

T d3 vxa2 C ab ( f Ma xa2 ).

ηab = T a

(A.29)

(A.33)

The above field particle operator guarantees the conservation of density, momentum, and energy. However, accumulation of numerical errors due to finite grid resolution and velocity space truncation may lead to unphysical evolutions of plasma profiles. To avoid such numerical issues, we introduce a modified field particle operator following Ref. [32],





F ˜ C˜ ab ( f b ) = f Ma ( v ) A ab I ab ( v ) + v  B ab J ab ( v ) + C ab K ab ( v ) ,

(A.34)

I ab ( v ) = 1 − K ab ( v ),

(A.35)

where

J ab ( v ) = K ab ( v ) =

3

√

π 2G (xb )

4

3

√ 4

xa

π

2

αab xb

αab (θab − 1) + , 2 3/2 (1 + αab )

(A.36)

 2α (θ − 1) ˆ xb ) − xb φˆ  (xb )(1 + α 2 ) + ab ab (xa2 − 3/2). φ( ab 2 3/2 (1 + αab )



(A.37)

J ab ( v ) and K ab ( v ) are equal, up to a constant factor, to the original field particle operator, while I ab ( v ) is chosen such that the energy conservation is not modified by a correction on the particle number. A ab , B ab , C ab are coefficients to be determined at each time step through the conservation laws:



F ˜ d3 v C˜ ab ( fb) = −







E T ˜ d3 v C ab ( f Ma , f Mb ) + C ab ( fa) ,

F ˜ d3 vma v  C˜ ab ( fb) = −



1 F ˜ d3 v ma v 2 C˜ ab ( fb) = − 2

d3 vmb v 





 E T ˜ C ba ( f Mb , f Ma ) + C ba ( fb) ,

  1 E T ˜ d3 v mb v 2 C ba ( f Mb , f Ma ) + C ba ( fb) , 2

which can be imposed by solving the following linear system at each grid:

(A.38) (A.39) (A.40)

530

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

 ab = Y ba , M ab X

(A.41)

where

M ab =

⎛

 ab X

⎛

I ab ( v ) d3 v f Ma ( v ) ⎝ ma v  I ab ( v ) 1 m v 2 I ab ( v ) 2 a

⎞

A ab = ⎝ B ab ⎠ , C ab

⎛

Y ba = −



⎜ ⎜ ⎝

⎞

J ab ( v ) ma v  J ab ( v ) 1 m v 2 J ab ( v ) 2 a

K ab ( v ) ma v  K ab ( v ) ⎠ , 1 m v 2 K ab ( v ) 2 a

(A.42)

(A.43)





E T ˜ C ab ( f Ma , f Mb ) + C ab ( fa)

⎞

 ⎟ ⎟ ⎟. ⎠  1 E T ˜ 2 m v C ba ( f Mb , f Ma ) + C ba ( f b ) 2 b 

E T ˜ d3 v ⎜ mb v  C ba ( f Mb , f Ma ) + C ba ( fb)

(A.44)

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]

[11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23]

[24] [25] [26] [27] [28] [29] [30] [31] [32]

A.J. Brizard, T.S. Hahm, Foundations of nonlinear gyrokinetic theory, Rev. Mod. Phys. 79 (2007) 421–468. W.W. Lee, Gyrokinetic particle simulation model, J. Comput. Phys. 72 (1987) 243–269. Y. Idomura, T.H. Watanabe, H. Sugama, Kinetic simulations of turbulent fusion plasmas, C. R. Phys. 7 (2006) 650–669. X. Garbet, Y. Idomura, L. Villard, T.H. Watanabe, Gyrokinetic simulations of turbulent transport, Nucl. Fusion 50 (2010) 043002. Y. Idomura, H. Urano, N. Aiba, S. Tokuda, Study of ion turbulent transport and profile formations using global gyrokinetic full- f Vlasov simulation, Nucl. Fusion 49 (2009) 065029. G. Dif-Pradalier, V. Grandgirard, Y. Sarazin, X. Garbet, Ph. Ghendrih, Interplay between gyrokinetic turbulence, flows, and collisions: perspectives on transport and poloidal rotation, Phys. Rev. Lett. 103 (2009) 065002. S. Ku, J. Abiteboul, P.H. Diamond, G. Dif-Pradalier, J.M. Kwon, Y. Sarazin, T.S. Hahm, X. Garbet, C.S. Chang, G. Latu, E.S. Yoon, Ph. Ghendrih, S. Yi, A. Strugarek, W. Solomon, V. Grandgirard, Physics of intrinsic rotation in flux-driven itg turbulences, Nucl. Fusion 52 (2012) 063013. P.H. Diamond, T.S. Hahm, On the dynamics of turbulent transport near marginal stability on the dynamics of turbulent transport near marginal stability on the dynamics of turbulent transport near marginal stability, Phys. Plasmas 2 (1995) 3640. Y. Idomura, Full- f gyrokinetic simulation over a confinement time, Phys. Plasmas 21 (2014) 022517. Y. Sarazin, V. Grandgirard, J. Abiteboul, S. Allfrey, X. Garbet, Ph. Ghendrih, G. Latu, A. Strugarek, G. Dif-Pradalier, P.H. Diamond, S. Ku, C.S. Chang, B.F. McMillan, T.M. Tran, L. Villard, S. Jolliet, A. Bottino, P. Angelino, Predictions on heat transport and plasma rotation from global gyrokinetic simulations, Nucl. Fusion 51 (2011) 103023. S. Jolliet, Y. Idomura, Plasma size scaling of avalanche-like heat transport in tokamaks, Nucl. Fusion 52 (2012) 023026. M. Nakata, Y. Idomura, Plasma size and collisionality scaling of ion temperature gradient driven turbulence, Nucl. Fusion 53 (2013) 113039. Y. Idomura, M. Nakata, Plasma size and power scaling of ion temperature gradient driven turbulence, Phys. Plasmas 21 (2014) 020706. C.S. Chang, S. Ku, P.H. Diamond, Z. Lin, S. Parker, T.S. Hahm, N. Samatova, Compressed ion temperature gradient turbulence in diverted tokamak edge, Phys. Plasmas 16 (2009) 056108. D. Zarzoso, Y. Sarazin, X. Garbet, R. Dumont, A. Strugarek, J. Abiteboul, T. Cartier-Michaud, G. Dif-Pradalier, Ph. Ghendrih, V. Grandgirard, G. Latu, C. Passeron, O. Thomine, Impact of energetic-particle-driven geodesic acoustic modes on turbulence, Phys. Rev. Lett. 110 (2013) 125002. K. Miki, Y. Idomura, Finite-orbit-width effects on energetic-particle-induced geodesic acoustic mode, J. Plasma Fusion Res. 10 (2015) 3403068. F. Jenko, Massively parallel Vlasov simulation of electromagnetic drift-wave turbulence, Comput. Phys. Commun. 125 (2000) 196–209. A. Bottino, A.G. Peeters, O. Sauter, J. Vaclavik, L. Villard, ASDEX Upgrade Team, Simulations of global electrostatic microinstabilities in ASDEX upgrade discharges, Phys. Plasmas 11 (2004) 198. Y. Idomura, S. Tokuda, Y. Kishimoto, Gyrokinetic simulations of tokamak micro-turbulence including kinetic electron effects, J. Plasma Fusion Res. 6 (2004) 17. T. Vernay, S. Brunner, L. Villard, B.F. McMillan, S. Jolliet, A. Bottino, T. Gorler, F. Jenko, Global gyrokinetic simulations of TEM microturbulence, Plasma Phys. Control. Fusion 55 (2013) 074016. Z. Lin, Y. Nishimura, Y. Xiao, I. Holod, W.L. Zhang, L. Chen, Global gyrokinetic particle simulations with kinetic electrons, Plasma Phys. Control. Fusion 49 (2007) B163–B172. H. Sugama, T.-H. Watanabe, M. Nunami, Linearized model collision operators for multiple ion species plasmas and gyrokinetic entropy balance equations, Phys. Plasmas 16 (2009) 112503. 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, J. Weiland, Comparisons and physics basis of tokamak transport models and turbulence simulations, Phys. Plasmas 7 (2000) 969–983. S. Maeyama, Y. Idomura, M. Nakata, M. Yagi, N. Miyato, T.-H. Watanabe, A. Ishizawa, M. Nunami, Multi-scale plasma turbulence simulations with real ion-to-electron mass ration and beta value, Phys. Rev. Lett. 114 (2015) 255002. G. Rewoldt, Z. Lin, Y. Idomura, Linear comparison of gyrokinetic codes with trapped electrons, Comput. Phys. Commun. 177 (2007) 775. Y. Camenen, Y. Idomura, S. Jolliet, A.G. Peeters, Consequences of profile shearing on toroidal momentum transport, Nucl. Fusion 51 (2011) 073039. Y. Idomura, Accuracy of momentum transport calculations in full-f gyrokinetic simulations, Comput. Sci. Discov. 5 (2012) 014018. Y. Idomura, M. Ida, S. Tokuda, L. Villard, New conservative gyrokinetic full- f Vlasov code and its comparison to gyrokinetic δ f particle-in-cell code, J. Comput. Phys. 226 (2007) 244–262. Y. Idomura, M. Ida, T. Kano, N. Aiba, S. Tokuda, Conservative global gyrokinetic toroidal full- f five-dimensional Vlasov simulation, Comput. Phys. Commun. 179 (2008) 391–403. Y. Morinishi, T.S. Lund, O.V. Vasilyev, P. Moin, Fully conservative higher order finite difference schemes for incompressible flow, J. Comput. Phys. 143 (1998) 90. S. Jolliet, B.F. McMillan, L. Villard, T. Vernay, P. Angelino, T.M. Tran, S. Brunner, A. Bottino, Y. Idomura, Parallel filtering in global gyrokinetic simulations, J. Comput. Phys. 231 (2012) 745–758. S. Satake, R. Kanno, H. Sugama, Development of a non-local neoclassical transport code for helical configurations, J. Plasma Fusion Res. 3 (2008) S1062.

Y. Idomura / Journal of Computational Physics 313 (2016) 511–531

531

[33] X. Zhong, Additive semi-implicit Runge–Kutta methods for computing high-speed nonequilibrium reactive flows, J. Comput. Phys. 128 (2007) 19–31. [34] S.C. Eisenstat, H.C. Elman, M. Schultz, Variational iterative methods for nonsymmetric systems of linear equations, SIAM J. Numer. Anal. 20 (1983) 345–357. [35] Y. Idomura, S. Jolliet, Performance evaluations of gyrokinetic Eulerian code GT5D on massively parallel multi-core platforms, in: Proceedings of SC11, 2011. [36] Y. Idomura, M. Nakata, S. Yamada, M. Machida, T. Imamura, T.-H. Watanabe, M. Nunami, H. Inoue, S. Tsutsumi, I. Miyoshi, N. Shida, Communicationoverlap techniques for improved strong scaling of gyrokinetic Eulerian code beyond 100k cores on the k-computer, Int. J. High Perform. Comput. Appl. 28 (2014) 73–86. [37] A. Bottino, Modelling of global electrostatic microinstabilities in tokamaks: effects of E × B flow and magnetic shear, PhD thesis, Ecole Polytechnique Fédérale de Lausanne, 2004. [38] S.P. Hirshman, D.J. Sigmar, Neoclassical transport of impurities in tokamak plasmas, Nucl. Fusion 21 (1981) 1079. [39] H. Sugama, T.-H. Watanabe, Collisionless damping of zonal flows in helical systems, Phys. Plasmas 13 (2006) 012501. [40] M.N. Rosenbluth, F.L. Hinton, Poloidal flow driven by ion-temperature-gradient turbulence in tokamaks, Phys. Rev. Lett. 80 (1998) 724. [41] Y. Idomura, S. Tokuda, Y. Kishimoto, Global gyrokinetic simulation of ion temperature gradient driven turbulence in plasmas using a canonical Maxwellian distribution, Nucl. Fusion 43 (2003) 60167. [42] I. Holod, Z. Lin, Verification of electromagnetic fluid-kinetic hybrid electron model in global gyrokinetic particle simulation, Phys. Plasmas 20 (2013) 032309. [43] G. Rewoldt, W.M. Tang, Toroidal microinstability studies of high-temperature tokamaks, Phys. Fluids B 2 (1990) 318. [44] M.B. Isichenko, A.V. Gruzinov, P.H. Diamond, Invariant measure and turbulent pinch in tokamaks, Phys. Rev. Lett. 74 (1995) 4436. [45] X. Garbet, L. Garzotti, P. Mantica, H. Nordman, M. Valovic, H. Weisen, C. Angioni, Turbulent particle transport in magnetized plasmas, Phys. Rev. Lett. 91 (2003) 035001. [46] C. Bourdelle, X. Garbet, F. Imbeaux, A. Casati, N. Dubuit, R. Guirlet, T. Parisot, A new gyrokinetic quasilinear transport model applied to particle transport in tokamak plasmas, Phys. Plasmas 14 (2007) 112501. [47] J. Dominski, S. Brunner, T. Gorler, F. Jenko, D. Told, L. Villard, How non-adiabatic passing electron layers of linear microinstabilities affect turbulent transport, Phys. Plasmas 22 (2015) 062303. [48] R.E. Waltz, M.E. Austin, K.H. Burrell, J. Candy, Gyrokinetic simulations of off-axis minimum-q profile corrugations, Phys. Plasmas 13 (2006) 052301. [49] J. Abiteboul, X. Garbet, V. Grandgirard, S.J. Allfrey, Ph. Ghendrih, G. Latu, Y. Sarazin, A. Strugarek, Conservation equations and calculation of mean flows in gyrokinetics, Phys. Plasmas 18 (2011) 082503.