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 Amorisco et al. (2024) and 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="../machine_configs/MAST-U/MAST-U_like_active_coils.pickle",
    passive_coils_path="../machine_configs/MAST-U/MAST-U_like_passive_coils.pickle",
    limiter_path="../machine_configs/MAST-U/MAST-U_like_limiter.pickle",
    wall_path="../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). Though it does not need to encompass any/all of the active coils or passive structures.

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.

The following table indicates which parameterisations are currently available in FreeGSNKE and for which types of equilibrium simulation they currently work for:

Name	Static	Evolutive (linear)	Evolutive (nonlinear)
ConstrainPaxisIp	✅	✅	✅
ConstrainBetapIp	✅	✅	✅
FiestaTopeol	✅	✅	✅
Lao85	✅	✅	✅
TensionSpline	✅	❌	❌
GeneralPprimeFFprime	✅	❌	❌

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 (\(R B_{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,    # pressure on the magnetic axis
    Ip=6e5,       # plasma current
    fvac=0.5,     # fvac = R B_{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 static solver object. 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.

Note: A random seed can also be set here. This comes into play for tricky forward/inverse solves where the internal Newton-Krylov solver gets stuck and needs to search in random Krylov direction. See code documentation for further details.

from freegsnke import GSstaticsolver
GSStaticSolver = GSstaticsolver.NKGSsolver(
    eq=eq,    # eq object
    # seed=42,  # random seed
    )    

Fixed coil currents

Firstly, we set the values for any currents in specific active poloidal field coils that we may know. 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

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) together with prescribed plasma current density profiles, in order to obtain a desired equilibrium shape.

FreeGSNKE uses a constrain object, which accepts the following types of constraints:

• Null points (null_points)
  A null point \((R_X, Z_X)\) (X-point or O-point) is defined by vanishing poloidal field, i.e. we can impose $\( \nabla \psi(R_X, Z_X) = \mathbf{0}. \)\( In particular, an X-point is a saddle point of \)\psi\( while an O-point (magnetic axis) is a local extremum (typically a minimum) of \)\psi\(. The above is equivalent to both: \)\( \frac{\partial \psi}{\partial R}(R_X, Z_X) = \frac{\partial \psi}{\partial Z}(R_X, Z_X) = 0, \)\( or \)\( B_r (R_X, Z_X) = B_z(R_X, Z_X) = 0. \)$ This means that for each null point considered, two constraints are imposed internall within the inverse solver.

• Isoflux sets (isoflux_set)
  An isoflux set is a set of points \(\{(R_i, Z_i)\}_{i=1}^k\) that lie on the same flux surface, i.e. we can impose $\( \psi(R_1, Z_1) = \psi(R_2, Z_2) = \cdots = \psi(R_k, Z_k). \)\( Internally, the inverse solver enforces each constraint such that \)\( \psi(R_i, Z_i) - \psi(R_j, Z_j) = 0, \qquad i \neq j, \)\( meaning that for \)k\( points in an isoflux set, there are \)k(k-1)/2$ constraints.

• Coil current limits (coil_current_limits)
  If we wnt to enforce bounds on the coil currents (for machine safety purposes), we can impose $\( I_m^{\min} \le I_m \le I_m^{\max}, \qquad m = 1, \dots, N. \)$

Note: see the next notebooks for some more advanced constraints!

At least one constraint (preferably more) is required to carry out an inverse solve. Although it may be tempting to utilise all of the above constraints, one must be careful not to overconstrain the inverse problem, else constraints will fight one another and the solver will not converge to an equilibrium.

In the following, we start simple and specify two X-point locations (to obtain a double-null plasma), a few isoflux locations, and some coil current limits. The isoflux constraints define the core plasma boundary and part of the divertor legs.

import numpy as np 
from freegsnke.inverse import Inverse_optimizer

Rx = 0.6      # X-point radius
Zx = 1.1      # X-point height
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.0, 1.0, .8,.8], [Zx, -Zx, 0.,0., 2.0, -2.0, 1.62, -1.62]]
    ])

