23P140-B
Gyrokinetic simulations in helical plasmas with the hybrid method ofsemi-Lagrangian and semi-implicit Runge-Kutta schemes
セミラグランジアン法と半陰的ルンゲクッタ法の複合解法を用いたヘリカル
プラズマのジャイロ運動論的シミュレーション
Shinya Maeyama1, Akihiro Ishizawa2, Tomohiko Watanabe2, Noriyoshi Nakajima2, Shunji
Tsuji-Iio1, Hiroaki Tsutsui1
1
2
前山伸也 ，石澤明宏 ，渡邉智彦2，中島徳嘉2，飯尾俊二1，筒井広明1
1

Tokyo Institute of Technology
2-12-1, Ookayama, Meguro-ku, Tokyo 152-8550, Japan
東京工業大学 〒152-8550 東京都目黒区大岡山2-12-1
2
National Institute for Fusion Science322-6, Oroshi-cho, Toki, Gifu 509-5292, Japan核融合科学研究所 〒509-5292 岐阜県土岐市下石町322-6
A hybrid method of semi-Lagrangian and additive semi-implicit Runge-Kutta schemes for gyrokinetic f
Vlasov simulations is presented. This method is free from the Courant-Friedrichs-Lewy condition for thelinear terms in the gyrokinetic equation. Simulations of parallel dynamics and of the ion-temperaturegradient instability in fusion plasmas confined by helical magnetic fields are carried out by means of thismethod. It is demonstrated that their results show good agreements with those obtained by using theexplicit Runge-Kutta-Gill scheme, while the new numerical method has no time-step restrictions for thelinear terms and substantially reduces the computational cost.

1. Introduction
Three-dimensional configurations of helicalplasmas bring numerical difficulties to gyrokinetic simulations of micro-instabilities. To
investigate the ion-temperature-gradient (ITG)driven turbulence in helical plasmas, Watanabe etal. employ a large number of grid points along amagnetic field line in the gyrokinetic Vlasov
simulation code GKV [1]. As a result, theCourant-Friedrichs-Lewy (CFL) condition on theparallel advection severely restricts time steps. To
address this problem, we propose an efficient
numerical method for linear analysis of microinstabilities in helical plasmas by means of gyrokinetic Vlasov simulations, employing semiLagrangian and additive semi-implicit RungeKutta schemes. The new scheme is free from theCFL restrictions for the linear terms.
2. Equations and schemes
Employing the flux tube model, the linearizedgyrokinetic Vlasov equation for the perturbed ion
gyrocenter distribution function fk(z,v||,) in theelectrostatic limit is given by || B    v||  ||  i d 
 fkmi v||  t(1)eF  M v||  ||  i d  i* J 0 (k   ti )k ,Tiwhere B, mi, e, FM, Ti, d, ∗ and ti are themagnetic field strength, ion mass, elementarycharge, Maxwell distribution function, ion
temperature, magnetic and diamagnetic drift
frequencies and ion Larmor radius, respectively.
The radial and poloidal wavenumber kx and ky, thefield-aligned coordinate z, the parallel velocity v||
and the magnetic moment  are employed as thephase space coordinates. The electrostatic potential
k is given by the quasi-neutrality condition with anadiabatic electron response.
By means of the operator splitting method [2], Eq.
(1) can be split into the parallel motions and theothers. The former consists of two linear advection
equations in z and v||. They are easily computed bya semi-Lagrangian scheme [3]. Practically, one has
to evaluate the value of fk by using one-dimensional
interpolations. For more details, see Ref. [4]. The
latter is regarded as an additive operator of theperpendicular drift, source and collision terms. To
solve this problem, we employ additive semiimplicit Runge-Kutta schemes (ASIRK) [5] andtreat the magnetic drift term implicitly. Since thecoefficient matrix of the magnetic drift term isdiagonal, one can easily compute its semi-implicit
time integration without using matrix solvers.

Fig 1. (a) Contour lines of the particle kinetic energy mv||2/2+B with the helical field, and (b)-(d)Snapshots of equi-contours of the distribution function f(z,v||) in parallel phase space (whereB0/Ti = 4.0). Horizontal and vertical axes are defined by z and v||, respectively.
3. Parallel dynamics in helical plasmas
Here we compute only the parallel dynamics by asemi-Lagrangian scheme with a second-order
operator splitting method. We employ 192×256 grid
points and a time step size t/ttr = 0.1, where ttr=Ln/vti with the density scale length Ln and the ion
thermal velocity vti.The initial profile is given by f(z,v||;t=0) = FM(mv||2/2+B)(1+cosz), and theperiodic boundary condition is employed in z.
Characteristic curves of the dynamics are given as
contour lines of the particle kinetic energy, as
shown in Fig. 1 (a). There are trajectories ofhelical-ripple-trapped particles as well as those ofpassing particles. Snapshots of the contour lines ofthe distribution functions are shown in Fig. 1 (b)-(d).
The distribution function is advected along thecontour lines of the particle kinetic energy. Whilepassing particles elongate the profile, trapped
particles stay in the trapped regions. Thus,
fine-scale structures appear at the trapped-passing
boundary.
4. Linear ion-temperature-gradient instabilities
in helical plasmas
Employing the hybrid method of semi-(a)Lagrangian and additive semi-implicit RungeKutta schemes (SLASIRK), we have carried out
linear ITG simulations of a helical plasma.
Physical and numerical settings are the same as
those for the inward-shifted LHD case shown in
Ref. [1].
The linear growth rates and the real frequencies
are plotted as a function of the poloidal wave
number ky in Fig. 2. The results agree well withthe results obtained by using the fourth-order
Runge-Kutta-Gill scheme (RKG), while taking thetime step size larger than that of RKG. The
presented numerical method gives sufficiently
accurate results for t/ttr<0.4, which is comparable
to the transit time of passing particles through one
helical ripple.
References[1] T.-H. Watanabe, H. Sugama, S. Ferrando-Margalet,
Phys. Rev. Lett. 100, 195002 (2008).
[2] G. Strang, SIAM J. Numer. Anal. 5, 506 (1968).
[3] A. Staniforth, J. Côté, Mon. Wea. Rev. 119, 2206
(1991).
[4] S. Maeyama, et al., Plasma Fus. Res. 6, 2401028(2011).
[5] X. Zhong, J. Comp. Phys. 128, 19 (1996).

(b)Fig. 2. (a) Linear growth rates l and (b) real frequencies r as a function of the poloidal wave
number kyti. The solid, dashed and dotted lines represent the results obtained by RKG
with t/ttr = 0.005 and SLASIRK with t/ttr = 0.1, 0.4, respectively.