Notes on Mapping Field-Line-Following Coordinates to Equivalent Annulus

G. W. Hammett, W. Dorland, E. Belli

Field-line-following coordinates are used in gyrokineitc microturbulence codes such as GS2. While we often think of the simulation domain as being a flux-tube, one can show that this flux-tube is exactly equivalent to a toroidal annular wedge, covering 1/N_{0} of the toroidal angle. Though we often use small gyroradius, large N_{0}, approximations, formally one can consider the N_{0} = 1 special case, and then be simulating the full toroidal extent. Because it is sometimes hard to explain or visualize the flux-tube domain, it is useful just to make plots and movies in the toroidal annulus wedge. These notes are based on Dorland and Kotschenreuther’s “Notes on Local Equilibrium Implementation”, http://gk.umd.edu/g short.pdf (called DK-geo here).

The magnetic field can be written as B⃗ = ∇α × ∇ψ. Many equilibria codes use general (not field-line following) toroidal coordinates (ψ, ϕ, θ). The coordinates used in GS2 are (ψ, α, θ), where ψ is the poloidal flux, α is a second perpendicular coordinate that labels fields line on a flux surface, and the poloidal angle θ is used as a parallel coordinate. α can be written as

α = ϕ − q(ψ)θ − ν(ψ, ϕ, θ)

where ν = ν(ψ, θ, ϕ) must be a periodic function in θ and the toroidal angle ϕ. In an axisymmetric system, simplies to ν = ν(ψ, θ). Furthermore, in the thin flux-tube approximation, we can simplify this to

ν = ν_{0}(θ) + ν_{1}(θ)(ψ − ψ_{0})

where ψ_{0} is ψ in the middle of the flux-tube. Just as dq/dψ is related to the global shear of the magnetic field (and thus to the shear in the flux-tube as one moves along the field line), dν/dψ = ν_{1}(θ) is related to the local shear in the magnetic field. ν_{0}(θ) corresponds to a shift of the whole flux tube in the toroidal direction.

GS2’s geometry and eikcoeffs routines calculate quantities of interest, handling the mapping between various coordinates. In particular, ν_{0}(θ) is determined by Eq.6 of DK-geo, though we point out that the ψ variation of the integrand in Eq. 6 is ignored in the thin flux-tube limit and there is an arbitrary constant of integration so that Eq.6 is really expressing α vs. θ at fixed ψ = ψ_{0} and fixed ϕ = 0, or α(θ, ψ = ψ_{0}, ϕ = 0) = −q_{0}θ − ν_{0}(θ). So Eq. 6 would read

−qθ − ν_{0}(θ) = ∫_{0}^{θ} dθ\frac{B_{0} · ∇ϕ}{∇θ × ∇Ψ · ∇ϕ}  (1)

This can be used to determine ν_{0}(θ).

To determine ν_{1}θ, we use

∇α = \frac{∂α}{∂ψ}∇ψ + \frac{∂α}{∂ϕ}∇ϕ + \frac{∂α}{∂θ}∇θ

Dotting this equation with ∇ϕ × ∇θ (related to a contravariant basis), we get

\frac{∂α}{∂ψ} = \frac{∇ϕ × ∇θ · ∇α}{∇ϕ × ∇θ · ∇ψ}

(the denominator is related to the Jacobian). We can now determine ν_{1} from

ν_{1}(θ) = \frac{∂ν}{∂ψ} = −\frac{∂α}{∂ψ} − \frac{∂q}{∂ψ} θ  (2)

GS2’s geometry routines can be used to evaluate ν_{0}(θ) and ν_{1}(θ) using the above equations. Once these quantities are written to a file, the mapping from flux-tube coordinates to toroidal annulus coordinates for visualization can be done in a postprocessor.

1 Post-processing mapping

We want to map quantities such as the temperature T (ψ, α, θ) from GS2’s flux-tube coordinates to annulus coordinates for visualizations, which is straightforwardly:

T_{annulus}(ψ, ϕ, θ) = T (ψ, α(ψ, ϕ, θ), θ)