# set the coil current limits (upper and lower)
# coil ordering in this case: PX,  D1,  D2,  D3,  Dp,  D5,  D6,  D7,  P4,  P5,  P6
coil_current_limits = [
    [5e3, 9e3, 9e3, 7e3, 7e3, 5e3, 4e3, 5e3, 0.0, 0.0, None],
    [-5e3, -9e3, -9e3, -7e3, -7e3, -5e3, -4e3, -5e3, -10e3, -10e3, None]
]
           
# instantiate the freegsnke constrain object
constrain = Inverse_optimizer(
    null_points=null_points,
    isoflux_set=isoflux_set,
    coil_current_limits=coil_current_limits
)

# if you find coil limits are being violated or an undesireable solution is being produced,
# you can try increasing the penalty factor for violating the coil limits
# (here we just set it to its default value of 1e5)
constrain.mu_coils = 1e5

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.

The coil limits are handled in a slightly different manner because they are box constraints that cannot be encoded in the minimisation problem. Instead, FreeGSNKE adds additional constraints to the solver that bound the solution \(\Delta \vec{I}^c\). These bounds have a tolerance (slack) which allows the coils to violate their limits, but penalises violations; this allows the coil limits to be violated ‘on the path’ to a solution that respects the prescribed limits.

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 l2_reg parameter defines the Tikonov regularisation (i.e. \(\gamma\) in the cell above) used by the optmiser. This can be set as a scalar or as a vector (equal to the number of coil currents being solved for). Larger values force coil currents to stay closer to their original values while lower values will encourage more freedom. This can be useful in particular for vertical control coils (e.g. the P6 coil in MAST-U), in which we don’t want the coil current to “jump around” during optimisation and cause vertical instability issues.

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 Newton-Krylov iterations).

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

   • check tolerances.

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. It is often a bit of an art when solving inverse problems so the best strategy is to start with some basic isoflux constraints before trying to add null points and others.

# 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-12]*10+[1e-6]), 
                     )

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

Successfully computed GS residual using initial plasma_psi and tokamak_psi guesses.
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): ['-7.66e+02', '3.01e+04', '1.51e+04', '1.98e+04', '-9.06e+03', '8.18e+03', '-1.02e+04', '-1.04e+04', '-2.63e+04', '6.08e+03', '2.98e-10']
Constraint losses = 7.05e-01
Relative update of tokamak psi (in plasma core): 2.49e+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): ['1.91e+03', '-7.46e+03', '-2.48e+03', '-5.86e+03', '1.76e+03', '-1.27e+03', '3.30e+03', '2.69e+03', '7.32e+03', '-2.99e+03', '-6.19e-10']
Constraint losses = 1.62e-01
Relative update of tokamak psi (in plasma core): 2.66e-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.40e+03', '-5.45e+03', '-2.30e+03', '-4.19e+03', '1.59e+03', '-8.16e+02', '2.21e+03', '1.89e+03', '5.25e+03', '-2.21e+03', '-6.09e-10']
Constraint losses = 1.14e-01
Relative update of tokamak psi (in plasma core): 2.67e-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): ['1.05e+03', '-4.26e+03', '-1.26e+03', '-3.31e+03', '5.19e+02', '-9.54e+02', '4.18e+03', '6.37e+02', '2.75e+03', '-1.44e+03', '-2.46e-10']
Constraint losses = 7.97e-02
Relative update of tokamak psi (in plasma core): 2.47e-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): ['1.06e+03', '-3.28e+03', '-3.05e+03', '-2.47e+03', '1.72e+03', '-3.50e+02', '1.62e+03', '2.25e+03', '8.27e+02', '-1.01e+03', '4.82e-11']
Constraint losses = 5.65e-02
Relative update of tokamak psi (in plasma core): 2.01e-02
Handing off to forward solve (with updated currents).
Relative error =  4.79e-01
-----
Iteration: 5
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['2.20e+02', '-1.70e+03', '-9.20e+02', '-3.61e+02', '5.89e+02', '-2.50e+03', '-7.34e+02', '2.03e+03', '3.38e+03', '-1.04e+03', '-1.36e-10']
Constraint losses = 4.05e-02
Relative update of tokamak psi (in plasma core): 2.33e-02
Handing off to forward solve (with updated currents).
Relative error =  2.46e-01
-----
Iteration: 6
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.49e+02', '-1.01e+03', '-3.68e+02', '-1.30e+02', '2.97e+02', '-1.90e+03', '-9.16e+02', '1.31e+03', '2.91e+03', '-8.80e+02', '-3.41e-10']
Constraint losses = 2.84e-02
Relative update of tokamak psi (in plasma core): 1.93e-02
Handing off to forward solve (with updated currents).
Relative error =  3.83e-02
-----
Iteration: 7
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-4.00e+02', '-5.12e+02', '-3.18e+02', '-6.33e+02', '1.58e+02', '-1.55e+02', '3.20e+02', '2.24e+02', '6.05e+02', '-4.73e+02', '2.62e-14']
Constraint losses = 1.54e-02
Relative update of tokamak psi (in plasma core): 9.82e-03
Handing off to forward solve (with updated currents).

