diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index 6471685831..23a62f6c19 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -13,7 +13,8 @@ import matplotlib.pyplot as plt import numpy as np from matplotlib import patches -from matplotlib.patches import Circle, Rectangle +from matplotlib.axes import Axes +from matplotlib.patches import Circle, FancyBboxPatch, Rectangle from matplotlib.path import Path as mplPath from scipy.interpolate import interp1d @@ -185,6 +186,34 @@ class RadialBuild: rtangle2 = 2 * rtangle +def _box_style(colour: str): + return {"boxstyle": "round", "facecolor": colour, "alpha": 1.0, "linewidth": 2} + + +def _text_layout(fig): + return { + "fontsize": 9, + "verticalalignment": "bottom", + "horizontalalignment": "left", + "transform": fig.transFigure, + } + + +def _setup_axis(axis, xmin, xmax, ymin, ymax): + axis.set_ylim(ymin, ymax) + axis.set_xlim(xmin, xmax) + axis.set_axis_off() + axis.set_autoscaley_on(False) + axis.set_autoscalex_on(False) + + +def _colourbar(contour_fill, axis, colourbar_axis): + # Use a dedicated colorbar axes when provided so the main axes width is unchanged. + if colourbar_axis is None: + return axis.figure.colorbar(contour_fill, ax=axis, pad=0.02) + return axis.figure.colorbar(contour_fill, cax=colourbar_axis) + + def plot_plasma( axis: plt.Axes, mfile: MFile, @@ -419,17 +448,12 @@ def poloidal_cross_section( plot_tf_coils(axis, mfile, scan, colour_scheme) plot_pf_coils(axis, mfile, scan, colour_scheme) - # Ranges - # --- - # DEMO : Fixed ranges for comparison if demo_ranges: - axis.set_ylim([-15, 15]) - axis.set_xlim([0, 20]) + axis.set_ylim(-15, 15) + axis.set_xlim(0, 20) - # Adaptive ranges else: - axis.set_xlim([0, axis.get_xlim()[1]]) - # --- + axis.set_xlim(0, axis.get_xlim()[1]) def plot_full_machine_poloidal_cross_section( @@ -522,7 +546,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur # Display the plasma image over the figure, not the axes new_ax = axis.inset_axes( - [-0.15, 0.6, 0.45, 0.45], transform=axis.transAxes, zorder=1 + (-0.15, 0.6, 0.45, 0.45), transform=axis.transAxes, zorder=1 ) new_ax.imshow(plasma) new_ax.axis("off") @@ -543,7 +567,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur neutron = mpimg.imread(img_path.open("rb")) new_ax = axis.inset_axes( - [0.2, 0.85, 0.03, 0.03], transform=axis.transAxes, zorder=10 + (0.2, 0.85, 0.03, 0.03), transform=axis.transAxes, zorder=10 ) new_ax.imshow(neutron) new_ax.axis("off") @@ -690,7 +714,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur # Display the alpha particle image over the figure, not the axes new_ax = axis.inset_axes( - [0.16, 0.95, 0.025, 0.025], transform=axis.transAxes, zorder=10 + (0.16, 0.95, 0.025, 0.025), transform=axis.transAxes, zorder=10 ) new_ax.imshow(alpha) new_ax.axis("off") @@ -718,10 +742,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.37, 0.775, f"$P_{{\\text{{neutron}}}}$:\n{mfile.get('p_neutron_total_mw', scan=scan):,.2f} MW", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "grey", @@ -741,16 +762,8 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.0725, 0.83, f"$P_{{\\text{{HCD,primary}}}}$: {mfile.get('p_hcd_primary_injected_mw', scan=scan) + mfile.get('p_hcd_primary_extra_heat_mw', scan=scan):.2f} MW", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightyellow", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("lightyellow"), ) # Add HCD secondary injected power @@ -758,16 +771,8 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.0725, 0.725, f"$P_{{\\text{{HCD,secondary}}}}$: {mfile.get('p_hcd_secondary_injected_mw', scan=scan) + mfile.get('p_hcd_secondary_extra_heat_mw', scan=scan):.2f} MW", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightyellow", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("lightyellow"), ) # Load the HCD injector image @@ -776,11 +781,11 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur # Display the injector image over the figure, not the axes new_ax = axis.inset_axes( - [-0.2, 0.8, 0.15, 0.15], transform=axis.transAxes, zorder=10 + (-0.2, 0.8, 0.15, 0.15), transform=axis.transAxes, zorder=10 ) new_ax.imshow(hcd_injector_1) new_ax.axis("off") - new_ax = axis.inset_axes([-0.2, 0.5, 0.15, 0.5], transform=axis.transAxes, zorder=10) + new_ax = axis.inset_axes((-0.2, 0.5, 0.15, 0.5), transform=axis.transAxes, zorder=10) new_ax.imshow(hcd_injector_2) new_ax.axis("off") @@ -816,16 +821,8 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.04, 0.45, "\n\nH&CD Power Supply\n\n", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightyellow", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("lightyellow"), zorder=4, ) @@ -850,17 +847,8 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.2, 0.435, f"$P_{{\\text{{secondary,loss}}}}$:\n{mfile.get('p_hcd_secondary_electric_mw', scan=scan) * (1.0 - mfile.get('eta_hcd_secondary_injector_wall_plug', scan=scan)):.2f} MWe", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightblue", - "alpha": 1.0, - "linewidth": 2, - "linestyle": "dashed", - }, + **_text_layout(fig), + bbox=_box_style("lightblue") | {"linestyle": "dashed"}, ) # Draw an arrow from HCD secondary losses to the total secondary heat power @@ -932,17 +920,8 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.2, 0.485, f"$P_{{\\text{{primary,loss}}}}$:\n{mfile.get('p_hcd_primary_electric_mw', scan=scan) * (1.0 - mfile.get('eta_hcd_primary_injector_wall_plug', scan=scan)):.2f} MWe", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightblue", - "alpha": 1.0, - "linewidth": 2, - "linestyle": "dashed", - }, + **_text_layout(fig), + bbox=_box_style("lightblue") | {"linestyle": "dashed"}, ) # Draw arrow from HCD primary electric box to HCD power supply box @@ -977,17 +956,9 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur axis.text( 0.12, 0.35, - f"$P_{{\\text{{secondary}}}}$:\n{mfile.get('p_hcd_secondary_electric_mw', scan=scan):.2f} MWe \n$\\eta$: {mfile.get('eta_hcd_secondary_injector_wall_plug', scan=scan):.2f}", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightyellow", - "alpha": 1.0, - "linewidth": 2, - }, + f"$P_{{\\text{{secondary}}}}$:\n{mfile.get('p_hcd_secondary_electric_mw', scan=scan):.2f} MWe\n$\\eta$: {mfile.get('eta_hcd_secondary_injector_wall_plug', scan=scan):.2f}", + **_text_layout(fig), + bbox=_box_style("lightyellow"), ) # Plot HCD primary electric box @@ -995,16 +966,8 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.025, 0.35, f"$P_{{\\text{{primary}}}}$:\n{mfile.get('p_hcd_primary_electric_mw', scan=scan):.2f} MWe\n$\\eta$: {mfile.get('eta_hcd_primary_injector_wall_plug', scan=scan):.2f}", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightyellow", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("lightyellow"), ) # ============================================= @@ -1022,13 +985,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur verticalalignment="bottom", horizontalalignment="center", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightblue", - "alpha": 1.0, - "linewidth": 2, - "linestyle": "dashed", - }, + bbox=_box_style("lightblue") | {"linestyle": "dashed"}, zorder=4, ) @@ -1043,7 +1000,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur turbine = mpimg.imread(img_path.open("rb")) # Display the turbine image over the figure, not the axes - new_ax = axis.inset_axes([1.1, 0.0, 0.15, 0.15], transform=axis.transAxes, zorder=10) + new_ax = axis.inset_axes((1.1, 0.0, 0.15, 0.15), transform=axis.transAxes, zorder=10) new_ax.imshow(turbine) new_ax.axis("off") @@ -1051,17 +1008,9 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur axis.text( 0.9, 0.25, - f"$P_{{\\text{{primary,thermal}}}}$:\n{mfile.get('p_plant_primary_heat_mw', scan=scan):,.2f} MW \n$\\eta_{{\\text{{turbine}}}}$: {mfile.get('eta_turbine', scan=scan):.3f}", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "orange", - "alpha": 1.0, - "linewidth": 2, - }, + f"$P_{{\\text{{primary,thermal}}}}$:\n{mfile.get('p_plant_primary_heat_mw', scan=scan):,.2f} MW\n$\\eta_{{\\text{{turbine}}}}$: {mfile.get('eta_turbine', scan=scan):.3f}", + **_text_layout(fig), + bbox=_box_style("orange"), ) # Draw arrow from bend to turbine inlet @@ -1100,7 +1049,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur # Display the generator image over the figure, not the axes new_ax = axis.inset_axes( - [0.96, 0.0, 0.15, 0.15], transform=axis.transAxes, zorder=10 + (0.96, 0.0, 0.15, 0.15), transform=axis.transAxes, zorder=10 ) new_ax.imshow(generator) new_ax.axis("off") @@ -1125,10 +1074,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.79, 0.16, "Generator", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), zorder=20, ) @@ -1169,7 +1115,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur # Display the pylon image over the figure, not the axes new_ax = axis.inset_axes( - [0.925, -0.1, 0.1, 0.1], transform=axis.transAxes, zorder=10 + (0.925, -0.1, 0.1, 0.1), transform=axis.transAxes, zorder=10 ) new_ax.imshow(pylon) new_ax.axis("off") @@ -1179,16 +1125,8 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.68, 0.15, f"$P_{{\\text{{gross}}}}$:\n{mfile.get('p_plant_electric_gross_mw', scan=scan):,.2f} MWe", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lime", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("lime"), ) # Gross to net electric power @@ -1211,17 +1149,8 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.875, 0.05, f"$P_{{\\text{{loss}}}}$:\n{mfile.get('p_turbine_loss_mw', scan=scan):,.2f} MWth", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "orange", - "alpha": 1.0, - "linewidth": 2, - "linestyle": "dashed", - }, + **_text_layout(fig), + bbox=_box_style("orange") | {"linestyle": "dashed"}, ) # Shield primary thermal to plant total primary thermal arrow @@ -1244,16 +1173,8 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.68, 0.05, f"$P_{{\\text{{net,electric}}}}$:\n{mfile.get('p_plant_electric_net_mw', scan=scan):,.2f} MWe", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lime", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("lime"), ) # Plot the recirculated electric power box @@ -1264,16 +1185,8 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur f"$P_{{\\text{{recirc,electric}}}}$:\n{mfile.get('p_plant_electric_recirc_mw', scan=scan):,.2f} MWe\n" f"$f_{{\\text{{recirc}}}}$:\n{mfile.get('f_p_plant_electric_recirc', scan=scan):,.2f}" ), - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lime", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("lime"), ) # Gross to recirculated power arrow @@ -1377,7 +1290,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur fw = mpimg.imread(img_path.open("rb")) # Display the first wall image over the figure, not the axes - new_ax = axis.inset_axes([0.4, 0.625, 0.4, 0.4], transform=axis.transAxes, zorder=10) + new_ax = axis.inset_axes((0.4, 0.625, 0.4, 0.4), transform=axis.transAxes, zorder=10) new_ax.imshow(fw) new_ax.axis("off") @@ -1397,16 +1310,8 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.46, 0.85, f"$P_{{\\text{{FW, }}\\alpha}}$:\n{mfile.get('p_fw_alpha_mw', scan=scan):.2f} MW", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "red", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("red"), ) # Neutron power incident on first wall box @@ -1414,10 +1319,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.46, 0.775, f"$P_{{\\text{{FW,nuclear}}}}$:\n{mfile.get('p_fw_nuclear_heat_total_mw', scan=scan):,.2f} MW", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "grey", @@ -1431,10 +1333,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.46, 0.71, f"$P_{{\\text{{FW,rad}}}}$:\n{mfile.get('p_fw_rad_total_mw', scan=scan):,.2f} MW", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "dodgerblue", @@ -1493,10 +1392,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.5, 0.555, f"Primary thermal\n(inc pump): {mfile.get('p_fw_heat_deposited_mw', scan=scan):,.2f} MWth", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "orange", @@ -1509,10 +1405,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.7, 0.555, f"Primary thermal\n(inc pump): {mfile.get('p_blkt_heat_deposited_mw', scan=scan):,.2f} MWth", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "orange", @@ -1525,10 +1418,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.875, 0.555, f"Primary thermal:\n{mfile.get('p_shld_heat_deposited_mw', scan=scan):.2f} MWth", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "orange", @@ -1616,16 +1506,8 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.6, 0.49, f"Primary thermal (inc pump): {mfile.get('p_fw_blkt_heat_deposited_mw', scan=scan):,.2f} MWth\n", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "orange", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("orange"), ) # Load the blanket image @@ -1636,7 +1518,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur # Display the blanket image over the figure, not the axes new_ax = axis.inset_axes( - [0.75, 0.625, 0.4, 0.4], transform=axis.transAxes, zorder=10 + (0.75, 0.625, 0.4, 0.4), transform=axis.transAxes, zorder=10 ) new_ax.imshow(blanket) new_ax.axis("off") @@ -1657,14 +1539,11 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.625, 0.775, ( - f"$P_{{\\text{{Blkt,nuclear}}}}$:\n{mfile.get('p_blkt_nuclear_heat_total_mw', scan=scan):,.2f} MW \n" + f"$P_{{\\text{{Blkt,nuclear}}}}$:\n{mfile.get('p_blkt_nuclear_heat_total_mw', scan=scan):,.2f} MW\n" f"$P_{{\\text{{Blkt,multiplication}}}}$:\n{mfile.get('p_blkt_multiplication_mw', scan=scan):,.2f} MW\n" f"$f_{{\\text{{multiplication}}}}$:\n{mfile.get('f_p_blkt_multiplication', scan=scan):,.2f}" ), - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "grey", @@ -1679,7 +1558,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur # Display the vacuum vessel image over the figure, not the axes new_ax = axis.inset_axes( - [0.975, 0.625, 0.4, 0.4], transform=axis.transAxes, zorder=10 + (0.975, 0.625, 0.4, 0.4), transform=axis.transAxes, zorder=10 ) new_ax.imshow(vv) new_ax.axis("off") @@ -1700,17 +1579,8 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.38, 0.375, f"$P_{{\\text{{shld,secondary}}}}$:\n{mfile.get('p_shld_secondary_heat_mw', scan=scan):,.2f} MWth", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightblue", - "alpha": 1.0, - "linewidth": 2, - "linestyle": "dashed", - }, + **_text_layout(fig), + bbox=_box_style("lightblue") | {"linestyle": "dashed"}, ) # Shield secondary power box to secondary heat total @@ -1782,7 +1652,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur divertor = mpimg.imread(img_path.open("rb")) # Display the divertor image over the figure, not the axes - new_ax = axis.inset_axes([0.1, 0.4, 0.3, 0.25], transform=axis.transAxes, zorder=10) + new_ax = axis.inset_axes((0.1, 0.4, 0.3, 0.25), transform=axis.transAxes, zorder=10) new_ax.imshow(divertor) new_ax.axis("off") @@ -1791,10 +1661,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.29, 0.57, f"$P_{{\\text{{div,rad}}}}$:\n{mfile.get('p_div_rad_total_mw', scan=scan):,.2f} MW", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "dodgerblue", @@ -1808,10 +1675,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.4, 0.58, f"$P_{{\\text{{div,nuclear}}}}$:\n{mfile.get('p_div_nuclear_heat_total_mw', scan=scan):,.2f} MW", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "grey", @@ -1829,10 +1693,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur f"Solid angle fraction: {mfile.get('f_ster_div_single', scan=scan):.3f}\n" f"Primary heat fraction: {mfile.get('f_p_div_primary_heat', scan=scan):.3f}" ), - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "orange", @@ -1846,17 +1707,8 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.3, 0.375, f"$P_{{\\text{{div,secondary}}}}$:\n{mfile.get('p_div_secondary_heat_mw', scan=scan):.2f} MWth", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightblue", - "alpha": 1.0, - "linewidth": 2, - "linestyle": "dashed", - }, + **_text_layout(fig), + bbox=_box_style("lightblue") | {"linestyle": "dashed"}, ) # Divertor to divertor secondary heat arrow @@ -1917,10 +1769,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.55, 0.33, f"$P_{{\\text{{div,pump}}}}$: {mfile.get('p_div_coolant_pump_mw', scan=scan):.2f} MW", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "wheat", @@ -1979,10 +1828,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.875, 0.325, f"$P_{{\\text{{shld,pump}}}}$:\n{mfile.get('p_shld_coolant_pump_mw', scan=scan):.2f} MW", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "wheat", @@ -1996,10 +1842,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.725, 0.4, f"$P_{{\\text{{FW + Blkt}}}}$:\n{mfile.get('p_fw_blkt_coolant_pump_mw', scan=scan):.2f} MW", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "wheat", @@ -2061,10 +1904,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur f"Coolant pumps electric:\n{mfile.get('p_coolant_pump_elec_total_mw', scan=scan):.3f} MWe\n" f"$\\eta$: {mfile.get('eta_coolant_pump_electric', scan=scan):.3f}" ), - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "lime", @@ -2078,10 +1918,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.7, 0.325, f"Coolant pumps total:\n{mfile.get('p_coolant_pump_total_mw', scan=scan):.3f} MW", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "wheat", @@ -2110,10 +1947,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.5, 0.235, f"Coolant pumps losses total:\n{mfile.get('p_coolant_pump_loss_total_mw', scan=scan):.3f} MWth", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "lightblue", @@ -2166,10 +2000,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.49, 0.05, f"Cryo Plant:\n{mfile.get('p_cryo_plant_electric_mw', scan=scan):.3f} MWe", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "burlywood", @@ -2198,10 +2029,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.4, 0.05, f"Tritium Plant:\n{mfile.get('p_tritium_plant_electric_mw', scan=scan):.3f} MWe", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "burlywood", @@ -2230,10 +2058,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.575, 0.05, f"Vacuum pumps:\n{mfile.get('vachtmw', scan=scan):.3f} MWe", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "burlywood", @@ -2266,10 +2091,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur f"Minimum base load:\n{mfile.get('p_plant_electric_base', scan=scan) * 1.0e-6:.3f} MWe\n" f"Plant floor power density:\n{mfile.get('pflux_plant_floor_electric', scan=scan) * 1.0e-3:.3f} kW$\\text{{m}}^{{-2}}$" ), - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "burlywood", @@ -2283,10 +2105,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.325, 0.075, f"TF coils:\n{mfile.get('p_tf_electric_supplies_mw', scan=scan):.3f} MWe", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "burlywood", @@ -2300,10 +2119,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.25, 0.05, f"PF coils:\n{mfile.get('p_pf_electric_supplies_mw', scan=scan):.3f} MWe", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "burlywood", @@ -2362,10 +2178,7 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.46, 0.285, f"$P_{{\\text{{HCD,loss}}}}$:\n{mfile.get('p_hcd_secondary_heat_mw', scan=scan):.2f} MWth", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, + **_text_layout(fig), bbox={ "boxstyle": "round", "facecolor": "lightblue", @@ -2412,17 +2225,8 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur 0.155, 0.25, f"$P_{{\\text{{TF,nuclear}}}}$:\n{mfile.get('p_tf_nuclear_heat_mw', scan=scan):.2f} MWth", - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightblue", - "alpha": 1.0, - "linewidth": 2, - "linestyle": "dashed", - }, + **_text_layout(fig), + bbox=_box_style("lightblue") | {"linestyle": "dashed"}, ) # TF nuclear heat to secondary heat total box arrow @@ -2444,67 +2248,45 @@ def plot_main_power_flow(axis: plt.Axes, mfile: MFile, scan: int, fig: plt.Figur # =========================================== -def plot_main_plasma_information( - axis: plt.Axes, - mfile: MFile, - scan: int, - colour_scheme: Literal[1, 2], - fig: plt.Figure, +def plasma_main_page( + axis, + fig, + rmajor, + rminor, + triang, + radius_plasma_core_norm, + kappa, + i_single_null, + plasma_square, + big_q_plasma, + p_fw_alpha_mw, + f_p_alpha_plasma_deposited, + p_neutron_total_mw, + pflux_plasma_surface_neutron_avg_mw, ): - """Plots the main plasma information including plasma shape, geometry, currents, heating, - confinement, and other relevant plasma parameters. - - Parameters - ---------- - axis : plt.Axes - The matplotlib axis object to plot on. - mfile : MFile - The MFILE data object containing plasma parameters. - scan : int - The scan number to use for extracting data. - colour_scheme : int - The colour scheme to use for plots. - fig : plt.Figure - The matplotlib figure object for additional annotations. - """ - # Import key variables - triang = mfile.get("triang", scan=scan) - kappa = mfile.get("kappa", scan=scan) - - # Remove the axes - axis.axis("off") + """Plot plasma from main page""" + white_box = {"boxstyle": "round", "facecolor": "white", "alpha": 1.0} - # Plot the main plasma shape - plot_plasma(axis, mfile, scan, colour_scheme) - - rmajor = mfile.get("rmajor", scan=scan) - rminor = mfile.get("rminor", scan=scan) - # Get the plasma permieter points for the core plasma region pg = plasma_geometry( rmajor=rmajor, - rminor=mfile.get("rminor", scan=scan) - * mfile.get("radius_plasma_core_norm", scan=scan), - triang=mfile.get("triang", scan=scan), - kappa=mfile.get("kappa", scan=scan), - i_single_null=mfile.get("i_single_null", scan=scan), + rminor=rminor * radius_plasma_core_norm, + triang=triang, + kappa=kappa, + i_single_null=i_single_null, i_plasma_shape=1, - square=mfile.get("plasma_square", scan=scan), + square=plasma_square, ) - # Plot the core plasma boundary line axis.plot(pg.rs, pg.zs, color="black", linestyle="--") - - # Plot the centre of the plasma axis.plot(rmajor, 0, "r+", markersize=20, markeredgewidth=2) - # Add Q plasma information box axis.text( 0.725, 0.175, - f"$Q_{{\\text{{plasma}}}}$: {mfile.get('big_q_plasma', scan=scan):.2f}", + f"$Q_{{\\text{{plasma}}}}$: {big_q_plasma:.2f}", fontsize=15, verticalalignment="center", horizontalalignment="center", - bbox={"boxstyle": "round", "facecolor": "white", "alpha": 1.0}, + bbox=white_box, transform=fig.transFigure, ) @@ -2526,7 +2308,7 @@ def plot_main_plasma_information( fontsize=9, color="black", ha="center", - bbox={"boxstyle": "round", "facecolor": "white", "alpha": 1.0}, + bbox=white_box, ) # ============================================ @@ -2547,7 +2329,7 @@ def plot_main_plasma_information( fontsize=9, color="black", ha="center", - bbox={"boxstyle": "round", "facecolor": "white", "alpha": 1.0}, + bbox=white_box, ) # ============================================ @@ -2564,13 +2346,13 @@ def plot_main_plasma_information( axis.text( 0.3, 0.75, - f"$\\kappa$: {mfile.get('kappa', scan=scan):.2f}", + f"$\\kappa$: {kappa:.2f}", fontsize=9, color="black", rotation=270, verticalalignment="center", transform=axis.transAxes, - bbox={"boxstyle": "round", "facecolor": "white", "alpha": 1.0}, + bbox=white_box, ) # ============================================= @@ -2587,23 +2369,18 @@ def plot_main_plasma_information( axis.text( rmajor - (rminor * triang * 0.75), kappa * rminor * 0.3, - f"$\\delta$: {mfile.get('triang', scan=scan):.2f}", + f"$\\delta$: {triang:.2f}", fontsize=9, color="black", - rotation=0, verticalalignment="center", - bbox={"boxstyle": "round", "facecolor": "white", "alpha": 1.0}, + bbox=white_box, ) - # ============================================= - - radius_plasma_core_norm = mfile.get("radius_plasma_core_norm", scan=scan) - # Draw a double-ended arrow for the plasma core region axis.annotate( "", - xy=(rmajor, -rminor * 0.1 * kappa), # Inner plasma edge - xytext=(rmajor + (rminor * radius_plasma_core_norm), -rminor * 0.1 * kappa), + xy=(rmajor, -rminor * 0.1 * kappa), + xytext=(rmajor + rminor * radius_plasma_core_norm, -rminor * 0.1 * kappa), arrowprops={"arrowstyle": "<->", "color": "black"}, ) # Add a label for core region @@ -2613,606 +2390,464 @@ def plot_main_plasma_information( f"$\\rho_{{\\text{{core}}}}$: {radius_plasma_core_norm:.2f}", fontsize=9, color="black", - rotation=0, verticalalignment="center", - bbox={"boxstyle": "round", "facecolor": "white", "alpha": 1.0}, + bbox=white_box, ) - # ================================================ - - # Add plasma volume, areas and shaping information - - geom_type = PlasmaGeometryModelType(mfile.get("i_plasma_geometry", scan=scan)) + with ( + resources.path( + "process.core.io.plot.images", "alpha_particle.png" + ) as image_path, + image_path.open("rb") as image_file, + ): + alpha_particle = mpimg.imread(image_file) - textstr_plasma = ( - f"$\\mathbf{{Shaping:}}$\n \n" - f"$\\kappa_{{95}}$: {mfile.get('kappa95', scan=scan):.2f} ({geom_type.kappa95_model.description}) | $\\delta_{{95}}$: {mfile.get('triang95', scan=scan):.2f} ({geom_type.triang95_model.description}) | $\\zeta$: {mfile.get('plasma_square', scan=scan):.2f}\n" - f"$\\kappa$: {mfile.get('kappa', scan=scan):.2f} ({geom_type.kappa_model.description}) | $\\delta$: {mfile.get('triang', scan=scan):.2f} ({geom_type.triang_model.description}) | A: {mfile.get('aspect', scan=scan):.2f}\n" - f"$ V_{{\\text{{p}}}}:$ {mfile.get('vol_plasma', scan=scan):,.2f}$ \\ \\text{{m}}^3$ | $ A_{{\\text{{p,surface}}}}:$ {mfile.get('a_plasma_surface', scan=scan):,.2f}$ \\ \\text{{m}}^2$ | $ A_{{\\text{{p,poloidal}}}}:$ {mfile.get('a_plasma_poloidal', scan=scan):,.3f}$ \\ \\text{{m}}^2$\n" - f"$ L_{{\\text{{p,poloidal}}}}:$ {mfile.get('len_plasma_poloidal', scan=scan):,.3f}$ \\ \\text{{m}}$" + image_axis = axis.inset_axes( + (0.975, 0.275, 0.075, 0.075), + transform=axis.transAxes, + zorder=10, ) + image_axis.imshow(alpha_particle) + image_axis.axis("off") + axis.annotate( + "", + xy=(rmajor + rminor, -rminor * kappa * 0.55), + xytext=(rmajor + 0.2 * rminor, -rminor * kappa * 0.25), + arrowprops={"facecolor": "red", "edgecolor": "grey", "lw": 1}, + ) axis.text( - 0.365, - 0.975, - textstr_plasma, + 1.0, + 0.275, + f"$P_{{\\alpha,\\text{{loss}}}}$ {p_fw_alpha_mw:.2f} MW\n" + f"$f_{{\\alpha,\\text{{coupled}}}}$ " + f"{f_p_alpha_plasma_deposited:.2f}", fontsize=9, + color="black", verticalalignment="top", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightyellow", - "alpha": 1.0, - "linewidth": 2, - }, + transform=axis.transAxes, + bbox={"boxstyle": "round", "facecolor": "red", "alpha": 1.0, "linewidth": 2}, ) - # ============================================ + with ( + resources.path("process.core.io.plot.images", "neutron.png") as image_path, + image_path.open("rb") as image_file, + ): + neutron = mpimg.imread(image_file) - # Draw a red arrow coming from the right and pointing at the plasma - axis.annotate( - "", - xy=(rmajor + (rminor * 0.8), kappa * rminor * 0.2), # Pointing at the plasma - xytext=( - rmajor + (rminor * 1.4), - kappa * rminor * 0.2, - ), # Starting point of the arrow - arrowprops={"facecolor": "red", "edgecolor": "red", "lw": 2}, + image_axis = axis.inset_axes( + (0.975, 0.75, 0.075, 0.075), + transform=axis.transAxes, + zorder=10, ) + image_axis.imshow(neutron) + image_axis.axis("off") - # Draw a red arrow coming from the right and pointing at the plasma axis.annotate( "", - xy=(rmajor + (rminor * 0.8), -kappa * rminor * 0.2), # Pointing at the plasma - xytext=( - rmajor + (rminor * 1.4), - -kappa * rminor * 0.2, - ), # Starting point of the arrow - arrowprops={"facecolor": "red", "edgecolor": "red", "lw": 2}, - ) - - i_hcd_primary = mfile.get("i_hcd_primary", scan=scan) - i_hcd_secondary = mfile.get("i_hcd_secondary", scan=scan) - - # Add heating and current drive information - textstr_hcd = ( - f"$\\mathbf{{Heating \\ & \\ current \\ drive:}}$\n \n" - f"Total injected heat: {mfile.get('p_hcd_injected_total_mw', scan=scan):.3f} MW \n" - f"Ohmic heating power: {mfile.get('p_plasma_ohmic_mw', scan=scan):.3f} MW \n\n" - f"$\\mathbf{{Primary \\ system: {CurrentDriveModel(i_hcd_primary).abbreviation}}}$ \n" - f"Current driving power {mfile.get('p_hcd_primary_injected_mw', scan=scan):.4f} MW\n" - f"Extra heat power: {mfile.get('p_hcd_primary_extra_heat_mw', scan=scan):.4f} MW\n" - f"$\\eta_{{\\text{{CD,prim}}}}$: {mfile.get('eta_cd_hcd_primary', scan=scan):.4f} A/W | $\\langle\\zeta_{{\\text{{CD,prim}}}}\\rangle$: {mfile.get('eta_cd_dimensionless_hcd_primary', scan=scan):.4f} \n" - f"$\\gamma_{{\\text{{CD,prim}}}}$: {mfile.get('eta_cd_norm_hcd_primary', scan=scan):.4f} $\\times 10^{{20}} \\mathrm{{A}} / \\mathrm{{Wm}}^2$\n" - f"Current driven by primary: {mfile.get('c_hcd_primary_driven', scan=scan) / 1e6:.3f} MA\n\n" - f"$\\mathbf{{Secondary \\ system: {CurrentDriveModel(i_hcd_secondary).abbreviation}}}$ \n" - f"Current driving power {mfile.get('p_hcd_secondary_injected_mw', scan=scan):.4f} MW\n" - f"Extra heat power: {mfile.get('p_hcd_secondary_extra_heat_mw', scan=scan):.4f} MW\n" - f"$\\eta_{{\\text{{CD,sec}}}}$: {mfile.get('eta_cd_hcd_secondary', scan=scan):.4f} A/W | $\\langle\\zeta_{{\\text{{CD,sec}}}}\\rangle$: {mfile.get('eta_cd_dimensionless_hcd_secondary', scan=scan):.4f} \n" - f"$\\gamma_{{\\text{{CD,sec}}}}$: {mfile.get('eta_cd_norm_hcd_secondary', scan=scan):.4f} $\\times 10^{{20}} \\mathrm{{A}} / \\mathrm{{Wm}}^2$\n" - f"Current driven by secondary: {mfile.get('c_hcd_secondary_driven', scan=scan) / 1e6:.3f} MA" + xy=(rmajor + rminor, rminor * kappa * 0.65), + xytext=(rmajor, rminor * kappa * 0.5), + arrowprops={"facecolor": "grey", "edgecolor": "grey", "lw": 1}, ) - axis.text( - 0.73, - 0.675, - textstr_hcd, + 0.775, + 0.875, + f"$P_{{\\text{{n,total}}}}$ {p_neutron_total_mw:.2f} MW\n" + f"$\\phi_{{\\text{{n,avg}}}}$ " + f"{pflux_plasma_surface_neutron_avg_mw:.3f} MW/m²", fontsize=9, + color="black", verticalalignment="top", - transform=plt.gcf().transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "paleturquoise", - "alpha": 1.0, - "linewidth": 2, - "edgecolor": "black", - }, - ) - - # Add injected power label - axis.text( - 0.92, - 0.625, - "$P_{\\text{inj}}$", - fontsize=23, - verticalalignment="top", - transform=fig.transFigure, + transform=axis.transAxes, + bbox={"boxstyle": "round", "facecolor": "grey", "alpha": 0.8, "linewidth": 2}, ) + for kap in (-kappa, kappa): + axis.annotate( + "", + xy=(rmajor + rminor * 0.8, kap * rminor * 0.2), + xytext=(rmajor + rminor * 1.4, kap * rminor * 0.2), + arrowprops={"facecolor": "red", "edgecolor": "red", "lw": 2}, + ) - # ================================================ - # Add beta information - textstr_beta = ( - f"$\\mathbf{{Beta \\ Information:}}$\n \n" - f"Total beta,$ \\ \\langle \\beta \\rangle$: {mfile.get('beta_total_vol_avg', scan=scan):.4f}\n" - f"Thermal beta,$ \\ \\langle \\beta_{{\\text{{thermal}}}} \\rangle$: {mfile.get('beta_thermal_vol_avg', scan=scan):.4f}\n" - f"Toroidal beta,$ \\ \\langle \\beta_{{\\text{{t}}}} \\rangle$: {mfile.get('beta_toroidal_vol_avg', scan=scan):.4f}\n" - f"Poloidal beta,$ \\ \\langle \\beta_{{\\text{{p}}}} \\rangle$: {mfile.get('beta_poloidal_vol_avg', scan=scan):.4f}\n" - f"Fast-alpha beta,$ \\ \\langle \\beta_{{\\alpha}} \\rangle$: {mfile.get('beta_fast_alpha', scan=scan):.4f}\n" - f"Upper limit on {BetaComponentLimits(int(mfile.get('i_beta_component', scan=scan))).full_name}: $ \\langle \\beta \\rangle$: {mfile.get('beta_vol_avg_max', scan=scan):.4f}\n" - f"Normalised total beta,$ \\ \\beta_{{\\text{{N}}}}$: {mfile.get('beta_norm_total', scan=scan):.4f}\n" - f"Normalised thermal beta,$ \\ \\beta_{{\\text{{N,thermal}}}}$: {mfile.get('beta_norm_thermal', scan=scan):.4f}\n" - f"Maximum normalised beta ({BetaNormMaxModel(int(mfile.get('i_beta_norm_max', scan=scan))).full_name}),$ \\ \\beta_{{\\text{{N,max}}}}$: {mfile.get('beta_norm_max', scan=scan):.4f}" - ) +def plot_main_plasma_information( + axis: plt.Axes, + mfile: MFile, + scan: int, + colour_scheme: Literal[1, 2], + fig: plt.Figure, +): + """Plots the main plasma information including plasma shape, geometry, currents, heating, + confinement, and other relevant plasma parameters. - axis.text( - 0.025, - 0.975, - textstr_beta, - fontsize=9, - verticalalignment="top", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightblue", - "alpha": 1.0, - "linewidth": 2, - }, - ) + Parameters + ---------- + axis : plt.Axes + The matplotlib axis object to plot on. + mfile : MFile + The MFILE data object containing plasma parameters. + scan : int + The scan number to use for extracting data. + colour_scheme : int + The colour scheme to use for plots. + fig : plt.Figure + The matplotlib figure object for additional annotations. + """ - # Add beta label - axis.text( - 0.27, - 0.94, - "$\\beta$", - fontsize=23, - verticalalignment="top", - transform=fig.transFigure, - ) + def value(name: str): + return mfile.get(name, scan=scan) + + def add_panel( + x: float, + top: float, + width: float, + height: float, + text: str, + style: dict, + text_x_offset: float = 0.004, + ): + """Add one fixed-size box and its unchanged original text.""" + box_expansion = 0.002 + left = x - box_expansion + bottom = top - height - box_expansion + expanded_width = width + 2 * box_expansion + expanded_height = height + 2 * box_expansion + + patch = FancyBboxPatch( + (left, bottom), + expanded_width, + expanded_height, + boxstyle="round,pad=0,rounding_size=0.003", + transform=fig.transFigure, + facecolor=style.get("facecolor", "white"), + edgecolor=style.get("edgecolor", "black"), + linewidth=style.get("linewidth", 1.0), + alpha=style.get("alpha", 1.0), + clip_on=False, + zorder=20, + ) + axis.add_patch(patch) - # ================================================ + return axis.text( + x + text_x_offset, + top - 0.004, + text, + fontsize=9, + color="black", + verticalalignment="top", + horizontalalignment="left", + transform=fig.transFigure, + clip_on=False, + zorder=21, + ) - # Add volt-second information - textstr_volt_second = ( - f"$\\mathbf{{Volt-second \\ requirements:}}$\n \n" - f"Total volt-second consumption: {mfile.get('vs_plasma_total_required', scan=scan):.4f} Vs \n" - f" - Internal volt-seconds: {mfile.get('vs_plasma_internal', scan=scan):.4f} Vs\n" - f" - Volt-seconds needed for burn: {mfile.get('vs_plasma_burn_required', scan=scan):.4f} Vs\n" - f" - Volt-seconds needed for ramp: {mfile.get('vs_plasma_ramp_required', scan=scan):.4f} Vs | $C_{{\\text{{ejima}}}}$: {mfile.get('ejima_coeff', scan=scan):.4f}\n" - f"$V_{{\\text{{loop}}}}$: {mfile.get('v_plasma_loop_burn', scan=scan):.4f} V\n" - f"$\\Omega_{{\\text{{p}}}}$: {mfile.get('res_plasma', scan=scan):.4e} $\\Omega$\n" - f"Plasma resistive diffusion time: {mfile.get('t_plasma_res_diffusion', scan=scan):,.4f} s\n" - f"Plasma inductance: {mfile.get('ind_plasma', scan=scan):.4e} H | ITER $l_i(3)$: {mfile.get('ind_plasma_internal_norm_iter_3', scan=scan):.4f}\n" - f"Plasma stored magnetic energy: {mfile.get('e_plasma_magnetic_stored', scan=scan) / 1e9:.4f} GJ\n" - f"Plasma normalised internal inductance, $l_i$ ({IndInternalNormModel(int(mfile.get('i_ind_plasma_internal_norm', scan=scan))).full_name}) :{mfile.get('ind_plasma_internal_norm', scan=scan):.3f}" - ) + def add_symbol(panel_name: str, x: float, top: float): + """Add a large symbol at its original figure-relative position.""" + return axis.text( + x, + top, + panel_name, + fontsize=23, + color="black", + verticalalignment="top", + transform=fig.transFigure, + clip_on=False, + zorder=22, + ) - axis.text( - 0.025, - 0.78, - textstr_volt_second, - fontsize=9, - verticalalignment="top", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightgreen", - "alpha": 1.0, - "linewidth": 2, - }, - ) + triang = value("triang") + kappa = value("kappa") + rmajor = value("rmajor") + rminor = value("rminor") + radius_plasma_core_norm = value("radius_plasma_core_norm") - # Add volt second label - axis.text( - 0.30, - 0.77, - "Vs", - fontsize=23, - verticalalignment="top", - transform=fig.transFigure, + axis.axis("off") + plot_plasma(axis, mfile, scan, colour_scheme) + plasma_main_page( + axis, + fig, + rmajor, + rminor, + triang, + radius_plasma_core_norm, + kappa, + value("i_single_null"), + value("plasma_square"), + value("big_q_plasma"), + value("p_fw_alpha_mw"), + value("f_p_alpha_plasma_deposited"), + value("p_neutron_total_mw"), + value("pflux_plasma_surface_neutron_avg_mw"), + ) + + geom_type = PlasmaGeometryModelType(value("i_plasma_geometry")) + textstr_plasma = ( + "$\\mathbf{Shaping:}$\n\n" + f"$\\kappa_{{95}}$: {value('kappa95'):.2f} ({geom_type.kappa95_model.description}) | " + f"$\\delta_{{95}}$: {value('triang95'):.2f} ({geom_type.triang95_model.description}) | " + f"$\\zeta$: {value('plasma_square'):.2f}\n" + f"$\\kappa$: {kappa:.2f} ({geom_type.kappa_model.description}) | " + f"$\\delta$: {triang:.2f} ({geom_type.triang_model.description}) | " + f"A: {value('aspect'):.2f}\n" + f"$V_{{\\text{{p}}}}$: {value('vol_plasma'):,.2f} $\\mathrm{{m}}^3$ | " + f"$A_{{\\text{{p,surface}}}}$: {value('a_plasma_surface'):,.2f} $\\mathrm{{m}}^2$ | " + f"$A_{{\\text{{p,poloidal}}}}$: {value('a_plasma_poloidal'):,.3f} $\\mathrm{{m}}^2$\n" + f"$L_{{\\text{{p,poloidal}}}}$: {value('len_plasma_poloidal'):,.3f} $\\mathrm{{m}}$" + ) + add_panel(0.365, 0.975, 0.340, 0.110, textstr_plasma, _box_style("lightyellow")) + + i_hcd_primary = value("i_hcd_primary") + i_hcd_secondary = value("i_hcd_secondary") + textstr_hcd = ( + "$\\mathbf{Heating\\ &\\ current\\ drive:}$\n\n" + f"Total injected heat: {value('p_hcd_injected_total_mw'):.3f} MW\n" + f"Ohmic heating power: {value('p_plasma_ohmic_mw'):.3f} MW\n\n" + f"$\\mathbf{{Primary\\ system: {CurrentDriveModel(i_hcd_primary).abbreviation}}}$\n" + f"Current driving power: {value('p_hcd_primary_injected_mw'):.4f} MW\n" + f"Extra heat power: {value('p_hcd_primary_extra_heat_mw'):.4f} MW\n" + f"$\\eta_{{\\text{{CD,prim}}}}$: {value('eta_cd_hcd_primary'):.4f} A/W | " + f"$\\langle\\zeta_{{\\text{{CD,prim}}}}\\rangle$: " + f"{value('eta_cd_dimensionless_hcd_primary'):.4f}\n" + f"$\\gamma_{{\\text{{CD,prim}}}}$: {value('eta_cd_norm_hcd_primary'):.4f} " + "$\\times 10^{20}\\ \\mathrm{A}/\\mathrm{Wm}^2$\n" + f"Current driven by primary: {value('c_hcd_primary_driven') / 1e6:.3f} MA\n\n" + f"$\\mathbf{{Secondary\\ system: {CurrentDriveModel(i_hcd_secondary).abbreviation}}}$\n" + f"Current driving power: {value('p_hcd_secondary_injected_mw'):.4f} MW\n" + f"Extra heat power: {value('p_hcd_secondary_extra_heat_mw'):.4f} MW\n" + f"$\\eta_{{\\text{{CD,sec}}}}$: {value('eta_cd_hcd_secondary'):.4f} A/W | " + f"$\\langle\\zeta_{{\\text{{CD,sec}}}}\\rangle$: " + f"{value('eta_cd_dimensionless_hcd_secondary'):.4f}\n" + f"$\\gamma_{{\\text{{CD,sec}}}}$: {value('eta_cd_norm_hcd_secondary'):.4f} " + "$\\times 10^{20}\\ \\mathrm{A}/\\mathrm{Wm}^2$\n" + f"Current driven by secondary: {value('c_hcd_secondary_driven') / 1e6:.3f} MA" + ) + add_panel( + 0.730, + 0.675, + 0.245, + 0.300, + textstr_hcd, + _box_style("paleturquoise") | {"edgecolor": "black"}, ) + add_symbol("$P_{\\text{inj}}$", 0.92, 0.625) - # ========================================= + beta_component = BetaComponentLimits(value("i_beta_component")).full_name + beta_norm_model = BetaNormMaxModel(value("i_beta_norm_max")).full_name + textstr_beta = ( + "$\\mathbf{Beta\\ Information:}$\n\n" + f"Total beta, $\\langle\\beta\\rangle$: {value('beta_total_vol_avg'):.4f}\n" + f"Thermal beta, $\\langle\\beta_{{\\text{{thermal}}}}\\rangle$: " + f"{value('beta_thermal_vol_avg'):.4f}\n" + f"Toroidal beta, $\\langle\\beta_{{\\text{{t}}}}\\rangle$: " + f"{value('beta_toroidal_vol_avg'):.4f}\n" + f"Poloidal beta, $\\langle\\beta_{{\\text{{p}}}}\\rangle$: " + f"{value('beta_poloidal_vol_avg'):.4f}\n" + f"Fast-alpha beta, $\\langle\\beta_{{\\alpha}}\\rangle$: " + f"{value('beta_fast_alpha'):.4f}\n" + f"Upper limit on {beta_component}, $\\langle\\beta\\rangle$: " + f"{value('beta_vol_avg_max'):.4f}\n" + f"Normalised total beta, $\\beta_{{\\text{{N}}}}$: " + f"{value('beta_norm_total'):.4f}\n" + f"Normalised thermal beta, $\\beta_{{\\text{{N,thermal}}}}$: " + f"{value('beta_norm_thermal'):.4f}\n" + f"Maximum normalised beta ({beta_norm_model}), $\\beta_{{\\text{{N,max}}}}$: " + f"{value('beta_norm_max'):.4f}" + ) + add_panel(0.025, 0.975, 0.325, 0.180, textstr_beta, _box_style("lightblue")) + add_symbol("$\\beta$", 0.27, 0.94) + + inductance_model = IndInternalNormModel( + value("i_ind_plasma_internal_norm") + ).full_name + textstr_volt_second = ( + "$\\mathbf{Volt-second\\ requirements:}$\n\n" + f"Total volt-second consumption: {value('vs_plasma_total_required'):.4f} Vs\n" + f" - Internal volt-seconds: {value('vs_plasma_internal'):.4f} Vs\n" + f" - Volt-seconds needed for burn: {value('vs_plasma_burn_required'):.4f} Vs\n" + f" - Volt-seconds needed for ramp: {value('vs_plasma_ramp_required'):.4f} Vs | " + f"$C_{{\\text{{ejima}}}}$: {value('ejima_coeff'):.4f}\n" + f"$V_{{\\text{{loop}}}}$: {value('v_plasma_loop_burn'):.4f} V\n" + f"$\\Omega_{{\\text{{p}}}}$: {value('res_plasma'):.4e} $\\Omega$\n" + f"Plasma resistive diffusion time: {value('t_plasma_res_diffusion'):,.4f} s\n" + f"Plasma inductance: {value('ind_plasma'):.4e} H | " + f"ITER $l_i(3)$: {value('ind_plasma_internal_norm_iter_3'):.4f}\n" + f"Plasma stored magnetic energy: {value('e_plasma_magnetic_stored') / 1e9:.4f} GJ\n" + f"Plasma normalised internal inductance, $l_i$ ({inductance_model}): " + f"{value('ind_plasma_internal_norm'):.3f}" + ) + add_panel(0.025, 0.780, 0.335, 0.195, textstr_volt_second, _box_style("lightgreen")) + add_symbol("Vs", 0.30, 0.77) - # Add divertor information textstr_div = ( - f"\n$P_{{\\text{{sep}}}}$: {mfile.get('p_plasma_separatrix_mw', scan=scan):.2f} MW \n" - f"$\\frac{{P_{{\\text{{sep}}}}}}{{R}}$: {mfile.get('p_plasma_separatrix_rmajor_mw', scan=scan):.2f} MW/m \n" - f"$\\frac{{P_{{\\text{{sep}}}}B_T}}{{q_{{95}} A R}}$: {mfile.get('p_div_bt_q_aspect_rmajor_mw', scan=scan):.2f} MW T/m " - ) - - axis.text( - 0.35, - 0.12, - textstr_div, - fontsize=9, - verticalalignment="top", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "orange", - "alpha": 1.0, - "linewidth": 2, - }, + f"$P_{{\\text{{sep}}}}$: {value('p_plasma_separatrix_mw'):.2f} MW\n" + f"$\\frac{{P_{{\\text{{sep}}}}}}{{R}}$: " + f"{value('p_plasma_separatrix_rmajor_mw'):.2f} MW/m\n" + f"$\\frac{{P_{{\\text{{sep}}}}B_T}}{{q_{{95}} A R}}$: " + f"{value('p_div_bt_q_aspect_rmajor_mw'):.2f} MW T/m" ) + add_panel(0.350, 0.120, 0.165, 0.095, textstr_div, _box_style("orange")) + add_symbol("$P_{\\text{div}}$", 0.45, 0.10) - # Add divertor label - axis.text( - 0.45, - 0.1, - "$P_{\\text{div}}$", - fontsize=23, - verticalalignment="top", - transform=fig.transFigure, - ) - - # ================================================ - - # Add confinement information textstr_confinement = ( - f"$\\mathbf{{Confinement:}}$\n \n" - f"Confinement scaling law: {mfile.get('tauelaw', scan=scan)}\n" - f"Confinement $H$ factor: {mfile.get('hfact', scan=scan):.4f}\n" - f"Energy confinement time from scaling: {mfile.get('t_energy_confinement', scan=scan):.4f} s\n" - f"Fusion double product: {mfile.get('ntau', scan=scan):.4e} s/m³\n" - f"Lawson Triple product: {mfile.get('nttau', scan=scan):.4e} keV·s/m³\n" - f"Transport loss power assumed in scaling law: {mfile.get('p_plasma_loss_mw', scan=scan):.4f} MW\n" - f"Plasma thermal energy (inc. $\\alpha$), $W$: {mfile.get('e_plasma_beta', scan=scan) / 1e9:.4f} GJ\n" - f"Alpha particle confinement time: {mfile.get('t_alpha_confinement', scan=scan):.4f} s | $\\tau_{{\\alpha}}/\\tau_{{e}}$: {mfile.get('f_t_alpha_energy_confinement', scan=scan):.4f}" - ) - - axis.text( + "$\\mathbf{Confinement:}$\n\n" + f"Confinement scaling law: {value('tauelaw')}\n" + f"Confinement $H$ factor: {value('hfact'):.4f}\n" + f"Energy confinement time from scaling: {value('t_energy_confinement'):.4f} s\n" + f"Fusion double product: {value('ntau'):.4e} s/m³\n" + f"Lawson Triple product: {value('nttau'):.4e} keV·s/m³\n" + f"Transport loss power assumed in scaling law: {value('p_plasma_loss_mw'):.4f} MW\n" + f"Plasma thermal energy (inc. $\\alpha$), $W$: " + f"{value('e_plasma_beta') / 1e9:.4f} GJ\n" + f"Alpha particle confinement time: {value('t_alpha_confinement'):.4f} s | " + f"$\\tau_{{\\alpha}}/\\tau_{{e}}$: {value('f_t_alpha_energy_confinement'):.4f}" + ) + add_panel( 0.025, - 0.57, + 0.570, + 0.325, + 0.155, textstr_confinement, - fontsize=9, - verticalalignment="top", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "gainsboro", # Changed to a not normal color (Aquamarine) - "alpha": 1.0, - "linewidth": 2, - "edgecolor": "black", - }, - ) - - # Add tau label - axis.text( - 0.3, - 0.55, - "$\\tau_{\\text{e}} $", - fontsize=23, - verticalalignment="top", - transform=fig.transFigure, - ) - - # ========================================= - - # Load the neutron image - with resources.path( - "process.core.io.plot.images", "alpha_particle.png" - ) as alpha_particle_image_path: - # Use importlib.resources to locate the image - alpha_particle = mpimg.imread(alpha_particle_image_path.open("rb")) - - # Display the neutron image over the figure, not the axes - new_ax = axis.inset_axes( - [0.975, 0.275, 0.075, 0.075], transform=axis.transAxes, zorder=10 - ) - new_ax.imshow(alpha_particle) - new_ax.axis("off") - - axis.annotate( - "", - xy=(rmajor + rminor, -rminor * kappa * 0.55), # Pointing at the plasma - xytext=(rmajor + 0.2 * rminor, -rminor * kappa * 0.25), - arrowprops={"facecolor": "red", "edgecolor": "grey", "lw": 1}, - ) - - textstr_alpha = ( - f"$P_{{\\alpha,\\text{{loss}}}}$ {mfile.get('p_fw_alpha_mw', scan=scan):.2f} MW \n" - f"$f_{{\\alpha,\\text{{coupled}}}}$ {mfile.get('f_p_alpha_plasma_deposited', scan=scan):.2f}" - ) - - axis.text( - 1.0, - 0.275, - textstr_alpha, - fontsize=9, - verticalalignment="top", - transform=axis.transAxes, - bbox={"boxstyle": "round", "facecolor": "red", "alpha": 1.0, "linewidth": 2}, - ) - - # ========================================= - with resources.path( - "process.core.io.plot.images", "neutron.png" - ) as neutron_image_path: - neutron = mpimg.imread(neutron_image_path.open("rb")) - new_ax = axis.inset_axes( - [0.975, 0.75, 0.075, 0.075], transform=axis.transAxes, zorder=10 - ) - new_ax.imshow(neutron) - new_ax.axis("off") - - # Draw a red arrow coming from the right and pointing at the plasma - axis.annotate( - "", - xy=(rmajor + rminor, rminor * kappa * 0.65), # Pointing at the plasma - xytext=(rmajor, rminor * kappa * 0.5), - arrowprops={"facecolor": "grey", "edgecolor": "grey", "lw": 1}, - ) - - textstr_neutron = ( - f"$P_{{\\text{{n,total}}}}$ {mfile.get('p_neutron_total_mw', scan=scan):.2f} MW \n" - f"$\\phi_{{\\text{{n,avg}}}}$ {mfile.get('pflux_plasma_surface_neutron_avg_mw', scan=scan):.3f} MW/m²" + _box_style("gainsboro") | {"edgecolor": "black"}, ) + add_symbol("$\\tau_{\\text{e}}$", 0.30, 0.55) - axis.text( - 0.775, - 0.875, - textstr_neutron, - fontsize=9, - verticalalignment="top", - transform=axis.transAxes, - bbox={"boxstyle": "round", "facecolor": "grey", "alpha": 0.8, "linewidth": 2}, - ) - - # =============================================== - - # Add fusion reaction information textstr_reactions = ( - f"$\\mathbf{{Fusion \\ Reactions:}}$\n \n" - f"Fuel mixture: \n" - f"| D: {mfile.get('f_plasma_fuel_deuterium', scan=scan):.2f} | T: {mfile.get('f_plasma_fuel_tritium', scan=scan):.2f} | 3He: {mfile.get('f_plasma_fuel_helium3', scan=scan):.2f} |\n\n" - f"Fusion Power, $P_{{\\text{{fus}}}}:$ {mfile.get('p_fusion_total_mw', scan=scan):,.2f} MW\n" - f"D-T Power, $P_{{\\text{{fus,DT}}}}:$ {mfile.get('p_dt_total_mw', scan=scan):,.2f} MW\n" - f"D-D Power, $P_{{\\text{{fus,DD}}}}:$ {mfile.get('p_dd_total_mw', scan=scan):,.2f} MW\n" - f"D-3He Power, $P_{{\\text{{fus,D3He}}}}:$ {mfile.get('p_dhe3_total_mw', scan=scan):,.2f} MW\n" - f"Alpha Power, $P_{{\\alpha}}:$ {mfile.get('p_alpha_total_mw', scan=scan):,.2f} MW" - ) - - axis.text( + "$\\mathbf{Fusion\\ Reactions:}$\n\n" + "Fuel mixture:\n" + f"| D: {value('f_plasma_fuel_deuterium'):.2f} | " + f"T: {value('f_plasma_fuel_tritium'):.2f} | " + f"3He: {value('f_plasma_fuel_helium3'):.2f} |\n\n" + f"Fusion Power, $P_{{\\text{{fus}}}}$: {value('p_fusion_total_mw'):,.2f} MW\n" + f"D-T Power, $P_{{\\text{{fus,DT}}}}$: {value('p_dt_total_mw'):,.2f} MW\n" + f"D-D Power, $P_{{\\text{{fus,DD}}}}$: {value('p_dd_total_mw'):,.2f} MW\n" + f"D-3He Power, $P_{{\\text{{fus,D3He}}}}$: {value('p_dhe3_total_mw'):,.2f} MW\n" + f"Alpha Power, $P_{{\\alpha}}$: {value('p_alpha_total_mw'):,.2f} MW" + ) + add_panel( 0.025, - 0.4, + 0.400, + 0.185, + 0.165, textstr_reactions, - fontsize=9, - verticalalignment="top", - transform=fig.transFigure, - bbox={"boxstyle": "round", "facecolor": "red", "alpha": 0.6, "linewidth": 2}, + {"boxstyle": "round", "facecolor": "red", "alpha": 0.6, "linewidth": 2}, ) - # ================================================ - - # Add fuelling information textstr_fuelling = ( - f"$\\mathbf{{Fuelling:}}$\n \n" - f"Plasma mass: {mfile.get('m_plasma', scan=scan) * 1000:.4f} g\n" - f" - Average mass of all plasma ions: {mfile.get('m_ions_total_amu', scan=scan):.3f} amu\n" - f"Fuel mass: {mfile.get('m_plasma_fuel_ions', scan=scan) * 1000:.4f} g\n" - f" - Average mass of all fuel ions: {mfile.get('m_fuel_amu', scan=scan):.3f} amu\n\n" - f"Fueling rate: {mfile.get('molflow_plasma_fuelling_required', scan=scan):.3e} nucleus-pairs/s\n" - f"Fuel burn-up rate: {mfile.get('rndfuel', scan=scan):.3e} reactions/s \n" - f"Burn-up fraction: {mfile.get('burnup', scan=scan):.4f} \n" - ) - - axis.text( + "$\\mathbf{Fuelling:}$\n\n" + f"Plasma mass: {value('m_plasma') * 1000:.4f} g\n" + f" - Average mass of all plasma ions: {value('m_ions_total_amu'):.3f} amu\n" + f"Fuel mass: {value('m_plasma_fuel_ions') * 1000:.4f} g\n" + f" - Average mass of all fuel ions: {value('m_fuel_amu'):.3f} amu\n\n" + f"Fueling rate: {value('molflow_plasma_fuelling_required'):.3e} nucleus-pairs/s\n" + f"Fuel burn-up rate: {value('rndfuel'):.3e} reactions/s\n" + f"Burn-up fraction: {value('burnup'):.4f}" + ) + add_panel( 0.025, - 0.22, + 0.220, + 0.255, + 0.175, textstr_fuelling, - fontsize=9, - verticalalignment="top", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "khaki", - "alpha": 1.0, - "linewidth": 2, - "edgecolor": "black", - }, + _box_style("khaki") | {"edgecolor": "black"}, ) - # ================================================ - - # Add ion density information impurity_data = ImpurityRadiationData() textstr_ions = ( - f" $\\mathbf{{Ion \\ to \\ electron}}$\n" - f" $\\mathbf{{relative \\ number}}$\n" - f" $\\mathbf{{densities:}}$\n \n" - f" Effective charge: {mfile.get('n_charge_plasma_effective_vol_avg', scan=scan):.3f}\n\n" + "$\\mathbf{Ion\\ to\\ electron}$\n" + "$\\mathbf{relative\\ number}$\n" + "$\\mathbf{densities:}$\n\n" + f"Effective charge: {value('n_charge_plasma_effective_vol_avg'):.3f}\n\n" + "\n".join( - f" {label.replace('_', '') + ':':<6}" - f"{mfile.get(f'f_nd_impurity_electrons({index:02d})', scan=scan):.4e}" + f"{label.replace('_', '') + ':':<6}" + f"{value(f'f_nd_impurity_electrons({index:02d})'):.4e}" for index, label in enumerate(impurity_data.imp_label[:14], start=1) ) ) - - axis.text( + add_panel( 0.805, 0.335, + 0.170, + 0.310, textstr_ions, - fontsize=9, - verticalalignment="top", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "olivedrab", - "alpha": 0.7, - "linewidth": 2, - }, + {"boxstyle": "round", "facecolor": "olivedrab", "alpha": 0.7, "linewidth": 2}, + text_x_offset=0.045, ) + add_symbol("$Z$", 0.815, 0.29) - # Add ion charge label - axis.text( - 0.815, - 0.29, - "$Z$", - fontsize=23, - verticalalignment="top", - transform=fig.transFigure, - ) - - # ================================================ - - # Add plasma current information + current_model = PlasmaCurrentModel(value("i_plasma_current")).full_name + bootstrap_model = BootstrapCurrentFractionModel( + value("i_bootstrap_current") + ).full_name + diamagnetic_model = PlasmaDiamagneticCurrentModel( + value("i_diamagnetic_current") + ).full_name textstr_currents = ( - f"$\\mathbf{{Plasma\\ currents:}}$\n\n" - f"Plasma current ({PlasmaCurrentModel(int(mfile.get('i_plasma_current', scan=scan))).full_name}): {mfile.get('plasma_current_ma', scan=scan):.4f} MA \n" - f" - Bootstrap fraction ({BootstrapCurrentFractionModel(int(mfile.get('i_bootstrap_current', scan=scan))).full_name}): {mfile.get('f_c_plasma_bootstrap', scan=scan):.4f}\n" - f" - Diamagnetic fraction ({PlasmaDiamagneticCurrentModel(int(mfile.get('i_diamagnetic_current', scan=scan))).full_name}): {mfile.get('f_c_plasma_diamagnetic', scan=scan):.4f}\n" - f" - Pfirsch-Schlüter fraction {mfile.get('f_c_plasma_pfirsch_schluter', scan=scan):.4f}\n" - f" - Auxiliary fraction {mfile.get('f_c_plasma_auxiliary', scan=scan):.4f}\n" - f" - Inductive fraction {mfile.get('f_c_plasma_inductive', scan=scan):.4f}" - ) + "$\\mathbf{Plasma\\ currents:}$\n\n" + f"Plasma current ({current_model}): {value('plasma_current_ma'):.4f} MA\n" + f" - Bootstrap fraction ({bootstrap_model}): " + f"{value('f_c_plasma_bootstrap'):.4f}\n" + f" - Diamagnetic fraction ({diamagnetic_model}): " + f"{value('f_c_plasma_diamagnetic'):.4f}\n" + f" - Pfirsch-Schlüter fraction: {value('f_c_plasma_pfirsch_schluter'):.4f}\n" + f" - Auxiliary fraction: {value('f_c_plasma_auxiliary'):.4f}\n" + f" - Inductive fraction: {value('f_c_plasma_inductive'):.4f}" + ) + add_panel(0.720, 0.975, 0.255, 0.130, textstr_currents, _box_style("#C8A2C8")) + add_symbol("$I_{\\text{p}}$", 0.93, 0.90) - axis.text( - 0.72, - 0.975, - textstr_currents, - fontsize=9, - verticalalignment="top", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "#C8A2C8", # Hex code for lilac color - "alpha": 1.0, - "linewidth": 2, - }, - ) - - # Add plasma current label - axis.text( - 0.93, - 0.9, - "$I_{\\text{p}} $", - fontsize=23, - verticalalignment="top", - transform=fig.transFigure, - ) - - # Add magnetic field information textstr_fields = ( - f"$\\mathbf{{Magnetic\\ fields:}}$\n\n" - f"Toroidal field at $R_0$, $B_{{T}}$: {mfile.get('b_plasma_toroidal_on_axis', scan=scan):.4f} T \n" - f" Ripple at outboard , $\\delta$: {mfile.get('ripple_b_tf_plasma_edge', scan=scan):.2f}% \n" - f"Surface average poloidal field, $\\langle B_{{p}}(a) \\rangle$: {mfile.get('b_plasma_surface_poloidal_average', scan=scan):.4f} T\n" - f"Total field, $B_{{tot}}$: {mfile.get('b_plasma_total', scan=scan):.4f} T \n" - f"Vertical field, $B_{{vert}}$: {mfile.get('b_plasma_vertical_required', scan=scan):.4f} T" - ) + "$\\mathbf{Magnetic\\ fields:}$\n\n" + f"Toroidal field at $R_0$, $B_T$: {value('b_plasma_toroidal_on_axis'):.4f} T\n" + f" Ripple at outboard, $\\delta$: {value('ripple_b_tf_plasma_edge'):.2f}%\n" + f"Surface average poloidal field, $\\langle B_p(a)\\rangle$: " + f"{value('b_plasma_surface_poloidal_average'):.4f} T\n" + f"Total field, $B_{{\\text{{tot}}}}$: {value('b_plasma_total'):.4f} T\n" + f"Vertical field, $B_{{\\text{{vert}}}}$: " + f"{value('b_plasma_vertical_required'):.4f} T" + ) + add_panel(0.5325, 0.140, 0.2575, 0.115, textstr_fields, _box_style("royalblue")) + add_symbol("$B$", 0.75, 0.12) - axis.text( - 0.5325, - 0.14, - textstr_fields, - fontsize=9, - verticalalignment="top", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "royalblue", - "alpha": 1.0, - "linewidth": 2, - }, - ) - - # Add magnetic field label - axis.text( - 0.75, - 0.12, - "$B$", - fontsize=23, - verticalalignment="top", - transform=fig.transFigure, - ) - - # Add radiation information textstr_radiation = ( - f" $\\mathbf{{Radiation:}}$\n\n" - f" Total radiation power {mfile.get('p_plasma_rad_mw', scan=scan):.4f} MW\n" - f" Separatrix radiation fraction {mfile.get('f_p_plasma_separatrix_rad', scan=scan):.4f}\n" - f" Core radiation power {mfile.get('p_plasma_inner_rad_mw', scan=scan):.4f} MW\n" - f" - $f_{{\\text{{core,reduce}}}}$ {mfile.get('f_p_plasma_core_rad_reduction', scan=scan):.4f}\n" - f" Edge radiation power {mfile.get('p_plasma_outer_rad_mw', scan=scan):.4f} MW\n" - f" Synchrotron radiation power {mfile.get('p_plasma_sync_mw', scan=scan):.4f} MW\n" - f" Synchrotron wall reflectivity {mfile.get('f_sync_reflect', scan=scan):.4f}" - ) - - axis.text( - 0.72, - 0.83, + "$\\mathbf{Radiation:}$\n\n" + f"Total radiation power: {value('p_plasma_rad_mw'):.4f} MW\n" + f"Separatrix radiation fraction: {value('f_p_plasma_separatrix_rad'):.4f}\n" + f"Core radiation power: {value('p_plasma_inner_rad_mw'):.4f} MW\n" + f" - $f_{{\\text{{core,reduce}}}}$: {value('f_p_plasma_core_rad_reduction'):.4f}\n" + f"Edge radiation power: {value('p_plasma_outer_rad_mw'):.4f} MW\n" + f"Synchrotron radiation power: {value('p_plasma_sync_mw'):.4f} MW\n" + f"Synchrotron wall reflectivity: {value('f_sync_reflect'):.4f}" + ) + add_panel( + 0.720, + 0.830, + 0.255, + 0.145, textstr_radiation, - fontsize=9, - verticalalignment="top", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lavender", - "alpha": 1.0, - "linewidth": 2, - "edgecolor": "black", - }, + _box_style("lavender") | {"edgecolor": "black"}, + text_x_offset=0.032, ) + add_symbol("$\\gamma$", 0.725, 0.78) - # Add radiation label - axis.text( - 0.725, - 0.78, - "$\\gamma$", - fontsize=23, - verticalalignment="top", - transform=fig.transFigure, - ) - - # Add L-H threshold information - textstr_lh = ( - f"$\\mathbf{{L-H \\ threshold:}}$\n" - f"({PlasmaConfinementTransitionModel(int(mfile.get('i_l_h_threshold', scan=scan))).full_name})\n\n" - f"$P_{{\\text{{L-H}}}}:$ {mfile.get('p_l_h_threshold_mw', scan=scan):.4f} MW\n" - ) - - # Wrap long model names to new line - model_name = PlasmaConfinementTransitionModel( - int(mfile.get("i_l_h_threshold", scan=scan)) - ).full_name + model_name = PlasmaConfinementTransitionModel(value("i_l_h_threshold")).full_name if len(model_name) > 20: model_name = "\n".join(textwrap.wrap(model_name, width=20)) - textstr_lh = ( - f"$\\mathbf{{L-H \\ threshold:}}$\n" + "$\\mathbf{L-H\\ threshold:}$\n" f"{model_name}\n\n" - f"$P_{{\\text{{L-H}}}}:$ {mfile.get('p_l_h_threshold_mw', scan=scan):.4f} MW" - ) - - axis.text( - 0.22, - 0.4, - textstr_lh, - fontsize=9, - verticalalignment="top", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "peachpuff", - "alpha": 1.0, - "linewidth": 2, - }, + f"$P_{{\\text{{L-H}}}}$: {value('p_l_h_threshold_mw'):.4f} MW" ) + add_panel(0.220, 0.400, 0.110, 0.080, textstr_lh, _box_style("peachpuff")) - # Add density limit information + density_model = DensityLimitModel(value("i_density_limit")).full_name textstr_density_limit = ( - f"$\\mathbf{{Density \\ limit:}}$\n" - f"({DensityLimitModel(int(mfile.get('i_density_limit', scan=scan))).full_name})\n" - f"$n_{{\\text{{e,limit}}}}: {mfile.get('nd_plasma_electrons_max', scan=scan):.3e} \\ m^{{-3}}$\n" - f"$f_{{\\text{{GW}}}}$: {mfile.get('f_nd_plasma_greenwald', scan=scan):.4f}" - ) - - axis.text( - 0.22, - 0.31, - textstr_density_limit, - fontsize=9, - verticalalignment="top", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "pink", - "alpha": 1.0, - "linewidth": 2, - }, + "$\\mathbf{Density\\ limit:}$\n" + f"({density_model})\n" + f"$n_{{\\text{{e,limit}}}}$: {value('nd_plasma_electrons_max'):.3e} " + "$\\mathrm{m}^{-3}$\n" + f"$f_{{\\text{{GW}}}}$: {value('f_nd_plasma_greenwald'):.4f}" ) + add_panel(0.220, 0.310, 0.130, 0.075, textstr_density_limit, _box_style("pink")) def plot_current_profiles_over_time(axis: plt.Axes, mfile: MFile, scan: int): @@ -3421,8 +3056,8 @@ def plot_system_power_profiles_over_time(axis: plt.Axes, mfile: MFile, scan: int # Add energy produced info textstr_energy = ( f"$\\mathbf{{Energy \\ Production:}}$\n\n" - f"Energy produced over whole pulse: {mfile.get('e_plant_net_electric_pulse_mj', scan=scan):,.4f} MJ \n" - f"Energy produced over whole pulse: {mfile.get('e_plant_net_electric_pulse_kwh', scan=scan):,.4f} kWh \n" + f"Energy produced over whole pulse: {mfile.get('e_plant_net_electric_pulse_mj', scan=scan):,.4f} MJ\n" + f"Energy produced over whole pulse: {mfile.get('e_plant_net_electric_pulse_kwh', scan=scan):,.4f} kWh\n" ) axis.text( @@ -3432,12 +3067,7 @@ def plot_system_power_profiles_over_time(axis: plt.Axes, mfile: MFile, scan: int fontsize=9, verticalalignment="top", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "grey", - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("grey"), ) # Add energy produced info @@ -3460,7 +3090,7 @@ def plot_system_power_profiles_over_time(axis: plt.Axes, mfile: MFile, scan: int fontsize=9, verticalalignment="top", transform=fig.transFigure, - bbox={"boxstyle": "round", "facecolor": "grey", "alpha": 1.0, "linewidth": 2}, + bbox=_box_style("grey"), ) @@ -3508,8 +3138,8 @@ def plot_cryostat( def color_key(axis: plt.Axes, mfile: MFile, scan: int, colour_scheme: Literal[1, 2]): """Function to plot the colour key""" - axis.set_ylim([0, 10]) - axis.set_xlim([0, 10]) + axis.set_ylim(0, 10) + axis.set_xlim(0, 10) axis.set_axis_off() axis.set_autoscaley_on(False) axis.set_autoscalex_on(False) @@ -3553,7 +3183,7 @@ def color_key(axis: plt.Axes, mfile: MFile, scan: int, colour_scheme: Literal[1, axis.text(x_pos, y_pos, text, ha="left", va="top", size="small") axis.add_patch( patches.Rectangle( - [x_pos + 1.5, y_pos - 0.35], + (x_pos + 1.5, y_pos - 0.35), 0.5, 0.4, lw=0 if color != "none" else 1, @@ -3796,13 +3426,13 @@ def calc_xy(rt, e=e): # --- # DEMO : Fixed ranges for comparison if demo_ranges: - axis.set_ylim([0, 20]) - axis.set_xlim([0, 20]) + axis.set_ylim(0, 20) + axis.set_xlim(0, 20) # Adaptive ranges else: - axis.set_ylim([0.0, axis.get_ylim()[1]]) - axis.set_xlim([0.0, axis.get_xlim()[1]]) + axis.set_ylim(0.0, axis.get_ylim()[1]) + axis.set_xlim(0.0, axis.get_xlim()[1]) # --- @@ -4051,14 +3681,14 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int): # Ranges # --- # DEMO : Fixed ranges for comparison - ax_main.set_xlim([0, 1]) - ax_impurity.set_xlim([0, 1]) + ax_main.set_xlim(0, 1) + ax_impurity.set_xlim(0, 1) if demo_ranges: - ax_main.set_ylim([0, 20]) + ax_main.set_ylim(0, 20) # Adaptive ranges else: - ax_main.set_ylim([0, ax_main.get_ylim()[1]]) + ax_main.set_ylim(0, ax_main.get_ylim()[1]) # Use logarithmic scale for impurity axis if any impurity values are very small impurity_data = [ imp_frac[i] * ne / 1e16 @@ -4068,7 +3698,7 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int): if impurity_data and np.min(impurity_data) / np.max(impurity_data) < 0.01: # If range spans more than 100x, use log scale ax_impurity.set_yscale("log") - ax_impurity.set_ylim([1e-3, ax_impurity.get_ylim()[1]]) + ax_impurity.set_ylim(1e-3, ax_impurity.get_ylim()[1]) if i_plasma_pedestal != 0: # Print pedestal lines @@ -4200,7 +3830,7 @@ def plot_jprofile(prof, mfile: MFile, scan: int): prof.set_ylabel(r"Current density $[kA/m^2]$") prof.set_title("$J$ profile") prof.minorticks_on() - prof.set_xlim([0, 1.0]) + prof.set_xlim(0, 1.0) rho = np.linspace(0, 1) y2 = (j_plasma_0 * (1 - rho**2) ** alphaj) / 1e3 @@ -4313,14 +3943,14 @@ def plot_t_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int): # Ranges # --- - prof.set_xlim([0, 1]) + prof.set_xlim(0, 1) # DEMO : Fixed ranges for comparison if demo_ranges: - prof.set_ylim([0, 50]) + prof.set_ylim(0, 50) # Adaptive ranges else: - prof.set_ylim([0, prof.get_ylim()[1]]) + prof.set_ylim(0, prof.get_ylim()[1]) if i_plasma_pedestal != 0: # Plot pedestal lines @@ -4422,14 +4052,14 @@ def plot_qprofile(prof, demo_ranges: bool, mfile: MFile, scan: int): # Ranges # --- - prof.set_xlim([0, 1]) + prof.set_xlim(0, 1) # DEMO : Fixed ranges for comparison if demo_ranges: - prof.set_ylim([0, 10]) + prof.set_ylim(0, 10) # Adaptive ranges else: - prof.set_ylim([0, q95 * 1.2]) + prof.set_ylim(0, q95 * 1.2) prof.text( 0.6, @@ -4490,7 +4120,7 @@ def read_imprad_data(_skiprows, data_path): if "infinite confinement" in header.content: zav = np.asarray(header.data, dtype=float) - lzdata[i] = np.column_stack((Te, lz, zav)) + lzdata[i] = np.column_stack([Te, lz, zav]) # then switch string to floats return np.array(lzdata, dtype=float) @@ -4501,12 +4131,10 @@ def profiles_with_pedestal(mfile, scan: int): alphan = mfile.get("alphan", scan=scan) alphat = mfile.get("alphat", scan=scan) nd_plasma_electron_on_axis = mfile.get("nd_plasma_electron_on_axis", scan=scan) - temp_plasma_electron_on_axis_kev = mfile.get( - "temp_plasma_electron_on_axis_kev", scan=scan - ) - radius_plasma_pedestal_density_norm = mfile.get( - "radius_plasma_pedestal_density_norm", scan=scan + temp_plasma_electron_on_axis_kev = mfile.get( + "temp_plasma_electron_on_axis_kev", scan=scan ) + radius_plasma_pedestal_temp_norm = mfile.get( "radius_plasma_pedestal_temp_norm", scan=scan ) @@ -4514,7 +4142,6 @@ def profiles_with_pedestal(mfile, scan: int): n_plasma_profile_elements = int(mfile.get("n_plasma_profile_elements", scan=scan)) i_plasma_pedestal = mfile.get("i_plasma_pedestal", scan=scan) nd_plasma_pedestal_electron = mfile.get("nd_plasma_pedestal_electron", scan=scan) - ne0 = mfile.get("nd_plasma_electron_on_axis", scan=scan) radius_plasma_pedestal_density_norm = mfile.get( "radius_plasma_pedestal_density_norm", scan=scan ) @@ -4685,16 +4312,16 @@ def plot_line_brem_power_density_profile( # Ranges # --- axis.legend(loc="upper left", bbox_to_anchor=(-0.1, -0.1), ncol=4) - axis.set_xlim([0, 1.0]) + axis.set_xlim(0, 1.0) axis.set_yscale("log") axis.yaxis.grid(True, which="both", alpha=0.2) # DEMO : Fixed ranges for comparison if demo_ranges: - axis.set_ylim([1e-4, 0.5]) + axis.set_ylim(1e-4, 0.5) # Adaptive ranges else: - axis.set_ylim([1e-4, axis.get_ylim()[1]]) + axis.set_ylim(1e-4, axis.get_ylim()[1]) # --- @@ -4847,7 +4474,7 @@ def plot_line_brem_power_profile( # Ranges # --- - axis.set_xlim([0, 1.0]) + axis.set_xlim(0, 1.0) axis.legend(loc="upper left", bbox_to_anchor=(1.0, 1.0), ncol=1) axis.set_yscale("log") axis.yaxis.grid(True, which="both", alpha=0.2) @@ -4926,12 +4553,12 @@ def plot_line_brem_loss_function_profile( axis.set_xlabel(r"$\rho \quad [r/a]$") axis.set_ylabel(r"$L_z$ $[\mathrm{W}\mathrm{m}^3]$") axis.set_title("Line & Bremsstrahlung Loss Function ($L_z$) Profiles") - axis.set_xlim([0, 1.0]) + axis.set_xlim(0, 1.0) axis.set_yscale("log") axis.yaxis.grid(True, which="both", alpha=0.2) -def plot_rad_density_contour(axis: "mpl.axes.Axes", mfile: "Any", scan: int, impp: str): +def plot_rad_density_contour(axis: Axes, mfile: MFile, scan: int, impp: str): """Plots the contour of line and bremsstrahlung radiation density [MW/m³] for a plasma cross-section. @@ -5629,8 +5256,8 @@ def plot_first_wall_top_down_cross_section(axis: plt.Axes, mfile: MFile, scan: i axis.set_xlabel("X [cm]") axis.set_ylabel("R [cm]") axis.set_title("First Wall Top-Down Cross Section") - axis.set_xlim([-1, 2 * dx_fw_module + 1]) - axis.set_ylim([-1, 2 * (dr_fw_wall + radius_fw_channel) + 1]) + axis.set_xlim(-1, 2 * dx_fw_module + 1) + axis.set_ylim(-1, 2 * (dr_fw_wall + radius_fw_channel) + 1) def plot_first_wall_poloidal_cross_section(axis: plt.Axes, mfile: MFile, scan: int): @@ -5751,8 +5378,8 @@ def plot_first_wall_poloidal_cross_section(axis: plt.Axes, mfile: MFile, scan: i axis.set_xlabel("R [m]") axis.set_ylabel("Z [m]") axis.set_title("First Wall Poloidal Cross Section") - axis.set_xlim([-0.01, (dx_fw_module + radius_fw_channel * 2) + 0.01]) - axis.set_ylim([-0.2, len_fw_channel + 0.2]) + axis.set_xlim(-0.01, (dx_fw_module + radius_fw_channel * 2) + 0.01) + axis.set_ylim(-0.2, len_fw_channel + 0.2) def plot_firstwall( @@ -6057,7 +5684,7 @@ def plot_superconducting_tf_wp(axis: plt.Axes, mfile: MFile, scan: int, fig): if i_tf_sup == 1: axis.add_patch( Circle( - [0, 0], + (0, 0), r_tf_inboard_in, facecolor="none", edgecolor="black", @@ -6068,7 +5695,7 @@ def plot_superconducting_tf_wp(axis: plt.Axes, mfile: MFile, scan: int, fig): if i_tf_case_geom == TFPlasmaCaseType.CIRCULAR: axis.add_patch( Circle( - [0, 0], + (0, 0), r_tf_inboard_out, facecolor="none", edgecolor="black", @@ -6576,7 +6203,7 @@ def plot_superconducting_tf_wp(axis: plt.Axes, mfile: MFile, scan: int, fig): # Add info about the steel casing surrounding the WP textstr_casing = ( - f"$\\mathbf{{Casing:}}$\n \n" + f"$\\mathbf{{Casing:}}$\n\n" f"Coil half angle: {mfile.get('rad_tf_coil_inboard_toroidal_half', scan=scan):.3f} radians\n\n" f"$\\text{{Full Coil Case:}}$\n" f"$r_{{start}} \\rightarrow r_{{end}}$: {mfile.get('r_tf_inboard_in', scan=scan):.3f} $\\rightarrow$ {mfile.get('r_tf_inboard_out', scan=scan):.3f} m\n" @@ -6613,7 +6240,7 @@ def plot_superconducting_tf_wp(axis: plt.Axes, mfile: MFile, scan: int, fig): # Add info about the steel casing surrounding the WP textstr_wp_insulation = ( - f"$\\mathbf{{Ground \\ Insulation:}}$\n \n" + f"$\\mathbf{{Ground \\ Insulation:}}$\n\n" f"Area of insulation around WP: {mfile.get('a_tf_wp_ground_insulation', scan=scan):.3f} $\\mathrm{{m}}^2$\n" f"$\\Delta r$: {mfile.get('dx_tf_wp_insulation', scan=scan):.4f} m\n\n" f"WP Insertion Gap:\n" @@ -6637,7 +6264,7 @@ def plot_superconducting_tf_wp(axis: plt.Axes, mfile: MFile, scan: int, fig): # Add info about the Winding Pack textstr_wp = ( - f"$\\mathbf{{Winding \\ Pack:}}$\n \n" + f"$\\mathbf{{Winding \\ Pack:}}$\n\n" f"$N_{{\\text{{turns}}}}$: " f"{int(mfile.get('n_tf_coil_turns', scan=scan))} turns\n" f"$r_{{start}} \\rightarrow r_{{end}}$: {mfile.get('r_tf_wp_inboard_inner', scan=scan):.3f} $\\rightarrow$ {mfile.get('r_tf_wp_inboard_outer', scan=scan):.3f} m\n" @@ -6674,7 +6301,7 @@ def plot_superconducting_tf_wp(axis: plt.Axes, mfile: MFile, scan: int, fig): # Add info about the Winding Pack textstr_general_info = ( - f"$\\mathbf{{General \\ info:}}$\n \n" + f"$\\mathbf{{General \\ info:}}$\n\n" f"$N_{{\\text{{TF,coil}}}}$: {mfile.get('n_tf_coils', scan=scan)}\n" f"Self inductance of single coil: {mfile.get('ind_tf_coil', scan=scan) * 1e6:.4f} $\\mu$H\n" f"Stored energy of all coils: {mfile.get('e_tf_magnetic_stored_total_gj', scan=scan):.4f} GJ\n" @@ -6736,7 +6363,7 @@ def plot_resistive_tf_wp(axis: plt.Axes, mfile: MFile, scan: int, fig): axis.add_patch( Circle( - [0, 0], + (0, 0), r_tf_inboard_in, facecolor="none", edgecolor="black", @@ -6747,7 +6374,7 @@ def plot_resistive_tf_wp(axis: plt.Axes, mfile: MFile, scan: int, fig): if i_tf_case_geom == TFPlasmaCaseType.CIRCULAR: axis.add_patch( Circle( - [0, 0], + (0, 0), r_tf_inboard_out, facecolor="none", edgecolor="black", @@ -7100,7 +6727,7 @@ def plot_resistive_tf_info(axis: plt.Axes, mfile: MFile, scan: int, fig): """Plot info about the resistive TF coils""" # Add info about the steel casing surrounding the WP textstr_casing = ( - f"$\\mathbf{{Casing:}}$\n \n" + f"$\\mathbf{{Casing:}}$\n\n" f"Coil half angle: {mfile.get('rad_tf_coil_inboard_toroidal_half', scan=scan):.3f} radians\n\n" f"$\\text{{Full Coil Case:}}$\n" f"$r_{{start}} \\rightarrow r_{{end}}$: {mfile.get('r_tf_inboard_in', scan=scan):.3f} $\\rightarrow$ {mfile.get('r_tf_inboard_out', scan=scan):.3f} m\n" @@ -7123,12 +6750,12 @@ def plot_resistive_tf_info(axis: plt.Axes, mfile: MFile, scan: int, fig): verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={"boxstyle": "round", "facecolor": "grey", "alpha": 1.0, "linewidth": 2}, + bbox=_box_style("grey"), ) # Add info about the steel casing surrounding the WP textstr_wp_insulation = ( - f"$\\mathbf{{Insulation:}}$\n \n" + f"$\\mathbf{{Insulation:}}$\n\n" f"Area of insulation around WP: {mfile.get('a_tf_wp_ground_insulation', scan=scan):.3f} $\\mathrm{{m}}^2$\n" f"$\\Delta r$: {mfile.get('dx_tf_wp_insulation', scan=scan):.4f} m\n\n" f"$\\text{{Turn Insulation:}}$\n" @@ -7147,7 +6774,7 @@ def plot_resistive_tf_info(axis: plt.Axes, mfile: MFile, scan: int, fig): # Add info about the Winding Pack textstr_wp = ( - f"$\\mathbf{{Winding Pack:}}$\n \n" + f"$\\mathbf{{Winding Pack:}}$\n\n" f"$N_{{\\text{{turns}}}}$: " f"{int(mfile.get('n_tf_coil_turns', scan=scan))} turns\n" f"$r_{{start}} \\rightarrow r_{{end}}$: {mfile.get('r_tf_wp_inboard_inner', scan=scan):.3f} $\\rightarrow$ {mfile.get('r_tf_wp_inboard_outer', scan=scan):.3f} m\n" @@ -7172,7 +6799,7 @@ def plot_resistive_tf_info(axis: plt.Axes, mfile: MFile, scan: int, fig): # Add info about the Winding Pack textstr_general_info = ( - f"$\\mathbf{{General \\ info:}}$\n \n" + f"$\\mathbf{{General \\ info:}}$\n\n" f"Self inductance: {mfile.get('ind_tf_coil', scan=scan) * 1e6:.4f} $\\mu$H\n" f"Stored energy of all coils: {mfile.get('e_tf_magnetic_stored_total_gj', scan=scan):.4f} GJ\n" ) @@ -7189,7 +6816,7 @@ def plot_resistive_tf_info(axis: plt.Axes, mfile: MFile, scan: int, fig): # Add info about the Winding Pack textstr_cooling = ( - f"$\\mathbf{{Cooling \\ info:}}$\n \n" + f"$\\mathbf{{Cooling \\ info:}}$\n\n" f"Coolant inlet temperature: {mfile.get('temp_cp_coolant_inlet', scan=scan):.2f} K\n" f"Coolant temperature rise: {mfile.get('dtemp_cp_coolant', scan=scan):.2f} K\n" f"Coolant velocity: {mfile.get('vel_cp_coolant_midplane', scan=scan):.2f} $\\mathrm{{ms^{{-1}}}}$\n\n" @@ -7209,12 +6836,7 @@ def plot_resistive_tf_info(axis: plt.Axes, mfile: MFile, scan: int, fig): verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "wheat", - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("wheat"), ) @@ -7439,7 +7061,7 @@ def _pack_strands_rectangular_with_obstacles( if TFWPIntegerTurnType(i_tf_turns_integer) == TFWPIntegerTurnType.NON_INTEGER: axis.add_patch( Rectangle( - [0, 0], + (0, 0), turn_width, turn_width, facecolor="red", @@ -7449,7 +7071,7 @@ def _pack_strands_rectangular_with_obstacles( # Plot the steel conduit axis.add_patch( Rectangle( - [insulation_thickness, insulation_thickness], + (insulation_thickness, insulation_thickness), (turn_width - 2 * insulation_thickness), (turn_width - 2 * insulation_thickness), facecolor="grey", @@ -7460,10 +7082,10 @@ def _pack_strands_rectangular_with_obstacles( # Plot the cable space with rounded corners axis.add_patch( patches.FancyBboxPatch( - [ + ( insulation_thickness + steel_thickness, insulation_thickness + steel_thickness, - ], + ), (turn_width - 2 * (insulation_thickness + steel_thickness)), (turn_width - 2 * (insulation_thickness + steel_thickness)), boxstyle=patches.BoxStyle( @@ -7477,10 +7099,10 @@ def _pack_strands_rectangular_with_obstacles( # Plot dashed line around the cable space axis.add_patch( Rectangle( - [ + ( insulation_thickness + steel_thickness, insulation_thickness + steel_thickness, - ], + ), (turn_width - 2 * (insulation_thickness + steel_thickness)), (turn_width - 2 * (insulation_thickness + steel_thickness)), facecolor="none", @@ -7493,7 +7115,7 @@ def _pack_strands_rectangular_with_obstacles( # Plot the coolant channel axis.add_patch( Circle( - [(turn_width / 2), (turn_width / 2)], + ((turn_width / 2), (turn_width / 2)), he_pipe_diameter / 2, facecolor="white", edgecolor="black", @@ -7544,7 +7166,7 @@ def _pack_strands_rectangular_with_obstacles( elif TFWPIntegerTurnType(i_tf_turns_integer) == TFWPIntegerTurnType.INTEGER: axis.add_patch( Rectangle( - [0, 0], + (0, 0), turn_width, turn_height, facecolor="red", @@ -7555,7 +7177,7 @@ def _pack_strands_rectangular_with_obstacles( # Plot the steel conduit axis.add_patch( Rectangle( - [insulation_thickness, insulation_thickness], + (insulation_thickness, insulation_thickness), (turn_width - 2 * insulation_thickness), (turn_height - 2 * insulation_thickness), facecolor="grey", @@ -7566,10 +7188,10 @@ def _pack_strands_rectangular_with_obstacles( # Plot the cable space with rounded corners axis.add_patch( patches.FancyBboxPatch( - [ + ( insulation_thickness + steel_thickness, insulation_thickness + steel_thickness, - ], + ), (turn_width - 2 * (insulation_thickness + steel_thickness)), (turn_height - 2 * (insulation_thickness + steel_thickness)), boxstyle=patches.BoxStyle( @@ -7582,10 +7204,10 @@ def _pack_strands_rectangular_with_obstacles( # Plot dashed line around the cable space axis.add_patch( Rectangle( - [ + ( insulation_thickness + steel_thickness, insulation_thickness + steel_thickness, - ], + ), (turn_width - 2 * (insulation_thickness + steel_thickness)), (turn_height - 2 * (insulation_thickness + steel_thickness)), facecolor="none", @@ -7597,7 +7219,7 @@ def _pack_strands_rectangular_with_obstacles( ) axis.add_patch( Circle( - [(turn_width / 2), (turn_height / 2)], + ((turn_width / 2), (turn_height / 2)), he_pipe_diameter / 2, facecolor="white", edgecolor="black", @@ -7656,12 +7278,7 @@ def _pack_strands_rectangular_with_obstacles( verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "red", - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("red"), ) # Add info about the steel casing surrounding the WP @@ -7678,12 +7295,7 @@ def _pack_strands_rectangular_with_obstacles( verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "grey", - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("grey"), ) if TFWPIntegerTurnType(i_tf_turns_integer) == TFWPIntegerTurnType.NON_INTEGER: @@ -7692,15 +7304,15 @@ def _pack_strands_rectangular_with_obstacles( f"$\\mathbf{{Cable \\ Space:}}$\n\n" f"$\\Delta r:$ {cable_space_width:.3e} m\n" f"Corner radius, $r$: {radius_tf_turn_cable_space_corners:.3e} m\n" - f"Cable area with no cooling \nchannel or gaps: {a_tf_turn_cable_space_no_void:.3e} m$^2$\n" + f"Cable area with no cooling\nchannel or gaps: {a_tf_turn_cable_space_no_void:.3e} m$^2$\n" f"Extra cable space area void fraction: {f_a_tf_turn_cable_space_extra_void}\n" f"True cable space area: {a_tf_turn_cable_space_effective:.3e} m$^2$" ) elif TFWPIntegerTurnType(i_tf_turns_integer) == TFWPIntegerTurnType.INTEGER: textstr_turn_cable_space = ( f"$\\mathbf{{Cable \\ Space:}}$\n\n" - f"Cable space: \n$\\Delta r$: {cable_space_width_radial:.3e} m \n" - f"$\\Delta x$: {cable_space_width_toroidal:.3e} m \n" + f"Cable space:\n$\\Delta r$: {cable_space_width_radial:.3e} m\n" + f"$\\Delta x$: {cable_space_width_toroidal:.3e} m\n" f"Corner radius, $r$: {radius_tf_turn_cable_space_corners:.3e} m\n" f"Cable area with no cooling channel or gaps: {a_tf_turn_cable_space_no_void:.3e} m$^2$\n" f"Extra cable space area void fraction: {f_a_tf_turn_cable_space_extra_void}\n" @@ -7715,12 +7327,7 @@ def _pack_strands_rectangular_with_obstacles( verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "royalblue", - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("royalblue"), ) if TFWPIntegerTurnType(i_tf_turns_integer) == TFWPIntegerTurnType.NON_INTEGER: @@ -7745,12 +7352,7 @@ def _pack_strands_rectangular_with_obstacles( verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "wheat", - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("wheat"), ) # Add info about the steel casing surrounding the WP @@ -7768,32 +7370,27 @@ def _pack_strands_rectangular_with_obstacles( verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "white", - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("white"), ) textstr_superconductor = ( - f"$\\mathbf{{Superconductor:}}$\n \n" - f"Superconductor used: \n" + f"$\\mathbf{{Superconductor:}}$\n\n" + f"Superconductor used:\n" f"{SuperconductorModel(mfile.get('i_tf_sc_mat', scan=scan)).full_name}\n" - f"Critical field at zero \ntemperature and strain: {mfile.get('b_tf_superconductor_critical_zero_temp_strain', scan=scan):.4f} T\n" - f"Critical temperature at \nzero field and strain: {mfile.get('temp_tf_superconductor_critical_zero_field_strain', scan=scan):.4f} K\n" + f"Critical field at zero\ntemperature and strain: {mfile.get('b_tf_superconductor_critical_zero_temp_strain', scan=scan):.4f} T\n" + f"Critical temperature at\nzero field and strain: {mfile.get('temp_tf_superconductor_critical_zero_field_strain', scan=scan):.4f} K\n" f"Temperature at conductor: {mfile.get('tftmp', scan=scan):.4f} K\n" f"Field at conductor: {mfile.get('b_tf_inboard_peak_with_ripple', scan=scan):.4f} T\n" - f"Superconductor critical current density at \noperating conditions: {mfile.get('j_tf_superconductor_critical', scan=scan):.2e} A/m$^2$\n" + f"Superconductor critical current density at\noperating conditions: {mfile.get('j_tf_superconductor_critical', scan=scan):.2e} A/m$^2$\n" f"$I_{{\\text{{TF,turn critical}}}}$: {mfile.get('c_turn_cables_critical', scan=scan):,.2f} A\n" f"$I_{{\\text{{TF,turn}}}}$: {mfile.get('c_tf_turn', scan=scan):,.2f} A\n" f"Critcal current ratio: {mfile.get('f_c_tf_turn_operating_critical', scan=scan):,.4f}\n" - f"Superconductor temperature \nmargin: {mfile.get('temp_tf_superconductor_margin', scan=scan):,.4f} K\n" - f"\n$\\mathbf{{Quench:}}$\n \n" + f"Superconductor temperature\nmargin: {mfile.get('temp_tf_superconductor_margin', scan=scan):,.4f} K\n" + f"\n$\\mathbf{{Quench:}}$\n\n" f"Quench dump time: {mfile.get('t_tf_superconductor_quench', scan=scan):.4f} s\n" f"Quench detection time: {mfile.get('t_tf_quench_detection', scan=scan):.4f} s\n" - f"User input max temperature \nduring quench: {mfile.get('temp_tf_conductor_quench_max', scan=scan):.2f} K\n" - f"Required maxium WP current \ndensity for heat protection:\n{mfile.get('j_tf_wp_quench_heat_max', scan=scan):.2e} A/m$^2$\n" + f"User input max temperature\nduring quench: {mfile.get('temp_tf_conductor_quench_max', scan=scan):.2f} K\n" + f"Required maxium WP current\ndensity for heat protection:\n{mfile.get('j_tf_wp_quench_heat_max', scan=scan):.2e} A/m$^2$\n" ) axis.text( 0.75, @@ -7803,12 +7400,7 @@ def _pack_strands_rectangular_with_obstacles( verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "#6dd3f7", # light blue for superconductors - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("#6dd3f7"), # light blue for superconductors ) @@ -7852,7 +7444,7 @@ def plot_tf_croco_turn(axis: plt.Axes, fig, mfile: MFile, scan: int): if TFWPIntegerTurnType(i_tf_turns_integer) == TFWPIntegerTurnType.NON_INTEGER: axis.add_patch( Rectangle( - [0, 0], + (0, 0), turn_width, turn_width, facecolor="red", @@ -7862,7 +7454,7 @@ def plot_tf_croco_turn(axis: plt.Axes, fig, mfile: MFile, scan: int): # Plot the steel conduit axis.add_patch( Rectangle( - [insulation_thickness, insulation_thickness], + (insulation_thickness, insulation_thickness), (turn_width - 2 * insulation_thickness), (turn_width - 2 * insulation_thickness), facecolor="grey", @@ -7870,27 +7462,20 @@ def plot_tf_croco_turn(axis: plt.Axes, fig, mfile: MFile, scan: int): ), ) - # Plot the central cable space - axis.add_patch( - Circle( - [(turn_width / 2), (turn_width / 2)], - 1.5 * dia_tf_turn_croco_cable, - facecolor="white", - edgecolor="black", - linewidth=1.2, - ), - ) - - # PLot the central copper cyclinder - axis.add_patch( - Circle( - [(turn_width / 2), (turn_width / 2)], - dia_tf_turn_croco_cable / 2, - facecolor="#B87333", - edgecolor="black", - linewidth=1.2, - ), - ) + # Plot the central cable space and copper cylinder + for rad, col in [ + (1.5 * dia_tf_turn_croco_cable, "white"), + (dia_tf_turn_croco_cable / 2, "#B87333"), + ]: + axis.add_patch( + Circle( + ((turn_width / 2), (turn_width / 2)), + rad, + facecolor=col, + edgecolor="black", + linewidth=1.2, + ), + ) # Plot six surrounding Croco cables in a hexagonal layout. center_x = turn_width / 2 @@ -7937,12 +7522,7 @@ def plot_tf_croco_turn(axis: plt.Axes, fig, mfile: MFile, scan: int): verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "red", - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("red"), ) # Add info about the steel casing surrounding the WP @@ -7959,12 +7539,7 @@ def plot_tf_croco_turn(axis: plt.Axes, fig, mfile: MFile, scan: int): verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "grey", - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("grey"), ) if TFWPIntegerTurnType(i_tf_turns_integer) == TFWPIntegerTurnType.NON_INTEGER: @@ -7973,15 +7548,15 @@ def plot_tf_croco_turn(axis: plt.Axes, fig, mfile: MFile, scan: int): f"$\\mathbf{{Cable \\ Space:}}$\n\n" f"$\\Delta r:$ {cable_space_width:.3e} m\n" f"Corner radius, $r$: {radius_tf_turn_cable_space_corners:.3e} m\n" - f"Cable area with no cooling \nchannel or gaps: {a_tf_turn_cable_space_no_void:.3e} m$^2$\n" + f"Cable area with no cooling\nchannel or gaps: {a_tf_turn_cable_space_no_void:.3e} m$^2$\n" f"Extra cable space area void fraction: {f_a_tf_turn_cable_space_extra_void}\n" f"True cable space area: {a_tf_turn_cable_space_effective:.3e} m$^2$" ) elif TFWPIntegerTurnType(i_tf_turns_integer) == TFWPIntegerTurnType.INTEGER: textstr_turn_cable_space = ( f"$\\mathbf{{Cable \\ Space:}}$\n\n" - f"Cable space: \n$\\Delta r$: {cable_space_width_radial:.3e} m \n" - f"$\\Delta x$: {cable_space_width_toroidal:.3e} m \n" + f"Cable space:\n$\\Delta r$: {cable_space_width_radial:.3e} m\n" + f"$\\Delta x$: {cable_space_width_toroidal:.3e} m\n" f"Corner radius, $r$: {radius_tf_turn_cable_space_corners:.3e} m\n" f"Cable area with no cooling channel or gaps: {a_tf_turn_cable_space_no_void:.3e} m$^2$\n" f"Extra cable space area void fraction: {f_a_tf_turn_cable_space_extra_void}\n" @@ -7996,12 +7571,7 @@ def plot_tf_croco_turn(axis: plt.Axes, fig, mfile: MFile, scan: int): verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "royalblue", - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("royalblue"), ) if TFWPIntegerTurnType(i_tf_turns_integer) == TFWPIntegerTurnType.NON_INTEGER: @@ -8026,12 +7596,7 @@ def plot_tf_croco_turn(axis: plt.Axes, fig, mfile: MFile, scan: int): verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "wheat", - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("wheat"), ) # Add info about the steel casing surrounding the WP @@ -8049,31 +7614,26 @@ def plot_tf_croco_turn(axis: plt.Axes, fig, mfile: MFile, scan: int): verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "white", - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("white"), ) textstr_superconductor = ( - f"$\\mathbf{{Superconductor:}}$\n \n" + f"$\\mathbf{{Superconductor:}}$\n\n" f"Superconductor used: {SuperconductorModel(mfile.get('i_tf_sc_mat', scan=scan)).full_name}\n" - f"Critical field at zero \ntemperature and strain: {mfile.get('b_tf_superconductor_critical_zero_temp_strain', scan=scan):.4f} T\n" - f"Critical temperature at \nzero field and strain: {mfile.get('temp_tf_superconductor_critical_zero_field_strain', scan=scan):.4f} K\n" + f"Critical field at zero\ntemperature and strain: {mfile.get('b_tf_superconductor_critical_zero_temp_strain', scan=scan):.4f} T\n" + f"Critical temperature at\nzero field and strain: {mfile.get('temp_tf_superconductor_critical_zero_field_strain', scan=scan):.4f} K\n" f"Temperature at conductor: {mfile.get('tftmp', scan=scan):.4f} K\n" f"Field at conductor: {mfile.get('b_tf_inboard_peak_with_ripple', scan=scan):.4f} T\n" - f"Superconductor critical current density at \noperating conditions: {mfile.get('j_tf_superconductor_critical', scan=scan):.2e} A/m$^2$\n" + f"Superconductor critical current density at\noperating conditions: {mfile.get('j_tf_superconductor_critical', scan=scan):.2e} A/m$^2$\n" f"$I_{{\\text{{TF,turn critical}}}}$: {mfile.get('c_turn_cables_critical', scan=scan):,.2f} A\n" f"$I_{{\\text{{TF,turn}}}}$: {mfile.get('c_tf_turn', scan=scan):,.2f} A\n" f"Critcal current ratio: {mfile.get('f_c_tf_turn_operating_critical', scan=scan):,.4f}\n" - f"Superconductor temperature \nmargin: {mfile.get('temp_tf_superconductor_margin', scan=scan):,.4f} K\n" - f"\n$\\mathbf{{Quench:}}$\n \n" + f"Superconductor temperature\nmargin: {mfile.get('temp_tf_superconductor_margin', scan=scan):,.4f} K\n" + f"\n$\\mathbf{{Quench:}}$\n\n" f"Quench dump time: {mfile.get('t_tf_superconductor_quench', scan=scan):.4e} s\n" f"Quench detection time: {mfile.get('t_tf_quench_detection', scan=scan):.4e} s\n" - f"User input max temperature \nduring quench: {mfile.get('temp_tf_conductor_quench_max', scan=scan):.2f} K\n" - f"Required maxium WP current \ndensity for heat protection:\n{mfile.get('j_tf_wp_quench_heat_max', scan=scan):.2e} A/m$^2$\n" + f"User input max temperature\nduring quench: {mfile.get('temp_tf_conductor_quench_max', scan=scan):.2f} K\n" + f"Required maxium WP current\ndensity for heat protection:\n{mfile.get('j_tf_wp_quench_heat_max', scan=scan):.2e} A/m$^2$\n" ) axis.text( 0.75, @@ -8083,12 +7643,7 @@ def plot_tf_croco_turn(axis: plt.Axes, fig, mfile: MFile, scan: int): verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "#6dd3f7", # light blue for superconductors - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("#6dd3f7"), ) @@ -8151,7 +7706,7 @@ def plot_cable_in_conduit_cable(axis: plt.Axes, fig, mfile: MFile, scan: int): ) textstr_cable = ( - f"$\\mathbf{{Cable:}}$\n \n" + f"$\\mathbf{{Cable:}}$\n\n" f"Cable diameter: {cable_diameter_mm:,.4f} mm\n" f"Copper area fraction: {mfile.get('f_a_tf_turn_cable_copper', scan=scan):.4f}\n" f"Number of strands per turn: {int(mfile.get('n_tf_turn_superconducting_cables', scan=scan)):,}\n" @@ -8166,12 +7721,7 @@ def plot_cable_in_conduit_cable(axis: plt.Axes, fig, mfile: MFile, scan: int): verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "#cccccc", # grayish color - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("#cccccc"), ) axis.set_aspect("equal") @@ -8264,8 +7814,12 @@ def plot_pf_coils( ) # Get axis height for fontsize scaling - bbox = axis.get_window_extent().transformed(axis.figure.dpi_scale_trans.inverted()) - axis_height = bbox.height + axis_height = ( + axis + .get_window_extent() + .transformed(axis.figure.dpi_scale_trans.inverted()) + .height + ) for i in range(len(coils_r)): mirrored_r_points = [x_scale * r for r in r_points[i]] @@ -8389,16 +7943,7 @@ def plot_header(axis: plt.Axes, mfile: MFile, scan: int): scan : scan number to use """ - xmin = 0 - xmax = 1 - ymin = -16 - ymax = 1 - - axis.set_ylim([ymin, ymax]) - axis.set_xlim([xmin, xmax]) - axis.set_axis_off() - axis.set_autoscaley_on(False) - axis.set_autoscalex_on(False) + _setup_axis(axis, xmin=0, xmax=1, ymin=-16, ymax=1) data2 = [ (f"!{mfile.get('runtitle', scan=-1)}", "Run title", ""), @@ -8511,17 +8056,8 @@ def plot_geometry_info(axis: plt.Axes, mfile: MFile, scan: int): scan : scan number to use """ - xmin = 0 - xmax = 1 - ymin = -16 - ymax = 1 - axis.text(-0.05, 1, "Geometry:", ha="left", va="center") - axis.set_ylim([ymin, ymax]) - axis.set_xlim([xmin, xmax]) - axis.set_axis_off() - axis.set_autoscaley_on(False) - axis.set_autoscalex_on(False) + _setup_axis(axis, xmin=0, xmax=1, ymin=-16, ymax=1) in_blanket_thk = mfile.get("dr_shld_inboard", scan=scan) + mfile.get( "dr_blkt_inboard", scan=scan @@ -8560,17 +8096,8 @@ def plot_physics_info(axis: plt.Axes, mfile: MFile, scan: int): scan : scan number to use """ - xmin = 0 - xmax = 1 - ymin = -16 - ymax = 1 - axis.text(-0.05, 1, "Physics:", ha="left", va="center") - axis.set_ylim([ymin, ymax]) - axis.set_xlim([xmin, xmax]) - axis.set_axis_off() - axis.set_autoscaley_on(False) - axis.set_autoscalex_on(False) + _setup_axis(axis, xmin=0, xmax=1, ymin=-16, ymax=1) nong = mfile.get("nd_plasma_electron_line", scan=scan) / mfile.get( "nd_plasma_electron_max_array(7)", scan=scan @@ -8640,17 +8167,8 @@ def plot_magnetics_info(axis: plt.Axes, mfile: MFile, scan: int): # Check for Copper magnets i_tf_sup = int(mfile.get("i_tf_sup", scan=scan)) if "i_tf_sup" in mfile.data else 1 - xmin = 0 - xmax = 1 - ymin = -16 - ymax = 1 - axis.text(-0.05, 1, "Coil currents etc:", ha="left", va="center") - axis.set_ylim([ymin, ymax]) - axis.set_xlim([xmin, xmax]) - axis.set_axis_off() - axis.set_autoscaley_on(False) - axis.set_autoscalex_on(False) + _setup_axis(axis, xmin=0, xmax=1, ymin=-16, ymax=1) # Number of coils (1 is OH coil) number_of_coils = 0 @@ -8676,22 +8194,21 @@ def plot_magnetics_info(axis: plt.Axes, mfile: MFile, scan: int): t_plant_pulse_burn = mfile.get("t_plant_pulse_burn", scan=scan) / 3600.0 - if "i_tf_bucking" in mfile.data: - i_tf_bucking = int(mfile.get("i_tf_bucking", scan=scan)) - else: - i_tf_bucking = 1 + i_tf_bucking = ( + int(mfile.get("i_tf_bucking", scan=scan)) if "i_tf_bucking" in mfile.data else 1 + ) # Get superconductor material (i_tf_sc_mat) # If i_tf_sc_mat not present, assume resistive - if "i_tf_sc_mat" in mfile.data: - i_tf_sc_mat = int(mfile.get("i_tf_sc_mat", scan=scan)) - else: - i_tf_sc_mat = 0 + i_tf_sc_mat = ( + int(mfile.get("i_tf_sc_mat", scan=scan)) if "i_tf_sc_mat" in mfile.data else 0 + ) - if i_tf_sc_mat > 0: - tftype = SuperconductorModel(int(mfile.get("i_tf_sc_mat", scan=scan))).full_name - else: - tftype = "Resistive Copper" + tftype = ( + SuperconductorModel(int(mfile.get("i_tf_sc_mat", scan=scan))).full_name + if i_tf_sc_mat > 0 + else "Resistive Copper" + ) vssoft = mfile.get("vs_plasma_res_ramp", scan=scan) + mfile.get( "vs_plasma_ind_ramp", scan=scan @@ -8766,17 +8283,8 @@ def plot_power_info(axis: plt.Axes, mfile: MFile, scan: int): scan : scan number to use """ - xmin = 0 - xmax = 1 - ymin = -16 - ymax = 1 - axis.text(-0.05, 1, "Power flows:", ha="left", va="center") - axis.set_ylim([ymin, ymax]) - axis.set_xlim([xmin, xmax]) - axis.set_axis_off() - axis.set_autoscaley_on(False) - axis.set_autoscalex_on(False) + _setup_axis(axis, xmin=0, xmax=1, ymin=-16, ymax=1) gross_eff = 100.0 * ( mfile.get("p_plant_electric_gross_mw", scan=scan) @@ -8857,28 +8365,17 @@ def plot_current_drive_info(axis: plt.Axes, mfile: MFile, scan: int): scan : scan number to use """ - xmin = 0 - xmax = 1 - ymin = -16 - ymax = 1 + _setup_axis(axis, xmin=0, xmax=1, ymin=-16, ymax=1) + i_hcd_primary = mfile.get("i_hcd_primary", scan=scan) - nbi = False - ecrh = False - ebw = False - lhcd = False - iccd = False - - if i_hcd_primary in {5, 8}: - nbi = True + + if nbi := (i_hcd_primary in {5, 8}): axis.text(-0.05, 1, "Neutral Beam Current Drive:", ha="left", va="center") - if i_hcd_primary in {3, 7, 10, 11, 13}: - ecrh = True + if ecrh := (i_hcd_primary in {3, 7, 10, 11, 13}): axis.text(-0.05, 1, "Electron Cyclotron Current Drive:", ha="left", va="center") - if i_hcd_primary == 12: - ebw = True + if ebw := (i_hcd_primary == 12): axis.text(-0.05, 1, "Electron Bernstein Wave Drive:", ha="left", va="center") - if i_hcd_primary in {1, 4, 6}: - lhcd = True + if lhcd := (i_hcd_primary in {1, 4, 6}): axis.text( -0.05, 1, @@ -8886,30 +8383,22 @@ def plot_current_drive_info(axis: plt.Axes, mfile: MFile, scan: int): ha="left", va="center", ) - if i_hcd_primary == 2: - iccd = True + if iccd := (i_hcd_primary == 2): axis.text(-0.05, 1, "Ion Cyclotron Current Drive:", ha="left", va="center") - if "i_hcd_secondary" in mfile.data: + i_hcd_secondary = mfile.get("i_hcd_secondary", scan=scan) or 0 + if i_hcd_secondary in {5, 8}: + secondary_heating = "NBI" + elif i_hcd_secondary in {3, 7, 10, 11, 13}: + secondary_heating = "ECH" + elif i_hcd_secondary == 12: + secondary_heating = "EBW" + elif i_hcd_secondary in {1, 4, 6}: + secondary_heating = "LHCD" + elif i_hcd_secondary == 2: + secondary_heating = "ICCD" + else: secondary_heating = "" - i_hcd_secondary = mfile.get("i_hcd_secondary", scan=scan) - - if i_hcd_secondary in {5, 8}: - secondary_heating = "NBI" - if i_hcd_secondary in {3, 7, 10, 11, 13}: - secondary_heating = "ECH" - if i_hcd_secondary == 12: - secondary_heating = "EBW" - if i_hcd_secondary in {1, 4, 6}: - secondary_heating = "LHCD" - if i_hcd_secondary == 2: - secondary_heating = "ICCD" - - axis.set_ylim([ymin, ymax]) - axis.set_xlim([xmin, xmax]) - axis.set_axis_off() - axis.set_autoscaley_on(False) - axis.set_autoscalex_on(False) pinjie = mfile.get("p_hcd_injected_total_mw", scan=scan) p_plasma_separatrix_mw = mfile.get("p_plasma_separatrix_mw", scan=scan) @@ -8929,166 +8418,88 @@ def plot_current_drive_info(axis: plt.Axes, mfile: MFile, scan: int): # Assume Martin scaling if pthresh is not printed # Accounts for pthresh not being written prior to issue #679 and #680 - if "p_l_h_threshold_mw" in mfile.data: - pthresh = mfile.get("p_l_h_threshold_mw", scan=scan) - else: - pthresh = mfile.get("l_h_threshold_powers(6)", scan=scan) + pthresh_name = ( + "p_l_h_threshold_mw" + if "p_l_h_threshold_mw" in mfile.data + else "l_h_threshold_powers(6)" + ) + pthresh = mfile.get(pthresh_name, scan=scan) flh = p_plasma_separatrix_mw / pthresh hstar = mfile.get("hstar", scan=scan) - if ecrh: - data = [ - (pinjie, "Steady state auxiliary power", "MW"), - ("p_hcd_primary_extra_heat_mw", "Power for heating only", "MW"), - ("f_c_plasma_bootstrap", "Bootstrap fraction", ""), - ("f_c_plasma_auxiliary", "Auxiliary fraction", ""), - ("f_c_plasma_inductive", "Inductive fraction", ""), - ("p_plasma_loss_mw", "Plasma heating used for H factor", "MW"), - ( - "eta_cd_hcd_primary", - "Current drive efficiency", - "A W$^{-1}$", - ), - (pdivr, r"$\frac{P_{\mathrm{div}}}{R_{0}}$", "MW m$^{-1}$"), - ( - pdivnr, - r"$\frac{P_{\mathrm{div}}}{\langle n \rangle R_{0}}$", - r"$\times 10^{-20}$ MW m$^{2}$", - ), - (flh, r"$\frac{P_{\mathrm{div}}}{P_{\mathrm{LH}}}$", ""), - (hstar, "H* (non-rad. corr.)", ""), - ] - # i_hcd_secondary is now always in the MFILE with = 0 meaning no fixed heating - if mfile.get("i_hcd_secondary", scan=scan) != 0: - data.insert( - 1, ("pinjmwfix", f"{secondary_heating} secondary auxiliary power", "MW") - ) - data[0] = ((pinjie - pinjmwfix), "Primary auxiliary power", "MW") - data.insert(2, (pinjie, "Total auxillary power", "MW")) - - if nbi: - data = [ - (pinjie, "Steady state auxiliary power", "MW"), - ("p_hcd_primary_extra_heat_mw", "Power for heating only", "MW"), - ("f_c_plasma_bootstrap", "Bootstrap fraction", ""), - ("f_c_plasma_auxiliary", "Auxiliary fraction", ""), - ("f_c_plasma_inductive", "Inductive fraction", ""), + data = [ + (pinjie, "Steady state auxiliary power", "MW"), + ("p_hcd_primary_extra_heat_mw", "Power for heating only", "MW"), + ("f_c_plasma_bootstrap", "Bootstrap fraction", ""), + ("f_c_plasma_auxiliary", "Auxiliary fraction", ""), + ("f_c_plasma_inductive", "Inductive fraction", ""), + ("p_plasma_loss_mw", "Plasma heating used for H factor", "MW"), + (pdivr, r"$\frac{P_{\mathrm{div}}}{R_{0}}$", "MW m$^{-1}$"), + ( + pdivnr, + r"$\frac{P_{\mathrm{div}}}{\langle n \rangle R_{0}}$", + r"$\times 10^{-20}$ MW m$^{2}$", + ), + (flh, r"$\frac{P_{\mathrm{div}}}{P_{\mathrm{LH}}}$", ""), + (hstar, "H* (non-rad. corr.)", ""), + ] + # Optional override based on condition + field_overrides = { + "ecrh": ( + "eta_cd_hcd_primary", + r"$\frac{P_{\mathrm{div}}}{R_{0}}$", + "A W$^{-1}$", + ), + "nbi": ( ("gamnb", "NB gamma", "$10^{20}$ A W$^{-1}$ m$^{-2}$"), ("e_beam_kev", "NB energy", "keV"), - ("p_plasma_loss_mw", "Plasma heating used for H factor", "MW"), - (pdivr, r"$\frac{P_{\mathrm{div}}}{R_{0}}$", "MW m$^{-1}$"), - ( - pdivnr, - r"$\frac{P_{\mathrm{div}}}{\langle n \rangle R_{0}}$", - r"$\times 10^{-20}$ MW m$^{2}$", - ), - (flh, r"$\frac{P_{\mathrm{div}}}{P_{\mathrm{LH}}}$", ""), - (hstar, "H* (non-rad. corr.)", ""), - ] - if mfile.get("i_hcd_secondary", scan=scan) != 0: - data.insert( - 1, ("pinjmwfix", f"{secondary_heating} secondary auxiliary power", "MW") - ) - data[0] = ((pinjie - pinjmwfix), "Primary auxiliary power", "MW") - data.insert(2, (pinjie, "Total auxillary power", "MW")) - - if ebw: - data = [ - (pinjie, "Steady state auxiliary power", "MW"), - ("p_hcd_primary_extra_heat_mw", "Power for heating only", "MW"), - ("f_c_plasma_bootstrap", "Bootstrap fraction", ""), - ("f_c_plasma_auxiliary", "Auxiliary fraction", ""), - ("f_c_plasma_inductive", "Inductive fraction", ""), - ("p_plasma_loss_mw", "Plasma heating used for H factor", "MW"), - ( - "eta_cd_norm_hcd_primary", - "Normalised current drive efficiency of primary HCD system", - "(10$^{20}$ A/(Wm$^{2}$))", - ), - (pdivr, r"$\frac{P_{\mathrm{div}}}{R_{0}}$", "MW m$^{-1}$"), - ( - pdivnr, - r"$\frac{P_{\mathrm{div}}}{\langle n \rangle R_{0}}$", - r"$\times 10^{-20}$ MW m$^{2}$", - ), - (flh, r"$\frac{P_{\mathrm{div}}}{P_{\mathrm{LH}}}$", ""), - (hstar, "H* (non-rad. corr.)", ""), - ] - if "i_hcd_secondary" in mfile.data: - data.insert( - 1, ("pinjmwfix", f"{secondary_heating} secondary auxiliary power", "MW") - ) - data[0] = ((pinjie - pinjmwfix), "Primary auxiliary power", "MW") - data.insert(2, (pinjie, "Total auxillary power", "MW")) - - if lhcd: - data = [ - (pinjie, "Steady state auxiliary power", "MW"), - ("p_hcd_primary_extra_heat_mw", "Power for heating only", "MW"), - ("f_c_plasma_bootstrap", "Bootstrap fraction", ""), - ("f_c_plasma_auxiliary", "Auxiliary fraction", ""), - ("f_c_plasma_inductive", "Inductive fraction", ""), - ("p_plasma_loss_mw", "Plasma heating used for H factor", "MW"), - ( - "eta_cd_norm_hcd_primary", - "Normalised current drive efficiency", - "(10$^{20}$ A/(Wm$^{2}$))", - ), - (pdivr, r"$\frac{P_{\mathrm{div}}}{R_{0}}$", "MW m$^{-1}$"), - ( - pdivnr, - r"$\frac{P_{\mathrm{div}}}{\langle n \rangle R_{0}}$", - r"$\times 10^{-20}$ MW m$^{2}$", - ), - (flh, r"$\frac{P_{\mathrm{div}}}{P_{\mathrm{LH}}}$", ""), - (hstar, "H* (non-rad. corr.)", ""), - ] - if "i_hcd_secondary" in mfile.data: - data.insert( - 1, ("pinjmwfix", f"{secondary_heating} secondary auxiliary power", "MW") - ) - data[0] = ((pinjie - pinjmwfix), "Primary auxiliary power", "MW") - data.insert(2, (pinjie, "Total auxillary power", "MW")) + ), + "ebw": ( + "eta_cd_norm_hcd_primary", + "Normalised current drive efficiency of primary HCD system", + "(10$^{20}$ A/(Wm$^{2}$))", + ), + "lhcd": ( + "eta_cd_norm_hcd_primary", + "Normalised current drive efficiency", + "(10$^{20}$ A/(Wm$^{2}$))", + ), + "iccd": ( + "eta_cd_norm_hcd_primary", + "Normalised current drive efficiency", + "(10$^{20}$ A/(Wm$^{2}$))", + ), + } - if iccd: - data = [ - (pinjie, "Steady state auxiliary power", "MW"), - ("p_hcd_primary_extra_heat_mw", "Power for heating only", "MW"), - ("f_c_plasma_bootstrap", "Bootstrap fraction", ""), - ("f_c_plasma_auxiliary", "Auxiliary fraction", ""), - ("f_c_plasma_inductive", "Inductive fraction", ""), - ("p_plasma_loss_mw", "Plasma heating used for H factor", "MW"), - ( - "eta_cd_norm_hcd_primary", - "Normalised current drive efficiency", - "(10$^{20}$ A/(Wm$^{2}$))", - ), - (pdivr, r"$\frac{P_{\mathrm{div}}}{R_{0}}$", "MW m$^{-1}$"), - ( - pdivnr, - r"$\frac{P_{\mathrm{div}}}{\langle n \rangle R_{0}}$", - r"$\times 10^{-20}$ MW m$^{2}$", - ), - (flh, r"$\frac{P_{\mathrm{div}}}{P_{\mathrm{LH}}}$", ""), - (hstar, "H* (non-rad. corr.)", ""), - ] - if "i_hcd_secondary" in mfile.data: - data.insert( - 1, ("pinjmwfix", f"{secondary_heating} secondary auxiliary power", "MW") - ) - data[0] = ((pinjie - pinjmwfix), "Primary auxiliary power", "MW") - data.insert(2, (pinjie, "Total auxillary power", "MW")) + if ecrh: + data.insert(6, field_overrides["ecrh"]) + elif nbi: + data.insert(6, field_overrides["nbi"][0]) + data.insert(7, field_overrides["nbi"][1]) + elif ebw: + data.insert(6, field_overrides["ebw"]) + elif lhcd: + data.insert(6, field_overrides["lhcd"]) + elif iccd: + data.insert(6, field_overrides["iccd"]) + + # Secondary heating logic — common across all cases + if mfile.get("i_hcd_secondary", scan=scan) != 0: + data.insert( + 1, ("pinjmwfix", f"{secondary_heating} secondary auxiliary power", "MW") + ) + data[0] = ((pinjie - pinjmwfix), "Primary auxiliary power", "MW") + data.insert(2, (pinjie, "Total auxillary power", "MW")) coe = mfile.get("coe", scan=scan) - if coe == 0.0: # noqa: RUF069 - data.append(("", "", "")) - data.append(("#Costs", "", "")) - data.append(("", "Cost output not selected", "")) - else: - data.append(("", "", "")) - data.append(("#Costs", "", "")) - data.append((coe, "Cost of electricity", r"\$/MWh")) + data.extend(( + ("", "", ""), + ("#Costs", "", ""), + ("", "Cost output not selected", "") + if coe == 0.0 # noqa: RUF069 + else (coe, "Cost of electricity", r"\$/MWh"), + )) plot_info(axis, data, mfile, scan) @@ -9105,43 +8516,28 @@ def plot_bootstrap_comparison(axis: plt.Axes, mfile: MFile, scan: int): scan : scan number to use """ - boot_ipdg = mfile.get("f_c_plasma_bootstrap_iter89", scan=scan) - boot_sauter = mfile.get("f_c_plasma_bootstrap_sauter", scan=scan) - boot_nenins = mfile.get("f_c_plasma_bootstrap_nevins", scan=scan) - boot_wilson = mfile.get("f_c_plasma_bootstrap_wilson", scan=scan) - boot_sakai = mfile.get("f_c_plasma_bootstrap_sakai", scan=scan) - boot_aries = mfile.get("f_c_plasma_bootstrap_aries", scan=scan) - boot_andrade = mfile.get("f_c_plasma_bootstrap_andrade", scan=scan) - boot_hoang = mfile.get("f_c_plasma_bootstrap_hoang", scan=scan) - boot_wong = mfile.get("f_c_plasma_bootstrap_wong", scan=scan) - boot_gi_I = mfile.get("bscf_gi_i", scan=scan) # noqa: N806 - boot_gi_II = mfile.get("bscf_gi_ii", scan=scan) # noqa: N806 - boot_sugiyama_l = mfile.get("f_c_plasma_bootstrap_sugiyama_l", scan=scan) - boot_sugiyama_h = mfile.get("f_c_plasma_bootstrap_sugiyama_h", scan=scan) - # Data for the box plot data = { - "IPDG": boot_ipdg, - "Sauter": boot_sauter, - "Nevins": boot_nenins, - "Wilson": boot_wilson, - "Sakai": boot_sakai, - "ARIES": boot_aries, - "Andrade": boot_andrade, - "Hoang": boot_hoang, - "Wong": boot_wong, - "Gi-I": boot_gi_I, - "Gi-II": boot_gi_II, - "Sugiyama (L-mode)": boot_sugiyama_l, - "Sugiyama (H-mode)": boot_sugiyama_h, + "IPDG": mfile.get("f_c_plasma_bootstrap_iter89", scan=scan), + "Sauter": mfile.get("f_c_plasma_bootstrap_sauter", scan=scan), + "Nevins": mfile.get("f_c_plasma_bootstrap_nevins", scan=scan), + "Wilson": mfile.get("f_c_plasma_bootstrap_wilson", scan=scan), + "Sakai": mfile.get("f_c_plasma_bootstrap_sakai", scan=scan), + "ARIES": mfile.get("f_c_plasma_bootstrap_aries", scan=scan), + "Andrade": mfile.get("f_c_plasma_bootstrap_andrade", scan=scan), + "Hoang": mfile.get("f_c_plasma_bootstrap_hoang", scan=scan), + "Wong": mfile.get("f_c_plasma_bootstrap_wong", scan=scan), + "Gi-I": mfile.get("bscf_gi_i", scan=scan), + "Gi-II": mfile.get("bscf_gi_ii", scan=scan), + "Sugiyama (L-mode)": mfile.get("f_c_plasma_bootstrap_sugiyama_l", scan=scan), + "Sugiyama (H-mode)": mfile.get("f_c_plasma_bootstrap_sugiyama_h", scan=scan), } # Create the violin plot - axis.violinplot(data.values(), showextrema=False) + data_values = list(data.values()) + axis.violinplot(data_values, showextrema=False) # Create the box plot - axis.boxplot( - data.values(), showfliers=True, showmeans=True, meanline=True, widths=0.3 - ) + axis.boxplot(data_values, showfliers=True, showmeans=True, meanline=True, widths=0.3) # Scatter plot for each data point colors = plt.cm.plasma(np.linspace(0, 1, len(data.values()))) @@ -9150,7 +8546,6 @@ def plot_bootstrap_comparison(axis: plt.Axes, mfile: MFile, scan: int): axis.legend(loc="upper left", bbox_to_anchor=(1, 1)) # Calculate average, standard deviation, and median - data_values = list(data.values()) avg_bootstrap = np.mean(data_values) std_bootstrap = np.std(data_values) median_bootstrap = np.median(data_values) @@ -9176,7 +8571,7 @@ def plot_bootstrap_comparison(axis: plt.Axes, mfile: MFile, scan: int): axis.set_title("Bootstrap Current Fraction ($f_\\text{BS}$) Comparison") axis.set_ylabel("Bootstrap Current Fraction") - axis.set_xlim([0.5, 1.5]) + axis.set_xlim(0.5, 1.5) axis.set_xticks([]) axis.set_xticklabels([]) axis.set_facecolor("#f0f0f0") @@ -9217,22 +8612,21 @@ def plot_sol_power_decay_length_comparison(axis: plt.Axes, mfile: MFile, scan: i f"{OutbordSOLPowerDecayLengthModel.EICH_2011_JET.description}": len_plasma_sol_eich11_jet_power_decay_mm, f"{OutbordSOLPowerDecayLengthModel.EICH_2011_JET_ASDEX.description}": len_plasma_sol_eich11_jet_asdex_power_decay_mm, } + data_values = list(data.values()) + # Create the violin plot - axis.violinplot(data.values(), showextrema=False) + axis.violinplot(data_values, showextrema=False) # Create the box plot - axis.boxplot( - data.values(), showfliers=True, showmeans=True, meanline=True, widths=0.3 - ) + axis.boxplot(data_values, showfliers=True, showmeans=True, meanline=True, widths=0.3) # Scatter plot for each data point - colors = plt.cm.plasma(np.linspace(0, 1, len(data.values()))) + colors = plt.cm.plasma(np.linspace(0, 1, len(data_values))) for index, (key, value) in enumerate(data.items()): axis.scatter(1, value, color=colors[index], label=key, alpha=1.0) axis.legend(loc="upper left", bbox_to_anchor=(1, 1)) # Calculate average, standard deviation, and median - data_values = list(data.values()) avg_decay_length = np.mean(data_values) std_decay_length = np.std(data_values) median_decay_length = np.median(data_values) @@ -9262,7 +8656,7 @@ def plot_sol_power_decay_length_comparison(axis: plt.Axes, mfile: MFile, scan: i axis.set_title("SOL Power Decay Length ($\\lambda_q$) Comparison") axis.set_ylabel("Power Decay Length [mm]") - axis.set_xlim([0.5, 1.5]) + axis.set_xlim(0.5, 1.5) axis.set_xticks([]) axis.set_xticklabels([]) axis.set_facecolor("#f0f0f0") @@ -9323,11 +8717,11 @@ def plot_brunner_divertor_power_split_comparison_stackplot( alpha=0.5, label="$\u0394 r_{\\mathrm{sep}}$", ) - axis.set_ylim([0.0, 1.0]) - axis.set_xlim([ + axis.set_ylim(0.0, 1.0) + axis.set_xlim( -5 * len_plasma_sol_outboard_pd, 5 * len_plasma_sol_outboard_pd, - ]) + ) axis.grid(True, which="both", linestyle="--", linewidth=0.5, alpha=0.35) axis.set_title("Brunner Divertor Power Split Fractions") axis.set_xlabel("$\\Delta r_{\\mathrm{sep}}$ [m]") @@ -9539,65 +8933,49 @@ def plot_h_threshold_comparison(axis: plt.Axes, mfile: MFile, scan: int, u_seed= u_seed : (Default value = None) """ - iter_nominal = mfile.get("l_h_threshold_powers(1)", scan=scan) - iter_upper = mfile.get("l_h_threshold_powers(2)", scan=scan) - iter_lower = mfile.get("l_h_threshold_powers(3)", scan=scan) - iter_1997_1 = mfile.get("l_h_threshold_powers(4)", scan=scan) - iter_1997_2 = mfile.get("l_h_threshold_powers(5)", scan=scan) - martin_nominal = mfile.get("l_h_threshold_powers(6)", scan=scan) - martin_upper = mfile.get("l_h_threshold_powers(7)", scan=scan) - martin_lower = mfile.get("l_h_threshold_powers(8)", scan=scan) - snipes_nominal = mfile.get("l_h_threshold_powers(9)", scan=scan) - snipes_upper = mfile.get("l_h_threshold_powers(10)", scan=scan) - snipes_lower = mfile.get("l_h_threshold_powers(11)", scan=scan) - snipes_closed_nominal = mfile.get("l_h_threshold_powers(12)", scan=scan) - snipes_closed_upper = mfile.get("l_h_threshold_powers(13)", scan=scan) - snipes_closed_lower = mfile.get("l_h_threshold_powers(14)", scan=scan) - hubbard_nominal = mfile.get("l_h_threshold_powers(15)", scan=scan) - hubbard_lower = mfile.get("l_h_threshold_powers(16)", scan=scan) - hubbard_upper = mfile.get("l_h_threshold_powers(17)", scan=scan) - hubbard_2017 = mfile.get("l_h_threshold_powers(18)", scan=scan) - martin_aspect_nominal = mfile.get("l_h_threshold_powers(19)", scan=scan) - martin_aspect_upper = mfile.get("l_h_threshold_powers(20)", scan=scan) - martin_aspect_lower = mfile.get("l_h_threshold_powers(21)", scan=scan) - # Data for the box plot data = { - "ITER 1996 Nominal": iter_nominal, - "ITER 1996 Upper": iter_upper, - "ITER 1996 Lower": iter_lower, - "ITER 1997 (1)": iter_1997_1, - "ITER 1997 (2)": iter_1997_2, - "Martin Nominal": martin_nominal, - "Martin Upper": martin_upper, - "Martin Lower": martin_lower, - "Snipes Nominal": snipes_nominal, - "Snipes Upper": snipes_upper, - "Snipes Lower": snipes_lower, - "Snipes Closed Divertor Nominal": snipes_closed_nominal, - "Snipes Closed Divertor Upper": snipes_closed_upper, - "Snipes Closed Divertor Lower": snipes_closed_lower, - "Hubbard Nominal (I-mode)": hubbard_nominal, - "Hubbard Lower (I-mode)": hubbard_lower, - "Hubbard Upper (I-mode)": hubbard_upper, - "Hubbard 2017 (I-mode)": hubbard_2017, - "Martin Aspect Corrected Nominal": martin_aspect_nominal, - "Martin Aspect Corrected Upper": martin_aspect_upper, - "Martin Aspect Corrected Lower": martin_aspect_lower, + "ITER 1996 Nominal": mfile.get("l_h_threshold_powers(1)", scan=scan), + "ITER 1996 Upper": mfile.get("l_h_threshold_powers(2)", scan=scan), + "ITER 1996 Lower": mfile.get("l_h_threshold_powers(3)", scan=scan), + "ITER 1997 (1)": mfile.get("l_h_threshold_powers(4)", scan=scan), + "ITER 1997 (2)": mfile.get("l_h_threshold_powers(5)", scan=scan), + "Martin Nominal": mfile.get("l_h_threshold_powers(6)", scan=scan), + "Martin Upper": mfile.get("l_h_threshold_powers(7)", scan=scan), + "Martin Lower": mfile.get("l_h_threshold_powers(8)", scan=scan), + "Snipes Nominal": mfile.get("l_h_threshold_powers(9)", scan=scan), + "Snipes Upper": mfile.get("l_h_threshold_powers(10)", scan=scan), + "Snipes Lower": mfile.get("l_h_threshold_powers(11)", scan=scan), + "Snipes Closed Divertor Nominal": mfile.get( + "l_h_threshold_powers(12)", scan=scan + ), + "Snipes Closed Divertor Upper": mfile.get("l_h_threshold_powers(13)", scan=scan), + "Snipes Closed Divertor Lower": mfile.get("l_h_threshold_powers(14)", scan=scan), + "Hubbard Nominal (I-mode)": mfile.get("l_h_threshold_powers(15)", scan=scan), + "Hubbard Lower (I-mode)": mfile.get("l_h_threshold_powers(16)", scan=scan), + "Hubbard Upper (I-mode)": mfile.get("l_h_threshold_powers(17)", scan=scan), + "Hubbard 2017 (I-mode)": mfile.get("l_h_threshold_powers(18)", scan=scan), + "Martin Aspect Corrected Nominal": mfile.get( + "l_h_threshold_powers(19)", scan=scan + ), + "Martin Aspect Corrected Upper": mfile.get( + "l_h_threshold_powers(20)", scan=scan + ), + "Martin Aspect Corrected Lower": mfile.get( + "l_h_threshold_powers(21)", scan=scan + ), } - + data_values = list(data.values()) # Create the violin plot - axis.violinplot(data.values(), showextrema=False) + axis.violinplot(data_values, showextrema=False) # Create the box plot - axis.boxplot( - data.values(), showfliers=True, showmeans=True, meanline=True, widths=0.3 - ) + axis.boxplot(data_values, showfliers=True, showmeans=True, meanline=True, widths=0.3) # Scatter plot for each data point - colors = plt.cm.plasma(np.linspace(0, 1, len(data.values()))) + colors = plt.cm.plasma(np.linspace(0, 1, len(data_values))) generator = np.random.default_rng(seed=u_seed) - x_values = generator.normal(loc=1, scale=0.01, size=len(data.values())) + x_values = generator.normal(loc=1, scale=0.01, size=len(data_values)) for index, (key, value) in enumerate(data.items()): if "ITER 1996" in key: color = "blue" @@ -9621,7 +8999,6 @@ def plot_h_threshold_comparison(axis: plt.Axes, mfile: MFile, scan: int, u_seed= axis.legend(loc="upper left", bbox_to_anchor=(-1.1, 1), ncol=2) # Calculate average, standard deviation, and median - data_values = list(data.values()) avg_threshold = np.mean(data_values) std_threshold = np.std(data_values) median_threshold = np.median(data_values) @@ -9651,7 +9028,7 @@ def plot_h_threshold_comparison(axis: plt.Axes, mfile: MFile, scan: int, u_seed= axis.set_title("L-H Threshold ($P_\\text{LH}$) Comparison") axis.set_ylabel("L-H threshold power [MW]") - axis.set_xlim([0.5, 1.5]) + axis.set_xlim(0.5, 1.5) axis.set_xticks([]) axis.set_xticklabels([]) @@ -9904,14 +9281,13 @@ def plot_confinement_time_comparison( rf"{ConfinementTimeModel.ITPA20.full_name}": itpa20, rf"{ConfinementTimeModel.ITPA20_IL.full_name}": itpa20_ilc, } + data_values = list(data.values()) # Create the violin plot - axis.violinplot(data.values(), showextrema=False) + axis.violinplot(data_values, showextrema=False) # Create the box plot - axis.boxplot( - data.values(), showfliers=True, showmeans=True, meanline=True, widths=0.3 - ) + axis.boxplot(data_values, showfliers=True, showmeans=True, meanline=True, widths=0.3) # Scatter plot for each data point # Use a set of distinct colors for better differentiation @@ -9956,7 +9332,6 @@ def plot_confinement_time_comparison( axis.legend(loc="upper left", bbox_to_anchor=(-1.3, 0.75), ncol=2) # Calculate average, standard deviation, and median - data_values = list(data.values()) avg_threshold = np.mean(data_values) std_threshold = np.std(data_values) median_threshold = np.median(data_values) @@ -9993,7 +9368,7 @@ def plot_confinement_time_comparison( axis.set_title("Confinement time ($\\tau_{\\text{E}}$) Comparison") axis.set_ylabel("Confinement time, $\\tau_{\\text{E}}$ [s]") - axis.set_xlim([0.5, 1.5]) + axis.set_xlim(0.5, 1.5) axis.set_xticks([]) axis.set_xticklabels([]) @@ -10481,14 +9856,13 @@ def plot_density_limit_comparison(axis: plt.Axes, mfile: MFile, scan: int): "Greenwald": greenwald, "ASDEX New": asdex_new, } + data_values = list(data.values()) # Create the violin plot - axis.violinplot(data.values(), showextrema=False) + axis.violinplot(data_values, showextrema=False) # Create the box plot - axis.boxplot( - data.values(), showfliers=True, showmeans=True, meanline=True, widths=0.3 - ) + axis.boxplot(data_values, showfliers=True, showmeans=True, meanline=True, widths=0.3) # Scatter plot for each data point colors = plt.cm.plasma(np.linspace(0, 1, len(data.values()))) @@ -10497,7 +9871,6 @@ def plot_density_limit_comparison(axis: plt.Axes, mfile: MFile, scan: int): axis.legend(loc="upper left", bbox_to_anchor=(1, 1)) # Calculate average, standard deviation, and median - data_values = list(data.values()) avg_density_limit = np.mean(data_values) std_density_limit = np.std(data_values) median_density_limit = np.median(data_values) @@ -10529,7 +9902,7 @@ def plot_density_limit_comparison(axis: plt.Axes, mfile: MFile, scan: int): axis.set_title("Density Limit Comparison") axis.set_ylabel(r"Density Limit [$10^{20}$ m$^{-3}$]") axis.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{x * 1e-20:.1f}")) - axis.set_xlim([0.5, 1.5]) + axis.set_xlim(0.5, 1.5) axis.set_xticks([]) axis.set_xticklabels([]) axis.set_facecolor("#f0f0f0") @@ -10687,39 +10060,25 @@ def plot_cs_coil_structure( ) textstr_cs = ( - f"$\\mathbf{{Coil \\ parameters:}}$\n \n" + f"$\\mathbf{{Coil \\ parameters:}}$\n\n" f"CS height vs TF internal height: {mfile.get('f_z_cs_tf_internal', scan=scan):.2f}\n" f"CS thickness: {mfile.get('dr_cs', scan=scan):.4f} m\n" f"CS radial middle: {mfile.get('r_cs_middle', scan=scan):.4f} m\n" f"CS full height: {mfile.get('dz_cs_full', scan=scan):.4f} m\n" f"CS full width: {mfile.get('dr_cs_full', scan=scan):.4f} m\n" f"CS poloidal area: {mfile.get('a_cs_poloidal', scan=scan):.4f} m$^2$\n" - f"CS top-down toroidal area: {mfile.get('a_cs_toroidal', scan=scan):.4f} m$^2$\n" - f"$N_{{\\text{{turns}}}}:$ {mfile.get('n_pf_coil_turns[n_cs_pf_coils-1]', scan=scan):,.2f}\n" - f"$I_{{\\text{{peak}}}}:$ {mfile.get('c_pf_cs_coils_peak_ma[n_cs_pf_coils-1]', scan=scan):.3f} MA\n" - f"$B_{{\\text{{peak}}}}:$ {mfile.get('b_pf_coil_peak[n_cs_pf_coils-1]', scan=scan):.3f} T\n" - f"$F_{{\\text{{z,self,peak}}}}:$ {mfile.get('forc_z_cs_self_peak_midplane', scan=scan) / 1e6:.3f} MN\n" - f"$\\sigma_{{\\text{{z,self,peak}}}}:$ {mfile.get('stress_z_cs_self_peak_midplane', scan=scan) / 1e6:.3f} MPa\n" - f"$\\sigma_{{\\text{{mises,peak}}}}:$ {mfile.get('stress_mises_cs_peak', scan=scan) / 1e6:.3f} MPa\n" - f"$\\tau_{{\\text{{shear,peak}}}}:$ {mfile.get('stress_shear_cs_peak', scan=scan) / 1e6:.3f} MPa " - ) - - axis.text( - 0.5, - 0.6, - textstr_cs, - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightyellow", - "alpha": 1.0, - "linewidth": 2, - }, + f"CS top-down toroidal area: {mfile.get('a_cs_toroidal', scan=scan):.4f} m$^2$\n" + f"$N_{{\\text{{turns}}}}:$ {mfile.get('n_pf_coil_turns[n_cs_pf_coils-1]', scan=scan):,.2f}\n" + f"$I_{{\\text{{peak}}}}:$ {mfile.get('c_pf_cs_coils_peak_ma[n_cs_pf_coils-1]', scan=scan):.3f} MA\n" + f"$B_{{\\text{{peak}}}}:$ {mfile.get('b_pf_coil_peak[n_cs_pf_coils-1]', scan=scan):.3f} T\n" + f"$F_{{\\text{{z,self,peak}}}}:$ {mfile.get('forc_z_cs_self_peak_midplane', scan=scan) / 1e6:.3f} MN\n" + f"$\\sigma_{{\\text{{z,self,peak}}}}:$ {mfile.get('stress_z_cs_self_peak_midplane', scan=scan) / 1e6:.3f} MPa\n" + f"$\\sigma_{{\\text{{mises,peak}}}}:$ {mfile.get('stress_mises_cs_peak', scan=scan) / 1e6:.3f} MPa\n" + f"$\\tau_{{\\text{{shear,peak}}}}:$ {mfile.get('stress_shear_cs_peak', scan=scan) / 1e6:.3f} MPa " ) + axis.text(0.5, 0.6, textstr_cs, **_text_layout(fig), bbox=_box_style("lightyellow")) + # Plot the current filament points as blue dots and label them axis.plot( @@ -10862,16 +10221,8 @@ def plot_cs_turn_structure(axis: plt.Axes, fig, mfile: MFile, scan: int): 0.7, 0.375, textstr_turn, - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightyellow", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("lightyellow"), ) axis.set_xlim(-dr_cs_turn * 0.2, dr_cs_turn * 1.2) @@ -11936,6 +11287,68 @@ def plot_tf_stress(axis: plt.Axes, mfile: MFile): plt.tight_layout() +def draw_bend( + ax: Axes, + elbow_radius: float, + theta_span: float, + radius_pipe: float, + title: str = "Bend", + alpha: float = 0.8, +): + """ + Draws a circular pipe bend with centerline and inner/outer boundaries. + + Parameters + ---------- + ax: + Target axes for plotting. + elbow_radius: + Radius of the elbow [m]. + theta_span: + Array of angles [0, θ] where θ is pi/2 or pi [rad] + radius_pipe: + Pipe radius (fallback to 0.1m if not provided) [m] + title: + Plot title string. + alpha: + fill opacity + """ + # Convert all inputs to mm + elbow_radius_mm = elbow_radius * 1000 + pipe_radius_mm = radius_pipe * 1000 + + theta = np.linspace(0, theta_span, 100) + x_center = elbow_radius_mm * np.cos(theta) + y_center = elbow_radius_mm * np.sin(theta) + + # Outer and inner walls (offset by ± pipe radius in mm) + x_outer = (elbow_radius_mm + pipe_radius_mm) * np.cos(theta) + y_outer = (elbow_radius_mm + pipe_radius_mm) * np.sin(theta) + x_inner = (elbow_radius_mm - pipe_radius_mm) * np.cos(theta) + y_inner = (elbow_radius_mm - pipe_radius_mm) * np.sin(theta) + + # Plot + ax.plot(x_center, y_center, color="black", linestyle="--", label="Centerline") + ax.plot(x_outer, y_outer, color="black") + ax.plot(x_inner, y_inner, color="black") + ax.fill( + np.concatenate([x_outer, x_inner[::-1]]), + np.concatenate([y_outer, y_inner[::-1]]), + color="lightgrey", + alpha=alpha, + ) + + ax.set_aspect("equal") + ax.set_xlabel("X [mm]") + ax.set_ylabel("Y [mm]") + ax.set_title(title) + ax.grid(True, linestyle="--", alpha=0.3) + + # Legend: Centerline + pipe radius info + legend_text = f"Centerline\nPipe radius: {pipe_radius_mm:.2f} mm\nElbow radius: {elbow_radius_mm:.2f} mm" + ax.legend([legend_text], loc="upper right") + + def plot_blkt_pipe_bends(fig, m_file, scan: int): """Plot the blanket pipe bends on the given axis, with axes in mm. @@ -11951,80 +11364,18 @@ def plot_blkt_pipe_bends(fig, m_file, scan: int): ax_90 = fig.add_subplot(341) ax_180 = fig.add_subplot(342) - # Get pipe radius from m_file, fallback to 0.1 m r = m_file.get("radius_blkt_channel", scan=scan) - elbow_radius_90 = m_file.get("radius_blkt_channel_90_bend", scan=scan) - - # --- 90 degree bend --- - theta_90 = np.linspace(0, np.pi / 2, 100) - # Convert coordinates from meters to millimeters - x_center_90 = elbow_radius_90 * np.cos(theta_90) * 1000 - y_center_90 = elbow_radius_90 * np.sin(theta_90) * 1000 - x_outer_90 = (elbow_radius_90 + r) * np.cos(theta_90) * 1000 - y_outer_90 = (elbow_radius_90 + r) * np.sin(theta_90) * 1000 - x_inner_90 = (elbow_radius_90 - r) * np.cos(theta_90) * 1000 - y_inner_90 = (elbow_radius_90 - r) * np.sin(theta_90) * 1000 - - ax_90.plot( - x_center_90, y_center_90, color="black", linestyle="--", label="Centerline" - ) - ax_90.plot(x_outer_90, y_outer_90, color="black") - ax_90.plot(x_inner_90, y_inner_90, color="black") - ax_90.fill( - np.concatenate([x_outer_90, x_inner_90[::-1]]), - np.concatenate([y_outer_90, y_inner_90[::-1]]), - color="lightgrey", - alpha=1.0, - ) - ax_90.set_aspect("equal") - ax_90.set_xlabel("X [mm]") - ax_90.set_ylabel("Y [mm]") - ax_90.set_title("Blanket Pipe 90° Bend") - ax_90.grid(True, linestyle="--", alpha=0.3) - # Add legend with radius values - ax_90.legend( - [ - f"Centerline\nPipe radius: {r * 1000:.2f} mm\nElbow radius: {elbow_radius_90 * 1000:.2f} mm" - ], - loc="upper right", - ) - - # --- 180 degree bend --- - - elbow_radius_180 = m_file.get("radius_blkt_channel_180_bend", scan=scan) + fallback_radius = 0.1 # meters - theta_180 = np.linspace(0, np.pi, 100) - x_center_180 = elbow_radius_180 * np.cos(theta_180) * 1000 - y_center_180 = elbow_radius_180 * np.sin(theta_180) * 1000 - x_outer_180 = (elbow_radius_180 + r) * np.cos(theta_180) * 1000 - y_outer_180 = (elbow_radius_180 + r) * np.sin(theta_180) * 1000 - x_inner_180 = (elbow_radius_180 - r) * np.cos(theta_180) * 1000 - y_inner_180 = (elbow_radius_180 - r) * np.sin(theta_180) * 1000 - - ax_180.plot( - x_center_180, y_center_180, color="black", linestyle="--", label="Centerline" + elbow_radius_90 = ( + m_file.get("radius_blkt_channel_90_bend", scan=scan) or fallback_radius ) - ax_180.plot(x_outer_180, y_outer_180, color="black") - ax_180.plot(x_inner_180, y_inner_180, color="black") - ax_180.fill( - np.concatenate([x_outer_180, x_inner_180[::-1]]), - np.concatenate([y_outer_180, y_inner_180[::-1]]), - color="lightgrey", - alpha=1.0, + elbow_radius_180 = ( + m_file.get("radius_blkt_channel_180_bend", scan=scan) or fallback_radius ) - ax_180.set_aspect("equal") - ax_180.set_xlabel("X [mm]") - ax_180.set_ylabel("Y [mm]") - ax_180.set_title("Blanket Pipe 180° Bend") - ax_180.grid(True, linestyle="--", alpha=0.3) - # Add legend with radius values - ax_180.legend( - [ - f"Centerline\nPipe radius: {r * 1000:.2f} mm\nElbow radius: {elbow_radius_180 * 1000:.2f} mm" - ], - loc="upper right", - ) + draw_bend(ax_90, elbow_radius_90, np.pi / 2, r, title="Blanket Pipe 90° Bend") + draw_bend(ax_180, elbow_radius_180, np.pi, r, title="Blanket Pipe 180° Bend") def plot_fw_90_deg_pipe_bend(ax, m_file, scan: int): @@ -12043,45 +11394,13 @@ def plot_fw_90_deg_pipe_bend(ax, m_file, scan: int): r = m_file.get("radius_fw_channel", scan=scan) elbow_radius = m_file.get("radius_fw_channel_90_bend", scan=scan) - # --- 90 degree bend --- - theta_90 = np.linspace(0, np.pi / 2, 100) - # Convert coordinates from meters to millimeters - x_center_90 = elbow_radius * np.cos(theta_90) * 1000 - y_center_90 = elbow_radius * np.sin(theta_90) * 1000 - x_outer_90 = (elbow_radius + r) * np.cos(theta_90) * 1000 - y_outer_90 = (elbow_radius + r) * np.sin(theta_90) * 1000 - x_inner_90 = (elbow_radius - r) * np.cos(theta_90) * 1000 - y_inner_90 = (elbow_radius - r) * np.sin(theta_90) * 1000 - - ax.plot(x_center_90, y_center_90, color="black", linestyle="--", label="Centerline") - ax.plot(x_outer_90, y_outer_90, color="black") - ax.plot(x_inner_90, y_inner_90, color="black") - ax.fill( - np.concatenate([x_outer_90, x_inner_90[::-1]]), - np.concatenate([y_outer_90, y_inner_90[::-1]]), - color="lightgrey", - alpha=1.0, - ) - ax.set_aspect("equal") - ax.set_xlabel("X [mm]") - ax.set_ylabel("Y [mm]") - ax.set_title("First Wall Pipe 90° Bend") - ax.grid(True, linestyle="--", alpha=0.3) - ax.legend( - [ - f"Centerline\nPipe radius: {r * 1000:.2f} mm\nElbow radius: {elbow_radius * 1000:.2f} mm" - ], - loc="upper right", + draw_bend( + ax, elbow_radius, np.pi / 2, r, title="First Wall Pipe 90° Bend", alpha=1.0 ) def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): """Plot the fusion rate density profiles on the given axis""" - fusden_plasma_dt_profile = [] - fusden_plasma_dd_triton_profile = [] - fusden_plasma_dd_helion_profile = [] - fusden_plasma_dhe3_profile = [] - n_plasma_profile_elements = int(mfile.get("n_plasma_profile_elements", scan=scan)) fusden_plasma_dt_profile = [ @@ -12236,9 +11555,9 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): ) axis.set_yscale("log") axis.grid(True, which="both", linestyle="--", alpha=0.5) - axis.set_xlim([0, 1.025]) + axis.set_xlim(0, 1.025) axis.minorticks_on() - axis.set_ylim([1e10, 1e23]) + axis.set_ylim(1e10, 1e23) axis.yaxis.set_major_locator(plt.LogLocator(base=10.0, numticks=10)) axis.yaxis.set_minor_locator( plt.LogLocator(base=10.0, subs=np.arange(1, 10) * 0.1, numticks=100) @@ -12268,23 +11587,15 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): 0.05, 0.85, textstr_general, - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightyellow", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("lightyellow"), ) # ============================================================================ textstr_dt = ( f"Total fusion power: {mfile.get('p_dt_total_mw', scan=scan):,.2f} MW\n" - f"Plasma fusion power: {mfile.get('p_plasma_dt_mw', scan=scan):,.2f} MW \n" + f"Plasma fusion power: {mfile.get('p_plasma_dt_mw', scan=scan):,.2f} MW\n" f"Volume-averaged fusion power density: plasma: {mfile.get('pden_plasma_dt_vol_avg_mw', scan=scan):,.3f} MW/m³\n" f"Beam fusion power: {mfile.get('p_beam_dt_mw', scan=scan):,.2f} MW\n" ) @@ -12293,16 +11604,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): 0.05, 0.75, textstr_dt, - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightyellow", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("lightyellow"), ) axis.text( @@ -12319,23 +11622,15 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): textstr_dd = ( f"Total fusion power: {mfile.get('p_dd_total_mw', scan=scan):,.2f} MW\n" f"Volume-averaged total power density: {mfile.get('pden_dd_total_vol_avg_mw', scan=scan):,.3e} MW/m³\n" - f"Tritium branching ratio: {mfile.get('f_dd_branching_trit', scan=scan):.4f} \n" + f"Tritium branching ratio: {mfile.get('f_dd_branching_trit', scan=scan):.4f}\n" ) axis.text( 0.05, 0.65, textstr_dd, - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightyellow", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("lightyellow"), ) axis.text( @@ -12350,7 +11645,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): # ================================================= textstr_dhe3 = ( - f"Total fusion power: {mfile.get('p_dhe3_total_mw', scan=scan):,.2f} MW \n" + f"Total fusion power: {mfile.get('p_dhe3_total_mw', scan=scan):,.2f} MW\n\n" f"Volume-averaged total power density: {mfile.get('pden_dhe3_total_vol_avg_mw', scan=scan):,.3e} MW/m³\n\n" ) @@ -12358,16 +11653,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): 0.05, 0.55, textstr_dhe3, - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "lightyellow", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("lightyellow"), ) axis.text( @@ -12397,16 +11684,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): 0.05, 0.25, textstr_alpha, - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "red", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("red"), ) axis.text( @@ -12432,16 +11711,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): 0.05, 0.1, textstr_neutron, - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "grey", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("grey"), ) axis.text( @@ -12536,12 +11807,7 @@ def plot_cover_page( ha="left", va="top", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "#e0f7fa", - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("#e0f7fa"), ) # Box 2: File/Branch Info @@ -12559,12 +11825,7 @@ def plot_cover_page( ha="left", va="top", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "#fffde7", - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("#fffde7"), ) # Box 3: Run Settings @@ -12588,12 +11849,7 @@ def plot_cover_page( ha="left", va="top", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "#f3e5f5", - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("#f3e5f5"), ) axis.text( @@ -12690,7 +11946,7 @@ def plot_plasma_pressure_profiles(axis: plt.Axes, mfile: MFile, scan: int): axis.grid(which="minor", linestyle=":", linewidth=0.5, alpha=0.5) axis.set_title("Plasma Thermal Pressure Profiles") axis.grid(True, linestyle="--", alpha=0.5) - axis.set_xlim([0, 1.025]) + axis.set_xlim(0, 1.025) axis.set_ylim(bottom=0) axis.legend() @@ -12770,12 +12026,11 @@ def plot_plasma_current_comparison(axis: plt.Axes, mfile: MFile, scan: int): } # Create the violin plot - axis.violinplot(data.values(), showextrema=False) + data_values = list(data.values()) + axis.violinplot(data_values, showextrema=False) # Create the box plot - axis.boxplot( - data.values(), showfliers=True, showmeans=True, meanline=True, widths=0.3 - ) + axis.boxplot(data_values, showfliers=True, showmeans=True, meanline=True, widths=0.3) # Scatter plot for each data point colors = plt.cm.plasma(np.linspace(0, 1, len(data.values()))) @@ -12815,7 +12070,7 @@ def plot_plasma_current_comparison(axis: plt.Axes, mfile: MFile, scan: int): axis.set_title("Plasma Current ($I_p$) Comparison") axis.set_ylabel(r"Plasma Current [MA]") axis.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{x * 1e-6:.1f}")) - axis.set_xlim([0.5, 1.5]) + axis.set_xlim(0.5, 1.5) axis.set_xticks([]) axis.set_xticklabels([]) axis.set_facecolor("#f0f0f0") @@ -12849,14 +12104,12 @@ def plot_max_normalised_beta_comparison(axis: plt.Axes, mfile: MFile, scan: int) f"{BetaNormMaxModel.THOLERUS.full_name}": beta_norm_max_tholerus, f"{BetaNormMaxModel.STAMBAUGH.full_name}": beta_norm_max_stambaugh, } - + data_values = list(data.values()) # Create the violin plot - axis.violinplot(data.values(), showextrema=False) + axis.violinplot(data_values, showextrema=False) # Create the box plot - axis.boxplot( - data.values(), showfliers=True, showmeans=True, meanline=True, widths=0.3 - ) + axis.boxplot(data_values, showfliers=True, showmeans=True, meanline=True, widths=0.3) # Scatter plot for each data point colors = plt.cm.plasma(np.linspace(0, 1, len(data.values()))) @@ -12895,7 +12148,7 @@ def plot_max_normalised_beta_comparison(axis: plt.Axes, mfile: MFile, scan: int) axis.set_title("Max Normalised Beta ($\\beta_N$) Comparison") axis.set_ylabel("Max Normalised Beta $\\beta_N$ [unitless]") - axis.set_xlim([0.5, 1.5]) + axis.set_xlim(0.5, 1.5) axis.set_xticks([]) axis.set_xticklabels([]) axis.set_facecolor("#f0f0f0") @@ -12946,7 +12199,7 @@ def plot_plasma_pressure_gradient_profiles(axis: plt.Axes, mfile: MFile, scan: i axis.grid(which="minor", linestyle=":", linewidth=0.5, alpha=0.5) axis.set_title("Plasma Thermal Pressure Gradient Profiles") axis.grid(True, linestyle="--", alpha=0.5) - axis.set_xlim([0, 1.025]) + axis.set_xlim(0, 1.025) axis.legend() @@ -12992,7 +12245,6 @@ def plot_plasma_poloidal_pressure_contours(axis: plt.Axes, mfile: MFile, scan: i # 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") # Add colorbar for pressure (now in kPa) @@ -13040,45 +12292,35 @@ def interp1d_profile(profile, mfile: MFile, scan: int): ) # 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 = np.linspace(0, 1, 500) + theta = np.linspace(0, 2 * np.pi, 720) 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 + bdry_r = pg.rs + bdry_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] + bdry_theta = np.arctan2(bdry_z - pg.zs.mean(), bdry_r - pg.rs.mean()) + # Ensure bdry_theta is monotonic and covers [0, 2pi] + bdry_theta = np.unwrap(bdry_theta) + # Sort bdry_theta and corresponding r/z for monotonic interpolation + sort_idx = np.argsort(bdry_theta) + bdry_theta = bdry_theta[sort_idx] + bdry_r = bdry_r[sort_idx] + bdry_z = bdry_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]])) + if bdry_theta[0] > 0 or bdry_theta[-1] < 2 * np.pi: + bdry_theta = np.concatenate(([0], bdry_theta, [2 * np.pi])) + bdry_r = np.concatenate(([bdry_r[0]], bdry_r, [bdry_r[-1]])) + bdry_z = np.concatenate(([bdry_z[0]], bdry_z, [bdry_z[-1]])) # Map theta to boundary r/z f_r = interp1d( - boundary_theta, - boundary_r, - kind="linear", - fill_value="extrapolate", - assume_sorted=True, + bdry_theta, bdry_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, + bdry_theta, bdry_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 @@ -13291,13 +12533,13 @@ def plot_hts_tape_geometry( def plot_tf_corc_cable_summary_box(axis, fig, mfile: MFile, scan: int): """Plot TF CORC cable summary box""" textstr_cable = ( - f"$\\mathbf{{CroCo \\ Cable:}}$\n \n" + f"$\\mathbf{{CroCo \\ Cable:}}$\n\n" f"Cable diameter: {mfile.get('dia_tf_turn_croco_cable', scan=scan) * 1e3:,.4f} mm\n" f"Copper width: {mfile.get('dx_tf_croco_strand_copper', scan=scan) * 1e3:,.4f} mm\n" f"Diameter of solder tape region: {mfile.get('dia_tf_croco_strand_tape_region', scan=scan) * 1e3:,.4f} mm\n" f"Height of tape stack: {mfile.get('dx_tf_croco_strand_tape_stack', scan=scan) * 1e3:,.4f} mm\n" f"Width of HTS tape / tape stack: {mfile.get('dr_tf_hts_tape', scan=scan) * 1e3:,.4f} mm\n" - f"Number of HTS tape layers: {int(mfile.get('n_tf_croco_strand_hts_tapes', scan=scan))}\n \n" + f"Number of HTS tape layers: {int(mfile.get('n_tf_croco_strand_hts_tapes', scan=scan))}\n\n" f"Total copper area: {mfile.get('a_tf_croco_strand_copper_total', scan=scan) * 1e6:,.4f} mm²\n" f"Total hastelloy area: {mfile.get('a_tf_croco_strand_hastelloy', scan=scan) * 1e6:,.4f} mm²\n" f"Total solder area: {mfile.get('a_tf_croco_strand_solder', scan=scan) * 1e6:,.4f} mm²\n" @@ -13313,12 +12555,7 @@ def plot_tf_corc_cable_summary_box(axis, fig, mfile: MFile, scan: int): verticalalignment="top", horizontalalignment="left", transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "#cccccc", # grayish color - "alpha": 1.0, - "linewidth": 2, - }, + bbox=_box_style("#cccccc"), # grayish color ) @@ -13350,10 +12587,7 @@ def reaction_plot_grid( ax.set_xlabel("R [m]") ax.set_xlim(rmajor - 1.2 * rminor, rmajor + 1.2 * rminor) - ax.set_ylim( - -1.2 * rminor * kappa, - 1.2 * kappa * rminor, - ) + ax.set_ylim(-1.2 * rminor * kappa, 1.2 * kappa * rminor) ax.set_ylabel("Z [m]") ax.plot( rmajor, @@ -13391,51 +12625,24 @@ def reaction_plot_grid( ax.legend(legend_handles, legend_labels, loc="upper right", fontsize=8) -def plot_fusion_rate_contours( - fig1, - fig2, - mfile: MFile, - scan: int, -): +def plot_fusion_rate_contours(fig1, fig2, mfile: MFile, scan: int): """Plot fusion rate density contours""" - fusden_plasma_dt_profile = [] - fusden_plasma_dd_triton_profile = [] - fusden_plasma_dd_helion_profile = [] - fusden_plasma_dhe3_profile = [] - rmajor = mfile.get("rmajor", scan=scan) rminor = mfile.get("rminor", scan=scan) kappa = mfile.get("kappa", scan=scan) n_plasma_profile_elements = int(mfile.get("n_plasma_profile_elements", scan=scan)) - fusden_plasma_dt_profile = [ - mfile.get(f"fusden_plasma_dt_profile{i}", scan=scan) - for i in range(n_plasma_profile_elements) - ] - - fusden_plasma_dd_triton_profile = [ - mfile.get(f"fusden_plasma_dd_triton_profile{i}", scan=scan) - for i in range(n_plasma_profile_elements) - ] - - fusden_plasma_dd_helion_profile = [ - mfile.get(f"fusden_plasma_dd_helion_profile{i}", scan=scan) - for i in range(n_plasma_profile_elements) - ] - fusden_plasma_dhe3_profile = [ - mfile.get(f"fusden_plasma_dhe3_profile{i}", scan=scan) - for i in range(n_plasma_profile_elements) - ] - - dt_grid, _r_grid, _z_grid = interp1d_profile(fusden_plasma_dt_profile, mfile, scan) + def fusrat(name): + fusrat_dat = [ + mfile.get(f"fusrat_plasma_{name}_profile{i}", scan=scan) + for i in range(n_plasma_profile_elements) + ] + return interp1d_profile(fusrat_dat, mfile, scan) - dd_triton_grid, _r_grid, _z_grid = interp1d_profile( - fusden_plasma_dd_triton_profile, mfile, scan - ) - dd_helion_grid, _r_grid, _z_grid = interp1d_profile( - fusden_plasma_dd_helion_profile, mfile, scan - ) - dhe3_grid, r_grid, z_grid = interp1d_profile(fusden_plasma_dhe3_profile, mfile, scan) + dt_grid, _r_grid, _z_grid = fusrat("dt") + dd_triton_grid, _r_grid, _z_grid = fusrat(" dd_triton ") + dd_helion_grid, _r_grid, _z_grid = fusrat(" dd_helion ") + dhe3_grid, r_grid, z_grid = fusrat(" dhe3") dt_axes = fig1.add_subplot(121, aspect="equal") dd_triton_axes = fig1.add_subplot(122, aspect="equal") @@ -13512,7 +12719,7 @@ def plot_magnetic_fields_in_plasma(axis: plt.Axes, mfile: MFile, scan: int): verticalalignment="center", horizontalalignment="center", transform=axis.transAxes, - bbox={"boxstyle": "round", "facecolor": "wheat", "alpha": 1.0, "linewidth": 2}, + bbox=_box_style("wheat"), ) # Text box for outboard toroidal field @@ -13524,7 +12731,7 @@ def plot_magnetic_fields_in_plasma(axis: plt.Axes, mfile: MFile, scan: int): verticalalignment="center", horizontalalignment="center", transform=axis.transAxes, - bbox={"boxstyle": "round", "facecolor": "wheat", "alpha": 1.0, "linewidth": 2}, + bbox=_box_style("wheat"), ) axis.set_xlabel("Radial Position [m]") @@ -14135,7 +13342,7 @@ def plot_debye_length_profile(axis: plt.Axes, mfile_data: MFile, scan: int): axis.set_xlabel("$\\rho \\ [r/a]$") axis.grid(True, which="both", linestyle="--", alpha=0.5) - axis.set_xlim([0, 1.025]) + axis.set_xlim(0, 1.025) axis.minorticks_on() axis.legend() @@ -14201,7 +13408,7 @@ def plot_velocity_profile(axis: plt.Axes, mfile_data: MFile, scan: int) -> None: axis.set_ylabel("Velocity [m/s]") axis.set_xlabel("$\\rho \\ [r/a]$") axis.grid(True, which="both", linestyle="--", alpha=0.5) - axis.set_xlim([0, 1.025]) + axis.set_xlim(0, 1.025) axis.minorticks_on() axis.legend() @@ -15000,20 +14207,14 @@ def plot_blkt_structure( r_fw_inboard_out = r_blkt_inboard_in + dr_blkt_inboard + dr_fw_inboard # Plot a horizontal line at dz_blkt_half (blanket half height) - ax.axhline( - dz_blkt_half, - color="purple", - linestyle="--", - linewidth=1.5, - label="Blanket Half Height", - ) - ax.axhline( - -dz_blkt_half, - color="purple", - linestyle="--", - linewidth=1.5, - label="Blanket Half Height", - ) + for dz_blkt in (dz_blkt_half, -dz_blkt_half): + ax.axhline( + dz_blkt, + color="purple", + linestyle="--", + linewidth=1.5, + label="Blanket Half Height", + ) if DivertorNumberModels(i_single_null) == DivertorNumberModels.DOUBLE_NULL: # Plot arrows for the outboard blanket angles @@ -15077,21 +14278,14 @@ def plot_blkt_structure( ) # Plot arrows for the inboard blanket angles - ax.annotate( - "", - xy=(rmajor, 0), - xytext=(r_fw_inboard_out, dz_blkt_half), - arrowprops={"arrowstyle": "<-", "color": "green"}, - zorder=5, - ) - - ax.annotate( - "", - xy=(rmajor, 0), - xytext=(r_fw_inboard_out, -dz_blkt_half), - arrowprops={"arrowstyle": "<-", "color": "green"}, - zorder=5, - ) + for dz_blkt in (dz_blkt_half, -dz_blkt_half): + ax.annotate( + "", + xy=(rmajor, 0), + xytext=(r_fw_inboard_out, dz_blkt), + arrowprops={"arrowstyle": "<-", "color": "green"}, + zorder=5, + ) # Plot arc showing the angle between the two inboard blanket arrows arc_radius = 1.0 @@ -15174,7 +14368,7 @@ def plot_blkt_structure( arc_radius = 1.5 # 3 to 6 o'clock is -90 degrees angle_start = -90.0 - angle_end = -90.0 - deg_div_poloidal_plasma + angle_end = angle_start - deg_div_poloidal_plasma theta = np.linspace(np.deg2rad(angle_start), np.deg2rad(angle_end), 50) arc_x = rmajor + arc_radius * np.cos(theta) @@ -15209,40 +14403,26 @@ def plot_blkt_structure( ) # Plot vertical lines at the inner and outer radial boundaries of the blanket - ax.axvline( - r_blkt_inboard_in, color="black", linestyle="--", linewidth=1.5, zorder=10 - ) - ax.axvline( - r_blkt_outboard_out, color="black", linestyle="--", linewidth=1.5, zorder=10 - ) - ax.axvline(r_fw_inboard_out, color="black", linestyle="--", linewidth=1.5, zorder=10) - - ax.axvline(r_fw_outboard_in, color="black", linestyle="--", linewidth=1.5, zorder=10) + linestyle = {"color": "black", "linestyle": "--", "linewidth": 1.5, "zorder": 10} + ax.axvline(r_blkt_inboard_in, **linestyle) + ax.axvline(r_blkt_outboard_out, **linestyle) + ax.axvline(r_fw_inboard_out, **linestyle) + ax.axvline(r_fw_outboard_in, **linestyle) ax.axvline( - rmajor, - color="black", - linestyle="--", - linewidth=1.5, - label="Major Radius $R_0$", + rmajor, color="black", linestyle="--", linewidth=1.5, label="Major Radius $R_0$" ) # Plot midplane line (horizontal dashed line at Z=0) - ax.axhline( - 0.0, - color="black", - linestyle="--", - linewidth=1.5, - label="Midplane", - ) + ax.axhline(0.0, color="black", linestyle="--", linewidth=1.5, label="Midplane") textstr_blkt_areas = ( f"$\\mathbf{{Blanket \\ Areas:}}$\n\n" - f"Inboard blanket, with holes and gaps: {m_file.get('a_blkt_inboard_surface', scan=scan):,.3f} $\\text{{m}}^2$ \n" - f"Outboard blanket, with holes and gaps: {m_file.get('a_blkt_outboard_surface', scan=scan):,.3f} $\\text{{m}}^2$ \n" - f"Total blanket, with holes and gaps: {m_file.get('a_blkt_total_surface', scan=scan):,.3f} $\\text{{m}}^2$ \n\n" - f"Inboard blanket, full coverage: {m_file.get('a_blkt_inboard_surface_full_coverage', scan=scan):,.3f} $\\text{{m}}^2$ \n" - f"Outboard blanket, full coverage: {m_file.get('a_blkt_outboard_surface_full_coverage', scan=scan):,.3f} $\\text{{m}}^2$ \n" + f"Inboard blanket, with holes and gaps: {m_file.get('a_blkt_inboard_surface', scan=scan):,.3f} $\\text{{m}}^2$\n" + f"Outboard blanket, with holes and gaps: {m_file.get('a_blkt_outboard_surface', scan=scan):,.3f} $\\text{{m}}^2$\n" + f"Total blanket, with holes and gaps: {m_file.get('a_blkt_total_surface', scan=scan):,.3f} $\\text{{m}}^2$\n\n" + f"Inboard blanket, full coverage: {m_file.get('a_blkt_inboard_surface_full_coverage', scan=scan):,.3f} $\\text{{m}}^2$\n" + f"Outboard blanket, full coverage: {m_file.get('a_blkt_outboard_surface_full_coverage', scan=scan):,.3f} $\\text{{m}}^2$\n" f"Total blanket, full coverage: {m_file.get('a_blkt_total_surface_full_coverage', scan=scan):,.3f} $\\text{{m}}^2$ " ) @@ -15250,25 +14430,17 @@ def plot_blkt_structure( 0.05, 0.3, textstr_blkt_areas, - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "wheat", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("wheat"), ) textstr_blkt_volumes = ( f"$\\mathbf{{Blanket \\ Volumes:}}$\n\n" - f"Inboard blanket, with holes and gaps: {m_file.get('vol_blkt_inboard', scan=scan):,.3f} $\\text{{m}}^3$ \n" - f"Outboard blanket, with holes and gaps: {m_file.get('vol_blkt_outboard', scan=scan):,.3f} $\\text{{m}}^3$ \n" - f"Total blanket, with holes and gaps: {m_file.get('vol_blkt_total', scan=scan):,.3f} $\\text{{m}}^3$ \n\n" - f"Inboard blanket, full coverage: {m_file.get('vol_blkt_inboard_full_coverage', scan=scan):,.3f} $\\text{{m}}^3$ \n" - f"Outboard blanket, full coverage: {m_file.get('vol_blkt_outboard_full_coverage', scan=scan):,.3f} $\\text{{m}}^3$ \n" + f"Inboard blanket, with holes and gaps: {m_file.get('vol_blkt_inboard', scan=scan):,.3f} $\\text{{m}}^3$\n" + f"Outboard blanket, with holes and gaps: {m_file.get('vol_blkt_outboard', scan=scan):,.3f} $\\text{{m}}^3$\n" + f"Total blanket, with holes and gaps: {m_file.get('vol_blkt_total', scan=scan):,.3f} $\\text{{m}}^3$\n\n" + f"Inboard blanket, full coverage: {m_file.get('vol_blkt_inboard_full_coverage', scan=scan):,.3f} $\\text{{m}}^3$\n" + f"Outboard blanket, full coverage: {m_file.get('vol_blkt_outboard_full_coverage', scan=scan):,.3f} $\\text{{m}}^3$\n" f"Total blanket, full coverage: {m_file.get('vol_blkt_total_full_coverage', scan=scan):,.3f} $\\text{{m}}^3$ " ) @@ -15276,16 +14448,8 @@ def plot_blkt_structure( 0.05, 0.05, textstr_blkt_volumes, - fontsize=9, - verticalalignment="bottom", - horizontalalignment="left", - transform=fig.transFigure, - bbox={ - "boxstyle": "round", - "facecolor": "wheat", - "alpha": 1.0, - "linewidth": 2, - }, + **_text_layout(fig), + bbox=_box_style("wheat"), ) @@ -16018,11 +15182,7 @@ def plot_cs_radial_stress_contour_profile( linewidth=2, ) - # Use a dedicated colorbar axes when provided so the main axes width is unchanged. - if colorbar_axis is None: - cbar = axis.figure.colorbar(contour_fill, ax=axis, pad=0.02) - else: - cbar = axis.figure.colorbar(contour_fill, cax=colorbar_axis) + cbar = _colourbar(contour_fill, axis, colorbar_axis) cbar.set_label("Radial Stress (MPa)") axis.set_xlabel("R [m]") @@ -16112,11 +15272,7 @@ def plot_cs_hoop_stress_contour_profile( linewidth=2, ) - # Use a dedicated colorbar axes when provided so the main axes width is unchanged. - if colorbar_axis is None: - cbar = axis.figure.colorbar(contour_fill, ax=axis, pad=0.02) - else: - cbar = axis.figure.colorbar(contour_fill, cax=colorbar_axis) + cbar = _colourbar(contour_fill, axis, colorbar_axis) cbar.set_label("Hoop Stress (MPa)") axis.set_xlabel("R [m]") @@ -16226,11 +15382,7 @@ def plot_vertical_stress_contour_profile( linewidth=2, ) - # Use a dedicated colorbar axes when provided so the main axes width is unchanged. - if colorbar_axis is None: - cbar = axis.figure.colorbar(contour_fill, ax=axis, pad=0.02) - else: - cbar = axis.figure.colorbar(contour_fill, cax=colorbar_axis) + cbar = _colourbar(contour_fill, axis, colorbar_axis) cbar.set_label("Vertical Stress (MPa)") axis.set_xlabel("R [m]") @@ -16299,7 +15451,6 @@ def plot_cs_tresca_2d_contour( ) # Plot filled contour of Tresca stress distribution - contour_fill = axis.contourf(r, z, tresca_data, levels=15, cmap="RdYlBu_r") contour_lines = axis.contour( r, z, @@ -16336,11 +15487,8 @@ def plot_cs_tresca_2d_contour( linewidth=2, ) - # Use a dedicated colorbar axes when provided so the main axes width is unchanged. - if colorbar_axis is None: - cbar = axis.figure.colorbar(contour_fill, ax=axis, pad=0.02) - else: - cbar = axis.figure.colorbar(contour_fill, cax=colorbar_axis) + contour_fill = axis.contourf(r, z, tresca_data, levels=15, cmap="RdYlBu_r") + cbar = _colourbar(contour_fill, axis, colorbar_axis) cbar.set_label("Tresca Stress (MPa)") axis.set_xlabel("R [m]") @@ -16413,7 +15561,6 @@ def plot_cs_von_mises_2d_contour( ) # Plot filled contour of Von Mises stress distribution - contour_fill = axis.contourf(r, z, von_mises_data, levels=15, cmap="RdYlBu_r") contour_lines = axis.contour( r, z, @@ -16450,11 +15597,8 @@ def plot_cs_von_mises_2d_contour( linewidth=2, ) - # Use a dedicated colorbar axes when provided so the main axes width is unchanged. - if colorbar_axis is None: - cbar = axis.figure.colorbar(contour_fill, ax=axis, pad=0.02) - else: - cbar = axis.figure.colorbar(contour_fill, cax=colorbar_axis) + contour_fill = axis.contourf(r, z, von_mises_data, levels=15, cmap="RdYlBu_r") + cbar = _colourbar(contour_fill, axis, colorbar_axis) cbar.set_label("Von Mises Stress (MPa)") axis.set_xlabel("R [m]") @@ -17220,14 +16364,14 @@ def _add_page(name: str | None = None): demo_ranges, colour_scheme, ) - ax_full_toroidal.set_ylim([ + ax_full_toroidal.set_ylim( -ax_full_toroidal.get_ylim()[1], ax_full_toroidal.get_ylim()[1], - ]) - ax_full_toroidal.set_xlim([ + ) + ax_full_toroidal.set_xlim( -ax_full_toroidal.get_xlim()[1], ax_full_toroidal.get_xlim()[1], - ]) + ) ax18 = _add_page().add_subplot(211) ax18.set_position([0.1, 0.33, 0.8, 0.6]) diff --git a/process/models/stellarator/neoclassics.py b/process/models/stellarator/neoclassics.py index b1e00bd16d..975c4ee1d9 100644 --- a/process/models/stellarator/neoclassics.py +++ b/process/models/stellarator/neoclassics.py @@ -4,6 +4,7 @@ from dataclasses import dataclass import numpy as np +from numpy.polynomial.polynomial import polyval from process.core import constants from process.core.model import Model @@ -30,9 +31,14 @@ class Neoclassics(Model): """Module containing neoclassics routines""" @property - def no_roots(self): - """Obtain number of Gauss Laguerre roots""" - return self.data.neoclassics.roots.shape[0] + def mass(self): + """Component mass array for electron, deuterium, tritium and Helium""" + return np.array([ + constants.ELECTRON_MASS, + constants.PROTON_MASS * 2.0, + constants.PROTON_MASS * 3.0, + constants.PROTON_MASS * 4.0, + ]) def output(self): """Neoclassics model doesn't have any output""" @@ -156,142 +162,65 @@ def init_neoclassics(self, r_effin, eps_effin, iotain): ) self.data.neoclassics.q_flux = self.neoclassics_calc_q_flux() - def init_profile_values_from_PROCESS(self, rho): + def init_profile_values_from_PROCESS(self, rho) -> tuple[np.ndarray, ...]: """Initialises the profile_values object from PROCESS' parabolic profiles Parameters ---------- rho : + Returns + ------- + dens: + density of electron, deuterium, tritum and alpha [m⁻³] + temp: + temperature of electron, deuterium, tritum and alpha [J] + dr_dens: + derivative density of electron, deuterium, tritum and alpha [m⁻³] + dr_temp: + derivative temperature of electron, deuterium, tritum and alpha [J] """ - tempe = ( - self.data.physics.temp_plasma_electron_on_axis_kev - * (1 - rho**2) ** self.data.physics.alphat - * KEV - ) - tempT = ( - self.data.physics.temp_plasma_ion_on_axis_kev - * (1 - rho**2) ** self.data.physics.alphat - * KEV - ) - tempD = ( - self.data.physics.temp_plasma_ion_on_axis_kev - * (1 - rho**2) ** self.data.physics.alphat - * KEV - ) - tempa = ( - self.data.physics.temp_plasma_ion_on_axis_kev - * (1 - rho**2) ** self.data.physics.alphat - * KEV - ) + t_suffix = (1 - rho**2) ** self.data.physics.alphat * KEV - dense = ( - self.data.physics.nd_plasma_electron_on_axis - * (1 - rho**2) ** self.data.physics.alphan - ) - densT = ( - (1 - self.data.physics.f_plasma_fuel_deuterium) - * self.data.physics.nd_plasma_ions_on_axis - * (1 - rho**2) ** self.data.physics.alphan - ) - densD = ( + dens_suffix = (1 - rho**2) ** self.data.physics.alphan + + deriv_prefix = -2.0 * 1.0 / self.data.physics.rminor * rho + + # array are set up as electron, deuterium, tritum, alpha + density_array = np.array([ + self.data.physics.nd_plasma_electron_on_axis, self.data.physics.f_plasma_fuel_deuterium - * self.data.physics.nd_plasma_ions_on_axis - * (1 - rho**2) ** self.data.physics.alphan - ) - densa = ( + * self.data.physics.nd_plasma_ions_on_axis, + (1 - self.data.physics.f_plasma_fuel_deuterium) + * self.data.physics.nd_plasma_ions_on_axis, self.data.physics.nd_plasma_alphas_thermal_vol_avg - * (1 + self.data.physics.alphan) - * (1 - rho**2) ** self.data.physics.alphan - ) + * (1 + self.data.physics.alphan), + ]) + + temperature_array = np.array([ + self.data.physics.temp_plasma_electron_on_axis_kev, + self.data.physics.temp_plasma_ion_on_axis_kev, + self.data.physics.temp_plasma_ion_on_axis_kev, + self.data.physics.temp_plasma_ion_on_axis_kev, + ]) # Derivatives in real space - dr_tempe = ( - -2.0 - * 1.0 - / self.data.physics.rminor - * self.data.physics.temp_plasma_electron_on_axis_kev - * rho - * (1.0 - rho**2) ** (self.data.physics.alphat - 1.0) - * self.data.physics.alphat - * KEV - ) - dr_tempT = ( - -2.0 - * 1.0 - / self.data.physics.rminor - * self.data.physics.temp_plasma_ion_on_axis_kev - * rho - * (1.0 - rho**2) ** (self.data.physics.alphat - 1.0) - * self.data.physics.alphat - * KEV - ) - dr_tempD = ( - -2.0 - * 1.0 - / self.data.physics.rminor - * self.data.physics.temp_plasma_ion_on_axis_kev - * rho - * (1.0 - rho**2) ** (self.data.physics.alphat - 1.0) - * self.data.physics.alphat - * KEV + dr_dens = ( + deriv_prefix + * density_array + * (1.0 - rho**2) ** (self.data.physics.alphan - 1.0) + * self.data.physics.alphan ) - dr_tempa = ( - -2.0 - * 1.0 - / self.data.physics.rminor - * self.data.physics.temp_plasma_ion_on_axis_kev - * rho + dr_temp = ( + deriv_prefix + * temperature_array * (1.0 - rho**2) ** (self.data.physics.alphat - 1.0) * self.data.physics.alphat * KEV ) - dr_dense = ( - -2.0 - * 1.0 - / self.data.physics.rminor - * rho - * self.data.physics.nd_plasma_electron_on_axis - * (1.0 - rho**2) ** (self.data.physics.alphan - 1.0) - * self.data.physics.alphan - ) - dr_densT = ( - -2.0 - * 1.0 - / self.data.physics.rminor - * rho - * (1 - self.data.physics.f_plasma_fuel_deuterium) - * self.data.physics.nd_plasma_ions_on_axis - * (1.0 - rho**2) ** (self.data.physics.alphan - 1.0) - * self.data.physics.alphan - ) - dr_densD = ( - -2.0 - * 1.0 - / self.data.physics.rminor - * rho - * self.data.physics.f_plasma_fuel_deuterium - * self.data.physics.nd_plasma_ions_on_axis - * (1.0 - rho**2) ** (self.data.physics.alphan - 1.0) - * self.data.physics.alphan - ) - dr_densa = ( - -2.0 - * 1.0 - / self.data.physics.rminor - * rho - * self.data.physics.nd_plasma_alphas_thermal_vol_avg - * (1 + self.data.physics.alphan) - * (1.0 - rho**2) ** (self.data.physics.alphan - 1.0) - * self.data.physics.alphan - ) - - dens = np.array([dense, densD, densT, densa]) - temp = np.array([tempe, tempD, tempT, tempa]) - dr_dens = np.array([dr_dense, dr_densD, dr_densT, dr_densa]) - dr_temp = np.array([dr_tempe, dr_tempD, dr_tempT, dr_tempa]) - + dens = density_array * dens_suffix + temp = temperature_array * t_suffix return dens, temp, dr_dens, dr_temp def calc_neoclassics(self): @@ -307,16 +236,6 @@ def calc_neoclassics(self): self.data.stellarator.iotabar, ) - q_PROCESS = ( - ( - self.data.physics.f_p_alpha_plasma_deposited - * self.data.physics.pden_alpha_total_vol_avg_mw - - self.data.physics.pden_plasma_core_rad_mw - ) - * self.data.physics.vol_plasma - / self.data.physics.a_plasma_surface - * self.data.impurity_radiation.radius_plasma_core_norm - ) q_PROCESS_r1 = ( ( self.data.physics.f_p_alpha_plasma_deposited @@ -327,81 +246,25 @@ def calc_neoclassics(self): / self.data.physics.a_plasma_surface ) - q_neo = sum(self.data.neoclassics.q_flux * 1e-6) - gamma_neo = sum( - self.data.neoclassics.gamma_flux * self.data.neoclassics.temperatures * 1e-6 - ) - - total_q_neo = sum( - self.data.neoclassics.q_flux * 1e-6 - + self.data.neoclassics.gamma_flux - * self.data.neoclassics.temperatures - * 1e-6 - ) - - total_q_neo_e = ( - 2 - * 2 - * ( - self.data.neoclassics.q_flux[0] * 1e-6 - + self.data.neoclassics.gamma_flux[0] - * self.data.neoclassics.temperatures[0] - * 1e-6 - ) - ) - - q_neo_e = self.data.neoclassics.q_flux[0] * 1e-6 - q_neo_D = self.data.neoclassics.q_flux[1] * 1e-6 - q_neo_a = self.data.neoclassics.q_flux[3] * 1e-6 - q_neo_T = self.data.neoclassics.q_flux[2] * 1e-6 - - g_neo_e = ( - self.data.neoclassics.gamma_flux[0] - * 1e-6 - * self.data.neoclassics.temperatures[0] - ) - g_neo_D = ( - self.data.neoclassics.gamma_flux[1] - * 1e-6 - * self.data.neoclassics.temperatures[1] - ) - g_neo_a = ( - self.data.neoclassics.gamma_flux[3] - * 1e-6 - * self.data.neoclassics.temperatures[3] - ) - g_neo_T = ( - self.data.neoclassics.gamma_flux[2] - * 1e-6 - * self.data.neoclassics.temperatures[2] - ) + q_PROCESS = q_PROCESS_r1 * self.data.impurity_radiation.radius_plasma_core_norm dndt_neo_e = self.data.neoclassics.gamma_flux[0] - dndt_neo_D = self.data.neoclassics.gamma_flux[1] - dndt_neo_a = self.data.neoclassics.gamma_flux[3] - dndt_neo_T = self.data.neoclassics.gamma_flux[2] + q_neo_e = self.data.neoclassics.q_flux[0] * 1e-6 + g_neo_e = 1e-6 * dndt_neo_e * self.data.neoclassics.temperatures[0] + total_q_neo_e = 4 * (q_neo_e + g_neo_e) - dndt_neo_fuel = ( - (dndt_neo_D + dndt_neo_T) - * self.data.physics.a_plasma_surface - * self.data.impurity_radiation.radius_plasma_core_norm - ) - dmdt_neo_fuel = ( - dndt_neo_fuel * self.data.physics.m_fuel_amu * constants.PROTON_MASS * 1.0e6 - ) # mg dmdt_neo_fuel_from_e = ( - 4 + 4e6 * dndt_neo_e * self.data.physics.a_plasma_surface * self.data.impurity_radiation.radius_plasma_core_norm * self.data.physics.m_fuel_amu * constants.PROTON_MASS - * 1.0e6 ) # kg chi_neo_e = -( self.data.neoclassics.q_flux[0] - + self.data.neoclassics.gamma_flux[0] * self.data.neoclassics.temperatures[0] + + dndt_neo_e * self.data.neoclassics.temperatures[0] ) / ( self.data.neoclassics.densities[0] * self.data.neoclassics.dr_temperatures[0] + self.data.neoclassics.temperatures[0] @@ -409,57 +272,52 @@ def calc_neoclassics(self): ) chi_PROCESS_e = self.st_calc_eff_chi() - - nu_star = NormalisedCollisionality( - e=self.data.neoclassics.nu_star_averaged[0], - D=self.data.neoclassics.nu_star_averaged[1], - T=self.data.neoclassics.nu_star_averaged[2], - He=self.data.neoclassics.nu_star_averaged[3], - ) - return ( q_PROCESS, q_PROCESS_r1, - q_neo, - gamma_neo, - total_q_neo, total_q_neo_e, q_neo_e, - q_neo_D, - q_neo_a, - q_neo_T, g_neo_e, - g_neo_D, - g_neo_a, - g_neo_T, dndt_neo_e, - dndt_neo_D, - dndt_neo_a, - dndt_neo_T, - dndt_neo_fuel, - dmdt_neo_fuel, dmdt_neo_fuel_from_e, chi_neo_e, chi_PROCESS_e, - nu_star, + NormalisedCollisionality( + e=self.data.neoclassics.nu_star_averaged[0], + D=self.data.neoclassics.nu_star_averaged[1], + T=self.data.neoclassics.nu_star_averaged[2], + He=self.data.neoclassics.nu_star_averaged[3], + ), ) def neoclassics_calc_KT(self): """Calculates the energy on the given grid which is given by the gauss laguerre roots. """ - k = np.repeat((self.data.neoclassics.roots / KEV)[:, np.newaxis], 4, axis=1) - - return (k * self.data.neoclassics.temperatures).T + return ( + self.data.neoclassics.roots[None] / KEV + ) * self.data.neoclassics.temperatures[:, None] + + @staticmethod + def _nu_erf(xk, expxk): + """Error function""" + # Rational approximation for erf + # t term + t = 1.0 / (1.0 + 0.3275911 * np.sqrt(xk)) + # Expand: t * (c0 + c1*t + c2*t^2 + ...) + coeffs = np.array([ + 0, + 0.254829592, + -0.284496736, + 1.421413741, + -1.453152027, + 1.061405429, + ]) + return 1.0 - polyval(t, coeffs) * expxk def neoclassics_calc_nu(self): """Calculates the collision frequency""" - mass = np.array([ - constants.ELECTRON_MASS, - constants.PROTON_MASS * 2.0, - constants.PROTON_MASS * 3.0, - constants.PROTON_MASS * 4.0, - ]) + mass = self.mass z = np.array([-1.0, 1.0, 1.0, 2.0]) * constants.ELECTRON_CHARGE # transform the temperature back in eV @@ -472,76 +330,41 @@ def neoclassics_calc_nu(self): * np.log10(self.data.neoclassics.temperatures[0] / constants.ELECTRON_CHARGE) ) - neoclassics_calc_nu = np.zeros((4, self.no_roots), order="F") - - for j in range(4): - for i in range(self.no_roots): - x = self.data.neoclassics.roots[i] - for k in range(4): - xk = ( - (mass[k] / mass[j]) - * ( - self.data.neoclassics.temperatures[j] - / self.data.neoclassics.temperatures[k] - ) - * x - ) - expxk = np.exp(-xk) - t = 1.0 / (1.0 + 0.3275911 * np.sqrt(xk)) - erfn = ( - 1.0 - - t - * ( - 0.254829592 - + t - * ( - -0.284496736 - + t - * (1.421413741 + t * (-1.453152027 + t * 1.061405429)) - ) - ) - * expxk - ) - phixmgx = (1.0 - 0.5 / xk) * erfn + expxk / np.sqrt(np.pi * xk) - v = np.sqrt( - 2.0 * x * self.data.neoclassics.temperatures[j] / mass[j] - ) - neoclassics_calc_nu[j, i] += ( - self.data.neoclassics.densities[k] - * (z[j] * z[k]) ** 2 - * lnlambda - * phixmgx - / (4.0 * np.pi * constants.EPSILON0**2 * mass[j] ** 2 * v**3) - ) - - return neoclassics_calc_nu + roots = self.data.neoclassics.roots + temp = self.data.neoclassics.temperatures - def neoclassics_calc_nu_star(self): - """Calculates the normalized collision frequency""" - k = np.repeat(self.data.neoclassics.roots[:, np.newaxis], 4, axis=1) - kk = (k * self.data.neoclassics.temperatures).T + inv_mass = 1 / mass + xk = np.einsum("k,j,j,k,r->jkr", mass, inv_mass, temp, 1 / temp, roots) - mass = np.array([ - constants.ELECTRON_MASS, - constants.PROTON_MASS * 2.0, - constants.PROTON_MASS * 3.0, - constants.PROTON_MASS * 4.0, - ]) + expxk = np.exp(-xk) - v = np.empty((4, self.no_roots)) - v[0, :] = constants.SPEED_LIGHT * np.sqrt( - 1.0 - (kk[0, :] / (mass[0] * constants.SPEED_LIGHT**2) + 1) ** (-1) - ) - v[1, :] = constants.SPEED_LIGHT * np.sqrt( - 1.0 - (kk[1, :] / (mass[1] * constants.SPEED_LIGHT**2) + 1) ** (-1) - ) - v[2, :] = constants.SPEED_LIGHT * np.sqrt( - 1.0 - (kk[2, :] / (mass[2] * constants.SPEED_LIGHT**2) + 1) ** (-1) + erfn = self._nu_erf(xk, expxk) + + phixmgx = (1.0 - 0.5 * 1 / xk) * erfn + expxk / np.sqrt(np.pi * xk) + + v = np.sqrt(2.0 * np.einsum("r,j,j->jr", roots, temp, inv_mass)) + denom = 1 / (4 * np.pi * constants.EPSILON0**2 * mass[:, None] ** 2 * v**3) + + return lnlambda * np.einsum( + "k,jk,jkr,jr->jr", + self.data.neoclassics.densities, + np.einsum("i,j->ij", z, z) ** 2, + phixmgx, + denom, ) - v[3, :] = constants.SPEED_LIGHT * np.sqrt( - 1.0 - (kk[3, :] / (mass[3] * constants.SPEED_LIGHT**2) + 1) ** (-1) + + def neoclassics_calc_nu_star(self): + """Calculates the normalized collision frequency""" + kk = ( + self.data.neoclassics.roots[None] + * self.data.neoclassics.temperatures[:, None] ) + mass = self.mass + + v = constants.SPEED_LIGHT * np.sqrt( + 1.0 - (kk / (mass[:, None] * constants.SPEED_LIGHT**2) + 1) ** -1 + ) return ( self.data.physics.rmajor * self.data.neoclassics.nu @@ -557,15 +380,13 @@ def neoclassics_calc_nu_star_fromT(self, iota): """ - temp = ( - np.array([ - self.data.physics.temp_plasma_electron_vol_avg_kev, - self.data.physics.temp_plasma_ion_vol_avg_kev, - self.data.physics.temp_plasma_ion_vol_avg_kev, - self.data.physics.temp_plasma_ion_vol_avg_kev, - ]) - * KEV - ) + temp = KEV * np.array([ + self.data.physics.temp_plasma_electron_vol_avg_kev, + self.data.physics.temp_plasma_ion_vol_avg_kev, + self.data.physics.temp_plasma_ion_vol_avg_kev, + self.data.physics.temp_plasma_ion_vol_avg_kev, + ]) + density = np.array([ self.data.physics.nd_plasma_electrons_vol_avg, self.data.physics.nd_plasma_fuel_ions_vol_avg @@ -575,12 +396,7 @@ def neoclassics_calc_nu_star_fromT(self, iota): self.data.physics.nd_plasma_alphas_thermal_vol_avg, ]) - mass = np.array([ - constants.ELECTRON_MASS, - constants.PROTON_MASS * 2.0, - constants.PROTON_MASS * 3.0, - constants.PROTON_MASS * 4.0, - ]) + mass = self.mass z = np.array([-1.0, 1.0, 1.0, 2.0]) * constants.ELECTRON_CHARGE # transform the temperature back in eV @@ -592,120 +408,57 @@ def neoclassics_calc_nu_star_fromT(self, iota): + 2.3 * np.log10(temp[0] / constants.ELECTRON_CHARGE) ) - neoclassics_calc_nu_star_fromT = np.zeros((4,)) - - for j in range(4): - v = np.sqrt(2.0 * temp[j] / mass[j]) - for k in range(4): - xk = (mass[k] / mass[j]) * (temp[j] / temp[k]) - - expxk = 0.0 - if xk < 200.0: - expxk = np.exp(-xk) - - t = 1.0 / (1.0 + 0.3275911 * np.sqrt(xk)) - erfn = ( - 1.0 - - t - * ( - 0.254829592 - + t - * ( - -0.284496736 - + t * (1.421413741 + t * (-1.453152027 + t * 1.061405429)) - ) - ) - * expxk - ) - phixmgx = (1.0 - 0.5 / xk) * erfn + expxk / np.sqrt(np.pi * xk) - neoclassics_calc_nu_star_fromT[j] += ( - density[k] - * (z[j] * z[k]) ** 2 - * lnlambda - * phixmgx - / (4.0 * np.pi * constants.EPSILON0**2 * mass[j] ** 2 * v**4) - * self.data.physics.rmajor - / iota - ) - return neoclassics_calc_nu_star_fromT + inv_mass = 1 / mass + v = np.sqrt(2.0 * temp * inv_mass) + xk = np.einsum("k,j,j,k->jk", mass, inv_mass, temp, 1 / temp) - def neoclassics_calc_vd(self): - """Calculates the drift velocity on GL roots""" - vde = ( - self.data.neoclassics.roots - * self.data.neoclassics.temperatures[0] - / ( - constants.ELECTRON_CHARGE - * self.data.physics.rmajor - * self.data.physics.b_plasma_toroidal_on_axis - ) - ) - vdD = ( - self.data.neoclassics.roots - * self.data.neoclassics.temperatures[1] - / ( - constants.ELECTRON_CHARGE - * self.data.physics.rmajor - * self.data.physics.b_plasma_toroidal_on_axis + # exp(-xk), clipped at xk >= 200 + mask_xk_lt_200 = xk < 200.0 + expxk = np.zeros_like(xk) + expxk[mask_xk_lt_200] = np.exp(-xk[mask_xk_lt_200]) + + erfn = self._nu_erf(xk, expxk) + + phixmgx = (1.0 - 0.5 * 1 / xk) * erfn + expxk / np.sqrt(np.pi * xk) + + denom = 1 / (4 * np.pi * constants.EPSILON0**2 * mass**2 * v**4) + + # sum over k dimension + return ( + lnlambda + * self.data.physics.rmajor + / iota + * np.einsum( + "k,jk,jk,j->j", density, np.einsum("i,j->ij", z, z) ** 2, phixmgx, denom ) ) - vdT = ( - self.data.neoclassics.roots - * self.data.neoclassics.temperatures[2] - / ( - constants.ELECTRON_CHARGE - * self.data.physics.rmajor - * self.data.physics.b_plasma_toroidal_on_axis - ) + + def neoclassics_calc_vd(self): + """Calculates the drift velocity on GL roots""" + # alpha denominator is 2*vd_suffix + vd_suffix = ( + constants.ELECTRON_CHARGE + * self.data.physics.rmajor + * self.data.physics.b_plasma_toroidal_on_axis + * np.array([1, 1, 1, 2]) ) - vda = ( - self.data.neoclassics.roots - * self.data.neoclassics.temperatures[3] - / ( - 2.0 - * constants.ELECTRON_CHARGE - * self.data.physics.rmajor - * self.data.physics.b_plasma_toroidal_on_axis - ) + return ( + self.data.neoclassics.roots[None] + * self.data.neoclassics.temperatures[:, None] + / vd_suffix[:, None] ) - vd = np.empty((4, self.no_roots)) - - vd[0, :] = vde - vd[1, :] = vdD - vd[2, :] = vdT - vd[3, :] = vda - - return vd - def neoclassics_calc_D11_plateau(self): """Calculates the plateau transport coefficients (D11_star sometimes)""" - mass = np.array([ - constants.ELECTRON_MASS, - constants.PROTON_MASS * 2.0, - constants.PROTON_MASS * 3.0, - constants.PROTON_MASS * 4.0, - ]) + mass = self.mass - v = np.empty((4, self.no_roots)) - v[0, :] = constants.SPEED_LIGHT * np.sqrt( + v = constants.SPEED_LIGHT * np.sqrt( 1.0 - - (self.data.neoclassics.kt[0, :] / (mass[0] * constants.SPEED_LIGHT**2) + 1) - ** (-1) - ) - v[1, :] = constants.SPEED_LIGHT * np.sqrt( - 1.0 - - (self.data.neoclassics.kt[1, :] / (mass[1] * constants.SPEED_LIGHT**2) + 1) - ** (-1) - ) - v[2, :] = constants.SPEED_LIGHT * np.sqrt( - 1.0 - - (self.data.neoclassics.kt[2, :] / (mass[2] * constants.SPEED_LIGHT**2) + 1) - ** (-1) - ) - v[3, :] = constants.SPEED_LIGHT * np.sqrt( - 1.0 - - (self.data.neoclassics.kt[3, :] / (mass[3] * constants.SPEED_LIGHT**2) + 1) + - ( + self.data.neoclassics.kt[None] + / (mass[:, None] * constants.SPEED_LIGHT**2) + + 1 + ) ** (-1) ) diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py index d9b58d781c..93d0076094 100644 --- a/process/models/stellarator/stellarator.py +++ b/process/models/stellarator/stellarator.py @@ -2406,24 +2406,10 @@ def st_phys(self, output): ( q_PROCESS, q_PROCESS_r1, - _q_neo, - _gamma_neo, - _total_q_neo, total_q_neo_e, q_neo_e, - _q_neo_D, - _q_neo_a, - _q_neo_T, g_neo_e, - _g_neo_D, - _g_neo_a, - _g_neo_T, dndt_neo_e, - _dndt_neo_D, - _dndt_neo_a, - _dndt_neo_T, - _dndt_neo_fuel, - _dmdt_neo_fuel, dmdt_neo_fuel_from_e, chi_neo_e, chi_PROCESS_e,