An exponential integrator for the drift-kinetic model

TL;DR

This paper introduces an exponential integrator for the drift-kinetic equation that removes the CFL constraint, conserves mass, reduces computational effort, and performs efficiently in plasma instability simulations.

Contribution

The paper presents a novel exponential integrator for the drift-kinetic model that improves efficiency, conservation, and allows higher order methods, outperforming traditional splitting approaches.

Findings

Removes CFL constraint in simulations

Achieves mass conservation up to machine precision

Enables larger time steps comparable to existing methods

Abstract

We propose an exponential integrator for the drift-kinetic equation in cylindrical geometry. This approach removes the CFL condition from the linear part of the system (which is often the most stringent requirement in practice) and treats the remainder explicitly using Arakawa's finite difference scheme. The present approach is mass conservative, up to machine precision, and significantly reduces the computational effort per time step. In addition, we demonstrate the efficiency of our method by performing numerical simulations in the context of the ion temperature gradient instability. In particular, we find that our numerical method can take time steps comparable to what has been reported in the literature for the (predominantly used) splitting approach. In addition, the proposed numerical method has significant advantages with respect to conservation of energy and efficient higher order methods can be obtained easily. We demonstrate this by investigating the performance of a fourth order implementation.

Tables

Table 1. Table 1: The run times in seconds for the semi-Lagrangian method (using an Euler solver to determine the feet of the characteristics for the two-dimensional interpolation) and the second order exponential integrator (for the perturbation and the direct formulation) are shown. The coarse problem uses 32 × 32 × 32 × 64 32 32 32 64 32\times 32\times 32\times 64 grid points and integrates up to a final time of T = 8000 𝑇 8000 T=8000 . The fine problem uses 64 × 64 × 64 × 128 64 64 64 128 64\times 64\times 64\times 128 grid points and integrates up to a final time of T = 6000 𝑇 6000 T=6000 .

	semi-Lagrangian	perturbation	direct
coarse ($\Delta t=10$)	2196	482	495
fine ($\Delta t=4$)	54122	13259	13116

Table 2. Table 2: The maximum stable time step Δ ​ t Δ 𝑡 \Delta t for the second and fourth order scheme is listed for both the coarse and fine grid.

	coarse			fine
	direct	perturb.		direct	perturb.
2nd order	15	11		4	4.5
4th order	51	38		4	10

Equations

$${\partial }_{t}f-\frac{{\partial }_{\theta }\varphi }{r}{\partial }_{r}f+\frac{{\partial }_r\varphi }{r}{\partial }_{\theta }f+v{\partial }_{z}f-{\partial }_z\varphi {\partial }_{v}f=0,$$

$$-\left[{\partial }_r^2\varphi +\left(\frac{1}{r}+\frac{{\partial }_r{n}_0(r)}{n_0(r)}\right){\partial }_r\varphi +\frac{1}{r^2}{\partial }_{\theta }^2\varphi \right]+\frac{1}{T_e(r)}\left(\varphi -⟨\varphi ⟩\right)$$

$$=\frac{1}{n_0(r)}{\int }_{\mathbb{R}}rfdv-1,$$

$${\partial }_{t}f+v{\partial }_{z}f=F(f),$$

$$F(f)=\frac{{\partial }_{\theta }\varphi }{r}{\partial }_{r}f-\frac{{\partial }_r\varphi }{r}{\partial }_{\theta }f+{\partial }_z\varphi {\partial }_{v}f.$$

$${\partial }_t\hat{f}+ikv\hat{f}=\mathcal{F}\left(F(f)\right),$$

$$\hat{f}(t)=e^{-ikvt}\hat{f}(0)+{\int }_0^t{e}^{ikv(s-t)}\mathcal{F}\left(F(f(s))\right)ds,$$

$$\hat{f}(t^{n+1})$$

$$+e^{-ikv\Delta t}{\int }_0^{\Delta t}{e}^{ikvs}\mathcal{F}\left(F(f(t^n+s))\right)ds.$$

$${\phi }_{\ell }(z)=z{\phi }_{\ell +1}(z)+\frac{1}{\ell !}$$

$${\phi }_1(z)=\frac{e^z-1}{z},{\phi }_2(z)=\frac{e^z-1-z}{z^2},{\phi }_3(z)=\frac{e^z-1-z-\frac{z^2}{2}}{z^3}.$$

$${\hat{f}}^{n+1}=e^z{\hat{f}}^n+\Delta t{\phi }_1(z)\mathcal{F}\left(F(f^n)\right)$$

$${\hat{f}}^{n+1}={\hat{f}}^n+\Delta t\mathcal{F}(F(f^n))$$

$${\hat{k}}_1$$

$${\hat{f}}^{n+1}$$

$${\hat{k}}_1$$

$${\hat{k}}_2$$

$${\hat{k}}_3$$

$${\hat{f}}^{n+1}$$

