Numerical aspects of drift kinetic turbulence: ill-posedness, regularization and a priori estimates of sub-grid-scale terms

Turbulent magnetized plasmas are encountered in a wide variety of astrophysical situations (e.g. the solar corona, accretion disks, etc) and in magnetic fusion devices such as tokamaks. Moreover, these plasmas are weakly collisional, requiring a kinetic description. The challenges of simulating plasma turbulence in a six-dimensional (6D) phase space is formidable. In the presence of a strong magnetic guide field, the physical system can be reduced to five dimensions by averaging over the gyroradius of charged particles—a formalism termed gyrokinetic (See for a review [1, 2]). Recently, an idea from under-resolved hydrodynamic turbulence simulations was implemented in the context of gyrokinetic microturbulence simulations [3, 4]. The idea is that the large scales of turbulence are resolved on a given mesh, and that turbulent scales smaller than a certain filter width (generally related to the mesh spacing) are modeled. Such simulations are dubbed 'large eddy simulation(s) (LES)' in hydrodynamic turbulence and the physical models for the small unresolved scales are the 'sub-grid scale (SGS)' models. To develop SGS models for plasma turbulence (4D drift kinetic, 5D gyrokinetic and 6D kinetic) is one of our objectives. One rationale for developing SGS models is that under-resolved simulations using an Eulerian (or semi-Lagrangian) method to solve the kinetic equations on a velocity-and-physical-space mesh can incorporate the contribution of sub-grid terms involving nonlinear correlations. Toward achieving our final objective of developing SGS models, we first perform 'direct numerical simulations (DNS)' with a view to quantifying the SGS terms—such quantifications may lead to the development of phenomenological SGS models. As a first step toward these objectives, we have developed a 4D drift-kinetic turbulence code using an Eulerian formulation before embarking on a gyrokinetic code. Such a phased development approach was also employed by Grandgirard et al [5] (4D drift kinetic simulations) and Grandgirard et al [6] (5D gyrokinetic simulations). In this work, we focus on electrostatic 4D drift kinetic turbulence governed by the Vlasov–Poisson system of equations.

Our numerical approach to the 4D drift kinetic equations is Eulerian and the differencing method is high-order finite differences. The Eulerian approach has previously been used in conservative gyrokinetic Vlasov simulations by Idomura et al [7]. One question we are frequently confronted with in performing resolved simulations is the following. What is the physical cut-off mechanism, i.e. what is the smallest scale that needs to be resolved in a DNS of plasma turbulence? Without a physical collisional model, some form of regularization is required so that the small-scale generation does not proceed ad infinitum. Stated differently, one would like to terminate the process of small-scale generation and/or phase space filamentation of the distribution function by incorporating either collisions or some regularization to attain mesh-independent solutions.

The outline of the remainder of the paper is as follows. In section 2, we present the system of equations and the numerical method for solving these. In section 3, we present the results of simulations with a discussion of ill-posedness and regularization, as well as a comparison between different choices of the order of accuracy. In section 4, we present a formal derivation of SGS terms and quantify these in our drift kinetic simulations, and we conclude in section 5 with a summary and future directions.

The 4D drift kinetic system of equations presented in this section assumes a uniform magnetic field, $\boldsymbol{B}= B \skew3\hat{e}_z$. The geometry under consideration is cylindrical. For ions, the finite Larmor radius effects are neglected, which means that the ion trajectories are governed by the guiding center trajectories. The radial and azimuthal components of the drift velocity are given by

$$\begin{eqnarray} {\boldsymbol{v}}_{GC}\equiv{\boldsymbol{v}}_\perp&=& \frac {\boldsymbol{E} \times {\boldsymbol B}}{B^2}, \end{eqnarray} \noindent\tag{1}$$

where $\boldsymbol{E}$ is the electric field, and the parallel, i.e. the component of the guiding center velocity in the direction of the magnetic field, is denoted by v_∥. The evolution of the distribution function f≡f(r, θ,z, v_∥, t) is described by the drift kinetic Vlasov equation

$$\begin{equation} \frac{\partial f}{\partial t} + {\boldsymbol{v}}_\perp \cdot \nabla_\perp{f} + v_\parallel\frac{\partial f}{\partial z}+ \dot{v_\parallel}\frac{\partial f}{\partial v_\parallel}=0, \end{equation} \noindent\tag{2}$$

where $\nabla_\perp\equiv(\frac{\partial }{\partial r},\frac{1}{r}\frac{\partial }{\partial \theta})$ is the gradient operator in the perpendicular plane (i.e. the (r, θ) plane). Within the electrostatic approximation, the electric field is expressed as ${\boldsymbol{E}}=-\nabla \Phi$, where Φ is the electrostatic potential. The quasi-neutrality equation relates the electrostatic potential Φ(r, θ,z) to the first moment of the distribution function. The quasi-neutrality condition is expressed as a Poisson equation governing Φ,

