1 Introduction

A routine in tokamak operation is to reconstruct magnetic field under the constraints of MHD force balance and various measurements. This kind of task can be done by various codes, e.g., EFIT. I developed a similar code, heq (github.com/youjunhu/heq). I will discuss details related to this code.

I also discuss the fixed-boundary equilibrium problem, where the boundary magnetic surface is given and one is asked to solve the force-balance equation within that boundray.

Another topic is how to construct spatial coordinates that conform to the shape of the magnetic surfaces and makes field line equation simple by defining suitable poloidal/toroidal angles, based on discrete numerical equilibrium data output by the equilibrium reconstructing codes.

Let us begin with the basics.

2 Axisymmetric magnetic field

Due to the divergence-free condition, i.e., ∇⋅ B = 0, magnetic field can be expressed as the curl of a vector field,

(1)

where A is called the vector potential of B. (The usefullness of this representation is that, once B is in this form, we do not need to worry about the divergence-free condition.) In cylindrical coordinates (R,ϕ,Z), the above expression is written

(2)

We consider axisymmetric magnetic field. The axial symmetry means that, when expressed in the cylindrical coordinate system (R,ϕ,Z), the components of B, namely B_R, B_Z, and B_ϕ, are all independent of ϕ. For this case, it can be proved that an axisymmetric vector potential A suffices for expressing the magnetic field, i.e., all the components of the vector potential A can also be taken independent of ϕ. Using this, Eq. (2) is written

(3)

In tokamak literature, direction is called the toroidal direction, and (R,Z) planes (i.e., ϕ = const planes) are called poloidal planes.

2.1 Poloidal magnetic field

Equation (3) indicates that the two poloidal components of B, namely B_R and B_Z, are determined by a single component of A, namely A_ϕ. This motivates us to define a function Ψ(R,Z):

(4)

Then Eq. (3) implies the poloidal components, B_R and B_Z, can be written as

(5)

(6)