Relative error =  2.08e-05
-----
Iteration: 8
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.38e+02', '1.46e+02', '-4.16e+01', '1.07e+01', '-1.21e+01', '-5.25e+00', '-1.01e+01', '-3.44e+01', '-1.15e+02', '-1.52e+02', '2.04e-12']
Constraint losses = 4.47e-02
Relative update of tokamak psi (in plasma core): 1.98e-02
Handing off to forward solve (with updated currents).
Relative error =  2.80e-05
-----
Iteration: 9
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.99e+02', '3.34e+02', '3.55e+00', '5.00e+01', '-4.28e+01', '2.01e+01', '-7.03e+01', '-7.07e+01', '-1.96e+02', '1.65e+02', '2.12e-10']
Constraint losses = 4.76e-03
Relative update of tokamak psi (in plasma core): 3.09e-03
Handing off to forward solve (with updated currents).

Relative error =  3.80e-06
-----
Iteration: 10
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): ['-1.18e+00', '-6.61e-01', '-1.30e-01', '2.33e-01', '-1.67e-01', '-2.86e+00', '-1.68e+00', '-2.80e+00', '-4.92e+00', '-1.52e+01', '3.75e-06']
Constraint losses = 5.48e-03
Relative update of tokamak psi (in plasma core): 1.40e-03
Handing off to forward solve (with updated currents).
Relative error =  7.77e-06
-----
Iteration: 11
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): ['-7.63e-01', '-1.35e-01', '3.01e-01', '6.73e-01', '1.48e+00', '9.39e-02', '1.09e+00', '5.67e-01', '1.43e-01', '-4.89e+00', '3.13e-06']
Constraint losses = 2.11e-03
Relative update of tokamak psi (in plasma core): 2.57e-04
Handing off to forward solve (with updated currents).
Relative error =  7.77e-07
-----
Inverse static solve SUCCESS. Tolerance 7.77e-07 (vs. requested 1e-06) reached in 12/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.

