11.5 Visualization of gridpoints in field aligned coordinate system

In this section, I try to visualize gridpoints in the field aligned coordinates. The directions of the covariant basis vectors of (ψ,𝜃,α) coordinates are as follows:

(311)

(312)

(313)

Here ∂r∕∂ψ|_𝜃,α is a combination of the usual radial and toroidal direction, which needs some clarification. Note that, ϕ is related to α by Eq. (299), i.e.,

(314)

(where the second equality becomes exact if 𝜃 is the straight-field-line poloidal angle defined in Sec. 6.4.), which indicates that, for q′(ψ)≠0 and 𝜃≠0, the usual toroidal angle ϕ is changing when changing ψ and holding 𝜃 and α fixed. Figure 16b shows how the usual toroidal angle ϕ changes when we change ψ and hold 𝜃 and α fixed.

Fig. 16: (a): A 𝜃 contour on ϕ = 0 plane. Here 𝜃 = 9×2π∕63.  This is the direction of ∂r∕∂ψ|_𝜃,ϕ. (b): ψ coordinate lines in (ψ,𝜃,α) coordinates ( ∂r∕∂ψ|_𝜃,α are the tangent lines to these curves) on the isosurface of 𝜃 = 9 × 2π∕63. Here different lines correspond to different values of α. (c): α coordinate lines (∂r∕∂α|_ψ,𝜃 are tangent lines to these curves), which are along the usual toroidal direction . (d): Grid on the isosurface of 𝜃 = 9×2π∕63, where the red lines are α coordinate lines while the blue lines are ψ coordinate lines. Magnetic field from EAST discharge #59954@3.03s.

The relation ϕ ≈ α + q(ψ)𝜃 given by Eq. (314) indicates that the toroidal shift along ∂r∕∂ψ|_α,𝜃 for a radial change form ψ₁ to ψ₂ is given by (q(ψ₂) − q(ψ₂))𝜃, which is larger on 𝜃 isosurface with larger value of 𝜃. An example for this is shown in Fig. 17 for 𝜃 = 19 × 2π∕63, where ∂r∕∂ψ|_α,𝜃 has larger toroidal shift than that in Fig. 16 for 𝜃 = 9 × 2π∕63.

The ∂r∕∂ψ|_𝜃,α curves can be understood from another perspective. Examine a family of magnetic field lines that start from 𝜃 = 0 and ϕ = ϕ₁ but different radial coordinates. These starting points all have the same value of α, which is equal to ϕ₁. When following these field lines to another isosurfce of 𝜃, the intersecting points of these field lines with the 𝜃 isosurface will trace out a ∂r∕∂ψ|_𝜃,α line with α = ϕ₁. Examine another family of magnetic field lines similar to the above but with the starting toroidal angle ϕ = ϕ₂. They will trace out another ∂r∕∂ψ|_𝜃,α line (with α = ϕ₂) on the 𝜃 isosurface. Continue the process, we finally get those curves in Fig. 16b and Fig. 17b.

Radial wavenumber in field-aligned coordinates (ψ,𝜃,α)

For a harmonic in (ψ,𝜃,ϕ) coordinates given by A(ψ,𝜃,ϕ) = exp(ik_ψψ + im𝜃 + inϕ), the radial wave number is k_ψ. Let us calculate the corresponding radial wavenumber k_ψ^⋆ in the new coordinates (ψ,𝜃,α), which is defined by

(315)

where the phase is given by phase = k_ψψ + m𝜃 + nϕ. Then the above expression is written as

Let us examine how many ψ grid points are needed to resolve the ψ dependence in (ψ,𝜃,α) coordinates on the high-field side (𝜃 = π). Assume k_ψ ≈ 0, then k_ψ^⋆ at 𝜃 = π is given by k_ψ^⋆ = nπq′. The corresponding wave-length is given by λ_ψ^⋆ = 2π∕k_ψ^⋆. The grid spacing Δψ should be less than half of this wave-length (sampling theorem). Then the grid number should satisfy that

(318)

where L_ψ is the radial width of the computational domain.

The number of Fourier harmonics that need to be included is given by

(319)

For DIII-D cyclone base case, choose the radial coordinate ψ as r. At the radial location ψ = r₀ = 0.24m, q₀ = 1.4, ŝ₀ = 0.78, then q′₀ = s₀q₀∕r₀ = 4.5m⁻¹. Then N_r^⋆ = n× 0.45 for the radial width L_r = 0.10m.

Figure 18 plots ∂r∕∂ψ|_𝜃,α lines on the 𝜃 = 0,2π isosurfaces, which are chosen to be on the low-field-side midplane. On 𝜃 = 0 surface, ∂r∕∂ψ|_𝜃,α lines are identical to ∂r∕∂ψ|_𝜃,ϕ lines. On 𝜃 = 2π surface, each ∂r∕∂ψ|_𝜃,α line has large ϕ shift. In old version of my code, 𝜃 = 0,2π surfaces are chosen as the 𝜃 cuts (in the new version 𝜃 = [−π,π]). A connection condition for the perturbations is needed between these two surfaces. This connection condition is discussed in Sec. 11.5.

Periodic conditions of physical quantity along 𝜃 and α in field-line-following coordinates (ψ,𝜃,α)

