Skip to content

Add EBW/ECRH coupling efficiency graph plotting functionality #3984

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
timothy-nunn merged 3 commits into from
Dec 5, 2025
Merged

Add EBW/ECRH coupling efficiency graph plotting functionality #3984

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
305ac71
Add EBW/ECRH coupling efficiency graph plotting functionality
chris-ashe Nov 14, 2025
7baeb31
Refactor EBW/ECRH coupling graph plotting: update B-field range and i…
chris-ashe Dec 2, 2025
c4eacba
Add EBW coupling efficiency to CurrentDrive output and update plot fu…
chris-ashe Dec 4, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 7 additions & 1 deletion process/current_drive.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
)
Expand Down Expand Up @@ -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,
Expand Down
160 changes: 127 additions & 33 deletions process/io/plot_proc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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(0.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=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),
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("Current drive efficiency [A/W]")
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,
Expand Down Expand Up @@ -12462,6 +12548,7 @@ def main_plot(
fig24,
fig25,
fig26,
fig27,
m_file_data,
scan,
imp="../data/lz_non_corona_14_elements/",
Expand Down Expand Up @@ -12516,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)
Expand Down Expand Up @@ -12604,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,
Expand All @@ -12615,85 +12709,81 @@ 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,
colour_scheme,
)

# 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)

# 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]")
Expand Down Expand Up @@ -12736,13 +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_system_power_profiles_over_time(ax24, m_file_data, scan, fig27)


def main(args=None):
Expand Down Expand Up @@ -13062,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(
Expand Down Expand Up @@ -13092,6 +13183,7 @@ def main(args=None):
page24,
page25,
page26,
page27,
m_file,
scan=scan,
demo_ranges=demo_ranges,
Expand Down Expand Up @@ -13127,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:
Expand Down Expand Up @@ -13159,6 +13252,7 @@ def main(args=None):
plt.close(page24)
plt.close(page25)
plt.close(page26)
plt.close(page27)


if __name__ == "__main__":
Expand Down
Loading