• (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_currents = eq.tokamak.getCurrents()
inverse_currents_arr = eq.tokamak.getCurrentsVec()[0:12]

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

Finally, we can check that the currents we have found respect the coil limits that we prescribed earlier:

# storage
coil_names = []
active_currents = []
min_limits = []
max_limits = []

for coil_name, coil_current, ul, ll in zip(
    eq.tokamak.coils_dict,
    eq.tokamak.getCurrentsVec()[:12],
    [None] + coil_current_limits[0],   # upper limits
    [None] + coil_current_limits[1],   # lower limits
):
    # skip solenoid (no limits)
    if ul is None or ll is None:
        continue

    coil_names.append(coil_name)
    active_currents.append(coil_current)
    min_limits.append(ll)
    max_limits.append(ul)
    
    if (coil_current >= ll) & (coil_current <= ul):
        within_limit = True
    else:
        within_limit = False

    # print
    print(f"{coil_name}: {coil_current:.2f} --> [{ll},{ul}] [A] --> limits met = {within_limit}")

# convert to numpy
active_currents = np.array(active_currents)
min_limits = np.array(min_limits)
max_limits = np.array(max_limits)

PX: 3696.81 --> [-5000.0,5000.0] [A] --> limits met = True
D1: 6941.24 --> [-9000.0,9000.0] [A] --> limits met = True
D2: 4390.86 --> [-9000.0,9000.0] [A] --> limits met = True
D3: 2924.36 --> [-7000.0,7000.0] [A] --> limits met = True
Dp: -2488.20 --> [-7000.0,7000.0] [A] --> limits met = True
D5: 249.19 --> [-5000.0,5000.0] [A] --> limits met = True
D6: -314.99 --> [-4000.0,4000.0] [A] --> limits met = True
D7: 511.08 --> [-5000.0,5000.0] [A] --> limits met = True
P4: -3564.38 --> [-10000.0,0.0] [A] --> limits met = True
P5: -3954.23 --> [-10000.0,0.0] [A] --> limits met = True

We can also visualise the currents with respect to their limits (normalising between -1 and 1). Note the Solenoid and P6 are missing as they have no limits in this case.

# VISUALISE HOW CLOSE THEY ARE TO THE LIMITS (normalised)

# normalized to [-1, 1] scale (relative to min/max)
norm_current = 2 * (active_currents - min_limits) / (max_limits - min_limits) - 1

# plot
plt.figure(figsize=(12, 5), dpi=70)
x = np.arange(len(coil_names))

# line plot
plt.plot(x, norm_current, marker='x', linestyle='-', color='red', label="Current")
plt.axhline(-1, color='k', linestyle='-', linewidth=1, label="Min. limit")
plt.axhline(1, color='k', linestyle='-', linewidth=1, label="Max. limit")

# formatting
plt.xticks(x, coil_names, rotation=45, ha='right')
plt.ylabel('Coil name [Amps]')
plt.ylabel('Normalised current')
plt.title('Coil current proximity to limits (normalised)')
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend()
plt.tight_layout()
plt.show()

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.

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

# initialise the profiles
profiles = ConstrainPaxisIp(
    eq=eq,        # equilibrium object
    paxis=6e3,    # pressure on the magnetic axis  <-- changed
    Ip=4e5,       # plasma current                 <-- changed
    fvac=0.5,     # fvac = R B_{tor}
    alpha_m=1.8,  # profile function parameter
    alpha_n=1.2   # profile function parameter
)

# first we specify some alternative constraints
Rout = 1.4    # outboard midplane radius
Rin = 0.255   # inboard midplane radius (less than the wall radius to make it limited)

# 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],
    [Zx, -Zx,  0,  0,  .1, -.1, .2, -.2, 1.6, -1.6, 1.7, -1.7]
    ]

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

# also make it perfectly up/down symmetric (as P6 is wired in anti-series vs. the other coils)
eq.tokamak.set_coil_current('P6', 0)
eq.tokamak['P6'].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),
                     )

-----
Inverse static solve starting...
Successfully computed GS residual using initial plasma_psi and tokamak_psi guesses.
Initial relative error = 9.73e-01
-----
Iteration: 0
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['1.04e+04', '1.80e+02', '1.99e+04', '4.19e+04', '-1.14e+04', '3.17e+03', '-1.98e+04', '-1.57e+04', '-3.84e+04', '1.83e+04']
Constraint losses = 8.36e-01
Relative update of tokamak psi (in plasma core): 1.43e+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): ['-3.55e+03', '-2.11e+03', '-6.36e+03', '-1.33e+04', '3.55e+03', '-9.29e+02', '6.18e+03', '5.02e+03', '1.25e+04', '-5.01e+03']
Constraint losses = 2.80e-01
Relative update of tokamak psi (in plasma core): 1.90e-02
Handing off to forward solve (with updated currents).
Relative error =  9.12e-01
-----
Iteration: 2
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.46e+03', '-1.51e+03', '-4.60e+03', '-9.64e+03', '2.55e+03', '-6.80e+02', '4.46e+03', '3.62e+03', '9.00e+03', '-3.69e+03']
Constraint losses = 1.99e-01
Relative update of tokamak psi (in plasma core): 2.15e-02
Handing off to forward solve (with updated currents).
Relative error =  8.53e-01
-----
Iteration: 3
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.80e+03', '-1.08e+03', '-3.28e+03', '-6.87e+03', '1.81e+03', '-4.89e+02', '3.18e+03', '2.58e+03', '6.41e+03', '-2.69e+03']
Constraint losses = 1.40e-01
Relative update of tokamak psi (in plasma core): 2.38e-02
Handing off to forward solve (with updated currents).
Relative error =  7.71e-01
-----
Iteration: 4
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.

