# AstroGK: Astrophysical Gyrokinetics Code

## Abstract

The gyrokinetic simulation code AstroGK is developed to study fundamental aspects of kinetic plasmas and for applications mainly to astrophysical problems. AstroGK is an Eulerian slab code that solves the electromagnetic gyrokinetic-Maxwell equations in five-dimensional phase space, and is derived from the existing gyrokinetics code GS2 by removing magnetic geometry effects. Algorithms used in the code are described. The code is benchmarked using linear and nonlinear problems. Serial and parallel performance scalings are also presented.

###### keywords:

Gyrokinetic simulation, Eulerian, Numerical Methods###### Pacs:

52.30.Gz, 52.65.Tt, 94.05.-a, 95.30.Qd^{1}

=20 \biboptionssort&compress

## 1 Introduction

Gyrokinetics is a limit of kinetic theory that describes the low-frequency dynamics of weakly collisional plasmas in a mean magnetic field. Developed for the study of magnetically confined fusion plasmas, it has proven to be a valuable tool in understanding the dynamics of drift-wave turbulence, a key cause of the enhanced transport measured in modern fusion experiments that leads to poor device performance. Tremendous theoretical, computational, and experimental efforts have been devoted to this problem, with steady progress over three decades.

It has been recently recognized that the gyrokinetic approach is also well-suited to the study of astrophysical plasmas, including galaxy clusters, accretion disks around compact objects, the interstellar medium, and the solar corona and solar wind HowesCowleyDorland_06 (); HowesCowleyDorland_08 (); SchekochihinCowleyDorland_09 (). Taking advantage of the knowledge and computational techniques developed in the simulation of turbulence in fusion plasmas, AstroGK sourceforge_gyrokinetics (), a gyrokinetic simulation code, is developed specifically for the study of astrophysical problems. In this paper, we describe the algorithms employed in AstroGK and present verification tests and performance results.

Gyrokinetics describes the low-frequency fluctuations of magnetized plasmas by exploiting the timescale separation between the low-frequency dynamics of interest and the fast cyclotron motion of particles, , where is the typical frequency of fluctuations and is the cyclotron frequency. By averaging the kinetic Vlasov–Landau (or Boltzmann) equation and Maxwell’s equations over the fast cyclotron motion, a self-consistent gyrokinetic-Maxwell (GK-M) system is defined in five-dimensional phase space. This system orders out the fast MHD waves and the cyclotron resonance, but retains finite Larmor radius (FLR) effects and collisionless wave-particle interactions via the Landau resonance Landau_46 (); Barnes_66 ().

The theoretical foundation of gyrokinetics has been developed extensively over the past four decades RutherfordFrieman_68 (); TaylorHastie_68 (); Catto_78 (); AntonsenLane_80 (); CattoTangBaldwin_81 (); FriemanChen_82 (); DubinKrommesOberman_83 (); HahmLeeBrizard_88 (); Brizard_92 (); Sugama_00 (); HowesCowleyDorland_06 (); BrizardHahm_07 (), and gyrokinetics is now broadly employed for numerical studies of turbulence driven by microinstabilities in laboratory plasmas. The reduction of the phase-space dimensionality and the relaxed timestep constraints under the gyrokinetic approximation have made possible gyrokinetic simulations of fusion devices, yet it is still computationally demanding. The simulation of the kinetic dynamics of magnetized plasmas in both laboratory and astrophysical settings presents a challenge to the scientific community requiring the most advanced computing technology available.

