On multiple solutions of the Grad–Shafranov equation

The Grad–Shafranov equation (GSE) for axisymmetric magnetohydromagnetic (MHD) equilibria is a nonlinear, scalar partial differential equation (PDE) which in principle can have no solution, one solution or multiple solutions. The topic of solutions of the GSE has not been widely explored in the literature so far even though calculating the correct equilibrium is vital to any stability or transport analysis that is carried out. In particular, an improved understanding of the existence of multiple equilibria may give insights into plasma stability and triggers for major disruptions, the avoidance of which is key on the route to a fusion power plant. It may also have implications for predictive scenario modelling or equilibrium reconstruction.

The topic has been discussed in Solano [1] at a conceptual level but no equilibria with multiple solutions were shown in that paper. Indeed, that paper called for developments to allow the calculation and understanding of multiple solutions of the GSE. This paper answers that call. Solano [1] also raises the important question of what happens when a solution branch disappears as transport in the plasma changes the profiles. We must assume that the initial plasma equilibrium is lost on an Alfvénic timescale and the plasma will evolve to a new state. It may transition to a different nontrivial equilibrium (albeit at a much lower plasma pressure i.e. a major disruption), a periodic orbit, or a chaotic trajectory. The work by Berestycki and Brezis [2] discusses existence and uniqueness but the results depend on restrictions to the profile shapes. A mathematically focussed treatment of the problem has also been given in Jeyakumar et al [3] but again without the calculation of examples. Results based on a force-free plasma may also provide insights to this problem, see Taylor [4] for example.

An analytic demonstration of multiple solutions of the GSE was given in Schnack [5], in the case of a free boundary equilibrium surrounded by a perfectly conducting wall. This problem was shown to have either two solutions, one deeply confined and the other shallow, or no solution. However, the equilibrium is assumed to be a thin, vertically extended plasma. This does not produce an equilibrium that it is easy to match with a numerical model. Our analysis overcomes this limitation, and provides solutions that match well with numerical simulations.

Our analytic model is still an idealisation and does not provide solutions to the complete GSE. We must understand the solutions in realistic geometry to be sure of the relevance to tokamak experiments. We therefore deploy a recently developed algorithm called deflated continuation [6, 7] to find multiple solutions in a full geometry, discretising the GSE with the finite element method using Firedrake [8].

This paper is organised as follows. We describe the GSE and give the weak form in section 2. We develop a new analytic model—assuming a large aspect ratio torus—which produces an example of no (non-trivial) solution, moving to two solutions as a control parameter is varied. In section 3, we describe deflated continuation, a method for finding multiple solutions of PDEs. We apply this technique to the test case that we developed and solved analytically in section 2 and find good agreement. In section 4, we use deflated continuation to understand the effect of aspect ratio and plasma shaping on multiple solutions of the GSE. In section 5 we discuss possible implications for this work and finally we give conclusions.

The GSE is very well known in magnetic confinement fusion. It is used ubiquitously, but it is important to understand when multiple solutions (or non-existence of non-trivial solutions) may occur and the possible consequences.

In a non-rotating plasma the GSE is derived from force balance, namely $\textbf{J} \times \textbf{B} = \nabla p$, where J is the plasma current density, B is the magnetic field and p is the isotropic plasma pressure (see, for example Goedbloed et al [9] or Freidberg [10] for more details on the derivation). Using cylindrical coordinates, see figure 1, the GSE becomes

$$\begin{align} & \nabla^2 \psi - \frac{2}{R}\frac{\partial \psi}{\partial R} \equiv R\frac{\partial}{\partial R}\left( \frac{1}{R}\frac{\partial \psi}{\partial R} \right)\nonumber\\ & \quad + \frac{\partial^2 \psi}{\partial Z} = -F\left(\psi\right)F^{^{\prime}}\left(\psi\right)-R^2p^{^{\prime}}\left(\psi\right) \end{align} \tag{ 1 }$$

where $\psi(R,Z)$ is the poloidal magnetic flux to be solved for, R is the coordinate in the major radius direction, Z is the vertical coordinate, and $F = RB_{\phi}$ is the poloidal current stream function, where B_φ is the component of the magnetic field in the toroidal direction.