Change in coil currents (being controlled): ['-1.33e+03', '-7.65e+02', '-2.34e+03', '-4.89e+03', '1.28e+03', '-3.53e+02', '2.27e+03', '1.84e+03', '4.56e+03', '-1.98e+03']
Constraint losses = 9.83e-02
Relative update of tokamak psi (in plasma core): 2.58e-02
Handing off to forward solve (with updated currents).
Relative error =  6.54e-01
-----
Iteration: 5
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.01e+03', '-5.41e+02', '-1.67e+03', '-3.47e+03', '9.07e+02', '-2.56e+02', '1.61e+03', '1.30e+03', '3.24e+03', '-1.48e+03']
Constraint losses = 6.85e-02
Relative update of tokamak psi (in plasma core): 2.71e-02
Handing off to forward solve (with updated currents).
Relative error =  4.91e-01
-----
Iteration: 6
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-8.03e+02', '-3.77e+02', '-1.18e+03', '-2.44e+03', '6.32e+02', '-1.87e+02', '1.14e+03', '9.16e+02', '2.28e+03', '-1.12e+03']
Constraint losses = 4.74e-02
Relative update of tokamak psi (in plasma core): 2.67e-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): ['-6.69e+02', '-2.45e+02', '-7.98e+02', '-1.63e+03', '4.12e+02', '-1.34e+02', '7.60e+02', '6.05e+02', '1.51e+03', '-8.32e+02']
Constraint losses = 3.15e-02
Relative update of tokamak psi (in plasma core): 2.23e-02
Handing off to forward solve (with updated currents).
Relative error =  6.90e-02
-----
Iteration: 8
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.64e+02', '-3.14e+01', '-1.71e+01', '-2.70e+01', '6.31e+01', '9.46e-01', '5.20e+01', '3.22e+01', '4.54e+01', '-1.34e+02']
Constraint losses = 2.04e-02
Relative update of tokamak psi (in plasma core): 7.65e-03
Handing off to forward solve (with updated currents).

Relative error =  2.28e-02
-----
Iteration: 9
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.54e+01', '-3.59e+00', '1.23e+00', '-1.36e-01', '8.77e+00', '-2.17e+01', '-7.62e+00', '-1.81e+01', '-3.21e+01', '-1.31e+02']
Constraint losses = 1.08e-01
Relative update of tokamak psi (in plasma core): 2.13e-02
Handing off to forward solve (with updated currents).

Relative error =  1.98e-03
-----
Iteration: 10
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-5.45e+01', '-9.65e+00', '-3.09e+00', '-6.38e+00', '1.59e+01', '-2.60e+01', '-2.91e+00', '-1.73e+01', '-3.06e+01', '-1.69e+02']
Constraint losses = 5.29e-02
Relative update of tokamak psi (in plasma core): 2.05e-02
Handing off to forward solve (with updated currents).
Relative error =  5.00e-05
-----
Iteration: 11
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.22e+02', '-3.61e+01', '-3.79e+01', '-4.68e+01', '1.02e+01', '2.26e+01', '4.11e+01', '4.44e+01', '8.56e+01', '9.53e+01']
Constraint losses = 1.91e-02
Relative update of tokamak psi (in plasma core): 1.30e-02
Handing off to forward solve (with updated currents).

Relative error =  7.67e-05
-----
Iteration: 12
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-4.50e+01', '-7.74e+00', '-4.61e+00', '-8.52e+00', '8.78e+00', '-2.35e+01', '-5.25e+00', '-1.65e+01', '-2.67e+01', '-1.39e+02']
Constraint losses = 4.11e-02
Relative update of tokamak psi (in plasma core): 1.68e-02
Handing off to forward solve (with updated currents).
Relative error =  5.73e-05
-----
Iteration: 13
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-9.00e+01', '-2.61e+01', '-3.32e+01', '-4.28e+01', '8.26e-01', '2.02e+01', '3.13e+01', '3.71e+01', '7.50e+01', '1.08e+02']
Constraint losses = 1.75e-02
Relative update of tokamak psi (in plasma core): 1.31e-02
Handing off to forward solve (with updated currents).

