Notes on tokamak equilibrium

Youjun Hu1
Institute of Plasma Physics, Chinese Academy of Sciences
Email: yjhu@ipp.cas.cn

Note: This document has been written using GNU TE Xmacs [?].
Abstract. This document discusses the free boundary equilibrium fitting problem, fixed-boundary equilibrium problem, and magnetic coordinates.

 

 1 Axisymmetric magnetic field
 2 Non-axisymmetric magnetic perturbations
 3 Plasma current density in terms of Ψ and g
 4 Constraint of force-balance on magnetic field
 5 Free-boundary fitting problem
 6 Magnetic control (to be continued)
 7 Circuit equation
 8 Coordinates based on magnetic field structure
 9 Constructing magnetic surface coordinate system from discrete Ψ(R,Z) data
 10 Magnetic surface averaging
 11 Flux Surface Functions
 12 Magnetic flux diffusion equation
 13 Fixed boundary equilibrium problem
 14 Magnetic coordinates with general toroidal angle
 15 Field-line-following coordinates

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

B = ∇ × A,
(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 constraint.) In cylindrical coordinates (R,ϕ,Z), the above expression is written

    (             )    (                  )     (           )
B =   1∂AZ- − ∂A-ϕ Rˆ +  1-∂(RA-ϕ)− -1∂AR-  ˆZ +  ∂AR- − ∂AZ-  ˆϕ.
      R ∂ϕ    ∂Z         R   ∂R     R  ∂ϕ         ∂Z    ∂R
(2)

We consider axisymmetric magnetic field, which means that, when expressed in the cylindrical coordinate system (R,ϕ,Z), the components of B, namely BR, BZ, 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

     ∂A       1∂(RA  )    (∂A     ∂A  )
B = −---ϕRˆ+ -------ϕ-ˆZ +  ---R − ---Z  ˆϕ.
      ∂Z     R   ∂R         ∂Z     ∂R
(3)

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

1.1 Poloidal magnetic field

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

Ψ (R, Z) ≡ RA ϕ(R,Z ).
(4)

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

       1∂ Ψ
BR = −R-∂Z-,
(5)

      1-∂Ψ-
BZ =  R ∂R.
(6)

(Note that it is the property of being axisymmetric and divergence-free that enables us to express the two components of B, namely BR and BZ, 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:

           (            )
             ∂Ψ-ˆ   ∂Ψ-ˆ
B ⋅∇ Ψ = B ⋅  ∂RR +  ∂ZZ
          1 ∂Ψ ∂Ψ   1 ∂Ψ ∂Ψ
      = − R-∂Z-∂R-+ R-∂R-∂Z-
      = 0.                                         (7)
]

(We note that Ψ is related to poloidal magnetic flux, as we will discuss in Sec. 1.7.)

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

Bp = BRRˆ +BZ ˆZ
       1∂ Ψ     1∂Ψ
   = −R-∂Z-Rˆ+ R-∂R-ˆZ
     1
   = R-∇Ψ × ˆϕ
   = ∇Ψ × ∇ ϕ                                   (8)

1.2 Toroidal magnetic field

Next, let’s examine the toroidal component Bϕ. Equation (3) indicates that Bϕ involves both AR and AZ. 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

            g
Bϕ = Bϕϕˆ=  Rˆϕ = g∇ ϕ.
(9)

1.3 General form of axisymmetric magnetic field

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

B = Bp + B ϕ
  = ∇ Ψ × ∇ϕ + g∇ϕ,                           (10)
which is a general expression of axisymmetric magnetic field. Expression (10) is a famous formula in tokamak physics.

1.4 Gauge freedom of Ψ

Let us discuss the gauge freedom of A in the axisymmetric case. Magnetic field remains the same under the following gauge transformation of the vector potential:

Anew = A +∇f,
(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

      ∂f     ∂f    1 ∂f
∇f = ---Rˆ+ ---ˆZ + ----ϕˆ.
     ∂R     ∂Z     R ∂ϕ
(12)

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

 new       -1∂f-
Aϕ  = A ϕ + R ∂ϕ
           C-
    = A ϕ + R ,                             (13)
(We see that the requirement of being axial symmetry greatly reduces gauge freedom of Aϕ.) Multiplying Eq. (13) with R, we obtain the corresponding gauge transformation for Ψ,
Ψ new = Ψ + C,
(14)

which indicates Ψ is similar to that of electrostatic 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 the equilibrium reconstruction problem.

1.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

dΨ = 0,
(15)

i.e.,

∂Ψ-dR + ∂Ψ-dZ = 0.
∂R      ∂Z
(16)

   1 ∂Ψ      1∂Ψ
⇒  ----dR + -----dZ = 0.
   R∂R      R ∂Z
(17)

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

BZdR  − BRdZ = 0,
(18)

i.e.,

dZ-=  BZ-,
dR    BR
(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).


pict

 

Fig. 1: Blue lines are contours of Ψ. The black line is the machine wall.

Points where the poloidal field is zero (i.e., Ψ=0) are called magnetic null points. There are two types of 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

         ∂2Ψ ∂2Ψ   (  ∂2Ψ )2
S(R,Z ) = ∂R2-∂Z2-−  ∂R∂Z--  ,
(20)

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

1.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. Some Ψ contours intersect the wall before they can form closed curves. These flux surfaces are called “open”. Otherwsie, they are called closed flux surfaces.

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 (the area between the X-point and the material divertor., i.e, the region between divertor legs that is unconnected to the plasma).]

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


pict

Fig. 2: Closed magnetic surfaces (blue) and various toroidal ribbons used to define the poloidal magnetic flux discussed in Sec. 1.7. The magnetic flux through the toroidal ribbons S2 and S3 is equal to each other. Also the magnetic flux through S4 and S5 is equal to each other; the magnetic flux through S0 and S1 is equal to each other.

The innermost magnetic surface reduces to a curve, which is called magnetic axis (in Fig. 2, Ψ0 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)). For closed flux surfaces, since Bp = Ψ ×∇ϕ, the condition ΨLCFS Ψaxis > 0 means Bp points in the anticlockwise direction (viewed along ϕ direction), and ΨLCFS Ψaxis < 0 means Bp points in the clockwise direction.

1.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.


pict

Fig. 3: The poloidal magnetic flux Ψp(12) between the two magnetic surfaces Ψ1 and Ψ2 is given by Ψp(12) = 2π(Ψ2 Ψ1), where the surface orientation for the poloidal flux definition is chosen to be +ˆZ.

In Fig. 3, there are two magnetic surfaces labeled, respectively, by Ψ = Ψ1 and Ψ = Ψ2. The poloidal magnetic flux through any toroidal ribbons between the two magnetic surfaces is equal to each other (due to the magnetic Gauss theorem). Denote this poloidal magnetic flux by Ψp(12).  Next, we calculate this flux. To make the calculation easy, we select a plane perpendicular to the Z axis, as is shown by the dash line in Fig. (3). In this case, only BZ contribute to the poloidal magnetic flux: (the positive direction of the plane is chosen to be ˆZ)

  (12)  ∫ R2
Ψ p  =     Bz (R, Z)2πRdR
       ∫R1R2
     =      1∂-Ψ2πRdR
        R1  R∂R
          ∫ R2 ∂Ψ
     = 2π     ∂R-dR
           R1
     = 2π[Ψ2 − Ψ1].                              (21)
Equation (21) indicates that the difference of Ψ between two magnetic surfaces is equal to the poloidal flux divided by 2π.

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

Ψp = 2π[Ψ (R, Z)− Ψ(R = 0,Z)],
(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 the equilibrium reconstruction problem. Then this flux is written as

Ψp = 2πΨ (R,Z).
(23)

Due to this relation, Ψ is often called the poloidal magnetic flux per radian (SI unit: web/rad), or simply “the poloidal flux”. This relation allows the poloidal flux measurements to be used to constrain the GS equation in the equilibrium reconstrunction (discussed in Sec. 5). Note that the positive normal direction of the surface (where the magnetic flux Ψp is defined) is chosen in the +Z direction. This sign choice along with the 2π factor appearing in Eq. (23) are often confusing people when they try to relate the experimentally measured flux with the Ψ appearing in the GS equation. The 2π factor also appears when we relate the Green function of Ψ to the mutual inductance, i.e., M = 2πG (discussed later).

A further compliction is that when one talks about  “the poloidal magnetic flux of a magnetic surface”, there are two possibilities of the defintion. One of them is defined relative to the R = 0 (as is discussed above), and another is relative to the magnetic axis. The first definition is the flux through the central hole of the magnetic surface, i.e., the poloidal flux through S0 in Fig. 2. In this case, as is discussed above, the poloidal magnetic flux is related to Ψ by

Ψp = 2π(Ψ − Ψ(R = 0,Z )) = 2πΨ,
(24)

where the postive direction is +ˆZ. The is the more often used defintion, and we will stick to this in the equilibrium reconstruction problem.

The second definition is the flux enclosed by the closed flux surface, i.e., the flux through the toroidal ribbon S2. Denote this flux by Ψp(axis), then

Ψ(axis) = 2π(Ψaxis − Ψ),
 p
(25)

where Ψaxis is the value of Ψ at the magnetic axis, the S2 orientation is in the clockwise direction when an observer looks along the direction of ˆϕ.

1.8 Measuring poloidal magnetic flux in experiments

By measuring the voltage around a toroidal wire loop, we can obtain the time derivative of the poloidal flux and, after integrating over time, the flux itself.

Suppose that there is a toroidal wire loop 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, hence this flux is called the poloidal flux). Then Faraday’s law

∮          ∫  ∂B
   E⋅dl = −   ∂t-⋅dS,
(26)

where the direction of the loop integration and the direction of dS is related by the right-hand rule, is written as

𝜀 = − dΨp-,
      dt
(27)

where 𝜀 = E dl is often called electromotive force (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𝜀. Then Eq. (27) is written as

V = − NdΨp-.
        dt
(28)

Integrating the above equation over time, we obtain

                       1 ∫ t
Ψp(R,Z,t) = Ψp(R,Z,0)− --   Vdt.
                       N  0
(29)

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 (29) tells us how to calcualte Ψp from the measured loop voltage V . Then Ψ is obtained by Ψ = Ψp(2π).

There are usually many flux loops (e.g. 35 on EAST[?]) at different locations in the poloidal plane (see Fig. 4). They are outside of the plasma region and thus are “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. 5.


pict

Fig. 4: Flux loops (circles) and magnetic probes (blue arows) on EAST tokamak. Arrows indicate the positive normal direction of the probes.

1.9 Safety factor

A magnetic field line on a closed magnetic surface travel a closed curve in the poloidal plane. For these 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.

q ≡ △ϕ-,
    2π
(30)

where ϕ as the change of the toroidal angle when a magnetic field line travels a full poloidal loop.

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

Rdϕ   B
----= --ϕ,
dℓp   Bp
(31)

where dℓp is the line element along the direction of Bp on the poloidal plane. Equation (31) can be arranged in the form

     1-B-ϕ
dϕ = R Bp dℓp,
(32)

which can be integrated over dℓp to give

     ∮  1 Bϕ
△ϕ =    R-B-dℓp,
           p
(33)

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

     1 ∮ 1 B
q = ---  ---ϕ-dℓp.
    2π   R Bp
(34)

Expression of safety factor in terms of magnetic flux

The safety factor given by Eq. (34) 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

Δ Ψ = 2πR ΔxB
   p          p
(35)

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 Bp). Note that Δx, as well as R and Bp, generally depends on the poloidal location whereas ΔΨp is independent of the poloidal location.

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

      1  Δ Ψp
Bp = 2πR-Δx--.
(36)

Substituting Eq. (36) into Eq. (34), we obtain

      ∮               ∮                 ∮
q = 1--  1-Bϕ-dℓp = -1-  1-2πR-ΔxB-ϕdℓp =   ΔxB-ϕdℓp.
    2π   R Bp      2π   R    ΔΨp           Δ Ψp
(37)

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

        ∮
   --1-
q = Δ Ψp  ΔxB ϕdℓp
(38)

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

    Δ Φ
q = ----
    ΔΨp
(39)

Equation (39) 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.


pict pict

Fig. 5: Left: A magnetic field line (blue) on a rational surface with q = 2.1 = 2110 (magnetic field is from EAST discharge #59954@3.1s). This field line closes itself after traveling 21 toroidal loops (meanwhile, it travels 10 poloidal loops). Right: The intersecting points of the magnetic field line with the ϕ = 0 plane when it is traveling toroidally. The sequence of the intersecting points is indicated by the number labels. The 22nd intersecting point coincides with the 1st point and then the intersecting points repeat themselves.

 

 

2 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 AR and AZ) will appear in the expression of the poloidal magnetic field. Let us re-examine Eq. (2) for a general magnetic perturbation:

     (1 ∂δA     ∂δA )     ( 1 ∂(R δA )   1 ∂δA  )
δB  =  ------Z−  ---ϕ- ˆR +   -------ϕ--− -----R- Zˆ
     (R  ∂ϕ      ∂Z)        R   ∂R      R  ∂ϕ
   +  ∂δAR- − ∂δAZ- ϕˆ.                                       (40)
       ∂Z      ∂R
When studying tearing modes and electromagnetic turbulence, most authors narrow the possible perturbations by setting δAR = δAZ = 0, i.e.,
        1∂ δΨ
δBR = − -----,
        R ∂Z
(41)

δBZ  = 1-∂δΨ,
       R ∂R
(42)

δBϕ = 0.
(43)

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

δB = ∇ δΨ × ∇ϕ.
(44)

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

B = ∇ (Ψ + δΨ )× ∇ϕ + g∇ϕ.
(45)

**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).

3 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):

μ0J = ∇ × B                (            )
       ∂B-ϕˆ   1-∂(RB-ϕ)ˆ    ∂BR-   ∂BZ-  ˆ
   = −  ∂Z R + R   ∂R   Z+    ∂Z −  ∂R   ϕ,               (46)
where μ0 is vacuum magnetic permeability.

3.1 Poloidal current density

Use Eq. (46) and the definition g RBϕ, the poloidal components of the current density, JZ and JR, can be written as

μ0JR = −-1∂g-,
        R ∂Z
(47)

and

μ0JZ =  1-∂g,
        R∂R
(48)

respectively.

3.2 Toroidal current density

Ampere’s law (46) indicates the toroidal current density Jϕ is given by

      ∂BR    ∂BZ
μ0Jϕ =-∂Z- − -∂R-
        1 ∂2Ψ    ∂ ( 1 ∂Ψ)
    = − ----2-− ---  ----- .                      (49)
        R ∂Z    ∂R   R ∂R
Define by
                (      )
△ ∗ ≡ ∂2--+ R-∂-  1--∂- ,
      ∂Z2    ∂R   R ∂R
(50)

then Eq. (49) is written as

Jϕ = − -1-△ ∗Ψ.
       μ0R
(51)

(Many authors incorrectly refer to Eq. (51) as the Grad-Shafranov (GS) equation. Eq. (51) is just Ampere’s law, which has nothing to do with the force-balance. Only after we express Jϕ in terms of the plasma pressure, can Eq. (51) be called the GS equation, as is discussed Sec. 4.3.)

The operator is different from the Laplacian operator △≡∇⋅∇. In terms of the nabla operator , the operator is writen as

△∗Ψ = R2∇ ⋅ 1-∇ Ψ.
            R2
(52)

4 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

  (           )
ρ  ∂u-+ u ⋅∇u   = ρqE + J × B − ∇ ⋅ℙ
    ∂t
(53)

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 ??) and the plasma pressure is isotropic, then the steady state momentum equation (force balance equation) is written as

J × B = ∇P,
(54)

where P is the scalar plasma pressure.

Is the force balance (54) always satisfied in a real toakamak discharge? To answer this question, we need to go back to the original momentum equation (53). 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 < Cs, where Cs is the the sound speed), the leading order of the momentum equation is the force balance (54).[?]. I.e. the inertial effect can be neglected (plasma mass is approximately zero).

4.1 Parallel force balance

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

0 = B ⋅∇P,
(55)

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. ??) 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 = dP-
dΨB ⋅∇Ψ = 0,

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

4.2 Toroidal force balance

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

JZBR  − JRBZ = -1∂P-.
               R  ∂ϕ
(56)

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

JZBR − JRBZ  = 0
(57)

Using the expressions of the poloidal current density (47) and (48) in the force balance equation (57) yields

∂g      ∂g
---BR + ---BZ = 0,
∂R      ∂Z
(58)

which can be further written

B ⋅∇g = 0.
(59)

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. ??.)

4.3 Force balance along the major radius

Consider the force balance in Rˆ direction. The ˆR component of Eq. (54) is written

JϕBZ − JZB ϕ = ∂P
              ∂R
(60)

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

Jϕ 1-∂Ψ-− -1--∂g-g-= ∂P-.
   R ∂R   μ0R ∂R R   ∂R
(61)

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

Jϕ 1-∂Ψ-−-1--dg-∂Ψ-g-= dP-∂Ψ-,
  R ∂R   μ0R dΨ ∂R R   dΨ ∂R
(62)

which can be simplified to

J = R dP-+ --1--dg-g,
 ϕ    dΨ   μ0R dΨ
(63)

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

            dP    dg
△∗Ψ = − μ0R2---− ---g,
            dΨ   dΨ
(64)

i.e.,

           (      )
-∂2Ψ- +R -∂-  1-∂Ψ-  = − μ R2 dP-− dgg.
∂Z2     ∂R   R ∂R       0   dΨ   dΨ
(65)

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

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

JRBϕ JϕBR = ∂∂PZ-
⇒−-∂g-
∂Z1-
R-g
R -1
RΨ1-
R∂Ψ-
∂Z = μ0dP-
dΨ∂Ψ-
∂Z
⇒−-dg
dΨ∂Ψ-
∂Z1-
Rg-
R 1-
RΨ-1
R∂Ψ-
∂Z = μ0dP
dΨ∂Ψ-
∂Z
⇒−-dg
dΨ1-
R-g
R -1
RΨ1-
R = μ0dP
dΨ
⇒△Ψ = μ0R2dP
dΨ dg
dΨ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, ˆ
ϕ, ˆ
R, and B direction) and thus it must be satisfied in all the directions.]

4.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.,

B = ∇Ψ × ∇ ϕ+  g∇ϕ ,
    ◟--◝◜--◞   ◟◝◜◞
     poloidal   toroidal
(66)

For the above axisymmetric magnetic field to be consistent with the force balance equation (54), 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 (65) 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 (65) 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 ??.

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, jis 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 implicitly involves Ψ). This indicates that iterations are needed when numerically solving the GS equation.

 

5 Free-boundary fitting problem

A routine in tokamak operation is to reconstruct magnetic field under the constraints of MHD force balance and expermental 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 on this code.

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

  ∗                 ∑
△  Ψ = − μ0RJ ϕ − μ0R  Iiδ(R − Ri)δ(Z − Zi),               (67)
                     i
where Jϕ is the plasma toroidal current density,  Ii is a toroidal curent filament located at (R = Ri,Z = Zi), δ() is Dirac’s delta function. The solution to the above equation is given by
          ∫                              ∑
Ψ (R,Z) =   G (R′,Z ′,R, Z)Jϕ(R′,Z ′)dR ′dZ′ +   G (Ri,Zi,R, Z)Ii,
           Ω                              i
(68)

where G is the fundamental solution given by

  G(R ′,Z′,R,Z)
  μ0R ∫ 2π              R′cosϕ′dϕ′
= -4π-    ∘-(Z-−-Z′)2-+-(R-−-R-′cosϕ′)2 +-(R′sin-ϕ′)2,           (69)
       0
