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.

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.

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

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.

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 I m ( t ), featuring all currents in both active and passive tokamak metal structures, I m ( t ) = ( I act ( t ) , I pas ( t ) );

  • the value of the total current in the plasma, I p ( t );

  • a set of parameters describing the plasma profile, θ ( t ), the details of which depend on the chosen parameterization.

The force balance equation relevant to toroidally symmetric, isotropic static plasma equilibria is the GS equation,27,28
(1)
Here:
  • ( R , z , ϕ ) are cylindrical coordinates and μ0 is the magnetic permeability.

  • ψ ( R , z ) is the poloidal flux function, i.e. the poloidal magnetic flux divided by 2 π. 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, ψ = ψ p + ψ m.

  • J p , ϕ ( R , z ) is the plasma toroidal current density distribution, J m , ϕ ( R , z ) is the toroidal current density distribution of metal conductors. In absence of ambiguity, we can drop the subscript ϕ.

Specifically, the plasma current density is given by
(2)
where p is the isotropic plasma pressure and F = R B ϕ , B ϕ being the azimuthal component of the magnetic field. Non-linearities in p ( ψ ) and F F ( ψ ), together with the peculiar non-linearity introduced by the unknown plasma boundary, make Eq. (1) a non-linear elliptic partial differential equation.
For this paper, we focus on profile functions parameterized as in Luxon and Brown29 and Jeon,30 and as adopted in FreeGS,31,
(3)
Here, ψn is the normalized flux function,
where ψa and ψb being the values of the flux function at the magnetic axis and at the plasma boundary, respectively. The radius r R / R 0 is a dimensionless cylindrical radius, where R0 is a scale radius, with no influence on the profile properties. We refer to Luxon and Brown29 and Jeon30 and FreeGS itself for a more detailed description of the physical meaning of the different free parameters ( α m , α n , β 0 , λ ) appearing in Eq. (3). Briefly, FreeGSNKE adheres to FreeGS on how these are used:
  • The user sets values for the pair ( α m , α n ) 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 ( α m , α n ), the product λ β 0 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 λ β 0 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.

In conclusion, there are two possible valid definitions of the plasma-tokamak state S in FreeGSNKE
or
We stress, however, that the computational tools used in FreeGSNKE are agnostic to the specific parameterization used for the profile function. In the following, we assume that one of the two modes above has been chosen, but use of different or more complex parameterizations would indeed be possible with little additional work. Finally, Eq. (3) is only valid inside the plasma domain, while Jp = 0 outside it. Identification of the unknown plasma boundary is part of the free-boundary GS problem itself. The plasma boundary may be defined by either a limiter structure, or, in the case of a diverted plasma, by saddle points in the flux function ψ ( R , z ).

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, ψ p ( R , z ), which satisfies the GS equation (1)

