Example: Advanced constraints in static inverse free-boundary equilibrium calculations

---

In this example notebook, we demonstrate some of the more advanced types of constraints and techniques that can be employed in the inverse solver.

Instatiate the objects

As before, we start by instantiating the machine, equilibrium, profiles, and solver objects.

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

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
)

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

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

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.

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

(Normalised) poloidal flux constraints at specific locations

In addition to the null points, isoflux set, and coil_current_limits constraints introduced in the previous notebook, additional methods are available for more advanced control of the inverse problem.

We can also constrain:

• Flux values (psi_vals)
  If the flux is known at a location \((R_j, Z_j)\), we can impose $\( \psi(R_j, Z_j) = \psi_j^{\text{target}}. \)\( More generally, one could prescribe \)\psi(R,Z)$ over a region (e.g. a full flux map). Flux values are typically difficult to estimate a priori to simulation and so the following constraint is often more useful.

• Normalised flux values (psi_norm_limits)
  At a location \((R_j, Z_j)\), we can impose upper (or lower) bound constraints on the normalised flux $\( \hat{\psi}(R,Z) = \frac{\psi(R,Z) - \psi_{axis}}{\psi_{boundary} - \psi_{axis}}, \)\( such that \)\( \hat{\psi}(R_j, Z_j) \leq \hat{\psi}_j^{\text{target}}. \)\( or \)\( \hat{\psi}(R_j, Z_j) \geq \hat{\psi}_j^{\text{target}}. \)\( Note that \)\psi_{axis}\( and \)\psi_{boundary}$ are the values of the flux on the magnetic axis and plasma boundary, respectively. As we’ll see, this can be useful for setting explicit constraints on flux behaviour near the wall if there are certain safety limits to adhere to.

Once again, we stress that users should be mindful of specifying too many constraints that may conflict with one another during an inverse solve. For example, a given normalised \(\psi\) constraint may be impossible/difficult to achieve given other isoflux constraints or coil current limits. More generally, specifying too many constraints could make the problem ill-posed and reduce the overall quality of the final solution.

Recall that if you find that the solver is violating the limit (inequality) constraints (coil limits and normalised \(\psi\)), you should look to reduce the number of constraints or try modifying the mu_* parameters. These parameters control how much a violation of the constraint is penalised in the solver, and increasing this penalty (by increasing mu_*) will increase the likelihood the constraint is adequately satisfied. Again, it may also be that the constraint is impossible to satisfy given your other constraints.

In the following, we will show how to use the normalised flux constraints to control the normalised flux on the divertor nose.

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, weights]
isoflux_set = np.array([
    [
        [Rx, Rx, Rin, Rout], 
        [Zx, -Zx, 0., 0.],
    ]
])

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

# normalised psi constraints set with format:
# [R, Z, psiN, 1] --> ψ̂(R,Z) ≥ psiN  (the +1 or -1 defines ≥ or ≤)
# [R, Z, psiN, -1] --> ψ̂(R,Z) ≤ psiN
psi_norm_limits = [
    [0.82, -1.55, 1.2, 1],
    [0.82, -1.55, 1.3, -1],
]

# instantiate the freegsnke constrain object
constrain = Inverse_optimizer(
    null_points=null_points,
    isoflux_set=isoflux_set,
    coil_current_limits=coil_current_limits,
    psi_norm_limits=psi_norm_limits,
    mu_psi_norm = 1e7, # penalise how much the normalised psi constraint is violated (default 1e6)
)

We solve, as in the previous example, by passing the equilibrium, profiles, and constraints to the static solver. Again, we apply regularisation to encourage a solution with low coil currents.

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): ['-4.16e+03', '4.08e+04', '2.02e+04', '-2.40e+03', '-1.02e+04', '-7.31e+03', '-5.53e+03', '-8.38e+03', '-1.83e+04', '6.42e+03', '9.15e-01']
Constraint losses = 6.07e-01
Relative update of tokamak psi (in plasma core): 2.91e+01
Handing off to forward solve (with updated currents).
Relative error =  9.32e-01
-----
Iteration: 1
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.