$$\begin{equation} -\nabla_\perp\cdot\left[\frac{n_0(r)}{B\Omega_0}\nabla_\perp\right] \Phi +\frac{e n_0(r)}{T_\mathrm{e}(r)}\left(\Phi-\langle \Phi\rangle\right) =n_\mathrm{i} -n_0, \end{equation} \noindent\tag{3}$$

where the first term is the linearized polarization term, and the second one corresponds to the adibatic response of electrons [5], and . is the average in the z-direction. Here Ω₀=q_iB/m_i is the ion-cyclotron frequency, and n₀(r) and T_e(r) are, respectively, the ion density and electron temperature radial profiles, m_i is the ion mass and e is the electron charge. n_i≡n_i(r, θ,z) is the number density of the ions and is defined as the zeroth moment of the distribution function, i.e. n_i=∫_v∥fdv_∥.

The equilibrium distribution function is given by

$$\begin{eqnarray} &&f_{\mathrm{eq}}(r,v_\parallel) = \frac{n_0(r)}{2\pi T_\mathrm{i}(r)/m_\mathrm{i}}\exp\kern-2pt\left(-\frac{m_\mathrm{i}v_\parallel^2}{2T_\mathrm{i}(r)}\right), \end{eqnarray} \noindent\tag{4}$$

where T_i(r) is radial profile of the ion temperature. The above Vlasov–Poisson system (equations (2) and (3)) can be non-dimensionalized by normalizing the various quantities as follows:

$$\begin{eqnarray*} &&\mathrm{Temperature{:}}\ \skew3\hat{T}= \frac{T}{T_{\mathrm{e}0}}, \quad \frac{T_\mathrm{e}(r_0)}{T_{\mathrm{e}0}}=1,\\ &&\mathrm{Time{:}}\ \skew3\hat{t}= \Omega_0 t,\\ &&\mathrm{Velocity{:}}\ \skew3\hat{v}=\frac{v}{c_\mathrm{s}}, \quad c_\mathrm{s}=\sqrt{\frac{T_{\mathrm{e}0}}{m_\mathrm{i}}},\\ &&\mathrm{Length{:}}\ \skew3\hat{l}= \frac{\Omega_0}{c_\mathrm{s}} l, \\ &&\mathrm{Potential{:}}\ \skew3\hat{\Phi}=\frac{q_i}{T_{\mathrm{e}0}} \Phi,\\ &&\mathrm{Electric~field{:}}\ \skew3\hat{E}=\frac{1}{c_\mathrm{s} B} E, \end{eqnarray*}$$

where T_e0 is the electron temperature at a particular radial location r₀, and c_s is the sound speed. Dropping the $\hat{}$  over the various quantities, the non-dimensional form of equations (2)–(4) can be written as

$$\begin{eqnarray} &&\frac{\partial f}{\partial t} + {\boldsymbol{v}}_\perp \cdot \nabla_\perp{f} + v_\parallel\frac{\partial f}{\partial z}+ \dot{v_\parallel}\frac{\partial f}{\partial v_\parallel}=0, \end{eqnarray}\tag{5}$$

$$\begin{eqnarray} &&-\nabla_\perp\cdot\left[\frac{n_0(r)}{B}\nabla_\perp\right] \Phi +\frac {n_0(r)}{T_\mathrm{e}(r)}\left(\Phi-\langle \Phi\rangle\right) =n_\mathrm{i} -n_0, \end{eqnarray}\tag{6}$$

$$\begin{eqnarray} &&f_{\mathrm{eq}}(r,v_\parallel) = \frac{n_0(r)}{2\pi T_\mathrm{i}(r)}\exp\kern-2pt\left(-\frac{v_\parallel^2}{2T_\mathrm{i}(r)}\right). \end{eqnarray} \noindent\tag{7}$$

Note that the equation governing the time evolution of f has the same form as before, and that no new non-dimensional parameters have been introduced.

For the uniform magnetic field case, one may use the Liouville theorem, i.e.

$$\begin{eqnarray} &&\nabla_\perp\cdot {\boldsymbol{v}}_\perp +\frac{\partial v_\parallel}{\partial z}+\frac{\partial \dot{v_\parallel}}{\partial v_\parallel}=0, \end{eqnarray} \noindent\tag{8}$$

so the Vlasov equation can be written in conservative form as follows:

