GYRO Technical Guide Manual

User Manual:

Open the PDF directly: View PDF .
Page Count: 74

Download
Open PDF In Browser	View PDF

GYRO Technical Guide
J. Candy and E. Belli
General Atomics, P.O. Box 85608, San Diego, CA 92186-5608, USA

December 3, 2018

Contents
1 Navigation of Key Results

1

2 Flux-surface Geometry
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Field-aligned coordinates and flux functions . . . . . . . . .
2.2.1 Clebsch representation . . . . . . . . . . . . . . . . .
2.2.2 Periodicity constraints . . . . . . . . . . . . . . . . .
2.2.3 Toroidal and poloidal flux . . . . . . . . . . . . . . .
2.2.4 Flux-surface averaging and the divergence theorem .
2.2.5 Additional flux functions . . . . . . . . . . . . . . .
2.3 Local Grad-Shafranov equilibria . . . . . . . . . . . . . . . .
2.3.1 Required operators . . . . . . . . . . . . . . . . . . .
2.3.2 Metric coefficients . . . . . . . . . . . . . . . . . . .
2.3.3 Mercier-Luc coordinates . . . . . . . . . . . . . . . .
2.3.4 Solution of the Grad-Shafranov equation . . . . . . .
2.3.5 Magnetic field derivatives . . . . . . . . . . . . . . .
2.3.6 Calculation of the eikonal function . . . . . . . . . .
2.3.7 Gyrokinetic drift velocity in Mercier-Luc coordinates
2.3.8 Perpendicular Laplacian in Mercier-Luc coordinates
2.3.9 Coriolis drift terms . . . . . . . . . . . . . . . . . . .
2.3.10 Detailed catalogue of shape functions . . . . . . . .
2.3.11 Required operators in terms of shape functions . . .
2.3.12 Pressure-gradient effects in the drift velocity . . . .
2.4 Specification of the plasma shape . . . . . . . . . . . . . . .
2.4.1 Model flux-surface shape . . . . . . . . . . . . . . . .
2.4.2 General flux-surface shape . . . . . . . . . . . . . . .
2.4.3 A measure of the error . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

2
2
3
3
4
4
5
5
6
7
7
8
10
11
11
12
13
13
13
15
15
16
16
16
17

3 The Gyrokinetic Model
3.1 Foundations and Notation . . . . . . . . .
3.2 Reduction of the Fokker-Planck Equation
3.2.1 Lowest-order constraints . . . . . .
3.2.2 Equilibrium equation and solution
3.2.3 The drift-kinetic equation . . . . .
3.2.4 The gyro-kinetic equation . . . . .
3.3 The Gyrokinetic Equation in Detail . . .
3.3.1 Ordering . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

21
21
21
22
23
23
23
24
25

iii

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

3.4
3.5
3.6
3.7

3.3.2 Rotation and rotation shear parameters . . . . . . .
3.3.3 Comment on the Hahm-Burrell shearing rate . . . .
Maxwell equations . . . . . . . . . . . . . . . . . . . . . . .
Transport Fluxes and Heating . . . . . . . . . . . . . . . . .
3.5.1 Ambipolarity and Exchange Symmetries . . . . . . .
Entropy production . . . . . . . . . . . . . . . . . . . . . . .
Simplified fluxes and field equations with operator notation
3.7.1 Operator notation . . . . . . . . . . . . . . . . . . .
3.7.2 Maxwell equations: Ha -form . . . . . . . . . . . . .
3.7.3 Maxwell equations: ha -form . . . . . . . . . . . . . .
3.7.4 Transport coefficients . . . . . . . . . . . . . . . . .

4 Normalization of Fields and Equations
4.1 Dimensionless fields and profiles . . . . . . . . .
4.2 Velocity space normalization . . . . . . . . . .
4.2.1 Velocity variables . . . . . . . . . . . . .
4.2.2 Dimensionless velocity-space integration
4.3 Dimensionless equations . . . . . . . . . . . . .
4.3.1 Normalized gyrokinetic equation . . . .
4.3.2 Normalized Maxwell equations . . . . .
4.3.3 Normalized Transport Fluxes . . . . . .
4.3.4 Diffusivities . . . . . . . . . . . . . . . .
4.3.5 GyroBohm normalization . . . . . . . .

.
.
.
.
.
.
.
.
.
.

5 Spatial Discretization
5.1 Foreword . . . . . . . . . . . . . . . . . . . . . .
5.2 Spectral Decomposition in Toroidal Direction . .
5.2.1 Expansion of fields . . . . . . . . . . . . .
5.2.2 Poloidal wavenumber . . . . . . . . . . . .
5.3 Operator Discretization Methods . . . . . . . . .
5.3.1 Finite-difference operators for derivatives
5.3.2 Upwind schemes . . . . . . . . . . . . . .
5.3.3 Banded pseudospectral gyro-orbit integral
5.3.4 Banded Approximations . . . . . . . . . .
5.4 Discretization of the gyrokinetic equation . . . .
5.4.1 Parallel motion on an orbit-time grid . . .
5.4.2 Er shear . . . . . . . . . . . . . . . . . . .
5.4.3 Drift motion . . . . . . . . . . . . . . . .
5.4.4 Diamagnetic effects . . . . . . . . . . . . .
5.4.5 Poisson bracket nonlinearity . . . . . . . .
5.5 Blending-function expansion for fields . . . . . .
5.6 Velocity-Space Discretization . . . . . . . . . . .
5.6.1 Decomposition of FV . . . . . . . . . . .
5.6.2 Energy Integration . . . . . . . . . . . . .
5.6.3 λ Integration . . . . . . . . . . . . . . . .
5.6.4 Discretization Summary . . . . . . . . . .
5.7 Connection to Ballooning Modes . . . . . . . . .

General Atomics Report GA-A26818

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
operators
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

25
25
26
27
27
27
28
28
29
29
29

.
.
.
.
.
.
.
.
.
.

30
30
31
31
31
32
32
32
33
33
33

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

34
34
34
34
35
35
35
35
36
38
38
39
40
40
41
41
42
45
45
45
46
47
47

iv

6 Temporal
6.0.1
6.0.2
6.0.3
6.0.4
6.0.5

Discretization
Reduction to canonical form . . . . . . . . . . .
IMEX-RK-SSP schemes of Pareschi and Russo
Implementation in GYRO . . . . . . . . . . . .
Procedural summary . . . . . . . . . . . . . . .
Time-Integration Considerations . . . . . . . .

7 Source and Boundary Conditions
7.1 Radial Domain and Boundary Conditions
7.1.1 Periodic . . . . . . . . . . . . . . .
7.1.2 Nonperiodic . . . . . . . . . . . . .
7.2 Long-wavelength Source . . . . . . . . . .
7.2.1 Formulation of the problem . . . .
7.2.2 Solution by damping . . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

49
49
50
51
53
53

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

54
54
54
54
56
56
56

8 Collisions
8.1 Pitch-angle Scattering Operator . . . . . . . . . . .
8.1.1 The Radial Basis Function (RBF) Method .
8.1.2 Basic RBF expansion . . . . . . . . . . . .
8.1.3 Influence of boundaries . . . . . . . . . . .
8.1.4 The method in detail . . . . . . . . . . . . .
8.1.5 Additional comments . . . . . . . . . . . .
8.2 Conservative Krook Operator . . . . . . . . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

58
58
59
60
60
60
61
62

9 Maxwell Dispersion Matrix Eigenvalue Solver
9.1 Motivation and Related Work . . . . . . . . . .
9.2 The Linearized Gyrokinetic Equation . . . . . .
9.3 Construction of the Dispersion Matrix . . . . .
9.4 Discretization and Implementation in LAPACK

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

63
63
64
64
65

General Atomics Report GA-A26818

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.

v

Chapter 1

Navigation of Key Results
The basic equations, simplified but not yet normalized:
Gyrokinetic equation: Eq. (3.37).
Maxwell equations: Eq. (3.76), (3.77) and (3.78)
Transport coefficients: Eq. (3.79), (3.80), (3.81) and (3.82)

1

Chapter 2

Flux-surface Geometry
2.1

Introduction

The goal of this chapter is to present a self-contained derivation and cataloguing of all geometrical
quantities required for gyrokinetic or neoclassical simulation of plasmas with an arbitrary crosssectional shape. The method represents a refinement and generalization of a class of approaches
which fall under the category of local equilibrium techniques [ML74]. These have a long history of
application to ballooning mode analysis [GC81, BKC+ 84, MCG+ 98], and are described in historical
detail by Miller [MCG+ 98]. The unified method we propose can be consistently applied to either
a local or global equilibrium, and to flux-surfaces of model or general shape. A global equilibrium,
in this context, refers to a preexisting numerical solution of the Grad-Shafranov equation over
the entire plasma volume, whereas a local equilibrium refers to Grad-Shafranov force-balance at a
single surface only, without reference to adjacent surfaces. The latter scenario has proven to be a
powerful tool for understanding the effect of plasma shape [KWC07, BHD08] on transport.
The unified approach requires that the major radius R(Ψt , θ) and elevation Z(Ψt , θ) along a
given flux-surface contour are tabulated as functions of the poloidal angle θ (or, more generally,
any parameter that labels position on the flux surface at constant toroidal angle). In the local case,
the toroidal flux, Ψt , need not be specified, and geometrical quantities can be evaluated up to an
unknown constant (the so-called effective field, Bunit ), with flux-surfaces labeled by the midplane
minor radius, r. Both Bunit and r will be precisely defined later. If a global plasma equilibrium
exists, such that the toroidal flux, Ψt , has been computed, then the functions r = r(Ψt ) and Bunit (r)
can be uniquely determined for all r (i.e., on each flux-surface). A special case, in which R and
Z are approximated by up-down symmetric contours with finite ellipticity and triangularity (the
so-called Miller geometry method [MCG+ 98]), is treated by numerous codes worldwide [DJKR00,
CW03, CPC+ 03]. Unfortunately, the implementations of this method are not standardized, and
documentation is typically unavailable. For example, in certain cases, the model shape is assumed
but Grad-Shafranov force balance is not enforced. This gives rise to implementation differences
which significantly complicate code benchmarking exercises. The challenges are further amplified
for the consistent treatment of general (numerically-generated) flux-surface shape corresponding to
global plasma equilibria [XMJ+ 08]. Then, not only is there an even greater liklihood of significant
implementation difference, but the connections to the local limit and to the limit of model shape
are obscured. The present work is an attempt to define a standard approach which unifies the
treatment of local and global equilibria, as well as model and general flux-surface shape. By using
a general Fourier-series expansion for the flux-surface shape, we show how to estimate the error in
approximating a general closed contour with a model shape (often referred to as the Miller shape

2

[MCG+ 98, WM99]). We also show how to consistently identify the local limit of a global equilibrium
on each flux surface. The results herein modify certain conventions in the local equilibrium theory
and, further, correct a variety of sign and typographical errors in an earlier paper by Waltz and
Miller [WM99].
In Sec. 2.2, we give the most general definitions of coordinates and associated flux functions. In
Sec. 2.3, we use the local equilibrium method to compute all geometrical quantities required for gyrokinetic or neoclassical simulation of plasmas, given (R, Z) coordinates and associated derivatives
on a flux-surface. We specialize these results, in Sec. 2.4, to the specific cases of model (Miller) and
general (Fourier series) flux-surface parameterizations, and provide error estimates for each method.
Finally, in the appendix, we give a large-aspect-ratio expansion of the geometry coefficients, useful
for purposes of code checking, and to make contact with the s-α model [CHT78].

2.2
2.2.1

Field-aligned coordinates and flux functions
Clebsch representation

In what follows we adopt the right-handed, field-aligned coordinate system (ψ, θ, α) together with
the Clebsch representation [KK58] for the magnetic field
B = ∇α × ∇ψ

such that

B · ∇α = B · ∇ψ = 0 .

(2.1)

The angle α is written in terms of the toroidal angle ϕ as
.
α = ϕ + ν(ψ, θ) .

(2.2)

In Eqs. (2.1) and (2.2), ψ (as we will show) is poloidal flux divided by 2π, and θ refers simultaneously
to (a) an angle in the poloidal plane (at fixed ϕ), or (b) a parameterization of distance along a field
line (at fixed α). In these coordinates, the Jacobian is
.
Jψ =

1
1
=
.
∇ψ × ∇θ · ∇α
∇ψ × ∇θ · ∇ϕ

(2.3)

Since the coordinates (ψ, θ, α) and (ψ, θ, ϕ) form right-handed systems, the Jacobian Jψ is positivedefinite. In the latter coordinates, the magnetic field becomes
B = ∇ϕ × ∇ψ +

∂ν
∇θ × ∇ψ
∂θ

Using the definition of the safety factor, q(ψ), we may deduce

Z 2π
Z 2π 
B · ∇ϕ
1
∂ν
ν(ψ, 0) − ν(ψ, 2π)
. 1
q(ψ) =
dθ =
−
dθ =
.
2π 0 B · ∇θ
2π 0
∂θ
2π

(2.4)

(2.5)

For concreteness, we choose the following boundary conditions for ν:
ν(ψ, 2π) = − 2π q(ψ) ,
ν(ψ, 0) = 0 .

(2.6)
(2.7)

By writing B in the standard form for up-down symmetric equilibria,
B = ∇ϕ × ∇ψ + I(ψ)∇ϕ ,
General Atomics Report GA-A26818

(2.8)
3

we can derive the following integral for ν:
ν(ψ, θ) = −I(ψ)

Z
0

θ

Jψ |∇ϕ|2 dθ .

(2.9)

We remark that in the case of concentric (unshifted) circular flux surfaces, one will obtain the
approximate result ν(ψ, θ) ∼ −q(ψ)θ.

2.2.2

Periodicity constraints

The periodicity requirements for physical functions are a potential source of confusion. Consider a
physical function – for example, potential or density – represented by the (real) field z(ψ, θ, α). In
terms of the physical angle ϕ, we have
z(ψ, θ, α) = z (ψ, θ, ϕ + ν[ψ, θ])

(2.10)

We impose the following topological requirements on the function z:
1. z is 2π/∆n-periodic in ϕ for fixed ψ and θ,
2. z is 2π-periodic in θ for fixed ψ and ϕ.
From this point onwards, we will refer to ψ and χt simply and the poloidal and toroidal fluxes,
respectively.

2.2.3

Toroidal and poloidal flux

We can start from the general forms of the toroidal and poloidal fluxes [DHCS91]
ZZ
ZZZ
1
.
Ψt =
B · dS =
B · ∇ϕ dV ,
2π
S
V
Z Zt Z
Z tZ
1
.
Ψp =
B · ∇θ dV .
B · dS =
2π
Sp

(2.11)

(2.12)

Vp

Explicitly inserting the field-aligned coordinate system of the previous section, and differentiating
these with respect to ψ, gives
Z 2π
Z 2π
dΨt
1
=
dϕ
dθ B · ∇ϕ Jψ ,
(2.13)
dψ
2π 0
0
Z 2π
Z 2π
1
B · ∇ϕ
=
dϕ
dθ
,
(2.14)
2π 0
B · ∇θ
0
= 2π q(ψ) ,
(2.15)
Z 2π
Z 2π
dΨp
1
=
dϕ
dθ B · ∇θ Jψ ,
dψ
2π 0
0
Z 2π
Z 2π
1
=
dϕ
dθ ,
2π 0
0
= 2π .
General Atomics Report GA-A26818

(2.16)
(2.17)
(2.18)
4

Thus, ψ is the poloidal flux divided by 2π. For this reason, it is useful to also define the toroidal
flux divided by 2π:
. 1
χt =
Ψt .
(2.19)
2π
According to these conventions,
dΨt = q dΨp

2.2.4

and dχt = q dψ .

(2.20)

Flux-surface averaging and the divergence theorem

It is also convenient to define a related Jacobian
Jr =

∂ψ
1
=
Jψ ,
∇r × ∇θ · ∇ϕ
∂r

(2.21)

where r is the midplane minor radius. The flux-surface average of a function f can be expressed as
I
I
. 1
0
F(f ) = 0 dθ dϕ Jr f where V (r) = dθ dϕJr .
(2.22)
V
The normalizing factor V 0 is evidently the derivative of the flux-surface volume with respect to
r. The element of flux-surface volume can be written as dV = dn dS = dr dθ dϕ Jr where dS is
the element of surface area on a flux-surface, and n is the length along the normal vector n. The
relation between n and r is given by dr = (∂r/∂n) dn = |∇r| dn. This allows one to write a form of
the divergence theorem useful for the derivation of transport equations. Integrating the divergence
of an arbitrary vector A over volume gives
Z
I
I
dS
dV ∇ · A = dS A · n =
A · ∇r = V 0 (r) F (A · ∇r) ,
(2.23)
|∇r|
where we have used
F(f ) =

2.2.5

1
V0

I

dS
f
|∇r|

and V 0 (r) =

I

dS
.
|∇r|

(2.24)

Additional flux functions

