Verification of Gyrokinetic codes: theoretical background and applications

TL;DR

This paper introduces a comprehensive theoretical and numerical framework for verifying gyrokinetic codes, crucial for accurate plasma turbulence modeling in fusion research, by combining variational formulations and benchmark tests.

Contribution

It presents a novel unified approach using variational methods for theoretical verification and benchmark testing for numerical validation of gyrokinetic codes.

Findings

Established a hierarchy of reduced GK Vlasov-Maxwell equations

Verified two major GK codes: ORB5 and GENE

Applied methodology to CYCLONE Base Case scenarios

Abstract

In fusion plasmas the strong magnetic field allows the fast gyro-motion to be systematically removed from the description of the dynamics, resulting in a considerable model simplification and gain of computational time. Nowadays, the gyrokinetic (GK) codes play a major role in the understanding of the development and the saturation of turbulence and in the prediction of the subsequent transport. Naturally, these codes require thorough verification and validation. Here we present a new and generic theoretical framework and specific numerical applications to test the faithfulness of the implemented models to theory and to verify the domain of applicability of existing GK codes. For a sound verification process, the underlying theoretical GK model and the numerical scheme must be considered at the same time, which has rarely been done and therefore makes this approach pioneering. At the analytical level, the main novelty consists in using advanced mathematical tools such as variational formulation of dynamics for systematization of basic GK code's equations to access the limits of their applicability. The verification of numerical scheme is proposed via the benchmark effort. In this work, specific examples of code verification are presented for two GK codes: the multi-species electromagnetic ORB5 (PIC) and the radially global version of GENE (Eulerian). The proposed methodology can be applied to any existing GK code. We establish a hierarchy of reduced GK Vlasov-Maxwell equations implemented in the ORB5 and GENE codes using the Lagrangian variational formulation. At the computational level, detailed verifications of global electromagnetic test cases developed from the CYCLONE Base Case are considered, including a parametric $\beta$-scan covering the transition from ITG to KBM and the spectral properties at the nominal $\beta$ value.

Tables

Table 1. Table 1: Reduced Hamiltonian dynamics comparison.

ORB5 (Hamiltonian, low $\beta$)	GENE (Modified parallel symplectic)
	local: $\alpha =1$; global $\alpha =0$
${𝐀}^{\ast }=𝐀+\left(c/e\right)p_z\widehat{𝐛}$	${𝐀}^{\ast \ast }=𝐀+\left(c/e\right)p_{\parallel }\widehat{𝐛}+{ϵ}_{\delta }⟨A_{1\parallel }\left(𝐗+{𝝆}_0\right)⟩\widehat{𝐛}$
$H_0^{\mathrm{ORB5}}=p_z^2/\left(2m\right)+\mu B$	$H_0^{\mathrm{GENE}}=p_{\parallel }^2/\left(2m\right)+\mu B$
$H_1^{\mathrm{ORB5}}=e⟨{\psi }_1^{\mathrm{ORB5}}⟩$,	$H_1^{\mathrm{GENE}}=e⟨{\psi }_1^{\mathrm{GENE}}⟩$,
with ${\psi }_1^{\mathrm{ORB5}}={\varphi }_1-A_{1\parallel }{p}_z/m$	with ${\psi }_1^{\mathrm{GENE}}={\varphi }_1-\alpha {𝐀}_{1⟂}\cdot {𝐯}_{⟂}$
$H_2^{\mathrm{ORB5}}=-\frac{e^2}{2B}\frac{\partial }{\partial \mu }⟨{\tilde{\varphi }}_1{\left(𝐗+{𝝆}_0\right)}^2⟩$	$H_2^{\mathrm{GENE}}=-\frac{e^2}{2B}\frac{\partial }{\partial \mu }⟨{\tilde{\psi }}_1^{\mathrm{GENE}}{\left(𝐗+{𝝆}_0\right)}^2⟩$
$+\frac{e^2}{2m{c}^2}\left(A_{1\parallel }{\left(𝐗\right)}^2+m{\left(\frac{c}{e}\right)}^2\frac{\mu }{B}{A}_{1\parallel }\left(𝐗\right){\nabla }_{⟂}^2{A}_{1\parallel }\left(𝐗\right)\right)$
$H_{2,\mathrm{FLR}}^{\mathrm{ORB5}}=-\frac{m{c}^2}{2B^2}{\left|\nabla {\varphi }_1\left(𝐗\right)\right|}^2$
$+\frac{e^2}{2m{c}^2}\left(A_{1\parallel }{\left(𝐗\right)}^2+m{\left(\frac{c}{e}\right)}^2\frac{\mu }{B}{A}_{1\parallel }{\nabla }_{⟂}^2{A}_{1\parallel }\left(𝐗\right)\right)$

Equations

$$\mu =m{v}_{\perp }^2/2B,$$

$$⟨\psi (\mathbf{X}+{\mathbf{\rho }}_0)⟩=\frac{1}{2\pi }{\int }_0^{2\pi }d\theta \psi (\mathbf{X}+{\mathbf{\rho }}_0(\theta )),$$

$$\psi =⟨\psi ⟩+\tilde{\psi }.$$

$$L_p=\left(\frac{e}{c}\mathbf{A}+\left(\frac{e}{c}{ϵ}_{\delta }{A}_{1\Vert }+m{v}_{\Vert }\right)\hat{\mathbf{b}}+{ϵ}_{\delta }\frac{e}{c}{\mathbf{A}}_{1\perp }\cdot m{\mathbf{v}}_{\perp }\right)\cdot \dot{\mathbf{X}}+\frac{mc}{e}\mu \dot{\theta }-H.$$

$$p_z=m{v}_{\Vert }+\frac{e}{c}{ϵ}_{\delta }{A}_{1\Vert }(\mathbf{X}+{\mathbf{\rho }}_0).$$