$$\begin{eqnarray} &&\frac{\partial f}{\partial t} + \nabla_\perp\cdot{{\boldsymbol{v}}_\perp f} + \frac{\partial v_\parallel f}{\partial z}+ \frac{\partial \dot{v_\parallel} f}{\partial v_\parallel}=0, \end{eqnarray}\tag{9}$$

$$\begin{eqnarray} &&\frac{\partial f}{\partial t} + \nabla\cdot{\boldsymbol{F}}=0. \end{eqnarray} \noindent\tag{10}$$

where $\boldsymbol{F}$ is defined as the flux of the distribution function and the operator ∇ is the divergence operator in 4D space.

In our Eulerian implementation, the 4D domain is discretized using finite volumes and all quantities (the distribution function, electrostatic potential, etc) are defined at the cell center of each finite volume. In equation (10), the fluxes are evaluated using higher-order upwind finite difference or a higher-order tuned centered finite difference (TCD). The derivative of flux F in the x-direction (where x may represent any of the radial, azimuthal, axial and parallel velocity directions and F the flux of the distribution function in x ) is written as

$$\begin{equation} \left(\frac{\partial F}{\partial x}\right)_i=\frac{ \skew3\hat F_{i+\frac{1}{2}}- \skew3\hat F_{i-\frac{1}{2}}}{h}, \end{equation} \noindent\tag{11}$$

where $\skew3\hat{F}$ is the numerical flux, h is the mesh spacing in the x-direction and $i\pm\frac{1}{2}$ indicates the boundary of the finite volume in the x-direction. Let c represent the wave velocity in the x-direction. Then the numerical flux is written as (for a seventh-order upwind scheme)

$$\begin{equation} \fl \skew3\hat F_{i+\frac{1}{2}}= \frac{1}{105}F_{i-2}-\frac{19}{210}F_{i-1}+\frac{107}{210}F_i +\frac{319}{210} F_{i+1}-\frac{101}{420}F_{i+2}+\frac{5}{84}F_{i+3}-\frac{1}{140}F_{i+4} \quad (c>0), \end{equation}\tag{12}$$

$$\begin{equation} \fl \skew3\hat F_{i+\frac{1}{2}}= \frac{1}{105}F_{i+3}-\frac{19}{210}F_{i+2}+\frac{107}{210}F_{i+1} +\frac{319}{210} F_{i}-\frac{101}{420}F_{i-1}+\frac{5}{84}F_{i-2}-\frac{1}{140}F_{i-3} \quad(c<0). \end{equation} \noindent\tag{13}$$

We omit details for the fifth- and third-order fluxes. The detailed coefficients of such high-order upwind biased schemes may be found in [8]. In addition, we have also implemented tuned central finite-difference schemes [9] for the derivatives in equation (10). For the fourth-order tuned central difference the numerical flux is given by

$$\begin{eqnarray} &&\fl \skew3\hat F_{i+\frac{1}{2}}= 0.0605F_{i-2}-0.264\,833\,3F_{i-1}+0.704\,333\,3F_i+0.704\,333\,3F_{i+1}-0.264\,833\,3F_{i+2}+0.0605F_{i+3}.\nonumber\\ \end{eqnarray}\tag{14}$$

A modified wave number of differencing schemes shows that the upwind-biased schemes have a measure of implicit dissipation, while central difference schemes have dispersion but no dissipation errors (not taking boundaries into account). In the absence of any physical cut-off mechanism in the governing equations, the implicit dissipation in upwind schemes provides a form of regularization. The modified wave number characteristics of tuned schemes and upwind schemes are given in [9] and [8], respectively. Unlike the semi-Lagrangian approach [5], there is no splitting along directions, and all derivatives in the flux divergence expression are computed at once. The time stepping scheme chosen is either a second- or a third-order TVD Runge–Kutta scheme.

The following modification to equation (10) is made for enhanced numerical stability. The distribution function is split into a perturbed and an equilibrium portion. This so-called δf formulation can be expressed as f'=f−f_eq. In terms of f' the Vlasov equation (equation (10)) becomes

$$\begin{equation} \frac{\partial f'}{\partial t} + \nabla_\perp\cdot{{\boldsymbol{v}}_\perp f'} + \frac{\partial v_\parallel f'}{\partial z}+ \frac{\partial \dot{v_\parallel} f'}{\partial v_\parallel}+S=0. \end{equation} \noindent\tag{15}$$

The source S is given in terms of the equilibrium distribution function, and the density and temperature profiles, and their derivatives as

$$\begin{equation} S=E_\theta f_{\mathrm{eq}}\frac{T_\mathrm{i}'(r)}{T_\mathrm{i}(r)}\left [ \frac{v_\parallel^2}{2 T_\mathrm{i}(r)^2} + \frac{n_0'(r)}{n_0(r)} - \frac{1}{2} \right] - E_z f_{\mathrm{eq}} \frac{v_\parallel}{T_\mathrm{i}(r)}. \end{equation} \noindent\tag{16}$$

