Example: Static inverse free-boundary equilibrium calculations

---

This example notebook shows how to use FreeGSNKE to solve static inverse free-boundary Grad-Shafranov (GS) problems.

In the inverse solve mode we seek to estimate the active poloidal field coil currents using user-defined constraints (e.g. on isoflux, x-point, o-point, and psi values) and plasma current density profiles for a desired equilibrium shape.

Note that during this solve, currents are not found in any specified passive structures.

Below, we illustrate how to use the solver for both diverted and limited plasma configurations in a MAST-U-like tokamak using stored pickle files containing the machine description. These machine description files partially come from the FreeGS repository and are not an exact replica of MAST-U.

The static free-boundary Grad-Shafranov problem

Here we will outline the static free-boundary GS problem that is solved within both the forward and inverse solvers, though we encourage you to see Pentland et al. (2024) for more details.

Using a cylindrical coordinate system \((R,\phi,Z)\), the aim is to solve the GS equation:

\[ \Delta^* \psi \equiv \left( \frac{\partial^2}{\partial R^2} - \frac{1}{R} \frac{\partial}{\partial R} + \frac{\partial^2}{\partial Z^2} \right) \psi = -\mu_0 R J_{\phi}(\psi, R, Z), \qquad (R,Z) \in \Omega, \]

for the poloidal flux \(\psi(R,Z)\) (which here has units Weber/\(2\pi\)) in the rectangular computational domain \(\Omega\). The flux has contributions from both the plasma and the coils (metals) such that \(\psi = \psi_p + \psi_c\). This flux defines the toroidal current density \(J_{\phi} = J_p(\psi,R,Z) + J_c(R,Z)\), also containing a contribution from both the plasma and coils, respectively. We have the plasma current density (only valid in the core plasma region \(\Omega_p\)):

\[ J_p(\psi,R,Z) = R \frac{\mathrm{d}p}{\mathrm{d}\psi} +\frac{1}{\mu_0 R} F \frac{\mathrm{d} F}{\mathrm{d} \psi}, \qquad (R,Z) \in \Omega_p, \]

where \(p(\psi)\) is the plasma pressure profile and \(F(\psi)\) is the toroidal magnetic field profile. The current density generated by \(N\) active coils and passive structures is given by:

\[\begin{split} J_c(R,Z) = \sum_{j=1}^{N} \frac{I^c_j(R,Z)}{A_j^c}, \qquad (R,Z) \in \Omega, \quad \text{where} \quad I_j^c(R,Z) = \begin{cases} I_j^c & \text{if } (R,Z) \in \Omega_j^c, \\ 0 & \text{elsewhere}. \end{cases}.\end{split}\]

This makes use of the current \(I^c_j\) in each metal and its cross-sectional area \(A^c_j\) (the domain of each metal is denoted by \(\Omega_j^c\)).

To complete the problem, we need the integral (Dirichlet) free-boundary condition:

\[ \psi(R,Z) = \int_{\Omega} G(R,Z;R',Z') J_{\phi}(\psi, R',Z') \ \mathrm{d}R' \mathrm{d}Z', \qquad (R,Z) \in \partial \Omega, \]

where \(G\) is the (known) Green’s function for the elliptic operator above.

We’ll now go through the steps required to solve the inverse problem in FreeGSNKE.

Create the machine object

First, we build the machine object from previously created pickle files in the “machine_configs/MAST-U” directory.

FreeGSNKE requires the following paths in order to build the machine:

• active_coils_path

• passive_coils_path

• limiter_path

• wall_path

• magnetic_probe_path (not required here)

# build machine
from freegsnke import build_machine
tokamak = build_machine.tokamak(
    active_coils_path=f"../machine_configs/MAST-U/MAST-U_like_active_coils.pickle",
    passive_coils_path=f"../machine_configs/MAST-U/MAST-U_like_passive_coils.pickle",
    limiter_path=f"../machine_configs/MAST-U/MAST-U_like_limiter.pickle",
    wall_path=f"../machine_configs/MAST-U/MAST-U_like_wall.pickle",
)

Active coils --> built from pickle file.
Passive structures --> built from pickle file.
Limiter --> built from pickle file.
Wall --> built from pickle file.

Magnetic probes --> none provided.

Resistance (R) and inductance (M) matrices --> built using actives (and passives if present).
Tokamak built.

# plot the machine
import matplotlib.pyplot as plt

fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)
plt.tight_layout()

tokamak.plot(axis=ax1, show=False)
ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")

ax1.grid(alpha=0.5)
ax1.set_aspect('equal')
ax1.set_xlim(0.1, 2.15)
ax1.set_ylim(-2.25, 2.25)
ax1.set_xlabel(r'Major radius, $R$ [m]')
ax1.set_ylabel(r'Height, $Z$ [m]')

Text(10.027777777777777, 0.5, 'Height, $Z$ [m]')

Instantiate an equilibrium

We are now ready to build a plasma equilibrium object for our tokamak. This is done using the freegs4e.Equilibrium class, which implicitly defines the rectangular domain of the solver as well as the grid resolution.

Equilibrium has sensible defaults, but it is recommended to define the radial and vertical domain of the grid using the Rmin, Rmax, Zmin and Zmax parameters, as well as the grid resolution in the radial and vertical directions with the nx and ny parameters. The grid will be initialised using fourth-order finite differences. Note that the computational grid must encompass the domain enclosed by the limiter object (as this is where the plasma will be confined to).