Zoom In Zoom Out Reset image size

The weak form of the GSE is to find $\psi \in H^1(\Omega)$ such that

$$\begin{align} \int \frac{1}{R} \nabla \psi \cdot \nabla \xi d \Omega = \int \left[\mu_0 R p^{^{\prime}} +\frac{FF^{^{\prime}}}{R} \right] \xi d\Omega \quad \text{for all } \xi \in H^1\left(\Omega\right), \end{align} \tag{ 2 }$$

where Ω is a bounded Lipschitz domain, ξ is the test function, $p^{^{\prime}} = \frac{dp(\psi)}{d\psi}$ and $F^{^{\prime}} = \frac{dF(\psi)}{d\psi}$.

The functions $p(\psi)$ and $F(\psi)$ are specified in this problem. These profiles could in principle be the output of a transport model. The pressure, p, is generally always a monotonically decreasing function. This means that pʹ will be negative across the whole plasma. There may be regions of large gradient and possibly very small gradient. The function FFʹ can be both positive and negative and so the right-hand side may change sign in the domain.

2.1. An analytic model with multiple solutions

We assume a large aspect ratio toroidal plasma. Specifically, we assume that the radius of the magnetic axis, R₀, is much larger than the minor radius, a, so that $R_0 \gg a$. In this case we see that the $\frac{2}{R} \frac{\partial \psi}{\partial R}$ term can be ordered small compared to the other terms in (1). This results in the Grad–Shafranov operator reducing to a Laplacian on the left hand side. If we now move to a cylindrical coordinate system with the axis along the axis of the plasma, $(R_0, Z_0)$, and we assume poloidal symmetry then we can reduce the problem to one dimension. We also assume that at sufficiently large aspect ratio the flux surfaces are concentric, nested tori which coincide with the surfaces of constant radius, and that the right-hand side is a pure function of ψ. This is reasonable since at large aspect ratio the variation in R over the minor radius will be small. The theory of cylindrical equilibrium is well known. This is equivalent to using the $\beta_p \approx 1$, $\epsilon \ll 1$ expansion of the GSE given in [10]. This simple model will allow us to explore multiple solutions of the GSE and provide test cases for comparison to numerical simulations.

In cylindrical coordinates, assuming no poloidal variation, (1) becomes

$$\begin{equation} \frac{\partial^2 \psi}{\partial r^2} + \frac{1}{r}\frac{\partial \psi}{\partial r} = G\left(\psi\right) \end{equation} \tag{ 3 }$$

where r is the minor radius variable and $G(\psi)$ is a function of ψ which mimics the behaviour of the right-hand side. This function is often monotonically decreasing in physical cases but this does not necessarily hold.

Clearly if G is a constant then the differential equation is linear, and a unique solution will exist. If $G(\psi) = \lambda \psi$, then a unique solution exists whenever λ is not an eigenvalue. However, if G is nonlinear then the differential equation becomes nonlinear and thus permits more interesting behaviour of solutions. As an example, we take a tanh function

$$\begin{equation} G\left(\psi\right) = \frac{C}{2}\left[ 1- \tanh\left(\frac{\psi-\psi_p}{w}\right) \right] \end{equation} \tag{ 4 }$$

where we take C and w as parameters that can be varied. This function allows us to use a continuous function to investigate the situation where we have a linear pressure profile which transitions to a flat pressure profile at a given value. In the limit that w → 0 we have a discontinuity which is the function that was used by Schnack [5]. Figure 2 shows this function for $C = -5$, $\psi_p = -1$ and w = 0.05. Other nonlinear functions may also produce multiple solutions.

Zoom In Zoom Out Reset image size

We use a shooting method to find the solutions of equation (3), fixing $\psi_p = -1$ and w = 0.05, enforcing the boundary condition that $\psi(2) = 0$. As we vary C, we find that for C < 2.7 there is only one solution, the trivial solution. However, for C > 2.7 we find that there are two additional solutions. The bifurcation diagram is shown in figure 3.