$$F(f^n)=\frac{1}{r}\{{\varphi }^n,f^n\}+{\partial }_z{\varphi }^n{\partial }_v{f}^n,$$

$$\{{\varphi }^n,f^n\}\approx D_r(f^n)D_{\theta }({\varphi }^n)-D_{\theta }(f^n)D_r({\varphi }^n),$$

$$\{{\varphi }^n,f^n\}$$

$$J^{++}$$

$$J^{+x}$$

$$J^{x+}$$

$$\sum _{i=1}^N\sum _{j=1}^N[D_{\theta }(\varphi )D_r(f){]}_{i,j}=\frac{1}{4h_r{h}_{\theta }}\sum _{i=1}^N\sum _{j=1}^N({\varphi }_{i,j+1}-{\varphi }_{i,j-1})(f_{i+1,j}-f_{i-1,j}),$$

$$4h_r{h}_{\theta }\sum _{i=1}^N\sum _{j=1}^N[D_{\theta }(\varphi )D_r(f){]}_{i,j}=\sum _{j=1}^N\left[\sum _{i=1}^N{f}_{i,j}({\varphi }_{i-1,j+1}-{\varphi }_{i-1,j-1})-\sum _{i=1}^N{f}_{i,j}({\varphi }_{i+1,j+1}-{\varphi }_{i+1,j-1})\right].$$

$$\sum _{i=1}^N\sum _{j=1}^N({\varphi }_{i+1,j}-{\varphi }_{i+1,j})(f_{i,j+1}-f_{i,j-1})=4h_r{h}_{\theta }\sum _{i=1}^N\sum _{j=1}^N[D_r(\varphi )D_{\theta }(f){]}_{i,j}.$$

$$\sum _{i=1}^N\sum _{j=1}^N[D_{\theta }(\varphi D_r(f)){]}_{i,j}=0.$$

$$\sum _{i=1}^N[D_r(\varphi D_{\theta }(f)){]}_{i,j}=\frac{1}{2h_r}\sum _{i=1}^N(g_{i+1,j}-g_{i-1,j})=\frac{1}{2h_r}(g_{N+1,j}+g_{N,j}-g_{1,j}-g_{0,j}),$$

$$g_{0,j}=-g_{1,j},g_{N+1,j}=-g_{N,j},\forall j=1,\dots ,N.$$

$$\sum _{i=1}^N\sum _{j=1}^N(J_{i,j}^{+x}+J_{i,j}^{x+})=\sum _{i=1}^N\sum _{j=1}^N([D_r(f{D}_{\theta }(\varphi )){]}_{i,j}-[D_r(\varphi D_{\theta }(f)){]}_{i,j}).$$

$$\sum _{i=1}^N\sum _{j=1}^N(J_{i,j}^{+x}+J_{i,j}^{x+})=\frac{1}{h_r}\sum _{j=1}^N(g_{1,j}-g_{N,j})$$

$${\hat{k}}_1$$

$${\hat{f}}^{n+1}$$

Full text

In addition, we demonstrate the efficiency of our method by performing numerical simulations in the context of the ion temperature gradient instability. In particular, we find that our numerical method can take time steps comparable to what has been reported in the literature for the (predominantly used) splitting approach. In addition, the proposed numerical method has significant advantages with respect to conservation of energy and efficient higher order methods can be obtained easily. We demonstrate this by investigating the performance of a fourth order implementation.

1 Introduction

The evolution of strongly magnetized plasmas, like those encountered in tokamak devices, are governed by gyrokinetic type equations, where the highly oscillatory motion of the charged particles around the magnetic field lines is averaged out. In a simplified slab geometry, gyrokinetic models reduce to the drift-kinetic equation, where the unknown $f$ depends on three cylindrical spatial coordinates and one velocity direction. This model is composed of a guiding-center type dynamics in the plane orthogonal to the magnetic field lines and of a Vlasov type dynamics in the parallel (to the magnetic field lines) direction. Moreover, the electric field is determined in a nonlinear way through a quasi-neutrality equation. This drift-kinetic model is helpful to investigate plasma turbulence with an acceptable computational effort, and consequently has received a lot of attention in recent decades.

The goal of this work is to develop accurate and efficient numerical methods for the drift-kinetic model. Several different numerical approaches have been proposed for the drift-kinetic equation or more generally for kinetic equations. Among them is the family of particle in cell (PIC) schemes which only discretizes the self-consistent electric and magnetic fields. Particles are, on the other hand, traced in phase space along the characteristics of the corresponding kinetic equation. Although these methods are relatively cheap from a computational point of view, they suffer from inaccurate results in low density regions and slow convergence as the number of macro-particles is increased (see, for example, [2, 24, 25]).