Change in coil currents (being controlled): ['6.21e+03', '-2.00e+04', '-6.29e+03', '1.01e+04', '9.01e+03', '8.82e+02', '-1.12e+03', '3.90e+02', '2.35e+03', '-1.51e+03', '-1.15e-01']
Constraint losses = 1.41e-01
Relative update of tokamak psi (in plasma core): 1.71e-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): ['9.31e+02', '-4.57e+03', '-2.16e+03', '-9.39e+02', '-1.53e+03', '1.25e+03', '3.43e+03', '2.57e+03', '4.47e+03', '-2.68e+03', '-1.33e-01']
Constraint losses = 1.02e-01
Relative update of tokamak psi (in plasma core): 2.09e-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.02e+03', '-3.80e+03', '-1.46e+03', '-4.74e+03', '-4.12e+02', '1.82e+03', '5.39e+03', '1.14e+03', '9.53e+02', '-1.65e+03', '-2.15e-01']
Constraint losses = 7.24e-02
Relative update of tokamak psi (in plasma core): 1.83e-02
Handing off to forward solve (with updated currents).
Relative error =  6.51e-01
-----
Iteration: 4
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['4.74e+02', '-2.47e+03', '-9.30e+02', '-2.64e+03', '-6.69e+02', '2.38e+03', '6.19e+02', '4.10e+03', '7.13e+02', '-1.70e+03', '-2.37e-01']
Constraint losses = 5.18e-02
Relative update of tokamak psi (in plasma core): 1.93e-02
Handing off to forward solve (with updated currents).
Relative error =  4.78e-01
-----
Iteration: 5
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.

Change in coil currents (being controlled): ['-7.53e+01', '-9.83e+02', '-1.27e+03', '-1.45e+02', '-4.04e+01', '5.81e+02', '-1.03e+03', '4.94e+02', '3.76e+03', '-1.43e+03', '-2.08e-01']
Constraint losses = 3.67e-02
Relative update of tokamak psi (in plasma core): 2.41e-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): ['-3.13e+02', '5.31e+00', '-2.55e+03', '-2.04e+02', '3.51e+02', '4.25e+02', '-5.74e+02', '4.50e+02', '2.23e+03', '-1.01e+03', '-1.52e-01']
Constraint losses = 2.49e-02
Relative update of tokamak psi (in plasma core): 1.75e-02
Handing off to forward solve (with updated currents).
Relative error =  3.89e-02
-----
Iteration: 7
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.

Change in coil currents (being controlled): ['-4.56e+02', '-5.69e+02', '-2.73e+02', '-9.73e+01', '-2.97e+01', '9.08e+01', '2.86e+02', '2.51e+02', '4.32e+02', '-4.75e+02', '-1.49e-05']
Constraint losses = 1.20e-02
Relative update of tokamak psi (in plasma core): 1.04e-02
Handing off to forward solve (with updated currents).
Relative error =  2.86e-05
-----
Iteration: 8
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.

Change in coil currents (being controlled): ['-2.45e+02', '1.51e+01', '3.09e+01', '1.69e+01', '6.47e+01', '-3.92e+01', '-3.81e+01', '-5.70e+01', '-1.19e+02', '-1.04e+02', '-2.00e-06']
Constraint losses = 3.35e-02
Relative update of tokamak psi (in plasma core): 1.77e-02
Handing off to forward solve (with updated currents).
Relative error =  7.18e-05
-----
Iteration: 9
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-4.51e+02', '1.52e+02', '1.24e+02', '5.88e+01', '1.40e+02', '-7.30e+01', '-1.35e+02', '-1.53e+02', '-3.09e+02', '7.61e+01', '1.90e-04']
Constraint losses = 5.67e-03
Relative update of tokamak psi (in plasma core): 1.44e-02
Handing off to forward solve (with updated currents).

