Time diffusion method for gyrokinetic simulation of electrostatic turbulence with kinetic electrons

Abstract

A new method based on a time diffusion technique is developed to mitigate the numerical instability induced by the high frequency electrostatic shear Alfvén wave (${\omega }_H$ mode) in the gyrokinetic simulation of electrostatic turbulence. It can preserve fully drift-kinetic electron effects and can be incorporated efficiently with explicit time integration scheme. This method has been implemented in the semi-Lagrangian gyrokinetic code NLT (Lei Ye, et al. (2016) [25]) and applied successfully to linear and nonlinear ITG/TEM turbulence simulation. Numerical results show that this time diffusion method can greatly enlarge time step size and improve numerical stability for electrostatic turbulence simulation in tokamaks.

Introduction

Large scale gyrokinetic simulation provides an indispensable alternative to study the micro-turbulence in magnetized plasma, which plays a dominant role in the anomalous transport of particle and heat in tokamaks. The drift wave turbulence in a tokamak usually has a characteristic frequency spectrum on the order of diamagnetic drift frequency, which is much lower than the ion cyclotron frequency. Based on this fact, the fast gyro-motion of a particle can be decoupled from the slow drift motion of the gyro-center in the gyrokinetic theory, which enables ones to focus on the low frequency phenomena in fusion plasma while keeping the small scale finite-Larmor-radius effects [1]. At the same time, the 6D particle phase space is reduced to 5D gyrocenter phase space in gyrokinetic theory by averaging out of the gyro-angle. This reduction can save a large amount of computational cost for numerical simulation. In recent years, quite a number of simulation codes have been developed based on gyrokinetic theoretical model and applied to investigate the properties of turbulence transport in tokamaks [2].

Usually the gyrokinetic model is employed to describe ions in a gyrokinetic code, while the physical models for electrons can be different. The adiabatic electron model is the simplest one for modeling electron which excludes all the kinetic effects of electrons. With this adiabatic electron approximation the typical time step size of gyrokinetic simulation is on the order of ion cyclotron period. However, if kinetic electrons are completely included, the limitation on the time step will become stringent. Basically speaking, this is due to the large mass ratio of $m_i/m_e$ and the fast steaming of the passing electron along the field line. Here $m_i$ and $m_e$ are the masses of ion and electron, respectively. On the one hand, the electron Courant-Friedrichs-Lewy (CFL) condition, $k_{\parallel }{v}_{te}\mathrm{\Delta }t\lesssim 1$, should be satisfied for passing electrons because of the accuracy requirement [3]. Here $k_{\parallel }$ is the wave vector parallel to the equilibrium magnetic field, $v_{te}=\sqrt{T_e/m_e}$ is the thermal velocity of the electron. On the other hand, the shear Alfvén wave will reduce to the electrostatic shear Alfvén wave (${\omega }_H$ mode) in the electrostatic limit or the low β limit [3], [4]. Here $\beta ={\mu }_{0}nT/B^2$ is the ratio between plasma pressure and magnetic pressure. This ${\omega }_H$ mode mainly comes from the response of the ion polarization density to the kinetic passing electron and has a frequency of ${\omega }_H=(k_{\parallel }/k_{\perp })\sqrt{m_i/m_e}{\mathrm{\Omega }}_i$. In the long wavelength situations, i.e. $k_{\parallel }/k_{\perp }\ll \sqrt{m_i/m_e}$, ${\omega }_H$ is comparable to or even greater than the ion gyro-frequency ${\mathrm{\Omega }}_i$. Here $k_{\perp }$ is the wave vector perpendicular to the equilibrium magnetic field. It is believed that such high frequency ${\omega }_H$ mode has little influence on the electrostatic turbulence. However, it does impose severe limitation, ${\omega }_H\mathrm{\Delta }t\lesssim 1$, on the time step of electrostatic (or low β) gyrokinetic simulation, especially with explicit integration scheme. Besides, it also can easily cause the numerical instability originates from the aliasing effects on a temporal grid [3].

There are generally two approaches for circumventing the time step constraint of simulation imposed by a high frequency oscillation, physically insignificant wave [5]. One way is to replace the original equation with a set of approximate equations or models that do not support the high frequency wave like ${\omega }_H$ mode. For example, the kinetic trapped electron model or the hybrid electron model [6], [7], [8], [9], [10], [11] has been developed to preserve part of the kinetic electron effects in gyrokinetic simulation. In these models, passing electrons are treated with the adiabatic approximation to avoid numerical problems mentioned above, while fully drift-kinetic effects of trapped electrons are preserved. Thus, the trapped electron mode (TEM) as well as the ion temperature gradient (ITG) turbulence, two of the most important candidates for anomalous transport, can be successfully simulated with these models. However, it has been pointed out that in some conditions kinetic passing electron can have significant influence on the turbulence saturation level [12], [10].

Another method relies on numerical techniques to avoid the numerical instability induced by the high frequency oscillation, such as the implicit scheme [13], [14], [15], [16], [17] or the split-weight scheme [18], [4], [19], [20], [21] for electron equations. The implicit scheme is generally unconditionally stable and allows for larger time step. The split-weight scheme can be used to suppress the unwanted high frequency oscillations and relax the CFL condition in slab geometry [18], [4] and tokamak [20], [21]. A generalization of the split-weight scheme has been developed for electromagnetic (EM) turbulence simulation [19]. Recently, a new pull-back mitigation method based on the mixed variables has been developed to mitigate the cancellation problem in EM simulation [22], [23]. It has been shown that this method is essentially equivalent to an implicit scheme so that can enlarge simulation step significantly [24].

As we know, the computational efficiency has always been a critical concern for gyrokinetic simulation even with state-of-the-art supercomputer. Therefore, it is still an open question to perform turbulence simulation with an acceptable time step in electrostatic or low-β limit by suppressing the numerical instability caused by ${\omega }_H$ mode, while retaining the kinetic effect of passing electrons. In this work, an explicit time diffusion method is developed to mitigate the numerical problem caused by the ${\omega }_H$ mode and relax the CFL condition for the gyrokinetic simulation in the electrostatic limit. This method has been successfully implemented into the semi-Lagrange gyrokinetic code NLT [25] to include fully drift-kinetic effects of the electron. By taking this new approach, NLT can be used to simulate the electrostatic turbulence transport in tokamak, including the ITG and TEM turbulence.

The remaining part of this paper is organized as follows. In Sec. 2, the time diffusion method and its effect on the dispersion relation of ${\omega }_H$ mode are introduced. In Sec. 3, the numerical implementation of the time diffusion method in NLT code is described. Benchmark tests for linear ITG, TEM and zonal flows are presented in Sec. 4. Nonlinear simulation results of electrostatic turbulence with kinetic electrons are shown in Sec. 5. Finally, a brief summary is given in Sec. 6.