which is often called Green’s function in this context (which is obtained by using a formula similar to the Biot-Savart Law, see ??).

If all the toroidal currents (plasma current + coil/vessel curents) are known, then Ψ can be calculated by using Eq. (68). Unfortunately, Jϕ is usually unknown. Meanwhile the force-balance indicates that Jϕ can be expressed in terms of Ψ via Eq. (63), i.e.,

      dP     1  dg
Jϕ = R dΨ-+ μ-R-gdΨ-.
             0
(70)

Substituting this into Eq. (68) 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 Ii, are known ). We can see that this is a Picard iteration. Will the iteration converge? We do not know for sure. Numerical experiments indicate it does in some cases. Let us discuss some specific cases. We assume that the 1D functions, P(Ψ) and g(Ψ), can be modeled as (EFIT’s model[?]):

        { ∑NP −1
dP-(Ψ-)=     j=0  αj(xj − xNp )forx ≤ 1
  dΨ      0    for x > 1
(71)

 dg   { ∑NF −1βj(xj − xNF )forx ≤ 1
gdΨ-=   0 j=0 forx > 1
(72)

where

    Ψ − Ψ
x ≡ -----A-,
    ΨB − ΨA
(73)

ΨA and ΨB are values of Ψ at the magnetic axis and LCFS, the coefficients αj and βj are to be determined, Np and NF are integers chosen by users. Expressions (71) and (72) guarantee that dP∕dΨ and gdg∕dΨ are zero at and outside the LCFS, and thus no plasma current there. Note that expressions (71) and (72) are nonlinear functions of Ψ even for Np = NF = 1. This is because the unknowns, ΨA and ΨB, appear in the denominator of expression (73). Another nonlinearity is related to the fact that the unkonw Ψ determines where the LCFS is and thus determines the region where the current desnsity is set to zero. This also gives rise to the name “free-boundary” for this kind of problems.

Choosing an initial guess of Ψ(R,Z) on gridpoints, we can get the values of αj, βj, and Ii 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. 5.1.)

After the coefficients αj and βj are obtained, Jϕ can be updated by using Eqs. (70)-(72). Then, we use the latest Jϕ and Ii in Eq. (68) to re-compute the values of Ψ on gridpoints by using Eq. (68). The procedure is repeted 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. 5.4. The values of Ψ on the bounary still have to be obtained by using Eq. (68).)

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

5.1 Design matrix of the least square problem

For x 1, plugging expressions (71) and (72) into (70), we obtain

    ′  ′    ′N∑p−1    j   Np    -1--NF∑− 1   j    NF
Jϕ(R ,Z ) = R     αj(x − x  )+  μ0R′     βj(x  − x  ),          (74)
              j=0                    j=0