Relative error =  5.31e-05
-----
Iteration: 10
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['2.53e+02', '-1.67e+01', '-3.53e+01', '-2.09e+01', '-7.70e+01', '4.51e+01', '4.13e+01', '6.38e+01', '1.34e+02', '1.40e+02', '6.16e-05']
Constraint losses = 1.79e-02
Relative update of tokamak psi (in plasma core): 1.59e-02
Handing off to forward solve (with updated currents).
Relative error =  5.50e-05
-----
Iteration: 11
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.

Change in coil currents (being controlled): ['-1.76e+02', '-3.35e+01', '-1.33e+00', '3.39e+00', '4.25e+01', '-2.35e+01', '4.78e-01', '-1.86e+01', '-4.61e+01', '-1.59e+02', '-2.08e-06']
Constraint losses = 1.73e-02
Relative update of tokamak psi (in plasma core): 1.51e-02
Handing off to forward solve (with updated currents).
Relative error =  8.13e-05
-----
Iteration: 12
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['1.80e+02', '1.79e+01', '-7.76e+00', '-6.96e+00', '-4.65e+01', '2.88e+01', '1.12e+01', '3.00e+01', '6.66e+01', '1.49e+02', '7.39e-05']
Constraint losses = 1.10e-02
Relative update of tokamak psi (in plasma core): 1.34e-02
Handing off to forward solve (with updated currents).

Relative error =  4.76e-05
-----
Iteration: 13
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.47e+02', '-2.04e+01', '3.01e+00', '4.33e+00', '3.65e+01', '-2.15e+01', '-4.72e+00', '-2.02e+01', '-4.67e+01', '-1.24e+02', '4.57e-05']
Constraint losses = 1.84e-02
Relative update of tokamak psi (in plasma core): 1.26e-02
Handing off to forward solve (with updated currents).
Relative error =  5.63e-05
-----
Iteration: 14
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.

Change in coil currents (being controlled): ['1.12e+02', '2.41e+01', '2.32e+00', '-1.35e+00', '-2.61e+01', '1.72e+01', '3.90e-01', '1.33e+01', '3.12e+01', '1.15e+02', '-8.92e-07']
Constraint losses = 5.24e-03
Relative update of tokamak psi (in plasma core): 9.62e-03
Handing off to forward solve (with updated currents).
Relative error =  6.78e-05
-----
Iteration: 15
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.30e+02', '-2.08e+01', '1.13e+00', '3.19e+00', '3.17e+01', '-1.89e+01', '-2.72e+00', '-1.67e+01', '-3.92e+01', '-1.16e+02', '6.16e-05']
Constraint losses = 1.55e-02
Relative update of tokamak psi (in plasma core): 1.13e-02
Handing off to forward solve (with updated currents).

Relative error =  1.95e-05
-----
Iteration: 16
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['1.37e+02', '2.00e+01', '-2.47e+00', '-4.00e+00', '-3.44e+01', '2.06e+01', '4.49e+00', '1.92e+01', '4.42e+01', '1.20e+02', '-1.24e-05']
Constraint losses = 6.09e-03
Relative update of tokamak psi (in plasma core): 1.06e-02
Handing off to forward solve (with updated currents).

Relative error =  2.51e-05
-----
Iteration: 17
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.12e+02', '-1.66e+01', '1.63e+00', '2.99e+00', '2.75e+01', '-1.65e+01', '-3.13e+00', '-1.51e+01', '-3.49e+01', '-9.79e+01', '-5.85e-06']
Constraint losses = 1.67e-02
Relative update of tokamak psi (in plasma core): 9.68e-03
Handing off to forward solve (with updated currents).
Relative error =  1.67e-05
-----
Iteration: 18
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.