(Note that it is the property of being axisymmetric and divergence-free that enables us to express the two components of B, namely B_R and B_Z, in terms of a single function Ψ(R,Z).) Furthermore, it is ready to prove that Ψ is constant along a magnetic field line, i.e. B ⋅∇Ψ = 0. [Proof:

We note that Ψ is related to the poloidal magnetic flux, as we will discuss in Sec. 2.7.

Using Eqs. (5) and (6), the poloidal magnetic field B_p is written as

2.2 Toroidal magnetic field

Next, let’s examine the toroidal component B_ϕ. Equation (3) indicates that B_ϕ involves both A_R and A_Z. This indicates that using the potential form does not enable useful simplification for B_ϕ. Therfore we will directly use B_ϕ. Define g ≡ RB_ϕ(R,Z) (the reason that we define this quantity will become clear when we discuss the forece balance), then the toroidal magnetic field is written as

(9)

2.3 General form of axisymmetric magnetic field

Combining Eqs. (8) and (9), we obtain

2.4 Gauge transform of Ψ

Let us discuss the gauge transformation of A in the axisymmetric case. As is well known, magnetic field remains the same under the following gauge transformation:

(11)

where f is an arbitrary scalar field. Here we require that ∇f be axisymmetric because, as mentioned above, an axisymmetric vector potential suffices for describing an axisymmetric magnetic field. In cylindrical coordinates, ∇f is given by

(12)

Since ∇f is axisymmetric, it follows that all the three components of ∇f are independent of ϕ, i.e., ∂²f∕∂R∂ϕ = 0, ∂²f∕∂Z∂ϕ = 0, and ∂²f∕∂ϕ² = 0, which implies that ∂f∕∂ϕ is independent of R, Z, and ϕ, i.e., ∂f∕∂ϕ  is actually a spatial constant. Using this, the ϕ component of the gauge transformation (11) is written

(14)

which indicates Ψ is similar to that of electric potential, i.e., adding a constant to it does not make difference. Note that the definition Ψ(R,Z) ≡ RA_ϕ does not imply Ψ(R = 0,Z) = 0 because A_ϕ can adopt 1∕R dependence under the gauge transformation (13). If A_ϕ is finite at R = 0, then Ψ is zero there. This is the case we encounter in equilibrium reconstruction.

2.5 Contours of Ψ in the poloidal plane

Because Ψ is constant along a magnetic field line and Ψ is independent of ϕ, the projection of a magnetic field line onto (R,Z) plane is a contour of Ψ. Inversely, a contour of Ψ is also projections of a magnetic field line onto the plane. [Proof. A contour of Ψ on (R,Z) plane satisfies

(15)

i.e.,

(16)

(17)

Using Eqs. (5) and (6), the above equation is written

(18)

i.e.,

(19)

which is the equation of the projection of a magnetic field line in (R,Z) plane. This proves that contours of Ψ are projections of magnetic field lines in (R,Z) plane.]

Figure 1 shows typical contours of Ψ in a tokamak (CFETR tokamak as an example).

Points where the poloidal field is zero (i.e., ∇Ψ=0) are called magnetic null points. There are two types of magnetgic null points : O-points and X-points, which can be visually identified by viewing contours of Ψ. Mathematically, O-points and X-points are distinguished by the sign of S defined by

(20)

where S ≥ 0 corresponds to O-points and S < 0 corresponds to X-points.

2.6 Magnetic surfaces

Surfaces of revolution generated by rotating Ψ contours around the axis of symmetry (Z axis) are called magnetic surfaces or flux surfaces. No field line intersects these surfaces. We are only interested in flux surfaces within the machine wall. We note that some Ψ contours may intersect the wall before they can form closed curves. Field lines on the corresponding flux surfaces are called open field lines.

The value of Ψ is constant on a magnetic surface. Meanwhile, the values of Ψ on different magnetic surfaces are usually different. These two properties enable Ψ to be used as labels of magnetic surfaces. In cases that there are multiple magnetic surfaces of the same value of Ψ, the ambiguity can be resolved by specifying which region the flux surfaces lie in, e.g., within or outside the last-closed-flux surface region, near the high-field side or low-field side, in the Scrape-Off Layer or the private flux region (PFR).

PFR:  the area between the X-point and the material divertor., i.e, the region between divertor legs that is unconnected to the plasma.

2.7 Relation of Ψ with the poloidal magnetic flux

Note that Ψ is defined by Ψ = RA_ϕ, which is just a component of the vector potential A, thereby having no obvious physical meaning. Next, we show that Ψ has a simple relation with the poloidal magnetic flux that can be be measured in experiments.

In Fig. 2, there are two magnetic surfaces labeled, respectively, by Ψ = Ψ₁ and Ψ = Ψ₂. The poloidal magnetic flux through any toroidal ribbons between the two magnetic surfaces is equal to each other (due to the magnetic Gauss theorem). Next, we calculate this poloidal magnetic flux, denoted by Ψ_p⁽¹²⁾. To make the calculation easy, we select a plane perpendicular to the Z axis, as is shown by the dash line in Fig. (2). In this case, only B_Z contribute to the poloidal magnetic flux: (the positive direction of the plane is chosen in the direction of )

In experiments, we measure the poloidal magnetic flux through toroidal loops around the central symmetric axis (discussed in Sec. 2.8). Consider one of the loops that located at point (R,Z), then, using (21), the flux through the loop can be written as

(22)

The central symmetric axis (R = 0) is a field line. If A_ϕ is finite at R = 0, then Ψ = A_ϕR is zero there. This is the case we encounter in equilibrium reconstruction. Then this flux is written as

(23)

Due to this relation, Ψ is often called the poloidal magnetic flux per radian (SI unit: web/rad). This relation allows the poloidal flux measurements to be used to constrain the GS equation in the equilibrium reconstrunction (discussed in Sec. 6). Note that the positive normal direction of the surface (where the magnetic flux Ψ_p is defined) is chosen in the +Z direction.

2.8 Measuring poloidal magnetic flux in experiments

Tokamaks usually have some toroidal wire loops around the central symmetirc axis, called flux loops, which are used to measure the poloidal magnetic flux through the loops. (By measuring the voltage around the loops, we can obtain the time derivative of the flux and, after integrating over time, the flux itself.)

Suppose that the loop is located at (R,Z) and denote the magnetic flux through the loop by Ψ_p(R,Z,t) (only the poloidal magnetic field contribute to this flux since the loop is in the toroidal direction). Then Faraday’s law gives

(24)

where 𝜀 is the emf. If the loop is a coil with N turns, the induced voltage V in the coil is N times the emf 𝜀, i.e., V = N𝜀. Using this, Eq. (24) is written as

(25)

Integrating the above equation over time, we obtain

(26)

The starting time t = 0 can be chosen as when Ψ_p(R,Z,0) is easy to know (e.g., when there is no plasma). Equation (26) tells us how to calcualte Ψ_p from the measured loop voltage V .

There are usually many flux loops (e.g. 35 on EAST[32]) at different locations in the poloidal plane (see Fig. 3). They are outside of the plasma region and thus are  called “external magnetic measurements”. The measured poloidal flux, along with the poloidal field measurement by magnetic probes, can be used as constraints in reconstructing the magnetic field within the plasma region. This is discussed in Sec. 6.

2.9 Closed magnetic surfaces

In most part of a tokamak plasma, contours of Ψ in (R,Z) plane are closed before they touch the machine wall. Figure 4 shows some examples of closed flux surfaces.

The innermost magnetic surface reduces to a curve, which is called magnetic axis (in Fig. 4, Ψ₀ labels the magnetic axis). ∇Ψ is zero at the magnetic axis since Ψ(R,Z) reach maximum/minimum there. As a result, the poloidal magnetic field is zero there (refer to Eq. (8)).

As is discussed in Sec. 2.7, the poloidal magnetic flux enclosed by a magnetic surface Ψ (the poloidal magnetic flux through the toroidal surfaces S₂) is given by

(27)

Here the positive direction of the surface S₂ is defined to be in the clockwise direction when an observer looks along the direction of . In practice, we need to pay attention to the positive direction chosen (there can be a sign difference when choosing different positive directions).

Also note, the poloidal magnetic flux enclosed by a closed magnetic surface can have another definitions: the poloidal magnetic flux through the central hole of the magnetic surface, i.e., the poloidal flux through S₁ in Fig. 4. In this case, as is discussed in Sec. 2.7, the poloidal magnetic flux is related to Ψ by

(28)

where the positive direction of the surface S₁ is in the clockwise direction.

Also note that, since B_p = ∇Ψ ×∇ϕ, the condition Ψ_LCFS − Ψ_axis > 0 means B_p points in the anticlockwise direction (viewed along ϕ direction), and Ψ_LCFS − Ψ_axis < 0 means B_p points in the clockwise direction.

3 Non-axisymmetric magnetic perturbations

In the above, the magnetic field is assumed to be axisymmetric. With this assumption, the poloidal magnetic field (having two components) can be expressed in terms of a single component of the vector potential A, A_ϕ (specifically via Ψ ≡ A_ϕR). This kind of simplification is not achievable if the axisymmetricity assumption is dropped, because other components of the vector potential (namely A_R and A_Z) will appear in the expression of the poloidal magnetic field. Let us re-examine Eq. (2) for a general magnetic perturbation:

(30)

(31)

(32)

where δΨ = RδA_ϕ. Therefore this kind of magnetic perturbation can still be written in the same form of the equilibrium poloidal magnetic field:

(33)

The above approximation is widely used in practice, e.g., in turbulence simulation, where δA_ϕ is replaced by δA_∥. (Do we miss some magnetic perturbations that is important for plasma transport when using the above specific form?)

The total magnetic field is then written as

(34)

**check**Can the projection of the total magnetic field line in the poloidal plane can be traced by tracing the contour of Ψ + δΨ? No. The contours of Ψ + δΨ will not show island structures in the poloidal plane. To show the expected island structures, we need to subtract non-reconnecting poloidal magnetic field from the total poloidal field? **check** The contours of the so-called helical flux will give the expected island structures near the resonant surfaces?**check**

Next, we go back to discuss the 2D case (i.e., assuming axisymmetry).

4 Plasma current density in terms of Ψ and g

When the displacement current term is neglectable (the case we consider here), the conductive current is just another representation of the magnetic field. Specifically, the current density is proportional to the curl of the magnetic field (Ampère’s law):

4.1 Poloidal current density

Use Eq. (35) and the definition g ≡ RB_ϕ, the poloidal components of the current density, J_Z and J_R, can be written as

(36)

and

(37)

respectively.

4.2 Toroidal current density

Ampere’s law (35) indicates the toroidal current density J_ϕ is given by

(39)

which is the Laplace operator in cylindrical coordinates for the axisymmertic case, then Eq. (38) is written as

(40)

(Many authors refer to Eq. (40) as the Grad-Shafranov (GS) equation. No, it is not. It is just Ampere’s law, which has nothing to do with the force-balance. Only after you express J_ϕ in terms of the plasma pressure, can Eq. (40) be called the GS equation, as is discussed Sec. 5.3.)

5 Constraint of force-balance on magnetic field

Let us consider what constraint the force balance imposes on the axisymmetric magnetic field discussed above. The MHD momentum equation is given by

(41)

where ρ, ρ_q, ℙ, J, E, and B are mass density, charge density, thermal pressure tensor, current density, electric field, and magnetic field, respectively. The electric field force ρ_qE is usually ignored due to either ρ_q = 0 or E = 0. Further assume that there is no plasma flow (u = 0, the flow effect is discussed in A.4) and the plasma pressure is isotropic, then the steady state momentum equation (force balance equation) is written as

(42)

where P is the scalar plasma pressure.

Is the force balance (42) always satisfied in a real toakamak discharge? To answer this question, we need to go back to the original momentum equation (41). The imbalance between J × B and ∇P will give rise to the compressional Alfven waves, the time-scale of which, τ_A, is much shorter than the time-scale τ we are interested in. Therefore, on the time scale τ (and for slow flow with u < C_s, where C_s is the the sound speed), the leading order of the momentum equation is the force balance (42).[26].

5.1 Parallel force balance

Consider the force balance in the direction of B. Dotting the equilibrium equation (42) by B, we obtain

(43)

which implies that P is constant along a magnetic field line. Since Ψ is also constant along a magnetic field line, P can be expressed in terms of only Ψ on a single magnetic line. Note that this does not necessarily mean P is a single-valued function of Ψ, (i.e. P = P(Ψ)). This is because P still has the freedom of taking different values on different magnetic field lines with the same value of Ψ while still satisfying B ⋅∇P = 0. This situation can appear when there are saddle points (X points) in Ψ contours (refer to Sec. A.13) and P takes different functions of Ψ in islands of Ψ sepearated by a X point. For pressure within a single island of Ψ, P = P(Ψ) can be safely assumed.

On the other hand, if P = P(Ψ), then we obtain

B ⋅∇P = B ⋅∇Ψ = 0,

i.e., Eq. (43) is satisfied, indicating P = P(Ψ) is a sufficient condition for the force balance in the parallel (to the magnetic field) direction.

5.2 Toroidal force balance

Consider the force balance in the toroidal direction. The ϕ component of Eq. (42) is written

(44)

Since P = P(Ψ), which implies ∂P∕∂ϕ = 0, equation (44) reduces to

(45)

Using the expressions of the poloidal current density (36) and (37) in the force balance equation (45) yields

(46)

which can be further written

(47)

According to the same reasoning for the pressure, we conclude that g = g(Ψ) is a sufficient condition for the toroidal force balance. (The function g defined here is usually called the “poloidal current function” in tokamak literature. The reason for this name is discussed in Sec. A.5.)

5.3 Force balance along the major radius

Consider the force balance in direction. The component of Eq. (42) is written

(48)

Using the expressions of the current density and magnetic field [Eqs. (6) and (37)], equation (48) is written

(49)

Assuming the sufficient condition discussed above, i.e., P and g are a function of only Ψ, i.e., P = P(Ψ) and g = g(Ψ), Eq. (49) is written

(50)

which can be simplified to

(51)

which is the requirement of force-balance along the major radius. On the other hand, we know J_ϕ can be expressed in Ψ via Eq. (40). Combining this with Eq. (51) yields

(52)

i.e.,

(53)

Equation (53) is known as Grad-Shafranov (GS) equation.

[Note that the Z component of the force balance equation is written

JRBϕ − JϕBR = ⇒− −△∗Ψ = μ0 ⇒− −△∗Ψ = μ0 ⇒− −△∗Ψ = μ0 ⇒△∗Ψ = −μ0R2 −g

which turns out to be identical with the Grad-Shafranov equation. This is not a coincidence. The reason is that the force balance equation has been satisfied in three different directions (namely, , , and B direction) and thus it must be satisfied in all the directions.]

5.4 Axisymmetric equilibrium magnetic field

A general axisymmetric magnetic field (which does not necessarily satisfy the force balance), is given by Eq. (10), i.e.,

(54)

For the above axisymmetric magnetic field to be consistent with the force balance equation (42), there are additional requirements for Ψ and g. Specifically, Ψ is restricted by the GS equation and g should be a function of only Ψ. Therefore an axisymmetric equilibrium magnetic field is fully determined by two functions, Ψ = Ψ(R,Z) and g = g(Ψ). The Ψ is determined by solving the GS equation with specified RHS source terms and boundary conditions.

The RHS source terms in the GS equation (53) are P(Ψ) and g(Ψ), both of which must be specified before the GS equation can be solved. For most cases, the source terms are nonlinear about Ψ and thus the GS equation is a two-dimensional (in R and Z) nonlinear partial differential equation for Ψ.

For most choices of P(Ψ) and g(Ψ), the GS equation (53) has to be solved numerically. For some particular choices of P and g profiles, analytical solutions are available, one of which is the Solovév equilibrium and is discussed in Appendix A.3.

Note that we solve the GS equation in order to obtain the poloidal magnetic flux Ψ and thus the poloidal magnetic field. The toroidal magnetic field must be specified in some way before we can solve the GS equation. There are several ways of specifying the toroidal magnetic field: (1) given g(Ψ), (2) given ⟨j_∥⟩, (3) given the safety factor q(Ψ). There are simple relations between g, ⟨j_∥⟩, and q, which allows translation form one to another (discussed later). In transport simulations, ⟨j_∥⟩ is obtained from current drive models and neoclassical bootstrap current models. Note that the specification of the source terms (P, g, q, and ⟨j_∥⟩) usually involve the unknown Ψ (via not only the explicit presence of Ψ, but also the flux-surface averaging which implicit involves Ψ). This indicates that iterations are needed when numerically solving the GS equation.

5.5 Axisymmetric equilibrium current density

Since J = μ₀⁻¹∇× B, the current density J can be inferred from a given magnetic field. The components of J (expressed in terms of g and Ψ) are given by Eqs. (36), (37) and (40) and these expressions can be further simplified by using the equilibrium constraints, such as g = g(Ψ) and △^∗Ψ = −μ₀R² −g. The simplified J expressions are given in A.11.

On the other hand, in the kinetic equilibrium reconstruction[20], it is the current that is first (partially) specified (e.g., by summing sources of current drive and bootstrap current) and then the current is used as constraints for the GS equation, i.e, constraints for the magnetic field.

5.6 Equilibrium scaling

The GS equation is given by Eq. (53), i.e.,

(55)

If a solution to the GS equation is obtained, the solution can be scaled to obtain a family of solutions. Given an equilibrium with Ψ(R,Z), P(Ψ), g(Ψ), then it is ready to prove that Ψ₂ = sΨ, P₂ = s²P(Ψ), and g₂ = ±sg(Ψ) is also a solution to the GS equation, where s is a constant. In this case, both the poloidal and toroidal magnetic fields are increased by a factor of s, and thus the safety factor remains unchanged. Also note that the pressure is increased by s² factor and thus the value of β (the ratio of the therm pressure to magnetic pressure) remains unchanged. Note that g₂ = ±sg(Ψ), which indicates that the direction of the toroidal magnetic field can be reversed without breaking the force balance. Also note that Ψ₂ = sΨ and s can be negative, which indicates that the direction of the toroidal current can also be reversed without breaking the force balance.

The second kind of scaling is to set Ψ₂ = Ψ, P₂ = P(Ψ), and g₂² = g²(Ψ) + c. It is ready to prove that the scaled expression is still a solution to the GS equation because g₂g₂′ = gg′. This scaling keep the pressure and the poloidal field unchanged and thus the poloidal beta β_p remains unchanged. This scaling scales the toroidal field and thus can be used to generate a family of equilibria with different profiles of safety factor.

Another scaling, which is trivial, is to set Ψ₂ = Ψ, P₂ = P(Ψ) + c, and g₂ = g(Ψ). This scaling can be used to test the effects of the pressure (not the pressure gradient) on various physical processes.

When a numerical equilibrium is obtained, one can use these scaling methods together to generate new equilibria that satisfy particular global conditions. Note that the shape of magnetic surfaces of the scaled equilibrium remains the same as the original one.

The above scaling is made under the constraint that the the GS equation (55), is satisfied. In practice, we may scale Ψ by a factor while fixing g and P. This does not satisfy the GS equation, but allows more flexibility in changing the the safety factor profile. Scaling Ψ by a factor corresponds to scaling the plasma current and the poloidal magnetic field.

6 Free-boundary fitting problem

Ampere’s law in Eq. (40) can be generalized to include discrete toroidal currents:

(57)

where G is the fundamental solution given by

If all the toroidal currents (plasma current + PF coil curents) are known, then Ψ can be calculated by using Eq. (57). Unfortunately J_ϕ is usually unknown. On the other hand, the force-balance indicates that J_ϕ can be expressed in terms of Ψ via Eq. (51), i.e.,

(59)

Substituting this into Eq. (57) gives an implicit formula for Ψ, which can be iterated (with an initial guess of Ψ(R,Z), and assuming the function forms, P(Ψ) and g(Ψ), as well as the coil currents I_i, are known ). Will the iteration converge? We do not know for sure. Numerical experiments indicate it does in some cases[18]. Let us discuss these cases. We assume that the 1D functions, P(Ψ) and g(Ψ) can be modeled as (following EFIT’s model):

(60)

(61)

where

(62)

Ψ_M and Ψ_B are values of Ψ at the magnetic axis and LCFS, the coefficients α_j and β_j are to be determined, N_p and N_F are integers chosen by users. One feature of expressions (60) and (61) is that they guarantee that dP∕dΨ and gdg∕dΨ are zero at and outside the LCFS, and thus no plasma current there. Also note that expressions (60) and (61) are nonlinear functions of Ψ even for N_p = N_F = 1. This is because the unknowns, Ψ_M and Ψ_B, appear in the denominator of expression (62). Another nonlinearity is related to the fact that the unkonw Ψ determines where is the LCFS and thus determine region where the current desnsity is set to zero. This also gives rise to the name “free-boundary” when referring to this kind of problem.

Choosing an initial guess of Ψ(R,Z) on gridpoints, we can get the values of α_j, β_j, and I_i by solving a least square problem that minimises the difference between quantities computed and the corresponding quantities measured in actual experiments. (Details are given in Sec. 6.1.)

After the coefficients α_j and β_j are obtained, J_ϕ can be updated by using Eqs. (59)-(62). Then, we use the latest J_ϕ and I_i in Eq. (57) to update the values of Ψ on gridpoints. We repeat the above procedure until convergecne in Ψ(R,Z). This is called Picard iteration. (Alternatively, the values of Ψ on the inner gridpoints can be updated by inverting the Laplace operator △^⋆, see Sec. 6.4. The values of Ψ on the bounary still have to be obtained by using Eq. (57).)

I implemented the above method in a Python code HEQ (https://github.com/youjunhu/heq). The following subsections discuss details in the code: Sec. 6.1 discusses the least square problem. Sec. 6.3 discusses the vertical displacement instability stabilizer. Without the stabilizer, the Picard iteration diverges for many elongated configurations, due to the vertical displacement instability.

6.1 Response matrix Γ

Plugging expressions (60) and (61) into (59), we obtain

(64)

and the corresponding expansion coefficients

(65)

then Eq. (63) is written as

(66)

where N_b = N_P + N_F. Plugging expression (66) into Eq. (57), we obtain

(68)

for j = 0,…,N_b − 1, and

(69)

for j = N_b,…,N − 1. Here (R_i^(FL),Z_i^(FL)) is the location where the ith flux measurement is made.

The matrix Γ , which is of shape (M,N), is called “response matrix”. (In least square problems, this matrix is called “design matrix’.) The value of M is equal to number of measurements included (measurements are merged to Γ by row stacking). We will include M_FL measurements of the poloidal flux, M_B measurements of the poloidal magnetic field, 1 measurement of plasma current, a constraint for the on-axis safety factor. Then M = M_FL + M_B + 1 + 1. (For EAST, M_FL = 35,M_B = 38.)

Next, we consider the poloidal field measurement:

(70)

where (R_i,Z_i) is the location of the No. i probe, _i is a unit vector denoting the positive normal direction of the No. i magnetic probe, cos𝜃_pi = _i ⋅ e_R, sin𝜃_pi = _i ⋅ e_Z.

B_p(R_i,Z_i) can be inferred from u via

where

(72)

and G_BR and G_BZ is given by Eq. (681) and (682). The rhs of Eq. (71) indicates that the matrix elements of Γ are given by:

(73)

for (j = 0,…,N_b − 1), and

(74)

for (j = N_b,…,N − 1) . I place the magnetic field equations after the flux equations, i.,e the row number i is in the range [M_FL,M_FL + M_B − 1].

Next, we consider the constraint of plasma current. The plasma current I_p can be inferred from u via

(75)

The rhs of Eq. (75) indicates that the matrix elements of Γ are given by:

(76)

for (j = 0,…,N_b − 1), and zero for (j = N_b,…,N − 1) . This equation is placed after the equations of flux and magnetic field measurement, i.e., its row number is  i = M_FL + M_B.

Let us consider the constraint of on-axis safety factor value. A formula for computing q value at the magnetic axis is given by [ref. P. M. Bellan’s paper]:

(77)

where

(78)

is the elongation at the magnetic axis. Using

expression (77) is written as

(79)

In HEQ, this constraint is not included in some cases.

6.2 Algorithm summary

Step 0. Initial guess of Ψ: Ψ^(k) with k = 0.

Step 1. Compute elements of response matrix Γ using Ψ^(k).

Step 2. Solve linear least-square problem to get (c₁,c₂,c₂,c₄,I₁,…,I_Nc).

Step 3. Update Ψ using the Picard iteration:

(80)

or in its discrete form:

(81)

Step 4: k = k + 1 and goto Step 1.

In Step 1, we first search for the magnetic axis and boundary surface of Ψ^(k) so that we can get the values of Ψ^(k) at those locations. These values are needed in computing x = (Ψ^(k) − Ψ_M)∕(Ψ_B − Ψ_M).

6.3 Vertical Displacement Event (VDE) stabilizer in free-boundary solver[17]

To stabilize the vertical displacement instability in the Picard iteration, we add two virtual coils that carry oppositive currents (denoted by I_vc and −I_vc), and are up-down symmetric, located at (R_vc,Z_vc) and (R_vc,−Z_vc), respectively. And set the current I_vc by

Ivc = −gz,

where (R_cur,Z_cur) is the location of plasma current center defined by

Rcur = ∫ ΩplRJϕ(R,Z)dRdZ,

Zcur = ∫ ΩplZJϕ(R,Z)dRdZ,

B_R,vacuum is the B_R generated by all the real PF coils, g_z is a positive constant chosen by users. Note that g = 1 corresponds to that the B_R generated by the virtual coils exactly cancles the B_R generted by all the real coils. Usually we need value larger than 2, i.e., over compensation. The virtual coils locatons, (R_vc,Z_vc) and (R_vc,−Z_vc), are usually chosen near the center of the top/bottom boundary of the computational box.

6.4 Solving by matrix inversion

If we want to invert the finite-difference version of the Laplace operator △^∗ to solve the GS equation for Ψ, we encounter two issues: (1) the RHS of the GS equation involves 2 unknown functions P′(Ψ) and gg′(Ψ), in addition to the main unknown function Ψ(R,Z); (2) the boundary condition needs to be given: i.e., the values of Ψ on a rectangular computational boundry need to be provided.

To address issue (1):  function forms of P′(Ψ) and g′g(Ψ) are chosen by users, often as polymials in Ψ_N with coefficients that are assumed known (which are later obtained in the least-square method to make the resulting solution approximately match some measurements, as is discussed in the previous section).

To address issue (2): The value of Ψ on the rectangular boundary is updated using the Green’s function method (assume that external coil currents are given), as is discussed in the previous section.

Guess values of Ψ on both the inner gridpoints and boundary points. After this, we can invert △^∗ to update Ψ within the computational box, and use the Green’s function method to update values of Ψ on the boundary. Then iterate to converge.

Two iterations: one for Ψ values on the inner gridpoints, one over Ψ values on the computational boundary.

External PF coil currents enter the problem via the boundary condition: specifically via its contribution to Ψ on the boundary.

The above procedure is just a minor modification from the pure Green’s function algorithm discussed in Sec. 6.2, i.e., replacing the Green’s function method for the inner gridpoints with the finte-difference Laplace operator inverting.

For reference easy, the finite-difference schemes for the Laplace operator is listed below.

(82)

Method 1: discretize Eq. (83) as

where b = −−, a = , c = , d = e = . Method 2: discetize Eq. (84)

− + = −μ0RiJϕi,j,

which can be organized as

where b = − − a = + , c = −, d = e =

My numerical experiments indicate this scheme is more stable than than the first one (for the same grid size, and without the VDE stabilizer, the first scheme blow up while the latter scheme works well).

7 Semi-free bounday equilibrium problem

Semi-free boundary equilibrium problems refer to the cases where the LCFS shape and the value of Ψ on it are given and one is asked to solve the currents in the PF coils. This is similar to the fitting problem discussed above: we consider the control points on the LCFS as virtual flux loops that measure the poloidal flux.

8 Magnetic control

An equilibrium reconstruction code (e.g., rtefit) can be used to “observe” the inner structure of Ψ. PF coil currents can provide actions on Ψ. With the sensors and actuators, we can feedback controll the magnetic configuration (flux map).

For example, consider controlling the shape of LCFS (often called iso-flux control). We choose a target surface. To make sure that the target surface is a magnetic surface, the value of Ψ on it should be made a spatial constant. Dentoe this constant by Ψ_t. To make sure that there is no closed flux surface outside the target surface, we select the target surface as a surface that is tangential to the first wall (for limiter configuration) or there is a X-point on it (for divertor configuration). Then discretize the target surface by a series of discrete coordinates (R_i,Z_i) with i = 1,2,…M. The actions (coil current changes) are chosen to minimise the following residual:

(87)

where

(88)

(89)

the superscritp (old) indicates values before the control actions, Ψ_i^(old) is obtained by a real-time equilibrium reconstruction code.

This assumes that the plasma current density remains the same when the actions are imposed (i.e., no plasma response is included), so that only the coil current changes contribute to the Ψ change.

If we need to feedback control X points to derised positions, then the above sum can be extended to include the following:

9 Fixed boundary equilibrium and choices of coordinates

Plasma physists love to consider a kind of simplified problem: the fixed boundary equilibrium problem, where the shape of the boundary flux surface is given (the value of Ψ is a constant on this boundary). In dealing with the fixed boundary problem, the curvilinear coordinate system is useful. Specifically, the convenience is that the coordinates can be adjusted to make one of the coordinate surfaces coincide with the given boundary flux surface, so that the boundary condition becomes trivial.

More generally, curvilinear coordinates can be ajusted to make coordinate surfaces coincide with magnetic surfaces (dicussed later). This often simplifies analysis of wave and transport problems, especially when we properly choose the poloidal/toroidal coordinate to make magnetic field lines look like straight lines in terms of these coordinates.

Next section discusses the basic theory of curvilinear coordinates system[4].

In many studies of tokamak plasmas, one need construct a curvilinear coordinate system based on a given magnetic cofiguration in order to make the problem amenable to analytical methods or numerical methods. Specifically, many theories and numerical codes use the curvilinear coordinate systems that are constructed with one coordinate surface coinciding with magnetic surfaces. In these coordinate systems, we need to choose a poloidal coordinate 𝜃 and a toroidal coordinate ζ. As metioned above, a particular choice for 𝜃 and ζ is one that makes the magnetic field lines be straight lines in (𝜃,ζ) plane. These kinds of coordinates are often called magnetic coordinates. That is, “magnetic coordinates  are defined so they conform to the shape of the magnetic surfaces and trivialize the equations for the field lines.”

A further tuned magnetic coordinate system is the so-called field aligned (or filed-line following) coordinate system, in which changing one of the three coordinates with the other two fixed would correspond to following a magnetic field line. The field aligned coordinates are discussed in Sec. 17.

Next, let us discuss some general properties about coordinates transformation.

9.1 Safety factor

A magnetic field line on a closed magnetic surface travel a closed curve in a poloidal plane. For these kinds of field lines, we can define the safety factor q: the number of toroidal loops a magnetic field line travels when it makes one poloidal loop, i.e.

(91)

where △ϕ as the change of the toroidal angle when a magnetic field line travels a full poloidal loop. The safety factor can also be understood as the average pitch angle of a magnetic field line in (𝜃,ϕ) plane of a closed magnetic surface.

For open field line region (where a field line touches the wall before its poloidal projection can close itself), the “connection length” is often used to characterize the magnetic field.

Expression of safety factor in terms of magnetic field

The equation of magnetic field lines is given by

(92)

where dℓ_p is the line element along the direction of B_p on the poloidal plane. Equation (92) can be arranged in the form

(93)

which can be integrated over dℓ_p to give

(94)

where the line integration is along the poloidal magnetic field (the contour of Ψ on the poloidal plane). Using this, Eq. (91) is written

(95)

Expression of safety factor in terms of magnetic flux

The safety factor given by Eq. (95) is expressed in terms of the components of the magnetic field. The safety factor can also be expressed in terms of the magnetic flux. Define ΔΨ_p as the poloidal magnetic flux enclosed by two neighboring magnetic surface, then ΔΨ_p is given by

(96)

where Δx is the length of a line segment in the poloidal plane between the two magnetic surfaces, which is perpendicular to the first magnetic surface (so perpendicular to the B_p). Note that Δx, as well as R and B_p, generally depends on the poloidal location whereas ΔΨ_p is independent of the poloidal location.

Using Eq. (96), the poloidal magnetic field is written as

(97)

Substituting Eq. (97) into Eq. (95), we obtain

(98)

We know ΔΨ_p is a constant independent of the poloidal location, so ΔΨ_p can be taken outside the integration to give

(99)

It is ready to realise that the integral appearing in Eq. (99) is the toroidal magnetic flux enclosed by the two magnetic surfaces, ΔΨ_t. Using this, Eq. (99) is written as

(100)

Equation (100) indicates that the safety factor of a magnetic surface is equal to the differential of the toroidal magnetic flux with respect to the poloidal magnetic flux enclosed by the magnetic surface.

Rational surfaces vs. irrational surfaces

If the safety factor of a magnetic surface is a rational number, i.e., q = m∕n, where m and n are integers, then this magnetic surface is called a rational surface, otherwise an irrational surface. A field line on a rational surface with q = m∕n closes itself after it travels n poloidal loops. An example of a field line on a rational surface is shown in Fig. 5.

10 Curvilinear coordinate system

10.1 Coordinates transformation

In the Cartesian coordinates, a point is described by its coordinates (x,y,z), which, in the vector form, is written as

(101)

where r is the location vector of the point; , , and are the basis vectors of the Cartesian coordinates, which are constant, independent of spactial location. The transformation between the Cartesian coordinates system, (x,y,z), and a general coordinates system, (x₁,x₂,x₃), can be expressed as

(102)

For example, cylindrical coordinates (R,ϕ,Z) can be considered as a general coordinate systems, which are defined by  r = R cosϕ + R sinϕ + Z.

The transformation function in Eq. (102) can be written as

(103)

10.2 Jacobian

A useful quality characterizing coordinate transformation is the Jacobian determinant (or simply called Jacobian), which, for the transformation in Eq. (103), is defined by

(104)

which can also be written as

(105)

It is easy to prove that the Jacobian 𝒥 in Eq. (105) can also be written (the derivation is given in my notes on Jacobian)

(106)

Conventionally, the Jacobian of the transformation from the Cartesian coordinates to a particular coordinate system σ is called the Jacobian of σ, without explitly mentioning that this transformation is with respect to the Cartesian coordinates.

Using the defintion in Eq. (104), the Jacobian 𝒥 of the Cartesian coordinates can be calculated, yielding 1. Likewise, the Jacobian of the cylindrical coordinates (R,ϕ,Z) can be calculated as follows:

10.3 Orthogonality relation between two sets of basis vectors

In a curvilinear coordinate system (x₁,x₂,x₃), there are two kinds of basis vectors: ∇x_i and ∂r∕∂x_i, with i = 1,2,3. These two kinds of basis vectors satisfy the following orthogonality relation:

(107)

where δ_ij is the Kronical delta function. [Proof: Working in a Cartesian coordinate system (x,y,z) with the corresponding basis vectors denoted by (,,), then the left-hand side of Eq. (107) can be written as

[The cylindrical coordinate system (R,ϕ,Z) is an example of general coordinates. As an exercise, we can verify that the cylindrical coordinates have the property given in Eq. (107). In this case, x = x₁ cosx₂, y = x₁ sinx₂, z = x₃, where x₁ ≡ R, x₂ ≡ ϕ, x₃ ≡ Z.]

It can be proved that ∇xⁱ is a contravariant vector while ∂r∕∂xⁱ is a covariant vector (I do not prove this and do not bother with the meaning of these names, just using this as a naming scheme for easy reference).

The orthogonality relation in Eq. (107) is fundamental to the theory of general coordinates. The orthogonality relation allows one to write the covariant basis vectors in terms of contravariant basis vectors and vice versa. For example, the orthogonality relation tells that ∂r∕∂x₁ is orthogonal to ∇x₂ and ∇x₃, thus, ∂r∕∂x₁ can be written as

(110)

where A is a unknown variable to be determined. To determine A, dotting Eq. (110) by ∇x₁, and using the orthogonality relation again, we obtain

(111)

which gives

(113)

Similarly, we obtain

(114)

and

(115)

Equations (113)-(115) can be generally written

(116)

where (i,j,k) represents the cyclic order in the variables (x₁,x₂,x₃). Equation (116) expresses the covariant basis vectors in terms of the contravariant basis vectors. On the other hand, from Eq. (113)-(115), we obtain

(117)

which expresses the contravariant basis vectors in terms of the covariant basis vectors.

10.4 An example: (ψ,𝜃,ζ) coordinates

Suppose (ψ,𝜃,ζ) is an arbitrary general coordinate system. Following Einstein’s notation, contravariant basis vectors are denoted with upper indices as

(120)

In term of the contravairant basis vectors, A is written

(121)

where the components are easily obtained by taking scalar product with e_ψ,e_𝜃,ande_ζ, yielding A_ψ = A ⋅ e_ψ, A_𝜃 = A ⋅ e_𝜃, and A_ζ = A ⋅ e_ζ. Similarly, in term of the covariant basis vectors, A is written

(122)

where A^ψ = A ⋅ e^ψ, A^𝜃 = A ⋅ e^𝜃, and A^ζ = A ⋅ e^ζ.

Using the above notation, the relation in Eq. (116) is written as

(123)

(124)

(125)

where  𝒥 = [(∇ψ ×∇𝜃) ⋅∇ζ]⁻¹. Similarly, the relation in Eq. (117) is written as

(126)

(127)

(128)

10.5 Gradient and directional derivative in general coordinates (ψ,𝜃,ζ)

The gradient of a scalar function f(ψ,𝜃,ζ) is readily calculated from the chain rule,

(129)

Note that the gradient of a scalar function is in the covariant representation. The inverse form of this expression is obtained by dotting the above equation respectively by the three contravariant basis vectors, yielding

(130)

(131)

(132)

Using Eq. (129), the directional derivative in the direction of ∇ψ is written as

(133)

10.6 Divergence operator in general coordinates (ψ,𝜃,ζ)

To calculate the divergence of a vector, it is desired that the vector should be in the contravariant form because we can make use of the fact:

(134)

for any scalar quantities α and β. Therefore we write vector A as

(135)

where A^(ψ) = A ⋅∇ψ, A^(𝜃) = A ⋅∇𝜃, A^(ζ) = A ⋅∇ζ.   Then the divergence of A is readily calculated as

10.7 Laplacian operator in general coordinates (ψ,𝜃,ζ)

The Laplacian operator is defined by ∇² ≡∇⋅∇. Then ∇²f is written as (f is an arbitrary function)

10.8 Curl operator in general coordinates (ψ,𝜃,ζ)

To take the curl of a vector, it should be in the covariant representation since we can make use of the fact that ∇×∇α = 0. Thus the curl of A is written as

10.9 Metric tensor for general coordinate system

Consider a general coordinate system (ψ,𝜃,ζ). I define the metric tensor as the transformation matrix between the covariant basis vectors and the contravariant ones. Equations (116) and (117) express the relation between the two sets of basis vectors using cross product. Next, let us express the relation in matrix from. To obtain the metric matrix, we write the contrariant basis vectors in terms of the covariant ones, such as

(144)

Taking the scalar product respectively with ∇ψ, ∇𝜃, and ∇ζ, Eq. (144) is written as

(145)

(146)

(147)

Similarly, we write

(148)

Taking the scalar product with ∇ψ, ∇𝜃, and ∇ζ, respectively, the above becomes

(149)

(150)

(151)

The same situation applies for the ∇ζ basis vector,

(152)

Taking the scalar product with ∇ψ, ∇𝜃, and ∇ζ, respectively, the above equation becomes

(153)

(154)

(155)

Summarizing the above results in matrix form, we obtain

(156)

Similarly, to convert contravariant basis vector to covariant one, we write

(157)

Taking the scalar product respectively with ∇𝜃 ×∇ζ𝒥 , ∇ζ ×∇ψ𝒥 , and ∇ψ ×∇𝜃𝒥 , the above equation becomes

(158)

(159)

(160)

For the second contravariant basis vector

(161)

(162)

(163)

(164)

For the third contravariant basis vector

(165)

(166)

(167)

(168)

Summarizing these results, we obtain

(169)

where

M = ,

This matrix and the matrix in Eqs. (156) should be the inverse of each other. It is ready to prove this by directly calculating the product of the two matrix.

Special case: metric tensor for (ψ,𝜃,ϕ) coordinate system

Suppose that (ψ,𝜃,ϕ) are arbitrary general coordinates except that ϕ is the usual toroidal angle in cylindrical coordinates. Then ∇ϕ = 1∕R is perpendicular to both ∇ψ and ∇𝜃. Using this, Eq. (156) is simplified to

(170)

Similarly, Eq. (169) is simplified to

(171)

[Note that the matrix in Eqs. (170) and (171) should be the inverse of each other. The product of the two matrix,

(172)

can be calculated to give

,

where

By using the definition of the Jacobian in Eq. (106), it is easy to verify that A = 1, i.e.,

(173)

]

11 Covariant/contravariant representation of equilibrium magnetic field

The axisymmetric equilibrium magnetic field is given by Eq. (54), i.e.,

(174)

In a general coordinate system (ψ,𝜃,ϕ) (not necessarily magnetic surface coordinates), the above expression can be written as

(175)

where the subscripts denote the partial derivatives with the corresponding subscripts. Note that Eq. (175) is a mixed representation, which involves both covariant and contravariant basis vectors. Equation (175) can be converted to the contravariant form by using the metric tensor, giving

(176)

Similarly, Eq. (175) can also be transformed to the covariant form, giving

(177)

For the convenience of notation, define

(178)

then Eq. (177) is written as

(179)

12 Magnetic surface coordinates (ψ,𝜃,ϕ)

A coordinate system (ψ,𝜃,ϕ), where ϕ is the usual cylindrical toroidal angle, is called a magnetic surface coordinate system if Ψ is a function of only ψ, i.e., ∂Ψ∕∂𝜃 = 0 (we also have ∂Ψ∕∂ϕ = 0 since we are considering axially symmetrical case). In terms of (ψ,𝜃,ϕ) coordinates, the contravariant form of the magnetic field, Eq. (176), is written as

(180)

where Ψ′≡ dΨ∕dψ. The covariant form of the magnetic field, Eq. (177), is written as

(181)

12.1 Local safety factor

The local safety factor is defined by

(182)

which characterizes the local pitch angle of a magnetic field line in (𝜃,ϕ) plane of a magnetic surface. Substituting the contravariant representation of the magnetic field, Eq. (180), into the above equation, the local safety factor is written

(183)

Note that the expression in Eq. (183) depends on the Jacobian 𝒥 . This is because the definition of depends on the definition of 𝜃, which in turn depends on the the Jacobian 𝒥 .

In terms of , the contravariant form of the magnetic field, Eq. (180), is written

(184)

and the parallel differential operator B₀ ⋅∇ is written as

(185)

If happens to be independent of 𝜃 (i.e., field lines are straight in (𝜃,ϕ) plane), then the above operator becomes a constant coefficient differential oprator (after divided by 𝒥⁻¹). This simplification is useful because different poloidal harmonics are decoupled in this case. We will discuss this issue futher in Sec. 16.

12.2 Global safety factor

The global safety factor defined in Eq. (95) is actually the poloidal average of the local safety factor, i.e.,

Next, let us transform the 𝜃 integration in expression (187) to a curve integral in the poloidal plane. Using the relation dℓ_p and d𝜃 [Eq. (195)], expression (187) is further written

(189)

We note that expression (189) involves magnetic surface averaging, which is unknown before we know Ψ. Therefore iteration is usually needed in solving the fixed-q equilibrium (i.e., we guess the unknown Ψ, so that the magnetic surface averaging in expression (189) can be performed, yielding the values of g.)

Using B_p = |∇Ψ|∕R and B_ϕ = g∕R,  the absolute value of q in expression (188) is written

12.3 Relation between Jacobian and poloidal angle 𝜃

Given the definition of a magnetic surface coordinate system (ψ,𝜃,ϕ), the Jacobian of this system is fully determined. On the other hand, given the definition of ψ, ϕ, and the Jacobian, the definition of 𝜃 is fully determined (can have some trivial shifting freedoms). Next, let us discuss how to calculate 𝜃 in this case. In (ψ,𝜃,ϕ) coordinates, a line element is written

(192)

The line element that lies on a magnetic surface (i.e., dψ = 0) and in a poloidal plane (i.e., dϕ = 0) is then written

(194)

Using the fact that ∇ψ and ∇ϕ are orthogonal and ∇ϕ = ∕R, the above equation is written as

(195)

Given |𝒥∇ψ|, Eq. (195) can be integrated to determine the 𝜃 coordinate of points on a magnetic surface.

—–

(196)

12.4 Calculating poloidal angle

Once |𝒥∇ψ| is known, the value of 𝜃 of a point can be obtained by integrating expression (195), i.e.,

(197)

where the curve integration is along the contour Ψ = Ψ_j, x_ref,j is a reference point on the contour, where value of the poloidal angle is chosen as 𝜃_ref,j. The choice of the positive direction of 𝜃 is up to users. Depending on the positive direction chosen, the sign of the Jacobian of the constructed coordinates can have a sign difference from the 𝒥 appearing in Eq. (197). Denote the Jacobian of the constructed coordinates by 𝒥′, then

(198)

This sign can be determined after the radial coordinate and the positive direction of the poloidal angle are chosen. In GTAW code, I choose the positive direction of 𝜃 to be in anticlockwise direction when observers look along the direction of . To achieve this, the line integration in Eq. (197) should be along the anticlockwise direction. (I use the determination of the direction matrix, a well known method in graphic theory, to determine the direction from a given set of discrete points on a magnetic surface.)

Normalized poloidal angle

The span of 𝜃 defined by Eq. (197) is usually not 2π in one poloidal loop. This poloidal angle can be scaled by s(ψ) to define a new poloidal coordinate 𝜃, whose span is 2π in one poloidal loop, where s(ψ) is a magnetic surface function given by

(199)

Then 𝜃 is written as

[The reference points x_ref,j and the values of poloidal angle at these points can be chosen by users. One choice of the reference points x_ref,j are those points on the horizontal ray in the midplane that starts from the magnetic axis and points to the low filed side of the device and 𝜃_ref,j at these points is chosen as zero (this is my choice in the GTAW code). In the TEK code, the reference points are chosen at the high-field side of the midplane and 𝜃_ref,j = −π at the reference points.]

Jacobian for equal-arc-length poloidal angle

If the Jacobian 𝒥 is chosen to be of the following form

(202)

Then 𝜃_i,j in Eq. (200) is written

(203)

and the Jacobian of new coordinates (ψ,𝜃,ϕ), 𝒥_new, which is given by Eq. (201), now takes the form

(204)

Equation (203) indicates a set of poloidal points with equal arc intervals corresponds to a set of uniform 𝜃_i points. Therefore this choice of the Jacobian is called the equal-arc-length Jacobian. Note that Eq. (203) does not involve the radial coordinate ψ. Therefore the values of 𝜃 of points on any magnetic surface can be determined before the radial coordinate is chosen.

Jacobian for equal-volume poloidal angle (Hamada coordinate)

The volume element in (ψ,𝜃,ϕ) coordinates is given by dV = |𝒥|d𝜃dϕdψ. If we choose a Jacobian that is independent of 𝜃, then uniform 𝜃 grids will correspond to grids with uniform volume interval. In this case, 𝒥 is written as

(205)

where h(ψ) is a function independent of 𝜃. Then 𝜃_i,j in Eq. (200) is written

(206)

and the Jacobian of the new coordinates (ψ,𝜃,ϕ), 𝒥_new,  is given by Eq. (201), which now takes the following form:

(207)

Note that both 𝜃_i,j and 𝒥_new are independent of the function h(ψ) introduced in Eq. (205). (h(ψ) is eliminated by the normalization procedure specified in Sec. 12.4 due to the fact that h(ψ) is constant on a magnetic surface.) The equal-volume poloidal angle is also called Hamada poloidal angle.

The equal-volume poloidal angle is useful in achieving loading balance for parallel particle simulations. Assume that markers are loaded uniform in space and the poloidal angle is domain decomposed and assigned to different MPI process. Then the equal-volume poloidal angle can make marker number in each MPI process be equal to each other and thus work loading to each process be equal. (**check**If the domain decomposition is also applied to the radial direction, to achieve loading balance, then the radial coordinate ψ should be chosen in a way that makes ∮ (R∕∇ψ)dl_p be independent of ψ, so that 𝒥_new in Eq. (207) is constant in space.**)

Jacobian of Boozer’s poloidal angle

If the Jacobian 𝒥 is chosen to be of the Boozer form:

(208)

then the poloidal angle in Eq. (200) is written as

(209)

The final Jacobian is given by

(210)

The usefullness of Boozer poloidal angle will be further discussed in Sec. 16.8 after we introduce a gneralized toroidal angle.

Jacobian for PEST poloidal angle (straight-field-line poloidal angle)

If the Jacobian 𝒥 is chosen to be of the following form

then Eq. (183) implies that the local safety factor, (ψ,𝜃) = −g∕Ψ′, is a magnetic surface function, i.e., the magnetic field lines are straight in (ψ,𝜃) plane. Then the poloidal angle in Eq. (200) is written

(211)

The Jacobian 𝒥_new given by Eq. (201) now takes the form

(212)

Let us denote an arbitrary poloidal angle by 𝜃 and the above straight-field-line poloidal angle by 𝜗, then it is ready to find the following relation between 𝜃 and 𝜗:

(213)

where is the local safety factor corresponding to the arbitrary poloidal angle 𝜃, i.e., = B⋅∇ϕ∕(B⋅∇𝜃). [Proof: Using d𝜃 = dl_p, the poloidal angle 𝜗 given in Eq. (211) is written as

Note that Boozer poloidal angle is very close to the poloidal angel disccused here because the two Jacobians are very similar:

(215)

Comparison between different types of poloidal angle

All Jacobians introduced above can be written in a general form:

(216)

The choice of (i = 2,j = k = 0) gives the PEST coordinate, (i = j = 0,k = 2) give the Boozer coordinate, (i = j = 1,k = 0) gives the equal-arc coordinate, (i = j = k = 0) gives the Hammada coordinate.

Figure 6 compares the equal-arc-poloidal angle and the straight-line poloidal angle, which shows that the resolution of the straight-line poloidal angle is not good near the low-field-side midplane. Since ballooning modes take larger amplitude near the low-field-side midplane, better resolution is desired there. This is one reason that I often avoid using the straight-line poloidal angle in my numerical codes.

Fig. 6: The equal-arc poloidal angle (left) and the straight-line poloidal angle (right) for EAST equilibrium #38300@3.9s. The blue lines correspond to the equal-𝜃 lines, with uniform 𝜃 interval between neighbour lines. The red lines correspond to the magnetic surfaces, which start from = 0.01714 (the innermost magnetic surface) and end at = 0.9851 (the boundary magnetic surface), and are equally spaced in .

Verification of Jacobian

After the magnetic coordinates are constructed, we can evaluate the Jacobian 𝒥_new by using directly the definition of the Jacobian, i.e.,

(217)

which can be further written as

(218)

where the partial differential can be evaluated by using numerical differential schemes. The results obtained by this way should agree with results obtained from the analytical form of the Jacobian. This consistency check provide a verification for the correctness of the theory derivation and numerical implementation. In evaluating the Jacobian by using the analytical form, we may need to evaluate ∇ψ, which finally reduces to evaluating ∇Ψ. The value of |∇Ψ| is obtained numerically based on the numerical data of Ψ given in cylindrical coordinate grids. Then the cubic spline interpolating formula is used to obtain the value of |∇Ψ| at desired points. (𝒥_new calculated by the second method (i.e. using analytic form) is used in the GTAW code; the first methods are also implemented in the code for the benchmark purpose.) In the following sections, for notation ease, the Jacobiban of the constructed coordinate system will be denoted by 𝒥 , rather than 𝒥_new.

Radial coordinate

The radial coordinate ψ can be chosen to be various surface function, e.g., volume, poloidal or toroidal magnetic flux within a magnetic surface. The frequently used radial coordinates include Ψ, and , where Ψ is defined by

(219)

where Ψ₀ and Ψ_a are the values of Ψ at the magnetic axis and LCFS, respectively. Other choices of the radial coordinates: the toroidal magnetic flux and its square root, Ψ_t, and , where Ψ_t and Ψ_t are defined by

(220)

and

(221)

respectively, where Ψ_t(0) and Ψ_t(1) are the values of Ψ_t at the magnetic axis and LCFS, respectively.

If ψ = , then

(222)

The cylindrical coordinates (R,ϕ,Z) is a right-hand system, with the positive direction of Z pointing vertically up. In GTAW code, the positive direction of 𝜃 is chosen in the anticlockwise direction when observers look along the direction of . Then the definition 𝒥⁻¹ = ∇ψ ×∇𝜃 ⋅∇ϕ indicates that (1) 𝒥 is negative if ∇ψ points from the magnetic axis to LCFS; (2) 𝒥 is positive if ∇ψ points from the LCFS to the magnetic axis. This can be used to determine the sign of Jacobian after using the analytical formula to obtain the absolute value of Jacobian.

If ψ = Ψ or ψ = , ∇ψ points from the magnetic axis to LCFS.

The volume between magnetic surfaces can also be used as a radial coordinate. The differential volume element is written as

(223)

Integrating over the toroidal angle, we obtain

(224)

Further integrating over the poloidal angle, we obtain

(225)

i.e.,

(226)

In codes I wrote, I stick to using Ψ as the radial coordinate when doing computation, and transform to other radial coordinates when presenting the results if needed.

13 Constructing magnetic surface coordinate system from discrete Ψ(R,Z) data

Given an axisymmetric tokamak equilibrium in (R,ϕ,Z) coordinates (e.g., 2D data Ψ(R,Z) on a rectangular grids (R,Z) in G-file), we can construct a magnetic surface coordinates (ψ,𝜃,ϕ) by the following two steps. (1) Find out a series of magnetic surfaces on (R,Z) plane and select radial coordinates for each magnetic surface (e.g. the poloidal flux within each magnetic surface). (2) Specify the Jacobian or some property that we want the poloidal angle to have. Then calculate the poloidal angle of each point on each flux surface (on the ϕ = const plane) by using Eq. (195) (if the Jacobian is specified) or some method specified by us to achieve some property we prefer for the poloidal angle (if a Jacobian is not directly specified). Then we obtain the magnetic surface coordinates system (ψ,𝜃,ϕ).

13.1 Finding magnetic surfaces

Two-dimensional data Ψ(R,Z) on a rectangular grids (R,Z) is read from the G_EQDSK file (G-file) of EFIT code. Based on the 2D array data, I use 2D cubic spline interpolation to construct a interpolating function Ψ = Ψ(R,Z). To construct a magnetic surface coordinate system, I need to find the contours of Ψ, i.e., magnetic surfaces. The values of Ψ on the magnetic axis, Ψ₀, and the value of Ψ on the last closed flux surface (LCFS), Ψ_b, are given in G-file. Using these two values, I construct a 1D array “psival” with value of elements changing uniform from Ψ₀ to Ψ_b. Then I try to find the contours of Ψ with contour level value ranging from Ψ₀ to Ψ_b. This is done in the following way: construct a series of straight line (in the poloidal plane) that starts from the location of the magnetic axis and ends at one of the points on the LCFS. Combine the straight line equation, Z = Z(R), with the interpolating function Ψ(R,Z), we obtain a one variable function h = Ψ(R,Z(R)). Then finding the location where Ψ is equal to a specified value Ψ_i, is reduced to finding the root of the equation Ψ(R,Z(R)) − Ψ_i = 0. Since this is a one variable equation, the root can be easily found by using simple root finding scheme, such as bisection method (bisection method is used in GTAW code). After finding the roots for each value in the array “psival” on each straight lines, the process of finding the contours of Ψ is finished. The contours of Ψ found this way are plotted in Fig. 7.

In the above, we mentioned that the point of magnetic axis and points on the LCFS are needed to construct the straight lines. In G-file, points on LCFS are given explicitly in an array. The location of magnetic axis is also explicitly given in G-file. It is obvious that some of the straight lines Z = Z(R) that pass through the location of magnetic axis and points on the LCFS will have very large or even infinite slope. On these lines, finding the accurate root of the equation Ψ(R,Z(R)) −Ψ_i = 0 is difficult or even impossible. The way to avoid this situation is obvious: switch to use function R = R(Z) instead of Z = Z(R) when the slope of Z = Z(R) is large (the switch condition I used is |dZ∕dR| > 1).

In constructing the flux surface coordinate with desired Jacobian, we will need the absolute value of the gradient of Ψ, |∇Ψ|, on some specified spatial points. To achieve this, we need to construct a interpolating function for |∇Ψ|. The |∇Ψ| can be written as

(227)

By using the center difference scheme to evaluate the partial derivatives with respect to R and Z in the above equation (using one side difference scheme for the points on the rectangular boundary), we can obtain an 2D array for the value of |∇Ψ| on the rectangular (R,Z) grids. Using this 2D array, we can construct an interpolating function for ∇Ψ by using the cubic spline interpolation scheme.

13.2 Expression of metric elements of magnetic coordinates (ψ,𝜃,ϕ)

Metric elements of the (ψ,𝜃,ϕ) coordinates, e.g., ∇ψ ⋅∇𝜃, are often needed in practical calculations. Next, we express these metric elements in terms of the cylindrical coordinates (R,Z) and their partial derivatives with respect to ψ and 𝜃. Note that, in this case, the coordinate system is (ψ,𝜃,ϕ) while R and Z are functions of ψ and 𝜃, i.e.,

(228)

(229)

Then ∇R and ∇Z are written as

(230)

(231)

wehre R_ψ ≡ ∂R∕∂ψ, etc. Equations (230) and (231) can be solved to give

(232)

(233)

Using the above expressions, the Jacobian of (ψ,𝜃,ϕ) coordinates, 𝒥 , is written as

(235)

Using this, Expressions (232) and (233) are written as

(236)

and

(237)

Then the elements of the metric matrix are written as

(238)

(239)

and

(240)

Equations (238), (239), and (240) are the expressions of the metric elements in terms of R, R_ψ, R_𝜃, Z_ψ, and Z_𝜃. [Combining the above results, we obtain

(241)

Equation (240) is used in GTAW code. Using the above results, h^αβ = ∇α ⋅∇β are written as

(242)

(243)

(244)

As a side product of the above results, we can calculate the arc length in the poloidal plane along a constant ψ surface, dℓ_p, which is expressed as

14 Constructing model tokamak magnetic field

In some cases (e.g., turbulence simulation), model tokamak magnetic field, which is not an exact solution to the GS equation, is often used. In the model, the safety factor profile q(ψ), toroidal field function g(ψ), and the magnetic surface shape (R(ψ,𝜃),Z(ψ,𝜃)) are given. To use this model in a simulation, we need to calculate its poloidal magnetic field, which is determined by the poloidal magnetic flux function Ψ. To determine Ψ, we need to use Eq. (187), i.e.,

(246)

and re-organize the formula as

(247)

which can be integrated to obtain Ψ and thus the poloidal magnetic field B_p = ∇Ψ ×∇ϕ, where 𝒥 = [(∇ψ ×∇𝜃) ⋅∇ϕ]⁻¹. The toroidal magnetic field can be obtained from g(ψ) by B_ϕ = g∕R. In most papers, g(ψ) is chosen to be a constant.

A typical magnetic surface shape used in simulations is the Miller shape, which is given by

(248)

where r = ψ is the radial coordinate, R₀(r), δ(r), and κ_r(r) are the Shafronov shift, triangularity, and elongation profiles, which can be arbitrarily specified. For the special case of R₀ being a constant, δ(r) = 0, and κ(r) = 1, this shape reduces to a concentric-circular magnetic field. An example of calculating the poloidal field of the model field is given in Sec. 18. This kind of model field can be called theoretical physicists’ tokamak. Computational physicists often use this model for code benchmarking purpose. The famous DIII-D cyclone base case is an example, which was extensively used for benchmarking gyrokinetic simulation of ion temperature driven turbulence.

15 Magnetic surface averaging

15.1 Definition

The magnetic surface average of a physical quantity G(ψ,𝜃,ϕ) is defined by

(249)

where the volume integration is over a small volume between two adjacent flux surfaces with Ψ differing by △Ψ. [This formula (with finite ΔΨ) is used in TEK code to calculate the radial heat flux.]

The above 3D volume integration can also be written as a 2D surface integration. The differential volume element is given by d³V = |𝒥|dψd𝜃dϕ, where 𝒥 is the Jacobian of (ψ,𝜃,ϕ) coordinates. Using this, equation (249) is written as

(251)

where G₀(𝜃) is defined by the following Fourier expansion:

(252)

_________________________________________________________________________________________________________________________

“Zonal” and “mean” components

⟨G⟩ is sometimes called the “zonal” component of G if the radial wavelength of ⟨G⟩ is much smaller than the equilibrium scale length. If the radial wavelength of ⟨G⟩ is comparable to the equilibrium scale length, ⟨G⟩ is usually called “mean” component in tokamak literature. For example, mean flows are of system space scale and thus are easy to be observed in experiments. On the other hand, the “zonal” flow, which usually refers to the turbulence generated secondary flow, is of much smaller radial scale (the radial wavelength of zonal flow is of several Larmor radius) and thus is difficult to observe in experiments.

_________________________________________________________________________________________________________________________

Sometimes, we do not want the Jacobian to explicitly appear in the formula. This can be achieved by writing the differential volume element as

(253)

Using B_p = |∇Ψ|∕R, the volume element is further written as

(254)

Using this, the averaging defined in Eq. (249) is written as

(256)

(Equation (256) is used in the GTAW code to calculate the magnetic surface averaging.) Using Eq. (195) and B_p = |∇Ψ|∕R, equation (256) can also be written as

(257)

Using the expression of the volume element dτ = |𝒥|d𝜃dϕdψ, the volume within a magnetic surface is written

(258)

Using this, the differential of V with respect to ψ is written as

(259)

Using this, Eq. (257) is written as

⟨G⟩ = ∫ 02πG|𝒥|d𝜃

15.2 Flux Surface Functions—to be deleted

Next, examine the meaning of the following volume integral

(260)

where the volume V = V (ψ) is the volume within the magnetic surface labeled by ψ. Using ∇⋅ B = 0, the quantity D can be further written as

(261)

Note that 𝜃 is not a single-value function of the spacial points. In order to evaluate the integration in Eq. (261), we need to select one branch of 𝜃, which can be chosen to be 0 ≤ 𝜃 < 2π. Note that function 𝜃 = 𝜃(R,Z) is not continuous in the vicinity of the contour of 𝜃 = 0. Next, we want to use the Gauss’s theorem to convert the above volume integration to surface integration. Noting the discontinuity of the integrand 𝜃B in the vicinity of the contour of 𝜃 = 0, the volume should be cut along the contour, thus, generating two surfaces. Denote these two surfaces by S₁ and S₂, then equation (261) is written as

(263)

Using the expression of the volume element dτ = |𝒥|d𝜃dϕdψ, Ψ_p can be further written in terms of flux surface averaged quantities.

Similarly, the toroidal flux within a flux surface is written as

(265)

the poloidal current within a flux surface is written as

(266)

and toroidal current within a flux surface is written as

(267)

(**check**)The toroidal magnetic flux is written as

⇒ Ψt′ = g

⇒ = g

(269)

Next, calculate the derivative of the toroidal flux with respect to the poloidal flux.

(271)

By using the contravariant representation of current density (502), the poloidal current within a magnetic surface is written as

−∇ψ ×∇𝜃 − g′∇ϕ ×∇ψ,

The toroidal current is written as

(274)

Next, calculate another useful surface-averaged quantity,

16 Magnetic coordinates (ψ,𝜃,ζ) with general toroidal angle ζ

16.1 General toroidal angle ζ

In Sec. 12.1, we introduced the local safety factor (ψ,𝜃). Equation (183) indicates that if the Jacobian is chosen to be of the particular form 𝒥 = h(ψ)R², then the local safety factor is independent of 𝜃, i.e., magnetic line is straight in (𝜃,ϕ) plane. On the other hand, if we want to make field line straight in (𝜃,ϕ) plane, the Jacobian must be chosen to be of the specific form 𝒥 = h(ψ)R². We note that, as mentioned in Sec. 12.3, the poloidal angle is fully determined by the choice of the Jacobian. The specific choice of 𝒥 = α(ψ)R² is usually too restrictive for achieve a desired poloidal resolution (for example, the equal-arc poloidal angle can not be achieved by this choice of Jacobian). Is there any way that we can make the field line straight in a coordinate system at the same time ensure that the Jacobian can be freely adjusted to obtain desired poloidal angle? The answer is yes. The obvious way to achieve this is to define a new toroidal angle ζ that generalizes the usual toroidal angle ϕ. Define a new toroidal angle ζ by[10]

(277)

where δ = δ(ψ,𝜃) is a unknown function to be determined by the constraint of field line being straight in (𝜃,ζ) plane. Using Eq. (184), the new local safety factor in (ψ,𝜃,ζ) coordinates is written as

(279)

where c(ψ) can be an arbitrary function of ψ. A convenient choice for c(ψ) is c(ψ) = q, i.e., making the new local safety factor be equal to the original global safety factor, i.e., _new = q. In this case, equation (279) is written as

(280)

which, on a magnetic surface labed by ψ, can be integrated over 𝜃 to give

(281)

where 𝜃_ref is an starting poloidal angle arbitrarily chosen for the integration, and δ(ψ,𝜃_ref) is the constant of integration. In the following, both 𝜃_ref and δ(ψ,𝜃_ref) will be chosen to be zero. Then the above equations is written

(282)

Substituting the above expression into the definition of ζ (Eq. 277), we obtain

(283)

which is the formula for calculating the general toroidal angle. If 𝜃 is a straight-field line poloidal angle, then ζ in Eq. (283) reduces to the usual toroidal angle ϕ.

In summary, magnetic field line is straight in (𝜃,ζ) plane with slope being q if ζ is defined by Eq. (283). In this method, we make the field line straight by defining a new toroidal angle, instead of requiring the Jacobian to take particular forms. Thus, the freedom of choosing the form of the Jacobian is still available to be used later to choose a good poloidal angle coordinate. Note that the Jacobian of the new coordinates (ψ,𝜃,ζ) is equal to that of (ψ,𝜃,ϕ). [Proof:

[In numerical implementation, the term ∫ ₀^𝜃d𝜃 appearing in δ is computed by using

For later use, from Eq. (282), we obtain

16.2 Contravariant form of magnetic field in (ψ,𝜃,ζ) coordinates

Recall that the contravariant form of the magnetic field in (ψ,𝜃,ϕ) coordinates is given by Eq. (184), i.e.,

(287)

Next, let us derive the corresponding form in (ψ,𝜃,ζ) coordinates. Using the definition of ζ, equation (287) is written as

The expression of the magnetic field in Eq. (289) can be rewritten in terms of the flux function Ψ_p and Ψ_t discussed in Sec. 15.2. Equation (289) is

(290)

which, by using Eq. (264), i.e., ∇Ψ = ∇Ψ_p∕(2π), is rewritten as

(291)

which, by using Eq. (271), i.e., q = dΨ_t∕dΨ_p, is further written as

(292)

16.3 Relation between the partial derivatives in (ψ,𝜃,ϕ) and (ψ,𝜃,ζ) coordinates

Noting the simple fact that

(293)

where c is a constant, we conclude that

(294)

(since ϕ = ζ + q(ψ)δ(ψ,𝜃), where the part q(ψ)δ(ψ,𝜃) acts as a constant when we hold ψ and 𝜃 constant), i.e., the symmetry property with respect to the new toroidal angle ζ is identical with the one with respect to the old toroidal angle ϕ. On the other hand, generally,

(295)

and

(296)

In the special case that f is axisymmetric (i.e., f is independent of ϕ in (ψ,𝜃,ϕ) coordinates), then two sides of Eqs. (295) and (296) are equal to each other. Note that the partial derivatives ∂∕∂ψ and ∂∕∂𝜃 in Sec. 16.1 and 16.2 are taken in (ψ,𝜃,ϕ) coordinates. Because the quantities involved in Sec. 16.1 and 16.2 are axisymmetric, these partial derivatives are equal to their counterparts in (ψ,𝜃,ζ) coordinates.

16.4 Steps to construct a straight-line magnetic coordinate system

In Sec. 13, we have provided the steps to construct the magnetic surface coordinate system (ψ,𝜃,ϕ). Only one additional step is needed to construct the straight-line flux coordinate system (ψ,𝜃,ζ). The additional step is to calculate the generalized toroidal angle ζ according to Eq. (277), where δ is obtained from Eq. (282). Also note that the Jacobian of (ψ,𝜃,ζ) coordinates happens to be equal to that of (ψ,𝜃,ϕ) coordinates.

16.5 Form of operator B ⋅∇ in (ψ,𝜃,ζ) coordinates

The usefulness of the contravariant form [Eq. (289] of the magnetic field lies in that it allows a simple form of B ⋅∇ operator in a coordinate system. (The operator B₀ ⋅∇ is usually called magnetic differential operator.) In (ψ,𝜃,ζ) coordinate system, by using the contravariant form Eq. (289), the operator is written as

(298)

where h = h(ψ,𝜃,ζ) is some known function. Using Eq. (297), the magnetic differential equation is written as

(299)

Note that the coefficients before the two partial derivatives of the above equation are all independent of 𝜃 and ζ. This indicates that different Fourier harmonics in 𝜃 and ζ are decoupled. As a result of this fact, if f is Fourier expanded as

(300)

(note that, following the convention adopted in tokamak literature[10], the Fourier harmonics are chosen to be e^{i(m𝜃−nζ)}, instead of e^i(m𝜃+nζ)), and the right-hand side is expanded as

(301)

then Eq. (299) can be readily solved to give

(302)

The usefulness of the straight line magnetic coordinates (ψ,𝜃,ζ) lies in that, as mentioned previously, it makes the coefficients before the two partial derivatives both independent of 𝜃 and ζ, thus, allowing a simple solution to the magnetic differential equation.

16.6 Resonant surface of a perturbation

Equation (302) indicates that, for the differential equation (298), there is a resonant response to a perturbation e^{i(m𝜃−nζ)} on a magnetic surface with m − nq = 0. Therefore, the magnetic surface with q = m∕n is called the “resonant surface” for the perturbation e^{i(m𝜃−nζ)}.

The phase change of the perturbation e^{i(m𝜃−nζ)} along a magnetic field is given by mΔ𝜃 − nΔζ, which can be written as Δ𝜃(m−nq). Since m−nq = 0 on a resonant surface, this indicates that there is no phase change along a magnetic field line on a resonant surface, i.e., the parallel wavenumber k_∥ is zero on a resonant surface.

16.7 Helical angle used in tearing mode theory

Next, we discuss a special poloidal angle, which is useful in describling a perturnbation of single harmonic (m,n). This poloidal angle is defined by

(303)

where (m,n) are the mode numbers of the perturbation. The poloidal angle χ is often called helical angle and is special in that its definition is associated with a perturbation (the mode numbers of the perturbation appear in the definition) while the definition of the poloidal angles discussed previously only involve the equilibrium quantities.

The poloidal angle χ is designed to make 3D perturbations of the form ∼ f(ψ,m𝜃 −nζ) reduce to 2D perturbations, i.e.,

(304)

It is ready to verify that the Jacobian of coordinates (ψ,χ,ζ) is equal to that of coordinates (ψ,𝜃,ζ) [proof: (𝒥′)⁻¹ = ∇ψ ×∇χ ⋅∇ζ = ∇ψ ×∇(𝜃 − nζ∕m) ⋅∇ζ = ∇ψ ×∇𝜃 ⋅∇ζ = 𝒥⁻¹].

The component of B along ∇χ direction (i.e., the covariant component) is written

On the other hand, the component of B along ∇𝜃 direction is written

(307)

16.8 Covariant form of magnetic field in (ψ,𝜃,ζ) coordinate system

In the above, we have obtained the covariant form of the magnetic field in (ψ,𝜃,ϕ) coordinates (i.e., Eq. (181)). Next, we derive the corresponding form in (ψ,𝜃,ζ) coordinate. In order to do this, we need to express the ∇ϕ basis vector in terms of ∇ψ, ∇𝜃, and ∇ζ basis vectors. Using the definition of the generalized toroidal angle, we obtain

(309)

This expression can be further simplified by using equation (280) to eliminate ∂δ∕∂𝜃, which gives

(312)

with I(ψ) = h(ψ) − gq. The magnetic field expression in Eq. (312) frequently appears in tokamak literature[30]. In this form, the coefficients before both ∇𝜃 and ∇ζ depends on only the radial coordinate. In terms of I(ψ), the Jabobian can also be written as

(313)

16.9 Form of operator (B ×∇ψ∕B²) ⋅∇ in (ψ,𝜃,ζ) coordinates

In solving the MHD eigenmode equations in toroidal geometries, besides the B ⋅∇ operator, we will also encounter another surface operator (B×∇ψ∕B²) ⋅∇. Next, we derive the form of the this operator in (ψ,𝜃,ζ) coordinate system. Using the covariant form of the equilibrium magnetic field [Eq.  (311)], we obtain

(314)

Using this, the (B ×∇ψ∕B²) ⋅∇ operator is written as

Examining Eq. (316), we find that the coefficients before the two partial derivatives will be independent of 𝜃 and ζ if the Jacobian 𝒥 is chosen to be of the form 𝒥 = h(ψ)∕B², where h is some magnetic surface function. It is obvious that the independence of the coefficients on 𝜃 and ζ will be advantageous to some applications. The coordinate system (ψ,𝜃,ζ) with the particular choice of 𝒥 = h(ψ)∕B² is called the Boozer coordinates. The usefulness of the new toroidal angle ζ is highlighted in Boozer’s choice of the Jacobian, which makes both B ⋅∇ and (B ×∇ψ∕B²) ⋅∇ be a constant-coefficient differential operator. For other choices of the Jacobian, only the B⋅∇ operator is a constant-coefficient differential operator.

16.10 Radial differential operator

In solving the MHD eigenmode equations in toroidal geometry, we also need the radial differential operator ∇ψ ⋅∇. Next, we derive the form of the operator in (ψ,𝜃,ζ) coordinates. Using

∇f = ∇ψ + ∇𝜃 + ∇ζ,

the radial differential operator is written as

(318)

This formula is used in GTAW code.

17 Field-line-following coordinates

17.1 Definition of the field-line-following coordinates (ψ,𝜃,α)

In (ψ,𝜃,ζ) coordinates, a magnetic field line is straight in (𝜃,ζ) plane with slope being q. Then the equation for a magnetic field line is written as

(319)

where α is a constant in (𝜃,ζ) plane and can be used to label magnetic field lines on a magnetic surface. This motivates us to use α, i.e.,

(320)

to replace ζ. Then the magnetic field in Eq. (289) is written as

(321)

which is called the Clebsch form. The direction

(322)

is parallel (or anti-parallel) to the magnetic field direction. Due to this fact, (ψ,𝜃,α) coordinates are usually called “field-line-following coordinates” or “field-aligned coordinates” [2, 6].

Equation (321) implies that

(323)

and

(324)

i.e., both α and ψ are constant along a magnetic field line. Taking scalar product of Eq. (321) with ∇𝜃, we obtain

(325)

which is nonzero, i.e., only 𝜃 among (ψ,𝜃,α) is changing along a magnetic field line. (Here 𝒥 = (∇ψ ×∇𝜃 ⋅∇α)⁻¹ is the Jacobian of the coordinate system (ψ,𝜃,α), which happens to be equal to the Jacobian of (ψ,𝜃,ζ) coordinates.)

Using Eq. (321), the magnetic differential operator B ⋅∇ in the new coordinate system (ψ,𝜃,α) is written

(326)

which is just a partial derivative over 𝜃, as is expected, since only 𝜃 is changing along a magnetic field line.

By the way, note that (B ⋅∇α)∕(B ⋅∇𝜃) = 0, i.e., the magnetic field lines are straight with zero slope on (𝜃,α) plane.

Using Eqs. (320) and (283), α can be written as

It is widely believed that turbulence responsible for energy transport in tokamak plasmas usually has k_∥≪ k_⊥, where k_∥ and k_⊥ are the parallel and perpendicular wavenumbers, respectively. Due to this elongated structure along the parallel direction, less grids can be used in the parallel direction than that in the perpendicular direction in turbulence simulation. In this case, the field-aligned coordinates (ψ,𝜃,α) provide suitable coordinates to be used, where less gridpoints can be used for 𝜃 coordinate in simulations and even some ∂∕∂𝜃 derivatives can be neglected (high-n approximation), which simplifies the equations that need to be solved. There is another reason why almost all gyrokinetic codes use field-aligned coordinates: the  stability of numerical algorthims is improved when we use coarse grids in the parallel direction because the parallel Courant condition (for explicit schemes) Δt ≤ ΔL_∥∕v_∥ can be more easily satisfied (especially for the cases with kinetic electrons), where ΔL_∥ is the parallel grid spacing, which is larger when coarse grids are used in the parallel direction. This is also mentioned in Ref. [23] and it seems to be right from my experiences of testing several algorithm but a strict test is needed to verify this. This can also be understood in the following way: the coarse parallel grid automatically filters out physically irrelevant but numerically problematic high-k_∥ modes, permitting much longer time steps for explicit time stepping, in both particle and fluid codes[12].

17.2 Some discussions

The fact B ⋅∇α = 0 implies that α is constant along a magnetic field line. At first glance, a magnetic line on an irrational surface seems to sample all the points on the surface. This seems to indicate that α is a flux surface label for irrational surface. However, α must be a non-flux-surface-function so that it can provide a suitable toroidal coordinate. I had once been confused by this conflict for a long time. The key point to resolve this confusion is to realize that it is wrong to say there is only one magnetic line on an irrational surface, i.e. it is wrong to say a magnetic line on an irrational surface samples all the points on the surface. There are still infinite number of magnetic field lines that can not be connected with each other on an irrational surface. Then the fact B ⋅∇α = 0 does not imply that α must be the same on these different magnetic field lines. In fact, although B ⋅∇α = 0, the gradient of α on a flux-surface along the perpendicular (to B) direction is nonzero, i.e., B ×∇Ψ ⋅∇α≠0. [Proof:

In (ψ,𝜃,ϕ) coordinates, ∇ϕ is perpendicular to ∇ψ. However, in field-line-following coordinates (ψ,𝜃,α), ∇α is not perpendicular to ∇ψ. Therefore ∇α is not along the binormal direction B ×∇ψ.

17.3 Expression of ∇α, ∇ψ ⋅∇α, and ∇α ⋅∇𝜃

The generalized toroidal angle α is defined by Eq. (327), i.e., α = ϕ −δ, where δ = ∫ ₀^𝜃 d𝜃. The gradient of α is then written as

(330)

(331)

Another method of computing the above quantities: Using Eqs. (236) and (237), expression (329) is written as

17.4 Field-aligned coordinates in GEM[7] and GENE[13] codes

In GEM[7] and GENE[13] codes, the field-aligned coordinates (x,y,z) are defined by

(335)

(336)

(337)

where r is an arbitrary flux surface label with length dimension, which is often chosen in GEM to be the minor radius of a magnetic surface in the midplane. Here r₀ and R₀ are constant quantities of length dimension, r₀ is the minor radius of a reference magnetic surface (usually corresponding to the center of the radial simulation box), R₀ is the major radius of the magnetic axis, q₀ is the safety factor value on the r = r₀ surface. The constant length q₀R₀ introduced in the definition of z is to make z approximately correspond to the length along the field line in the large-aspect ratio limit. The constant length r₀∕q₀ introduced in the definition of y is to make y corresponds to the arc-length in the poloidal plane traced by a field line when its usual toroidal angle increment Δϕ is α. This explanation makes y look like a poloidal coordinate whereas y is actually a toroidal coordinate.

Next, let us calculate the wave number along the y direction, k_y, for a mode with toroidal wave number n. The wavelength along the y direction, λ_y, is given by

(338)

Then the wavenumber k_y is written as

(339)

On the other hand, the poloidal wavenumber k_𝜃 is given by

(340)

where m is the poloidal mode number of the mode in (ψ,𝜃,ϕ) coordinates. If the mode has the property k_∥≈ 0, i.e., m ≈ nq₀, then k_𝜃 is equal to the k_y. The motivation of introducing the constant length r₀∕q₀ in the definition of y is to make k_y ≈ k_𝜃 for a mode with k_∥≈ 0.

Some authors call k_y or k_𝜃 by the name “binormal wavenumber”, which is not an appropriate name in my opinion. Some authors call y the binormal direction, which is also an inappropriate name since neither ∇y nor ∂r∕∂y is along the binormal direction B₀ ×∇ψ.

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

(341)

(342)

(343)

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

(344)

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

(345)

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

(348)

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

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

(349)

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

(350)

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

(351)

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

(352)

However, generally there is no periodic condition along 𝜃,

(353)

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

(354)

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

(355)

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:

(356)

or equivalently

(357)

Numerical implementation of periodic condition (356)

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

(358)

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.