13.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:

(314)

(315)

(316)

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

(317)

(where the second equality becomes exact if 𝜃 is the straight-field-line poloidal angle defined in Sec. 8.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. (317) 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

(318)

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

(321)

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

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

(322)

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. 13.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 𝜃:

(323)

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

(324)

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 α:

(325)

However, generally there is no periodic condition along 𝜃,

(326)

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

(327)

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

(328)

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:

(329)

or equivalently

(330)

Numerical implementation of periodic condition (329)

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

(331)

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 (329).

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

(332)

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. (329) requires that

(333)

To satisfy the above condition, we can choose

(334)

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

(335)

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 (335) is to choose N = −n × NINT((q_max + q_min)∕2), so that m′ is given by

(336)

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. (336) 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. (332) with m′∼ 1 and n ≫ 1 has large poloidal mode number m when expressed in (ψ,𝜃,ϕ) coordinates. [Proof: Expression (332) can be written as

(337)

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

(338)

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: Values of 𝜃 and δ ≡∫ ₀^𝜃 d𝜃 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.

The above grid is often called “flux-tube”, since it is created by following field lines, and it looks like a tube if the field lines are very near to each other (e.g., when we choose a very narrow radial and toroidal region to start from). Since no magnetic field line passes through the sides of the tube, the flux through any cross section of the tube is equal. The term “flux tube” is often used in astrophysics.