A number of computational codes for gyrokinetics have been developed and are actively refined worldwide by the fusion community KotschenreutherRewoldtTang_95 (); LinHahmLee_98 (); DorlandJenkoKotschenreuther_00 (); JenkoDorlandKotschenreuther_00 (); CandyWaltz_03 (); ParkerChenWan_04 (); WatanabeSugama_04 (); GrandgirardSarazinAngelino_07 (); IdomuraIdaKano_09 (). These codes are typically classified by the following characteristics: Eulerian continuum vs. Lagrangian particle (PIC), local flux tube vs. global toroidal, vs. full, and electromagnetic vs. electrostatic. AstroGK has been derived from the Eulerian continuum, local flux tube, , electromagnetic code GS2 KotschenreutherRewoldtTang_95 (); DorlandJenkoKotschenreuther_00 (). GS2 was designed to simulate plasma dynamics in fusion devices where the magnetic geometry plays a central role. The handling in the code of the toroidal magnetic geometry and of particle trapping requires a rather complicated numerical implementation. In contrast, the study of the fundamental properties of kinetic plasmas for application to astrophysical situations demands the simulation of the dynamics at scales of order the particle Larmor radii on which the magnetic field is well approximated as straight or gently curved, so the coding to handle complicated magnetic geometries is unnecessary. Therefore AstroGK was created by stripping GS2 of the cumbersome coding necessary to describe the magnetic geometry effects, leading to a simplified, faster code ideally suited to study weakly collisional astrophysical plasmas. In addition to this, we believe the code is also an ideal testbed for new ideas of additional physics, diagnostics, numerical algorithms, and optimizations which may be exported back to GS2.

AstroGK has already proven its usefulness in a number of studies. For example, it has produced the first kinetic simulations of turbulence describing the transition from Alfvén to kinetic Alfvén wave turbulence at the scale of the ion Larmor radius in an attempt to understand solar wind turbulence HowesDorlandCowley_08 (), revealed nonlinear phase-mixing properties of turbulence TatsunoDorlandSchekochihin_09 (), enabled the study of the statistical properties of phase-space structures of plasma turbulence TatsunoBarnesCowley_10 (), explored the transition from collisional to collisionless tearing instabilities NumataDorlandHowes_10 (), and described the Alfvén wave dynamics in the LAPD experiment NielsonHowesTatsuno_10 ().

This paper is organized as follows: In Section 2, the set of GK-M equations solved in the code is given. Section 3 describes algorithms employed in the code, including the velocity-space discretization and integration, the finite-difference formalism of the GK-M system along the mean field direction and its solution by the special technique developed by Kotschenreuther et al. KotschenreutherRewoldtTang_95 (), the parallelization scheme, and some additional features of the code. Section 4 presents tests of AstroGK, ranging from linear electrostatic problems to nonlinear fully electromagnetic problems. Comparisons to analytic solutions of the linear problems enable the verification of the code and nonlinear examples show potential ability to apply the code to complicated problems. Serial and parallel performance measurements of the code on cutting edge supercomputers are given in Section 5. In Section 6, we present a summary of the paper.

## 2 Gyrokinetic-Maxwell equations

In this section, we present the gyrokinetic-Maxwell (GK-M) system of equations solved in AstroGK. For notational simplicity, we summarize all the symbols and their definitions in A.

We first assume that scale separations in space and time are well satisfied such that small fluctuations are locally embedded in a background plasma which is slowly varying spatially and temporally. We consider a temporally constant mean magnetic field . The mean field is pointing almost parallel to , but is allowed to have curvature . We also assume the amplitude of the mean field is constant along its direction, , but has finite gradient perpendicular to its direction, . Under this assumption, particle trapping do not take place. In the presence of a mean magnetic field, we can adopt the gyrokinetic ordering and average over the fast cyclotron motion to reduce the Vlasov–Maxwell equations to the GK-M equations; see Howes et al. HowesCowleyDorland_06 () and Schekochihin et al. SchekochihinCowleyDorland_09 () for derivations of these equations expressly intended for the study of astrophysical plasmas.

Under the gyrokinetic ordering, the distribution function of particles up to the first order is given by

(1) |

where is the zeroth-order, equilibrium
Maxwellian distribution function. The first-order part of the
distribution function is composed of the Boltzmann response term, a
term due to gradients of the equilibrium,
and the gyro-center distribution function defined in the
gyro-center coordinate .^{2}

(2) |

where parallel and perpendicular subscripts refer to directions with respect to the mean magnetic field. The perpendicular drift velocity is given by

(3) |

