From d63fa511d2f1fe94d99bac9c349c69aee405c46c Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 28 Aug 2025 11:04:56 +0100 Subject: [PATCH 01/12] :sparkle: Add fusion reaction rate profile variable for D-T plasma, fusrat_plasma_dt_profile --- process/data_structure/physics_variables.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 46d2571608..191ff128ba 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -511,6 +511,8 @@ fusrat_total: float = None """fusion reaction rate, from beams and plasma (reactions/sec)""" +fusrat_plasma_dt_profile: list[float] = None +"""Profile of D-T fusion reaction rate in plasma, (reactions/sec)""" fusden_plasma: float = None """fusion reaction rate, just from plasma (reactions/m3/sec)""" @@ -1420,6 +1422,7 @@ def init_physics_variables(): global f_tritium global fusden_total global fusrat_total + global fusrat_plasma_dt_profile global fusden_plasma global f_c_plasma_non_inductive global ejima_coeff @@ -1664,6 +1667,7 @@ def init_physics_variables(): f_tritium = 0.5 fusden_total = 0.0 fusrat_total = 0.0 + fusrat_plasma_dt_profile = [] fusden_plasma = 0.0 f_c_plasma_non_inductive = 1.0 ejima_coeff = 0.4 From a7a74e73c7311ad9d3032c373238e91461e7bcea Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 28 Aug 2025 14:13:38 +0100 Subject: [PATCH 02/12] Add fusion rate profile plotting and data handling for D-T plasma --- process/fusion_reactions.py | 5 +++++ process/io/plot_proc.py | 25 +++++++++++++++++++++++++ process/physics.py | 8 ++++++++ 3 files changed, 38 insertions(+) diff --git a/process/fusion_reactions.py b/process/fusion_reactions.py index 88720b0b9c..1add9e8d14 100644 --- a/process/fusion_reactions.py +++ b/process/fusion_reactions.py @@ -216,6 +216,11 @@ def dt_reaction(self) -> None: x=self.plasma_profile.neprofile.profile_x, dx=self.plasma_profile.neprofile.profile_dx, ) + physics_variables.fusrat_plasma_dt_profile = ( + fusion_rate_integral(self.plasma_profile, dt) + * (physics_variables.f_deuterium * physics_variables.nd_fuel_ions) + * (physics_variables.f_tritium * physics_variables.nd_fuel_ions) + ) # Store the average fusion reaction rate self.sigmav_dt_average = sigmav diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 07d1072626..258a37cb11 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -9835,6 +9835,25 @@ def plot_fw_90_deg_pipe_bend(ax, m_file_data, scan): ) +def plot_fusion_rate_profiles(axis, mfile_data, scan): + # Plot the fusion rate profiles on the given axis + fusrat_plasma_dt_profile = [] + fusrat_plasma_dt_profile = [ + mfile_data.data[f"fusrat_plasma_dt_profile{i}"].get_scan(scan) + for i in range(500) + ] + # Change from 0 to 1 index to align with poloidal cross-section plot numbering + axis.plot( + np.linspace(0, 1, len(fusrat_plasma_dt_profile)), + fusrat_plasma_dt_profile, + linestyle="--", + ) + + axis.set_xlabel("$\\rho$") + axis.set_ylabel("Fusion Rate [MW/m^3]") + axis.legend() + + def main_plot( fig1, fig2, @@ -10035,6 +10054,12 @@ def main_plot( plot_33 = fig17.add_subplot(111, aspect="equal") plot_main_power_flow(plot_33, m_file_data, scan, fig17) + plot_33 = fig17.add_subplot(122) + plot_fusion_rate_profiles(plot_33, m_file_data, scan) + + plot_33 = fig17.add_subplot(122) + plot_fusion_rate_profiles(plot_33, m_file_data, scan) + def main(args=None): # TODO The use of globals here isn't ideal, but is required to get main() diff --git a/process/physics.py b/process/physics.py index 5cc4b403b8..c81bd49ed1 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1528,6 +1528,7 @@ def _trapped_particle_fraction_sauter( class Physics: def __init__(self, plasma_profile, current_drive): self.outfile = constants.NOUT + self.mfile = constants.mfile self.plasma_profile = plasma_profile self.current_drive = current_drive @@ -4958,6 +4959,13 @@ def outplas(self): physics_variables.p_dt_total_mw, "OP ", ) + for i in range(len(physics_variables.fusrat_plasma_dt_profile)): + po.ovarre( + self.mfile, + f"DT fusion rate at point {i}", + f"fusrat_plasma_dt_profile{i}", + physics_variables.fusrat_plasma_dt_profile[i], + ) po.ovarre( self.outfile, "D-T fusion power: plasma (MW)", From 837efbb7fd543786ead9a859948496142d9eca40 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 28 Aug 2025 14:18:55 +0100 Subject: [PATCH 03/12] :sparkle: Add fusion reaction rate profile variable for D-T plasma, fusrat_plasma_dd_triton_profile --- process/data_structure/physics_variables.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 191ff128ba..c28b6596b2 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -514,6 +514,9 @@ fusrat_plasma_dt_profile: list[float] = None """Profile of D-T fusion reaction rate in plasma, (reactions/sec)""" +fusrat_plasma_dd_triton_profile: list[float] = None +"""Profile of D-D fusion reaction rate (triton branch) in plasma, (reactions/sec)""" + fusden_plasma: float = None """fusion reaction rate, just from plasma (reactions/m3/sec)""" @@ -1423,6 +1426,7 @@ def init_physics_variables(): global fusden_total global fusrat_total global fusrat_plasma_dt_profile + global fusrat_plasma_dd_triton_profile global fusden_plasma global f_c_plasma_non_inductive global ejima_coeff @@ -1668,6 +1672,7 @@ def init_physics_variables(): fusden_total = 0.0 fusrat_total = 0.0 fusrat_plasma_dt_profile = [] + fusrat_plasma_dd_triton_profile = [] fusden_plasma = 0.0 f_c_plasma_non_inductive = 1.0 ejima_coeff = 0.4 From 2479a7ced6da2c2d03dc87336e3fff709174a26b Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 28 Aug 2025 14:19:43 +0100 Subject: [PATCH 04/12] :sparkle: Add fusion reaction rate profile variable for D-D plasma, fusrat_plasma_dd_helion_profile --- process/data_structure/physics_variables.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index c28b6596b2..f468793fad 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -514,6 +514,9 @@ fusrat_plasma_dt_profile: list[float] = None """Profile of D-T fusion reaction rate in plasma, (reactions/sec)""" +fusrat_plasma_dd_helion_profile: list[float] = None +"""Profile of D-D fusion reaction rate (helium branch) in plasma, (reactions/sec)""" + fusrat_plasma_dd_triton_profile: list[float] = None """Profile of D-D fusion reaction rate (triton branch) in plasma, (reactions/sec)""" @@ -1427,6 +1430,7 @@ def init_physics_variables(): global fusrat_total global fusrat_plasma_dt_profile global fusrat_plasma_dd_triton_profile + global fusrat_plasma_dd_helion_profile global fusden_plasma global f_c_plasma_non_inductive global ejima_coeff @@ -1673,6 +1677,7 @@ def init_physics_variables(): fusrat_total = 0.0 fusrat_plasma_dt_profile = [] fusrat_plasma_dd_triton_profile = [] + fusrat_plasma_dd_helion_profile = [] fusden_plasma = 0.0 f_c_plasma_non_inductive = 1.0 ejima_coeff = 0.4 From de4ef0d9b30023ed06ceab4d31c39601d592d563 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 28 Aug 2025 15:08:35 +0100 Subject: [PATCH 05/12] Add D-D and D-T fusion rate profile calculations and plotting --- process/fusion_reactions.py | 35 ++++++++++++-- process/io/plot_proc.py | 91 +++++++++++++++++++++++++++++++++++-- process/physics.py | 15 ++++++ 3 files changed, 131 insertions(+), 10 deletions(-) diff --git a/process/fusion_reactions.py b/process/fusion_reactions.py index 1add9e8d14..0134207efd 100644 --- a/process/fusion_reactions.py +++ b/process/fusion_reactions.py @@ -210,17 +210,22 @@ def dt_reaction(self) -> None: # Initialize Bosch-Hale constants for the D-T reaction dt = BoschHaleConstants(**REACTION_CONSTANTS_DT) + physics_variables.fusrat_plasma_dt_profile = ( + bosch_hale_reactivity( + (physics_variables.ti / physics_variables.te) + * self.plasma_profile.teprofile.profile_y, + dt, + ) + * (physics_variables.f_deuterium * physics_variables.nd_fuel_ions) + * (physics_variables.f_tritium * physics_variables.nd_fuel_ions) + ) + # Calculate the fusion reaction rate integral using Simpson's rule sigmav = integrate.simpson( fusion_rate_integral(self.plasma_profile, dt), x=self.plasma_profile.neprofile.profile_x, dx=self.plasma_profile.neprofile.profile_dx, ) - physics_variables.fusrat_plasma_dt_profile = ( - fusion_rate_integral(self.plasma_profile, dt) - * (physics_variables.f_deuterium * physics_variables.nd_fuel_ions) - * (physics_variables.f_tritium * physics_variables.nd_fuel_ions) - ) # Store the average fusion reaction rate self.sigmav_dt_average = sigmav @@ -364,6 +369,16 @@ def dd_helion_reaction(self) -> None: dx=self.plasma_profile.neprofile.profile_dx, ) + physics_variables.fusrat_plasma_dd_helion_profile = ( + bosch_hale_reactivity( + (physics_variables.ti / physics_variables.te) + * self.plasma_profile.teprofile.profile_y, + dd1, + ) + * (physics_variables.f_deuterium * physics_variables.nd_fuel_ions) + * (physics_variables.f_deuterium * physics_variables.nd_fuel_ions) + ) + # Reaction energy in MegaJoules [MJ] reaction_energy = constants.DD_HELIUM_ENERGY / 1.0e6 @@ -437,6 +452,16 @@ def dd_triton_reaction(self) -> None: dx=self.plasma_profile.neprofile.profile_dx, ) + physics_variables.fusrat_plasma_dd_triton_profile = ( + bosch_hale_reactivity( + (physics_variables.ti / physics_variables.te) + * self.plasma_profile.teprofile.profile_y, + dd2, + ) + * (physics_variables.f_deuterium * physics_variables.nd_fuel_ions) + * (physics_variables.f_deuterium * physics_variables.nd_fuel_ions) + ) + # Reaction energy in MegaJoules [MJ] reaction_energy = constants.DD_TRITON_ENERGY / 1.0e6 diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 258a37cb11..1d1c4550db 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -32,6 +32,7 @@ import process.io.mfile as mf import process.superconducting_tf_coil as sctf from process.data_structure import physics_variables +from process.fortran import constants from process.geometry.blanket_geometry import ( blanket_geometry_double_null, blanket_geometry_single_null, @@ -9838,20 +9839,100 @@ def plot_fw_90_deg_pipe_bend(ax, m_file_data, scan): def plot_fusion_rate_profiles(axis, mfile_data, scan): # 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_dt_profile = [ mfile_data.data[f"fusrat_plasma_dt_profile{i}"].get_scan(scan) for i in range(500) ] - # Change from 0 to 1 index to align with poloidal cross-section plot numbering + + fusrat_plasma_dd_triton_profile = [ + mfile_data.data[f"fusrat_plasma_dd_triton_profile{i}"].get_scan(scan) + for i in range(500) + ] + + fusrat_plasma_dd_helion_profile = [ + mfile_data.data[f"fusrat_plasma_dd_helion_profile{i}"].get_scan(scan) + for i in range(500) + ] + + 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, - linestyle="--", + 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}$", + ) + + # 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="-.", ) - axis.set_xlabel("$\\rho$") - axis.set_ylabel("Fusion Rate [MW/m^3]") - axis.legend() + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.set_ylabel("Fusion Rate [reactions/second]") + axis.legend( + loc="lower left", edgecolor="black", facecolor="white", labelcolor="black" + ) + axis.set_yscale("log") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.set_xlim([0, 1.01]) + axis.minorticks_on() + axis.set_ylim([1e10, 1e19]) + 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 main_plot( diff --git a/process/physics.py b/process/physics.py index c81bd49ed1..80de017553 100644 --- a/process/physics.py +++ b/process/physics.py @@ -4966,6 +4966,21 @@ def outplas(self): f"fusrat_plasma_dt_profile{i}", physics_variables.fusrat_plasma_dt_profile[i], ) + + for i in range(len(physics_variables.fusrat_plasma_dd_triton_profile)): + po.ovarre( + self.mfile, + f"DT fusion rate at point {i}", + f"fusrat_plasma_dd_triton_profile{i}", + physics_variables.fusrat_plasma_dd_triton_profile[i], + ) + for i in range(len(physics_variables.fusrat_plasma_dd_helion_profile)): + po.ovarre( + self.mfile, + f"DT fusion rate at point {i}", + f"fusrat_plasma_dd_helion_profile{i}", + physics_variables.fusrat_plasma_dd_helion_profile[i], + ) po.ovarre( self.outfile, "D-T fusion power: plasma (MW)", From 710878f85bcbed32620190e85a8a34b5f861e270 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 29 Aug 2025 09:56:26 +0100 Subject: [PATCH 06/12] :sparkle: Add D-3He fusion reaction rate profile variable and update initialization --- process/data_structure/physics_variables.py | 9 +++++++-- process/io/plot_proc.py | 1 + 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index f468793fad..ec826aa88c 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -517,8 +517,11 @@ fusrat_plasma_dd_helion_profile: list[float] = None """Profile of D-D fusion reaction rate (helium branch) in plasma, (reactions/sec)""" -fusrat_plasma_dd_triton_profile: list[float] = None -"""Profile of D-D fusion reaction rate (triton branch) in plasma, (reactions/sec)""" +fusrat_plasma_dd_helion_profile: list[float] = None +"""Profile of D-D fusion reaction rate (helium branch) in plasma, (reactions/sec)""" + +fusrat_plasma_dhe3_profile: list[float] = None +"""Profile of D-3He fusion reaction rate in plasma, (reactions/sec)""" fusden_plasma: float = None """fusion reaction rate, just from plasma (reactions/m3/sec)""" @@ -1431,6 +1434,7 @@ def init_physics_variables(): global fusrat_plasma_dt_profile global fusrat_plasma_dd_triton_profile global fusrat_plasma_dd_helion_profile + global fusrat_plasma_dhe3_profile global fusden_plasma global f_c_plasma_non_inductive global ejima_coeff @@ -1678,6 +1682,7 @@ def init_physics_variables(): fusrat_plasma_dt_profile = [] fusrat_plasma_dd_triton_profile = [] fusrat_plasma_dd_helion_profile = [] + fusrat_plasma_dhe3_profile = [] fusden_plasma = 0.0 f_c_plasma_non_inductive = 1.0 ejima_coeff = 0.4 diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 1d1c4550db..0af2ab4e64 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -9895,6 +9895,7 @@ def plot_fusion_rate_profiles(axis, mfile_data, scan): 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, From 5a7de053bac6fd4e53fe76968875e16c958e10e4 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 29 Aug 2025 10:07:49 +0100 Subject: [PATCH 07/12] :sparkle: Add D-3He fusion reaction rate profile calculations and plotting --- process/fusion_reactions.py | 10 ++++++++++ process/io/plot_proc.py | 20 +++++++++++++++++++- process/physics.py | 11 +++++++++-- 3 files changed, 38 insertions(+), 3 deletions(-) diff --git a/process/fusion_reactions.py b/process/fusion_reactions.py index 0134207efd..6ab8ce1f85 100644 --- a/process/fusion_reactions.py +++ b/process/fusion_reactions.py @@ -298,6 +298,16 @@ def dhe3_reaction(self) -> None: x=self.plasma_profile.neprofile.profile_x, dx=self.plasma_profile.neprofile.profile_dx, ) + + physics_variables.fusrat_plasma_dhe3_profile = ( + bosch_hale_reactivity( + (physics_variables.ti / physics_variables.te) + * self.plasma_profile.teprofile.profile_y, + dhe3, + ) + * (physics_variables.f_deuterium * physics_variables.nd_fuel_ions) + * (physics_variables.f_helium3 * physics_variables.nd_fuel_ions) + ) # Reaction energy in MegaJoules [MJ] reaction_energy = constants.D_HELIUM_ENERGY / 1.0e6 diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 0af2ab4e64..4370c77d2a 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -9841,7 +9841,8 @@ def plot_fusion_rate_profiles(axis, mfile_data, scan): fusrat_plasma_dt_profile = [] fusrat_plasma_dd_triton_profile = [] fusrat_plasma_dd_helion_profile = [] - + fusrat_plasma_dhe3_profile = [] + fusrat_plasma_dt_profile = [ mfile_data.data[f"fusrat_plasma_dt_profile{i}"].get_scan(scan) for i in range(500) @@ -9856,6 +9857,10 @@ def plot_fusion_rate_profiles(axis, mfile_data, scan): mfile_data.data[f"fusrat_plasma_dd_helion_profile{i}"].get_scan(scan) for i in range(500) ] + fusrat_plasma_dhe3_profile = [ + mfile_data.data[f"fusrat_plasma_dhe3_profile{i}"].get_scan(scan) + for i in range(500) + ] axis.spines["left"].set_color("red") axis.yaxis.label.set_color("black") @@ -9883,6 +9888,13 @@ def plot_fusion_rate_profiles(axis, mfile_data, scan): 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}$", + ) # Plot fusion power (solid lines, right axis) with axis color and different linestyles ax2 = axis.twinx() @@ -9908,6 +9920,12 @@ def plot_fusion_rate_profiles(axis, mfile_data, scan): 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="--", + ) axis.set_xlabel("$\\rho \\ [r/a]$") axis.set_ylabel("Fusion Rate [reactions/second]") diff --git a/process/physics.py b/process/physics.py index 80de017553..b62350ced9 100644 --- a/process/physics.py +++ b/process/physics.py @@ -4970,17 +4970,24 @@ def outplas(self): for i in range(len(physics_variables.fusrat_plasma_dd_triton_profile)): po.ovarre( self.mfile, - f"DT fusion rate at point {i}", + f"D-D -> T fusion rate at point {i}", f"fusrat_plasma_dd_triton_profile{i}", physics_variables.fusrat_plasma_dd_triton_profile[i], ) for i in range(len(physics_variables.fusrat_plasma_dd_helion_profile)): po.ovarre( self.mfile, - f"DT fusion rate at point {i}", + f"D-D -> 3He fusion rate at point {i}", f"fusrat_plasma_dd_helion_profile{i}", physics_variables.fusrat_plasma_dd_helion_profile[i], ) + for i in range(len(physics_variables.fusrat_plasma_dhe3_profile)): + po.ovarre( + self.mfile, + f"D-3He fusion rate at point {i}", + f"fusrat_plasma_dhe3_profile{i}", + physics_variables.fusrat_plasma_dhe3_profile[i], + ) po.ovarre( self.outfile, "D-T fusion power: plasma (MW)", From 844a4ec2cf8522e4f1b0b7fed6cc7c3171adf92e Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 8 Sep 2025 11:03:53 +0100 Subject: [PATCH 08/12] :sparkle: Update plotting functions to include fusion rate profiles and adjust subplot configurations --- process/io/plot_proc.py | 73 +++++++++++++++++++++-------------------- 1 file changed, 37 insertions(+), 36 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 4370c77d2a..f670c6a1d8 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -10064,86 +10064,89 @@ def main_plot( plot_13 = fig4.add_subplot(4, 3, 12) plot_13.set_position([0.7, 0.125, 0.25, 0.15]) plot_qprofile(plot_13, demo_ranges, m_file_data, scan) + + plot_14 = fig5.add_subplot(122) + plot_fusion_rate_profiles(plot_14, m_file_data, scan) # Plot poloidal cross-section - plot_14 = fig5.add_subplot(121, aspect="equal") - poloidal_cross_section(plot_14, m_file_data, scan, demo_ranges, colour_scheme) + plot_15 = fig6.add_subplot(121, aspect="equal") + poloidal_cross_section(plot_15, m_file_data, scan, demo_ranges, colour_scheme) # Plot toroidal cross-section - plot_15 = fig5.add_subplot(122, aspect="equal") - toroidal_cross_section(plot_15, m_file_data, scan, demo_ranges, colour_scheme) + plot_16 = fig6.add_subplot(122, aspect="equal") + toroidal_cross_section(plot_16, m_file_data, scan, demo_ranges, colour_scheme) # fig4.subplots_adjust(bottom=-0.2, top = 0.9, left = 0.1, right = 0.9) # Plot color key - plot_16 = fig5.add_subplot(222) - plot_16.set_position([0.5, 0.5, 0.5, 0.5]) - color_key(plot_16, m_file_data, scan, colour_scheme) + plot_17 = fig6.add_subplot(222) + plot_17.set_position([0.5, 0.5, 0.5, 0.5]) + color_key(plot_17, m_file_data, scan, colour_scheme) - plot_17 = fig6.add_subplot(211) - plot_17.set_position([0.1, 0.33, 0.8, 0.6]) # x0, y0, width, height (2/3 vertical) - plot_radial_build(plot_17, m_file_data, colour_scheme) + plot_18 = fig7.add_subplot(211) + plot_18.set_position([0.1, 0.33, 0.8, 0.6]) # x0, y0, width, height (2/3 vertical) + plot_radial_build(plot_18, m_file_data, colour_scheme) # Make each axes smaller vertically to leave room for the legend - plot_175 = fig7.add_subplot(211) - plot_175.set_position([0.1, 0.61, 0.8, 0.32]) # x0, y0, width, height + plot_185 = fig8.add_subplot(211) + plot_185.set_position([0.1, 0.61, 0.8, 0.32]) # x0, y0, width, height - plot_17 = fig7.add_subplot(212) - plot_17.set_position([0.1, 0.13, 0.8, 0.32]) # x0, y0, width, height - plot_upper_vertical_build(plot_175, m_file_data, colour_scheme) - plot_lower_vertical_build(plot_17, m_file_data, colour_scheme) + plot_18 = fig8.add_subplot(212) + plot_18.set_position([0.1, 0.13, 0.8, 0.32]) # x0, y0, width, height + plot_upper_vertical_build(plot_185, m_file_data, colour_scheme) + plot_lower_vertical_build(plot_18, m_file_data, colour_scheme) # Can only plot WP and turn structure if superconducting coil at the moment if m_file_data.data["i_tf_sup"].get_scan(scan) == 1: # TF coil with WP - plot_19 = fig8.add_subplot(231, aspect="equal") + plot_19 = fig9.add_subplot(231, aspect="equal") plot_19.set_position([ 0.025, 0.5, 0.45, 0.45, ]) # Half height, a bit wider, top left - plot_superconducting_tf_wp(plot_19, m_file_data, scan, fig8) + plot_superconducting_tf_wp(plot_19, m_file_data, scan, fig9) # TF coil turn structure - plot_20 = fig8.add_subplot(325, aspect="equal") + plot_20 = fig9.add_subplot(325, aspect="equal") plot_20.set_position([0.025, 0.1, 0.3, 0.3]) - plot_tf_turn(plot_20, fig8, m_file_data, scan) + plot_tf_turn(plot_20, fig9, m_file_data, scan) else: - plot_19 = fig8.add_subplot(211, aspect="equal") + plot_19 = fig9.add_subplot(211, aspect="equal") plot_19.set_position([0.06, 0.55, 0.675, 0.4]) plot_resistive_tf_wp(plot_19, m_file_data, scan, fig8) - plot_21 = fig9.add_subplot(111, aspect="equal") + plot_21 = fig10.add_subplot(111, aspect="equal") plot_tf_coil_structure(plot_21, m_file_data, scan, colour_scheme) - axes = fig10.subplots(nrows=3, ncols=1, sharex=True).flatten() + axes = fig11.subplots(nrows=3, ncols=1, sharex=True).flatten() plot_tf_stress(axes) - plot_23 = fig11.add_subplot(221) + plot_23 = fig12.add_subplot(221) plot_bootstrap_comparison(plot_23, m_file_data, scan) - plot_24 = fig11.add_subplot(224) + plot_24 = fig12.add_subplot(224) plot_h_threshold_comparison(plot_24, m_file_data, scan) - plot_25 = fig12.add_subplot(221) + plot_25 = fig13.add_subplot(221) plot_density_limit_comparison(plot_25, m_file_data, scan) - plot_26 = fig12.add_subplot(224) + plot_26 = fig13.add_subplot(224) plot_confinement_time_comparison(plot_26, m_file_data, scan) - plot_27 = fig13.add_subplot(111) + plot_27 = fig14.add_subplot(111) plot_current_profiles_over_time(plot_27, m_file_data, scan) - plot_28 = fig14.add_subplot(121, aspect="equal") - plot_cs_coil_structure(plot_28, fig14, m_file_data, scan) + plot_28 = fig15.add_subplot(121, aspect="equal") + plot_cs_coil_structure(plot_28, fig15, m_file_data, scan) - plot_29 = fig14.add_subplot(224, aspect="equal") - plot_cs_turn_structure(plot_29, fig14, m_file_data, scan) + plot_29 = fig15.add_subplot(224, aspect="equal") + plot_cs_turn_structure(plot_29, fig15, m_file_data, scan) - plot_30 = fig15.add_subplot(221, aspect="equal") + plot_30 = fig16.add_subplot(221, aspect="equal") plot_first_wall_top_down_cross_section(plot_30, m_file_data, scan) - plot_31 = fig15.add_subplot(122) + plot_31 = fig16.add_subplot(122) plot_first_wall_poloidal_cross_section(plot_31, m_file_data, scan) plot_32 = fig15.add_subplot(337) @@ -10157,8 +10160,6 @@ def main_plot( plot_33 = fig17.add_subplot(122) plot_fusion_rate_profiles(plot_33, m_file_data, scan) - plot_33 = fig17.add_subplot(122) - plot_fusion_rate_profiles(plot_33, m_file_data, scan) def main(args=None): From 0818fbae9ea0ed0eb50b0b0ead68212f148dd685 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 8 Sep 2025 13:25:59 +0100 Subject: [PATCH 09/12] Enhance fusion rate profile plotting by adding detailed power and density information for D-T, D-D, D-3He, alpha, and neutron reactions --- process/fusion_reactions.py | 2 +- process/io/plot_proc.py | 220 +++++++++++++++++++++++++++++++++--- process/physics.py | 2 +- 3 files changed, 209 insertions(+), 15 deletions(-) diff --git a/process/fusion_reactions.py b/process/fusion_reactions.py index 6ab8ce1f85..f6e6fdbf0d 100644 --- a/process/fusion_reactions.py +++ b/process/fusion_reactions.py @@ -298,7 +298,7 @@ def dhe3_reaction(self) -> None: x=self.plasma_profile.neprofile.profile_x, dx=self.plasma_profile.neprofile.profile_dx, ) - + physics_variables.fusrat_plasma_dhe3_profile = ( bosch_hale_reactivity( (physics_variables.ti / physics_variables.te) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index f670c6a1d8..1075203134 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -9836,13 +9836,13 @@ def plot_fw_90_deg_pipe_bend(ax, m_file_data, scan): ) -def plot_fusion_rate_profiles(axis, mfile_data, scan): +def plot_fusion_rate_profiles(axis, fig, mfile_data, scan): # 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 = [] - + fusrat_plasma_dt_profile = [ mfile_data.data[f"fusrat_plasma_dt_profile{i}"].get_scan(scan) for i in range(500) @@ -9907,7 +9907,7 @@ def plot_fusion_rate_profiles(axis, mfile_data, scan): 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, @@ -9953,6 +9953,199 @@ def plot_fusion_rate_profiles(axis, mfile_data, scan): ) ax2.tick_params(axis="y", which="minor", colors="blue") + # ================================================= + + # Add plasma volume, areas and shaping information + textstr_general = ( + f"Total fusion rate: {mfile_data.data['fusrat_total'].get_scan(scan):.4e} reactions/s\n" + f"Total fusion rate density: {mfile_data.data['fusden_total'].get_scan(scan):.4e} reactions/m3/s\n" + f"Plasma fusion rate density: {mfile_data.data['fusden_plasma'].get_scan(scan):.4e} reactions/m3/s\n" + ) + + axis.text( + 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, + }, + ) + + # ============================================================================ + + textstr_dt = ( + f"Total fusion power: {mfile_data.data['p_dt_total_mw'].get_scan(scan):.2f} MW\n" + f"Plasma fusion power: {mfile_data.data['p_plasma_dt_mw'].get_scan(scan):.2f} MW \n" + f"Beam fusion power: {mfile_data.data['p_beam_dt_mw'].get_scan(scan):.2f} MW\n" + ) + + axis.text( + 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, + }, + ) + + axis.text( + 0.24, + 0.8, + "$\\text{D - T}$", + fontsize=20, + verticalalignment="top", + transform=fig.transFigure, + ) + + # ================================================= + + textstr_dd = ( + f"Total fusion power: {mfile_data.data['p_dd_total_mw'].get_scan(scan):.2f} MW\n" + f"Tritium branching ratio: {mfile_data.data['f_dd_branching_trit'].get_scan(scan):.2f} \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, + }, + ) + + axis.text( + 0.22, + 0.68, + "$\\text{D - D}$", + fontsize=20, + verticalalignment="top", + transform=fig.transFigure, + ) + + # ================================================= + + textstr_dhe3 = f"Total fusion power: {mfile_data.data['p_dhe3_total_mw'].get_scan(scan):.2f} MW \n" + + axis.text( + 0.05, + 0.6, + textstr_dhe3, + fontsize=9, + verticalalignment="bottom", + horizontalalignment="left", + transform=fig.transFigure, + bbox={ + "boxstyle": "round", + "facecolor": "lightyellow", + "alpha": 1.0, + "linewidth": 2, + }, + ) + + axis.text( + 0.225, + 0.6, + "$\\text{D - 3He}$", + fontsize=20, + verticalalignment="top", + transform=fig.transFigure, + ) + + # ================================================= + + textstr_alpha = ( + f"Total power: {mfile_data.data['p_alpha_total_mw'].get_scan(scan):.2f} MW\n" + f"Plasma power: {mfile_data.data['p_plasma_alpha_mw'].get_scan(scan):.2f} MW\n" + f"Beam power: {mfile_data.data['p_beam_alpha_mw'].get_scan(scan):.2f} MW\n\n" + f"Rate density total: {mfile_data.data['fusden_alpha_total'].get_scan(scan):.4e} particles/m3/sec\n" + f"Rate density, plasma: {mfile_data.data['fusden_plasma_alpha'].get_scan(scan):.4e} particles/m3/sec\n\n" + f"Total power density: {mfile_data.data['pden_alpha_total_mw'].get_scan(scan):.4e} MW/m3\n" + f"Plasma power density: {mfile_data.data['pden_plasma_alpha_mw'].get_scan(scan):.4e} MW/m3\n\n" + f"Power per unit volume transferred to electrons: {mfile_data.data['f_pden_alpha_electron_mw'].get_scan(scan):.4e} MW/m3\n" + f"Power per unit volume transferred to ions: {mfile_data.data['f_pden_alpha_ions_mw'].get_scan(scan):.4e} MW/m3\n\n" + ) + + axis.text( + 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, + }, + ) + + axis.text( + 0.35, + 0.45, + "$\\alpha$", + fontsize=20, + verticalalignment="top", + transform=fig.transFigure, + ) + + # ================================================= + + textstr_neutron = ( + f"Total power: {mfile_data.data['p_neutron_total_mw'].get_scan(scan):.2f} MW\n" + f"Plasma power: {mfile_data.data['p_plasma_neutron_mw'].get_scan(scan):.2f} MW\n" + f"Beam power: {mfile_data.data['p_beam_neutron_mw'].get_scan(scan):.2f} MW\n\n" + f"Total power density: {mfile_data.data['pden_neutron_total_mw'].get_scan(scan):.4e} MW/m3\n" + f"Plasma power density: {mfile_data.data['pden_plasma_neutron_mw'].get_scan(scan):.4e} MW/m3\n" + ) + + axis.text( + 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, + }, + ) + + axis.text( + 0.25, + 0.2, + "$n$", + fontsize=20, + verticalalignment="top", + transform=fig.transFigure, + ) + def main_plot( fig1, @@ -9972,6 +10165,7 @@ def main_plot( fig15, fig16, fig17, + fig18, m_file_data, scan, imp="../data/lz_non_corona_14_elements/", @@ -10064,9 +10258,9 @@ def main_plot( plot_13 = fig4.add_subplot(4, 3, 12) plot_13.set_position([0.7, 0.125, 0.25, 0.15]) plot_qprofile(plot_13, demo_ranges, m_file_data, scan) - + plot_14 = fig5.add_subplot(122) - plot_fusion_rate_profiles(plot_14, m_file_data, scan) + plot_fusion_rate_profiles(plot_14, fig5, m_file_data, scan) # Plot poloidal cross-section plot_15 = fig6.add_subplot(121, aspect="equal") @@ -10149,17 +10343,13 @@ def main_plot( plot_31 = fig16.add_subplot(122) plot_first_wall_poloidal_cross_section(plot_31, m_file_data, scan) - plot_32 = fig15.add_subplot(337) + plot_32 = fig16.add_subplot(337) plot_fw_90_deg_pipe_bend(plot_32, m_file_data, scan) - plot_blkt_pipe_bends(fig16, m_file_data, scan) - - plot_33 = fig17.add_subplot(111, aspect="equal") - plot_main_power_flow(plot_33, m_file_data, scan, fig17) - - plot_33 = fig17.add_subplot(122) - plot_fusion_rate_profiles(plot_33, m_file_data, scan) + plot_blkt_pipe_bends(fig17, m_file_data, scan) + plot_33 = fig18.add_subplot(111, aspect="equal") + plot_main_power_flow(plot_33, m_file_data, scan, fig18) def main(args=None): @@ -10453,6 +10643,7 @@ def main(args=None): page15 = plt.figure(figsize=(12, 9), dpi=80) page16 = plt.figure(figsize=(12, 9), dpi=80) page17 = plt.figure(figsize=(12, 9), dpi=80) + page18 = plt.figure(figsize=(12, 9), dpi=80) # run main_plot main_plot( @@ -10473,6 +10664,7 @@ def main(args=None): page15, page16, page17, + page18, m_file, scan=scan, demo_ranges=demo_ranges, @@ -10498,6 +10690,7 @@ def main(args=None): pdf.savefig(page15) pdf.savefig(page16) pdf.savefig(page17) + pdf.savefig(page18) # show fig if option used if args.show: @@ -10520,6 +10713,7 @@ def main(args=None): plt.close(page15) plt.close(page16) plt.close(page17) + plt.close(page18) if __name__ == "__main__": diff --git a/process/physics.py b/process/physics.py index b62350ced9..08f9e0a5cf 100644 --- a/process/physics.py +++ b/process/physics.py @@ -4987,7 +4987,7 @@ def outplas(self): f"D-3He fusion rate at point {i}", f"fusrat_plasma_dhe3_profile{i}", physics_variables.fusrat_plasma_dhe3_profile[i], - ) + ) po.ovarre( self.outfile, "D-T fusion power: plasma (MW)", From 3d61ceb8dca08b9fc0f6f8d065d2f8655912eb3a Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 8 Sep 2025 14:37:49 +0100 Subject: [PATCH 10/12] Refactor constants usage in plotting and physics modules for consistency --- process/data_structure/physics_variables.py | 4 ++-- process/io/plot_proc.py | 24 ++++++++++----------- process/physics.py | 2 +- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index ec826aa88c..7d207f031a 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -514,8 +514,8 @@ fusrat_plasma_dt_profile: list[float] = None """Profile of D-T fusion reaction rate in plasma, (reactions/sec)""" -fusrat_plasma_dd_helion_profile: list[float] = None -"""Profile of D-D fusion reaction rate (helium branch) in plasma, (reactions/sec)""" +fusrat_plasma_dd_triton_profile: list[float] = None`` +"""Profile of D-D fusion reaction rate (tritium branch) in plasma, (reactions/sec)""" fusrat_plasma_dd_helion_profile: list[float] = None """Profile of D-D fusion reaction rate (helium branch) in plasma, (reactions/sec)""" diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 1075203134..4dd4c6d0a8 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -29,10 +29,10 @@ from matplotlib.path import Path import process.confinement_time as confine +import process.constants as constants import process.io.mfile as mf import process.superconducting_tf_coil as sctf from process.data_structure import physics_variables -from process.fortran import constants from process.geometry.blanket_geometry import ( blanket_geometry_double_null, blanket_geometry_single_null, @@ -9903,26 +9903,26 @@ def plot_fusion_rate_profiles(axis, fig, mfile_data, scan): 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, + 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, + 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, + 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, + np.array(fusrat_plasma_dhe3_profile) * constants.D_HELIUM_ENERGY, color=ax2.spines["right"].get_edgecolor(), linestyle="--", ) @@ -10015,7 +10015,7 @@ def plot_fusion_rate_profiles(axis, fig, mfile_data, scan): textstr_dd = ( f"Total fusion power: {mfile_data.data['p_dd_total_mw'].get_scan(scan):.2f} MW\n" - f"Tritium branching ratio: {mfile_data.data['f_dd_branching_trit'].get_scan(scan):.2f} \n" + f"Tritium branching ratio: {mfile_data.data['f_dd_branching_trit'].get_scan(scan):.4f} \n" ) axis.text( @@ -10036,7 +10036,7 @@ def plot_fusion_rate_profiles(axis, fig, mfile_data, scan): axis.text( 0.22, - 0.68, + 0.685, "$\\text{D - D}$", fontsize=20, verticalalignment="top", @@ -10045,11 +10045,11 @@ def plot_fusion_rate_profiles(axis, fig, mfile_data, scan): # ================================================= - textstr_dhe3 = f"Total fusion power: {mfile_data.data['p_dhe3_total_mw'].get_scan(scan):.2f} MW \n" + textstr_dhe3 = f"Total fusion power: {mfile_data.data['p_dhe3_total_mw'].get_scan(scan):.2f} MW \n\n" axis.text( 0.05, - 0.6, + 0.55, textstr_dhe3, fontsize=9, verticalalignment="bottom", @@ -10064,8 +10064,8 @@ def plot_fusion_rate_profiles(axis, fig, mfile_data, scan): ) axis.text( - 0.225, - 0.6, + 0.21, + 0.59, "$\\text{D - 3He}$", fontsize=20, verticalalignment="top", @@ -10106,7 +10106,7 @@ def plot_fusion_rate_profiles(axis, fig, mfile_data, scan): 0.35, 0.45, "$\\alpha$", - fontsize=20, + fontsize=22, verticalalignment="top", transform=fig.transFigure, ) diff --git a/process/physics.py b/process/physics.py index 08f9e0a5cf..99bf252e84 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1528,7 +1528,7 @@ def _trapped_particle_fraction_sauter( class Physics: def __init__(self, plasma_profile, current_drive): self.outfile = constants.NOUT - self.mfile = constants.mfile + self.mfile = constants.MFILE self.plasma_profile = plasma_profile self.current_drive = current_drive From d019de5ecf30a8b49e83a603c41670b9d7901227 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 18 Sep 2025 13:24:27 +0100 Subject: [PATCH 11/12] :recycle: Update fusion reaction calculations to use density profile --- process/data_structure/physics_variables.py | 2 +- process/fusion_reactions.py | 36 ++++++++++++++++----- process/io/plot_proc.py | 2 +- 3 files changed, 30 insertions(+), 10 deletions(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 7d207f031a..62d393b131 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -514,7 +514,7 @@ fusrat_plasma_dt_profile: list[float] = None """Profile of D-T fusion reaction rate in plasma, (reactions/sec)""" -fusrat_plasma_dd_triton_profile: list[float] = None`` +fusrat_plasma_dd_triton_profile: list[float] = None """Profile of D-D fusion reaction rate (tritium branch) in plasma, (reactions/sec)""" fusrat_plasma_dd_helion_profile: list[float] = None diff --git a/process/fusion_reactions.py b/process/fusion_reactions.py index f6e6fdbf0d..eb652971f7 100644 --- a/process/fusion_reactions.py +++ b/process/fusion_reactions.py @@ -216,8 +216,13 @@ def dt_reaction(self) -> None: * self.plasma_profile.teprofile.profile_y, dt, ) - * (physics_variables.f_deuterium * physics_variables.nd_fuel_ions) - * (physics_variables.f_tritium * physics_variables.nd_fuel_ions) + * physics_variables.f_deuterium + * physics_variables.f_tritium + * ( + self.plasma_profile.neprofile.profile_y + * (physics_variables.nd_fuel_ions / physics_variables.dene) + ) + ** 2 ) # Calculate the fusion reaction rate integral using Simpson's rule @@ -305,8 +310,13 @@ def dhe3_reaction(self) -> None: * self.plasma_profile.teprofile.profile_y, dhe3, ) - * (physics_variables.f_deuterium * physics_variables.nd_fuel_ions) - * (physics_variables.f_helium3 * physics_variables.nd_fuel_ions) + * physics_variables.f_deuterium + * physics_variables.f_helium3 + * ( + self.plasma_profile.neprofile.profile_y + * (physics_variables.nd_fuel_ions / physics_variables.dene) + ) + ** 2 ) # Reaction energy in MegaJoules [MJ] @@ -385,8 +395,13 @@ def dd_helion_reaction(self) -> None: * self.plasma_profile.teprofile.profile_y, dd1, ) - * (physics_variables.f_deuterium * physics_variables.nd_fuel_ions) - * (physics_variables.f_deuterium * physics_variables.nd_fuel_ions) + * physics_variables.f_deuterium + * physics_variables.f_deuterium + * ( + self.plasma_profile.neprofile.profile_y + * (physics_variables.nd_fuel_ions / physics_variables.dene) + ) + ** 2 ) # Reaction energy in MegaJoules [MJ] @@ -468,8 +483,13 @@ def dd_triton_reaction(self) -> None: * self.plasma_profile.teprofile.profile_y, dd2, ) - * (physics_variables.f_deuterium * physics_variables.nd_fuel_ions) - * (physics_variables.f_deuterium * physics_variables.nd_fuel_ions) + * physics_variables.f_deuterium + * physics_variables.f_deuterium + * ( + self.plasma_profile.neprofile.profile_y + * (physics_variables.nd_fuel_ions / physics_variables.dene) + ) + ** 2 ) # Reaction energy in MegaJoules [MJ] diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 4dd4c6d0a8..596e0a0269 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -9934,7 +9934,7 @@ def plot_fusion_rate_profiles(axis, fig, mfile_data, scan): ) axis.set_yscale("log") axis.grid(True, which="both", linestyle="--", alpha=0.5) - axis.set_xlim([0, 1.01]) + axis.set_xlim([0, 1.025]) axis.minorticks_on() axis.set_ylim([1e10, 1e19]) axis.yaxis.set_major_locator(plt.LogLocator(base=10.0, numticks=10)) From 2a0fef55e0ea9e13feee276fd8ad63072e852a5a Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 19 Sep 2025 12:39:01 +0100 Subject: [PATCH 12/12] Add total fusion rate profile calculation and plotting to fusion rate profiles --- process/io/plot_proc.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 596e0a0269..7aab35fad2 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -9862,6 +9862,14 @@ def plot_fusion_rate_profiles(axis, fig, mfile_data, scan): for i in range(500) ] + fusrat_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)) + ] + axis.spines["left"].set_color("red") axis.yaxis.label.set_color("black") axis.tick_params(axis="y", colors="red") @@ -9895,6 +9903,15 @@ def plot_fusion_rate_profiles(axis, fig, mfile_data, scan): linestyle="--", label=r"$\mathrm{D-3He}$", ) + axis.plot( + np.linspace(0, 1, len(fusrat_plasma_total_profile)), + fusrat_plasma_total_profile, + 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() @@ -9926,6 +9943,20 @@ def plot_fusion_rate_profiles(axis, fig, mfile_data, scan): 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]")