On Multiple Solutions of the Grad-Shafranov Equation PDF

LETTER OPEN ACCESS You may also like - Intuition for the radial penetration of flux On multiple solutions of the Grad–Sha...

LETTER OPEN ACCESS You may also like - Intuition for the radial penetration of flux On multiple solutions of the Grad–Shafranov surface shaping in tokamaks Justin Ball and Felix I Parra equation - Magnetics only real-time equilibrium reconstruction on ASDEX Upgrade L Giannone, M Weiland, R Fischer et al. To cite this article: C.J. Ham and P.E. Farrell 2024 Nucl. Fusion 64 034001 - Shafranov shift correction to the Furth–Yoshikawa scaling of tokamak adiabatic compression T Nicolas, H Lütjens, O Sauter et al. View the article online for updates and enhancements. This content was downloaded from IP address 178.211.159.10 on 14/11/2024 at 08:33  International Atomic Energy Agency Nuclear Fusion Nucl. Fusion 64 (2024) 034001 (11pp) https://doi.org/10.1088/1741-4326/ad1d77 Letter On multiple solutions of the Grad–Shafranov equation C.J. Ham1,∗ and P.E. Farrell2 1 UKAEA, Culham Science Centre, Abingdon OX14 3DB, United Kingdom of Great Britain and Northern Ireland 2 Mathematical Institute, University of Oxford, Oxford OX2 6GG, United Kingdom of Great Britain and Northern Ireland E-mail: [email protected] Received 16 October 2023, revised 11 December 2023 Accepted for publication 11 January 2024 Published 2 February 2024 Abstract The Grad–Shafranov equation (GSE) for axisymmetric MHD equilibria is a nonlinear, scalar PDE which in principle can have zero, one or more non-trivial solutions. The conditions for the existence of multiple solutions has been little explored in the literature so far. We develop a simple analytic model to calculate multiple solutions in the large aspect ratio limit. We compare the results to the recently developed deflated continuation method to find multiple solutions in a realistic geometry and right-hand side of the GSE using the finite element method. The analytic model is surprisingly accurate in calculating multiple solutions of the GSE for given boundary conditions and the two methods agree well in limiting cases. We examine the effect of plasma shaping and aspect ratio on the multiple solutions and show that shaping generally does not alter the number of solutions. We discuss implications for predictive modelling, equilibrium reconstruction, plasma stability and disruptions. Keywords: MHD, Grad–Shafranov equation, nonlinear PDEs, multiple solutions 1. Introduction of multiple equilibria may give insights into plasma stability and triggers for major disruptions, the avoidance of which is The Grad–Shafranov equation (GSE) for axisymmetric mag- key on the route to a fusion power plant. It may also have netohydromagnetic (MHD) equilibria is a nonlinear, scalar implications for predictive scenario modelling or equilibrium partial differential equation (PDE) which in principle can have reconstruction. no solution, one solution or multiple solutions. The topic of The topic has been discussed in Solano at a conceptual solutions of the GSE has not been widely explored in the liter- level but no equilibria with multiple solutions were shown in ature so far even though calculating the correct equilibrium that paper. Indeed, that paper called for developments to allow is vital to any stability or transport analysis that is carried the calculation and understanding of multiple solutions of the out. In particular, an improved understanding of the existence GSE. This paper answers that call. Solano also raises the important question of what happens when a solution branch ∗ disappears as transport in the plasma changes the profiles. Author to whom any correspondence should be addressed. We must assume that the initial plasma equilibrium is lost on an Alfvénic timescale and the plasma will evolve to a new Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any state. It may transition to a different nontrivial equilibrium further distribution of this work must maintain attribution to the author(s) and (albeit at a much lower plasma pressure i.e. a major disrup- the title of the work, journal citation and DOI. tion), a periodic orbit, or a chaotic trajectory. The work by © 2024 Crown copyright, United Kindom Atomic Energy Authority & 1741-4326/24/034001+11$33.00 Printed in the UK 1 The Author(s).Published by IOP Publishing Ltd on behalf of the IAEA Nucl. Fusion 64 (2024) 034001 C.J. Ham and P.E. Farrell Berestycki and Brezis 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 but again without the cal- culation of examples. Results based on a force-free plasma may also provide insights to this problem, see Taylor for example. An analytic demonstration of multiple solutions of the GSE was given in Schnack , in the case of a free boundary equi- librium surrounded by a perfectly conducting wall. This prob- lem was shown to have either two solutions, one deeply con- fined and the other shallow, or no solution. However, the equi- librium 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 limita- tion, 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 relev- ance 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 Figure 1. Cartoon of the coordinate system used for the with the finite element method using Firedrake. Grad–Shafranov equation. This paper is organised as follows. We describe the GSE and give the weak form in section 2. We develop a new ana- lytic model—assuming a large aspect ratio torus—which pro- coordinate, and F = RBϕ is the poloidal current stream func- duces an example of no (non-trivial) solution, moving to two tion, where Bϕ is the component of the magnetic field in the solutions as a control parameter is varied. In section 3, we toroidal direction. describe deflated continuation, a method for finding multiple The weak form of the GSE is to find ψ ∈ H1 (Ω) such that solutions of PDEs. We apply this technique to the test case ˆ ˆ [ ] that we developed and solved analytically in section 2 and find 1 FF ′ ∇ψ · ∇ξ dΩ = µ0 Rp ′ + ξ dΩ for all ξ ∈ H1 (Ω) , good agreement. In section 4, we use deflated continuation to R R understand the effect of aspect ratio and plasma shaping on (2) multiple solutions of the GSE. In section 5 we discuss possible implications for this work and finally we give conclusions. where Ω is a bounded Lipschitz domain, ξ is the test function, p ′ = dpd(ψ) ′ dF(ψ) ψ and F = dψ. The functions p(ψ) and F(ψ) are specified in this problem. 2. GSE These profiles could in principle be the output of a transport The GSE is very well known in magnetic confinement fusion. model. The pressure, p, is generally always a monotonically It is used ubiquitously, but it is important to understand when decreasing function. This means that p′ will be negative across multiple solutions (or non-existence of non-trivial solutions) the whole plasma. There may be regions of large gradient and may occur and the possible consequences. possibly very small gradient. The function FF ′ can be both In a non-rotating plasma the GSE is derived from force bal- positive and negative and so the right-hand side may change ance, namely J × B = ∇p, where J is the plasma current dens- sign in the domain. ity, B is the magnetic field and p is the isotropic plasma pres- sure (see, for example Goedbloed et al or Freidberg for 2.1. An analytic model with multiple solutions more details on the derivation). Using cylindrical coordinates, see figure 1, the GSE becomes We assume a large aspect ratio toroidal plasma. Specifically, ( ) we assume that the radius of the magnetic axis, R0 , is much 2 ∂ψ ∂ 1 ∂ψ larger than the minor radius, a, so that R0 ≫ a. In this case we ∇ ψ− 2 ≡R R ∂R ∂R R ∂R see that the R2 ∂ψ ∂ R term can be ordered small compared to the ∂2ψ other terms in (1). This results in the Grad–Shafranov oper- + = −F (ψ) F ′ (ψ) − R2 p ′ (ψ) (1) ator reducing to a Laplacian on the left hand side. If we now ∂Z move to a cylindrical coordinate system with the axis along where ψ(R, Z) is the poloidal magnetic flux to be solved for, R the axis of the plasma, (R0 , Z0 ), and we assume poloidal sym- is the coordinate in the major radius direction, Z is the vertical metry then we can reduce the problem to one dimension. We 2 Nucl. Fusion 64 (2024) 034001 C.J. Ham and P.E. Farrell also assume that at sufficiently large aspect ratio the flux sur- faces are concentric, nested tori which coincide with the sur- faces 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 βp ≈ 1, ϵ ≪ 1 expansion of the GSE given in. 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 ∂ 2 ψ 1 ∂ψ + = G (ψ) (3) ∂r2 r ∂r where r is the minor radius variable and G(ψ) is a function of ψ which mimics the behaviour of the right-hand side. This Figure 2. Plot of the tanh function with C = −5, ψp = −1 and function is often monotonically decreasing in physical cases w = 0.05. 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(ψ) = λψ, then a unique solution exists whenever λ is not an eigen- value. However, if G is nonlinear then the differential equation becomes nonlinear and thus permits more interesting beha- viour of solutions. As an example, we take a tanh function [ ( )] C ψ − ψp G (ψ) = 1 − tanh (4) 2 w 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 func- tion that was used by Schnack. Figure 2 shows this function for C = −5, ψp = −1 and w = 0.05. Other nonlinear functions may also produce multiple solutions. We use a shooting method to find the solutions of Figure 3. Bifurcation diagram for our toy model varying the equation (3), fixing ψp = −1 and w = 0.05, enforcing the parameter C with w = 0.05 and ψp = −1. The trivial solution boundary condition that ψ(2) = 0. As we vary C, we find always exists but at C ≈ 2.7 two new solutions appear at a fold that for C < 2.7 there is only one solution, the trivial solu- bifurcation and persist as C increases. tion. However, for C > 2.7 we find that there are two additional solutions. The bifurcation diagram is shown in figure 3. We can also produce a bifurcation diagram by varying w 3. Deflated continuation 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. Our analytic model is useful for insight, but the GSE cannot We can also cast the problem of finding multiple solutions generally be reduced to a one dimensional equation for solu- in a different way. We can ask that we have a fixed value of tion in this way. It is fundamentally a two dimensional prob- central poloidal flux and investigate the equilibria at differ- lem. In this section, we will apply an algorithm of bifurcation ent values of another control parameter. We have varied the analysis, deflated continuation [6, 7], to find multiple solutions value of w with C = 2.9 and fixed values of ψ 0 to confirm that of the GSE in full geometry. The algorithm is applied to a we can find two solutions, one with high w and one with low finite element discretisation implemented in Firedrake [8, 11]. w. This is a situation more akin to equilibrium reconstruction, Firedrake lets users write the variational formulation of their especially with magnetics only, where measurements of the problem in a domain-specific language. It also allows users to internal plasma profiles are not used. flexibly construct finite element spaces for use in Galerkin’s 3 Nucl. Fusion 64 (2024) 034001 C.J. Ham and P.E. Farrell agrees between the two methods to 0.9% or better. This excel- lent 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 solu- tions of the reduced GSE arrived at by completely different numerical methods and not numerical artifacts. An issue to note is that deflated continuation is not guar- anteed 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 with a dis- continuity, 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 ⩾ 0.1 there is excellent agreement between deflated con- tinuation 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. Figure 4. Bifurcation diagram for our toy model varying the parameter w with C = 2.9 and ψp = −1.0. The most negative solution is always there but at w ≈ 0.88 two new solutions appear at 4. Multiple solutions of Grad–Shafranov a fold bifurcation and persist as w decreases. In section 3 we showed excellent agreement between our ana- lytic model and finite element simulations for a plasma with method applied to this variational formulation 3. We will com- circular cross-section. We now apply the same computational pare the results from deflated continuation with our analytic techniques to solve the problem with a more complete geo- model in appropriate limiting cases. metry to understand the effects of aspect ratio and plasma The aim of the deflated continuation algorithm is to com- shaping on the multiple equilibria. In this section we will pute the solutions of an equation assume that the profile of F is a constant so that FF ′ = 0 and f (u, λ) = 0, (5) that the pressure has the profile we used previously, [ ( )] where u is the solution of a PDE and λ ∈ R is a parameter ′ C ψ − ψp p (ψ) = 1 − tanh. (6) on which the PDE depends. For each parameter value, each 2 w known solution branch is continued to that parameter, and This represents a plasma that has a constant gradient of pres- then the nonlinear problem is modified so that those solutions sure with flux up to the value ψ p and then transitions to zero are excluded (the solutions are deflated, just as known roots after that. The weak formulation of the problem that we now of polynomials can be deflated by dividing by the appropri- solve then becomes: find ψ ∈ H1 (Ω) such that ate factor). The algorithm then seeks solutions of the modified nonlinear problem; if any are found, they are new solutions of ˆ ˆ [ [ ( )]] 1 Cµ0 R ψ − ψp the original problem. Importantly, this approach can compute ∇ψ · ∇ξ dΩ = 1 − tanh ξ dΩ R 2 w disconnected bifurcation diagrams, like those already presen- ted in figures 3 and 4. for all ξ ∈ H1 (Ω) , (7) We now apply deflated continuation to our model since FF ′ = 0 if F is taken as a constant. The boundary con- equation (3). We take the weak form of (3) and approx- dition is that ψ = 0. imate its solutions with a piecewise linear finite element The difference between equation (7) and (the weak form of) method on a disk (radius, a = 2 and origin (0, 0) in two (3) is the 1/R term that appears in the left hand side of (7) and dimensions). We then use deflated continuation (as imple- also the factor of R that appears on the right hand side of (7). mented in the Defcon library4 ) to explore the solutions of These extra geometric terms break the poloidal symmetry of the problem. We have solved the cylindrical variant of the the plasma which exists in the cylindrical case but not in the problem, i.e. without toroidal effects, in these cases so R0 has toroidal case. no effect. We first study the effect of aspect ratio. We take a circular Figures 5–7 show three solutions that have been calculated cross section plasma of minor radius, a = 1, and change the for C = 2.9 with this procedure. The central value of the flux major radius R0 with a fixed value of w = 0.18 and R20 C = 6. The results are plotted in figure 8. We see that the aspect ratio 3 https://firedrakeproject.org. makes little difference to the number of solutions in this case 4 https://bitbucket.org/pefarrell/defcon. and also to their values. It has the strongest effect at small 4 Nucl. Fusion 64 (2024) 034001 C.J. Ham and P.E. Farrell Figure 5. First of the non-trivial solutions of (3) with C = 2.9 and w = 0.7. Figure 6. Second of the non-trivial solutions of (3) with C = 2.9 and w = 0.7. aspect ratio when the variation of R across the plasma cross Z = aκ sin (θ) , (9) section is relatively largest. The results of the two approaches converge at large aspect ratio. where θ is the poloidal angle. We used Gmsh to We next study the effect of plasma shaping at low aspect generate appropriate triangular meshes of these geomet- ratio. We parametrise the plasma shape with the standard ries. We show an example of the mesh and solution for a elongation, κ, and triangularity, δ as defined by plasma with triangularity δ = 0.3 and elongation κ = 2 in figure 9. We fix R0 = 3 for this computation and vary elonga- R = R0 + a cos (θ + δ sin (θ)) , (8) tion, figure 10 and triangularity, figure 11. We see that 5 Nucl. Fusion 64 (2024) 034001 C.J. Ham and P.E. Farrell Figure 7. Third of the non-trivial solutions of (3) with C = 2.9 and w = 0.7. Figure 8. Solutions of (7) with circular cross-section as the aspect ratio of the torus is varied. The central value of poloidal magnetic flux (ψ 0 ) for the three solutions is shown. plasma triangularity does not have a strong influence on as high elongation has other benefits. The strong effect of the number of solutions. However, for small enough elonga- increasing elongation on the lower solution is partly due to tion the two solution branches coalesce and disappear. This the plasma cross sectional area changing significantly with is not in a regime where tokamaks are normally operated elongation. 6 Nucl. Fusion 64 (2024) 034001 C.J. Ham and P.E. Farrell Figure 9. Example of solution mesh, generated with Gmsh, and solution, computed by Firedrake/Defcon, for a plasma with elongation κ = 2, triangularity δ = 0.3, and major radius R0 = 3. Figure 10. Scan of plasma elongation, κ, with fixed major radius R0 = 3 with fixed minor radius a = 1. The central value of poloidal magnetic flux (ψ 0 ) for the three solutions is shown. Two solution branches coalesce at around κ = 0.9 and only the trivial solution exists for smaller values of elongation. Figure 11. Scan of plasma triangularity major radius R0 with fixed minor radius a = 1 with circular cross section which amounts to a scan in aspect ratio of the torus. The central value of poloidal magnetic flux (ψ 0 ) for the three solutions is shown. 7 Nucl. Fusion 64 (2024) 034001 C.J. Ham and P.E. Farrell 5. Discussion and conclusions A further observation is that our analytic model is more accurate and robust than we might expect. This means that We have presented methods to calculate multiple tokamak we can use this method to explore a wide range of potential equilibria solutions of the GSE and demonstrate their use in profile shapes to find out when multiple solutions might be relatively simple cases. It should be noted that if there is non- important. We may also use simple transport models along linearity in the profiles used in the GSE then there are likely to with this analytic model to understand when equilibria may be multiple solutions. The results are generic in the mathemat- be lost. This could improve our understanding of some types ical sense. We have shown that as system parameters are varied of major disruption or how we can plot a path to higher per- we can lose equilibrium solutions. One immediate question is formance plasmas. what will happen to the plasma in this case. If the plasma is no We have assumed that the plasma here is axisymmetric, longer in an equilibrium state then it will evolve on an Alfvénic indeed this is required as we are using the GSE. However, if timescale until it reaches a new equilibrium, or transitions to we have a non-axisymmetric model of the plasma we could a periodic orbit, or to a chaotic orbit. In principle, this may compute not only additional axisymmetric equilibria but also result in a complete loss of the plasma, but this is not seen in non-axisymmetric states that could be thought of as saturated experiments where the initial stage of a disruption is the loss instabilities such as the helical core mode. This may also of the thermal energy, the thermal quench. The plasma that is be the topic of further work. left is much cooler and broadly force-free. It may be that some In this work we have answered a call to develop techniques types of major tokamak disruption are caused by such a loss to calculate multiple solutions to the GSE. We have provided of equilibrium. For example, if we follow the highest confine- both analytical and numerical techniques and shown that they ment equilibrium branch in figure 3 and we say the parameter agree in the correct limits. We have shown that even in a rel- C is continuously dropping (as may occur in experiment as a atively simple model that bifurcation points can occur. If we result of plasma transport) the equilibrium will be lost once have more complicated profiles of FF ′ and p′ then we would C < 2.7 and this would be experienced as a major disruption expect more of this behaviour. These tools can be exploited in experimentally. We will attempt to make more detailed com- future work to potentially improve our understanding of dis- parisons to experimental results in future work. ruptions and to to find routes to improved performance in toka- A key part of understanding tokamak plasmas is equilib- mak plasmas. rium reconstruction where diagnostic measurements are used to constrain an equilibrium model to infer the experiment- Acknowledgments ally realised equilibrium. There are various ways this is done, for various purposes. If only magnetic diagnostic data is used We are grateful to Wayne Arter, Andrew Kirk and Lucy Kogan then there are no measurements of the internal profiles. In for their comments on this work. C J H dedicates this work to this case highly simplified profiles are used for FF ′ and p′ Alan Ham and Kathleen Ham. which strongly constrains the potential shapes of the recon- This work has been funded by the EPSRC Energy structed equilibria. We have shown that two very different Programme [grant number EP/W006839/1]. To obtain fur- profile shapes can produce the same flux on axis. If internal ther information on the data and models underlying this paper measurements of the plasma are used, for example from MSE please contact [email protected] or Thomson Scattering diagnostics, in what is often called For the purpose of open access, the author(s) has applied a kinetic equilibrium reconstruction, then the profile shapes are Creative Commons Attribution (CC BY) licence to any Author given more freedom by using more basis functions. However, Accepted Manuscript version arising. there are more measurements to constrain the final equilib- rium. 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 Appendix that the choice of basis functions for equilibrium reconstruc- tion does not determine the conclusions of studies of the sta- In this appendix we demonstrate that the β and plasma cur- bility and transport of a given plasma. rent are relevant to practical situations for the large aspect ratio The effect of multiple equilibria should be evaluated in pre- case. The theory of cylindrical equilibrium is well known. We dictive modelling where an equilibrium model is coupled to a reproduce some results here for reference based on the βp ≈ 1, transport model. These two are evolved in time to simulate a ϵ ≪ 1 expansion of the GSE given in Freidberg. We start tokamak pulse, therefore ramping up the current to flat top and with the equations of static ideal MHD then ramping down. It may be imagined that in this process bifurcation points will be passed where an equilibrium model J × B = ∇p, (10) only finds one of the available equilibrium branches. It may be that an improved shot evolution can be found (or a deleterious ∇ · B = 0, (11) one avoided) if all bifurcation points are properly understood. The work done on theoretically predicting super-H mode may give indications of how this should be done. 5 See also https://doi.org/10.14468/yk04-h275. 8 Nucl. Fusion 64 (2024) 034001 C.J. Ham and P.E. Farrell 1 1 dψ0 J= ∇ × B. (12) Bθ0 (r) = , (21) µ0 R0 dr We assume cylindrical symmetry so all variables are inde- pendent of θ and z. We can then calculate the current density −F (ψ0 ) Bz0 (r) =. (22) R0 [ ] dBz (r) 1 d We have assumed that Bz is a constant so µ0 J = 0, − , (rBθ (r)) , (13) dr r dr rBz q (r) =. (23) and so the radial component of the static momentum equation ψ0′ becomes [ ( 2 )] d Bθ + B2z B2 p+ =− θ , (14) A.2. Toroidal current profile dr 2µ0 µ0 r We need an expression for the toroidal current J ϕ which can where we may specify two of the profiles, i.e. Bθ and p(r) with be obtained from the derivation of the GSE (and in the βp 1, the third, Bz (r) in this case, being determined by these two. ϵ ≪ 1 expansion of the GSE given in Freidberg ) We may multiply equation (14) by π r2 and integrate from 0 to a to get dp 1 dF 2 − µ0 R0 Jz (r) = ∇2 ψ = −µ0 R20 − , (24) ⟨ ⟩ dψ 2 dψ B2z B2 (a) B2θ (a) ⟨p⟩ + − z − = 0, (15) 2µ0 2µ0 2µ0 and so in our case, given we have taken F as a constant, where [ ( )] dp R0 C ψ − ψp ˆ a Jz (r) = R0 = 1 − tanh. (25) 1 dψ 2 w ⟨X⟩ = 2π rX (r) dr (16) π a2 0 is the average value of X(r) over a cross-section of radius a. A.3. Equilibrium scaling The ratio of the thermal pressure to the magnetic pressure is the plasma β: The equilibrium can be scaled following the rules given in Lutjens et al. These allow us first to specify the value 2µ0 ⟨p⟩ of the q profile at one radial location, therefore allowing us to ⟨β⟩ =. (17) B2 (a) 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 β. The poloidal and toroidal parameters are defined as We cannot alter the poloidal β or the internal inductance with these scalings. The first rescaling is 2µ0 ⟨p⟩ 2µ0 ⟨p⟩ βp = , βt =. (18) B2θ (a) B2z (a) ψnew = α1 ψold , Fnew = α1 Fold , pnew = α12 pold , (26) We have assumed for simplicity in this paper that Bz is a con- stant so (15) results in βp = 1. and the second rescaling is F 2 new = F 2 old + α2. (27) A.1. Safety factor profile The twist of the magnetic field lines is captured by the safety We use (27) first to pick the q-profile and then (26) is used factor which in a cylindrical screw pinch is to set the β t and toroidal current. We have used the scalings to produce a set of equilibria with pressure profiles shown in rBz (r) figure 12 and safety factor profiles shown in figure 13. These q (r) =. (19) show that the pressures and safety factor profiles for these R0 Bθ (r) simple cases are not pathological in any sense. The equilib- In this approximation we have rium scaling shows that suitable equilibria can always be found and that when a bifurcation point exists the two equilibria can Br0 (r) = 0, (20) always be scaled to be sensible. 9 Nucl. Fusion 64 (2024) 034001 C.J. Ham and P.E. Farrell Figure 12. Plot of the pressure profiles for the three solutions of the large aspect ratio GSE with C = 2.9. Figure 13. Plot of the safety factor profiles for the three solutions of the large aspect ratio GSE with C = 2.9. ORCID iD Schnack D.D. 2009 Lectures in Magnetohydrodynamics (Springer) C.J. Ham  https://orcid.org/0000-0001-9190-8310 Farrell P.E., Birkisson Á. and Funke S.W. 2015 Deflation techniques for finding distinct solutions of nonlinear partial differential equations SIAM J. Sci. Comput.3 37 A2026–45 References Farrell P.E., Beentjes C.H.L. and Birkisson Á. 2016 The computation of disconnected bifurcation diagrams Solano E.R. 2004 Criticality of the Grad–Shafranov equation: (arXiv:1603.00809) transport barriers and fragile equilibria Plasma Phys. Rathgeber F., Ham D.A., Mitchell L., Lange M., Luporini F., Control. Fusion 46 L7–L13 Mcrae A.T.T., Bercea G.-T., Markall G.R. and Kelly P.H.J. Berestycki H. and Brezis H. 1980 On a free boundary problem 2016 Firedrake: automating the finite element method by arising in plasma physics Nonlinear Anal. Theory Methods composing abstractions ACM Trans. Math. Softw. 43 1–24 Appl. 4 415–36 Goedbloed J.P., Keppens R. and Poedts S. 2010 Advanced Jeyakumar S., Pfefferlé D., Hole M.J. and Qu Z.S. 2021 Magnetohydrodynamics (Cambridge University Press) Analysis of the isotropic and anisotropic Grad-Shafranov Freidberg J. 2014 Ideal MHD (Cambridge University Press) equation J. Plasma Phys. 87 905870506 Ham D.A. et al 2023 Firedrake User Manual 1st edn (Imperial Taylor J.B. 1986 Relaxation and magnetic reconnection in College London and University of Oxford and Baylor plasmas Rev. Mod. Phys. 58 741–63 University and University of Washington) 10 Nucl. Fusion 64 (2024) 034001 C.J. Ham and P.E. Farrell Geuzaine C. and Remacle J.-F. 2009 Gmsh: a Cooper W.A., Graves J.P., Pochelon A., Sauter O. and three-dimensional finite element mesh generator with Villard L. 2010 Tokamak magnetohydrodynamic built-in pre- and post-processing facilities Int. J. Numer. equilibrium states with axisymmetric boundary and a 3D Methods Eng. 79 1309–31 helical core Phys. Rev. Lett. 105 035003 Snyder P.B. 2019 High fusion performance in super H-mode Lutjens H., Bondeson A. and Sauter O. 1996 The CHEASE experiments on Alcator C-Mod and DIII-D Nucl. Fusion code for toroidal MHD equilibria Comput. Phys. Commun. 59 086017 97 219–60 11