Relative error =  2.85e-05
-----
Iteration: 14
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.02e+01', '-3.98e+00', '-2.85e+00', '-6.24e+00', '4.98e+00', '-1.93e+01', '-6.19e+00', '-1.46e+01', '-2.33e+01', '-1.08e+02']
Constraint losses = 3.83e-02
Relative update of tokamak psi (in plasma core): 1.34e-02
Handing off to forward solve (with updated currents).
Relative error =  6.77e-05
-----
Iteration: 15
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-7.06e+01', '-1.59e+01', '-2.44e+01', '-3.43e+01', '-1.02e+00', '5.09e+00', '1.66e+01', '1.80e+01', '4.16e+01', '3.64e+01']
Constraint losses = 1.07e-02
Relative update of tokamak psi (in plasma core): 5.22e-03
Handing off to forward solve (with updated currents).
Relative error =  4.49e-05
-----
Iteration: 16
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-5.59e+01', '-6.74e+00', '-1.16e+01', '-2.02e+01', '2.96e+00', '-2.04e+01', '-2.83e+00', '-1.12e+01', '-1.22e+01', '-1.04e+02']
Constraint losses = 1.53e-02
Relative update of tokamak psi (in plasma core): 1.14e-02
Handing off to forward solve (with updated currents).

Relative error =  7.97e-04
-----
Iteration: 17
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.54e+01', '-5.60e+00', '-1.14e+01', '-1.61e+01', '2.34e+00', '1.73e+01', '1.85e+01', '2.36e+01', '4.38e+01', '8.77e+01']
Constraint losses = 2.14e-02
Relative update of tokamak psi (in plasma core): 9.95e-03
Handing off to forward solve (with updated currents).
Relative error =  5.30e-05
-----
Iteration: 18
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.61e+01', '-2.66e+00', '-7.73e+00', '-1.45e+01', '-1.17e-01', '-1.71e+01', '-5.13e+00', '-1.12e+01', '-1.33e+01', '-8.09e+01']
Constraint losses = 1.46e-02
Relative update of tokamak psi (in plasma core): 9.33e-03
Handing off to forward solve (with updated currents).

Relative error =  1.85e-04
-----
Iteration: 19
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.53e+01', '-3.81e+00', '-1.06e+01', '-1.61e+01', '2.04e+00', '1.31e+01', '1.56e+01', '1.94e+01', '3.71e+01', '6.83e+01']
Constraint losses = 1.50e-02
Relative update of tokamak psi (in plasma core): 8.17e-03
Handing off to forward solve (with updated currents).
Relative error =  8.03e-05
-----
Iteration: 20
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.44e+01', '-4.08e-01', '-5.25e+00', '-1.07e+01', '-1.06e+00', '-1.37e+01', '-5.21e+00', '-9.71e+00', '-1.19e+01', '-6.23e+01']
Constraint losses = 1.33e-02
Relative update of tokamak psi (in plasma core): 7.37e-03
Handing off to forward solve (with updated currents).

Relative error =  2.52e-04
-----
Iteration: 21
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.59e+01', '-2.05e+00', '-1.00e+01', '-1.65e+01', '1.99e+00', '9.47e+00', '1.32e+01', '1.59e+01', '3.15e+01', '5.10e+01']
Constraint losses = 1.03e-02
Relative update of tokamak psi (in plasma core): 6.46e-03
Handing off to forward solve (with updated currents).
Relative error =  4.69e-05
-----
Iteration: 22
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.86e+01', '8.15e-01', '-3.84e+00', '-8.66e+00', '-1.09e+00', '-1.21e+01', '-5.11e+00', '-8.99e+00', '-1.14e+01', '-5.45e+01']
Constraint losses = 1.23e-02
Relative update of tokamak psi (in plasma core): 6.49e-03
Handing off to forward solve (with updated currents).

Relative error =  1.72e-04
-----
Iteration: 23
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.41e+01', '-3.65e-01', '-9.11e+00', '-1.60e+01', '2.40e+00', '8.16e+00', '1.24e+01', '1.45e+01', '2.90e+01', '4.41e+01']
Constraint losses = 8.30e-03
Relative update of tokamak psi (in plasma core): 5.76e-03
Handing off to forward solve (with updated currents).
Relative error =  2.90e-05
-----
Iteration: 24
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.27e+01', '1.29e+00', '-2.51e+00', '-6.26e+00', '-8.99e-01', '-9.58e+00', '-4.34e+00', '-7.37e+00', '-9.61e+00', '-4.30e+01']
Constraint losses = 1.16e-02
Relative update of tokamak psi (in plasma core): 5.15e-03
Handing off to forward solve (with updated currents).