A few preliminary calculations using the original form of equation (10) resulted in a large error stemming from taking finite differences in regions where the temperature (or density) profile exhibited large gradients. Using equation (16) ensures that no numerical derivative of the imposed temperature (or density) profile is undertaken, and this leads to an increased robustness of the code. Note that no linearization is made of the original Vlasov equation and that we are not ignoring any of the terms.

The Poisson equation (6) is solved by averaging over the axial direction, to first obtain the following:

$$\begin{eqnarray} &&-\nabla_\perp\cdot\left[\frac{n_0(r)}{B}\nabla_\perp\right] \langle\Phi\rangle =\langle n_\mathrm{i} -n_0\rangle. \end{eqnarray} \noindent\tag{17}$$

The above is a 2D equation in the radial and azimuthal directions. The correction, denoted by δΦ≡Φ−Φ, to the above average potential is obtained by subtracting equation (17) from (6), given by

$$\begin{equation} -\nabla_\perp\cdot\left[\frac{n_0(r)}{B}\nabla_\perp\right] \delta \Phi +\frac {n_0(r)}{T_\mathrm{e}(r)}\delta\Phi =n_\mathrm{i} -n_0 - \langle n_\mathrm{i}-n_0\rangle. \end{equation} \noindent\tag{18}$$

The linear operator on the left-hand side of equation (18) is not a function of z, whereas the right-hand side is indeed a function of z. This implies that independent 2D solves can be performed for each z-plane. In our simulations, no filtering of any kind of the distribution function or the electrostatic potential is employed for numerical robustness or stability reasons.

The boundary conditions are f'=0 in the radial direction, periodic in θ, zero gradient in z and f'=0 in the v_∥-direction. The initial conditions are similar to those used by Grandgirard et al [5] with the initial perturbation given by

$$\begin{equation} f'(r,\theta,z,v_\parallel,t=0)=f_{\mathrm{eq}}(r,v_\parallel)\ g(r) \ h(v_\parallel) \ \delta p(z,\theta), \end{equation}\tag{19}$$

where g(r) and h(v_∥) are exponential functions, centered at the mid-radial position and v_∥=0, respectively, and decay to zero well before the boundaries in r and v_∥. The perturbation δp(z, θ) is

$$\begin{equation} \delta p(z,\theta)=\sum_{m,n} \epsilon_{mn} \cos\left(\frac{2\pi n}{L_z}z + m\theta + \phi_{mn}\right), \end{equation}\tag{20}$$

where _mn is a small random amplitude (10⁻⁴) and _mn is a randomized phase.

For the electrostatic potential, Φ, periodic boundary conditions are enforced in the axial and azimuthal directions. It is convenient to expand both the average axial potential and the correction δΦ in Fourier modes. We express the radial derivatives in the linear variable coefficient Poisson operator in equations (17) and (18) using fourth-order central finite differences. The Fourier modes are then governed by a penta-diagonal linear equation which can be efficiently solved. We choose the radial boundary condition for Φ to be homogeneous Neumann, which implies no poloidal rotation on both radial boundaries. Grandgirard et al [5] also suggest that homogeneous Neumann in the radial direction is the correct boundary condition for Φ, but they actually implement homogeneous Dirichlet boundary conditions instead.

Results of our simulations of ITG drift kinetic turbulence (see section 3) show that the nonlinear regime is characterized by the rapid generation of small scales of the distribution function. Furthermore, there are clear signatures of small-scale generation in the energy spectra (also presented in section 3). A natural question that arises is then: what is the physical cut-off mechanism to terminate this cascade to fine scales? Consider the following transformations:

$$\begin{equation} \begin{array}{l} x\rightarrow\alpha x', \quad t\rightarrow\alpha t', \quad B\rightarrow B'/\alpha,\quad \Phi\rightarrow\Phi'/\alpha^2,\\ E\rightarrow E'/\alpha, \quad f\rightarrow f', \quad n\rightarrow n', \quad v_\parallel\rightarrow v'_\parallel, \end{array} \end{equation}\tag{21}$$