$${\mathbf{A}}^{\ast }$$

$${\mathbf{A}}^{\ast \ast }$$

$$L_p^{\ast }$$

$$L_p^{\ast \ast }$$

$$\mathcal{L}$$

$$\int dV\frac{|{\mathbf{E}}_1{|}^2}{8\pi }+\int dWdV{F}_0\frac{m{c}^2}{2B^2}|{\nabla }_{\perp }{\varphi }_1{|}^2$$

$${\mathbf{B}}_{1\perp }=\hat{\mathbf{b}}\times \nabla A_{1\Vert },$$

$${\mathcal{A}}^{ORB5}={\int }_{t_0}^{t_1}dt{\mathcal{L}}^{ORB5}$$

$$0=\frac{\delta {\mathcal{L}}^{ORB5}}{\delta {\varphi }_1}\circ {\hat{\varphi }}_1\Rightarrow 0$$

$$0$$

$${\mathcal{A}}^{\mathrm{GENE}}={\int }_{t_0}^{t_1}dt{\mathcal{L}}^{\mathrm{GENE}}$$

$${\mathbf{B}}^{\ast \ast }\approx \mathbf{B}+\frac{c}{e}{p}_{\Vert }\nabla \times \hat{\mathbf{b}}-\hat{\mathbf{b}}\times \nabla ⟨A_{1\Vert }⟩,$$

$$B_{\Vert }^{\ast \ast }=\hat{\mathbf{b}}\cdot {\mathbf{B}}^{\ast \ast }\approx B_{\Vert }^{\ast }.$$

$$0=\frac{\delta {\mathcal{L}}^{\mathrm{GENE}}}{\delta {\varphi }_1}\circ {\hat{\varphi }}_1$$

$$0=\frac{\partial {\mathcal{L}}^{\mathrm{GENE}}}{\partial A_{1\Vert }}\circ {\hat{A}}_{1\Vert }$$

$$0$$

$$\frac{\delta \mathcal{A}}{\delta \mathbf{Z}}=0\Rightarrow \frac{\delta \mathcal{L}}{\delta \mathbf{Z}}=0$$

$$\frac{d}{dt}\frac{\partial L_p}{\partial \dot{\mathbf{Z}}}=\frac{\partial L_p}{\partial \mathbf{Z}}.$$

$$H=H_0+{ϵ}_{\delta }{H}_1,$$

$$\dot{\mathbf{X}}$$

$${\dot{p}}_z$$

$$\dot{\mathbf{X}}$$

$${\dot{p}}_{\Vert }$$

$$0=\frac{d(F_0+{ϵ}_{\delta }{F}_1)}{dt}=\frac{\partial (F_0+{ϵ}_{\delta }{F}_1)}{\partial t}+\dot{\mathbf{X}}\cdot \nabla (F_0+{ϵ}_{\delta }{F}_1)+\dot{p}\frac{\partial (F_0+{ϵ}_{\delta }{F}_1)}{\partial p},$$

Full text

Here we present a new and generic theoretical framework and specific numerical applications to test the faithfulness of the implemented models to theory and to verify the domain of applicability of existing GK codes. For a sound verification process, the underlying theoretical GK model and the numerical scheme must be considered at the same time, which has rarely been done and therefore makes this approach pioneering. At the analytical level, the main novelty consists in using advanced mathematical tools such as variational formulation of dynamics for systematization of basic GK code’s equations to access the limits of their applicability. The verification of numerical scheme is proposed via the benchmark effort.

In this work, specific examples of code verification are presented for two GK codes: the multi-species electromagnetic ORB5 (PIC) and the radially global version of GENE (Eulerian). The proposed methodology can be applied to any existing GK code. We establish a hierarchy of reduced GK Vlasov-Maxwell equations implemented in the ORB5 and GENE codes using the Lagrangian variational formulation.At the computational level, detailed verifications of global electromagnetic test cases developed from the CYCLONE Base Case are considered, including a parametric $\beta$-scan covering the transition from ITG to KBM and the spectral properties at the nominal $\beta$ value.

1 Introduction

The gyrokinetic theory represents one of the most important theoretical frameworks for theoretical and numerical modeling of magnetised plasmas. Historically, it has been considered in the context of fusion plasmas as an accurate tool for assessment of turbulent transport, which in its turn, represents a serious issue for plasma confinement.

At the same time, more recently, gyrokinetic simulations have also been applied to astrophysics [21], [23] in order to access turbulence at the large spectra of scales, and especially at the small scales, where the MHD approximations fail.

Different numerical implementations of the gyrokinetic equations have been actively developed during the last decades. Since the pioneering implementation of the gyrokinetic Vlasov-Poisson system realised by W. W. Lee in $1983$ [17] in the framework of the Particle-In-Cell method, that approach has undergone important developments and is very popular nowadays, see for example [13], [15]. At the same time, the Eulerian numerical realisation of the gyrokinetic Vlasov-Maxwell equations, which appeared in $1995$ [18] attracted a similar popularity, for example [14], [10]. For a detailed review on gyrokinetic simulation see for example [9].

Access to High Performance Computing (HPC) facilities allowed the theoretical plasma physicists to bring gyrokinetic codes to a significantly more advanced level of realism. However, the question about the validity of those advanced numerical tools employed for the investigation of new physics has not been sufficiently investigated yet.

Therefore, a two fold verification framework, which allows one to simultaneously verify the implemented model together with the numerical scheme for the gyrokinetic codes needs to be established. For this purpose, one should deal with two main types of difficulties. The first difficulty consists in the understanding of the gyrokinetic models implemented in a given gyrokinetic code. This is mainly related to use of the different nomenclature and different orderings for the derivation of the equations of the model. The second difficulty comes from the fact that, the implementation typically uses different discretizations and approximations which may further alter the results.