Change in coil currents (being controlled): ['3.48e+01', '8.54e+00', '1.25e+00', '-2.50e-01', '-7.90e+00', '5.07e+00', '-3.43e-01', '3.62e+00', '8.74e+00', '3.59e+01', '-2.23e-05']
Constraint losses = 1.82e-03
Relative update of tokamak psi (in plasma core): 2.98e-03
Handing off to forward solve (with updated currents).
Relative error =  2.86e-05
-----
Iteration: 19
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-9.15e+01', '-1.16e+01', '2.37e+00', '2.85e+00', '2.28e+01', '-1.35e+01', '-3.41e+00', '-1.30e+01', '-2.99e+01', '-7.67e+01', '-9.69e-05']
Constraint losses = 4.50e-03
Relative update of tokamak psi (in plasma core): 7.27e-03
Handing off to forward solve (with updated currents).

Relative error =  2.64e-05
-----
Iteration: 20
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['8.97e+01', '9.68e+00', '-3.48e+00', '-3.43e+00', '-2.33e+01', '1.34e+01', '4.45e+00', '1.37e+01', '3.11e+01', '7.19e+01', '-6.30e-06']
Constraint losses = 9.94e-03
Relative update of tokamak psi (in plasma core): 6.46e-03
Handing off to forward solve (with updated currents).
Relative error =  1.43e-05
-----
Iteration: 21
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.

Change in coil currents (being controlled): ['-7.73e+01', '-7.23e+00', '3.43e+00', '3.00e+00', '1.99e+01', '-1.16e+01', '-4.19e+00', '-1.21e+01', '-2.73e+01', '-6.06e+01', '-3.58e-05']
Constraint losses = 3.92e-03
Relative update of tokamak psi (in plasma core): 5.94e-03
Handing off to forward solve (with updated currents).
Relative error =  2.72e-05
-----
Iteration: 22
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['7.10e+01', '7.55e+00', '-2.80e+00', '-2.72e+00', '-1.84e+01', '1.06e+01', '3.61e+00', '1.09e+01', '2.48e+01', '5.68e+01', '-1.30e-05']
Constraint losses = 7.93e-03
Relative update of tokamak psi (in plasma core): 5.18e-03
Handing off to forward solve (with updated currents).

Relative error =  2.43e-05
-----
Iteration: 23
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-5.82e+01', '-5.52e+00', '2.53e+00', '2.23e+00', '1.50e+01', '-8.71e+00', '-3.02e+00', '-8.99e+00', '-2.04e+01', '-4.58e+01', '-7.46e-06']
Constraint losses = 3.08e-03
Relative update of tokamak psi (in plasma core): 4.43e-03
Handing off to forward solve (with updated currents).
Relative error =  1.02e-05
-----
Iteration: 24
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.

Change in coil currents (being controlled): ['5.04e+01', '5.68e+00', '-1.81e+00', '-1.85e+00', '-1.30e+01', '7.56e+00', '2.43e+00', '7.64e+00', '1.74e+01', '4.10e+01', '-5.19e-06']
Constraint losses = 5.88e-03
Relative update of tokamak psi (in plasma core): 3.74e-03
Handing off to forward solve (with updated currents).
Relative error =  7.52e-06
-----
Iteration: 25
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.40e-01', '-2.37e-01', '-2.54e-01', '-3.54e-01', '-8.30e-01', '-1.33e+00', '-1.25e+00', '-1.46e+00', '-2.12e+00', '-4.55e+00', '3.67e-01']
Constraint losses = 2.10e-03
Relative update of tokamak psi (in plasma core): 5.12e-04
Handing off to forward solve (with updated currents).

Relative error =  1.06e-06
-----
Iteration: 26
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): ['-9.27e-02', '-9.88e-02', '-9.04e-02', '-1.41e-01', '-3.02e-01', '-5.96e-01', '-5.10e-01', '-6.33e-01', '-9.30e-01', '-2.22e+00', '-2.19e-01']
Constraint losses = 1.08e-03
Relative update of tokamak psi (in plasma core): 2.37e-04
Handing off to forward solve (with updated currents).
Relative error =  3.06e-07
-----
Inverse static solve SUCCESS. Tolerance 3.06e-07 (vs. requested 1e-06) reached in 27/100 iterations.

import matplotlib.pyplot as plt

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

We can now check that the normalised \(\psi\) constraint has indeed been satisfied!