Zoom In Zoom Out Reset image size

We can also produce a bifurcation diagram by varying w and keeping C = 2.9 fixed. Figure 4 shows that there are three solutions when w < 0.88 and only one when w > 0.88.

Zoom In Zoom Out Reset image size

We can also cast the problem of finding multiple solutions in a different way. We can ask that we have a fixed value of central poloidal flux and investigate the equilibria at different values of another control parameter. We have varied the value of w with C = 2.9 and fixed values of ψ₀ to confirm that we can find two solutions, one with high w and one with low w. This is a situation more akin to equilibrium reconstruction, especially with magnetics only, where measurements of the internal plasma profiles are not used.

Our analytic model is useful for insight, but the GSE cannot generally be reduced to a one dimensional equation for solution in this way. It is fundamentally a two dimensional problem. In this section, we will apply an algorithm of bifurcation analysis, deflated continuation [6, 7], to find multiple solutions of the GSE in full geometry. The algorithm is applied to a finite element discretisation implemented in Firedrake [8, 11]. Firedrake lets users write the variational formulation of their problem in a domain-specific language. It also allows users to flexibly construct finite element spaces for use in Galerkin's method applied to this variational formulation ³ . We will compare the results from deflated continuation with our analytic model in appropriate limiting cases.

The aim of the deflated continuation algorithm is to compute the solutions of an equation

$$\begin{equation} f\left(u, \lambda\right) = 0, \end{equation} \tag{ 5 }$$

where u is the solution of a PDE and $\lambda \in \mathbb{R}$ is a parameter on which the PDE depends. For each parameter value, each known solution branch is continued to that parameter, and then the nonlinear problem is modified so that those solutions are excluded (the solutions are deflated, just as known roots of polynomials can be deflated by dividing by the appropriate factor). The algorithm then seeks solutions of the modified nonlinear problem; if any are found, they are new solutions of the original problem. Importantly, this approach can compute disconnected bifurcation diagrams, like those already presented in figures 3 and 4.

We now apply deflated continuation to our model equation (3). We take the weak form of (3) and approximate its solutions with a piecewise linear finite element method on a disk (radius, a = 2 and origin $(0,0$) in two dimensions). We then use deflated continuation (as implemented in the Defcon library ⁴ ) to explore the solutions of the problem. We have solved the cylindrical variant of the problem, i.e. without toroidal effects, in these cases so R₀ has no effect.

Figures 5–7 show three solutions that have been calculated for C = 2.9 with this procedure. The central value of the flux agrees between the two methods to $0.9\%$ or better. This excellent agreement gives confidence to move forward to applying the same computational methods to the full GSE. It should also be noted that the finite element solutions and the shooting method are both well converged. These are real multiple solutions of the reduced GSE arrived at by completely different numerical methods and not numerical artifacts.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An issue to note is that deflated continuation is not guaranteed to find all the solutions of the equation. We observed failure of the algorithm if there is a discontinuity or too sharp a gradient in the function G. For the case in [5] with a discontinuity, our code only finds one non-trivial solution unless a very specific initial guess is given. In the example above if w < 0.1 then again our code does not find the roots. However, for $w \unicode{x2A7E} 0.1$ there is excellent agreement between deflated continuation and the results from the shooting method. We bear this caveat in mind as we proceed to a more realistic case. It should be noted that the shooting method is very robust.

In section 3 we showed excellent agreement between our analytic model and finite element simulations for a plasma with circular cross-section. We now apply the same computational techniques to solve the problem with a more complete geometry to understand the effects of aspect ratio and plasma shaping on the multiple equilibria. In this section we will assume that the profile of F is a constant so that $FF^{^{\prime}} = 0$ and that the pressure has the profile we used previously,

$$\begin{equation} p^{^{\prime}}\left(\psi\right) = \frac{C}{2}\left[ 1- \tanh\left(\frac{\psi-\psi_p}{w}\right) \right]. \end{equation} \tag{ 6 }$$

This represents a plasma that has a constant gradient of pressure with flux up to the value ψ_p and then transitions to zero after that. The weak formulation of the problem that we now solve then becomes: find $\psi \in H^1(\Omega)$ such that

