Example: “Anamak” tokamak
Here, we generate a few different equilibria using the “Anamak” example tokamak from the FGE code paper.
The equilbirium\profile parameters are completely made up based on some equilbiria in the above paper - please experiment on your own and change them to more realistic values as you please!
Create the machine object
The Anamak machine is a useful toy tokamak because the geometry can be varied in a simple and controlled way. In particular, the number and placement of active coils, passive structures, and the vessel can all be adjusted using a small set of geometric parameters.
Following the parameterisation used in the reference above, the \(N_{\mathrm{coils}}\) active coils are placed on an ellipse centred at \((R_0, 0)\) in the poloidal \((R, Z)\) plane. Their coordinates are defined by
\[
R_{\mathrm{coil}}^{(i)} = R_0 + \hat{R}_{\mathrm{coil}} \cos(\theta_i),
\]
\[
Z_{\mathrm{coil}}^{(i)} = \kappa \, \hat{R}_{\mathrm{coil}} \sin(\theta_i),
\]
where
\(R_0\) is the major radius at which the ellipse is centred.
\(\hat{R}_{\mathrm{coil}}\) is the minor radius of the coils.
\(\bold{\theta}\) is a vector (length \(N_{\mathrm{coils}}\)) of poloidal angles \(\theta \in [0, 2\pi)\).
\(\kappa\) is the “elongation” of the ellipse.
Similarly, the \(N_{\mathrm{passives}}\) passive structures are placed on an ellipse at the same centre:
\[
R_{\mathrm{passives}}^{(i)} = R_0 + \hat{R}_{\mathrm{passive}} \cos(\theta_i),
\]
\[
Z_{\mathrm{passives}}^{(i)} = \kappa \, \hat{R}_{\mathrm{passive}} \sin(\theta_i),
\]
but with a potentially different minor radius \( \hat{R}_{\mathrm{passive}}\).
Finally the \(N_{\mathrm{limiter}}\) limiter/wall points are set up similarly:
\[
R_{\mathrm{limiter}}^{(i)} = R_0 + \hat{R}_{\mathrm{limiter}} \cos(\theta_i),
\]
\[
Z_{\mathrm{limiter}}^{(i)} = \kappa \, \hat{R}_{\mathrm{limiter}} \sin(\theta_i).
\]
We should note that given we are only doing static solves below, the passive structures are not really used here.
In the next cell, we can choose the parameters for building the Anamak.
import matplotlib.pyplot as plt
import numpy as np
# major radius of the machine
R0 = 1.0
# "elongation" of the machine (same for actives, passives, and limiter here)
kappa = 1.0
# for each of the coils, passives, and limiter, define:
# -> N = the number of each you desire
# -> R = the minor radius of its ellipse
# actives parameters
N_actives = 8
R_coil = 0.7
# passives parameters
N_passives = 50
R_passive = 0.5
# limiter parameters
N_limiter = 101
R_limiter = 0.45
The following functions are used to build the Anamak.
# builds an ellipse of (R,Z) points
def build_ellipse(
R0, # major radius
N, # number of points in [0,2\pi)
R_minor, # minor radius
kappa, # "elongation"
):
# poloidal locations
theta = np.arange(0, 2*np.pi, 2*np.pi/N)
# locations
r_locs = R0 + R_minor*np.cos(theta)
z_locs = R_minor*kappa*np.sin(theta)
return r_locs, z_locs
# builds active coils using (R,Z) locations of filaments
def build_actives(
r_locs,
z_locs,
):
# populate active coils dictionary with single filament coils (point-sources)
active_coils_data = {}
for i, (r, z) in enumerate(zip(r_locs, z_locs)):
active_coils_data[i] = {
"R": [r],
"Z": [z],
"dR": 0.05, # size only used for plotting in the simulations below
"dZ": 0.05, # size only used for plotting in the simulations below
"resistivity": 1.55e-8, # not used below
"polarity": 1, # default
"multiplier": 1, # single turn coils
}
return active_coils_data
# builds passive structures using (R,Z) locations of filaments
def build_passives(
r_locs,
z_locs,
):
# populate passive coils list with single filament passives (point-sources)
passive_coils_data = []
for r, z in zip(r_locs, z_locs):
passive_coils_data.append({
"R": r,
"Z": z,
"dR": 0.02, # size only used for plotting in the simulations below
"dZ": 0.02, # size only used for plotting in the simulations below
"resistivity": 5.5e-7 # resistivity of material (steel)
})
return passive_coils_data
# builds limiter geometry using (R,Z) locations
def build_limiter(
r_locs,
z_locs,
):
# populate limiter list
limiter_data = []
for r, z in zip(r_locs, z_locs):
limiter_data.append({"R": r, "Z": z})
return limiter_data
# build the active coils data
r_active_locs, z_active_locs = build_ellipse(
R0=R0,
N=N_actives,
R_minor=1.4*R_passive,
kappa=kappa,
)
active_coils_data = build_actives(
r_locs=r_active_locs,
z_locs=z_active_locs
)
# build the passive structures data
r_passive_locs, z_passive_locs = build_ellipse(
R0=R0,
N=N_passives,
R_minor=R_passive,
kappa=kappa,
)
passive_coils_data = build_passives(
r_locs=r_passive_locs,
z_locs=z_passive_locs
)
# build the limiter data
r_limiter_locs, z_limiter_locs = build_ellipse(
R0=R0,
N=N_limiter,
R_minor=R_limiter,
kappa=kappa,
)
limiter_data = build_limiter(
r_locs=r_limiter_locs,
z_locs=z_limiter_locs
)
# build the Anamak
from freegsnke import build_machine
tokamak = build_machine.tokamak(
active_coils_data=active_coils_data,
passive_coils_data=passive_coils_data,
limiter_data=limiter_data,
wall_data=limiter_data,
)
# plot the machine
fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)
ax1.grid(zorder=0, alpha=0.75)
ax1.set_aspect('equal')
tokamak.plot(axis=ax1,show=False)
ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)
ax1.set_xlabel("Major radius [m]")
ax1.set_ylabel("Height [m]")
Active coils --> built from user-provided data.
Passive structures --> built from user-provided data.
Limiter --> built from user-provided data.
Wall --> built from user-provided data.
Magnetic probes --> none provided.
Resistance (R) and inductance (M) matrices --> built using actives (and passives if present).
Tokamak built.
Text(0, 0.5, 'Height [m]')
Instantiate an equilibrium
from freegsnke import equilibrium_update
eq = equilibrium_update.Equilibrium(
tokamak=tokamak, # provide tokamak object
Rmin=0.1, Rmax=1.9, # radial range
Zmin=-0.8, Zmax=0.8, # 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)
)
Instantiate a profile object
We use the Lao profile object here, providing constant pressure and toroidal field profiles for simplicity.
from freegsnke.jtor_update import Lao85
profiles = Lao85(
eq=eq, # eq object
Ip=2e5, # plasma current
fvac=0.5*R0, # toroidal field * major radius
alpha=[58213.6], # linear profile
beta=[0.582136] # linear profile
)
Load the static nonlinear solver
from freegsnke import GSstaticsolver
GSStaticSolver = GSstaticsolver.NKGSsolver(eq)
Inverse solver constraints
In this first example, we’ll set up some simple constraints to generate an up/down symmetric double null plasma (with relatively high triangularity).
# null points constraints
Rx = 0.8 # X-point radius
Zx = 0.35 # X-point height
null_points = [[Rx, Rx], [Zx, -Zx]]
# isoflux constraints (including null points here)
Rout = 1.22 # outboard midplane radius
isoflux_set = np.array([
[
[Rx, Rx, Rout],
[Zx, -Zx, 0.0]
]]
)
# instantiate the freegsnke constrain object
from freegsnke.inverse import Inverse_optimizer
constrain = Inverse_optimizer(
null_points=null_points,
isoflux_set=isoflux_set,
weight_isoflux=1.0,
weight_nulls=0.8, # slightly lower weight on null points to get correct shape
)
The inverse solve
Solve using the inverse solver with forced up/down symmetry.
Setting force_up_down_symmetric=True projects the initial plasma flux onto its even-in-Z component and keeps the plasma response even in every forward solve, including those used to construct the full inverse Jacobian. The grid, machine, and requested constraints should also be up/down symmetric. This flag constrains the plasma response, not the coil-current optimisation: if exactly equal currents are required in separate upper and lower coils, those currents must additionally be linked or constrained by the machine/control model.
GSStaticSolver.inverse_solve(eq=eq,
profiles=profiles,
constrain=constrain,
target_relative_tolerance=1e-5,
target_relative_psit_update=1e-3,
verbose=True,
l2_reg=np.array([1e-12]*N_actives),
force_up_down_symmetric=True,
)
plasma_flux_parity_error = np.linalg.norm(eq.plasma_psi - eq.plasma_psi[:, ::-1]) / np.linalg.norm(eq.plasma_psi)
print(f"Relative plasma-flux parity error = {plasma_flux_parity_error:.3e}")
-----
Inverse static solve starting...
Successfully computed GS residual using initial plasma_psi and tokamak_psi guesses.
Initial relative error = 9.74e-01
-----
Iteration: 0
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['6.54e+05', '-5.57e+05', '3.92e+05', '1.55e+06', '-9.89e+05', '1.55e+06', '3.92e+05', '-5.57e+05']
Constraint losses = 5.03e-01
Relative update of tokamak psi (in plasma core): 1.73e+07
Handing off to forward solve (with updated currents).
Relative error = 9.34e-01
-----
Iteration: 1
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.40e+05', '1.06e+05', '-9.66e+04', '-3.20e+05', '2.07e+05', '-3.20e+05', '-9.66e+04', '1.06e+05']
Constraint losses = 1.01e-01
Relative update of tokamak psi (in plasma core): 4.65e-02
Handing off to forward solve (with updated currents).
Relative error = 8.81e-01
-----
Iteration: 2
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.19e+05', '8.23e+04', '-6.85e+04', '-2.54e+05', '1.63e+05', '-2.54e+05', '-6.85e+04', '8.23e+04']
Constraint losses = 7.57e-02
Relative update of tokamak psi (in plasma core): 5.67e-02
Handing off to forward solve (with updated currents).
Relative error = 8.06e-01
-----
Iteration: 3
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-8.48e+04', '5.92e+04', '-4.93e+04', '-1.82e+05', '1.17e+05', '-1.82e+05', '-4.93e+04', '5.92e+04']
Constraint losses = 5.45e-02
Relative update of tokamak psi (in plasma core): 7.94e-02
Handing off to forward solve (with updated currents).
Relative error = 7.03e-01
-----
Iteration: 4
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-6.00e+04', '4.25e+04', '-3.55e+04', '-1.30e+05', '8.32e+04', '-1.30e+05', '-3.55e+04', '4.25e+04']
Constraint losses = 3.91e-02
Relative update of tokamak psi (in plasma core): 5.57e-02
Handing off to forward solve (with updated currents).
Relative error = 5.61e-01
-----
Iteration: 5
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-4.20e+04', '3.04e+04', '-2.55e+04', '-9.21e+04', '5.91e+04', '-9.21e+04', '-2.55e+04', '3.04e+04']
Constraint losses = 2.81e-02
Relative update of tokamak psi (in plasma core): 5.02e-02
Handing off to forward solve (with updated currents).
Relative error = 3.64e-01
-----
Iteration: 6
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.86e+04', '2.15e+04', '-1.82e+04', '-6.43e+04', '4.13e+04', '-6.43e+04', '-1.82e+04', '2.15e+04']
Constraint losses = 2.00e-02
Relative update of tokamak psi (in plasma core): 4.28e-02
Handing off to forward solve (with updated currents).
Relative error = 8.86e-02
-----
Iteration: 7
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['1.66e+03', '1.10e+03', '-1.15e+03', '-7.48e+02', '5.70e+02', '-7.48e+02', '-1.15e+03', '1.10e+03']
Constraint losses = 1.40e-02
Relative update of tokamak psi (in plasma core): 2.28e-03
Handing off to forward solve (with updated currents).
Relative error = 3.43e-06
-----
Iteration: 8
Using full Jacobian (of constraints wrt coil currents) to optimsise currents.
- calculating derivatives for coil 1/8
- calculating derivatives for coil 2/8
- calculating derivatives for coil 3/8
- calculating derivatives for coil 4/8
- calculating derivatives for coil 5/8
- calculating derivatives for coil 6/8
- calculating derivatives for coil 7/8
- calculating derivatives for coil 8/8
Change in coil currents (being controlled): ['-1.81e+00', '-1.28e+00', '-2.15e-01', '-9.48e-02', '2.13e-01', '-9.48e-02', '-2.15e-01', '-1.16e+00']
Constraint losses = 1.50e-02
Relative update of tokamak psi (in plasma core): 5.27e-06
Handing off to forward solve (with updated currents).
Relative error = 8.18e-06
-----
Inverse static solve SUCCESS. Tolerance 8.18e-06 (vs. requested 1e-05) reached in 9/100 iterations.
Relative plasma-flux parity error = 0.000e+00
# refine with a forward solve (now the currents are known)
GSStaticSolver.solve(eq=eq,
profiles=profiles,
constrain=None,
target_relative_tolerance=1e-9,
verbose=False, # print output
)
Forward static solve SUCCESS. Tolerance 8.78e-10 (vs. requested 1.00e-09) reached in 8/100 iterations.
Using only a few constraints, we can see that the inverse solver has found the desired plasma. It is almost perfectly up/down symmetric and has high triangularity.
print(f"Triangularity = {eq.triangularity()}")
print(f"Avg. Z (Jtor) position = {eq.Zcurrent()} [m]")
Triangularity = 0.8328932260406068
Avg. Z (Jtor) position = 1.492261613633297e-07 [m]
# plot the equilibrium
fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), 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
eq.plot(axis=ax1,show=False) # plots the equilibrium
constrain.plot(axis=ax1, show=False) # plots the contraints
<Axes: xlabel='Major radius [m]', ylabel='Height [m]'>
Next, we’ll use the same machine description to make a negtive triangularity plasma (again up/down symmetric).
# null points constraints
Rx = 1.1 # X-point radius
Zx = 0.35 # X-point height
null_points = [[Rx, Rx], [Zx, -Zx]]
# isoflux constraints (including null points here)
Rin = 0.8 # inboard midplane radius
isoflux_set = np.array([
[
[Rx, Rx, Rin],
[Zx, -Zx, 0.0]
]]
)
# instantiate the freegsnke constrain object
from freegsnke.inverse import Inverse_optimizer
constrain = Inverse_optimizer(
null_points=null_points,
isoflux_set=isoflux_set,
weight_isoflux=1.0,
weight_nulls=0.8, # slightly lower weight on null points to get correct shape
)
# solve
GSStaticSolver.inverse_solve(eq=eq,
profiles=profiles,
constrain=constrain,
target_relative_tolerance=1e-5,
target_relative_psit_update=1e-3,
verbose=True,
l2_reg=np.array([1e-14]*N_actives),
force_up_down_symmetric=True,
)
-----
Inverse static solve starting...
Successfully computed GS residual using initial plasma_psi and tokamak_psi guesses.
Initial relative error = 1.63e-07
-----
Iteration: 0
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-8.22e+03', '9.93e+03', '7.02e+03', '-2.90e+03', '7.14e+02', '-2.90e+03', '7.02e+03', '9.93e+03']
Constraint losses = 2.13e-01
Relative update of tokamak psi (in plasma core): 1.99e-02
Handing off to forward solve (with updated currents).
Relative error = 4.01e-03
-----
Iteration: 1
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.14e+04', '1.03e+04', '7.12e+03', '-3.52e+03', '1.14e+03', '-3.52e+03', '7.12e+03', '1.03e+04']
Constraint losses = 1.93e-01
Relative update of tokamak psi (in plasma core): 1.98e-02
Handing off to forward solve (with updated currents).
Relative error = 1.81e-02
-----
Iteration: 2
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.19e+04', '1.05e+04', '6.93e+03', '-3.66e+03', '1.25e+03', '-3.66e+03', '6.93e+03', '1.05e+04']
Constraint losses = 1.83e-01
Relative update of tokamak psi (in plasma core): 1.96e-02
Handing off to forward solve (with updated currents).
Relative error = 3.09e-02
-----
Iteration: 3
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.24e+04', '1.08e+04', '6.76e+03', '-3.84e+03', '1.39e+03', '-3.84e+03', '6.76e+03', '1.08e+04']
Constraint losses = 1.73e-01
Relative update of tokamak psi (in plasma core): 1.94e-02
Handing off to forward solve (with updated currents).
Relative error = 4.22e-02
-----
Iteration: 4
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.33e+04', '1.14e+04', '6.70e+03', '-4.13e+03', '1.58e+03', '-4.13e+03', '6.70e+03', '1.14e+04']
Constraint losses = 1.64e-01
Relative update of tokamak psi (in plasma core): 1.92e-02
Handing off to forward solve (with updated currents).
Relative error = 5.15e-02
-----
Iteration: 5
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.85e+04', '1.42e+04', '7.41e+03', '-5.56e+03', '2.39e+03', '-5.56e+03', '7.41e+03', '1.42e+04']
Constraint losses = 1.58e-01
Relative update of tokamak psi (in plasma core): 1.91e-02
Handing off to forward solve (with updated currents).
Relative error = 5.56e-02
-----
Iteration: 6
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.11e+04', '1.59e+04', '7.46e+03', '-6.38e+03', '2.89e+03', '-6.38e+03', '7.46e+03', '1.59e+04']
Constraint losses = 1.47e-01
Relative update of tokamak psi (in plasma core): 1.88e-02
Handing off to forward solve (with updated currents).
Relative error = 6.49e-02
-----
Iteration: 7
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.53e+04', '1.86e+04', '7.63e+03', '-7.71e+03', '3.70e+03', '-7.71e+03', '7.63e+03', '1.86e+04']
Constraint losses = 1.35e-01
Relative update of tokamak psi (in plasma core): 1.86e-02
Handing off to forward solve (with updated currents).
Relative error = 5.37e-02
-----
Iteration: 8
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-5.95e+04', '3.52e+04', '1.47e+04', '-1.65e+04', '8.37e+03', '-1.65e+04', '1.47e+04', '3.52e+04']
Constraint losses = 1.25e-01
Relative update of tokamak psi (in plasma core): 1.83e-02
Handing off to forward solve (with updated currents).
Relative error = 2.27e-03
-----
Iteration: 9
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.98e+04', '1.73e+04', '4.86e+03', '-1.01e+04', '5.90e+03', '-1.01e+04', '4.86e+03', '1.73e+04']
Constraint losses = 8.88e-02
Relative update of tokamak psi (in plasma core): 1.80e-02
Handing off to forward solve (with updated currents).
Relative error = 2.92e-05
-----
Iteration: 10
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.37e+05', '7.54e+04', '1.34e+04', '-3.89e+04', '2.28e+04', '-3.89e+04', '1.34e+04', '7.54e+04']
Constraint losses = 9.04e-02
Relative update of tokamak psi (in plasma core): 1.77e-02
Handing off to forward solve (with updated currents).
Relative error = 5.28e-04
-----
Iteration: 11
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.05e+05', '7.25e+04', '-2.30e+04', '-3.73e+04', '2.62e+04', '-3.73e+04', '-2.30e+04', '7.25e+04']
Constraint losses = 6.61e-02
Relative update of tokamak psi (in plasma core): 1.73e-02
Handing off to forward solve (with updated currents).
Relative error = 1.79e-05
-----
Iteration: 12
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.99e+04', '2.41e+04', '-1.43e+04', '-1.40e+04', '1.07e+04', '-1.40e+04', '-1.43e+04', '2.41e+04']
Constraint losses = 4.68e-02
Relative update of tokamak psi (in plasma core): 1.70e-02
Handing off to forward solve (with updated currents).
Relative error = 1.54e-04
-----
Iteration: 13
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-4.28e+04', '2.91e+04', '-1.91e+04', '-1.63e+04', '1.27e+04', '-1.63e+04', '-1.91e+04', '2.91e+04']
Constraint losses = 4.39e-02
Relative update of tokamak psi (in plasma core): 1.66e-02
Handing off to forward solve (with updated currents).
Relative error = 6.85e-05
-----
Iteration: 14
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-4.05e+04', '2.79e+04', '-1.95e+04', '-1.57e+04', '1.23e+04', '-1.57e+04', '-1.95e+04', '2.79e+04']
Constraint losses = 3.91e-02
Relative update of tokamak psi (in plasma core): 1.62e-02
Handing off to forward solve (with updated currents).
Relative error = 5.41e-05
-----
Iteration: 15
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.96e+04', '2.81e+04', '-2.08e+04', '-1.57e+04', '1.25e+04', '-1.57e+04', '-2.08e+04', '2.81e+04']
Constraint losses = 3.42e-02
Relative update of tokamak psi (in plasma core): 1.58e-02
Handing off to forward solve (with updated currents).
Relative error = 1.82e-05
-----
Iteration: 16
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.66e+04', '2.55e+04', '-1.99e+04', '-1.45e+04', '1.16e+04', '-1.45e+04', '-1.99e+04', '2.55e+04']
Constraint losses = 2.83e-02
Relative update of tokamak psi (in plasma core): 1.54e-02
Handing off to forward solve (with updated currents).
Relative error = 8.47e-05
-----
Iteration: 17
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.57e+04', '2.58e+04', '-2.12e+04', '-1.46e+04', '1.18e+04', '-1.46e+04', '-2.12e+04', '2.58e+04']
Constraint losses = 2.34e-02
Relative update of tokamak psi (in plasma core): 1.50e-02
Handing off to forward solve (with updated currents).
Relative error = 3.89e-03
-----
Iteration: 18
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.06e+04', '2.08e+04', '-1.73e+04', '-1.21e+04', '9.80e+03', '-1.21e+04', '-1.73e+04', '2.08e+04']
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 = 1.24e-04
-----
Iteration: 19
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.62e+04', '1.47e+04', '-1.47e+04', '-7.95e+03', '6.79e+03', '-7.95e+03', '-1.47e+04', '1.47e+04']
Constraint losses = 1.40e-02
Relative update of tokamak psi (in plasma core): 7.17e-03
Handing off to forward solve (with updated currents).
Relative error = 3.59e-05
-----
Iteration: 20
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.53e+04', '9.84e+03', '-8.76e+03', '-5.91e+03', '4.86e+03', '-5.91e+03', '-8.76e+03', '9.84e+03']
Constraint losses = 8.43e-03
Relative update of tokamak psi (in plasma core): 6.92e-03
Handing off to forward solve (with updated currents).
Relative error = 3.79e-05
-----
Iteration: 21
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-5.19e+03', '6.63e+03', '-7.58e+03', '-3.35e+03', '3.00e+03', '-3.35e+03', '-7.58e+03', '6.63e+03']
Constraint losses = 7.32e-03
Relative update of tokamak psi (in plasma core): 2.55e-03
Handing off to forward solve (with updated currents).
Relative error = 6.40e-06
-----
Iteration: 22
Using full Jacobian (of constraints wrt coil currents) to optimsise currents.
- calculating derivatives for coil 1/8
- calculating derivatives for coil 2/8
- calculating derivatives for coil 3/8
- calculating derivatives for coil 4/8
- calculating derivatives for coil 5/8
- calculating derivatives for coil 6/8
- calculating derivatives for coil 7/8
- calculating derivatives for coil 8/8
Change in coil currents (being controlled): ['-1.52e-01', '-9.54e-03', '-7.53e-03', '-2.55e-03', '3.41e-03', '-6.87e-03', '-7.53e-03', '-9.54e-03']
Constraint losses = 3.83e-03
Relative update of tokamak psi (in plasma core): 1.12e-07
Handing off to forward solve (with updated currents).
Relative error = 6.28e-06
-----
Inverse static solve SUCCESS. Tolerance 6.28e-06 (vs. requested 1e-05) reached in 23/100 iterations.
# refine with a forward solve (now the currents are known)
GSStaticSolver.solve(eq=eq,
profiles=profiles,
constrain=None,
target_relative_tolerance=1e-9,
verbose=False, # print output
)
Forward static solve SUCCESS. Tolerance 4.81e-10 (vs. requested 1.00e-09) reached in 6/100 iterations.
print(f"Triangularity = {eq.triangularity()}")
print(f"Avg. Z (Jtor) position = {eq.Zcurrent()} [m]")
Triangularity = -0.784878022381521
Avg. Z (Jtor) position = -5.21864031802719e-10 [m]
# plot the equilibrium
fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), 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
eq.plot(axis=ax1,show=False) # plots the equilibrium
constrain.plot(axis=ax1, show=False) # plots the contraints
<Axes: xlabel='Major radius [m]', ylabel='Height [m]'>
Now lets look at changing the shape of the Anamak.
We’ll now increase the number of active coils and increase the “elongation” of the device.
# "elongation" of the machine (same for actives, passives, and limiter here)
kappa = 1.8
# for each of the coils, passives, and limiter, define:
# -> N = the number of each you desire
# -> R = the minor radius of its ellipse
# actives parameters
N_actives = 24
R_coil = 0.7
# build the active coils data
r_active_locs, z_active_locs = build_ellipse(
R0=R0,
N=N_actives,
R_minor=1.4*R_passive,
kappa=kappa,
)
active_coils_data = build_actives(
r_locs=r_active_locs,
z_locs=z_active_locs
)
# build the passive structures data
r_passive_locs, z_passive_locs = build_ellipse(
R0=R0,
N=N_passives,
R_minor=R_passive,
kappa=kappa,
)
passive_coils_data = build_passives(
r_locs=r_passive_locs,
z_locs=z_passive_locs
)
# build the limiter data
r_limiter_locs, z_limiter_locs = build_ellipse(
R0=R0,
N=N_limiter,
R_minor=R_limiter,
kappa=kappa,
)
limiter_data = build_limiter(
r_locs=r_limiter_locs,
z_locs=z_limiter_locs
)
# build the Anamak
from freegsnke import build_machine
tokamak = build_machine.tokamak(
active_coils_data=active_coils_data,
passive_coils_data=passive_coils_data,
limiter_data=limiter_data,
wall_data=limiter_data,
)
# plot the machine
fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)
ax1.grid(zorder=0, alpha=0.75)
ax1.set_aspect('equal')
tokamak.plot(axis=ax1,show=False)
ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)
ax1.set_xlabel("Major radius [m]")
ax1.set_ylabel("Height [m]")
Active coils --> built from user-provided data.
Passive structures --> built from user-provided data.
Limiter --> built from user-provided data.
Wall --> built from user-provided data.
Magnetic probes --> none provided.
Resistance (R) and inductance (M) matrices --> built using actives (and passives if present).
Tokamak built.
Text(0, 0.5, 'Height [m]')
Re-initialise the equilibrium and profiles objects: here we’ve increased the vertical size of the computational grid and increased the plasma current.
from freegsnke import equilibrium_update
eq = equilibrium_update.Equilibrium(
tokamak=tokamak, # provide tokamak object
Rmin=0.1, Rmax=1.9, # radial range
Zmin=-1.5, Zmax=1.5, # 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)
)
from freegsnke.jtor_update import Lao85
profiles = Lao85(
eq=eq, # eq object
Ip=5e5, # plasma current
fvac=0.5*R0, # toroidal field * major radius
alpha=[58213.6], # linear profile
beta=[0.582136] # linear profile
)
from freegsnke import GSstaticsolver
GSStaticSolver = GSstaticsolver.NKGSsolver(eq)
Let’s try to create a simple up/down symmetric double null plasma in this new machine.
# null points constraints
Rx = 0.8 # X-point radius
Zx = 0.65 # X-point height
null_points = [[Rx, Rx], [Zx, -Zx]]
# isoflux constraints (including null points here)
Rout = 1.3 # outboard midplane radius
isoflux_set = np.array([
[
[Rx, Rx, Rout],
[Zx, -Zx, 0.0]
]]
)
# instantiate the freegsnke constrain object
from freegsnke.inverse import Inverse_optimizer
constrain = Inverse_optimizer(
null_points=null_points,
isoflux_set=isoflux_set,
weight_isoflux=1.0,
weight_nulls=0.7, # slightly lower weight on null points to get correct shape
)
# solve
GSStaticSolver.inverse_solve(eq=eq,
profiles=profiles,
constrain=constrain,
target_relative_tolerance=1e-5,
target_relative_psit_update=1e-3,
verbose=True,
l2_reg=np.array([1e-12]*N_actives),
force_up_down_symmetric=True,
)
-----
Inverse static solve starting...
Successfully computed GS residual using initial plasma_psi and tokamak_psi guesses.
Initial relative error = 9.81e-01
-----
Iteration: 0
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['2.39e+05', '4.33e+04', '-1.81e+05', '-2.09e+05', '-1.16e+05', '-8.78e+03', '8.27e+04', '1.76e+05', '2.90e+05', '3.26e+05', '2.54e+04', '-1.24e+05', '-1.20e+05', '-1.24e+05', '2.54e+04', '3.26e+05', '2.90e+05', '1.76e+05', '8.27e+04', '-8.78e+03', '-1.16e+05', '-2.09e+05', '-1.81e+05', '4.33e+04']
Constraint losses = 3.19e-01
Relative update of tokamak psi (in plasma core): 7.40e+06
Handing off to forward solve (with updated currents).
Relative error = 9.20e-01
-----
Iteration: 1
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-6.98e+04', '-2.86e+04', '2.23e+04', '3.61e+04', '2.50e+04', '8.65e+03', '-6.31e+03', '-2.20e+04', '-4.20e+04', '-5.15e+04', '-5.59e+03', '1.91e+04', '1.86e+04', '1.91e+04', '-5.59e+03', '-5.15e+04', '-4.20e+04', '-2.20e+04', '-6.31e+03', '8.65e+03', '2.50e+04', '3.61e+04', '2.23e+04', '-2.86e+04']
Constraint losses = 3.33e-02
Relative update of tokamak psi (in plasma core): 4.24e-02
Handing off to forward solve (with updated currents).
Relative error = 8.19e-01
-----
Iteration: 2
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-5.34e+04', '-2.16e+04', '1.74e+04', '2.79e+04', '1.92e+04', '6.49e+03', '-5.09e+03', '-1.72e+04', '-3.27e+04', '-4.00e+04', '-4.29e+03', '1.48e+04', '1.44e+04', '1.48e+04', '-4.29e+03', '-4.00e+04', '-3.27e+04', '-1.72e+04', '-5.09e+03', '6.49e+03', '1.92e+04', '2.79e+04', '1.74e+04', '-2.16e+04']
Constraint losses = 2.57e-02
Relative update of tokamak psi (in plasma core): 3.97e-02
Handing off to forward solve (with updated currents).
Relative error = 6.92e-01
-----
Iteration: 3
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.92e+04', '-1.57e+04', '1.31e+04', '2.07e+04', '1.42e+04', '4.73e+03', '-3.89e+03', '-1.29e+04', '-2.44e+04', '-2.98e+04', '-3.17e+03', '1.11e+04', '1.08e+04', '1.11e+04', '-3.17e+03', '-2.98e+04', '-2.44e+04', '-1.29e+04', '-3.89e+03', '4.73e+03', '1.42e+04', '2.07e+04', '1.31e+04', '-1.57e+04']
Constraint losses = 1.91e-02
Relative update of tokamak psi (in plasma core): 2.78e-02
Handing off to forward solve (with updated currents).
Relative error = 5.24e-01
-----
Iteration: 4
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.85e+04', '-1.13e+04', '9.80e+03', '1.53e+04', '1.04e+04', '3.41e+03', '-2.97e+03', '-9.61e+03', '-1.81e+04', '-2.20e+04', '-2.33e+03', '8.18e+03', '7.96e+03', '8.18e+03', '-2.33e+03', '-2.20e+04', '-1.81e+04', '-9.61e+03', '-2.97e+03', '3.41e+03', '1.04e+04', '1.53e+04', '9.80e+03', '-1.13e+04']
Constraint losses = 1.41e-02
Relative update of tokamak psi (in plasma core): 1.92e-02
Handing off to forward solve (with updated currents).
Relative error = 2.99e-01
-----
Iteration: 5
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.02e+04', '-7.82e+03', '7.29e+03', '1.11e+04', '7.51e+03', '2.39e+03', '-2.25e+03', '-7.08e+03', '-1.32e+04', '-1.60e+04', '-1.68e+03', '5.97e+03', '5.80e+03', '5.97e+03', '-1.68e+03', '-1.60e+04', '-1.32e+04', '-7.08e+03', '-2.25e+03', '2.39e+03', '7.51e+03', '1.11e+04', '7.29e+03', '-7.82e+03']
Constraint losses = 1.03e-02
Relative update of tokamak psi (in plasma core): 1.26e-02
Handing off to forward solve (with updated currents).
Relative error = 6.50e-03
-----
Iteration: 6
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.63e+02', '-1.02e+02', '2.15e+02', '2.89e+02', '2.00e+02', '7.93e+01', '-2.91e+01', '-1.41e+02', '-2.82e+02', '-3.54e+02', '-4.42e+01', '1.29e+02', '1.27e+02', '1.29e+02', '-4.42e+01', '-3.54e+02', '-2.82e+02', '-1.41e+02', '-2.91e+01', '7.93e+01', '2.00e+02', '2.89e+02', '2.15e+02', '-1.02e+02']
Constraint losses = 7.44e-03
Relative update of tokamak psi (in plasma core): 1.04e-04
Handing off to forward solve (with updated currents).
Relative error = 1.91e-06
-----
Inverse static solve SUCCESS. Tolerance 1.91e-06 (vs. requested 1e-05) reached in 7/100 iterations.
# refine with a forward solve (now the currents are known)
GSStaticSolver.solve(eq=eq,
profiles=profiles,
constrain=None,
target_relative_tolerance=1e-9,
verbose=False, # print output
)
Forward static solve SUCCESS. Tolerance 6.45e-10 (vs. requested 1.00e-09) reached in 5/100 iterations.
# plot the equilibrium
fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), 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
eq.plot(axis=ax1,show=False) # plots the equilibrium
constrain.plot(axis=ax1, show=False) # plots the contraints
<Axes: xlabel='Major radius [m]', ylabel='Height [m]'>
Now let us try to use the snowflake (i.e. second order null) constraints to make something a bit more exotic.
# null points constraints
Rx = 0.8 # X-point radius
Zx = 0.65 # X-point height
null_points_2nd_order = [[Rx, Rx], [Zx, -Zx]]
# isoflux constraints (including null points here, not really needed but should be fine)
isoflux_set = np.array([
[
[Rx, Rx],
[Zx, -Zx]
]]
)
# instantiate the freegsnke constrain object
from freegsnke.inverse import Inverse_optimizer
constrain = Inverse_optimizer(
null_points_2nd_order=null_points_2nd_order,
# isoflux_set=isoflux_set,
)
# solve
GSStaticSolver.inverse_solve(eq=eq,
profiles=profiles,
constrain=constrain,
target_relative_tolerance=1e-5,
target_relative_psit_update=1e-3,
verbose=True,
Picard_handover=1e-5,
l2_reg=np.array([1e-12]*N_actives),
force_up_down_symmetric=True,
)
-----
Inverse static solve starting...
Successfully computed GS residual using initial plasma_psi and tokamak_psi guesses.
Initial relative error = 6.60e-10
-----
Iteration: 0
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.07e+04', '-1.26e+04', '7.55e+03', '2.24e+04', '1.86e+04', '4.30e+03', '-1.18e+04', '-3.01e+04', '-5.04e+04', '-2.64e+04', '3.43e+04', '3.17e+03', '-6.45e+03', '3.17e+03', '3.43e+04', '-2.64e+04', '-5.04e+04', '-3.01e+04', '-1.18e+04', '4.30e+03', '1.86e+04', '2.24e+04', '7.55e+03', '-1.26e+04']
Constraint losses = 1.38e+00
Relative update of tokamak psi (in plasma core): 1.99e-02
Handing off to forward solve (with updated currents).
Relative error = 1.27e-02
-----
Iteration: 1
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.07e+04', '-1.76e+04', '1.44e+04', '3.65e+04', '2.86e+04', '5.05e+03', '-2.03e+04', '-4.82e+04', '-7.69e+04', '-3.22e+04', '5.61e+04', '1.83e+03', '-1.29e+04', '1.83e+03', '5.61e+04', '-3.22e+04', '-7.69e+04', '-4.82e+04', '-2.03e+04', '5.05e+03', '2.86e+04', '3.65e+04', '1.44e+04', '-1.76e+04']
Constraint losses = 1.20e+00
Relative update of tokamak psi (in plasma core): 1.98e-02
Handing off to forward solve (with updated currents).
Relative error = 1.45e-02
-----
Iteration: 2
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-5.34e+04', '-2.69e+04', '3.52e+04', '7.41e+04', '5.40e+04', '6.63e+03', '-4.15e+04', '-9.15e+04', '-1.36e+05', '-3.02e+04', '1.12e+05', '-7.27e+03', '-3.36e+04', '-7.27e+03', '1.12e+05', '-3.02e+04', '-1.36e+05', '-9.15e+04', '-4.15e+04', '6.63e+03', '5.40e+04', '7.41e+04', '3.52e+04', '-2.69e+04']
Constraint losses = 8.37e-01
Relative update of tokamak psi (in plasma core): 5.68e-03
Handing off to forward solve (with updated currents).
Relative error = 2.79e-03
-----
Iteration: 3
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.67e+04', '-1.65e+04', '7.58e+03', '2.51e+04', '2.34e+04', '1.22e+04', '4.14e+02', '-1.09e+04', '-1.77e+04', '1.75e+04', '3.32e+04', '-1.26e+04', '-1.89e+04', '-1.26e+04', '3.32e+04', '1.75e+04', '-1.77e+04', '-1.09e+04', '4.14e+02', '1.22e+04', '2.34e+04', '2.51e+04', '7.58e+03', '-1.65e+04']
Constraint losses = 2.69e-01
Relative update of tokamak psi (in plasma core): 1.60e-02
Handing off to forward solve (with updated currents).
Relative error = 1.02e-03
-----
Iteration: 4
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.93e+04', '-1.95e+04', '3.95e+03', '2.17e+04', '2.15e+04', '1.22e+04', '2.13e+03', '-7.43e+03', '-1.23e+04', '2.10e+04', '3.10e+04', '-1.35e+04', '-1.91e+04', '-1.35e+04', '3.10e+04', '2.10e+04', '-1.23e+04', '-7.43e+03', '2.13e+03', '1.22e+04', '2.15e+04', '2.17e+04', '3.95e+03', '-1.95e+04']
Constraint losses = 2.32e-01
Relative update of tokamak psi (in plasma core): 9.96e-03
Handing off to forward solve (with updated currents).
Relative error = 1.75e-03
-----
Iteration: 5
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.39e+04', '-1.58e+04', '3.30e+03', '1.71e+04', '1.64e+04', '9.02e+03', '1.45e+03', '-5.15e+03', '-6.47e+03', '2.41e+04', '2.56e+04', '-1.32e+04', '-1.73e+04', '-1.32e+04', '2.56e+04', '2.41e+04', '-6.47e+03', '-5.15e+03', '1.45e+03', '9.02e+03', '1.64e+04', '1.71e+04', '3.30e+03', '-1.58e+04']
Constraint losses = 1.99e-01
Relative update of tokamak psi (in plasma core): 8.20e-03
Handing off to forward solve (with updated currents).
Relative error = 1.23e-03
-----
Iteration: 6
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.04e+04', '-1.35e+04', '2.70e+03', '1.43e+04', '1.37e+04', '7.59e+03', '1.43e+03', '-3.75e+03', '-4.06e+03', '2.26e+04', '2.16e+04', '-1.19e+04', '-1.52e+04', '-1.19e+04', '2.16e+04', '2.26e+04', '-4.06e+03', '-3.75e+03', '1.43e+03', '7.59e+03', '1.37e+04', '1.43e+04', '2.70e+03', '-1.35e+04']
Constraint losses = 1.74e-01
Relative update of tokamak psi (in plasma core): 6.96e-03
Handing off to forward solve (with updated currents).
Relative error = 1.00e-03
-----
Iteration: 7
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.78e+04', '-1.18e+04', '2.18e+03', '1.23e+04', '1.19e+04', '6.66e+03', '1.44e+03', '-2.93e+03', '-3.07e+03', '1.98e+04', '1.86e+04', '-1.04e+04', '-1.32e+04', '-1.04e+04', '1.86e+04', '1.98e+04', '-3.07e+03', '-2.93e+03', '1.44e+03', '6.66e+03', '1.19e+04', '1.23e+04', '2.18e+03', '-1.18e+04']
Constraint losses = 1.51e-01
Relative update of tokamak psi (in plasma core): 5.62e-03
Handing off to forward solve (with updated currents).
Relative error = 5.11e-04
-----
Iteration: 8
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.55e+04', '-1.05e+04', '1.39e+03', '1.02e+04', '1.02e+04', '6.10e+03', '1.88e+03', '-1.68e+03', '-1.82e+03', '1.67e+04', '1.53e+04', '-8.81e+03', '-1.11e+04', '-8.81e+03', '1.53e+04', '1.67e+04', '-1.82e+03', '-1.68e+03', '1.88e+03', '6.10e+03', '1.02e+04', '1.02e+04', '1.39e+03', '-1.05e+04']
Constraint losses = 1.29e-01
Relative update of tokamak psi (in plasma core): 4.38e-03
Handing off to forward solve (with updated currents).
Relative error = 6.80e-04
-----
Iteration: 9
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.27e+04', '-8.61e+03', '1.03e+03', '8.18e+03', '8.31e+03', '5.14e+03', '1.89e+03', '-7.69e+02', '-5.32e+02', '1.46e+04', '1.22e+04', '-7.50e+03', '-9.26e+03', '-7.50e+03', '1.22e+04', '1.46e+04', '-5.32e+02', '-7.69e+02', '1.89e+03', '5.14e+03', '8.31e+03', '8.18e+03', '1.03e+03', '-8.61e+03']
Constraint losses = 1.07e-01
Relative update of tokamak psi (in plasma core): 3.59e-03
Handing off to forward solve (with updated currents).
Relative error = 7.99e-04
-----
Iteration: 10
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-9.87e+03', '-6.82e+03', '4.03e+02', '5.93e+03', '6.41e+03', '4.44e+03', '2.38e+03', '8.10e+02', '1.48e+03', '1.24e+04', '8.53e+03', '-6.09e+03', '-7.20e+03', '-6.09e+03', '8.53e+03', '1.24e+04', '1.48e+03', '8.10e+02', '2.38e+03', '4.44e+03', '6.41e+03', '5.93e+03', '4.03e+02', '-6.82e+03']
Constraint losses = 8.51e-02
Relative update of tokamak psi (in plasma core): 3.07e-03
Handing off to forward solve (with updated currents).
Relative error = 5.28e-04
-----
Iteration: 11
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-8.16e+03', '-5.85e+03', '-2.35e+02', '4.36e+03', '5.21e+03', '4.03e+03', '2.62e+03', '1.47e+03', '1.76e+03', '8.99e+03', '6.05e+03', '-4.53e+03', '-5.37e+03', '-4.53e+03', '6.05e+03', '8.99e+03', '1.76e+03', '1.47e+03', '2.62e+03', '4.03e+03', '5.21e+03', '4.36e+03', '-2.35e+02', '-5.85e+03']
Constraint losses = 6.64e-02
Relative update of tokamak psi (in plasma core): 2.04e-03
Handing off to forward solve (with updated currents).
Relative error = 6.78e-04
-----
Iteration: 12
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-5.13e+03', '-3.91e+03', '-8.93e+02', '1.90e+03', '3.08e+03', '3.25e+03', '3.21e+03', '3.39e+03', '4.46e+03', '7.43e+03', '2.10e+03', '-3.31e+03', '-3.37e+03', '-3.31e+03', '2.10e+03', '7.43e+03', '4.46e+03', '3.39e+03', '3.21e+03', '3.25e+03', '3.08e+03', '1.90e+03', '-8.93e+02', '-3.91e+03']
Constraint losses = 5.05e-02
Relative update of tokamak psi (in plasma core): 1.80e-03
Handing off to forward solve (with updated currents).
Relative error = 4.85e-04
-----
Iteration: 13
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.22e+03', '-2.78e+03', '-1.49e+03', '2.47e+02', '1.78e+03', '2.79e+03', '3.50e+03', '4.24e+03', '5.14e+03', '4.45e+03', '-5.61e+02', '-1.85e+03', '-1.55e+03', '-1.85e+03', '-5.61e+02', '4.45e+03', '5.14e+03', '4.24e+03', '3.50e+03', '2.79e+03', '1.78e+03', '2.47e+02', '-1.49e+03', '-2.78e+03']
Constraint losses = 4.00e-02
Relative update of tokamak psi (in plasma core): 1.11e-03
Handing off to forward solve (with updated currents).
Relative error = 3.40e-04
-----
Iteration: 14
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.05e+03', '-1.22e+03', '-1.46e+03', '-9.96e+02', '3.15e+02', '1.74e+03', '2.97e+03', '4.22e+03', '5.43e+03', '3.08e+03', '-2.27e+03', '-9.72e+02', '-3.57e+02', '-9.72e+02', '-2.27e+03', '3.08e+03', '5.43e+03', '4.22e+03', '2.97e+03', '1.74e+03', '3.15e+02', '-9.96e+02', '-1.46e+03', '-1.22e+03']
Constraint losses = 3.43e-02
Relative update of tokamak psi (in plasma core): 8.28e-04
Handing off to forward solve (with updated currents).
Relative error = 2.58e-04
-----
Iteration: 15
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['8.85e+01', '-3.44e+02', '-1.24e+03', '-1.41e+03', '-3.90e+02', '9.75e+02', '2.20e+03', '3.41e+03', '4.39e+03', '1.48e+03', '-2.69e+03', '-2.53e+02', '3.71e+02', '-2.53e+02', '-2.69e+03', '1.48e+03', '4.39e+03', '3.41e+03', '2.20e+03', '9.75e+02', '-3.90e+02', '-1.41e+03', '-1.24e+03', '-3.44e+02']
Constraint losses = 2.72e-02
Relative update of tokamak psi (in plasma core): 4.22e-04
Handing off to forward solve (with updated currents).
Relative error = 8.29e-06
-----
Inverse static solve SUCCESS. Tolerance 8.29e-06 (vs. requested 1e-05) reached in 16/100 iterations.
# plot the equilibrium
fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), 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
eq.plot(axis=ax1,show=False) # plots the equilibrium
constrain.plot(axis=ax1, show=False) # plots the contraints
<Axes: xlabel='Major radius [m]', ylabel='Height [m]'>
And now perhaps a single null plasma - note we change some solver settings here too.
# null points constraints
Rx = 0.8 # X-point radius
Zx = 0.6 # X-point height
null_points = [[Rx], [-Zx]]
# isoflux constraints (including null points here)
Rout = 1.32 # outboard midplane radius
isoflux_set = np.array([
[
[Rx, Rout, 0.85, 1.17, 0.65],
[-Zx, 0.0, 0.57, 0.45, 0.0]
]]
)
# instantiate the freegsnke constrain object
from freegsnke.inverse import Inverse_optimizer
constrain = Inverse_optimizer(
null_points=null_points,
isoflux_set=isoflux_set,
weight_isoflux=1.0,
weight_nulls=0.7, # slightly lower weight on null points to get correct shape
)
# solve
eq.set_plasma_psi(eq.create_psi_plasma_default())
eq.tokamak.set_all_coil_currents(current_vec=np.zeros(N_actives))
GSStaticSolver.inverse_solve(eq=eq,
profiles=profiles,
constrain=constrain,
target_relative_tolerance=1e-5,
target_relative_psit_update=1e-3,
verbose=True,
l2_reg=np.array([1e-14]*N_actives),
Picard_handover=1e-5,
full_jacobian_handover=[2e-2,1e-2],
force_up_down_symmetric=False,
)
-----
Inverse static solve starting...
Successfully computed GS residual using initial plasma_psi and tokamak_psi guesses.
Initial relative error = 9.81e-01
-----
Iteration: 0
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['2.23e+05', '8.98e+04', '-7.68e+04', '-7.78e+04', '-3.67e+04', '-7.18e+03', '1.18e+04', '2.71e+04', '4.32e+04', '5.76e+04', '5.34e+04', '3.09e+04', '-2.67e+04', '-1.86e+05', '-7.90e+04', '2.93e+05', '3.01e+05', '2.08e+05', '1.21e+05', '2.69e+04', '-1.00e+05', '-2.44e+05', '-2.78e+05', '-5.05e+04']
Constraint losses = 2.19e-01
Relative update of tokamak psi (in plasma core): 7.54e+06
Handing off to forward solve (with updated currents).
Relative error = 9.45e-01
-----
Iteration: 1
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-7.11e+04', '-5.10e+04', '-7.42e+03', '9.46e+03', '1.24e+04', '1.16e+04', '1.03e+04', '8.90e+03', '7.04e+03', '2.62e+03', '-9.60e+03', '-3.24e+04', '-3.41e+04', '5.09e+04', '7.65e+04', '-2.91e+04', '-5.51e+04', '-4.23e+04', '-2.56e+04', '-5.12e+03', '2.44e+04', '5.97e+04', '6.78e+04', '7.14e+03']
Constraint losses = 4.22e-02
Relative update of tokamak psi (in plasma core): 3.30e-02
Handing off to forward solve (with updated currents).
Relative error = 8.46e-01
-----
Iteration: 2
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-5.27e+04', '-3.74e+04', '-5.04e+03', '7.54e+03', '9.73e+03', '9.16e+03', '8.13e+03', '7.10e+03', '5.75e+03', '2.54e+03', '-6.43e+03', '-2.31e+04', '-2.46e+04', '3.62e+04', '5.44e+04', '-2.12e+04', '-3.96e+04', '-3.02e+04', '-1.82e+04', '-3.46e+03', '1.77e+04', '4.28e+04', '4.82e+04', '3.84e+03']
Constraint losses = 3.11e-02
Relative update of tokamak psi (in plasma core): 3.70e-02
Handing off to forward solve (with updated currents).
Relative error = 7.15e-01
-----
Iteration: 3
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-3.93e+04', '-2.72e+04', '-3.16e+03', '6.13e+03', '7.69e+03', '7.20e+03', '6.39e+03', '5.61e+03', '4.61e+03', '2.30e+03', '-4.13e+03', '-1.60e+04', '-1.73e+04', '2.53e+04', '3.75e+04', '-1.56e+04', '-2.81e+04', '-2.12e+04', '-1.26e+04', '-2.15e+03', '1.28e+04', '3.04e+04', '3.35e+04', '1.22e+03']
Constraint losses = 2.32e-02
Relative update of tokamak psi (in plasma core): 2.83e-02
Handing off to forward solve (with updated currents).
Relative error = 5.47e-01
-----
Iteration: 4
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.92e+04', '-1.95e+04', '-1.69e+03', '5.08e+03', '6.13e+03', '5.69e+03', '5.04e+03', '4.43e+03', '3.71e+03', '2.08e+03', '-2.39e+03', '-1.06e+04', '-1.15e+04', '1.69e+04', '2.45e+04', '-1.15e+04', '-1.94e+04', '-1.44e+04', '-8.35e+03', '-1.12e+03', '9.12e+03', '2.09e+04', '2.24e+04', '-8.83e+02']
Constraint losses = 1.75e-02
Relative update of tokamak psi (in plasma core): 2.11e-02
Handing off to forward solve (with updated currents).
Relative error = 3.30e-01
-----
Iteration: 5
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-2.13e+04', '-1.35e+04', '-4.55e+02', '4.35e+03', '4.99e+03', '4.58e+03', '4.04e+03', '3.57e+03', '3.05e+03', '1.98e+03', '-9.16e+02', '-6.14e+03', '-6.88e+03', '1.01e+04', '1.38e+04', '-8.25e+03', '-1.23e+04', '-8.85e+03', '-4.88e+03', '-2.30e+02', '6.19e+03', '1.33e+04', '1.33e+04', '-2.81e+03']
Constraint losses = 1.37e-02
Relative update of tokamak psi (in plasma core): 1.67e-02
Handing off to forward solve (with updated currents).
Relative error = 5.51e-02
-----
Iteration: 6
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.20e+04', '-6.43e+03', '4.86e+02', '3.03e+03', '3.36e+03', '3.12e+03', '2.82e+03', '2.59e+03', '2.42e+03', '2.16e+03', '1.41e+03', '2.12e+02', '-2.46e+02', '8.66e+02', '-1.53e+03', '-4.95e+03', '-3.40e+03', '-1.53e+03', '-1.13e+02', '1.25e+03', '2.72e+03', '3.55e+03', '1.33e+03', '-5.85e+03']
Constraint losses = 1.21e-02
Relative update of tokamak psi (in plasma core): 1.59e-02
Handing off to forward solve (with updated currents).
Relative error = 2.66e-02
-----
Iteration: 7
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-1.11e+04', '-4.54e+03', '2.08e+03', '3.33e+03', '2.70e+03', '2.00e+03', '1.48e+03', '1.09e+03', '7.53e+02', '4.91e+02', '5.71e+02', '1.16e+03', '1.68e+03', '4.92e+02', '-3.49e+03', '-5.18e+03', '-2.91e+03', '-1.13e+03', '5.04e+01', '1.07e+03', '2.04e+03', '2.22e+03', '-2.59e+02', '-6.60e+03']
Constraint losses = 1.19e-02
Relative update of tokamak psi (in plasma core): 1.57e-02
Handing off to forward solve (with updated currents).
Relative error = 2.13e-02
-----
Iteration: 8
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['1.00e+04', '6.26e+03', '1.92e+03', '1.63e+03', '2.19e+03', '2.42e+03', '2.47e+03', '2.53e+03', '2.66e+03', '2.56e+03', '1.20e+03', '-1.71e+03', '-4.24e+03', '-3.05e+03', '3.76e+03', '9.08e+03', '6.64e+03', '3.86e+03', '1.81e+03', '-1.55e+02', '-2.46e+03', '-4.45e+03', '-3.41e+03', '3.24e+03']
Constraint losses = 1.57e-02
Relative update of tokamak psi (in plasma core): 2.21e-02
Handing off to forward solve (with updated currents).
Relative error = 1.20e-02
-----
Iteration: 9
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['-5.00e+03', '-2.79e+03', '-1.05e+03', '-8.25e+02', '-9.22e+02', '-9.24e+02', '-8.74e+02', '-8.28e+02', '-7.73e+02', '-5.10e+02', '5.67e+02', '2.77e+03', '3.86e+03', '-1.24e+03', '-6.65e+03', '-4.24e+03', '-1.19e+03', '-8.72e+01', '3.07e+02', '4.18e+02', '1.91e+02', '-7.05e+02', '-2.40e+03', '-4.32e+03']
Constraint losses = 1.81e-02
Relative update of tokamak psi (in plasma core): 1.96e-02
Handing off to forward solve (with updated currents).
Relative error = 1.12e-02
-----
Iteration: 10
Using simplified Green's Jacobian (of constraints wrt coil currents) to optimise the currents.
Change in coil currents (being controlled): ['4.72e+03', '3.10e+03', '6.97e+02', '6.47e+02', '1.09e+03', '1.32e+03', '1.42e+03', '1.53e+03', '1.70e+03', '1.87e+03', '1.63e+03', '1.03e+03', '-2.22e+02', '-3.73e+03', '-2.89e+03', '3.78e+03', '4.59e+03', '3.36e+03', '2.09e+03', '6.56e+02', '-1.34e+03', '-3.65e+03', '-4.17e+03', '-2.38e+02']
Constraint losses = 5.92e-03
Relative update of tokamak psi (in plasma core): 7.65e-03
Handing off to forward solve (with updated currents).
Relative error = 1.44e-03
-----
Iteration: 11
Using full Jacobian (of constraints wrt coil currents) to optimsise currents.
- calculating derivatives for coil 1/24
- calculating derivatives for coil 2/24
- calculating derivatives for coil 3/24
- calculating derivatives for coil 4/24
- calculating derivatives for coil 5/24
- calculating derivatives for coil 6/24
- calculating derivatives for coil 7/24
- calculating derivatives for coil 8/24
- calculating derivatives for coil 9/24
- calculating derivatives for coil 10/24
- calculating derivatives for coil 11/24
- calculating derivatives for coil 12/24
- calculating derivatives for coil 13/24
- calculating derivatives for coil 14/24
- calculating derivatives for coil 15/24
- calculating derivatives for coil 16/24
- calculating derivatives for coil 17/24
- calculating derivatives for coil 18/24
- calculating derivatives for coil 19/24
- calculating derivatives for coil 20/24
- calculating derivatives for coil 21/24
- calculating derivatives for coil 22/24
- calculating derivatives for coil 23/24
- calculating derivatives for coil 24/24
Change in coil currents (being controlled): ['-1.48e+00', '-1.45e+00', '-1.17e+00', '-8.54e-01', '-5.93e-01', '-4.18e-01', '-2.98e-01', '-2.13e-01', '-8.49e-02', '-9.49e-03', '-1.15e-03', '6.16e-02', '9.92e-02', '9.76e-03', '3.87e-02', '-5.28e-02', '-1.11e-01', '-1.63e-01', '-2.29e-01', '-3.12e-01', '-4.30e-01', '-5.91e-01', '-8.28e-01', '-1.17e+00']
Constraint losses = 5.41e-03
Relative update of tokamak psi (in plasma core): 7.58e-06
Handing off to forward solve (with updated currents).
Relative error = 8.02e-06
-----
Inverse static solve SUCCESS. Tolerance 8.02e-06 (vs. requested 1e-05) reached in 12/100 iterations.
# plot the equilibrium
fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), 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
eq.plot(axis=ax1,show=False) # plots the equilibrium
constrain.plot(axis=ax1, show=False) # plots the contraints
<Axes: xlabel='Major radius [m]', ylabel='Height [m]'>