for psi_con in psi_norm_limits:
    sign = ">=" if psi_con[3] >= 0 else "<="
    print(f"Psi norm at ({psi_con[0]}, {psi_con[1]}) = {eq.psiNRZ(psi_con[0], psi_con[1]):.3f} ({sign} {psi_con[2]})")

Psi norm at (0.82, -1.55) = 1.216 (>= 1.2)
Psi norm at (0.82, -1.55) = 1.216 (<= 1.3)

Weighting constraints within the solver

It is possible to specify weights on the different types of constraints within the inverse solver, enabling the solver to prioritise certain constraints over others.

There are two ways of doing this:

1. Weighting different types of constraints relative to one another:

   Different classes of constraint can be weighted relative to each other (e.g. null points, isoflux, and psi constraints). This allows the solver to be instructed that one type of constraint is more important than others. For example, a user may want null point constraints to be strictly satisfied, but may be less concerned about the isoflux set being met with the same strictness.

2. Weighting individual isoflux constraints within a set:

   Individual isoflux constraints within a set can optionally be weighted relative to one another, instructing the solver to prioritise satisfying some subset over others.

   For example, consider a set of three isoflux constraints: one at the X-point, one at the outer midplane, and one at a desired strike point. The strike point constraint could be assigned a lower weight to allow the solver to find a solution with a strike point close to — but not necessarily exactly on — the target location.

Option 1: weighting different types of constraints

Here, we’ll show how to weight the null point constraints more than the isoflux constraints. This solver then prioritises fitting the null points exactly and the isoflux points less so.

# core constraints
Rx = 0.55     # X-point radius
Zx = 1.2      # 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, 0.75, 1.0, 0.75, 1.0], 
        [Zx, -Zx, 0.0, 0.0, -1.6, -2.0, 1.6, 2.0],
    ]
])

# instantiate the freegsnke constrain object
constrain = Inverse_optimizer(
    null_points=null_points,
    isoflux_set=isoflux_set,
    coil_current_limits=coil_current_limits,
    weight_isoflux=0.2, # <-- weight the isofluxes less than the null points (default 1.0)
    weight_nulls=1.0,   # <-- weight the null points more than the null points (default 1.0)
    # weight_psi=1.0    # <-- weight for psi values not used here
)