$$\begin{align} & \int \frac{1}{R} \nabla \psi \cdot \nabla \xi d \Omega = \int \left[ \frac{C \mu_0 R }{2}\left[ 1- \tanh\left(\frac{\psi-\psi_p}{w}\right) \right] \right] \xi d\Omega \nonumber\\ & \quad \text{for all } \xi \in H^1\left(\Omega\right), \end{align} \tag{ 7 }$$

since $FF^{^{\prime}} = 0$ if F is taken as a constant. The boundary condition is that ψ = 0.

The difference between equation (7) and (the weak form of) (3) is the $1/R$ term that appears in the left hand side of (7) and also the factor of R that appears on the right hand side of (7). These extra geometric terms break the poloidal symmetry of the plasma which exists in the cylindrical case but not in the toroidal case.

We first study the effect of aspect ratio. We take a circular cross section plasma of minor radius, a = 1, and change the major radius R₀ with a fixed value of w = 0.18 and $R_0^2 C = 6$. The results are plotted in figure 8. We see that the aspect ratio makes little difference to the number of solutions in this case and also to their values. It has the strongest effect at small aspect ratio when the variation of R across the plasma cross section is relatively largest. The results of the two approaches converge at large aspect ratio.

Zoom In Zoom Out Reset image size

We next study the effect of plasma shaping at low aspect ratio. We parametrise the plasma shape with the standard elongation, κ, and triangularity, δ as defined by

$$\begin{align} R & = R_0 + a \cos\left(\theta +\delta \sin\left(\theta\right) \right), \end{align} \tag{ 8 }$$

$$\begin{align} Z & = a \kappa \sin\left(\theta\right), \end{align} \tag{ 9 }$$

where θ is the poloidal angle. We used Gmsh [12] to generate appropriate triangular meshes of these geometries. We show an example of the mesh and solution for a plasma with triangularity δ = 0.3 and elongation κ = 2 in figure 9.

Zoom In Zoom Out Reset image size

We fix $R_0 = 3$ for this computation and vary elongation, figure 10 and triangularity, figure 11. We see that plasma triangularity does not have a strong influence on the number of solutions. However, for small enough elongation the two solution branches coalesce and disappear. This is not in a regime where tokamaks are normally operated as high elongation has other benefits. The strong effect of increasing elongation on the lower solution is partly due to the plasma cross sectional area changing significantly with elongation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have presented methods to calculate multiple tokamak equilibria solutions of the GSE and demonstrate their use in relatively simple cases. It should be noted that if there is nonlinearity in the profiles used in the GSE then there are likely to be multiple solutions. The results are generic in the mathematical sense. We have shown that as system parameters are varied we can lose equilibrium solutions. One immediate question is what will happen to the plasma in this case. If the plasma is no longer in an equilibrium state then it will evolve on an Alfvénic timescale until it reaches a new equilibrium, or transitions to a periodic orbit, or to a chaotic orbit. In principle, this may result in a complete loss of the plasma, but this is not seen in experiments where the initial stage of a disruption is the loss of the thermal energy, the thermal quench. The plasma that is left is much cooler and broadly force-free. It may be that some types of major tokamak disruption are caused by such a loss of equilibrium. For example, if we follow the highest confinement equilibrium branch in figure 3 and we say the parameter C is continuously dropping (as may occur in experiment as a result of plasma transport) the equilibrium will be lost once C < 2.7 and this would be experienced as a major disruption experimentally. We will attempt to make more detailed comparisons to experimental results in future work.

A key part of understanding tokamak plasmas is equilibrium reconstruction where diagnostic measurements are used to constrain an equilibrium model to infer the experimentally realised equilibrium. There are various ways this is done, for various purposes. If only magnetic diagnostic data is used then there are no measurements of the internal profiles. In this case highly simplified profiles are used for FFʹ and pʹ which strongly constrains the potential shapes of the reconstructed equilibria. We have shown that two very different profile shapes can produce the same flux on axis. If internal measurements of the plasma are used, for example from MSE or Thomson Scattering diagnostics, in what is often called kinetic equilibrium reconstruction, then the profile shapes are given more freedom by using more basis functions. However, there are more measurements to constrain the final equilibrium. Extreme care still needs to be taken here to ensure that the basis functions used do not overly constrain the potential profiles or generate spurious detail. It is important to ensure that the choice of basis functions for equilibrium reconstruction does not determine the conclusions of studies of the stability and transport of a given plasma.