A systematic derivation, which guarantees the energetic consistency of gyrokinetic models requires advanced mathematical tools such as differential geometry and variational calculus on functional spaces. In the electrostatic limit, i.e. in the case of the Vlasov-Poisson equations, a systematic theoretical derivation from the first principles of dynamics has been presented in [1] . For the numerical schemes, numerous verification studies can be found – amongst others, for instance, the cross-code benchmarks in Refs. [7], [19] and [16].

Regarding the full gyrokinetic Vlasov-Maxwell system including electromagnetic fluctuations, several cross-code comparisons exist in the flux tube limit [3, 20, 2]. However, detailed analytic comparisons and radially global benchmarks are hardly to be found which has therefore been identified as one of the goals for the European VeriGyro project launched in 2014. In very recent works the theoretical foundations for the Particle-In-Cell codes [24] and the result of cross-code benchmark for global electromagnetic gyrokinetic codes,[11] have been summarized.

In this work we the present theoretical framework for the systematic derivation of the ORB5 and GENE codes models as two generic representers of the PIC and Eulerian implementations of the gyrokinetic Vlasov-Maxwell equations. We also briefly present two inter-code benchmark test cases in order to explicitly illustrate the differences between the implemented models.

In Sec. 2 we review the main ideas behind the gyrokinetic dynamical reduction, paying special attention to comparison between the theoretical assumptions and the code implementations, for example, the difference between the theoretical and code orderings is explained. In Sec. 3 we start with a comparison between the reduced particle (gyrocenter) models for the ORB5 and GENE codes followed in Sec. 4 by the presentation of Lagrangian variational formulations for the reduced Maxwell-Vlasov models for both codes. Section 5 contains the intercode benchmark test cases specific for the code models comparison. The conclusions are summarized into the Sec. 6.

2 Gyrokinetic dynamical reduction

Gyrokinetic theory is based on the procedure of an asymptotic dynamical reduction for a multi-scaled dynamical system. Such a procedure is aiming at consrtucting a new set of phase space variables in which the dynamics of a system is restricted on a hypersurface, i.e. an invariant of the motion is serving as one of the phase space variables. To develop such an asymptotic dynamical reduction, first of all one needs to put in evidence the fact that there exists a dynamical invariant in the considered system and then to define orderings assumed during the asymptotic procedure construction.

The idea of the gyrokinetic dynamical reduction is based on the physical property of the considered system: the presence of a strong background magnetic field induces a scales of motion separation on the dynamics of the charged particles moving in the superposition of that strong background guide field and some additional fluctuating perturbative electromagnetic fields.

In fact, particle dynamics is decomposed into the fast rotation around the magnetic field lines and slow drift motion. The cyclotron frequency $\Omega=eB/mc$, where $e$ and $m$ are, respectively, the charge and mass of particles, $B$ is the magnetic field amplitude and $c$ is the speed of light, sets the scale of gyromotion.

The gyromotion is described by a fast gyroangle variable $\theta$ to which corresponds to a canonically conjugated slowly varying magnetic moment $\mu$. At the lowest order:

[TABLE]

where $v_{\perp}$ is the perpendicular velocity of particles with respect to the magnetic field lines. In the case of a constant and uniform background magnetic field, $\mu$ is an exact dynamical invariant.

The sources for the violation of the magnetic moment invariance can be attributed to two different reasons. First, the spatial variation of background quantities such as the magnetic field non-uniformity and curvature and, second, the presence of electromagnetic fluctuations. The gyrokinetic dynamical reduction uses the fact that averaged over the long times magnetic moment is still conserved, i.e. $\left\langle\dot{\mu}\right\rangle_{t}=0$ even for the perturbed system.

2.1 Gyrocenter phase space variables

The goal of the gyrokinetic dynamical reduction consists in building up a new set of phase space variables, such that the $\theta$ dependence is completely uncoupled and the magnetic moment $\mu$ has a trivial dynamics, i.e. $\dot{\mu}=0$. Therefore, the reduced particle dynamics is described in the $5$-dimensional phase space with variables $\left(\mathbf{X},p;\mu\right)$, where $\mathbf{X}$ represents the reduced particle position and $p$ is a scalar moment coordinate and the magnetic moment $\mu$ has a trivial dynamics. This change of coordinate is constructed via a perturbative series of near-identity phase space transformations. These transformations are invertible at each step of the perturbative procedure. The particle dynamics on the new reduced phase space is derived within the same near-identity phase space transformation procedure, which is performed on the corresponding phase space Lagrangian. We give the expression in Sec. (3).

The reduced position $\bf X$ has a simple geometrical meaning: it is the instantaneous center of the fast particle’s rotation around the magnetic field line. Therefore, from the space coordinate viewpoint the gyrokinetic transformation is a shift between the initial particle coordinate $\bf x$ and the instantaneous center of its rotation $\bf X$. The difference between both positions is the polarization displacement, defined at the lowest order of the dynamical reduction procedure as the Larmor vector of the particle $\bm{\rho}_{0}=mc/e\ (2\mu/mB)^{1/2}\widehat{\bm{\rho}}$, where $\widehat{\bm{\rho}}$ represents a fast rotating unit vector, perpendicular to the direction of the background magnetic field $\widehat{\mathbf{b}}$ and explicitly depending on the fast variable $\theta$. In this work, the exact gyrocenter spatial coordinate transformation is not considered, but instead the lowest order polarization displacement $\bm{\rho}_{0}$ will be taken in account: it will be shown to be sufficient in particular to expose the main differences between the ORB5 and GENE code models. Performing numerical simulations on the $5$-dimensional phase space instead of the $6$-dimensional one results in a drastic reduction of the computational costs.