# solve!
GSStaticSolver.solve(eq=eq, 
                     profiles=profiles, 
                     constrain=constrain, 
                     target_relative_tolerance=1e-6,
                     target_relative_psit_update=1e-3,
                     verbose=True,
                     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 = 3.06e-07
-----
Iteration: 0
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['1.10e+03', '-2.50e+03', '1.55e+02', '9.46e+00', '5.15e+02', '-7.41e+01', '2.33e+02', '1.65e+02', '3.53e+02', '-2.27e+02', '1.53e-07']
Constraint losses = 8.92e-02
Relative update of tokamak psi (in plasma core): 9.80e-03
Handing off to forward solve (with updated currents).

Relative error =  7.99e-05
-----
Iteration: 1
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['1.53e+03', '-3.43e+03', '3.47e+02', '3.26e+02', '5.65e+02', '-2.31e+01', '1.97e+02', '1.55e+02', '2.84e+02', '-2.43e+02', '2.96e-06']
Constraint losses = 7.02e-02
Relative update of tokamak psi (in plasma core): 8.87e-03
Handing off to forward solve (with updated currents).
Relative error =  3.45e-05
-----
Iteration: 2
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.

Change in coil currents (being controlled): ['-1.13e+02', '4.49e+01', '1.40e+02', '5.24e+02', '-2.67e+02', '1.77e+02', '-1.59e+02', '-7.28e+01', '-2.82e+02', '1.05e+02', '-1.85e-06']
Constraint losses = 4.96e-02
Relative update of tokamak psi (in plasma core): 8.99e-03
Handing off to forward solve (with updated currents).
Relative error =  3.04e-05
-----
Iteration: 3
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.86e+02', '-2.41e+01', '-1.82e+02', '1.83e+02', '-1.41e+02', '1.80e+02', '9.96e+01', '1.18e+02', '7.33e+01', '-8.33e+01', '-8.51e-06']
Constraint losses = 4.06e-02
Relative update of tokamak psi (in plasma core): 4.57e-03
Handing off to forward solve (with updated currents).

Relative error =  2.70e-05
-----
Iteration: 4
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.22e+02', '7.84e+01', '6.85e+01', '3.45e+02', '-1.91e+02', '1.28e+02', '-9.66e+01', '-3.95e+01', '-1.79e+02', '6.50e+01', '-9.01e-06']
Constraint losses = 3.32e-02
Relative update of tokamak psi (in plasma core): 5.42e-03
Handing off to forward solve (with updated currents).
Relative error =  1.05e-05
-----
Iteration: 5
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.86e+02', '-5.74e+00', '-1.11e+02', '1.25e+02', '-9.48e+01', '1.16e+02', '5.73e+01', '7.06e+01', '3.77e+01', '-4.69e+01', '-6.13e-07']
Constraint losses = 2.72e-02
Relative update of tokamak psi (in plasma core): 2.78e-03
Handing off to forward solve (with updated currents).

Relative error =  8.14e-06
-----
Iteration: 6
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.83e-01', '3.99e-02', '-3.83e-01', '4.45e-01', '-2.77e+00', '-1.89e+00', '-2.92e+00', '-2.94e+00', '-4.50e+00', '-5.76e+00', '2.69e-05']
Constraint losses = 2.24e-02
Relative update of tokamak psi (in plasma core): 9.17e-04
Handing off to forward solve (with updated currents).
Relative error =  2.84e-07
-----
Inverse static solve SUCCESS. Tolerance 2.84e-07 (vs. requested 1e-06) reached in 7/100 iterations.

Notice how the null points constraints are satisfied pefectly while the isoflux constraints are “less so” (especially in the divertor region).

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

Option 2: weighting different isoflux constraints

Here, we’ll show how to weight individual isoflux constraints differently from one another. For example, here we want to prioritise the core shape over the exact location of the strikepoint.

# core constraints
Rx = 0.55     # X-point radius
Zx = 1.2      # 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, weights]
isoflux_set = np.array([
    [
        [Rx, Rx, Rin, Rout, 1.0, 1.0, 0.75, 1.0, 0.75, 1.0], 
        [Zx, -Zx, 0.0, 0.0, -0.9, 0.9,  -1.6, -2.0, 1.6, 2.0],
        [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.1, 1.0, 0.1],
    ]
])

# instantiate the freegsnke constrain object
constrain = Inverse_optimizer(
    null_points=null_points,
    isoflux_set=isoflux_set,
    coil_current_limits=coil_current_limits,
)

# 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 = 2.84e-07
-----
Iteration: 0
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-9.38e+01', '-6.35e+01', '-4.35e+02', '7.11e+01', '4.51e+02', '5.64e+01', '6.13e+01', '-1.15e+02', '-6.86e+02', '2.30e+02', '2.42e-07']
Constraint losses = 2.80e-02
Relative update of tokamak psi (in plasma core): 9.80e-03
Handing off to forward solve (with updated currents).

Relative error =  2.19e-05
-----
Iteration: 1
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['1.76e+02', '-2.52e+02', '-1.36e+03', '1.09e+02', '1.19e+03', '2.41e+02', '3.58e+02', '-1.31e+02', '-1.55e+03', '5.41e+02', '-8.06e-05']
Constraint losses = 1.70e-02
Relative update of tokamak psi (in plasma core): 8.37e-03
Handing off to forward solve (with updated currents).
Relative error =  1.94e-05
-----
Iteration: 2
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.

Change in coil currents (being controlled): ['3.25e+02', '-1.64e+02', '-5.12e+02', '-8.75e+01', '2.99e+02', '8.60e+01', '2.61e+02', '1.04e+02', '-1.44e+02', '-2.17e+02', '-2.43e-05']
Constraint losses = 1.21e-02
Relative update of tokamak psi (in plasma core): 8.84e-03
Handing off to forward solve (with updated currents).
Relative error =  1.27e-05
-----
Iteration: 3
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['1.54e+02', '1.83e+01', '-1.32e+02', '5.15e+01', '1.92e+01', '7.29e+01', '5.05e+01', '3.07e+01', '-6.35e+01', '1.05e+02', '-5.80e-06']
Constraint losses = 1.75e-02
Relative update of tokamak psi (in plasma core): 8.53e-03
Handing off to forward solve (with updated currents).

Relative error =  1.11e-05
-----
Iteration: 4
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-4.77e+01', '3.27e+01', '-5.16e+01', '4.24e+01', '-8.61e+00', '1.78e+01', '1.23e+00', '-1.02e+01', '-6.98e+01', '-7.72e+01', '-4.73e-06']
Constraint losses = 8.19e-03
Relative update of tokamak psi (in plasma core): 8.23e-03
Handing off to forward solve (with updated currents).
Relative error =  1.75e-05
-----
Iteration: 5
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['1.28e+02', '5.24e+00', '-1.22e+02', '4.34e+01', '1.67e+01', '6.77e+01', '5.20e+01', '3.43e+01', '-4.35e+01', '8.56e+01', '-4.57e-05']
Constraint losses = 1.41e-02
Relative update of tokamak psi (in plasma core): 7.94e-03
Handing off to forward solve (with updated currents).

Relative error =  4.13e-05
-----
Iteration: 6
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-7.67e+01', '1.34e+01', '-6.94e+00', '2.06e+01', '-1.67e+01', '9.43e-02', '-6.71e+00', '-8.26e+00', '-2.25e+01', '-8.16e+01', '1.12e-05']
Constraint losses = 9.29e-03
Relative update of tokamak psi (in plasma core): 7.48e-03
Handing off to forward solve (with updated currents).
Relative error =  2.29e-05
-----
Iteration: 7
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['1.02e+02', '7.90e+00', '-1.02e+02', '4.35e+01', '7.36e+00', '6.08e+01', '4.24e+01', '2.95e+01', '-3.68e+01', '7.58e+01', '1.59e-04']
Constraint losses = 1.13e-02
Relative update of tokamak psi (in plasma core): 7.00e-03
Handing off to forward solve (with updated currents).

Relative error =  5.42e-05
-----
Iteration: 8
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-7.91e+01', '1.20e+01', '9.99e+00', '1.64e+01', '-2.00e+01', '-5.07e+00', '-1.24e+01', '-1.00e+01', '-1.26e+01', '-7.12e+01', '-3.87e-05']
Constraint losses = 9.27e-03
Relative update of tokamak psi (in plasma core): 6.67e-03
Handing off to forward solve (with updated currents).
Relative error =  3.30e-05
-----
Iteration: 9
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['8.56e+01', '8.74e+00', '-8.98e+01', '4.24e+01', '2.72e+00', '5.57e+01', '3.70e+01', '2.64e+01', '-3.28e+01', '6.65e+01', '-3.33e-05']
Constraint losses = 9.08e-03
Relative update of tokamak psi (in plasma core): 6.25e-03
Handing off to forward solve (with updated currents).

Relative error =  1.88e-05
-----
Iteration: 10
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-7.17e+01', '9.32e+00', '1.49e+01', '1.21e+01', '-1.81e+01', '-7.12e+00', '-1.28e+01', '-9.67e+00', '-7.56e+00', '-6.19e+01', '6.86e-06']
Constraint losses = 9.13e-03
Relative update of tokamak psi (in plasma core): 5.82e-03
Handing off to forward solve (with updated currents).
Relative error =  2.29e-05
-----
Iteration: 11
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['7.11e+01', '9.74e+00', '-8.11e+01', '4.18e+01', '-7.76e-01', '5.17e+01', '3.27e+01', '2.38e+01', '-3.07e+01', '5.68e+01', '1.65e-04']
Constraint losses = 6.94e-03
Relative update of tokamak psi (in plasma core): 5.38e-03
Handing off to forward solve (with updated currents).

Relative error =  2.40e-05
-----
Iteration: 12
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-6.18e+01', '7.30e+00', '1.57e+01', '9.04e+00', '-1.55e+01', '-7.48e+00', '-1.19e+01', '-8.76e+00', '-4.92e+00', '-5.25e+01', '2.32e-06']
Constraint losses = 8.76e-03
Relative update of tokamak psi (in plasma core): 4.95e-03
Handing off to forward solve (with updated currents).
Relative error =  2.01e-05
-----
Iteration: 13
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['5.78e+01', '1.10e+01', '-7.51e+01', '4.24e+01', '-3.99e+00', '4.89e+01', '2.94e+01', '2.15e+01', '-3.03e+01', '4.74e+01', '4.70e-05']
Constraint losses = 5.05e-03
Relative update of tokamak psi (in plasma core): 4.53e-03
Handing off to forward solve (with updated currents).

Relative error =  3.78e-05
-----
Iteration: 14
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-5.16e+01', '5.53e+00', '1.45e+01', '6.66e+00', '-1.27e+01', '-7.07e+00', '-1.03e+01', '-7.53e+00', '-3.29e+00', '-4.39e+01', '3.07e-05']
Constraint losses = 8.22e-03
Relative update of tokamak psi (in plasma core): 4.14e-03
Handing off to forward solve (with updated currents).
Relative error =  2.43e-05
-----
Iteration: 15
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.

Change in coil currents (being controlled): ['3.50e+01', '9.61e+00', '-5.50e+01', '3.42e+01', '-5.66e+00', '3.65e+01', '2.07e+01', '1.52e+01', '-2.40e+01', '2.90e+01', '2.94e-05']
Constraint losses = 3.50e-03
Relative update of tokamak psi (in plasma core): 2.83e-03
Handing off to forward solve (with updated currents).
Relative error =  2.76e-05
-----
Iteration: 16
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-4.21e+01', '6.70e+00', '4.71e+00', '1.17e+01', '-1.25e+01', '-7.55e-01', '-6.37e+00', '-4.70e+00', '-7.31e+00', '-3.51e+01', '3.23e-05']
Constraint losses = 5.10e-03
Relative update of tokamak psi (in plasma core): 3.27e-03
Handing off to forward solve (with updated currents).

Relative error =  1.00e-05
-----
Iteration: 17
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['3.75e+01', '5.79e+00', '-4.57e+01', '2.51e+01', '-1.95e+00', '2.97e+01', '1.84e+01', '1.36e+01', '-1.75e+01', '2.99e+01', '6.05e-06']
Constraint losses = 4.34e-03
Relative update of tokamak psi (in plasma core): 2.91e-03
Handing off to forward solve (with updated currents).
Relative error =  1.40e-05
-----
Iteration: 18
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.35e+01', '5.71e+00', '3.08e+00', '9.99e+00', '-1.03e+01', '-5.41e-02', '-4.86e+00', '-3.58e+00', '-6.29e+00', '-2.76e+01', '7.91e-05']
Constraint losses = 4.42e-03
Relative update of tokamak psi (in plasma core): 2.57e-03
Handing off to forward solve (with updated currents).

Relative error =  6.18e-06
-----
Iteration: 19
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): ['-2.97e-01', '2.70e-01', '2.07e-01', '4.39e-01', '7.49e-01', '1.54e+00', '1.32e+00', '1.62e+00', '2.33e+00', '5.38e+00', '5.78e-05']
Constraint losses = 5.27e-03
Relative update of tokamak psi (in plasma core): 5.82e-04
Handing off to forward solve (with updated currents).
Relative error =  1.50e-07
-----
Inverse static solve SUCCESS. Tolerance 1.50e-07 (vs. requested 1e-06) reached in 20/100 iterations.

Notice how the strikepoint isoflux constraint is no longer satisfied.

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