Multiconstrained Variational Problems in Magnetohydrodynamics, I: Equilibrium

Bruce Turkington Alexander Eydeland Alexander Lifschitz

Joel Spruck

Department of Mathematics University of Massachusetts Amherst, MA 01003

September 20, 1990

1 Introduction

2 Variational formulation

In this section we pose the variational problem whose solutions represent axisymmetric equilibrium configurations of a plasma confined in a tokomak or some other toroidal device. The geometry of the configurations we consider is depicted schematically in Figure 1. In the usual cylindrical coordinates the toroidal region D = {x = (r, ϕ, z) ∈ ^{3} : (r, z) ∈ Ω, 0 ≤ ϕ < 2π} is determined by its cross-section Ω, a domain in the half-plane _{+}^{2} = {(r, z) : r > 0}. In D the plasma-vacuum system is governed by the ideal magnetohydrodynamic equilibrium equations

J × B = ∇p, ∇ × B = J, ∇ · B = 0 . (2.1)

The plasma lies within the free-boundary surface S on which the natural interface conditions hold—namely, p = 0 , n · B = 0 on S and n × B is continuous across S, where n is the unit normal on S. On the fixed boundary ∂D, which is taken to be a perfectly conducting shell, the normal component of B vanishes. These equations and all those to follow are expressed in non-dimensional variables. The physical quantities representing magnetic field, current density and plasma pressure are retrieved as B_{0}B, (μ_{0}L_{0})^{−1}B_{0}J, μ_{0}^{−1}B_{0}^{2}p, respectively, in terms of a characteristic length L_{0}, a characteristic field strength B_{0} and a magnetic permeability μ_{0}.

Under the assumption of axisymmetry with respect to the toroidal angle ϕ, the fields in (2.1) admit the representation

B = ∇Ψ × ∇ϕ + f∇ϕ (2.2)

J = ∇f × ∇ϕ + (LΨ)∇ϕ (2.3)

in terms of the scalar unknowns Ψ = Ψ(r, z), f = f(r, z) and the linear elliptic partial differential operator

L: = −r\frac{∂}{∂r} (\frac{1}{r} \frac{∂}{∂r} ) − \frac{∂^{2}}{∂z^{2}} . (2.4)

The standard interpretation [8, 18, 1] of 2πΨ and 2πf as poloidal flux and poloidal current, respectively, is evident upon separating each of the expressions for B and J into a poloidal part (perpendicular to ∇ϕ) and a toroidal part (parallel to ∇ϕ). In view of (2.2) and (2.3), the force balance equation in (2.1) separates into its poloidal and toroidal parts

− r^{−2}f∇f + r^{−2}(LΨ)∇Ψ = ∇p, ∇f × ∇Ψ = 0 . (2.5)

Then the resulting equilibrium relations f = f(Ψ) and p = p(Ψ) reduce (2.5) to the Grad-Shafranov equation [14, 22, 20]

LΨ = f(Ψ)f′(Ψ) + r^{2}p′(Ψ) . (2.6)

Here and throughout the sequel the prime denotes differentiation with respect to Ψ. While (2.6) holds in the plasma, the vacuum field is governed by

∇f = 0, LΨ = r\overline{J}_{ϕ} . (2.7)

where \overline{J}_{ϕ} is the prescribed toroidal current density determined by the external coils. The boundary condition on the shell is normalized to be

Ψ = 0 on ∂Ω . (2.8)

Also, the plasma-vacuum interface S is necessarily a magnetic surface, Ψ = σ_{0} (say), and the interface conditions are simply that p = 0 on S and ∇Ψ , f are continuous across S.

We can make the sign convention that the total plasma current is positive, and so we can characterize the plasma region as the set {Ψ > σ_{0}}, except under certain exceptional circumstances. In the most typical case, the magnetic surface {Ψ = σ_{0}} encloses a single connected component, although our formulation is not restricted to such configurations only. Instead, we can consider the more general case when the plasma region consists of several components whose union is the set {Ψ > σ_{0}}, provided that all of the components are disjoint from the external current coils. However, in the exceptional case when some external coils have (sufficiently large) positive current there may be some components of {Ψ > σ_{0}} that are neighborhoods of those coils. Then it is necessary to exclude those spurious components from the plasma region. This can be accomplished by introducing a limiting subdomain Ω_{0} ⊂ Ω with respect to which the external coils are exterior, and then defining the plasma region to be the (truncated) set {Ψ > σ_{0}} ∩ Ω_{0}. Of course, Ω_{0} must be chosen so that in equilibrium the actual plasma region is situated a finite distance from ∂Ω_{0}. The freedom in the choice of Ω_{0} is therefore balanced by the necessity of satisfying this a posteriori condition.

Since the poloidal magnetic field is generated both by currents induced in the plasma and by external coils in the vacuum region, it is useful to split the total poloidal flux function into the sum Ψ = ψ + \overline{ψ} accordingly. The external poloidal field ∇\overline{ψ} × ∇ϕ is determined by solving

L\overline{ψ} = r\overline{J}_{ϕ} in Ω, \overline{ψ} = 0 on ∂Ω . (2.9)

For instance, \overline{ψ} is a finite sum of Green functions when \overline{J}_{ϕ} = Σ I_{ℓ}δ(r − R_{ℓ})δ(z − Z_{ℓ}), corresponding to a finite collection of elementary coils located at (R_{ℓ}, Z_{ℓ}) carrying currents I_{ℓ}. Hence, \overline{ψ} can be considered as known, and ψ can be used in place of the unknown Ψ. In order that this splitting be consistent, the plasma region, which carries the current density J_{ϕ} = r^{−1}Lψ, is always assumed to be disjoint from the support of \overline{J}_{ϕ}.

The nonuniqueness of equilibrium configurations governed by the above equations resulting from the presence of arbitrary profile functions f(Ψ) and p(Ψ) in (2.6) has been stressed in the literature [17, 10, 12]. This difficulty may be overcome either by specifying these profile functions within some class or by determining them along with the solution Ψ in response to some subsidiary constraints derived from the underlying physical dynamics. The former approach, although it is rather artificial, has been often adopted ([3, 16]) to expedite the solution of the equilibrium problem, thereby reducing it to a tractable free-boundary problem for a single nonlinear elliptic equation. On the other hand, the latter approach has become accepted as the physically correct formulation of the relevant equilibrium problem. Following the basic principles established by Kruskal and Kulsrud [17] (which are valid even without a restriction to axisymmetry), two natural constraint families are distinguished—namely, the flux and mass within each magnetic surface {Ψ = σ} in the plasma. The physical significance of these quantities stems from the fact that they are conserved under evolution governed by ideal magnetohydrodynamics. Moreover, as is also shown in [17], these constraint families characterize the class of admissible functions (ψ, f, p) (or equivalently the trial states B, p) in a natural variational principle for equilibria. This variational principle, which is the basis for our computational method, is constructed as follows.

The total potential energy is given by the volume integral

∫_{D} [\frac{1}{2}|B|^{2} + \frac{1}{γ − 1}p] dV ,

where the form of the internal energy supposes an equation of state p = e^{s}ρ^{γ}, in which s is entropy per unit mass, ρ is mass density, and γ is the adiabatic index (1 < γ < ∞). The mass within any magnetic surface {Ψ = σ} is defined by

∫_{\{Ψ>σ\}} p^{\frac{1}{γ}} dV ;

strictly speaking, this integral also involves the entropy within the magnetic surface and so is identified with the mass only in a loose sense ([17]). The toroidal flux within the magnetic surface can also be expressed as a volume integral

∫_{{Ψ>σ}} B · ∇ϕ dV ;

and the vacuum toroidal flux is similarly expressible. An equilibrium configuration can be characterized as a minimizer of energy over the class of trial states (B, p) satisfying ∇ · B = 0 in D, appropriate fixed-and free-boundary conditions, and having prescribed values of the mass and flux constraints displayed above. A formal verification that such a constrained minimizer satisfies the equilibrium equations (2.1) is given in [17].

The precise formulation of this variational principle in the present context of axisymmetric toroidal equilibria makes use of the functionals