where the terms correspond, from left to right, to the gradient-
drift (), the curvature drift,
and a nonlinear drift.^{3}

(4) |

where . (The gyro-average at fixed particle coordinate can also be defined by switching roles of and .)

In the GK-M system, the electromagnetic fields are specified by the
three scalar functions , ,
and
^{4}

(5) |

Maxwell’s equations in the gyrokinetic limit reduce to the quasi-neutrality condition, and the parallel and perpendicular components of Ampère’s law:

(6) | |||

(7) | |||

(8) |

To lowest order, the quasi-neutrality condition and the perpendicular Ampère’s law imply the constraints on the background plasma of quasi-neutrality and total pressure balance:

(9) |

It is advantageous, both analytically and in the code, to Fourier transform the equations only in the plane perpendicular to the mean magnetic field. The distribution function and the fields are decomposed as

(10) | ||||

(11) |

with . The virtue of expressing variables in Fourier space is that the gyro-averaging operation becomes a multiplication by a Bessel function. For example, the gyro-average of the gyrokinetic potential is given by

(12) |

where is the Bessel function of the first kind with the argument , taking . The Fourier coefficients of the fields are now functions of and , for example .

In the large-scale limit , Alfvénic perturbations have a gyro-center distribution function that is largely canceled by the Boltzmann response term (see, for example, Section 5 of SchekochihinCowleyDorland_09 ()). To avoid numerical error arising from this near cancellation, a complementary distribution function is introduced, given by

(13) |

After normalization described in A.2, we finally obtain a normalized version of the GK-M equations as follows. The normalized gyrokinetic equation is

(14) |

where and the normalized gyrokinetic potential is

(15) |

The nonlinear term, given by the Poisson bracket , is evaluated in real space and then transformed into Fourier space (denoted by ). The terms in the gyrokinetic equation due to gradients of the background plasma, given in square brackets, contribute at the same order as the other terms, and are characterized by the parameters , , , and the curvature . The normalized Maxwell’s equations are

(16) | ||||

(17) | ||||

(18) |

The plasma beta of the reference species, defined by the ratio of the thermal pressure of the reference species to the magnetic pressure of the mean magnetic field, is given by . The operators to take the th order moment of the distribution function are:

(19) | ||||

(20) | ||||

(21) |

The function arises from the integration over perpendicular velocity space of products of two Bessel functions. For , it is given by

(22) |

where is the modified Bessel function of the first kind, and the argument is HowesCowleyDorland_06 ().

The background plasma must also satisfy the normalized quasi-neutrality and total pressure balance constraints:

(23) |

We defer description of the explicit form of the collision operator used in the code to B, as it has a rather cumbersome form and is fully documented in BarnesAbelDorland_09 (). We mention here the basic properties of the operator. The operator is based on the linearized Landau collision operator transformed into the gyro-center coordinate. It has second-order velocity derivatives providing diffusion in velocity space and conserving terms which include integrations over velocity space. It is constructed to satisfy Boltzmann’s H-theorem and the conservation of particles, momentum, and energy. It contains both like-species collisions and inter-species collisions, but the inter-species collisions account only for the collisions of electrons with one species of ions with large mass. Note that the linearized collision operator for a given species can be made independent of the first-order evolution of any other species. The theoretical basis of the collision operator is discussed in detail in AbelBarnesCowley_08 ().

## 3 Algorithm description

This section describes the numerical algorithms used in AstroGK to evolve the GK-M system of Eqs. (14), and (16)–(18). The gyrokinetic equation combined with the field equations together comprise a set of integro-differential equations for the evolution of the distribution function defined in the five-dimensional phase space . In this section, the species subscript is omitted unless necessary. Periodic boundary conditions are assumed for the spatial dimensions , and the derivatives in the plane perpendicular to the mean field, , are handled using a Fourier-spectral method. Except for the nonlinear term, each Fourier mode is independent of the others in the gyrokinetic equation, so we omit the subscript for simplicity. In the numerical implementation described here, because fields are calculated separately from the gyrokinetic equations, the gyrokinetic equation for each species is essentially independent of that for the other species. The mean field parallel direction and the time are discretized by () and , where is fixed and may vary to satisfy the Courant–Friedrichs–Lewy (CFL) condition CourantFriedrichsLewy_28 () for the nonlinear term. Velocity space is discretized with grid points chosen by Gaussian quadrature rules for optimal integration, generating nonuniform meshes.