For notation ease, define the following basis functions:
          {
      ′    R ′(xj − xNp)         for  0 ≤ j < Np
bj(x,R ) =  (xj−Np − xNF )∕(μ0R ′) for  Np ≤ j < Np + NF
(75)

and the corresponding expansion coefficients

    {
     αj   for  0 ≤ j < Np
cj =  βj−Np   for  Np ≤ j < Np + NF
(76)

then Eq. (74) is written as

    Nb∑−1        ′
Jϕ =    cjbj(x,R ),
     j=0
(77)

where Nb = NP + NF. Plugging expression (77) into Eq. (68), we obtain

         Nb∑−1  ∫      ′ ′           ′   ′  ′
Ψ(R,Z ) =    cj Ω G(R ,Z ,R,Z )bj(x,R )dR dZ
          j=0
         Nc∑−1    (c)  (c)
       +     G (Rj ,Zj ,R,Z )Ij.                           (78)
          j=0
Collect all the free parameters as a column vector: u = (c0,cNb1,I0,,INc1)T. Denote the size of vector u by N, which is the number of free parameters, i.e., N = Nb + Nc, then the rhs of Eq. (78) can be written as matrix-vector product, Γu. Denote the total number of measurements of the poloidal flux by MFL. Then matrix elements of Γ for i = 0,,MFL 1 are given by:
     ∫
Γ ij =  G (R′,Z ′,R (FiL),Z(iFL))bj(x,R ′)dR ′dZ ′,
      Ω
(79)

for j = 0,,Nb 1, and

Γ ij = G (R (cj)−Nb,Z(jc)−Nb,R(iFL),Zi(FL)),
(80)

for j = Nb,,N 1. Here (Ri(FL),Zi(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/constraints included (measurements/constraints are merged to Γ by row stacking). We will include MFL measurements of the poloidal flux, MB measurements of the poloidal magnetic field, 1 measurement of plasma current, a constraint for the on-axis safety factor. Then M = MFL + MB + 1 + 1. (For EAST, MFL = 35,MB = 38.)

Next, we consider the poloidal field measurement:

Bp (Ri,Zi) = Bp ⋅nˆi = BR cos𝜃pi + BZ sin𝜃pi,
(81)

where (Ri,Zi) is the location of the No. i probe, ˆni is a unit vector denoting the positive normal direction of the No. i magnetic probe, cos𝜃pi = ˆni eR, sin𝜃pi = nˆi eZ.

Bp(Ri,Zi) can be inferred from u via

           N −1  ∫
            b∑            ′  ′           ′   ′  ′
Bp(Ri,Zi) =    cj Ω GB (R ,Z ,Ri,Zi)bj(x,R )dR dZ
            j=0
           Nc∑−1     (c)  (c)
         +     GB (Rj ,Zj ,Ri,Zi)Ij.                         (82)
            j=0

where

GB (R′,Z′,Ri,Zi) = GBR(R ′,Z′,Ri,Zi)cos𝜃pi + GBZ (R′,Z′,Ri,Zi)sin𝜃pi,
(83)

and GBR and GBZ is given by Eq. (??) and (??). The rhs of Eq. (82) indicates that the matrix elements of Γ are given by:

     ∫
Γ ij =  GB (R′,Z′,Ri,Zi)bj(x,R ′)dR ′dZ ′,
      Ω
(84)

for (j = 0,,Nb 1), and

          (c)    (c)
Γ ij = GB (Rj−Nb,Zj−Nb,Ri,Zi),
(85)

for (j = Nb,,N 1) . I place the magnetic field equations after the flux equations, i.,e the row number i is in the range [MFL,MFL + MB 1].

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

    ∫            Nb−1  ∫
I =    J dR′dZ′ = ∑  c    b (x,R ′)dR ′dZ′.
 p   Ω  ϕ         j=0 j Ω  j
(86)

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

     ∫        ′  ′  ′
Γ ij =  bj(x,R )dR dZ ,
      Ω
(87)

for (j = 0,,Nb 1), and zero for (j = Nb,,N 1) . This equation is placed after the equations of flux and magnetic field measurement, i.e., its row number is  i = MFL + MB.

Let us consider the constraint of the on-axis safety factor qaxis. A formula for computing qaxis is given by:

      𝜖2 + 1 B
qaxis =----- --ϕaxis-,
      𝜖Raxis μ0Jϕaxis
(88)

where

   ∘ --------------
     ( ∂2Ψ∕∂R2 )
𝜖 =    -2-----2    ,
       ∂ Ψ∕∂Z   axis
(89)

is the ellipticity (elongation) at the magnetic axis. Organizing Eq. (88) as

         2
--1- = 𝜖R2axis-μ0-Jϕaxis,
qaxis   𝜖 + 1g(0)
(90)

and using

        N∑b−1
Jϕ axis =    cjbj(x = 0,Raxis),
        j=0
(91)

equation (90) is written as

   2      Nb∑−1
-𝜖Raxisμ0--    cjbj(x = 0,Raxis) =--1-.
(𝜖2 +1)g(0) j=0                  qaxis
(92)

This provides one linear equation for cj, with the matrix elements given by

         2
Γ ij = -𝜖Raxisμ0--bj(x = 0,Raxis).
      (𝜖2 + 1)g(0)
(93)

for (j = 0,,Nb 1), and zero for (j = Nb,,N 1) .

In HEQ, the ellipticity 𝜖 is calculated by computing the elongation of the innermost magnetic surface found by the contour program.

5.2 Algorithm summary

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

Step 1. Search for the magnetic axis and LCFS of Ψ(k) so that we get the values of Ψ(k) at those locations. These values are needed in computing x = (Ψ(k) ΨM)(ΨB ΨM). The location of the LCFS are also needed because we need to set the plasma current to zero outside the LCFS.

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

Step 3. Solve linear least-square problem to get (c1,,cNb,I1,,INc).

Step 4. Update Ψ using the latest plasma and coil currents:

Ψ(k+1) = GJϕ(R(k),c 1,c2,c3,c4)dRdZ+ n=1Nc GIn,

or in its discrete form:

Ψij(k+1) = i,jGijij(          )
 ∑Nb
     cnbn,i′j′
 n=1dSij + n=1Nc Gn,ij(c)I n.

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

 

 

5.3 Vertical displacement instability stabilizer in free-boundary solver[?]

 

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

Ivc = gz--------------BR,vacuum(Rcur,Zcur)--------------
GBR (Rvc,Zvc,Rcur,Zcur) − GBR (Rvc,− Zvc,Rcur,Zcur),

where (Rcur,Zcur) is the location of plasma current center defined by

Rcur = 1-
Ip ΩplRJϕ(R,Z)dRdZ,

Zcur = 1
I-
 p ΩplZJϕ(R,Z)dRdZ,

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

5.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 gg(Ψ) 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. 5.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.

△ ∗Ψ = − μ0RJϕ
(94)

          (     )     2
△∗ = R-∂-  1--∂-  + -∂--                        (95)
      ∂R   R ∂R     ∂Z2
     -∂--  1--∂-   ∂2--
   = ∂R2 − R ∂R +  ∂Z2                          (96)

Method 1: discretize Eq. (95) as

        ( (     )       (     )     )
     -1-   1-∂Ψ-         -1∂Ψ-
  Ri ΔR    R ∂R  i+1∕2 − R ∂R  i−1∕2

+  Ψi,j+1-−-2Ψi,2j +-Ψi,j−1-= − μ0RiJϕi,j,                   (97)
          Δ Z
which can be further discretized as
       (                                    )
  R -1-  --1---Ψi+1,j −-Ψi,j-−--1---Ψi,j −-Ψi−1,j
   iΔR   Ri+1∕2    ΔR       Ri−1∕2    ΔR
  Ψi,j+1 −-2Ψi,j +-Ψi,j−1
+         Δ2         = − μ0RiJϕi,j,                          (98)
           Z
which can be organized as
i1,j + i,j + i+1,j + i,j1 + i,j+1 = μ0RiJϕi,j,

where b = -1-
Δ2R(--Ri-   -Ri--)
 Ri+1∕2 + Ri−1∕2-2-
Δ2Z, a = -1-
Δ2R-Ri--
Ri−1∕2, c = -1-
Δ2R--Ri-
Ri+1∕2, d = e = 1--
Δ2Z. Method 2: discetize Eq. (96)

Ψi+1,j-− 2Ψij +-Ψi−1,j − 1-Ψi+1,j −-Ψi−1,j-+ Ψi,j+1 −-2Ψi,j +-Ψi,j−1
        Δ2R           R     2ΔR                Δ2Z
= − μ0RiJ ϕi,j,
which can be organized as
i1,j + i,j + i+1,j + i,j1 + i,j+1 = μ0RiJϕi,j,

where b = Δ22-
 R Δ22-
 Z a = Δ12-
 R + 2R1ΔR-, c = Δ12-
  R 2R1ΔR-, d = e = Δ12-
  Z

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

5.5 Semi-free bounday equilibrium problem

Semi-free boundary equilibrium problems refer to the problems 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. This mode can be used to set the feedforward PF coil currents when a target plasma state time evolution is given.

6 Magnetic control (to be continued)

Using magnetic measurements as input, a realtime equilibrium reconstruction code (e.g., rtefit) can infer the inner structure of Ψ (observer). On the other hand, PF coil currents can provide actions on Ψ (actuactor). With the observer and actuator, 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 (Ri,Zi) with i = 1,2,M. The actions (coil current changes) are chosen to minimise the following residual:

M∑  (      N∑c                  ) 2
   (Δ Ψi −   G (R ′j,Z ′j,Ri,Zi)ΔIj)  ,
 i         j
(99)

where

ΔΨi = Ψ (target) − Ψ (oild),
(100)

      (new)   (old)
ΔIj = Ij   − Ij  ,
(101)

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 desired positions, then the above sum can be extended to include the following:

     (                                           )2
  M∑X       (X)  (X)   N∑c       ′  ′  (X)  (X)
     (BR (Ri  ,Zi  )+    GBR (Rj,Zj,Ri  ,Zi  )ΔIj )
   i                   j
  M∑X (                N∑c                         )2
+    (BZ (R(iX),Z(iX))+    GBZ (R′j,Z′j,R(iX),Z(iX))ΔIj ) ,         (102)
   i                   j
where (Ri(X),Zi(X)) are desired poistions of X points.

7 Circuit equation

Faraday’s law

          ∂B-
∇ × E = − ∂t ,
(103)

can be written in the integral form:

∮            ∫
  E ⋅dl = − ∂   B⋅dS.
           ∂t
(104)

(The direction of the loop integration is related to the surface orientation by the right-hand rule.) I.e.,

∮
           ∂Ψp-
  E ⋅dl = − ∂t ,
(105)

where Ψp is the flux linked with the loop.

7.1 Circuit equation for PF coil currents

For the No. k coil, the time derivative of the flux through it can be expressed as

    |   (          )
∂Ψp-||   (∑Nc    dIj)   -d∫
 ∂t |k =     Mjk  dt  + dt   Mk(r)JϕdS,
         j=1
(106)

where Mjk is the mutual indctance between the No. j coil and No. k coil, Mk(r) is the mutual inductance between the No. k coil and a plasma filament located at r, and JϕdS is the current carried by the filament. The integration is over a fixed region (time independent) that is large enough to include all the plasma and small enough to not include any coil.

Define the effective mutual inductance between the plasma and the No. k coil by

       1-∫
Mk,p = Ip  Mk (r)J ϕdS,
(107)

where Ip is the plasma current, then expression (106) is written as

   |   ( N∑c        )
∂Ψp|| = (    Mjk dIj) + d-(IpMk,p)
∂t |k    j=1    dt     dt
       (           )
       ( N∑c     dIj)       dIp
     ≈      Mjk dt   + Mk,pdt ,                    (108)
         j=1
where the term proportional to the time derivative of Mk,p is ignored (is this a good approximation?). Including the above emf in the circuit equation, we get
     (          )
       N∑c    dIj        dIp
Vk − (    Mjk-dt) − Mk,p-dt = Ikℛk,
       j=1
(109)

where k is the resistence of the coil, Ik is the current in the present coil,  V k is the voltage imposed on the coil by the power supply. (Here the positive direction of V k and Ik must be chosen consistent with the loop integration discussed above i.e., anti-clockwise in top view.)

The VDE controlling coils are made of two coils (one upper and one lower, with updown symmetry), which are anti-conneted in series and powered by one voltage source. For the upper (U) coil, the circuit equation is

    (  Nc       )
VU −( ∑  MjU dIj) − MU,p dIp-= IUℛU .
      j=1     dt         dt
(110)

Similarly, for the lower (L) coil, the circuit equation is

     ( Nc       )
V  − (∑  M   dIj) − M    dIp-= I ℛ  .
 L    j=1  jL dt      L,p dt    L  L
(111)

The anti-series connection implies that

IL = − IU,
(112)

and

VU − VL = VVDE,
(113)

where V VDE is the source voltage (assuming the voltage source positive terminal is connected to the upper coil that goes anti-clockwise in top view). Combining these equations, we obtain

       ( N               )
       (∑ c            dIj)                dIp
VVDE −     (MjU − MjL )dt   − (MU,p − ML,p)dt = IU(ℛU +ℛL ).
        j=1
(114)

This is used as an equation for the upper VDE coil current. Eq. (112) is used as an equation for the lower VDE coil current.

7.2 Circuit equation for total plasma current

Because the mutual inductance between two objects is symmetric (i.e., MAB = MBA), the mutual inductance Mk,p given by expression (107) can also be used to calculate the effective emf induced in the plasma by the coils:

  (         )
d- ∑Nc            N∑c     dIk
dt    Mk,pIk  ≈ −    Mk,p dt .
   k=1            k=1
(Here Mk,p can be understood as Jϕ weighted average of Mk(r). Why to use Jϕ weighted, rather than just equal-weighted spatial average? This is to ensure energy conservation: the power is E J, so the current density should serve as a weight.)

Similarly, the effective emf self-induced in the plasma is given by

  d ( 1 ∫ ∫                      )
−dt  I-     Jϕi′j′Mi ′j′ijJϕijdSijdSi′j′ ,
      p
(115)

where Mijij is the mutual inductance between ij filament and ijfilament; both the 2D integrations are over a fixed region that includes the plasma but does not include any coil; one of the two 2D integrations can be regarded as Jϕ weighted average and the other is the integration over the source. Expresion (115) motivates us to define the effective self-inductance of plasma by

       ∫ ∫
Lp = 12     Jϕi′j′Mi′j′ijJ ϕijdSijdSi′j′,
     Ip
(116)

Then the effective emf in expression (115) is written as

  d            dI
−--(LpIp) ≈ − Lp-p,
 dt             dt
where the term proportional to the time derivative of Lp is ignored. Then the circuit equation for the plasma current is written as
         ∑Nc
− LpdIp−    Mk,p dIk= Ipℛp,
     dt  k=1     dt
(117)

where p is the effective plasma resistivity.

7.3 Time discrete form of coupled coil and plasma currents

Using the Euler implicit scheme, the coil current evolution equation (109) is written as

V k(n) j=1Nc Mjk (n+1)   (n)
Ij----−-Ij-
     Δt Mk,p (n+1)   (n)
Ip----−-Ip--
    Δt = Ik(n+1) k.

Moving all the unknown (coil currents Ij(n+1) and plasma current Ip(n+1)) to the left-hand side, we get

            N       (n+1)
     (n+1)  ∑c     Ij-----      Ip(n+1)-
  ℛkIk    + j=1Mjk  Δt  + Mk,p  Δt
                (n)
   (n)  N∑c     Ij--      I(np)-
= Vk  +    Mjk Δt  +Mk,p Δt .                        (118)
        j=1
For each coil, we get a linear equation like the above. For Nc coils, we get Nc equations. (Note that Mk,p depends on Jϕ, and thus is time dependent. If we use value of Mk,p at tn+1 in the above equation, then we need to iterate because Mk,p(n+1) is unkown. If we use the value at tn, then no iteration is needed, whch is the case used in HEQ).

Similarly, the plasma current evolution equation (117) is discretized as

   I(pn+1)− I(pn)   N∑c     I(n+1) − I(n)
− Lp----Δt---- −    Mk,p-k--Δt--k--= I(pn+1)ℛp,
                 k=1
(119)

i.e.,

 N             (        )                 N
∑ cMk,p- (n+1)         Lp-  (n+1)     I(pn)-  ∑c     I(kn)
    Δt  Ij   +   ℛp + Δt  Ip    = LpΔt  +    Mk,p Δt .
k=1                                       k=1
(120)

 

7.4 Reinforcement learning –to be continued

Define the plasma state by the position of the magnetic axis (Ra,Za), total plasma current Ip, and the shape of the LCFS. Discretize the LCFS in the cylindrical coordinates by Ns points. Then the plasma state can be represented by 2Ns + 3 real numbers: Ra, Za, Ip, and (Ri,Zi) for i = 1,2,Ns.

Inputs to the policy network: the magnetic measurements at the present time step tn (which serves as a proxy to the actual plasma state), and target plasma states (provided by the designer) at tn,tn+1,,tn+k (where k may be a small integer, say k = 1).

Output of the policy network: voltage command to the coils imposed at tn which can hopefully make the plasma follow the target state in future time steps.

How to train the policy network on a plasma simulator?

 

8 Coordinates based on magnetic field structure

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 ζ. 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. 15.

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

8.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

r = xˆx+ yˆy + zˆz,
(121)

where r is the location vector of the point; ˆx, ˆy, and ˆz 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, (x1,x2,x3), can be expressed as

r = x(x1,x2,x3)ˆx + y(x1,x2,x3)ˆy+ z(x1,x2,x3)ˆz.
(122)

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

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

x = x(x1,x2,x3)
y = y(x1,x2,x3)
z = z(x1,x2,x3)
(123)

8.2 Jacobian

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

    |           |
    || ∂x-∂x- ∂x-||
𝒥 = || ∂x∂1y-∂∂xy2 ∂x∂3y-||,
    || ∂x∂1z-∂∂xz2 ∂x∂3z-||
      ∂x1 ∂x2 ∂x3
(124)

which can also be written as

     ∂r-- -∂r- -∂r-
𝒥 =  ∂x1 × ∂x2 ⋅∂x3 .
(125)

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

                    −1
𝒥  = (∇x1 × ∇x2 ⋅∇x3 ) .
(126)

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. (124), 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:

    |        |  |                    |
    ||∂∂xR-∂∂xϕ ∂∂xZ||  || ∂R∂cRosϕ-∂Rc∂oϕsϕ ∂R∂cZosϕ-||
𝒥 = ||∂∂yR-∂∂yϕ ∂∂yZ||= || ∂R∂sRinϕ-∂R∂siϕnϕ ∂R∂sZinϕ-||
    ||∂z-∂z ∂z||  ||  ∂Z-    ∂Z-   ∂Z-  ||
    |∂R ∂ϕ ∂Z     |∂R     ∂ϕ    ∂Z
    ||cosϕ − R sinϕ 0||
  = ||sinϕ  R cos ϕ 0||= R
    | 0     0    1|
If the Jacobian of a coordinate system is greater than zero, it is called a right-handed coordinate system. Otherwise it is called a left-handed system.

8.3 Orthogonality relation between two sets of basis vectors

In a curvilinear coordinate system (x1,x2,x3), there are two kinds of basis vectors: xi and r∕∂xi, with i = 1,2,3. These two kinds of basis vectors satisfy the following orthogonality relation:

     ∂r--
∇xi ⋅∂xj = δij,
(127)

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

          [(    )    (    )    (    )  ] [                                       ]
∇x  ⋅ ∂r-=   ∂xi  ˆx+   ∂xi  ˆy+   ∂xi  ˆz ⋅ -∂x-ˆx+ x-∂ˆx-+ -∂yyˆ+ y ∂-ˆy-+ ∂z-ˆz +z ∂ˆz-
  i  ∂xj     ∂x        ∂y        ∂z       ∂xj     ∂xj   ∂xj     ∂xj   ∂xj     ∂xj
          [( ∂xi)    ( ∂xi)    ( ∂xi)  ] [ ∂x        ∂y        ∂z     ]
        =    ∂x-  ˆx+   ∂y-  ˆy+   ∂z-  ˆz ⋅ ∂xj-ˆx+ 0+  ∂xj ˆy+ 0 + ∂xjˆz+ 0
           ∂x ∂x    ∂x  ∂y   ∂x  ∂z
        =  --i----+ --i----+ --i----
           ∂x ∂xj   ∂y ∂xj   ∂z ∂xj
        =  ∂xi-                                                               (128)
           ∂xj
        = δij,
where the second equality is due to ˆx∕∂xj = 0,∂ˆy∕∂xj = 0,∂ˆz∕∂xj = 0 since ˆx,ˆy,ˆz are constant vectors independent of spatial location; the chain rule has been used in obtaining Eq. (128)]

[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. (127). In this case, x = x1 cosx2, y = x1 sinx2, z = x3, where x1 R, x2 ϕ, x3 Z.]

It can be proved that xi is a contravariant vector while r∕∂xi 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. (127) 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∕∂x1 is orthogonal to x2 and x3, thus, r∕∂x1 can be written as

-∂r-
∂x1 = A∇x2 × ∇x3,
(129)

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

1 = A (∇x2 × ∇x3)⋅∇x1,
(130)

which gives

A = -------1--------
    (∇x2 × ∇x3 )⋅∇x1
 = 𝒥                                         (131)
Thus r∕∂x1 is written, in terms of x1, x2, and x3, as
 ∂r
∂x1-= 𝒥 ∇x2 × ∇x3.
(132)

Similarly, we obtain

-∂r-= 𝒥 ∇x  × ∇x
∂x2       3     1
(133)

and

-∂r-= 𝒥 ∇x1 × ∇x2.
∂x3
(134)

Equations (132)-(134) can be generally written

 ∂r
--- = 𝒥∇xj × ∇xk,
∂xi
(135)

where (i,j,k) represents the cyclic order in the variables (x1,x2,x3). Equation (135) expresses the covariant basis vectors in terms of the contravariant basis vectors. On the other hand, from Eq. (132)-(134), we obtain

∇x  = 𝒥− 1 ∂r-×-∂r-,
   i      ∂xj  ∂xk
(136)

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

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

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

eψ ≡ ∇ ψ; e𝜃 ≡ ∇𝜃; eζ ≡ ∇ζ.                    (137)
and the covariant basis vectors are denoted with low indices as
     ∂r              ∂r
eψ ≡ --; e𝜃 ≡ ∂∂r𝜃; eζ ≡-.                       (138)
     ∂ψ              ∂ζ
Then the orthogonality relation, Eq. (127), is written as
eα ⋅eβ = δαβ.
(139)

In term of the contravairant basis vectors, A is written

        ψ      𝜃     ζ
A  = Aψe  + A𝜃e + A ζe ,
(140)

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

A  = Aψeψ + A𝜃e𝜃 + A ζeζ,
(141)

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

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

        𝜃   ζ
eψ = 𝒥 e × e
(142)

e𝜃 = 𝒥eζ × eψ
(143)

eζ = 𝒥eψ × e𝜃
(144)

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

eψ = 𝒥 −1e𝜃 × eζ
(145)

 𝜃    − 1
e  = 𝒥  eζ × eψ
(146)

 ζ    − 1
e  = 𝒥  eψ × e𝜃
(147)

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

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

      ∂f      ∂f     ∂f
∇f =  --∇ ψ + --∇ 𝜃+ ---∇ζ.
      ∂ψ      ∂𝜃     ∂ ζ
(148)

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

∂f-= (𝒥 ∇𝜃 ×∇ ζ)⋅∇f =  ∂r-⋅∇f
∂ψ                     ∂ψ
(149)

∂f-= (𝒥 ∇ζ ×∇ ψ)⋅∇f  = ∂r-⋅∇f
∂𝜃                     ∂𝜃
(150)

∂f-= (𝒥 ∇ψ × ∇𝜃)⋅∇f  = ∂r-⋅∇f
∂ζ                     ∂ζ
(151)

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

∇ ψ ⋅∇f = |∇ ψ|2∂f-+ (∇𝜃 ⋅∇ψ)∂f-+ (∇ ζ ⋅∇ ψ)∂f.
               ∂ψ            ∂𝜃           ∂ζ
(152)

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

∇ ⋅(∇ α× ∇ β) = 0,
(153)

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

     (ψ)            (𝜃)             (ζ)
A = A   𝒥∇ 𝜃× ∇ ζ + A 𝒥 ∇ ζ × ∇ ψ+ A  𝒥∇ ψ × ∇𝜃,
(154)

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

∇  ⋅A  = (∇ 𝜃× ∇ ζ)⋅∇ (A (ψ)𝒥 )+ (∇ ζ × ∇ ψ)⋅∇ (A(𝜃)𝒥 )+ (∇ψ × ∇𝜃) ⋅∇(A(ζ)𝒥 ) (155)
        1 ( ∂A(ψ)𝒥   ∂A (𝜃)𝒥   ∂A (ζ)𝒥 )
      = 𝒥-  --∂ψ---+ --∂𝜃---+ --∂ζ--- ,                               (156)
where the second equality is obtained by using Eqs. (149), (150), and (151).

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

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

∇2f = ∇ ⋅∇f
         ( ∂f      ∂f     ∂f   )
    = ∇ ⋅ ∂-ψ∇ ψ+  ∂𝜃∇ 𝜃+ ∂ζ-∇ζ  .                  (157)
To proceed, we can use the divergence formula (156) to express the divergence in the above expression. However, the vector in the above (blue term) is not in the covariant form desired by the divergence formula (156). If we want to directly use the formula (156), we need to transform the vector (blue term in expression (157)) to the covariant form. This process seems to be a little complicated. Therefore, I choose not to use this method. Instead, I try to simplify expression (157) by using basic vector identities:
  2     ( ∂f)         ( ∂f)         (∂f )
∇ f = ∇   ∂ψ- ⋅∇ ψ+ ∇   ∂𝜃- ⋅∇ 𝜃+ ∇  -∂ζ  ⋅∇ζ
        ∂f       ∂f      ∂f
      + --∇2 ψ + --∇2 𝜃+ ---∇2ζ.                          (158)
        ∂ψ       ∂𝜃      ∂ζ
Using the gradient formula, the above expression is further written as
∇2f = ∂2f-|∇ ψ|2 +-∂2f-∇ ψ⋅∇ 𝜃+ -∂2f-∇ ψ⋅∇ ζ
      ∂ψ2        ∂ψ ∂𝜃         ∂ψ ∂ζ
         ∂2f          ∂2f    2   ∂2f
      + ∂𝜃∂ψ-∇𝜃 ⋅∇ ψ + ∂𝜃2|∇𝜃| + ∂𝜃∂ζ-∇𝜃 ⋅∇ζ
         ∂2f           ∂2f          ∂2f
      + -----∇ζ ⋅∇ψ + ----∇ ζ ⋅∇ 𝜃+---2|∇ ζ|2
        ∂ζ∂ψ          ∂ζ∂𝜃         ∂ ζ
      + ∂f-∇2ψ + ∂f∇2 𝜃+ ∂f-∇2ζ,                         (159)
        ∂ψ       ∂𝜃      ∂ζ
and can be simplified as
        2           2              2
∇2f =  ∂-f2|∇ψ |2 + 2-∂-f-∇ ψ⋅∇ 𝜃+ 2-∂-f-∇ ψ⋅∇ ζ
       ∂ψ         ∂ ψ∂𝜃          ∂ ψ∂ζ
      + ∂2f|∇ 𝜃|2 + 2-∂2f-∇𝜃 ⋅∇ζ
        ∂𝜃2        ∂𝜃∂ ζ
        ∂2f    2
      + ∂ζ2|∇ ζ|
        ∂f       ∂f      ∂f
      + ---∇2ψ + --∇2 𝜃+ ---∇2ζ.                          (160)
        ∂ψ       ∂𝜃       ∂ζ
Assume (ψ,𝜃,ζ) are field-line following coordinates with r∕∂𝜃 along the field line direction, then neglect all the parallel derivatives, i.e., derivative over 𝜃, then the above expression is reduced to
∇2f = ∂2f|∇ ψ|2 + 2 ∂2f-∇ ψ ⋅∇ζ + ∂2f|∇ζ|2
      ∂ψ2         ∂ψ∂ζ          ∂ζ2
      ∂f- 2    ∂f- 2
    + ∂ψ ∇ ψ + ∂ζ∇  ζ.                                 (161)
This approximation reduces the Laplacian operator from being three-dimensional to being two-dimensional. This approximation is often called the high-n approximation, where n is the toroidal mode number (mode number along ζ direction).

8.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

∇ × A = ∇ × (A1∇ψ + A2∇ 𝜃+ A3∇ ζ)
      = ∇A1 × ∇ ψ+ ∇A2  × ∇𝜃+ ∇A3  × ∇ζ
        ( ∂A      ∂A    )        ( ∂A       ∂A    )        ( ∂A       ∂A    )
      =   --1∇ 𝜃+ ---1∇ζ  × ∇ ψ+   ---2∇ψ + ---2∇ζ  × ∇ 𝜃+   --3∇ ψ+  --3∇ 𝜃  ×∇ ζ
          ∂(𝜃       ∂ζ )             ∂ψ(      ∂ζ  )           ∂ψ  (    ∂𝜃    )
      = 1-  ∂A2-− ∂A1- 𝒥 ∇ ψ× ∇ 𝜃+ -1  ∂A1-− ∂A3-  𝒥∇ ζ × ∇ ψ + 1  ∂A3-− ∂A2- 𝒥 ∇𝜃 ×(1∇6ζ2.)
        𝒥   ∂ ψ    ∂𝜃              𝒥    ∂ζ    ∂ψ               𝒥   ∂𝜃    ∂ ζ
Note that taking the curl of a vector in the covariant form leaves the vector in the contravariant form.

8.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 (135) and (136) 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

      1            2             3
∇ψ = a 𝒥 ∇𝜃× ∇ ζ + a 𝒥 ∇ζ × ∇ψ + a 𝒥∇ ψ × ∇𝜃.
(163)

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

 1      2
a = |∇ψ |,
(164)

a2 = ∇ ψ⋅∇ 𝜃,
(165)

 3
a = ∇ ψ⋅∇ ζ.
(166)

Similarly, we write

∇ 𝜃 = b1𝒥∇ 𝜃× ∇ ζ + b2𝒥∇ ζ × ∇ ψ + b3𝒥 ∇ψ × ∇ 𝜃,
(167)

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

b1 = ∇ 𝜃⋅∇ ψ
(168)

 2      2
b = |∇𝜃| ,
(169)

b3 = ∇ 𝜃⋅∇ ζ.
(170)

The same situation applies for the ζ basis vector,

∇ ζ = c1∇𝜃 × ∇ζ𝒥 + c2∇ζ × ∇ψ 𝒥 + c3∇ ψ × ∇𝜃𝒥 ,
(171)

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

 1
c = ∇ ζ ⋅∇ ψ
(172)

c2 = ∇ ζ ⋅∇ 𝜃
(173)

c3 = |∇ ζ|2
(174)

Summarizing the above results in matrix form, we obtain

(     )   (                      ) (           )
  ∇ ψ        |∇ ψ|2 ∇ ψ⋅∇ 𝜃 ∇ ψ ⋅∇ζ    ∇𝜃 ×∇ ζ𝒥
(  ∇𝜃 ) = ( ∇𝜃 ⋅∇ψ  |∇𝜃|2  ∇𝜃 ⋅∇ ζ) ( ∇ ζ × ∇ ψ𝒥 )
   ∇ζ       ∇ζ ⋅∇ψ ∇ ζ ⋅∇ 𝜃 |∇ ζ|2     ∇ ψ × ∇𝜃𝒥
(175)

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

∇𝜃 × ∇ζ𝒥 = d1∇ ψ+ d2∇ 𝜃+ d3∇ ζ
(176)

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

d1 = |∇ 𝜃× ∇ ζ|2𝒥 2
(177)

d2 = (∇ 𝜃× ∇ ζ𝒥)⋅(∇ ζ × ∇ ψ𝒥)
(178)

d3 = (∇ 𝜃× ∇ ζ𝒥)⋅(∇ ψ× ∇ 𝜃𝒥)
(179)

For the second contravariant basis vector

∇ϕ × ∇ψ 𝒥 = e1∇ ψ + e2∇ 𝜃+ e3∇ ζ
(180)

                          2
e1 = (∇ ζ × ∇ ψ)⋅(∇ 𝜃× ∇ ζ)𝒥
(181)

e2 = (∇ ζ × ∇ ψ)⋅(∇ ζ × ∇ψ )𝒥 2
(182)

e3 = (∇ ζ × ∇ ψ)⋅(∇ ψ× ∇ 𝜃)𝒥 2
(183)

For the third contravariant basis vector

∇ψ × ∇𝜃𝒥  = f1∇ ψ + f2∇ 𝜃+ f3∇ ζ
(184)

f1 = (∇ψ × ∇𝜃)⋅(∇ 𝜃× ∇ ζ)𝒥 2
(185)

f2 = (∇ ψ × ∇𝜃)⋅(∇ ζ × ∇ψ )𝒥 2
(186)

f3 = (∇ ψ × ∇𝜃)⋅(∇ ψ× ∇ 𝜃)𝒥 2
(187)

Summarizing these results, we obtain

(          )      (    )
  ∇𝜃 × ∇ζ𝒥          ∇ψ
( ∇ζ × ∇ψ𝒥 )  = M ( ∇ 𝜃) ,
  ∇ψ × ∇𝜃𝒥          ∇ ζ
(188)

where

M = (               2 2                           2                      2)
(      |∇𝜃 × ∇ζ| 𝒥     2 (∇𝜃 × ∇ζ)⋅(∇ ζ × ∇ψ )𝒥 2 (∇𝜃 ×∇ ζ)⋅(∇ψ × ∇ 𝜃)𝒥 2)
  (∇ζ × ∇ψ )⋅(∇𝜃 ×∇ ζ)𝒥2 (∇ ζ × ∇ ψ)⋅(∇ ζ × ∇ ψ)𝒥2 (∇ ζ × ∇ ψ)⋅(∇ ψ× ∇ 𝜃)𝒥2
  (∇ ψ× ∇ 𝜃)⋅(∇𝜃× ∇ ζ)𝒥  (∇ ψ × ∇𝜃)⋅(∇ ζ × ∇ ψ)𝒥 (∇ ψ × ∇𝜃)⋅(∇ ψ× ∇ 𝜃)𝒥,

This matrix and the matrix in Eqs. (175) 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. (175) is simplified to

( ∇ψ )   (  |∇ ψ|2 ∇ ψ ⋅∇𝜃  0  ) ( ∇ 𝜃× ∇ ϕ𝒥 )
( ∇𝜃 ) = ( ∇ ψ⋅∇ 𝜃 |∇ 𝜃|2   0  ) ( ∇ ϕ× ∇ ψ𝒥 )
  ∇ϕ          0      0    1∕R2    ∇ ψ× ∇ 𝜃𝒥
(189)

Similarly, Eq. (188) is simplified to

(  ∇𝜃× ∇ ϕ𝒥 )   (   |∇𝜃|2𝒥2∕R2   − ∇ 𝜃⋅∇ ψ𝒥 2∕R2 0 ) ( ∇ ψ)
( ∇ ϕ× ∇ ψ𝒥 ) = ( − ∇ψ ⋅∇ 𝜃𝒥2∕R2   |∇ ψ|2𝒥 2∕R2   0 ) ( ∇ 𝜃)
  ∇ ψ × ∇𝜃𝒥              0             0       R2     ∇ ϕ
(190)

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

(                                 ) (                    )
    |∇ 𝜃|2𝒥 2∕R2  − ∇ 𝜃⋅∇ψ 𝒥2∕R2  0      |∇ ψ|2 ∇ ψ⋅∇ 𝜃  0
( − ∇ ψ ⋅∇𝜃𝒥 2∕R2   𝒥R22|∇ψ |2     0 ) ( ∇ψ ⋅∇𝜃  |∇𝜃|2   0  ) ,
        0              0       R2        0      0   1∕R2
(191)

can be calculated to give

(       )
  A  0 0
(  0 A 0 )
   0 0 1,

where

A = |∇𝜃|2|∇ψ|2𝒥2∕R2 (𝜃 ⋅∇ψ)2𝒥2.

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

                         2
|∇𝜃|2|∇ ψ|2 − (∇𝜃 ⋅∇ ψ)2 = R-,
                        𝒥2
(192)

]

8.10 Covariant/contravariant representation of equilibrium magnetic field

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

B = ∇ Ψ × ∇ϕ + g∇ϕ.
(193)

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

B  = − Ψψ∇ ϕ× ∇ ψ − Ψ𝜃∇ϕ × ∇𝜃 + g∇ϕ,
(194)

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

                               𝒥
B = − Ψψ∇ ϕ ×∇ ψ − Ψ𝜃∇ϕ × ∇𝜃 + g-2∇ ψ × ∇𝜃.
                               R
(195)

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

    (                        )      (                          )
B =  Ψ  𝒥-∇ ψ ⋅∇𝜃 + Ψ 𝒥-|∇ 𝜃|2  ∇ψ +  − Ψ -𝒥-|∇ ψ|2 − Ψ-𝒥-∇𝜃 ⋅∇ψ  ∇ 𝜃+ g∇ ϕ.
       ψR2           𝜃R2                ψR2         𝜃R2
(196)

For the convenience of notation, define

hαβ = 𝒥-∇ α⋅∇ β,
      R2
(197)

then Eq. (196) is written as

B  = (Ψψh ψ𝜃 + Ψ𝜃h𝜃𝜃)∇ ψ + (− Ψψhψψ − Ψ𝜃hψ𝜃)∇𝜃 + g∇ ϕ.
(198)

8.11 Flux coordinates (ψ,𝜃,ϕ)

A coordinate system (ψ,𝜃,ϕ), where ϕ is the usual cylindrical toroidal angle, is called a magnetic/flux 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. (195), is written as

                   𝒥
B = − Ψ ′∇ ϕ× ∇ ψ +g-2∇ ψ ×∇ 𝜃,
                   R
(199)

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

    (            )      (           )
B =   Ψ′ 𝒥-∇ ψ ⋅∇𝜃 ∇ψ +  − Ψ′ 𝒥-|∇ ψ|2 ∇ 𝜃+ g∇ ϕ.
        R2                   R2
(200)

8.12 Local safety factor

The local safety factor ˆq is defined by

   B ⋅∇ ϕ
ˆq =------,
   B ⋅∇ 𝜃
(201)

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. (199), into the above equation, the local safety factor is written

ˆq(ψ,𝜃) = − g2-𝒥′.
          R  Ψ
(202)

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

In terms of ˆq , the contravariant form of the magnetic field, Eq. (199), is written

B = − Ψ′(∇ϕ × ∇ψ + ˆq∇ψ × ∇ 𝜃).
(203)

and the parallel differential operator B0 ⋅∇ is written as

                                             (         )
B  ⋅∇ = − Ψ′(∇ ϕ× ∇ ψ + ˆq∇ ψ × ∇𝜃)⋅∇ = − Ψ′𝒥 −1 -∂-+ ˆq-∂- .
  0                                           ∂ 𝜃   ∂ϕ
(204)

If ˆq 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 𝒥1). This simplification is useful because different poloidal harmonics are decoupled in this case. We will discuss this issue futher in Sec. 14.

 

8.13 Global safety factor

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

          ∫ 2π
q(ψ) ≡ -1-    ˆqd𝜃                              (205)
       2π  0
        -1-g-∫ 2π 𝒥--
     = −2π Ψ′  0  R2d𝜃.                        (206)
Note that q and ˆq defined this way can be negative, which depends on the choice of the positive direction of ϕ and 𝜃 coordinates (note that the safety factor given in G-eqdsk file is always positive, i.e. it is the absolute value of the safety factor defined here).

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

             ∮
q(ψ) = − 1-g′  sign(𝒥 )-dℓp--
        2πΨ        ∮  R|∇ψ |
    = − 1-g-sign(𝒥)   -dℓp--.                     (207)
        2π sign(Ψ ′)   R |∇ Ψ|
Expression (207) is used in the GTAW code to numerically calculate the value of q on magnetic surfaces (as a benchmarking of the q profile specified in the G-eqdsk file). Expression (207) can also be considered as a relation between q and g. In the equilibrium problem where q is given (fixed-q equilibrium), we can use expression (207) to obtain the corresponding g (which explicitly appears in the GS equation):
        ( ∮ --dℓp--)−1 sign(Ψ′)
g = − 2πq   R |∇ Ψ|    sign(𝒥 ).
(208)

We note that expression (208) 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 (208) can be performed, yielding the values of g.)

