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:

(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

(2)

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

(3)

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

1.1 Poloidal magnetic field

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

(4)

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

(5)

(6)

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

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

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

1.2 Toroidal magnetic field

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

(9)

1.3 General form of axisymmetric magnetic field

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

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:

(11)

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

(12)

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

(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

(15)

i.e.,

(16)

(17)

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

(18)

i.e.,

(19)

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

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

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

(20)

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

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.

The innermost magnetic surface reduces to a curve, which is called magnetic axis (in Fig. 2, Ψ₀ 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 B_p = ∇Ψ ×∇ϕ, the condition Ψ_LCFS − Ψ_axis > 0 means B_p points in the anticlockwise direction (viewed along ϕ direction), and Ψ_LCFS − Ψ_axis < 0 means B_p points in the clockwise direction.

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.

Fig. 3: The poloidal magnetic flux Ψ_p⁽¹²⁾ between the two magnetic surfaces Ψ₁ and Ψ₂ is given by Ψ_p⁽¹²⁾ = 2π(Ψ₂ − Ψ₁), where the surface orientation for the poloidal flux definition is chosen to be +.

In Fig. 3, there are two magnetic surfaces labeled, respectively, by Ψ = Ψ₁ and Ψ = Ψ₂. The poloidal magnetic flux through any toroidal ribbons between the two magnetic surfaces is equal to each other (due to the magnetic Gauss theorem). Denote this poloidal magnetic flux by Ψ_p⁽¹²⁾.  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 B_Z contribute to the poloidal magnetic flux: (the positive direction of the plane is chosen to be )

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

(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

(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 S₀ in Fig. 2. In this case, as is discussed above, the poloidal magnetic flux is related to Ψ by

(24)

where the postive direction is +. 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 S₂. Denote this flux by Ψ_p^(axis), then

(25)

where Ψ_axis is the value of Ψ at the magnetic axis, the S₂ 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

(26)

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

(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

(28)

Integrating the above equation over time, we obtain

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

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.

(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

(31)

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

(32)

which can be integrated over dℓ_p to give

(33)

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

(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

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

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

(36)

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

(37)

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

(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

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

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

(41)

(42)

(43)

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

(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

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

3.1 Poloidal current density

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

(47)

and

(48)

respectively.

3.2 Toroidal current density

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

(50)

then Eq. (49) is written as

(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

(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

(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

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

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

(56)

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

(57)

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

(58)

which can be further written

(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 direction. The component of Eq. (54) is written

(60)

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

(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

(62)

which can be simplified to

(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

(64)

i.e.,

(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 = ⇒− −△∗Ψ = μ0 ⇒− −△∗Ψ = μ0 ⇒− −△∗Ψ = μ0 ⇒△∗Ψ = −μ0R2 −g

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

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

(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, ⟨j_∥⟩ is obtained from current drive models and neoclassical bootstrap current models. Note that the specification of the source terms (P, g, q, and ⟨j_∥⟩) usually involve the unknown Ψ (via not only the explicit presence of Ψ, but also the flux-surface averaging which 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:

(68)

where G is the fundamental solution given by

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

(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 I_i, 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[?]):

(71)

(72)

where

(73)

Ψ_A and Ψ_B are values of Ψ at the magnetic axis and LCFS, the coefficients α_j and β_j are to be determined, N_p and N_F are integers chosen by users. 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 N_p = N_F = 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 I_i by solving a least square problem that minimises the difference between quantities computed and the corresponding quantities measured in actual experiments. (Details are given in Sec. 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 I_i 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

(75)

and the corresponding expansion coefficients

(76)

then Eq. (74) is written as

(77)

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

(79)

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

(80)

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

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

Next, we consider the poloidal field measurement:

(81)

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

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

where

(83)

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

(84)

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

(85)

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

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

(86)

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

(87)

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

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

(88)

where

(89)

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

(90)

and using

(91)

equation (90) is written as

(92)

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

(93)

for (j = 0,…,N_b − 1), and zero for (j = N_b,…,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 (c₁,…,c_Nb,I₁,…,I_Nc).

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

or in its discrete form:

Ψij(k+1) = ∑ i′,j′Gi′j′ijdSi′j′ + ∑ 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 I_vc and −I_vc), and are up-down symmetric, located at (R_vc,Z_vc) and (R_vc,−Z_vc), respectively. And set the current I_vc by

Ivc = −gz,

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

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

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

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

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 g′g(Ψ) are chosen by users, often as polymials in Ψ_N with coefficients that are assumed known (which are later obtained in the least-square method to make the resulting solution approximately match some measurements, as is discussed in the previous section).

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

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

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

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

The above procedure is just a minor modification from the pure Green’s function algorithm discussed in Sec. 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.

(94)

Method 1: discretize Eq. (95) as

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

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

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

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 (R_i,Z_i) with i = 1,2,…M. The actions (coil current changes) are chosen to minimise the following residual:

(99)

where

(100)

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

7 Circuit equation

Faraday’s law

(103)

can be written in the integral form:

(104)

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

(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

(106)

where M_jk is the mutual indctance between the No. j coil and No. k coil, M_k(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

(107)

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

(109)

where ℛ_k is the resistence of the coil, I_k 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 I_k 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

(110)

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

(111)

The anti-series connection implies that

(112)

and

(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

(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., M_AB = M_BA), the mutual inductance M_k,p given by expression (107) can also be used to calculate the effective emf induced in the plasma by the coils:

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

(115)

where M_i′j′ij is the mutual inductance between ij filament and i′j′ filament; 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

(116)

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

(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 − Mk,p = Ik(n+1)ℛ k.

Moving all the unknown (coil currents I_j⁽ⁿ⁺¹⁾ and plasma current I_p⁽ⁿ⁺¹⁾) to the left-hand side, we get

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

(119)

i.e.,

(120)

7.4 Reinforcement learning –to be continued

Define the plasma state by the position of the magnetic axis (R_a,Z_a), total plasma current I_p, and the shape of the LCFS. Discretize the LCFS in the cylindrical coordinates by N_s points. Then the plasma state can be represented by 2N_s + 3 real numbers: R_a, Z_a, I_p, and (R_i,Z_i) for i = 1,2,…N_s.

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

Output of the policy network: voltage command to the coils imposed at t_n 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

(121)

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

(122)

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

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

(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

(124)

which can also be written as

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

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

8.3 Orthogonality relation between two sets of basis vectors

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

(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 (,,), then the left-hand side of Eq. (127) can be written as

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

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

The orthogonality relation in Eq. (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∕∂x₁ is orthogonal to ∇x₂ and ∇x₃, thus, ∂r∕∂x₁ can be written as

(129)

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

(130)

which gives

(132)

Similarly, we obtain

(133)

and

(134)

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

(135)

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

(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

(139)

In term of the contravairant basis vectors, A is written

(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

(141)

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

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

(142)

(143)

(144)

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

(145)

(146)

(147)

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

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

(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

(149)

(150)

(151)

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

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

(153)

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

(154)

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

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

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

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

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

(163)

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

(164)

(165)

(166)

Similarly, we write

(167)

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

(168)

(169)

(170)

The same situation applies for the ∇ζ basis vector,

(171)

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

(172)

(173)

(174)

Summarizing the above results in matrix form, we obtain

(175)

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

(176)

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

(177)

(178)

(179)

For the second contravariant basis vector

(180)

(181)

(182)

(183)

For the third contravariant basis vector

(184)

(185)

(186)

(187)

Summarizing these results, we obtain

(188)

where

M = ,

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

(189)

Similarly, Eq. (188) is simplified to

(190)

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

(191)

can be calculated to give

,

where

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

(192)

]

8.10 Covariant/contravariant representation of equilibrium magnetic field

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

(193)

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

(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

(195)

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

(196)

For the convenience of notation, define

(197)

then Eq. (196) is written as

(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

(199)

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

(200)

8.12 Local safety factor

The local safety factor is defined by

(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

(202)

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

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

(203)

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

(204)

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

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

(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 B_p = |∇Ψ|∕R and B_ϕ = g∕R,  the absolute value of q in expression (207) is written

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

(211)

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

(213)

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

(214)

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

—–

(215)

8.15 Calculating poloidal angle

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

(216)

where the curve integration is along the contour Ψ = Ψ_j, x_ref,j is a reference point on the contour, where value of the poloidal angle is chosen as 𝜃_ref,j. The choice of the positive direction of 𝜃 is up to users. Depending on the positive direction chosen, the sign of the Jacobian of the constructed coordinates can have a sign difference from the 𝒥 appearing in Eq. (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

(218)

Then 𝜃 is written as

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

Jacobian for equal-arc-length poloidal angle

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

(221)

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

(222)

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

(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

(224)

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

(225)

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

(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∕∇ψ)dl_p 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:

(227)

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

(228)

The final Jacobian is given by

(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

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

(230)

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

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

(232)

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

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

(234)

Comparison between different types of poloidal angle

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

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

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

Verification of Jacobian

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

(236)

which can be further written as

(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

(238)

where Ψ₀ 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

(239)

Integrating over the toroidal angle, we obtain

(240)

Further integrating over the poloidal angle, we obtain

(241)

i.e.,

(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 𝒥⁻¹ = ∇ψ ×∇𝜃 ⋅∇ϕ 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, Ψ₀, and the value of Ψ on the last closed flux surface (LCFS), Ψ_b, are given in G-file. Using these two values, I construct a 1D array “psival” with value of elements changing uniform from Ψ₀ to Ψ_b. Then I try to find the contours of Ψ with contour level value ranging from Ψ₀ to Ψ_b. This is done in the following way: construct a series of straight line (in the poloidal plane) that starts from the location of the magnetic axis and ends at one of the points on the LCFS. Combine the straight line equation, Z = Z(R), with the interpolating function Ψ(R,Z), we obtain a one variable function h = Ψ(R,Z(R)). Then finding the location where Ψ is equal to a specified value Ψ_i, is reduced to finding the root of the equation Ψ(R,Z(R)) − Ψ_i = 0. Since this is a one variable equation, the root can be easily found by using simple root finding scheme, such as bisection method (bisection method is used in GTAW code). After finding the roots for each value in the array “psival” on each straight lines, the process of finding the contours of Ψ is finished. The contours of Ψ found this way are plotted in Fig. 7.

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

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

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

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

(244)

(245)

Then ∇R and ∇Z are written as

(246)

(247)

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

(248)

(249)

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

(251)

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

(252)

and

(253)

Then the elements of the metric matrix are written as

(254)

(255)

and

(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

(257)

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

(258)

(259)

(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

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

(262)

and re-organize the formula as

(263)

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

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

(264)

where r = ψ is the radial coordinate, R₀(r), δ(r), and κ_r(r) are the Shafronov shift, triangularity, and elongation profiles, which can be arbitrarily specified. For the special case of R₀ being a constant, δ(r) = 0, and κ(r) = 1, this shape reduces to a concentric-circular magnetic field. An example of calculating the poloidal field of the model field is given in Sec. ??. 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

(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

(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 d³V = |𝒥|dψd𝜃dϕ, where 𝒥 is the Jacobian of (ψ,𝜃,ϕ) coordinates. equation (265) is written as

(268)

where G₀(ψ,𝜃) ≡ (2π)⁻¹ ∫ ₀^2πG(ψ,𝜃,ϕ)dϕ, which is the n = 0 harmonic of G, where n is the toroidal mode number. Note that the surface averaging of any n≠0 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

The volume within a magnetic surface is written

(270)

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

(271)

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

(272)

Then expression (268) is written as

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 μ₀J = ∇× B, along with magnetic field expression (200) and the curl formula (162), we obtain

⟨.⟩ = ∫ 02π(.)|𝒥|d𝜃

we obtain

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:

(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

(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 S₁ and S₂, then equation (278) is written as

(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ϕdψ,

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

(285)

and

(286)

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

(287)

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

(288)

Let us express I(ψ) in flux coordinates:

(291)

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

(292)

then J ⋅∇ϕ is written as

12 Magnetic flux diffusion equation

12.1 In laboratory static coordinate system

Faraday’s law

(295)

can be written in the integral form:

(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

= −2πREϕ

i.e.,

= −2πR2E ⋅∇ϕ,

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

(297)

Use Ohm’s law

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

(298)

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

(299)

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

(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

(301)

i.e.,

(302)

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

(303)

i.e.,

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

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

(306)

is the coordinate velocity. Further define

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

(308)

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

(309)

Noting u = u_c + u_r, the above equation is simplified to

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

(311)

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

(312)

and performing the magnetic surface average on it, we obtain

(313)

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

(314)

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

(315)

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

𝒥x = 𝒥∇⋅,

i.e.,

x −x = 𝒥∇⋅.

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

(318)

then Eq. (317) simplifies to

(319)

Integrating Eq. (319) over 𝜃:

∫ 02πd𝜃 = ∫ 02π𝒥∇⋅d𝜃,

we obtain

(320)

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

⟨.⟩ = ∫ 02π(.)𝒥d𝜃

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

(321)

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

(322)

Eq. (321) is written as

(V ′g⟨R−2⟩) =

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

(323)

12.3 Choice of coordinates

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

Φ = ∫ 0ψdψ.

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

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

(326)

Combining this with Eqs. (325), we get

(327)

This is a condition that the coordinate velocity u_c (and thus u_r) 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⟩ + −⟨R ⋅ g∇ϕ⟩− =

which simplifies to (u_r terms cancel out)

−⟨B ⋅ R⟩ = ,

i.e.,

(328)

Then the poloidal magnetic flux diffusion equation is written as

(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 ψ = h₂(Ψ) 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 ψ = h₂(Ψ).

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 − J_ni), then Eq. (329) is written as

(330)

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

(331)

Define ρ by

(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

(333)

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

(334)

where use has been made of Eq. (286), i.e., dΦ∕dV = g⟨R⁻²⟩∕(2π).

12.4 Boundary conditions for poloidal magnetic flux diffusion equation

(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

(336)

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

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

(337)

where

G = V ′.

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

(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

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

(341)

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

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

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

(343)

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

(344)

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

(346)

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

(348)

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

(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

(350)

where

(351)

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

(353)

(354)

and

(355)

where the Jacobian

(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

(357)

(358)

with and ĝ chosen to be of polynomial form:

(359)

(360)

where

(361)

with Ψ_b the value of Ψ on the boundary, Ψ_a the value of Ψ on the magnetic axis, α, β, γ, P₀, P_b, and g₀ are free parameters. Using the profiles of P and g given by Eqs. (357) and (358), we obtain

(362)

where Δ = Ψ_b − Ψ_a, and

(363)

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

(364)

The value of parameters P₀, P_b, and g₀ 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.,

(365)

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

(367)

Eq. (366) is written as

(368)

from which we solve for γ, giving

(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

(370)

(371)

with 𝜃 changing from 0 to 2π. According to the definition in Eqs. (??), (??), and (??) we can readily verify that the parameters a, R₀, κ 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

(372)

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

(373)

(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

(375)

(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

(377)

(378)

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

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.

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 c₂ = 0, c₁ = B₀(κ₀² + 1)∕(R₀²κ₀q₀), κ₀ = 1.5, and q₀ = 1.5. The comparison of the analytic and numerical results are shown in Fig. 16.

Note that the parameter c₀ in the Solovev equilibrium seems to be not needed in the numerical calculation. In fact this impression is wrong: the c₀ parameter is actually needed in determining the boundary magnetic surface of the numerical equilibrium (in the case considered here c₀ is chosen as c₀ = B₀∕(R₀²κ₀q₀)).

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.

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

(382)

(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 R⁻² gives

(384)

Surface-averaging the above equation, we obtain

(385)

(386)

(387)

(388)

(389)

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

(390)

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

(391)

(392)

where

(393)

The GS equation is

(394)

(395)

(396)

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

A = −

Define

α = ,

β = .

A = −(2π)4β