### 3.1 Velocity-space integration

The velocity grid in AstroGK is specified by ,
where the pitch
angle^{5}

(24) |

Gaussian quadrature evaluates an integral of a function with weight by

(25) |

where is the th root of the th order polynomial, and is the corresponding discretized weight. The weights for the Gauss–Legendre and Gauss–Laguerre rules are given by

(26) |

where and are the th order Legendre and Laguerre polynomials, respectively AbramowitzStegun_72a (). The superscripts ‘Leg’ and ‘Lag’ to and explicitly denote that they are associated with the Gauss–Legendre and Gauss–Laguerre rules. The weight function and the integration range are chosen according to the polynomial.

For the pitch-angle integration, the Gauss–Legendre quadrature (, , ) is immediately applied by defining :

(27) |

The integration range is changed using a linear transformation to fit the range of the Gauss–Legendre rule. is the number of grid points describing the grid and is specified by the user input, ngauss: .

A rather careful treatment of the energy integral is necessary because there are singularities of the integrand at and which may prevent a simple approximation from achieving spectral convergence. To avoid the problem, we split the integration range at into a lower and an upper range and change the variable to from for the lower range integration: the Gauss–Legendre scheme is used for the lower range; and the Gauss–Laguerre (, , ) is used for the upper range. Therefore, the integration is approximated by

(28) |

where

(29) |

and the integration ranges are shifted appropriately. We allow users to specify , , , and in the code input.

We also have another mode to evaluate the energy integral in the code called the ’egrid mode’, whereas the above method is called the ’vgrid mode’. In the egrid mode, the energy integral is calculated by the method suggested by Candy and Waltz CandyWaltz_03 (), which is not exponentially accurate. We note, however, that if is very small (), we find empirically that the egrid mode may give better results than the vgrid mode. Therefore, the optimal choice of energy grid mode is governed by the simulation parameters.

Further discussion of our velocity-space coordinates can be found in BarnesDorlandTatsuno_10 ().

### 3.2 Time integration

The gyrokinetic equation is symbolically denoted by

(30) |

where is the linear term except the collision term, is the collision term which is also linear, and is the nonlinear term. We consider the time derivative using first-order finite differentiation for the linear term, treating the collision term by the implicit Euler method, and handling the nonlinear term explicitly by the third-order Adams–Bashforth method (AB3):

(31) |

where . The time-centering parameter (fexp in the code input) may be chosen within the range , where () represents a fully explicit (implicit) scheme. If , the scheme is stable for any . We mainly use an implicit trapezoidal rule which is second-order accurate and free from time step restrictions due to the linear term. Hereafter, we fix .

Note that Eq. (31) is linear with respect to due to the explicitness of the nonlinear term. We then employ a Godunov splitting technique Godunov_59 (), which is first-order accurate in , to separate the collision term:

(32) | ||||

(33) |

which greatly reduces the size of the matrix to be inverted. Solving for , we obtain:

(34) |

The method is first-order accurate in time. The use of AB3 for the nonlinear term is to make nonlinear runs stable. Note that the first time step for the nonlinear term is evaluated by the (explicit) Euler method (which is also first order in ), and the second timestep is evaluated using the second-order Adams–Bashforth scheme (AB2).

A second-order accurate method may, in principle, be derived by applying a higher-order scheme for the collision term as well, and by using a Strang splitting Strang_68 () for the operator splitting. The first two steps for the nonlinear term could also be computed using a higher-order method. These ideas are not implemented in the current version.

### 3.3 Gyrokinetic solver

We now describe the implicit advance of the linear terms in the gyrokinetic equation. Basically, the fields in the gyrokinetic equation can be obtained by a separate procedure, and the gyrokinetic equation becomes a differential equation with the given fields. Thus the collisionless gyrokinetic equation is written as