Using Bp = |∇Ψ|∕R and Bϕ = g∕R,  the absolute value of q in expression (207) is written

        ∮
|q| =-1-g  --dℓp-.                           (209)
    2π    R |∇ Ψ|
     1 ∮  1|Bϕ|
  = 2π-  R--Bp-dℓp                          (210)
which is the familiar formula we see in textbooks.

8.14 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

dl =-∂rdψ + ∂rd𝜃 + ∂rdϕ
    ∂ψ      ∂𝜃     ∂ϕ
(211)

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

     ∂r-
dℓp = ∂ 𝜃d𝜃
   = 𝒥 ∇ϕ × ∇ ψd𝜃.                          (212)
We use the convention that dℓp and d𝜃 take the same sign, i.e.,
dℓp = |𝒥 ∇ϕ × ∇ψ |d𝜃.
(213)

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

       R
d𝜃 = ------dlp
     |𝒥 ∇ ψ|
(214)

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

—–

-  ∫ 𝜃       ∫ 𝜃 g--𝒥
δ = 0 ˆqd𝜃 = − 0  R2Ψ ′d𝜃
(215)

8.15 Calculating poloidal angle

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

           ∫ xi,j   R
𝜃i,j = 𝜃ref,j +     -----dlp,
            xref,j |𝒥 ∇ψ |
(216)

where the curve integration is along the contour Ψ = Ψj, xref,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. (216). Denote the Jacobian of the constructed coordinates by 𝒥′, then

𝒥′ = ± 𝒥
(217)

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

s(ψ) = ∮--2π----.
        |𝒥R∇ψ|dlp
(218)

Then 𝜃 is written as

-
𝜃i,j = s(ψj)𝜃i,j
          (       ∫ xi,j        )
   = s(ψj) 𝜃ref,j +     --R---dlp
                   xref,j|𝒥∇ ψ|
     -          ∫ xi,j --R---
   = 𝜃ref,j + s(ψj) x    |𝒥 ∇ψ |dlp,                     (219)
                  ref,j
where 𝜃ref,j = s(ψj)𝜃ref,j. Sine we have modified the definition of the poloidal angle, the Jacobian of the new coordinates (ψ,𝜃) is different from that of (ψ,𝜃,ϕ). The Jacobian 𝒥new of the new coordinates (ψ,𝜃) is written as
       ------1-------
𝒥new ≡ (∇ ψ × ∇𝜃)⋅∇ ϕ
               1
    =  {∇ψ-×-∇[s(ψ)𝜃]}⋅∇-ϕ-

    =  -1--------1-------
       s(ψ )(∇ψ × ∇𝜃)⋅∇ ϕ
    =  -1--𝒥′
       s(ψ )
        -1--
    = ± s(ψ)𝒥                                   (220)
The Jacobian 𝒥 can be set to various forms to achieve various poloidal coordinates, which will be discussed in the next section. After the Jacobian and the radial coordinate ψ is chosen, all the quantities on the right-hand side of Eq. (219) are known and the integration can be performed to obtain the value of 𝜃i,j of each point on each flux surface.

[The reference points xref,j and the values of poloidal angle at these points can be chosen by users. One choice of the reference points xref,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

      R
𝒥 = ---- .
    |∇ ψ|
(221)

Then 𝜃i,j in Eq. (219) is written

                 ∫
𝜃   = 𝜃   + ∮2π--  xi,jdl.
 i,j    ref,j    dlp  xref,j p
(222)

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

        ∮ dl  R      ∮ dl  R  dΨ
𝒥new = ±---p-----= ± ---p--------.
         2π  |∇ ψ|     2π  |∇ Ψ|dψ
(223)

Equation (222) 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. (222) 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

𝒥 = h (ψ ),
(224)

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

-    -        2π    ∫ xi,j  R
𝜃i,j = 𝜃ref,j + ∮-R--dl      |∇ψ|dlp
             |∇ψ| p  xref,j
(225)

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

         1 ∮   R
𝒥new = ± 2π-  |∇-ψ|dlp.
(226)

Note that both 𝜃i,j and 𝒥new are independent of the function h(ψ) introduced in Eq. (224). (h(ψ) is eliminated by the normalization procedure specified in Sec. 8.15 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∕ψ)dlp be independent of ψ, so that 𝒥new in Eq. (226) is constant in space.**)

Jacobian of Boozer’s poloidal angle

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

         h(ψ)
𝒥(ψ,𝜃) = -B2-.
(227)

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

                    ∫ xi,j  2
𝜃i,j = 𝜃ref,j + ∮-B22πR--      B-R-dlp.
              |∇-ψ|dlp xref,j |∇ ψ|
(228)

The final Jacobian is given by

        ∮  2
          B|∇-Rψ|dlp 1
𝒥new = ±---2π---B2-.
(229)

The usefullness of Boozer poloidal angle will be further discussed in Sec. 14.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

𝒥 (ψ,𝜃) = R2,

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

                     ∫ xi,j
𝜃i,j = 𝜃ref,j + ∮-2π---      --1--dlp,
              R1|∇-ψ|dlp xref,j R|∇ψ |
(230)

The Jacobian 𝒥new given by Eq. (220) now takes the form

          ∮
         2  R|1∇ψ|dlp
𝒥new = ±R  ---2π----.
(231)

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

            ∫ 𝜃
𝜗 = 𝜗ref,j + 1    qˆd𝜃,
           q 𝜃ref,j
(232)

where ˆq is the local safety factor corresponding to the arbitrary poloidal angle 𝜃, i.e., ˆq = B⋅∇ϕ∕(B⋅∇𝜃). [Proof: Using d𝜃 = |𝒥R∇ψ-|dlp, the poloidal angle 𝜗 given in Eq. (230) is written as

                    ∫ x
𝜗 = 𝜗   + ∮---2π---    i,j--1---dl
     ref,j    R|1∇ψ|dlp  xref,jR |∇ ψ| p
                2π      ∫ xi,j   1   |𝒥 ∇ψ |
  = 𝜗ref,j + ∮--1--|𝒥∇ψ|--      -----------d𝜃
            R|∇ψ|--R--d𝜃 xref,j R|∇ψ |  R
          --2π---∫ xi,j|𝒥-|
  = 𝜗ref,j + ∮ |𝒥|d𝜃 x   R2 d𝜃.                            (233)
            R2     ref,j
Using ˆq = g-
R2𝒥-
Ψ′, where 𝒥 is the Jacobian of (ψ,𝜃,ϕ) coordinates, then the above 𝜗 is reduced to expression (232).]

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

𝒥  = -12-∼ R2.
     B
(234)

Comparison between different types of poloidal angle

 

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

         i
𝒥  = --R-j-k,
     |∇ ψ| B
(235)

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.


pict pict

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 ∘ Ψ--
    t = 0.01714 (the innermost magnetic surface) and end at ∘ Ψ--
    t = 0.9851 (the boundary magnetic surface), and are equally spaced in ∘ ---
  Ψt.

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.,

       ------1-------
𝒥new = (∇ψ × ∇ 𝜃) ⋅∇ϕ ,
(236)

which can be further written as

𝒥new = R (R 𝜃Zψ − RψZ𝜃),
(237)

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

--  Ψ − Ψ0
Ψ = Ψa-−-Ψ0,
(238)

where Ψ0 and Ψa are the values of Ψ at the magnetic axis and LCFS, respectively. Other choices of the radial coordinates: the normalized toroidal magnetic flux and its square root, Φ, and √ --
  Φ.

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

 3
d V = 𝒥 dψd𝜃dϕ.
(239)

Integrating over the toroidal angle, we obtain

 2
d V = 2π𝒥dψd 𝜃
(240)

Further integrating over the poloidal angle, we obtain

           ∮
d1V = 2πdψ   𝒥 d𝜃,
(241)

i.e.,

d1V      ∮
----= 2π   𝒥 d𝜃.
 dψ
(242)

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 𝒥1 = ψ ×∇𝜃 ⋅∇ϕ 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.