A tokamak object should be supplied to the tokamak parameter to assign the desired machine to the equilibrium.

If available, an initial guess for the plasma flux \(\psi_p\) (dimensions nx x ny) can be provided via the psi parameter (commented out in the following code). If not, the default initialisation will be used.

The eq object will store a lot of important information and derived quantites once the equilibrium has been calculated (see future notebook on this).

from freegsnke import equilibrium_update

eq = equilibrium_update.Equilibrium(
    tokamak=tokamak,      # provide tokamak object
    Rmin=0.1, Rmax=2.0,   # radial range
    Zmin=-2.2, Zmax=2.2,  # vertical range
    nx=65,                # number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)
    ny=129,               # number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)
    # psi=plasma_psi
)

Instantiate a profile object

We can now instatiate a profile object that contains the chosen parameterisation of the toroidal plasma current density \(J_p\) (i.e. on right hand side of the GS equation). We can then set the paramters for the chosen current density profiles.

A number of commonly used profile parameterisations exist in FreeGSNKE, including:

• ConstrainPaxisIp

• ConstrainBetapIp

• Fiesta_Topeol

• Lao85

• TensionSpline

In this notebook, we will make use of the ConstrainPaxisIp (and ConstrainBetapIp) profiles (see Jeon (2015)). Others will be utilised in later notebooks. If there is a profile parameterisation you require that does not exist, please do create an issue.

Both ConstrainPaxisIp and ConstrainBetapIp are parameterised as follows: $\(J_{p}(\psi, R, Z) = \lambda\big[ \beta_{0} \frac{R}{R_{0}} \left( 1-\tilde{\psi}^{\alpha_m} \right)^{\alpha_n} + (1-\beta_{0}) \frac{R_0}{R} \left( 1-\tilde{\psi}^{\alpha_m} \right)^{\alpha_n} \big] \quad (R,Z) \in \Omega_p,\)$

where the first term is the pressure profile and the second is the toroidal current profile. Here, \(\tilde{\psi}\) denotes the normalised flux: $\( \tilde{\psi} = \frac{\psi - \psi_a}{\psi_b - \psi_a}, \)\( where \)\psi_a\( and \)\psi_b$ are the values of the flux on the magnetic axis and plasma boundary, respectively.

The parameters required to define this particular profile are:

• Ip (total plasma current).

• fvac (\(rB_{tor}\), vacuum toroidal field strength).

• alpha_m>0, and alpha_n>0 (that define the shape/peakedness of the profiles).

• If ConstrainPaxisIp is used, then paxis (pressure on the magnetic axis) is required.

• If ConstrainBetapIp is used, then betap (proxy of the poloidal beta) is required.

The values of \(\lambda\) and \(\beta_0\) are found using the above parameters as constraints (\(R_0\) is a fixed scaling constant) in the following:

• For ConstrainPaxisIp, we can re-arrange the following equations to solve for the unknowns:

\[ p_{\text{axis}} = \lambda \beta_{0} \frac{R}{R_{0}} \int^{\psi_b}_{\psi_a} \left( 1-\tilde{\psi}^{\alpha_m} \right)^{\alpha_n} \mathrm{d} \tilde{\psi} \]

and

\[ I_p = \int^{Z_{\text{max}}}_{Z_{\text{min}}} \int^{R_{\text{max}}}_{R_{\text{min}}} J_p(\psi, R, Z) \ \mathrm{d}R \mathrm{d}Z. \]

• For ConstrainBetapIp, we can instead re-arrange and solve the following:

\[ \beta_p = \frac{8 \pi}{\mu_0 I_p^2} \int^{Z_{\text{max}}}_{Z_{\text{min}}} \int^{R_{\text{max}}}_{R_{\text{min}}} p(\psi) \ \mathrm{d}R \mathrm{d}Z. \]

and

\[I_p = \int^{Z_{\text{max}}}_{Z_{\text{min}}} \int^{R_{\text{max}}}_{R_{\text{min}}} J_p(\psi, R, Z) \ \mathrm{d}R \mathrm{d}Z. \]

In what follows, we use ConstrainPaxisIp. Note that the equilibrium (eq) object is passed to the profile to inform calculations relating to the machine description.

# initialise the profiles
from freegsnke.jtor_update import ConstrainPaxisIp
profiles = ConstrainPaxisIp(
    eq=eq,        # equilibrium object
    paxis=8e3,    # profile object
    Ip=6e5,       # plasma current
    fvac=0.5,     # fvac = rB_{tor}
    alpha_m=1.8,  # profile function parameter
    alpha_n=1.2   # profile function parameter
)

Load the static nonlinear solver

We can now load FreeGSNKE’s Grad-Shafranov static solver. The equilibrium is used to inform the solver of the computational domain and of the tokamak properties. The solver below can be used for both inverse and forward solves.

Note: It’s not necessary to instantiate a new solver when aiming to use it on new or different equilibria, as long as the integration domain, mesh grid, and tokamak are consistent across solves.

from freegsnke import GSstaticsolver
GSStaticSolver = GSstaticsolver.NKGSsolver(eq)    

Constraints

Recall that in the inverse solve mode we seek to estimate the active poloidal field coil currents using user-defined constraints (e.g. on isoflux, null points, and psi values) and plasma current density profiles for a desired equilibrium shape.

FreeGSNKE uses a constrain object, which accepts constraints on the location of:

• any null points, which are either X-points or O-points (null_points).