for arbitrary nonzero alpha. Under this transformation, the Vlasov–Poisson system of equations is unchanged, indicative of the fact that there is no physical cut-off mechanism in these equations. Such a system of equations will lead to numerical solutions that are non-convergent upon mesh refinement. Mathematically, the system of equations is then considered to be not well-posed. Stated differently, the distribution function will exhibit smaller and smaller features when the mesh size is increased. This phenomenon has been observed in the form of filamentation of the distribution function in phase space (see [10] and references therein). The most physically relevant manner in restoring well-posedness to the Vlasov–Poisson system is to include models for collisions. Schekochihin et al [11] point out that the termination of the net forward energy cascade to fine scales, and the satisfaction of the Boltzmann H-theorem, is by collisions, no matter how small their contribution. Collisions provide the only physically relevant mechanism to finally dissipate the kinetic energy to heat. Within the context of gyrokinetic simulations, Garbet et al [2] provide a good discussion of the need for physically relevant collision models or numerical dissipation mechanisms, because without these the gyrokinetic equation also leads to the continuous production of small-scale structures by linear and nonlinear mixing effects. In the present work on drift kinetic simulations, we do not implement a physically correct collision model. Instead, we rely on a few different regularization mechanisms (enumerated below) to provide a cut-off for the simulations. Similar ad-hoc regularizations have been used, for example, in recent works on gyrokinetic simulations [3] and those on LES of gyrokinetic turbulence. The regularization methods employed in the simulations are as follows.

(i)	Physically relevant collision model: We have extended the 1D collision model of Rathmann and Denavit [12]. We dub this model as 'physically relevant' because it is motivated by an actual collision model in 1D and in a sense is different from other regularizations (e.g. hyperviscosity) discussed below, but we do not expect this model to be exactly physically correct. The Vlasov equation (equation (5)) is then modified as $$\begin{equation} \frac{\partial f}{\partial t} + {\boldsymbol{v}}_\perp \cdot \nabla_\perp{f} + v_\parallel\frac{\partial f}{\partial z}+ \dot{v_\parallel}\frac{\partial f}{\partial v_\parallel}=\left(\frac{\partial f}{\partial t}\right)_{\mathrm{coll}}, \end{equation}\tag{22}$$ where the collision term is given by $$\begin{eqnarray} \left(\frac{\partial f}{\partial t}\right)_{\mathrm{coll}}&=&\frac{\partial }{\partial v_\parallel}\left[\gamma v_\parallel f + D \frac{\partial }{\partial v_\parallel}\left (\gamma f\right)\right],\nonumber\\ &=&\frac{\partial }{\partial v_\parallel}\left (\gamma v_\parallel f\right) + D\frac{\partial ^2}{\partial v_\parallel^2}\left(\gamma \right), \end{eqnarray}\tag{23}$$ consisting of a friction-like term and a diffusion term. The coefficients γ and D are given by $$\begin{equation} \gamma=\frac{ \left (3 \langle v_\parallel^2\rangle \right )^{\frac{3}{2}}\gamma_0}{ \left (2 \langle v_\parallel^2\rangle+v_\parallel^2 \right )^{\frac{3}{2}}}, \quad D=\frac{\langle \gamma v_\parallel^2\rangle}{\langle \gamma\rangle}, \end{equation}\tag{24}$$ where $$\begin{eqnarray*} &&\langle \cdot \rangle = \int \cdot f \mathrm{d}V \mathrm{d}v_\parallel. \end{eqnarray*}$$ Numerically, the derivative terms in the collision operator are computed using standard fourth-order centered finite differences.
(ii)	Hyperviscosity: A fourth-order hyperviscosity term is added to the right-hand side of the Vlasov equation. In our simulations, the hyperviscosity term has a constant coefficient. The choice of this term is common in both hydrodynamic turbulence [13] and gyrokinetic turbulence simulations [3].
(iii)	Artificial viscosity: A second-order Laplacian term is added to the right-hand side of the Vlasov equation. This practice is somewhat similar to shock hydrodynamics simulations with artificial viscosity. In our simulations, the artificial viscosity term is ν∇2f, where the viscosity coefficient is computed using local gradient information as ν=ν0|∇f|/∥∇f∥∞, with ∥∇f∥∞ being the max norm of the gradient of the distribution function.
(iv)	Implicit dissipation: Upwind schemes have a measure of implicit dissipation which is a form of regularization.

In this section we present the results on ion temperature gradient (ITG) drift kinetic turbulence as computed by the Eulerian method discussed above (see Krommes [14] for a review of plasma microturbulence and ITG turbulence, and Grandgirard et al [5] on drift-kinetic ITG turbulence). Before reporting on the main simulation results, we briefly digress to mention a linear benchmark case. We performed a simulation with a single mode perturbation (m,n)=(11, 3) similar to that discussed in [5] (see section 4.2.2 of [5]). The growth rate quantified in our code was γ=6.2184×10⁻³, which compares well with Grandgirard et al [5], and the analytical growth rate of 6.259×10⁻³.