9 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. (214) (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 (ψ,𝜃,ϕ).

9.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, Ψ0, 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 Ψ0 to Ψb. Then I try to find the contours of Ψ with contour level value ranging from Ψ0 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.

 


pict

Fig. 7: Verification of the numerical code that calculates the contours of the poloidal flux Ψ. The bold line in the figure indicates the LCFS. The contour lines (solid lines) given by the gnuplot program agrees well with the results I calculate by using interpolation and root-finding method (the two sets of contours are indistinguishable in this scale). My code only calculate the contour lines within the LCFS, while those given by gnuplot contains additional contour lines below the X point and on the left top in the figure. Eqdisk file of the equilibrium was provided by Dr. Guoqiang Li (filename: g013606.07104).

 

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

       ∘ ----------------
         (∂ Ψ)2   (∂ Ψ)2
|∇Ψ | =   ---   +  ---  ,
          ∂R       ∂Z
(243)

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.


pict pict

Fig. 8: Grid points (the intersecting points of two curves in the figure) corresponding to uniform poloidal flux and uniform poloidal arc length for EAST equilibrium shot 13606 at 7.1s (left) (G-file name: g013606.07104) and shot 38300 at 3.9s (right) (G-file name: g038300.03900).


pict

Fig. 9: |∇Ψ| as a function of the poloidal angle. The different lines corresponds to the values of |∇Ψ| on different magnetic surfaces. The stars correspond to the values on the boundary magnetic surface while the plus signs correspond to the value on the innermost magnetic surface (the magnetic surface adjacent to the magnetic axis). The equilibrium is for EAST shot 38300 at 3.9s.


pict pict

Fig. 10: The Poloidal magnetic field Bp = |∇Ψ|∕R (left) and toroidal magnetic field Bϕ = g∕R (right) as a function of the poloidal angle. The different lines corresponds to the values on different magnetic surfaces. The stars correspond to the values on the boundary magnetic surface while the plus signs correspond to the value on the innermost magnetic surface (the magnetic surface adjacent to the magnetic axis). The equilibrium is for EAST shot 38300 at 3.9s.

 

 


pict

Fig. 11: Equal-arc Jacobian as a function of the poloidal angle on different magnetic surfaces. The dotted line corresponds to the values of Jacobian on the boundary magnetic surface. The equilibrium is for EAST shot 38300 at 3.9s.


pict

Fig. 12: |∇Ψ| as a function of the poloidal angle. The different lines corresponds to values of |∇Ψ| on different magnetic surfaces. The stars correspond to the values of |∇Ψ| on the boundary magnetic surface while the plus signs correspond to the value on the innermost magnetic surface (the magnetic surface adjacent to the magnetic axis). The equilibrium is a Solovev equilibrium.

 


pict

Fig. 13: Jacobian on different magnetic surfaces as a function of the poloidal angle. The equilibrium is a Solovev equilibrium and the Jacobian is an equal-arc Jacobian. The stars correspond to the values of Jacobian on the boundary magnetic surface while the plus signs correspond to the value on the innermost magnetic surface (the magnetic surface adjacent to the magnetic axis).

 

9.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.,

R = R (ψ,𝜃),
(244)

Z = R (ψ,𝜃).
(245)

Then R and Z are written as

      ˆ
∇R  = R = Rψ∇ ψ + R𝜃∇ 𝜃,
(246)

∇Z  = ˆZ = Zψ∇ ψ + Z𝜃∇𝜃,
(247)

wehre Rψ ∂R∕∂ψ, etc. Equations (246) and (247) can be solved to give

           1
∇ψ =  R-Z--−-Z-R--(Z 𝜃Rˆ − R 𝜃ˆZ),
       ψ 𝜃    ψ 𝜃
(248)

     -----1------    ˆ     ˆ
∇𝜃 = ZψR 𝜃 − R ψZ𝜃(ZψR − R ψZ).
(249)

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

  −1
𝒥    = ∇ψ × ∇𝜃 ⋅∇ϕ
     = −-------1-------(Z  R ˆϕ − R Z  ˆϕ)× ϕˆ
        (RψZ 𝜃 − ZψR 𝜃)2 𝜃 ψ     𝜃 ψ     R
               1
     = −(Z𝜃R-ψ −-R-𝜃Z-ψ)R,                                (250)
i.e.,
𝒥 = R (R 𝜃Zψ − R ψZ𝜃).
(251)

Using this, Expressions (248) and (249) are written as

        R
∇ ψ = − 𝒥 (Z 𝜃Rˆ− R 𝜃ˆZ)
(252)

and

∇𝜃 = R-(ZψRˆ− R ψˆZ).
     𝒥
(253)

Then the elements of the metric matrix are written as

        R2
|∇ψ|2 = 𝒥2(Z2𝜃 + R2𝜃),
(254)

|∇ 𝜃|2 = R2-(Z2 +R2 ),
       𝒥 2  ψ    ψ
(255)

and

            2
∇ψ ⋅∇ 𝜃 = − R-(Z 𝜃Z ψ + R 𝜃Rψ).
           𝒥2
(256)

Equations (254), (255), and (256) are the expressions of the metric elements in terms of R, Rψ, R𝜃, Zψ, and Z𝜃. [Combining the above results, we obtain

∇ψ ⋅∇ 𝜃    Z Z  + R R
-----2-= − -𝜃--ψ2---𝜃2-ψ.
 |∇ ψ|         Z𝜃 + R 𝜃
(257)

Equation (256) is used in GTAW code. Using the above results, hαβ = 𝒥-
R2α ⋅∇β are written as

       𝒥         1
hψ ψ =--2|∇ψ|2 = -(Z2𝜃 + R2𝜃)
      R          𝒥
(258)

      𝒥         1
h𝜃𝜃 = -2|∇ 𝜃|2 = -(Z2ψ + R2ψ ),
      R         𝒥
(259)

      𝒥             1
hψ𝜃 = R2-∇ψ ⋅∇ 𝜃 = − 𝒥-(Z 𝜃Z ψ + R 𝜃Rψ)
(260)

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

     ∘ ------------
dℓp =  (dR)2 + (dZ)2
     ∘ -------------2---------------2
   =   (Rψdψ + R𝜃d𝜃) + (Zψdψ + Z𝜃d𝜃) .
Note that = 0 since we are considering the arc length along a constant ψ surface in (R,Z) plane. Then the above expression is reduced to
      ∘----------------
dℓp =  (R-𝜃d𝜃)2 +(Z𝜃d𝜃)2
      ∘  2    2
    =  R 𝜃 + Z𝜃d𝜃
    = |𝒥-∇-ψ|d𝜃,                                (261)
        R
which agrees with Eq. (214).]

9.3 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. (206), i.e.,

         1 g ∫ 2π 𝒥
q(ψ) = −2π-Ψ′     R2d𝜃,
               0
(262)

and re-organize the formula as

              ∫
dΨ-= − -1-g(ψ-)  2π-𝒥-d𝜃,
dψ     2π q(ψ ) 0  R2
(263)

which can be integrated to obtain Ψ and thus the poloidal magnetic field Bp = Ψ ×∇ϕ, where 𝒥 = [(ψ ×∇𝜃) ⋅∇ϕ]1. 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

                          −1
R (r,𝜃) = R0 (r)+ r cos[𝜃+ (sin δ(r))sin𝜃]
Z (r,𝜃) = κ(r)rsin 𝜃
(264)

where r = ψ is the radial coordinate, R0(r), δ(r), and κr(r) are the Shafronov shift, triangularity, and elongation profiles, which can be arbitrarily specified. For the special case of R0 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. ??. 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.

10 Magnetic surface averaging

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

           ( ∫ ∫ ∫      )
                 ΔψGd3V
⟨G ⟩ ≡ Δliψm→0  ∫-∫ ∫--d3V-- ,
                 Δ ψ
(265)

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 definition can be written as

       d ∫
⟨G ⟩ =---   Gd3V,
      dV  V
(266)

where the volume integration is over the volume enclosed by the magnetic surface labeled by ψ.

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

          (                         )
            ∫2π∫ 2π ∫ψ+Δψ G|𝒥 |dψd𝜃dϕ
⟨G ⟩ = Δlψim→0 (-0∫2π0∫2π∫ψψ+Δψ-----------)
              0  0  ψ     |𝒥 |dψd 𝜃dϕ
     ∫ 2π ∫2πG |𝒥 |d𝜃dϕ
   = -0∫2π∫02π---------,                                (267)
       0  0  |𝒥 |d𝜃dϕ
which is a 2D averaging over a magnetic surface. Noting that 𝒥 is independent of ϕ, expression (267) can be written as
      ∫
      -20πG0-|𝒥-|d𝜃-
⟨G⟩ =  ∫2π|𝒥 |d𝜃  ,
        0
(268)

where G0(ψ,𝜃) (2π)1 02πG(ψ,𝜃,ϕ), which is the n = 0 harmonic of G, where n is the toroidal mode number. Note that the surface averaging of any n0 harmonic is zero. I.e., only the n = 0 harmonic can contribute to the magnetic surface average.

In flux coordinates, Eq. (266) can be written as

      dψ d ∫ ψ ∫ 2π ∫ 2π
⟨G⟩ = ------           G|𝒥|dψd𝜃dϕ
      dV d∫ψ2π0∫ 2π0   0
   =  dψ-        G|𝒥|d𝜃dϕ
      dV  0   0
      2π∫ 2π
   =  V′-   G0 |𝒥 |d𝜃,                                (269)
         0
where V = dV∕dψ. (For the surface average, in some parts of these notes, I assume that the Jacobian is postive and thus the absolute value operator is dropped. In most part of this document, axisymmetry is assumed for G, i.e., G = G0.)

The volume within a magnetic surface is written

       ∫ ∫  ∫               ∫ ψ∫ 2π
V (ψ) =       |𝒥 |d𝜃dϕdψ = 2π        |𝒥 |d𝜃dψ.
                             ψ0 0
(270)

Then the differential of V with respect to ψ is written as

     dV     ∫ 2π
V′ ≡ ---= 2π    |𝒥|d𝜃.
     dψ      0
(271)

 

Sometimes, we do not want the Jacobian to explicitly appear in the formula. This can be achieved by using  Eq. (214), i.e.,

|𝒥 |d𝜃 = -R--dl,
        |∇ ψ| p
(272)

Then expression (268) is written as

     ∫ 2π
     -0--G0|R∇ψ|dlp
⟨G ⟩ =  ∫2π-R--dl
     ∫  0 |∇ψ| p
      02π G0B1p dlp
   = --∫2π-1----.                           (273)
        0 Bpdlp
(Equation (273) is used in GTAW code[?].)

10.1 Expression of surface-averaged parallel current density

The parallel current density appears in the flux diffusion equation. Let us derive a formula for it in the flux coordinates. Using μ0J = ∇× B, along with magnetic field expression (200) and the curl formula (162), we obtain

        [ ∂ (   𝒥      )   ∂ (   𝒥        )]
μ0J = −  ∂ψ- Ψ′R2-|∇ ψ|2 +  ∂𝜃- Ψ′R2∇ ψ ⋅∇𝜃   ∇ψ × ∇ 𝜃− g′∇ ϕ ×∇ ψ,   (274)
where Ψ= dΨ∕dψ, g= dg∕dψ. Expression (274) is the contravariant form of J. Taking dot product between Eqs. (200) and (274), we obtain
           [ ∂-( ′-𝒥-    2)   ∂-(  ′ 𝒥-      )]
μ0J ⋅B = −  ∂ψ  Ψ R2 |∇ ψ|  +  ∂𝜃  Ψ R2∇ ψ ⋅∇𝜃   g∇ψ × ∇ 𝜃⋅∇ϕ
            (    𝒥      )
         − g′ − Ψ′-2|∇ ψ|2 ∇ ϕ× ∇ ψ⋅∇ 𝜃
           [   ( R        )     (            )]
       = −  ∂-- Ψ′-𝒥-|∇ ψ|2 +  ∂-- Ψ′ 𝒥-∇ ψ ⋅∇𝜃  g𝒥− 1
            ∂ψ    R2          ∂𝜃    R2
           ′ − 1 ′ 𝒥-   2
         +g 𝒥  Ψ  R2|∇ψ|
           2  −1{ 1 ∂ (  ′ 𝒥    2)   1 ∂ (  ′ 𝒥       )   g′ ′ 𝒥     2}
       = − g 𝒥    g∂ψ- Ψ  R2|∇ψ |  + g∂𝜃- Ψ  R2∇ ψ⋅∇ 𝜃  − g2Ψ R2-|∇ ψ|
                [ ∂ ( Ψ′ 𝒥      )   ∂ ( Ψ′𝒥         )]
       = − g2𝒥 −1---  ----2|∇ ψ|2 +  --- ----2∇ψ ⋅∇ 𝜃  .                (275)
                 ∂ψ   g R           ∂𝜃  g R
Using the defintion of the magnetic surface average, Eq. (269), i.e.,
.= 2π-
V ′ 02π(.)|𝒥|d𝜃

we obtain

            ⟨       d ( Ψ′𝒥      ) ⟩
⟨μ0J ⋅B ⟩ = −  g2𝒥 −1---  ----2|∇ ψ|2
               ∫ 2π dψ   g R   (          )
        = − 2π-    |𝒥 |g2𝒥 −1-d-  Ψ′𝒥--|∇ ψ|2  d𝜃
            V ′ 0          dψ   g R2
            2π 2 d ( Ψ′∫ 2π |𝒥 |   2  )
        = − V-′g dψ-  g-     R2-|∇ ψ|d𝜃
             2   (  ′   0 ∫ 2π          )
        = − g--d-- Ψ-V ′2π-    |𝒥-||∇ψ |2d𝜃
            V ′dψ(  g  V⟨′ 0  ⟩R2)
            g2-d-- Ψ′  ′ |∇ψ-|2
        = − V ′dψ   g V   R2     ,                        (276)
which agrees with Eq. (5) in Ref. [?].

11 Flux Surface Functions

This is an interesting way of defining flux within a magnetic surface using volume integration, instend of surface integration. Examine the meaning of the following volume integral:

          ∫
D (ψ) ≡ 1--  B ⋅∇𝜃dτ,
        2π V
(277)

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

       ∫
D = -1-   ∇ ⋅(𝜃B )dτ.
    2π  V
(278)

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

     1 ∫           1 ∫            1 ∫
D = ---   𝜃B ⋅dS+  ---  𝜃B ⋅dS + ---   𝜃B ⋅dS,
    2π  S1         2π S2         2π  S3
where the orientation of surface S1 is in the negative direction of 𝜃, the orientation of S2 is in the positive direction of 𝜃, and the surface S3 is the toroidal magnetic surface ψ = ψ0. The surface integration through S3 is obviously zero since B lies in this surface. Therefore, we have
       ∫             ∫
    -1-           -1-
D = 2π  S1 𝜃B ⋅dS+ 2π S2 𝜃B ⋅dS + 0
     1 ∫           1 ∫
  = 2π-   0B ⋅dS+ 2π-   2πB ⋅dS
    ∫   S1            S2
  =    B ⋅dS.                                       (279)
     S2
Eq. (279) indicates that D is the magnetic flux through the S2 surface. This is the flux Ψp(axis) defined by Eq. (25), i.e.,
           ∫
 (axis)  -1-
Ψp    = 2π  V B ⋅∇ 𝜃dτ.
(280)

[We can also directly verify Ψp(axis) defined by Eq. (280) is the poloidal flux, without using the divergence theorem. Using the expression of the volume element = |𝒥|d𝜃dϕdψ,

          ∫
Ψ(axis)=  1--  B ⋅∇𝜃|𝒥 |d𝜃dϕdψ
 p      2π  V
       ∫ ψ   ∫ 2π
     =  0  dψ 0  B ⋅∇ 𝜃|𝒥 |d𝜃
       ∫ ψ   ∫ 2π
     =     dψ    Ψ ′∇ ψ ×∇ ϕ ⋅∇𝜃|𝒥|d𝜃
        0     0 ∫     ∫
                  ψ     2π  ′
     = − sign(𝒥 ) 0 dψ  0  Ψ (ψ )d𝜃
                  ∫ ψ
     = − 2πsign (𝒥 )   Ψ ′(ψ)dψ
                   0
     = − 2πsign (𝒥 )[Ψ (ψ )− Ψ(0)].                       (281)
Note that the sign of the Jacobian appears in Eq. (281), which is due to the positive direction of surface S2 is determined by the positive direction of 𝜃, which in turn is determined by the sign of the Jacobian (In my code, however, the positive direction of 𝜃 is chosen by me and the sign of the Jacobian is determined by the positive direction of 𝜃). We can verify the sign of Eq. (281) is consistent with that in Eq. (25).]

 

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

         ∫
Φ(ψ) = 1--  B ⋅∇ϕd τ                         (282)
       2π  V
       1-∫  -g-
    =  2π  V R2 dτ.                           (283)
In flux coordinates, Φ is written as
          ∫
Φ(ψ) = -1-  -g-|𝒥 |d𝜃dϕdψ
       2π   R2
       ∫ ψ   ∫ 2π -g-
     =  0 dψ  0  R2 |𝒥 |d𝜃
       ∫ ψ[   ′⟨   ⟩ ]
     =     gV--  12-  dψ,                       (284)
        0   2π   R
from which, we can also obtain
dΦ    V ′⟨ 1 ⟩
dψ-= g2π-  R2- ,
(285)

and

dΦ   dΦ dψ     1 ⟨  1⟩
dV-= dψ-dV-= g 2π- R2- .
(286)

Using dΦ∕dΨ = 2πq, we obtain

                    ⟨    ⟩
-dΨ = dΨ-dΦ-= -g--1-- -12 .
dV    dΦ dV   2πq 2π  R
(287)

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

       1 ∫
I(ψ) =---   J⋅∇ ϕdτ.
      2π  V
(288)

Let us express I(ψ) in flux coordinates:

       1 ∫
I(ψ) = 2π-  J ⋅∇ϕdτ
         ∫ Vρ ∫ 2π ∫ 2π
     = 1--           J⋅∇ ϕ|𝒥 |dψd𝜃dϕ
       2∫π ∫0  0   0
         ψ  2π
     =  0  0  J ⋅∇ϕ|𝒥|dψd𝜃                           (289)
       ∫ ψV ′
     =    ---⟨J ⋅∇ϕ ⟩dψ.                               (290)
        0 2π
Use Eq. (51), i.e.,
        1   ∗
Jϕ = − μ0R△  Ψ.
(291)

and Eq. (341), i.e.,

         {    [          ]     [              ]}
△⋆Ψ = R2-  -∂- dΨ-𝒥-|∇ ψ|2  + ∂-- dΨ-𝒥-(∇ψ ⋅∇ 𝜃)  ,
       𝒥   ∂ψ  dψ R2         ∂𝜃  dψR2
(292)

then J ⋅∇ϕ is written as

        Jϕ
J⋅∇ ϕ = R--
         1    ∗
      = μ0R2△  Ψ
         1  { ∂  [dΨ 𝒥      ]   ∂ [dΨ  𝒥         ]}
      = ----  --- ----2|∇ψ |2 +  --------2(∇ ψ ⋅∇𝜃)
        μ0𝒥   ∂ψ  dψ R          ∂𝜃 dψ R
Then the surface average is written as
          2π∫ 2π
⟨J ⋅∇ ϕ⟩ = V′-   J ⋅∇ ϕ|𝒥 |d𝜃
            ∫02π   {    [          ]     [              ]}
       =  2π-    1-- ∂-- dΨ-𝒥-|∇ ψ|2  + ∂-- dΨ-𝒥-(∇ψ ⋅∇ 𝜃)  d𝜃
          V′∫0   μ0{ ∂ψ [dψ R2     ]}  ∂𝜃  dψR2
          2π- 2π 1-- ∂-- dΨ-𝒥--   2
       =  V′ 0   μ0  ∂ψ  dψ R2|∇ ψ|   d𝜃
          1 2π ∂  (dΨ ∫ 2π |∇ ψ|2    )
       =  ----′--- ---     --2--𝒥d𝜃
          μ0V  ∂ψ [dψ  0⟨   R  ⟩]
       =  1--1-∂-- dΨ-V′  |∇-ψ|2  .                                (293)
          μ0V ′∂ψ  dψ      R2
Using this in Eq. (290), we obtain
       ∫ ψV ′1  1  ∂ [dΨ   ⟨ |∇ ψ|2⟩]
I(ψ) =    -------′--- ---V′  --2--   dψ
        0 2π∫μ0ψV  ∂ψ[  dψ ⟨    R ⟩]
     = 1--1-    ∂-- dΨV ′  |∇-ψ|2-  dψ
       μ02π  0  ∂ψ  dψ      R2
       1  1 dΨ  ′⟨ |∇ ψ|2⟩
     = μ02π-dψ-V   -R2--  − 0.                         (294)

12 Magnetic flux diffusion equation

12.1 In laboratory static coordinate system

Faraday’s law

∇ × E = − ∂B-,
          ∂t
(295)

can be written in the integral form:

∮
           ∂Ψp-
  E ⋅dl = − ∂t ,
(296)

where Ψp is the flux linked with the loop. Choosing an axisymmetric loop along ϕ, and assuming E is axisymmetric, Eq. (296) is written as

∂-Ψp
 ∂t = 2πREϕ

i.e.,

∂Ψp-
 ∂t = 2πR2E ⋅∇ϕ,

Using ∂Ψp∕∂t = 2π∂Ψ∕∂t, the above equation is written as

∂Ψ-     2
 ∂t = − R E ⋅∇ ϕ.
(297)

Use Ohm’s law

E + u × B = R,

(where R stands for all the terms on the rhs of Ohm’s law, typically R = η(J Jni), where J is the total current density and Jni is the non-inductive part), to eliminate E, then Eq. (297) is written as

∂Ψ
---= − R2 (R − u × B) ⋅∇ϕ.
∂t
(298)

Using B = Ψ ×∇ϕ + gϕ to eliminate B on the rhs, we obtain

∂ Ψ
--- = − R2R ⋅∇ϕ + R2((∇ Ψ × ∇ϕ) ×∇ ϕ)⋅u.
 ∂t
(299)

Noting that Ψ ⋅∇ϕ = 0, the above equation is simplified to

∂Ψ-= − R2R ⋅∇ ϕ− u ⋅∇Ψ.
 ∂t
(300)

This is an evolution equation for the poloidal flux function Ψ in the lab frame.

Next, let us derive an evolution equation for the toroidal field function g BϕR. The ϕ projection of Faraday’s law (295) is written as

∂-
∂t(B ⋅∇ϕ ) = − ∇ ϕ⋅∇ × E,
(301)

i.e.,

∂ ( g )
∂t  R2- = ∇ ⋅(∇ϕ × E).
(302)

Use Ohm’s law to eliminate E, then Eq. (302) is written as

   (   )
-∂  -g- = ∇ ⋅(∇ ϕ× (R − u ×B )),
∂t  R2
(303)

i.e.,

∂ ( g )      (                     g  )
∂t  R2- = ∇ ⋅ ∇ϕ × R + (u⋅∇ ϕ)B − R2u  .
(304)

12.2 Transform to time-dependent magnetic coordinates

Next, let us write Eqs. (300) and (304) in a general magnetic coordinate system (ψ,𝜃,ϕ) that is time-dependent (relative to the laboratory static coordinate system). Why do we need a magnetic coordinate system? Because it makes it easy to perform the magnetic surface average when we derive 1D transport equations. Then why do we require the magnetic coordinates to be time-dependent? Because the magnetic surface shapes are evolving, and thus a static transform can not remain conform to the evolving shapes.

For a time-dependent coordinate transform:

x(ψ,𝜃,ϕ,t),

where x stands for a laboratory static Cartesian coordinate system (x,y,z), it follows that

  ||         ||          ||
∂f||    =  ∂f|| + ∂f-⋅ ∂x||
∂t ψ,𝜃,ϕ   ∂t|x  ∂x   ∂t ψ,𝜃,ϕ
          ∂f||
       =  ∂t|x + uc ⋅∇f,                          (305)
where
     ∂x ||
uc ≡ ∂t-||   ,
        ψ,𝜃,ϕ
(306)

is the coordinate velocity. Further define

ur ≡ u − uc,
(307)

which is the fluid velocity observed in (ψ,𝜃,ϕ) coordinate system.

Note that the flux evolution equation (300) is written in a laboratory static coordinate system, i.e.,

   |
∂-Ψ||      2
 ∂t|x = − R R ⋅∇ ϕ− u ⋅∇Ψ.
(308)

Then in time-dependent (ψ,𝜃,ϕ) coordinate system, the above equation is written as

∂Ψ||                  2
∂t||    − uc ⋅∇Ψ = − R R ⋅∇ϕ − u ⋅∇Ψ,
   ψ,𝜃,ϕ
(309)

Noting u = uc + ur, the above equation is simplified to

  |
∂Ψ||         2
∂t|ψ,𝜃,ϕ = − R R ⋅∇ ϕ− ur ⋅∇Ψ.
(310)

Now we restrict (ψ,𝜃,ϕ) to be flux coordinates (i.e., contours of ψ are contours of Ψ) and remain flux coordinates when both ψ and Ψ evolve. Then the rhs of Eq. (310) must be chosen to be independent of 𝜃, i.e.,

   2
− R R ⋅∇ ϕ− ur ⋅∇Ψ = f(ψ,t),
(311)

where f(ψ,t) is a function to be determined. For later use, multiplying the above equation by g∕R2:

           g           g
− R ⋅g∇ ϕ − R2ur ⋅∇Ψ = R2-f(ψ),
(312)

and performing the magnetic surface average on it, we obtain

            ⟨         ⟩   ⟨       ⟩
− ⟨R ⋅g∇ϕ ⟩− -g-ur ⋅∇ Ψ =  -g-f(ψ) .
             R2            R2
(313)

Next, consider the evolution equation for g, Eq. (304), which can be written as

       |
∂-( g-)||      (                           -g- )
∂t  R2 |x = ∇ ⋅ ∇ϕ × R + ((ur +uc) ⋅∇ϕ)B − R2 u ,
(314)

We note that uc ⋅∇ϕ = 0. Then the above equation is written as

∂ ( g )||      (                      g  )
∂t R2- || = ∇ ⋅ ∇ ϕ ×R  +(ur ⋅∇ϕ)B − R2-u ,
        x
(315)

To facilitate later use, we multiply the above equation by the Jacobian J of (ψ,𝜃,ϕ) coordinates:

𝒥       |
∂-( g-)||
∂t  R2 |x = 𝒥∇⋅(                     g-- )
 ∇ ϕ× R + (ur ⋅∇ ϕ)B − R2u,

i.e.,

∂ (   g )||
-- 𝒥 --2 ||
∂t   Rx g
-2-
R∂𝒥 ||
---||
∂tx = 𝒥∇⋅(                      g  )
 ∇ ϕ ×R  +(ur ⋅∇ϕ)B − --2u
                      R.

Tansforming to time-dependent (ψ,𝜃,ϕ) coordinate system, the above equation is written as

                                   (                  )
   ∂ (  g )||           (   g )   g   ∂𝒥 ||
  ∂t  𝒥 R2-||    − uc ⋅∇ 𝒥 R2- −  R2- -∂t||    − uc ⋅∇ 𝒥 .
      (     ψ,𝜃,ϕ                )        ψ,𝜃,ϕ
= 𝒥∇ ⋅  ∇ϕ × R + (ur ⋅∇ϕ)B − -g-u  .                            (316)
                            R2
For notation ease, the subscripts (ψ,𝜃,ϕ) of the time derivative is dropped in the following, with the understanding that the time derivatives are taken with (ψ,𝜃,ϕ) held fixed. Equation (316) simplifies to
  ∂ (   g )         ( g )    g ∂𝒥
  ∂t 𝒥 R2- − uc ⋅𝒥∇   R2- − R2-∂t-
      (                     -g- )
= 𝒥∇ ⋅ ∇ ϕ × R + (ur ⋅∇ϕ)B − R2 u .
Noting that u = uc + ur, the above equation is written as
   ∂ (  g )         ( g )   g ∂ 𝒥
  --  𝒥 -2- − uc ⋅𝒥 ∇ --2 −  -2---
  ∂t   (R            R      Rg ∂t)       ( g   )
= 𝒥 ∇ ⋅ ∇ϕ × R + (ur ⋅∇ ϕ)B − R2ur  − 𝒥 ∇ ⋅ R2-uc
which simplifies to
    (     )
  ∂-  𝒥-g-  − g--∂𝒥-
  ∂t   R(2     R2 ∂t              )
= 𝒥 ∇ ⋅ ∇ϕ × R + (ur ⋅∇ϕ )B − -g2ur − 𝒥 -g2∇ ⋅uc.            (317)
                            R         R
Use the well known relation (to be proved):
∂-𝒥 = 𝒥 ∇ ⋅u ,
 ∂t         c
(318)

then Eq. (317) simplifies to

  (     )       (                          )
∂- 𝒥 -g- = 𝒥 ∇ ⋅ ∇ ϕ× R + (ur ⋅∇ ϕ)B −-gur
∂t   R2                               R2
(319)

Integrating Eq. (319) over 𝜃:

02π∂
∂t(   g )
 𝒥 R2-d𝜃 = 02π𝒥∇⋅(                      g   )
 ∇ ϕ ×R + (ur ⋅∇ϕ)B − R2-urd𝜃,

we obtain

   (∫          )  ∫
-∂    2π  -g-       2π     (                     -g-  )
∂t   0  𝒥 R2 d𝜃 =   0  𝒥∇ ⋅ ∇ ϕ× R + (ur ⋅∇ ϕ)B − R2 ur d𝜃
(320)

Using the defintion of the surface average, Eq. (269), i.e.,

.= 2π
-′-
V 02π(.)𝒥d𝜃

(where V = dV∕dψ, V is the volume enclosed by surface ψ), Eq. (320) is written as

 1 ∂             V′ ⟨   (                     g   )⟩
-----(V ′g⟨R −2⟩) = --- ∇ ⋅ ∇ϕ × R + (ur ⋅∇ ϕ)B −-2ur   .
2π ∂t            2π                           R
(321)

Noting the relation (proof is given in Sec. ??):

         1--∂-  ′
⟨∇ ⋅A⟩ = V′∂ψ (V ⟨A  ⋅∇ ψ⟩),
(322)

Eq. (321) is written as

-1-
2π∂-
∂t(V gR2) = V-′
2π1--
V′-∂-
∂ψ(  ′⟨(                     g-- )     ⟩)
 V    ∇ ϕ× R + (ur ⋅∇ ϕ)B − R2ur ⋅∇ ψ

Noting that B ⋅∇ψ = 0, the above equation is simplified to

 ∂             ∂  [  ⟨                g       ⟩]
∂t(V ′g⟨R− 2⟩) = ∂ψ- V′ (∇ ψ× ∇ ϕ)⋅R − R2-ur ⋅∇ ψ .
(323)

12.3 Choice of coordinates

The toroidal magnetic flux is given by Eq. (284), i.e.,

Φ = 0ψ[   ′⟨   ⟩]
 gV-- -1-
  2π  R2dψ.

Then the time partial derivative in (ψ,𝜃,ϕ,t) coordinates is written as

        ( ∫  [    ⟨   ⟩ ]  )
∂-Φ = ∂-    ψ gV-′  1--  dψ
 ∂t   ∂t   0   2π   R2
        ∫ ψ   [    ⟨   ⟩]
    = 1--   ∂- gV ′ -1-   dψ.                     (324)
      2π  0 ∂t      R2
Use Eq. (323), then Eq. (324) is written as
        ∫
∂Φ-  -1-  ψ ∂--[ ′⟨               -g-      ⟩]
∂t = 2π  0  ∂ψ V   (∇ψ × ∇ ϕ)⋅R − R2 ur ⋅∇ ψ dψ
      1   ′⟨               g        ⟩||ψ= ψ
   = 2π-V   (∇ ψ × ∇ϕ) ⋅R − R2ur ⋅∇ψ  |ψ=0
      1   ⟨                g       ⟩
   = 2π-V′ (∇ ψ× ∇ ϕ)⋅R − R2-ur ⋅∇ ψ ,                     (325)
where use has been made of the fact that ψ = 0 at ψ = 0.

If we choose ψ as the toroidal magnetic flux, i.e., ψ = Φ (or more general ψ = h(Φ), where h(.) is an arbitrary time-independent monotonical function), then, by defintion

∂ Φ
--- = 0.
 ∂t
(326)

Combining this with Eqs. (325), we get

                ⟨ g       ⟩
⟨(∇ ψ ×∇ ϕ)⋅R ⟩−  --2ur ⋅∇ ψ = 0.
                 R
(327)

This is a condition that the coordinate velocity uc (and thus ur) satisfies when you choose ψ as h(Φ).

Meanwhile, as is discussed above, the requirement of (ψ,𝜃,ϕ) being flux coordinates imposes a constraint given by Eq. (311). Next, we combine this constraint with Eq. (327) to uniquely determine f(ψ,t) appearing in Eq. (311). Mutiplying Eq. (327) by Ψand then adding it to Eq. (313), we obtain

−⟨(Ψ ×∇ϕ) R+ ⟨          ⟩
  g-ur ⋅∇Ψ
  R2−⟨R gϕ⟩−⟨         ⟩
 -g-ur ⋅∇ Ψ
 R2 = ⟨         ⟩
  g-f(ψ,t)
  R2

which simplifies to (ur terms cancel out)

−⟨B R= ⟨ g      ⟩
 R2-f(ψ, t),

i.e.,

          ⟨B ⋅R⟩
f(ψ,t) = −---−2-.
          g⟨R  ⟩
(328)

Then the poloidal magnetic flux diffusion equation is written as

∂Ψ     ⟨B ⋅R ⟩
∂t-= − g⟨R−2⟩.
(329)

Note that it is in the moving flux coordinates (Φ,𝜃,ϕ) that the poloidal magnetic flux diffusion can be written in the form given by Eq. (329). Before I did the above derivation, I considered Eq. (329) (which appears in many papers) as an approximate equation in static lab frame (approximating the toroidal component by the parallel component). It turns out that it is an exact equation in the moving coordinate system.

We can switch the role of Ψ and Φ and do similar derivation as the above. I.e., we choose ψ = h2(Ψ) and derive an evolution equation for the toroidal flux Φ. (It would be interesting to do this derivation. to be done.) But most researchers have chosen ψ = h(Φ). Some reasons behind this choice are as follows. Like the poloidal flux, the toroidal flux is contributed by two parts: coils and plasma. Different from the poloidal flux, the toroidal flux is mainly contributed by the coils. This is because the plasma poloidal current (which contributes to the toroidal flux) is much smaller than the TF coil currents. Because TF coil currents usually do not change during one discharge, the toroidal flux Φ(R,Z,t) remains almost constant in a discharge (small variation is contributed by the varying plasma poloidal current, which is much smaller than the TF coil current). We usually prefer to our coordinate not moving significantly with respect to the lab frame. So we prefer ψ = h(Φ) over ψ = h2(Ψ).

A confusion I once had: if you choose ψ = h(Φ), then ψ can also be considered as a function of Ψ because Φ is a function of Ψ. Then why does it matter which one is chosen in the above derivation? The reason why we can not easily switch between Φ and Ψ is that thier relation is time-dependent, i.e., Φ = Φ(Ψ,t). This time-dependence will show up when we want to switch from one to another.

Assume R = η(J Jni), then Eq. (329) is written as

∂Ψ-    η⟨B-⋅J⟩    ⟨B-⋅Jni⟩
 ∂t = − g⟨R −2⟩ + η g⟨R−2⟩ .
(330)

Using Eq. (276), the above equation is written as

                  (      ⟨     ⟩ )
∂Ψ-=  ---η----g-d-- Ψ′V ′ |∇ψ-|2    +η ⟨B-⋅Jni⟩,
∂t    μ0⟨R −2⟩V ′dψ   g     R2         g⟨R− 2⟩
(331)

Define ρ by

    ∘ -----
      -Φ---
ρ =   πBϕ0,
(332)

where Bϕ0 is typical value of the toroidal field (a space-time constant). Since ρ is a time-independent function of Φ, so ρ can be used as the radial coordinate, i.e., ψ = ρ. Then Eq. (331) is written as

                  (  ′⟨    2⟩    )
∂Ψ-= ---η−-2 g′-∂- V--  |∇-ρ2|- ∂-Ψ  + η⟨B-⋅J−n2i⟩,
∂t   μ0⟨R  ⟩ V ∂ρ   g    R     ∂ρ      g⟨R   ⟩
(333)

which agrees with Eq. (8) in Ref. [?]. The radial derivative of volume is written as

                                    2
V′ = ∂V = ∂V-∂Φ-=  -2π---2πρBϕ0 = 4π-ρBϕ0,
     ∂ρ   ∂Φ ∂ρ    g⟨R −2⟩         g⟨R −2⟩
(334)

where use has been made of Eq. (286), i.e., dΦ∕dV = gR2(2π).

 

12.4 Boundary conditions for poloidal magnetic flux diffusion equation

 

∂Ψ (Φ,t)   ∂Ψ ∂Φ    1           ρ
------- = ------= ----2ρπBϕ0 = -Bϕ0
   ∂ρ     ∂Φ ∂ρ   2πq          q
(335)

We know that Φ = 0 and ρ = 0 at the magnetic axis. We also know that the safety factor q is nonzero at the magnetic axis. Then Eq. (335) implies

  |
∂Ψ||   = 0,
∂ρ|ρ=0
(336)

which serves as an inner boundary for the diffusion equation (333).

Using Eq. (294), the outer boundary condition can be written as

   |
∂Ψ-||    = 2πμ0I(ρ = ρb)∕G(ρ = ρb),
 ∂ρ|ρ= ρb
(337)

where

G = V ⟨     ⟩
  |∇-ρ|2-
   R2.

13 Fixed boundary equilibrium problem

The fixed boundary equilibrium problem (also called the “inverse equilibrium problem” by some authors) refers to the case where the shape of a boundary magnetic surface is given and one is asked to solve the equilibrium within this magnetic surface. To make it convenient to deal with the shape of the boundary, one usually uses a general coordinates system which has one coordinate surface coinciding with the given magnetic surface. This makes it trivial to deal with the irregular boundary. To obtain the equilibrium, one needs to solve the GS equation in the general coordinate system.

13.1 Toroidal elliptic operator in general coordinates

Next we derive the form of the toroidal elliptic operator in a general curvilinear coordinate system. Specifically, we consider a curvilinear coordinate system (ψ,𝜃,ϕ), which is an arbitrary coordinate system except that ϕ is perpendicular to both ψ and 𝜃.

The toroidal elliptic operator takes the form given by Eq. (52), i.e.,

            ( 1    )
△ ∗Ψ ≡ R2∇ ⋅ R2-∇ Ψ  .
(338)

The gradient of Ψ is written as (note that Ψ is independent of ϕ)

     ∂ Ψ      ∂Ψ
∇Ψ = ---∇ ψ + --∇ 𝜃.
      ∂ψ      ∂𝜃
(339)

Using this expression and the divergence formula (156), the elliptic operator in Eq. (338) is written

            (  1 ∂Ψ       1 ∂Ψ   )
△ ∗Ψ = R2∇ ⋅  -2---∇ ψ + -2---∇ 𝜃
              R( ∂ψ      R  ∂𝜃                 )
     = R2 1-∂-- 𝒥 -1-∂Ψ∇ ψ ⋅∇ψ + 𝒥 1-∂-Ψ∇ 𝜃⋅∇ ψ
          𝒥 ∂ψ(   R2 ∂ψ            R2 ∂𝜃      )
        2 1-∂--   1-∂Ψ-           -1-∂Ψ-
     + R  𝒥 ∂𝜃  𝒥 R2∂ψ ∇ ψ⋅∇ 𝜃+ 𝒥 R2 ∂𝜃∇ 𝜃⋅∇ 𝜃
         2[(          )    (            )  ]
     = R--   Ψψ 𝒥-|∇ ψ|2  +   Ψ𝜃 𝒥-∇ψ ⋅∇𝜃
        𝒥      R2       ψ      R2         ψ
       R2 [(   𝒥        )    (   𝒥      ) ]
     + -𝒥-  Ψ ψR2-∇ψ ⋅∇𝜃   +  Ψ 𝜃R2|∇𝜃|2   .                (340)
                          𝜃              𝜃
where the subscripts denotes partial derivatives, 𝒥 is the Jacobian of the coordinate system (ψ,𝜃,ϕ).

In a flux coordinate system, Ψ is independent of 𝜃, i.e., Ψ𝜃 = 0. Then the toroidal elliptic operator is written as

          {                                }
  ⋆    R2  [   𝒥     2]    [  𝒥          ]
△  Ψ = 𝒥--  Ψψ R2|∇ψ |   +  ΨψR2-(∇ ψ ⋅∇𝜃)    ,
                       ψ                  𝜃
(341)

Using Eq. (340), the GS equation (65) is written

  2[(           )    (            ) ]
 R--  Ψψ-𝒥-|∇ψ|2   +  Ψ𝜃-𝒥-∇ψ ⋅∇ 𝜃
 𝒥      R2       ψ      R2         ψ
  R2 [(   𝒥        )    (   𝒥     )  ]
+ ---  Ψψ -2∇ ψ ⋅∇𝜃   +  Ψ𝜃--2|∇ 𝜃|2
  𝒥       R          𝜃     R        𝜃
= − μ0R2dP-−  dgg,                                    (342)
        dΨ    dΨ
which is the form of the GS equation in (ψ,𝜃,ϕ) coordinate system.

13.2 Finite difference form of toroidal elliptic operator in general coordinate system

The toroidal elliptic operator in Eq. (340) can be written

       R2
△ ∗Ψ = ---[(hψψΨψ)ψ + (h 𝜃𝜃Ψ𝜃)𝜃 + (hψ𝜃Ψ𝜃)ψ + (hψ𝜃Ψψ)𝜃],
       𝒥
(343)

where h is defined by Eq. (197), i.e.,

hαβ = 𝒥-∇ α⋅∇ β.
      R2
(344)

Next, we derive the finite difference form of the toroidal elliptic operator. The finite difference form of the term (hψψΨψ)ψ is written

                [       (           )         (           ) ]
(hψψΨψ)ψ|i,j =-1- hψψ      Ψi,j+1-−-Ψi,j- − hψψ     Ψi,j −-Ψi,j−1
             δψ   i,j+1∕2      δψ          i,j− 1∕2      δψ
           = H ψψ   (Ψi,j+1 − Ψi,j)− H ψψ  (Ψi,j − Ψi,j−1),           (345)
               i,j+1∕2                i,j− 1∕2
where
       hψψ
Hψψ = ----2.
      (δψ)
(346)

The finite difference form of (h𝜃𝜃Ψ𝜃)𝜃 is written

              [       (           )          (           )]
(h𝜃𝜃Ψ𝜃)𝜃|i,j = 1- h 𝜃𝜃i+1∕2,j  Ψi+1,j −-Ψi,j − h𝜃i𝜃−1∕2,j  Ψi,j-−-Ψi−1,j
            δ𝜃              δ𝜃                    δ𝜃
          = Hi𝜃𝜃+1∕2,j(Ψi+1,j − Ψi,j)− H𝜃i𝜃−1∕2,j(Ψi,j − Ψi−1,j),             (347)
where
H 𝜃𝜃 =-h𝜃𝜃-.
      (δ𝜃)2
(348)

The finite difference form of (hψ𝜃Ψ𝜃)ψ is written as

        |     1 [       (Ψi+1,j+1∕2 − Ψi−1,j+1∕2 )        ( Ψi+1,j− 1∕2 − Ψi−1,j−1∕2)]
(Ψ𝜃hψ𝜃)ψ|i,j =--- hψi𝜃,j+1∕2  -------------------- − hψi𝜃,j−1∕2  --------------------  .
             δψ                  2δ𝜃                             2δ𝜃
(349)

Approximating the value of Ψ at the grid centers by the average of the value of Ψ at the neighbor grid points, Eq. (349) is written as

    ψ𝜃        ψ𝜃                                    ψ𝜃
(Ψ𝜃h  )ψ|i,j = H i,j+1∕2(Ψi+1,j+Ψi+1,j+1− Ψi−1,j− Ψi−1,j+1)− H i,j−1∕2(Ψi+1,j− 1+ Ψi+1,j− Ψi−1,j−1− Ψi− 1,j).
(350)

where

  ψ𝜃   -hψ𝜃--
H    = 4δψd𝜃.
(351)

Similarly, the finite difference form of (hψ𝜃Ψψ)𝜃 is written as

               [       (                    )         (                    )]
(Ψ hψ𝜃)|  =  1- hψ𝜃     Ψi+1∕2,j+1 −-Ψi+1∕2,j−1 − hψ𝜃      Ψi−-1∕2,j+1 −-Ψi−1∕2,j−1
  ψ    𝜃i,j   δ𝜃   i+1∕2,j         2δψ             i−1∕2,j          2δψ
          = H ψ𝜃   (Ψ      + Ψ    − Ψ       − Ψ   ) − H ψ𝜃  (Ψ    + Ψ       − Ψ    − Ψ  (352).)
              i+1∕2,j  i+1,j+1   i,j+1   i+1,j−1    i,j−1     i−1∕2,j  i,j+1   i−1,j+1    i,j−1    i− 1,j− 1
Using the above results, the finite difference form of the operator 𝒥△Ψ∕R2 is written as
 𝒥     ||
R2-△∗Ψ ||  = (h ψ𝜃Ψ 𝜃)ψ + (h ψψΨψ)ψ + (h𝜃𝜃Ψ𝜃)𝜃 + (h ψ𝜃Ψ ψ)𝜃
       i,j
          = Hψiψ,j+1∕2(Ψi,j+1 − Ψi,j) − Hψiψ,j−1∕2(Ψi,j − Ψi,j−1)
          + H𝜃𝜃    (Ψ     − Ψ ) − H𝜃𝜃   (Ψ   − Ψ   )
             i+1∕2,j i+1,j    i,j    i−1∕2,j i,j    i−1,j
          + Hψi𝜃,j+1∕2(Ψi+1,j + Ψi+1,j+1 − Ψi− 1,j − Ψi−1,j+1)− Hiψ,j𝜃−1∕2(Ψi+1,j−1 + Ψi+1,j − Ψi− 1,j−1 − Ψi−1,j)
             ψ𝜃                                        ψ𝜃
          + Hi+1∕2,j(Ψi+1,j+1 +Ψi,j+1 − Ψi+1,j− 1 − Ψi,j−1)− Hi−1∕2,j(Ψi,j+1 + Ψi− 1,j+1 − Ψi− 1,j−1 − Ψi,j−1)
          = Ψi,j(− Hiψ,jψ+1∕2 − H ψi,ψj− 1∕2 − H 𝜃i+𝜃1∕2,j − H 𝜃i𝜃−1∕2,j)
                     ψ𝜃       ψ𝜃              ψψ       ψ𝜃       ψ𝜃
          + Ψi− 1,j−1(H i,j−1∕2 +H i−1∕2,j) +Ψi,j−1(H i,j− 1∕2 − H i+1∕2,j +H i−1∕2,j)
          + Ψi+1,j− 1(− H ψ𝜃   − H ψ𝜃   )+ Ψi−1,j(H𝜃𝜃    − Hψ 𝜃   + H ψ𝜃   )
                      i,j−1∕2ψ𝜃  i+1∕2,jψ𝜃        i−1∕2,j   iψ,j𝜃+1∕2   i,jψ−𝜃1∕2
          + Ψi+1,j(H 𝜃i+𝜃1∕2,j + Hi,j+1∕2 − Hi,j−1∕2)+ Ψi− 1,j+1(− Hi,j+1∕2 − Hi−1∕2,j)
          + Ψ    (H ψψ    + Hψ𝜃    − Hψ𝜃    )+ Ψ     (H ψ𝜃   + H ψ𝜃   )
             i,j+1  i,j+1∕2    i+1∕2,j   i−1∕2,j    i+1,j+1   i,j+1∕2    i+1∕2,j
The coefficients are given by
       𝒥         1
hψψ = --2|∇ ψ|2 =--(R2𝜃 + Z2𝜃),
      R         𝒥
(353)

h𝜃𝜃 = 𝒥-|∇ 𝜃|2 = 1(R2 + Z2 ),
      R2        𝒥   ψ   ψ
(354)

and

hψ𝜃 = 𝒥-∇ ψ ⋅∇𝜃 = −-1(R R  + Z Z  ),
      R2           𝒥   𝜃  ψ    𝜃 ψ
(355)

where the Jacobian

𝒥 = R (R 𝜃Zψ − R ψZ𝜃).
(356)

The partial derivatives, R𝜃, Rψ, Z𝜃, and Zψ, appearing in Eqs. (353)-(356) are calculated by using the central difference scheme. The values of hψψ, h𝜃𝜃, hψ𝜃 and 𝒥 at the middle points are approximated by the linear average of their values at the neighbor grid points.

 

13.3 Pressure and toroidal field function profile

The function p(Ψ) and g(Ψ) in the GS equation are free functions which must be specified by users before solving the GS equation. Next, we discuss one way to specify the free functions. Following Ref. [?],  we take P(Ψ) and g(Ψ) to be of the forms

                     --
P(Ψ) = P0 − (P0 − Pb)ˆp(Ψ),
(357)

1       1         --
2g2(Ψ) = 2 g20(1 − γ ˆg(Ψ )),
(358)

with ˆp and ĝ chosen to be of polynomial form:

  --  --α
ˆp(Ψ) = Ψ ,
(359)

  --  --β
ˆg(Ψ) = Ψ ,
(360)

where

--  Ψ − Ψa
Ψ ≡ Ψ-−-Ψ--,
     b    a
(361)

with Ψb the value of Ψ on the boundary, Ψa the value of Ψ on the magnetic axis, α, β, γ, P0, Pb, and g0 are free parameters. Using the profiles of P and g given by Eqs. (357) and (358), we obtain

                     --
dP-(Ψ-)= − 1-(P0 − Pb)αΨα−1
  dΨ      Δ
(362)

where Δ = Ψb Ψa, and

g dg-= − 1g2γ 1-βΨβ−1.
  dΨ     2 0 Δ
(363)

Then the term on the r.h.s (nonlinear source term) of the GS equation is written

                      [            --  ]          --
− μ0R2dP- − dgg = μ0R2  1-(P0 − Pb)α Ψα−1 +  1g20γ 1-βΨβ−1.
       dΨ   dΨ          Δ                  2   Δ
(364)

The value of parameters P0, Pb, and g0 in Eqs. (357) and (358), and the value of α and β in Eqs. (359) and (360) can be chosen arbitrarily. The parameter γ is used to set the value of the total toroidal current. The toroidal current density is given by Eq. (??), i.e.,

J  = RdP- + 1-dg-g-,
 ϕ    dΨ    μ0dΨ R
(365)

which can be integrated over the poloidal cross section within the boundary magnetic surface to give the total toroidal current,

    ∫
Iϕ =   JϕdS
    ∫  (              )
          dP-  1--dgg-
  =     R dΨ + μ0dΨ R   dRdZ
    ∫  [ (   1)           --    1  1   1   --]
  =     R  −--  (P0 − Pb)ˆp′(Ψ) −-----g20γ--ˆg′(Ψ)  dRdZ          (366)
            Δ                  μ0R 2   Δ
Using
        1
dRdZ = R-𝒥 dψd𝜃,
(367)

Eq. (366) is written as

    ∫ [(    )                               ]
          -1           ′--   -1--1  2 1- ′--
Iϕ =      −Δ   (P0 − Pb)ˆp(Ψ)− μ0R22 g0γ Δ ˆg(Ψ) 𝒥 dψd𝜃,
(368)

from which we solve for γ, giving

                    ∫   --
    −-ΔIϕ-− (P0-−-Pb)-[pˆ′(Ψ)]𝒥-dψd𝜃
γ =     -1g2 ∫ [-12ˆg′(Ψ)]𝒥dψd 𝜃  .
        2μ00   R
(369)

If the total toroidal current Iϕ is given, Eq. (369) can be used to determine the value of γ.

13.4 Boundary magnetic surface and initial coordinates

In the fixed boundary equilibrium problem, the shape of the boundary magnetic surface (it is also the boundary of the computational region) is given while the shape of the inner flux surface is to be solved. A simple analytical expression for a D-shaped magnetic surface takes the form

R = R0 + acos(𝜃 +Δ sin𝜃),
(370)

Z = κa sin𝜃,
(371)

with 𝜃 changing from 0 to 2π. According to the definition in Eqs. (??), (??), and (??) we can readily verify that the parameters a, R0, κ appearing in Eqs. (370) and (371) are indeed the minor radius, major radius, and ellipticity, respectively. According to the definition of triangularity Eq. (??), the triangularity δ for the magnetic surface defined by Eqs. (370) and (371) is written as

δ = sin Δ.
(372)

Another common expression for the shape of a magnetic surface was given by Miller[??], which is written as

R = R0 + acos[𝜃+ arcsin(Δ sin𝜃)],
(373)

Z = κa sin𝜃.
(374)

Note that Miller’s formula is only slightly different from the formula (370). For Miller’s formula, it is easy to prove that the triangularity δ is equal to Δ (instead of δ = sinΔ as given in Eq. (372)).

In the iterative metric method[?] for solving the fixed boundary equilibrium problem, we need to provide an initial guess of the shape of the inner flux surface (this initial guess is used to construct a initial generalized coordinates system). A common guess of the inner flux surfaces is given by

R = ψα(RLCFS − R0)+ R0,
(375)

Z = ψαZLCFS,
(376)

where α is a parameter, ψ is a label parameter of flux surface. If the shape of the LCFS is given by Eqs. (370) and (371), then Eqs. (375) and (376) are written as

R = R0 + ψ αacos(𝜃 + Δ sin𝜃),
(377)

      a
Z = ψ κa sin 𝜃.
(378)

Fig. 14 plots the shape given by Eqs. (377) and (378) for a =0.4, R0 = 1.7, κ = 1.7, Δ = arcsin(0.6), and α = 1 with ψ varying from zero to one.


pict pict

Fig. 14: Shape of flux surface given by Eqs. (377) and (378). Left figure: points with the same value of ψ are connected to show ψ coordinate surface; Right figure: dot plot. Parameters are a =0.4m, R0 = 1.7m, κ = 1.7, Δ = arcsin(0.6), α = 1 with ψ varying from zero (center) to one (boundary). The shape parameters of LCFS are chosen according to the parameters of EAST tokamak.

13.5 Fixed boundary equilibrium numerical code

The tokamak equilibrium problem where the shape of the LCFS is given is called fixed boundary equilibrium problem. I wrote a numerical code that uses the iterative metric method[?] to solve this kind of equilibrium problem. Figure 15 describes the steps involved in the iterative metric method.


pict pict

pict pict

pict

Fig. 15: Upper left figure plots the initial coordinate surfaces. After solving the GS equation to obtain the location of the magnetic axis, I shift the origin point of the initial coordinate system to the location of the magnetic axis (upper right figure). Then, reshape the coordinate surface so that the coordinate surfaces ψ = const lies on magnetic surfaces (middle left figure). Recalculate the radial coordinate ψ that is consistent with the Jacobian constraint and interpolate flux surface to uniform ψ coordinates (middle right figure). Recalculate the poloidal coordinate 𝜃 that is consistent with the Jacobian constraint and interpolate poloidal points on every flux surface to uniform 𝜃 coordinates (bottom left figure).

13.6 Benchmark of the code

To benchmark the numerical code, I set the profile of P and g according to Eqs. (??) and (??) with the parameters c2 = 0, c1 = B0(κ02 + 1)(R02κ0q0), κ0 = 1.5, and q0 = 1.5. The comparison of the analytic and numerical results are shown in Fig. 16.


pict

Fig. 16: Agreement between the magnetic surfaces (contours of Ψ) given by the Solovev analytic formula and those calculated by the numerical code. The two sets of magnetic surfaces are indistinguishable at this scale. The boundary magnetic surface is selected by the requirement that Ψ = 0.11022 in the Solovev formula (??). Parameters: κ0 = 1.5, q0 = 1.5.

Note that the parameter c0 in the Solovev equilibrium seems to be not needed in the numerical calculation. In fact this impression is wrong: the c0 parameter is actually needed in determining the boundary magnetic surface of the numerical equilibrium (in the case considered here c0 is chosen as c0 = B0(R02κ0q0)).

13.7 Low-beta equilibrium vs. high-beta equilibrium

With the pressure increasing, the magnetic axis usually shifts to the low-field-side of the device, as is shown in Fig. 17.


pict pict

Fig. 17: Comparison of two equilibria obtained with P0 = 104 Pasca (left) and P0 = 105 Pasca (right), respectively, where P0 is the pressure at the magnetic axis. All the other parameters are the same for the two equilibria, α = 1.0, β = 1.0, Pb = 101 Pasca, g0 = 1.0Tm, Iϕ = 500kA, and the LCFS is given by miller’s formulas (373) and (374) with R0 = 1.7, a = 0.45, κ = 1.7, and δ = 0.6.

 

 

 

13.8 Grad-Shafranov equation with prescribed safety factor profile (to be finished)

In the GS equation, g(Ψ) is one of the two free functions which can be prescribed by users. In some cases, we want to specify the safety factor profile q(Ψ), instead of g(Ψ), in solving the GS equation. Next, we derive the form of the GS equation that contains q(Ψ), rather than g(Ψ), as a free function. The safety factor defined in Eq. (206) can be written

        1  g ∫ 2π𝒥
q(ψ) = − 2πΨ-′   R2-d𝜃                             (379)
             ∫02π  −2    ∫ 2π
    = − 1--g -0∫R---𝒥d𝜃-    𝒥d𝜃
        2πΨ ′  02π𝒥 d𝜃   0
        1  g      ∫ 2π
    = − ----′⟨R −2⟩    𝒥d𝜃                          (380)
        2πΨ        0
    = − sgn(𝒥 )-V-′-⟨R−2⟩g.                         (381)
              (2π)2 Ψ ′
Equation (381) gives the relation between the safety factor q and the toroidal field function g. This relation can be used in the GS equation to eliminate g in favor of q, which gives
            (2π)2
g = − sgn(𝒥 )-′-−2-Ψ ′q.
           V ⟨R  ⟩
(382)

                 [           ]
  -dg    -(2π)2-- --(2π-)2--  ′
⇒ dΨ g = V′⟨R −2⟩q V ′⟨R −2⟩qΨ  ψ.
(383)

Note that this expression involves the flux surface average, which depends on the flux surface shape and the shape is unknown before Ψ is determined.

Multiplying Eq. (??) by R2 gives

  [                                 ]
1  (   𝒥      )    (   𝒥         )        dp       dg
𝒥-   Ψ′R2|∇ψ |2   +  Ψ ′R2(∇ψ ⋅∇ 𝜃)   + μ0dΨ-+ R −2dΨ-g = 0.
                ψ                  𝜃
(384)

Surface-averaging the above equation, we obtain

       [                                ]
2π∫     (  ′ 𝒥    2)    ( ′ 𝒥         )       dp     −2 dg
V′-  d𝜃  Ψ R2-|∇ψ|    +  Ψ R2-(∇ψ ⋅∇𝜃)    + μ0dΨ-+ ⟨R  ⟩dΨ-g = 0,
                    ψ                  𝜃
(385)

     ∫    (          )
⇒  2π-- d𝜃  Ψ′-𝒥-|∇ ψ|2   +μ0 dp-+ ⟨R −2⟩dgg = 0,
   V′        R2       ψ     dΨ        dΨ
(386)

     [       (       ) ]
   2π-  ′∫      |∇ψ-|2        -dp     −2 dg-
⇒  V′ Ψ    d𝜃 𝒥   R2    ψ + μ0dΨ + ⟨R  ⟩dΨg = 0,
(387)

     [     ⟨     ⟩ ]
   1--  ′ ′ |∇ψ-|2        -dp     −2 dg-
⇒  V′ V Ψ     R2    ψ + μ0dΨ + ⟨R  ⟩dΨg = 0.
(388)

   [    ⟨     2⟩]
⇒   V′Ψ′  |∇-ψ2|-    + μ0V′ dp-+ ⟨R−2⟩V′ dg-g = 0,
           R     ψ       dΨ          dΨ
(389)

Substitute Eq. (383) into the above equation to eliminate gdg∕dΨ, we obtain

  [    ⟨      ⟩ ]                  [       ]
     ′ ′ |∇ψ|2          ′dp-       4 --qΨ-′--
⇒  V Ψ    R2     ψ + μ0V dΨ + q(2π) V ′⟨R− 2⟩ ψ = 0,
(390)

Eq. (390) agrees with Eq. (5.55) in Ref. [?].

     ⟨     ⟩     [  ⟨      ⟩]                   [        ]          [       ]
⇒ V ′ |∇ψ-|2  Ψ′′+  V′  |∇-ψ|2    Ψ′+μ V ′ dp+q (2π)4 ---q----Ψ ′′+q (2π)4  ---q----  Ψ′ = 0
        R2             R2    ψ     0   dΨ         V′⟨R −2⟩            V′⟨R−2⟩ ψ
(391)

              {                                              }
    ′′     1   [  ′⟨|∇ψ |2⟩ ]   ′      4[   q    ]   ′     ′ dp
⇒  Ψ  = −V-′D   V   --R2-    Ψ  +q(2π)  V-′⟨R-−2⟩  Ψ + μ0V  dΨ-
                            ψ                    ψ
(392)

where

    ⟨|∇ ψ|2⟩    (2π )4 ( q )2
D =  ---2-  + ---−2  --′
       R      ⟨R   ⟩ V
(393)

The GS equation is

△ ⋆Ψ = − μ R2 dp-− g dg
         0   dΨ     dΨ
(394)

                          4 [     ′  ]
⇒  △ ⋆Ψ = − μ0R2-dp−  q(2π)-----qΨ----
               dΨ    V′⟨R −2⟩ V ′⟨R −2⟩ ψ
(395)

                         4 { [       ]      [        ]  }
⇒ △ ⋆Ψ = − μ0R2-dp− q(2π)---  ---q----  Ψ′ + ---q---- Ψ′′
              dΨ    V′⟨R −2⟩   V ′⟨R−2⟩ ψ      V ′⟨R− 2⟩
(396)

Using Eq. (392) to eliminate Ψ′′ in the above equation, the coefficients before (μ0dp∕dΨ) is written as

                { [       ]       }
B = R2 − -q(2π)4--  ---q---- -1--V′
         V′⟨R −2⟩   V′⟨R−2⟩  V′D
    1 {  2    ( q )2 (2π)4 }
  = D-  R D −  V-′  ⟨R−-2⟩2  .                        (397)
Substituting the expression of D into the above equation, we obtain
     1 { 2 ⟨|∇ψ |2⟩     2 (2π)4 ( q )2 ( q )2 (2π)4}
B = D-  R   --R2-  + R  ⟨R-−2⟩  V′- −   V′- ⟨R-−2⟩2
       {   ⟨    2⟩   (   )2    4 (          ) }
  = -1  R2  |∇ψ2|  +  -q′  -(2π−)2  R2 − --1−2-   ,             (398)
    D         R       V    ⟨R   ⟩      ⟨R  ⟩
which is equal to the expression (5.58) in Ref. [?]. The coefficient before Ψis written as
A =  q(2π)4
--′--−2-
V ⟨R   ⟩{[    q   ]   [    q   ]  1  [(   ⟨|∇ψ |2⟩ )         [    q   ] ]}
   -′--−2-- −   -′--−2----′-   V ′ ---2-     + q(2π)4  -′--−2--
   V ⟨R   ⟩ ψ    V ⟨R   ⟩ V  D         R     ψ          V ⟨R   ⟩ ψ

Define

α = ⟨     ⟩
 |∇ψ-|2
   R2,

β = ---q----
V ′⟨R− 2⟩.

A = (2π)4β{                           }
 β  − β--1-[(V′α) + q(2π)4β  ]
  ψ    V ′D       ψ        ψ

      D                1   ′          4
=⇒ −-(2π)4β-A = Dβψ − βV-′[(V α)ψ +q(2π) βψ]
(399)

= ⇒ ---D----A = D β − β-1[V′′α + V ′α  + q(2π )4β ]
    − β(2π)4       ψ    V ′         ψ         ψ
(400)

Using

D = α + (2π )4⟨R− 2⟩β2
(401)

Eq. () is written as

                                                                    ]