(2.10)  E(ψ, f, g) = \frac{1}{2} ∫_{Ω} [r^{−1}|∇ψ|^{2} + r^{−1}f^{2} + rg^{2}] drdz
(2.11)  F_{0}(f) = ∫_{Ω} r^{−1}f drdz
(2.12)  F_{σ}(ψ, f) = ∫_{Ω} r^{−1}f(ψ + \overline{ψ} − σ)_{+}drdz (σ ≥ σ_{0})
(2.13)  G_{σ}(ψ, g) = ∫_{Ω} rg^{2/γ}(ψ + \overline{ψ} − σ)_{+}drdz (σ ≥ σ_{0}) ,

where the notation s_{+} = max(s, 0) is used. Here and in the sequel the pressure p is replaced by the unknown

g: = (\frac{2p}{γ − 1})^{\frac{1}{2}} , (2.14)

as this substitution is convenient for technical reasons in what follows. The functional E clearly represents total potential energy. It is worth noting however that E involves the flux function ψ induced by the plasma current density and not the external poloidal field determined by \overline{ψ}. When rewritten in terms of the total flux function Ψ, this energy expression takes the form

E = \frac{1}{2} ∫_{Ω} [r^{−1}|∇Ψ|^{2} − 2Ψ\overline{J}_{ϕ} + r^{−1}f^{2} + rg^{2}] drdz

+ \frac{1}{2} ∫_{Ω} r^{−1}|∇\overline{ψ}|^{2} drdz

in which the contribution due to the external current density \overline{J}_{ϕ} enters explicitly. Nevertheless, given that \overline{ψ} is fixed according to (2.9), the form (2.10) furnishes the most concise expression for E needed in the further development. The functional F_{0} obviously is the total toroidal flux. The functionals F_{σ} and G_{σ}, on the other hand, require some interpretation. These two families of functionals are parametrized by the flux variable σ which runs through the range of Ψ = ψ + \overline{ψ} in the plasma. Differentiating with respect to this variable yields the formula

−\frac{d}{dσ} F_{σ}(ψ, f) = ∫_{\{ψ+\overline{ψ}>σ\}} r^{−1}f drdz ,

the right-hand side of which is recognized as the toroidal flux within the magnetic surface \{ψ + \overline{ψ} = σ\}. Similarly,

−\frac{d}{dσ} G_{σ}(ψ, g) = ∫_{\{ψ+\overline{ψ}>σ\}} rg^{2/γ} drdz

is identified with the mass within \{ψ +\overline{ψ} = σ\}. Thus, the σ-parametrized constraints F_{σ} and G_{σ} are simply the σ-antiderivatives of the classical integral invariants discussed above. (In the exceptional case mentioned above in which a limiting subdomain Ω_{0} is needed, F_{σ} and G_{σ} are defined by integrals extended over Ω_{0} rather than Ω.) Consequently, the functionals given in (2.11)-(2.13) constitute the complete family of conserved quantities associated with ideal MHD evolution, the conservation of poloidal flux being implied by the σ-parametrization of the magnetic surfaces.

The basic variational problem governing the equilibrium configurations we study can now be stated as