(35) |

where , coefficients and are functions of and , and contains the AB3 nonlinear term. The collision term is always consecutively applied after (35), but separately (Section 3.3.4).

We use a compact finite-difference method for evaluation of to achieve up to second-order accuracy with keeping the matrix bi-diagonal. To take into account global information, compact finite differencing schemes use the derivatives at neighboring grids to evaluate the derivatives. A general formula of the compact finite-difference scheme using two neighboring grids is

(36) |

for . ( is a coefficient of . See (14).) Information at instead of should be used for . (In the following discussions, we show equations for only.) is the space-centering parameter specified by in the code. We fix in the following discussion, and therefore second-order accuracy is achieved.

Combining the trapezoidal rule for the time derivatives with the compact
scheme for the space derivatives^{6}

(37) |

where denotes the average value of the variables at and grids, and similarly for . Then, the gyrokinetic equation is cast into the following symbolic form:

(38) |

#### Kotschenreuther’s Green’s function approach

Kotschenreuther et al. developed an efficient way to solve the gyrokinetic equation by breaking the large matrix to be inverted into a number of small matrices KotschenreutherRewoldtTang_95 (). In the method, fields at the future timestep are obtained using a Green’s function formalism to decouple parts of the matrix related to velocity-space integrations in the gyrokinetic equation. Consequently, the matrix to be inverted becomes at each velocity grid for each species. (Due to the periodic boundary condition, modes are not independent.) The actual scheme Belli_06 () that was originally implemented in GS2, and has been inherited by AstroGK, differs slightly from that described in the original paper by Kotschenreuther et al., as described below.

Because (38) is linear with respect to variables at timestep , the solution to the equation may consist of any linear combination of solutions to parts of the equation. Thus we may split the solution at timestep into two pieces, , each satisfying the following equations (here we ignore the spatial index in the interest of clarity):

(39) | ||||

(40) |

where . An inhomogeneous piece depends only on the known quantities at timestep , thus is immediately solved, while a homogeneous piece depends on the fields at timestep . The unknown portion of the fields may be solved as a separate step using a Green’s function approach.

A formal homogeneous solution is given by

(41) |

where , , and are called the plasma response matrices. Using the response matrices, the field Eqs. (16)–(18) can be written as

(42) |

where

(43) | ||||

(44) | ||||

(45) | ||||

(46) | ||||

(47) | ||||

(48) | ||||

(49) | ||||

(50) | ||||

(51) |

with being the identity matrix and

(52) | ||||

(53) | ||||

(54) |

By solving the field Eq. (42), we obtain , and ultimately . Given , we solve (38) for . This is equivalent to solving (39) with replaced by . Finally, the full procedure to solve the gyrokinetic equation becomes:

Up to this point in the subsection, we have ignored the index for notational simplicity. In fact, the gyrokinetic equation and the field equations should be solved at all grids simultaneously. The field vector in (42) is actually a vector of length , where is the number of fields evolved. An electrostatic simulation requires only the evolution of and so has , whereas a fully electromagnetic simulation evolves , , and and so has . Each of the elements of the matrix in (42) represents an matrix, so the entire matrix is of size . Similarly, the vector is of length .

To evaluate the computational efficiency of Kotschenreuther’s approach, let us compare it to a brute-force approach to solving the GK-M system. A brute-force approach requires the inversion of a dense -size square matrix (where is the number of species), which generally takes operations. On the other hand, Kotschenreuther’s method requires the following: (a) for the gyrokinetic equation, inversions of bi-diagonal -size square matrices, which costs ; and, (b) for the fields, the inversion of the matrix in (42), a dense -size square matrix. For a fixed timestep , however, the matrix does not change, so this matrix inversion need only be performed once during an initialization stage. Therefore, during each timestep the field solver requires only an matrix multiplication operation. Ignoring the factor since , Kotschenreuther’s approach requires operations per timestep, much more efficient than the brute-force approach.