As an alternative, Eulerian methods, such as finite volumes or finite differences, have been introduced. Since they perform a discretization of the full phase space they are able to overcome many of the difficulties inherent in PIC schemes. It should be duly noted, however, that such schemes are expensive (due to the high dimensionality of the problem) and the time step size is often limited by a stringent CFL condition. The latter deficiency can be overcome by employing semi-Lagrangian schemes. These methods use characteristic curves (to remove the stability constraint) but still maintain a discretization of the full phase space. In a sense, semi-Lagrangian methods can be considered a good compromise between the Lagrangian approach (for example, PIC) and the Eulerian approach. Nevertheless, such methods often require large computational resources. Consequently, it is vital to develop more efficient numerical methods to reduce the computational cost associated with Eulerian or semi-Lagrangian schemes.

Within the semi-Lagrangian approach, splitting methods are very popular. The reason for this is that (as has been realized in the seminal paper by Cheng & Knorr [5]) splitting schemes often decouple different parts of the equation and allow us to perform computations on lower-dimensional slices of phase space. A significant part of the literature has been devoted to the construction and analysis of splitting methods (see, for example, [4, 5, 12, 13, 22]). In addition, a range of space discretization schemes has been proposed. Among these interpolation by cubic splines is perhaps the most popular (see, for example, [22, 28]). Nevertheless, spectral methods and methods based on discontinuous Galerkin type approximations have been considered as well (see, [10, 27]).

In particular, for the drift-kinetic equation, which we consider in this paper, splitting methods have been employed extensively and production-level computer codes are available that implement such schemes (see [16, 17, 23]). However, while for the Vlasov–Poisson equation splitting allows us to reduce the problem to a sequence of one-dimensional advections, this is not true in the present case. In addition, we have to perform a three term splitting compared to only two terms in the Vlasov–Poisson equation. Some high-order numerical methods that avoid splitting have also been proposed for gyrokinetic type models. The time integration is mostly based on explicit Runge-Kutta methods (see [3, 9]) which suffer from a CFL condition.

In the present paper we propose an alternative to what has been considered in the literature. In particular, we will use the observation that in many problems the most stringent CFL condition is associated with the linear part of the drift-kinetic equation. We propose to solve this linear part using Fourier techniques as part of an exponential integrator. Note that high-order exponential integrators are easily available (see, for example, [20]). A somewhat related idea has been used in [21] for a more complex gyrokinetic model. There the linear part (which involves the stiffest terms) is solved implicitly whereas the nonlinear part is treated explicitly. In [23], the linear part is solved using a semi-Lagrangian approach where the function value at the feet of the characteristics is computed using a $4D$ interpolation. In our approach, we exploit the structure of the drift-kinetic model by solving the linear part exactly. This can be done efficiently by using Fourier techniques. It should also be noted that contrary to some of the methods described above, our proposed scheme requires no multi-dimensional interpolation.

The rest of this work is structured as follows. First, the drift-kinetic model is introduced in section 1.1. Then, in section 2, our proposed numerical method is described and its (theoretical) properties are compared to that of the splitting approach. In section 3, we perform numerical simulations that allow us to conduct an extensive comparison of our exponential integrator with a splitting scheme. Finally, we conclude in section 4.

1.1 Drift-kinetic equations

In this work, we focus on the numerical approximation of the drift-kinetic model. More specifically, our goal is to compute an approximation of $f=f(t,r,\theta,z,v)$ using the initial condition $f_{0}(r,\theta,z,v)=f(t=0,r,\theta,z,v)$ and satisfying the following $4D$ slab drift-kinetic equation (for more details see [8, 17])

[TABLE]

for $(r,\theta,z,v)\in\Omega\times[0,L]\times\mathbb{R}$, $\Omega=[r_{\rm min},r_{\rm max}]\times[0,2\pi]$, where $r_{max}>r_{min}>0$. The self-consistent potential $\phi=\phi(r,\theta,z)$ is determined from the quasi-neutrality equation

[TABLE]

with $\langle\phi\rangle=\frac{1}{L}\int_{0}^{L}\phi(r,\theta,z)\,dz$. For the potential, periodic boundary conditions are considered in the $\theta$ and $z$ directions whereas in the radial direction, homogeneous Dirichlet boundary conditions are imposed. See [8], [23] and [26] for more details. For the solution $f$, periodic boundary conditions are imposed in the $\theta,z,v$ directions, whereas in the radial $r$ directions, Dirichlet boundary conditions are imposed. The problem dependent functions $n_{0}$ and $T_{e}$ only depend on $r$.

2 Numerical methods

In this section we will first propose our exponential integrator based numerical method (section 2.1) and then the splitting/semi-Lagrangian based approach (section 2.2). Finally, we discuss the computational cost of both methods (section 2.3).

2.1 Exponential integrator based scheme

2.1.1 Semi-discretization in time

Let us note that the periodic variable $z$ is usually responsible for the most stringent CFL condition. Thus, our goal is to derive a numerical method which eliminates this CFL condition. To that end we write equation (1) as follows

