diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index 15ca2beb11..3ad0496bd8 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 @@ -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) @@ -11557,7 +11557,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)) @@ -11580,7 +11580,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): 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] @@ -11622,8 +11622,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): 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", @@ -11633,8 +11633,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): # 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, @@ -11685,7 +11685,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): 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 @@ -11737,8 +11737,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): # 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_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): 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_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): ) 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_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"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_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): ) axis.text( - 0.22, + 0.36, 0.685, "$\\text{D - D}$", fontsize=20, @@ -11804,7 +11804,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\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_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): ) axis.text( - 0.21, + 0.36, 0.59, "$\\text{D - 3He}$", fontsize=20, @@ -11830,12 +11830,12 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): 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_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): axis.text( 0.35, - 0.45, + 0.425, "$\\alpha$", fontsize=22, verticalalignment="top", @@ -11861,8 +11861,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): 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_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): ) axis.text( - 0.25, - 0.2, + 0.3, + 0.175, "$n$", fontsize=20, verticalalignment="top", @@ -11883,6 +11883,172 @@ 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, 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 = [] + fusden_plasma_dhe3_profile = [] + + vol_plasma = mfile.get("vol_plasma", 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) + ] + + 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) + ] + + fusden_plasma_total_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)) + ] + + # Convert each reaction-rate density into shell contributions, then integrate + # cumulatively to obtain the fusion-rate profile for that reaction. + reaction_profiles = { + "$\\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( + 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["$\\Sigma\\ \\Gamma_{\\mathrm{Total}}$"] = np.cumsum( + cum_fusrat_profiles["$\\Gamma_{\\mathrm{Total}}$"] + ) + + line_colors = { + "$\\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=(":" 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=":", + label="Total Fusion Rate", + ) + + half_total_rate = fusrat_total / 2 + half_total_rate_position = np.interp( + half_total_rate, + cum_fusrat_profiles["$\\Sigma\\ \\Gamma_{\\mathrm{Total}}$"], + rho, + ) + + axis.axvline( + x=half_total_rate_position, + color="tab:pink", + linestyle="-", + label="50% Total Fusion Rate", + ) + + # ================================================= + + 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.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="upper 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( axis: plt.Axes, mfile: MFile, @@ -12403,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) @@ -12792,15 +12959,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") @@ -16337,8 +16504,19 @@ def _add_page(name: str | None = None): _add_page("line_brem_power").add_subplot(121), m_file, scan, imp ) + plot_fusion_rate_density_profiles( + _add_page("fusion_rate_density").add_subplot(122), + pages["fusion_rate_density"], + m_file, + scan, + ) + + fusion_rate_page = _add_page("fusion_rate") plot_fusion_rate_profiles( - _add_page("fusion_rate").add_subplot(122), pages["fusion_rate"], m_file, scan + 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") @@ -16408,9 +16586,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(