• any pairs of points that lie on the same flux surface (isoflux_set). The flux at these locations will be constrained to be the same. Multiple sets can be defined.

• any known flux points (psi_vals), i.e. one can specify \(\psi\) at any location \((R,Z)\), possibly an entire flux map (not used here).

At least one constraint (preferably many more) is required to carry out an inverse solve.

Here, we specify two X-point locations (we want a double null plasma) and a number of isoflux locations. The isofluxes here will define the core plasma shape and part of the divertor legs.

import numpy as np 

Rx = 0.6      # X-point radius
Zx = 1.1      # X-point height
Ra = .85
Rout = 1.4    # outboard midplane radius
Rin = 0.34    # inboard midplane radius

# set desired null_points locations
# this can include X-point and O-point locations
null_points = [[Rx, Rx], [Zx, -Zx]]

# set desired isoflux constraints with format 
# isoflux_set = [isoflux_0, isoflux_1 ... ] 
# with each isoflux_i = [R_coords, Z_coords]
isoflux_set = np.array([[[Rx, Rx, Rin, Rout, 1.3, 1.3, .8,.8], [Zx, -Zx, 0.,0., 2.1, -2.1,1.62,-1.62]]])
           
# instantiate the freegsnke constrain object
from freegsnke.inverse import Inverse_optimizer
constrain = Inverse_optimizer(null_points=null_points,
                              isoflux_set=isoflux_set)

Given that there may be more or less constraints than unknown parameters (i.e. coil currents), the inverse problem may be over- or under-constrained. This means that there may be zero, one, or many solutions to the problem.

A number of quadratic (i.e.Tikhonov) regularization parameters will be used to combat this. Larger values will encourage lower absolute coil current values. It is sometimes useful to experiment with different values to explore whether the converged solution departs from the desired constraints.

The linear system of constraints

During an inverse solve, a minmisation problem involving the responses, changes in coil currents, and constraints, is repeatedly solved:

\[ \min_{x} \| A\vec{x} - \vec{b}\|^2 + \| \gamma \vec{x} \|^2, \]

where

• \(A\) is the fixed response matrix (that determines how a change in coil currents \(x\), affects constraint values).

• \(\vec{x} = \Delta \vec{I}^c\) is the step change in active coil currents required to match the constraints.

• \(\vec{b}\) is the vector of constraint values being enforced.

• \(\gamma > 0\) is the regularisation parameter/vector.

We solve for \(\vec{x}\) using a gradient-based optimiser.

Song et al. (2024) provide a nice overview of the inverse problem.

Fixed coil currents

It is also possible to set values for the current in specific active poloidal field coils. The inverse solver will not allow the current in these coils to vary if the control parameter is set to False (i.e. it will be excluded during the optimisation).

Note: any passive structures in the tokamak automatically have their control parameter set to False and are therefore not included in an inverse solve.

As an example, we will fix the Solenoid current and seek a solution in which this value is fixed, rather than estimated by the inverse solve.

eq.tokamak.set_coil_current('Solenoid', 5000)
eq.tokamak['Solenoid'].control = False  # ensures the current in the Solenoid is fixed

The inverse solve

The following cell will execute the solve. Since a constrain object is provided, this is interpreted as a call to the inverse solver, if constrain=None, then the forward solver will be called (see next notebook). The target_relative_tolerance is the maximum relative error on the plasma flux function allowed for convergence and target_relative_psit_update ensures that the relative update to the plasma flux (caused by the update in the control currents) is lower than this target value. Both are required to be met for the inverse problem to be considered successfully solved.

The verbose=True option will provide an indication of the progression of the solve.

The solver also needs stabilisation for any up-down symmetric coils used for vertical control of the plasma (in this case the P6 coil). To do this, one can set the l2_reg parameter: the Tikonov regularisation used by the optmiser of the coil currents. Here, P6 is the last coil in the list of those available for control, and so we use a larger regularisation for it specifically. This is the length of the coils being controlled (so it excludes Solenoid).

The solver steps are (roughly):

1. Solve the linear system to find initial coil currents that approximately satisfy the constraints for the initial equilibrium.

2. While tolerance is not met:

   • use the coil currents to solve the forward GS problem (with NK iterations).

   • solve the linear system to update coil currents to satisfy constraints for current equilibrium (from forward solve).

   • check tolerances.

GSStaticSolver.solve(eq=eq, 
                     profiles=profiles, 
                     constrain=constrain, 
                     target_relative_tolerance=1e-6,
                     target_relative_psit_update=1e-3,
                     verbose=True, # print output
                     l2_reg=np.array([1e-12]*10+[1e-6]), 
                     )

-----
Inverse static solve starting...

Initial guess for plasma_psi successful, residual found.
Initial relative error = 1.02e+00
-----
Iteration: 0
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-5.76e+01', '2.98e+04', '1.11e+04', '1.30e+04', '-7.91e+03', '2.73e+04', '-1.35e+04', '-6.50e+03', '-3.55e+04', '4.80e+03', '-7.74e-10']
Constraint losses = 7.38e-01
Relative update of tokamak psi (in plasma core): 2.53e+01
Handing off to forward solve (with updated currents).