Relative error =  2.14e-05
-----
Iteration: 25
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.97e+01', '1.10e+00', '-6.97e+00', '-1.34e+01', '1.70e+00', '3.27e+00', '7.81e+00', '8.53e+00', '1.84e+01', '2.01e+01']
Constraint losses = 5.45e-03
Relative update of tokamak psi (in plasma core): 2.95e-03
Handing off to forward solve (with updated currents).
Relative error =  1.18e-05
-----
Iteration: 26
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.64e+01', '2.19e+00', '-4.30e+00', '-9.93e+00', '-6.53e-01', '-8.27e+00', '-2.13e+00', '-4.73e+00', '-3.98e+00', '-3.47e+01']
Constraint losses = 6.26e-03
Relative update of tokamak psi (in plasma core): 3.88e-03
Handing off to forward solve (with updated currents).
Relative error =  2.18e-05
-----
Iteration: 27
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.06e+01', '1.23e+00', '-4.49e+00', '-8.67e+00', '1.45e+00', '4.79e+00', '6.88e+00', '8.21e+00', '1.63e+01', '2.61e+01']
Constraint losses = 6.73e-03
Relative update of tokamak psi (in plasma core): 3.42e-03
Handing off to forward solve (with updated currents).

Relative error =  2.10e-05
-----
Iteration: 28
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.23e+01', '2.40e+00', '-3.38e+00', '-8.19e+00', '-5.49e-01', '-6.60e+00', '-1.67e+00', '-3.71e+00', '-2.94e+00', '-2.72e+01']
Constraint losses = 5.38e-03
Relative update of tokamak psi (in plasma core): 3.01e-03
Handing off to forward solve (with updated currents).
Relative error =  1.91e-05
-----
Iteration: 29
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-8.99e+00', '1.71e+00', '-3.85e+00', '-7.86e+00', '1.16e+00', '3.32e+00', '5.47e+00', '6.38e+00', '1.31e+01', '1.91e+01']
Constraint losses = 5.25e-03
Relative update of tokamak psi (in plasma core): 2.58e-03
Handing off to forward solve (with updated currents).
Relative error =  1.16e-05
-----
Iteration: 30
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.

Change in coil currents (being controlled): ['-9.29e+00', '2.40e+00', '-2.76e+00', '-6.89e+00', '-4.15e-01', '-5.01e+00', '-1.07e+00', '-2.60e+00', '-1.62e+00', '-2.01e+01']
Constraint losses = 4.46e-03
Relative update of tokamak psi (in plasma core): 2.20e-03
Handing off to forward solve (with updated currents).
Relative error =  1.67e-05
-----
Iteration: 31
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-7.78e+00', '2.09e+00', '-3.41e+00', '-7.34e+00', '8.46e-01', '2.04e+00', '4.27e+00', '4.83e+00', '1.04e+01', '1.31e+01']
Constraint losses = 4.06e-03
Relative update of tokamak psi (in plasma core): 1.85e-03
Handing off to forward solve (with updated currents).
Relative error =  9.49e-06
-----
Iteration: 32
Using full Jacobian (of constraints wrt coil currents) to optimsise currents.

 - calculating derivatives for coil 1/10
 - calculating derivatives for coil 2/10

 - calculating derivatives for coil 3/10
 - calculating derivatives for coil 4/10

 - calculating derivatives for coil 5/10
 - calculating derivatives for coil 6/10

 - calculating derivatives for coil 7/10
 - calculating derivatives for coil 8/10

 - calculating derivatives for coil 9/10
 - calculating derivatives for coil 10/10

Change in coil currents (being controlled): ['-6.06e-01', '3.77e-01', '-3.75e-01', '-9.38e-01', '6.19e-02', '-8.90e-01', '-5.84e-03', '-3.87e-01', '-1.38e+00', '-4.28e+00']
Constraint losses = 2.71e-03
Relative update of tokamak psi (in plasma core): 4.26e-04
Handing off to forward solve (with updated currents).
Relative error =  7.72e-07
-----
Inverse static solve SUCCESS. Tolerance 7.72e-07 (vs. requested 1e-06) reached in 33/100 iterations.

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

# save coil currents to file
with open('data/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>