The effect of multiple equilibria should be evaluated in predictive modelling where an equilibrium model is coupled to a transport model. These two are evolved in time to simulate a tokamak pulse, therefore ramping up the current to flat top and then ramping down. It may be imagined that in this process bifurcation points will be passed where an equilibrium model only finds one of the available equilibrium branches. It may be that an improved shot evolution can be found (or a deleterious one avoided) if all bifurcation points are properly understood. The work done on theoretically predicting super-H mode [13] may give indications of how this should be done.

A further observation is that our analytic model is more accurate and robust than we might expect. This means that we can use this method to explore a wide range of potential profile shapes to find out when multiple solutions might be important. We may also use simple transport models along with this analytic model to understand when equilibria may be lost. This could improve our understanding of some types of major disruption or how we can plot a path to higher performance plasmas.

We have assumed that the plasma here is axisymmetric, indeed this is required as we are using the GSE. However, if we have a non-axisymmetric model of the plasma we could compute not only additional axisymmetric equilibria but also non-axisymmetric states that could be thought of as saturated instabilities such as the helical core mode [14]. This may also be the topic of further work.

In this work we have answered a call to develop techniques to calculate multiple solutions to the GSE. We have provided both analytical and numerical techniques and shown that they agree in the correct limits. We have shown that even in a relatively simple model that bifurcation points can occur. If we have more complicated profiles of FFʹ and pʹ then we would expect more of this behaviour. These tools can be exploited in future work to potentially improve our understanding of disruptions and to to find routes to improved performance in tokamak plasmas.

We are grateful to Wayne Arter, Andrew Kirk and Lucy Kogan for their comments on this work. C J H dedicates this work to Alan Ham and Kathleen Ham.

This work has been funded by the EPSRC Energy Programme [grant number EP/W006839/1]. To obtain further information on the data and models underlying this paper please contact PublicationsManager@ukaea.uk . ⁵

For the purpose of open access, the author(s) has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

In this appendix we demonstrate that the β and plasma current are relevant to practical situations for the large aspect ratio case. The theory of cylindrical equilibrium is well known. We reproduce some results here for reference based on the $\beta_p \approx 1$, $\epsilon \ll 1$ expansion of the GSE given in Freidberg [10]. We start with the equations of static ideal MHD

$$\begin{eqnarray} &&\textbf{J} \times \textbf{B} = \nabla p, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} &&\nabla \cdot \textbf{B} = 0, \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} &&\textbf{J} = \frac{1}{\mu_0}\nabla \times \textbf{B}. \end{eqnarray} \tag{ 12 }$$

We assume cylindrical symmetry so all variables are independent of θ and z. We can then calculate the current density

$$\begin{equation} \mu_0 \textbf{J} = \left[0, -\frac{dB_z\left(r\right)}{dr}, \frac{1}{r}\frac{d}{dr}\left(rB_{\theta}\left(r\right)\right)\right], \end{equation} \tag{ 13 }$$

and so the radial component of the static momentum equation becomes

$$\begin{equation} \frac{d}{dr}\left[p + \frac{\left(B_{\theta}^2+B_{z}^2\right)}{2\mu_0} \right] = -\frac{B_{\theta}^2}{\mu_0 r}, \end{equation} \tag{ 14 }$$

where we may specify two of the profiles, i.e. B_θ and p(r) with the third, $B_z(r)$ in this case, being determined by these two.

We may multiply equation (14) by $\pi r^2$ and integrate from 0 to a to get