Relative error =  9.31e-01
-----
Iteration: 1
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['2.01e+03', '-7.50e+03', '-2.81e+03', '-3.41e+03', '2.30e+03', '-7.39e+03', '3.88e+03', '1.88e+03', '9.76e+03', '-2.45e+03', '-5.64e-10']
Constraint losses = 1.70e-01
Relative update of tokamak psi (in plasma core): 2.62e-02
Handing off to forward solve (with updated currents).
Relative error =  8.67e-01
-----
Iteration: 2
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['1.42e+03', '-5.47e+03', '-2.05e+03', '-2.48e+03', '1.68e+03', '-5.37e+03', '2.82e+03', '1.37e+03', '7.09e+03', '-1.85e+03', '-2.06e-10']
Constraint losses = 1.23e-01
Relative update of tokamak psi (in plasma core): 2.61e-02
Handing off to forward solve (with updated currents).
Relative error =  7.77e-01
-----
Iteration: 3
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['9.70e+02', '-3.99e+03', '-1.50e+03', '-1.81e+03', '1.22e+03', '-3.91e+03', '2.06e+03', '9.94e+02', '5.15e+03', '-1.41e+03', '5.14e-10']
Constraint losses = 8.82e-02
Relative update of tokamak psi (in plasma core): 2.52e-02
Handing off to forward solve (with updated currents).
Relative error =  6.52e-01
-----
Iteration: 4
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['6.28e+02', '-2.90e+03', '-1.10e+03', '-1.33e+03', '8.95e+02', '-2.84e+03', '1.50e+03', '7.19e+02', '3.73e+03', '-1.09e+03', '4.07e-10']
Constraint losses = 6.33e-02
Relative update of tokamak psi (in plasma core): 2.32e-02
Handing off to forward solve (with updated currents).
Relative error =  4.80e-01
-----
Iteration: 5
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['3.51e+02', '-2.09e+03', '-8.13e+02', '-9.87e+02', '6.57e+02', '-2.06e+03', '1.09e+03', '5.15e+02', '2.68e+03', '-8.48e+02', '5.22e-10']
Constraint losses = 4.53e-02
Relative update of tokamak psi (in plasma core): 1.96e-02
Handing off to forward solve (with updated currents).
Relative error =  2.48e-01
-----
Iteration: 6
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['8.95e+01', '-1.46e+03', '-6.02e+02', '-7.45e+02', '4.83e+02', '-1.48e+03', '7.77e+02', '3.58e+02', '1.88e+03', '-6.78e+02', '5.04e-10']
Constraint losses = 3.22e-02
Relative update of tokamak psi (in plasma core): 1.43e-02
Handing off to forward solve (with updated currents).
Relative error =  3.66e-02
-----
Iteration: 7
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.04e+02', '-6.33e+02', '-3.50e+02', '-4.77e+02', '2.75e+02', '-7.56e+02', '3.95e+02', '1.60e+02', '8.84e+02', '-4.44e+02', '3.12e-13']
Constraint losses = 1.83e-02
Relative update of tokamak psi (in plasma core): 9.07e-03
Handing off to forward solve (with updated currents).

Relative error =  5.49e-05
-----
Iteration: 8
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.01e+02', '1.02e+02', '-2.87e+01', '-5.33e+01', '7.63e+00', '4.84e+01', '-1.15e+01', '-2.18e+01', '-1.15e+02', '-1.70e+02', '2.61e-12']
Constraint losses = 4.19e-02
Relative update of tokamak psi (in plasma core): 1.88e-02
Handing off to forward solve (with updated currents).
Relative error =  2.08e-05
-----
Iteration: 9
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-9.83e+01', '1.85e+02', '3.79e+01', '-3.53e+00', '-1.75e+01', '6.61e+01', '-6.91e+01', '-4.49e+01', '-1.41e+02', '1.92e+02', '2.11e-10']
Constraint losses = 6.93e-03
Relative update of tokamak psi (in plasma core): 6.83e-03
Handing off to forward solve (with updated currents).
Relative error =  5.46e-05
-----
Iteration: 10
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.90e+02', '1.14e+01', '-3.92e+01', '-4.67e+01', '1.40e+01', '7.64e+00', '2.14e+01', '3.47e+00', '-2.33e+01', '-1.87e+02', '1.04e-09']
Constraint losses = 1.85e-02
Relative update of tokamak psi (in plasma core): 1.41e-02
Handing off to forward solve (with updated currents).

Relative error =  9.56e-05
-----
Iteration: 11
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['2.21e+02', '-5.95e+01', '2.91e+01', '4.44e+01', '-8.91e+00', '-3.01e+01', '-2.76e+00', '9.24e+00', '6.87e+01', '1.64e+02', '8.28e-09']
Constraint losses = 2.21e-02
Relative update of tokamak psi (in plasma core): 1.25e-02
Handing off to forward solve (with updated currents).
Relative error =  3.66e-05
-----
Iteration: 12
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.66e+02', '3.57e+01', '-2.50e+01', '-3.53e+01', '8.12e+00', '1.84e+01', '6.29e+00', '-4.39e+00', '-4.33e+01', '-1.32e+02', '1.09e-07']
Constraint losses = 1.83e-02
Relative update of tokamak psi (in plasma core): 1.12e-02
Handing off to forward solve (with updated currents).