(for visualization packages, one would also provide a separate set of arrays giving the x, y, z cartesion grid locations for each grid point in (ψ, ϕ, θ), but that is relatively straightforward and I won’t describe it here.) For visualizations, we will have a regular grid in (ψ, ϕ, θ), with ϕ spanning the range (0, 2π/N_{0}), and the one little trick is that to determine the corresponding α we make use of periodicity:

α = (ϕ − q(ψ)θ − ν(ψ, θ)) Mod (2π/N_{0})

where N_{0} is the fraction of the toroidal angle covered by the flux tube. Formally, the N_{0} parameter is arbitrarily big in the small ρ_{∗}, thin flux-tube limit. One can either determine N_{0} from the actual ρ_{∗} = ρ/a of the experiment, or one can pick a “typical” value of N_{0} = 5, at least to start with.

An interesting point to note is that the data mapped to an annulus requires many more points to represent it than the original flux-tube data. This is because the original flux-tube data T (ρ, α, θ) is known on an N x*N y*N theta grid, where typical parameters for a moderate resolution GS2 run are 64*32*16). But the data mapped on an annular grid needs to be on a grid of order:

N r = N x ∼ 64

N phi = N y ∼ 32

N theta = N y*N 0*q ∼ 32*5*2 ∼ 320

I.e., as θ goes from zero to 2*pi, the alpha coordinate wraps around it’s full range qN_{0} times. I suppose technically one should choose N theta=max(N y*N 0*q, N z), but I’m assuming that N z will always be negligible by comparison.

2 Details Specific to GS2 Calculation

The key reference to bring together the preceding and what is actually calculated in gs2 is C. M. Bishop, et al., NF, 24, 1579 (1984). To make headway easier, it is useful to note a couple of identities. We denote the major radius by R; Bishop uses

R ≡ X_{0}h_{0}.

Bishop defines R to be the poloidal curvature of the flux surface in the poloidal plane. In the geometry modules, I call this R_{pol}. Note that R_{pol} = R_{pol}(θ).

The key result for the present purposes is Eq. (9) of Bishop, which in my notation, and with a couple of typos corrected (according to me) reads:

α = ϕ − ∫_{0}^{θ} \frac{dθ}{∇Ψ × ∇θ · ∇ϕ} \{\frac{I}{R^{2}B_{p}} \}

+ρB_{p}IR ∫_{0}^{θ} \frac{dθ}{∇Ψ × ∇θ · ∇ϕ} \frac{1}{R^{3}B_{p}} \{\frac{B_{p}RI′}{I} + \frac{2 sin u}{R} + \frac{2}{R_{pol}} + \frac{Rp′}{B_{p}} + \frac{II′}{RB_{p}} \}

From this expression, we can make a close correspondence with Greg’s expressions. For example, the term after ϕ is exactly ν_{0}. I’m confused by a detail, though: Bishop does not explicitly write down the secular-in-θ part of α. That is, there is no term in Bishop’s Eq. (9) like −qθ. Pushing on, though, the next term is in fact Greg’s ν_{1}, again leaving out the secular part:

ν_{1} = \frac{∂α}{∂ρ} = B_{p}IR ∫_{0}^{θ} \frac{dθ}{∇Ψ × ∇θ · ∇ϕ} \frac{1}{R^{3}B_{p}} \{\frac{B_{p}RI′}{I} + \frac{2 sin u}{R} + \frac{2}{R_{pol}} + \frac{Rp′}{B_{p}} + \frac{II′}{RB_{p}} \} .

The point of this addendum is point out that this expression offers an alternative way to find ν_{1}. I just have to add back in the secular part.

These particular expressions are already explicitly calculated in the geometry module. The remaining issues in my mind are relatively easy to get straight, and include getting the output to make sense for any choice of radial coordinate. The notes above are rigorously correct only for when the poloidal flux is used, while most gs2 users use other choices.

In principle, there is also the problem of calculating the grid quantities when the Bishop relations aren’t being used. However, I have recommended against using bishop = 0 for several years, because the expressions that we got from this paper are much more accurate and well-behaved than the alternative, in every case. So, I’m going to ignore the implementation of grid-writing routines for this choice.