[TABLE]

where

[TABLE]

Next, we perform a Fourier transform in the $z$-direction which gives

[TABLE]

where $\hat{g}=\mathcal{F}(g)=\hat{g}(t,r,\theta,k,v)$ denotes the Fourier transform of an arbitrary function $g(t,r,\theta,z,v)$ in $z$. The associated variable in frequency space is denoted by $k$. Using the variation of constants formula, the exact solution can then be written as

[TABLE]

or equivalently for a time step of size $\Delta t$ from $t^{n}=n\Delta t$ to $t^{n+1}=(n+1)\Delta t$

[TABLE]

To obtain a first order exponential integrator we approximate $f(t^{n}+s)$ by $f(t^{n})$ and introduce the $\varphi_{\ell}$-functions defined by the recurrence relation

[TABLE]

with $\varphi_{0}(z)=e^{z}$. The first three $\varphi_{\ell}$-functions are

[TABLE]

Consequently, we can write the first order exponential Euler method in the following form (denoting by $f^{n}$ an approximation to $f(t^{n})$)

[TABLE]

for $z=-ikv\Delta t$. Note that for $v\to 0$ we recover the explicit Euler method

[TABLE]

which makes it clear that, while we have eliminated the CFL condition in $z$, we still have to satisfy the CFL condition corresponding to $F(f^{n})$.

A significant literature has been dedicated to exponential integrators (see for example [7, 20]). In particular, order conditions have been derived that allow for the construction of higher order methods. In 1987 Strehmel and Weiner [29] introduced an exponential integrator of second order by using rational functions rather than $\varphi$-functions. To the best of our knowledge, the first systematic approach to exponential integrators using $\varphi$-functions was done in [19] by Hochbruck & Ostermann. In the present work we will mostly consider the two-stage second order scheme from [19]

[TABLE]

and the four-stage fourth order scheme, namely the Cox-Matthews scheme [7]

[TABLE]

2.1.2 Phase space approximation

To complete the numerical scheme, one has to detail the phase space approximation. That is, we have to specify how to discretize the explicit part of the exponential integrator in space. In the proposed algorithms, one has to approximate $F(f(t^{n}))$ with respect to the $r$, $\theta$, and $v$ directions. Using the notation of the advection formulation (1), the right hand side is of the form

[TABLE]

where the Poisson bracket is given by $\{\phi^{n},f^{n}\}=\partial_{r}f^{n}\partial_{\theta}\phi^{n}-\partial_{\theta}f^{n}\partial_{r}\phi^{n}$ and is discretized by the second order Arakawa finite difference scheme [1]. The remaining part $\partial_{z}\phi^{n}\partial_{v}f^{n}$ is discretized by a second order upwind scheme. Consequently, we obtain a CFL type condition for the $r,\theta,$ and $v$ directions, but not in the stiffest $z$-direction.

Let us discuss the choice of the finite difference scheme to discretize the Poisson bracket in more detail. The first idea could be to use

[TABLE]

where $D_{r}$ and $D_{\theta}$ are the standard centered difference approximations of the first derivative. The resulting scheme is second order accurate and conserves mass exactly. Unfortunately, the present scheme is not as favorable with respect to momentum and energy conservation. In [1], Arakawa introduced a finite difference method that can simultaneously conserve mass, momentum, and energy. This numerical scheme considers an equal superposition of three parts

[TABLE]

Each of the stencils $J^{++}$, $J^{+x}$, and $J^{x+}$ are mass conservative on their own but only their sum conserves momentum and energy as well (in the case of periodic boundary conditions). Thus, in our case this method is attractive because the only error made in these quantities is due to the time discretization (i.e. due to the exponential integrator employed). Since it is known that exponential integrators do conserve mass (see the proof in [11]) the same is true for the numerical method proposed here. Let us also note that this scheme has been generalized to discontinuous Galerkin type methods and thus higher order variants can be easily constructed [14]. However, since we will focus on time integration in this paper we will not consider such methods further.

Although this is not commonly acknowledged in the literature, even for homogeneous Dirichlet boundary conditions mass is not conserved up to machine precision. Let us investigate this in more detail. We will proceed as follows. First we show that $J^{++}$ conserves mass to machine precision. This is a straightforward but rather tedious calculation. Then we will show why mass conservation is violated for $J^{+x}$. This result will also allow us to suggest a procedure to fix this deficiency. To do so, we introduce a mesh in $r$ and $\theta$: $r_{i}=r_{\min}+ih_{r}$ for $i=0,\dots,N+1$ ($h_{r}=(r_{\max}-r_{\min})/(N+1)$) and $\theta_{j}=jh_{\theta}$ for $j=0,\dots,N$ ($h_{\theta}=2\pi/N$); moreover, we define $\phi_{i,j}$ (resp. $f_{i,j}$) as an approximation of $\phi(r_{i},\theta_{j})$ (resp. of $f(r_{i},\theta_{j})$).

