From 305ac71f0d7bab25ab08b712f5e178896475fc88 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 14 Nov 2025 17:12:46 +0000 Subject: [PATCH 1/3] Add EBW/ECRH coupling efficiency graph plotting functionality --- process/current_drive.py | 2 +- process/io/plot_proc.py | 88 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 89 insertions(+), 1 deletion(-) diff --git a/process/current_drive.py b/process/current_drive.py index fb341cbdfc..19ac4c92b0 100644 --- a/process/current_drive.py +++ b/process/current_drive.py @@ -803,7 +803,7 @@ def electron_cyclotron_freethy( 1 / (2 * np.pi) * np.sqrt( - (nd_plasma_electrons_vol_avg / 1.0e19) + (nd_plasma_electrons_vol_avg) * constants.ELECTRON_CHARGE**2 / (constants.ELECTRON_MASS * constants.EPSILON0) ) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 22612945c9..da03e377b4 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -36,6 +36,7 @@ import process.io.mfile as mf import process.superconducting_tf_coil as sctf from process.build import Build +from process.current_drive import ElectronBernstein, ElectronCyclotron from process.data_structure import impurity_radiation_module, physics_variables from process.geometry.blanket_geometry import ( blanket_geometry_double_null, @@ -12434,6 +12435,91 @@ def plot_ion_charge_profile(axis, mfile_data, scan): axis.grid(which="both", linestyle="--", alpha=0.5) +def plot_ebw_ecrh_coupling_graph(axis, mfile_data, scan): + # Plot EBW and ECRH coupling efficiency graph + ebw = ElectronBernstein(plasma_profile=0) + ecrg = ElectronCyclotron(plasma_profile=0) + b_on_axis = mfile_data.data["b_plasma_toroidal_on_axis"].get_scan(scan) + bs = np.linspace(b_on_axis - 2.0, b_on_axis + 2.0, 500) + # Use a color map for harmonics + colors = ["red", "green", "blue"] + linestyles = ["-", "--"] # EBW: solid, ECRH: dashed + + for idx, n_harmonic in enumerate(range(1, 4)): + eta_ebw_vals = [] + # For ECRH, store results for both wave modes (0: O-mode, 1: X-mode) + eta_ecrh_vals_omode = [] + eta_ecrh_vals_xmode = [] + for b in bs: + eta_ebw = ebw.electron_bernstein_freethy( + te=mfile_data.data["temp_plasma_electron_vol_avg_kev"].get_scan(scan), + rmajor=mfile_data.data["rmajor"].get_scan(scan), + dene20=mfile_data.data["nd_plasma_electrons_vol_avg"].get_scan(scan) + / 1e20, + b_plasma_toroidal_on_axis=b, + n_ecrh_harmonic=n_harmonic, + xi_ebw=0.8, + ) + eta_ecrh_omode = ecrg.electron_cyclotron_freethy( + te=mfile_data.data["temp_plasma_electron_vol_avg_kev"].get_scan(scan), + zeff=mfile_data.data["zeff"].get_scan(scan), + rmajor=mfile_data.data["rmajor"].get_scan(scan), + nd_plasma_electrons_vol_avg=mfile_data.data[ + "nd_plasma_electrons_vol_avg" + ].get_scan(scan), + b_plasma_toroidal_on_axis=b, + n_ecrh_harmonic=n_harmonic, + i_ecrh_wave_mode=0, # O-mode + ) + eta_ecrh_xmode = ecrg.electron_cyclotron_freethy( + te=mfile_data.data["temp_plasma_electron_vol_avg_kev"].get_scan(scan), + zeff=mfile_data.data["zeff"].get_scan(scan), + rmajor=mfile_data.data["rmajor"].get_scan(scan), + nd_plasma_electrons_vol_avg=mfile_data.data[ + "nd_plasma_electrons_vol_avg" + ].get_scan(scan), + b_plasma_toroidal_on_axis=b, + n_ecrh_harmonic=n_harmonic, + i_ecrh_wave_mode=1, # X-mode + ) + eta_ebw_vals.append(eta_ebw) + eta_ecrh_vals_omode.append(eta_ecrh_omode) + eta_ecrh_vals_xmode.append(eta_ecrh_xmode) + # EBW: solid, ECRH O-mode: dashed, ECRH X-mode: dotted, same color for same harmonic + axis.plot( + bs, + eta_ebw_vals, + label=f"EBW (harmonic {n_harmonic})", + color=colors[idx], + linestyle=linestyles[0], + ) + axis.plot( + bs, + eta_ecrh_vals_omode, + label=f"ECRH O-mode (harmonic {n_harmonic})", + color=colors[idx], + linestyle="--", + ) + axis.plot( + bs, + eta_ecrh_vals_xmode, + label=f"ECRH X-mode (harmonic {n_harmonic})", + color=colors[idx], + linestyle=":", + ) + axis.set_xlabel("On axis toroidal B-field [T]") + axis.set_ylabel("Coupling Efficiency") + axis.set_title("EBW/ECRH Coupling Efficiency vs Toroidal B-field") + axis.legend() + axis.grid(True) + # Plot a vertical line at the on-axis value of the toroidal B-field + b_on_axis = mfile_data.data["b_plasma_toroidal_on_axis"].get_scan(scan) + axis.axvline( + b_on_axis, color="black", linestyle="-", linewidth=2.5, label="On-axis $B_T$" + ) + axis.minorticks_on() + + def main_plot( fig0, fig1, @@ -12744,6 +12830,8 @@ def main_plot( ax24.set_position([0.08, 0.35, 0.84, 0.57]) plot_system_power_profiles_over_time(ax24, m_file_data, scan, fig26) + plot_ebw_ecrh_coupling_graph(fig26.add_subplot(111), m_file_data, scan) + def main(args=None): # TODO The use of globals here isn't ideal, but is required to get main() From 7baeb31e3ae56a8192262725d77897497dd17955 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 2 Dec 2025 09:22:40 +0000 Subject: [PATCH 2/3] Refactor EBW/ECRH coupling graph plotting: update B-field range and improve axis labels --- process/io/plot_proc.py | 78 ++++++++++++++++++++++------------------- 1 file changed, 42 insertions(+), 36 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index da03e377b4..faedd82c28 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -12440,7 +12440,7 @@ def plot_ebw_ecrh_coupling_graph(axis, mfile_data, scan): ebw = ElectronBernstein(plasma_profile=0) ecrg = ElectronCyclotron(plasma_profile=0) b_on_axis = mfile_data.data["b_plasma_toroidal_on_axis"].get_scan(scan) - bs = np.linspace(b_on_axis - 2.0, b_on_axis + 2.0, 500) + bs = np.linspace(0.0, b_on_axis + 2.0, 500) # Use a color map for harmonics colors = ["red", "green", "blue"] linestyles = ["-", "--"] # EBW: solid, ECRH: dashed @@ -12508,7 +12508,7 @@ def plot_ebw_ecrh_coupling_graph(axis, mfile_data, scan): linestyle=":", ) axis.set_xlabel("On axis toroidal B-field [T]") - axis.set_ylabel("Coupling Efficiency") + axis.set_ylabel("Current drive efficiency [A/W]") axis.set_title("EBW/ECRH Coupling Efficiency vs Toroidal B-field") axis.legend() axis.grid(True) @@ -12548,6 +12548,7 @@ def main_plot( fig24, fig25, fig26, + fig27, m_file_data, scan, imp="../data/lz_non_corona_14_elements/", @@ -12690,9 +12691,16 @@ def main_plot( plot_magnetic_fields_in_plasma(fig10.add_subplot(122), m_file_data, scan) plot_beta_profiles(fig10.add_subplot(221), m_file_data, scan) + plot_ebw_ecrh_coupling_graph(fig11.add_subplot(111), m_file_data, scan) + + plot_bootstrap_comparison(fig12.add_subplot(221), m_file_data, scan) + plot_h_threshold_comparison(fig12.add_subplot(224), m_file_data, scan) + plot_density_limit_comparison(fig13.add_subplot(221), m_file_data, scan) + plot_confinement_time_comparison(fig13.add_subplot(224), m_file_data, scan) + # Plot poloidal cross-section poloidal_cross_section( - fig11.add_subplot(121, aspect="equal"), + fig14.add_subplot(121, aspect="equal"), m_file_data, scan, demo_ranges, @@ -12701,7 +12709,7 @@ def main_plot( # Plot toroidal cross-section toroidal_cross_section( - fig11.add_subplot(122, aspect="equal"), + fig14.add_subplot(122, aspect="equal"), m_file_data, scan, demo_ranges, @@ -12709,19 +12717,19 @@ def main_plot( ) # Plot color key - ax17 = fig11.add_subplot(222) + ax17 = fig14.add_subplot(222) ax17.set_position([0.5, 0.5, 0.5, 0.5]) color_key(ax17, m_file_data, scan, colour_scheme) - ax18 = fig12.add_subplot(211) + ax18 = fig15.add_subplot(211) ax18.set_position([0.1, 0.33, 0.8, 0.6]) plot_radial_build(ax18, m_file_data, colour_scheme) # Make each axes smaller vertically to leave room for the legend - ax185 = fig13.add_subplot(211) + ax185 = fig16.add_subplot(211) ax185.set_position([0.1, 0.61, 0.8, 0.32]) - ax18b = fig13.add_subplot(212) + ax18b = fig16.add_subplot(212) ax18b.set_position([0.1, 0.13, 0.8, 0.32]) plot_upper_vertical_build(ax185, m_file_data, colour_scheme) plot_lower_vertical_build(ax18b, m_file_data, colour_scheme) @@ -12729,57 +12737,53 @@ def main_plot( # 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 - ax19 = fig14.add_subplot(221, aspect="equal") + ax19 = fig17.add_subplot(221, aspect="equal") ax19.set_position([ 0.025, 0.45, 0.5, 0.5, ]) # Half height, a bit wider, top left - plot_superconducting_tf_wp(ax19, m_file_data, scan, fig14) + plot_superconducting_tf_wp(ax19, m_file_data, scan, fig17) # TF coil turn structure - ax20 = fig15.add_subplot(325, aspect="equal") + ax20 = fig18.add_subplot(325, aspect="equal") ax20.set_position([0.025, 0.5, 0.4, 0.4]) - plot_tf_cable_in_conduit_turn(ax20, fig15, m_file_data, scan) - plot_205 = fig15.add_subplot(223, aspect="equal") + plot_tf_cable_in_conduit_turn(ax20, fig18, m_file_data, scan) + plot_205 = fig18.add_subplot(223, aspect="equal") plot_205.set_position([0.075, 0.1, 0.3, 0.3]) - plot_cable_in_conduit_cable(plot_205, fig15, m_file_data, scan) + plot_cable_in_conduit_cable(plot_205, fig18, m_file_data, scan) else: - ax19 = fig14.add_subplot(211, aspect="equal") + ax19 = fig17.add_subplot(211, aspect="equal") ax19.set_position([0.06, 0.55, 0.675, 0.4]) - plot_resistive_tf_wp(ax19, m_file_data, scan, fig14) + plot_resistive_tf_wp(ax19, m_file_data, scan, fig17) plot_tf_coil_structure( - fig16.add_subplot(111, aspect="equal"), m_file_data, scan, colour_scheme + fig19.add_subplot(111, aspect="equal"), m_file_data, scan, colour_scheme ) - plot_plasma_outboard_toroidal_ripple_map(fig17, m_file_data, scan) + plot_plasma_outboard_toroidal_ripple_map(fig20, m_file_data, scan) - axes = fig18.subplots(nrows=3, ncols=1, sharex=True).flatten() + axes = fig21.subplots(nrows=3, ncols=1, sharex=True).flatten() plot_tf_stress(axes) - plot_bootstrap_comparison(fig19.add_subplot(221), m_file_data, scan) - plot_h_threshold_comparison(fig19.add_subplot(224), m_file_data, scan) - plot_density_limit_comparison(fig20.add_subplot(221), m_file_data, scan) - plot_confinement_time_comparison(fig20.add_subplot(224), m_file_data, scan) - plot_current_profiles_over_time(fig21.add_subplot(111), m_file_data, scan) + plot_current_profiles_over_time(fig22.add_subplot(111), m_file_data, scan) plot_cs_coil_structure( - fig22.add_subplot(121, aspect="equal"), fig22, m_file_data, scan + fig23.add_subplot(121, aspect="equal"), fig23, m_file_data, scan ) plot_cs_turn_structure( - fig22.add_subplot(224, aspect="equal"), fig22, m_file_data, scan + fig23.add_subplot(224, aspect="equal"), fig23, m_file_data, scan ) plot_first_wall_top_down_cross_section( - fig23.add_subplot(221, aspect="equal"), m_file_data, scan + fig24.add_subplot(221, aspect="equal"), m_file_data, scan ) - plot_first_wall_poloidal_cross_section(fig23.add_subplot(122), m_file_data, scan) - plot_fw_90_deg_pipe_bend(fig23.add_subplot(337), m_file_data, scan) + plot_first_wall_poloidal_cross_section(fig24.add_subplot(122), m_file_data, scan) + plot_fw_90_deg_pipe_bend(fig24.add_subplot(337), m_file_data, scan) - plot_blkt_pipe_bends(fig24, m_file_data, scan) - ax_blanket = fig24.add_subplot(122, aspect="equal") + plot_blkt_pipe_bends(fig25, m_file_data, scan) + ax_blanket = fig25.add_subplot(122, aspect="equal") plot_blanket(ax_blanket, m_file_data, scan, colour_scheme) plot_firstwall(ax_blanket, m_file_data, scan, colour_scheme) ax_blanket.set_xlabel("Radial position [m]") @@ -12822,15 +12826,13 @@ def main_plot( ) plot_main_power_flow( - fig25.add_subplot(111, aspect="equal"), m_file_data, scan, fig25 + fig26.add_subplot(111, aspect="equal"), m_file_data, scan, fig26 ) - ax24 = fig26.add_subplot(111) + ax24 = fig27.add_subplot(111) # set_position([left, bottom, width, height]) -> height ~ 0.66 => ~2/3 of page height ax24.set_position([0.08, 0.35, 0.84, 0.57]) - plot_system_power_profiles_over_time(ax24, m_file_data, scan, fig26) - - plot_ebw_ecrh_coupling_graph(fig26.add_subplot(111), m_file_data, scan) + plot_system_power_profiles_over_time(ax24, m_file_data, scan, fig27) def main(args=None): @@ -13150,6 +13152,7 @@ def main(args=None): page24 = plt.figure(figsize=(12, 9), dpi=80) page25 = plt.figure(figsize=(12, 9), dpi=80) page26 = plt.figure(figsize=(12, 9), dpi=80) + page27 = plt.figure(figsize=(12, 9), dpi=80) # run main_plot main_plot( @@ -13180,6 +13183,7 @@ def main(args=None): page24, page25, page26, + page27, m_file, scan=scan, demo_ranges=demo_ranges, @@ -13215,6 +13219,7 @@ def main(args=None): pdf.savefig(page24) pdf.savefig(page25) pdf.savefig(page26) + pdf.savefig(page27) # show fig if option used if args.show: @@ -13247,6 +13252,7 @@ def main(args=None): plt.close(page24) plt.close(page25) plt.close(page26) + plt.close(page27) if __name__ == "__main__": From c4eacba348e4ad972fa09ba88c040d406eb69f6d Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 4 Dec 2025 09:09:54 +0000 Subject: [PATCH 3/3] Add EBW coupling efficiency to CurrentDrive output and update plot function --- process/current_drive.py | 6 ++++++ process/io/plot_proc.py | 4 ++-- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/process/current_drive.py b/process/current_drive.py index 19ac4c92b0..a2addcf001 100644 --- a/process/current_drive.py +++ b/process/current_drive.py @@ -2136,6 +2136,12 @@ def output_current_drive(self): "(n_ecrh_harmonic)", current_drive_variables.n_ecrh_harmonic, ) + po.ovarre( + self.outfile, + "EBW coupling efficiency", + "(xi_ebw)", + current_drive_variables.xi_ebw, + ) if current_drive_variables.i_hcd_primary == 13: po.ovarin( self.outfile, diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index faedd82c28..eab42a11dd 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -12458,7 +12458,7 @@ def plot_ebw_ecrh_coupling_graph(axis, mfile_data, scan): / 1e20, b_plasma_toroidal_on_axis=b, n_ecrh_harmonic=n_harmonic, - xi_ebw=0.8, + xi_ebw=mfile_data.data["xi_ebw"].get_scan(scan), ) eta_ecrh_omode = ecrg.electron_cyclotron_freethy( te=mfile_data.data["temp_plasma_electron_vol_avg_kev"].get_scan(scan), @@ -12603,7 +12603,7 @@ def main_plot( plot_power_info(fig1.add_subplot(235), m_file_data, scan) # Current drive - plot_current_drive_info(fig1.add_subplot(236), m_file_data, scan) + # plot_current_drive_info(fig1.add_subplot(236), m_file_data, scan) fig1.subplots_adjust(wspace=0.25, hspace=0.25) ax7 = fig2.add_subplot(121)