Example: Obtaining results from the eq and profiles objects
Here we will take a look at how to access different results stored (or calculated using methods) in the eq and profiles objects.
To do this we first need to run a static forward simulation in the MAST-U-like tokamak using previously saved coil currents to generate our results. Of course, the same methods below can be applied after an inverse solve!
Static forward simulation
import matplotlib.pyplot as plt
import numpy as np
# build machine
from freegsnke import build_machine
tokamak = build_machine.tokamak(
active_coils_path=f"../machine_configs/MAST-U/MAST-U_like_active_coils.pickle",
passive_coils_path=f"../machine_configs/MAST-U/MAST-U_like_passive_coils.pickle",
limiter_path=f"../machine_configs/MAST-U/MAST-U_like_limiter.pickle",
wall_path=f"../machine_configs/MAST-U/MAST-U_like_wall.pickle",
)
# initialise equilibrium object
from freegsnke import equilibrium_update
eq = equilibrium_update.Equilibrium(
tokamak=tokamak,
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 profile object
from freegsnke.jtor_update import ConstrainPaxisIp
profiles = ConstrainPaxisIp(
eq=eq,
paxis=8.1e3,
Ip=6.2e5,
fvac=0.5,
alpha_m=1.8,
alpha_n=1.2
)
# initialise solver
from freegsnke import GSstaticsolver
GSStaticSolver = GSstaticsolver.NKGSsolver(eq)
# set coil currents
import pickle
with open('data/simple_diverted_currents_PaxisIp.pk', 'rb') as f:
current_values = pickle.load(f)
for key in current_values.keys():
eq.tokamak[key].current = current_values[key]
eq.tokamak["P6"].current += 100
# carry out forward solve
GSStaticSolver.solve(eq=eq,
profiles=profiles,
constrain=None,
target_relative_tolerance=1e-9)
# plot the resulting equilbrium
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)
ax1.set_xlim(0.1, 2.15)
ax1.set_ylim(-2.25, 2.25)
plt.tight_layout()
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.
Forward static solve SUCCESS. Tolerance 1.93e-10 (vs. requested 1.00e-09) reached in 24/100 iterations.
Dashboard of results
It is worth playing around with the fields/methods in the eq and profiles objects yourself to see which quantities can be calculated from the plasma equilibrium using built-in functionality. Below, we provide a mini dashboard of different quantities and how to generate them.
If there are certain quantites that you wish to be added that don’t exist within FreeGSNKE at the moment, please do submit a feature request to the Github repository.
General equillibrium quantites
Here, we display a number of used equilibrium quantites.
# plasma quantities
print(rf"Plasma current: {eq.plasmaCurrent()} [A]")
print("---")
print(rf"Poloidal beta (definition 1): {eq.poloidalBeta1()}")
print(rf"Poloidal beta (definition 2): {eq.poloidalBeta2()}")
print(rf"Poloidal beta (definition 3): {eq.poloidalBeta3()}")
print(rf"Poloidal beta (definition 4): {eq.poloidalBeta4()}")
print("---")
print(rf"Toroidal beta (definition 1): {eq.toroidalBeta1()}")
print(rf"Toroidal beta (definition 2): {eq.toroidalBeta2()}")
print(rf"Toroidal beta (definition 3): {eq.toroidalBeta3()}")
print(rf"Toroidal beta (definition 4): {eq.toroidalBeta4()}")
print("---")
print(rf"Normalised total beta: {eq.normalised_total_Beta()}")
print("---")
print(fr"Plasma internal inductance (definition 1): {eq.internalInductance1()}")
print(fr"Plasma internal inductance (definition 2): {eq.internalInductance2()}")
print(fr"Plasma internal inductance (definition 3): {eq.internalInductance3()}")
print(fr"Plasma internal inductance (definition 4): {eq.internalInductance4()}")
print(fr"Plasma internal inductance (definition 5): {eq.internalInductance5()}")
Plasma current: 620000.0 [A]
---
Dicrepancy between psi_func and plasma_psi detected. psi_func has been re-set.
Poloidal beta (definition 1): 0.1936817345907993
Poloidal beta (definition 2): 0.2821178924022873
Poloidal beta (definition 3): 104220.19909336207
Poloidal beta (definition 4): 0.19180164253089887
---
Toroidal beta (definition 1): 0.011019794406383024
Toroidal beta (definition 2): 0.014914930674373989
Toroidal beta (definition 3): 0.014166006179373084
Toroidal beta (definition 4): 0.018988713030348738
---
Normalised total beta: 0.01042656059247916
---
Plasma internal inductance (definition 1): 0.7846804064742476
Plasma internal inductance (definition 2): 0.9284555730637306
Plasma internal inductance (definition 3): 0.7139912365432
Plasma internal inductance (definition 4): 0.6557807890367433
Plasma internal inductance (definition 5): 0.995269997254816
# plasma geometry quantities
print(rf"Minor radius: {eq.minorRadius()} [m]")
print("---")
print(rf"Magnetic axis location: {eq.magneticAxis()[0:2]} [m]")
print(rf"Geometric axis location: {eq.geometricAxis()} [m]")
print(rf"Shafranov shift: {eq.shafranov_shift()} [m]")
print(rf"Inner and outer mipdplane (Z=0) radii: {eq.innerOuterSeparatrix(Z=0)} [m]")
print(rf"Inner and outer mipdplane (Z=0) radii (alternative): {eq.innerOuterSeparatrix2(Z=0)} [m]")
print("---")
print(rf"LCFS circumference: {eq.separatrix_length()} [m]")
print(rf"LCFS area: {eq.separatrix_area()} [m^2]")
print(rf"Plasma volume: {eq.plasmaVolume()} [m^3]")
print("---")
print(rf"Aspect ratio: {eq.aspectRatio()}")
print("---")
print(rf"Geometric elongation: {eq.geometricElongation()}")
print(rf"Geometric elongation (upper): {eq.geometricElongation_upper()}")
print(rf"Geometric elongation (lower): {eq.geometricElongation_lower()}")
print(rf"Effective elongation: {eq.effectiveElongation()}")
print("---")
print(rf"Triangularity: {eq.triangularity()}")
print(rf"Triangularity (upper): {eq.triangularity_upper()}")
print(rf"Triangularity (lower): {eq.triangularity_lower()}")
print("---")
s_uo, s_ui, s_lo, s_li = eq.squareness()
print(rf"Squareness (upper outer): {s_uo}")
print(rf"Squareness (upper inner): {s_ui}")
print(rf"Squareness (lower outer): {s_lo}")
print(rf"Squareness (lower inner): {s_li}")
print("---")
L = eq.closest_wall_point()
print(rf"Point on wall that is closest to the LCFS: {L[0]} [m]")
print(rf"Corresponding distance from point to LCFS: {L[1]} [m]")
print("---")
print(rf"Is this a limited plasma?: {eq.flag_limiter}")
print(rf"Does the core plasma boundary intersect the wall?: {eq.intersectsWall()}")
Minor radius: 0.5469603178958806 [m]
---
Magnetic axis location: [ 0.98048366 -0.00971724] [m]
Geometric axis location: [ 0.90054656 -0.03772983] [m]
Shafranov shift: [0.07993711 0.02801259] [m]
Inner and outer mipdplane (Z=0) radii: (0.3532642706221092, 1.4476337685436893) [m]
Inner and outer mipdplane (Z=0) radii (alternative): [(0.3536175527042067, 0.0), (1.447486838852618, 0.0)] [m]
---
LCFS circumference: 5.200320958063151 [m]
LCFS area: 1.7313991012071586 [m^2]
Plasma volume: 9.405450330147868 [m^3]
---
Aspect ratio: 1.646456841041555
---
Geometric elongation: 1.9826035747592203
Geometric elongation (upper): 1.9164845021516368
Geometric elongation (lower): 2.0487226473668034
Effective elongation: 1.7686098701478838
---
Triangularity: 0.46117433164616467
Triangularity (upper): 0.40766884702212436
Triangularity (lower): 0.5146798162702049
---
Squareness (upper outer): -0.07021099683353893
Squareness (upper inner): -0.15820810643352018
Squareness (lower outer): -0.09716133549580039
Squareness (lower inner): -0.30981331781450683
---
Point on wall that is closest to the LCFS: [ 1.4099508 -0.22000329] [m]
Corresponding distance from point to LCFS: 0.11502924090715125 [m]
---
Is this a limited plasma?: False
Does the core plasma boundary intersect the wall?: False
print(fr"Poloidal flux on magnetic axis: {eq.psi_axis} [Webers/2\pi]")
print(fr"Poloidal flux on plasma boundary: {eq.psi_bndry} [Webers/2\pi]")
Poloidal flux on magnetic axis: 0.0978919222684091 [Webers/2\pi]
Poloidal flux on plasma boundary: 0.03027176752441405 [Webers/2\pi]
# extract the plasma core boundary (specify the number of points you want)
eq.separatrix(ntheta=20)
array([[ 0.98048366, 0.92249757],
[ 1.21387596, 0.70859038],
[ 1.33567521, 0.47916198],
[ 1.38835407, 0.28661796],
[ 1.42927248, 0.13610308],
[ 1.44747551, -0.00971724],
[ 1.42690582, -0.15476859],
[ 1.38735701, -0.30532804],
[ 1.3363373 , -0.49950775],
[ 1.21290871, -0.72504799],
[ 0.98048366, -0.94157961],
[ 0.61903712, -1.12213531],
[ 0.46337981, -0.72144963],
[ 0.39201072, -0.43726786],
[ 0.36176471, -0.21075122],
[ 0.35358624, -0.00971724],
[ 0.36218169, 0.19118125],
[ 0.39312694, 0.4170224 ],
[ 0.46655168, 0.69764945],
[ 0.63998828, 1.03821978]])
# extract the X points (R, Z, psi)
eq.xpt
array([[ 6.10134679e-01, -1.12513012e+00, 3.02717675e-02],
[ 6.10388837e-01, 1.11325556e+00, 2.96505117e-02],
[ 1.79006012e+00, -2.04057217e+00, 9.08956348e-03],
[ 1.79065390e+00, 2.04121340e+00, 8.92577674e-03],
[ 1.06978437e+00, -1.45302603e+00, 6.11273360e-03],
[ 1.06929267e+00, 1.45175065e+00, 5.72822488e-03],
[ 1.50382342e+00, -1.51407196e+00, 2.23478535e-03],
[ 1.50411773e+00, 1.51406254e+00, 1.97038970e-03],
[ 1.83128199e+00, -8.33987223e-01, -2.19613143e-02],
[ 1.82831778e+00, 8.33759568e-01, -2.23205025e-02]])
# extract the strike points where the separatric intersects the wall, if any (R, Z)
eq.strikepoints()
array([[ 0.99648299, -1.97607357],
[ 0.332797 , -1.27798218],
[ 0.332797 , 1.2831153 ],
[ 0.99479494, 1.9671284 ]])
# here we plot some of the features above on different axes to show where they are located (see code documentation for more details)
sep = eq.separatrix(ntheta=360)
strikes = eq.strikepoints()
fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 8), dpi=80)
plt.subplots_adjust(wspace=0.25) # adjust the horizontal space between subplots
ax1.grid(True, which='both')
eq.tokamak.plot(axis=ax1,show=False)
ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)
ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], linestyles='solid', colors='r') #, label="Separatrix")
ax1.scatter(strikes[:,0], strikes[:,1], color='k', marker='x', s=40, zorder=2, label="Strikepoints")
ax1.set_aspect('equal')
ax1.set_xlim(0.1, 2.15)
ax1.set_ylim(-2.25, 2.25)
ax1.set_xlabel(r'Major radius, $R$ [m]')
ax1.set_ylabel(r'Height, $Z$ [m]')
ax1.legend(loc="upper right")
ax2.grid(True, which='both')
eq.tokamak.plot(axis=ax2,show=False)
ax2.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)
ax2.scatter(sep[:,0], sep[:,1], color='k', marker='o', s=2, zorder=2, label="LCFS (plasma boundary)")
ax2.scatter(eq.xpt[:,0], eq.xpt[:,1], color='r', marker='x', s=40, zorder=2, label="X-points")
ax2.scatter(eq.magneticAxis()[0], eq.magneticAxis()[1], color='g', marker='x', s=40, zorder=2, label="Magnetic axis")
ax2.scatter(eq.geometricAxis()[0], eq.geometricAxis()[1], color='b', marker='x', s=40, zorder=2, label="Geometric axis")
ax2.set_aspect('equal')
ax2.set_xlim(0.1, 2.15)
ax2.set_ylim(-2.25, 2.25)
ax2.set_xlabel(r'Major radius, $R$ [m]')
ax2.set_ylabel(r'Height, $Z$ [m]')
ax2.legend(loc="upper right")
ax3.grid(True, which='both')
eq.tokamak.plot(axis=ax3,show=False)
ax3.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)
ax3.scatter(sep[:,0], sep[:,1], color='k', marker='o', s=2, zorder=2, label="LCFS (plasma boundary)")
ax3.scatter(eq._sep_Rmin, eq._sep_ZRmin, color='r', marker='x', s=40, zorder=2, label="(Rmin, ZRmin)")
ax3.scatter(eq._sep_Rmax, eq._sep_ZRmax, color='m', marker='x', s=40, zorder=2, label="(Rmax, ZRmax)")
ax3.scatter(eq._sep_RZmin, eq._sep_Zmin, color='g', marker='x', s=40, zorder=2, label="(RZmin, Zmin)")
ax3.scatter(eq._sep_RZmax, eq._sep_Zmax, color='b', marker='x', s=40, zorder=2, label="(ZRmax, Zmax)")
ax3.set_aspect('equal')
ax3.set_xlim(0.1, 2.15)
ax3.set_ylim(-2.25, 2.25)
ax3.set_xlabel(r'Major radius, $R$ [m]')
ax3.set_ylabel(r'Height, $Z$ [m]')
ax3.legend(loc="upper right")
plt.tight_layout()
Here we visualise the dr sep quantities.
These two distances are defined as the radial separation (on the inboard and outboard sides) at the vertical position Z_level between flux surfaces passing through the lower and upper X-points. The default Z_level value is at the geometric axis (i.e. eq.Zgeometric()), however, here we have chosen the midplance (at \(Z=0\)).
dr_seps = eq.dr_sep(Z_level=0)
print(f"dr_sep_in = {dr_seps[0]} [m].")
print(f"dr_sep_out = {dr_seps[1]} [m].")
dr_sep_in = 0.005883440979607579 [m].
dr_sep_out = -0.0027507969968811885 [m].
In the plot below, dr_sep_in is the difference between the radial position of the blue and the red flux surfaces at \(Z=0\).
Similarly, dr_sep_out is the difference between the radial position of the red and the blue flux surfaces at \(Z=0\).
import shapely as sh
# visualising these quantities
psi_bndry = eq.psi_bndry
xpts = eq.xpt
# define your polygon (example)
polygon = sh.Polygon(np.array([eq.tokamak.wall.R, eq.tokamak.wall.Z]).T)
# vectorized-ish approach: boolean mask for points inside
mask = np.array([polygon.contains(sh.Point(x, y)) for x, y in xpts[:,0:2]])
# select points inside as a NumPy array
xpts_inside_wall = xpts[mask,:]
# get indices of the two cloest to magnetic axis
closest_xpts_idxs = np.argsort(np.linalg.norm(xpts_inside_wall[:,0:2] - eq.opt[0,0:2], axis=1))
closest_xpts_sorted = xpts_inside_wall[closest_xpts_idxs[0:2],:]
# plot the flux contours for the values of psi_boundary at each X-point
fig1, axes = plt.subplots(1, 3, figsize=(15, 8), dpi=80)
plt.subplots_adjust(wspace=0.25) # adjust the horizontal space between subplots
axes[0].grid(True, which='both')
eq.tokamak.plot(axis=axes[0],show=False)
axes[0].fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)
axes[0].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[0, 2]], linestyles='solid', colors='r')
axes[0].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[1, 2]], linestyles='dashed', colors='b')
axes[0].scatter(closest_xpts_sorted[0,0], closest_xpts_sorted[0,1], color='r', marker='x', s=40, zorder=2, label="Separatrix from lower X-point")
axes[0].scatter(closest_xpts_sorted[1,0], closest_xpts_sorted[1,1], color='b', marker='x', s=40, zorder=2, label="Separatrix from upper X-point")
axes[0].set_aspect('equal')
axes[0].set_xlim(0.1, 2.15)
axes[0].set_ylim(-2.25, 2.25)
axes[0].set_xlabel(r'Major radius, $R$ [m]')
axes[0].set_ylabel(r'Height, $Z$ [m]')
axes[0].legend(loc="upper right")
axes[1].grid(True, which='both')
eq.tokamak.plot(axis=axes[1],show=False)
axes[1].fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)
axes[1].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[0, 2]], linestyles='solid', colors='r')
axes[1].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[1, 2]], linestyles='dashed', colors='b')
axes[1].scatter(closest_xpts_sorted[0,0], closest_xpts_sorted[0,1], color='r', marker='x', s=40, zorder=2, label="Separatrix from lower X-point")
axes[1].scatter(closest_xpts_sorted[1,0], closest_xpts_sorted[1,1], color='b', marker='x', s=40, zorder=2, label="Separatrix from upper X-point")
axes[1].set_aspect('equal')
axes[1].set_xlim(0.18, 0.38)
axes[1].set_ylim(-0.05, 0.05)
axes[1].set_xlabel(r'Major radius, $R$ [m]')
axes[1].set_ylabel(r'Height, $Z$ [m]')
axes[1].set_title("Inboard side: dr_sep_in")
axes[2].grid(True, which='both')
eq.tokamak.plot(axis=axes[2],show=False)
axes[2].fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)
axes[2].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[0, 2]], linestyles='solid', colors='r')
axes[2].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[1, 2]], linestyles='dashed', colors='b')
axes[2].scatter(closest_xpts_sorted[0,0], closest_xpts_sorted[0,1], color='r', marker='x', s=40, zorder=2, label="Separatrix from lower X-point")
axes[2].scatter(closest_xpts_sorted[1,0], closest_xpts_sorted[1,1], color='b', marker='x', s=40, zorder=2, label="Separatrix from upper X-point")
axes[2].set_aspect('equal')
axes[2].set_xlim(1.3, 1.5)
axes[2].set_ylim(-0.05, 0.05)
axes[2].set_xlabel(r'Major radius, $R$ [m]')
axes[2].set_ylabel(r'Height, $Z$ [m]')
axes[2].set_title("Outboard side: dr_sep_out")
Text(0.5, 1.0, 'Outboard side: dr_sep_out')
Flux quantites
Here, we plot the total poloidal flux \(\psi\) [Webers / \(2\pi\)] and its two components \(\psi_p\) (the plasma flux) and \(\psi_c\) (the coil flux): \(\psi = \psi_p + \psi_c\).
fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 8), dpi=80)
plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots
ax1.grid(True, which='both')
eq.tokamak.plot(axis=ax1,show=False)
# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
im1 = ax1.contour(eq.R, eq.Z, eq.psi(), levels=50) # total psi
ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], colors='r') # psi boundary contour
ax1.set_aspect('equal')
ax1.set_xlim(0.1, 2.15)
ax1.set_ylim(-2.25, 2.25)
ax1.set_title(r"Total flux, $\psi$")
ax1.set_xlabel(r'Major radius, $R$ [m]')
ax1.set_ylabel(r'Height, $Z$ [m]')
cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)
ax2.grid(True, which='both')
eq.tokamak.plot(axis=ax2,show=False)
# ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
im2 = ax2.contour(eq.R, eq.Z, eq.plasma_psi, levels=50) # plasma psi
ax2.set_aspect('equal')
ax2.set_xlim(0.1, 2.15)
ax2.set_ylim(-2.25, 2.25)
ax2.set_title(r"Plasma flux, $\psi_p$")
ax2.set_xlabel(r'Major radius, $R$ [m]')
ax2.set_ylabel(r'Height, $Z$ [m]')
cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)
ax3.grid(True, which='both')
eq.tokamak.plot(axis=ax3,show=False)
# ax3.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax3.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
im3 = ax3.contour(eq.R, eq.Z, eq.tokamak_psi, levels=50) # coil psi
ax3.set_aspect('equal')
ax3.set_xlim(0.1, 2.15)
ax3.set_ylim(-2.25, 2.25)
ax3.set_title(r"Coil flux, $\psi_c$")
ax3.set_xlabel(r'Major radius, $R$ [m]')
ax3.set_ylabel(r'Height, $Z$ [m]')
cbar = plt.colorbar(im3, ax=ax3, fraction=0.09)
plt.tight_layout()
We can also extract the flux produced by individual coils (or pasive structures) if required. We do this my using the pre-calculated Green’s functions, the coil currents, multipliers, and polarities.
currents = eq.tokamak.getCurrents()
# calculate the active/passive coil fluxes
psi_coils = dict()
for i, name in enumerate(currents.keys()):
coil = eq.tokamak.coils_dict[name]
scaling = coil["multiplier"]*coil["polarity"]
greens_matrix = 0.0
if type(eq._pgreen[name]) is dict:
num_coils = len(eq._pgreen[name])
for i, ind in enumerate(eq._pgreen[name]):
greens_matrix += eq._pgreen[name][ind]*scaling[i*(len(scaling)//num_coils)]
else:
num_coils = 1
greens_matrix = eq._pgreen[name]*scaling[0]
psi_coils[name] = greens_matrix*currents[name]
Here, we plot a few of the fluxes from the coils. You can change the name parameters to visualise different coil fluxes.
fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 8), dpi=80)
plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots
name = "Solenoid"
ax1.grid(True, which='both')
eq.tokamak.plot(axis=ax1,show=False)
# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
im1 = ax1.contour(eq.R, eq.Z, psi_coils[name], levels=50)
ax1.set_aspect('equal')
ax1.set_xlim(0.1, 2.15)
ax1.set_ylim(-2.25, 2.25)
ax1.set_title(rf"$\psi_c$ ({name})")
ax1.set_xlabel(r'Major radius, $R$ [m]')
ax1.set_ylabel(r'Height, $Z$ [m]')
cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)
name = "D1"
ax2.grid(True, which='both')
eq.tokamak.plot(axis=ax2,show=False)
# ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
im2 = ax2.contour(eq.R, eq.Z, psi_coils[name], levels=50)
ax2.set_aspect('equal')
ax2.set_xlim(0.1, 2.15)
ax2.set_ylim(-2.25, 2.25)
ax2.set_title(rf"$\psi_c$ ({name})")
ax2.set_xlabel(r'Major radius, $R$ [m]')
ax2.set_ylabel(r'Height, $Z$ [m]')
cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)
name = "P6"
ax3.grid(True, which='both')
eq.tokamak.plot(axis=ax3,show=False)
# ax3.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax3.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
im3 = ax3.contour(eq.R, eq.Z, psi_coils[name], levels=50)
ax3.set_aspect('equal')
ax3.set_xlim(0.1, 2.15)
ax3.set_ylim(-2.25, 2.25)
ax3.set_title(rf"$\psi_c$ ({name})")
ax3.set_xlabel(r'Major radius, $R$ [m]')
ax3.set_ylabel(r'Height, $Z$ [m]')
cbar = plt.colorbar(im3, ax=ax3, fraction=0.09)
plt.tight_layout()
Magnetic fields
Here, we plot the different magnetic field components over the domain.
fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 8), dpi=80)
plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots
ax1.grid(True, which='both')
eq.tokamak.plot(axis=ax1,show=False)
# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
im1 = ax1.contour(eq.R, eq.Z, eq.Br(eq.R, eq.Z), levels=200)
ax1.set_aspect('equal')
ax1.set_xlim(0.1, 2.15)
ax1.set_ylim(-2.25, 2.25)
ax1.set_title(r"$B_R$ [T]")
ax1.set_xlabel(r'Major radius, $R$ [m]')
ax1.set_ylabel(r'Height, $Z$ [m]')
cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)
ax2.grid(True, which='both')
eq.tokamak.plot(axis=ax2,show=False)
# ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
im2 = ax2.contour(eq.R, eq.Z, eq.Btor(eq.R, eq.Z), levels=100)
ax2.set_aspect('equal')
ax2.set_xlim(0.1, 2.15)
ax2.set_ylim(-2.25, 2.25)
ax2.set_title(r"$B_{\phi}$ [T]")
ax2.set_xlabel(r'Major radius, $R$ [m]')
ax2.set_ylabel(r'Height, $Z$ [m]')
cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)
ax3.grid(True, which='both')
eq.tokamak.plot(axis=ax3,show=False)
# ax3.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax3.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
im3 = ax3.contour(eq.R, eq.Z, eq.Bz(eq.R, eq.Z), levels=100)
ax3.set_aspect('equal')
ax3.set_xlim(0.1, 2.15)
ax3.set_ylim(-2.25, 2.25)
ax3.set_title(r"$B_Z$ [T]")
ax3.set_xlabel(r'Major radius, $R$ [m]')
ax3.set_ylabel(r'Height, $Z$ [m]')
cbar = plt.colorbar(im3, ax=ax3, fraction=0.09)
plt.tight_layout()
# this is just the contribution from the plasma
fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 8), dpi=80)
plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots
ax1.grid(True, which='both')
eq.tokamak.plot(axis=ax1,show=False)
# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
im1 = ax1.contour(eq.R, eq.Z, eq.plasmaBr(eq.R, eq.Z), levels=30)
ax1.set_aspect('equal')
ax1.set_xlim(0.1, 2.15)
ax1.set_ylim(-2.25, 2.25)
ax1.set_title(r"$B_R$ (from plasma) [T]")
ax1.set_xlabel(r'Major radius, $R$ [m]')
ax1.set_ylabel(r'Height, $Z$ [m]')
cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)
ax2.grid(True, which='both')
eq.tokamak.plot(axis=ax2,show=False)
# ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
im2 = ax2.contour(eq.R, eq.Z, eq.plasmaBz(eq.R, eq.Z), levels=30)
ax2.set_aspect('equal')
ax2.set_xlim(0.1, 2.15)
ax2.set_ylim(-2.25, 2.25)
ax2.set_title(r"$B_Z$ (from plasma) [T]")
ax2.set_xlabel(r'Major radius, $R$ [m]')
ax2.set_ylabel(r'Height, $Z$ [m]')
cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)
plt.tight_layout()
# this is just the contribution from the active coils
Br_actives = 0.0
Bz_actives = 0.0
for name in eq.tokamak.coils_list[0:12]:
Br_actives += eq.tokamak[name].Br(eq.R, eq.Z)
Bz_actives += eq.tokamak[name].Bz(eq.R, eq.Z)
fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 8), dpi=80)
plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots
ax1.grid(True, which='both')
eq.tokamak.plot(axis=ax1,show=False)
# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
im1 = ax1.contour(eq.R, eq.Z, Br_actives, levels=200)
ax1.set_aspect('equal')
ax1.set_xlim(0.1, 2.15)
ax1.set_ylim(-2.25, 2.25)
ax1.set_title(r"$B_R$ (from coils) [T]")
ax1.set_xlabel(r'Major radius, $R$ [m]')
ax1.set_ylabel(r'Height, $Z$ [m]')
cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)
ax2.grid(True, which='both')
eq.tokamak.plot(axis=ax2,show=False)
# ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
im2 = ax2.contour(eq.R, eq.Z, Bz_actives, levels=200)
ax2.set_aspect('equal')
ax2.set_xlim(0.1, 2.15)
ax2.set_ylim(-2.25, 2.25)
ax2.set_title(r"$B_Z$ (from coils) [T]")
ax2.set_xlabel(r'Major radius, $R$ [m]')
ax2.set_ylabel(r'Height, $Z$ [m]')
cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)
plt.tight_layout()
# # this is just the contribution from the passive coils (note this is zero because there are no currents in the passives!)
# Br_passives = 0.0
# Bz_passives = 0.0
# for name in eq.tokamak.coils_list[12:]:
# Br_passives += eq.tokamak[name].Br(eq.R, eq.Z)
# Bz_passives += eq.tokamak[name].Bz(eq.R, eq.Z)
# fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 8), dpi=80)
# plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots
# ax1.grid(True, which='both')
# eq.tokamak.plot(axis=ax1,show=False)
# # ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
# ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
# im1 = ax1.contour(eq.R, eq.Z, Br_passives, levels=200)
# ax1.set_aspect('equal')
# ax1.set_xlim(0.1, 2.15)
# ax1.set_ylim(-2.25, 2.25)
# ax1.set_title(r"$B_R$ (from coils) [T]")
# ax1.set_xlabel(r'Major radius, $R$ [m]')
# ax1.set_ylabel(r'Height, $Z$ [m]')
# cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)
# ax2.grid(True, which='both')
# eq.tokamak.plot(axis=ax2,show=False)
# # ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
# ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
# im2 = ax2.contour(eq.R, eq.Z, Bz_passives, levels=200)
# ax2.set_aspect('equal')
# ax2.set_xlim(0.1, 2.15)
# ax2.set_ylim(-2.25, 2.25)
# ax2.set_title(r"$B_Z$ (from coils) [T]")
# ax2.set_xlabel(r'Major radius, $R$ [m]')
# ax2.set_ylabel(r'Height, $Z$ [m]')
# cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)
# plt.tight_layout()
# total poloidal magnetic field
fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=100)
ax1.grid(True, which='both')
eq.tokamak.plot(axis=ax1,show=False)
# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
im1 = ax1.contour(eq.R, eq.Z, eq.Bpol(eq.R, eq.Z), levels=200)
ax1.set_aspect('equal')
ax1.set_xlim(0.1, 2.15)
ax1.set_ylim(-2.25, 2.25)
ax1.set_title(r"$B_{pol}$ [T]")
ax1.set_xlabel(r'Major radius, $R$ [m]')
ax1.set_ylabel(r'Height, $Z$ [m]')
cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)
plt.tight_layout()
Coil currents
Here, we visualise the size of the currents in the poloidal field coils.
# listing the currents may not be that informative
eq.tokamak.getCurrents()
{'Solenoid': 5000.0,
'PX': 3696.8090944766755,
'D1': 6941.236553652974,
'D2': 4390.861230885664,
'D3': 2924.355186006588,
'Dp': -2488.1988314333917,
'D5': 249.19273211994647,
'D6': -314.9949422096607,
'D7': 511.07691868766733,
'P4': -3564.379440164539,
'P5': -3954.2343617676433,
'P6': 100.00000687925065,
'vessel_1': 0.0,
'vessel_2': 0.0,
'vessel_3': 0.0,
'vessel_4': 0.0,
'vessel_5': 0.0,
'vessel_6': 0.0,
'vessel_7': 0.0,
'vessel_8': 0.0,
'vessel_9': 0.0,
'vessel_10': 0.0,
'vessel_11': 0.0,
'vessel_12': 0.0,
'vessel_13': 0.0,
'vessel_14': 0.0,
'vessel_15': 0.0,
'vessel_16': 0.0,
'vessel_17': 0.0,
'vessel_18': 0.0,
'vessel_19': 0.0,
'vessel_20': 0.0,
'vessel_21': 0.0,
'vessel_22': 0.0,
'centrecolumn_1': 0.0,
'centrecolumn_2': 0.0,
'centrecolumn_3': 0.0,
'centrecolumn_4': 0.0,
'centrecolumn_5': 0.0,
'centrecolumn_6': 0.0,
'centrecolumn_7': 0.0,
'centrecolumn_8': 0.0,
'centrecolumn_9': 0.0,
'centrecolumn_10': 0.0,
'colosseum_upper_1': 0.0,
'colosseum_upper_2': 0.0,
'colosseum_upper_3': 0.0,
'colosseum_lower_1': 0.0,
'colosseum_lower_2': 0.0,
'colosseum_lower_3': 0.0,
'colosseum_outer_upper_1': 0.0,
'colosseum_outer_upper_2': 0.0,
'colosseum_outer_lower_1': 0.0,
'colosseum_outer_lower_2': 0.0,
'ring_plate_minor_upper_1': 0.0,
'ring_plate_minor_upper_2': 0.0,
'ring_plate_minor_lower_1': 0.0,
'ring_plate_minor_lower_2': 0.0,
'ring_plate_major_upper_1': 0.0,
'ring_plate_major_lower_1': 0.0,
'psp_upper_1': 0.0,
'psp_lower_1': 0.0,
'gas_baffle_upper_1': 0.0,
'gas_baffle_upper_2': 0.0,
'gas_baffle_upper_3': 0.0,
'gas_baffle_lower_1': 0.0,
'gas_baffle_lower_2': 0.0,
'gas_baffle_lower_3': 0.0,
'd1_case_upper_0': 0.0,
'd1_case_upper_1': 0.0,
'd1_case_upper_2': 0.0,
'd1_case_upper_3': 0.0,
'd1_case_lower_0': 0.0,
'd1_case_lower_1': 0.0,
'd1_case_lower_2': 0.0,
'd1_case_lower_3': 0.0,
'd2_case_upper_0': 0.0,
'd2_case_upper_1': 0.0,
'd2_case_upper_2': 0.0,
'd2_case_upper_3': 0.0,
'd2_case_lower_0': 0.0,
'd2_case_lower_1': 0.0,
'd2_case_lower_2': 0.0,
'd2_case_lower_3': 0.0,
'd3_case_upper_0': 0.0,
'd3_case_upper_1': 0.0,
'd3_case_upper_2': 0.0,
'd3_case_upper_3': 0.0,
'd3_case_lower_0': 0.0,
'd3_case_lower_1': 0.0,
'd3_case_lower_2': 0.0,
'd3_case_lower_3': 0.0,
'd5_case_upper_0': 0.0,
'd5_case_upper_1': 0.0,
'd5_case_upper_2': 0.0,
'd5_case_upper_3': 0.0,
'd5_case_lower_0': 0.0,
'd5_case_lower_1': 0.0,
'd5_case_lower_2': 0.0,
'd5_case_lower_3': 0.0,
'd6_case_upper_0': 0.0,
'd6_case_upper_1': 0.0,
'd6_case_upper_2': 0.0,
'd6_case_upper_3': 0.0,
'd6_case_lower_0': 0.0,
'd6_case_lower_1': 0.0,
'd6_case_lower_2': 0.0,
'd6_case_lower_3': 0.0,
'd7_case_upper_0': 0.0,
'd7_case_upper_1': 0.0,
'd7_case_upper_2': 0.0,
'd7_case_upper_3': 0.0,
'd7_case_lower_0': 0.0,
'd7_case_lower_1': 0.0,
'd7_case_lower_2': 0.0,
'd7_case_lower_3': 0.0,
'dp_case_upper_0': 0.0,
'dp_case_upper_1': 0.0,
'dp_case_upper_2': 0.0,
'dp_case_upper_3': 0.0,
'dp_case_lower_0': 0.0,
'dp_case_lower_1': 0.0,
'dp_case_lower_2': 0.0,
'dp_case_lower_3': 0.0,
'p4_case_upper_0': 0.0,
'p4_case_upper_1': 0.0,
'p4_case_upper_2': 0.0,
'p4_case_upper_3': 0.0,
'p4_case_lower_0': 0.0,
'p4_case_lower_1': 0.0,
'p4_case_lower_2': 0.0,
'p4_case_lower_3': 0.0,
'p5_case_upper_0': 0.0,
'p5_case_upper_1': 0.0,
'p5_case_upper_2': 0.0,
'p5_case_upper_3': 0.0,
'p5_case_lower_0': 0.0,
'p5_case_lower_1': 0.0,
'p5_case_lower_2': 0.0,
'p5_case_lower_3': 0.0,
'p6_case_upper_0': 0.0,
'p6_case_upper_1': 0.0,
'p6_case_upper_2': 0.0,
'p6_case_upper_3': 0.0,
'p6_case_upper_4': 0.0,
'p6_case_lower_0': 0.0,
'p6_case_lower_1': 0.0,
'p6_case_lower_2': 0.0,
'p6_case_lower_3': 0.0,
'p6_case_lower_4': 0.0}
You can play around with the colour schemes below to visualise how the current is distributed around the machine. Note here that there is no current in the passive structures but they can be included if present.
from matplotlib.patches import Rectangle, Polygon
from matplotlib.colors import Normalize, TwoSlopeNorm
import matplotlib.cm as cm
# create colormap based on magnitude of currents
currents_array = []
for key in list(eq.tokamak.coils_dict.keys()):
currents_array.append(eq.tokamak[key].current)
max_curr = np.max(np.abs(currents_array))
# norm = Normalize(vmin=-max_curr, vmax=max_curr) # alternative colorbar
norm = TwoSlopeNorm(vmin=np.min(currents_array), vcenter=0, vmax=np.max(currents_array))
cmap = cm.bwr
# plot
fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=120)
plt.tight_layout()
ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
ax1.grid(alpha=0.5)
ax1.set_aspect('equal')
ax1.set_xlim(0.1, 2.15)
ax1.set_ylim(-2.25, 2.25)
ax1.set_xlabel(r'Major radius, $R$ [m]')
ax1.set_ylabel(r'Height, $Z$ [m]')
for name in list(eq.tokamak.coils_dict.keys()):
coil = eq.tokamak.coils_dict[name]
current = eq.tokamak[name].current
color = cmap(norm(current)) # map the current to a color
# plot active coils (and currents)
if coil["active"]:
for i in range(0, len(coil["coords"][0,:])):
patch = Rectangle(
(coil["coords"][0,i] - coil["dR"] / 2, coil["coords"][1,i] - coil["dZ"] / 2),
width=coil["dR"],
height=coil["dZ"],
facecolor=color,
edgecolor='k',
linewidth=0.2,
)
ax1.add_patch(patch)
# plot passive structures (currents are zero here)
else:
patch = Polygon(
coil["vertices"].T,
facecolor=color,
edgecolor="k",
linewidth=0.2
)
ax1.add_patch(patch)
# add a colorbar
sm = cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([])
cbar = fig1.colorbar(sm, ax=ax1, orientation='vertical', fraction=0.09)
cbar.set_label('Current [A]')
Forces on the coils
We can also extract the radial and vertical forces on the coils.
eq.printForces()
Forces on coils
Solenoid (circuit)
Solenoid : R = 7967.14 kN , Z = -1.52 kN
PX (circuit)
PX1 : R = 152.82 kN , Z = 13.30 kN
PX2 : R = 153.51 kN , Z = -13.12 kN
D1 (circuit)
D11 : R = 440.20 kN , Z = -38.24 kN
D12 : R = 440.88 kN , Z = 38.59 kN
D2 (circuit)
D21 : R = 46.58 kN , Z = -34.21 kN
D22 : R = 46.80 kN , Z = 34.38 kN
D3 (circuit)
D31 : R = 17.53 kN , Z = -11.58 kN
D32 : R = 17.61 kN , Z = 11.67 kN
Dp (circuit)
Dp1 : R = 45.96 kN , Z = -8.71 kN
Dp2 : R = 45.91 kN , Z = 8.43 kN
D5 (circuit)
D51 : R = -0.27 kN , Z = 0.47 kN
D52 : R = -0.27 kN , Z = -0.47 kN
D6 (circuit)
D61 : R = 1.57 kN , Z = -1.81 kN
D62 : R = 1.58 kN , Z = 1.80 kN
D7 (circuit)
D71 : R = 0.09 kN , Z = 4.31 kN
D72 : R = 0.07 kN , Z = -4.30 kN
P4 (circuit)
P41 : R = 64.10 kN , Z = 8.21 kN
P42 : R = 63.91 kN , Z = -8.04 kN
P5 (circuit)
P51 : R = 115.64 kN , Z = 15.18 kN
P52 : R = 116.44 kN , Z = -14.69 kN
P6 (circuit)
P61 : R = -0.83 kN , Z = -0.55 kN
P62 : R = 0.86 kN , Z = -0.57 kN
vessel_1 : R = 0.00 kN , Z = -0.00 kN
vessel_2 : R = 0.00 kN , Z = -0.00 kN
vessel_3 : R = 0.00 kN , Z = -0.00 kN
vessel_4 : R = 0.00 kN , Z = 0.00 kN
vessel_5 : R = 0.00 kN , Z = -0.00 kN
vessel_6 : R = 0.00 kN , Z = 0.00 kN
vessel_7 : R = 0.00 kN , Z = 0.00 kN
vessel_8 : R = 0.00 kN , Z = 0.00 kN
vessel_9 : R = 0.00 kN , Z = -0.00 kN
vessel_10 : R = 0.00 kN , Z = 0.00 kN
vessel_11 : R = 0.00 kN , Z = -0.00 kN
vessel_12 : R = 0.00 kN , Z = -0.00 kN
vessel_13 : R = 0.00 kN , Z = -0.00 kN
vessel_14 : R = 0.00 kN , Z = -0.00 kN
vessel_15 : R = 0.00 kN , Z = 0.00 kN
vessel_16 : R = 0.00 kN , Z = -0.00 kN
vessel_17 : R = 0.00 kN , Z = 0.00 kN
vessel_18 : R = 0.00 kN , Z = 0.00 kN
vessel_19 : R = 0.00 kN , Z = 0.00 kN
vessel_20 : R = 0.00 kN , Z = 0.00 kN
vessel_21 : R = 0.00 kN , Z = -0.00 kN
vessel_22 : R = 0.00 kN , Z = 0.00 kN
centrecolumn_1 : R = 0.00 kN , Z = 0.00 kN
centrecolumn_2 : R = 0.00 kN , Z = -0.00 kN
centrecolumn_3 : R = 0.00 kN , Z = 0.00 kN
centrecolumn_4 : R = 0.00 kN , Z = -0.00 kN
centrecolumn_5 : R = 0.00 kN , Z = 0.00 kN
centrecolumn_6 : R = 0.00 kN , Z = -0.00 kN
centrecolumn_7 : R = 0.00 kN , Z = 0.00 kN
centrecolumn_8 : R = 0.00 kN , Z = -0.00 kN
centrecolumn_9 : R = 0.00 kN , Z = 0.00 kN
centrecolumn_10 : R = 0.00 kN , Z = -0.00 kN
colosseum_upper_1 : R = 0.00 kN , Z = -0.00 kN
colosseum_upper_2 : R = 0.00 kN , Z = -0.00 kN
colosseum_upper_3 : R = 0.00 kN , Z = -0.00 kN
colosseum_lower_1 : R = 0.00 kN , Z = 0.00 kN
colosseum_lower_2 : R = 0.00 kN , Z = 0.00 kN
colosseum_lower_3 : R = 0.00 kN , Z = 0.00 kN
colosseum_outer_upper_1 : R = 0.00 kN , Z = -0.00 kN
colosseum_outer_upper_2 : R = 0.00 kN , Z = -0.00 kN
colosseum_outer_lower_1 : R = 0.00 kN , Z = 0.00 kN
colosseum_outer_lower_2 : R = 0.00 kN , Z = 0.00 kN
ring_plate_minor_upper_1 : R = 0.00 kN , Z = -0.00 kN
ring_plate_minor_upper_2 : R = 0.00 kN , Z = -0.00 kN
ring_plate_minor_lower_1 : R = 0.00 kN , Z = 0.00 kN
ring_plate_minor_lower_2 : R = 0.00 kN , Z = 0.00 kN
ring_plate_major_upper_1 : R = 0.00 kN , Z = -0.00 kN
ring_plate_major_lower_1 : R = 0.00 kN , Z = 0.00 kN
psp_upper_1 : R = 0.00 kN , Z = -0.00 kN
psp_lower_1 : R = 0.00 kN , Z = 0.00 kN
gas_baffle_upper_1 : R = 0.00 kN , Z = 0.00 kN
gas_baffle_upper_2 : R = 0.00 kN , Z = 0.00 kN
gas_baffle_upper_3 : R = 0.00 kN , Z = 0.00 kN
gas_baffle_lower_1 : R = 0.00 kN , Z = -0.00 kN
gas_baffle_lower_2 : R = 0.00 kN , Z = -0.00 kN
gas_baffle_lower_3 : R = 0.00 kN , Z = -0.00 kN
d1_case_upper_0 : R = 0.00 kN , Z = -0.00 kN
d1_case_upper_1 : R = 0.00 kN , Z = 0.00 kN
d1_case_upper_2 : R = 0.00 kN , Z = -0.00 kN
d1_case_upper_3 : R = 0.00 kN , Z = 0.00 kN
d1_case_lower_0 : R = 0.00 kN , Z = 0.00 kN
d1_case_lower_1 : R = 0.00 kN , Z = 0.00 kN
d1_case_lower_2 : R = 0.00 kN , Z = -0.00 kN
d1_case_lower_3 : R = 0.00 kN , Z = 0.00 kN
d2_case_upper_0 : R = 0.00 kN , Z = -0.00 kN
d2_case_upper_1 : R = 0.00 kN , Z = 0.00 kN
d2_case_upper_2 : R = 0.00 kN , Z = -0.00 kN
d2_case_upper_3 : R = 0.00 kN , Z = -0.00 kN
d2_case_lower_0 : R = 0.00 kN , Z = 0.00 kN
d2_case_lower_1 : R = 0.00 kN , Z = 0.00 kN
d2_case_lower_2 : R = 0.00 kN , Z = -0.00 kN
d2_case_lower_3 : R = 0.00 kN , Z = 0.00 kN
d3_case_upper_0 : R = 0.00 kN , Z = -0.00 kN
d3_case_upper_1 : R = 0.00 kN , Z = 0.00 kN
d3_case_upper_2 : R = 0.00 kN , Z = -0.00 kN
d3_case_upper_3 : R = 0.00 kN , Z = -0.00 kN
d3_case_lower_0 : R = 0.00 kN , Z = 0.00 kN
d3_case_lower_1 : R = 0.00 kN , Z = 0.00 kN
d3_case_lower_2 : R = 0.00 kN , Z = -0.00 kN
d3_case_lower_3 : R = 0.00 kN , Z = 0.00 kN
d5_case_upper_0 : R = 0.00 kN , Z = 0.00 kN
d5_case_upper_1 : R = 0.00 kN , Z = 0.00 kN
d5_case_upper_2 : R = 0.00 kN , Z = -0.00 kN
d5_case_upper_3 : R = 0.00 kN , Z = 0.00 kN
d5_case_lower_0 : R = 0.00 kN , Z = -0.00 kN
d5_case_lower_1 : R = 0.00 kN , Z = 0.00 kN
d5_case_lower_2 : R = 0.00 kN , Z = -0.00 kN
d5_case_lower_3 : R = 0.00 kN , Z = -0.00 kN
d6_case_upper_0 : R = 0.00 kN , Z = 0.00 kN
d6_case_upper_1 : R = 0.00 kN , Z = 0.00 kN
d6_case_upper_2 : R = 0.00 kN , Z = 0.00 kN
d6_case_upper_3 : R = 0.00 kN , Z = 0.00 kN
d6_case_lower_0 : R = 0.00 kN , Z = -0.00 kN
d6_case_lower_1 : R = 0.00 kN , Z = -0.00 kN
d6_case_lower_2 : R = 0.00 kN , Z = 0.00 kN
d6_case_lower_3 : R = 0.00 kN , Z = -0.00 kN
d7_case_upper_0 : R = 0.00 kN , Z = 0.00 kN
d7_case_upper_1 : R = 0.00 kN , Z = -0.00 kN
d7_case_upper_2 : R = 0.00 kN , Z = 0.00 kN
d7_case_upper_3 : R = 0.00 kN , Z = 0.00 kN
d7_case_lower_0 : R = 0.00 kN , Z = -0.00 kN
d7_case_lower_1 : R = 0.00 kN , Z = 0.00 kN
d7_case_lower_2 : R = 0.00 kN , Z = -0.00 kN
d7_case_lower_3 : R = 0.00 kN , Z = -0.00 kN
dp_case_upper_0 : R = 0.00 kN , Z = 0.00 kN
dp_case_upper_1 : R = 0.00 kN , Z = 0.00 kN
dp_case_upper_2 : R = 0.00 kN , Z = -0.00 kN
dp_case_upper_3 : R = 0.00 kN , Z = 0.00 kN
dp_case_lower_0 : R = 0.00 kN , Z = -0.00 kN
dp_case_lower_1 : R = 0.00 kN , Z = -0.00 kN
dp_case_lower_2 : R = 0.00 kN , Z = 0.00 kN
dp_case_lower_3 : R = 0.00 kN , Z = -0.00 kN
p4_case_upper_0 : R = 0.00 kN , Z = -0.00 kN
p4_case_upper_1 : R = 0.00 kN , Z = 0.00 kN
p4_case_upper_2 : R = 0.00 kN , Z = -0.00 kN
p4_case_upper_3 : R = 0.00 kN , Z = -0.00 kN
p4_case_lower_0 : R = 0.00 kN , Z = -0.00 kN
p4_case_lower_1 : R = 0.00 kN , Z = -0.00 kN
p4_case_lower_2 : R = 0.00 kN , Z = 0.00 kN
p4_case_lower_3 : R = 0.00 kN , Z = -0.00 kN
p5_case_upper_0 : R = 0.00 kN , Z = -0.00 kN
p5_case_upper_1 : R = 0.00 kN , Z = 0.00 kN
p5_case_upper_2 : R = 0.00 kN , Z = -0.00 kN
p5_case_upper_3 : R = 0.00 kN , Z = -0.00 kN
p5_case_lower_0 : R = 0.00 kN , Z = -0.00 kN
p5_case_lower_1 : R = 0.00 kN , Z = 0.00 kN
p5_case_lower_2 : R = 0.00 kN , Z = -0.00 kN
p5_case_lower_3 : R = 0.00 kN , Z = 0.00 kN
p6_case_upper_0 : R = 0.00 kN , Z = -0.00 kN
p6_case_upper_1 : R = 0.00 kN , Z = -0.00 kN
p6_case_upper_2 : R = 0.00 kN , Z = -0.00 kN
p6_case_upper_3 : R = 0.00 kN , Z = -0.00 kN
p6_case_upper_4 : R = 0.00 kN , Z = -0.00 kN
p6_case_lower_0 : R = 0.00 kN , Z = 0.00 kN
p6_case_lower_1 : R = 0.00 kN , Z = 0.00 kN
p6_case_lower_2 : R = 0.00 kN , Z = 0.00 kN
p6_case_lower_3 : R = 0.00 kN , Z = 0.00 kN
p6_case_lower_4 : R = 0.00 kN , Z = 0.00 kN
1D plasma current density profiles (and others)
Here, we visualise the \(p'\) and \(FF'\) profiles used in our equilbirium solve.
# plot the p' and FF' profiles
psi_n = eq.psiN_1D(N=65)
fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,6), dpi=80)
ax1.grid(zorder=0, alpha=0.75)
ax1.plot(psi_n, profiles.pprime(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)
ax1.set_xlabel(r'$\hat{\psi}$')
ax1.set_ylabel(r"$p'(\hat{\psi})$")
ax1.ticklabel_format(axis='y', scilimits=(0,0))
ax2.grid(zorder=0, alpha=0.75)
ax2.plot(psi_n, profiles.ffprime(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)
ax2.set_xlabel(r'$\hat{\psi}$')
ax2.set_ylabel(r"$FF'(\hat{\psi})$")
Text(0, 0.5, "$FF'(\\hat{\\psi})$")
# plot the p and F profiles
psi_n = eq.psiN_1D(N=65)
fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,6), dpi=80)
ax1.grid(zorder=0, alpha=0.75)
ax1.plot(psi_n, profiles.pressure(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)
ax1.set_xlabel(r'$\hat{\psi}$')
ax1.set_ylabel(r"$p(\hat{\psi})$")
ax1.ticklabel_format(axis='y', scilimits=(0,0))
ax2.grid(zorder=0, alpha=0.75)
ax2.plot(psi_n, profiles.fpol(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)
ax2.set_xlabel(r'$\hat{\psi}$')
ax2.set_ylabel(r"$F(\hat{\psi})$")
Text(0, 0.5, '$F(\\hat{\\psi})$')
# plot the p and F profiles
psi_n = eq.psiN_1D(N=65)
fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,6), dpi=80)
ax1.grid(zorder=0, alpha=0.75)
ax1.plot(psi_n, profiles.pressure(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)
ax1.set_xlabel(r'$\hat{\psi}$')
ax1.set_ylabel(r"$p(\hat{\psi})$")
ax1.ticklabel_format(axis='y', scilimits=(0,0))
ax2.grid(zorder=0, alpha=0.75)
ax2.plot(psi_n, profiles.fpol(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)
ax2.set_xlabel(r'$\hat{\psi}$')
ax2.set_ylabel(r"$F(\hat{\psi})$")
Text(0, 0.5, '$F(\\hat{\\psi})$')
# plot q profile
psi_n = np.linspace(0.01,0.99,65) # values of q at 0 and 1 can be problematic
fig1, ax1 = plt.subplots(1, 1, figsize=(6,6), dpi=80)
ax1.grid(zorder=0, alpha=0.75)
ax1.plot(psi_n, eq.q(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)
ax1.set_xlabel(r'$\hat{\psi}$')
ax1.set_ylabel(r"$q(\hat{\psi})$")
Text(0, 0.5, '$q(\\hat{\\psi})$')
Toroidal current density
Here, we visualise the toroidal current density inside the core of the plasma.
fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=100)
ax1.grid(True, which='both')
eq.tokamak.plot(axis=ax1,show=False)
ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
im1 = ax1.contourf(eq.R, eq.Z, eq._profiles.jtor, levels=np.linspace(0.01, np.max( eq._profiles.jtor,), 20))
ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], colors='r')
ax1.set_aspect('equal')
ax1.set_xlim(0.1, 2.15)
ax1.set_ylim(-2.25, 2.25)
ax1.set_title(r"Toroidal current density, $J_p$ [A]")
ax1.set_xlabel(r'Major radius, $R$ [m]')
ax1.set_ylabel(r'Height, $Z$ [m]')
cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)
plt.tight_layout()
Masking arrays
There are a number of masking arrays that are built during the equilibrium solve that may be useful for plotting or other calculations you may carry out.
fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 8), dpi=80)
plt.subplots_adjust(wspace=0.25) # adjust the horizontal space between subplots
ax1.grid(True, which='both')
# eq.tokamak.plot(axis=ax1,show=False)
# ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor=None, zorder=0)
# ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], linestyles='solid', colors='r') #, label="Separatrix")
ax1.pcolormesh(eq.R, eq.Z, eq.mask_inside_limiter, cmap="Greys")
ax1.set_aspect('equal')
ax1.set_xlim(0.1, 2.15)
ax1.set_ylim(-2.25, 2.25)
ax1.set_xlabel(r'Major radius, $R$ [m]')
ax1.set_ylabel(r'Height, $Z$ [m]')
ax1.set_title("Mask of the region allowed to the plasma")
ax2.grid(True, which='both')
# eq.tokamak.plot(axis=ax2,show=False)
# ax2.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor=None, zorder=0)
# ax2.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], linestyles='solid', colors='r') #, label="Separatrix")
ax2.pcolormesh(eq.R, eq.Z, eq._profiles.diverted_core_mask, cmap="Greys")
ax2.set_aspect('equal')
ax2.set_xlim(0.1, 2.15)
ax2.set_ylim(-2.25, 2.25)
ax2.set_xlabel(r'Major radius, $R$ [m]')
ax2.set_ylabel(r'Height, $Z$ [m]')
ax2.set_title("Mask within the plasma core boundary")
plt.tight_layout()
Flux averaged quantities
The equilibrium class provides a method to calculate the “flux averaged” value of a user-defined 2D scalar field \(f(R,Z)\), using line integrals, on a given (normalised) flux surface of \(\psi_n\) (within the last closed flux surface). The flux average \(\langle f \rangle\) is given by
\[
\langle f \rangle (\psi_n) = \frac{ \int_{C(\psi_n)} \frac{f(R,Z)}{B_{\text{pol}}(R,Z)}\, ds}{ \int_{C(\psi_n)} \frac{1}{B_{\text{pol}}(R,Z)} \, ds },
\]
where:
\(f(R,Z)\) = 2D scalar field function (e.g. the plasma current density function \(J_p(R,Z)\)).
\(\psi_n\) = value of normalised flux at which to evaluate line integrals.
\(C(\psi_n)\) = curve of \((R, Z)\) points satisfying \(\psi_n = \text{const}\) (i.e. a flux contour).
\(B_{\text{pol}}(R,Z)\) = 2D scalar poloidal magnetic field function.
\(ds\) = the arc length element (where \(ds = \sqrt{(dR/dl)^2 + (dZ/dl)^2} dl\) and \(l \in [0,L]\) is a parameterised length going from the beginning to the end of the contour).
The definition of the flux average was taken from Song et al. (2024).
Here we will calculate the flux averaged values of the plasma current density by first defining a function that returns the current density at arbitrary \((R,Z)\) locations.
from scipy.interpolate import RectBivariateSpline
def f(R,Z):
jtor = RectBivariateSpline(eq.R_1D, eq.Z_1D, eq._profiles.jtor)
return jtor(R, Z, grid=False)
Next we can call the method in the equilibrium object and plot the results. Given the notation above, we note that \(\psi_n\) and \(\hat{\psi}\) are equivalent.
# call the method
flux_averaged_jtor, psi_n = eq.flux_averaged_function(
f=f,
psi_n=np.linspace(0.0,1.0,101)
)
# plot
fig1, ax1 = plt.subplots(1, 1, figsize=(6,6), dpi=80)
ax1.grid(zorder=0, alpha=0.75)
ax1.plot(psi_n, flux_averaged_jtor, color='k', linewidth=1, marker='x', markersize=2, zorder=10)
ax1.set_xlabel(r'$\hat{\psi}$')
ax1.set_ylabel(r"$\langle J_p \rangle (\hat{\psi})$")
Text(0, 0.5, '$\\langle J_p \\rangle (\\hat{\\psi})$')
# # for example you could flux average other quantities of interest
# # 1/R
# def f(R,Z):
# g = RectBivariateSpline(eq.R_1D, eq.Z_1D, 1/eq.R)
# return g(R, Z, grid=False)