Relative error =  5.56e-05
-----
Iteration: 13
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['1.69e+02', '-3.41e+01', '2.51e+01', '3.44e+01', '-8.10e+00', '-1.90e+01', '-6.07e+00', '4.81e+00', '4.53e+01', '1.40e+02', '1.59e-06']
Constraint losses = 1.36e-02
Relative update of tokamak psi (in plasma core): 1.05e-02
Handing off to forward solve (with updated currents).
Relative error =  5.69e-05
-----
Iteration: 14
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.52e+02', '3.07e+01', '-2.31e+01', '-3.18e+01', '7.51e+00', '1.64e+01', '6.08e+00', '-3.87e+00', '-3.92e+01', '-1.24e+02', '1.23e-05']
Constraint losses = 2.01e-02
Relative update of tokamak psi (in plasma core): 1.05e-02
Handing off to forward solve (with updated currents).
Relative error =  5.74e-05
-----
Iteration: 15
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['1.41e+02', '-2.42e+01', '2.21e+01', '2.90e+01', '-7.26e+00', '-1.43e+01', '-6.55e+00', '3.19e+00', '3.51e+01', '1.22e+02', '8.48e-05']
Constraint losses = 9.89e-03
Relative update of tokamak psi (in plasma core): 9.11e-03
Handing off to forward solve (with updated currents).

Relative error =  2.09e-05
-----
Iteration: 16
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.34e+02', '2.63e+01', '-2.04e+01', '-2.79e+01', '6.65e+00', '1.43e+01', '5.50e+00', '-3.34e+00', '-3.43e+01', '-1.10e+02', '4.48e-06']
Constraint losses = 1.92e-02
Relative update of tokamak psi (in plasma core): 9.24e-03
Handing off to forward solve (with updated currents).
Relative error =  3.67e-05
-----
Iteration: 17
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['1.29e+02', '-2.26e+01', '2.01e+01', '2.66e+01', '-6.62e+00', '-1.32e+01', '-5.88e+00', '2.97e+00', '3.22e+01', '1.11e+02', '1.63e-04']
Constraint losses = 7.41e-03
Relative update of tokamak psi (in plasma core): 8.40e-03
Handing off to forward solve (with updated currents).

Relative error =  2.11e-05
-----
Iteration: 18
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.10e+02', '2.10e+01', '-1.70e+01', '-2.31e+01', '5.57e+00', '1.16e+01', '4.75e+00', '-2.66e+00', '-2.80e+01', '-9.23e+01', '4.80e-06']
Constraint losses = 1.91e-02
Relative update of tokamak psi (in plasma core): 7.70e-03
Handing off to forward solve (with updated currents).
Relative error =  3.95e-05
-----
Iteration: 19
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['5.40e+01', '-6.23e+00', '9.12e+00', '1.11e+01', '-3.08e+00', '-4.46e+00', '-3.39e+00', '7.36e-01', '1.18e+01', '5.03e+01', '1.04e-04']
Constraint losses = 3.21e-03
Relative update of tokamak psi (in plasma core): 3.75e-03
Handing off to forward solve (with updated currents).

Relative error =  1.45e-06
-----
Iteration: 20
Using full Jacobian (of constraints wrt coil currents) to optimsise currents.
 - calculating derivatives for coil 1/11

 - calculating derivatives for coil 2/11

 - calculating derivatives for coil 3/11

 - calculating derivatives for coil 4/11

 - calculating derivatives for coil 5/11

 - calculating derivatives for coil 6/11

 - calculating derivatives for coil 7/11

 - calculating derivatives for coil 8/11

 - calculating derivatives for coil 9/11

 - calculating derivatives for coil 10/11

 - calculating derivatives for coil 11/11

Change in coil currents (being controlled): ['-8.57e-01', '-8.93e-01', '-6.63e-01', '-6.29e-01', '-2.36e+00', '-3.92e+00', '-3.88e+00', '-4.84e+00', '-7.69e+00', '-1.79e+01', '1.86e-04']
Constraint losses = 1.34e-03
Relative update of tokamak psi (in plasma core): 1.86e-03
Handing off to forward solve (with updated currents).
Relative error =  1.28e-06
-----
Iteration: 21
Using full Jacobian (of constraints wrt coil currents) to optimsise currents.
 - calculating derivatives for coil 1/11

 - calculating derivatives for coil 2/11

 - calculating derivatives for coil 3/11

 - calculating derivatives for coil 4/11

 - calculating derivatives for coil 5/11

 - calculating derivatives for coil 6/11

 - calculating derivatives for coil 7/11

 - calculating derivatives for coil 8/11

 - calculating derivatives for coil 9/11

 - calculating derivatives for coil 10/11

 - calculating derivatives for coil 11/11

Change in coil currents (being controlled): ['-4.39e-01', '-2.45e-01', '-5.63e-02', '1.43e-01', '-2.76e-01', '-6.49e-01', '-7.43e-01', '-1.08e+00', '-2.26e+00', '-7.02e+00', '5.62e-06']
Constraint losses = 7.96e-04
Relative update of tokamak psi (in plasma core): 6.10e-04
Handing off to forward solve (with updated currents).
Relative error =  6.32e-07
-----
Inverse static solve SUCCESS. Tolerance 6.32e-07 (vs. requested 1e-06) reached in 22/100 iterations.

The following plots show how to display:

1. the tokamak with:

   • active coil filaments (rectangles with blue interior)

   • passive structures (blue circles if defined as filaments or thin black outline/grey interior if defined as parallelograms)

   • limiter/wall (solid black)

2. the tokamak + the equilibrium with:

   • separatrix/last closed flux surface (solid red lines)

   • poloidal flux (yellow/green/blue contours, colours indicates magnitude)

   • X-points (red circles)

   • O-points (green crosses)

