From 634db6bdedade9ad90ad0c77d25d262210e7a901 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 28 Sep 2026 10:00:40 +0100 Subject: [PATCH 1/8] Add fusion rate density profile plotting functionality --- process/core/io/plot/summary.py | 178 +++++++++++++++++++++++++++++++- 1 file changed, 177 insertions(+), 1 deletion(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index 15ca2beb11..b97d2a1596 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -16,6 +16,7 @@ from matplotlib.axes import Axes from matplotlib.patches import Circle, FancyBboxPatch, Rectangle from matplotlib.path import Path as mplPath +from scipy.integrate import cumulative_simpson from scipy.interpolate import interp1d from process.core import constants @@ -11557,7 +11558,7 @@ def plot_fw_90_deg_pipe_bend(ax, m_file, scan: int): ) -def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): +def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): """Plot the fusion rate density profiles on the given axis""" n_plasma_profile_elements = int(mfile.get("n_plasma_profile_elements", scan=scan)) @@ -11882,6 +11883,181 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): transform=fig.transFigure, ) +def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): + """Plot the fusion rate profiles on the given axis""" + fusrat_plasma_dt_profile = [] + fusrat_plasma_dd_triton_profile = [] + fusrat_plasma_dd_helion_profile = [] + fusrat_plasma_dhe3_profile = [] + + n_plasma_profile_elements = int(mfile.get("n_plasma_profile_elements", scan=scan)) + vol_plasma = mfile.get("vol_plasma", scan=scan) + + rho, ne, te = profiles_with_pedestal(mfile, scan) + + + + fusrat_plasma_dt_profile = [ + mfile.get(f"fusrat_plasma_dt_profile{i}", scan=scan) + for i in range(n_plasma_profile_elements) + ] + + fusrat_plasma_dd_triton_profile = [ + mfile.get(f"fusrat_plasma_dd_triton_profile{i}", scan=scan) + for i in range(n_plasma_profile_elements) + ] + + fusrat_plasma_dd_helion_profile = [ + mfile.get(f"fusrat_plasma_dd_helion_profile{i}", scan=scan) + for i in range(n_plasma_profile_elements) + ] + fusrat_plasma_dhe3_profile = [ + mfile.get(f"fusrat_plasma_dhe3_profile{i}", scan=scan) + for i in range(n_plasma_profile_elements) + ] + + fusden_plasma_total_profile = [ + fusrat_plasma_dt_profile[i] + + fusrat_plasma_dd_triton_profile[i] + + fusrat_plasma_dd_helion_profile[i] + + fusrat_plasma_dhe3_profile[i] + for i in range(len(fusrat_plasma_dt_profile)) + ] + + # Integrand for the volume integral N = vol_plasma * 2 * integral(fusden(rho) * rho, drho) + fusrat_plasma_total_profile = np.array([ + fusden_plasma_total_profile[i] * vol_plasma * 2 * rho[i] + for i in range(len(fusden_plasma_total_profile)) + ]) + + # Simpson's rule matches the quadrature used internally for fusrat_total, so the + # cumulative integral naturally lands on the reported total at rho=1 + cum_fusrat = cumulative_simpson(fusrat_plasma_total_profile, x=rho, initial=0.0) + print(cum_fusrat[-1]) + axis.spines["left"].set_color("red") + axis.yaxis.label.set_color("black") + axis.tick_params(axis="y", colors="red") + + # Plot fusion rates (dashed lines, left axis) with axis color and different linestyles + # axis.plot( + # np.linspace(0, 1, len(fusrat_plasma_dt_profile)), + # fusrat_plasma_dt_profile, + # color=axis.spines["left"].get_edgecolor(), + # linestyle="-", + # label=r"$\mathrm{D-T}$", + # ) + # axis.plot( + # np.linspace(0, 1, len(fusrat_plasma_dd_triton_profile)), + # fusrat_plasma_dd_triton_profile, + # color=axis.spines["left"].get_edgecolor(), + # linestyle=":", + # label=r"$\mathrm{D-D \ Triton}$", + # ) + # axis.plot( + # np.linspace(0, 1, len(fusrat_plasma_dd_helion_profile)), + # fusrat_plasma_dd_helion_profile, + # color=axis.spines["left"].get_edgecolor(), + # linestyle="-.", + # label=r"$\mathrm{D-D \ Helion}$", + # ) + # axis.plot( + # np.linspace(0, 1, len(fusrat_plasma_dhe3_profile)), + # fusrat_plasma_dhe3_profile, + # color=axis.spines["left"].get_edgecolor(), + # linestyle="--", + # label=r"$\mathrm{D-3He}$", + # ) + axis.plot( + np.linspace(0, 1, len(fusrat_plasma_total_profile)), + cum_fusrat, + color=axis.spines["left"].get_edgecolor(), + linestyle="None", + marker="d", + markersize=1, + label=r"Total", + ) + + # # Plot fusion power (solid lines, right axis) with axis color and different linestyles + # ax2 = axis.twinx() + # ax2.spines["right"].set_color("blue") + # ax2.yaxis.label.set_color("black") + # ax2.tick_params(axis="y", colors="blue") + # ax2.plot( + # np.linspace(0, 1, len(fusrat_plasma_dt_profile)), + # np.array(fusrat_plasma_dt_profile) * constants.D_T_ENERGY, + # color=ax2.spines["right"].get_edgecolor(), + # linestyle="-", + # ) + + # ax2.plot( + # np.linspace(0, 1, len(fusrat_plasma_dd_triton_profile)), + # np.array(fusrat_plasma_dd_triton_profile) * constants.DD_TRITON_ENERGY, + # color=ax2.spines["right"].get_edgecolor(), + # linestyle=":", + # ) + # ax2.plot( + # np.linspace(0, 1, len(fusrat_plasma_dd_helion_profile)), + # np.array(fusrat_plasma_dd_helion_profile) * constants.DD_HELIUM_ENERGY, + # color=ax2.spines["right"].get_edgecolor(), + # linestyle="-.", + # ) + # ax2.plot( + # np.linspace(0, 1, len(fusrat_plasma_dhe3_profile)), + # np.array(fusrat_plasma_dhe3_profile) * constants.D_HELIUM_ENERGY, + # color=ax2.spines["right"].get_edgecolor(), + # linestyle="--", + # ) + # ax2.plot( + # np.linspace(0, 1, len(fusrat_plasma_total_profile)), + # ( + # np.array(fusrat_plasma_dhe3_profile) * constants.D_HELIUM_ENERGY + # + np.array(fusrat_plasma_dd_helion_profile) * constants.DD_HELIUM_ENERGY + # + np.array(fusrat_plasma_dd_triton_profile) * constants.DD_TRITON_ENERGY + # + np.array(fusrat_plasma_dt_profile) * constants.D_T_ENERGY + # ), + # color=ax2.spines["right"].get_edgecolor(), + # linestyle="None", + # marker="d", + # markersize=1, + # label=r"Total", + # ) + + # ================================================= + + + + # ================================================= + + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.set_ylabel("Fusion Rate [reactions/second]") + axis.legend( + loc="lower left", + edgecolor="black", + facecolor="white", + labelcolor="black", + framealpha=1.0, + frameon=True, + ) + #axis.set_yscale("log") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.set_xlim([0, 1.025]) + axis.minorticks_on() + #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) + # ) + # axis.tick_params(axis="y", which="minor", colors="red") + + # ax2.set_title("Fusion Rate and Fusion Power Profiles") + # ax2.set_ylabel("Fusion Power [W]") + # ax2.set_yscale("log") + # ax2.minorticks_on() + # ax2.yaxis.set_major_locator(plt.LogLocator(base=10.0, numticks=10)) + # ax2.yaxis.set_minor_locator( + # plt.LogLocator(base=10.0, subs=np.arange(1, 10) * 0.1, numticks=100) + # ) + # ax2.tick_params(axis="y", which="minor", colors="blue") def plot_cover_page( axis: plt.Axes, From b1cd9caf8fdda2b6aba4ceff57d21d4b1957c20c Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 1 Oct 2026 11:19:27 +0100 Subject: [PATCH 2/8] Refactor fusion rate profile plotting to use density profiles and improve integration calculations --- process/core/io/plot/summary.py | 238 +++++++++++++------------------- 1 file changed, 94 insertions(+), 144 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index b97d2a1596..b3b45a3610 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -16,7 +16,6 @@ from matplotlib.axes import Axes from matplotlib.patches import Circle, FancyBboxPatch, Rectangle from matplotlib.path import Path as mplPath -from scipy.integrate import cumulative_simpson from scipy.interpolate import interp1d from process.core import constants @@ -11581,7 +11580,7 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i for i in range(n_plasma_profile_elements) ] - fusrat_plasma_total_profile = [ + fusden_plasma_total_profile = [ fusden_plasma_dt_profile[i] + fusden_plasma_dd_triton_profile[i] + fusden_plasma_dd_helion_profile[i] @@ -11623,8 +11622,8 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i label=r"$\mathrm{D-3He}$", ) axis.plot( - np.linspace(0, 1, len(fusrat_plasma_total_profile)), - fusrat_plasma_total_profile, + np.linspace(0, 1, len(fusden_plasma_total_profile)), + fusden_plasma_total_profile, color=axis.spines["left"].get_edgecolor(), linestyle="None", marker="d", @@ -11634,8 +11633,8 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i # Show the plasma volume-averaged rate density and its position on the # profile. - profile_positions = np.linspace(0, 1, len(fusrat_plasma_total_profile)) - profile_rates = np.asarray(fusrat_plasma_total_profile) + profile_positions = np.linspace(0, 1, len(fusden_plasma_total_profile)) + profile_rates = np.asarray(fusden_plasma_total_profile) average_rate = mfile.get("fusden_plasma_vol_avg", scan=scan) axis.axhline( average_rate, @@ -11686,7 +11685,7 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i linestyle="--", ) ax2.plot( - np.linspace(0, 1, len(fusrat_plasma_total_profile)), + np.linspace(0, 1, len(fusden_plasma_total_profile)), ( np.array(fusden_plasma_dhe3_profile) * constants.D_HELIUM_ENERGY + np.array(fusden_plasma_dd_helion_profile) * constants.DD_HELIUM_ENERGY @@ -11883,148 +11882,107 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i transform=fig.transFigure, ) -def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): + +def plot_fusion_rate_profiles(axis: plt.Axes, mfile: MFile, scan: int): """Plot the fusion rate profiles on the given axis""" - fusrat_plasma_dt_profile = [] - fusrat_plasma_dd_triton_profile = [] - fusrat_plasma_dd_helion_profile = [] - fusrat_plasma_dhe3_profile = [] + 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)) vol_plasma = mfile.get("vol_plasma", scan=scan) - - rho, ne, te = profiles_with_pedestal(mfile, scan) - - - - fusrat_plasma_dt_profile = [ - mfile.get(f"fusrat_plasma_dt_profile{i}", scan=scan) + n_plasma_profile_elements = int(mfile.get("n_plasma_profile_elements", scan=scan)) + fusrat_total = mfile.get("fusrat_total", scan=scan) + + rho, _, _ = profiles_with_pedestal(mfile, scan) + + fusden_plasma_dt_profile = [ + mfile.get(f"fusden_plasma_dt_profile{i}", scan=scan) for i in range(n_plasma_profile_elements) ] - fusrat_plasma_dd_triton_profile = [ - mfile.get(f"fusrat_plasma_dd_triton_profile{i}", scan=scan) + fusden_plasma_dd_triton_profile = [ + mfile.get(f"fusden_plasma_dd_triton_profile{i}", scan=scan) for i in range(n_plasma_profile_elements) ] - fusrat_plasma_dd_helion_profile = [ - mfile.get(f"fusrat_plasma_dd_helion_profile{i}", scan=scan) + fusden_plasma_dd_helion_profile = [ + mfile.get(f"fusden_plasma_dd_helion_profile{i}", scan=scan) for i in range(n_plasma_profile_elements) ] - fusrat_plasma_dhe3_profile = [ - mfile.get(f"fusrat_plasma_dhe3_profile{i}", scan=scan) + fusden_plasma_dhe3_profile = [ + mfile.get(f"fusden_plasma_dhe3_profile{i}", scan=scan) for i in range(n_plasma_profile_elements) ] fusden_plasma_total_profile = [ - fusrat_plasma_dt_profile[i] - + fusrat_plasma_dd_triton_profile[i] - + fusrat_plasma_dd_helion_profile[i] - + fusrat_plasma_dhe3_profile[i] - for i in range(len(fusrat_plasma_dt_profile)) + fusden_plasma_dt_profile[i] + + fusden_plasma_dd_triton_profile[i] + + fusden_plasma_dd_helion_profile[i] + + fusden_plasma_dhe3_profile[i] + for i in range(len(fusden_plasma_dt_profile)) ] - - # Integrand for the volume integral N = vol_plasma * 2 * integral(fusden(rho) * rho, drho) - fusrat_plasma_total_profile = np.array([ - fusden_plasma_total_profile[i] * vol_plasma * 2 * rho[i] - for i in range(len(fusden_plasma_total_profile)) - ]) - # Simpson's rule matches the quadrature used internally for fusrat_total, so the - # cumulative integral naturally lands on the reported total at rho=1 - cum_fusrat = cumulative_simpson(fusrat_plasma_total_profile, x=rho, initial=0.0) - print(cum_fusrat[-1]) + # Convert each reaction-rate density into shell contributions, then integrate + # cumulatively to obtain the fusion-rate profile for that reaction. + reaction_profiles = { + r"D-T": fusden_plasma_dt_profile, + r"D-D Triton": fusden_plasma_dd_triton_profile, + r"D-D Helion": fusden_plasma_dd_helion_profile, + r"D-$^3$He": fusden_plasma_dhe3_profile, + "Total": fusden_plasma_total_profile, + } + cum_fusrat_profiles = { + label: calculate_profile_shell_contributions( + profile_x=rho, + profile_y=np.array(profile), + vol_plasma=vol_plasma, + profile_dx=rho[1] - rho[0], + ) + for label, profile in reaction_profiles.items() + } + cum_fusrat_profiles["Total cumulative"] = np.cumsum(cum_fusrat_profiles["Total"]) axis.spines["left"].set_color("red") axis.yaxis.label.set_color("black") axis.tick_params(axis="y", colors="red") - # Plot fusion rates (dashed lines, left axis) with axis color and different linestyles - # axis.plot( - # np.linspace(0, 1, len(fusrat_plasma_dt_profile)), - # fusrat_plasma_dt_profile, - # color=axis.spines["left"].get_edgecolor(), - # linestyle="-", - # label=r"$\mathrm{D-T}$", - # ) - # axis.plot( - # np.linspace(0, 1, len(fusrat_plasma_dd_triton_profile)), - # fusrat_plasma_dd_triton_profile, - # color=axis.spines["left"].get_edgecolor(), - # linestyle=":", - # label=r"$\mathrm{D-D \ Triton}$", - # ) - # axis.plot( - # np.linspace(0, 1, len(fusrat_plasma_dd_helion_profile)), - # fusrat_plasma_dd_helion_profile, - # color=axis.spines["left"].get_edgecolor(), - # linestyle="-.", - # label=r"$\mathrm{D-D \ Helion}$", - # ) - # axis.plot( - # np.linspace(0, 1, len(fusrat_plasma_dhe3_profile)), - # fusrat_plasma_dhe3_profile, - # color=axis.spines["left"].get_edgecolor(), - # linestyle="--", - # label=r"$\mathrm{D-3He}$", - # ) - axis.plot( - np.linspace(0, 1, len(fusrat_plasma_total_profile)), - cum_fusrat, - color=axis.spines["left"].get_edgecolor(), - linestyle="None", - marker="d", - markersize=1, - label=r"Total", - ) + line_colors = { + "D-T": "red", + "D-D Triton": "tab:blue", + "D-D Helion": "tab:green", + "D-$^3$He": "tab:orange", + "Total": "tab:purple", + "Total cumulative": "black", + } + for label, cum_fusrat in cum_fusrat_profiles.items(): + axis.plot( + rho, + cum_fusrat, + color=line_colors[label], + linestyle="-", + label=label, + ) - # # Plot fusion power (solid lines, right axis) with axis color and different linestyles - # ax2 = axis.twinx() - # ax2.spines["right"].set_color("blue") - # ax2.yaxis.label.set_color("black") - # ax2.tick_params(axis="y", colors="blue") - # ax2.plot( - # np.linspace(0, 1, len(fusrat_plasma_dt_profile)), - # np.array(fusrat_plasma_dt_profile) * constants.D_T_ENERGY, - # color=ax2.spines["right"].get_edgecolor(), - # linestyle="-", - # ) - - # ax2.plot( - # np.linspace(0, 1, len(fusrat_plasma_dd_triton_profile)), - # np.array(fusrat_plasma_dd_triton_profile) * constants.DD_TRITON_ENERGY, - # color=ax2.spines["right"].get_edgecolor(), - # linestyle=":", - # ) - # ax2.plot( - # np.linspace(0, 1, len(fusrat_plasma_dd_helion_profile)), - # np.array(fusrat_plasma_dd_helion_profile) * constants.DD_HELIUM_ENERGY, - # color=ax2.spines["right"].get_edgecolor(), - # linestyle="-.", - # ) - # ax2.plot( - # np.linspace(0, 1, len(fusrat_plasma_dhe3_profile)), - # np.array(fusrat_plasma_dhe3_profile) * constants.D_HELIUM_ENERGY, - # color=ax2.spines["right"].get_edgecolor(), - # linestyle="--", - # ) - # ax2.plot( - # np.linspace(0, 1, len(fusrat_plasma_total_profile)), - # ( - # np.array(fusrat_plasma_dhe3_profile) * constants.D_HELIUM_ENERGY - # + np.array(fusrat_plasma_dd_helion_profile) * constants.DD_HELIUM_ENERGY - # + np.array(fusrat_plasma_dd_triton_profile) * constants.DD_TRITON_ENERGY - # + np.array(fusrat_plasma_dt_profile) * constants.D_T_ENERGY - # ), - # color=ax2.spines["right"].get_edgecolor(), - # linestyle="None", - # marker="d", - # markersize=1, - # label=r"Total", - # ) + axis.axhline( + y=fusrat_total, + color="tab:brown", + linestyle="-", + label="Total Fusion Rate", + ) - # ================================================= + half_total_rate = fusrat_total / 2 + half_total_rate_position = np.interp( + half_total_rate, + cum_fusrat_profiles["Total cumulative"], + rho, + ) - + axis.axvline( + x=half_total_rate_position, + color="tab:pink", + linestyle="-", + label="Half Total Fusion Rate Position", + ) # ================================================= @@ -12038,26 +11996,13 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): framealpha=1.0, frameon=True, ) - #axis.set_yscale("log") + axis.set_yscale("log") axis.grid(True, which="both", linestyle="--", alpha=0.5) axis.set_xlim([0, 1.025]) + axis.set_ylim(bottom=1e12) + axis.yaxis.minorticks_on() axis.minorticks_on() - #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) - # ) - # axis.tick_params(axis="y", which="minor", colors="red") - - # ax2.set_title("Fusion Rate and Fusion Power Profiles") - # ax2.set_ylabel("Fusion Power [W]") - # ax2.set_yscale("log") - # ax2.minorticks_on() - # ax2.yaxis.set_major_locator(plt.LogLocator(base=10.0, numticks=10)) - # ax2.yaxis.set_minor_locator( - # plt.LogLocator(base=10.0, subs=np.arange(1, 10) * 0.1, numticks=100) - # ) - # ax2.tick_params(axis="y", which="minor", colors="blue") + def plot_cover_page( axis: plt.Axes, @@ -16513,10 +16458,15 @@ def _add_page(name: str | None = None): _add_page("line_brem_power").add_subplot(121), m_file, scan, imp ) - plot_fusion_rate_profiles( - _add_page("fusion_rate").add_subplot(122), pages["fusion_rate"], m_file, scan + plot_fusion_rate_density_profiles( + _add_page("fusion_rate_density").add_subplot(122), + pages["fusion_rate_density"], + m_file, + scan, ) + plot_fusion_rate_profiles(_add_page("fusion_rate").add_subplot(121), m_file, scan) + _add_page("rx_1_2"), _add_page("rx_3_4") if m_file.get("i_plasma_shape", scan=scan) == PlasmaShapeModelType.SAUTER: plot_fusion_rate_contours(pages["rx_1_2"], pages["rx_3_4"], m_file, scan) From c5bcd009139db7f03c140fa1484be18a5fdd00d7 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 5 Oct 2026 14:09:18 +0100 Subject: [PATCH 3/8] Enhance fusion rate profile plotting with zoom functionality and improved labeling --- process/core/io/plot/summary.py | 97 +++++++++++++++++++++++++-------- 1 file changed, 74 insertions(+), 23 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index b3b45a3610..6dcb49d467 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -4474,7 +4474,7 @@ def plot_line_brem_power_profile( # Ranges # --- - axis.set_xlim(0, 1.0) + axis.set_xlim(0, 1.01) 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) @@ -11883,8 +11883,23 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i ) -def plot_fusion_rate_profiles(axis: plt.Axes, mfile: MFile, scan: int): - """Plot the fusion rate profiles on the given axis""" +def plot_fusion_rate_profiles( + axis: plt.Axes, mfile: MFile, scan: int, zoom_axis: plt.Axes | None = None +): + """Plot the fusion rate profiles on the given axis + + Parameters + ---------- + axis : plt.Axes + The axis for the full (log scale) fusion rate profiles. + mfile : MFile + The MFILE data object. + scan : int + The scan number to use for extracting data. + zoom_axis : plt.Axes | None + If given, an axis on which only the cumulative total fusion rate is plotted + on a linear scale. + """ fusden_plasma_dt_profile = [] fusden_plasma_dd_triton_profile = [] fusden_plasma_dd_helion_profile = [] @@ -11926,11 +11941,11 @@ def plot_fusion_rate_profiles(axis: plt.Axes, mfile: MFile, scan: int): # Convert each reaction-rate density into shell contributions, then integrate # cumulatively to obtain the fusion-rate profile for that reaction. reaction_profiles = { - r"D-T": fusden_plasma_dt_profile, - r"D-D Triton": fusden_plasma_dd_triton_profile, - r"D-D Helion": fusden_plasma_dd_helion_profile, - r"D-$^3$He": fusden_plasma_dhe3_profile, - "Total": fusden_plasma_total_profile, + "$\\Gamma_{D-T}$": fusden_plasma_dt_profile, + "$\\Gamma_{\\mathrm{D-D\\ Triton}}$": fusden_plasma_dd_triton_profile, + "$\\Gamma_{\\mathrm{D-D\\ Helion}}$": fusden_plasma_dd_helion_profile, + "$\\Gamma_{\\mathrm{D-}^3\\mathrm{He}}$": fusden_plasma_dhe3_profile, + "$\\Gamma_{\\mathrm{Total}}$": fusden_plasma_total_profile, } cum_fusrat_profiles = { label: calculate_profile_shell_contributions( @@ -11941,39 +11956,39 @@ def plot_fusion_rate_profiles(axis: plt.Axes, mfile: MFile, scan: int): ) for label, profile in reaction_profiles.items() } - cum_fusrat_profiles["Total cumulative"] = np.cumsum(cum_fusrat_profiles["Total"]) - axis.spines["left"].set_color("red") - axis.yaxis.label.set_color("black") - axis.tick_params(axis="y", colors="red") + cum_fusrat_profiles["$\\Sigma\\ \\Gamma_{\\mathrm{Total}}$"] = np.cumsum( + cum_fusrat_profiles["$\\Gamma_{\\mathrm{Total}}$"] + ) line_colors = { - "D-T": "red", - "D-D Triton": "tab:blue", - "D-D Helion": "tab:green", - "D-$^3$He": "tab:orange", - "Total": "tab:purple", - "Total cumulative": "black", + "$\\Gamma_{D-T}$": "red", + "$\\Gamma_{\\mathrm{D-D\\ Triton}}$": "tab:blue", + "$\\Gamma_{\\mathrm{D-D\\ Helion}}$": "tab:green", + "$\\Gamma_{\\mathrm{D-}^3\\mathrm{He}}$": "tab:orange", + "$\\Gamma_{\\mathrm{Total}}$": "tab:purple", + "$\\Sigma\\ \\Gamma_{\\mathrm{Total}}$": "black", } for label, cum_fusrat in cum_fusrat_profiles.items(): axis.plot( rho, cum_fusrat, color=line_colors[label], - linestyle="-", + linestyle=(":" if label == "$\\Gamma_{\\mathrm{Total}}$" else "-"), label=label, + linewidth=2 if label == "$\\Gamma_{\\mathrm{Total}}$" else 1, ) axis.axhline( y=fusrat_total, color="tab:brown", - linestyle="-", + linestyle=":", label="Total Fusion Rate", ) half_total_rate = fusrat_total / 2 half_total_rate_position = np.interp( half_total_rate, - cum_fusrat_profiles["Total cumulative"], + cum_fusrat_profiles["$\\Sigma\\ \\Gamma_{\\mathrm{Total}}$"], rho, ) @@ -11981,7 +11996,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, mfile: MFile, scan: int): x=half_total_rate_position, color="tab:pink", linestyle="-", - label="Half Total Fusion Rate Position", + label="50% Total Fusion Rate", ) # ================================================= @@ -12002,6 +12017,36 @@ def plot_fusion_rate_profiles(axis: plt.Axes, mfile: MFile, scan: int): axis.set_ylim(bottom=1e12) axis.yaxis.minorticks_on() axis.minorticks_on() + axis.set_title("Thermal Fusion Rate Profiles (log scale)") + + if zoom_axis is not None: + zoom_axis.plot( + rho, + cum_fusrat_profiles["$\\Gamma_{\\mathrm{Total}}$"], + color="tab:purple", + linestyle=":", + label="$\\Gamma_{\\mathrm{Total}}$", + ) + zoom_axis.axvline( + x=half_total_rate_position, + color="tab:pink", + linestyle="-", + label="50% Total Fusion Rate", + ) + zoom_axis.set_xlabel("$\\rho \\ [r/a]$") + zoom_axis.set_ylabel("Fusion Rate [reactions/second]") + zoom_axis.set_title("Thermal Fusion Rate Profiles (linear scale)") + zoom_axis.set_xlim([0, 1.01]) + zoom_axis.set_ylim(bottom=0.0) + zoom_axis.legend( + loc="lower right", + edgecolor="black", + facecolor="white", + framealpha=1.0, + frameon=True, + ) + zoom_axis.grid(True, which="both", linestyle="--", alpha=0.5) + zoom_axis.minorticks_on() def plot_cover_page( @@ -16465,7 +16510,13 @@ def _add_page(name: str | None = None): scan, ) - plot_fusion_rate_profiles(_add_page("fusion_rate").add_subplot(121), m_file, scan) + fusion_rate_page = _add_page("fusion_rate") + plot_fusion_rate_profiles( + fusion_rate_page.add_subplot(121), + m_file, + scan, + zoom_axis=fusion_rate_page.add_subplot(122), + ) _add_page("rx_1_2"), _add_page("rx_3_4") if m_file.get("i_plasma_shape", scan=scan) == PlasmaShapeModelType.SAUTER: From 5255aef5eca3ba99623141faeca45861294a547e Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 5 Oct 2026 14:23:02 +0100 Subject: [PATCH 4/8] Update fusion rate density profile labels and axis limits for improved clarity --- process/core/io/plot/summary.py | 46 ++++++++++++++++----------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index 6dcb49d467..d9c403f615 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -4312,7 +4312,7 @@ 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.01) axis.set_yscale("log") axis.yaxis.grid(True, which="both", alpha=0.2) # DEMO : Fixed ranges for comparison @@ -11737,8 +11737,8 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i # Add plasma volume, areas and shaping information textstr_general = ( f"Total fusion rate: {mfile.get('fusrat_total', scan=scan):.4e} reactions/s\n" - f"Total volume averaged fusion rate density: {mfile.get('fusden_total_vol_avg', scan=scan):.4e} reactions/m3/s\n" - f"Plasma volume averaged fusion rate density: {mfile.get('fusden_plasma_vol_avg', scan=scan):.4e} reactions/m3/s\n" + f"Total volume averaged fusion rate density: {mfile.get('fusden_total_vol_avg', scan=scan):.4e} reactions/m³/s\n" + f"Plasma volume averaged fusion rate density: {mfile.get('fusden_plasma_vol_avg', scan=scan):.4e} reactions/m³/s" ) axis.text( @@ -11754,8 +11754,8 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i 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"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" + 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" ) axis.text( @@ -11767,8 +11767,8 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i ) axis.text( - 0.24, - 0.8, + 0.39, + 0.79, "$\\text{D - T}$", fontsize=20, verticalalignment="top", @@ -11779,8 +11779,8 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i 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"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}" ) axis.text( @@ -11792,7 +11792,7 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i ) axis.text( - 0.22, + 0.36, 0.685, "$\\text{D - D}$", fontsize=20, @@ -11804,7 +11804,7 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i textstr_dhe3 = ( 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" + f"Volume-averaged total power density: {mfile.get('pden_dhe3_total_vol_avg_mw', scan=scan):,.3e} MW/m³ \n" ) axis.text( @@ -11816,7 +11816,7 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i ) axis.text( - 0.21, + 0.36, 0.59, "$\\text{D - 3He}$", fontsize=20, @@ -11830,12 +11830,12 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i f"Total power: {mfile.get('p_alpha_total_mw', scan=scan):.2f} MW\n" f"Plasma power: {mfile.get('p_plasma_alpha_mw', scan=scan):.2f} MW\n" f"Beam power: {mfile.get('p_beam_alpha_mw', scan=scan):.2f} MW\n\n" - f"Volume-averaged rate density total: {mfile.get('fusden_alpha_total_vol_avg', scan=scan):.4e} particles/m3/sec\n" - f"Volume-averaged rate density, plasma: {mfile.get('fusden_plasma_alpha_vol_avg', scan=scan):.4e} particles/m3/sec\n\n" - f"Volume-averaged total power density: {mfile.get('pden_alpha_total_vol_avg_mw', scan=scan):.4e} MW/m3\n" - f"Volume-averaged plasma power density: {mfile.get('pden_plasma_alpha_vol_avg_mw', scan=scan):.4e} MW/m3\n\n" - f"Power per unit volume transferred to electrons: {mfile.get('f_pden_alpha_electron_mw', scan=scan):.4e} MW/m3\n" - f"Power per unit volume transferred to ions: {mfile.get('f_pden_alpha_ions_mw', scan=scan):.4e} MW/m3\n\n" + f"Volume-averaged rate density total: {mfile.get('fusden_alpha_total_vol_avg', scan=scan):.4e} particles/m³/sec\n" + f"Volume-averaged rate density, plasma: {mfile.get('fusden_plasma_alpha_vol_avg', scan=scan):.4e} particles/m³/sec\n\n" + f"Volume-averaged total power density: {mfile.get('pden_alpha_total_vol_avg_mw', scan=scan):.4e} MW/m³\n" + f"Volume-averaged plasma power density: {mfile.get('pden_plasma_alpha_vol_avg_mw', scan=scan):.4e} MW/m³\n\n" + f"Power per unit volume transferred to electrons: {mfile.get('f_pden_alpha_electron_mw', scan=scan):.4e} MW/m³\n" + f"Power per unit volume transferred to ions: {mfile.get('f_pden_alpha_ions_mw', scan=scan):.4e} MW/m³" ) axis.text( @@ -11848,7 +11848,7 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i axis.text( 0.35, - 0.45, + 0.425, "$\\alpha$", fontsize=22, verticalalignment="top", @@ -11861,8 +11861,8 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i f"Total power: {mfile.get('p_neutron_total_mw', scan=scan):,.2f} MW\n" f"Plasma power: {mfile.get('p_plasma_neutron_mw', scan=scan):,.2f} MW\n" f"Beam power: {mfile.get('p_beam_neutron_mw', scan=scan):,.2f} MW\n\n" - f"Volume-averaged total power density: {mfile.get('pden_neutron_total_vol_avg_mw', scan=scan):,.4e} MW/m3\n" - f"Volume-averaged plasma power density: {mfile.get('pden_plasma_neutron_vol_avg_mw', scan=scan):,.4e} MW/m3\n" + f"Volume-averaged total power density: {mfile.get('pden_neutron_total_vol_avg_mw', scan=scan):,.4e} MW/m³\n" + f"Volume-averaged plasma power density: {mfile.get('pden_plasma_neutron_vol_avg_mw', scan=scan):,.4e} MW/m³" ) axis.text( @@ -11874,8 +11874,8 @@ def plot_fusion_rate_density_profiles(axis: plt.Axes, fig, mfile: MFile, scan: i ) axis.text( - 0.25, - 0.2, + 0.3, + 0.175, "$n$", fontsize=20, verticalalignment="top", From 45882140d923af2acf1333a874a773d02693aa4b Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 12 Oct 2026 15:48:19 +0100 Subject: [PATCH 5/8] Update legend position in fusion rate profile plot for better visibility --- process/core/io/plot/summary.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index d9c403f615..e888941bfb 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -12039,7 +12039,7 @@ def plot_fusion_rate_profiles( zoom_axis.set_xlim([0, 1.01]) zoom_axis.set_ylim(bottom=0.0) zoom_axis.legend( - loc="lower right", + loc="upper right", edgecolor="black", facecolor="white", framealpha=1.0, From a58947c0e0d788629a29c5999184c75ef01f4c70 Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Mon, 12 Oct 2026 16:09:07 +0100 Subject: [PATCH 6/8] fix fusrat function --- process/core/io/plot/summary.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index e888941bfb..3f7b59a787 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -12958,15 +12958,15 @@ def plot_fusion_rate_contours(fig1, fig2, mfile: MFile, scan: int): def fusrat(name): fusrat_dat = [ - mfile.get(f"fusrat_plasma_{name}_profile{i}", scan=scan) + mfile.get(f"fusden_plasma_{name}_profile{i}", scan=scan) for i in range(n_plasma_profile_elements) ] return interp1d_profile(fusrat_dat, 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") + 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") From ec9c630f7f0b663df81c127275dc438668e9b6bd Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 12 Oct 2026 16:17:09 +0100 Subject: [PATCH 7/8] Conditionalize EBW ECRH coupling graph plotting based on current drive model --- process/core/io/plot/summary.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index 3f7b59a787..9fedfb1c56 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -16585,9 +16585,11 @@ def _add_page(name: str | None = None): plot_cumulative_plasma_thermal_energy_profiles( ax_thermal_energy_cumulative, m_file, scan ) - - plot_ebw_ecrh_coupling_graph(_add_page().add_subplot(111), m_file, scan) - + if CurrentDriveModel(m_file.get("i_hcd_primary", scan=scan)) in { + CurrentDriveModel.USER_INPUT_ELECTRON_BERNSTEIN, + CurrentDriveModel.FREETHY_ELECTRON_CYCLOTRON, + }: + plot_ebw_ecrh_coupling_graph(_add_page().add_subplot(111), m_file, scan) plot_bootstrap_comparison(_add_page("current").add_subplot(221), m_file, scan) plot_plasma_current_comparison(pages["current"].add_subplot(224), m_file, scan) plot_h_threshold_comparison( From 75631879f127729915bbf379d4127d0befe7ddd3 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 12 Oct 2026 16:21:51 +0100 Subject: [PATCH 8/8] Add filled contour plot for plasma poloidal pressure in summary plot --- process/core/io/plot/summary.py | 1 + 1 file changed, 1 insertion(+) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index 9fedfb1c56..3ad0496bd8 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -12569,6 +12569,7 @@ 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 + 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)