Example : Using the magnetic probes
This example will show how to use the probes class.
Generate a starting equilibrum (via a forward solve)
First, we need a tokamak and equilibrium since the probes take properties from the equilibrium as inputs.
We will copy the code from example_1 to generate a sample equilibrium.
import matplotlib.pyplot as plt
# 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",
magnetic_probe_path=f"../machine_configs/MAST-U/MAST-U_like_magnetic_probes.pickle",
)
# initialise the equilibrium
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 the profiles
from freegsnke.jtor_update import ConstrainPaxisIp
profiles = ConstrainPaxisIp(
eq=eq, # equilibrium object
paxis=8e3, # profile object
Ip=6e5, # plasma current
fvac=0.5, # fvac = rB_{tor}
alpha_m=1.8, # profile function parameter
alpha_n=1.2 # profile function parameter
)
# load the nonlinear solver
from freegsnke import GSstaticsolver
GSStaticSolver = GSstaticsolver.NKGSsolver(eq)
# set the coil currents
import pickle
with open('simple_diverted_currents_PaxisIp.pk', 'rb') as f:
current_values = pickle.load(f)
for key in current_values.keys():
eq.tokamak.set_coil_current(key, current_values[key])
# carry out the foward solve to find the equilibrium
GSStaticSolver.solve(eq=eq,
profiles=profiles,
constrain=None,
target_relative_tolerance=1e-9)
# updates the plasma_psi (for later on)
eq._updatePlasmaPsi(eq.plasma_psi)
# plot the resulting equilbria
fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)
ax1.grid(True, which='both')
eq.plot(axis=ax1, show=False)
eq.tokamak.plot(axis=ax1, show=False)
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 --> built from pickle file.
Resistance (R) and inductance (M) matrices --> built using actives (and passives if present).
Tokamak built.
Forward static solve SUCCESS. Tolerance 4.87e-10 (vs. requested 1.00e-09) reached in 25/100 iterations.
Using the probe objects
The tokamak object has a probe object attribute pre-initialised, however, if we wanted to we could create a stand alone one by importing freegsnke.magnetic_probes and then setting up with probes = magnetic_probes.Probe().
In the tokamak, the probes object is located in tokamak.probes. When initialised it reads the information on the names, positions and orientations of the different probes from the pickle file.
We need to initialise the Greens functions that are used in the calculations. This is done by running tokamak.probes.initialise_setup(eq). This takes an input equilibrium object and saves probe attributes for each current source (i.e. the coils and the plasma) and evaluates the Greens functions at the positions of the probes. The purpose of the equilibrium here is to provide the size/shape of the grid that is used when determining the plasma Greens functions.
tokamak.probes.initialise_setup(eq)
Once initialised we can also now access information from the different probes:
flux loops measure \(\psi(r,z)\).
pickup coils measure \(\vec B(r,z)\cdot \hat n\) where \(\hat n\) is the orientation vector of the pickup coil.
# display the first five flux loop info
tokamak.probes.floops[0:5]
#print(tokamak.probes.floop_order[:5])
#print(tokamak.probes.floop_pos[:5])
[{'name': 'f_nu_01', 'position': array([0.901 , 1.3582], dtype=float32)},
{'name': 'f_nu_02', 'position': array([0.9544 , 1.3041999], dtype=float32)},
{'name': 'f_nu_03', 'position': array([1.0445 , 1.2150999], dtype=float32)},
{'name': 'f_nu_04', 'position': array([1.1239, 1.1366], dtype=float32)},
{'name': 'f_nu_a05', 'position': array([1.1505, 1.1987], dtype=float32)}]
# display the first five pickup coils info
tokamak.probes.pickups[0:5]
#print(tokamak.probes.floop_order[:5])
#print(tokamak.probes.floop_pos[:5])
[{'name': 'b_c1_p01',
'position': array([2.76900e-01, 3.00000e+02, 1.26203e+00], dtype=float32),
'orientation': 'PARALLEL',
'orientation_vector': array([0., 0., 1.], dtype=float32)},
{'name': 'b_c1_t02',
'position': array([2.7689108e-01, 2.9999680e+02, 1.2245095e+00], dtype=float32),
'orientation': 'TOROIDAL',
'orientation_vector': array([0., 1., 0.], dtype=float32)},
{'name': 'b_c1_p03',
'position': array([2.7689999e-01, 3.0000000e+02, 1.1870301e+00], dtype=float32),
'orientation': 'PARALLEL',
'orientation_vector': array([0., 0., 1.], dtype=float32)},
{'name': 'b_c1_t04',
'position': array([2.7689108e-01, 2.9999680e+02, 1.1495094e+00], dtype=float32),
'orientation': 'TOROIDAL',
'orientation_vector': array([0., 1., 0.], dtype=float32)},
{'name': 'b_c1_p05',
'position': array([2.76900e-01, 3.00000e+02, 1.11203e+00], dtype=float32),
'orientation': 'PARALLEL',
'orientation_vector': array([0., 0., 1.], dtype=float32)}]
Now in principle we could update and re-solve for the equilibrium. Doing this will not require any changes to the probe object, assuming the machine setup doesn’t change.
Once we have the equilibrium we want to analyse with the probes, we can call the probe functions calculate_fluxloop_value(eq) and calculate_pickup_value(eq), the outputs of which are arrays with probe values for each probe.
If one is interested in certain probes, then tokamak.probes.floop_order and tokamak.probes.pickup_order contain a list of the probe names which can be used to find the appropriate element of the output list. Alternatively the could be combined into a dictionary.
# compute probe values from equilibrium
floops_vals = tokamak.probes.calculate_fluxloop_value(eq)
pickups_vals = tokamak.probes.calculate_pickup_value(eq)
# create dictionary to access specific values (show here for the flux loops)
dict = {}
for i, el in enumerate(tokamak.probes.floop_order):
dict[el] = floops_vals[i]
dict
{'f_nu_01': 0.019508025863511077,
'f_nu_02': 0.016389267393185993,
'f_nu_03': 0.010517957501648606,
'f_nu_04': 0.004385540352841895,
'f_nu_a05': 0.0007045269159035722,
'f_nu_b05': 0.00016332281280780764,
'f_nu_06': -0.003421468914553631,
'f_nl_01': 0.01950022938981641,
'f_nl_02': 0.016372744616121655,
'f_nl_03': 0.01046834261776642,
'f_nl_04': 0.004290838141890491,
'f_nl_a05': 0.0006751691687334679,
'f_nl_b05': 0.00013643292544146535,
'f_nl_06': -0.0034825015085262895,
'f_p6u_01': -0.0038238514031687515,
'f_p6l_01': -0.0047346024057185465,
'f_bu_01': 0.015279126013471485,
'f_bu_02': 0.006688823459354964,
'f_bu_03': -0.00029955594545197015,
'f_bu_04': 0.0025362916706833695,
'f_bl_01': 0.015453431199033734,
'f_bl_02': 0.007009050674338595,
'f_bl_03': 0.00019506834465508188,
'f_bl_04': 0.003208642389300239,
'f_p5u_01': -0.029312355345764915,
'f_p5u_02': -0.033567535368408674,
'f_p5u_03': -0.03169037112097868,
'f_p5u_04': -0.023722660275558013,
'f_p5l_04': -0.025320947639182786,
'f_p5l_03': -0.03247944695592157,
'f_p5l_02': -0.03323452115618375,
'f_p5l_01': -0.028284085800739087,
'f_xu_01': 0.03503783782340049,
'f_xu_02': 0.0319585366310251,
'f_xu_03': 0.03214410781055678,
'f_xu_04': 0.02971911585010223,
'f_xu_05': 0.026676067716392513,
'f_xu_06': 0.026280983305174352,
'f_xu_07': 0.03196680693104395,
'f_xl_01': 0.03488006710738594,
'f_xl_02': 0.032038989091360055,
'f_xl_03': 0.03236492776456409,
'f_xl_04': 0.029572970476739503,
'f_xl_05': 0.026552573157495745,
'f_xl_06': 0.02621351733970166,
'f_xl_07': 0.0319961966261132,
'f_dpu_02': 0.022192511619817907,
'f_dpu_01': 0.021263909270750447,
'f_dpu_03': 0.01593546486978594,
'f_dpu_04': 0.014703129391506605,
'f_dpl_04': 0.01473294944496102,
'f_dpl_01': 0.02127759812670344,
'f_dpl_03': 0.01612368703328661,
'f_dpl_02': 0.022343856039771307,
'f_d3u_02': 0.03289853664373493,
'f_d3u_03': 0.03429489769898719,
'f_d3u_01': 0.03445142591612482,
'f_d3u_04': 0.034878303279153686,
'f_d3l_01': 0.03457521572338677,
'f_d3l_02': 0.035150125476364896,
'f_d3l_04': 0.032240269249439964,
'f_d3l_03': 0.03371180773191461,
'f_d2u_03': 0.03780633953207976,
'f_d2u_02': 0.035952632154454406,
'f_d2u_04': 0.039659877446915584,
'f_d2u_01': 0.04009138428591663,
'f_d2l_02': 0.039710079500917365,
'f_d2l_01': 0.03964557498882192,
'f_d2l_03': 0.03702123375680339,
'f_d2l_04': 0.034865222025403785,
'f_d1u_03': 0.046648012415156426,
'f_d1u_02': 0.03879872664795093,
'f_d1u_04': 0.04923450859781968,
'f_d1u_01': 0.043071468751650674,
'f_d1l_04': 0.0518039581863362,
'f_d1l_01': 0.045079475304175,
'f_d1l_03': 0.04425442052951883,
'f_d1l_02': 0.03654158328579665,
'f_c_a01': 0.029908808583542248,
'f_c_b01': 0.029897203786296348,
'f_c_a02': 0.026517114102207216,
'f_c_b02': 0.026466050119017228,
'f_c_a03': 0.0195791118666565,
'f_c_b03': 0.019559218093552933,
'f_c_a04': 0.018081942240838467,
'f_c_b04': 0.01809028158127342,
'f_c_a05': 0.018951618736368463,
'f_c_b05': 0.018961331313067126,
'f_c_a06': 0.019840130690278256,
'f_c_b06': 0.01984779951600365,
'f_c_a07': 0.01978945134710404,
'f_c_b07': 0.019781526118287547,
'f_c_a08': 0.01888842296993975,
'f_c_b08': 0.018878687676683888,
'f_c_a09': 0.01802876135213343,
'f_c_b09': 0.01802070914806749,
'f_c_a10': 0.019718196738282828,
'f_c_b10': 0.019741170845805116,
'f_c_a11': 0.026832787239055328,
'f_c_b11': 0.02687902122742693,
'f_c_a12': 0.02997855858211923,
'f_c_b12': 0.02998835182001161}
Suppose we want to compute a new equilibrium with a different grid spacing or shape. We don’t need to update the probe objects, we simply pass the new equilibrium to the ‘calculate’ functions. The first time a new grid is encountered it will create a new set of greens functions and save them to a dictionary so that they can be reused in the future if the same grid is used again.
Below is a new equilibrium with a modified grid shape and spacing. When a new grid is encountered, a message is displayed to tell that new greens are computed. Note it only does it the first time.
Note that computing on a different grid but with same plasma setup should give the same values at the probes (which it does).
# initialise the equilibrium
eq_new = equilibrium_update.Equilibrium(
tokamak=tokamak,
Rmin=0.1, Rmax=2.0, # Radial range
Zmin=-2.0, Zmax=2.0, # 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
)
from freegsnke.jtor_update import ConstrainPaxisIp
profiles = ConstrainPaxisIp(
eq=eq_new,
paxis=8e3,
Ip=6e5,
fvac=0.5,
alpha_m=1.8,
alpha_n=1.2
)
from freegsnke import GSstaticsolver
GSStaticSolver = GSstaticsolver.NKGSsolver(eq_new)
import pickle
with open('simple_diverted_currents_PaxisIp.pk', 'rb') as f:
current_values = pickle.load(f)
for key in current_values.keys():
eq_new.tokamak.set_coil_current(key, current_values[key])
GSStaticSolver.solve(eq=eq_new,
profiles=profiles,
constrain=None,
target_relative_tolerance=1e-9)
Forward static solve SUCCESS. Tolerance 7.11e-10 (vs. requested 1.00e-09) reached in 25/100 iterations.
floops_vals_new = tokamak.probes.calculate_fluxloop_value(eq_new)
pickups_vals_new = tokamak.probes.calculate_pickup_value(eq_new)
new equilibrium grid - computed new greens functions
new equilibrium grid - computed new greens functions
# compare values
print(floops_vals[:5])
print(floops_vals_new[:5])
[0.01950803 0.01638927 0.01051796 0.00438554 0.00070453]
[0.01950725 0.01638834 0.01051674 0.00438404 0.00070317]
# compare values
print(pickups_vals[:5])
print(pickups_vals_new[:5])
[ 2.22001732e-02 1.80576420e+00 -6.77851406e-04 1.80576420e+00
-1.02435028e-02]
[ 2.21979295e-02 1.80576420e+00 -6.80208110e-04 1.80576420e+00
-1.02459068e-02]
If we re-run this same line of code, we don’t get the message that the greens functions have been recalculated. They are stored in a dictionary labeled by a key containing the grid specification in the form key = (Rmin,Rmax,Zmin,Zmax,nx,ny).
pickup_vals_new2 = tokamak.probes.calculate_pickup_value(eq_new)
# show greens function keys
tokamak.probes.greens_B_plasma_oriented.keys()
dict_keys([(0.1, 2.0, -2.2, 2.2, 65, 129), (0.1, 2.0, -2.0, 2.0, 65, 129)])
We can also plot the fluxloop and the pickup coil locations (and orientations) on our machine model.
# plot the resulting equilbria
fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)
ax1.grid(True, which='both')
# eq.plot(axis=ax1, show=False)
eq.tokamak.plot(axis=ax1, show=False)
eq.tokamak.probes.plot(axis=ax1, show=False, floops=True, pickups=True, pickups_scale=0.05)
ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle="--")
ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle="-")
ax1.set_xlim(0.1, 2.15)
ax1.set_ylim(-2.25, 2.25)
plt.tight_layout()