3. the tokamak + the equilibrum + the contraints with:

   • null-point constraints (purple Y-shaped crosses)

   • isoflux contour constraints (usual crosses)

Note: Setting ‘show=True’ can toggle the legend on/off.

fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 8), dpi=80)

ax1.grid(zorder=0, alpha=0.75)
ax1.set_aspect('equal')
eq.tokamak.plot(axis=ax1,show=False)                                                          # plots the active coils and passive structures
ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)   # plots the limiter
ax1.set_xlim(0.1, 2.15)
ax1.set_ylim(-2.25, 2.25)

ax2.grid(zorder=0, alpha=0.75)
ax2.set_aspect('equal')
eq.tokamak.plot(axis=ax2,show=False)                                                          # plots the active coils and passive structures
ax2.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)   # plots the limiter
eq.plot(axis=ax2,show=False)                                                                  # plots the equilibrium
ax2.set_xlim(0.1, 2.15)
ax2.set_ylim(-2.25, 2.25)


ax3.grid(zorder=0, alpha=0.75)
ax3.set_aspect('equal')
eq.tokamak.plot(axis=ax3,show=False)                                                          # plots the active coils and passive structures
ax3.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)   # plots the limiter
eq.plot(axis=ax3,show=False)                                                                  # plots the equilibrium
constrain.plot(axis=ax3, show=True)                                                          # plots the contraints
ax3.set_xlim(0.1, 2.15)
ax3.set_ylim(-2.25, 2.25)

plt.tight_layout()

<Figure size 640x480 with 0 Axes>

A solve call will modify the equilibrium object in place. That means that certain quantities within the object will be updated as a result of the solve.

Various different quantities and functions can be accessed via the ‘eq’ and ‘profile’ objects. For example:

• the total flux can be accessed with eq.psi().

• the plasma flux with eq.plasma_psi.

• the active coil + passive structure flux with eq.tokamak.psi(eq.R, eq.Z).

• (Total flux = plasma flux + coil flux)

Explore eq. to see more (also profiles., e.g. the plasma current distribution over the domain can be found with profiles.jtor). Also see the example3 notebook for how to access many different equilibrium-derived quantities of interest.

The set of optimised coil currents can be retrieved using eq.tokamak.getCurrents() having been assigned to the equilibrium object during the inverse solve.

The following lines will save the calculated currents to a pickle file (we will use these in future notebooks).

inverse_current_values = eq.tokamak.getCurrents()

# save coil currents to file
import pickle
with open('simple_diverted_currents_PaxisIp.pk', 'wb') as f:
    pickle.dump(obj=inverse_current_values, file=f)

Second inverse solve: limited plasma example

Below we carry out an inverse solve seeking coil current values for a limited plasma configuration (rather than a diverted one).

In a limiter configuration the plasma “touches” the limiter of the tokamak and is confined by the solid structures of the vessel. The last closed flux surface (LCFS) is the closed contour that is furthest from the magnetic axis that just barely touches (and is tangent to) the limiter.

We instantiate a new equilibrium but use the same profile (illustrating how to vary the total plasma current and pressure on axis within it) and solver objects.

# a new eq object resets the coil currents and equilibrium
eq = equilibrium_update.Equilibrium(
    tokamak=tokamak,      # provide tokamak object
    Rmin=0.1, Rmax=2.0,   # radial range
    Zmin=-2.2, Zmax=2.2,  # vertical range
    nx=65,                # number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)
    ny=129,               # number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)
    # psi=plasma_psi
)

# modify the total plasma current
profiles.Ip = 4e5

# modify the pressure on the magnetic axis
profiles.paxis = 6e3

# first we specify some alternative constraints
Rout = 1.4   # outboard midplane radius
Rin = 0.24   # inboard midplane radius

# locations of X-points
Rx = 0.45
Zx = 1.18
null_points = [[Rx, Rx], [Zx, -Zx]]

# isoflux constraints
isoflux_set = [[Rx, Rx, Rout, Rin, Rin, Rin, Rin, Rin, .75, .75, .85, .85, 1.3, 1.3],
           [Zx, -Zx,  0,  0,  .1, -.1, .2, -.2, 1.6, -1.6, 1.7, -1.7, 2.1, -2.1]]

# let's assume we're also seeking an equilibrium with no solenoid current
eq.tokamak.set_coil_current('Solenoid', 0)
eq.tokamak['Solenoid'].control = False # fixes the current

# pass the magnetic constraints to a new constrain object
constrain = Inverse_optimizer(null_points=null_points,
                              isoflux_set=isoflux_set)

# carry out the solve
GSStaticSolver.solve(eq=eq, 
                     profiles=profiles, 
                     constrain=constrain, 
                     target_relative_tolerance=1e-6,
                     target_relative_psit_update=1e-3,
                     verbose=True, # print output
                     l2_reg=np.array([1e-9]*10+[1e-5]), 
                     )

-----
Inverse static solve starting...
Initial guess for plasma_psi successful, residual found.
Initial relative error = 9.72e-01
-----
Iteration: 0
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.