Therefore, the static forward GS problem defines the functional link as follows:
(4)
For clarity, we only seek a solution for the plasma flux ψ p ( R , z ) as given Im, the contribution to the total flux function associated with the tokamak conductors, ψm, is easily found using the relevant free space Green function for toroidally symmetric sources: ψ m ( R , z ) = ψ m ( I m ). We refer the reader to, for example, Jeon.30 It is also worth noticing that our problem statement above does not include an explicit characterization of the toroidal magnetic field Btor. In FreeGS, for example, this is provided by prescribing the value of the toroidal field in the vacuum Fvac = RBtor. Specifying a value for Fvac is necessary to many applications, as for example to calculating the plasma safety factor q, the toroidal and total beta, or the connection length. However, Fvac does not directly enter the GS poloidal force balance or the evolutive problem we focus on. We can therefore safely omit it here.
Jeon30 provides a comprehensive description of how to tackle the forward GS problem above using Picard iterations.32 This involves transforming Eq. (1) to an algebraic equation through spatial discretization of the elliptic differential operator. A solution for ψ p ( R , z ) is then sought iteratively, by solving the so-called “linear GS problem” repeatedly, until convergence. In the kth iteration, an estimate of the source term on the RHS of Eq. (1) is constructed based on the flux function at the previous iteration, ψ p k 1 ( R , z ). The resulting equation for ψ p k ( R , z ) is a linear algebraic equation, which can be solved trivially. One such iterations defines the following functional link:
(5)
between the two subsequent estimates of the plasma flux function ψ p k and ψ p k + 1. Note that, despite this being referred to as solving the linear GS problem, the operation of constructing the source term of the GS equation (1) from the latest estimate of the flux function is intrinsically non-linear, as it requires identifying the plasma boundary.
FreeGS implements the solution method outlined above on a rectangular spatial domain, using a fourth order finite difference representation of the GS elliptic operator, which we also adopt. The iterative method is effective at solving “inverse” (or interpretative) static GS problems.20,30 However, it is not efficient and often unsuccessful when presented with forward problems like those we are interested in, due to a well-known numerical instability characterized by a vertical displacement of the plasma during the numerical iterations.33 Tackling the non-linear dynamic problem requires a stable, efficient forward solver, which is why in FreeGSNKE we replace Picard iterations with an implementation of the Newton–Krylov (NK) method, which we describe in Sec. 1 of  Appendix. An analogous strategy is adopted in Carpanese.21 We use FreeGS functionalities to solve the linear GS problem and formulate the static forward GS problem as a root problem in ψp by using Eq. (5),
(6)
The latter can then be solved using the NK method. We use a relative criterion for convergence, and consider solutions ψp as converged when these satisfy the following condition:
(7)
where the maximum is over the discretization of the spatial domain and we adopt x max ( x ) min ( x ). The dimensionless constant εGS can be adjusted as needed. In FreeGSNKE, it is customary to adopt values for log 10 ε G S between −10 and −5, depending on the requirements of the specific application. In the rest of this paper, we use an upper case F and lower case f to refer, respectively, to absolute and relative residuals of any relevant root problems.
We proceed to validate FreeGSNKE's forward static solver by setting up a comparison with Solov'ev equilibria, following Jeon.30 Solov'ev solutions are fixed boundary solutions of the GS equation (1) defined by a plasma current density (2) of the following form:
(8)
inside the plasma domain, while the current density is zero elsewhere. The total flux function which solves the GS equation (1) for such a current density distribution is
(9)
inside the same plasma domain,34 where ψ0 satisfies
As such, Solov'ev solutions are set by the values of the constants (A1, A2) and by the geometry of the plasma domain. The constants (A1, A2) are in a one-to-one relation with the pair ( I p , β p ): fixing the latter fixes the pair (A1, A2). The plasma domain can be adjusted by selecting an appropriate form for the contribution ψ 0 ( R , z ). We follow Zheng et al.,34 in considering solutions for ψ0 written in the form of an infinite series [see their Eq. (9)], which we also truncate at the fourth term [see their Eq. (14)]. The four coefficients of the non-zero terms in ψ 0 ( R , z ) set the geometry of the plasma domain.
In terms of the profile parameterization (3) used in FreeGSNKE, Eq. (8) is equivalent to setting α n = 0, which removes the dependence of the current distribution on the flux itself. We, therefore, set α n = 0 and use our static NK solver to build an up-down symmetric free-boundary equilibrium for a plasma in a limiter configuration in a simple tokamak geometry, featuring a total of five poloidal field coils. For reference, we use ( I p , β p ) = ( 0.39 M A , 0.27 ). Isoflux contours for our numerical free-boundary solution, ψNK, are displayed as dashed orange lines in the left panel of Fig. 1, together with the location of the poloidal field coils and with the geometry of the limiter wall. We then construct the analytic solution ψS which describes a plasma with the same values for ( I p , β p ), which, as just mentioned, fix the values of (A1, A2) in the RHS of the GS equation. We constrain the separatrix in analytic case to be as close as possible to the one in the numerical case as allowed by the adopted truncation of ψ0. We find the best values of the four independent coefficients which define ψ0 by imposing the set of four constraints
(10)
(11)
FIG. 1.

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.

FIG. 1.

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.

Close modal

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.

The right panel of Fig. 1 further demonstrates and quantifies the accuracy of our numerical solution ψNK by displaying the residue δRHS defined as follows:
(12)
By construction, we expect δ RHS ( R , z ) = 0, since the GS equations solved by ψS and ψNK have identical right hand sides. Therefore, the only contributions to the difference ψ N K ψ S can be originated from a mismatch in ψ0 (for example, due to having adopted a finite number of terms in its series expansion rather than an infinite one). As ψ0 is defined by Δ * ψ 0 = 0, we should find that δ RHS ( R , z ) = 0. Figure 1 shows that this is indeed the case in our numerical solution, where we find | δ RHS ( R , z ) | 3 × 10 6 in the entire domain of validity of the Solov'ev solution. For reference, our solution ψNK is built on a 129 × 129 discretization of the spatial domain and with a target relative accuracy of ε G S = 10 5 for the relevant Eq. (6).

Finally, it is worth mentioning that we have proceeded to verifying that FreeGSNKE's forward static solver agrees with FreeGS' inverse static solver.

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, θ ( t + d t );

  • the set of voltage values V act ( t ) applied to the active poloidal field coils between time t and t + dt.

FreeGSNKE solves for (i) all metal current values I m ( t + d t ) and (ii) the total plasma current I p ( t + d t ). Together with the prescribed set θ ( t + d t ), these complete the tokamak-plasma state S ( t + d t ), thereby providing a full characterization of the equilibrium at time t + dt. The following sections detail the conditions that the state S ( t + d t ) have to satisfy and describe the coupling between plasma and metal structures.

1. The metals

We start with the currents in the tokamak's metal structures, including both poloidal field coils and the system of toroidal coils used to model the tokamak's passive elements. The equation which describes the currents' evolution is the combination of Ohm's law and Faraday's law,
(13)
Here, we use subscripts as labels for the spaces the matrices act on and map to, as opposed to actual matrix indices. Specifically,
  • M m , m is the matrix of mutual inductances, describing the coupling between each pair of metal circuits;

  • R m , m 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 V m = ( V act , 0 ), where Vact is the set of active voltage values applied to the active poloidal field coils;

  • d Ψ p / d t 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 Ψ p should not be confused with the plasma flux ψp.