The stencil $J^{++}$ discretizes the Poisson bracket in advection form. We start with the sum over the first term

[TABLE]

where, for simplicity, we have assumed that both the $r$ and the $\theta$-directions are discretized using $N$ grid points. In addition, we have dropped the superscript $n$, which denotes the value at the previous time step, from both $\phi$ and $f$. Let us further note that we do not consider the volume element (arising from cylindrical coordinates) in the calculation. This is valid since the Poisson bracket is multiplied by $1/r$ which immediately annihilates the volume element. Therefore, all results obtained here apply equally to cylindrical and cartesian coordinates.

The homogeneous Dirichlet boundary condition in the $r$-direction is then enforced by setting $f_{0,j}$, $f_{N+1,j}$, $\phi_{0,j}$, and $\phi_{N+1,j}$ (the ghost cells) to zero for all $j$. Using these relations and applying summation by parts in the $r$-direction we obtain

[TABLE]

Now, integration by parts in the $\theta$-direction is straightforward as periodic boundary conditions are imposed in that direction. We obtain

[TABLE]

Since this is exactly the sum over the second term in $J^{++}$, we conclude at once that mass is conserved up to machine precision.

Now, let us consider $J^{+x}$. Since we have periodic boundary conditions in the $\theta$-direction we get

[TABLE]

For the second term we have (focusing only on the summation in $i$)

[TABLE]

where $g=\phi D_{\theta}(f)$. The only local choice to ensure conservation of mass is

[TABLE]

This should not be that surprising. In fact, showing mass conservation for the standard second order finite volume scheme proceeds very similar to the previous argument. In this case the condition above has the natural interpretation that $g$ is zero at the cell interface (in our notation $g_{1/2,j}=\frac{1}{2}(g_{0,j}+g_{1,j})=0,\,\forall j=1,\dots,N$). The above calculation gives the same result if applied to $J^{x+}$.

Note that strictly speaking we have assumed here that $J^{+x}$ and $J^{x+}$ conserve mass separately. This is, of course, not necessarily true as the error in mass committed by $J^{+x}$ could be compensated for by $J^{x+}$. To check this possibility we consider (there is no need to take $J^{++}$ into account as that term conserves mass exactly)

[TABLE]

Enforcing the Dirichlet boundary conditions and performing summation by parts gives (note that, as before, we define $g=\phi D_{\theta}(f)$)

[TABLE]

which is not zero in general.

The previous discussion leaves us with the realization that no matter how we enforce the boundary conditions in at least one of the stencils, mass conservation will be violated at the boundary. We can, however, remedy this deficiency by enforcing different boundary conditions for the different stencils. This will introduce a first order error at the boundary but will ensure mass conservation up to machine precision. We will investigate this issue further in section 3.2, where we will find that this approach also has a positive effect on energy conservation but a negative effect on $L^{2}$ norm conservation.

2.2 Splitting/semi-Lagrangian method

In view of the numerical comparison conducted in the next section, we recall the basic steps of the algorithm used to approximate the drift-kinetic model (1) coupled with the quasi-neutrality equation (1.1) by means of the splitting/semi-Lagrangian method (see [8, 17]).

Assuming $f^{n}$ and $\phi^{n}$ some known approximations to $f(t^{n})$ and $\phi(t^{n})$ (with $t^{n}=n\Delta t$), the main steps of the algorithm to compute $f^{n+1}$ and $\phi^{n+1}$ are

• •

solve the $1D$ advection $\partial_{t}f+v\partial_{z}f=0$ with step size $\Delta t/2$ using a semi-Lagrangian method,

• •

solve the $1D$ advection $\partial_{t}f-\partial_{z}\phi^{n}\partial_{v}f=0$ with step size $\Delta t/2$ using a semi-Lagrangian method to get $f^{\star}$,

• •

solve the quasi-neutrality equation (1.1) with right hand side $\int_{\mathbb{R}}rf^{\star}dv$ to compute the potential $\phi^{\star}$, with FFT in $\theta,z$ and second order finite difference in $r$,

• •

compute the derivatives $(\partial_{r}\phi^{\star},\partial_{\theta}\phi^{\star},\partial_{z}\phi^{\star})$ of $\phi^{\star}$ with cubic splines,

• •

solve the $2D$ advection in $(r,\theta)$ with step size $\Delta t$ using a semi-Lagrangian method,

• •

solve the $1D$ advection $\partial_{t}f-\partial_{z}\phi^{\star}\partial_{v}f=0$ with step size $\Delta t/2$ using a semi-Lagrangian method,

• •

solve the $1D$ advection $\partial_{t}f+v\partial_{z}f=0$ with step size $\Delta t/2$ using a semi-Lagrangian method, to get $f^{n+1}$

• •