$$\begin{equation} \langle p \rangle + \left \langle \frac{B_z^2}{2 \mu_0} \right \rangle - \frac{B_z^2\left(a\right)}{2 \mu_0}- \frac{B_{\theta}^2\left(a\right)}{2 \mu_0} = 0, \end{equation} \tag{ 15 }$$

where

$$\begin{equation} \langle X \rangle = \frac{1}{\pi a^2}\int_0^a 2 \pi r X\left(r\right) dr \end{equation} \tag{ 16 }$$

is the average value of X(r) over a cross-section of radius a.

The ratio of the thermal pressure to the magnetic pressure is the plasma β:

$$\begin{equation} \langle \beta \rangle = \frac{2 \mu_0 \langle p \rangle}{B^2\left(a\right)}. \end{equation} \tag{ 17 }$$

The poloidal and toroidal parameters are defined as

$$\begin{equation} \beta_p = \frac{2 \mu_0 \langle p \rangle}{B_{\theta}^2\left(a\right)}, \qquad \beta_t = \frac{2 \mu_0 \langle p \rangle}{B_z^2\left(a\right)}. \end{equation} \tag{ 18 }$$

We have assumed for simplicity in this paper that B_z is a constant so (15) results in $\beta_p = 1$.

A.1. Safety factor profile

The twist of the magnetic field lines is captured by the safety factor which in a cylindrical screw pinch is

$$\begin{equation} q\left(r\right) = \frac{r B_z\left(r\right)}{R_0 B_{\theta}\left(r\right)}. \end{equation} \tag{ 19 }$$

In this approximation we have [10]

$$\begin{eqnarray} B_{r0}\left(r\right) & = & 0, \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} B_{\theta 0}\left(r\right) & = & \frac{1}{R_0}\frac{d\psi_0}{dr}, \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} B_{z 0}\left(r\right) & = & \frac{-F\left(\psi_0\right)}{R_0}. \end{eqnarray} \tag{ 22 }$$

We have assumed that B_z is a constant so

$$\begin{equation} q\left(r\right) = \frac{r B_z}{\psi_0^{^{\prime}}}. \end{equation} \tag{ 23 }$$

A.2. Toroidal current profile

We need an expression for the toroidal current J_φ which can be obtained from the derivation of the GSE (and in the $\beta_p ~1$, $\epsilon \ll 1$ expansion of the GSE given in Freidberg [10])

$$\begin{equation} -\mu_0 R_0 J_z\left(r\right) = \nabla^2 \psi = -\mu_0 R_0^2 \frac{dp}{d\psi}- \frac{1}{2}\frac{dF^2}{d\psi}, \end{equation} \tag{ 24 }$$

and so in our case, given we have taken F as a constant,

$$\begin{equation} J_z\left(r\right) = R_0 \frac{dp}{d\psi} = \frac{R_0 C}{2}\left[ 1- \tanh\left(\frac{\psi-\psi_p}{w}\right) \right]. \end{equation} \tag{ 25 }$$

A.3. Equilibrium scaling

The equilibrium can be scaled following the rules given in Lutjens et al [15]. These allow us first to specify the value of the q profile at one radial location, therefore allowing us to move the q profile up and down but not to change its shape. We can also scale the value of ψ and so we can alter the toroidal β. We cannot alter the poloidal β or the internal inductance with these scalings. The first rescaling is

$$\begin{equation} \psi_\textrm{new} = \alpha_1 \psi_\textrm{old}, \qquad F_\textrm{new} = \alpha_1 F_\textrm{old}, \qquad p_{new} = \alpha^2_1 p_\textrm{old}, \end{equation} \tag{ 26 }$$

and the second rescaling is

$$\begin{equation} F^2_\textrm{new} = F^2_\textrm{old} + \alpha_2. \end{equation} \tag{ 27 }$$

We use (27) first to pick the q-profile and then (26) is used to set the β_t and toroidal current. We have used the scalings to produce a set of equilibria with pressure profiles shown in figure 12 and safety factor profiles shown in figure 13. These show that the pressures and safety factor profiles for these simple cases are not pathological in any sense. The equilibrium scaling shows that suitable equilibria can always be found and that when a bifurcation point exists the two equilibria can always be scaled to be sensible.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size