(P_{∞}) \{ \begin{matrix}E(ψ, f, g) → min subject to \\ F_{0}(f) = F_{0}^{∗}, F_{σ}(ψ, f) = F_{σ}^{∗}, G_{σ}(ψ, g) = G_{σ}^{∗} (σ ≥ σ_{0}) , \end{matrix}

where F_{0}^{∗}, F_{σ}^{∗}, G_{σ}^{∗} denote given constraint values. This minimization problem is novel in the sense that it imposes two (continuously) infinite families of constraints (as σ runs through the magnetic surfaces within the plasma). Consequently, the calculation of the variational equation (a Lagrange multiplier rule) satisfied by a minimizer (ψ, f, g) is not straightforward, since certain regularity properties are required to justify such a calculation. In order to arrive at a variational problem that is tractable to analysis and computation we therefore prefer to formulate a natural discretization of the constraint families in (P_{∞}). The resulting multiconstrained variational problem that we call (P_{n}) is constructed next.

Let σ_{0} < σ_{1} < σ_{2} < ⋯ < σ_{n} < +∞ be a partition of the interval σ_{0} ≤ σ < +∞, and let Δσ_{i}: = σ_{i} −σ_{i−1}. We then define what we refer to as “basis functions” relative to this partition by

Φ_{i}(s) : = \frac{1}{Δσ_{i}} ∫_{σ_{i−1}}^{σ_{i}} (s − σ)_{+} dσ − \frac{1}{Δσ_{i+1}} ∫_{σ_{i}}^{σ_{i+1}} (s − σ)_{+} dσ , (2.15)

or explicitly,

Φ_{i}(s) = \{ \begin{matrix}0 & s ≤ σ_{i−1} \\ (s − σ_{i−1})^{2}/2Δσ_{i} & σ_{i−1} ≤ s ≤ σ_{i} \\ s − (s − σ_{i})^{2}/2Δσ_{i+1} − (σ_{i−1} + σ_{i})/2 & σ_{i} ≤ s ≤ σ_{i+1} \\ (σ_{i+1} − σ_{i−1})/2 & s ≥ σ_{i+1} \end{matrix}

In fact, these functions are precisely the s-antiderivatives of the usual finite element functions Φ′_{i}(s); namely, Φ′_{i}(1 ≤ i ≤ n−1) is piecewise-linear on σ_{i−1} ≤ s ≤ σ_{i+1}, zero elseswhere, and Φ′_{i}(σ_{i}) = 1, while Φ′_{n} is linear on σ_{n−1} ≤ s ≤ σ_{n}, constant on s ≥ σ_{n}, and Φ′_{n}(σ_{n}) = 1. With these basis functions in hand we introduce the constraint functionals

(2.16)  F_{i}(ψ, f) = ∫_{Ω} r^{−1}fΦ_{i}(ψ + \overline{ψ}) drdz
(2.17)  G_{i}(ψ, g) = ∫_{Ω} rg^{2/γ}Φ_{i}(ψ + \overline{ψ}) drdz .

Now we can pose the multiconstrained variational problem

(P_{n}) \{ \begin{matrix}E(ψ, f, g) → min subject to \\ F_{0}(f) = F_{0}^{∗}, F_{i}(ψ, f) = F_{i}^{∗}, G_{i}(ψ, g) = G_{i}^{∗} (i = 1, …, n) , \end{matrix}

where these constraint values are derived from those in (P_{∞}) according to

F_{i}^{∗} = \frac{1}{Δσ_{i}} ∫_{σ_{i−1}}^{σ_{i}} F_{σ}^{∗}dσ − \frac{1}{Δσ_{i+1}} ∫_{σ_{i}}^{σ_{i+1}} F_{σ}^{∗}dσ

and similarly for G_{i}^{∗}.

The (first-order) variational equations associated with (P_{n}) are readily derived. Let F′ = (D_{ψ}F, D_{f}F, D_{g}F) denote the functional (or Fréchet) derivative of any (differentiable) functional F in the triple (ψ, f, g), employing the D-notational for the partial functional derivatives. The standard Lagrange multiplier rule ([15]) states that a solution (ψ, f, g) of (P_{n}) satisfies

E′(ψ, f, g) = λ_{0}F′_{0}(f) + ∑_{i=1}^{n} λ_{i}F′_{i}(ψ, f) + μ_{i}G′_{i}(ψ, g) (2.18)

for some multipliers λ_{0}, λ_{i}, μ_{i} ∈ . (Technically, these multipliers exist and are uniquely determined by the minimizer provided the constraint derivatives are linearly independent; this condition can be checked in the present context, but the necessary analysis is rather tedious and so will be omitted.) An explicit calculation of the expressions entering into (2.18) yields the triple of equations

(2.19)  Lψ = fλ · Φ′(ψ + \overline{ψ}) + r^{2}g^{2/γ} μ · Φ′(ψ + \overline{ψ})
(2.20)  f = λ_{0} + λ · Φ(ψ + \overline{ψ})
(2.21)  g = (2/γ)g^{(2−γ)/γ}μ · Φ(ψ + \overline{ψ}) .

Here the abbreviated notation λ · Φ = ∑_{i=1}^{n} λ_{i}Φ_{i} is used. These equations are merely a version of the Grad-Shafranov equation (2.6). Indeed, it is immediate that (2.20) defines the poloidal current profile f = f(Ψ), and it follows after some manipulation that (2.21) defines the pressure profile p = p(Ψ) = (γ − 1)g(Ψ)^{2}/2; thus, the (otherwise arbitrary) profile functions in (2.6) are determined by the multipliers λ_{0}, λ_{i}, μ_{i} according to

f(Ψ) = λ_{0} + λ · Φ(Ψ), p(Ψ) = \frac{(γ − 1)}{2} [\frac{2}{γ} μ · Φ(Ψ)]^{γ/(γ−1)} . (2.22)

In view of these relations, the equations (2.19)–(2.21) are seen to be equivalent to the classical equilibrium equations (2.6) and (2.7) in the plasma and vacuum regions. Moreover, the free-boundary conditions on S = {Ψ = σ_{0}} are also valid, because the continuity (and continuous differentiability) of ∇ψ, f and g are implied by these equations, by virtue of standard regularity theory ([9]). In summary, any minimizer (ψ, f, g) for (P_{n}) produces an exact solution of the complete equilibrium problem for the plasma-vacuum system.

The variational problem (P_{n}) can be viewed as a particular formulation of the general variational principle for ideal magnetohydrodynamic equilibria due to Woltjer [23, 24]. In his principle the constraints, which are derived as invariants of the associated evolution equations, take exactly the same form as our functionals F_{i} and G_{i}, except that the functions Φ_{i}(Ψ) are left as arbitrary. A derivation of the associated variational equations (analogous to (2.19)-(2.21)) is included in his work [23, 24]. Our presentation differs from his only in two ways: the choice of specific basic functions, which are tailored to the numerical analysis; and the complete treatment of the plasma-vacuum free boundary. It is important to notice that the development we give in this section connects the Kruskal-Kulsrud principle with the Woltjer principle by means of the idea that the constraints for (P_{∞}) (which are parametrized by the magnetic surface variable σ) are naturally discretized into the constraints for (P_{n}) (which are determined by the finite basis Φ_{i}). In a sense this discretization of constraints is achieved by averaging the σ-parametrized constraints over subintervals of the partition {σ_{i}}. Therefore, as the class of admissible triples (ψ, f, g) is enlarged in going from (P_{∞}) to (P_{n}), the solutions of (P_{n}) are themselves exact solutions of the governing equilibrium equations. Their associated profile functions f(Ψ) and p(Ψ), however, are special in that they are constructed from linear combinations of the basis functions Φ_{i}(Ψ), i = 1, …, n. As n tends to infinity, these (interpolating) profile functions tend to the general profile functions appropriate to (P_{∞}), at least formally. In this light we may view the formulation of (P_{n}) as the first step in the computation of the equilibrium configurations governed by (P_{∞}).

Some comments are in order concerning the relation between our formulation and the approach pioneered by Grad et. al. [10, 11, 12, 13], which relies on averaging over magnetic surfaces. For this purpose we recall the interpretation of our constraint functionals F_{σ} and G_{σ} as σ-antiderivatives of the (toroidal) flux and mass within {Ψ = σ}, and we observe therefore that

\frac{d^{2}}{dσ^{2}} F_{σ}(ψ, f) = ∫_{\{Ψ=σ\}} \frac{f}{r|∇Ψ|} dℓ = 2πq(σ)

\frac{d^{2}}{dσ^{2}} G_{σ}(ψ, g) = ∫_{\{Ψ=σ\}} \frac{rg^{2/γ}}{|∇Ψ|} dℓ = 2πm(σ) ,

where these identities define q and m, respectively, the specific toroidal flux (or safety factor) and specific mass (in the loose sense mentioned above) per unit poloidal flux. These integrals with respect to the length element dℓ on {Ψ = σ} ⊆ Ω (actually the one-dimensional Hausdorff measure on the level set) are obtained by applying the socalled coarea formula ([7]). Clearly, q(σ) and m(σ) (σ ≥ σ_{0}) constitute an equivalent form of the constraint data F_{σ}^{∗} and G_{σ}^{∗} in (P_{∞}); indeed, it is easy to check that the latter can be determined from the former by means of the formulas

F_{σ}^{∗} = 2π ∫_{σ_{0}}^{+∞}(s − σ)_{+}q(s)ds, G_{σ}^{∗} = 2π ∫_{σ_{0}}^{+∞}(s − σ)_{+}m(s)ds .

In the approach based on “generalized differential equations” (GDE) initiated by Grad, the prescribed data in the equilibrium problem consist precisely of the two functions q(σ) and m(σ). The GDE is then derived from the Grad-Shafranov equation (with unspecified poloidal current and pressure profile functions) using an averaging procedure, which replaces a quantity a = a(r, z) by its magnetic surface average

⟨a⟩: = ∫_{\{Ψ=σ\}} \frac{r a}{|∇Ψ|} dℓ;

such an averaged quantity can be considered as a function either of the flux σ or of the enclosed volume V = vol {Ψ > σ}. The GDE resulting from this procedure has the novel feature that it involves both the second-order elliptic operator L in (r, z) and second-order differentiation with respect to V (say). It therefore poses a nonstandard problem from the standpoint of either analysis ([19]) or computation ([12, 13, 11, 4]). On the other hand, our approach is based on the observation that the second-order antiderivatives of the prescribed data q(σ) and m(σ) furnish the σparametrized constraint data for the natural variational principle (P_{∞}) governing the same equilibrium problem, and further that a simple relaxation of those constraints leads to the variational problem (P_{n}) whose standard form makes it amenable to analytical and computational methods. Hence, our formulation obviates the necessity of employing the GDE formalism.

3 Iterative algorithm

We now construct an iterative procedure that solves the variational problem (P_{n}) formulated in the preceding section. In this section we develop the abstract form of the algorithm and derive its general convergence properties. In the succeeding section we state the algorithm in its concrete form and discuss its implementation for toroidal equilibrium computations.

A preliminimary transformation of (P_{n}) into an equivalent variational problem is required before the algorithm given below can be applied. For this reason we consider the modified variational problem

(P̃_{n}) \{ \begin{matrix}E(ψ, f, g) → min subject to \\ F̃_{0}(ψ, f) = F̃_{0}^{∗}, F_{i}(ψ, f) = F_{i}^{∗}, G_{i}(ψ, g) = G_{i}^{∗} (i = 1, …, n) , \end{matrix}

where the total toroidal flux constraint is replaced by a constraint imposed on the functional

F̃_{0}(ψ, f) = F_{0}(f) − ∑_{i=1}^{n} ω_{i}F_{i}(ψ, f) (3.1)

defined by some constants ω_{i} ≥ 0. This minor change in (P_{n}) is necessary because the structure of the algorithm rests on the property that the (Lagrange) multipliers corresponding to a minimizer are strictly positive. However, under typical circumstances, the multipliers for (P_{n}) are expected to have the signs λ_{0} > 0, λ_{i} < 0, μ_{i} > 0. Indeed, apart from the sign of λ_{0} which is set by the direction of the toroidal magnetic field (in the vacuum), the signs of these multipliers can be inferred from the monotonicity of the profile functions in (2.22)—namely, from the inequalities f′(Ψ) < 0, p′(Ψ) > 0, which are anticipated on physical grounds. This is so by virtue of the choice of the basis functions Φ_{i} as antiderivatives of the finite element functions, combined with the formulas

(3.2) f′(Ψ) = λ · Φ′(Ψ), p′(Ψ) = [\frac{2}{γ} μ · Φ(Ψ)]^{1/(γ−1)} μ · Φ′(Ψ) .

Now, if λ̃_{0}, λ̃_{i}, μ̃_{i} denote the multipliers for (P̃_{n}), then they are related to those for (P_{n}) by λ_{0} = λ̃_{0}, λ_{i} = λ̃_{i} − λ_{0} ω_{i} , μ_{i} = μ̃_{i} (i = 1, …, n) , (3.3)

as is immediate from the comparison of the variational equations (2.18)–(2.21) for (P_{n}) with their analogues for (P̃_{n}). Consequently, the constants ω_{1}, …, ω_{n} can always be fixed large enough (depending upon the solution) to ensure that all of the multipliers for (P̃_{n}) are strictly positive. Throughout the remainder of the section we assume that this is done.

For the sake of simplicity of exposition we let u = (ψ, f, g) denote the unknown triple, and we write the variational problem (P̃_{n}) in the abstract form

(3.4) E(u) → min subject to F_{i}(u) = γ_{i} (i = 0, …, 2n) .

In an obvious change of notation, the constraint family for (P̃_{n}) is now rewritten with the functionals F_{0}, …, F_{2n} replacing F̃_{0}, F_{1}, …, F_{n}, G_{1}, … G_{n}, and with corresponding constraint values γ_{0}, …, γ_{2n}. The admissible triple u is assumed to belong to the space H = H_{0}^{1}(Ω) × L^{2}(Ω) × L^{2}(Ω), for which a norm ∥ · ∥_{H} is given by

∥u∥_{H}^{2} = ∫_{Ω}[|∇ψ|^{2} + f^{2} + g^{2}] drdz .

Let ⟨ ·, ·⟩ denote the standard L^{2}-pairing in the sense of distribution theory— namely,

⟨u_{1}, u_{2}⟩: = ∫_{Ω}[ψ_{1}ψ_{2} + f_{1}f_{2} + g_{1}g_{2}] drdz ,

where f_{1}, f_{2}, g_{1}, g_{2} ∈ L^{2}(Ω), and ψ_{1} ∈ H_{0}^{1}(Ω), ψ_{2} ∈ H^{−1}(Ω) or vice versa. In terms of this pairing the (objective) functional E, which is differentiable at any u ∈ H, has the derivative

E′(u) = (D_{ψ}E, D_{f}E, D_{g}E) = (r^{−1}Lψ, r^{−1}f, rg) ∈ H^{−1}(Ω) × L^{2}(Ω) × L^{2}(Ω) ,

meaning that E(u + δu) = E(u) + ⟨E′(u), δu⟩ + o(∥δu∥_{H}) as δu tends to zero in H. Similar remarks pertain to the (constraint) functionals F_{i}. The variational equations satisfied by a minimizer u ∈ H, interpreted in this sense, are simply

E′(u) = ∑_{j=0}^{2n} λ_{j}F′_{j}(u) (3.5)

with (positive) multipliers λ_{0}, …, λ_{2n} (which are identified with λ_{0}, λ̃_{1}, …, λ̃_{n}, μ_{1}, …, μ_{n} in (P̃_{n})). The linear independence of the constraint gradients is needed here; this can be verified, but the necessary analysis will be omitted.

Convexity properties of the objective and constraint functionals are fundamental to the construction of the iterative algorithm. The objective functional E is both strictly convex and quadratic. Thus it admits the useful expansion

E(v) = E(u) + ⟨E′(u), v − u⟩ + E(v − u) , (3.6)

in which the positive second-order term is explicit. On the other hand, the constraint functionals F_{i} are not convex. Therefore, the algorithm involves a certain “convexification” of them, and requires that constants α_{i} ≥ 0 be fixed such that

(3.7) F_{i} + α_{i}E is convex (i = 0, …, 2n) .

Clearly, α_{i} can be chosen large enough to guarantee this property. With ω_{1}, …, ω_{n} and α_{0}, …, α_{2n} fixed appropriately the algorithm can be stated.

Let u^{0} ∈ H satisfying F_{i}(u^{0}) ≥ γ_{i} (i = 0, …, 2n) be an (otherwise arbitrary) initialization. The iterative sequence u^{k}, k = 0, 1, …, is defined inductively by solving the following quadratic programming subproblem at the iterative step: u = u^{k+1} solves

\{ \begin{matrix}E(u) → min subject to \\ F_{i}(u^{k}) + ⟨F′_{i}(u^{k}) + α_{i}E′(u^{k}), u − u^{k}⟩ ≥ γ_{i} & (i = 0, …, 2n) . \end{matrix} (3.8)

This subproblem is a convex optimization problem having linear (actually affine) inequality constraints. Therefore, it has a unique solution which is characterized by the corresponding Kuhn-Tucker conditions (the analogue of the Lagrange multiplier rule) [15]; namely, u^{k+1} ∈H and its associated multiplier (Kuhn-Tucker) vector ξ^{k+1} ∈ 2n+1 are uniquely determined by

E′(u^{k+1}) = ∑_{j=0}^{2n+1} ξ_{j}^{k+1} [F′_{j}(u^{k}) + α_{j} E′(u^{k})] (3.9)

ξ_{i}^{k+1} ≥ 0 (3.10) F_{i}(u^{k}) + ⟨F′(u^{k}) + α_{i}E′(u^{k}), u^{k+1} − u^{k}⟩ − γ_{i} ≥ 0 (3.11) ξ_{i}^{k+1} [F_{i}(u^{k}) + ⟨F′_{i}(u^{k}) + α_{i} E′(u^{k}), u^{k+1} − u^{k}⟩ − γ_{i}] = 0 . (3.12)

An explicit description of the algorithm defined by (3.8) can be given in terms of the so-called dual subproblem ([15, 21]). Let a matrix (a_{ij}^{k}) and a vector (c_{i}^{k}) be defined by the iterate u^{k} according to

a_{ij}^{k} = ⟨F′_{i}(u^{k}) + α_{i} E′(u^{k}), M^{−1}[F′_{j}(u^{k}) + α_{j} E′(u^{k})]⟩ (3.13) c_{i}^{k} = γ_{i} − F_{i}(u^{k}) + ⟨F′_{i}(u^{k}) + α_{i} E′(u^{k}), u^{k}⟩ , (3.14)

where M^{−1} denotes the (bounded) inverse of the (unbounded) linear operator M = E′′; specifically, M is that operator which corresponds to the positive definite quadratic form

E(u) = \frac{1}{2}⟨Mu, u⟩ (u ∈ H) . (3.15)

The subproblem in ξ that is dual to (3.8) is then expressible as

\frac{1}{2} ∑_{i,j=0}^{2n+1} a_{ij}^{k}ξ_{i}ξ_{j} − ∑_{i=1}^{2n+1} c_{i}^{k}ξ_{i} → min over ξ_{i} ≥ 0 . (3.16)

Now the algorithm can be described in a form that leads directly to a concrete numerical implementation. The iterative step is defined by

u^{k+1} = ∑_{j=0}^{2n+1} ξ_{j}^{k+1}M^{−1}[F′_{j}(u^{k}) + α_{j} E′(u^{k})] , (3.17)

where ξ^{k+1} is the (unique) solution of (3.16). The equivalence of (3.16), (3.17) with (3.8) is standard. The variational inequalities satisfied by ξ^{k+1} are simply

∑_{j=1}^{2n+1} a_{ij}^{k}ξ_{j}^{k+1} − c_{i}^{k} \{ \begin{matrix}≥ 0 & if ξ_{i}^{k+1} = 0 \\ = 0 & if ξ_{i}^{k+1} > 0 . \end{matrix} (3.18)

Noticing that (3.17) coincides with (3.9) by virtue of the identity E′(u) = Mu, and substituting this expression for u^{k+1} into the complementarity conditions (3.10)– (3.12), it is evident that (3.18) is just a restatement of (3.10)–(3.12).

We now turn to a discussion of the convergence properties of the algorithm. These special properties dictate the structure of the algorithm, and hence justify the particular construction described above. It is essential to explain them in conjunction with the definition of the iterative algorithm.

In short, the algorithm defined in (3.8) possesses the monotonicity properties:

E(u^{k+1}) ≤ E(u^{k}) , (3.19) F_{i}(u^{k}) ≥ γ_{i} (i = 0, …, 2n + 1) (3.20)

for every k and any admissible initialization u^{0}. An inductive argument on k proves this claim. Suppose that F_{i}(u^{k}) ≥ γ_{i} (i = 0, …, 2n + 1). Then u^{k} is admissible in the subproblem (3.8) which defines u^{k+1}, and hence (3.19) follows. Since F_{i} + α_{i}E is convex by hypothesis, there holds

F_{i}(u^{k+1}) + α_{i}E(u^{k+1}) ≥ F_{i}(u^{k} + α_{i} E(u^{k})

+ ⟨F′_{i}(u^{k}) + α_{i} E′(u^{k}), u^{k+1} − u^{k}⟩ (3.21)

≥ γ_{i} + α_{i} E(u^{k}) ,

invoking the linear inequality constraints in (3.8). Now (3.21) and (3.19) combine to give F_{i}(u^{k+1}) ≥ γ_{i} (i = 0, …, 2n + 1), completing the induction.

A stronger version of the monotonicity property (3.19) is true; namely,

E(u^{k} − u^{k+1}) ≤ E(u^{k}) − E(u^{k+1}) , (3.22)

for every k. This inequality can be derived from the Kuhn-Tucker conditions (3.9)– (3.12) along with the expansion (3.6). The required calculation is

(3.23)  E(u^{k}) − E(u^{k+1}) = ⟨E′(u^{k+1}), u^{k} − u^{k+1}⟩ + E(u^{k} − u^{k+1})
= ∑ ξ_{j}^{k+1}⟨F′_{j}(u^{k}) + α_{j} E′(u^{k}), u^{k} − u^{k+1}⟩
+ E(u^{k} − u^{k+1})
≥ ∑ ξ_{j}^{k+1}[F_{j}(u^{k}) − γ_{j}] + E(u^{k} − u^{k+1})
≥ E(u^{k} − u^{k+1}) .

From (3.22) follows both the monotonic convergence of the objective functional values

E(u^{k}) ↓ E^{∗} (say) as k → ∞ ,

and the convergence property of the iterative sequence:

∥u^{k} − u^{k+1}∥_{H} → 0 as k → ∞ . (3.24)

Without a rate for this convergence it cannot be asserted that the (entire) sequence {u^{k}} converges to a (single) limit point u^{∗} ∈ H. Nevertheless, the bound E(u^{k}) ≤ E(u^{0}) ensures at least that every subsequence has a further subsequence converging weakly in H. Now, if u^{∗} denotes such a weak limit point, then a straightforward argument (which is omitted here) furnishes the corresponding weak convergence of the Kuhn-Tucker conditions; namely, the associated subsequence of multiplier vectors ξ_{j}^{k} tends to a limit ξ^{∗}, and the pair u^{∗}, ξ^{∗} satisfies

⟨E′(u^{∗}), v⟩ = ∑_{j=0}^{2n+1} ξ_{j}^{∗}⟨F′_{j}(u^{∗}) + α_{j} E′(u^{∗}), v⟩ for all v ∈ H , (3.25)

ξ_{i}^{∗} ≥ 0 (3.26) F_{i}(u^{∗}) − γ_{i} ≥ 0 (3.27) ξ_{i}^{∗}[F_{i}(u^{∗}) − γ_{i}] = 0 . (3.28)

But (3.25) is just the weak form of the equation

E′(u^{∗}) = ∑_{j=0}^{2n+1} ξ_{j}^{∗}[F′_{j}(u^{∗}) + α_{j} E′(u^{∗})] ,

which in turn is equivalent to the governing variational equations (3.5) when the multipliers λ_{i}^{∗} associated with u^{∗} are defined by

λ_{i}^{∗} = ξ_{i}^{∗}/(1 − α · ξ^{∗}) . (3.29)

Thus, u^{∗} is a critical point (presumably a minimum point) for the variational problem (3.4) whenever (a) ξ_{i}^{∗} > 0 for all i, and (b) α · ξ^{∗} < 1. As remarked above, condition (a), which in view of condition (b) is the same as the condition that λ_{i}^{∗} > 0 for all i, is guaranteed by the choice of ω (and some physical considerations). Hence the complementarity conditions (3.28) enforce the equality constraints in (3.4). Therefore, it suffices to verify condition (b). This can be inferred from the identity

(1 − α · ξ^{∗})⟨E′(u^{∗}), u^{∗}⟩ = ∑ ξ_{j}^{∗}⟨F′_{j}(u^{∗}), u^{∗}⟩ ,

which results from taking v = u^{∗} in (3.25). Noting that the form of the objective and constraint functionals generally implies that

⟨E′(u), u⟩ > 0 , ⟨F′_{j}(u), u⟩ > 0 ,

for all admissible functions u, it is then clear that the latter identity can be satisfied only when condition (b) holds. (The stated inequality for E is always valid, while the inequalities for F_{j} are valid under some very mild conditions on the solution u^{∗} = (ψ^{∗}, f^{∗}, g^{∗}) of (P̃_{n}), which hold in all cases of real interest.)

The convergence properties demonstrated above refer to subsequences of the iterative sequence {u^{k}}. An alternative statement of these makes use of the distance in the space H between a point u∈H and a set S^{∗} ⊆ H:

dist_{H}(u, S^{∗}) : = inf{∥u − v∥_{H} : v ∈ S^{∗}} .

If S^{∗} is taken to be the set of solutions (actually critical points) of the variational problem (3.4), then it can be shown that

dist_{H}(u^{k}, S^{∗}) → 0 as k → ∞ . (3.30)

This generalized form of the convergence statement is required because of the possibility that solutions of (3.4) are not unique. The proof of (3.30) depends upon standard elliptic regularity theory (applied to the operator L) which can be used to demonstrate that the weakly convergent sequences of u^{k} converge strongly in H. As this proof is of technical interest only, it will not be included in the present discussion. (A further analysis based on the same techniques yields the analogue of (3.30) in the strongest norm compatible with the optimal regularity of the solution u^{∗}. For instance, an analysis of (P̃_{n}) shows that the solution u^{∗} = (ψ^{∗}, f^{∗}, g^{∗}) belongs to C^{2}(\overline{Ω}) × C^{1}(\overline{Ω}) × C^{1}(\overline{Ω}), hence furnishing at least the regularity of classical solutions; this analysis depends upon the particular form of the piecewise-quadratic basis functions Φ_{i}.) The reader is referred to our earlier work [6] which contains detailed convergence proofs in the case of a prototype problem derived from a simple version of (P_{n}).

In implemented computations such as those documented in §5 the iterative sequence is observed to converge to a single limit point with a linear rate of convergence. Thus, the theoretical issues raised above can be ignored in practice. Heuristically, this convergence behavior is expected whenever the algorithm is applied to compute a (local) energy minimizing solution subject to the given flux and mass constraints. Indeed, the convergence rate can be related to the stability of the equilibrium to axisymmetric perturbations; however, this relationship will not be pursued here.

4 Numerical implementation

The concrete form of the iterative algorithm defined in §3 can now be presented. For this purpose we first record in component form some of the expressions used in the abstract statement of the algorithm. For the sake of clarity in the presentation we choose γ = 2 here and in the implemented computations discussed in §5. The objective functional E can be represented as the quadratic form (3.15) associated with the (symmetric) linear operator

Mu = E′(u) = (r^{−1}Lψ, r^{−1}f, rg) . (4.1)

The constraint functionals F_{0}, …, F_{2n}, which are identified with F̃_{0}, F_{1}, …, F_{n}, G_{1}, …, G_{n}, have the derivatives

F′_{0}(u) = (r^{−1}f Φ̃′_{0}(ψ + \overline{ψ}), r^{−1}Φ̃_{0}(ψ + \overline{ψ}), 0) (4.2)

F′_{i}(u) = (r^{−1}f Φ′_{i}(ψ + \overline{ψ}), r^{−1}Φ_{i}(ψ + \overline{ψ}), 0)
F′_{n+i}(u) = (rg Φ′_{i}(ψ + \overline{ψ}), 0, rΦ_{i}(ψ + \overline{ψ}))

for i = 1, …, n. The basis function Φ̃_{0}(s) : = 1 − ω · Φ(s) is included so that the functional F̃_{0} may be expressed in the form

F̃_{0}(u) = ∫_{Ω} r^{−1}f Φ̃_{0}(ψ + \overline{ψ}) drdz .

Some functionals derived from these expressions are introduced for convenience:

P_{ij}(u) : = ⟨F′_{i}(u), M^{−1}F′_{j}(u)⟩ (4.3) Q_{i}(u) : = ⟨F′_{i}(u), u⟩ , (4.4)

for i = 0, …, 2n. According to (4.1) the terms involving M^{−1} are calculated to be

M^{−1}F′_{0}(u) = (w_{0}, Φ̃_{0}(ψ + \overline{ψ}), 0) (4.5)

M^{−1}F′_{i}(u) = (w_{i}, Φ_{i}(ψ + \overline{ψ}), 0)
M^{−1}F′_{n+i}(u) = (w_{n+i}, 0, Φ_{i}(ψ + \overline{ψ})) ,

where w_{0}, …, w_{2n} are determined by solving the linear elliptic boundary value problems:

L w_{0} = f Φ̃′_{0}(ψ + \overline{ψ})  in Ω,  w_{0} = 0  on ∂Ω
L w_{i} = f Φ′_{i}(ψ + \overline{ψ})  in Ω,  w_{i} = 0  on ∂Ω
L w_{n+i} = r^{2}g Φ′_{i}(ψ + \overline{ψ})  in Ω,  w_{n+i} = 0  on ∂Ω .

(4.6)

The iterative step that generates u^{k+1} = (ψ^{k+1}, f^{k+1}, g^{k+1}) from u^{k} = (ψ^{k}, f^{k}, g^{k}) can be described as the three stage process:

(Stage 1) Find the solutions w_{o}^{k}, …, w_{2n}^{k} of (4.6) corresponding to (ψ^{k}, f^{k}, g^{k}).

(Stage 2) Evaluate the integrals E(u^{k}), F_{i}(u^{k}), P_{ij}(u^{k}), Q_{i}(u^{k}) (i = 0, …, 2n); then assemble the coefficients

a_{ij}^{k} = P_{ij}(u^{k}) + α_{i} Q_{j}(u^{k}) + α_{j} Q_{i}(u^{k}) + 2α_{i}α_{j} E(u^{k}) (4.7) c_{i}^{k} = γ_{i} − F_{i}(u^{k}) + Q_{i}(u^{k}) + 2α_{i} E(u^{k}) , (4.8) and solve the (dual) quadratic programming subproblem (3.16) for ξ^{k+1}.

(Stage 3) Set

ψ^{k+1} = ∑_{j=0}^{2n+1} ξ_{j}^{k+1} w_{j} + (α · ξ^{k+1}) ψ^{k} (4.9)

f^{k+1} = ξ_{0}^{k+1} Φ̃_{0}(ψ^{k} + \overline{ψ}) + ∑_{j=1}^{n} ξ_{j}^{k+1} Φ_{j}(ψ^{k} + \overline{ψ}) + (α · ξ^{k+1}) f^{k} (4.10)

g^{k+1} = ∑_{j=1}^{n} ξ_{n+j}^{k+1} Φ_{j}(ψ^{k} + \overline{ψ}) + (α · ξ^{k+1}) g^{k} . (4.11)

Stage 1 involves the solution of 2n + 1 Dirichlet problems in Ω. These standard linear problems can be treated by any available numerical method, since the structure of the algorithm is independent of the particular discretization used.

Stage 2 constructs and then solves the dual subproblem for ξ^{k+1}. The several numerical integrations required to compute the coefficients given in (4.7) and (4.8) can be accomplished by any numerical quadrature appropriate to the discretization used. The determination of ξ^{k+1}, a convex quadratic minimization problem with simple inequality constraints, can be achieved by a variety of known optimization methods. In the event that none of the constraints holds as an equality (ξ_{i}^{k+1} > 0 for all i), the multiplier vector ξ^{k+1} is just the solution of the linear system

∑_{j=0}^{2n+1} a_{ij}^{k} ξ_{j}^{k+1} = c_{i}^{k} (i = 0, …, 2n) .

This event is typical in implemented computations, since the basic variational problem (P̃_{n}) is formulated in such a way that its solution u^{∗} corresponds to positive multipliers ξ_{i}^{∗}. Accordingly, the determination of ξ_{i}^{k+1} is substantially simplified. (An effective strategy is as follows: if the solution ξ^{k+1} of the linear system satisfies ξ_{i}^{k+1} > 0 for all i, then it is accepted as the multiplier vector; otherwise, the multiplier vector is defined to be the solution ξ^{k+1} of the quadratic programming problem.)

Stage 3 defines u^{k+1} in terms of u^{k} and ξ^{k+1} according to the explicit formula (3.17) for the solution of the primal subproblem (3.8).

The algorithm exhibited above is completely specified once the “shifting” parameters ω_{1}, …, ω_{n} and the “convexifying” parameters α_{0}, …, α_{2n} are fixed. As explained in §3 the parameters ω_{i} depend upon the solution through f(Ψ), and consequently they must be chosen adaptively to ensure the positivity of the multipliers λ̃_{i}^{∗} associated with the solution (ψ^{∗}, f^{∗}, g^{∗}) of (P̃_{n}). The choice of the parameters α_{i} may also be made adaptively on the basis of a posteriori information obtained in practice, but certain a priori estimates for these parameters are available in theory. In particular, lower bounds for the parameters α_{i} which constitute sufficient conditions for the convexity properties (3.7) can be inferred from the positivity of the second-order derivatives of the functionals F_{i} + α_{i}E. For example, in order to give an a priori lower bound for α_{i} with 1 ≤ i ≤ n, the following calculations suffice:

⟨F′′_{i}(u)δu, δu⟩ = ∫_{Ω} r^{−1}[f Φ′′_{i}(ψ + \overline{ψ}) (δψ)^{2} + 2Φ′_{i}(ψ + \overline{ψ}) δψδf] drdz
≥ − ∫_{Ω} r^{−1}[{1 + (Δσ_{i+1})^{−1}|f|} (δψ)^{2} + (δf)^{2}] drdz
⟨E′′(u)δu, δu⟩ = ∫_{Ω}[r^{−1}|∇δψ|^{2} + r^{−1}(δf)^{2} + r(δg)^{2}] drdz
≥ ∫_{Ω}[Λ_{Ω}r^{−1}(δψ)^{2} + r^{−1}(δf)^{2} + r(δg)^{2}] drdz;

here the Poincaré inequality ([9]) is utilized with a constant Λ_{Ω} (the smallest eigenvalue for L in Ω). Then (3.7) holds whenever

Λ_{Ω} α_{i} ≥ 1 + (Δσ_{i+1})^{−1} max_{Ω} |f| (1 ≤ i ≤ n) .

Similarly, lower bounds for α_{i} with either i = 0 or n + 1 ≤ i ≤ 2n can be found. Such estimates are useful as guides to the appropriate choice of the parameters α_{i}. Nevertheless, they are too crude to be employed effectively in the implemented algorithm, in which the resulting monotonicity and convergence properties are observed to persist for values of the parameters α_{i} that are considerably smaller than the estimated values. Moreover, the performance of the algorithm improves as the parameters α_{i} are decreased (within allowable limits), both with respect to the conditioning of the dual subproblem (3.16) and the rate of convergence of the iterative sequence. In this regard, our computations suggest that the optimal choice of α_{i} is roughly 1/10 of its estimated value, and that nearby choices (within a factor of 3) result in nearly the same performance; on the other hand, our tests indicate that the algorithm performs poorly when the parameters α_{i} are chosen much too small or much too large.

We now address the question of how to select the constraint values F_{0}^{∗}, F_{i}^{∗}, G_{i}^{∗} as appropriate given data for the variational problem (P_{n}). From the standpoint of modelling and computation the most natural answer is first to derive this data from some known equilibrium configuration, and then to use the iterative algorithm to determine a family of solutions to (P_{n}) corresponding to constraint data that are varied (incrementally) in some prescribed manner. For instance, a family of fluxconserving equilibria are generated if the flux constraints F_{0}^{∗}, F_{i}^{∗} are fixed while the mass constraints G_{i}^{∗} are varied. Alternatively, a family of adiabatically compressed equilibria are obtained if all of the constraints F_{0}^{∗}, F_{i}^{∗}, G_{i}^{∗} are fixed and the external poloidal field (defined by \overline{ψ}) is varied by changing the currents in the external coils. The computation of such families of equilibrium configurations is indeed the main goal of the present paper. In addition, an obvious initialization of the algorithm is available when solutions are computed incrementally.

In view of these remarks, we need only supply a method of solving the standard equilibrium problem with prescribed toroidal current profile and total toroidal current, in order to start the variational problem (P_{n}) in the sense just explained. Below we describe a simple and efficient method of this kind, which is also based on a variational formulation. Here we only sketch the method since it is fully discussed in another paper [5], where the analogous hydrodynamic problems are treated.

The equilibrium problem for the plasma-vacuum system is to be solved with a given total toroidal current I_{0} > 0, and given profile functions in the Grad-Shafranov equation (2.6) having the form

f(Ψ) = f_{0}((Ψ − σ_{0})_{+}), p(Ψ) = p_{0}((Ψ − σ_{0})_{+}) ;

the prescribed functions f_{0}(s) and p_{0}(s) are assumed to be smooth for s ≥ 0 and to satisfy f_{0}(0) > 0, f′_{0}(0) = 0 and p_{0}(0) = 0, p′_{0}(0) = 0. The problem is then to find ψ and σ_{0} such that

Lψ = rj(r, ψ + \overline{ψ} − σ_{0}) in Ω, ψ = 0 on ∂Ω , (4.12)

∫_{Ω} j(r, ψ + \overline{ψ} − σ_{0}) drdz = I_{0} , (4.13) where the toroidal current density in the plasma is written as

j(r, s) = r^{−1}f_{0}(s_{+})f′_{0}(s_{+}) + rp′_{0}(s_{+}) .

In contrast to (P_{n}) the flux constant σ_{0} is an unknown in this problem and determines the location of the plasma-vacuum interface \{ψ + \overline{ψ} = σ_{0}\} in response to the total current constraint.

The iterative algorithm appropriate to this problem can be stated as follows: given ψ^{k}, define ψ^{k+1} and σ_{0}^{k+1} by

Lψ^{k+1} = rj(r, ψ^{k} + \overline{ψ} − σ_{0}^{k+1}) in Ω, ψ^{k+1} = 0 on ∂Ω (4.14)

∫_{Ω} j(r, ψ^{k} + \overline{ψ} − σ_{0}^{k+1}) drdz = I_{0} . (4.15)

This iterative step is easily implemented by first finding σ_{0}^{k+1} to satisfy (4.15) and then solving the linear elliptic boundary value problem (4.14) for ψ^{k+1}. The reader is referred to the abovementioned paper [5] for a discussion of the variational structure of free-boundary problem (4.12), (4.13) and the concomitant global convergence theory for the algorithm (4.14), (4.15).

5 Computed examples

In this section we present the results of some computations made with the general algorithm described above. For the sake of brevity, we limit our discussion to two representative cases, each of which is chosen to exhibit an aspect of the performance of the algorithm. In the first case, we generate a family of flux-conserving equilibria by increasing the mass constraints while fixing the flux constraints. In the second case, we generate a family of adiabatically compressed equilibria by varying the external poloidal field (determined by toroidal current coils) while fixing all of the constraints. These two cases suffice to demonstrate the effectiveness of our method under conditions that model those encountered in real confinement devices.

As is explained in §4, such sequences of equilibria are initiated by computing an equilibrium with prescribed toroidal current profile j(r, s) and total toroidal current I_{0} from which (initial) constraint values F_{0}^{∗}, F_{i}^{∗}, G_{i}^{∗} are derived. In the cases discussed below the prescribed profiles (which define the Grad-Shafranov equation) are taken to have the form

\frac{1}{2}f_{0}(s)^{2} = a_{0} + a_{1}s_{+}^{κ+1}/(κ + 1)
p_{0}(s) = a_{2}s_{+}^{κ+1}/(κ + 1)
j(r, s) = (a_{1}/r + a_{2}r) s_{+}^{κ} ,

where a_{0} > 0, a_{1} < 0, a_{2} > 0 and κ > 0. For appropriate values of the given parameters I_{0}, a_{0}, a_{1}, a_{2} and κ, solutions of the standard equilibrium problem (4.12), (4.13) are furnished by the simple iterative algorithm (4.14), (4.15). The above parameters can be varied in order to bring the physical properties of the solutions within acceptable bounds, as is common practice (see [3], for instance).

For the sake of simplicity, the cross-sectional domain is taken to be a rectangle Ω = {(r, z): r_{0} < r < r_{1}, |z| < z_{1}} in each case, although the general algorithm is not restricted to this geometry. The operator L in Ω is discretized by a standard finite difference method, and the linear elliptic boundary-value problems constituting Stage 1 of the algorithm are treated with a corresponding fast Poisson solver. All of the integrals required in Stage 2 of the algorithms are computed using appropriate quadrature formulas consistent with the discretization. The (dual) quadratic programming subproblem for ξ^{k+1} is solved by a standard routine based on active set method (see [?]). This method, which employs the direct linear solution of reduced systems, is precise and efficient; moreover, it exploits the special feature of the algorithm that, by construction, all of the constraints on ξ^{k+1} are inactive for sufficiently large k. (In fact, only a few constraints on ξ^{k+1} are active at any iteration k under most circumstances.) The explicit formulas for (ψ^{k+1}, f^{k+1}, g^{k+1}) which form Stage 3 of the algorithm are evaluated at each grid node.

In each case the iterations are terminated when the following stopping criteria are achieved:

max \{\frac{∥ ψ^{k+1} − ψ^{k} ∥_{2}}{∥ ψ^{k} ∥_{2}} , \frac{∥ f^{k+1} − f^{k} ∥_{2}}{∥ f^{k} ∥_{2}} , \frac{∥ g^{k+1} − g^{k} ∥_{2}}{∥ g^{k} ∥_{2}} \} < 0.003

max \{\frac{F̃_{0}^{k+1} − F̃_{0}^{∗}}{F_{0}^{∗}} , \frac{F_{i}^{k+1} − F_{i}^{∗}}{F_{i}^{∗}} , \frac{G_{i}^{k+1} − G_{i}^{∗}}{G_{i}^{∗}} \} < 0.003 ,

where F_{0}^{k+1}, F_{i}^{k+1}, G_{i}^{k+1} denote the constraints evaluated at (ψ^{k+1}, f^{k+1}, g^{k+1}). Among the examples given below (which are representative) between 10 and 50 iterations are required to satisfy the above criteria. At the beginning of the iterations the performance of the algorithm is chiefly influenced by the choice of initialization. The global convergence properties demanded in the construction of the algorithm are then crucial, since they ensure computational robustness even when a rough guess is used to initialize the iterative sequence. At the end of the iterations the algorithm always exhibits a linear rate of convergence, which depends strongly on the choice of the parameters ω_{1}, …, ω_{n} and (especially) α_{0}, …, α_{2n}. The best convergence rates are achieved by varying these parameters adaptively along the iterative sequence, taking them no larger than necessary to obtain the required monotonicity and convergence properties of the algorithm. A full description of the adaptive method employed in our implemented code is given in our technical report [?].

The first family of computed equilibria illustrates the transition from low-β to high-β, the volume-averaged β being defined by

β: = ∫ pdV/ ∫ \frac{1}{2}B^{2}dV .

The computational domain Ω = {0.5 < r < 1.5, |z| < .5} is discretized with a grid having 101 × 51 nodes. There are no external field coils. The multiconstrained variational problem (P_{n}) is considered with n = 3 corresponding to equal flux increments Δσ_{i} = (max ψ − σ_{0})/3, where (ψ, σ_{0}) is the solution of (4.12), (4.13) used to initiate the family in the sense described above; the parameters determining this solution are taken to be

a_{0} = 50, a_{1} = −5, a_{2} = 50, κ = 1, I_{0} = 3.5 .

The constraint values furnished by this solution are

F_{0} = 0.0 and σ_{0} = 0.0 ,

F_{1} = 0.0 , F_{2} = 0.0 , F_{3} = 0.0 (5.1)

G_{1} = 0.0 , G_{2} = 0.0 , G_{3} = 0.0 .

A family of four solutions of (P_{3}) is then computed by fixing F_{0}^{∗}, F_{1}^{∗}, F_{2}^{∗}, F_{3}^{∗} and σ_{0} at these values while increasing G_{1}^{∗}, G_{2}^{∗}, G_{3}^{∗} incrementally starting from these values. We refer to the resulting solutions as equilibria 1.1 through 1.4. Table 1 lists the prescribed values of G_{i}^{∗} used, and the corresponding computed values of β, I (total toroidal plasma current) and max ψ (total poloidal flux between the magnetic axis and the conducting shell). Figure 2 displays the magnetic surface plots associated with these equilibria; several level curves {ψ = σ} are shown in the plasma and the vacuum regions, and the free-boundary {ψ = σ_{0}} is indicated by the bolder curve. Table 2 contains the multipliers λ_{0}, λ_{i}, μ_{i} associated with the solutions, and Figure 3 displays the corresponding cross-sectional (z = 0) profiles of f, p and J_{ϕ}.

Along this first family β increases from % to %, and the outward shift of the magnetic surfaces (including the plasma-vacuum interface) grows correspondingly; also, I increases from 0.0 to 0.0, and the toroidal current density profile becomes peaked on the outward side of the plasma. Thus, the family exhibits the expected behavior of flux-conserving tokomak equilibria, as anticipated by other treatments based on the GDE method [?, 4]. For the purpose of comparison, however, it is important to note that our results pertain to the full free-boundary problem formulated in §1 and, as such, therefore appear to be new.

The tabulation of max ψ is included in Table 1 to demonstrate the degree to which (poloidal) flux is conserved by this family. Namely, while flux-conservation is imposed only in the (volume-averaged) sense that F_{i}(i = 0, …, 3) are constant along the family, the variation of max ψ, which is not fixed in (P_{3}), is remarkably small. Similar results are observed over a wide range of conditions whenever (P_{n}) is solved with n ≥ 3. On the other hand, max ψ can vary greatly when n = 1. This evidence strongly supports the claim that exact (infinitely-constrained) flux-conservation can be effectively approximated by a few integral constraints (say, 3 ≤ n ≤ 5). Since the computational effort of each iteration grows rapidly with n, this claimed property of the variational formulation (P_{n}) is very significant whenever the algorithm is applied in practice.

Table 1. Some computed quantities for equilibria 1.1–1.4 corresponding to F_{i}^{∗} as in (5.1) and G_{i}^{∗} as shown.

Equilibrium G_{1}^{∗} G_{2}^{∗} G_{3}^{∗} β I max ψ
1.1  * * * * *  *
1.2  * * * * *  *
1.3  * * * * *  *
1.4  * * * * *  *

Table 2. Multipliers associated with equilibria 1.1–1.4

Equilibrium λ_{0} λ_{1} λ_{2} λ_{3} μ_{1} μ_{2} μ_{3}
1.1  * * * * * * *
1.2  * * * * * * *
1.3  * * * * * * *
1.4  * * * * * * *

The second family computed equilibria illustrates the change in shape of the magnetic surfaces during compression due to varying the external poloidal field. The computational domain Ω = {2.5 < r < 3.5, |z| < 1.5} is used with a grid having 51 × 51 nodes. The external field is induced by two elementary coils located at (R_{1}, Z_{1}) = (2.7, 0) and (R_{2}, Z_{2}) = (3.3, 0), carrying currents I_{1} = I_{2} > 0. The flux function \overline{ψ} for this field is the sum of two Green functions (for L) corresponding to the coil locations (R_{ℓ}, Z_{ℓ}); their singular parts can be evaluated precisely by an explicit formula (in terms of elliptic integrals), and their regular parts can then be computed numerically. As in the construction of the first family, the multiconstrained variational problem is posed with n = 3 and equal flux increments Δσ_{i} derived from a solution of (4.12), (4.13) with given parameters

a_{0} = 5, a_{1} = −0.5, a_{2} = 5, κ = 1, I_{0} = 1.5 .

The constraint values fixed by this solution are

F_{0} = 0.0 and σ_{0} = 0.0

F_{1} = 0.0 , F_{2} = 0.0 , F_{3} = 0.0 , (5.2)

G_{1} = 0.0 , G_{2} = 0.0 , G_{3} = 0.0 .

A family of four solutions of (P_{n}) is then generated by successively increasing the external coil currents I_{1} = I_{2}, while maintaining all of these constraint values. We refer to the resulting solutions as equilibria 2.1–2.4. Table 3 records the prescribed values of I_{1} and I_{2}, and the corresponding values of β, I and max Ψ. Figure 3 shows the magnetic surface plots for these equilibria.

The most striking feature evidenced by this family is the large z-shift of the plasma from equilibrium 2.3 to 2.4, and the associated loss of z-symmetry. This behavior can be interpreted as the result of a two-dimensional (axisymmetric) instability due to increasing compression, whereupon an energetically more favorable equilibrium is found by the (energy-decreasing) algorithm. By symmetry, the z-reflection of equilibrium 2.4 has identical properties, and so under these circumstances the algorithm may converge to either equilibrium 2.4 or its reflection, depending only on the (very small) numerical errors that initiate a departure from z-symmetry. Consequently, many iterations may be needed to develop such a shift, if no other pertubations are added. It is noteworthy that again max Ψ shows only a small variation in the course of this large change in the magnetic surfaces {Ψ = σ}.

Table 3. Some computed quantities for equilibria 2.1–2.4 corresponding to F_{0}, F_{i}, G_{i} as in (5.2) and I_{1}, I_{2} as shown.

Equilibrium I_{1} I_{2} β I max Ψ
2.1  * * * *  *
2.2  * * * *  *
2.3  * * * *  *
2.4  * * * *  *

It is possible to devise other geometries and external coil configurations which will exhibit bifurcation of the plasma region rather than the above z-shift. The variational formulation (P_{n}) remains unchanged through such a bifurcation, even though the interior of magnetic surfaces {ψ > σ} may not be connected. (The integrals defining the constraints then extend over all of the components of these sets.) Thus, the above algorithm can be applied (without any modifications) to situations where magnetic island formation is allowed, without the necessity of providing a priori information about the topology of the magnetic surfaces. In these situations our formulation and algorithm appear to have a distinct advantage over the GDE method, which encounters difficulties if the magnetic surfaces are not (regularly) nested.

References

[1] G. Bateman, MHD Instabilities, MIT press, Cambridge, Massachusetts, 1978.

[2] F. Bauer, O. Betancourt and P. Garabedian, A Computational Method in Plasma Physics, Springer-Verlag, Berlin/Heidelberg/New York, 1978.

[3] J. Blum, Numerical Simulation and Optimal Control in Plasma Physics, with Applications to Tokomaks, Wiley/Gauthiers-Villars, Chichester/New York, 1989.

[4] R.A. Dory and Y.-K. M. Peng, High pressure flux-preserving tokomak equilibria, Nucl. Fusion 17 (1977), 21-31.

[5] A. Eydeland and B. Turkington, A computational method of solving freeboundary problems in vortex dynamics, J. Comp. Phys. 78 (1988), 194-214.

[6] A. Eydeland, J. Spruck and B. Turkington, Multiconstrained variational problems of nonlinear eigenvalue type: new formulations and algorithms, Math. of Comp., in press.

[7] H. Federer, Geometric Measure Theory, Springer-Verlag, Berlin/Heidelberg/New York, 1969.

[8] J. P. Freidberg, Ideal Magnetohydrodynamics, Plenum Press, New York, 1987.

[9] D. Gilbarg and N. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer-Verlag, Berlin/Heidelberg/New York, 1977.

[10] H. Grad, Magnetic confinement fusion energy research, Proceeding of Symposia in Applied Math. 21 (1977), 3-40.

[11] H. Grad, Survey of 1\frac{1}{2} D transport codes, Courant Institute of Mathematical Sciences Report MF-93, New York University, 1978 (unpublished).

[12] H. Grad, P.N. Hu and D.C. Stevens, Adiabatic evolution of plasma equilibria, Proc. Nat. Acad. Sci. 72 (1975), 3789-3793.

[13] H. Grad, P.N. Hu, D.C. Stevens and E. Turkel, TITLE?, Plasma Phys. Contr. Nucl. Fusion Res. 2 (1976), 355-365.

[14] H. Grad and H. Rubin, Hydromagnetic equilibria and force-free fields, Proc. of 2nd U. N. International Conf. on the Peaceful Uses of Atomic Energy, U.N., Geneva (1958), vol. 31, 190-197.

[15] A.D. Ioffe and V.M. Tihomirov, Theory of Extremal Problems, Elsevier North-Holland, New York, 1979.

[16] J.L. Johnson et. al., Numerical determination of axisymmetric toroidal magnetohydrodynamic equilibria, J. Comp. Phys. 32 (1979), 212-234.

[17] M. D. Kruskal and R. M. Kulsrud, Equilibrium of a magnetically confined plasma in a toroid, Physics of Fluids 1(1958), 265–274.

[18] A. E. Lifschitz, Magnetohydrodynamics and Spectral Theory, Kluwer Academic Publishers, Dordrecht, 1989.

[19] P. Laurence and E. Stredulinsky, A new approach to queer differential equations, Commun. Pure Appl. Math. 38 (1985), 333-355.

[20] R. Lust and A. Schluter, Axial symmetrische magnetohydrodynamische Gleichgewichtskonfigurationen, Z. Naturforsch. 12a (1957), 850-854.

[21] R.T. Rockafellar, Convex Analysis, Princeton Univ. Press, Princeton, 1970.

[22] V.D. Shafranov, On equilibrium magnetohydrodynamical configurations, Sov. Phys. - JETP 6 (1958), 545-554.

[23] L. Woltjer, Hydromagnetic equilibrium,III: axisymmetric incompressible media, Astrophysical J. 130 (1959), 400-404.

[24] ———-, Hydromagnetic equilibrium,IV: axisymmetric compressible media, Astrophysical J. 130 (1959), 405-413.