solve the quasi-neutrality equation (1.1) with right hand side $\int_{\mathbb{R}}rf^{n+1}dv$, to compute the potential $\phi^{n+1}$, with FFT in $\theta,z$ and second order finite difference in $r$.

Contrary to our method, this approach is not restricted by any CFL condition. However, as we shall see, our method has some advantages from a computational point of view. Moreover, for this drift-kinetic model, the semi-Lagrangian method is not conservative (due to the $2D$ spline interpolation in the $(r,\theta)$ plane used in [8, 17]). Thus the corresponding long time behavior of the numerical solution may become unsuitable.

2.3 Computational cost

For the actual implementation it is advantageous to rewrite the second order exponential integrator given in (4) in the following form

[TABLE]

In this format we see that the method requires four FFTs (compute $f^{n}$ and $k_{1}$ from $\hat{f}^{n}$ and $\hat{k_{1}}$, compute $\mathcal{F}({F(f^{n})})$ and $\mathcal{F}({F(k_{1})})$ from $F(f^{n})$ and $F(k_{1})$).

In addition, we require the evaluation of two right hand sides which consists of a stencil code and the quasi-neutrality equation (1.1). The latter in turn requires two FFTs and the Thomas algorithm (which scales as $\mathcal{O}(N_{r})$, with $N_{r}$ the number of points in the $r$ direction) but only takes place in a three-dimensional subspace. Consequently the main numerical effort is determined by the four one-dimensional FFTs that have to be applied to the entire four-dimensional phase space.

For comparison, the semi-Lagrangian approach requires four one-dimensional semi-Lagrangian steps in one dimension and one semi-Lagrangian step in two dimensions. The solution of the quasi-neutrality equation incurs the same computational effort in both schemes. In [8, 17] the semi-Lagrangian steps are done using cubic spline interpolation. For the constant advection case in the $z$-direction two FFTs could be performed instead of the cubic spline interpolation (note, however, that this is not possible for the remaining semi-Lagrangian steps). It is not entirely straight forward to compare the performance of cubic spline interpolation to FFT based techniques. However, since FFT libraries are available that operate close to the theoretical hardware limit with respect to performance, we might expect that constructing and then evaluating the cubic spline interpolation is at least (which would require that both the construction as well as the evaluation step operate close to the theoretical bandwidth limit) as expensive as two FFTs. Under this assumption the remaining cost of the semi-Lagrangian method results from $2$ one-dimensional advections and a single two-dimensional advection, while the remaining cost of the exponential integrator is due to $3$ evaluations of $F$. We note that the evaluation of $F$ is a standard stencil code for which it is well known that performance comparable to the theoretical peak performance of the hardware can be achieved on virtually all architectures. Consequently, the exponential integrator approach introduced in this paper has a clear advantage from a computational point of view.

In addition, the splitting approach still has to tackle the two-dimensional semi-Lagrangian step which is a further drawback of this method (both with respect to efficiency but also with respect to ease of implementation).

In Table 1 we directly compare the runtime of the semi-Lagrangian approach to the second order exponential integrator proposed in this paper (the differences between the perturbation formulation and the direct formulation as well as further details about our numerical method are discussed in the next section). Even though we use the same computer (dual socket Intel Intel Xeon E5-2630 v3 with 32GB of main memory) and a sequential run for all the simulations, these numbers have to be considered with care, since the two codes are written in different languages (the semi-Lagrangian code uses the library SeLaLib and is written in Fortran while the integrator proposed here is implemented in C). However, we observe that the proposed method is more than four times faster compared to the semi-Lagrangian approach.

3 Numerical simulations

In this section, we detail the physical parameters of the considered test case. First, the initial value is given by

[TABLE]

where the equilibrium function is

[TABLE]

The radial profiles $\{T_{i},T_{e},n_{0}\}$ have the analytic expressions

[TABLE]

where the constants are

[TABLE]

Finally, we consider the parameters of [6] (MEDIUM case)

[TABLE]

and use a $v$-range of

[TABLE]

For the direct formulation, we use the boundary conditions

[TABLE]

Note that these are not homogeneous Dirichlet boundary conditions. It is well known (and supported by our numerical experiments) that the Arakawa scheme works better for homogeneous boundary conditions. In addition to the direct formulation, we therefore also introduce a so-called perturbation formulation. First, we note that the equilibrium function $f_{\rm eq}$ defined in (7) is a steady state for our problem. We therefore divide $f$ into

[TABLE]

With this formulation, our problem (1) becomes

[TABLE]

where $E_{\theta}=-\partial_{\theta}\phi$, $E_{r}=-\partial_{r}\phi$ and $E_{z}=-\partial_{z}\phi$. Expanding the various terms we obtain

[TABLE]

which can be written as

[TABLE]

Note that the equation is very similar to equation (1). We, however, have obtained two additional source terms, which depend on the equilibrium distribution $f_{eq}$ as well as on the electric field. Furthermore, the right hand side of the quasi-neutrality equation (1.1) becomes

