The ability to model computationally the behavior of magnetically confined plasma in the edge region of a tokamak fusion reactor is a key component in the development of a predictive simulation capability for the entire device. As shown schematically in Figure 1, the edge region spans both sides of the magnetic separatrix, encompassing part of the core region, part or all of the region between the plasma and the reactor wall, and the divertor plates. In addition to geometric complexity, an important feature that distinguishes the edge from the core is the development of a region of steep gradients in the density and temperature profiles called the pedestal (Figure 2), the height of which determines the quality of plasma confinement, and hence fusion gain. A kinetic plasma model is needed in this region, because the radial width of the pedestal observed in experiments is comparable to the radial width of individual particle orbits (leading to large distortions of the local distribution functions from a Maxwellian), while the mean free path can be comparable to the scale length for temperature variations along the magnetic field (violating the assumptions underlying a collisional fluid model).
Because of the large number of independent variables in a fully kinetic model, as well as the fast time scale represented by the ion gyrofrequency, gyrokinetic models have been developed that remove the fast gyromotion asymptotically and thereby facilitate numerical treatments. Continuum models consist of a Vlasov operator describing the evolution of plasma species distribution functions in a particular coordinate system combined with some variant of Maxwell’s equations and collision terms. Codes such as GENE[21, 22], GS2[15, 26] and GYRO[4, 6, 5] have been successfully employed to model core plasmas for many years. In addition to requiring simpler geometries, these codes exploit the fact that in the core, distribution functions are typically small perturbations about a known Maxwellian distribution , providing a simpler, and even sometimes linear, model. To model the edge plasma all the way to the reactor walls, a method to solve nonlinear gyrokinetic models for the entire distribution function (so-called full-) in edge-relevant geometries is needed.
In this paper, we describe an approach for the discretization of a continuum full-, gyrokinetic Vlasov system in axisymmetric edge geometries. The system is treated as a conservation law describing the advection of plasma species distribution functions in a four dimensional phase space (2 configuration space and 2 velocity space coordinates). A (mostly) flux-surface-aligned, mapped-multiblock (MMB) coordinate system (Figure 3) is used to accommodate strong anisotropy in the direction of magnetic field lines, similar to fluid edge codes such as UEDGE [31, 30, 29] and SOLPS-ITER . We formulate a high-order spatial discretization to reduce the phase space grid size needed for a specific level of accuracy relative to a lower-order approximation, as well as to reduce numerical dissipation in long-time integrations. In particular, we employ a semi-discretization based on a general formalism [7, 27] for the creation of arbitrarily high-order finite-volume spatial discretizations in mapped coordinates. This formalism is summarized in Sections 3.1 and 3.3. A unique feature of the edge geometry application is that a key requirement of the formalism, namely, that the block mappings are smooth up to and some distance beyond block boundaries, cannot be satisfied by a completely flux-surface-aligned mapping due to the singularity of the metric factors at the X point, where the magnetic separatrix intersects itself and the poloidal field component vanishes. As discussed in Section 4, some dealignment is therefore necessary in a neighborhood of the X point to enable high-order differentiation there.
Since the phase space velocity in gyrokinetic models is developed via
Hamiltonian equations of motion, the zero velocity divergence implied by the
area-preserving property of such dynamical systems is an important
constraint that must be maintained during discretization to avoid the
accumulation of truncation error in long-time integrations. For the
gyrokinetic model  addressed here and presented in
Section 2, we demonstrate in
Section 3.2 how the divergence-free property can be
satisfied to machine precision in the context of our high-order, MMB,
finite-volume discretization. Specifically, in calculating the phase
space cell averages of the gyrokinetic flux divergence, the velocity
enters through integrals of the normal components over the cell faces,
which in turn are used to compute flux normals via a high-order
product formula. By exploiting the fact that the phase space velocity can be
written as a sum of vectors including a certain skew-symmetric,
second-order tensor divergence, Stokes’ theorem can be invoked
to reduce the cell face velocity
integrals to face boundary integrals that telescopically sum to zero
in the cell-averaged finite volume divergence calculation. Additional
benefits accrue from the fact that the face boundary integrals
corresponding to three of the four velocity terms (
how the divergence-free property can be satisfied to machine precision in the context of our high-order, MMB, finite-volume discretization. Specifically, in calculating the phase space cell averages of the gyrokinetic flux divergence, the velocity enters through integrals of the normal components over the cell faces, which in turn are used to compute flux normals via a high-order product formula. By exploiting the fact that the phase space velocity can be written as a sum of vectors including a certain skew-symmetric, second-order tensor divergence, Stokes’ theorem can be invoked to reduce the cell face velocity integrals to face boundary integrals that telescopically sum to zero in the cell-averaged finite volume divergence calculation. Additional benefits accrue from the fact that the face boundary integrals corresponding to three of the four velocity terms (i.e., parallel streaming, curvature and drifts) are computed exactly, and no metric factors appear.
The algorithms described herein provide the foundation of the COGENT (COntinuum Gyrokinetic Edge New Technology) code, which we have been developing for the solution of continuum gyrokinetic models in multiblock geometries, including those describing the tokamak edge. An overview of COGENT is the topic of Section 5, and in Section 6, we utilize COGENT to demonstrate the accuracy of the spatial discretization described in the preceding sections. The efficient discretization of the gyrokinetic Vlasov operator is necessary but not sufficient for the complete simulation of the edge plasma problem. Given its fundamental importance, however, we have limited the scope of this paper to only that part. As described in Section 5, COGENT includes a number of the other components (e.g., collision operators, self-consistent electrostatic fields and high-order semi-explicit time integration methods) needed to address edge-relevant problems, which have been used in several verification studies [16, 11, 13, 9, 14, 12, 17]. To our knowledge, COGENT’s ability to solve a continuum gyrokinetic model in edge geometries spanning both sides of the magnetic separatrix is unique.
2 The gyrokinetic Vlasov system
We target a reduced version of the full- gyrokinetic model of :
The unknown quantity is the plasma species distribution function in gyrocenter phase space coordinates , which are described further below and whose equations of motion are given by (2)-(3). and are the magnitude and direction of the magnetic field , respectively. and are the species charge state and mass, respectively. is the electric potential. For our present purposes, we assume the long wavelength limit in which the Larmor radius is much smaller than the characteristic length scales for electrostatic potential variations. We utilize a particular normalization described in Appendix 9.1 that nondimensionalizes all quantities, including the Larmor number, , defined in Table 5.
Gyrocenter coordinates play a key role in gyrokinetic models in two
important ways. First, they reduce what would otherwise be a
six-dimensional phase space to five dimensions: is the
three-dimensional configuration space coordinate, is the
velocity space component along field lines, and the magnetic moment
is the velocity space component along field lines, and the magnetic momentis related to the velocity perpendicular to field lines. Through the use of asymptotic orderings, gyrocenter coordinates are specifically constructed so as to make the distribution function symmetric with respect to gyrophase. The latter component, which would have been the third velocity component, can then be ignored. The magnetic moment , an adiabatic invariant, is assumed to be constant in the development of gyrokinetic theories, which is why no evolution equation appears for it. The second benefit of gyrocenter coordinates is that the gyrofrequency is eliminated, which would otherwise represent a fast time scale that would need to be resolved.
The gyrokinetic phase space velocity (2) models strong flow along magnetic field lines, represented by the quantity, together with curvature, and drift terms containing the factor, which is assumed small in the gyrokinetic asymptotic ordering. For typical tokamak parameters, the difference in the magnitude of the streaming and drift terms can be a few orders of magnitude. The resulting strong anisotropy resulting from the disparity in parallel and perpendicular quantities has many implications in the development of numerical algorithms. Because gyrocenter coordinates are developed as a Hamiltonian dynamical system, the velocity (2) satisfies the area preserving property
where is the Jacobian of the mapping between lab frame and gyrocenter coordinates. The analytic verification of (4) is contained in Appendix 9.2. As noted in , the gyrokinetic Vlasov equation can therefore be expressed in either convective or conservative form. We choose the latter with the objective of achieving a correspondingly conservative numerical discretization. The potential in (3c) is evaluated by solving some form of Maxwell’s equations. For the purposes of this paper, we assume that is known, although the COGENT code used for our numerical example includes additional options, as discussed in Section 5.
The gyrokinetic system is posed in a domain defined by the tokamak magnetic geometry, which is comprised of field lines lying on concentric flux surfaces (Figure 1). In this paper, we assume an axisymmetric geometry, where all quantities are assumed to be constant in the toroidal direction. The configuration space domain therefore consists of a single poloidal slice. Due to large variations of plasma parameters along and across field lines (and therefore along and across sliced flux surfaces in the poloidal plane), there is strong motivation to discretize in coordinates where at least one of the coordinate directions is defined by the flux surfaces. Since a single, smooth, flux-surface-aligned coordinate mapping cannot be constructed over the entire domain, we consider a multiblock decomposition such as the one depicted in the left-hand side of Figure 3. This decomposition is constructed by first recognizing the natural partitioning defined by the magnetic separatrix into the core, scrape-off and private flux regions. Each region is further decomposed into blocks such that (i) each block can be mapped from a logically rectangular computational domain and (ii) adjacent blocks must abut along entire boundaries. The core region is therefore subdivided into the left (LCORE), middle (MCORE) and right (RCORE) blocks; the scrape-off layer is decomposed into the left (LSOL), left-central (LCSOL), middle-central (MCSOL), right-central (RCSOL) and right (RSOL) blocks, and the private flux region is decomposed into the left (LPF) and right (RPF) blocks. The right-hand side of Figure 3 depicts the inter-block connectivity in the mapped coordinate domain. Within each mapped block, rectangular grids are introduced, resulting in a block rectilinear gridding of the physical domain. We require that the grids be conformal across inter-block boundaries, but otherwise no additional grid smoothness across block boundaries is assumed. For parallelization purposes, the rectangular block grids are further decomposed in 4D, load balanced and assigned to processors in a fully general manner.
We comment that the multiblock decomposition rules imposed above could have been satisfied using only six blocks by combining the LCORE, MCORE and RCORE blocks into a single block and the LCSOL, MCSOL and RCSOL blocks into another block. However, the generation of smooth mappings on the resulting “long and skinny” merged blocks is non-trivial. The advantage of the ten-block decomposition shown in Figure 3 is that eight of the blocks are curvilinear quadrilaterals that are relatively modest deformations of the unit square; the MCORE and MCSOL blocks can be mapped to the unit square using nearly polar coordinate transformations.
On this mapped multiblock grid, we consider the discretization of the
system (1)-(3). Among our requirements
is the discrete enforcement of the zero velocity divergence condition
(4) assumed in the conservative formulation. A second
requirement is high-order (i.e. , greater than second-order)
accuracy, which reduces the number of phase space degrees of freedom
to achieve a given level of accuracy. A high-order method is
also important for reducing numerical dissipation in long-time
integrations. We therefore begin with a review of a general
approach ) is satisfied to
machine precision. In our finite volume approach, the mapped
grid normal velocities are used in flux calculations on cell faces
together with discretized distribution functions. On cell faces near
interblock boundaries, calculation of the latter is accomplished using suitably
, greater than second-order) accuracy, which reduces the number of phase space degrees of freedom to achieve a given level of accuracy. A high-order method is also important for reducing numerical dissipation in long-time integrations. We therefore begin with a review of a general approach for constructing fourth-order, finite-volume discretizations of hyperbolic conservation laws in mapped coordinates on a single block. In this context, we then describe the calculation of the gyrokinetic velocities in the mapped coordinate system such that the divergence free condition (4
) is satisfied to machine precision. In our finite volume approach, the mapped grid normal velocities are used in flux calculations on cell faces together with discretized distribution functions. On cell faces near interblock boundaries, calculation of the latter is accomplished using suitably high-order interpolation, enabling the extension of the discretization to multiblock geometries such as the edge geometry in Figure 3.
3 Phase space discretization
Fundamental to our approach is a general strategy for the systematic development of high-order, finite-volume discretizations in mapped coordinates. Additional details are contained in .
3.1 Fourth-order, finite-volume discretization in a mapped block
Consider a smooth mapping from the unit cube onto the spatial domain :
Given this mapping, the divergence of a vector field on can be written in terms of derivatives in , which will serve as our computational domain. That is,
where denotes the matrix obtained by replacing the row of the matrix by the vector , and denotes the unit vector in the coordinate direction.
In a finite volume approach, is discretized as a union of control volumes. For Cartesian grid finite volume methods, a control volume takes the form
where is the grid spacing. When using mapped coordinates, we define control volumes in as the images of the cubic control volumes . Then, by changing variables and applying the divergence theorem, we obtain the flux divergence integral over a physical control volume by
where the and are lower and upper faces of cell in the -th direction, respectively. As described in , the integrals on the cell faces can be approximated using the following formula for the average of a product in terms of fourth-order accurate face averages of each factor:
Here, is the second-order accurate central difference approximation to the component of the gradient operator orthogonal to the -th direction: , and the operator denotes a fourth-order accurate average over the cell face centered at :
Alternative expressions to (7) are obtained by replacing the averages and/or used in the transverse gradients by the corresponding face-centered pointwise values and/or , respectively.
When applied in the discretization of the phase space divergence operator in the gyrokinetic system (1), we consider the flux , where
and the species subscript is dropped for brevity. We therefore obtain
where, taking and as the factors in the product formula (7),
We note that an alternative discretization can be obtained by preserving the flux as one of the product formula factors, yielding
where is the -th component of and is the -th element of the matrix . In , it is demonstrated that the computation of the face averages can be reduced to integrals over cell edges. Moreover, assuming that the edge integrals are performed with the same quadratures wherever they appear,
which guarantees the freestream property that the divergence of a constant vector field computed by (6) is identically zero. Free-stream preservation is an extremely important property in the simulation of flows using mapped coordinate systems, since it represents a constraint on the approximation of the metric terms that reduces the dependence of computed solutions on the choice of mapping and metric discretization . However, as will be demonstrated in the next section, the factorization (10) enables a discretization in which the divergence of the gyrokinetic velocity vanishes to machine precision, thereby enforcing the assumption (4) made in the conservative formulation (1). Moreover, the formulation is free of discretized metric terms, which naturally achieves one of the goals addressed in free-stream preserving approaches.
Calculation of the face-averaged fluxes (10) is therefore reduced to the calculation of face-averaged distribution functions and mapped normal velocity components. One choice for the former is obtained from the fourth-order, centered-difference formula
where denotes the average of on cell . This results in a dissipationless scheme to which a limiter can also be added. Alternatively, an upwind method with order higher than four, such as the fifth-order WENO scheme , may be employed. Boundary conditions are also implemented here through the setting of inflow conditions in physical boundary ghost cells. Fourth-order face averages of the mapped normal velocity components can be computed directly from (2) using the product formula (7) and metric factor face averages, but as shown in the next section, a more careful exploitation of the specific structure of the gyrokinetic velocity results in a discretization that is also discretely divergence free.
3.2 Gyrokinetic velocity discretization
As shown in Appendix 9.2, the divergence of the gyrokinetic velocity (2) is zero. In this section, we demonstrate how to preserve this property when computing the divergence using the mapped-grid, finite-volume discretization described in the preceding section. Specifically, we show how the mapped normal velocity components can be computed such that the divergence integral (6) with is zero to machine precision.
Let denote a magnetic potential (i.e., ) and let
Letting , we define
where is the vector with entries
is divergence-free, and therefore so is . Using the components of , we define the 3-form
The mapped normal velocity components of can be expressed (see Appendix 9.3) in terms of surface integrals of using the parameterization given by our mapping (which implies a particular surface orientation that has already been accounted for in (6)) as
Since is divergence-free, its exterior derivative is zero. Poincaré’s Lemma (applied on the contractible manifolds ) therefore guarantees the existence of a 2-form
such that . It therefore follows that
where and are the low and high side faces of in direction , respectively. The pullback form is defined by
In terms of coordinates, the pullback integrals are then 
For example, consider a cylindrical coordinate mapping of the form:
where for and . The evaluation of (25) yields
and , , and . Noting that the quantities are independent of for all , we have from (23) that
which, together with (18), yields
and therefore is discretely divergence free.
The assumption of axisymmetric fields allows additional simplifications of the quantities , , and . Assuming that has no toroidal component, we observe that
and similarly, since has no toroidal component,
Hence, since is constant,
Letting denote the magnetic flux, we may take , yielding
Summarizing the above, the face integrals of the mapped normal velocities are
We note that no metric factors appear in , , and , which are computed exactly from the evaluation of magnetic field data at configuration space cell vertices, except for the cell face integrals of the transverse derivatives divided by . These latter integrals can be computed using the fourth-order product formula (7), in which the resulting integrals can be evaluated exactly in terms of nodal values. This fact plays a critical role in avoiding a severe instability that can arise in the high-order, finite-volume modeling of drift waves .
By tracing the individual terms of the gyrokinetic velocity (2) through the above development, the physical meaning of the terms in (34), (35), (36) and (38) can be identified. The and acceleration terms in (38) are clear. The first and second terms in (34) and (35) are the and drifts, respectively. The second term in (36) is the curvature drift. The first term in (36) is the parallel streaming contribution, whose magnitude, as mentioned in Section 2, dominates those of the drift terms by a few orders of magnitude for typical tokamak parameters. In addition to the fact that the parallel streaming, drift and curvature drift terms are all computed exactly, an important consequence of the velocity discretization described here is that, on radial cell faces contained in a flux surface, parallel streaming makes zero (to machine precision) contribution to the normal velocity component (37b), due to the subtraction of the uniform values effected by the sum. This fact ensures that parallel velocity discretization error cannot dominate the drift terms, thereby masking the contribution of the latter via what is sometimes referred to as “numerical pollution” .
3.3 Multiblock extension
The single block discretization just described can be extended to the multiblock edge plasma geometry (Figure 3) using the approach detailed in . In order to apply at block interfaces the same fourth-order reconstruction (13) used in the block interiors, a halo of “extrablock” ghost cells is filled by interpolation from surrounding cell-averaged data. Assuming that the smooth mappings on each block possess smooth extensions beyond their respective boundaries, these extrablock ghost cells are generated by simply applying the mapping to an extended computational grid. Figure 4 shows an example of filling an extrablock ghost cell near the X point. The only aspect the calculation of fluxes at interblock boundaries that is special is a post-processing step performed to restore strict conservation. Since there is no guarantee that the fluxes calculated on each side of a multiblock interface using the above procedure will agree to better than fourth order, the two fluxes are averaged to define a consistent value for use on both adjacent blocks.
The key element is therefore the interpolation of valid neighbor block data to extrablock ghost cells, which we briefly summarize here. Full details are contained in . We use a least-squares approach that allows us to obtain high-order accuracy independent of the degree of smoothness of the grid. We compute a polynomial interpolant in the neighborhood of a ghost cell of the form
We assume that we know the conserved quantities in a collection of control volumes . In that case, we impose the conditions
The integrals on the left-hand side can be computed to any order from the known integrals of the conserved quantities , and the integrals of can be computed directly from the grid mapping. Thus, this constitutes a system of linear equations for the interpolation coefficients . Generally, we choose the number of equations to be greater than the number of unknowns in such a way that the resulting overdetermined system has maximal rank, so that it can be solved using least squares. In the case where we are computing an interpolant onto a finer grid from a coarser one in a locally-refined grid calculation, we impose the conservation condition as a linear constraint.
4 Mapping the edge geometry
The use of mapped coordinates to accommodate strong anisotropy along magnetic field lines requires the construction of block mappings from computational to physical coordinates in which one of the computational coordinates parameterizes flux surfaces. In this section, we describe an approach based on the assumed availability of the magnetic flux in cylindrical coordinates, from which the (assumed axisymmetric) magnetic field is obtained as