For the simulation results presented here, the domain of investigation is [0:15]×[0:2π]×[0:1500]×[−8:8] in our (r, θ,z, v_∥) space. The density and temperature profiles are given in figure 1, and these are fixed for all the results presented. The first set of results is for a mesh size of 256×256×32×256, and fluxes are computed with a fourth-order TCD. Moreover, for these results, the non-dimensional time step is fixed at dt=0.1. The initial perturbations δp (equation (20)) consisted of eight azimuthal and eight axial modes with random amplitudes and phases.

Zoom In Zoom Out Reset image size

Figure 2 shows a time sequence of the distribution function f visualized at v_∥=0 and presented in physical space on constant coordinate surfaces r=2, θ=π,z=0. The ITG drives the instability and the perturbations grow linearly until about t=1500, following which we observe a nonlinear saturation, zonal flow and 'mixing' in phase space, as well as the generation of small-scale structures in the distribution function. One diagnostic computed in the simulations is the following moment of the distribution function:

$$\begin{eqnarray} &&Q(t) = \frac{1}{2}\int f ~v_\parallel^2 ~r ~v_{\mathrm{GC},r}~\mathrm{d}v_\parallel~\frac{\mathrm{d}\theta}{2\pi}~\frac{\mathrm{d}z}{L_z}~\frac{\mathrm{d}r}{L_r}, \end{eqnarray}\tag{25}$$

where v_{GC, r} is the radial component of the guiding center drift velocity, and L_r, L_z are the domain lengths in the r, z directions, respectively. This diagnostic may be considered as a measure of the average heat flux in the radial direction. Q(t) for the present simulation is shown in figure 3(a). An examination of the electrostatic potential Φ (figure 4) shows the early linear development of the mode structure followed by the nonlinear phase during which Φ evolves into two dominant bands of opposite sign.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For this simulation, the code conserves mass and energy to within 0.001 and 0.02%, respectively. The relative energy conservation is plotted in figure 5(a). Defining the total kinetic (K) and potential (P) energies as

$$\begin{equation} \fl K = \frac{1}{2}\int (f-f_{\mathrm{eq}})~ v_\parallel^2 ~r \mathrm{d}r~\mathrm{d}\theta~\mathrm{d}z~\mathrm{d}v_\parallel,\quad P = \frac{1}{2}\int \phi~ (n_\mathrm{i}-n_0) ~r \,\mathrm{d}r~\mathrm{d}\theta~ \mathrm{d}z, \end{equation}\tag{26}$$

one may derive (see [5]) the total energy conservation as K+P=constant, and moreover we obtain the following:

$$\begin{equation} \frac{\mathrm{d}K}{\mathrm{d}t}+\frac{\mathrm{d}P}{\mathrm{d}t}=0. \end{equation}\tag{27}$$

As stated above, there is net exchange between kinetic and potential energies while maintaining the total energy constant. In figure 5(b), we observe that the code accurately preserves the energy conversion rate. The potential and kinetic energy spectra at various times are plotted in figure 6. The spectra continue to show an increase in energy at small scales during the early part of the simulation. We observe little difference between times t=3000 and t=4000, indicating that we have attained a statistical steady state.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The above simulation results were obtained with an artificial viscosity Laplacian term in the Vlasov equation. The artificial viscosity coefficient was computed as a local function of the gradient of the distribution function and we choose ν₀=10⁻⁴. Similar results were obtained with the hyperviscosity term. Both the artificial viscosity and hyperviscosity regularization terms did not have any measurable effect until sufficiently late into the nonlinear regime. The same simulation was repeated with the physically relevant collision model, and in this case the onset of the nonlinear regime was delayed for γ₀=10⁻³, while for γ₀=10⁻⁴, the results were similar to the Laplacian artificial viscosity case at least until about t=3200. The γ₀=10⁻⁴ still blows up, albeit later than t=3500, indicating that the level of collisions was not sufficient, while the γ₀=10⁻³ case does not blow up at all. These comparisons are shown in figure 3(b) for a 128×128×32×128 mesh. Also shown in the same figure is the case with no regularization, for which the simulation exhibits numerical instability shortly after t=3000, because there was no physical cut-off and the very small scales are not adequately resolved.