[TABLE]

Due to the similarity of the direct formulation and the perturbation formulation, the same code can be used for both by simply exchanging the right hand side of the quasi-neutrality equation, changing the boundary conditions, and adding the appropriate source terms. Thus, to implement the exponential integrator we consider the following equation

[TABLE]

with

[TABLE]

and proceed as in section 2.1.1 (with $F$ replaced by $F_{pert}$). The space discretization of the source terms can be done either analytically or using a numerical approximation. In our implementation we have used standard centered differences. The Arakawa scheme that is used to discretize $F(\delta f)$ now employs homogeneous Dirichlet boundary conditions for $\delta f$ in the $r$-direction.

3.1 Numerical results

In this section, we compare the standard BSL method introduced in section 2.2 to the newly proposed method. To do that, we consider the $4D$ drift-kinetic model in cylindrical geometry (1)-(1.1).

We are interested in the time history of the electric energy

[TABLE]

where $r_{p}=(r_{\rm max}+r_{\rm min})/2$, the total mass

[TABLE]

the $L^{2}$ norm

[TABLE]

and the total energy

[TABLE]

The quantity in (8) is known to exponentially increase with time (see, for example, [6, 17]). The growth rate can be computed by solving the linearized drift-kinetic model (1)-(1.1). This is done by performing a Fourier transform in $(z,\theta)$ and a Laplace transform in the time variable $t$. The resulting dispersion relation is derived in [6]. We have solved this relation using a computer algebra system and obtained the root with the largest imaginary part. This is precisely the (analytically derived) growth rate.

The total mass, $L^{2}$ norm, and energy are physically conserved quantities in the analytic model. They are therefore an important measure of the physical accuracy of the numerical method.

The following numerical results are shown for two different grid sizes for $(r,\theta,z,v)$, where the discretization of $32\times 32\times 32\times 64$ is referred to as the coarse grid and the discretization of $64\times 64\times 64\times 128$ is referred to as the fine grid. For the coarse grid, we integrate up to a final time $T=8000$, while the fine grid is integrated up to $T=6000$. Furthermore, we compare the perturbation formulation with the direct one.

3.1.1 Exponential integrator of order 2

In this section, we discuss the numerical results of the second order exponential integrator described by (4). In Fig. 1, we plot the time history of (8) obtained by the perturbation formulation as well as the direct formulation of the exponential integrator for both the coarse and the fine grid. The results are plotted using the maximum possible time step $\Delta t$ for each version. We can see, that all versions reproduce the linear part very well and are close to the analytic growth rate derived in [6].

We observe, as expected, that the fine grid requires a smaller time step to maintain stability. However, if we examine the CFL condition a bit more thoroughly we can see that we still can take a relatively large step size compared to the CFL condition in the $z$-direction. The CFL condition is given by

[TABLE]

where the speeds associated with these terms are

[TABLE]

Due to the non-linearity of the problem, the CFL condition depends on the numerical solution. It is therefore not straightforward to predict the CFL condition based on theoretical considerations or by numerical data gathered on a coarser grid (see, for example, [15]). Ideally, an adaptive time stepping scheme would be employed that automatically selects an appropriate step size (such as using the strategy in [18] with an embedded exponential integrator [20]). As we will see in this section this would be extremely beneficial in the linear regime (i.e. up to approximately $t=3000$) as the nonlinear part of the simulation dictates the time step size. However, since the semi-Lagrangian implementation uses a constant step size we will do the same here. This enables us to perform a more direct comparison with the semi-Lagrangian scheme (in particular, concerning the invariants). However, we consider developing an adaptive version of the present algorithm as future work.

We did evaluate the CFL number during the simulation. For the coarse grid (henceforth denoted with a subscript $c$), the simulation ran until the final time $T=8000$. The computed minimum CFL values for the second order exponential scheme (since these values differ only slightly between the direct and perturbation formulation and change only slightly for reasonable $\Delta t$, we exclusively report the overall minimum) are

[TABLE]

Using the maximum value $v=7.32$ and the domain length $L\approx 1506.8$ of our simulation, the CFL condition in the $z$-direction is

[TABLE]

which as the smallest value would be the most restrictive CFL condition for the method. However, with the exponential integrator proposed in this paper we are able to take time steps with $\Delta t=11>\text{CFL}_{z,c}$ for the perturbation formulation, while still obtaining physically reasonable results. It should be duly noted, however, that we still have to respect the CFL condition in the remaining directions. That is, the exponential integrators have to satisfy $\Delta t<\min(\text{CFL}_{r,c},\text{CFL}_{\theta,c},\text{CFL}_{v,c})$. The time history of the CFL number, defined as $\min\left(\text{CFL}_{r,c},\text{CFL}_{\theta,c},\text{CFL}_{v,c}\right)$, computed for both the coarse and the fine grid is shown in Fig. 2.