Recall that we have assumed that the coordinates (R, Z) of a each flux-surface are known via a onedimensional parameterization (the arc length, for example). As a robust measure of the flux-surface
elevation, we use the elevation of the centroid, defined as
I
dZ RZ
.
Z0 = I
.
(2.25)
dZ R
In terms of the centroid elevation, the effective minor radius, r, is defined as the half-width of the
flux surface at the elevation, Z0 :
. R+ − R−
r=
.
(2.26)
2
The quantities R+ and R− are the points of intersection of the flux surface with the line Z = Z0 ,
as illustrated in Fig. 2.1. We further define the effective major radius, R0 = (R+ + R− )/2, and the
effective field strength, Bunit , as
q dψ
. 1 dχt
Bunit =
=
.
(2.27)
r dr
r dr
General Atomics Report GA-A26818

5

120
90
60

Z (cm)

30
0
−30
−60

(R− , Z0 )

r (R0 , Z0 ) r

(R+ , Z0)

Z
ϕ

−90
−120

100

200
R (cm)

300

Figure 2.1: Illustration of the centroid elevation, Z0 , the effective major radius, R0 , and the effective
minor radius, r, for a near-separatrix flux surface in DIII-D. The orientation of the toroidal angle,
ϕ, is also indicated, such that (R, Z, ϕ) is a right-handed system.

The concept of an effective field was introduced by Waltz [WM99], and represents the equivalent
field that would be obtained if the flux surface was deformed to a circle with the penetrating flux
held fixed. It is emphasized that one must have access to a global equilibrium, not just the local
flux-surface shape, to determine Bunit . An important feature of the unified approach is that the
definitions of r and Bunit are the same in both cases; that is, for both the model Grad-Shafranov
(i.e., Miller) and general Grad-Shafranov equilibria to be defined shortly.
As an alternative to working with the toroidal flux directly, we introduce a function ρ, with
units of length, which parameterizes the toroidal flux:
Z
1
1
χt =
B · dS = Bref ρ(r)2 ,
(2.28)
2π
2
In the expression above, Bref is some reference field. While this would normally be the vacuum
toroidal field, or something similar, it is important to keep in mind that the choice is ultimately
arbitrary, but by knowing Bref and ρ, one may calculate χt . In terms of ρ, the area enclosed by a
flux surface is approximately given by A ' πρ2 so long as Bref is approximately equal to the on-axis
toroidal field strength. In this sense, ρ is more intuitive, and so in certain cases more convenient,
than χt .

2.3

Local Grad-Shafranov equilibria

The aim of this section is to compute convenient, standardized expressions for the common differential and integral operators appearing in the drift-kinetic and gyrokinetic equations.

General Atomics Report GA-A26818

6

2.3.1

Required operators

Below we collect, in coordinate-free form, the full set of differential and integral operators which
are to be evaluated. In what follows, Ωca = za eB/(ma c) is the gyroradius of species a and p =
P
a na kB Ta is the total plasma pressure.
Perpendicular drift:
vd · ∇ ⊥ =

vk2

2 /2
+ v⊥
Ωca B 2

B × ∇B · ∇⊥ +

4πvk2
b
Ωca B 2

Perpendicular Laplacian:
Radial gradient of eikonal:
Binormal gradient of eikonal:
Parallel gradient:
Flux-surface average:

× ∇p · ∇⊥ +

2vk ω0
b × s · ∇⊥ ,
Ωca

∇2⊥ = (∇ − bb · ∇)2 ,

(2.29)
(2.30)

∇ψ
∇ψ · ∇ν
· ∇α =
,
|∇ψ|
|∇ψ|

(2.31)

b×

(2.32)

∇ψ
B
· ∇α = −
,
|∇ψ|
|∇ψ|

1 ∂
,
Jψ B ∂θ
 I

dV −1
dθ dϕ Jψ f .
hf i =
dψ
b·∇=

(2.33)
(2.34)

The perpendicular drift operator, Eq. (2.29), is a sum of curvature and grad-B terms. This particular form, which neglects the parallel part of the drift response (see, for example, Eq. (54) of
Ref. [HW06] or Eq. (12) of Ref. [BC08]), is sufficient to recover all standard results of gyrokinetic
and neoclassical theory. In Eq. (2.34), hf i denotes the flux-surface average of the function f , and
I
dV
= dθ dϕ Jψ ,
(2.35)
dψ
where V is the volume enclosed by the flux surface. In the drift velocity, s is the following dimensionless vector
1 ∂R
I
s=
eϕ −
∇R .
(2.36)
Jψ B ∂θ
RB
which is discussed in more detail later.

2.3.2

Metric coefficients

Define a right-handed Cartesian coordinate system (x, y, z) via intermediate (R, Z, ϕ) coordinates
as
x = R(r, θ) cos(−ϕ) ,

(2.37)

y = R(r, θ) sin(−ϕ) ,

(2.38)

z = Z(r, θ) .

(2.39)

In this geometry the covariant basis vectors are
. ∂r
êr =
∂r

,

General Atomics Report GA-A26818

. ∂r
êθ =
∂θ

and

. ∂r
êϕ =
,
∂ϕ

(2.40)

7

where r = (x, y, z). The corresponding contravariant basis vectors are
.
êr = ∇r

,

.
êθ = ∇θ

and

.
êϕ = ∇ϕ .

(2.41)

With this information we can write the covariant and contravariant components of the metric tensor
g as gij = êi · êj and g ij = êi · êj . It is illustrative to write this in matrix form as


grr grθ
0
(2.42)
gij =  grθ gθθ 0  ,
0
0 gϕϕ
where



∂R 2
∂Z 2
=
+
,
∂r
∂r
∂R ∂R ∂Z ∂Z
=
+
,
∂r ∂θ
∂r ∂θ




∂Z 2
∂R 2
+
.
=
∂θ
∂θ


grr
grθ
gθθ

(2.43)
(2.44)
(2.45)

gϕϕ = R2

(2.46)

Since we also know that gij · g ij = I, where I is the identity matrix, we can easily determine the
contravariant counterpart by calculating the inverse of gij . This yields


gθθ gϕϕ −grθ gϕϕ
0
1

0
(2.47)
g ij = 2  −grθ gϕϕ grr gϕϕ
Jr
2
0
0
grr gθθ − grθ


(∇r)2 ∇r · ∇θ
0
 .
0
=  ∇r · ∇θ (∇θ)2
(2.48)
2
0
0
(∇ϕ)
Here, the Jacobians are
∂ψ
. ∂(x, y, z)
Jr =
=
Jψ
∂(r, θ, ϕ)
∂r

. ∂(x, y, z)
and Jψ =
,
∂(ψ, θ, ϕ)

(2.49)

where explicitly

Jr = det gij = R

∂R ∂Z
∂R ∂Z
−
∂r ∂θ
∂θ ∂r


>0.

Metric coefficients can then be computed straightforwardly. For example,


∂R ∂Z −1
1/2 ∂R ∂Z
.
|∇r| = gθθ
−
∂r ∂θ
∂θ ∂r

2.3.3

(2.50)

(2.51)

Mercier-Luc coordinates

Consider a reference flux surface ψ = ψs and the corresponding one-dimensional curve x(`) =
(Rs , Zs ) defined by the intersection of that surface and the plane ϕ = 0. If we choose ` to be the
arc length along x, then the tangent vector


dRs (`) dZs (`)
. dx
t=
=
,
(2.52)
d`
d`
d`
General Atomics Report GA-A26818

8

Z

(R, Z)
̺
u

(Rs , Zs )
ℓ

R

Figure 2.2: Mercier-Luc coordinate system.

is a unit vector. We can further define the unit tangent and unit (inward) normal vectors in terms
of a frame angle u [Gug77] as
.
t = (− sin u, cos u)

.
and n = (− cos u, − sin u) .

(2.53)

This definition of the frame angle is different than that used in Ref. [MCG+ 98] and [WM99]. We
believe the present convention is the natural choice, since in the circular limit we have u = θ and
` = rθ. The radius of curvature, rc , satisfies the equation
dt
n
=
d`
rc

so that

1
du
=
.
rc (`)
d`

Some algebra shows that the curvature can be written as

−1
2
∂Z ∂ 2 R
3/2 ∂R ∂ Z
rc (θ) = gθθ
−
.
∂θ ∂θ2
∂θ ∂θ2

(2.54)

(2.55)

Following Mercier and Luc [ML74], we introduce a right-handed, orthogonal coordinate system
(%, `, ϕ) which is defined in relation to the reference flux surface ψ = ψs through
R(r, θ) = Rs (`) + % cos u ,

(2.56)

Z(r, θ) = Zs (`) + % sin u .

(2.57)

The orientations of ` and u are shown in Fig. 2.2. The metric tensor in this case is diagonal, with
g%% = 1 ,


% 2
,
g`` = 1 +
rc
gϕϕ = R2 .
The Jacobian is



∂(x, y, z)
1
%
J% =
=
=R 1+
>0.
∂(%, `, ϕ)
∇% × ∇` · ∇ϕ
rc

General Atomics Report GA-A26818

(2.58)
(2.59)
(2.60)

(2.61)
9

By differentiating Eqs. (2.56) and (2.57) with respect to r, and evaluating the result at % = 0, we
obtain the following identities valid on the surface ψ = ψs :
∂ρ
∂r
∂`
∂r

= cos u

∂R
∂Z
+ sin u
,
∂r
∂r

(2.62)

= cos u

∂Z
∂R
− sin u
.
∂r
∂r

(2.63)

ψ=ψs

ψ=ψs

Differentiation with respect to θ yields
∂ρ
∂θ
∂`
∂θ

=0,

(2.64)

ψ=ψs

s
=
ψ=ψs

∂Z
∂θ

2


+

∂R
∂θ

2
=

√

gθθ .

(2.65)

On the surface, the Jacobian is equal to
Jr =

2.3.4

R ∂`
.
|∇r| ∂θ

(2.66)

Solution of the Grad-Shafranov equation

The solution of Grad-Shafranov equation determines the poloidal flux, ψ, as a function of sources
f and p:


∇ψ
2
R ∇·
= −4πR2 p0 (ψ) − II 0 (ψ) .
(2.67)
R2
The lefthand side can be written as
R2 ∇ ·



∇ψ
R2



= ∇2 ψ −

2
∇R · ∇ψ .
R

(2.68)

We wish to obtain a solution valid in a neighborhood of % = 0 and so expand
ψ = ψs + ψ1 (`) % + ψ2 (`) %2 + · · · .

(2.69)

The required operators can then be evaluated to leading order in % as
∇R =

1 ∂R
1 ∂R
ê% +
ê` ∼ ê% cos u − ê` sin u + O(%) ,
h% ∂%
h` ∂`

(2.70)

1 ∂ψ
1 ∂ψ
ê% +
ê` ∼ ψ1 ê% + O(%) ,
h% ∂%
h` ∂`

(2.71)

∇ψ =

 



1
∂ h` hϕ ∂ψ
∂ h% hϕ ∂ψ
∇ ψ=
+
,
h% h` hϕ ∂%
h% ∂%
∂`
h` ∂%
∂ 2 ψ ∂ψ 1 ∂
∼
+
(h` hϕ ) + O(%) ,
∂%2
∂% h` hϕ ∂ρ
2

. √
where hi = gii . Combining these results gives


∇ψ
ψ1
ψ1
2
R ∇·
∼
−
cos u + 2ψ2 + O(%) .
2
R
rc
Rs
General Atomics Report GA-A26818

(2.72)
(2.73)

(2.74)
10

The solution for ψ1 is obtained from
Bp2

1
= |∇ϕ| |∇ψ| ∼ 2
Rs
2

2



∂ψ
∂ρ

2
+ O(%) ,

(2.75)

from which it is evident that
ψ1 = Rs Bps ,

(2.76)

where Bps = Bp (ψs ). The solution of the Grad-Shafranov equation then gives an explicit expression
for ψ2 :


Bp R
1
2 0
0
Bp cos u −
− 4πR p − II
.
(2.77)
ψ2 =
2
rc
ψ=ψs

2.3.5

Magnetic field derivatives

We can start with the exact expressions
" 

2 #
2
1
∂ψ
1
∂ψ
Bp2 = 2
+
,
R
∂%
h` ∂`


I(ψ) 2
2
.
Bt =
R

(2.78)
(2.79)

Expanding the poloidal flux gives
R2 Bt2 ∼ Is2 + 2%ψ1 Is Is0 + O(%2 ) ,

(2.80)

R2 Bp2 ∼ ψ12 + 4%ψ1 ψ2 + O(%2 ) ,

(2.81)

where Is = I(ψs ). Expanding R and taking derivatives gives


∂ 2
1
Is2
0
B ∼ 2 −2
cos u + 2Is Is ψ1 ,
∂% t
Rs
Rs


∂ 2
1
ψ2
Bp ∼ 2 −2 1 cos u + 4ψ1 ψ2 .
∂%
Rs
Rs

(2.82)
(2.83)

Adding these contributions together gives, finally:
∂ 2
B ∼
∂%

2.3.6

Bp2
I2
−2 3 cos u − 2
− 8πRBp p0
R
rc

!
.

(2.84)

ψ=ψs

Calculation of the eikonal function

In order to compute radial derivatives required for the evaluation of the drift and Laplacian, we
need to determine the radial dependence of ν. Using the expansion, Eq. (2.69), for the poloidal
flux, together with the eikonal expansion
ν = νs + ν1 (`) % + · · · ,

(2.85)

∇ν × ∇ψ = f ∇ϕ

(2.86)

we solve the equation

General Atomics Report GA-A26818

11

order-by-order in %. First, an expansion of the gradients yields
∂ψ1
∇ψ ∼ (ψ1 + 2%ψ2 ) ∇% + %
∇` + O(%2 ) ,
∂`


∂ν1
∂νs
∇` + O(%2 ) .
+%
∇ν ∼ ν1 ∇% +
∂`
∂`
Replacing these formulae into Eq. (2.86) and dotting with ∇ϕ yields




∂νs
Is + %Is0 ψ1
%
∂ψ1 ∂ν1
∂νs
−
1+
.
ψ1 + % ν1
−
ψ1 − 2
ψ2 ∼
∂`
∂`
∂`
∂`
Rs + % cos u
rc

(2.87)
(2.88)

(2.89)

At each order in %, we have
O(1) :
O(%) :

∂νs
Is
,
ψ1 =
∂`
Rs
 
∂ ν1
I0
Is
Is cos u
∂νs ψ2
− ψ1
+ s +
−
.
=2
∂` ψ1
∂` ψ1 Rs Rs rc ψ1
Rs2 ψ1

−