Discarding 'primary' Xpoint! Please check final result
Change in coil currents (being controlled): ['7.00e+03', '8.04e+02', '2.84e+04', '3.23e+04', '-1.77e+04', '2.89e+04', '-2.57e+04', '-1.29e+04', '-4.43e+04', '1.54e+04', '-5.66e-05']
Constraint losses = 8.85e-01
Relative update of tokamak psi (in plasma core): 1.58e+00
Handing off to forward solve (with updated currents).
Relative error =  9.53e-01
-----
Iteration: 1
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.89e+03', '-1.94e+03', '-8.93e+03', '-1.01e+04', '5.62e+03', '-9.48e+03', '7.97e+03', '3.99e+03', '1.42e+04', '-3.89e+03', '-5.20e-05']
Constraint losses = 3.03e-01
Relative update of tokamak psi (in plasma core): 2.14e-02
Handing off to forward solve (with updated currents).
Relative error =  9.13e-01
-----
Iteration: 2
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.98e+03', '-1.34e+03', '-6.56e+03', '-7.23e+03', '3.98e+03', '-6.88e+03', '5.81e+03', '2.91e+03', '1.03e+04', '-2.90e+03', '-4.78e-05']
Constraint losses = 2.15e-01
Relative update of tokamak psi (in plasma core): 2.38e-02
Handing off to forward solve (with updated currents).
Relative error =  8.55e-01
-----
Iteration: 3
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.45e+03', '-9.54e+02', '-4.68e+03', '-5.15e+03', '2.83e+03', '-4.90e+03', '4.14e+03', '2.07e+03', '7.33e+03', '-2.13e+03', '-4.42e-05']
Constraint losses = 1.52e-01
Relative update of tokamak psi (in plasma core): 2.59e-02
Handing off to forward solve (with updated currents).
Relative error =  7.72e-01
-----
Iteration: 4
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.07e+03', '-6.80e+02', '-3.34e+03', '-3.67e+03', '2.02e+03', '-3.50e+03', '2.96e+03', '1.48e+03', '5.22e+03', '-1.59e+03', '-4.12e-05']
Constraint losses = 1.06e-01
Relative update of tokamak psi (in plasma core): 2.73e-02
Handing off to forward solve (with updated currents).
Relative error =  6.56e-01
-----
Iteration: 5
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-8.17e+02', '-4.83e+02', '-2.38e+03', '-2.61e+03', '1.43e+03', '-2.50e+03', '2.11e+03', '1.05e+03', '3.72e+03', '-1.20e+03', '-3.89e-05']
Constraint losses = 7.43e-02
Relative update of tokamak psi (in plasma core): 2.75e-02
Handing off to forward solve (with updated currents).

Relative error =  4.94e-01
-----
Iteration: 6
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-6.46e+02', '-3.40e+02', '-1.70e+03', '-1.85e+03', '1.01e+03', '-1.78e+03', '1.50e+03', '7.44e+02', '2.63e+03', '-9.33e+02', '-3.74e-05']
Constraint losses = 5.17e-02
Relative update of tokamak psi (in plasma core): 2.59e-02
Handing off to forward solve (with updated currents).
Relative error =  2.80e-01
-----
Iteration: 7
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-5.41e+02', '-2.29e+02', '-1.18e+03', '-1.27e+03', '6.94e+02', '-1.24e+03', '1.04e+03', '5.08e+02', '1.81e+03', '-7.36e+02', '-3.76e-05']
Constraint losses = 3.54e-02
Relative update of tokamak psi (in plasma core): 2.07e-02
Handing off to forward solve (with updated currents).
Relative error =  6.97e-02
-----
Iteration: 8
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.71e+02', '-5.78e+00', '1.69e+01', '-4.28e+01', '1.31e+02', '-1.43e+02', '3.36e+01', '-1.78e+01', '4.23e+01', '-1.11e+02', '-4.57e-09']
Constraint losses = 2.42e-02
Relative update of tokamak psi (in plasma core): 9.57e-03
Handing off to forward solve (with updated currents).

Relative error =  1.31e-02
-----
Iteration: 9
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.37e+01', '-4.48e+00', '2.32e+00', '6.05e+00', '1.60e+01', '-8.74e+00', '2.41e-01', '-1.19e+01', '-3.52e+01', '-1.75e+02', '-1.17e-03']
Constraint losses = 1.13e-01
Relative update of tokamak psi (in plasma core): 2.19e-02
Handing off to forward solve (with updated currents).

Relative error =  1.90e-02
-----
Iteration: 10
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.04e+02', '-1.66e+01', '-2.90e+00', '-9.98e+00', '5.01e+01', '-3.30e+01', '1.79e+01', '-7.90e+00', '-1.85e+01', '-2.38e+02', '-1.02e-03']
Constraint losses = 4.03e-02
Relative update of tokamak psi (in plasma core): 2.08e-02
Handing off to forward solve (with updated currents).

Relative error =  4.89e-04
-----
Iteration: 11
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-6.79e+01', '-2.13e+01', '-2.00e+01', '-4.94e+01', '3.02e+01', '-3.23e+01', '3.34e+01', '3.21e+01', '1.09e+02', '2.49e+02', '-1.48e-04']
Constraint losses = 4.41e-02
Relative update of tokamak psi (in plasma core): 1.94e-02
Handing off to forward solve (with updated currents).

Relative error =  7.26e-03
-----
Iteration: 12
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-4.71e+01', '-9.60e+00', '-4.35e+00', '-3.28e+00', '1.34e+01', '-1.45e+01', '1.07e-01', '-1.26e+01', '-3.02e+01', '-1.74e+02', '-1.68e-04']
Constraint losses = 5.28e-02
Relative update of tokamak psi (in plasma core): 1.78e-02
Handing off to forward solve (with updated currents).
Relative error =  3.53e-05
-----
Iteration: 13
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-9.41e+01', '-2.88e+01', '-2.77e+01', '-5.40e+01', '2.69e+01', '-3.97e+01', '2.71e+01', '1.93e+01', '8.02e+01', '1.16e+02', '4.74e-05']
Constraint losses = 1.92e-02
Relative update of tokamak psi (in plasma core): 1.08e-02
Handing off to forward solve (with updated currents).