=⇒  ---D---A = αβ  + (2π)4⟨R −2⟩β2β  − β-1-V′′α − β-1V ′α  − β 1-q(2π )4β
    − β(2π)4     ψ               ψ    V ′       V ′   ψ    V′      ψ
(402)

       D                 4  −2  2      1  ′′           1     4  ]
=⇒  −-β(2π-)4A = α βψ + (2π) ⟨R  ⟩β βψ − β V′V α − βαψ − βV-′q(2π) βψ
(403)

   A                                1
−-β(2π)4D = αβψ + (2π )4⟨R− 2⟩β2βψ − βV-′V ′′α − β αψ − β2⟨R −2⟩(2π)4βψ
(404)

       A               1   ′′
= ⇒ −-β(2π)4-D = αβψ − βV′V  α− βα ψ
(405)

                   (   )
=⇒  ---A---D  = − β2 α-   − β 1-V ′′α
    − β(2π)4          β  ψ    V′
(406)

      A            (α )       1   α
=⇒  − β(2π)4D = − β2 β   − β2V-′V ′′β-
                       ψ
(407)

                   [               ]
      A             (α )     1   α
=⇒ −-β(2π)4D = − β2  -β   + V-′V ′′β-
                        ψ
(408)

                      [(  )           ]
    ---A----      2-1-   α-    ′   ′′α-