2.2 Numerical schemes

First of all, let us clarify which system of equations from the mathematical point of view we are about to solve. In the gyrokinetic model, we have an equation of evolution for the particle distribution function $F$ on the $5$-dimensional phase space, i.e. we have to solve the new partial differential equation on the $5$-dimensional phase space $({\mathbf{X}},p,\mu)$. This can be solved either by advancing $F$ in time on the $5$D grid (Eulerian method) or by advancing in time the corresponding characteristics, which are $5$ non-linearly coupled ODEs (PIC). In what concerns the field equations, we are dealing with solving the $3$D integro-differential equations: Poisson (quasi-neutrality) and Ampère, which are using the information about the distribution function $F$ in order to evaluate charges and currents. To summarize: solving the system of the gyrokinetic Vlasov-Maxwell equations is a rather challenging numerical task. We provide details of the derivation of these equations below.

In this paper we compare the gyrokinetic Vlasov-Maxwell models implemented in two gyrokinetic codes, which are using two different numerical schemes.

On the one hand, we consider the ORB5 code, which is using the Lagrangian approach for solving the gyrokinetic equations of motion. That approach consists in sampling initial positions in phase space (loading of markers), then following marker orbits in $5$D (pushing) and obtaining the source terms (charge and density) for the field equations at the each time step. On the other hand, we consider the GENE code, which is using an Eulerian numerical scheme in order to solve the system of the gyrokinetic Vlasov-Maxwell equations. It consists in discretizing the phase space on a fixed grid, and applying finite differences or finite volume schemes for the differential and integral operators. The Eulerian approach is sometimes also called the ”Vlasov” approach.

The GENE code has two versions from the geometrical point of view: the global and the local version. The local version is also called the flux-tube code, in which the domain considered is a vicinity of a magnetic field line. The equations of motion are expanded in that vicinity such that all coefficients, which define the density and temperature profiles are constant: $\bm{\nabla}n=\mathrm{const}$, $\bm{\nabla}T=\mathrm{const}$ as well as the geometrical coefficients, i.e. the elements of the metric tensor and in particular $q=\mathrm{const}$ with a constant magnetic shear $\hat{s}=\mathrm{const}$ .

The ORB5 code has only the global version. It means that it takes the geometry of the whole plasma domain into account with consistent plasma profiles and gradients, as well as full metric of the background axisymmetric magnetic geometry.

2.3 Gyrokinetic theoretical ordering

From the point of view of the two-step derivation procedure, in which the effects of the magnetic moment non-invariance induced by the background fields non-uniformities are considered separately from those related to the presence of the electromagnetic fluctuating fields, the small parameters can be organised in two groups. In the first group of parameters, i.e. guiding-center transformation related parameters, we include those related to the variations of the background related quantities, i.e. $\epsilon_{B}=\rho_{\mathrm{th}}/L_{B}$, where $\rho_{\mathrm{th}}$ is the thermal Larmor radius of the particle and $L_{B}=|\bm{\nabla}B/B|^{-1}$ sets up the length scale of the background magnetic field variation.

The second group of the small parameters is related to the gyrocenter coordinate transformation, aiming to reestablish the invariance of magnetic moment $\mu$, destructed by the presence of fluctuating electromagnetic fields. The associated small parameter can be defined as $\epsilon_{\delta}\sim|\mathbf{B}_{1}|/B\sim|\mathbf{E}_{1\perp}|/(Bv_{\mathrm{th}})\sim\left(k_{\perp}\rho_{\mathrm{th}}\right)\ e\phi_{1}/T_{i}\equiv\epsilon_{\perp}e\phi_{1}/T_{i}$, where $v_{\mathrm{th}}$ is the thermal velocity, $B$ the amplitude of the background magnetic field, $T_{i}$ is the ion temperature and $\phi_{1}$ represents the amplitude of the fluctuating electrostatic potential. The parameter $\epsilon_{\perp}$ allows the distinction between the gyrokinetic theory with $\epsilon_{\perp}\sim\mathcal{O}(1)$ and the drift-kinetic theory with $\epsilon_{\perp}\ll 1$. In particular, the model which is truncated up to the second order in $\epsilon_{\perp}$ is called the long-wavelength approximation of the gyrokinetic theory. The long-wavelength approximation is rather popular for the numerical implementations. We will come to that model later in Sec. 3.

The parallel fluctuations of the electric field are pushed on the next level of smallness according to the anisotropy of the turbulence $|E_{1\|}|/|E_{1\perp}|\sim|k_{\|}|/|k_{\perp}|\ll 1$.

In addition to that, for parallel fluctuations of the magnetic field it is assumed that $|B_{1\|}|/|B|\sim\beta\epsilon_{\delta}$, where $\beta=(8\pi p_{\mathrm{th}})/B^{2}$, the ratio between the kinetic thermal pressure $p_{\mathrm{th}}$ and the magnetic pressure $B^{2}$, which means that they are only considered when the magnetic beta is close to $1$, i.e. $\beta\sim\mathcal{O}(1)$.

In the case of the maximal ordering, the effects from the variations of background quantitites and the electromagnetic fields fluctuations should be considered at the same order $\epsilon_{B}\sim\epsilon_{\delta}$. However, it is important to note, that each set of small parameters, related to a specific choice of physical configuration, will define a setup for the derivation of the corresponding gyrokinetic theory. In other words, it is important to emphasize that there is not only one specific set of reduced Vlasov-Maxwell equations, which is defined as the gyrokinetic Vlasov-Maxwell system, but rather different sets of reduced equations, which need to be put inside the common systematic framework.

2.4 Gyrokinetic code ordering

