Numerically stable method for kinetic electrons in ... - Semantic Scholar
Journal of Computational Physics 239 (2013) 22–29
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Numerically stable method for kinetic electrons in gyrokinetic particle-in-cell simulation of toroidal plasmas T. Korpilo a,⇑, J.A. Heikkinen b, S.J. Janhunen a, T.P. Kiviniemi a, S. Leerink a, F. Ogando c a
Euratom-Tekes Association, Aalto University, P.O. Box 14100, FI-00076 Aalto, Finland Euratom-Tekes Association, VTT, P.O. Box 1000, FI-02044 VTT, Finland c Universidad Nacional de Educacion a Distancia, C/Juan del Rosal, 12, 28040 Madrid, Spain b
a r t i c l e
i n f o
Article history: Received 21 March 2012 Received in revised form 5 October 2012 Accepted 22 December 2012 Available online 10 January 2013 Keywords: Kinetic electrons Gyrokinetic particle simulation Fusion plasma physics
a b s t r a c t The direct implicit method with a second-order implicit integration scheme is formulated for and applied to the electron parallel nonlinearity in global electrostatic gyrokinetic particle-in-cell simulations of toroidal fusion plasmas. The method shows improved numerical accuracy and stability properties compared to the direct implicit method with a first-order integration scheme. The conservation of total energy and toroidal angular momentum are analyzed by both techniques and the results are presented. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction The computational power provided by modern supercomputers is making predictive particle simulations of global plasma transport possible for magnetized fusion devices like tokamaks. Although the full picture of particle and heat transport is still beyond the reach of such simulations, many of the mechanisms responsible for electrostatic drift-wave turbulence are included by treating ions gyrokinetically and electrons drift-kinetically. Compared to full Lorentz force ions and electrons, these approximations reduce numerical noise and relax significantly the strict limitations for temporal and spatial resolution while retaining finite ion Larmor radius effects relevant for plasma turbulence. Kinetic effects, the electron parallel flow in particular, are important also for the formation and evolution of the parallel current and for electromagnetic fluctuations. Thus, in first-principle simulations of low-frequency turbulence, the full ion distribution is often solved gyrokinetically while electrons are modelled drift-kinetically. In addition to the full ion distribution, the proper description of turbulent transport necessitates also resolving the full electron distribution. The simplest solution is to push drift-kinetic electron motion explicitly. In electrostatic simulations, the time step is then constrained by both the accuracy condition set by the electron transit time along the magnetic field and the stability condition referred in the literature to being imposed by the unstable xh mode [1]. A way of avoiding numerical instabilities related to the explicit integration of electrons is the direct implicit method for the electron parallel nonlinearity [2]. The method makes no assumptions about the electron distribution function and is based on treating the electron acceleration by the parallel electric field (the electron parallel nonlinearity) implicitly while the remainder of the driftkinetic motion is pushed explicitly. The time step is limited only by the accuracy conditions. The direct implicit method is employed in an approximate fashion in Ref. [3] where the implicit response to the electron parallel nonlinearity is considered in the form of an analytic differential operator, thus including assumptions about the distribution function. There are ⇑ Corresponding author. Tel.: +358 50 4327580. E-mail address: tuomas.korpilo@aalto.fi (T. Korpilo). 0021-9991/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.12.033
T. Korpilo et al. / Journal of Computational Physics 239 (2013) 22–29
23
also other approximate models such as the fluid-kinetic hybrid scheme [4], the split-weight scheme [5,6] and the kinetic closure scheme [7,8] that are introduced for adding electron kinetic effects and electromagnetism into turbulence simulations. The main motivation of such schemes is to ease the strict time step constraint imposed by explicit electrons and thus reduce computational burden. In the electrostatic limit, electrons can also be modelled adiabatically, i.e. as a Maxwell–Boltzmann response, without resolving the distribution function. However, this heavy approximation eliminates net turbulent particle transport (due to the phasing between electron density and electric potential) and electron thermal transport but circumvents any numerical restriction set by the electrons. In the present work, the direct implicit method for the electron parallel nonlinearity is presented in detail in the electrostatic limit for tokamak plasma configurations. The value of the scheme lies in its inherent numerical stability for resolving the full electron distribution kinetically but it also coheres with the implicit construction of the gyrokinetic Poisson equation from the ion polarization drift motion that is used in the ELMFIRE code. The direct implicit method is originally based on the works [9–11] where the method is introduced and employed in particle-in-cell (PIC) simulations. Here, the implementation of the method for the electron parallel nonlinearity is shown for global gyrokinetic PIC drift-wave turbulence simulations. The present work will also extend the implicit method already demonstrated in Ref. [2] by considering a second-order integration scheme for the electron parallel nonlinearity. This more accurate integration scheme leads to improved energy conservation, through reduced numerical heating, and enables stable simulations when a continuous finite difference method is used for the electric field appearing in the parallel acceleration. The paper is organized as follows. In Section 2, the ELMFIRE code is briefly described to the extent required for presenting the implementation of the implicit method for electrons. In Section 3, the actual implementation as well as the new and the old integrator for the parallel nonlinearity are shown. In Section 4, the integrators are compared to each other with respect to the accuracy of the conservation of total energy and toroidal angular momentum observed in a turbulent simulation. A short summary is drawn in Section 5. 2. Code description The implicit method for simulating kinetic electrons is described here for ELMFIRE [2], for a global PIC code that advances full 5D distributions of gyrokinetic ions and drift-kinetic electrons under 3D electrostatic potential in a quasi-toroidal (axisymmetric) magnetic equilibrium of tokamak plasmas. In the code, in accordance with the gyrokinetic formalism in Ref. [12], the polarization drift is included in the gyrokinetic equations of motion, leaving only short-wavelength trembling to affect the gyroradius, and the gyrokinetic Poisson equation is reduced to the quasi-neutrality condition. In order to obtain the quasi-neutrality, the ion polarization drift and the electron acceleration by the parallel electric field are integrated implicitly in time, separately from the rest of the terms in the equations of motion that are integrated explicitly. The charge density responses to the ion polarization drift and the electron parallel acceleration are constructed by the direct implicit method and are separated from the quasi-neutral charge densities. The remaining charge densities comprise the charge separation caused by the explicit particle motion. The electrostatic potential is solved from the resulting linear system, after which the actual implicit motion is performed, neutralizing the plasma. As electrons are drift-kinetic in ELMFIRE, their motion is described by the standard guiding center motion. The equations for advancing the guiding center coordinates of a particle in time are given for three position coordinates, R, and parallel velocity, v k , by
! ml v k Bb þ rbþb rB þ r/ q q mv B q v_ k ¼ b þ k r b lrB þ r/ B m qB 1 R_ ¼ B
mv 2k
ð1Þ
together with l_ ¼ 0. Here, the magnetic moment, l, is defined as a function of the perpendicular velocity as v 2? ¼ 2Bl, and B ¼ B þ ðmv k =qÞb r b. In addition, / denotes the electric potential, B the strength of the magnetic field and b the unit vector of the magnetic field, and q and m are the charge and mass of the particle, respectively. The term in the electron motion used for obtaining the quasi-neutrality is now the parallel acceleration ak ¼ ðqB=mB Þb r/ in the equation of v_ k . This part of the electron motion is therefore integrated implicitly, while the rest of the motion is followed explicitly. For clarity and later purposes, we separate the explicit and implicit accelerations in Eq. (1) and rewrite the equation for the electron ~_ k þ ak . parallel acceleration as v_ k ¼ v The quasi-neutrality of the plasma and the concomitant electric potential are solved in a structured uniform grid set up on the quasi-ballooning coordinates, which are again built on the Boozer coordinates of the quasi-toroidal magnetic equilibrium. This equilibrium has concentric circular flux surfaces and thus the poloidal flux, Wp , can be substituted by the radial coordinate, r, in these coordinate systems. The charge density is sampled from electron guiding center positions and from four gyro-ring points of ions to the quasi-ballooning grid by trilinear interpolation by neglecting the non-orthogonality but retaining the radial shear of the coordinate system. Consequently, the charge sampling equals the linear cloud-in-cell (CIC) technique; the center of the cloud extends over radial shells by constant poloidal Boozer angle. The gyrocenter motion of ions and the guiding center motion of electrons, as expressed in the Boozer coordinates, are followed by the classical fourth-order Runge–Kutta scheme, except for the electron and ion acceleration by the parallel electric field and the ion polarization drift. The electron parallel acceleration ak is integrated in time by a second-order implicit scheme (to be discussed in
24
T. Korpilo et al. / Journal of Computational Physics 239 (2013) 22–29
the next Section) and, for symmetry reasons, the ion parallel acceleration is integrated in the same way. The integrator for the ion polarization drift is the first-order backward Euler method. Collisions among ions and electrons are modelled in the drift-kinetic limit, thus eliminating classical diffusion across the magnetic field, and are evaluated with the momentum and energy conserving binary collision model [13]. The primitive velocity vectors, in terms of which the collision model is cast, are obtained from the guiding center (gyrocenter) velocity components v k and v ? by decomposing v ? into its two orthogonal components with a random gyroangle. The covariant components of the electric field present in the equations of motion are interpolated from the quasi-ballooning grid to electron guiding center positions and via gyroaveraging to ion gyrocenter positions by two methods. First, the electric field for the E B drift motion is evaluated by a momentum conserving scheme [14]. Second, the electric field for the parallel acceleration (ion and electron) and ion polarization drift is calculated by an energy conserving scheme which P is expressed generally as EðxÞ ¼ r/ðxÞ ¼ /k rSðx xk Þ in the electrostatic case. Here, the subscript k indexes the grid points surrounding the particle, and S is a 3D shape function determining the interpolation scheme between the particle at x and the grid points at xk . As in the charge sampling, the shape function applied here is the trilinear interpolation, except that a second-order quadratic spline [15] is used in the direction of the derivative. The derivative of this spline equals the derivative of the second-order Lagrange interpolation. However, it is also possible to use the fully trilinear shape function in the electric field for the parallel acceleration. The spline utilizes three potential values from the grid (the three-point method) and produces a continuous electric field, whereas the linear interpolation uses two potential values from the grid (the two-point method) and gives a piecewise constant electric field.
3. Direct implicit method for kinetic electrons Electron parallel dynamics determines normally the largest time step for a gyrokinetic simulation. In order to extend the time step limitation near the electron parallel transit time in a stable way, the electron acceleration by the parallel electric field needs to be treated implicitly. In that case, the time step has to satisfy the accuracy condition given as kk v te Dt K 1, or Dt K Dz=3v te . Here, kk denotes the largest wavenumber of the parallel electric field, v te the electron thermal velocity, Dt the time step and Dz the metric grid spacing in the parallel direction. Another restriction for the time step is imposed by the particle motion perpendicular to the magnetic field, and in particular by the E B drift motion. The concomitant accuracy condition k? v EB Dt K 1 is, however, usually less stringent than the accuracy condition set by the electron parallel motion. An important factor in the direct implicit method is the choice for the integrator by which the implicit parallel acceleration, ak , is to be pushed. The only restrictions for the integrator are, first, that the parallel displacement resulting from the acceleration is linear with respect to the electric potential. This condition guarantees the solvability of the matrix equation arising in the implicit method. Second, to make physical sense, the linear integrator should advance the parallel coordinate and velocity in a consistent way. It is also worth noticing that the implicit parallel motion, together with the implicit ion polarization drift motion, gives a new potential that is used in the explicit particle push in the next time step. Therefore, the hybrid implicit-explicit integration scheme divides the time integration of the equations of motion on two consecutive time steps making the implicit particle motion an a priori estimation for the next time step in the sense of numerical accuracy of the integration. In this work, we consider the Velocity Verlet algorithm, a variant of the Leapfrog method, in an implicit form for the time integration of the electron parallel acceleration ak . The algorithm is the same as the undamped limit of the adjustable implicit particle mover proposed for implicit PIC simulations in Ref. [16]. This algorithm accurate to the second order in the time step fulfils the requirements for the integrator and is described in the guiding center coordinates as
1 Rnþ1 ¼ Rn þ R_ ðzn ; t n ÞDt þ ak ðznþ1 ; t nþ1 ÞDt 2 bðRnþ1 Þ 2 a ðz ; t Þ þ ak ðzn ; t n Þ v k;nþ1 ¼ v k;n þ v~_ k ðzn ; tn ÞDt þ k nþ1 nþ1 Dt 2
ð2Þ
~_ k , are included for completeness. The time where the guiding center equations of motion integrated explicitly in time, R_ and v level of the integration is given by the subscript n, and z represents the coordinates ðR; v k Þ. In the algorithm, the implicit acceleration ak is considered not only to change the particle parallel velocity but also to cause a small parallel shift, dR ¼ ðak Dt 2 =2Þb, to the particle position through the term ðB=B Þv k b present in the guiding center equations of motion (the contribution of B=B is omitted). The direct influence of the implicit acceleration is neglected on the other v k -dependent terms in the equations of motion. In ELMFIRE, electrons are not pushed exactly as described above but, as already mentioned, the explicit motion is integrated by the classical Runge–Kutta scheme unlike suggested by the algorithm. This is because the E B drift motion is not integrated accurately in time by a low-order scheme. In addition, the acceleration ak ðznþ1 ; tnþ1 Þ and ~ nþ1 ¼ Rn þ Rðz _ n ; t n ÞDt, instead of the magnetic field bðRnþ1 Þ are evaluated at the end position of the explicit motion, i.e. at R Rnþ1 since dR is short compared to any resulting parameter change and the size of a grid cell. Similarly, the parallel velocity v~ k;nþ1 ¼ v k;n þ v~_ k ðzn ; tn ÞDt is employed in the integration of the implicit motion, as ak has only a weak dependence on v k through B . Consequently, a costly iteration procedure is avoided in the algorithm. We will refer to the algorithm in Eq. (2) as the Velocity Verlet algorithm.
T. Korpilo et al. / Journal of Computational Physics 239 (2013) 22–29
25
After the selection of the integrator for the implicit particle movement, the direct implicit method is straightforward. The starting point for the method is the requirement for the plasma quasi-neutrality and, following the ideas presented in Ref. [9], the quasi-neutrality is enforced by the implicit particle movement. This means, with the notation of Eq. (2), that the implicit electron parallel acceleration is integrated at the end of a time step, in order to solve the unknown electric potential /ðtnþ1 Þ. Let us denote the quasi-neutral electron charge density discretized on the quasi-ballooning grid at the time level tnþ1 P P by qk ðxnþ1 Þ. The charge density is evaluated on the grid as qk ¼ qk;p ¼ qp Sðxp xk Þ=V k , where qp denotes the total charge of the pth electron cloud, V k the volume of the kth grid cell and x the charge sampling coordinates, i.e. the radial coordinate r, the Boozer poloidal angle h and the field-aligned (quasi-ballooning) angle f. Now, Taylor expanding the shape function S up to the first order with respect to the charge sampling coordinates of the pth electron at the end point of the explicit push, ~nþ1 , one obtains x
qk ðxnþ1 Þ ¼ qk ðx~nþ1 Þ þ
X @ qk;p ðx ~nþ1 Þ p
@fp
dfp ð~znþ1 ; t nþ1 Þ:
ð3Þ
The second term denotes here the linear increment to the former charge density resulting from the implicit electron parallel acceleration and thus dfp ¼ dRp rf (the small discrepancy between the orientation of b and the quasi-ballooning parallel vector ef is neglected in this paper for convenience). With the trilinear shape function, this linear increment to the initial electron charge density is exact for all the electrons not crossing any grid points in the parallel direction during the implicit shift. The increment is now depending linearly on the electric potential /k ðt nþ1 Þ through the shift dfp , and together with the contribution from the similarly constructed implicit ion polarization charge density, the quasi-neutrality condition reduces to a linear matrix equation for the electric potential, from which the potential is solved by conventional means. With the solution for /k ðt nþ1 Þ, the parallel shift of each electron, dfp , is solved, and the electrons are moved accordingly to ~ nþ1 þ dR. At the same time, the electron parallel velocity is updated to the value of their quasi-neutral positions Rnþ1 ¼ R v k;nþ1 . At the next time step, the time integration of the electron explicit motion is started from the point Rnþ1 ; v k;nþ1 (in the absence of binary collisions) by using the solved electric potential /k ðt nþ1 Þ. The explicit push reproduces a charge separation between ions and electrons, and the implicit step is repeated. The Velocity Verlet algorithm considers the electron parallel acceleration ak up to the second order in the time step for both the position and the velocity coordinate [16]. In Ref. [2], the implicit integrator for the kinetic electrons differs from the presented algorithm in that the change of the parallel velocity is integrated from the parallel acceleration by the first-order backward Euler method, i.e. v k;nþ1 ¼ v k;n þ ak ðznþ1 ; tnþ1 ÞDt. The concomitant parallel shift is obtained from the Taylor series expansion Rk;n ¼ Rk;nþ1 v k;nþ1 Dt þ ak ðznþ1 ; tnþ1 ÞDt2 =2 by substituting the expression for v k;nþ1 in it (the term B=B is again ignored). The resulting algorithm, with the inclusion of the explicitly integrated guiding center motion, is the same as the Velocity Verlet algorithm in Eq. (2), except that the velocity change due to the implicit parallel acceleration is evaluated as dv k ¼ ak ðznþ1 ; tnþ1 ÞDt instead of dv k ¼ ½ak ðznþ1 ; tnþ1 Þ þ ak ðzn ; tn ÞDt=2. In the next Section, the difference between the first-order and the second-order accurate dv k is considered relative to a numerical simulation. 4. Simulation results In the standard version of ELMFIRE, the change of the parallel velocity due to the implicit electron acceleration has been integrated in time by the backward Euler method. Moreover, the covariant component of the electric field in the parallel acceleration ak has been evaluated by the two-point method. To compare this Euler scheme to the Velocity Verlet algorithm by using both the two-point (discontinuous) and the three-point (continuous) method for the electric field, a series of simulations are performed with identical simulation parameters. In the simulations, the main parameters of deuterium plasma are similar to FT-2 tokamak discharges: minor radius of 8 cm, major radius of 55 cm, magnetic field of 2.2 T (on the axis) and total plasma current of 55 kA. The radial current density profile is linear. Further, the initial radial density profile is parabolic with the maximum value of 5:0 1019 m3 at the magnetic axis and the minimum value of 0:1 1019 m3 at the separatrix. The initial radial electron (and ion) temperature profile is also parabolic and has a maximum value of 120 eV at the axis and a minimum value of 10 eV at the separatrix. Together with the time step of 30 ns, the parallel accuracy condition is satisfied in the simulations, i.e. the parameter v te Dt=Dz lies between 0.45 and 0.15. The simulation domain is limited between the normalized minor radius of 0.25 and 1.0 and is divided into 50 radial shells each of which consists of 200 grid cells in poloidal direction and 8 grid cells in toroidal direction. The average number of electrons (and ions) per grid cell is 1100. The comparison between the integration schemes employed in the implicit electron acceleration is shown in Figs. 1–4, including the particle number and the time step scan. The comparison is based on the accuracy of the conservation of total energy and toroidal angular momentum, as electron parallel dynamics is vulnerable to the numerical heating/cooling and as the integration of the parallel acceleration ak can cause a numerical source term for toroidal angular momentum. In the simulations, total energy and toroidal angular momentum are evaluated in a diagnostic domain that is composed of the simulation domain excluding the five first and five last radial shells. This exclusion eliminates boundary effects from contributing to the conserving variables, and the only physical sinks and sources in the diagnostic domain to be considered are the energy flux and the toroidal angular momentum transport, carried by the particles, through the domain boundaries. The expressions for the conserved total energy K and toroidal angular momentum L, consistent with the gyrokinetic formalism used in ELMFIRE, are given in Refs. [12,14]. In the figures, the initial value of total energy, Kðt ¼ 0Þ, and toroidal angular momentum,
26
T. Korpilo et al. / Journal of Computational Physics 239 (2013) 22–29
Lðt ¼ 0Þ, in the diagnostic domain are 28 J and 2:4 106 m2 kg/s, respectively. The latter value comprises the ion banana current, arisen as ions are initialized on their neoclassical (collisionless) orbits. In the following, the Velocity Verlet algorithm is shown to provide a numerically stable method to obtain acceptable toroidal angular momentum and energy conservation. As shown in Fig. 1, both the Velocity Verlet algorithm and the three-point method for the electric field improve the conservation of total energy, compared to the Euler scheme with the two-point method for the electric field, while an opposite trend is seen in the conservation of total toroidal angular momentum. The improved energy conservation is explained by reduced numerical heating/cooling of electrons. The increased error in the conservation of toroidal angular momentum, on the other hand, is explained by increased asymmetry between the integration of the parallel acceleration ak of ions and electrons, enabled by the inhomogeneous charge density in the space between the grid points. The inhomogeneity exists even if the charge density sampled at the grid points is neutral. These results indicate that the numerical heating/cooling of electrons plays a dominant role in the time evolution of total energy and that the asymmetric integration of the electron and ion parallel acceleration affects notably the time evolution of total toroidal angular momentum. Fig. 2 illustrates a severe numerical instability arising when the three-point method for the electric field is used together with the Euler scheme. In other words, electrons experience violent numerical heating which sets in after 10th time step and as a result of which electron energy content doubles in a very short time period of 1 ls. The instability has been observed in all simulations, regardless of the simulation parameters. A similar type of instability results both with the Euler scheme and with the Velocity Verlet algorithm if a momentum conserving interpolation scheme is utilized in evaluating the electric field P for the parallel acceleration. In other words, if the electric field is calculated from the expression EðxÞ ¼ Ek Sðx xk Þ by employing the trilinear shape function and the central difference method for the covariant components of the electric field at the grid points.
1.06
1.08 Euler, 2−point
Verlet, 3−point
Verlet, 2−point
Verlet, 2−point
Verlet, 3−point
Euler, 2−point
1.06
L(t) / L(t=0)
K(t) / K(t=0)
1.04
1.04
1.02 1.02
1
500
1000
1500
2000
2500
1
3000
t / 30 ns
500
1000
1500
2000
2500
3000
t / 30 ns
Fig. 1. Accuracy of energy (left) and toroidal angular momentum (right) conservation resulting from integrating the electron response to the parallel electric field by Velocity Verlet and Euler schemes. The dashed line is reflected around x-axis for illustrative purposes in the picture on the left.
1.5 Euler, 3−point Euler, 2−point
1.4
K(t) / K(t=0)
1.3 1.2 1.1 1 0.9
10
20
30
40
50
t / 30 ns Fig. 2. The Euler scheme with the 3-point method for the electric field is unstable.
27
T. Korpilo et al. / Journal of Computational Physics 239 (2013) 22–29
The accuracy of the conservation of total energy and toroidal angular momentum is presented in Fig. 3 in terms of the number of simulated particles. For this, the simulation with the Velocity Verlet algorithm and with the three-point method for the electric field is repeated with twice (2 N), four times (4 N) and eight times (8 N) the particle number per grid cell. The particle numbers of 1 N, 2 N, 4 N and 8 N translate into the noise fluctuation levels dn=n of 0.84%, 0.57%, 0.40% and 0.28%, respectively, while the observed turbulence fluctuation level is 1–5% depending on the radial position. The simulation results reveal that the energy conservation improves according to 1=N and that the conservation of toroidal angular momentum improves even faster. However, the fluctuating time evolution of toroidal angular momentum in the 8 N case indicates that the error in the conservation of toroidal angular momentum is not vanishing only by increasing the number of particles. A potential reason for this, namely the numerical accuracy of the diagnostics, is discussed in Ref. [14]. Nevertheless, 0.25% error in the energy conservation and under 0.4% error in the toroidal angular momentum conservation is observed in the 8 N case during 90 ls or 3000 time steps. The Euler scheme and the Velocity Verlet algorithm with the two-point method for the electric field, in turn, improve the conservation of toroidal angular momentum in terms of the particle number as the Velocity Verlet algorithm with the three-point method, while in energy conservation the picture is different. In other words, the Velocity Verlet algorithm with the two-point method improves energy conservation slower than 1=N, and the Euler scheme shows practically no improvement. In the simulations, the growth of turbulent modes starts around the 500th time step after which turbulence is a strong candidate in explaining the changes in the time evolutions of total energy and toroidal angular momentum. The effect of the time step size on the accuracy of the conservation of total energy and toroidal angular momentum is presented in Fig. 4 for two simulation cases, one employing the Euler scheme and the other the Velocity Verlet algorithm with the two-point method for the electric field. In addition to the base step size, one-half and twice the time step size is
1.08 Verlet, 3−point, 1N
1.04
Verlet, 3−point, 1N
Verlet, 3−point, 2N
Verlet, 3−point, 2N
Verlet, 3−point, 4N Verlet, 3−point, 8N Verlet, 2−point, 4N
L(t) / L(t=0)
K(t) / K(t=0)
1.03
Verlet, 3−point, 4N
1.06
Euler, 2−point, 4N 1.02
Verlet, 3−point, 8N Verlet, 2−point, 4N Euler, 2−point, 4N
1.04
1.02 1.01 1 1
500
1000
1500
2000
2500
3000
500
1000
t / 30 ns
1500
2000
2500
3000
t / 30 ns
Fig. 3. Accuracy of energy (left) and toroidal angular momentum (right) conservation with respect to particle number. The crossed line is reflected around x-axis for illustrative purposes in the picture on the left. N equals 1100 particles per grid cell.
1.06
Verlet, 2−point, 2Δt
Verlet, 2−point, Δt
1.06
Verlet, 2−point, Δt/2
1.02 1 0.98 0.96
Verlet, 2−point, Δt Verlet, 2−point, Δt/2
1.05
L(t) / L(t=0)
K(t) / K(t=0)
1.04
Verlet, 2−point, 2Δt
Euler, 2−point, 2Δt Euler, 2−point, Δt
1.04
Euler, 2−point, Δt/2
1.03 1.02
Euler, 2−point, 2Δt 1.01
Euler, 2−point, Δt Euler, 2−point, Δt/2
0.94
15
30
45
t (μs)
60
75
90
1
15
30
45
60
75
t (μs)
Fig. 4. Accuracy of energy (left) and toroidal angular momentum (right) conservation with respect to the time step size. Dt equals 30 ns.
90
28
T. Korpilo et al. / Journal of Computational Physics 239 (2013) 22–29
considered in the scan in order to have a time step clearly shorter and clearly longer than the time step limit set by the parallel accuracy condition. The scan shows that shortening the time step improves the conservation of toroidal angular momentum regardless of whether the Velocity Verlet or Euler scheme is applied in the simulation. The observed convergence rate is, however, slower than 1=Dt for both schemes, thus making the particle number convergence clearly more effective. Energy conservation, on the other hand, improves steadily with respect to the time step size with the convergence rate nearly as good as with respect to the particle number when the Velocity Verlet scheme is considered. In the case of the Euler scheme, the picture is quite the opposite, concluding that energy conservation cannot be improved with this scheme by reducing the time step size within the given limits. The simulation case employing the Velocity Verlet algorithm with the three-point method for the electric field is omitted here from the scan, since the scheme becomes numerically unstable in the simulation with half the time step size. The instability appears as strong electron numerical heating in the low-temperature plasma edge region and is limited to the domain satisfying v te Dt=Dz K 0:15. The core plasma, instead, is stable explaining why the instability is not seen in the simulation with the base step size. The reason for the instability can be understood by comparing the schemes for evaluating the parallel velocity v k;nþ1 in the Velocity Verlet algorithm and in the backward Euler method. As the difference between the schemes is of the second order, effectively a_ k Dt2 =2, the Velocity Verlet algorithm converges towards the Euler method when shortening the time step, in consequence of which the numerical properties of the Euler method become more and more dominant. As shown already in Fig. 2, one of these numerical properties is violent electron numerical heating in the simulations employing the three-point method for the electric field. Although the simulation with half the time step size is found unstable, the simulations with twice the time step size as well as with the base step size are fully stable and, in the 2Dt simulation, the observed errors in the conservation of total energy and toroidal angular momentum are Kðt ¼ 90 lsÞ=Kðt ¼ 0Þ ¼ 1:041 and Lðt ¼ 90 lsÞ=Lðt ¼ 0Þ ¼ 1:090, respectively.
5. Summary A numerically stable method for simulating the full 5D distribution function of drift-kinetic electrons is presented for electrostatic gyrokinetic particle-in-cell simulations of toroidal plasmas. The feasibility of the method for gyrokinetic electrons and electromagnetism is also expected but, in the absence of verification, the discussion of this matter will be addressed elsewhere. In the method, to allow a time step size close to the accuracy limit set by the Courant condition of the grid for the electron parallel motion, the electron parallel nonlinearity is solved by the direct implicit method. This means that the electron parallel acceleration (bE part of it) and the concomitant parallel displacement are advanced implicitly in terms of electric potential by requiring quasi-neutrality after the push. The accuracy of the parallel push is determined by the integration scheme used for the push, in addition to the time step size. The requirement of the linearity of the scheme (with respect to electric potential) and the fact that the push is made prior to the rest of the electron motion in the sense of numerical accuracy of the hybrid implicit-explicit particle mover restrict the number of integration schemes suitable for the problem. The two stable integration schemes, known so far by the authors, are the Velocity Verlet algorithm presented in Eq. (2) and the Euler scheme, in which the parallel velocity due to the implicit parallel acceleration is integrated in time by the backward Euler method, contrary to the former scheme. The Velocity Verlet algorithm and the Euler scheme are compared to each other using both a continuous and a piecewise constant electric field interpolation in the implicit parallel acceleration. The electric field is evaluated by the socalled energy conserving interpolation scheme, and depending on the shape function in the parallel direction, the resulting electric field is continuous or piecewise constant. The simulation results show that the Velocity Verlet algorithm with the continuous electric field interpolation conserves the total energy best and, as the scheme shows the best convergence rate with respect to the particle number, the advantage gets bigger when the number of simulation particles is increased. The time step size scan, however, reveals that the scheme becomes numerically unstable in the simulation with the time step size clearly shorter than that given by the parallel accuracy condition. Shortening the time step size affects the stability of the scheme because the Velocity Verlet algorithm converges towards the Euler scheme with decreasing step size, and the Euler scheme with the continuous electric field interpolation is found severely unstable regardless of the simulation parameters. In the simulation case considered in this paper, the Velocity Verlet algorithm with the continuous electric field interpolation becomes unstable in the plasma region of v te Dt=Dz K 0:15 and, together with the parallel accuracy condition v te Dt=Dz K 0:3, a narrow operation window is suggested for the scheme. In situations where the scheme cannot be employed in a stable way with a reasonable time step size, it is possible to use the piecewise constant electric field interpolation. The Velocity Verlet algorithm with the piecewise constant electric field interpolation shows to improve energy conservation steadily with respect to the time step size, and the convergence rate with respect to the particle number is nearly as good as with the continuous electric field interpolation. The Euler scheme with the piecewise constant electric field interpolation, instead, fails to conserve total energy within a reasonable level, in consequence of poor convergence rates with respect to both the time step size and the particle number. Although the accuracy of the energy conservation depends greatly on the numerical schemes applied to the implicit parallel acceleration, the accuracy of the conservation of total toroidal angular momentum is within a factor of two between the schemes and, in addition, the convergence rates with respect to the time step size and particle number are similar between the schemes.
T. Korpilo et al. / Journal of Computational Physics 239 (2013) 22–29
29
Acknowledgements This work was supported by the Academy of Finland (Grant Nos. 122435 and 134977), EFDA Topical Group Activities, High Level Support Team, and the supercomputing facilities of CSC, HPC-FF and PRACE. References [1] W.W. Lee, Gyrokinetic particle simulation model, J. Comput. Phys. 72 (1987) 243–269. [2] J.A. Heikkinen, S.J. Janhunen, T.P. Kiviniemi, F. Ogando, Full f gyrokinetic method for particle simulation of tokamak transport, J. Comput. Phys. 227 (2008) 5582–5609. [3] B.I. Cohen, A.M. Dimits, Implicit, partially linearized, electromagnetic particle simulation of plasma drift-wave turbulence, Phys. Rev. E 56 (1997) 2151– 2160. [4] Z. Lin, Y. Nishimura, Y. Xiao, I. Holod, W.L. Zhang, L. Chen, Global gyrokinetic particle simulations with kinetic electrons, Plasma Phys. Control. Fus. 49 (2007) B163–B172. [5] I. Manuilskiy, W.W. Lee, The split-weight particle simulation scheme for plasmas, Phys. Plasmas 7 (2000) 1381–1385. [6] Y. Chen, S. Parker, Gyrokinetic turbulence simulations with kinetic electrons, Phys. Plasmas 8 (2001) 2095–2100. [7] B.I. Cohen, A.M. Dimits, W.M. Nevins, Y. Chen, S. Parker, Kinetic electron closures for electromagnetic simulation of drift and shear-Alfvén waves. I., Phys. Plasmas 9 (2002) 251–262. [8] B.I. Cohen, A.M. Dimits, W.M. Nevins, Y. Chen, S. Parker, Kinetic electron closures for electromagnetic simulation of drift and shear-Alfvén waves. II, Phys. Plasmas 9 (2002) 1915–1924. [9] A. Friedman, A.B. Langdon, B.I. Cohen, A direct method for implicit particle-in-cell simulation, Comments Plasma Phys. Control. Fus. 6 (1981) 225–236. [10] B.I. Cohen, A.B. Langdon, A. Friedman, Implicit time integration for plasma simulation, J. Comput. Phys. 46 (1982) 15–38. [11] A.B. Langdon, B.I. Cohen, A. Friedman, Direct implicit large time-step particle simulation of plasmas, J. Comput. Phys. 51 (1983) 107–138. [12] J.A. Heikkinen, M. Nora, Gyrokinetic equations and full f solution method based on Dirac’s constrained Hamiltonian and inverse Kruskal iteration, Phys. Plasmas 18 (2011) 022310–022313. [13] T. Takizuka, H. Abe, A binary collision model for plasma simulation with a particle code, J. Comput. Phys. 25 (1977) 205–219. [14] J.A. Heikkinen, T. Korpilo, S.J. Janhunen, T.P. Kiviniemi, S. Leerink, F. Ogando, Interpolation for momentum conservation in 3D toroidal gyrokinetic particle simulation of plasmas, Comp. Phys. Commun. 183 (2012) 1719–1727. [15] C.K. Birdsall, A.B. Langdon, Plasma Physics via Computer Simulation, IOP Publishing, Bristol, 1995. [16] A. Friedman, A second-order implicit particle mover with adjustable damping, J. Comput. Phys. 90 (1990) 292–312.
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Numerically stable method for kinetic electrons in gyrokinetic particle-in-cell simulation of toroidal plasmas T. Korpilo a,⇑, J.A. Heikkinen b, S.J. Janhunen a, T.P. Kiviniemi a, S. Leerink a, F. Ogando c a
Euratom-Tekes Association, Aalto University, P.O. Box 14100, FI-00076 Aalto, Finland Euratom-Tekes Association, VTT, P.O. Box 1000, FI-02044 VTT, Finland c Universidad Nacional de Educacion a Distancia, C/Juan del Rosal, 12, 28040 Madrid, Spain b
a r t i c l e
i n f o
Article history: Received 21 March 2012 Received in revised form 5 October 2012 Accepted 22 December 2012 Available online 10 January 2013 Keywords: Kinetic electrons Gyrokinetic particle simulation Fusion plasma physics
a b s t r a c t The direct implicit method with a second-order implicit integration scheme is formulated for and applied to the electron parallel nonlinearity in global electrostatic gyrokinetic particle-in-cell simulations of toroidal fusion plasmas. The method shows improved numerical accuracy and stability properties compared to the direct implicit method with a first-order integration scheme. The conservation of total energy and toroidal angular momentum are analyzed by both techniques and the results are presented. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction The computational power provided by modern supercomputers is making predictive particle simulations of global plasma transport possible for magnetized fusion devices like tokamaks. Although the full picture of particle and heat transport is still beyond the reach of such simulations, many of the mechanisms responsible for electrostatic drift-wave turbulence are included by treating ions gyrokinetically and electrons drift-kinetically. Compared to full Lorentz force ions and electrons, these approximations reduce numerical noise and relax significantly the strict limitations for temporal and spatial resolution while retaining finite ion Larmor radius effects relevant for plasma turbulence. Kinetic effects, the electron parallel flow in particular, are important also for the formation and evolution of the parallel current and for electromagnetic fluctuations. Thus, in first-principle simulations of low-frequency turbulence, the full ion distribution is often solved gyrokinetically while electrons are modelled drift-kinetically. In addition to the full ion distribution, the proper description of turbulent transport necessitates also resolving the full electron distribution. The simplest solution is to push drift-kinetic electron motion explicitly. In electrostatic simulations, the time step is then constrained by both the accuracy condition set by the electron transit time along the magnetic field and the stability condition referred in the literature to being imposed by the unstable xh mode [1]. A way of avoiding numerical instabilities related to the explicit integration of electrons is the direct implicit method for the electron parallel nonlinearity [2]. The method makes no assumptions about the electron distribution function and is based on treating the electron acceleration by the parallel electric field (the electron parallel nonlinearity) implicitly while the remainder of the driftkinetic motion is pushed explicitly. The time step is limited only by the accuracy conditions. The direct implicit method is employed in an approximate fashion in Ref. [3] where the implicit response to the electron parallel nonlinearity is considered in the form of an analytic differential operator, thus including assumptions about the distribution function. There are ⇑ Corresponding author. Tel.: +358 50 4327580. E-mail address: tuomas.korpilo@aalto.fi (T. Korpilo). 0021-9991/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.12.033
T. Korpilo et al. / Journal of Computational Physics 239 (2013) 22–29
23
also other approximate models such as the fluid-kinetic hybrid scheme [4], the split-weight scheme [5,6] and the kinetic closure scheme [7,8] that are introduced for adding electron kinetic effects and electromagnetism into turbulence simulations. The main motivation of such schemes is to ease the strict time step constraint imposed by explicit electrons and thus reduce computational burden. In the electrostatic limit, electrons can also be modelled adiabatically, i.e. as a Maxwell–Boltzmann response, without resolving the distribution function. However, this heavy approximation eliminates net turbulent particle transport (due to the phasing between electron density and electric potential) and electron thermal transport but circumvents any numerical restriction set by the electrons. In the present work, the direct implicit method for the electron parallel nonlinearity is presented in detail in the electrostatic limit for tokamak plasma configurations. The value of the scheme lies in its inherent numerical stability for resolving the full electron distribution kinetically but it also coheres with the implicit construction of the gyrokinetic Poisson equation from the ion polarization drift motion that is used in the ELMFIRE code. The direct implicit method is originally based on the works [9–11] where the method is introduced and employed in particle-in-cell (PIC) simulations. Here, the implementation of the method for the electron parallel nonlinearity is shown for global gyrokinetic PIC drift-wave turbulence simulations. The present work will also extend the implicit method already demonstrated in Ref. [2] by considering a second-order integration scheme for the electron parallel nonlinearity. This more accurate integration scheme leads to improved energy conservation, through reduced numerical heating, and enables stable simulations when a continuous finite difference method is used for the electric field appearing in the parallel acceleration. The paper is organized as follows. In Section 2, the ELMFIRE code is briefly described to the extent required for presenting the implementation of the implicit method for electrons. In Section 3, the actual implementation as well as the new and the old integrator for the parallel nonlinearity are shown. In Section 4, the integrators are compared to each other with respect to the accuracy of the conservation of total energy and toroidal angular momentum observed in a turbulent simulation. A short summary is drawn in Section 5. 2. Code description The implicit method for simulating kinetic electrons is described here for ELMFIRE [2], for a global PIC code that advances full 5D distributions of gyrokinetic ions and drift-kinetic electrons under 3D electrostatic potential in a quasi-toroidal (axisymmetric) magnetic equilibrium of tokamak plasmas. In the code, in accordance with the gyrokinetic formalism in Ref. [12], the polarization drift is included in the gyrokinetic equations of motion, leaving only short-wavelength trembling to affect the gyroradius, and the gyrokinetic Poisson equation is reduced to the quasi-neutrality condition. In order to obtain the quasi-neutrality, the ion polarization drift and the electron acceleration by the parallel electric field are integrated implicitly in time, separately from the rest of the terms in the equations of motion that are integrated explicitly. The charge density responses to the ion polarization drift and the electron parallel acceleration are constructed by the direct implicit method and are separated from the quasi-neutral charge densities. The remaining charge densities comprise the charge separation caused by the explicit particle motion. The electrostatic potential is solved from the resulting linear system, after which the actual implicit motion is performed, neutralizing the plasma. As electrons are drift-kinetic in ELMFIRE, their motion is described by the standard guiding center motion. The equations for advancing the guiding center coordinates of a particle in time are given for three position coordinates, R, and parallel velocity, v k , by
! ml v k Bb þ rbþb rB þ r/ q q mv B q v_ k ¼ b þ k r b lrB þ r/ B m qB 1 R_ ¼ B
mv 2k
ð1Þ
together with l_ ¼ 0. Here, the magnetic moment, l, is defined as a function of the perpendicular velocity as v 2? ¼ 2Bl, and B ¼ B þ ðmv k =qÞb r b. In addition, / denotes the electric potential, B the strength of the magnetic field and b the unit vector of the magnetic field, and q and m are the charge and mass of the particle, respectively. The term in the electron motion used for obtaining the quasi-neutrality is now the parallel acceleration ak ¼ ðqB=mB Þb r/ in the equation of v_ k . This part of the electron motion is therefore integrated implicitly, while the rest of the motion is followed explicitly. For clarity and later purposes, we separate the explicit and implicit accelerations in Eq. (1) and rewrite the equation for the electron ~_ k þ ak . parallel acceleration as v_ k ¼ v The quasi-neutrality of the plasma and the concomitant electric potential are solved in a structured uniform grid set up on the quasi-ballooning coordinates, which are again built on the Boozer coordinates of the quasi-toroidal magnetic equilibrium. This equilibrium has concentric circular flux surfaces and thus the poloidal flux, Wp , can be substituted by the radial coordinate, r, in these coordinate systems. The charge density is sampled from electron guiding center positions and from four gyro-ring points of ions to the quasi-ballooning grid by trilinear interpolation by neglecting the non-orthogonality but retaining the radial shear of the coordinate system. Consequently, the charge sampling equals the linear cloud-in-cell (CIC) technique; the center of the cloud extends over radial shells by constant poloidal Boozer angle. The gyrocenter motion of ions and the guiding center motion of electrons, as expressed in the Boozer coordinates, are followed by the classical fourth-order Runge–Kutta scheme, except for the electron and ion acceleration by the parallel electric field and the ion polarization drift. The electron parallel acceleration ak is integrated in time by a second-order implicit scheme (to be discussed in
24
T. Korpilo et al. / Journal of Computational Physics 239 (2013) 22–29
the next Section) and, for symmetry reasons, the ion parallel acceleration is integrated in the same way. The integrator for the ion polarization drift is the first-order backward Euler method. Collisions among ions and electrons are modelled in the drift-kinetic limit, thus eliminating classical diffusion across the magnetic field, and are evaluated with the momentum and energy conserving binary collision model [13]. The primitive velocity vectors, in terms of which the collision model is cast, are obtained from the guiding center (gyrocenter) velocity components v k and v ? by decomposing v ? into its two orthogonal components with a random gyroangle. The covariant components of the electric field present in the equations of motion are interpolated from the quasi-ballooning grid to electron guiding center positions and via gyroaveraging to ion gyrocenter positions by two methods. First, the electric field for the E B drift motion is evaluated by a momentum conserving scheme [14]. Second, the electric field for the parallel acceleration (ion and electron) and ion polarization drift is calculated by an energy conserving scheme which P is expressed generally as EðxÞ ¼ r/ðxÞ ¼ /k rSðx xk Þ in the electrostatic case. Here, the subscript k indexes the grid points surrounding the particle, and S is a 3D shape function determining the interpolation scheme between the particle at x and the grid points at xk . As in the charge sampling, the shape function applied here is the trilinear interpolation, except that a second-order quadratic spline [15] is used in the direction of the derivative. The derivative of this spline equals the derivative of the second-order Lagrange interpolation. However, it is also possible to use the fully trilinear shape function in the electric field for the parallel acceleration. The spline utilizes three potential values from the grid (the three-point method) and produces a continuous electric field, whereas the linear interpolation uses two potential values from the grid (the two-point method) and gives a piecewise constant electric field.
3. Direct implicit method for kinetic electrons Electron parallel dynamics determines normally the largest time step for a gyrokinetic simulation. In order to extend the time step limitation near the electron parallel transit time in a stable way, the electron acceleration by the parallel electric field needs to be treated implicitly. In that case, the time step has to satisfy the accuracy condition given as kk v te Dt K 1, or Dt K Dz=3v te . Here, kk denotes the largest wavenumber of the parallel electric field, v te the electron thermal velocity, Dt the time step and Dz the metric grid spacing in the parallel direction. Another restriction for the time step is imposed by the particle motion perpendicular to the magnetic field, and in particular by the E B drift motion. The concomitant accuracy condition k? v EB Dt K 1 is, however, usually less stringent than the accuracy condition set by the electron parallel motion. An important factor in the direct implicit method is the choice for the integrator by which the implicit parallel acceleration, ak , is to be pushed. The only restrictions for the integrator are, first, that the parallel displacement resulting from the acceleration is linear with respect to the electric potential. This condition guarantees the solvability of the matrix equation arising in the implicit method. Second, to make physical sense, the linear integrator should advance the parallel coordinate and velocity in a consistent way. It is also worth noticing that the implicit parallel motion, together with the implicit ion polarization drift motion, gives a new potential that is used in the explicit particle push in the next time step. Therefore, the hybrid implicit-explicit integration scheme divides the time integration of the equations of motion on two consecutive time steps making the implicit particle motion an a priori estimation for the next time step in the sense of numerical accuracy of the integration. In this work, we consider the Velocity Verlet algorithm, a variant of the Leapfrog method, in an implicit form for the time integration of the electron parallel acceleration ak . The algorithm is the same as the undamped limit of the adjustable implicit particle mover proposed for implicit PIC simulations in Ref. [16]. This algorithm accurate to the second order in the time step fulfils the requirements for the integrator and is described in the guiding center coordinates as
1 Rnþ1 ¼ Rn þ R_ ðzn ; t n ÞDt þ ak ðznþ1 ; t nþ1 ÞDt 2 bðRnþ1 Þ 2 a ðz ; t Þ þ ak ðzn ; t n Þ v k;nþ1 ¼ v k;n þ v~_ k ðzn ; tn ÞDt þ k nþ1 nþ1 Dt 2
ð2Þ
~_ k , are included for completeness. The time where the guiding center equations of motion integrated explicitly in time, R_ and v level of the integration is given by the subscript n, and z represents the coordinates ðR; v k Þ. In the algorithm, the implicit acceleration ak is considered not only to change the particle parallel velocity but also to cause a small parallel shift, dR ¼ ðak Dt 2 =2Þb, to the particle position through the term ðB=B Þv k b present in the guiding center equations of motion (the contribution of B=B is omitted). The direct influence of the implicit acceleration is neglected on the other v k -dependent terms in the equations of motion. In ELMFIRE, electrons are not pushed exactly as described above but, as already mentioned, the explicit motion is integrated by the classical Runge–Kutta scheme unlike suggested by the algorithm. This is because the E B drift motion is not integrated accurately in time by a low-order scheme. In addition, the acceleration ak ðznþ1 ; tnþ1 Þ and ~ nþ1 ¼ Rn þ Rðz _ n ; t n ÞDt, instead of the magnetic field bðRnþ1 Þ are evaluated at the end position of the explicit motion, i.e. at R Rnþ1 since dR is short compared to any resulting parameter change and the size of a grid cell. Similarly, the parallel velocity v~ k;nþ1 ¼ v k;n þ v~_ k ðzn ; tn ÞDt is employed in the integration of the implicit motion, as ak has only a weak dependence on v k through B . Consequently, a costly iteration procedure is avoided in the algorithm. We will refer to the algorithm in Eq. (2) as the Velocity Verlet algorithm.
T. Korpilo et al. / Journal of Computational Physics 239 (2013) 22–29
25
After the selection of the integrator for the implicit particle movement, the direct implicit method is straightforward. The starting point for the method is the requirement for the plasma quasi-neutrality and, following the ideas presented in Ref. [9], the quasi-neutrality is enforced by the implicit particle movement. This means, with the notation of Eq. (2), that the implicit electron parallel acceleration is integrated at the end of a time step, in order to solve the unknown electric potential /ðtnþ1 Þ. Let us denote the quasi-neutral electron charge density discretized on the quasi-ballooning grid at the time level tnþ1 P P by qk ðxnþ1 Þ. The charge density is evaluated on the grid as qk ¼ qk;p ¼ qp Sðxp xk Þ=V k , where qp denotes the total charge of the pth electron cloud, V k the volume of the kth grid cell and x the charge sampling coordinates, i.e. the radial coordinate r, the Boozer poloidal angle h and the field-aligned (quasi-ballooning) angle f. Now, Taylor expanding the shape function S up to the first order with respect to the charge sampling coordinates of the pth electron at the end point of the explicit push, ~nþ1 , one obtains x
qk ðxnþ1 Þ ¼ qk ðx~nþ1 Þ þ
X @ qk;p ðx ~nþ1 Þ p
@fp
dfp ð~znþ1 ; t nþ1 Þ:
ð3Þ
The second term denotes here the linear increment to the former charge density resulting from the implicit electron parallel acceleration and thus dfp ¼ dRp rf (the small discrepancy between the orientation of b and the quasi-ballooning parallel vector ef is neglected in this paper for convenience). With the trilinear shape function, this linear increment to the initial electron charge density is exact for all the electrons not crossing any grid points in the parallel direction during the implicit shift. The increment is now depending linearly on the electric potential /k ðt nþ1 Þ through the shift dfp , and together with the contribution from the similarly constructed implicit ion polarization charge density, the quasi-neutrality condition reduces to a linear matrix equation for the electric potential, from which the potential is solved by conventional means. With the solution for /k ðt nþ1 Þ, the parallel shift of each electron, dfp , is solved, and the electrons are moved accordingly to ~ nþ1 þ dR. At the same time, the electron parallel velocity is updated to the value of their quasi-neutral positions Rnþ1 ¼ R v k;nþ1 . At the next time step, the time integration of the electron explicit motion is started from the point Rnþ1 ; v k;nþ1 (in the absence of binary collisions) by using the solved electric potential /k ðt nþ1 Þ. The explicit push reproduces a charge separation between ions and electrons, and the implicit step is repeated. The Velocity Verlet algorithm considers the electron parallel acceleration ak up to the second order in the time step for both the position and the velocity coordinate [16]. In Ref. [2], the implicit integrator for the kinetic electrons differs from the presented algorithm in that the change of the parallel velocity is integrated from the parallel acceleration by the first-order backward Euler method, i.e. v k;nþ1 ¼ v k;n þ ak ðznþ1 ; tnþ1 ÞDt. The concomitant parallel shift is obtained from the Taylor series expansion Rk;n ¼ Rk;nþ1 v k;nþ1 Dt þ ak ðznþ1 ; tnþ1 ÞDt2 =2 by substituting the expression for v k;nþ1 in it (the term B=B is again ignored). The resulting algorithm, with the inclusion of the explicitly integrated guiding center motion, is the same as the Velocity Verlet algorithm in Eq. (2), except that the velocity change due to the implicit parallel acceleration is evaluated as dv k ¼ ak ðznþ1 ; tnþ1 ÞDt instead of dv k ¼ ½ak ðznþ1 ; tnþ1 Þ þ ak ðzn ; tn ÞDt=2. In the next Section, the difference between the first-order and the second-order accurate dv k is considered relative to a numerical simulation. 4. Simulation results In the standard version of ELMFIRE, the change of the parallel velocity due to the implicit electron acceleration has been integrated in time by the backward Euler method. Moreover, the covariant component of the electric field in the parallel acceleration ak has been evaluated by the two-point method. To compare this Euler scheme to the Velocity Verlet algorithm by using both the two-point (discontinuous) and the three-point (continuous) method for the electric field, a series of simulations are performed with identical simulation parameters. In the simulations, the main parameters of deuterium plasma are similar to FT-2 tokamak discharges: minor radius of 8 cm, major radius of 55 cm, magnetic field of 2.2 T (on the axis) and total plasma current of 55 kA. The radial current density profile is linear. Further, the initial radial density profile is parabolic with the maximum value of 5:0 1019 m3 at the magnetic axis and the minimum value of 0:1 1019 m3 at the separatrix. The initial radial electron (and ion) temperature profile is also parabolic and has a maximum value of 120 eV at the axis and a minimum value of 10 eV at the separatrix. Together with the time step of 30 ns, the parallel accuracy condition is satisfied in the simulations, i.e. the parameter v te Dt=Dz lies between 0.45 and 0.15. The simulation domain is limited between the normalized minor radius of 0.25 and 1.0 and is divided into 50 radial shells each of which consists of 200 grid cells in poloidal direction and 8 grid cells in toroidal direction. The average number of electrons (and ions) per grid cell is 1100. The comparison between the integration schemes employed in the implicit electron acceleration is shown in Figs. 1–4, including the particle number and the time step scan. The comparison is based on the accuracy of the conservation of total energy and toroidal angular momentum, as electron parallel dynamics is vulnerable to the numerical heating/cooling and as the integration of the parallel acceleration ak can cause a numerical source term for toroidal angular momentum. In the simulations, total energy and toroidal angular momentum are evaluated in a diagnostic domain that is composed of the simulation domain excluding the five first and five last radial shells. This exclusion eliminates boundary effects from contributing to the conserving variables, and the only physical sinks and sources in the diagnostic domain to be considered are the energy flux and the toroidal angular momentum transport, carried by the particles, through the domain boundaries. The expressions for the conserved total energy K and toroidal angular momentum L, consistent with the gyrokinetic formalism used in ELMFIRE, are given in Refs. [12,14]. In the figures, the initial value of total energy, Kðt ¼ 0Þ, and toroidal angular momentum,
26
T. Korpilo et al. / Journal of Computational Physics 239 (2013) 22–29
Lðt ¼ 0Þ, in the diagnostic domain are 28 J and 2:4 106 m2 kg/s, respectively. The latter value comprises the ion banana current, arisen as ions are initialized on their neoclassical (collisionless) orbits. In the following, the Velocity Verlet algorithm is shown to provide a numerically stable method to obtain acceptable toroidal angular momentum and energy conservation. As shown in Fig. 1, both the Velocity Verlet algorithm and the three-point method for the electric field improve the conservation of total energy, compared to the Euler scheme with the two-point method for the electric field, while an opposite trend is seen in the conservation of total toroidal angular momentum. The improved energy conservation is explained by reduced numerical heating/cooling of electrons. The increased error in the conservation of toroidal angular momentum, on the other hand, is explained by increased asymmetry between the integration of the parallel acceleration ak of ions and electrons, enabled by the inhomogeneous charge density in the space between the grid points. The inhomogeneity exists even if the charge density sampled at the grid points is neutral. These results indicate that the numerical heating/cooling of electrons plays a dominant role in the time evolution of total energy and that the asymmetric integration of the electron and ion parallel acceleration affects notably the time evolution of total toroidal angular momentum. Fig. 2 illustrates a severe numerical instability arising when the three-point method for the electric field is used together with the Euler scheme. In other words, electrons experience violent numerical heating which sets in after 10th time step and as a result of which electron energy content doubles in a very short time period of 1 ls. The instability has been observed in all simulations, regardless of the simulation parameters. A similar type of instability results both with the Euler scheme and with the Velocity Verlet algorithm if a momentum conserving interpolation scheme is utilized in evaluating the electric field P for the parallel acceleration. In other words, if the electric field is calculated from the expression EðxÞ ¼ Ek Sðx xk Þ by employing the trilinear shape function and the central difference method for the covariant components of the electric field at the grid points.
1.06
1.08 Euler, 2−point
Verlet, 3−point
Verlet, 2−point
Verlet, 2−point
Verlet, 3−point
Euler, 2−point
1.06
L(t) / L(t=0)
K(t) / K(t=0)
1.04
1.04
1.02 1.02
1
500
1000
1500
2000
2500
1
3000
t / 30 ns
500
1000
1500
2000
2500
3000
t / 30 ns
Fig. 1. Accuracy of energy (left) and toroidal angular momentum (right) conservation resulting from integrating the electron response to the parallel electric field by Velocity Verlet and Euler schemes. The dashed line is reflected around x-axis for illustrative purposes in the picture on the left.
1.5 Euler, 3−point Euler, 2−point
1.4
K(t) / K(t=0)
1.3 1.2 1.1 1 0.9
10
20
30
40
50
t / 30 ns Fig. 2. The Euler scheme with the 3-point method for the electric field is unstable.
27
T. Korpilo et al. / Journal of Computational Physics 239 (2013) 22–29
The accuracy of the conservation of total energy and toroidal angular momentum is presented in Fig. 3 in terms of the number of simulated particles. For this, the simulation with the Velocity Verlet algorithm and with the three-point method for the electric field is repeated with twice (2 N), four times (4 N) and eight times (8 N) the particle number per grid cell. The particle numbers of 1 N, 2 N, 4 N and 8 N translate into the noise fluctuation levels dn=n of 0.84%, 0.57%, 0.40% and 0.28%, respectively, while the observed turbulence fluctuation level is 1–5% depending on the radial position. The simulation results reveal that the energy conservation improves according to 1=N and that the conservation of toroidal angular momentum improves even faster. However, the fluctuating time evolution of toroidal angular momentum in the 8 N case indicates that the error in the conservation of toroidal angular momentum is not vanishing only by increasing the number of particles. A potential reason for this, namely the numerical accuracy of the diagnostics, is discussed in Ref. [14]. Nevertheless, 0.25% error in the energy conservation and under 0.4% error in the toroidal angular momentum conservation is observed in the 8 N case during 90 ls or 3000 time steps. The Euler scheme and the Velocity Verlet algorithm with the two-point method for the electric field, in turn, improve the conservation of toroidal angular momentum in terms of the particle number as the Velocity Verlet algorithm with the three-point method, while in energy conservation the picture is different. In other words, the Velocity Verlet algorithm with the two-point method improves energy conservation slower than 1=N, and the Euler scheme shows practically no improvement. In the simulations, the growth of turbulent modes starts around the 500th time step after which turbulence is a strong candidate in explaining the changes in the time evolutions of total energy and toroidal angular momentum. The effect of the time step size on the accuracy of the conservation of total energy and toroidal angular momentum is presented in Fig. 4 for two simulation cases, one employing the Euler scheme and the other the Velocity Verlet algorithm with the two-point method for the electric field. In addition to the base step size, one-half and twice the time step size is
1.08 Verlet, 3−point, 1N
1.04
Verlet, 3−point, 1N
Verlet, 3−point, 2N
Verlet, 3−point, 2N
Verlet, 3−point, 4N Verlet, 3−point, 8N Verlet, 2−point, 4N
L(t) / L(t=0)
K(t) / K(t=0)
1.03
Verlet, 3−point, 4N
1.06
Euler, 2−point, 4N 1.02
Verlet, 3−point, 8N Verlet, 2−point, 4N Euler, 2−point, 4N
1.04
1.02 1.01 1 1
500
1000
1500
2000
2500
3000
500
1000
t / 30 ns
1500
2000
2500
3000
t / 30 ns
Fig. 3. Accuracy of energy (left) and toroidal angular momentum (right) conservation with respect to particle number. The crossed line is reflected around x-axis for illustrative purposes in the picture on the left. N equals 1100 particles per grid cell.
1.06
Verlet, 2−point, 2Δt
Verlet, 2−point, Δt
1.06
Verlet, 2−point, Δt/2
1.02 1 0.98 0.96
Verlet, 2−point, Δt Verlet, 2−point, Δt/2
1.05
L(t) / L(t=0)
K(t) / K(t=0)
1.04
Verlet, 2−point, 2Δt
Euler, 2−point, 2Δt Euler, 2−point, Δt
1.04
Euler, 2−point, Δt/2
1.03 1.02
Euler, 2−point, 2Δt 1.01
Euler, 2−point, Δt Euler, 2−point, Δt/2
0.94
15
30
45
t (μs)
60
75
90
1
15
30
45
60
75
t (μs)
Fig. 4. Accuracy of energy (left) and toroidal angular momentum (right) conservation with respect to the time step size. Dt equals 30 ns.
90
28
T. Korpilo et al. / Journal of Computational Physics 239 (2013) 22–29
considered in the scan in order to have a time step clearly shorter and clearly longer than the time step limit set by the parallel accuracy condition. The scan shows that shortening the time step improves the conservation of toroidal angular momentum regardless of whether the Velocity Verlet or Euler scheme is applied in the simulation. The observed convergence rate is, however, slower than 1=Dt for both schemes, thus making the particle number convergence clearly more effective. Energy conservation, on the other hand, improves steadily with respect to the time step size with the convergence rate nearly as good as with respect to the particle number when the Velocity Verlet scheme is considered. In the case of the Euler scheme, the picture is quite the opposite, concluding that energy conservation cannot be improved with this scheme by reducing the time step size within the given limits. The simulation case employing the Velocity Verlet algorithm with the three-point method for the electric field is omitted here from the scan, since the scheme becomes numerically unstable in the simulation with half the time step size. The instability appears as strong electron numerical heating in the low-temperature plasma edge region and is limited to the domain satisfying v te Dt=Dz K 0:15. The core plasma, instead, is stable explaining why the instability is not seen in the simulation with the base step size. The reason for the instability can be understood by comparing the schemes for evaluating the parallel velocity v k;nþ1 in the Velocity Verlet algorithm and in the backward Euler method. As the difference between the schemes is of the second order, effectively a_ k Dt2 =2, the Velocity Verlet algorithm converges towards the Euler method when shortening the time step, in consequence of which the numerical properties of the Euler method become more and more dominant. As shown already in Fig. 2, one of these numerical properties is violent electron numerical heating in the simulations employing the three-point method for the electric field. Although the simulation with half the time step size is found unstable, the simulations with twice the time step size as well as with the base step size are fully stable and, in the 2Dt simulation, the observed errors in the conservation of total energy and toroidal angular momentum are Kðt ¼ 90 lsÞ=Kðt ¼ 0Þ ¼ 1:041 and Lðt ¼ 90 lsÞ=Lðt ¼ 0Þ ¼ 1:090, respectively.
5. Summary A numerically stable method for simulating the full 5D distribution function of drift-kinetic electrons is presented for electrostatic gyrokinetic particle-in-cell simulations of toroidal plasmas. The feasibility of the method for gyrokinetic electrons and electromagnetism is also expected but, in the absence of verification, the discussion of this matter will be addressed elsewhere. In the method, to allow a time step size close to the accuracy limit set by the Courant condition of the grid for the electron parallel motion, the electron parallel nonlinearity is solved by the direct implicit method. This means that the electron parallel acceleration (bE part of it) and the concomitant parallel displacement are advanced implicitly in terms of electric potential by requiring quasi-neutrality after the push. The accuracy of the parallel push is determined by the integration scheme used for the push, in addition to the time step size. The requirement of the linearity of the scheme (with respect to electric potential) and the fact that the push is made prior to the rest of the electron motion in the sense of numerical accuracy of the hybrid implicit-explicit particle mover restrict the number of integration schemes suitable for the problem. The two stable integration schemes, known so far by the authors, are the Velocity Verlet algorithm presented in Eq. (2) and the Euler scheme, in which the parallel velocity due to the implicit parallel acceleration is integrated in time by the backward Euler method, contrary to the former scheme. The Velocity Verlet algorithm and the Euler scheme are compared to each other using both a continuous and a piecewise constant electric field interpolation in the implicit parallel acceleration. The electric field is evaluated by the socalled energy conserving interpolation scheme, and depending on the shape function in the parallel direction, the resulting electric field is continuous or piecewise constant. The simulation results show that the Velocity Verlet algorithm with the continuous electric field interpolation conserves the total energy best and, as the scheme shows the best convergence rate with respect to the particle number, the advantage gets bigger when the number of simulation particles is increased. The time step size scan, however, reveals that the scheme becomes numerically unstable in the simulation with the time step size clearly shorter than that given by the parallel accuracy condition. Shortening the time step size affects the stability of the scheme because the Velocity Verlet algorithm converges towards the Euler scheme with decreasing step size, and the Euler scheme with the continuous electric field interpolation is found severely unstable regardless of the simulation parameters. In the simulation case considered in this paper, the Velocity Verlet algorithm with the continuous electric field interpolation becomes unstable in the plasma region of v te Dt=Dz K 0:15 and, together with the parallel accuracy condition v te Dt=Dz K 0:3, a narrow operation window is suggested for the scheme. In situations where the scheme cannot be employed in a stable way with a reasonable time step size, it is possible to use the piecewise constant electric field interpolation. The Velocity Verlet algorithm with the piecewise constant electric field interpolation shows to improve energy conservation steadily with respect to the time step size, and the convergence rate with respect to the particle number is nearly as good as with the continuous electric field interpolation. The Euler scheme with the piecewise constant electric field interpolation, instead, fails to conserve total energy within a reasonable level, in consequence of poor convergence rates with respect to both the time step size and the particle number. Although the accuracy of the energy conservation depends greatly on the numerical schemes applied to the implicit parallel acceleration, the accuracy of the conservation of total toroidal angular momentum is within a factor of two between the schemes and, in addition, the convergence rates with respect to the time step size and particle number are similar between the schemes.
T. Korpilo et al. / Journal of Computational Physics 239 (2013) 22–29
29
Acknowledgements This work was supported by the Academy of Finland (Grant Nos. 122435 and 134977), EFDA Topical Group Activities, High Level Support Team, and the supercomputing facilities of CSC, HPC-FF and PRACE. References [1] W.W. Lee, Gyrokinetic particle simulation model, J. Comput. Phys. 72 (1987) 243–269. [2] J.A. Heikkinen, S.J. Janhunen, T.P. Kiviniemi, F. Ogando, Full f gyrokinetic method for particle simulation of tokamak transport, J. Comput. Phys. 227 (2008) 5582–5609. [3] B.I. Cohen, A.M. Dimits, Implicit, partially linearized, electromagnetic particle simulation of plasma drift-wave turbulence, Phys. Rev. E 56 (1997) 2151– 2160. [4] Z. Lin, Y. Nishimura, Y. Xiao, I. Holod, W.L. Zhang, L. Chen, Global gyrokinetic particle simulations with kinetic electrons, Plasma Phys. Control. Fus. 49 (2007) B163–B172. [5] I. Manuilskiy, W.W. Lee, The split-weight particle simulation scheme for plasmas, Phys. Plasmas 7 (2000) 1381–1385. [6] Y. Chen, S. Parker, Gyrokinetic turbulence simulations with kinetic electrons, Phys. Plasmas 8 (2001) 2095–2100. [7] B.I. Cohen, A.M. Dimits, W.M. Nevins, Y. Chen, S. Parker, Kinetic electron closures for electromagnetic simulation of drift and shear-Alfvén waves. I., Phys. Plasmas 9 (2002) 251–262. [8] B.I. Cohen, A.M. Dimits, W.M. Nevins, Y. Chen, S. Parker, Kinetic electron closures for electromagnetic simulation of drift and shear-Alfvén waves. II, Phys. Plasmas 9 (2002) 1915–1924. [9] A. Friedman, A.B. Langdon, B.I. Cohen, A direct method for implicit particle-in-cell simulation, Comments Plasma Phys. Control. Fus. 6 (1981) 225–236. [10] B.I. Cohen, A.B. Langdon, A. Friedman, Implicit time integration for plasma simulation, J. Comput. Phys. 46 (1982) 15–38. [11] A.B. Langdon, B.I. Cohen, A. Friedman, Direct implicit large time-step particle simulation of plasmas, J. Comput. Phys. 51 (1983) 107–138. [12] J.A. Heikkinen, M. Nora, Gyrokinetic equations and full f solution method based on Dirac’s constrained Hamiltonian and inverse Kruskal iteration, Phys. Plasmas 18 (2011) 022310–022313. [13] T. Takizuka, H. Abe, A binary collision model for plasma simulation with a particle code, J. Comput. Phys. 25 (1977) 205–219. [14] J.A. Heikkinen, T. Korpilo, S.J. Janhunen, T.P. Kiviniemi, S. Leerink, F. Ogando, Interpolation for momentum conservation in 3D toroidal gyrokinetic particle simulation of plasmas, Comp. Phys. Commun. 183 (2012) 1719–1727. [15] C.K. Birdsall, A.B. Langdon, Plasma Physics via Computer Simulation, IOP Publishing, Bristol, 1995. [16] A. Friedman, A second-order implicit particle mover with adjustable damping, J. Comput. Phys. 90 (1990) 292–312.