Equation (13) is a system of N m = N act + N pas circuit equations, where Nact is the number of active poloidal field coils, and Npas is the number of passive coils used to model the tokamak's passive structures. We assume that the matrices M m , m and R m , m are both assigned, for example, as defined by the machine geometry and by the properties of the tokamak's power supply. Introducing the matrix of mutual inductances between each of the tokamak conductors and the grid points of the discretized spatial domain, M m , y, we can rewrite Eq. (13) as
(14)
where Iy is the plasma current density distribution discretized over the spatial domain.

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 I d , pas 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 I d , pas are the eigenvectors of the matrix R pas 1 / 2 M pas , pas 1 R pas 1 / 2, where Rpas and M pas , pas are the resistance and mutual inductance matrices relevant to the passive structures. The characteristic frequencies ω d , pas are the associated eigenvalues. Normal modes with fast dynamics are easily dropped on the basis of the frequencies ω d , pas. To identify modes with a weak coupling with the plasma, we examine the response matrix I y / I d, 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.

Using our normal mode decomposition, the circuit equations (14) can be written as
(15)
where the matrices Λ 1 and P are defined in Sec. 2 of  Appendix, as is the new set of unknown currents, Id.

2. The plasma

Different approaches have been used in the literature to model the plasma itself when coupled with the metals.13,17,25,35–37 Most are united by the assumption of a massless plasma, followed by the use of resistive MHD equations, aimed at deriving Poynting's Theorem for the balance of magnetic energy and Ohmic dissipation—while the work done by the plasma is often neglected. Here, we follow Albanese17 and Carpanese,21 and use the following evolutive equation:
(16)
in which
  • M y , y is the matrix of mutual inductances between pairs of grid points in the discretized spatial domain on which Iy is defined;

  • R y , y is the diagonal matrix of plasma resistance values;

  • M y , m = M m , y T, with M m , y defined in Sec. III A 1.

As for M m , m and R m , m, we assume that M y , y , M m , y and R y , y are all known matrices. In particular, we model the plasma resistance as
(17)
where σp is the assigned plasma resistivity, R and dA are the cylindrical radius and the area of each element of the Discretized spatial domain. Before the contraction at the left hand side, Eq. (16) describes the system of circuit equations defined by each infinitesimal toroidal plasma element. This system is then contracted with the normalized plasma current distribution I y T / I p, making it a one-dimensional constraint, sometimes referred to as a lumped parameter model. In the following, we use
For completeness, we rewrite Eq. (16) in terms of the vector of normal mode currents Id,
(18)
The latter has also been divided by the lumped plasma resistance
(19)
so to make it dimensionally coherent with Eq. (15).

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, R y , y , M y , y and M m , y. 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.

Let us focus on a specific time step, and assume that S(t) and θ ( t + d t ) are given, as specified in Sec. III. It is useful here to introduce a symbol for the full set of extensive currents, collecting together all metal currents and the total plasma current,
The plasma current distribution Iy is a function of the extensive currents and of the plasma flux ψp: through Eqs. (2) and (3), we have that
(20)
Therefore, the pair of the extensive currents I e ( t + d t ) and the plasma flux ψ p ( t + d t ) represent a complete set of unknowns of the evolutive problem. A pair [ I e ( t + d t ) , ψ p ( t + d t ) ] is a solution if it satisfies both
  • the root problem defined by the dynamics equations
    (21)

    where we have brought Eqs. (15) and (18) together;

  • the root problem associated to the GS equation (6)
    (22)
For completeness, we also report the relative convergence criteria our solutions have to satisfy in FreeGSNKE. These criteria are formulated in terms of the evolution between time t and t + dt of both currents Ie and plasma flux ψp,
(23)
(24)
Solutions are considered converged if
(25)
(26)
where both divisions are intended as element-wise operations and both maxima are over the dimensions of the relevant vectors. As for Eq. (7), we use x max ( x ) min ( x ). Both εI and ε G S dyn are adjustable dimensionless constants. Depending on the requirements of the specific applications, values between 10 6 and 10 2 may be used. Note that the way in which the residual (26) is made relative is different from the one used for the static GS solver, Eq. (7), thereby the additional superscript dyn. In particular, since δ ψ p is usually considerably smaller that ψ p , caution should be taken when comparing values of εGS and ε G S dyn. More, in general, the adoption of relative convergence criteria in FreeGSNKE ensures a coherent level of accuracy is achieved at each time step. Note, however, that the criteria (25) may not be directly comparable with the absolute thresholds used sometimes in other codes.

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 [ I e ( t + d t ) , ψ p ( t + d t ) ] or

  • by restricting the search to the manifold of GS solutions, i.e. only considering trial solutions of the following form:
    (27)

    where ψ p G S ( I e ) 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 I e ( t + d t ).

FreeGSNKE adopts solution methods belonging to both of these two classes, which we detail in Sec. IV.