Relative error =  2.44e-04
-----
Iteration: 14
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-4.66e+01', '-1.09e+01', '-7.99e+00', '-7.97e+00', '9.54e+00', '-1.37e+01', '6.29e-01', '-1.03e+01', '-2.29e+01', '-1.51e+02', '-9.00e-05']
Constraint losses = 3.44e-02
Relative update of tokamak psi (in plasma core): 1.45e-02
Handing off to forward solve (with updated currents).

Relative error =  4.45e-05
-----
Iteration: 15
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-6.52e+01', '-2.13e+01', '-2.44e+01', '-4.66e+01', '1.99e+01', '-2.36e+01', '2.76e+01', '2.52e+01', '8.32e+01', '1.52e+02', '1.99e-06']
Constraint losses = 2.13e-02
Relative update of tokamak psi (in plasma core): 1.38e-02
Handing off to forward solve (with updated currents).
Relative error =  2.15e-04
-----
Iteration: 16
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.76e+01', '-5.83e+00', '-4.38e+00', '-1.94e+00', '4.46e+00', '-7.08e+00', '-2.10e+00', '-1.01e+01', '-2.58e+01', '-1.33e+02', '1.01e-04']
Constraint losses = 4.18e-02
Relative update of tokamak psi (in plasma core): 1.35e-02
Handing off to forward solve (with updated currents).

Relative error =  1.45e-04
-----
Iteration: 17
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-5.21e+01', '-1.57e+01', '-1.97e+01', '-3.35e+01', '1.27e+01', '-1.77e+01', '1.71e+01', '1.26e+01', '4.45e+01', '4.33e+01', '1.16e-04']
Constraint losses = 1.14e-02
Relative update of tokamak psi (in plasma core): 4.84e-03
Handing off to forward solve (with updated currents).
Relative error =  9.34e-06
-----
Iteration: 18
Using full Jacobian (of constraints wrt coil currents) to optimsise currents.

 - calculating derivatives for coil 1/11

 - calculating derivatives for coil 2/11

 - calculating derivatives for coil 3/11

 - calculating derivatives for coil 4/11

 - calculating derivatives for coil 5/11

 - calculating derivatives for coil 6/11

 - calculating derivatives for coil 7/11

 - calculating derivatives for coil 8/11

 - calculating derivatives for coil 9/11

 - calculating derivatives for coil 10/11

 - calculating derivatives for coil 11/11

Change in coil currents (being controlled): ['-4.59e+00', '-8.53e-01', '-1.35e+00', '-2.57e+00', '1.55e+00', '-6.58e+00', '-1.08e+00', '-4.63e+00', '-1.56e+01', '-2.39e+01', '-2.56e-03']
Constraint losses = 8.93e-03
Relative update of tokamak psi (in plasma core): 2.90e-03
Handing off to forward solve (with updated currents).
Relative error =  1.86e-05
-----
Iteration: 19
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-4.49e+01', '-1.17e+01', '-1.53e+01', '-2.24e+01', '1.12e+01', '-1.07e+01', '1.28e+01', '6.45e+00', '1.95e+01', '-3.40e+01', '1.72e-07']
Constraint losses = 9.13e-03
Relative update of tokamak psi (in plasma core): 1.84e-03
Handing off to forward solve (with updated currents).

Relative error =  2.58e-06
-----
Iteration: 20
Using full Jacobian (of constraints wrt coil currents) to optimsise currents.
 - calculating derivatives for coil 1/11

 - calculating derivatives for coil 2/11

 - calculating derivatives for coil 3/11

 - calculating derivatives for coil 4/11

 - calculating derivatives for coil 5/11

 - calculating derivatives for coil 6/11

 - calculating derivatives for coil 7/11

 - calculating derivatives for coil 8/11

 - calculating derivatives for coil 9/11

 - calculating derivatives for coil 10/11

 - calculating derivatives for coil 11/11

Change in coil currents (being controlled): ['-3.59e+00', '-8.77e-01', '-1.53e+00', '-2.44e+00', '2.27e+00', '5.03e-01', '3.82e+00', '3.44e+00', '8.00e+00', '6.43e+00', '-1.56e-04']
Constraint losses = 7.76e-03
Relative update of tokamak psi (in plasma core): 9.89e-04
Handing off to forward solve (with updated currents).
Relative error =  8.08e-07
-----
Inverse static solve SUCCESS. Tolerance 8.08e-07 (vs. requested 1e-06) reached in 21/100 iterations.

# save the currents for later use
inverse_current_values = eq.tokamak.getCurrents()

# save coil currents to file
with open('simple_limited_currents_PaxisIp.pk', 'wb') as f:
    pickle.dump(obj=inverse_current_values, file=f)

# plot the resulting equilbria 
fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)
ax1.grid(True, which='both')
eq.plot(axis=ax1, show=False)
eq.tokamak.plot(axis=ax1, show=False)
constrain.plot(axis=ax1,show=True)
ax1.set_xlim(0.1, 2.15)
ax1.set_ylim(-2.25, 2.25)
plt.tight_layout()

<Figure size 640x480 with 0 Axes>