The majority of the gyrokinetic models implemented in the codes are derived within the assumption that $\epsilon_{B}\ll\epsilon_{\delta}$. Typically, all the background gradient corrections are taken into account at the first order, while the contributions from the fluctuating fields are considered at the second order. Such a treatment allows one to eliminate a significant number of terms, (see for example derivations performed within the maximal ordering in [5] and [26]) and simplified numerical implementation.

In what concerns the FLR or the $\epsilon_{\perp}$- ordering, both models: derived in the limit with full FLR corrections as well as the models truncated up to the second order in $\epsilon_{\perp}$ are implemented. Below, we will discuss differences between the long-wavelength approximated Vlasov-Maxwell models and those containing the full FLR corrections. In addition to that, we note that the ORB5 and the radially global version of the GENE code considers the perpendicular fluctuations of the magnetic field only i.e. are restricted to the low-beta ordering only, while the local version of the GENE code follows the full $\beta$ ordering.

2.5 Gyroaveraging

Implementing a gyroaveraging operator is an important and challenging task for the gyrokinetic simulations.We do not focus on the numerical details here, but just provide the theoretical definition, which is necessary for the models derivation.

For each function $\psi$, defined in the position $\mathbf{x}=\mathbf{X}+\bm{\rho}_{0}$ we define the gyroaverage operator as:

[TABLE]

such that each function on the reduced phase space can be decomposed in the gyroaveraged and fluctuating parts:

[TABLE]

3 GENE and ORB5 gyrocenter models

The $6$-dimensional reduced phase space Lagrangian represents a central object of the gyrokinetic dynamical reduction:

[TABLE]

The last term in the expression above is referred to as the Hamiltonian of the system. The remaining terms represent the symplectic part of the phase space Lagrangian. This phase space Lagrangian contains the fluctuating electromagnetic fields $(\phi_{1},A_{1\|},\mathbf{A}_{1\perp})$, which are explicitly time-dependent. The perturbed fields are evaluated at the spatial position ${\mathbf{x}}={\mathbf{X}}+\bm{\rho}_{0}$, i.e. they possess the fast gyroangle $\theta$ dependency, which is removed according to the near-identity phase space transformation. That coordinate transformation at the code relevant order is presented in [25], for the procedure at all orders, see for example [6]. Here we skip the detailed description of the reduction procedure and rather focus on its conceptual explanation and comparison between the reduced particle models, issued from the various representations of that procedure.

From the conceptual point of view, the fluctuating electric field $\phi_{1}$, as a scalar field, is included into the Hamiltonian $H$, whereas, for the magnetic field perturbations $A_{1\|}$ and ${\mathbf{A}}_{1\perp}$ different options are possible. The choice of including the perturbed magnetic fields $A_{1\|}$ and $\mathbf{A}_{1\perp}$ in the symplectic part or in the Hamiltonian defines its dynamical representation. The right choice of the representation is important for the corresponding numerical scheme realisation. Different dynamical representations are possible, the complete list with the corresponding nomenclature, which we are following below, is available in [6].

Consideration of the perpendicular part of the perturbed electromagnetic potential $\mathbf{A}_{1\perp}$ depends on the choice of the ordering. In the case of the low$\beta$ ordering, when parallel perturbations of magnetic field are neglected, the $\mathbf{A}_{1\perp}$ related term does not appear in the Eq.  (4).

It is possible to identify the near-identity phase space transformations, which would affect the Hamiltonian part of the phase space Lagrangian only. For realisation of such a transformation, one should always keep the fluctuating parts of all the perturbed electromagnetic potentials in the Hamiltonian part of $L_{p}$ and leave the symplectic part free of the explicit fast angle $\theta$ dependency. Here we propose two different options for the definition of such a transformation.

The first option is the Hamiltonian representation (adopted for the ORB5 code) which leaves the symplectic part of the phase space Lagrangian Eq. (4) completely free of the perturbative electromagnetic potentials, by including them into the Hamiltonian. That manipulation is possible when using the canonical parallel gyrocenter moment as one of the phase space variables:

[TABLE]

The second option is to redistribute the fluctuating potentials between the Hamiltonian and the symplectic parts of the phase space Lagrangian, so that it leaves the symplectic part free of the $\theta$-dependency. This corresponds to the parallel-symplectic representation. In that case, only the gyroaveraged component of the parallel part of the perturbed magnetic potential $\left\langle A_{1\|}\right\rangle$ is contained in the symplectic part of the phase space Lagrangian and its fluctuating part together with the perpendicular part of the perturbed magnetic potential is accounted in the Hamiltonian. This parallel-symplectic representation is used for the GENE code, however with some further approximations as explained below.

The parallel-symplectic representation, from the theoretical point of view is using a modified parallel moment consisting of the sum of the parallel kinetic moment and the fluctuating part of the parallel magnetic potential $\bar{p}_{\|}=mv_{\|}+\frac{e}{c}\widetilde{A}_{1\|}(\mathbf{X}+\bm{\rho}_{0})$ of the particle as the phase space variable. However, the fluctuating component $\widetilde{A}_{1\|}$ is neglected in the GENE code model and therefore the phase space variable is $p_{\|}=mv_{\|}$, i.e. the kinetic moment.

We define the symplectic magnetic potentials in the following form:

[TABLE]

therefore the reduced particle phase space Lagrangians are

[TABLE]

where $H^{*}$ and $H^{**}$ are the Hamiltonians, corresponding to the hamiltonian or parallel-symplectic representation of the dynamics, respectively.

