We present a Python-based numerical solver for the two-dimensional dynamic plasma equilibrium problem. We model the time evolution of toroidally symmetric free-boundary tokamak plasma equilibria in the presence of the non-linear magnetohydrodynamic coupling with both currents in the “active” poloidal field coils, with assigned applied voltages, and eddy currents in the tokamak passive structures. FreeGSNKE (FreeGS Newton–Krylov Evolutive) builds and expands on the framework provided by the Python package FreeGS (Free boundary Grad–Shafranov). FreeGS solves the static free-boundary Grad–Shafranov (GS) problem, discretized in space using finite differences, by means of Picard iterations. FreeGSNKE introduces: (i) a solver for the static free-boundary GS problem based on the Newton–Krylov (NK) method, with improved stability and convergence properties; (ii) a solver for the linearized dynamic plasma equilibrium problem; and (iii) a solver for the non-linear dynamic problem, based on the NK method. We propose a novel “staggered” solution strategy for the non-linear problem, in which we make use of a set of equivalent formulations of the non-linear dynamic problem we derive. The alternation of NK solution steps in the currents and in the plasma flux lends this strategy an increased resilience to co-linearity and stagnation problems, resulting in favorable convergence properties. FreeGSNKE can be used for any user-defined tokamak geometry and coil configuration. FreeGSNKE's flexibility and ease of use make it a suitably robust control-oriented simulator of plasma magnetic equilibria. FreeGSNKE is entirely written in Python and easily interfaced with Python libraries, which facilitates machine learning based approaches to plasma control.
I. INTRODUCTION
Magnetohydrodynamic (MHD) equilibrium calculations play a foundational role in the analysis of tokamak plasmas, both in terms of the interpretation of experiments on the current fleet of devices1 and for the design of the next generation of prototype reactors. A plethora of codes exist to address the varied equilibrium related tasks, ranging from codes like EFIT++,2 which allow inference of possible magnetic equilibria consistent with measurements from diagnostics like magnetic probes, motional Stark effect, and Thomson scattering systems, to codes like FIESTA, which are routinely used for power plant design studies.3 As well as static equilibrium simulations, codes exist which couple the Grad–Shafranov (GS) equation to a system of circuit equations describing the passive components and active coils in the tokamak and allow a full evolutive prediction of the trajectory of the plasma equilibrium over time.4 One key application of evolutive simulations is in the design of robust, performant control schemes capable of guiding the plasma through ramp-up, to the desired steady state operating point and then to a safe termination.
Recent work has demonstrated how machine learning is a useful tool to designing control schemes, with contributions to plasma equilibrium reconstruction and emulation,5–7 plasma detachment studies,8–10 disruption prediction,11,12 and magnetic stabilization and nonlinear control of the plasma position and shape.13 These examples highlight fresh challenges faced by the magnetic-confinement community, across different machines, as well as indicate a set of ideal criteria that a plasma evolution codebase should satisfy to be conducive to contemporary research in real-time control and integrated modeling. First, code to model the plasma dynamics should allow for a stateful representation of the plasma, i.e., a time-dependent description in terms of a given number of parameters, allowing a restart of any simulation from any intermediate state with ease. Second, such codebase should allow for a seamless co-integration between equilibrium calculations, dynamics, evolution of internal profiles due to transport, etc., and machine learning (ML) frameworks. These two criteria are of particular importance to nonlinear control studies based on emulation and/or reinforcement learning, where exploration and experience accumulation are based on numerous, repeated interactions between the ML framework itself and the simulation one. Third, any algorithmic and numerical implementation choices should be fully and easily inspectable, to facilitate contribution and customization by the user as well as to avoid duplication of efforts. Close-sourced languages do not appear as the prime candidate to promote this latter but fundamental point. A Python codebase which can make the most of open-source libraries would be better suited to satisfy all criteria above, but it is currently missing in the literature. This is becoming ever more pressing, for example, as transport modules for integrated modeling are adopting machine-learning emulators14 or just-in-time (JIT) compiling modules, the largest majority of which are written in Python.
The optimization of scenario and control design requires fast and accurate simulation capabilities. Without aiming to compile a complete list of all of the evolutive codes available in the literature, the diversity of codes being currently used for similar purpose is exemplified by DINA,15 TokSys/GS-Evolve,16 CREATE-NL,17 SPIDER,18 NICE,19 CEDRES++,20 and FGE.21 Each of these codes have specific features, they have different levels of complexity and computational costs, and were originally aimed at tackling different applications. For example, while CREATE-NL and FGE adopt solution strategies based on numerically approximated gradients, DINA solves the free boundary equilibrium problem using an inverse variable technique,22 and CEDRES++ uses a finite element method for which it is able to derive analytic gradients.20 Additional differences involve the coding languages themselves: TokSys, CREATE-NL, and FGE are Matlab/Simulink based codes, while NICE and CEDRES++ are in C++. While each of the codes above has been used in the literature with success, many of them use licensed languages. For some, it may be difficult to inspect and customize algorithmic choices to the needs of specific applications. Others do not actually feature a stateful representation of the plasma. For instance, interfacing the codes above with open-source reinforcement learning libraries or with transport modules based on emulation or JIT compilation would require redundant read-write operations and/or use of purpose-made socket interfaces for multi-language interaction.23 In addition to the latency and memory usage that these workarounds may cause whenever repeated calls to an evolutive solver are needed, the need for cross-language binders also introduces an unnecessary barrier to widening the community effort toward novel approaches, e.g., in integrated modeling and nonlinear control.
A. This work
In this work, we introduce a new evolutive code, FreeGSNKE (pronounced free-ge-snake and standing for FreeGS Newton–Krylov Evolve), which aims to satisfy the criteria outlined above. FreeGSNKE models the toroidal dynamics of tokamak plasmas and is entirely written in Python. This facilitates any user inspection, customization, and contribution as well as allows for seemless integration with emulatated transport modules and with any open-source Python-based ML libraries and frameworks. FreeGSNKE is also a stateful code, which benefits control design and reinforcement learning studies. Furthermore, this paper introduces and explores novel algorithmic strategies to solve the non-linear coupling of currents in the plasma and in the confining structures across a sequences of Grad–Shafranov equilibria.
One notable Python-based tool in equilibrium modeling is FreeGS,24 a freely available package which solves the GS equation for free boundary equilibria. FreeGS focuses on “interpretative” tasks, in which magnetic measurements are used to constrain the equilibrium itself. As such, FreeGS does not feature a “forward” GS solver (see Sec. II B) and it does not tackle the dynamic evolution problem (see Sec. III). The code we present here, FreeGSNKE, provides these added capabilities, by building on the FreeGS framework.
Like the evolutive codes mentioned above, FreeGSNKE is a quasi-static plasma evolution code, in which the plasma evolves across solutions of the GS equation. This working hypothesis stems from the realization that both (i) the Alfvén timescale, on which the plasma relaxes toward the lowest-order force balance and (ii) the parallel transport timescale, on which density and temperature homogenize on flux surfaces, are much shorter than both (a) the energy confinement timescale and (b) the characteristic timescales of the coupled circuit equations which describe the tokamak conductors. FreeGSNKE couples the plasma with these conductors, by pairing the GS equation with the system of circuit equations describing both active poloidal field coils and eddy currents in the tokamak passive structures. Given the time dependence of the applied voltages on the active poloidal field coils (including any Ohmic drive), FreeGSNKE calculates the evolution of the plasma equilibrium and the evolution of active and passive currents. At the present time, FreeGSNKE does not account for ferromagnetic structures in the tokamak.2 It is also worth pointing out that FreeGSNKE in and by itself does not include a transport code, to model self-consistently the transport of density or plasma energy across flux surfaces. The evolution of the parameters defining the dimensionless profile is prescribed, and FreeGSNKE treats it as given input. The absolute scaling of the internal profiles is given by the total plasma current, which is solved inside FreeGSNKE by enforcing Poynting's theorem on the plasma together with the circuit equations on the conducting structures and with Grad–Shafranov equilibrium sequences. Fully self-consistent evolution may be achieved by coupling FreeGSNKE with a transport code. Already in its current form, FreeGSNKE allows for a reparameterization of the internal profiles, so that it can be easily coupled with transport modules in a 1.5D modeling fashion,25,26 an approach which is common to the codes referenced above. This is beyond the scope of this paper, which is devoted to detailing the algorithms developed for the coupled evolution of currents and equilibria, and we defer it to future work.
This paper presents a description of the main computational methods adopted in FreeGSNKE and provides an illustration of its functionalities. Section II focuses on the free-boundary static GS problem: it describes FreeGSNKE's plasma profile conventions as well as similarities with and differences with FreeGS. It then validates FreeGSNKE's static solver on analytical Solov'ev equilibria. Section III introduces the dynamic problem, and describes how we model the coupling between metal structures and plasma. Section IV details our solution methods and implementation choices and exemplifies the use of both our linear and nonlinear solvers. In particular, it describes both our reformulations of the dynamic problem and our novel solution strategy. The Appendix provide additional details on: the NK method (Sec. 1 of Appendix), which we use in both our static GS solver and in the solver of the non-linear dynamic problem; the normal mode decomposition we adopt for the vessel structures (Sec. 2 of Appendix); the inner workings of our solution strategy (Sec. 3 of Appendix).
II. STATIC GS EQUILIBRIA
Here, we assume that FreeGSNKE has been provided with a suitable description of the tokamak's geometry, properties of the active poloidal field coils, and an axisymmetric toroidal model of its passive conducting structures.
A. Instantaneous states of the plasma-tokamak system
FreeGSNKE is a stateful simulator: at each time t, the plasma-tokamak system is described by its state S(t), which include the following quantities:
-
a set of current values , featuring all currents in both active and passive tokamak metal structures, ;
-
the value of the total current in the plasma, ;
-
a set of parameters describing the plasma profile, , the details of which depend on the chosen parameterization.
-
are cylindrical coordinates and μ0 is the magnetic permeability.
-
is the poloidal flux function, i.e. the poloidal magnetic flux divided by . This includes a contribution by the plasma, ψp, and a contribution associated to currents in any external metal conductors, ψm, accounting for both active circuits and passive structures. Therefore, .
-
is the plasma toroidal current density distribution, is the toroidal current density distribution of metal conductors. In absence of ambiguity, we can drop the subscript .
-
The user sets values for the pair directly, for example, to qualitatively regulate the peaked-ness of the plasma current distribution or to adjust the value of the plasma safety factor q at some desired set of locations in the plasma.
-
For a fixed pair , the product is in a one-to-one relationship with both (i) the ratio between the plasma pressure and the poloidal magnetic pressure, βp and (ii) the plasma pressure on the magnetic axis, pa. The product can, therefore, be set indirectly, by requiring a specific value for either βp or for pa.
-
Finally, the value of the total plasma current Ip provides the outstanding constraint to disentangle the individual parameters λ and β0.
B. The forward GS problem
We focus on the forward static GS problem, which can be conveniently stated in terms of the plasma-tokamak state S: Given the state S, including values for all currents in the conductors Im, the total plasma current Ip, and a valid set of plasma profile parameters , find the plasma contribution to the flux function, , which satisfies the GS equation (1)
C. Benchmark with an analytic solution of the static problem
A comparison between a numerical solution of the static GS problem and the corresponding analytic Solov'ev solution for a simple plasma in a limiter configuration (the dotted black line shows the limiter while the black × markers display the positions of the active coils). The left panel displays isoflux contours for both numerical solution (orange dashed lines) and analytic Solov'ev solution (black curves), for the same set of flux values. The right panel demonstrates and quantifies the accuracy of the numerical solution by displaying values of the residue (12) inside the plasma domain.
A comparison between a numerical solution of the static GS problem and the corresponding analytic Solov'ev solution for a simple plasma in a limiter configuration (the dotted black line shows the limiter while the black × markers display the positions of the active coils). The left panel displays isoflux contours for both numerical solution (orange dashed lines) and analytic Solov'ev solution (black curves), for the same set of flux values. The right panel demonstrates and quantifies the accuracy of the numerical solution by displaying values of the residue (12) inside the plasma domain.
Here, Ri and Ro are the radial coordinates of the equatorial innermost and outermost points; (Rt, Zt) is the highest point on the separatrix; and ψb is the value of the flux function at the boundary according to our numerical solution ψNK. Black lines in the left panel of Fig. 1 show isoflux contours for the analytic solution ψS in the plasma domain, for the same set of flux values used for ψNK. Curves overlap for the entire set of values. Within the domain of validity of the Solov'ev analytic solution, ψS and ψNK are hardly distinguishable.
Finally, it is worth mentioning that we have proceeded to verifying that FreeGSNKE's forward static solver agrees with FreeGS' inverse static solver.
III. THE EVOLUTIVE PROBLEM
We can now formulate the evolutive equilibrium problem, which we first define in terms of the states of the plasma-tokamak system S. Assuming a time step dt and given
-
a plasma-tokamak state S(t) like described in Sec. II A;
-
a prescription for the parameters of the dimensionless profile function at time t + dt, ;
-
the set of voltage values applied to the active poloidal field coils between time t and t + dt.
FreeGSNKE solves for (i) all metal current values and (ii) the total plasma current . Together with the prescribed set , these complete the tokamak-plasma state , thereby providing a full characterization of the equilibrium at time t + dt. The following sections detail the conditions that the state have to satisfy and describe the coupling between plasma and metal structures.
A. Defining the dynamic plasma equilibrium problem
1. The metals
-
is the matrix of mutual inductances, describing the coupling between each pair of metal circuits;
-
is the diagonal resistance matrix, including resistance values for all tokamak conductors;
-
Vm is the vector of applied voltages; by definition, the voltage applied on passive tokamak structures is zero, so that , where Vact is the set of active voltage values applied to the active poloidal field coils;
-
is the vector of electromotive forces on each tokamak conductor due to the change in plasma magnetic flux, which represents the nonlinear coupling between plasma and metals. The integrated fluxes should not be confused with the plasma flux ψp.
In FreeGSNKE, it is possible to simplify the description of the passive vessel structures using normal modes: rather than solving directly for the evolution of the eddy currents Ipas, we introduce the basis of currents which diagonalizes the dynamics of the passive structure currents. The advantage of doing so is that normal modes which either (i) evolve on exceedingly fast timescales or (ii) have limited effect on the plasma, may be dropped as desired as their dynamics is near-decoupled. Section 2 of Appendix lays down the details of the decomposition in normal modes we adopt in our implementation, here, we just report that are the eigenvectors of the matrix , where Rpas and are the resistance and mutual inductance matrices relevant to the passive structures. The characteristic frequencies are the associated eigenvalues. Normal modes with fast dynamics are easily dropped on the basis of the frequencies . To identify modes with a weak coupling with the plasma, we examine the response matrix , and in particular, we consider the norm of its columns—response matrices are defined in Sec. IV A, while the vector of currents Id is defined in Sec. 2 of Appendix. An illustration of the use of these two criteria is provided in Sec. IV A 1.
2. The plasma
-
is the matrix of mutual inductances between pairs of grid points in the discretized spatial domain on which Iy is defined;
-
is the diagonal matrix of plasma resistance values;
-
, with defined in Sec. III A 1.
We note that without specific adjustments, the dimensionality of Iy is the same as the dimensionality of the Discretized spatial domain, which may be superfluous. Depending on the tokamak's geometry, some regions may not be accessible to the plasma, which encourages the use of a suitably reduced spatial domain when defining Iy—and consequently, and . In FreeGSNKE, we limit Iy to the region defined by the limiter.
3. Basic problem formulation
It is now possible to define the quasi-static evolutive problem. This consists of satisfying at the same time the circuit equations (15) and the lumped Pointing equation (18) while forcing the plasma on solutions of the GS equation (1) FreeGSNKE approaches this problem by discretising it in time and by approximating time derivatives with forward finite differences.
It is also worth mentioning explicitly that one may tackle the problem described by the system of Eqs. (21) and (22) in two different ways
-
by seeking unrestricted trial solutions over the entire vector space of or
- by restricting the search to the manifold of GS solutions, i.e. only considering trial solutions of the following form:
where is the solution of the GS equation for the extensive currents Ie, as defined in Eq. (4). In this second case, Eq. (22) is satisfied by construction and the only actual unknown is the vector of extensive currents .
FreeGSNKE adopts solution methods belonging to both of these two classes, which we detail in Sec. IV.
IV. STRATEGIES, IMPLEMENTATION, AND EXAMPLES
A. Linearized problem
The linearized dynamic problem obtained by replacing Eq. (29) in Eq. (21) is a readily solvable linear system in the extensive currents. In FreeGSNKE, this linear system is solved using an implicit Euler method, as the coupled dynamics is notoriously stiff. If desired, the evolution of the plasma flux can be obtained retroactively by solving the associated set of forward GS problems once the evolution of the currents has been calculated. Clearly, use of the response matrices calculated for S = S0 remains accurate only for equilibria S(t) that remain suitably close to the expansion point itself. Caution should be exercised whenever calculating a “linearized time evolution” using fixed values for the response matrices, as results may substantially degrade when any differences between S0 and S(t) grow as part of the evolution. (see, for example, Sec. IV B 3).
1. Example: The growth-rate of a vertical displacement event
Next, we illustrate using FreeGSNKE's linearized evolutive solver capabilities by calculating the growth rates of a vertical displacement event (VDE) in a diverted plasma, focusing in particular, on the effect of removing fast evolving and weakly coupled normal modes.
As in Sec. II C, we adopt a 129 × 129 spatial grid. We take inspiration from the MAST Upgrade (MAST-U) tokamak and adopt a simplified geometry that qualitatively reproduces some of its main characteristics. This is displayed in the left-most panel of Fig. 2. Black crosses display the location of the poloidal field coils. Blue points show the location of the passive filaments we use to model the tokamak passive structures, for a total of 259 toroidally symmetric elements. The geometry of the limiter wall is also shown. We stress that we do not aim to provide an accurate model of MAST-U here and that both tokamak and initial equilibria should be intended as part of a toy problem setting.
A diverted plasma in a toy tokamak inspired by the geometry of MAST-Upgrade. The left-most panel shows the tokamak geometry: black crosses indicate poloidal field coils and blue dots represent the toroidal filaments we use to model the vessel passive structures. Orange lines and the cyan separatrix are isoflux contours for two up-down symmetric plasmas. Black plus signs in the middle panel show the number of vessel normal modes Nd with characteristic frequency smaller than . Orange and cyan symbols show the number of modes with and with non-negligible coupling with the plasma (see Sec. IV A 1 for more details). Black symbols in the right-most panel display the e-folding time of the vertical instability for both plasmas displayed in the right-most panel, as a function of the truncation frequency . Colored symbols refer to only including modes with non-negligible coupling with the plasma (see Sec. IV A 1 for more details).
A diverted plasma in a toy tokamak inspired by the geometry of MAST-Upgrade. The left-most panel shows the tokamak geometry: black crosses indicate poloidal field coils and blue dots represent the toroidal filaments we use to model the vessel passive structures. Orange lines and the cyan separatrix are isoflux contours for two up-down symmetric plasmas. Black plus signs in the middle panel show the number of vessel normal modes Nd with characteristic frequency smaller than . Orange and cyan symbols show the number of modes with and with non-negligible coupling with the plasma (see Sec. IV A 1 for more details). Black symbols in the right-most panel display the e-folding time of the vertical instability for both plasmas displayed in the right-most panel, as a function of the truncation frequency . Colored symbols refer to only including modes with non-negligible coupling with the plasma (see Sec. IV A 1 for more details).
We build the matrix of mutual inductances relevant to metal elements, , based on the geometry of our simple tokamak (and some simple assumption on the size of the poloidal section of each coil). For the resistance values , we take that that poloidal field coils are made of copper, while passive elements are made of steel. As described in Sec. III A 1 and Sec. 2 of Appendix, the matrices and define the normal modes Id of the vessel. Black plus signs in the middle panel of Fig. 2 display values of the characteristic frequencies ωd of these modes, by showing the number of modes Nd with characteristic frequency smaller than the frequency , on the abscissa. As it is well known in this setting and mentioned in Sec. III A 1, the frequencies ωd span several orders of magnitude, and reach down to timescales that are considerably shorter than those we wish to resolve here. Depending on the application, a suitable threshold frequency may be selected and used to drop modes with .
We also illustrate here the effect of selecting modes based on the strength of their coupling with the plasma. As touched upon in Sec. III A 1, we base this selection on the norm of columns in the response matrix . Orange crosses in both middle and right-most panel of Fig. 2 refer to including an additional criterion based on how strongly each mode actually affects the plasma. As touched upon in Sec. III A 1, we base this on the norm of columns in the response matrix . For reference, here, we use . As shown by the orange crosses in the central panel, this additional criterion substantially reduces the number of selected modes, with at most 48 independent modes actually being retained, independently of the value of . Orange crosses in the right-most panel of Fig. 2 show how the estimate of the VDE timescale τ is affected by this simplification. Despite more than 200 normal modes having been dropped the estimate remains effective. In particular, for KHz, the VDE growth rate remains consistent at better than . Henceforth, this is the selection criterion we adopt in the following.
We perform an analogous parallel analysis on a second equilibrium with the same set of properties as in Eq. (30), but with a different plasma shape, as displayed by the cyan line in the left-most panel of Fig. 2. This equilibrium features a considerably faster vertical instability, with an e-folding timescale shorter by approximately two orders of magnitude with respect to the equilibrium considered above, as displayed by the black crosses in the right-most panel of the same figure. This may be at least partly justified by a poorer coupling with the vessel passive structure, as a consequence of the more limited volume this second plasma fills within the limiter. Cyan diamonds in both middle and rightmost panel refer to only including modes with .
B. Non-linear problem
FreeGSNKE implements two different solution strategies for the fully non-linear evolution problem. Both approaches use NK solution methods, but these are applied to different, though equivalent, reformulations of the non-linear evolution problem, as well as using different algorithmic choices. Next, we examine these mathematical reformulations, then, in Sec. IV B 2, we focus on our algorithmic solution strategies, illustrating how they lead to good convergence properties. To the best of our knowledge, these reformulations and strategies are original contributions of this paper.
1. Equivalent reformulations of the evolutive problem
We formulate here two root problems that are equivalent to the combined system of circuit Eq. (21).
Figure 3 illustrates the equivalence between the root problems (21), (33), and (35) by showing the relative residuals corresponding to each of these three reformulations as they evolve toward convergence, during the solution of an individual time step of plasma evolution. As done for Eq. (21), we define the relative residual corresponding to Eq. (33), and the relative residual corresponding to Eq. (35), by using the current and plasma flux updates defined in Eq. (23), respectively. We also use such relative residuals to define relative convergence criteria analogous to Eq. (25). In Fig. 3, the same set of trial solutions is used to calculate each of the displayed relative residuals at each iteration step. All three relative residuals evolve coherently: an improvement in one corresponds to improvement in the others. For completeness, trial solutions at each iteration in Fig. 3 are obtained using our staggered solution strategy, presented in Sec. IV B 2. Figure 3 displays data from a time step in the VDE presented in both Secs. IV A 1 and IV B 3.
A demonstration of the equivalence between the standard circuit equations (21) leading to the relative residual , and our reformulations (33) and (35), leading to the relative residuals and , respectively. The figure displays the evolution of such residuals toward convergence, as a function of the iteration step.
A demonstration of the equivalence between the standard circuit equations (21) leading to the relative residual , and our reformulations (33) and (35), leading to the relative residuals and , respectively. The figure displays the evolution of such residuals toward convergence, as a function of the iteration step.
Having at our disposal a set of equivalent formulations of the evolutive problem is advantageous when solving the nonlinear problem numerically, especially as it is well known that NK methods can be prone to co-linearity and stagnation.38,39 When the NK method stagnates or fails due to collinearity on one of the root problems, the shift to an alternative formulation is likely to be advantageous. In the following, we show how use of the iterated root problem (33) leads to the reduction of problematic co-linear updates to the Krylov subspace (see Sec. 1 of the Appendix) when compared to using the standard circuit equations (21) We also show how a novel, “staggered” solution strategy adopting both our root problem reformulations (33) and (35) leads to a further improvement in the convergence properties.
2. Solution strategies
We compare here three solution strategies. Two of them are of the restricted kind described in Sec. III A 3, seeking solutions parameterized by the extensive currents alone and located on the manifold of the GS solutions in the space of all trial solutions . The third, and our best in terms of convergence properties, is of the unrestricted kind.
The first two solution strategies are elementary: we simply apply the NK method to the root problems (21) and (33), respectively, starting from a guess of the extensive currents provided by a solution of the linearized problem.40 As such, a comparison between the performances of these two is a comparison between adopting the standard circuit equations (21) leading to what we refer to with Strategy Ia (SIa) or adopting our equivalent formulation (33) leading to Strategy Ib (SIb).
We use both SIa and SIb to solve for the nonlinear VDE set out in Sec. IV B 3, based on the plasma described in Sec. IV A 1. The top-left panel in Fig. 4 displays the evolution of the relevant relative residuals for an individual example time step of such nonlinear evolution, as a function of the solution step. For reference, all time steps in Fig. 4 amount to of the e-folding time of the VDE. The relative residuals relevant to SIa are displayed in orange, while those relevant to SIb are shown in purple. The behavior displayed in the referenced panel is representative of a majority of individual time steps in the VDE: we find it challenging to reach good levels of convergence when using SIa. We find that SIa suffers from frequent co-linear updates to the Krylov subspace which seriously inhibit the progress of the NK iterations. In turn, and despite an otherwise identical setting, SIb is a successful strategy. We find that use of the reformulation (33) substantially decreases the co-linearity and stagnation properties of SIa, thereby leading to systematically superior convergence properties. We note that, as mentioned in Sec. 1 of Appendix, in both SIa and SIb we orthogonalise the Krylov basis. It would be interesting to explore whether the reformulation (33) is somehow “less nonlinear” than the standard root problem (21) and, therefore, easier to tackle using a NK approach, although we do not attempt this here. All in all, we find that SIb is substantially less prone to stagnation and capable of delivering high accuracy in the largest majority of cases. For these reasons, we do not consider SIa further in this paper and, while we implement SIb in FreeGSNKE, we do not do so for SIa.
A comparison of the convergence properties of our solution strategies for the nonlinear evolutive problem. Each panel displays the relevant relative residual as a function of the solution step for an individual time step, from a selection of time steps in the evolution of the VDE described in Sec. IV B 3. Our novel staggered solution strategy is SII. See the text for details.
A comparison of the convergence properties of our solution strategies for the nonlinear evolutive problem. Each panel displays the relevant relative residual as a function of the solution step for an individual time step, from a selection of time steps in the evolution of the VDE described in Sec. IV B 3. Our novel staggered solution strategy is SII. See the text for details.
We now introduce our Strategy II, which make use of both our reformulations (33) and (35) of the dynamic problem. We set up SII as a staggered approach, by applying the NK method separately and sequentially to the extensive currents and to the plasma flux , by using, respectively, Eqs. (33) and (35). This is different, for example, from the implementation used by Carpanese et al.,37 which appears to use the NK method directly on the composed root problem , thereby building updates for the entire vector on unknowns at once. Intuitively, it can be expected that alternating NK solution steps in and should lead to a reduction in stagnation problems: each of the two different and subsequent applications of the NK method are viable to propose a useful step toward the solution if/when the other has been unsuccessful. At the same time, if/when one of the two NK applications should get stuck in a co-linear cycle,39 a step by the other NK method would likely solve this, thereby improving issues caused by possible co-linear updates.
Finally, we define SII as the following algorithm:
-
Use the solution of the linearized problem (as defined in Sec. IV A) and the associated GS plasma flux, , as initial trial solution: .
-
If the current trial solution fails the desired relative convergence criterion for the GS residual equation (26) update using Eq. (36).
-
If is such that the relative convergence criterion set for eq. (35) fails, update by applying the NK method to the same root problem.
-
If is such that the convergence criterion set for Eq. (33) fails, update by applying the NK method to the same root problem.
-
Check the current trial solution against the relative convergence criteria set for both Eqs. (33) and (6). If either are not satisfied, restart from point 2.
As done with SIa and SIb, we use SII to solve for the same VDE as above. Data displayed in black in Fig. 4 pertains to SII. Each panel in that figure refers to a different time step, selected across the entire VDE, and displays the evolution of the relative residual associated with the root problem (26) as a function of the solution step. For SII, a solution step refers to a cycle in the algorithm outlined above. For completeness, we report that in this evolution we have used . This means that the residual on each independent component of the iterated currents is less than one part in ten thousand of the corresponding update in the same time step [see Eqs. (23) and (25)]. Similarly, the GS residual is smaller than one part in ten thousand of the update to the plasma flux due the time evolution in that time step [see Eqs. (23) and (26)].
Since SIb and SII both rely on the same iterated root problem (33) the residuals displayed for the two strategies are directly comparable in Fig. 4. Such comparison shows that SII has improved convergence properties when compared to SIb, which, at times, appears to include updates to the extensive currents that in fact lead to an increase in the maximum and (sometimes, though less often) the average relative residual. The most plausible cause of these erroneous updates is, again, co-linearity between Krylov projections vl (see Sec. 1 of Appendix), resulting in difficulties in the least squares problem (A4). In turn, the evolution of residuals obtained using SII is considerably better behaved, thanks to the alternation between NK updates to the extensive currents and the plasma flux. This gain is mirrored in the number of calls to our GS solver that are necessary to reach convergence in each individual time step. These numbers are displayed in the upper-right of each panel in Fig. 4, in black for SII and in purple for SIb. SII requires a systematically reduced number of calls to the GS solver. For these reasons, SII is our preferred solution strategy in FreeGSNKE.
3. Example: Linear vs nonlinear VDE
We provide here some additional details on the solutions produced for the VDE which we used in Sec. IV B 2. This starts with the first of the two equilibria described in Sec. IV A 1. In order to define a unique time evolution, we need to define a few more parameter choices.
-
We fix the resistivity of the plasma σp to the value . For reference, this results in a lumped plasma resistance of for our chosen equilibrium.
-
In absence of a self-consistent evolution of the profile parameters and of the value of the plasma pressure on axis pa, we prescribe constant values for these quantities, and keep them fixed as in Eq. (30) during the evolution.
-
Voltages applied to the active poloidal field coils are not driven by a control policy during the evolution. Instead, we keep voltage values constant at the values set by the initial equilibrium: , where is the resistance matrix of the poloidal field coils and describes the current values in the initial conditions.
With the choices above, we simulate the dynamical evolution of our reference equilibrium during the entire VDE, up to the time at which the plasma contacts the wall, which we find takes a total of approximately 45 ms.
We solve both the linearized problem and the complete non-linear problem, for which we use our SII. For the former, we use the response matrices (29) defined by our initial equilibrium and keep them constant during the evolution. As mentioned in Sec. IV A, it can be expected that this approximation grows less and less accurate with time. For the nonlinear solution we use time steps of size dt = 0.35 ms and we use relative convergence criteria as outlined in Sec. IV B 2. A selection of the results from the two different analyses is displayed in Fig. 5. Dashed lines refer to the linearized evolution, full lines pertain to the solution of the complete non-linear problem. The upper set of three panels displays snapshots of the separatrix during the VDE. The time evolution of the z coordinate of the plasma magnetic axis, ZO, is also shown in one of the lower panels in the same Figure. The approximation provided by the linearized evolution appears accurate while ms, corresponding to ZO of a few mm. After that, the speed of the VDE accelerates, but this is not captured by the solution of the linearized problem, which, by construction, continues to track the linear e-folding time τ. This is evidenced by the black dash-dotted line in the same panel, which shows the slope of the corresponding exponential growth as derived in the previous Section. Similarly, the linearized evolution underestimates the speed of the decay of the plasma current, which, as shown in the bottom panel of Fig. 5 accelerates in the later phases of the VDE, driven by the increased coupling with the tokamak passive structures.
The time evolution of a vertical displacement event in a toy tokamak inspired by the geometry of MAST Upgrade. See Sec. IV B 3 and Fig. 2 for further details and the initial equilibrium. Dashed lines in all panels pertain to the solution of the linearized problem. Full lines refer to the nonlinear solution. The upper set of panels displays snapshots of the separatrix during the evolution. The bottom set of panels shows the evolution of the z coordinate of the plasma magnetic axis, ZO, and of the total plasma current Ip. The black dotted–dashed line shows the slope of the exponential growth with e-folding time τ as predicted by the linear analysis of Sec. IV A 1.
The time evolution of a vertical displacement event in a toy tokamak inspired by the geometry of MAST Upgrade. See Sec. IV B 3 and Fig. 2 for further details and the initial equilibrium. Dashed lines in all panels pertain to the solution of the linearized problem. Full lines refer to the nonlinear solution. The upper set of panels displays snapshots of the separatrix during the evolution. The bottom set of panels shows the evolution of the z coordinate of the plasma magnetic axis, ZO, and of the total plasma current Ip. The black dotted–dashed line shows the slope of the exponential growth with e-folding time τ as predicted by the linear analysis of Sec. IV A 1.
V. DISCUSSION
The axisymmetric evolution of a plasma in a tokamak can be modeled as a sequence of Grad–Shafranov (GS) equilibria parameterized by the currents in the active poloidal field coils and passive vessel structures, and the total poloidal current in the plasma. In turn, these currents are constrained by the circuit equations and Pointing's theorem. This “lumped parameter” approach is justified by the ordering of the relevant characteristic timescales, which allows one to bypass the direct solution of resistive-MHD equations. Despite this simplification, the lumped-parameter model remains inherently nonlinear, with a system of circuit equations coupled to the GS equation itself. Due to the peculiar non-linearity introduced by the unknown plasma boundary, iterative methods are often insufficient. Therefore, modeling the coupled time evolutive problem requires stable and accurate solvers.
In this paper, we introduce a fully Python-based magnetic plasma simulator, FreeGSNKE. FreeGSNKE relies on the public FreeGS code release, which is imported as an external library and from which we mutuate the object framework of tokamak machines, equilibria, profiles, etc. In FreeGSNKE, we implement a Netwon–Krylov (NK) algorithm to solve both non-linear problems at hand: the static GS problem introduced in Sec. II B and the dynamic problem presented in Sec. III A 3. In the absence of analytic gradients, a Jacobian-free NK algorithm greatly reduces the computational expense when compared to a Newton method based on numerically estimated gradients. As part of this paper, we present two reformulations of the nonlinear evolutive problem, both equivalent to the standard set of circuit equations. However, while the first is also phrased as a root problem in the extensive currents, the second is phrased as a root problem in the space of the plasma flux function. Having a range of equivalent problem formulations is a useful asset, as these can then be used to resolve co-linearity and stagnation problems that often affect NK methods. We compare the performance of a set of solution strategies, and present a novel staggered algorithm (strategy II in Sec. IV B 2) leading to improved convergence properties. This strategy makes use of both our equivalent reformulations, by applying the NK solution method to both, separately and consecutively. We show that this approach substantially improves the issues caused by cycles of co-linear vector updates to the Krylov subspace. This novel strategy represents our go-to implementation in FreeGSNKE.
In Sec. II C, we validate the static GS solver against a family of analytic equilibria arising from the Solov'ev special case of the GS equation. For the dynamic problem, in Sec. IV A 1, we focus on a vertical displacement event. We first examine how the number of vessels eigenmodes included in the analysis affects the estimate of the linear growth rate of the vertical instability and then examine how the time evolution resulting from evolving the linearized problem compares to the one obtained by solving the full non-linear problem.
FreeGSNKE is a stateful code and, given that it is written entirely in Python, it is seamlessly integrated with machine-learning libraries and algorithms. This entirely overcomes the need for setting up cross-language binders or socket interfaces for multi-language interaction when using other Python-based resources, as it is often the case, for example, in control design and reinforcement learning studies. As such, FreeGSNKE is well positioned to lower the entry barrier to studies which seek to adopt machine-learning based or other novel approaches to plasma control. FreeGSNKE has already been used to generate large libraries of synthetic plasma equilibria, for the fast emulation of real-time shape control targets,7 and efforts are under way to train reinforcement-learning agents to discover robust, nonlinear shape control strategies.41 Work is also in progress to furnish FreeGSNKE with an accurate machine description for the MAST-U tokamak, including a calibration of the latter based on vacuum shots. This will enable a detailed comparison study between diagnostic data from MAST-U plasma shots and synthetic data from tailored simulations in FreeGSNKE.
ACKNOWLEDGMENTS
M.M. was supported by the UCL Centre for Doctoral Training in Data Intensive Science (STFC Training Grant No. ST/P006736/1). It is a pleasure to thank Ben Dudson for multiple helpful conversations and for providing the community with the open source FreeGS framework. We thank Geof Cunningham for reading an earlier version of this draft and for his expert comments. N.A., A.A., and G.H. also thank Federico Felici for a number of enlightening and motivating discussions and for his clear and invaluable lectures at the “Control & Operations of Tokamaks” intensive course in February 2023, and Dominic Richards for endless encouragement during the development of this work. We also thank Rob Akers and Vassil Alexandrov, the PIs of the STFC Hartree—UKAEA collaboration which funded this work, as well as all of the members on both STFC and UKAEA sides.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Nicola C. Amorisco: Conceptualization (lead); Formal analysis (equal); Investigation (lead); Methodology (lead); Software (lead); Writing—original draft (lead). Adriano Agnello: Conceptualization (supporting); Formal analysis (equal); Investigation (supporting); Methodology (supporting); Software (supporting); Writing—original draft (supporting); Writing—review & editing (supporting). George Holt: Conceptualization (supporting); Methodology (supporting); Software (supporting); Writing—review & editing (supporting). Matthijs Mars: Software (supporting); Writing—review & editing (supporting). James Christopher Buchanan: Conceptualization (supporting); Writing—original draft (supporting); Writing—review & editing (supporting). Stanislas J. P. Pamela: Conceptualization (supporting); Writing—review & editing (supporting).
DATA AVAILABILITY
FreeGSNKE is available upon reasonable request to the corresponding author for collaborative academic use.
APPENDIX: FURTHER IMPLEMENTATION DETAILS
1. Newton–Krylov method
We aim to solve a non-linear problem of the form f(x) = 0 in which f(x) has the same dimensions of the independent variable x. Both static and dynamic plasma equilibrium problems described in this paper can be written in this form. The NK method is a Jacobian-free method. Its stability and convergence properties are analogous to those of the Newton method. However, the NK method does not require estimating the gradient of the non-linear operator f, which may be intractable in some cases, and often exceedingly expensive. For a rigorous presentation of NK methods, we refer the reader to Knoll and Keyes.38