Since (ψ,𝜃,ϕ) and (ψ,𝜃 + 2π,ϕ) correspond to the same spatial point, a real space continuous quantity f expressed in terms of coordinates (ψ,𝜃,ϕ), i.e., f = f(ψ,𝜃,ϕ), must satisfy the following periodic conditions along 𝜃:

(320)

Since (ψ,𝜃,ϕ) and (ψ,𝜃,ϕ + 2π) correspond to the same spatial point, f must satisfy the following periodic conditions along ϕ:

(321)

Since (ψ,𝜃,α) and (ψ,𝜃,α + 2π) correspond to the same spatial point, a real space continuous quantity g expressed in terms of field-line-following coordinates (ψ,𝜃,α), i.e., g = g(ψ,𝜃,α), must satisfy the following periodic condition along α:

(322)

However, generally there is no periodic condition along 𝜃,

(323)

because P₁ = (ψ,𝜃,α) and P₂ = (ψ,𝜃 + 2π,α) are generally not the same spatial point. In fact, equation (299) implies, for point P₁, its toroidal angle ϕ₁ is given by

(324)

while for point P₂, its toroidal angle ϕ₂ is given by

(325)

i.e., ϕ₁ and ϕ₂ are different by 2πq. From this, we know that (ψ,𝜃,α) and (ψ,𝜃 + 2π,α − 2πq) correspond to the same spatial point. Therefore we have the following periodic condition:

(326)

or equivalently

(327)

Numerical implementation of periodic condition (326)

For the fully kinetic ion module of GEM code that I am developing, 𝜃 is chosen in the range [0 : 2π]. The condition (327) imposes the following boundary condition:

(328)

If α is on a grid, α + 2πq is usually not on a grid. Therefore, to get the value of g(ψ,0,α + 2πq), an interpolation of the discrete date over the generalized toroidal angle α (or equivalently ϕ) is needed, as is shown in Fig. 19.

α contours on a magnetic surface

Figure 20 compares a small number of ϕ contours and α contours on a magnetic surface.

As is shown in the left panel of Fig. 20, with ϕ fixed, an 𝜃 curve reaches its starting point when 𝜃 changes from zero to 2π. However, as shown in the right panel of Fig. 20, with α fixed, an 𝜃 curve (i.e. a magnetic field line) does not necessarily reach its starting point when 𝜃 changes from zero to 2π. There is a toroidal shift, 2πq, between the starting point and ending point. Therefore there is generally no periodic condition along 𝜃 since q is not always an integer. A mixed periodic condition involves both 𝜃 and α is given in (326).

In field-line-following coordinates (ψ,𝜃,α), a toroidal harmonic of a physical perturbation can be written as

(329)

where n is the toroidal mode number, m′, which is not necessary an integer, is introduced to describe the variation along a field line. The periodic condition given by Eq. (326) requires that

(330)

To satisfy the above condition, we can choose

(331)

where N is an arbitrary integer, i.e.,

(332)

We are interested in perturbation with a slow variation along the field line direction (i.e., along ∂r∕∂𝜃|_ψ,α) and thus we want the value of m′ to be small. One of the possible small values given by expression (332) is to choose N = −n × NINT((q_max + q_min)∕2), so that m′ is given by

(333)

where NINT is a function that return the nearest integer of its argument, q_max and q_min is the maximal and minimal value of the safety factor in the radial region in which we are interested. [In the past, I choose m′(ψ) = nq − NINT(nq). However, m′(ψ) in this case is not a continuous function of ψ and thus is not physical.] This form is used to set the initial density perturbation in the fully kinetic code I am developing. Note that m′ in Eq. (333) depends on the radial coordinate ψ through q(ψ). Also note that m′ here is different from the poloidal mode number m in (ψ,𝜃,ϕ) coordinate system. It is ready to show that the perturbation given by Eq. (329) with m′∼ 1 and n ≫ 1 has large poloidal mode number m when expressed in (ψ,𝜃,ϕ) coordinates. [Proof: Expression (329) can be written as

(334)

If 𝜃 is the straight-field-line poloidal angle in (ψ,𝜃,ϕ) coordinate system, then δ(ψ,𝜃) = q𝜃 and the above equation is written as

(335)

which indicates the poloidal mode number m in (ψ,𝜃,ϕ) coordinates is given by m = m′−nq. For the case with m′∼ 1 and n ≫ 1, m is much larger than one.]

Since α contours on a magnetic surface are magnetic field lines, they span out the 3D shape of magnetic surface when there are many α contours on a magnetic surface, as is shown by the right-panel of Fig. 21.

α contours in a toroidal annulus

Figure 22 compares the ϕ coordinate surface of (ψ,𝜃,ϕ) coordinates with the α coordinate surface of (ψ,𝜃,α) coordinates.

 

Fig. 24: The value of 𝜃 and δ on an annulus in the poloidal plane (R,Z). Upper panel: using magnetic coordinates, Lower panel: comparison with the results computed in cylindricall coordinates. Both 𝜃 and δ are discontinuous at the 𝜃 cut, 𝜃 = ±π, which is chosen on the high-field-side in this case. Here δ = ∫ ₀^𝜃 d𝜃.