We now compare the differences between the fourth-order TCD (central differences) with the seventh-order upwind scheme. The distribution function at times t=2000 and t=4000 is shown in figure 7. Comparing these snapshots with the previous time sequence in figures 2(c) and (d), we observe very little difference at the early nonlinear stage (t=2000), while the difference between the central difference and upwind methods is larger at late times (t=4000). The upwind method shows a smoother distribution function although the gross features (or large-scale structures) are quite similar between the two methods. In figure 8, we compare the regularized Vlasov equation computed with a fourth-order TCD method against the seventh-order upwind method applied to the Vlasov equation with no explicit regularization. The difference in Q(t) is less than 1% until t=2300. Even in the late nonlinear stage of the simulation the difference between these two methods appears to be quite small. The kinetic energy spectrum at various times (shell averaged kinetic energy spectrum after computing the kinetic energy in physical space) is shown in figure 8(b). The central difference and upwind methods agree well up to k=70. The upwind method has a steeper spectrum for high wavenumbers, whereas the central scheme shows some evidence for energy pileup at high wavenumbers. The seventh-order upwind method has implicit numerical dissipation and that is sufficient to regularize the Vlasov–Poisson system. Simulations at different mesh sizes (r, θ, v_∥): 256, 128, 64, 32 were undertaken fixing the z-direction mesh at 32 with the seventh-order upwind method. The effect of resolution on the upwind scheme is shown in figure 9(a). It is evident that the upwind scheme will result in a more robust method when under-resolved simulations are performed while one may conjecture that under-resolved computations with central schemes will probably not be very robust. In figure 9(b), we compare Q(t) computed with third-, fifth- and seventh-order upwind methods as well as the fourth-order TCD method on a 128×128×32×128 mesh. All methods reproduce the same result until about t=2000, i.e. the linear phase is very well captured by all the methods. Differences between the fifth- and seventh-order upwind methods are the least throughout the duration of the simulation. The maximum radial electric field is plotted in figure 9(c). It is evident that on this log scale plot the differences in the peak radial electric field are small (except for the mesh with 32 points in all directions), and the growth rate (numerical value of 0.0036) during the linear phase is comparable between all the methods shown.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In hydrodynamic turbulence, the idea of a kinetic energy cascade is generally considered well established and follows the Kolmogorov (K41) theory. Filtering the incompressible Navier–Stokes equation leads to nonlinear correlations, dubbed the SGS terms which must be modeled in terms of the resolved scales. Simulations which do not resolve the flow scales all the way to the dissipation scale are known as LES. For a discussion of filtering and LES, see [15–17]. A similar filtering process may be applied to the drift kinetic Vlasov equation. Generalizing the filtering introduced by Leonard [18] to include the velocity space, we may filter the distribution function as follows:

$$\begin{eqnarray} &&\skew4\bar f({\boldsymbol{x}},{\boldsymbol{v}},t)=\int G({\boldsymbol{r}},{\boldsymbol{w}}) ~f({\boldsymbol{x}}-{\boldsymbol{r}},{\boldsymbol{v}}-{\boldsymbol{w}},t)~\mathrm{d}{\boldsymbol{r}}~\mathrm{d}{\boldsymbol{w}}~\equiv~ G*f, \end{eqnarray}\tag{28}$$

where the filter obeys the normalization condition

$$\begin{eqnarray} &&\int G({\boldsymbol{r}},{\boldsymbol{w}})~\mathrm{d}{\boldsymbol{r}}~ \mathrm{d}{\boldsymbol{w}}=1. \end{eqnarray}\tag{29}$$

The filtered velocity of the guiding centers is defined using the same filter as $\bar {\boldsymbol{v}} = G*{\boldsymbol{v}}$. The parallel velocity coordinate is an independent variable and hence $\bar v_\parallel=v_\parallel$. The filtered Vlasov equation is

$$\begin{equation} \frac{\partial \skew4\bar f}{\partial t} + \nabla_\perp\cdot{\overline {{\boldsymbol{v}}_\perp} \skew4\bar f} + \frac{\partial v_\parallel \skew4\bar f}{\partial z}+ \frac{\partial \overline{\dot{v_\parallel}} \skew4\bar f}{\partial v_\parallel} + \nabla\cdot {\cal S}=0. \end{equation}\tag{30}$$

where ${\cal S}$ is the SGS vector. The divergence of this SGS vector is given by

$$\begin{eqnarray} &&\fl \nabla\cdot {\cal S}=\frac{1}{r}\frac{\partial }{\partial r}\left[r\left(\overline{v_{\mathrm{GC},r} f}-\overline{v_{\mathrm{GC},r}}{\skew4\bar f}\right)\right]~+~\frac{1}{r}\frac{\partial }{\partial \theta}\left[\overline{v_{\mathrm{GC},\theta} f}-\overline{v_{\mathrm{GC},\theta}}{\skew4\bar f}\right] + \frac{\partial }{\partial v_\parallel}\left[\overline{\dot v_\parallel f}-\overline{\dot v_\parallel}{\skew4\bar f}\right]. \end{eqnarray}\tag{31}$$

