Field-aligned interpolation for semi-Lagrangian gyrokinetic simulations (Q1743433)

The authors study the numerical solution of the gyrokinetic Vlasov equation by means of a semi-Lagrangian method using field aligned interpolation. The equation under consideration is \[ \Big (\frac{\partial}{\partial t} +\mathbf{u}(t,\mathbf{x},v_{\parallel},\mu) \cdot \nabla + a_{\parallel}(t,\mathbf{x},v_{\parallel},\mu) \frac{\partial}{\partial v_{\parallel}} \Big )f(t,\mathbf{x},v_{\parallel},\mu) = 0.\tag{1} \] Here, \(f\) is the gyro-center distribution function for a single ion species which is immersed in a static magnetic field \(\mathbf{B}(\mathbf{x})\). The phase space coordinates are \((\mathbf{x},v_{\parallel},\mu)\) in which \(\mathbf{x}\) is the gyro-center position, \(v_{\parallel}\) is velocity parallel to \(\mathbf{B}\), and \(\mu\) is the modified magnetic moment. The fields \(\mathbf{u}\) and \(a_{\parallel}\) are given in terms of \(\mathbf{B}\) and the gyro-center Hamiltonian \(H(t,\mathbf{x},v_{\parallel},\mu) = \frac12 mv_{\parallel}^2 + \mu \mathbf{B}(\mathbf{x}) + q\langle \phi \rangle_{\alpha} (t,\mathbf{x})\) in which \(\phi(t,\mathbf{x})\) is the electrostatic potential, and \(\langle \cdot \rangle_{\alpha}\) the gyro-average operator. The function \(\phi(t,\mathbf{x})\) is obtained as the solution to \[ -\Big (\nabla_{\perp} \cdot \frac{\rho_{th,i}^2}{\lambda_{D,i}^2} \nabla_{\perp} \Big ) \phi + \frac{1}{\lambda_{D,e}^2}(\phi - \langle \phi \rangle_f) = \frac{q_i}{\epsilon_0}(n_i - n_0). \tag{2} \] In (2), \(\nabla_{\perp}\) is a gradient operator perpendicular to the magnetic field, \(\langle \cdot \rangle_f\) an averaging operator, \(\rho_{th,i}(\mathbf{x})\) is the thermal ion Larmor radius, \(\lambda_{D,i}(\mathbf{x}),\; \lambda_{D,e} (\mathbf{x})\) Debye lengths for ions and electrons, \(n_i(t,\mathbf{x})\) is the instantaneous ion number density, \(n_0(\mathbf{x})\) is the equilibrium ion density, \(q_i\) the ion charge, and \(\epsilon_0\) the vacuum permittivity. In the paper, the authors consider two model problems. The first is a simplified gyrokinetic simulation in a cylindrical geometry, the second is a solution to (1) in a torus. Referring to the toroidal geometry, the equation (1) is solved in terms of variables \((r, \theta, \varphi)\) in which \(\theta\) is a poloidal angle, \(\varphi\) is a toroidal angle, and the radius \(r\) is defined through \(R = R_0 + r\cos(\theta)\) where \(R_0\) is the major radius of the torus. The numerical method employed to solve (1) is a semi-Lagrangian method in which characteristic trajectories are computed backward in time from grid points to points off the grid. Function values at the points off the grid at time \(t\) are interpolated from points on the grid. These values are then used to update function values on the grid at time \(t + \Delta t\). A field aligned interpolation is used to get the function values off the grid and also to get derivatives with respect to \(\varphi\) of the potential function \(\phi\). Such derivatives of \(\phi\) are needed for the problem in the torus. The motivation for using field aligned interpolation is that there is a relatively small and slow variation of relevant function values in the direction parallel to the magnetic field lines whereas variations in the perpendicular direction are large. By interpolating from grid points that align closely to the path of a magnetic field line, one can get a high accuracy in interpolated values using relatively fewer grid points. This, therefore, improves the efficiency of the numerical method. In the toroidal geometry, this 1-D interpolation along field lines is combined with a 2-D interpolation on the poloidal plane at points of intersection with the field lines. There are two components to this paper, one analytical, the other numerical. For the analytical part, the authors consider a simplified 2-D linear, constant coefficient advection equation in the variables \((t,\theta,\varphi)\). Using semi-discrete Fourier transform methods, it is proved that the backward, field-aligned semi-Lagrangian method is unconditionally stable as applied to the simplified advection equation. This is followed by a convergence analysis which assesses the error in the 2-D field-aligned method and compares it to the error of the classical non-aligned scheme. The last part of the paper describes numerical results. The numerical study compares results using the field-aligned semi-Lagrangian method with standard methods using non-aligned interpolation schemes. It is concluded that the field aligned method can achieve comparable accuracy to non-aligned methods but with considerable less computer memory. With regard to the toroidal geometry, the memory reduction is accomplished by using fewer grid points in the \(\varphi\) direction. Also, there is a moderate (but improvable) speed-up in simulation time compared with non-aligned methods.