As we can see, in the first case the symplectic part of the phase space Lagrangian is time independent and in the second case the symplectic part of the phase-space Lagrangian is explicitly time-dependent, because of the presence of $\left\langle A_{1\|}\right\rangle$. That has a direct impact on the equations for the corresponding characteristics of the dynamical variables on the reduced phase space. The Hamiltonian representation allows one to avoid explicit time derivative of the gyroaveraged parallel magnetic potential $\left\langle A_{1\|}\right\rangle$ on the r.h.s. of the reduced phase space $(\mathbf{X},p_{\|})$ characteristics, while the parallel-symplectic representation makes it appear explicitly in the characteristic for the kinetic parallel moment $p_{\|}$, see for example [24].

In the table below we summarize the reduced particle models implemented in both of the codes: first, we explicit the symplectic and then the Hamiltonian parts of the phase-space Lagrangian.

In table I, the modified and parallel-symplectic model means that the fluctuating component of the magnetic potential $\widetilde{A}_{1\|}$ has been ignored in the definition of the phase space variable $\bar{p}_{\|}$, which in its turn leads to neglect the $\left\langle\widetilde{A}_{1\|}(\mathbf{X}+\bm{\rho}_{0})^{2}\right\rangle$ term into the expression for the second order Hamiltonian; additionally the square of the perpendicular component of the electromagnetic potential $A_{1\perp}\left({\mathbf{X}}+\bm{\rho}_{0}\right)^{2}$ has also been omitted.

4 Gyrokinetic Vlasov-Maxwell equations

There exist two significantly different approaches for the derivation of the reduced gyrokinetic Vlasov-Maxwell equations. First of all, one can proceed with the direct calculation of the moments for the reduced gyrokinetic Vlasov distribution function in order to evaluate the charge and the current densities in the gyrokinetic Poisson and Ampère equations (zeroth order moment corresponding to gyrokinetic Poisson and the first moment to gyrokinetic Ampère equation). This approach is called the ”pull-back” transformation. Another possible approach is to get the reduced gyrokinetic field equations from the variational formulation in which the interaction between the reduced particle dynamics, described by the particle Lagrangian $L_{p}$ and the dynamics of the electromagnetic fields is included inside the field-particle Lagrangian.

The first approach is an intuitive approach, it has been introduced, for instance, in Ref. [12] within the electrostatic approximation. On the other hand, the variational approach is more formal and it has been formulated almost two decades later in Refs. [22] and [4]. For historical reasons, the GENE code equations have been derived within this intuitive approach, while the ORB5 model is already obtained within the variational framework.

To make a formal comparison between the models of both codes in this work we choose to derive the theoretical code models from the variational formulation of the reduced dynamics [22].

4.1 Gyrokinetic field formulation

In this section we present the variational framework for the derivation of consistently coupled gyrokinetic Vlasov-Maxwell equations suitable for a code implementation. The first term of the field-particle Lagrangian includes the reduced particle dynamics represented by the gyrocenter Lagrangian $L_{p}$ coupled to the Vlasov distribution function $F$; the second term contains the electromagnetic fields:

[TABLE]