Before tackling the full non-linear problem, it is useful to consider its linearized version. A linearization of the evolutive problem can be obtained within the framework of the restricted trial solutions defined in Sec. III A 3, of the form [ I e , ψ p G S ( I e ) ]. As the profile parameters θ are given, under these assumptions the plasma current distribution is a function of the extensive currents alone,
(28)
By expanding the latter around a given plasma-tokamak state S0, it is possible to rewrite the terms I ̇ y appearing in Eqs. (15) and (18) as
(29)
We stress that the superscript GS highlights that all derivatives in Eq. (29), sometimes called response matrices, are calculated on the manifold of GS solutions, i.e., assume use of Eqs. (27) and (28). In FreeGSNKE, we approximate the response matrices using finite differences, by solving the associated static GS problems around a user-defined equilibrium S0.

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 ψ p ( t ) 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.

FIG. 2.

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 ω d max. Orange and cyan symbols show the number of modes with ω d < ω d max 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 ω d max. Colored symbols refer to only including modes with non-negligible coupling with the plasma (see Sec. IV A 1 for more details).

FIG. 2.

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 ω d max. Orange and cyan symbols show the number of modes with ω d < ω d max 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 ω d max. Colored symbols refer to only including modes with non-negligible coupling with the plasma (see Sec. IV A 1 for more details).

Close modal

We build the matrix of mutual inductances relevant to metal elements, M m , m, 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 R m , m, 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 M m , m and R m , m 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 ω d max, 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 ω d max may be selected and used to drop modes with ω d > ω d max.

Orange lines in the left-most panel of Fig. 2 display isoflux contours for an up-down symmetric diverted plasma with
(30)
We linearize the dynamics around this equilibrium as described in Sec. IV A, which involves numerically approximating the response matrices (29) through finite differences. We can then estimate the e-folding time τ of the vertically unstable plasma by calculating the eigenvalues of the state matrix of the linearized circuit equations. Note that this linear growth rate does not assume a rigid plasma. Black plus signs in the right-most panel of Fig. 2 show such e-folding time τ as a function of the truncation frequency ω d max, assuming modes with ω d > ω d max have been dropped. As expected, a less accurate model of the vessel, i.e., with a smaller number of retained modes Nd, corresponds to a diminished coupling with the plasma, resulting in a shortened e-folding time τ and, therefore, an artificially faster vertical instability.

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 I y / I d. 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 I y / I d. For reference, here, we use | I y / I d | > 0.1. 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 ω d max. 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 ω d max = 10 2 KHz, the VDE growth rate remains consistent at better than 0.005. 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 | I y / I d | > 0.1.

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