Thus, the solution for νs is
νs (`) = −

Z
0

`

d`0
Is .
Rs2 Bps

(2.90)
(2.91)

(2.92)

Some additional algebra shows that


ν1 (`) = Rs (`)Bps (`) D0 (`) + D1 (`) II 0 + D2 (`) p0 ,

(2.93)

where the Di integrals are:
D0 = −
D1 = −
D2 = −

Z

`

0

Z

`

0

Z
0

`

d`0
Rs2 Bps
d`0
Rs2 Bps
d`0
Rs2 Bps



2
2 cos u
−
rc Rs Bps Rs2 Bps
 2
B
1
,
2
Bps
Is


4π
Is .
2
Bps


Is

(2.94)
(2.95)
(2.96)

.
We want to eliminate II 0 in favour of q 0 . Writing `(2π) = L, and expanding Eq. (2.6), we find
− 2π(qs + qs0 %ψ1 ) ∼ νs (L) + %ν1 (L) + O(%2 ) .

(2.97)

Therefore, we have the connection
− 2π qs0 = D0 (L) + D1 (L) II 0 + D2 (L) p0 .

2.3.7

(2.98)

Gyrokinetic drift velocity in Mercier-Luc coordinates

In this case, we work out the dominant part of the drift operator acting on a function f (ψ, θ, α).
This limit is appropriate for standard gyrokinetics. We proceed by simplifying Eq. (2.29) in the
Miller equilibrium case by beginning from the Mercier-Luc representation. On the surface ψ = ψs
we have
Bs = Bps ∇` + Is ∇ϕ .
(2.99)
General Atomics Report GA-A26818

12

Note that Eq. (2.99) cannot be used to obtain ∂B/∂ψ because the limit % → 0 has already been
taken. To evaluate ∇B, we must write
∇B =

∂B
∂B
∇% +
∇` ,
∂%
∂`

(2.100)

and then make use of Eq. (2.84). Taking the cross product of the above two quantities gives
Bs × ∇B = −Bps

∂B
∂B
∂B
∇% × ∇` + Is
∇ϕ × ∇% + Is
∇ϕ × ∇` .
∂%
∂%
∂`

Then, the perpendicular gradient operator, accurate to O(%), can be written as


∂
∂νs
∂
∇⊥ ∼ ∇%
+ ∇ϕ +
∇` + ν1 ∇% + O(%)
.
∂%
∂`
∂α

(2.101)

(2.102)

The final result, valid on the surface ψ = ψs , is
B × ∇B · ∇ = −Is Bps

2.3.8



B 2 ∂B Is ν1 ∂B ∂
∂B ∂
+ −
−
.
∂` ∂ψ
Rs Bps ∂%
Rs ∂` ∂α

(2.103)

Perpendicular Laplacian in Mercier-Luc coordinates

To evaluate the perpendicular Laplacian, we square Eq. (2.102) and simplify, where it is sufficiently
accurate to ignore variation of coefficients
"

 #
1
∂2
∂ ∂
∂ν0 2 ∂ 2
2
2
+
+ ν1 +
.
(2.104)
∇⊥ ∼ 2 + 2ν1
∂%
∂% ∂α
R2
∂`
∂α2

2.3.9

Coriolis drift terms

Some algebra show that
∂
b × s · ∇⊥ = − sin u
+
∂ρ




I
∂
cos u − ν1 sin u
,
2
R Bp
∂α

where we have used the relation
∂R
sin u = −
∂θ

2.3.10



∂`
∂θ

−1

.

(2.105)

(2.106)

Detailed catalogue of shape functions

For convenience, we give a complete summary of the equations needed to compute all the relevant
shape functions. In the expressions below, R and Z are taken to be functions of r and θ, and the
subscript s can be safely omitted.
s



∂`
∂R 2
∂Z 2
=
+
,
(2.107)
∂θ
∂θ
∂θ

Jr = R

General Atomics Report GA-A26818

∂R ∂Z
∂R ∂Z
−
∂r ∂θ
∂θ ∂r


,

(2.108)

13

|∇r| =

R ∂`
,
Jr ∂θ

−1
∂R ∂ 2 Z
∂Z ∂ 2 R
−
,
rc (θ) =
∂θ ∂θ2
∂θ ∂θ2
Z L
−1
I(r)
d`0
= 2π r
,
Bunit
0 R|∇r|


∂`
∂θ

3 

I
Bt (r, θ)
,
=
Bunit (r)
RBunit

(2.109)
(2.110)
(2.111)
(2.112)

Bp (r, θ)
r |∇r|
,
=
Bunit (r)
R q
s



Bp 2
Bt 2
B(r, θ)
= sgn(Bunit )
+
,
Bunit (r)
Bunit
Bunit


. Bt R0 ∂B
gsin(r, θ) =
,
B
B ∂`

(2.114)

. R0 ∂B
gcos(r, θ) = −
= gcos1 + gcos2 ,
B ∂%

(2.116)

gcos1 (r, θ) =

Bp2 R0
Bt2 R0
cos
u
+
,
B2 R
B 2 rc

2
1 Bunit
R0 |∇r| β ∗ ,
2 B2
usin(r, θ) = sin u ,

gcos2 (r, θ) = −

Bt
cos u ,
B


Z `(θ)
d`0 Bt r
r
E1 (r, θ) = 2
− cos u
R|∇r| Bp rc R
0
 2
Z `(θ)
d`0
B
E2 (r, θ) =
,
R|∇r| Bp2
0
 2 
Z
1 `(θ) d`0 Bt Bunit
,
E3 (r, θ) =
2 0
R Bp
Bp2


qs 1
1
∗
∗
f (r) =
2π − E1 (r, 2π) + β E3 (r, 2π)
E2 (r, 2π)
r
r




RBp
RBp
1
.
∗
∗
Θ(r, θ) = −
ν1 =
|∇r|
E1 + f E2 − β E3
B
B
r


. 1 rB
Gq (r, θ) =
,
q RBp
ucos(r, θ) =

B R 1 ∂`
. Jψ B
Gθ (r, θ) =
=
.
qR0
Bunit R0 r|∇r| ∂θ

General Atomics Report GA-A26818

(2.113)

(2.115)

(2.117)
(2.118)
(2.119)
(2.120)
(2.121)
(2.122)
(2.123)
(2.124)
(2.125)
(2.126)
(2.127)

14

2.3.11

Required operators in terms of shape functions
Perpendicular drift:
vk2

vd · ∇⊥ = −ikθ Gq
+ikθ Gq

−ikθ Gq

2vk ω0 R0
Ωca R0

vk2

+ µB

!
[gcos1 + gcos2 + Θ gsin]

Ωca R0

vk2 + µB

!

2vk ω0 R0
Ωca R0



!

Ωca R0

gcos2 −



Ωca R0


[ucos + Θ usin] −

|∇r| gsin

∂
,
∂r

|∇r| usin

∂
,
∂r

Perpendicular Laplacian:

∂
− kθ2 G2q 1 + Θ2 ,
∇2⊥ = |∇r|2 2 + 2iΘkθ Gq |∇r|
∂r
∂r
∂2

Radial gradient of eikonal:
Binormal gradient of eikonal:
Parallel gradient:

n

∇ψ
· ∇α = −kθ Gq Θ ,
|∇ψ|

nb ×
b·∇=
Z

∇ψ
· ∇α = −kθ Gq ,
|∇ψ|

1 1 ∂
,
Gθ qR0 ∂θ

(2.128)

(2.129)
(2.130)
(2.131)
(2.132)

2π

Gθ
f (θ)
B
0
hf i = Z 2π
.
Gθ
dθ
B
0
dθ

Flux-surface average:

(2.133)

Above, we have introduced the binormal wavenumber, kθ = i(q/r)∂/∂α.

2.3.12

Pressure-gradient effects in the drift velocity

In cases where compressional magnetic perturbations are neglected, one may also wish to ignore
finite-pressure effects in the drift velocity [WM99]. This is done by setting gcos2 = 0 in Eq. (2.128).

General Atomics Report GA-A26818

15

2.4

Specification of the plasma shape

It is assumed throughout that we have access to a global equilibrium for which the poloidal and
toroidal fluxes are known over the entire range of r. Then, we consider two general ways in which
to specify the plasma shape via a closed contour in the (R, Z) plane: (1) a model parameterization,
and (2) a general Fourier expansion.
In addition to the plasma shape, we also need to know the plasma pressure gradient, which we
will define via the effective inverse beta-gradient scale length:
8π ∂p
.
β ∗ (r) = − 2
>0.
Bunit ∂r

(2.134)

The safety factor q = dψ/dχt and the magnetic shear s = (r/q)dq/dr can be obtained directly via
numerical differentiation.

2.4.1

Model flux-surface shape

A popular model for the flux-surface shape in (R, Z) coordinates was introduced by Miller [MCG+ 98].
We generalize that model to include the effects of elevation and squareness. In particular, finite
elevation is required very close to r = 0, where the flux surface may lie entirely above the midplane
at Z = 0. Specifically, let
R(r, θ) = R0 (r) + r cos(θ + arcsin δ sin θ) ,

(2.135)

Z(r, θ) = Z0 (r) + κ(r)r sin(θ + ζ sin 2θ) ,

(2.136)

where κ is the elongation, δ is the triangularity, ζ is the squareness and Z0 is the elevation
[TLLM+ 99]. In order to evaluate the required quantities ∂R/∂r and ∂Z/∂r, we also need to
know the radial derivatives of R0 , Z0 , δ, κ and ζ. To this end, we define the associated shape
functions
. r ∂κ
. ∂δ
. ∂ζ
sκ =
, sδ = r
, sζ = r
.
(2.137)
κ ∂r
∂r
∂r
Note that this definition of sδ differs slightly from Refs. [MCG+ 98] and [WM99]. Thus, to compute
the shape functions for the flux-surface at r according to this model, we require the 10 parameters


dR0
dZ0
R0 ,
, Z0 ,
, κ, sκ , δ, sδ , ζ, sζ .
(2.138)
dr
dr

2.4.2

General flux-surface shape

In this case the flux-surface shape is an expansion of the form
N
X
 R

1 R
R(r, θ) = a0 (r) +
an (r) cos(nθ) + bR
n (r) sin(nθ) ,
2

Z(r, θ) =

1 Z
a (r) +
2 0

n=1
N
X

 Z

an (r) cos(nθ) + bZ
n (r) sin(nθ) .

(2.139)

(2.140)

n=1

Here, N is an integer which evidently controls the accuracy of the expansion, and is in principle
arbitrary. If (R, Z) are known as functions of arc length, then one can simply set θ = 2π`/L. To
compute the shape functions for the flux-surface at r, we require the 8(N + 1) parameters


R
R
Z
Z
R R Z Z dan dbn dan dbn
an , bn , an , bn ,
,
,
,
.
(2.141)
dr dr dr dr
General Atomics Report GA-A26818

16

2.4.3

A measure of the error

It is beyond the scope of the present paper to give a detailed assessement of the error in neoclassical
or turbulent transport calculations caused by errors in the flux-surface shape. Nevertheless, we can
d
define the error in the flux-surface approximation of the discrete data {Ri , Zi }ni=1
as
q
nd
. 1 X
min [R(θ) − Ri ]2 + [Z(θ) − Zi ]2 .
ε=
θ
nd r

(2.142)

i=1

Figure 2.3 shows the flux-surface fitting error for DIII-D discharge 132010. This discharge is part
of a series of experiments designed to accurately measure profiles across the entire confined plasma,
including the edge barrier region. Discharge 132010 is a standard DIII-D H-mode plasma, with
typical values of the magnetic field (Bt = 2.1T ), current (Ip = 1.2M A), elongation (κ = 1.8)
and triangularity (δ = 0.3). The width of the edge barrier (∼ 3% in r/a), height of the edge
barrier (∼ 11kP a), and global pressure (βN ∼ 2) are all fairly typical of DIII-D H-Modes. The
Grad-Shafranov equillbrium was computed by EFIT using a 129 × 129 mesh, and then mapped to
very-high-resolution (R, Z)-contours (400 flux surfaces, each with 512 points equally-spaced in arc
length along the surface) using the ELITE code. The Miller-type model equilibrium is seen to be
significantly less accurate, according to the metric defined in Eq. 2.142, at all radii than even the
N = 4 Fourier expansion (red curve). The N = 8 (solid black curve) and N = 12 (green curve)
results are also shown. In the region r/a < 0.4, the accuracy of the Fourier expansion probably
exceeds the accuracy of the original equilibrium solution, and is therefore limited by it. This result
shows that N > 12 is probably a good default choice for N . A further comparison of actual flux
surface shapes at r/a = 0.99 (i.e., very close to the separatrix) is shown in Fig. 2.4.
For core plasma parameters, our limited experience is that the difference, with respect to linear
growth rates, between general and model flux-surfaces is insignificant in moderately-shaped DIII-D
plasmas (roughly 5% at r/a = 0.9, and less than that deeper in the core). However, this relies on
the accuracy of the fitting procedure used to obtain the model parameters κ, δ, etc. If the model
flux-surface shape is poorly calculated, the absolute error can be significant.

General Atomics Report GA-A26818

17

100
Parameterized
N =4
N =8
N = 12

10−1

10−2
ε 10−3

10−4

10−5

10−6

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
r/a

1

Figure 2.3: Error, ε, in the flux-surface fits for DIII-D discharge 132010. The Miller-type model
equilibrium is seen to be significantly less accurate at all radii than even the N = 4 Fourier expansion
(red curve). The N = 8 (solid black curve) and N = 12 (green curve) results are also shown. In
the region r/a < 0.4, the accuracy of the Fourier expansion probably exceeds the accuracy of the
original equilibrium solution, and is therefore limited by it.

120

90

60

60

30

30
Z (cm)

Z (cm)

90

120
Parameterized
Exact

0

0

−30

−30

−60

−60

−90

−90

−120

100

200
R (cm)

N =4
N = 12
Exact

−120

100

200
R (cm)

Figure 2.4: Comparison of exact flux-surface (dashed curve) to parameterized (left) and general
Fourier expansion (right) at r/a = 0.99.

General Atomics Report GA-A26818

18

Appendix A: Large-aspect-ratio expansion
Consider a shifted-circle (Shafranov) flux-surface shape:
R(r, θ) = R0 + ∆(r) + r cos θ ,

(2.143)

Z(r, θ) = r sin θ .

(2.144)

Let ∆(r) = ar2 /(2R0 ) be the Shafranov shift, with a an O(1) constant, so that ∂∆/∂r = a (r/R0 ).
It is instructive to calculate the shape functions as Taylor series in the small parameter r/R0 . To
first order in r/R0 , we find
r
|∇r| ∼ 1 − a
cos θ ,
(2.145)
R0
Bt (r, θ)
r
cos θ ,
∼1−
Bunit (r)
R0


Bp (r, θ)
r
r
1 − (a + 1)
∼
cos θ ,
Bunit (r)
qR0
R0

(2.146)
(2.147)

r
B(r, θ)
∼1−
cos θ ,
Bunit (r)
R0


r
cos θ ,
gsin(r, θ) ∼ sin θ 1 −
R0


r
1
2
gcos1 (r, θ) ∼ cos θ −
cos θ − 2 ,
R0
q
r
Θ(r, θ) ∼ sθ − q 2 R0 β ∗ sin θ + Θ1
,
R0
r
cos θ ,
Gq (r, θ) ∼ 1 + (a − 1)
R0
r
cos θ .
Gθ (r, θ) ∼ 1 + a
R0

(2.148)
(2.149)
(2.150)
(2.151)
(2.152)
(2.153)

Here, Θ1 is the function


1
Θ1 = (1 − 2a)sθ sin θ + sin θ (3a − 1)s − 2(1 + a) + (a − 3)q 2 R0 β ∗ cos θ
2


.

(2.154)

Since the formulae do not represent an exact Grad-Shafranov equilibrium, but rather an equilibrium
accurate only first order in r/R0 , they should not be used for simulation purposes. However, they
can be used to provide a convenient asymptotic check of the numerical routines used to compute
the shape functions in the general case.
Finally, for completeness, we can make contact with the popular s-α equilibrium model [CHT78].
First, recall that
8π dp
.
αMHD = q 2 R0 β ∗ = −q 2 R0 2
.
(2.155)
Bunit dr
To recover the usual s-α formulae, we must take the limit r/R0 → 0 in all the expressions above,
except for B, which is taken to be B/Bunit = R0 /R. Evidently, the result is not a Grad-Shafranov
equilibrium (not even to first order in r/R0 ). In practice, simulation results obtained using a
model circular equilibrium (sometimes called a “Miller circle”) can differ significantly from the
corresponding s-α result. For example, Kinsey [KWC07] has shown that the nonlinear electron
energy flux increases by more than a factor of 1.5 when using a Miller circle instead of s-α geometry.
General Atomics Report GA-A26818

19

Appendix B: Translation to GS2 Geometry Variables
The normalizing magnetic field in GS2, BT 0 , is defined as
.
I(ψ) = RBt = RGEO BT 0 ,

(2.156)

where RGEO is a reference radius (input). For simplicity, let’s assume that RGEO = R0 (this appears
to be a typical choice), in which case we can then relate BT 0 to Bunit according to
Bunit
= R0
BT 0



I

−1

Bunit

,

where I/Bunit is given by Eq. (2.111). Explicitly, we can write this as


Z 2π
Bunit
R0
dθ ∂R ∂Z
∂R ∂Z
=
,
−
BT 0
2πr 0 R ∂r ∂θ
∂θ ∂r


sκ 1 dR0 r
' κ 1+
−
.
2
2 dr R0

(2.157)

(2.158)
(2.159)

Whereas Eq. (2.158) is exact, Eq. (2.159) is approximate and is given here only as an intuitive aid.
In practice the required integral is computed to high accuracy by GYRO and the result available
in the normalized GEO interface variable GEO f. Specifically,
Bunit
R0 1
=
.
BT 0
a GEO f

General Atomics Report GA-A26818

(2.160)

20

Chapter 3

The Gyrokinetic Model
3.1

Foundations and Notation

The definitive works in terms of the complete derivation of the gyrokinetic-Maxwell equations for
evolution of fluctuations, and the neoclassical equations for evaluation of collisional transport, are
due to Sugama and coworkers [SH97, SH98]. For historical reasons, the notation and presentation
here differ in some respects from these papers. Nevertheless, what is documented here should be
completely consistent with the formulae given in Ref. [SH98]. We prefer when applicable to use
Gaussian CGS units1 . Roughly speaking, the gyrokinetic approach [AL80, CTB81, FC82] is based
on the assumption that equilibrium quantities are slowly varying, while perturbations are smaller
but more rapidly varying in space. To describe this ordering, it is convenient to introduce the
parameters
r
Te
cs =
ion-sound speed
(3.1)
mi
cs
ρs =
ion-sound gyroradius
(3.2)
Ωci
where Ωci is the ion cyclotron frequency. The derivation of the gyrokinetic equation proceeds by
.
expanding the primitive Fokker-Planck equation in the small parameter, ρ∗ = ρs /a, where a is the
plasma minor radius.

3.2

Reduction of the Fokker-Planck Equation

The details of the derivation of the gyrokinetic (and neoclassical) equations is beyond the scope of
this manual. Still, we attempt to sketch the essential details. The Fokker-Planck equation provides
the fundamental theory for plasma equilibrium, fluctuations, transport. In this section, we use
Sugama’s notation. The FP equations is written as

 ∂ 
ea 
v
∂
+v·∇+
(E + Ê) + × (B + B̂) ·
(fa + fˆa ) = Ca (fa + fˆa ) + Sa
(3.3)
∂t
ma
c
∂v
where fa is the ensemble-averaged distribution, fˆa is the fluctuating distribution, Sa are sources
(beams, RF, etc), and
X
Ca =
Cab (fa + fˆa , fb + fˆb )
(3.4)
b
1

See http://wwwppd.nrl.navy.mil/nrlformulary/NRL FORMULARY 06.pdf

21

is the nonlinear collision operator. The general approach is to separate the FP equation into
ensemble-averaged, A, and fluctuating, F, components:
d
dt
d
F=
dt
A=

where

ens

fa − hCa iens − Da − Sa ,

(3.5)

 ∂
ea 
v
fˆa +
Ê + × B̂ ·
(fa + fˆa ) − Ca + hCa iens + Da ,
ma
c
∂v
ens
d
dt

 ∂
v
ea 
. ∂
=
E+ ×B ·
+v·∇+
,
∂t
ma
c
∂v
ens
*
+
 ∂ fˆ
ea 
v
.
a
Da = −
Ê + × B̂ ·
.
ma
c
∂v

(3.6)

(3.7)

(3.8)

ens

such that Da is the fluctuation-particle interaction operator. Ensemble averages are expanded in
powers of ρ∗ as
fa = fa0 + fa1 + fa2 + . . . ,
Sa =

(3.9)

Sa2 + . . . (transport ordering),

(3.10)

E = E0 + E1 + E2 + . . . ,

(3.11)

B = B0 .

(3.12)

Fluctuations are also expanded in powers of ρ∗ as

3.2.1

fˆa = fˆa1 + fˆa2 + . . . ,

(3.13)

Ê = Ê1 + Ê2 + . . . ,

(3.14)

B̂ = B̂1 + B̂2 + . . . .

(3.15)

Lowest-order constraints

The lowest-order ensemble-averaged equation gives the constraints
A−1 = 0 :

1
E0 + V0 × B = 0
c

and

∂fa0
=0
∂ξ

(3.16)

where ξ is the gyroangle. The only zeroth-order flow (i.e., sonic) flow that persists on the fluctuation
timescale is a purely toroidal flow [HW85], which we write as
V0 = V0 eϕ = R ω0 (ψ)eϕ

where

∂φ−1
.
ω0 = −c
.
∂ψ

(3.17)

The first-order flow V1 contains both toroidal and poloidal components which are self-consistenly
calculated by as moments of the ensemble averaged first-order distribution, fa1 (that is, they are
computed within the context of neoclassical theory).

General Atomics Report GA-A26818

22

3.2.2

Equilibrium equation and solution

The gyrophase average of the zeroth order ensemble-averaged equation gives the collisional equilibrium equation:
Z 2π


dξ
(3.18)
A0 = 0 :
V0 + vk0 b · ∇fa0 = Ca (fa0 )
2π
0
where v0 = v − V0 is the deviation from the mean flow velocity. The exact solution for fa0 is a
Maxwellian in the rotating frame, such that the centrifugal force causes the density to vary on the
flux surface:


ma (v 0 )2
na (ψ, θ)
exp −
fa0 =
= na FM a .
(3.19)
2Ta
(2πTa /ma )3/2
It is important to note here that to account for sonic rotation in the equilibrium, we do not
use the shifted Maxwellian approach [WSCH07]. Although the use of a shifted Maxwellian gives
rise to errors which are probably not significant for typical operating parameters, the approach is
conceptually incorrect. Instead, we work in the shifted velocity frame v0 = v − V0 . The latter
approach was first used in the context of neoclassical transport by Hinton and Wong [HW85].
Beyond this point, we limit our attention to the moderate flow regime, such that ρ∗  V0 /cs  1.
Operationally, this means that we will ignore terms quadratic in V0 /cs , whilst retaining all terms
which are linear in V0 /cs . Physically, this approach will not capture centrifugal effects like the
poloidal variation of density, since
na (ψ, θ) ∼ na (ψ) + O(V02 /c2s ) ,

(3.20)

but will correctly retain the symmetry-breaking effects of radial electric field shear, rotation shear
drive, and the Coriolis drift.

3.2.3

The drift-kinetic equation

Taking a gyroaverage of the first-order ensemble-averaged component A1 gives expressions for the
gyroangle-dependent and independent distributions, f˜a1 and f¯a1 :
Z

2π

dξ
A1 = 0 :
2π

0

fa1

= f˜a1 + f¯a1 ,

1
f˜a1 =
Ωa

Z

ξ

]
dξ Lf
a0

(3.21)

The function f¯a1 is determined by the solution of the drift kinetic equation.

3.2.4

The gyro-kinetic equation

The gyroaverage of first-order F1 gives an expression for first-order fluctuating distribution, fˆa1 , in
terms of the distribution of the gyrocenters, Ha (R):
Z
0

2π

dξ
F1 = 0 :
2π

ea δφ(x)
fˆa1 (x) = −
fa0 + Ha (R) ,
Ta

(3.22)

where x = R + ρ is the particle position, ρ = b × v0 /Ωca is the gyroradius vector, Ωca = ea B/(ma c)
is the cyclotron frequency, and R (X in Ref. [SH98]) is the guiding-center position. The function
Ha (R) (ha (R) in Ref. [SH98]) is determined by solution of the nonlinear gyrokinetic equation. Also
note that the perturbed potential δφ appears as φ̂ in Ref. [SH98].

General Atomics Report GA-A26818

23

3.3

The Gyrokinetic Equation in Detail

In what follows, we will drop the prime notation in reference to the velocity coordinates. To bring the gyrokinetic equation into a form convenient for numerical integration, we
introduce the function ha (R):
Ha (R) =

ea fa0
Ψa (R) + ha (R)
Ta

(3.23)

where Ψa (ψ̂a in Ref. [SH98]) is the following gyrophase average (or, more simply, gyroaverage) at
fixed R:


1
.
Ψa (R) = δφ(R + ρ) − (V0 + v) · δA(R + ρ)
.
(3.24)
c
R
The gyroaverage can be defined formally as
.
hz(R, ξ)iR =

I

dξ
z(R, ξ) ,
2π

(3.25)

for any function, z. Also, we note the following identities
b · ∇V0 = ω0 s
I
b · ∇V0 + ∇V0 · b = ∇ω0 ,
B

(3.26)
(3.27)

where s is the dimensionless vector
s=

1 ∂R
I
eϕ −
∇R .
Jψ B ∂θ
RB

(3.28)

We use a form of the gyrokinetic equation which can be obtained, after some rearrangement, from
Eq. (46) of Ref. [SH98].

∂ha
+ vk b + vd · ∇Ha + vE0 · ∇ha + δva · ∇ha
∂t


ma vk fa0 I
∇ω0
+ δva · ∇fa0 +
Ta
B



= CaGL [Ha ] .

(3.29)

The velocities are
2
4πvk2
2vk ω0
. vk + µB
vd =
b × ∇B +
b×s+
b × ∇p
Ωca B
Ωca
Ωca B 2
. c
vE0 =
b × ∇φ−1 ,
B
. c
δva =
b × ∇Ψa .
B

(3.30)
(3.31)
(3.32)

In this result, s is a dimensionless vector which is not exactly in the grad-B direction. The correct
form of s cannot be obtained using the shifted Maxwellian model. This form of the drift matches
Brizard’s result, and the resulting expression for vd · ∇ψ matches the familiar result from Hinton
and Wong [HW85]. By setting ω0 = 0, we recover the usual diamagnetic rotation ordering. For
normalization purposes it is useful to note that
Z
d 3v FM a = 1 .
(3.33)
General Atomics Report GA-A26818

24

Various terms can be simplified without refering to the geometry model. Within the accuracy of
the gyrokinetic ordering, we have
∂ha
c
b × ∇φ−1 · ∇ha ∼ ω0
,
B
∂α
c
∂fa0 ∂Ψa
δva · ∇fa0 =
b × ∇Ψa · ∇fa0 ∼ c
,
B
∂ψ ∂α
c
∂ha ∂Ψa
∂ha ∂Ψa
δva · ∇ha =
b × ∇Ψa · ∇ha ∼ c
−c
.
B
∂ψ ∂α
∂α ∂ψ
vE0 · ∇ha =

(3.34)
(3.35)
(3.36)

Using these expression, we find
vk ∂Ha
∂ha
∂ha
+
+ vd · ∇Ha + ω0
+ c [ha , Ψa ]ψ,α
∂t
Jψ B ∂θ
∂α


∂Ψa
∂fa0 ma vk I ∂ω0
+c
+
fa0
= CaGL [Ha ] .
∂ψ
Ta B ∂ψ
∂α

(3.37)

The expansion and simplification of the remaining operators has been treated in Chap. 2.

3.3.1

Ordering

We remark that the gyrokinetic ordering requires that
ea Ψa
ha
ω − k⊥ · V 0
∼
∼
∼ kk ρs ∼ ρ∗ .
Ta
fa0
Ωca

3.3.2

(3.38)

Rotation and rotation shear parameters

Recalling the definition of the rotation frequency:
∂φ−1
.
ω0 = −c
,
∂ψ

(3.39)

we define the Mach number, the Er shearing rate and the rotation shearing rate respectively as
. ω0 R0
M=
,
cs
r ∂ω0
.
γE = −
,
q ∂r
∂ω0
.
γp = − R0
.
∂r

(3.40)
(3.41)
(3.42)

These parameters are defined in this way for legacy reasons and are not independent; rather, we
have the constraint
R0
γp =
γE .
(3.43)
qr

3.3.3

Comment on the Hahm-Burrell shearing rate

Note that the shearing rate defined in Eq. (3.41) is not in general equal to the familiar Hahm-Burrell
shearing rate [Bur97]


2
Er
HB . (RBp ) ∂
γE =
,
(3.44)
B ∂ψ RBp
General Atomics Report GA-A26818

25

where Er is the radial electric field
∂φ−1
.
Er = −êr · ∇φ−1 = −|∇r|
.
∂r
In terms of the Miller geometry coefficients, γEHB can be written as


c q ∂φ−1
|∇r|
|∇r| r ∂
HB
=
γE .
γE =
Gq q ∂r Bunit r ∂r
Gq

3.4

(3.45)

(3.46)

Maxwell equations

.
.
Defining the scalar electromagnetic fields δAk = b · δA and δBk = b · ∇ × δA, we can write an
equation for each of the fields (δφ, δAk , δBk ) (see Appendix A, [SH98]). In each case, the species
summation runs over all species a (ions and electrons).
Poisson equation
−

∇2⊥ δφ(x)

= 4π

X
a

eza δna = 4π

X

Z
ea

d 3v fˆa1 (x) .

(3.47)

a

Parallel Ampère’s Law
Z
4π X
4π X
δjk,a =
ea d 3v vk fˆa1 (x) .
c a
c a

(3.48)

Z
4π X
4π X
∇⊥ δBk (x) × b =
δj⊥,a =
ea d 3v v⊥ fˆa1 (x)
c a
c a

(3.49)

− ∇2⊥ δAk (x) =
Perpendicular Ampère’s Law

The righthand sides can be written in terms of Ha according to
Z
Z
na ea
δφ(x) + d 3v Ha (x − ρ) ,
d 3v fˆa1 (x) = −
Ta
Z
Z
d 3v vk fˆa1 (x) =
d 3v vk Ha (x − ρ)
Z
Z
d 3v v⊥ fˆa1 (x) =
d 3v v⊥ Ha (x − ρ)

General Atomics Report GA-A26818

(3.50)
(3.51)
(3.52)

26

3.5

Transport Fluxes and Heating

For each species separately, we define a particle flux, a toroidal angular momentum flux, an energy
flux, and an exchange power density:
Z
Γa (r) = F d 3v Ha∗ (R) δva · ∇r ,
(3.53)
Z
1
Qa (r) = F d 3v Ha∗ (R) δva · ∇r ma v 2 ,
(3.54)
2
Z
Πa (r) = F d 3v Ha∗ (R)




c
1
[ma R(V0 + v) · eϕ ] b × ∇ δφ(x) − (V0 + v) · δA(x) · ∇r
(3.55)
B
c
R



Z
∂
∂
1
3
∗
+ V0 (x) ·
δφ(x) − (V0 + v) · δA(x)
Sa (r) = F d v Ha (R) ea
(3.56)
∂t
∂x
c
R

3.5.1

Ambipolarity and Exchange Symmetries

By summing the particle fluxes over species and using the Maxwell equations, one can prove the
exact ambipolarity property
X
ea Γ̄a = 0 ,
(3.57)
a

where an overbar denotes a perpendicular spatial and time average taken in the flux-tube limit.
Similarly, summing the exchange power density over species and using the Maxwell equations, one
can prove the net heating is zero:
X
S̄a = 0 .
(3.58)
a

Note that the time-average is only required to prove the exchange property, not the ambipolarity
property. Both these conditions are in general violated if profile variation is allowed.

3.6

Entropy production

The balance equation for entropy production is given by
Z
Z
Z
Z
∗
∗
∗
H∗
3 Ha ∂Ha
3 Ha
GL
3 Ha
σa − F d v
+F d v
Ca + F d v
Dτ + F d 3v a Dr → 0 ,
fa0 ∂t
fa0
fa0
fa0

(3.59)

where Dτ and Dr represent the (artificial) upwind dissipation terms added in the numerical discretization. The function σa is


1
3 1
1 Qa ∂ω0 Πa
1
.
σa =
−
Γa +
+
+ Sa
.
(3.60)
Lna 2 LT a
LT a Ta
∂r Ta
Ta
On taking a radial average, and the time-average over sufficiently long times, the sum of terms
should approach zero.

General Atomics Report GA-A26818

27

3.7
3.7.1

Simplified fluxes and field equations with operator notation
Operator notation

It is useful at this point to discuss the general spectral representation of fields and operators. First,
expanding an arbitrary field in a spectral (Fourier) basis gives
. X iS(R)
e
z̃(k⊥ ) ,
(3.61)
z(R) =
k⊥

where k⊥ = ∇⊥ S. In this case, the gyrophase dependence in Fourier space is harmonic
X
eiS(R) eik⊥ ·ρ z̃(k⊥ ) ,
z(R + ρ) =

(3.62)

k⊥

and the gyroaverage becomes
hz(R + ρ)iR =

X
k⊥

eiS(R) J0 (k⊥ ρa ) z̃(k⊥ ) ,

(3.63)

.
where ρa = v⊥ /Ωca . Thus, in real space, the gyroaverage can be represented as a linear operator
G0a whose spectral representation is J0 (k⊥ ρa ).
Now, we can write the field defined in Eq. 3.24 using operator notation. Moreover, we also
define a new field which is useful for calculation of momentum transport coefficients. These are
h
i
vk
.
Ψa (R) = G0a δφ(R) − δAk (R) +
c
h
i
v
k
.
Xa (R) = G2a δφ(R) − δAk (R) +
c

2
v⊥
G1a δBk (R) ,
Ωca c
2
v⊥
G3a δBk (R) .
Ωca c

(3.64)
(3.65)

Although we will construct explicit discrete approximations to the operators G0a , G1a and G2a in
the next chapter, it can be shown that they have the following spectral representations:
G0a → J0 (γa ) ,
1
G1a → [J0 (γa ) + J2 (γa )] ,
2
kx ρa
G2a → − i
[J0 (γa ) + J2 (γa )] ,
2
kx ρa
G3a → i 2 [J0 (γa ) − J1 (γa )/γa ] ,
γa

(3.66)
(3.67)
(3.68)
(3.69)

.
where γa = k⊥ ρa , and kx is defined explicitly in Sec. 5.3. The expressions above are adapted
directly from Sugama [SH98]. However, to make use of Sugama’s results, we use the following
identities, where ϕ = eϕ :
k⊥ k⊥ : (Rϕ̂)(∇ψ) = [k⊥ · (Rϕ̂)] [k⊥ · ∇ψ] ,


∂
q
∂
k⊥ · ∇ψ = − iRBp |∇r|
− Gq Θ
= RBp kx ,
∂r r
∂α
∂
k⊥ · (Rϕ̂) = − i
,
∂α

(3.70)
(3.71)
(3.72)

where Eq. (2.102) has been used to expand ∇⊥ .
General Atomics Report GA-A26818

28

3.7.2

Maxwell equations: Ha -form

The Maxwell equations are simplest when written in terms of Ha :
Poisson equation

Z
X  −na ea
1 2
3
δφ + d v G0a Ha
−
∇ δφ =
ea
4π ⊥
Ta
a

(3.73)

Parallel Ampère’s Law
X
1 2
−
∇⊥ δAk =
ea
4π
a

Z

d 3v

vk
G0a Ha
c

(3.74)

Perpendicular Ampère’s Law
−

3.7.3

X
1
ea
δBk =
4π
a

Z

d 3v

2
v⊥
G1a Ha
Ωca c

(3.75)

Maxwell equations: ha -form

For time-integration purposes, we will need to write the field equations in terms of ha :
Poisson equation
X e2 Z
1 2
2
na a
− ∇⊥ δφ +
d 3v FM a (1 − G0a
)δφ
4π
T
a
a
X Z
X e2 Z
v2
d 3v FM a G0a G1a ⊥ δBk =
ea d 3v G0a ha
−
na a
T
Ω
c
a
ca
a
a

(3.76)

Parallel Ampère’s Law
−

X e2
1 2
na a
∇⊥ δAk +
4π
Ta
a

Z

d 3v

vk2

2
FM a G0a
δAk =
c2

X
a

Z
ea

d 3v

vk
G0a ha
c

Perpendicular Ampère’s Law
 2
2
X e2 Z
v⊥
1
a
3
δB +
d v FM a
G1a δBk
na
4π k
Ta
Ωca c
a
2
X Z
X e2 Z
v⊥
v2
a
3
d v FM a
G1a G0a δφ = −
ea d 3v G1a ⊥ ha
+
na
Ta
Ωca c
Ωca c
a
a

3.7.4

(3.77)

(3.78)

Transport coefficients

Some algebra yields the simplifications:
Z
c
∂Ψa
Γa (r) = 0 F d 3v Ha∗ (R)
,
ψ
∂α
Z
c
1
∂Ψa
Qa (r) = 0 F d 3v Ha∗ (R) ma v 2
,
ψ
2
∂α



Z
Bp ∂Xa
c
Bt ∂Ψa
Πa (r) = 0 F d 3v Ha∗ (R)ma R V0 + vk
+ v⊥
,
ψ
B
∂α
B ∂α


Z
c
∂
∂
3
∗
Sa (r) = 0 F d v Ha (R) ea
+ ω0
Ψa .
ψ
∂t
∂α
General Atomics Report GA-A26818

(3.79)
(3.80)
(3.81)
(3.82)

29

Chapter 4

Normalization of Fields and Equations
4.1

Dimensionless fields and profiles

For consistency, we will use an overbar to denote reference quantities; that is, quantities which are
evaluated at the reference radius, r̄. Explicitly,

T̄e = Te (r̄)

(4.1)

n̄e = ne (r̄)
s
T̄e
c̄s =
mi

(4.2)
(4.3)

Here, mi is the mass of the main ion species (in practice, this will often be deuterium). Next, we
introduce the normalized fields
.
ĥa =
.
δ φ̂ =
.
δ Âk =
.
δ B̂k =

ha
n̄e FM a (r)
eδφ
T̄e
c̄s eδAk
,
c T̄e
δBk
,
Bunit (r)

(4.4)
(4.5)
(4.6)
(4.7)

and the normalized profiles
. na (r)
n̂a (r) =
,
n̄e
. Ta (r)
T̂a (r) =
,
T̄e
. a
ω̂0 (r) = ω0 (r) ,
c̄s
. a
γ̂E (r) = γE (r) ,
c̄s
. a
γ̂p (r) = γp (r) .
c̄s
30

(4.8)
(4.9)
(4.10)
(4.11)
(4.12)

We also have the additional normalized quantities
. B(r, θ)
B̂ =
Bunit (r)
. vk
v̂ka =
c̄s

(4.13)
(4.14)

For quantities which depend on the magnetic field strength, it is necessary to define unit quantities:
c̄s
ρs,unit (r) =
,
(4.15)
eBunit (r)/(mi c)
8πne Te
.
(4.16)
βe,unit (r) =
2
Bunit
At the reference radius, these are written as ρ̄s,unit and β̄e,unit . To measure the radial variation of
Bunit , we introduce the parameter
Bunit (r)
Gr (r) =
.
(4.17)
Bunit (r̄)

4.2
4.2.1

Velocity space normalization
Velocity variables

We also use the normalized velocity-space coordinates (ε, λ, ς), defined as
ma v 2
,
2Ta
v2
λ= ⊥ ,
v 2 B̂
ς = sgn(v̂ka ) .
ε=

(4.18)
(4.19)
(4.20)

Let us also note the identities


vk2 = v 2 1 − λB̂ ,

(4.21)

2
v⊥
= v 2 λB̂ = 2µB .

(4.22)

The pair (ε, λ) are unperturbed constants of motion. The sign of the parallel velocity, ς, is required
to separate two populations of trapped particles for each value of λ. With these definitions, the
normalized parallel velocity becomes
r
q
mi
v̂ka = ±
2εT̂a (1 − λB̂) .
(4.23)
ma

4.2.2

Dimensionless velocity-space integration

At this point, we must introduce the dimensionless velocity-space integration operator V [·]
Z ∞
Z 1
d(λB̂)
. X 1
−ε √
√
p
V [z] =
dε e
ε
z(R, λ, ε, ς) ,
(4.24)
2 π 0
0
1 − λB̂
ς=±1
where ς = sgn (v̂ka ) = ±1. It can be verified that V [1] = 1. In writing Eq. (4.24), we explicitly rule
out consideration of non-Maxwellian particle distributions.
General Atomics Report GA-A26818

31

4.3

Dimensionless equations

4.3.1

Normalized gyrokinetic equation

The normalized gyrokinetic equation is
v̂ka
∂ ĥa q ρ̄s,unit
∂ ĥa
∂ Ĥa vd ˆ
+
· ∇Ĥa + ω̂0
a[ĥa , Ψ̂a ]r,α
+
+
G
q(R
/a)
∂θ
c̄
∂α
rGr
∂ t̂
0
s
θ


h i
ma v̂ka Bt R
q ρ̄s,unit a
∂ Ψ̂a
a
+ (ε − 3/2)
+
γ̂p
− n̂a
= ĈaGL Ĥa ,
rGr
Lna
LT a
∂α
mi T̂a BR0

(4.25)

where
Ĥa = ĥa + za αa Ψ̂a ,

 2ελT̂
a
G1a δB̂k ,
Ψ̂a = G0a δφ̂ − v̂ka δÂk +
za
αa = n̂a /T̂a ,

(4.26)
(4.27)
(4.28)

and ea = eza . The inverse gradient scale lengths are defined as
1
1 ∂na
= −
,
Lna
na ∂r
1
1 ∂Ta
= −
.
LT a
Ta ∂r

4.3.2

(4.29)
(4.30)

Normalized Maxwell equations

Poisson equation:
h
h
i X
X
X
 i
2
− λ̄2D ∇2⊥ δφ̂ +
αa za2 V 1 − G0a
δφ̂ − 2
za n̂a V G0a G1a ελδB̂k =
za V [G0a ĥa ] .
a

a

(4.31)

a

Above, λ̄D is the Debye length at the reference radius

λ̄D =

T̄e
4πn̄e e2

1/2
.

(4.32)

Parallel Ampère’s Law
−

2ρ̄2s,unit
β̄e,unit

∇2⊥ δÂk +

X
a

2 2
αa za2 V [v̂ka
G0a δÂk ] =

X
a

za V [v̂ka G0a ĥa ] .

(4.33)

We remind the reader that the Ampère cancellation problem [CW03] will occur if one attempts to
2 ] = 1 rather than evaluate it numerically.
set V [v̂ka
Perperdicular Ampère’s Law
G2r

h
i X
h
i
h
i
X
X
δB̂k
2 2 2
+2
n̂a T̂a V G1a
ε λ δB̂k +
za n̂a V G1a G0a ελδφ̂ = −
T̂a V G1a ελĥa . (4.34)
β̄e,unit
a
a
a

General Atomics Report GA-A26818

32

4.3.3

Normalized Transport Fluxes

The normalized particle flux is
"
#
q ρ̄s,unit
Γa
∗ ∂ Ψ̂a
=
FV Ĥa
.
Γ̂a (r) =
n̄e c̄s
r Gr
∂α

(4.35)

The normalized energy flux is
"
#
q ρ̄s,unit
Qa
∂
Ψ̂
a
= T̂a
Q̂a (r) =
FV Ĥa∗
ε .
r Gr
∂α
n̄e T̄e c̄s
The normalized toroidal momentum flux is
Πa
Π̂a (r) =
,
n̄e mi c̄2s a
"
(
)#

B
V
B
∂
Ψ̂
∂
X̂
ma q ρ̄s,unit
R
p
0
t
a
a
FV Ĥa∗
+ v̂ka
+ v̂⊥
=
mi r Gr
a
c̄s
B
∂α
B ∂α

(4.36)

(4.37)
(4.38)

Finally, the normalized anomalous energy exchange is
 


Sa
1 q ρ̄s,unit
∂
∂
∗
Ŝa (r) =
= za
+ ω̂0
Ψ̂a .
FV Ĥa
Gr r G r
∂α
n̄e T̄e c̄s /a
∂ t̂

(4.39)

Above, FV represents the flux-surface average of the dimensionless velocity-space integration operator. The discrete respresentation of the product FV will be described in detail in the next
chapter.

4.3.4

Diffusivities

In terms of the fluxes, we further define a particle diffusivity Da according to
Γa = −Da

∂na
,
∂r

(4.40)

and an energy diffusivity χa according to
Qa = −na χa

4.3.5

∂Ta
.
∂r

(4.41)

GyroBohm normalization

In GYRO, the output fluxes and diffusivites also carry the so-called gyroBohm normalization. That
is, for output, we use

χa
χGB

Γa
ΓGB
Πa
ΠGB
Qa
QGB
Sa
SGB
Da
,
χGB

where

.
ΓGB = n̄e c̄s (ρ̄s,unit /a)2 ,

(4.42)

where

.
ΠGB = n̄e aT̄e (ρ̄s,unit /a)2 ,

(4.43)

where

.
QGB = n̄e c̄s T̄e (ρ̄s,unit /a)2 ,

(4.44)

where

.
SGB = n̄e (c̄s /a)T̄e (ρ̄s,unit /a)2 ,

(4.45)

where

.
χGB , = ρ̄2s,unit c̄s /a .

(4.46)

General Atomics Report GA-A26818

33

Chapter 5

Spatial Discretization
5.1

Foreword

The original explicit version of GYRO has been partially documented in a previous article [CW03].
This report supercedes that document in all respects. Among many other changes and improvements over [CW03], the current version of GYRO includes
1. the option to treat the fast electron parallel motion implicitly,
2. an improved and simplified treatment of boundary conditions,
3. a fully Arakawa-like nonlinear discretization scheme,
4. a generalization to an arbitrary number of kinetic impurities.
For an exhaustive, chronological list of changes, always refer to the CHANGES file with each
GYRO release.

5.2
5.2.1

Spectral Decomposition in Toroidal Direction
Expansion of fields

We expand the perturbed quantities (δφ̂, δÂk , δB̂k , ĥa ) as Fourier series in α. For example, the
potential is written as
δφ̂(r, θ, α) =

NX
n −1

δφn (r, θ) e−inα einω0 t

where

n = j∆n .

(5.1)

j=−Nn +1

In GYRO,
Nn → TOROIDAL GRID

∆n → TOROIDAL SEP

The hat is omitted on n-space quantities for brevity. Here, ω0 is a suitably-averaged rotation
frequency. In GYRO, this is taken to be the rotation frequency at the domain center. The θperiodicity condition (see condition 2, Sec. 2.2.2) requires that
δφ̂ (r, 0, ϕ + ν[ψ, 0]) = δφ̂ (r, 2π, ϕ + ν[ψ, 2π]) .
34

(5.2)

So, although the physical field, δφ̂, is 2π-periodic in θ, the Fourier representation has the implication
that the coefficients, δφn , are nonperiodic, and satisfy the phase condition
δφn (r, 0) = e2πinq(r) δφn (r, 2π) .

(5.3)

Since δφ̂ is real, the Fourier coefficients satisfy the relation δφ∗n = δφ−n . The spectral form given in
Eq. (5.1) is
1. (2π/∆n)-periodic in α at fixed (r, θ)
2. (2π/∆n)-periodic in ϕ at fixed (r, θ).

5.2.2

Poloidal wavenumber

We choose to define the poloidal wavenumber so that it is a proper flux-surface function:
kθ =

5.3
5.3.1

nq(r)
.
r

(5.4)

Operator Discretization Methods
Finite-difference operators for derivatives

The differential band width in the radial direction is denoted by the parameter id . First and second
.
derivatives can be discretized using nd -point centered differences, where nd = 2id + 1. These are
0

0

D1ii (nd , ∆x)f i =
0

0

D2ii (nd , ∆x)f i =

id
1 X
c1ν f i+ν ,
∆x

(5.5)

id
X
1
c2ν f i+ν ,
(∆x)2 ν=−i

(5.6)

ν=−id

d

where
c1ν

=

X
p6=ν

c2ν

=

X
p6=ν

Y (−j)
1
,
ν−p
ν−j

(5.7)

Y (−j)
1 X 1
.
ν−p
ν−q
ν−j

(5.8)

j6=ν,p

q6=ν,p

j6=ν,q,p

The argument nd in the operators refers to the number of points in the stencil,
 not the
 order or
accuracy of the stencil. The formal truncation error for both D1 and D2 is O (∆x)nd −1 ; in other
words, these stencils are said to be order-(nd − 1) accuracte. A typical case would be nd = 5; that
is, 5-point, 4th-order

5.3.2

Upwind schemes

To construct an arbitrary-order upwind scheme, we begin by writing the centered (nd − 1)th derivative as


id
X
1
ii0
i0
ν nd − 1
D∗ (nd , ∆x)f = −
(−1)
f i+ν .
(5.9)
ν + id
(∆x)nd −1
ν=−id

General Atomics Report GA-A26818

35

The (nd − 2)th-order upwind discretization (we will lose one order in accuracy because of the added
dissipation) of an advective derivative is then written as
v

∂
0
0
→ v D1ii (nd , ∆x) − γ|v|D∗ii (nd , ∆x)
∂x

where

.
γ = |c1id | (∆r)nd −2 .

(5.10)

The choice above for the dissipation, γ, recovers the usual first, third and higher-order upwind
schemes. For a more complete discussion of the discretization given in Eq. (5.10), see [CW03]. So,
we define a smoothing stencil
0
0
.
S ii (nd , ∆x) = |c1id | (∆r)nd −2 D∗ii (nd , ∆x) .
(5.11)
In GYRO, we add an adjustable parameter c to the upwind scheme:
v

∂
0
0
→ v D1ii (nd , ∆x) − c|v|S ii (nd , ∆x) ,
∂x

(5.12)

such that c = 1 gives the standard 1st-order and 3rd-order upwind schemes in the case nd = 3 and
nd = 5, respectively.

5.3.3

Banded pseudospectral gyro-orbit integral operators

In this section, we derive explicit forms for the gyroaverage and associated operators which were
introduced in Sec. 3.7.1. The gyroaverage of a function z(x) is defined as
Z 2π
dξ
G0a z(R) =
z [R + ρa (ξ)] .
(5.13)
2π
0
To perform the loop integral, we write the velocity and gyrovector as
v⊥
ρa =
(ex cos ξ + ey sin ξ)
Ωca
v⊥ = v⊥ (ex sin ξ − ey cos ξ)

(5.14)
(5.15)

where ex = ∇r/|∇r| and ey = b × ex . Then, we write z in spectral form
nr /2−1

z(R) =

X

X

znp (θ)e2πipr/L e−inα ,

(5.16)

znp (θ)e2πipr/L e−inα e2πipρa ·∇r/L e−inρa ·∇α .

(5.17)

n p=−nr /2
nr /2−1

z(R + ρa ) =

X

X

n p=−nr /2

We have neglected the θ-variation of the integrand, consistent with the gyrokinetic ordering (i.e.,
low parallel wavenumber). Some algebra shows
v⊥
|∇r| cos ξ ,
Ωca
v⊥
ρa · ∇α =
(ex · ∇α cos ξ + ey · ∇α sin ξ) .
Ωca
ρa · ∇r =

(5.18)
(5.19)

In terms of local equilibrium functions, we have
q
ex · ∇α = − Gq Θ ,
r
q
ey · ∇α = − Gq .
r
General Atomics Report GA-A26818

(5.20)
(5.21)

36

So, if we define
|∇r| nq
kx = 2πp
+
Gq Θ ,
L
r
nq
ky =
Gq
r

(5.22)
(5.23)

with ρa = v⊥ /Ωca , then the gyroaveraged potential for a single harmonic becomes
nr /2−1

G0a,n zn =

X

znp (θ)e

2πipr/L −inα

Z

e

0

p=−nr /2

2π

dξ i(kx ρa cos ξ+ky ρa sin ξ)
e
,
2π

(5.24)

nr /2−1

=

X
p=−nr /2

znp (θ)e2πipr/L e−inα J0 (k⊥ ρa ) ,

(5.25)

q
where k⊥ = kx2 + ky2 . To evaluate the gyroaverage explicitly, we assume that z is known on a
uniform mesh rj = j∆r:
J−1
X
j
znp e2πiprj /L ,
(5.26)
(zn ) =
p=−J

where L is the radial domain size. The radial domain and boundary conditions are described in
more detail in Sec. 7.1. The Fourier decomposition is easily inverted to yield
znp

1
=
nr

so that
0

0

jj
znj =
G0a,n

nr /2−1

X

(zn )j e−2πiprj /L ,

(5.27)

j=−nr /2

J−1
X
p=−J

e2πiprj /L J0 (k⊥ ρa )r=rj znp ,

where J0 is a Bessel function of the first kind, and
q
k⊥ = (2πp|∇r|/L + kθ Gq Θ)2 + (kθ Gq )2 .

(5.28)

(5.29)

The result gives the pseudospectral forms of the gyroaverage operator:
0

jj
G0a,n
=

1
nr

nr /2−1
0

X
p=−nr /2

wpj−j J0 (k⊥ ρa )r=rj ,

(5.30)

where indices {j, j 0 } run from −nr /2 to nr /2 − 1, and
.
wp = exp(2πip/nr ) .

(5.31)

In terms of normalized quantities,
r
ρa = ρ̄s,unit

General Atomics Report GA-A26818

p
ma 2εT̂a λB̂
mi za B̂Gr

(5.32)

37

Additional required operators are
0
(G 2 )jj
0a,n

jj 0
G1a,n
0

1
=
nr
1
=
nr

jj
G2a,n
=

5.3.4

1
nr

nr /2−1

X
p=−nr /2

0

wpj−j J02 (k⊥ ρa ) ,

nr /2−1

X

wpj−j

0

wpj−j

0

p=−nr /2
nr /2−1

X
p=−nr /2

(5.33)

1
[J0 (k⊥ ρa ) + J2 (k⊥ ρa )] ,
2

(5.34)

i
kx ρa [J0 (k⊥ ρa ) + J2 (k⊥ ρa )] ,
2

(5.35)

Banded Approximations
0

jj
The matrix G0a,n
is diagonally-dominant, such that elements {j, j 0 } for which |j − j 0 | > jg , where
jg is some sufficiently large integer, can be neglected. Thus, to convert the operators in Eqs. (5.30)
to banded form, it is enough to set
( 0
jj
0
if |j − j 0 | ≤ jg
G0a,n
jj
b
(5.36)
G0a,n =
0
if |j − j 0 | > jg

For n = 0, an additional correction is required. One must ensure that the long-wavelength limit is
asymptotically correct:
G0a,0 · 1 = 1 .
(5.37)
This is accomplished by making a small correction to the diagonal term
j+jg
jj
jj
Gb0a,0
→ Gb0a,0
+1−

0

X
j 0 =j−j

g

jj
G0a,0

(5.38)

The size of this correction (the sum above) decreases rapidly as ig is increased. We have observed
excellent results using these banded approximations for the various averaging operators, even when
the radial domain is nonperiodic. Physically, the validity of this method relies on the observation
that the numerical contribution to the averages decays rapidly at distances beyond a few gyroradii
from the gyrocenter. For linear benchmarks, when the radial domain is only one period of the
ballooning mode, we generally use full pseudospectral representations.

5.4

Discretization of the gyrokinetic equation

Upon spectral decomposition, the gyrokinetic equation finally takes the form which is solved in
GYRO. We write this symbolically as
∂ha,n
q ρ̄s,unit
− ic
ωθ Ha,n − ic
ωd Ha,n − iωc
ω∗ Ψa,n +
a{ĥa , Ψ̂a } = ĈaGL [Ha,n ] ,
E ha,n − ic
rGr
∂ t̂
where the various toroidal harmonics are given by

(5.39)

Ha,n = ha,n + za αa Ψa,n ,

(5.40)

 2ελT̂a
Ψa,n = G0a,n δφn − v̂ka δAkn +
G1a,n δBkn .
za

(5.41)

The Poisson bracket {·, ·} is defined in Eq. (5.60) of Sec. 5.4.5, the geometry functions Gθ , Gq , gsin,
gcos and Θ are defined in Sec. 2.3.10, and the profile function Gr is defined in Eq. (4.17).
General Atomics Report GA-A26818

38

5.4.1

Parallel motion on an orbit-time grid

The parallel advection term is written
− ic
ωθ Ha,n =

v̂ka
Ha,n
.
Gθ q(R0 /a) ∂θ

(5.42)

The presence of field quantities in the definition of Ha,n complicates the matter of poloidal discretization as we will discuss shortly. Equation (5.42) is subject to short-wavelength instability in
regions where the variation of v̂ka (θ) is sufficiently strong. This property is well-known [Dur99],
and time-explicit schemes must normally include dissipative smoothing if a solution is sought on an
equally-spaced θ-grid. Moreover, at bounce points θb , where vk (θb ) = 0, the distribution function
may develop cusps, bringing into question the accuracy of any finite-difference scheme on such a
grid.
This leads us to the observation that the poloidal angle, θ, is a potentially poor variable for
numerical solution of the GK equation. The obvious solution is to remove these cusps analytically
using the normalized orbit time for discretization in the poloidal direction. To this end, define
Z θ
G (θ) dθ0
1



q θ
if λ ≤
(trapped) ,


 −θb 1 − λB̂(θ0 )
B̂(π)
.
τ0 (λ, θ) = Z θ
(5.43)

Gθ (θ) dθ0
1


q
if λ >
(passing) ,


 −π 1 − λB̂(θ0 )
B̂(π)
where θb is the solution of B̂(θb ) = 1/λ. τ0 (λ, θ) must be computed numerically for general plasma
equilibria — a tedious but straightforward exercise in numerical analysis. Note that in Eq. (5.43),
and in the rest of this section, we suppress the radial dependence of τ0 (λ, r, θ), B̂(r, θ) and Gθ (r, θ)
for brevity.
Next, we introduce a normalized orbit time, τ , which runs from 0 to 2 for a given λ and
describes both signs of velocity. In this way τ parameterizes the solution on both Riemann sheets
by subsuming the two signs of velocity:
(
for 0 ≤ τ ≤ 1 (ς = 1) ,
. τ0 (λ, θ)/τ̄
τ (λ, θ) =
(5.44)
2 − τ0 (λ, θ)/τ̄ for 1 < τ ≤ 2 (ς = −1) ,
where τ̄ = τ0 (λ, θb ) for trapped particles, and τ̄ = τ0 (λ, π) for passing particles. The parallel
advection operator is reduced to one with constant velocity
p
r
v̂ka
Ha,n
∂Ha,n
mi
1
2T̂a ε
.
= Ωa (ε, λ)
, where Ωa (ε, λ) =
.
(5.45)
Gθ q(R0 /a) ∂θ
∂τ
ma q(R0 /a) τ̄
High-accuracy numerical discretization schemes for this form of the equation are well-documented.
Boundary conditions on physical distributions can be stated very simply:
• trapped are periodic on the interval [0, 2);
• co-passing are periodic on [0, 1);
• counter-passing are periodic on [1, 2).

General Atomics Report GA-A26818

39

However, the Fourier representation of Eq. (5.3)] requires that the functions h describing the passing
population are not periodic, but subject to phase conditions: h(1) = P h(0) for co-passing, and
h(2) = P h(1) for counter-passing. The important result is that as a function of τ , the trapped
distribution will be not only continuous, but also smooth, across bounce points no matter what
order difference scheme is used. The obvious physical interpretation of the location of orbit-time
gridpoints is that they are equally spaced in time, not space, along an orbit. In particular, because
they are highly stagnant near θ = ±π, particles close to the trapped-passing boundary benefit from
equal-time spacing.
For finite-n modes, we write the continuous derivative as the sum of a centered 4th-order 1st
derivative plus a 4th derivative smoother:


0
∂Ha,n
0
0
j0
Ωa
= Ωa D1jj (5, ∆τ )Ha,n
− c|Ωa |S jj (5, ∆τ )hja,n .
(5.46)
∂τ
j
.
The grid spacing for the parallel motion is taken to be ∆τ = 2/nτ , and c is a parameter which
measures the amount of numerical disipation. When c = 1 (and Ψa,n = 0) we are left with the
usual third-order upwind method. We have shown previously that in this scheme finite dissipation
is required to give proper Landau damping of n = 0 GAMs [CW03]. Conversely, dissipation-free
schemes suffer from recurrence problems. Note that the dissipation, above, acts on ha,n but not
Ha,n . In fact, attempting to add upwind diffusion to Ha,n (that is, also to fields) will generate
numerical modes with growth rate proportional to the square of the radial box size. The difference
scheme given in Eq. (5.46) is applied only if the species is to be integrated in time explicitly. When
implicit integration is used, we set c = 0.

5.4.2

Er shear

Since we have already shifted to a frame rotating with the mean frequency ω¯0 , the remaining
shearing term is
− iωc
(5.47)
E ha,n = −in [ω0 (r) − ω0 ] ha,n .
In a simulation with fixed shearing rate, this can be approximated as
− iωc
E ha,n → −in

5.4.3

∂ω0
(r − r̄)ha,n = ikθ γE (r − r̄)ha,n .
∂r

(5.48)

Drift motion

Using the geometry formulae presented in Sec. 2.3.11, we write the drift operator as
−ic
ωd Ha,n = − iGq kθ (Rk + R⊥ ) (gcos1 + gcos2 + Θgsin) Ha,n
+ iGq kθ Rk gcos2 Ha,n

− iGq kθ Rc (ucos + Θusin) Ha,n
∂Ha,n
− |∇r|(Rk + R⊥ )gsin
∂r
∂Ha,n
− |∇r|Rc usin
,
∂r

General Atomics Report GA-A26818

(5.49)
(5.50)
(5.51)
(5.52)
(5.53)

40

where the dimensionless quantities Rk /a, R⊥ /a and Rc /a are defined as

vk2
Rk .
ρ̄s,unit 2εT̂a 
=
=
1 − λB̂ ,
a
c̄a Ωca R0
R0 za Gr B̂
!
ρ̄s,unit 2εT̂a
λB̂
R⊥ .
µB
=
=
,
a
c̄a Ωca R0
R0 za Gr B̂
2
ρ̄s,unit ma 2v̂ka (ω0 R0 /c̄s )
Rc . 2vk ω0 R0
.
=
=
a
c̄a Ωca R0
R0 mi
za Gr B̂
The derivative operator is discretized according to:


∂Ha,n
0
0
0
i0
gsin a
= gsin D1ii (nd , ∆r/a)Ha,n
− c|gsin|S ii (nd , ∆r/a)hia,n ,
∂r
i

(5.54)
(5.55)
(5.56)

(5.57)

where
nd = 2 × RADIAL DERIVATIVE BAND + 1

5.4.4

Diamagnetic effects

The spectral form of the diamagnetic term is simply


ma v̂ka Bt R
kθ ρ̄s,unit a
a
+
− ic
ω∗ Ψa,n = in̂a
+ (ε − 3/2)
γ̂p Ψa,n .
Gr
Lna
LT a
mi T̂a BR0

5.4.5

(5.58)

(5.59)

Poisson bracket nonlinearity

Numerical discretizations of the nonlinear E × B motion (including electrstatic, flutter and compressional components) are subject to short-wavelength instabilities. A numerical scheme which
preserves continuous conservation laws for domain-integrated number, energy and enstrophy (and
thereby prevents both instability and a cascade to short wavelengths) was proposed by Arakawa in
1966 [Ara66]. We adapt this scheme to a semi-spectral variant suitable for use with GYRO. First,
let us begin by writing the continuous form of the bracket appearing in the gyrokinetic equation:
{F, G} =

∂F ∂G ∂G ∂F
−
.
∂α ∂r
∂α ∂r

Eq. (5.60) can be recast into the following alternative but equivalent forms




∂
∂G
∂
∂G
{F, G} =
F
−
F
,
∂α
∂r
∂r
∂α
and
{F, G} = −

∂
∂α





∂F
∂
∂F
G
+
G
.
∂r
∂r
∂α

The trick, first proposed by Arakawa, is to sum these and divide by 3 to yield






1 ∂
∂G
∂F
∂
∂F
∂G
∂F ∂G ∂G ∂F
{F, G} =
F
−G
+
G
−F
+
−
.
3 ∂α
∂r
∂r
∂r
∂α
∂α
∂α ∂r
∂α ∂r

General Atomics Report GA-A26818

(5.60)

(5.61)

(5.62)

(5.63)

41

Using the spectral decomposition of the bracket together with the discrete form of the derivative
ii0 (∆r/a), gives the final discrete form
operator, DN
{F, G}in =



1X
0
0
0
i0
(n + n0 ) Fni 0 D1ii Gin−n0 − Gin0 D1ii Fn−n
0
3 0
n

1 X 0 ii0  i0
0
i0
i0
+
n D1 Gn−n0 Fni 0 − Fn−n
. (5.64)
0 Gn0
3 0
n

This is the expression used in GYRO. Now, let us assume, in what follows, that the radial domain
is periodic. Deviations from the results below will as a consequence be solely limited to boundary
effects.
Lowest order invariant
The lowest integral invariant, which measures the rate of change of the distribution along the
nonlinear flow, is
Z
Z
X
dα dr {F, G} =
{F, G}i0 = 0 .
(5.65)
i

P ii0 i0
Using Eq. (5.64), the sum above vanishes for all DN which satisfy i DN
a = 0 for all vectors a.
Centered-difference formulae of all orders for DN will satsify this condition.
Higher-order invariants
Some algebra shows that the additional quantities vanish:
Z
Z
XX
dα dr G {F, G} =
Gin {F, G}i−n = 0 ,
n

Z

Z
dα

5.5

dr F {F, G} =

XX
n

(5.66)

i

i

Fni {F, G}i−n = 0 .

(5.67)

Blending-function expansion for fields

Since the distribution function ha,n is computed at a different set of points {θj } for each discrete
value of λ, there is no natural way to solve the Maxwell equations using finite-difference methods
on a fixed poloidal grid. Instead, we adopt a function-space approach, and expand the fields
(δφn , δAkn , δBkn ) in series of uniform polynomial blending functions. We will show that basis
functions which incorporate the complex phase conditions [see Eq. (5.3)] at −π and π can be
constructed from pairs of these blending functions. Equations for the expansion coefficients are
then obtained using the well-known Galerkin method.
In computing the poloidal dependence of the fields, there are a number of separate discretization
effects to consider. Since the fields are sums (integrals in the continuum limit) of distribution
functions, low velocity-space resolution will lead to poor poloidal accuracy — even if the there are
no other sources of discretization error. Conversely, even when a very large number of velocity-space
grids are used, discretization error in the orbit-time integration of ha,n will lead to poor poloidal
accuracy — even if the function-space method of the present section was exact. Thus, achieving an
absolute level of accuracy in the θ-dependence of the fields requires sufficient convergence in both
velocity-space and in orbit-time.
General Atomics Report GA-A26818

42

Examples of blending functions
A finite number of blending functions, Nm (s), are used to provide a basis for the field expansion.
An important feature of the Nm , however, is that they are all translates of a single function N (s),
such that N (s − m) = Nm (s). In practise, we use one of the following three types of blending
functions:
Piecewise linear:
N

(2)

.
(s) =

(
s
2−s

if 0 ≤ s ≤ 1 ,
if 1 ≤ s ≤ 2 .

(5.68)

Piecewise quadratic:

2

s /2
.
N (3) (s) = −(3/2) + 3s − s2


(3 − s)2 /2

if 0 ≤ s ≤ 1 ,
if 1 ≤ s ≤ 2 ,
if 2 ≤ s ≤ 3 .

(5.69)

Piecewise cubic:

3

(1/6)s



. (2/3) − 2s + 2s2 − (1/2)s3
N (4) (s) =

−(22/3) + 10s − 4s2 + (1/2)s3



(1/6)(4 − s)3

if
if
if
if

0≤s≤1,
1≤s≤2,
2≤s≤3,
3≤s≤4.

(5.70)

Representation of a quasi-periodic function
It remains to construct meaningful set of basis vectors from the prototype blending functions given
above. Let us begin by considering a function z(θ) which is not periodic in θ, but satisfies a phase
condition z(θ + 2π) = Pz(θ). On the infinite domain (−∞, ∞) we can write
∞
X

z(θ) =

=

cm Nm (θ/∆θ) ,

m=−∞
nX
blend
m=1

[· · · + cm−nblend Nm−nblend + cm Nm + cm+nblend Nm+nblend + · · · ] ,

(5.71)

(5.72)

.
where Nm (θ/∆θ) = N (θ/∆θ − m), ∆θ = 2π/nblend , and nblend is the number of basis functions
used to represent one 2π-segment of the function z. The phase condition in θ translates into an
equivalent phase condition for the blending coefficients themselves; namely cm+nblend = P cm . The
function z is therefore completely described by the coefficients {c1 , . . . , cnblend } and basis functions
{F1 , . . . , Fnblend } according to
z(θ) =

nX
blend
m=1

cm Fm (θ)

where

Fm (θ) = N (θ/∆θ − m) + PN (θ/∆θ − m − nblend ) .

(5.73)

In choosing a value of the parameter nblend , one is constrained by the number of points used in
the orbit time discretization, nτ . Given that a passing particle is described by nτ /2 orbit points,
choosing values of nblend larger than nτ /2 will lead to a breakdown of the method since there will
be less than one orbit point per blending function segment. For the simulations presented in this
General Atomics Report GA-A26818

43

F6
F5
F4
F3
F2
F1
-π

0
θ

π

Figure 5.1: Quadratic basis functions, Fm (θ), for nblend = 6 and phase P = −1.
paper, we have found that nblend = 6 and nτ = 20 are efficient choices. However, we expect that
the simulation of plasmas with strong equilibrium shaping, will require a simultaneous increase of
both nblend and nτ .
Figure 5.1 shows the nblend = 6 quadratic basis functions {F1 , . . . , F6 } for P = −1. We remark
that since P is generally a complex number, and is a function of radius, the basis functions are also
complex and functions of radius. Knowing this, we can write the expansion of fields in terms of
the basis functions, with explicit radial dependence indicated, as
δφn (ri , θ) =
δAkn (ri , θ) =
δBkn (ri , θ) =

nb
X
m=1
nb
X
m=1
nb
X

i
φ̃im Fm
(θ)

(5.74)

i
Ãim Fm
(θ)

(5.75)

i
i
B̃m
Fm
(θ)

(5.76)

m=1

Although the blending expansion coefficients φ̃ and Ã depend on the toroidal mode number n, we
supress this dependence for brevity. A few well-known but important observations concerning the
approximation properties of the linear, quadratic and cubic functions should be stated. In a region
of extremely rapid field variation, the cubic approximation is the least robust — suffering from
overshoot/undershoot oscillations. Conversely, as the fields become progressively smoother, the
linear approximation will give the poorest interpolation accuracy. A further undesireable property
of the linear approximation is its discontinuous first derivative. Although there is no general rule or
obviously best interpolation order for all situations, we find that the piecewise quadratic functions
are sufficiently accurate and robust for all cases studied to date.

General Atomics Report GA-A26818

44

5.6

Velocity-Space Discretization

To solve the Maxwell equations, and to compute the transport fluxes, we will need to develop a
discrete representation of the operator FV . In the present section, we discuss quadrature methods
for evaluation of these operators. Since the techniques which follow apply independently to all
species, we omit species indices for brevity.

5.6.1

Decomposition of FV

Since V [1] = 1, it follows automatically that FV [1] = 1. An equivalent statement of this result is
embodied in the identity
Z ∞
Z λ∗
I
1
−ε √
√
dε e
ε
dλ τ̄ dτ = J0 ,
(5.77)
2 π 0
0
with λ∗ = 1/B̂(r, 0) the maximum possible value of λ, and J0 defined in Sec. 2.2.4. Above, we have
used the integral identity
XZ
σ

π

−π

1/B̂(r,θ)

Z
dθ Gθ (r, θ)
0

dλ

Z

q
=
1 − λB̂(r, θ)

λ∗

I
dλ τ̄

dτ .

(5.78)

0

Since the numerical evaluation of integrals of this type will involve separate integration weights in
each of the variables ε, λ and τ , it is useful to make the decomposition
FV = Vε × Vλ × Vτ ,

(5.79)

with factors
.
Vε [h] =
.
Vλ [h] =
.
Vτ [h] =

2
√
π

Z

∞

√
dε e−ε ε h ,

(5.80)

0

Z λ∗
1
dλ τ̄ h ,
2 J0 0
I
1
dτ h .
2

(5.81)
(5.82)

These are normalized so that Vε [1] = Vλ [1] = Vτ [1] = 1. The task at hand, now, is the construction
of discrete forms of these operators, and therefore of FV .

5.6.2

Energy Integration

To develop a quadrature method for the energy integral, we split the interval of integration in Vε
— defined in Eq. (5.80) — into two regions: [0, ε∗ ) and [ε∗ , ∞), where ε∗ is the maximum energy
gridpoint (input). Integration over the first interval is done by changing variables according to
Z ε
√
. 2
x(ε) = √
dε e−ε ε .
(5.83)
π 0
.
We let x0 = x(ε∗ ), and evaluate the integral using Gauss-Legendre integration [BF85] over nε − 1
points:
Z x0
nX
ε −1
dx h[ε(x)] '
wi h[ε(xi )] .
(5.84)
0

General Atomics Report GA-A26818

i=1

45

i
1
2
3
4
5
6

nε = 4,
εi
0.2924572707
1.0404041729
2.2531263847
3.0000000000

ε∗ = 3.0
wi
0.2467749375
0.3948399000
0.2467749375
0.1116102251

nε = 6,
εi
0.1625303849
0.5442466206
1.1228428641
1.9784353093
3.2361246631
4.0000000000

ε∗ = 4.0
wi
0.1130127375
0.2283030745
0.2713566704
0.2283030745
0.1130127375
0.0460117057

Table 5.1: Sample energy abscissae and weights

The abscissae and weights (xi , wi ) are the usual Gauss-Legendre ones. Note that we must solve the
nonlinear equations xi = x(εi ) for εi to obtain the energy gridpoints. With the dominant part of
the energy integration done, we evaluate the remaining, infinite integral according to
Z ∞
Z ∞
√
2
2
−ε √
∗
√
√
dε e
dε e−ε ε = (1 − x0 )h(ε∗ ) .
ε h(ε) ' h(ε )
(5.85)
π ε∗
π ε∗
This gives the final weight 1 − x0 at the energy gridpoint ε∗ , for a total of nε gridpoints. We remark
that this method has the desireable property
nε
X

wi = 1

such that

i=1

∀i, wi > 0 .

(5.86)

Some sample abscissae and weights are given in Table 5.1 to limited precision (10 significant
digits). Note that it is straightforward to generate these, and thus to enforce the sum in Eq. (5.86),
to machine precision. The abcissae and weights are unique for a given nε and ε∗ . Outside of this
(ε)
section, to avoid ambiguity, we will use index iε and weight wiε to refer to energy integration.

5.6.3

λ Integration

The λ integration follows essentially the same strategy as the energy integration, with only minor
differences. First, introduce the integration variable
Z λ
. 1
dλ0 τ̄ (λ0 ) .
(5.87)
x(λ) =
J0 0
Then, then Vλ can be expressed as
Z
Vλ [h] =

1

dx h[x(λ)] .

(5.88)

0

Because h is rapidly-varying across the trapped-passing boundary at xt = x(λt ), where λt = 1/B(π),
it is wise to split the previous integral into two regions:
Z 1
Z xt
Z 1
dx h[x(λ)] =
dx h[x(λ)] +
dx h[x(λ)] .
(5.89)
0

0

xt

Gauss-Legendre rules are then applied to each integral separately to determine abscissae and weights
(xi , wi ). Then, as before, the equations xi = x(λi ) must be inverted numerically to obtain the
gridpoints λi .
General Atomics Report GA-A26818

46

i
1
2
3
4
5
6

xi
0.0701916516
0.3114046777
0.5526177038
0.6653193692
0.8114046777
0.9574899862

λi
0.1353809983
0.5237287429
0.7871934230
0.8581589333
0.9778736670
1.1217887702

wi
0.1730025988
0.2768041580
0.1730025988
0.1047751790
0.1676402866
0.1047751790

Table 5.2: Sample λ weights for npass = ntrap = 3.

We emphasize that there is no assumption of continuity across the trapped passing boundary.
In fact, in the collisionless limit, we do not expect the distribution to be continuous there. With
weak collisions, we expect the formation of of a boundary layer around x = xt . In the latter case,
we expect good layer resolution because the Gauss-Legendre scheme puts integration abcissae very
close to (but not on) xt . When combined with the orbit-time grid, the λ integration abcissae
provide the (θ, vk /v) gridpoint distribution shown in Fig. 8.1.
Sample values of abscissae and weights are given in Table 3 for a circular equilibrium with
(λ)
r/R0 = 1/6. Outside of this section, to avoid ambiguity, we will use index iλ and weight wiλ to
refer to pitch-angle integration.

5.6.4

Discretization Summary

Our methods for solution of the Maxwell equations have the implication that it is not V but FV for
which a discrete form is required. But we have already done enough to show that the discretization
takes the form
nλ
nε
nτ
X
X
X
(ε)
(λ)
(τ )
FV [h] →
wiε
wiλ
wiτ hiε iλ iτ ,
(5.90)
iε =1

iλ =1

iτ =1

(τ )

where, so far, no radial discretization has been employed. The τ -weights are simply wiτ =
(1/2)∆τ = 1/nτ . The numerical representation of the weights is such that when h = 1, the
sum in Eq. (5.90) is unity to machine precision.

5.7

Connection to Ballooning Modes

In what follows, we limit our discussion to s − α geometry. Recall the form of a toroidal harmonic
when expanded in a radial Fourier series
X
δφn (r, θ)e−inα = e−in[ϕ−q(r)θ]
e2πipr/L δφnp (θ)
(5.91)
p

where L, the radial domain size, is written as a multiple, M , of the natural quantization length,
1/(nq00 ), of the ballooning mode:
1
r0
.
L=M 0 =M
.
nq0
nq0 s0

General Atomics Report GA-A26818

(5.92)

47

So we can write
δφn (r, θ)e−inα = e−in[ϕ−q(r)θ]

X

−in[ϕ−q(r)θ]

X

0

e2πipr(nq0 /M ) δφnp (θ) ,

(5.93)

e2πip(nq/M ) Φ`B (θ + 2πp0 ) .

(5.94)

p

=e

p

where we have introduced the ballooning potential Φ`B via
δφnp (θ) = Φ`B (θ + 2πp0 )e2πinq0 p/M .

(5.95)

This transformation is motivated by the fact that the function δφn (r, θ) is not periodic but instead
satisfies the phase condition in Eq. (5.3). The meanings of p0 and ` are not yet specified. Writing
p = M p0 + `, and breaking the sum into parts yields
δφn (r, θ)e−inα = e−in[ϕ−q(r)θ]

M
−1 X
X
`=0

= e−inϕ

M
−1 X
X
`=0

0

`

einq(θB +2πp ) Φ`B (θ + 2πp0 ) ,

(5.96)

p0
0

`

einq(θ+θB +2πp ) Φ`B (θ + 2πp0 ) ,

(5.97)

p0

where the function Φ`B (θ) is continuous and bounded over the domain −∞ < θ < ∞. The usual
equation for Φ`B can be obtained by substituting Eq. (5.97) into the gyrokinetic equation, and then
` = 2π(`/M ) is the ballooning
solving for Φ`B for a given value of ` on the infinite θ-domain. Here, θB
angle. In GYRO, we have
M = BOX MULTIPLIER .
(5.98)
On a discrete grid, the indices have the following ranges of variation:
nr
nr
−1,
p = − ,...,
2
2
nr
nr
p0 = −
,...,
−1,
2M
2M
` = 0, . . . , M − 1 .

(5.99)
(5.100)
(5.101)

such that the discrete Fourier transform of the potential is
δφnp (θ) =

1
nr

nr /2−1

X

δφn (rj , θ)e−2πipj/nr .

(5.102)

j=−nr /2

In conclusion, over the interval
−π

n


n

r
+1 ≤θ <π
−1 ,
M
M
r

(5.103)

we have the reconstruction algorithm
0

`=0:

Φ0B (θ + 2πp0 ) = δφn,M p0 (θ)e−2πinq0 p ,

`=1:

Φ1B (θ
Φ2B (θ

`=2:

(5.104)

0

−2πinq0 p0 −2πinq0 (1/M )

,

(5.105)

0

−2πinq0 p0 −2πinq0 (2/M )

,

(5.106)

+ 2πp ) = δφn,M p0 +1 (θ)e
+ 2πp ) = δφn,M p0 +2 (θ)e

e
e

..
.
`=M −1:

(5.107)
p0

M −1
ΦB
(θ + 2πp0 ) = δφn,M p0 +M −1 (θ)e−2πinq0 e−2πinq0 (M −1)/M .

General Atomics Report GA-A26818

(5.108)

48

Chapter 6

Temporal Discretization
The gyrokinetic treatment of electrons is problematic in simulations because of the emergence of
a host of numerical instabilities connected with the discretization of the electron parallel motion.
The presence of the electrostatic and magnetic potentials in the advection term,

v̂ka
∂ ĥa
∂ 
+
ĥa + za αa Ψ̂a + · · ·
Gθ q(R0 /a) ∂θ
∂ t̂

(6.1)

for a = e give rise to stability considerations much more troublesome than the simple electron
parallel Courant limit. Physically, the parallel term gives rise to the high-frequency electrostatic
Alfvén wave [LLHL01] at β = 0. In the slab limit, one can show
r
kk mi
ωH =
Ωci .
(6.2)
kr me
This mode is physically pathological, because the frequency increases indefinitely as me decreases,
and also as kr decreases – that is, as the radial box size grows. While nonzero β provides a
cutoff, the stable numerical treatment of this term is nontrivial.1 To overcome this severe and
intrinsic simulation difficulty, we were persuaded to treat the electron advection implicitly. We
experimented with a variety of implicit-explicit (IMEX) Runge-Kutta (RK) schemes that were
popular in the period 2001-2003 [ARW95, ARS97, PR00, PR02, KC03], settling on an IMEX-SSP
scheme of Pareschi and Russo [PR02].

6.0.1

Reduction to canonical form

For brevity, we suppress the toroidal harmonic index and abbreviate ha,n as ha . With all spatial
operators discretized, the GKM system has the form
e a ({ha }, he ) , a = 1, . . . , nion
ḣa = H
e e ({ha }, he ) + He ({ĥa }, he ) ,
ḣe = H

(6.3)
(6.4)

where {ha } = h1 , . . . , hnion . In Eqs. (6.3) and (6.4), the use of a tilde on the RHS represents a
nonstiff term, while the absence of a tilde indicates a stiff term. The IMEX-RK schemes are written
for equations in the canonical form
ẏ = Ye (y) + Y (y) ,
(6.5)
1

The ability of Eulerian codes to treat this term both accurately and stably has been a major contributor to the
relative success of Eulerian solvers in comparison to their PIC counterparts.

49

with Ye the explicit RHS and Y the implicit RHS. Thus, we have the connection
#
"

 

ea}
{
H
0
{ĥa }
e
Y =
y=
Y =
.
e
He
he
He

(6.6)

Now, there is some hidden complexity in how we have written the GKM equations, as the functions
H depend on the distribution function h both directly, and indirectly through the fields. We
emphasize that it is only the fast electron advection term, He , which we wish to (or are able to)
treat implicitly:

i
a
∂ h
.
He = − v̂ke (r, θ, ε)
ĥe − αe δ φ̂ − v̂ke δ Âk − ελT̂e δ B̂k
(6.7)
qR0 Gθ ∂θ
h

i
= − Ω(r, θ, ε)∂τ ĥe − αe δ φ̂ − v̂ke δ Âk − ελT̂e δ B̂k
,
(6.8)
where from this point forward we suggest that the reader consider the object ∂τ as a matrix.

6.0.2

IMEX-RK-SSP schemes of Pareschi and Russo

The second-order IMEX schemes we consider can be summarized in Butcher tableau form as
Implicit
Explicit
0
a11
ã21
0
a21 a22
ã31 ã32 0
a31 a32 a33
w̃1 w̃2 w̃3
w1 w2 w3
This notation corresponds to the following explicit evolution equations:
y1 = y + ∆t a11 Y1
y2 = y + ∆t(ã21 Ye1 + a21 Y1 ) + ∆t a22 Y2

(6.10)

y3 = y + ∆t(ã31 Ye1 + ã32 Ye2 + a31 Y1 + a32 Y2 ) + ∆t a33 Y3

(6.11)

ȳ = y + ∆t

3
X
k=1

w̃k Yek + ∆t

3
X

(6.9)

wk Yk .

(6.12)

k=1

Above, y is the old field vector (time t) and ȳ is the new field vector (time t + ∆t). In GYRO, we
currently use the SSP2(3,2,2) scheme:
Implicit
Explicit
0
1/2
0
−1/2 1/2
SSP2(3,2,2) 0
0
1
0
0
1/2 1/2
0 1/2 1/2
0
1/2 1/2
Since the diagonal coefficients of the implicit stages are equal, we need to solve the same implicit
system at each stage. It would be good at some point to implement the SSP2(3,3,2) scheme, which
requires the solution of two different implicit problems, but has half the advective Courant limit in
explicit part.
Explicit
Implicit
0
1/4
1/2
0
0
1/4
SSP2(3,3,2)
1/2 1/2
0
1/3 1/3 1/3
1/3 1/3 1/3
1/3 1/3 1/3
General Atomics Report GA-A26818

50

6.0.3

Implementation in GYRO

The complication which arisese when applying these methods to GYRO are connected with the
field solve. For all the so-called SDIRK methods, the implicit equations at each stage have a generic
form. In the present case, at each integrator stage (k = 1, 2, 3), we are faced with:
X
e a )p ,
(ha )k = ha + ∆t
ãkp (H
(6.13)
p s0 , where s0 is the maximum linear growth rate. The inversion formula is given by the
Bromwich integral
Z c+i∞
1
. −1 e
f (t) = L f =
fe(s)est ds , 0 < t < ∞ .
(9.8)
2πi c−i∞
It will be convenient to use the variable ω = is in subsequent formulae. Now, it is then easy to
show that


1
ω − ω∗
ea f0a e
e a (ω) = LHa,n =
H
iha (0) +
Ψa ,
(9.9)
ω + ωθ + ωd + ωC
ω + ωθ + ωd + ωC
Ta
e a = LΨa,n . Upon substitution into the Laplace transform of the Maxwell equations, we
where Ψ
are left with a system of the form
0

0

Mσσ (ω)Φσ (ω) = S σ (ω)

General Atomics Report GA-A26818

(9.10)

64

where σ and σ 0 are field indices which run from 1 to 3. The field matrix and source are given by
X e2 Z
ω − ω∗
0
a
σσ 0
σ
(9.11)
M (ω) = δσσ0 L −
Gσ a ,
d 3v f0a Gσa
Ta
ω + ωθ + ωd + ωC
a
X Z
1
σ
S (ω) =
ea d 3v Gσa
iha (0) .
(9.12)
ω + ωθ + ωd + ωC
a
We have defined the additional 3-index vectors:

 . 
e δA
gk , δB
gk ,
Φ1 , Φ2 , Φ3 = δφ,
 .
L ,L ,L =
1

1a

2

2a

3

3a

G ,G ,G



X e2
1
a
− ∇2⊥ +
4π
T
a
a

(9.13)
Z

d 3v f0a ,

1 2
1
∇⊥ , −
4π
4π

!



2
vk
v⊥
= G0a , − G0a ,
G1a .
c
Ωca c

,

(9.14)
(9.15)

In what follows, we will refer to M(ω) as the Maxwell dispersion matrix. The roots of the equation
det M(ω) = 0 correspond to the normal modes of the system. When the velocity integrals are
taken along the real velocity axes, the integrals used to compute M define a function of s analytic
in region s > 0, which corresponds to the upper-half ω-plane. This means that unstable modes can
be readily computed using the same (real) velocity discretization as in the initial-value problem.
Calculation of stable normal modes, on the other hand, would require analytic continuation of
M(ω) into the lower half-plane, Im(s) ≤ 0, by deformation of the contour in velocity space, or
perhaps by numerical analytic continuation.
The numerical implementation of the eigenvalue solver, which employs the existing GYRO
spatial discretization methods, is described in detail in the Appendix. Here, we remark that the
size of the final matrix problem is small, with rank(M) = nσ nr nb , with nσ the number of fields,
nr the number of radial gridpoints, and nb the number of poloidal finite elements. For a basic
electrostatic case, this can be as small as rank(M) = 24. The dominant cost is therefore not
computing det M(ω), but rather computing the inverse P −1 , where P is the matrix representation
of the propagator P = ω + ωθ + ωd + ωC .

9.4

Discretization and Implementation in LAPACK

The spatial discretization of the differential operators required for construction of the Maxwell field
matrix follows precisely the same approach as for the GYRO initial-value solver, with the various
stencils and quadrature methods discussed in detail in Ref. [CW03]. The velocity-space variables
in GYRO are
. Bunit
. ma v 2
λ = 2 Bv 2 and ε =
(9.16)
2Ta
v⊥

r
Fluctuating quantities are evaluated on a species-independent mesh with radial nodes {ri }ni=1
,
nλ
nε
nτ
pitch-angle nodes {λk }k=1 , energy nodes {εµ }µ=1 and orbit-time nodes {τm }m=1 . The gyroaverage
e a , for species a, has the discrete representation
of the effective potential Ψ
X σaµ
e a )µ =
(Ψ
Gii0 km Φσi0 km ,
(9.17)
ikm

i0 σ

General Atomics Report GA-A26818

65

where the three-potential Φσ at the same point in configuration space is given by the complex
Galerkin representation
X
Φσikm =
cσij Fij (θkm ) .
(9.18)
j

cσij

Here,
are the so-called blending coefficients, and Fij the basis functions defined in Sections 5.2
and 5.3 of Ref. [CW03]. The propagator has the matrix form
aµ
aµ
0
0
0
0
0
Piiaµ0 kk0 mm0 = ωδii0 kk0 mm0 + (ωθ )aµ
ikmm0 δii kk + (ωd )ii0 km δmm kk + (ωC )ikk0 mm0 δii .

(9.19)

e a in terms of the inverse of the propagator as
We can write the nonadiabatic distribution H
!aµ
X σaµ
X
ea
ea na −1 aµ
H
aµ
=
(P )ii0 kk0 mm0 (ω − ω∗i
Gi0 i00 k0 m0
cσi00 j Fi00 j (θk0 m0 ) .
(9.20)
0 k 0 m0 )
fM a
Ta
σ
j

ikm

Constructing the Galerkin projections of all three Maxwell equations, using the technique described
in Section 5.3 of Ref. [CW03], yields the matrix equation
h
i 0
0
0
0
Miiσσ0 jj 0 cσi0 j 0 = Aσii0 jj 0 δσσ0 − Biiσσ0 jj 0 (ω) cσi0 j 0 = 0 ,
(9.21)
where
Aσii0 jj 0 =

X

∗
Fijkm
Lσii0 km Fi∗0 j 0 km

(9.22)

km

and
0

Biiσσ0 jj 0 (ω) =

X X X e 2 na µ
aµ
σ 0 aµ
a
∗
−1 aµ
wkm Fijkm
Gσaµ
)i000 i00 kk0 mm0 (ω − ω∗i
000 km (P
00 k 0 m0 )Gi00 i0 k 0 m0 Fi0 j 0 k 0 m0
ii
Ta
00 0 0
aµ 000
i km i k m

(9.23)
=

XX
aµ

aµ
aµ
Uqp
(P −1 )aµ
pp0 Vp0 q 0

(9.24)

p,p0
0

σσ
= Bqq
0 .

(9.25)

µ
Here, the weights wkm
are the products of the energy, pitch-angle and orbit-time weights defined
in Eq. (72) of Ref. [CW03]. In terms of these weights, the flux-surface average of the velocity-space
integration is written as
Z
X µ
X µ
F d 3vfM a Ψa →
wkm (Ψa )µikm with
wkm = 1 .
(9.26)
kmµ

kmµ

In addition, we have defined the lumped indices p = (i000 , k, m), p0 = (i00 , k 0 m0 ), q = (i, j, σ) and
q 0 = (i0 , j 0 , σ 0 ). High performance is achieved by computing the inverse P −1 using the LAPACK
routines ZGETRF (LU decomposition) followed by ZGETRI (inverse), with the subsequent matrix
triple-product U P −1 V evaluated using two sequential calls to the BLAS ZGEMM kernel. Finally,
det(M ) is computed by the formula
Y
det(M ) = ±
Lqq
(9.27)
q

where L is the lower-triangular matrix returned by the ZGETRF factorization. The upper (lower)
sign is taken if an even (odd) number of row permutations were made in the factorization.
General Atomics Report GA-A26818

66

Acknowledgements
This work was supported by the U.S. Department of Energy under DE-FG02-95ER54309.

67

Bibliography
[AL80]

T. Antonsen and B. Lane. Kinetic equations for low frequency instabilities in inhomogeneous plasmas. Phys. Fluids, 23:1205, 1980.

[Ara66]

A. Arakawa. Computational design for long term numerical integration of the equations
of fluid motion: Two-dimensional incompressible flow. J. Comput. Phys., 1:119, 1966.

[ARS97]

U.M. Ascher, S.J. Ruuth, and R.J. Spiteri. Implicit-explicit runge-kutta methods
for time-dependent partial differential equations. Applied Numerical Mathematics,
25:151–167, 1997.

[ARW95]

U.M. Ascher, S.J. Ruuth, and B.T. Wetton. Implicit-explicit methods for timedependent pdes. SIAM J. Numer. Anal., 32:797–823, 1995.

[Bas09]

E.M. Bass. Bull. Am. Phys. Soc., 54:344, 2009.

[BC08]

E.A. Belli and J. Candy. Kinetic calculation of neoclassical transport including selfconsistent electron and impurity dynamics. Plasma Phys. Control. Fusion, 50:095010,
2008.

[BF85]

R.L. Burden and J.D. Faires. Numerical Analysis. Prindle, Weber and Schmidt,
Boston, 1985.

[BHD08]

E.A. Belli, G.W. Hammett, and W. Dorland. Effects of plasma shaping on nonlinear
gyrokinetic turbulence. Phys. Plasmas, 15:092303, 2008.

[BKC+ 84]

C.M. Bishop, P. Kirby, J.W. Connor, R.J. Hastie, and J.B. Taylor. Ideal MHD ballooning stability in the vicinity of a separatrix. NF, 24:1579, 1984.

[Bur97]

K.H. Burrell. Effects of E × B velocity shear and magnetic shear on turbulence and
transport in magnetic confinement devices. Phys. Plasmas, 4:1499, 1997.

[CHT78]

J.W. Connor, R.J. Hastie, and J.B. Taylor. Shear, periodicity and plasma balloning
modes. Phys. Rev. Lett., 40:396, 1978.

[CPC+ 03]

Y. Chen, S.E. Parker, B.I. Cohen, A.M. Dimits, W.M. Nevins, D. Schumaker, V.K.
Decyk, and J.N. Leboeuf. Simulations of turbulent transport with kinetic electrons
and electromagnetic effects. Nucl. Fusion, 43:1121, 2003.

[CTB81]

P.J. Catto, W.M. Tang, and D.E. Baldwin. Ion transport in toroidally rotating tokamak plasmas. Plasma Phys., 23:639, 1981.

[CW03]

J. Candy and R.E. Waltz. An Eulerian gyrokinetic-Maxwell solver. J. Comput. Phys.,
186:545, 2003.
68

[DHCS91]

W.D. D’haeseleer, W.N.G. Hitchon, J.D. Callen, and J.L. Shohet. Flux coordinates
and magnetic field structure. Springer-Verlag, Berlin, 1991.

[DJKR00]

W. Dorland, F. Jenko, M. Kotschenreuther, and B.N. Rogers. Electron temperature
gradient turbulence. Phys. Rev. Lett., 85:5579, 2000.

[Dur99]

D.R. Durran. Numerical methods for wave equations in geophysical fluid dynamics.
Springer-Verlag, New York, 1999.

[FC82]

E.A. Frieman and L. Chen. Nonlinear gyrokinetic equations for low-frequency electromagnetic waves in general plasma equilibria. Phys. Fluids, 25:502, 1982.

[FVV03]

G.L. Falchetto, J. Vaclavik, and L. Villard. Global-gyrokinetic study of finite β effects
on linear microinstabilities. Phys. Plasmas, 10:1424, 2003.

[GC81]

J.M. Greene and M.S. Chance. Nucl. Fusion, 21:453, 1981.

[GST]

S. Gottlieb, C.-W. Shu, and E. Tadmor. Strong stability preserving high order time
discretization methods. SIAM Review, 43:97.

[Gug77]

H.W. Guggenheimer. Differential Geometry. Dover, New York, 1977.

[HW85]

F.L. Hinton and S.K. Wong. Neoclassical ion transport in rotating axisymmetric
plasmas. Phys. Fluids, 28:3082, 1985.

[HW06]

F.L. Hinton and R.E. Waltz. Gyrokinetic turbulent heating. Phys. Plasmas, 13:102301,
2006.

[KC03]

C.A. Kennedy and M.H. Carpenter. Additive runge-kutta schemes for convectiondiffusion-reaction equations. Applied Num. Math, 44:139, 2003.

[KK58]

M.D. Kruskal and R.M. Kulsrud. Equilibrium of a magnetically confined plasma in a
toroid. Phys. Fluids, 1:265, 1958.

[KMJ08]

M. Kammerer, F. Merz, and F. Jenko. Exceptional points in linear gyrokinetics. Phys.
Plasmas, 15:052102, 2008.

[KWC07]

J.E. Kinsey, R.E. Waltz, and J. Candy. The effect of plasma shaping on turbulent
transport and E×B shear quenching in nonlinear gyrokinetic simulations. Phys. Plasmas, 14:102306, 2007.

[LLHL01]

W.W. Lee, J.L.V. Lewandowski, T.S. Hahm, and Z. Lin. Shear Alfvén waves in
gyrokinetic plasmas. Phys. Plasmas, 8:4435, 2001.

[MCG+ 98] R.L. Miller, M.S. Chu, J.M. Greene, Y.R. Lin-liu, and R.E. Waltz. Noncircular, finite
aspect ratio, local equilibrium model. Phys. Plasmas, 5:973, 1998.
[ML74]

C. Mercier and N. Luc. Technical report, Commission of the European Communities,
Report No. EUR-5127e 140, Brussels, 1974.

[OKP+ 00]

M. Ono, S.M. Kaye, Y.-K.M. Peng, G. Barnes, W. Blanchard, M.D. Carter,
J. Chrzanowski, L. Dudek, R. Ewig, D. Gates, R.E. Hatcher, T. Jarboe, S.C. Jardin,
D. Johnson, R. Kaita, M. Kalish, C.E. Kessel, H.W. Kugel, R. Maingi, R. Majeski,
J. Manickam, B. McCormack, J. Menard, D. Mueller, B.A. Nelson, B.E. Nelson,

General Atomics Report GA-A26818

69

C. Neumeyer, G. Oliaro, F. Paoletti, R. Parsells, E. Perry, N. Pomphrey, S. Ramakrishnan, R. Raman, G. Rewoldt, J. Robinson, A.L. Roquemore, P. Ryan, S. Sabbagh,
D. Swain, E.J. Synakowski, M. Viola, M. Williams, J.R. Wilson, and NSTX Team.
Exploration of spherical torus physics in the nstx device. Nucl. Fusion, 40:557, 2000.
[PR00]

L. Pareschi and G. Russo. Implicit-explicit runge-kutta schemes for stiff systems of
differential equations. In L. Brugnano and D. Trigante, editors, Recent Trends in
Numerical Analysis, pages 269–289. Nova Science, 2000.

[PR02]

L. Pareschi and G. Russo. High order asymptotically strong-stability-preserving methods for hyperbolic systems with relaxation. In Hyperbolic problems: Theory, Numerics,
Applications. Springer, 2002.

[RTC82]

G. Rewoldt, W.M. Tang, and M.S. Chance. Electromagnetic kinetic toroidal eigenmodes for general magnetohydrodynamic equilibria. Phys. Fluids, 25:480, 1982.

[SH97]

H. Sugama and W. Horton. Neoclassical electron and ion transport in toroidally
rotating plasmas. Phys. Plasmas, 4:2215, 1997.

[SH98]

H. Sugama and W. Horton. Nonlinear electromagnetic gyrokinetic equation for plasmas with large mean flows. Phys. Plasmas, 5:2560, 1998.

[Sug99]

H. Sugama. Damping of toroidal ion temperature gradient modes. Phys. Plasmas,
6:3527, 1999.

[TLLM+ 99] A.D. Turnbull, Y.R. Lin-Liu, R.L. Miller, T.S. Taylor, and T.N. Todd. Improved magnetohydrodynamic stability through optimization of higher order moments in crosssection shape of tokamaks. Phys. Plasmas, 6:1113, 1999.
[WM99]

R.E. Waltz and R.L. Miller. Ion temperature gradient turbulence simulations and
plasma flux surface shape. Phys. Plasmas, 6:4265, 1999.

[WSCH07] R.E. Waltz, G.M. Staebler, J. Candy, and F.L. Hinton. Gyrokinetic theory and simulation of angular momentum transport. Phys. Plasmas, 14:122507, 2007.
[XMJ+ 08]

P. Xanthopoulos, D. Mikkelsen, F. Jenko, W. Dorland, and O. Kalentev. Verification
and application of numerically generated magnetic coordinate systems in gyrokinetics.
Phys. Plasmas, 15:122108, 2008.

[Zay96]

A. Zayed. Handbook of Function and Generalized Function Transformations. CRC
Press, Boca Raton, 1996.

General Atomics Report GA-A26818

70

Source Exif Data:

File Type                       : PDF
File Type Extension             : pdf
MIME Type                       : application/pdf
PDF Version                     : 1.5
Linearized                      : No
Page Count                      : 74
Page Mode                       : UseOutlines
Author                          : J. Candy
Title                           : GYRO Technical Guide
Subject                         : 
Creator                         : 
Producer                        : 
Create Date                     : 2018:12:03 19:32:54-08:00
Modify Date                     : 2018:12:03 19:32:54-08:00
Trapped                         : False
PTEX Fullbanner                 : This is pdfTeX, Version 3.14159265-2.6-1.40.16 (TeX Live 2015/Debian) kpathsea version 6.2.1

EXIF Metadata provided by EXIF.tools

GYRO Technical Guide Manual

Navigation menu

Versions of this User Manual:

Views

Navigation