Consider the radial component of ${\cal S}$ below, where we have decomposed every field into a filtered (with an overbar) and a residual field (denoted with a prime),

$$\begin{eqnarray} &&\fl {\cal S}_r=\overline{v_{\mathrm{GC},r} f}-\overline{v_{\mathrm{GC},r}}{\skew4\bar f}\nonumber\\ &&=\underbrace{\overline{\overline {v_{\mathrm{GC},r}} \overline{f}}-\overline{\overline {v_{\mathrm{GC},r}}}\overline{ \overline{f}}}_{\mathrm{Leonard}}\,+\,\underbrace{ \overline{\overline {v_{\mathrm{GC},r}} {f'}}-\overline{\overline {v_{\mathrm{GC},r}}}\overline{ {f'}} \,+\,\overline{{v_{\mathrm{GC},r}'} \overline{f}}-\overline{{v_{\mathrm{GC},r}'}}\overline{ \overline{f}} }_{\mathrm{Cross}}\,+\, \underbrace{ \overline{{v_{\mathrm{GC},r}'} {f'}}-\overline{{v_{\mathrm{GC},r}'}}~\overline{ {f'}} }_{\mathrm{Reynolds}},\nonumber\\ \end{eqnarray}\tag{32}$$

where the terms labeled Leonard, Cross and Reynolds are analogous to the Leonard stress, Cross stress and SGS Reynolds stresses in hydrodynamic turbulence, respectively. Other terms with v_{GC, r} replaced above by v_{GC, θ} or $\dot v_\parallel$ can be similarly written.

In this work, we do not attempt to quantify each of these terms separately, instead focusing on the overall SGS term ${\cal S}$. Borrowing the terminology from hydrodynamic turbulence, we distinguish between a posteriori tests and a priori tests [19] of SGS models. A posteriori testing refers to simulating the turbulence with an SGS model and testing against DNS, turbulence theory or experiments. On the other hand, a priori testing refers to performing DNS, then filtering the data to extract SGS terms and comparing these with any proposed SGS models. In the spirit of a priori testing, we quantify SGS terms from DNS results. Note, however, that we have not proposed an SGS model, and leave that for future work. We apply a sharp filter in Fourier space to the DNS data and quantify each of the SGS terms, i.e. the radial, azimuthal and axial components of the SGS vector ${\cal S}$. The filtering is at present done only in physical space with the cut-off wavenumber set at 16. In figure 10, we plot the time history of the root mean square (rms) of each of the SGS terms by filtering the DNS results obtained from the fourth-order TCD and seventh-order upwind methods. It is evident that for drift kinetic turbulence the rms axial component of the SGS vector is more than two orders of magnitude smaller than the radial and azimuthal components, both of the latter being nearly the same. This implies that modeling efforts may be concentrated on the radial and azimuthal terms. Another observation is that the SGS terms during the linear phase are quite negligible but grow by several orders of magnitude by the onset of the nonlinear phase and then slowly change during the nonlinear saturation phase of the simulation. After t=2500 the magnitude of the largest SGS terms compared with $\overline{v_{\mathrm{GC},r}f}$ and $\overline{v_{\mathrm{GC},\theta} f}$ ranges between 20 and 30%. Hence, under-resolved computations that do not model or take SGS terms into account are missing a significant contribution.

Zoom In Zoom Out Reset image size

A new 4D drift kinetic plasma turbulence code has been developed using an Eulerian formulation to solve the Vlasov–Poisson system of equations in cylindrical geometry. High-order finite differences are used to compute the fluxes, whereas the Poisson equation is solved with a combination of Fourier decomposition and fourth-order finite differences in the radial direction. The system of equations must be regularized to provide a cut-off mechanism. In this paper, we used simple regularization methods such as the explicit use of artificial viscosity, hyperviscosity, a simple collision model or implicit regularization provided by upwind methods. The regularized system is robust and stable, with the upwind method resulting in the most robust code. A comparison between different numerical methods was made. All methods (order of accuracy and regularization method chosen) showed very little difference during the linear growth phases with the differences increasing as the nonlinear saturation sets in, although, in a statistical sense, the differences even in the nonlinear phase were not very significant. We presented the formalism for filtering the Vlasov equation and presented a priori quantifications of the SGS terms for the 4D drift-kinetic turbulence simulations. These SGS terms can be a significant portion during the nonlinear phase and must not be ignored in under-resolved simulations. Extension of the present numerical method, the filtering formalism, quantifications of SGS terms applicable to 5D gyrokinetic simulations and developing SGS models are some of the future directions.

RS was supported for this work by the King Abdullah University of Science and Technology. Simulations in this paper were conducted using Shaheen, IBM Blue Gene P at KAUST.