= ⇒ − β(2π)4 D = − β V ′  β  ψV  + V β
(409)

       A            1 (α   )
=⇒  − β(2π)4D = − β2V-′ β-V′
                            ψ
(410)

           4 (        )3   [⟨     2⟩   −2    ]
=⇒  A = (2π)- ----q---   1--  |∇-ψ|-  ⟨R--⟩V ′2
         D    V ′⟨R −2⟩   V′    R2      q      ψ
(411)

But the expression of A is slightly different from that given in Ref. [?] [Eq. (5.57)]. Using the above coefficients, the GS equation with the q-profile held fixed is written as

(          )
   ⋆    -∂-           dp-
  Δ  − A∂ψ   Ψ = − B μ0dΨ.
(412)

 

14 Magnetic coordinates with general toroidal angle

14.1 General toroidal angle ζ

In Sec. 8.12, we introduced the local safety factor ˆq (ψ,𝜃). Equation (202) indicates that if the Jacobian is chosen to be of the particular form 𝒥 = h(ψ)R2, 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(ψ)R2. We note that, as mentioned in Sec. 8.14, the poloidal angle is fully determined by the choice of the Jacobian. The specific choice of 𝒥 = α(ψ)R2 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 obvious: Yes, we can achieve this by defining a new toroidal angle ζ that generalizes the usual (cylindrical) toroidal angle ϕ. Define a new toroidal angle ζ by[?]