where the reduced phase space variables are $\mathbf{Z}=\left({\mathbf{X}},p,\mu,\theta\right)$ with $d{\mathbf{X}}=dV$ and the reduced velocity phase space volume is chosen according to the representation of the reduced particle dynamics $dW=B_{\|}^{\#}dp\ d\mu$, i.e. $B_{\|}^{\#}=(B^{*}_{\|},B^{**}_{\|})$ and $p=(p_{z},p_{\|})$.

The main important idea about building up a consistent gyrokinetic Vlasov-Maxwell model implementable into the given code consists in the fact that all the approximations should be performed on the field-particle Lagrangian before the derivation of the gyrokinetic Vlasov-Maxwell equations. The corresponding equations of motion should be obtained according to the first principle of dynamics together with the corresponding conservation laws following the Noether method.

4.2 ORB5 and GENE codes models derivation

Here, we focus on the derivation of the ORB5 and GENE theoretical models. By theoretical we mean that the derived equations are written down in the form, which follows directly from the analytical calculation. Note that each numerical scheme has its own requirement for model rewriting before discretization, aiming to simplify the numerical resolution. We do not focus on the detailed comparison of the discretized models here. In order to perform the verification of the numerical implementations we realise a detailed intercode benchmark, linear results can be found in [11]. Two test cases connected to the theoretical models verification are presented in Sec. 5.

We provide a detailed list of approximations performed on the field particles Lagrangian (10) in order to obtain the models corresponding to the ORB5 and GENE codes. We start with presenting the approximations common to both codes and and subsequently focus on the specific details of each model.

The first common approximation considering the derivation of the ORB5 and GENE code models consists in the fact that the field-particle Lagrangian (10) is truncated up to the second order i.e. contains up to the $\mathcal{O}(\epsilon_{\delta}^{2})$ electromagnetic field terms. It means that the second order Lagrangian couples nonlinear terms related to the reduced particle dynamics $H_{2}$, to the background (non-dynamical) Vlasov distribution function $F_{0}$ only. Therefore, the corresponding gyrokinetic Vlasov equation contains exclusively the linear (i.e. $\sim\mathcal{O}(\epsilon_{\delta})$) terms.

The second common approximation, directly applied to the field particle Lagrangian (10) is the quasi-neutrality approximation. We compare the field term with the nonlinear second order particle contribution. According to the Tab. 1 both Hamiltonian second order reduced particle models $H_{2}^{\mathrm{ORB5}}$ and $H_{2}^{\mathrm{GENE}}$ coincide in the electrostatic limit . Therefore, in the long-wavelength approximation $H_{2,\mathrm{FLR}}^{\mathrm{es}}=-mc^{2}/2B^{2}|\bm{\nabla}_{\perp}\phi_{1}|^{2}$. With taking into account the gyrokinetic ordering for the electrostatic field with $E_{1\|}/E_{1\perp}=\epsilon_{\delta}$, we neglect the parallel contribution $E_{1\|}=\nabla_{||}\phi_{1}-1/c\ \partial_{t}A_{1\|}$ to the electric field $\mathbf{E}_{1}$ in the second term of the Eq. (10). Therefore,

[TABLE]

$\alpha=0$ corresponds to the low-beta approximation of the reduced dynamics and $\alpha=1$, to the finite $\beta$ approximation. In the strongly magnetised plasma the ratio of the sound Larmor radius and the Debye length $\rho_{s}^{2}/\lambda_{d}^{2}=c^{2}/v_{A}^{2}\gg 1$, i.e. we can systematically neglect the term $1/8\pi\int dV\left|\bm{\nabla}_{\perp}\phi_{1}\right|^{2}$ in Eq. (4.2).

The last common code model approximation is performed on the perturbed part of the magnetic field, i.e.

[TABLE]

which means that we have neglected the $A_{1\|}\bm{\nabla}\times\widehat{\mathbf{b}}\sim\mathcal{O}(\epsilon_{\delta}\epsilon_{B})$ term according to the general codes ordering, which we have discussed in Sec. 2.4.

4.3 Gyrokinetic Maxwell equations for ORB5

The second order linearised field-particle Lagrangian used in ORB5, including the second order Hamiltonian $H_{2}^{\mathrm{ORB5}}$ as written in Tab. 1 reads:

[TABLE]

where we have assumed the approximation on the magnetic field given by Eq. (12). With using the first principle of dynamics, we derive the gyrokinetic Vlasov-Maxwell equations corresponding to field-particle Lagrangian given by Eq. (13). Here we limit our derivation to the weak form of the equations of motion, it means, including the arbitrary test function $\widehat{A}_{1\|}$ and $\widehat{\phi}_{1}$. The weak form is suitable for the finite element discretisation implemented in the ORB5 code. We start our derivation with the gyrokinetic quasineutrality equation:

[TABLE]

The Ampère equation is:

[TABLE]

4.4 Gyrokinetic Maxwell equations for GENE

The second order (i.e. containing terms up to $\mathcal{O}(\epsilon_{\delta}^{2})$) linearised field-particle Lagrangian action with the second order Hamiltonian $H_{2}^{\mathrm{GENE}}$ for the GENE code is given by:

[TABLE]

where we have taken into account the fact that the term $A_{1\|}\bm{\nabla}\times\widehat{\mathbf{b}}\sim\mathcal{O}(\epsilon_{\delta}\epsilon_{B})$ is neglected, so that only the part of $\mathbf{B}_{1\perp}$ given by Eq. (12) is taken into account, such that ${\mathbf{B}}_{1}=\alpha\bm{\nabla}\times{\mathbf{A}}_{1\perp}+\widehat{\mathbf{b}}\times\bm{\nabla}A_{1\|}$, with $\alpha=0$ for the global and $\alpha=1$ for the local (flux-tube) code. This approximation further affects the field-particle Lagrangian given by Eq. (16): the symplectic magnetic field $\mathbf{B}^{**}=\bm{\nabla}\times\mathbf{A}^{**}$ is approximated up to:

[TABLE]

which means that the parallel component of the symplectic magnetic field coincides with the one from the Hamiltonian representation of dynamics, i.e.

[TABLE]

In addition to that the volume element of the velocity phase space $dW$ is further approximated up to $dW=B\ dp_{\|}\ d\mu\ d\theta$.

Following the same procedure as for the ORB5 code, we derive the corresponding gyrokinetic Vlasov-Maxwell equations from the first principle of dynamics in the weak form. The weak form is considered for the purpose of comparison of the models implemented in ORB5 and GENE codes in the same form.

We start with the gyrokinetic quasineutrality equation:

[TABLE]

The parallel component of the gyrokinetic Ampère equation:

[TABLE]

Finally, for the perpendicular component of the gyrokinetic Ampère equation (for the local GENE code only, i.e. with $\alpha=1$):

[TABLE]

4.5 Discussion

Let us now consider the main differences between the gyrokinetic Maxwell equations for both codes. The quasineutrality and the linear parallel Ampère equations for the ORB5 code and the ones for the global version of the GENE code differs only with respect to the phase space volume element: ORB5 code uses $B_{\|}^{*}$ and GENE implementation uses $B$. Terms proportional to the perturbed parallel magnetic potential $A_{1\|}$ appearing in the last term of the parallel Ampère equation are available in the ORB5 code up to the second order FLR decomposition. We note that those terms are not present in the GENE code. Finally, the ORB5 code and the global GENE code do not consider the perpendicular Ampère equation yet.

4.6 Gyrokinetic Vlasov equations for ORB5 and GENE codes

The gyrokinetic Vlasov equations for the ORB5 and GENE codes are reconstructed from the characteristics (i.e. dynamical equations for the phase space variables $\mathbf{X}$ and $p$). The characteristics, in their turn, are obtained from the fields-particles Lagrangian action, corresponding to the codes given by Eqs. (13) and (16).

[TABLE]

and therefore, the Euler-Lagrange equations for the characteristics are given by:

[TABLE]

Both codes are using the first order in $\epsilon_{\delta}$, linearized characteristics, i.e. corresponding to the linearized Hamiltonian

[TABLE]

where $H_{0}$ and $H_{1}$ are chosen correspondingly to the code model from the Tab. 1.

In the case of the ORB5 code, which corresponds to the Hamiltonian representation of dynamics, we have:

[TABLE]

The detailed derivation of the equations of motion for the ORB5 code in the case of $H_{2\ \mathrm{FLR}}^{\mathrm{ORB5}}$ and their comparison with the theoretical model, containing up to the second order terms in $\epsilon_{\delta}$ can be found in [24].

While in the case of the GENE code, which corresponds to the modified parallel-symplectic representation of dynamics, we have:

[TABLE]

where the numerical implementation into GENE takes into account the approximation on the symplectic magnetic field $\mathbf{B}^{**}$ given by Eqs. (17) and (18).

For both of the codes, the gyrokinetic Vlasov equation is reconstructed from the characteristics in the following way:

[TABLE]

where $p=(p_{z},p_{\|})$ depending on the code.

Let us now analyze the terms in the Vlasov equation. First of all, both codes take into account that the background distribution function is time independent, i.e. $dF_{0}/dt=0$.

In what concerns the ORB5 code all the nonlinear terms are implemented and no approximation has been made for the symplectic magnetic field $\mathbf{B}^{*}$, which corresponds to the Hamiltonian representation. The latter guarantees that the phase space volume is preserved following the gyrokinetic coordinate transformation. For the GENE code, however, the symplectic magnetic field is approximated up to the term $\left\langle A_{1\|}\right\rangle\bm{\nabla}\times\widehat{\mathbf{b}}\sim\mathcal{O}(\epsilon_{B}\epsilon_{\delta})$, an additional ordering is implemented, which allows one to eliminate the $\partial_{p_{\|}}$-derivatives in non-linear terms of the Vlasov equation, given by Eq. (27) is implemented.

5 Numerical verification: ORB5/GENE Benchmark

Establishing a framework for the detailed comparison of the basic equations naturally represents just one pillar in a code comparison. The second is to verify the correct numerical implementation of the underlying model by either comparisons with analytic results or via benchmarks with codes based on entirely different numerical approaches. In this section, we present examples for the latter which have been realised within the framework of the EUROfusion project VeriGyro as well. The detailed setups and results involving up to five different codes can be found in [11]. Here, however, we focus on the results from the PIC based ORB5 and the Eulerian code GENE in order to complement the theoretical model comparison introduced in the previous sections. In particular, we will focus on two specific test cases. The first one aims at cross-checking the stabilization of ion temperature gradient (ITG) driven modes with $\beta$ and the eventual onset of kinetic ballooning modes (KBMs) above a certain threshold value. Here, we employ linear growth rates and real frequencies as observables. Since comparisons of linear modes are not limited to these two quantities, a second test case at a given toroidal mode number $n$ and for a fixed finite $\beta$ is presented. Here, full poloidal and radial profiles of the electrostatic $\phi_{1}$ and the parallel electromagnetic $A_{1\|}$ potentials are considered. These extended comparisons furthermore reveal and illustrate the impact of including the full FLR corrections into the model for the reduced particles dynamics. The latter is achieved by comparing the full FLR solver for the gyrokinetic Poisson equation implemented in GENE code with the long-wavelength approximation solver implemented in the ORB5 code version, which has been used during the benchmark campaign.

5.1 Electromagnetic $\beta$ scan at fixed wave number

In Fig. 1, we present the result of the electromag- netic $\beta$-scan at the fixed toroidal wave number $n=19$. As mentioned in Sec.2.4 both ORB5 and GENE global codes are using the same low $\beta$ ordering (both codes neglect the $A_{1\perp}$ part of the electromagnetic potential). As one can see, the linear growth rates and real frequencies from both codes - marked by red and blue lines - demonstrate good agreement. The plot furthermore contains the results from the local (flux-tube) GENE version maximized over radius and ballooning angle (black lines) in order to quantify finite size effects and emphasize the need for global electromagnetic models in the scenario at hand. The flux tube growth rates are indeed found to be generally higher. More strikingly, the mode transition is observed at a different $\beta$ value since the finite size (so-called $\rho^{*}$) effects seem to be depend on the microinstability type as well. Differences between the local (flux-tube) and global results have been found in previous linear comparisons without electromagnetic effects and have furthermore been confirmed in full electrostatic $\beta=0$ (nonlinear) turbulence simulations, see e.g. [19] and references therein. For the case shown here, we have $\rho^{*}\approx 1/182$ such that mild but visible deviations as found in Fig. 1 are indeed expected.

5.2 Radial and poloidal mode structures at fixed $\beta$

Another example for a direct comparison of linear electromagnetic microinstabilities is found in Fig. 2. Here we compare radial and poloidal mode structures of the electrostatic and magnetic potential. While generally agreement between the two codes at hand is confirmed, they particularly differ with regard to fine-scale structures in the radial profile of the electrostatic potential which can be related to mode rational surfaces and which are absent in the ORB5 profiles. The reason is the choice of a long-wavelength approximation Poisson solver – the full FLR solver, which has been recently implemented in ORB5 [27], should recover the peaks. Considering the good agreement in Fig. 1, missing the radial fine-scale structures does not introduce too much harm to visibly alter the growth rates and frequencies. However, as discussed in [11] in more detail, deviations will become intolerable at the latest if higher toroidal mode numbers $n\gtrsim 40$ are considered. This problem is also remedied by the new solver option in ORB5.

6 Conclusions

The two fold verification framework for gyrokinetic codes has been established. First of all, the models currently implemented in ORB5 (PIC representative code) and GENE (Eulerian representative code) have been derived within the Lagrangian variational framework, which permitted to perform a close comparison of both models and to identify the corresponding approximations. From the theoretical point of view (i.e. before performing all the approximations for further numerical implementations), it has been identified that the linear models for both codes are identical. In addition to that, the approximations performed on the theoretical models before the discretization have been also stated.

On the other hand, the numerical part of the verification framework has shown an excellent agreement in the linear electromagnetic $\beta$ scan test case. Further comparison of the nonlinear models at the theoretical level as well as extension of the benchmark to nonlinear simulations is the part of our future work.

7 Acknowledgment

Authors would like to acknowledge Cristel Chandre for the advice and help during preparation of this manuscript. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training Programme No. 2014- 2018 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.