diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 5c13687ee0..c30155c8a4 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -806,10 +806,23 @@
pres_plasma_on_axis: float = None
"""central total plasma pressure (Pa)"""
+pres_plasma_total_profile: list[float] = None
+"""Profile of total pressure in plasma (Pa)"""
+
+pres_plasma_electron_profile: list[float] = None
+"""Profile of electron pressure in plasma (Pa)"""
+
+pres_plasma_ion_total_profile: list[float] = None
+"""Profile of ion pressure in plasma (Pa)"""
+
+pres_plasma_fuel_profile: list[float] = None
+"""Profile of fuel pressure in plasma (Pa)"""
j_plasma_on_axis: float = None
"""Central plasma current density (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)"""
@@ -1476,7 +1489,12 @@ def init_physics_variables():
global nd_plasma_ions_on_axis
global m_s_limit
global pres_plasma_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 n_plasma_profile_elements
global f_dd_branching_trit
global pden_plasma_alpha_mw
global pden_alpha_total_mw
@@ -1723,7 +1741,12 @@ def init_physics_variables():
nd_plasma_ions_on_axis = 0.0
m_s_limit = 0.3
pres_plasma_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
+ n_plasma_profile_elements = 500
f_dd_branching_trit = 0.0
pden_plasma_alpha_mw = 0.0
pden_alpha_total_mw = 0.0
diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
index 58573b8347..abd7533324 100644
--- a/process/io/plot_proc.py
+++ b/process/io/plot_proc.py
@@ -27,6 +27,7 @@
import numpy as np
from matplotlib.patches import Circle, Rectangle
from matplotlib.path import Path
+from scipy.interpolate import interp1d
import process.confinement_time as confine
import process.constants as constants
@@ -3671,6 +3672,35 @@ def plot_n_profiles(prof, demo_ranges, mfile_data, scan):
verticalalignment="top",
bbox=props_density,
)
+
+ textstr_ions = "\n".join((
+ r"$\langle n_{\text{ions-total}} \rangle $: "
+ f"{mfile_data.data['nd_ions_total'].get_scan(scan):.3e} m$^{{-3}}$",
+ r"$\langle n_{\text{fuel}} \rangle $: "
+ f"{mfile_data.data['nd_fuel_ions'].get_scan(scan):.3e} m$^{{-3}}$",
+ r"$\langle n_{\text{alpha}} \rangle $: "
+ f"{mfile_data.data['nd_alphas'].get_scan(scan):.3e} m$^{{-3}}$",
+ r"$\langle n_{\text{impurities}} \rangle $: "
+ f"{mfile_data.data['nd_impurities'].get_scan(scan):.3e} m$^{{-3}}$",
+ r"$\langle n_{\text{protons}} \rangle $:"
+ f"{mfile_data.data['nd_protons'].get_scan(scan):.3e} m$^{{-3}}$",
+ ))
+
+ prof.text(
+ 1.2,
+ -0.725,
+ textstr_ions,
+ fontsize=9,
+ verticalalignment="bottom",
+ horizontalalignment="left",
+ transform=prof.transAxes,
+ bbox={
+ "boxstyle": "round",
+ "facecolor": "wheat",
+ "alpha": 0.5,
+ },
+ )
+
prof.grid(True, which="both", linestyle="--", linewidth=0.5, alpha=0.2)
# ---
@@ -9928,23 +9958,27 @@ def plot_fusion_rate_profiles(axis, fig, mfile_data, scan):
fusrat_plasma_dd_helion_profile = []
fusrat_plasma_dhe3_profile = []
+ n_plasma_profile_elements = int(
+ mfile_data.data["n_plasma_profile_elements"].get_scan(scan)
+ )
+
fusrat_plasma_dt_profile = [
mfile_data.data[f"fusrat_plasma_dt_profile{i}"].get_scan(scan)
- for i in range(500)
+ for i in range(n_plasma_profile_elements)
]
fusrat_plasma_dd_triton_profile = [
mfile_data.data[f"fusrat_plasma_dd_triton_profile{i}"].get_scan(scan)
- for i in range(500)
+ for i in range(n_plasma_profile_elements)
]
fusrat_plasma_dd_helion_profile = [
mfile_data.data[f"fusrat_plasma_dd_helion_profile{i}"].get_scan(scan)
- for i in range(500)
+ for i in range(n_plasma_profile_elements)
]
fusrat_plasma_dhe3_profile = [
mfile_data.data[f"fusrat_plasma_dhe3_profile{i}"].get_scan(scan)
- for i in range(500)
+ for i in range(n_plasma_profile_elements)
]
fusrat_plasma_total_profile = [
@@ -10409,6 +10443,254 @@ def plot_cover_page(axis, mfile_data, scan, fig, colour_scheme):
inset_ax.axis("off")
+def plot_plasma_pressure_profiles(axis, mfile_data, scan):
+ # Plot the plasma pressure profiles on the given axis
+ n_plasma_profile_elements = int(
+ mfile_data.data["n_plasma_profile_elements"].get_scan(scan)
+ )
+
+ pres_plasma_profile = [
+ mfile_data.data[f"pres_plasma_electron_profile{i}"].get_scan(scan)
+ for i in range(n_plasma_profile_elements)
+ ]
+ pres_plasma_profile_ion = [
+ mfile_data.data[f"pres_plasma_ion_total_profile{i}"].get_scan(scan)
+ for i in range(n_plasma_profile_elements)
+ ]
+ pres_plasma_total_profile = [
+ mfile_data.data[f"pres_plasma_total_profile{i}"].get_scan(scan)
+ for i in range(n_plasma_profile_elements)
+ ]
+ pres_plasma_profile_fuel = [
+ mfile_data.data[f"pres_plasma_fuel_profile{i}"].get_scan(scan)
+ for i in range(n_plasma_profile_elements)
+ ]
+ pres_plasma_profile_kpa = [p / 1000.0 for p in pres_plasma_profile]
+ pres_plasma_profile_ion_kpa = [p / 1000.0 for p in pres_plasma_profile_ion]
+ pres_plasma_profile_fuel_kpa = [p / 1000.0 for p in pres_plasma_profile_fuel]
+ pres_plasma_profile_total_kpa = [p / 1000.0 for p in pres_plasma_total_profile]
+
+ axis.plot(
+ np.linspace(0, 1, len(pres_plasma_profile_kpa)),
+ pres_plasma_profile_kpa,
+ color="blue",
+ label="Electron",
+ )
+ axis.plot(
+ np.linspace(0, 1, len(pres_plasma_profile_ion_kpa)),
+ pres_plasma_profile_ion_kpa,
+ color="Red",
+ label="Ion-total",
+ )
+ axis.plot(
+ np.linspace(0, 1, len(pres_plasma_profile_fuel_kpa)),
+ pres_plasma_profile_fuel_kpa,
+ color="orange",
+ label="Fuel",
+ )
+ axis.plot(
+ np.linspace(0, 1, len(pres_plasma_profile_total_kpa)),
+ pres_plasma_profile_total_kpa,
+ color="green",
+ label="Total",
+ )
+ axis.set_xlabel("$\\rho$ [r/a]")
+ axis.set_ylabel("Pressure [kPa]")
+ axis.minorticks_on()
+ axis.grid(which="minor", linestyle=":", linewidth=0.5, alpha=0.5)
+ axis.set_title("Plasma Pressure Profiles")
+ axis.grid(True, linestyle="--", alpha=0.5)
+ axis.set_xlim([0, 1.025])
+ axis.set_ylim(bottom=0)
+ axis.legend()
+
+
+def plot_plasma_pressure_gradient_profiles(axis, mfile_data, scan):
+ # Get the plasma pressure profiles
+ n_plasma_profile_elements = int(
+ mfile_data.data["n_plasma_profile_elements"].get_scan(scan)
+ )
+
+ pres_plasma_profile = [
+ mfile_data.data[f"pres_plasma_electron_profile{i}"].get_scan(scan)
+ for i in range(n_plasma_profile_elements)
+ ]
+ pres_plasma_profile_ion = [
+ mfile_data.data[f"pres_plasma_ion_total_profile{i}"].get_scan(scan)
+ for i in range(n_plasma_profile_elements)
+ ]
+ pres_plasma_profile_total = [
+ mfile_data.data[f"pres_plasma_total_profile{i}"].get_scan(scan)
+ for i in range(n_plasma_profile_elements)
+ ]
+ pres_plasma_profile_fuel = [
+ mfile_data.data[f"pres_plasma_fuel_profile{i}"].get_scan(scan)
+ for i in range(n_plasma_profile_elements)
+ ]
+ pres_plasma_profile_kpa = np.array(pres_plasma_profile) / 1000.0
+ pres_plasma_profile_ion_kpa = np.array(pres_plasma_profile_ion) / 1000.0
+ pres_plasma_profile_fuel_kpa = np.array(pres_plasma_profile_fuel) / 1000.0
+ pres_plasma_profile_total_kpa = np.array(pres_plasma_profile_total) / 1000.0
+
+ # Calculate the normalized radius
+ rho = np.linspace(0, 1, len(pres_plasma_profile_kpa))
+
+ # Compute gradients using numpy.gradient
+ grad_electron = np.gradient(pres_plasma_profile_kpa, rho)
+ grad_ion = np.gradient(pres_plasma_profile_ion_kpa, rho)
+ grad_total = np.gradient(pres_plasma_profile_total_kpa, rho)
+ grad_fuel = np.gradient(pres_plasma_profile_fuel_kpa, rho)
+
+ axis.plot(rho, grad_electron, color="blue", label="Electron")
+ axis.plot(rho, grad_ion, color="red", label="Ion")
+ axis.plot(rho, grad_total, color="green", label="Total")
+ axis.plot(rho, grad_fuel, color="orange", label="Fuel")
+ axis.set_xlabel("$\\rho$ [r/a]")
+ axis.set_ylabel("$dP/dr$ [kPa / m]")
+ axis.minorticks_on()
+ axis.grid(which="minor", linestyle=":", linewidth=0.5, alpha=0.5)
+ axis.set_title("Plasma Pressure Gradient Profiles")
+ axis.grid(True, linestyle="--", alpha=0.5)
+ axis.set_xlim([0, 1.025])
+ axis.legend()
+
+
+def plot_plasma_poloidal_pressure_contours(
+ axis: plt.Axes, mfile_data: mf.MFile, scan: int
+) -> None:
+ """
+ Plot plasma poloidal pressure contours inside the plasma boundary.
+
+ :param axis: Matplotlib axis object to plot on.
+ :type axis: matplotlib.axes.Axes
+ :param mfile_data: MFILE data object containing plasma and geometry data.
+ :type mfile_data: mf.MFile
+ :param scan: Scan number to use for extracting data.
+ :type scan: int
+
+ This function visualizes the poloidal pressure distribution inside the plasma boundary
+ by interpolating the pressure profile onto a grid defined by the plasma geometry.
+ The pressure is shown as filled contours, with the plasma boundary overlaid.
+ """
+
+ n_plasma_profile_elements = int(
+ mfile_data.data["n_plasma_profile_elements"].get_scan(scan)
+ )
+
+ # Get pressure profile (function of normalized radius rho, 0..1)
+ pres_plasma_electron_profile = [
+ mfile_data.data[f"pres_plasma_electron_profile{i}"].get_scan(scan)
+ for i in range(n_plasma_profile_elements)
+ ]
+ pres_plasma_profile_ion = [
+ mfile_data.data[f"pres_plasma_ion_total_profile{i}"].get_scan(scan)
+ for i in range(n_plasma_profile_elements)
+ ]
+
+ # Convert pressure to kPa
+ pres_plasma_electron_profile_kpa = [
+ p / 1000.0 for p in pres_plasma_electron_profile
+ ]
+ pres_plasma_profile_ion_kpa = [p / 1000.0 for p in pres_plasma_profile_ion]
+ pres_plasma_profile = [
+ e + i
+ for e, i in zip(
+ pres_plasma_electron_profile_kpa, pres_plasma_profile_ion_kpa, strict=False
+ )
+ ]
+
+ # Get plasma geometry and boundary
+ pg = plasma_geometry(
+ rmajor=mfile_data.data["rmajor"].get_scan(scan),
+ rminor=mfile_data.data["rminor"].get_scan(scan),
+ triang=mfile_data.data["triang"].get_scan(scan),
+ kappa=mfile_data.data["kappa"].get_scan(scan),
+ i_single_null=mfile_data.data["i_single_null"].get_scan(scan),
+ i_plasma_shape=mfile_data.data["i_plasma_shape"].get_scan(scan),
+ square=mfile_data.data["plasma_square"].get_scan(scan),
+ )
+ rminor = mfile_data.data["rminor"].get_scan(scan)
+ rmajor = mfile_data.data["rmajor"].get_scan(scan)
+ # Create a grid of (R, Z) points inside the plasma boundary
+ n_rho = 500
+ n_theta = 720
+ rho = np.linspace(0, 1, n_rho)
+ theta = np.linspace(0, 2 * np.pi, n_theta)
+ rho_grid, theta_grid = np.meshgrid(rho, theta)
+
+ # Map (rho, theta) to (R, Z) using plasma boundary shape
+ # For each theta, get boundary (R, Z), then scale by rho
+ boundary_r = pg.rs
+ boundary_z = pg.zs
+ # Interpolate boundary for all theta
+ boundary_theta = np.arctan2(boundary_z - pg.zs.mean(), boundary_r - pg.rs.mean())
+ # Ensure boundary_theta is monotonic and covers [0, 2pi]
+ boundary_theta = np.unwrap(boundary_theta)
+ # Sort boundary_theta and corresponding r/z for monotonic interpolation
+ sort_idx = np.argsort(boundary_theta)
+ boundary_theta = boundary_theta[sort_idx]
+ boundary_r = boundary_r[sort_idx]
+ boundary_z = boundary_z[sort_idx]
+ # Extend boundary to cover full [0, 2pi] if needed
+ if boundary_theta[0] > 0 or boundary_theta[-1] < 2 * np.pi:
+ boundary_theta = np.concatenate(([0], boundary_theta, [2 * np.pi]))
+ boundary_r = np.concatenate(([boundary_r[0]], boundary_r, [boundary_r[-1]]))
+ boundary_z = np.concatenate(([boundary_z[0]], boundary_z, [boundary_z[-1]]))
+ # Map theta to boundary r/z
+ f_r = interp1d(
+ boundary_theta,
+ boundary_r,
+ kind="linear",
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )
+ # Map theta to boundary z
+ f_z = interp1d(
+ boundary_theta,
+ boundary_z,
+ kind="linear",
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )
+ # For each (theta, rho), get boundary (R, Z), then scale by rho
+ # Use the boundary center for scaling, not mean, to avoid vertical offset
+ r_center = rmajor
+ z_center = pg.zs.mean()
+ r_grid = r_center + (f_r(theta_grid) - r_center) * rho_grid
+ z_grid = z_center + (f_z(theta_grid) - z_center) * rho_grid
+
+ # Interpolate pressure profile for each rho
+ pressure_grid = np.interp(
+ rho_grid, np.linspace(0, 1, len(pres_plasma_profile)), pres_plasma_profile
+ )
+
+ # Mask points outside the plasma boundary (optional, but grid is inside by construction)
+ # Plot filled contour
+ c = axis.contourf(r_grid, z_grid, pressure_grid, levels=50, cmap="plasma")
+ c = axis.contourf(r_grid, -z_grid, pressure_grid, levels=50, cmap="plasma")
+
+ # Overlay plasma boundary
+ axis.plot(pg.rs, pg.zs, color="black", linewidth=2, label="Plasma Boundary")
+
+ # Add colorbar for pressure (now in kPa)
+ # You can control the location using the 'location' argument ('left', 'right', 'top', 'bottom')
+ # For more control, use 'ax' or 'fraction', 'pad', etc.
+ # Example: location="right", pad=0.05, fraction=0.05
+ axis.figure.colorbar(
+ c, ax=axis, label="Pressure [kPa]", location="left", anchor=(-0.25, 0.5)
+ )
+
+ axis.set_aspect("equal")
+ axis.set_xlabel("R [m]")
+ axis.set_xlim(rmajor - 1.2 * rminor, rmajor + 1.2 * rminor)
+ axis.set_ylim(
+ -1.2 * rminor * mfile_data.data["kappa"].get_scan(scan),
+ 1.2 * mfile_data.data["kappa"].get_scan(scan) * rminor,
+ )
+ axis.set_ylabel("Z [m]")
+ axis.set_title("Plasma Poloidal Pressure Contours")
+
+
def main_plot(
fig0,
fig1,
@@ -10429,6 +10711,7 @@ def main_plot(
fig16,
fig17,
fig18,
+ fig19,
m_file_data,
scan,
imp="../data/lz_non_corona_14_elements/",
@@ -10468,154 +10751,151 @@ def main_plot(
plot_cover_page(plot_0, m_file_data, scan, fig0, colour_scheme)
# Plot header info
- plot_1 = fig1.add_subplot(231)
- plot_header(plot_1, m_file_data, scan)
+ plot_header(fig1.add_subplot(231), m_file_data, scan)
# Geometry
- plot_2 = fig1.add_subplot(232)
- plot_geometry_info(plot_2, m_file_data, scan)
+ plot_geometry_info(fig1.add_subplot(232), m_file_data, scan)
# Physics
- plot_3 = fig1.add_subplot(233)
- plot_physics_info(plot_3, m_file_data, scan)
+ plot_physics_info(fig1.add_subplot(233), m_file_data, scan)
# Magnetics
- plot_4 = fig1.add_subplot(234)
- plot_magnetics_info(plot_4, m_file_data, scan)
+ plot_magnetics_info(fig1.add_subplot(234), m_file_data, scan)
# power/flow economics
- plot_5 = fig1.add_subplot(235)
- plot_power_info(plot_5, m_file_data, scan)
+ plot_power_info(fig1.add_subplot(235), m_file_data, scan)
# Current drive
- plot_6 = fig1.add_subplot(236)
- plot_current_drive_info(plot_6, m_file_data, scan)
+ plot_current_drive_info(fig1.add_subplot(236), m_file_data, scan)
fig1.subplots_adjust(wspace=0.25, hspace=0.25)
- plot_7 = fig2.add_subplot(121)
- plot_7.set_position([0.175, 0.1, 0.35, 0.8]) # Move plot slightly to the right
- plot_iteration_variables(plot_7, m_file_data, scan)
+ ax7 = fig2.add_subplot(121)
+ ax7.set_position([0.175, 0.1, 0.35, 0.8]) # Move plot slightly to the right
+ plot_iteration_variables(ax7, m_file_data, scan)
# Plot main plasma information
- plot_8 = fig3.add_subplot(111, aspect="equal")
- plot_main_plasma_information(plot_8, m_file_data, scan, colour_scheme, fig3)
+ plot_main_plasma_information(
+ fig3.add_subplot(111, aspect="equal"), m_file_data, scan, colour_scheme, fig3
+ )
# Plot density profiles
- plot_9 = fig4.add_subplot(231) # , aspect= 0.05)
- plot_9.set_position([0.075, 0.55, 0.25, 0.4])
- plot_n_profiles(plot_9, demo_ranges, m_file_data, scan)
+ ax9 = fig4.add_subplot(231)
+ ax9.set_position([0.075, 0.55, 0.25, 0.4])
+ plot_n_profiles(ax9, demo_ranges, m_file_data, scan)
# Plot temperature profiles
- plot_10 = fig4.add_subplot(232)
- plot_10.set_position([0.375, 0.55, 0.25, 0.4])
- plot_t_profiles(plot_10, demo_ranges, m_file_data, scan)
+ ax10 = fig4.add_subplot(232)
+ ax10.set_position([0.375, 0.55, 0.25, 0.4])
+ plot_t_profiles(ax10, demo_ranges, m_file_data, scan)
# Plot impurity profiles
- plot_11 = fig4.add_subplot(233)
- plot_11.set_position([0.7, 0.45, 0.25, 0.5])
- plot_radprofile(plot_11, m_file_data, scan, imp, demo_ranges)
+ ax11 = fig4.add_subplot(233)
+ ax11.set_position([0.7, 0.45, 0.25, 0.5])
+ plot_radprofile(ax11, m_file_data, scan, imp, demo_ranges)
# Plot current density profile
- plot_12 = fig4.add_subplot(4, 3, 10)
- plot_12.set_position([0.075, 0.125, 0.25, 0.15])
- plot_jprofile(plot_12)
+ ax12 = fig4.add_subplot(4, 3, 10)
+ ax12.set_position([0.075, 0.125, 0.25, 0.15])
+ plot_jprofile(ax12)
# Plot q profile
- plot_13 = fig4.add_subplot(4, 3, 12)
- plot_13.set_position([0.7, 0.125, 0.25, 0.15])
- plot_qprofile(plot_13, demo_ranges, m_file_data, scan)
+ ax13 = fig4.add_subplot(4, 3, 12)
+ ax13.set_position([0.7, 0.125, 0.25, 0.15])
+ plot_qprofile(ax13, demo_ranges, m_file_data, scan)
- plot_14 = fig5.add_subplot(122)
- plot_fusion_rate_profiles(plot_14, fig5, m_file_data, scan)
+ plot_fusion_rate_profiles(fig5.add_subplot(122), fig5, m_file_data, scan)
+
+ plot_plasma_pressure_profiles(fig6.add_subplot(222), m_file_data, scan)
+ plot_plasma_pressure_gradient_profiles(fig6.add_subplot(224), m_file_data, scan)
+ # Currently only works with Sauter geometry as plasma has a closed surface
+ if m_file_data.data["i_plasma_shape"].get_scan(scan) == 1:
+ plot_plasma_poloidal_pressure_contours(
+ fig6.add_subplot(121, aspect="equal"), m_file_data, scan
+ )
# Plot poloidal cross-section
- plot_15 = fig6.add_subplot(121, aspect="equal")
- poloidal_cross_section(plot_15, m_file_data, scan, demo_ranges, colour_scheme)
+ poloidal_cross_section(
+ fig7.add_subplot(121, aspect="equal"),
+ m_file_data,
+ scan,
+ demo_ranges,
+ colour_scheme,
+ )
# Plot toroidal cross-section
- plot_16 = fig6.add_subplot(122, aspect="equal")
- toroidal_cross_section(plot_16, m_file_data, scan, demo_ranges, colour_scheme)
- # fig4.subplots_adjust(bottom=-0.2, top = 0.9, left = 0.1, right = 0.9)
+ toroidal_cross_section(
+ fig7.add_subplot(122, aspect="equal"),
+ m_file_data,
+ scan,
+ demo_ranges,
+ colour_scheme,
+ )
# Plot color key
- plot_17 = fig6.add_subplot(222)
- plot_17.set_position([0.5, 0.5, 0.5, 0.5])
- color_key(plot_17, m_file_data, scan, colour_scheme)
+ ax17 = fig7.add_subplot(222)
+ ax17.set_position([0.5, 0.5, 0.5, 0.5])
+ color_key(ax17, m_file_data, scan, colour_scheme)
- plot_18 = fig7.add_subplot(211)
- plot_18.set_position([0.1, 0.33, 0.8, 0.6]) # x0, y0, width, height (2/3 vertical)
- plot_radial_build(plot_18, m_file_data, colour_scheme)
+ ax18 = fig8.add_subplot(211)
+ ax18.set_position([0.1, 0.33, 0.8, 0.6])
+ plot_radial_build(ax18, m_file_data, colour_scheme)
# Make each axes smaller vertically to leave room for the legend
- plot_185 = fig8.add_subplot(211)
- plot_185.set_position([0.1, 0.61, 0.8, 0.32]) # x0, y0, width, height
+ ax185 = fig9.add_subplot(211)
+ ax185.set_position([0.1, 0.61, 0.8, 0.32])
- plot_18 = fig8.add_subplot(212)
- plot_18.set_position([0.1, 0.13, 0.8, 0.32]) # x0, y0, width, height
- plot_upper_vertical_build(plot_185, m_file_data, colour_scheme)
- plot_lower_vertical_build(plot_18, m_file_data, colour_scheme)
+ ax18b = fig9.add_subplot(212)
+ ax18b.set_position([0.1, 0.13, 0.8, 0.32])
+ plot_upper_vertical_build(ax185, m_file_data, colour_scheme)
+ plot_lower_vertical_build(ax18b, m_file_data, colour_scheme)
# Can only plot WP and turn structure if superconducting coil at the moment
if m_file_data.data["i_tf_sup"].get_scan(scan) == 1:
# TF coil with WP
- plot_19 = fig9.add_subplot(231, aspect="equal")
- plot_19.set_position([
- 0.025,
- 0.5,
- 0.45,
- 0.45,
- ]) # Half height, a bit wider, top left
- plot_superconducting_tf_wp(plot_19, m_file_data, scan, fig9)
+ ax19 = fig10.add_subplot(231, aspect="equal")
+ ax19.set_position([0.025, 0.5, 0.45, 0.45])
+ plot_superconducting_tf_wp(ax19, m_file_data, scan, fig10)
# TF coil turn structure
- plot_20 = fig9.add_subplot(325, aspect="equal")
- plot_20.set_position([0.025, 0.1, 0.3, 0.3])
- plot_tf_turn(plot_20, fig9, m_file_data, scan)
+ ax20 = fig10.add_subplot(325, aspect="equal")
+ ax20.set_position([0.025, 0.1, 0.3, 0.3])
+ plot_tf_turn(ax20, fig10, m_file_data, scan)
else:
- plot_19 = fig9.add_subplot(211, aspect="equal")
- plot_19.set_position([0.06, 0.55, 0.675, 0.4])
- plot_resistive_tf_wp(plot_19, m_file_data, scan, fig8)
+ ax19 = fig10.add_subplot(211, aspect="equal")
+ ax19.set_position([0.06, 0.55, 0.675, 0.4])
+ plot_resistive_tf_wp(ax19, m_file_data, scan, fig9)
- plot_21 = fig10.add_subplot(111, aspect="equal")
- plot_tf_coil_structure(plot_21, m_file_data, scan, colour_scheme)
+ plot_tf_coil_structure(
+ fig11.add_subplot(111, aspect="equal"), m_file_data, scan, colour_scheme
+ )
- axes = fig11.subplots(nrows=3, ncols=1, sharex=True).flatten()
+ axes = fig12.subplots(nrows=3, ncols=1, sharex=True).flatten()
plot_tf_stress(axes)
- plot_23 = fig12.add_subplot(221)
- plot_bootstrap_comparison(plot_23, m_file_data, scan)
-
- plot_24 = fig12.add_subplot(224)
- plot_h_threshold_comparison(plot_24, m_file_data, scan)
-
- plot_25 = fig13.add_subplot(221)
- plot_density_limit_comparison(plot_25, m_file_data, scan)
-
- plot_26 = fig13.add_subplot(224)
- plot_confinement_time_comparison(plot_26, m_file_data, scan)
-
- plot_27 = fig14.add_subplot(111)
- plot_current_profiles_over_time(plot_27, m_file_data, scan)
-
- plot_28 = fig15.add_subplot(121, aspect="equal")
- plot_cs_coil_structure(plot_28, fig15, m_file_data, scan)
+ plot_bootstrap_comparison(fig13.add_subplot(221), m_file_data, scan)
+ plot_h_threshold_comparison(fig13.add_subplot(224), m_file_data, scan)
+ plot_density_limit_comparison(fig14.add_subplot(221), m_file_data, scan)
+ plot_confinement_time_comparison(fig14.add_subplot(224), m_file_data, scan)
+ plot_current_profiles_over_time(fig15.add_subplot(111), m_file_data, scan)
- plot_29 = fig15.add_subplot(224, aspect="equal")
- plot_cs_turn_structure(plot_29, fig15, m_file_data, scan)
-
- plot_30 = fig16.add_subplot(221, aspect="equal")
- plot_first_wall_top_down_cross_section(plot_30, m_file_data, scan)
-
- plot_31 = fig16.add_subplot(122)
- plot_first_wall_poloidal_cross_section(plot_31, m_file_data, scan)
+ plot_cs_coil_structure(
+ fig16.add_subplot(121, aspect="equal"), fig16, m_file_data, scan
+ )
+ plot_cs_turn_structure(
+ fig16.add_subplot(224, aspect="equal"), fig16, m_file_data, scan
+ )
- plot_32 = fig16.add_subplot(337)
- plot_fw_90_deg_pipe_bend(plot_32, m_file_data, scan)
+ plot_first_wall_top_down_cross_section(
+ fig17.add_subplot(221, aspect="equal"), m_file_data, scan
+ )
+ plot_first_wall_poloidal_cross_section(fig17.add_subplot(122), m_file_data, scan)
+ plot_fw_90_deg_pipe_bend(fig17.add_subplot(337), m_file_data, scan)
- plot_blkt_pipe_bends(fig17, m_file_data, scan)
+ plot_blkt_pipe_bends(fig18, m_file_data, scan)
- plot_33 = fig18.add_subplot(111, aspect="equal")
- plot_main_power_flow(plot_33, m_file_data, scan, fig18)
+ plot_main_power_flow(
+ fig19.add_subplot(111, aspect="equal"), m_file_data, scan, fig19
+ )
def main(args=None):
@@ -10925,6 +11205,7 @@ def main(args=None):
page16 = plt.figure(figsize=(12, 9), dpi=80)
page17 = plt.figure(figsize=(12, 9), dpi=80)
page18 = plt.figure(figsize=(12, 9), dpi=80)
+ page19 = plt.figure(figsize=(12, 9), dpi=80)
# run main_plot
main_plot(
@@ -10947,6 +11228,7 @@ def main(args=None):
page16,
page17,
page18,
+ page19,
m_file,
scan=scan,
demo_ranges=demo_ranges,
@@ -10974,6 +11256,7 @@ def main(args=None):
pdf.savefig(page16)
pdf.savefig(page17)
pdf.savefig(page18)
+ pdf.savefig(page19)
# show fig if option used
if args.show:
@@ -10998,6 +11281,7 @@ def main(args=None):
plt.close(page16)
plt.close(page17)
plt.close(page18)
+ plt.close(page19)
if __name__ == "__main__":
diff --git a/process/physics.py b/process/physics.py
index 5f32eb4509..ba0df63b78 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -4536,6 +4536,12 @@ def outplas(self):
po.oblnkl(self.outfile)
po.osubhd(self.outfile, "Temperature and Density (volume averaged) :")
+ po.ovarrf(
+ self.outfile,
+ "Number of radial points in plasma profiles",
+ "(n_plasma_profile_elements)",
+ physics_variables.n_plasma_profile_elements,
+ )
po.ovarrf(
self.outfile,
"Volume averaged electron temperature (keV)",
@@ -4610,6 +4616,24 @@ def outplas(self):
physics_variables.pres_plasma_vol_avg,
"OP ",
)
+ # As array output is not currently supported, each element is output as a float instance
+ # Output plasma pressure profiles to mfile
+ for name, arr in [
+ (
+ "Plasma electron pressure",
+ physics_variables.pres_plasma_electron_profile,
+ ),
+ ("Plasma ion pressure", physics_variables.pres_plasma_ion_total_profile),
+ ("Plasma total pressure", physics_variables.pres_plasma_total_profile),
+ ("Fuel pressure", physics_variables.pres_plasma_fuel_profile),
+ ]:
+ for i, val in enumerate(arr):
+ po.ovarre(
+ self.mfile,
+ f"{name} at point {i}",
+ f"{name.lower().replace(' ', '_')}_profile{i}",
+ val,
+ )
if stellarator_variables.istell == 0:
po.ovarre(
diff --git a/process/plasma_profiles.py b/process/plasma_profiles.py
index 596273a701..2d154fa4eb 100644
--- a/process/plasma_profiles.py
+++ b/process/plasma_profiles.py
@@ -43,6 +43,7 @@ def __init__(self) -> None:
"""
# Default profile_size = 501, 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
self.neprofile = profiles.NeProfile(self.profile_size)
self.teprofile = profiles.TeProfile(self.profile_size)
@@ -279,6 +280,37 @@ def calculate_profile_factors(self) -> None:
* constants.ELECTRON_CHARGE
)
+ # Electron pressure profile (Pa)
+ physics_variables.pres_plasma_electron_profile = self.neprofile.profile_y * (
+ self.teprofile.profile_y * constants.KILOELECTRON_VOLT
+ )
+
+ # Total ion pressure profile (Pa)
+ physics_variables.pres_plasma_ion_total_profile = (
+ physics_variables.nd_ions_total
+ * (self.neprofile.profile_y / physics_variables.nd_plasma_electrons_vol_avg)
+ ) * (
+ self.teprofile.profile_y
+ * constants.KILOELECTRON_VOLT
+ * physics_variables.f_temp_plasma_ion_electron
+ )
+
+ # Total pressure profile (Pa)
+ physics_variables.pres_plasma_total_profile = (
+ physics_variables.pres_plasma_electron_profile
+ + physics_variables.pres_plasma_ion_total_profile
+ )
+
+ # Fuel ion pressure profile (Pa)
+ physics_variables.pres_plasma_fuel_profile = (
+ physics_variables.nd_fuel_ions
+ * (self.neprofile.profile_y / physics_variables.nd_plasma_electrons_vol_avg)
+ ) * (
+ self.teprofile.profile_y
+ * constants.KILOELECTRON_VOLT
+ * physics_variables.f_temp_plasma_ion_electron
+ )
+
# 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)
# and
= .T_n where <...> denotes volume-averages and T_n is the
diff --git a/tests/integration/data/large_tokamak_MFILE.DAT b/tests/integration/data/large_tokamak_MFILE.DAT
index 21406bab9c..57106a9d99 100644
--- a/tests/integration/data/large_tokamak_MFILE.DAT
+++ b/tests/integration/data/large_tokamak_MFILE.DAT
@@ -455,6 +455,7 @@
R._Stambaugh_normalised_beta_upper_limit_________________________________ (beta_norm_max_stambaugh)______ 2.06087953258754553e+00 OP
Plasma_thermal_energy_derived_from_thermal_beta_(J)______________________ (e_plasma_beta_thermal)________ 9.66355464245535374e+08 OP
Plasma_thermal_energy_derived_from_the_total_beta_(J)____________________ (e_plasma_beta)________________ 1.10779607068200541e+09 OP
+ Number_of_radial_points_in_plasma_profiles_______________________________ (n_plasma_profile_elements)____ 500
Volume_averaged_electron_temperature_(kev)_______________________________ (temp_plasma_electron_vol_avg_kev)___________________________ 1.27326078294372866e+01
Ratio_of_ion_to_electron_volume-averaged_temperature_____________________ (f_temp_plasma_ion_electron)_______________________ 1.00000000000000000e+00 IP
Electron_temperature_on_axis_(keV)_______________________________________ (temp_plasma_electron_on_axis_kev)__________________________ 2.64346007040757769e+01 OP