First, we note that, if we take the normalized plasma current distribution I ̂ y ( t + d t ) to be an assigned quantity, the circuit equations (21) become a readily solvable linear system for the extensive currents I e ( t + d t ). Such a linear system is reported in Sec. 3 of  Appendix, where we have used an implicit Euler method. We refer to the set of extensive currents obtained by solving this system with
(31)
where the superscript c.eq. highlights that this mapping is provided by the circuit equations themselves. We now fold in that the normalized plasma current distribution I ̂ y is itself a function of Ie and ψp through Eq. (20), I ̂ y = I ̂ y ( I e , ψ p ). Therefore, for a given trial solution of the evolutive problem [ I e ( t + d t ) , ψ p ( t + d t ) ], we can calculate the set of extensive currents
(32)
by using Eq. (20) first, and then by solving the resulting linear system with Eq. (31) (see also Sec. 3 of  Appendix). Somewhat similarly to the formulation of the linear GS problem (5), the function (32) essentially “iterates” the trial set of extensive currents I e ( t + d t ) through the circuit equations (21) while assuming a trial plasma flux ψ p ( t + d t ). Pairs [ I e ( t + d t ) , ψ p ( t + d t ) ] for which I e c . e q . = I e also result in F c . e q . = 0. In other words, the root problem,
(33)
is equivalent to the standard formulation of the evolutive problem [Eq. (21)]. To differentiate our alternative formulation, we use the subscript it.c.eq, referencing its iterative nature. It is worth noticing that, computing the residuals F i t . c . e q . for a given trial solution [ I e ( t + d t ) , ψ p ( t + d t ) ] does not require, at any stage, solving a GS problem, as Eq. (20) does not do so.
Continuing from the above, we construct a second equivalent root problem, having, however, a residual in the space of the plasma function ψp rather than a residual in the space of the extensive currents Ie, as both Eqs. (21) and (33). By solving the GS problem for the set of extensive currents obtained through Eq. (32), we get
(34)
Like the extensive currents I c . e q ., we can interpret the plasma flux ψ p c . e q . as having been “iterated” thought the circuit equations—for some value of the trial extensive currents I e ( t + d t ). It is then natural to define the root problem
(35)
which, as Eq. (33), is also equivalent to the standard circuit equations (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 f i t . c . e q . corresponding to Eq. (33), and the relative residual f i t . ψ 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 [ I e ( t + d t ) , ψ p ( t + d t ) ] 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 [ I e ( t + d t ) , ψ p ( t + d t ) ] 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.

FIG. 3.

A demonstration of the equivalence between the standard circuit equations (21) leading to the relative residual f c . e q ., and our reformulations (33) and (35), leading to the relative residuals f i t . c . e q . and f i t . ψ, respectively. The figure displays the evolution of such residuals toward convergence, as a function of the iteration step.

FIG. 3.

A demonstration of the equivalence between the standard circuit equations (21) leading to the relative residual f c . e q ., and our reformulations (33) and (35), leading to the relative residuals f i t . c . e q . and f i t . ψ, respectively. The figure displays the evolution of such residuals toward convergence, as a function of the iteration step.

Close modal

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 [ I e , ψ p G S ( I e ) ] parameterized by the extensive currents alone and located on the manifold of the GS solutions in the space of all trial solutions [ I e ( t + d t ) , ψ p ( t + d t ) ]. 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 I e ( t + d t ) 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 1 / 10 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.

FIG. 4.

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.

FIG. 4.

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.

Close modal

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 I e ( t + d t ) and to the plasma flux ψ p ( t + d t ), 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 ( F c . e q . , F G S ) = 0, thereby building updates for the entire vector on unknowns [ I e ( t + d t ) , ψ p ( t + d t ) ] at once. Intuitively, it can be expected that alternating NK solution steps in I e ( t + d t ) and ψ p ( t + d t ) 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.

As SII is an unrestricted strategy, the root problem (22) associated with the GS equation is not satisfied by construction and needs enforcing. However, application of the NK method to Eqs. (33) and (35) does not lead to an improvement in the GS residual (22) which means that this strategy also requires a dedicated algorithmic step aimed at this. We use a convex update of the plasma flux: given the trial solution ( I e , ψ p ), the following blended combination:
(36)
corresponds to a reduced GS residual. Here, ψ p G S ( I e ) is the GS solution corresponding to set of extensive currents Ie and b is the dimensionless coefficient which sets the degree of blending in the convex combination. For completeness, our implementation uses b = 0.5, which we find to be advantageous in terms of performance.

Finally, we define SII as the following algorithm:

  1. Use the solution of the linearized problem I e ( t + d t ) (as defined in Sec. IV A) and the associated GS plasma flux, ψ p G S ( I e ( t + d t ) ), as initial trial solution: [ I e ( t + d t ) , ψ p G S ( I e ( t + d t ) ].

  2. If the current trial solution [ I e ( t + d t ) , ψ p ( t + d t ) ] fails the desired relative convergence criterion for the GS residual equation (26) update ψ p ( t + d t ) using Eq. (36).

  3. If [ I e ( t + d t ) , ψ p ( t + d t ) ] is such that the relative convergence criterion set for eq. (35) fails, update ψ p ( t + d t ) by applying the NK method to the same root problem.

  4. If [ I e ( t + d t ) , ψ p ( t + d t ) ] is such that the convergence criterion set for Eq. (33) fails, update I e ( t + d t ) by applying the NK method to the same root problem.

  5. Check the current trial solution [ I e ( t + d t ) , ψ p ( t + d t ) ] 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 f i t . c . e q . 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 ε I = ε G S dyn = 10 4. 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 f i t . c . e q . 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 σ p = 10 6 Ω m. For reference, this results in a lumped plasma resistance of R p = I ̂ y T R y , y I ̂ y = 4.2 × 10 6 Ω for our chosen equilibrium.

  • In absence of a self-consistent evolution of the profile parameters ( α m , α n ) 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: V act = R act , act I act ( t = 0 ), where R act , act is the resistance matrix of the poloidal field coils and I act ( t = 0 ) 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 t 10 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.

FIG. 5.

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.

FIG. 5.

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.

Close modal

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.

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.

The authors have no conflicts to disclose.

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

FreeGSNKE is available upon reasonable request to the corresponding author for collaborative academic use.

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 

Like the Newton method, the NK method is also iterative: given a guess for the solution xk and the corresponding residual f ( x k ), a suitable δx is sought so that
(A1)
In the Newton method, δx is found based on the knowledge of the full Jacobian f / x at x = xk, so to minimize the residual
(A2)
The NK method, instead, only uses projections of the Jacobian itself, vl,
(A3)
where the vectors ul are elements of a suitable basis set U. With the above-mentioned projections, the least squares problem (A2) becomes
(A4)
where ϖ l are the linear coefficients. Correspondingly, the update δx is obtained by combining the basis vectors ul so that
(A5)
This technique is generally referred to as generalized minimal residual method.42 Use of the Krylov subspace
(A6)
as a set of basis vectors gives the NK method the second part of its name. It is often useful to orthogonalize the basis U while this is constructed, which is a procedure we adopt in FreeGSNKE. Also, care is put in FreeGSNKE to identifying and avoiding cases of collinearity between the Jacobian projections vl, which may interfere with the solution of the least squares problem (A4).
The NK method is, therefore, composed of two nested loops. The outer loop updates the trial solution as in Eq. (A1) and according to Eqs. (A5) and (A4). The inner loop is used to sequentially build new basis vectors ul and projections vl and is set to truncate the Krylov subspace (A6) at an appropriate value of L as needed. As such, the NK method is advantageous when a closed analytical form of the Jacobian is not available, and numerical approximation would have to be used to estimate it, for example using finite differences. The number of projections L used at a given iteration of the outer loop is often substantially smaller than the number of dimensions of the independent variable x. Clearly, the projections defined in Eq. (A3) can be also be approximated using finite differences,
(A7)
A range of criteria can be used to identify a suitable value L where to interrupt the inner loop. Often, a compromise is to be sought, as using fewer projections is bound to result in the need for an increased number of iterations of the outer NK loop. In FreeGSNKE we determine L using a relative criterion on Eq. (A4): we interrupt the inner loop when the available projections vl satisfy
(A8)
We find that different values of the tolerance εLS are best suited when applying the NK method to different problems within FreeGSNKE. Similarly, different absolute and relative tolerance criteria are used to interrupt the outer NK loop on a converged solution.
2. Normal mode decomposition
We start by identifying the basis of normal modes which diagonalizes the dynamics of the isolated passive structures
(A9)
where Rpas and M pas , pas are the resistance and mutual inductance matrices relevant to the passive structures alone. The matrix Λpas is diagonal and collects the characteristic frequencies ω d , pas of each of the vessel normal modes. The matrix Ppas is an orthogonal matrix, P pas 1 = P pas T and each of its columns identifies one of such modes.
Using the matrix Ppas we can define a global change of basis matrix P, applicable to the entire set of metal currents Im. The orthogonal matrix P leaves the active coil currents unchanged, while it is equal to Ppas over the passive structures
(A10)
where we have used the symbol 11 ( n act ) for the identity matrix with dimension nact, the number of independent active coils. With the latter definition, we can introduce the vector of normal modes Id, which we use to formulate the evolutive problem,
(A11)
Individual normal modes can be dropped by simply dropping the corresponding columns of the matrix P (rows of the matrix PT). For simplicity, and despite the slight abuse of notation, in the following, we keep the symbols P (and P 1), although it is intended that some of their columns (rows) may have been dropped with respect to the definition of Eq. (A10).
With the above, we can compose the matrix Λ 1 which first appears in the main text in the circuit equations (15),
(A12)
Note that, differently from the case of Λpas, Λ, and Λ 1 are not diagonal matrices.
3. Linear system for the iterated circuit equations
We assume that the normalized plasma distribution I ̂ y ( t + d t ) is a given vector, so that I y ( t + d t ) = I p ( t + d t ) I ̂ y ( t + d t ), where only I p ( t + d t ) is unknown. With this, Eq. (21) becomes a linear system for the extensive currents I e ( t + d t ) = ( I d ( t + d t ) , I p ( t + d t ) ). Here, we have discretized such system using an implicit Euler setting
(A13)
where
(A14)
(A15)
with
(A16)
With these definitions, Eq. (A13) can be readily solved for I e ( t + d t ). In the main text, we refer to this solution with Eqs. (31) and (32).
1.
J. W.
Berkery
,
S. A.
Sabbagh
,
L.
Kogan
,
D.
Ryan
,
J. M.
Bialek
,
Y.
Jiang
,
D. J.
Battaglia
,
S.
Gibson
, and
C.
Ham
, “
Kinetic equilibrium reconstructions of plasmas in the mast database and preparation for reconstruction of the first plasmas in mast upgrade
,”
Plasma Phys. Controlled Fusion
63
,
055014
(
2021
).
2.
L.
Appel
and
I.
Lupelli
, and
JET Contributors
, “
Equilibrium reconstruction in an iron core tokamak using a deterministic magnetisation model
,”
Comput. Phys. Commun.
223
,
1
17
(
2018
).
3.
A.
Hudoba
,
S.
Newton
,
G.
Voss
,
G.
Cunningham
, and
S.
Henderson
, “
Divertor optimisation and power handling in spherical tokamak reactors
,”
Nucl. Mater. Energy
35
,
101410
(
2023
).
4.
H.
Anand
,
O.
Bardsley
,
D.
Humphreys
,
M.
Lennholm
,
A.
Welander
,
Z.
Xing
,
J.
Barr
,
M.
Walker
,
J.
Mitchell
, and
H.
Meyer
, “
Modelling, design and simulation of plasma magnetic control for the spherical tokamak for energy production (step)
,”
Fusion Eng. Des.
194
,
113724
(
2023
).
5.
S.
Joung
,
J.
Kim
,
S.
Kwak
,
J. G.
Bak
,
S. G.
Lee
,
H. S.
Han
,
H. S.
Kim
,
G.
Lee
,
D.
Kwon
, and
Y. C.
Ghim
, “
Deep neural network Grad–Shafranov solver constrained with measured magnetic signals
,”
Nucl. Fusion
60
,
016034
(
2020
).
6.
L. L.
Lao
,
S.
Kruger
,
C.
Akcay
,
P.
Balaprakash
,
T. A.
Bechtel
,
E.
Howell
,
J.
Koo
,
J.
Leddy
,
M.
Leinhauser
,
Y. Q.
Liu
,
S.
Madireddy
,
J.
McClenaghan
,
D.
Orozco
,
A.
Pankin
,
D.
Schissel
,
S.
Smith
,
X.
Sun
, and
S.
Williams
, “
Application of machine learning and artificial intelligence to extend EFIT equilibrium reconstruction
,”
Plasma Phys. Controlled Fusion
64
,
074001
(
2022
).
7.
A.
Agnello
,
N.
Amorisco
,
A.
Keats
,
G.
Holt
,
J.
Buchanan
, and
S.
Pamela
, “
Emulation for scenario design and classical control of tokamak plasmas
,”
Phys. Plasmas
31
,
043901
(
2024
).
8.
C. M.
Samuell
,
A. G.
Mclean
,
C. A.
Johnson
,
F.
Glass
, and
A. E.
Jaervinen
, “
Measuring the electron temperature and identifying plasma detachment using machine learning and spectroscopy
,”
Rev. Sci. Instrum.
92
,
043520
(
2021
).
9.
B.
Zhu
,
M.
Zhao
,
H.
Bhatia
,
X.-q.
Xu
,
P.-T.
Bremer
,
W.
Meyer
,
N.
Li
, and
T.
Rognlien
, “
Data-driven model for divertor plasma detachment prediction
,”
J. Plasma Phys.
88
,
895880504
(
2022
).
10.
G. a.
Holt
, “
Tokamak divertor plasma emulation with machine learning
,” unpublished (
2023
).
11.
J.
Kates-Harbeck
,
A.
Svyatkovskiy
, and
W.
Tang
, “
Predicting disruptive instabilities in controlled fusion plasmas through deep learning
,”
Nature
568
,
526
531
(
2019
).
12.
J.
Vega
,
A.
Murari
,
S.
Dormido-Canto
,
G.
Rattá
, and
M.
Gelfusa
, “
Disruption prediction with artificial intelligence techniques in tokamak plasmas
,”
Nat. Phys.
18
,
741
750
(
2022
).
13.
J.
Degrave
,
F.
Felici
,
J.
Buchli
,
M.
Neunert
,
B.
Tracey
,
F.
Carpanese
,
T.
Ewalds
,
R.
Hafner
,
A.
Abdolmaleki
,
D.
de las Casas
,
C.
Donner
,
L.
Fritz
,
C.
Galperti
,
A.
Huber
,
J.
Keeling
,
M.
Tsimpoukelli
,
J.
Kay
,
A.
Merle
,
J.-M.
Moret
,
S.
Noury
,
F.
Pesamosca
,
D.
Pfau
,
O.
Sauter
,
C.
Sommariva
,
S.
Coda
,
B.
Duval
,
A.
Fasoli
,
P.
Kohli
,
K.
Kavukcuoglu
,
D.
Hassabis
, and
M.
Riedmiller
, “
Magnetic control of tokamak plasmas through deep reinforcement learning
,”
Nature
602
,
414
419
(
2022
).
14.
A.
Ho
,
J.
Citrin
,
C.
Bourdelle
,
Y.
Camenen
,
F. J.
Casson
,
K. L.
van de Plassche
, and
H.
Weisen
, and
JET Contributors
, “
Neural network surrogate of QuaLiKiz using JET experimental data to populate training space
,”
Phys. Plasmas
28
,
032305
(
2021
).
15.
R.
Khayrutdinov
and
V.
Lukash
, “
Studies of plasma equilibrium and transport in a tokamak fusion device with the inverse-variable technique
,”
J. Comput. Phys.
109
,
193
201
(
1993
).
16.
D.
Humphreys
,
J.
Ferron
,
M.
Bakhtiari
,
J.
Blair
,
Y.
In
,
G.
Jackson
,
H.
Jhang
,
R.
Johnson
,
J.
Kim
,
R.
LaHaye
,
J.
Leuer
,
B.
Penaflor
,
E.
Schuster
,
M.
Walker
,
H.
Wang
,
A.
Welander
, and
D.
Whyte
, “
Development of ITER-relevant plasma control solutions at DIII-D
,”
Nucl. Fusion
47
,
943
(
2007
).
17.
R.
Albanese
, “
CREATE-NL+: A robust control-oriented free boundary dynamic plasma equilibrium solver
,”
Fusion Eng. Des.
96
,
664
667
(
2015
).
18.
A.
Ivanov
,
R.
Khayrutdinov
,
S.
Medvedev
, and
Y.
Poshekhonov
, “
New adaptive grid plasma evolution code SPIDER
,” in
32nd EPS Conference on Plasma Physics, Tarragona, Spain
(
International Atomic Energy Agency
,
2005
), Vol.
29C
, p.
P–5.063
.
19.
B.
Faugeras
, “
An overview of the numerical methods for tokamak plasma equilibrium computation implemented in the nice code
,”
Fusion Eng. Des.
160
,
112020
(
2020
).
20.
H.
Heumann
,
J.
Blum
,
C.
Boulbe
,
B.
Faugeras
,
G.
Selig
,
J.-M.
Ané
,
S.
Brémond
,
V.
Grandgirard
,
P.
Hertout
,
E.
Nardon
et al, “
Quasi-static free-boundary equilibrium of toroidal plasma with CEDRES: Computational methods and applications
,”
J. Plasma Phys.
81
,
905810301
(
2015
).
21.
F.
Carpanese
, “
Development of free-boundary equilibrium and transport solvers for simulation and real-time interpretation of tokamak experiments
,” Technical Report,
EPFL
,
2021
.
22.
L.
Degtyarev
and
V.
Drozdov
, “
An inverse variable technique in the MHD-equilibrium problem
,”
Comput. Phys. Rep.
2
,
341
387
(
1985
).
23.
S.
Kerboua-Benlarbi
,
R.
Nouailletas
,
B.
Faugeras
, and
P.
Moreau
, “
Deep reinforcement learning for magnetic control on WEST
,” in
Workshop on AI for Accelerating Fusion and Plasma Science
(
International Atomic Energy Agency
,
2023
).
24.
See https://github.com/freegs-plasma/freegs for the documentation of the latest public release.
25.
J. P.
Freidberg
,
Plasma Physics and Fusion Energy
(
Cambridge University Press
,
2008
).
26.
F.
Felici
,
O.
Sauter
,
T.
Goodman
, and
J.
Paley
, “
RAPTOR: Optimization, real-time simulation and control of the tokamak q profile evolution using a simplified transport model
,” in
APS Division of Plasma Physics Meeting Abstracts, APS Meeting Abstracts
(
American Physical Society
,
2010
), Vol.
52
, p.
BP9.090
.
27.
V.
Shafranov
, “
On equilibrium magnetohydrodynamic configurations
,”
Zh. Eksp. Teor. Fiz.
33
,
710
722
(
1957
).
28.
V. D.
Shafranov
, “
On magnetohydrodynamical equilibrium configurations
,”
Sov. J. Exp. Theor. Phys.
6
,
545
(
1958
).
29.
J. L.
Luxon
and
B. B.
Brown
, “
Magnetic analysis of non-circular cross-section tokamaks
,”
Nucl. Fusion
22
,
813
821
(
1982
).
30.
Y. M.
Jeon
, “
Development of a free-boundary tokamak equilibrium solver for advanced study of tokamak equilibria
,”
J. Korean Phys. Soc.
67
,
843
853
(
2015
).
31.
Although we note that FreeGS also allows for more generic functions p ( ψ ) and F 2 ( ψ ), for example as specified by spline functions.
32.
H.
Kreyszig
and
E.
Kreyszig
,
Advanced Engineering Mathematics
, Student Solutions Manual and Study Guide Vol.
1
(
John Wiley & Sons
,
2012
), Chap. 1–12.
33.
J. L.
Johnson
,
H.
Dalhed
,
J.
Greene
,
R.
Grimm
,
Y.
Hsieh
,
S.
Jardin
,
J.
Manickam
,
M.
Okabayashi
,
R.
Storer
,
A.
Todd
et al, “
Numerical determination of axisymmetric toroidal magnetohydrodynamic equilibria
,”
J. Comput. Phys.
32
,
212
234
(
1979
).
34.
S. B.
Zheng
,
A. J.
Wootton
, and
E. R.
Solano
, “
Analytical tokamak equilibrium for shaped plasmas
,”
Phys. Plasmas
3
,
1176
1178
(
1996
).
35.
J. B.
Lister
,
A.
Sharma
,
D. J. N.
Limebeer
,
Y.
Nakamura
,
J. P.
Wainwright
, and
R.
Yoshino
, “
Plasma equilibrium response modelling and validation on JT-60U
,”
Nucl. Fusion
42
,
708
724
(
2002
).
36.
M.
Ariola
,
G. D.
Tommasi
,
D.
Mazon
,
D.
Moreau
,
F.
Piccolo
,
A.
Pironti
,
F.
Sartori
, and
L.
Zabeo
, “
Integrated plasma shape and boundary flux control on JET tokamak
,”
Fusion Sci. Technol.
53
,
789
805
(
2008
).
37.
F.
Carpanese
,
F.
Felici
,
C.
Galperti
,
A.
Merle
,
J. M.
Moret
, and
O.
Sauter
, and
the TCV Team
. “
First demonstration of real-time kinetic equilibrium reconstruction on TCV by coupling LIUQE and RAPTOR
,”
Nucl. Fusion
60
,
066020
(
2020
).
38.
D. A.
Knoll
and
D. E.
Keyes
, “
Jacobian-free Newton–Krylov methods: A survey of approaches and applications
,”
J. Comput. Phys.
193
,
357
397
(
2004
).
39.
M. A.
Gomes-Ruggiero
,
V. L. R.
Lopes
, and
J. V.
Toledo-Benavides
, “
A safeguard approach to detect stagnation of GMRES(m) with applications in Newton–Krylov methods
,”
Comput. Appl. Math.
27
,
175
199
(
2008
).
40.
The associated response matrices need not be recalculated at each time step.
41.
N.
Amorisco
,
A.
Agnello
,
G.
Holt
,
M.
Mars
,
A.
Keats
,
S.
Pamela
,
J.
Buchanan
,
G.
McArdle
,
A.
Gray
,
D.
Richards
,
V.
Alexandrov
, and
R.
Akers
, “
Machine learning for plasma shape control on MAST-U
,” in
IAEA FEC23
(
2023
).
42.
Y.
Saad
and
M. H.
Schultz
, “
GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems
,”
SIAM J. Sci. Statist. Comput.
7
,
856
869
(
1986
).