The coarse grid corresponds exactly to the problem considered in [8]. In that work a time step size of $\Delta t=8$ is used. That is, we are able to run the simulation with similar (or even larger) time steps while using a significantly cheaper method (see the discussion in section 2.3).

For the fine grid, the simulation ran until $T=6000$ and resulted in the values

[TABLE]

while the computed CFL condition in the $z$-direction is

[TABLE]

which again is the most restrictive condition. However, in Fig. 1, we see that once again the time step can be chosen larger than that value and our method still produces a stable run. In Fig. 3, we plot the error of the total mass defined in (9) and compare our method with the results from the semi-Lagrangian method. We see that the mass in both the perturbation method as well as the direct formulation is much better preserved than for the semi-Lagrangian approach.

Similarly, the error of the $L^{2}$ norm defined in (10) is plotted in Fig. 4. Here, we can observe, that the error is smaller or equal compared to the semi-Lagrangian method.

In Fig. 5, the error in the total energy (as defined in (3.1)) is shown. Our method performs approximately an order of magnitude better compared to the semi-Lagrangian scheme.

3.1.2 Methods of higher order

In this section, we discuss the numerical results for the method of order 4 described in (5). In Fig. 6, we can see that all variants of our method, i.e, the perturbation formulation and the direct formulation for both the coarse and the fine grid reproduce the linear growth rate of (8) (as compared to the analytic result given in [6]). This behavior is similar compared to the second order method.

Using numerical simulations to determine the maximal allowed step size, we observe that the time step size for the fourth order method is significantly less stringent compared to the method of order 2. The maximum allowed time steps for all configurations are listed in Table 2. In particular, for the fourth order method in the perturbation formulation we can take significantly larger time steps compared to the second order method. Since all methods are able to resolve the growth phase very well, the ability to take larger time steps is the major advantage of the fourth order method in the linear regime. In addition, employing higher order schemes is useful as they introduce less diffusion and generally provide more accurate results.

In Fig. 7, we show the error in mass for the fourth order method. The behavior is similar to that of the second order method. In particular, we see a multiple orders of magnitude improvement compared to the semi-Lagrangian approach. This is true even for the significantly larger time steps we are now able to take.

The $L^{2}$ error shown in Fig. 8 and the energy error shown in Fig. 9 depict similar behavior compared to the second order exponential integrator. In particular, for the $L^{2}$ error we see a small improvement compared to the semi-Lagrangian approach, while for the energy error we observe an improvement in the error by approximately one order of magnitude.

3.2 Mass conservation up to machine precision

As discussed in section 2.1.2 and observed in the previous section, the Arakawa space discretization is not mass conservative up to machine precision (even for homogeneous Dirichlet boundary conditions). Although it should be noted that conservation of mass is still very good (up to $10^{-8}$ in the simulations conducted in the previous section) and, in particular, significantly better than conservation of both momentum and energy. Nevertheless, in some situations it is still beneficial to consider a numerical method that conserves mass up to machine precision.

Thus, in the following we will consider the performance of the numerical method suggested in section 2.1.2 which conserves mass to machine precision. The numerical results are shown in Fig. 10. As expected, we observe conservation of mass up to machine precision which gives a significant advantage compared to standard homogeneous boundary conditions in Arakawa’s finite difference method. In addition, we observe a small decrease in the energy error at the expense of a significant increase in the error of the $L^{2}$ norm. We have also conducted the corresponding simulation with the fourth order exponential integrator. However, since, as expected, the results are almost identical we do not show them here.

4 Conclusion & Outlook

In this paper we have demonstrated that using exponential integrators is a viable approach for numerically integrating the drift-kinetic equations in time. Combining this with a suitable space discretization (such as Arakawa’s method) yields a numerical method with improved conservation properties and significantly reduced computational effort that is still able to take large time steps compared to the most stringent CFL condition. Furthermore, the numerical method proposed in this paper could conceivably be applied to more complicated gyrokinetic models. We consider this as future work.

In this paper we have exclusively used Fourier techniques in order to solve the linear part of the drift-kinetic equation. However, in principle, any semi-Lagrangian scheme can be used in its place (in the spirit of [23], for example). For example, using a (one-dimensional) cubic spline or discontinuous Galerkin approach would still result in a numerical scheme with mass conservation up to machine precision. In particular the latter, i.e. using a semi-Lagrangian discontinuous Galerkin approach, might be beneficial in the context of high performance computing systems, where parallelization to a larger number of cores is essential to obtain good performance. We consider this as future work.

Acknowledgements

This paper is based upon work supported by the VSC Research Center funded by the Austrian Federal Ministry of Science, Research and Economy (bmwfw). N.C. is supported by the Enabling Research EUROFusion project CfP-WP14-ER- 01/IPP-03 and the IPL FRATRES.