ζ = ϕ− q(ψ)δ(ψ,𝜃),
(413)

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

ˆqnew ≡ B-⋅∇-ζ
      B ⋅∇ 𝜃
    = (∇-ϕ×-∇-ψ+-ˆq∇-ψ-×∇-𝜃)⋅∇-ζ
      (∇ ϕ× ∇ ψ+ ˆq∇ ψ ×∇ 𝜃)⋅∇ 𝜃
      ∇-ϕ×-∇-ψ-⋅∇ζ-+qˆ∇-ψ-×-∇𝜃-⋅∇ζ
    =        (∇ϕ × ∇ψ )⋅∇𝜃
      ∇ ϕ× ∇ ψ ⋅∇(− qδ)+ ˆq∇ ψ ×∇ 𝜃⋅∇ ϕ
    = ---------(∇ϕ-×-∇ψ-)⋅∇𝜃---------

    = − q∂δ+ ˆq.                                       (414)
         ∂𝜃
To make the new local safety factor be independent of 𝜃, the right-hand side of Eq. (414) should be independent of 𝜃, i.e.,
− q∂δ-+qˆ= c(ψ ),
   ∂𝜃
(415)

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., ˆq new = q. In this case, equation (415) is written as

∂δ-= ˆq− 1,
∂𝜃   q
(416)

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

                            1 ∫ 𝜃
δ(ψ,𝜃) = δ(ψ,𝜃ref)− (𝜃− 𝜃ref)+-     ˆqd𝜃,
                            q  𝜃ref
(417)

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

         1∫ 𝜃
δ = − 𝜃 + q 0 qˆd𝜃.
(418)

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

           ∫ 𝜃
ζ = ϕ + q𝜃−   ˆqd𝜃,
            0
(419)

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

In summary, magnetic field line is straight in (𝜃,ζ) plane with slope being q if ζ is defined by Eq. (419). 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:

𝒥−ne1w = ∇ ψ × ∇𝜃⋅∇ ζ
    = ∇ ψ × ∇𝜃⋅∇ (ϕ− qδ)

    = (∇ ψ× ∇ 𝜃)⋅∇ ϕ− (∇ψ × ∇ 𝜃) ⋅∇(qδ)
    = ∇ ψ × ∇𝜃⋅∇ ϕ − 0
    = 𝒥 −1.
] Also note that δ(ψ,𝜃) defined by Eq. (418) is a periodic function of 𝜃. [Proof: Equation (418) implies that
              ∫
             1  𝜃+2π
δ(ψ,𝜃+ 2π) = q 0    ˆqd𝜃− (𝜃− 𝜃ref)− 2π
             1∫ 𝜃                    1∫ 𝜃+2π
          =  -   ˆqd𝜃 − (𝜃 − 𝜃ref)− 2π +      qˆd𝜃
             q 0          ∫          q  𝜃
          = δ(ψ,𝜃)− 2π + 1  𝜃+2πˆqd𝜃
                         q 𝜃
          = δ(ψ,𝜃).                                        (420)
]

[In numerical implementation, the term 0𝜃B-⋅∇ϕ
B ⋅∇𝜃d𝜃 appearing in δ is computed by using

∫             ∫           ∫
  𝜃B-⋅∇-ϕ      𝜃 -g-𝒥-      𝜃 g--1 -R--
 0 B ⋅∇ 𝜃d𝜃 =  0 R2 Ψ′d𝜃 = 0  R2Ψ ′|∇ ψ|dlp
              ∫ 𝜃 g 1
            =    ------dlp
              ∫0 𝜃 R |∇ Ψ|
            =    1-Bϕdℓ .                               (421)
               0 R Bp  p

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

                 ∫               ∫
∂(δq)    -d-(-g )  𝜃-𝒥-    -g ∂--  𝜃𝒥--    dq-
 ∂ψ   = −dψ  Ψ ′  0 R2 d𝜃− Ψ ′∂ψ  0 R2 d𝜃− dψ 𝜃.           (422)
This formula is used in GTAW code, where the derivative (g∕Ψ)∕∂ψ is calculated numerically by using the central difference scheme.]

 

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

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

B = − Ψ′(∇ϕ × ∇ψ + ˆq∇ψ × ∇ 𝜃).
(423)

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

       ′                ′
B = − Ψ ∇′ (ζ + qδ)× ∇ψ′ − Ψ ˆq∇ψ × ∇ 𝜃′
  = − Ψ ∇ ζ × ∇ ψ − Ψ ∇ (qδ) × ∇ψ − Ψ ˆq∇ψ × ∇𝜃
       ′           ′∂-δ           ′
  = − Ψ ∇ ζ × ∇ ψ − Ψ q∂ 𝜃∇𝜃 × ∇ψ − Ψ ˆq∇ψ × ∇𝜃           (424)
Using Eq. (416), the above equation is simplified as
B  = − Ψ′(∇ζ × ∇ψ + q∇ψ × ∇ 𝜃).                    (425)
Equation (425) is the contravariant form of the magnetic field in (ψ,𝜃,ζ) coordinates.

The expression of the magnetic field in Eq. (425) can be rewritten in terms of the flux function Ψp and Ψt discussed in Sec. (). Equation (425) is

B = ∇ Ψ × ∇ ζ + q∇ 𝜃× ∇ Ψ,
(426)

which, by using Eq. (), i.e., Ψ = Ψp(2π), is rewritten as

B = -1-(∇ζ ×∇ Ψ + q∇ Ψ × ∇ 𝜃),
    2π         p      p
(427)

which, by using Eq. (), i.e., q = dΦ∕dΨp, is further written as

B = -1-(∇ ζ × ∇ Ψ + ∇ Φ ×∇ 𝜃).
    2π         p
(428)

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

Noting the simple fact that

-d-=  --d---,
dx    d(x + c)
(429)

where c is a constant, we conclude that

(   )     (    )
  ∂f-       ∂f-
  ∂ζ  ψ,𝜃 =  ∂ϕ  ψ,𝜃 ,
(430)

(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,

( ∂f )     ( ∂f)
  ∂ψ-    ⁄=  ∂-ψ
      𝜃,ζ         𝜃,ϕ
(431)

and

(   )     (    )
  ∂f-       ∂f-
  ∂𝜃 ψ,ζ ⁄=  ∂𝜃  ψ,ϕ .
(432)

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

14.4 Steps to construct a straight-line magnetic coordinate system

In Sec. 9, 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. (413), where δ is obtained from Eq. (418). Also note that the Jacobian of (ψ,𝜃,ζ) coordinates happens to be equal to that of (ψ,𝜃,ϕ) coordinates.

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

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

B ⋅∇f  = − Ψ ′(∇ζ × ∇ψ )⋅∇f(ψ,𝜃,ζ)− Ψ ′q(∇ψ × ∇ 𝜃)⋅∇f(ψ,𝜃,ζ)
                ( ∂     ∂ )
       = − Ψ ′𝒥 −1 --+ q---  f.                                  (433)
                  ∂𝜃   ∂ ζ
Next, consider the solution of the following magnetic differential equation:
B ⋅∇f = h.
(434)

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

( ∂        ∂)       1
  --+ q(ψ)--- f = − -′𝒥h(ψ,𝜃,ζ).
  ∂𝜃      ∂ζ        Ψ
(435)

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

          ∑         i(m𝜃−nζ)
f(ψ,𝜃,ζ) =   fmn(ψ)e       ,
          m,n
(436)

(note that, following the convention adopted in tokamak literature[?], the Fourier harmonics are chosen to be ei(m𝜃), instead of ei(m𝜃+)), and the right-hand side is expanded as

  1             ∑
− -′𝒥 h(ψ, 𝜃,ζ) =    γmn(ψ)ei(m𝜃−nζ),
  Ψ             m,n
(437)

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

         γ
fmn = ----mn--.
      i[m − nq]
(438)

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.

14.6 Resonant surface of a perturbation

Equation (438) indicates that, for the differential equation (434), there is a resonant response to a perturbation ei(m𝜃) on a magnetic surface with m nq = 0. Therefore, the magnetic surface with q = m∕n is called the “resonant surface” for the perturbation ei(m𝜃).

The phase change of the perturbation ei(m𝜃) along a magnetic field is given by mΔ𝜃 nΔζ, which can be written as Δ𝜃(mnq). Since mnq = 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.

14.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

        n
χ = 𝜃−  mζ,
(439)

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𝜃 ) reduce to 2D perturbations, i.e.,

  |
∂f||   = 0.
∂ζ|ψ,χ
(440)

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

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

B(χ) ≡ B ⋅∇ χ
    = − Ψ′(∇ ζ × ∇ ψ + q∇ ψ × ∇𝜃)⋅∇ (𝜃− nζ∕m )
         ′                 ′ n
    = − Ψ (∇ ζ × ∇ ψ ⋅∇𝜃)+ Ψ m-(q∇ψ × ∇ 𝜃⋅∇ ζ)
         ′ −1    ′ n  −1
    = − Ψ 𝒥   + Ψ m q𝒥
    = Ψ ′𝒥 −1( nq− 1) .                                   (441)
              m
At the resonant surface q = m∕n, equation (441) implies B(χ) = 0. The direction χ defines the reconnecting component of the magnetic field?

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

B(𝜃) ≡ B ⋅∇𝜃
    = − Ψ ′(∇ ζ × ∇ ψ+ q∇ ψ× ∇ 𝜃)⋅∇ 𝜃
    = − Ψ ′∇ ζ × ∇ ψ⋅∇ 𝜃
         ′ −1
    = − Ψ 𝒥  .                                      (442)
Using (442) and (441), the relation between B(𝜃) and B(χ) is written as
          (    nq)
B (χ) = B (𝜃) 1 − m- .
(443)

14.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. (200)). 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

g∇ ϕ = g∇ (ζ + qδ)
    = g∇ ζ + gq∇ δ+ gδ∇q
               (∂δ      ∂δ   )      ′
    = g∇ ζ + gq ∂ψ-∇ψ + ∂𝜃-∇𝜃  + gδq∇ ψ
      (           )
    =   gq∂δ-+ gδq′ ∇ ψ+ gq∂δ-∇𝜃 + g∇ ζ
          ∂ψ               ∂𝜃
    = g ∂(qδ)∇ ψ + gq ∂δ∇ 𝜃+ g∇ζ.                        (444)
         ∂ψ         ∂𝜃
Using Eq. (444), the covariant form of the magnetic field, Eq. (200), is written as
     (  𝒥            ∂(qδ))      (  ∂δ     𝒥      )
B =   Ψ′R2-∇ψ ⋅∇𝜃 + g-∂ψ-- ∇ ψ +  gq∂𝜃-− Ψ′R2|∇ψ |2  ∇𝜃 + g∇ζ.
(445)

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

     (  𝒥            ∂(qδ))      (  g2    R2          ) 𝒥
B  =  Ψ′--2∇ψ ⋅∇ 𝜃+ g----- ∇ ψ +  − -′ − gq--− Ψ ′|∇ψ |2  -2∇ 𝜃+ g∇ ζ
     (  R             ∂ψ  )      (  Ψ      𝒥        )   R
   =  Ψ′ 𝒥-∇ψ ⋅∇ 𝜃+ g∂(qδ) ∇ ψ +  − g2 +-|∇-Ψ|2− gqR2  𝒥-∇ 𝜃+ g∇ ζ.   (446)
        R2            ∂ψ               Ψ ′        𝒥   R2
Using B2 = (|∇Ψ|2 + g2)∕R2, the above equation is written as
    (                    )      (              )
B =  Ψ′-𝒥-∇ψ ⋅∇ 𝜃+ g∂(qδ) ∇ ψ +  − B2R2-− gqR2-  𝒥-∇ 𝜃+ g∇ ζ
       R2            ∂ψ             Ψ ′      𝒥   R2
    ( ′ 𝒥           ∂(qδ))      (  B2      )
  =  Ψ R2-∇ψ ⋅∇ 𝜃+ g-∂ψ-- ∇ ψ +  − Ψ′ 𝒥 − gq ∇ 𝜃+ g∇ζ.           (447)
Equation (447) is the covariant form of the magnetic field in (ψ,𝜃,ζ) coordinate system. For the particular choice of the radial coordinate ψ = Ψ and the Jacobian 𝒥 = h(ψ)∕B2 (i.e., Boozer’s Jacobian, discussed in Sec. 14.9), equation (447) reduces to
    (   𝒥           ∂(qδ))
B =   −--2∇ψ ⋅∇ 𝜃+ g----- ∇ ψ + I(ψ )∇ 𝜃+ g(ψ)∇ ζ,
       R             ∂ψ
(448)

with I(ψ) = h(ψ) gq. The magnetic field expression in Eq. (448) frequently appears in tokamak literature[?]. 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

𝒥  = gq+-I.
      B2
(449)

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

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

B × ∇ ψ    1 (  B2      )            g
----2-- = -2- − --′ 𝒥 − gq ∇ 𝜃× ∇ ψ+ --2∇ ζ × ∇ ψ.
  B       B     Ψ                   B
(450)

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

                (         )
B-×-∇-ψ ⋅∇ = -1- B2-𝒥 + gq  𝒥− 1 ∂-+-g-𝒥 −1 ∂-            (451)
  B2         B2   Ψ′            ∂ζ  B2     ∂𝜃
             ( 1    𝒥 −1 ) ∂     𝒥−1 ∂
           =  Ψ-′ + g-B2-q ∂ζ-+g B2--∂𝜃,                  (452)
which is the form of the operator in (ψ,𝜃,ζ) coordinate system.

Examining Eq. (452), 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(ψ)∕B2, 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(ψ)∕B2 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 ×∇ψ∕B2) ⋅∇ be a constant-coefficient differential operator. For other choices of the Jacobian, only the B⋅∇ operator is a constant-coefficient differential operator.

14.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 = ∂f
---
∂ψψ + ∂f
---
∂ 𝜃𝜃 + ∂f
---
∂ζζ,

the radial differential operator is written as

∇ ψ ⋅∇f = |∇ ψ|2∂f-+ (∇𝜃 ⋅∇ ψ)∂f-+ (∇ζ ⋅∇ψ )∂f-
               ∂ψ           ∂ 𝜃           ∂ζ
              2∂f-          ∂f-                     ∂f-
        = |∇ ψ| ∂ψ + (∇𝜃 ⋅∇ ψ)∂ 𝜃 + {∇ [ϕ − qδ(ψ,𝜃)]⋅∇ ψ}∂ζ
              2∂f           ∂f            ∂f
        = |∇ ψ| ∂ψ-+ (∇𝜃 ⋅∇ ψ)∂-𝜃 − ∇ [qδ]⋅∇ ψ∂ζ
               ∂f           ∂f                  ∂f
        = |∇ ψ|2---+ (∇𝜃 ⋅∇ ψ)---− [q∇δ +δ∇q ]⋅∇ψ ---
               ∂ψ           ∂ 𝜃  [ (            ∂ζ)        ]
        = |∇ ψ|2∂f-+ (∇𝜃 ⋅∇ ψ)∂f-−  q  ∂δ-∇ψ + ∂δ∇ 𝜃  + δq′∇ ψ ⋅∇ ψ∂f-
               ∂ψ           ∂ 𝜃  [   ∂ψ      ∂𝜃         ]        ∂ζ
              2∂f-          ∂f-   ∂-(qδ)    2   ∂δ-        ∂f-
        = |∇ ψ| ∂ψ + (∇𝜃 ⋅∇ ψ)∂ 𝜃 −  ∂ ψ |∇ψ| + q∂𝜃 ∇𝜃 ⋅∇ψ  ∂ζ,        (453)
where ()∕∂ψ and q∂δ∕∂𝜃 are given respectively by Eqs. (422) and (416). Using the above formula, ψ ⋅∇ζ is written as
            [∂(qδ)   2    ∂δ       ]
∇ ψ ⋅∇ζ = −  -∂ψ--|∇ ψ| + q∂𝜃∇ 𝜃⋅∇ ψ .
(454)

This formula is used in GTAW code.

15 Field-line-following coordinates

15.1 Definition

In (ψ,𝜃,ζ) coordinates, a magnetic field line is straight in (𝜃,ζ) plane with the slope being q. So the equation for magnetic field lines is given by ζ = q𝜃 + α, where α is a constant in (𝜃,ζ) plane. Note that α can be used to label different magnetic field lines on a magnetic surface. This motivates us to use α, i.e.,

α ≡ ζ − q𝜃,
(455)

to as a new coordinate to replace ζ. Then we obtain a new set of coordinate: (ψ,𝜃,α), which are called field-line-following coordinates. Using Eqs. (455) and (419), α can be written as

       ∫ 𝜃
α = ϕ −   ˆqd𝜃,                            (456)
        0
where ˆq = B ⋅∇ϕ∕B ⋅∇𝜃 is the local safety factor. (If we choose the straight-field-line 𝜃, then α is written as α = ϕ q𝜃.) Define δ = 0𝜃qˆ d𝜃, which is called tor_shift in TEK code, then α = ϕ δ.

In TEK, I choose 𝜃 [π,π) with 𝜃 = ±π corresponding to the high-field-side midplane, and 𝜃 is increasing along the counter-clockwise direction wehn viewed along ϕ.  There are two reasons why 𝜃 is chosen this way: 1. The toroidal shift δ remain small near the low-field-side 𝜃 = 0. This may be good for spatial resolution in this important region where ballooning modes often have larger amplitude. 2. The 𝜃 cut, i.e., 𝜃 = ±π, is far away from the important low-field-side. The gyrokinetic field equations are not solved at 𝜃 = +π. Perturbations there need to obtained by using interpolations at the 𝜃 cut, i.e., infer values of perturbations on α gridpoints at 𝜃 = +π from the corresponding values at 𝜃 = π. Numerical errors more likely appear there. So we prefer that the 𝜃 cut is located in less important area (area where mode amplitude is small).