diff --git a/documentation/physics-models/plasma_beta/plasma_beta.md b/documentation/physics-models/plasma_beta/plasma_beta.md
index d49bff6c36..1ecc214ae0 100644
--- a/documentation/physics-models/plasma_beta/plasma_beta.md
+++ b/documentation/physics-models/plasma_beta/plasma_beta.md
@@ -5,9 +5,19 @@
The efficiency of confinement of plasma pressure by the magnetic field is represented by the ratio:
$$
-\beta = \frac{2\mu_0p}{B^2}
+\beta(\rho) = \frac{2\mu_0p(\rho)}{\left(B(\rho)\right)^2}
$$
+Where $\beta$ is a function of normalised minor radius across the plasma ($\rho$), due to the change in pressure and magnetic field strength.
+
+The standard $\beta$ term used for comparison and to represent the plasma as a whole in many calculations is the volume averaged value given by:
+
+$$
+\langle \beta \rangle = \frac{2\mu_0 \langle p \rangle}{\langle B \rangle^2}
+$$
+
+Where $\langle p \rangle$ is the volume averaged plasma pressure and $\langle B \rangle$ is the average field.
+
There are several different measures of this type, arising from different choices of definition and from the need to quantify different equilibrium properties.
In its expanded form of magnetic field components:
@@ -92,6 +102,20 @@ $$
------------------------
+## Definitions
+
+### Volume averaged thermal toroidal beta
+
+$$
+\overbrace{\langle \beta_t \rangle}^{\texttt{beta_toroidal_thermal_vol_avg}} = \frac{2\mu_0 \langle p_{\text{thermal}} \rangle}{\underbrace{B_{\text{T,on-axis}}^2}_{\texttt{b_plasma_toroidal_on_axis}}}
+$$
+
+### Volume averaged thermal poloidal beta
+
+$$
+\langle \beta_p \rangle = \frac{2\mu_0 \langle p_{\text{thermal}} \rangle}{\langle B_{\text{P,average}} \rangle^2}
+$$
+
## Troyon Beta Limit
The Troyon plasma beta limit is given by[^0][^1]:
diff --git a/examples/data/large_tokamak_nof_2_MFILE.DAT b/examples/data/large_tokamak_nof_2_MFILE.DAT
index e1d77b2528..8920191321 100644
--- a/examples/data/large_tokamak_nof_2_MFILE.DAT
+++ b/examples/data/large_tokamak_nof_2_MFILE.DAT
@@ -231,8 +231,8 @@
Volume_averaged_electron_number_density_(/m3)____________________________ (nd_plasma_electrons_vol_avg)_________________________ 7.79622390002983731e+19
Electron_number_density_on_axis_(/m3)____________________________________ (nd_plasma_electron_on_axis)__________________________ 1.02466460344915771e+20 OP
Line-averaged_electron_number_density_(/m3)______________________________ (nd_electron_line)_________________________ 8.63127495812675338e+19 OP
- Plasma_pressure_on_axis_(Pa)_____________________________________________ (pres_plasma_on_axis)___________________________ 7.82668759018785902e+05 OP
- Volume_averaged_plasma_pressure_(Pa)_____________________________________ (pres_plasma_vol_avg)_____________ 2.26860509860517632e+05 OP
+ Plasma_pressure_on_axis_(Pa)_____________________________________________ (pres_plasma_thermal_on_axis)___________________________ 7.82668759018785902e+05 OP
+ Volume_averaged_plasma_pressure_(Pa)_____________________________________ (pres_plasma_thermal_vol_avg)_____________ 2.26860509860517632e+05 OP
Line-averaged_electron_density_/_Greenwald_density_______________________ (dnla_gw)______________________ 1.19833068191158598e+00 OP
Total_Ion_number_density_(/m3)___________________________________________ (nd_ions_total)________________ 7.06084473006886912e+19 OP
Fuel_ion_number_density_(/m3)____________________________________________ (nd_fuel_ions)_________________ 6.58453241268011663e+19 OP
diff --git a/examples/data/large_tokamak_nof_MFILE.DAT b/examples/data/large_tokamak_nof_MFILE.DAT
index cdd794dd89..0e3791a63a 100644
--- a/examples/data/large_tokamak_nof_MFILE.DAT
+++ b/examples/data/large_tokamak_nof_MFILE.DAT
@@ -303,8 +303,8 @@
Volume_averaged_electron_number_density_(/m3)____________________________ (nd_plasma_electrons_vol_avg)_________________________ 8.50078961768999158e+19
Electron_number_density_on_axis_(/m3)____________________________________ (nd_plasma_electron_on_axis)__________________________ 1.11814263037413868e+20 OP
Line-averaged_electron_number_density_(/m3)______________________________ (nd_electron_line)_________________________ 9.41418090734859715e+19 OP
- Plasma_pressure_on_axis_(Pa)_____________________________________________ (pres_plasma_on_axis)___________________________ 9.10652371594472905e+05 OP
- Volume_averaged_plasma_pressure_(Pa)_____________________________________ (pres_plasma_vol_avg)_____________ 2.63957209157818230e+05 OP
+ Plasma_pressure_on_axis_(Pa)_____________________________________________ (pres_plasma_thermal_on_axis)___________________________ 9.10652371594472905e+05 OP
+ Volume_averaged_plasma_pressure_(Pa)_____________________________________ (pres_plasma_thermal_vol_avg)_____________ 2.63957209157818230e+05 OP
Line-averaged_electron_density_/_Greenwald_density_______________________ (dnla_gw)______________________ 1.19999999999999951e+00 OP
Total_Ion_number_density_(/m3)___________________________________________ (nd_ions_total)________________ 7.34892735750798868e+19 OP
Fuel_ion_number_density_(/m3)____________________________________________ (nd_fuel_ions)_________________ 6.49714589130413015e+19 OP
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index c30155c8a4..a2ab4440a6 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -244,6 +244,14 @@
beta_toroidal: float = None
"""toroidal beta"""
+beta_toroidal_profile: list[float] = None
+"""toroidal beta profile"""
+
+beta_poloidal_profile: list[float] = None
+"""poloidal beta profile"""
+
+beta_total_profile: list[float] = None
+"""total beta profile"""
beta_thermal: float = None
"""thermal beta"""
@@ -284,10 +292,24 @@
b_plasma_poloidal_average: float = None
"""Plasma average poloidal field (T)"""
+b_plasma_poloidal_edge: float = None
+"""Plasma poloidal field at plasma edge (T)"""
b_plasma_toroidal_on_axis: float = None
"""Plasma toroidal field on axis (T) (`iteration variable 2`)"""
+b_plasma_toroidal_profile: list[float] = None
+"""toroidal field profile in plasma (T)"""
+
+b_plasma_poloidal_profile: list[float] = None
+"""poloidal field profile in plasma (T)"""
+
+b_plasma_circular_poloidal_profile: list[float] = None
+"""poloidal field profile in circular plasma (T)"""
+
+b_plasma_total_profile: list[float] = None
+"""total field profile in plasma (T)"""
+
b_plasma_total: float = None
"""Sum of plasma total toroidal + poloidal field (T)"""
@@ -803,7 +825,7 @@
"""margin to vertical stability"""
-pres_plasma_on_axis: float = None
+pres_plasma_thermal_on_axis: float = None
"""central total plasma pressure (Pa)"""
pres_plasma_total_profile: list[float] = None
@@ -821,11 +843,14 @@
j_plasma_on_axis: float = None
"""Central plasma current density (A/m2)"""
+j_plasma_circular_on_axis: float = None
+"""Central plasma current density for a circular plasma (A/m2)"""
+
n_plasma_profile_elements: int = None
"""Number of elements in plasma profile"""
-pres_plasma_vol_avg: float = None
-"""Volume averaged plasma pressure (Pa)"""
+pres_plasma_thermal_vol_avg: float = None
+"""Volume averaged thermal plasma pressure (Pa)"""
f_dd_branching_trit: float = None
@@ -951,6 +976,9 @@
plasma_current: float = None
"""plasma current (A)"""
+c_plasma_circular: float = None
+"""Plasma current for circular plasma (A)"""
+
p_plasma_neutron_mw: float = None
"""Neutron fusion power from just the plasma [MW]"""
@@ -1097,6 +1125,11 @@
qstar: float = None
"""cylindrical safety factor"""
+q_95_circular: float = None
+"""Safety factor at 95% flux surface for circular plasma"""
+
+q_circular_profile: list[float] = None
+"""safety factor profile for circular plasma"""
rad_fraction_sol: float = None
"""SoL radiation fraction"""
@@ -1380,6 +1413,9 @@ def init_physics_variables():
global beta_poloidal
global beta_poloidal_eps
global beta_toroidal
+ global beta_toroidal_profile
+ global beta_total_profile
+ global beta_poloidal_profile
global beta_thermal
global beta_thermal_poloidal
global beta_thermal_toroidal
@@ -1390,7 +1426,12 @@ def init_physics_variables():
global beta_norm_toroidal
global betbm0
global b_plasma_poloidal_average
+ global b_plasma_poloidal_edge
global b_plasma_toroidal_on_axis
+ global b_plasma_toroidal_profile
+ global b_plasma_poloidal_profile
+ global b_plasma_circular_poloidal_profile
+ global b_plasma_total_profile
global b_plasma_total
global burnup
global burnup_in
@@ -1488,12 +1529,13 @@ def init_physics_variables():
global nd_plasma_electron_on_axis
global nd_plasma_ions_on_axis
global m_s_limit
- global pres_plasma_on_axis
+ global pres_plasma_thermal_on_axis
global pres_plasma_total_profile
global pres_plasma_electron_profile
global pres_plasma_ion_total_profile
global pres_plasma_fuel_profile
global j_plasma_on_axis
+ global j_plasma_circular_on_axis
global n_plasma_profile_elements
global f_dd_branching_trit
global pden_plasma_alpha_mw
@@ -1526,6 +1568,7 @@ def init_physics_variables():
global pflux_fw_rad_mw
global pden_ion_electron_equilibration_mw
global plasma_current
+ global c_plasma_circular
global p_plasma_neutron_mw
global p_neutron_total_mw
global pden_neutron_total_mw
@@ -1556,6 +1599,8 @@ def init_physics_variables():
global tauratio
global q95_min
global qstar
+ global q_95_circular
+ global q_circular_profile
global rad_fraction_sol
global rad_fraction_total
global f_nd_alpha_electron
@@ -1632,6 +1677,9 @@ def init_physics_variables():
beta_poloidal = 0.0
beta_poloidal_eps = 0.0
beta_toroidal = 0.0
+ beta_toroidal_profile = []
+ beta_poloidal_profile = []
+ beta_total_profile = []
beta_thermal = 0.0
beta_thermal_poloidal = 0.0
beta_thermal_toroidal = 0.0
@@ -1642,7 +1690,12 @@ def init_physics_variables():
beta_norm_toroidal = 0.0
betbm0 = 1.5
b_plasma_poloidal_average = 0.0
+ b_plasma_poloidal_edge = 0.0
b_plasma_toroidal_on_axis = 5.68
+ b_plasma_toroidal_profile = []
+ b_plasma_poloidal_profile = []
+ b_plasma_circular_poloidal_profile = []
+ b_plasma_total_profile = []
b_plasma_total = 0.0
burnup = 0.0
burnup_in = 0.0
@@ -1740,12 +1793,13 @@ def init_physics_variables():
nd_plasma_electron_on_axis = 0.0
nd_plasma_ions_on_axis = 0.0
m_s_limit = 0.3
- pres_plasma_on_axis = 0.0
+ pres_plasma_thermal_on_axis = 0.0
pres_plasma_total_profile = []
pres_plasma_electron_profile = []
pres_plasma_ion_total_profile = []
pres_plasma_fuel_profile = []
j_plasma_on_axis = 0.0
+ j_plasma_circular_on_axis = 0.0
n_plasma_profile_elements = 500
f_dd_branching_trit = 0.0
pden_plasma_alpha_mw = 0.0
@@ -1778,6 +1832,7 @@ def init_physics_variables():
pflux_fw_rad_mw = 0.0
pden_ion_electron_equilibration_mw = 0.0
plasma_current = 0.0
+ c_plasma_circular = 0.0
p_plasma_neutron_mw = 0.0
p_neutron_total_mw = 0.0
pden_neutron_total_mw = 0.0
@@ -1808,6 +1863,8 @@ def init_physics_variables():
tauratio = 1.0
q95_min = 0.0
qstar = 0.0
+ q_95_circular = 0.0
+ q_circular_profile = []
rad_fraction_sol = 0.8
rad_fraction_total = 0.0
f_nd_alpha_electron = 0.1
diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
index e149229d3d..920fd79fd1 100644
--- a/process/io/plot_proc.py
+++ b/process/io/plot_proc.py
@@ -3706,40 +3706,48 @@ def plot_n_profiles(prof, demo_ranges, mfile_data, scan):
# ---
-def plot_jprofile(prof):
+def plot_jprofile(axis, mfile_data, scan):
"""Function to plot density profile
Arguments:
prof --> axis object to add plot to
"""
- prof.set_xlabel(r"$\rho \quad [r/a]$")
- prof.set_ylabel(r"Current density $[kA/m^2]$")
- prof.set_title("$J$ profile")
- prof.minorticks_on()
- prof.set_xlim([0, 1.0])
+ axis.set_xlabel(r"$\rho \quad [r/a]$")
+ axis.set_ylabel(r"Current density $[kA/m^2]$")
+ axis.set_title("$J$ profile")
+ axis.minorticks_on()
+ axis.set_xlim([0, 1.0])
+
+ j_plasma_circular_on_axis = mfile_data.data["j_plasma_circular_on_axis"].get_scan(
+ scan
+ )
rho = np.linspace(0, 1)
y2 = (j_plasma_0 * (1 - rho**2) ** alphaj) / 1e3
+ y_circular = (j_plasma_circular_on_axis * (1 - rho**2) ** alphaj) / 1e3
- prof.plot(rho, y2, label="$n_{i}$", color="red")
+ axis.plot(rho, y2, label="Real", color="red")
+ axis.plot(rho, y_circular, label="Circular", color="blue", linestyle="--")
+ axis.legend()
textstr_j = "\n".join((
- r"$j_0$: " + f"{y2[0]:.3f} kA m$^{{-2}}$\n",
+ r"$j_0$: " + f"{j_plasma_0 / 1000:.3f} kA m$^{{-2}}$\n",
+ r"$j_0,circular$: " + f"{j_plasma_circular_on_axis / 1000:.3f} kA m$^{{-2}}$\n",
r"$\alpha_J$: " + f"{alphaj:.3f}",
))
props_j = {"boxstyle": "round", "facecolor": "wheat", "alpha": 0.5}
- prof.text(
+ axis.text(
1.1,
0.75,
textstr_j,
- transform=prof.transAxes,
+ transform=axis.transAxes,
fontsize=9,
verticalalignment="top",
bbox=props_j,
)
- prof.text(
+ axis.text(
0.05,
0.04,
"*Current profile is assumed to be parabolic",
@@ -3747,7 +3755,7 @@ def plot_jprofile(prof):
ha="left",
transform=plt.gcf().transFigure,
)
- prof.grid(True, which="both", linestyle="--", linewidth=0.5, alpha=0.2)
+ axis.grid(True, which="both", linestyle="--", linewidth=0.5, alpha=0.2)
def plot_t_profiles(prof, demo_ranges, mfile_data, scan):
@@ -3860,9 +3868,27 @@ def plot_qprofile(prof, demo_ranges, mfile_data, scan):
q_r_nevin = q0 + (q95 - q0) * (rho + rho * rho + rho**3) / (3.0)
q_r_sauter = q0 + (q95 - q0) * (rho * rho)
+ n_plasma_profile_elements = int(
+ mfile_data.data["n_plasma_profile_elements"].get_scan(scan)
+ )
+
+ q_circular_profile = [
+ mfile_data.data[f"q_circular_profile{i}"].get_scan(scan)
+ for i in range(n_plasma_profile_elements)
+ ]
+
+ prof.plot(
+ np.linspace(0, 1, n_plasma_profile_elements),
+ q_circular_profile,
+ label="Circular",
+ linestyle="--",
+ color="black",
+ )
+
prof.plot(rho, q_r_nevin, label="Nevins")
prof.plot(rho, q_r_sauter, label="Sauter")
prof.legend()
+ prof.figure.tight_layout()
# Ranges
# ---
@@ -3876,9 +3902,9 @@ def plot_qprofile(prof, demo_ranges, mfile_data, scan):
prof.set_ylim([0, q95 * 1.2])
prof.text(
- 0.6,
+ 0.45,
0.04,
- "*Profile is not calculated, only $q_0$ and $q_{95}$ are known.",
+ "*Only circular profile is calculated, for others only $q_0$ and $q_{95}$ are known.",
fontsize=10,
ha="left",
transform=plt.gcf().transFigure,
@@ -3889,12 +3915,14 @@ def plot_qprofile(prof, demo_ranges, mfile_data, scan):
textstr_q = "\n".join((
r"$q_0$: " + f"{q0:.3f}\n",
r"$q_{95}$: " + f"{q95:.3f}\n",
- r"$q_{\text{cyl}}$: " + f"{mfile_data.data['qstar'].get_scan(scan):.3f}",
+ r"$q_{\text{cyl}}$: " + f"{mfile_data.data['qstar'].get_scan(scan):.3f}\n",
+ r"$q_{95,\text{circular}}$: "
+ + f"{mfile_data.data['q_95_circular'].get_scan(scan):.3f}",
))
props_q = {"boxstyle": "round", "facecolor": "wheat", "alpha": 0.5}
prof.text(
- -0.4,
+ -0.5,
0.75,
textstr_q,
transform=prof.transAxes,
@@ -10695,6 +10723,173 @@ def plot_plasma_poloidal_pressure_contours(
axis.set_title("Plasma Poloidal Pressure Contours")
+def plot_magnetic_fields_in_plasma(axis, mfile_data, scan):
+ # Plot magnetic field profiles inside the plasma boundary
+
+ n_plasma_profile_elements = int(
+ mfile_data.data["n_plasma_profile_elements"].get_scan(scan)
+ )
+
+ # Get toroidal magnetic field profile (in Tesla)
+ b_plasma_toroidal_profile = [
+ mfile_data.data[f"b_plasma_toroidal_profile{i}"].get_scan(scan)
+ for i in range(2 * n_plasma_profile_elements)
+ ]
+
+ b_plasma_poloidal_profile = [
+ mfile_data.data[f"b_plasma_poloidal_profile{i}"].get_scan(scan)
+ for i in range(2 * n_plasma_profile_elements)
+ ]
+
+ b_plasma_circular_poloidal_profile = [
+ mfile_data.data[f"b_plasma_circular_poloidal_profile{i}"].get_scan(scan)
+ for i in range(2 * n_plasma_profile_elements)
+ ]
+
+ b_plasma_total_profile = [
+ mfile_data.data[f"b_plasma_total_profile{i}"].get_scan(scan)
+ for i in range(2 * n_plasma_profile_elements)
+ ]
+
+ # Get major and minor radius for x-axis in metres
+ rmajor = mfile_data.data["rmajor"].get_scan(scan)
+ rminor = mfile_data.data["rminor"].get_scan(scan)
+
+ # Plot magnetic field first (background)
+ axis.plot(
+ np.linspace(rmajor - rminor, rmajor + rminor, len(b_plasma_toroidal_profile)),
+ b_plasma_toroidal_profile,
+ color="blue",
+ label="Toroidal B-field [T]",
+ )
+
+ axis.plot(
+ np.linspace(rmajor - rminor, rmajor + rminor, len(b_plasma_poloidal_profile)),
+ b_plasma_poloidal_profile,
+ color="red",
+ label="Poloidal B-field [T]",
+ )
+
+ axis.plot(
+ np.linspace(
+ rmajor - rminor, rmajor + rminor, len(b_plasma_circular_poloidal_profile)
+ ),
+ b_plasma_circular_poloidal_profile,
+ color="orange",
+ label="Circular Poloidal B-field [T]",
+ )
+
+ axis.plot(
+ np.linspace(rmajor - rminor, rmajor + rminor, len(b_plasma_total_profile)),
+ b_plasma_total_profile,
+ color="green",
+ label="Total B-field [T]",
+ )
+
+ # Plot plasma minor radius as a circle (ensure it's visible and not covered)
+ circle = plt.Circle(
+ (rmajor, 0),
+ rminor,
+ color="purple",
+ fill=False,
+ linewidth=2,
+ linestyle="--",
+ label="Plasma Outer Circlular Flux Surface",
+ zorder=10,
+ )
+ axis.add_patch(circle)
+
+ # Plot plasma on top of magnetic field, displaced vertically by bt
+ plot_plasma(axis, mfile_data, scan, colour_scheme=1)
+
+ # Plot plasma centre dot
+ axis.plot(rmajor, 0, marker="o", color="red", markersize=8, label="Plasma Centre")
+
+ # Plot vertical lines at plasma edge
+ axis.axvline(
+ rmajor - rminor, color="green", linestyle="--", linewidth=1, label="Plasma Edge"
+ )
+ axis.axvline(rmajor + rminor, color="green", linestyle="--", linewidth=1)
+ axis.set_xlabel("Radial Position [m]")
+ axis.set_ylabel("Toroidal Magnetic Field [T]")
+ axis.set_title("Toroidal Magnetic Field Profile in Plasma")
+ axis.grid(True, linestyle="--", alpha=0.5)
+ axis.minorticks_on()
+ axis.legend()
+ axis.set_xlim(rmajor - 1.25 * rminor, rmajor + 1.25 * rminor)
+
+
+def plot_beta_profiles(axis, mfile_data, scan):
+ # Plot the beta profiles on the given axis
+
+ n_plasma_profile_elements = int(
+ mfile_data.data["n_plasma_profile_elements"].get_scan(scan)
+ )
+
+ beta_plasma_toroidal_profile = [
+ mfile_data.data[f"beta_toroidal_profile{i}"].get_scan(scan)
+ for i in range(2 * n_plasma_profile_elements)
+ ]
+
+ beta_total_profile = [
+ mfile_data.data[f"beta_total_profile{i}"].get_scan(scan)
+ for i in range(2 * n_plasma_profile_elements)
+ ]
+
+ axis.plot(
+ np.linspace(-1, 1, 2 * n_plasma_profile_elements),
+ beta_plasma_toroidal_profile,
+ color="blue",
+ label="Beta Toroidal",
+ )
+
+ axis.plot(
+ np.linspace(-1, 1, 2 * n_plasma_profile_elements),
+ beta_total_profile,
+ color="green",
+ label="Beta Total",
+ )
+
+ axis.set_xlabel("$\\rho$ [r/a]")
+ axis.set_ylabel("$\\beta$ [%]")
+ axis.minorticks_on()
+ axis.grid(which="minor", linestyle=":", linewidth=0.5, alpha=0.5)
+ axis.set_title("Beta Profile")
+ axis.legend()
+ axis.axvline(x=0, color="black", linestyle="--", linewidth=1)
+ axis.grid(True, linestyle="--", alpha=0.5)
+
+
+def plot_beta_profiles_poloidal(axis, mfile_data, scan):
+ # Plot the beta profiles on the given axis
+
+ n_plasma_profile_elements = int(
+ mfile_data.data["n_plasma_profile_elements"].get_scan(scan)
+ )
+
+ beta_poloidal_profile = [
+ mfile_data.data[f"beta_poloidal_profile{i}"].get_scan(scan)
+ for i in range(2 * n_plasma_profile_elements)
+ ]
+
+ axis.plot(
+ np.linspace(-1, 1, 2 * n_plasma_profile_elements),
+ beta_poloidal_profile,
+ color="red",
+ label="Beta Poloidal",
+ )
+
+ axis.set_xlabel("$\\rho$ [r/a]")
+ axis.set_ylabel("$\\beta$ [%]")
+ axis.minorticks_on()
+ axis.grid(which="minor", linestyle=":", linewidth=0.5, alpha=0.5)
+ axis.set_title("Beta Profile")
+ axis.legend()
+ axis.set_yscale("log")
+ axis.axvline(x=0, color="black", linestyle="--", linewidth=1)
+ axis.grid(True, linestyle="--", alpha=0.5)
+
+
def main_plot(
fig0,
fig1,
@@ -10716,6 +10911,7 @@ def main_plot(
fig17,
fig18,
fig19,
+ fig20,
m_file_data,
scan,
imp="../data/lz_non_corona_14_elements/",
@@ -10800,7 +10996,7 @@ def main_plot(
# Plot current density profile
ax12 = fig4.add_subplot(4, 3, 10)
ax12.set_position([0.075, 0.125, 0.25, 0.15])
- plot_jprofile(ax12)
+ plot_jprofile(ax12, mfile_data=m_file_data, scan=scan)
# Plot q profile
ax13 = fig4.add_subplot(4, 3, 12)
@@ -10901,6 +11097,13 @@ def main_plot(
fig19.add_subplot(111, aspect="equal"), m_file_data, scan, fig19
)
+ plot_magnetic_fields_in_plasma(
+ fig20.add_subplot(122, aspect="equal"), m_file_data, scan
+ )
+
+ plot_beta_profiles(fig20.add_subplot(221), m_file_data, scan)
+ plot_beta_profiles_poloidal(fig20.add_subplot(223), m_file_data, scan)
+
def main(args=None):
# TODO The use of globals here isn't ideal, but is required to get main()
@@ -11210,6 +11413,7 @@ def main(args=None):
page17 = plt.figure(figsize=(12, 9), dpi=80)
page18 = plt.figure(figsize=(12, 9), dpi=80)
page19 = plt.figure(figsize=(12, 9), dpi=80)
+ page20 = plt.figure(figsize=(12, 9), dpi=80)
# run main_plot
main_plot(
@@ -11233,6 +11437,7 @@ def main(args=None):
page17,
page18,
page19,
+ page20,
m_file,
scan=scan,
demo_ranges=demo_ranges,
@@ -11261,6 +11466,7 @@ def main(args=None):
pdf.savefig(page17)
pdf.savefig(page18)
pdf.savefig(page19)
+ pdf.savefig(page20)
# show fig if option used
if args.show:
@@ -11286,6 +11492,7 @@ def main(args=None):
plt.close(page17)
plt.close(page18)
plt.close(page19)
+ plt.close(page20)
if __name__ == "__main__":
diff --git a/process/physics.py b/process/physics.py
index ba0df63b78..97bac88907 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -1617,7 +1617,7 @@ def physics(self):
# Calculate plasma current
(
- physics_variables.b_plasma_poloidal_average,
+ _,
physics_variables.qstar,
physics_variables.plasma_current,
) = self.calculate_plasma_current(
@@ -1628,7 +1628,7 @@ def physics(self):
physics_variables.i_plasma_current,
physics_variables.kappa,
physics_variables.kappa95,
- physics_variables.pres_plasma_on_axis,
+ physics_variables.pres_plasma_thermal_on_axis,
physics_variables.len_plasma_poloidal,
physics_variables.q95,
physics_variables.rmajor,
@@ -1735,6 +1735,110 @@ def physics(self):
+ physics_variables.b_plasma_poloidal_average**2
)
+ # Calculate the toroidal field across the plasma
+ # Calculate the toroidal field profile across the plasma (1/R dependence)
+ # Double element size to include both sides of the plasma
+ rho = np.linspace(
+ physics_variables.rmajor - physics_variables.rminor,
+ physics_variables.rmajor + physics_variables.rminor,
+ 2 * physics_variables.n_plasma_profile_elements,
+ )
+ # Avoid division by zero at the magnetic axis
+ rho = np.where(rho == 0, 1e-10, rho)
+ physics_variables.b_plasma_toroidal_profile = (
+ physics_variables.rmajor * physics_variables.b_plasma_toroidal_on_axis / rho
+ )
+
+ physics_variables.b_plasma_poloidal_profile = (
+ self.calculate_poloidal_field_profile(
+ j_plasma_on_axis=physics_variables.j_plasma_on_axis,
+ alphaj=physics_variables.alphaj,
+ n_profile_elements=physics_variables.n_plasma_profile_elements * 2,
+ rminor=physics_variables.rminor,
+ )
+ )
+ # Calculate the average poloidal field from the profile
+ physics_variables.b_plasma_poloidal_average = np.mean(
+ physics_variables.b_plasma_poloidal_profile
+ )
+
+ physics_variables.b_plasma_circular_poloidal_profile = (
+ self.calculate_poloidal_field_profile(
+ j_plasma_on_axis=physics_variables.j_plasma_circular_on_axis,
+ alphaj=physics_variables.alphaj,
+ n_profile_elements=physics_variables.n_plasma_profile_elements * 2,
+ rminor=physics_variables.rminor,
+ )
+ )
+
+ physics_variables.b_plasma_circular_poloidal_profile = (
+ (
+ constants.RMU0
+ * physics_variables.j_plasma_circular_on_axis
+ * physics_variables.rminor**2
+ * (
+ 1
+ - (
+ 1
+ - (
+ np.linspace(
+ -1, 1, 2 * physics_variables.n_plasma_profile_elements
+ )
+ )
+ ** 2
+ )
+ ** (physics_variables.alphaj + 1)
+ )
+ )
+ / (2 * (physics_variables.alphaj + 1))
+ * np.linspace(-1, 1, 2 * physics_variables.n_plasma_profile_elements)
+ * physics_variables.rminor
+ )
+
+ # physics_variables.q_circular_profile = (
+ # (
+ # np.linspace(
+ # 0,
+ # physics_variables.rminor,
+ # physics_variables.n_plasma_profile_elements,
+ # )
+ # / physics_variables.rmajor
+ # )
+ # * physics_variables.b_plasma_toroidal_profile[
+ # physics_variables.n_plasma_profile_elements :
+ # ]
+ # / physics_variables.b_plasma_circular_poloidal_profile[
+ # physics_variables.n_plasma_profile_elements :
+ # ]
+ # )
+ rho = np.linspace(
+ 0,
+ 1.0,
+ physics_variables.n_plasma_profile_elements,
+ )
+ physics_variables.q_circular_profile = (
+ (
+ (
+ 2
+ * (physics_variables.alphaj + 1)
+ * physics_variables.b_plasma_toroidal_on_axis
+ )
+ / (
+ constants.RMU0
+ * physics_variables.rmajor
+ * physics_variables.j_plasma_circular_on_axis
+ )
+ )
+ * rho**2
+ / (1 - (1 - rho**2) ** (physics_variables.alphaj + 1))
+ )
+
+ # Calculate the total magnetic field profile at each point by summing the squares of the toroidal and poloidal components
+ physics_variables.b_plasma_total_profile = np.sqrt(
+ np.square(physics_variables.b_plasma_toroidal_profile)
+ + np.square(physics_variables.b_plasma_poloidal_profile)
+ )
+
# ============================================
# -----------------------------------------------------
@@ -1747,6 +1851,36 @@ def physics(self):
/ physics_variables.b_plasma_toroidal_on_axis**2
)
+ # Mirror the pressure profiles to match the doubled toroidal field profile
+ pres_profile_total = np.concatenate([
+ physics_variables.pres_plasma_total_profile[::-1],
+ physics_variables.pres_plasma_total_profile,
+ ])
+
+ physics_variables.beta_toroidal_profile = np.array([
+ self.calculate_plasma_beta(
+ pres_plasma=pres_profile_total[i],
+ b_field=physics_variables.b_plasma_toroidal_profile[i],
+ )
+ for i in range(len(physics_variables.b_plasma_toroidal_profile))
+ ])
+
+ physics_variables.beta_poloidal_profile = np.array([
+ self.calculate_plasma_beta(
+ pres_plasma=pres_profile_total[i],
+ b_field=physics_variables.b_plasma_poloidal_profile[i],
+ )
+ for i in range(len(physics_variables.b_plasma_poloidal_profile))
+ ])
+
+ physics_variables.beta_total_profile = np.array([
+ self.calculate_plasma_beta(
+ pres_plasma=pres_profile_total[i],
+ b_field=physics_variables.b_plasma_total_profile[i],
+ )
+ for i in range(len(physics_variables.b_plasma_total_profile))
+ ])
+
# Calculate physics_variables.beta poloidal [-]
physics_variables.beta_poloidal = calculate_poloidal_beta(
physics_variables.b_plasma_total,
@@ -2011,8 +2145,8 @@ def physics(self):
current_drive_variables.cboot
* self.bootstrap_fraction_andrade(
beta_poloidal=physics_variables.beta_poloidal,
- core_pressure=physics_variables.pres_plasma_on_axis,
- average_pressure=physics_variables.pres_plasma_vol_avg,
+ core_pressure=physics_variables.pres_plasma_thermal_on_axis,
+ average_pressure=physics_variables.pres_plasma_thermal_vol_avg,
inverse_aspect=physics_variables.eps,
)
)
@@ -2602,8 +2736,8 @@ def physics(self):
physics_variables.beta_norm_max_thloreus = (
self.calculate_beta_norm_max_thloreus(
c_beta=physics_variables.c_beta,
- pres_plasma_on_axis=physics_variables.pres_plasma_on_axis,
- pres_plasma_vol_avg=physics_variables.pres_plasma_vol_avg,
+ pres_plasma_on_axis=physics_variables.pres_plasma_thermal_on_axis,
+ pres_plasma_vol_avg=physics_variables.pres_plasma_thermal_vol_avg,
)
)
@@ -3794,6 +3928,12 @@ def calculate_plasma_current(
* b_plasma_toroidal_on_axis
)
# i_plasma_current == 2 case covered above
+ physics_variables.c_plasma_circular = (
+ (constants.TWOPI / constants.RMU0)
+ * rminor**2
+ / (rmajor * q95)
+ * b_plasma_toroidal_on_axis
+ )
# Calculate cyclindrical safety factor from IPDG89
qstar = (
@@ -3804,6 +3944,14 @@ def calculate_plasma_current(
* (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3))
)
+ physics_variables.q_95_circular = (
+ 2
+ * np.pi
+ * rminor**2
+ * b_plasma_toroidal_on_axis
+ / (constants.RMU0 * rmajor * plasma_current)
+ )
+
# Normalised beta from Troyon beta limit
physics_variables.beta_norm_total = (
1.0e8
@@ -3829,6 +3977,63 @@ def calculate_plasma_current(
return b_plasma_poloidal_average, qstar, plasma_current
+ def calculate_plasma_beta(self, pres_plasma: float, b_field: float) -> float:
+ """
+ Calculate the plasma beta for a given pressure and field.
+
+ :param pres_plasma: Plasma pressure (in Pascals).
+ :type pres_plasma: float
+ :param b_field: Magnetic field strength (in Tesla).
+ :type b_field: float
+ :returns: The plasma beta (dimensionless).
+ :rtype: float
+
+ Plasma beta is the ratio of plasma pressure to magnetic pressure.
+ """
+ if b_field == 0:
+ return 0.0
+ return 2 * constants.RMU0 * pres_plasma / (b_field**2)
+
+ def calculate_poloidal_field_profile(
+ self,
+ j_plasma_on_axis: float,
+ alphaj: float,
+ n_profile_elements: int,
+ rminor: float,
+ ) -> np.ndarray:
+ """
+ This method uses Ampère's law and the circular plasma approximation with a parabolic current profile
+ to compute the poloidal magnetic field profile.
+
+ :param float j_plasma_on_axis: On-axis plasma current density (A/m^2).
+ :param float alphaj: Current profile index (dimensionless).
+ :param int n_profile_elements: Number of radial profile elements (half the total number of points across the diameter).
+ :param float rminor: Plasma minor radius (m).
+
+ :returns: Poloidal magnetic field profile (T) across the full plasma diameter (length 2 * n_profile_elements).
+ :rtype: numpy.ndarray
+
+ """
+ n_profile_elements = n_profile_elements // 2
+ rho = np.linspace(
+ 0, rminor, n_profile_elements
+ ) # Normalized minor radius (0 to 1)
+ r = rho * rminor
+
+ j_plasma_profile = j_plasma_on_axis * (1 - (rho / rminor) ** 2) ** alphaj
+
+ integrand = 2 * np.pi * (r * j_plasma_profile) # dI/dr'
+ dr = r[1] - r[0]
+ c_enclosed = np.cumsum(integrand * dr) # cumulative current
+
+ # Poloidal field (avoid division by zero at axis)
+ b_poloidal = np.zeros_like(r)
+ mask = rho > 0
+ b_poloidal[mask] = constants.RMU0 / (2 * np.pi * r[mask]) * c_enclosed[mask]
+
+ # Mirror the b_poloidal profile and concatenate to match the doubled toroidal field profile
+ return np.concatenate([b_poloidal[::-1], b_poloidal])
+
def outtim(self):
po.oheadr(self.outfile, "Times")
@@ -4166,6 +4371,13 @@ def outplas(self):
physics_variables.plasma_current / 1.0e6,
"OP ",
)
+ po.ovarrf(
+ self.outfile,
+ "Plasma current for a circular plasma (A)",
+ "(c_plasma_circular)",
+ physics_variables.c_plasma_circular,
+ "OP ",
+ )
if physics_variables.i_alphaj == 1:
po.ovarrf(
@@ -4200,6 +4412,13 @@ def outplas(self):
physics_variables.j_plasma_on_axis,
"OP ",
)
+ po.ovarrf(
+ self.outfile,
+ "On-axis circular plasma current density (A/m2)",
+ "(j_plasma_circular_on_axis)",
+ physics_variables.j_plasma_circular_on_axis,
+ "OP ",
+ )
po.ovarrf(
self.outfile,
"Vertical field at plasma (T)",
@@ -4214,6 +4433,41 @@ def outplas(self):
"(b_plasma_toroidal_on_axis)",
physics_variables.b_plasma_toroidal_on_axis,
)
+ for i in range(len(physics_variables.b_plasma_toroidal_profile)):
+ po.ovarre(
+ self.mfile,
+ f"Toroidal field in plasma at point {i}",
+ f"b_plasma_toroidal_profile{i}",
+ physics_variables.b_plasma_toroidal_profile[i],
+ )
+ for i in range(len(physics_variables.b_plasma_poloidal_profile)):
+ po.ovarre(
+ self.mfile,
+ f"Poloidal field in plasma at point {i}",
+ f"b_plasma_poloidal_profile{i}",
+ physics_variables.b_plasma_poloidal_profile[i],
+ )
+ for i in range(len(physics_variables.b_plasma_circular_poloidal_profile)):
+ po.ovarre(
+ self.mfile,
+ f"Poloidal field in circular plasma at point {i}",
+ f"b_plasma_circular_poloidal_profile{i}",
+ physics_variables.b_plasma_circular_poloidal_profile[i],
+ )
+ for i in range(len(physics_variables.b_plasma_total_profile)):
+ po.ovarre(
+ self.mfile,
+ f"Total field in plasma at point {i}",
+ f"b_plasma_total_profile{i}",
+ physics_variables.b_plasma_total_profile[i],
+ )
+ for i in range(len(physics_variables.q_circular_profile)):
+ po.ovarre(
+ self.mfile,
+ f"Safety factor for circular plasma at point {i}",
+ f"q_circular_profile{i}",
+ physics_variables.q_circular_profile[i],
+ )
po.ovarrf(
self.outfile,
"Average poloidal field (T)",
@@ -4257,6 +4511,12 @@ def outplas(self):
physics_variables.qstar,
"OP ",
)
+ po.ovarrf(
+ self.outfile,
+ "Edge safety factor for a circular plasma (q_circular)",
+ "(q_95_circular)",
+ physics_variables.q_95_circular,
+ )
if physics_variables.i_plasma_geometry == 1:
po.ovarrf(
@@ -4372,6 +4632,27 @@ def outplas(self):
physics_variables.beta_toroidal,
"OP ",
)
+ for i in range(len(physics_variables.beta_toroidal_profile)):
+ po.ovarre(
+ self.mfile,
+ f"Beta toroidal profile at point {i}",
+ f"beta_toroidal_profile{i}",
+ physics_variables.beta_toroidal_profile[i],
+ )
+ for i in range(len(physics_variables.beta_poloidal_profile)):
+ po.ovarre(
+ self.mfile,
+ f"Beta poloidal profile at point {i}",
+ f"beta_poloidal_profile{i}",
+ physics_variables.beta_poloidal_profile[i],
+ )
+ for i in range(len(physics_variables.beta_total_profile)):
+ po.ovarre(
+ self.mfile,
+ f"Beta total profile at point {i}",
+ f"beta_total_profile{i}",
+ physics_variables.beta_total_profile[i],
+ )
po.ovarre(
self.outfile,
"Fast alpha beta",
@@ -4605,15 +4886,15 @@ def outplas(self):
po.ovarre(
self.outfile,
"Plasma pressure on axis (Pa)",
- "(pres_plasma_on_axis)",
- physics_variables.pres_plasma_on_axis,
+ "(pres_plasma_thermal_on_axis)",
+ physics_variables.pres_plasma_thermal_on_axis,
"OP ",
)
po.ovarre(
self.outfile,
"Volume averaged plasma pressure (Pa)",
- "(pres_plasma_vol_avg)",
- physics_variables.pres_plasma_vol_avg,
+ "(pres_plasma_thermal_vol_avg)",
+ physics_variables.pres_plasma_thermal_vol_avg,
"OP ",
)
# As array output is not currently supported, each element is output as a float instance
diff --git a/process/plasma_profiles.py b/process/plasma_profiles.py
index 2d154fa4eb..aded19f186 100644
--- a/process/plasma_profiles.py
+++ b/process/plasma_profiles.py
@@ -27,7 +27,7 @@ class PlasmaProfile:
parameterise_plasma(): Initializes the density and temperature profile averages and peak values.
parabolic_paramterisation(): Parameterizes plasma profiles in the case where i_plasma_pedestal=0.
pedestal_parameterisation(): Instance temperature and density profiles then integrate them, setting physics variables temp_plasma_electron_density_weighted_kev and temp_plasma_ion_density_weighted_kev.
- calculate_profile_factors(): Calculate and set the central pressure (pres_plasma_on_axis) using the ideal gas law and the pressure profile index (alphap).
+ calculate_profile_factors(): Calculate and set the central pressure (pres_plasma_thermal_on_axis) using the ideal gas law and the pressure profile index (alphap).
calculate_parabolic_profile_factors(): The gradient information for i_plasma_pedestal = 0.
"""
@@ -41,7 +41,7 @@ def __init__(self) -> None:
neprofile (NeProfile): An instance of the NeProfile class.
teprofile (TeProfile): An instance of the TeProfile class.
"""
- # Default profile_size = 501, but it's possible to experiment with this value.
+ # Default profile_size = 500, but it's possible to experiment with this value.
self.profile_size = 501
physics_variables.n_plasma_profile_elements = self.profile_size
self.outfile = constants.NOUT
@@ -255,10 +255,10 @@ def pedestal_parameterisation(self) -> None:
def calculate_profile_factors(self) -> None:
"""
- Calculate and set the central pressure (pres_plasma_on_axis) using the ideal gas law and the pressure profile index (alphap).
+ Calculate and set the central pressure (pres_plasma_thermal_on_axis) using the ideal gas law and the pressure profile index (alphap).
- This method calculates the central pressure (pres_plasma_on_axis) using the ideal gas law and the pressure profile index (alphap).
- It sets the value of the physics variable `pres_plasma_on_axis`.
+ This method calculates the central pressure (pres_plasma_thermal_on_axis) using the ideal gas law and the pressure profile index (alphap).
+ It sets the value of the physics variable `pres_plasma_thermal_on_axis`.
Args:
None
@@ -269,7 +269,7 @@ def calculate_profile_factors(self) -> None:
# Central pressure (Pa), from ideal gas law : p = nkT
- physics_variables.pres_plasma_on_axis = (
+ physics_variables.pres_plasma_thermal_on_axis = (
(
physics_variables.nd_plasma_electron_on_axis
* physics_variables.temp_plasma_electron_on_axis_kev
@@ -312,7 +312,7 @@ def calculate_profile_factors(self) -> None:
)
# Pressure profile index (N.B. no pedestal effects included here)
- # N.B. pres_plasma_on_axis is NOT equal to * (1 + alphap), but p(rho) = n(rho)*T(rho)
+ # N.B. pres_plasma_thermal_on_axis is NOT equal to
* (1 + alphap), but p(rho) = n(rho)*T(rho)
# and
= .T_n where <...> denotes volume-averages and T_n is the
# density-weighted temperature
@@ -320,8 +320,9 @@ def calculate_profile_factors(self) -> None:
# Shall assume that the pressure profile is parabolic. Can find volume average from
# profile index and core value the same as for density and temperature
- physics_variables.pres_plasma_vol_avg = (
- physics_variables.pres_plasma_on_axis / (physics_variables.alphap + 1)
+ physics_variables.pres_plasma_thermal_vol_avg = (
+ physics_variables.pres_plasma_thermal_on_axis
+ / (physics_variables.alphap + 1)
)
# Central plasma current density (A/m^2)
@@ -335,6 +336,17 @@ def calculate_profile_factors(self) -> None:
)
)
+ physics_variables.j_plasma_circular_on_axis = (
+ (physics_variables.c_plasma_circular)
+ * 2
+ / (
+ sp.special.beta(0.5, physics_variables.alphaj + 1)
+ * 2
+ * np.pi
+ * physics_variables.rminor
+ )
+ )
+
@staticmethod
def calculate_parabolic_profile_factors() -> None:
"""
diff --git a/tests/integration/data/large_tokamak_MFILE.DAT b/tests/integration/data/large_tokamak_MFILE.DAT
index 1cf005fc1a..259813c18a 100644
--- a/tests/integration/data/large_tokamak_MFILE.DAT
+++ b/tests/integration/data/large_tokamak_MFILE.DAT
@@ -466,7 +466,7 @@
Electron_number_density_on_axis_(/m3)____________________________________ (nd_plasma_electron_on_axis)__________________________ 1.05329687553494565e+20 OP
Line-averaged_electron_number_density_(/m3)______________________________ (nd_electron_line)_____________ 8.86821331867437138e+19 OP
Plasma_pressure_on_axis_(Pa)_____________________________________________ (pres_plasma_on_axis)___________________________ 8.42056568743616925e+05 OP
- Volume_averaged_plasma_pressure_(Pa)_____________________________________ (pres_plasma_vol_avg)_____________ 2.44074367751773010e+05 OP
+ Volume_averaged_plasma_pressure_(Pa)_____________________________________ (pres_plasma_thermal_vol_avg)_____________ 2.44074367751773010e+05 OP
Line-averaged_electron_density_/_Greenwald_density_______________________ (dnla_gw)______________________ 1.19999998892791648e+00 OP
Total_Ion_number_density_(/m3)___________________________________________ (nd_ions_total)________________ 7.10762537507155067e+19 OP
Fuel_ion_number_density_(/m3)____________________________________________ (nd_fuel_ions)_________________ 6.43726858733289964e+19 OP
diff --git a/tests/unit/test_plasma_profiles.py b/tests/unit/test_plasma_profiles.py
index 4c5f633116..37dbe8e27a 100644
--- a/tests/unit/test_plasma_profiles.py
+++ b/tests/unit/test_plasma_profiles.py
@@ -141,7 +141,7 @@ class PlasmaProfilesParam(NamedTuple):
temp_plasma_electron_on_axis_kev: float = 0.0
- pres_plasma_on_axis: float = 0.0
+ pres_plasma_thermal_on_axis: float = 0.0
nd_plasma_separatrix_electron: float = 0.0
@@ -225,7 +225,7 @@ class PlasmaProfilesParam(NamedTuple):
alphap=0.0,
tbeta=2,
temp_plasma_electron_on_axis_kev=0.0,
- pres_plasma_on_axis=0.0,
+ pres_plasma_thermal_on_axis=0.0,
nd_plasma_separatrix_electron=3.6421334486704804e19,
temp_plasma_separatrix_kev=0.10000000000000001,
pcoef=0.0,
@@ -270,7 +270,7 @@ class PlasmaProfilesParam(NamedTuple):
alphap=2.4500000000000002,
tbeta=2,
temp_plasma_electron_on_axis_kev=27.369013322953624,
- pres_plasma_on_axis=868071.46874220832,
+ pres_plasma_thermal_on_axis=868071.46874220832,
nd_plasma_separatrix_electron=3.6421334486704804e19,
temp_plasma_separatrix_kev=0.10000000000000001,
pcoef=1.1110842637642833,
@@ -355,8 +355,8 @@ def test_plasma_profiles(plasmaprofilesparam, monkeypatch):
monkeypatch.setattr(
physics_variables,
- "pres_plasma_on_axis",
- plasmaprofilesparam.pres_plasma_on_axis,
+ "pres_plasma_thermal_on_axis",
+ plasmaprofilesparam.pres_plasma_thermal_on_axis,
)
monkeypatch.setattr(
@@ -482,7 +482,7 @@ def test_plasma_profiles(plasmaprofilesparam, monkeypatch):
plasmaprofilesparam.expected_te0
)
- assert physics_variables.pres_plasma_on_axis == pytest.approx(
+ assert physics_variables.pres_plasma_thermal_on_axis == pytest.approx(
plasmaprofilesparam.expected_p0
)