diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 7f6adaa12..f84d8e94b 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1308,6 +1308,12 @@ len_plasma_debye_electron_profile: list[float] = None """Profile of electron Debye length in plasma (m)""" +radius_plasma_deuteron_toroidal_larmor_isotropic_profile: list[float] = None +"""Profile of deuteron toroidal Larmor radius in plasma, assuming equal speeds in all directions (m)""" + +radius_plasma_triton_toroidal_larmor_isotropic_profile: list[float] = None +"""Profile of triton toroidal Larmor radius in plasma, assuming equal speeds in all directions (m)""" + len_plasma_debye_electron_vol_avg: float = None """Volume averaged electron Debye length in plasma (m)""" @@ -1320,6 +1326,12 @@ vel_plasma_triton_profile: list[float] = None """Profile of triton thermal velocity in plasma (m/s)""" +vel_plasma_alpha_thermal_profile: list[float] = None +"""Profile of thermal alpha particle velocity in plasma (m/s)""" + +vel_plasma_alpha_birth: float = None +"""Birth velocity of alpha particles in plasma (m/s)""" + plasma_coulomb_log_electron_electron_profile: list[float] = None """Profile of electron-electron Coulomb logarithm in plasma""" @@ -1329,9 +1341,18 @@ plasma_coulomb_log_electron_triton_profile: list[float] = None """Profile of electron-triton Coulomb logarithm in plasma""" +plasma_coulomb_log_deuteron_triton_profile: list[float] = None +"""Profile of deuteron-triton Coulomb logarithm in plasma""" + +plasma_coulomb_log_electron_alpha_thermal_profile: list[float] = None +"""Profile of electron-alpha Coulomb logarithm in plasma""" + freq_plasma_electron_profile: list[float] = None """Electron plasma frequency profile (Hz)""" +freq_plasma_deuteron_profile: list[float] = None +"""Deuteron plasma frequency profile (Hz)""" + freq_plasma_larmor_toroidal_electron_profile: list[float] = None """Profile of electron Larmor frequency in plasma due to toroidal magnetic field (Hz)""" @@ -1341,6 +1362,9 @@ freq_plasma_larmor_toroidal_triton_profile: list[float] = None """Profile of triton Larmor frequency in plasma due to toroidal magnetic field (Hz)""" +freq_plasma_upper_hybrid_profile: list[float] = None +"""Profile of upper hybrid frequency in plasma (Hz)""" + def init_physics_module(): """Initialise the physics module""" @@ -1659,17 +1683,25 @@ def init_physics_variables(): n_charge_plasma_effective_mass_weighted_vol_avg, \ j_plasma_bootstrap_sauter_profile, \ len_plasma_debye_electron_profile, \ + radius_plasma_deuteron_toroidal_larmor_isotropic_profile, \ + radius_plasma_triton_toroidal_larmor_isotropic_profile, \ len_plasma_debye_electron_vol_avg, \ vel_plasma_electron_profile, \ vel_plasma_deuteron_profile, \ vel_plasma_triton_profile, \ + vel_plasma_alpha_thermal_profile, \ + vel_plasma_alpha_birth, \ plasma_coulomb_log_electron_electron_profile, \ plasma_coulomb_log_electron_deuteron_profile, \ plasma_coulomb_log_electron_triton_profile, \ + plasma_coulomb_log_deuteron_triton_profile, \ + plasma_coulomb_log_electron_alpha_thermal_profile, \ freq_plasma_electron_profile, \ + freq_plasma_deuteron_profile, \ freq_plasma_larmor_toroidal_electron_profile, \ freq_plasma_larmor_toroidal_deuteron_profile, \ - freq_plasma_larmor_toroidal_triton_profile + freq_plasma_larmor_toroidal_triton_profile, \ + freq_plasma_upper_hybrid_profile m_beam_amu = 0.0 m_fuel_amu = 0.0 @@ -1928,14 +1960,22 @@ def init_physics_variables(): n_charge_plasma_effective_profile = [] n_charge_plasma_effective_mass_weighted_vol_avg = 0.0 len_plasma_debye_electron_profile = [] + radius_plasma_deuteron_toroidal_larmor_isotropic_profile = [] + radius_plasma_triton_toroidal_larmor_isotropic_profile = [] len_plasma_debye_electron_vol_avg = 0.0 vel_plasma_electron_profile = [] vel_plasma_deuteron_profile = [] vel_plasma_triton_profile = [] + vel_plasma_alpha_thermal_profile = [] + vel_plasma_alpha_birth = 0.0 plasma_coulomb_log_electron_electron_profile = [] plasma_coulomb_log_electron_deuteron_profile = [] plasma_coulomb_log_electron_triton_profile = [] + plasma_coulomb_log_deuteron_triton_profile = [] + plasma_coulomb_log_electron_alpha_thermal_profile = [] freq_plasma_electron_profile = [] + freq_plasma_deuteron_profile = [] freq_plasma_larmor_toroidal_electron_profile = [] freq_plasma_larmor_toroidal_deuteron_profile = [] freq_plasma_larmor_toroidal_triton_profile = [] + freq_plasma_upper_hybrid_profile = [] diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index d97d0bbcb..288e63405 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -12560,6 +12560,58 @@ def plot_ebw_ecrh_coupling_graph(axis: plt.Axes, mfile: mf.MFile, scan: int): axis.minorticks_on() +def plot_larmor_radius_profile(axis: plt.Axes, mfile_data: mf.MFile, scan: int): + """Plot the Larmor radius profile on the given axis.""" + radius_plasma_deuteron_larmor_profile = [ + mfile_data.data[ + f"radius_plasma_deuteron_toroidal_larmor_isotropic_profile{i}" + ].get_scan(scan) + for i in range( + 2 * int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + radius_plasma_triton_larmor_profile = [ + mfile_data.data[ + f"radius_plasma_triton_toroidal_larmor_isotropic_profile{i}" + ].get_scan(scan) + for i in range( + 2 * int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + radius_plasma_deuteron_larmor_profile_mm = [ + radius * 1e3 for radius in radius_plasma_deuteron_larmor_profile + ] + + radius_plasma_triton_larmor_profile_mm = [ + radius * 1e3 for radius in radius_plasma_triton_larmor_profile + ] + + axis.plot( + np.linspace(-1, 1, len(radius_plasma_deuteron_larmor_profile_mm)), + radius_plasma_deuteron_larmor_profile_mm, + color="red", + linestyle="-", + label=r"$\rho_{Larmor,toroidal,D}$", + ) + + axis.plot( + np.linspace(-1, 1, len(radius_plasma_triton_larmor_profile_mm)), + radius_plasma_triton_larmor_profile_mm, + color="green", + linestyle="-", + label=r"$\rho_{Larmor,toroidal,T}$", + ) + + axis.set_ylabel(r"Larmor Radii [mm]") + axis.set_title(r" Toroidal Larmor Radii ($v_{\perp}^2 = 2v_{th}^2$)") + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.minorticks_on() + axis.legend() + + def plot_debye_length_profile(axis: plt.Axes, mfile_data: mf.MFile, scan: int): """Plot the Debye length profile on the given axis. @@ -12623,6 +12675,12 @@ def plot_velocity_profile(axis, mfile_data, scan): mfile_data.data[f"vel_plasma_triton_profile{i}"].get_scan(scan) for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) ] + vel_plasma_alpha_thermal_profile = [ + mfile_data.data[f"vel_plasma_alpha_thermal_profile{i}"].get_scan(scan) + for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) + ] + + vel_plasma_alpha_birth = mfile_data.data["vel_plasma_alpha_birth"].get_scan(scan) axis.plot( np.linspace(0, 1, len(vel_plasma_electron_profile)), @@ -12634,7 +12692,7 @@ def plot_velocity_profile(axis, mfile_data, scan): axis.plot( np.linspace(0, 1, len(vel_plasma_deuteron_profile)), vel_plasma_deuteron_profile, - color="red", + color="pink", linestyle="-", label=r"$v_{D}$", ) @@ -12645,6 +12703,20 @@ def plot_velocity_profile(axis, mfile_data, scan): linestyle="-", label=r"$v_{T}$", ) + axis.plot( + np.linspace(0, 1, len(vel_plasma_alpha_thermal_profile)), + vel_plasma_alpha_thermal_profile, + color="red", + linestyle="-", + label=r"$v_{\alpha,thermal}$", + ) + axis.axhline( + vel_plasma_alpha_birth, + color="red", + linestyle="--", + linewidth=1.5, + label=r"$v_{\alpha,birth}$", + ) axis.set_yscale("log") axis.set_ylabel("Velocity [m/s]") @@ -12680,24 +12752,82 @@ def plot_electron_frequency_profile(axis, mfile_data, scan): ) ] + freq_plasma_upper_hybrid_electron_profile = [ + mfile_data.data[f"freq_plasma_upper_hybrid_profile{i}"].get_scan(scan) + for i in range( + 2 * int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + axis.plot( np.linspace(-1, 1, len(freq_plasma_larmor_toroidal_electron_profile)), np.array(freq_plasma_larmor_toroidal_electron_profile) / 1e9, color="red", linestyle="-", - label=r"$f_{Larmor,toroidal,e}$", + label=r"$f_{Larmor,toroidal,e}$ | Fundamental", ) + + axis.plot( + np.linspace(-1, 1, len(freq_plasma_larmor_toroidal_electron_profile)), + 2 * np.array(freq_plasma_larmor_toroidal_electron_profile) / 1e9, + color="red", + linestyle="--", + label=r"$f_{Larmor,toroidal,e}$ | 2nd harmonic", + ) + + axis.plot( + np.linspace(-1, 1, len(freq_plasma_larmor_toroidal_electron_profile)), + 3 * np.array(freq_plasma_larmor_toroidal_electron_profile) / 1e9, + color="red", + linestyle=":", + label=r"$f_{Larmor,toroidal,e}$ | 3rd harmonic", + ) + x = np.linspace(0, 1, len(freq_plasma_electron_profile)) y = np.array(freq_plasma_electron_profile) / 1e9 # original curve - axis.plot(x, y, color="blue", linestyle="-", label=r"$\omega_{p,e}$") + axis.plot( + x, y, color="blue", linestyle="-", label=r"$\omega_{p,e}$ | Plasma Frequency" + ) # mirrored across the y-axis (drawn at negative rho) axis.plot(-x, y, color="blue", linestyle="-", label="_nolegend_") + + axis.plot( + np.linspace(-1, 1, len(freq_plasma_upper_hybrid_electron_profile)), + np.array(freq_plasma_upper_hybrid_electron_profile) / 1e9, + color="purple", + linestyle="-", + label=r"$\omega_{UH,e}$ | Upper Hybrid", + ) + axis.set_xlim(-1.025, 1.025) + axis.set_ylim(None, max(freq_plasma_larmor_toroidal_electron_profile) / 1e9 * 1.6) + axis.set_xlabel("$\\rho$ [r/a]") axis.set_ylabel("Frequency [GHz]") axis.grid(True, which="both", linestyle="--", alpha=0.5) + # Add secondary x-axis showing radius in metres below the primary axis + ax2 = axis.twiny() + rmajor = mfile_data.get("rmajor", scan=scan) + rminor = mfile_data.get("rminor", scan=scan) + + # Convert normalized radius to actual radius + # rho ranges from -1 to 1, which corresponds to r = rmajor - rminor to rmajor + rminor + rho_ticks = np.array([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1]) + r_ticks = rmajor + rho_ticks * rminor + + ax2.set_xticks(rho_ticks) + ax2.set_xticklabels([f"{r:.2f}" for r in r_ticks]) + ax2.set_xlabel("Radius [m]") + ax2.minorticks_on() + ax2.set_xlim(axis.get_xlim()) + + # Move secondary axis to the bottom + ax2.xaxis.set_ticks_position("bottom") + ax2.xaxis.set_label_position("bottom") + ax2.spines["bottom"].set_position(("outward", 30)) + axis.legend() @@ -12771,6 +12901,18 @@ def plot_plasma_coloumb_logarithms(axis, mfile_data, scan): for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) ] + plasma_coulomb_log_deuteron_triton_profile = [ + mfile_data.data[f"plasma_coulomb_log_deuteron_triton_profile{i}"].get_scan(scan) + for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) + ] + + plasma_coulomb_log_electron_alpha_thermal_profile = [ + mfile_data.data[ + f"plasma_coulomb_log_electron_alpha_thermal_profile{i}" + ].get_scan(scan) + for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) + ] + axis.plot( np.linspace(0, 1, len(plasma_coulomb_log_electron_electron_profile)), plasma_coulomb_log_electron_electron_profile, @@ -12782,7 +12924,7 @@ def plot_plasma_coloumb_logarithms(axis, mfile_data, scan): axis.plot( np.linspace(0, 1, len(plasma_coulomb_log_electron_deuteron_profile)), plasma_coulomb_log_electron_deuteron_profile, - color="red", + color="pink", linestyle="-", label=r"$ln \Lambda_{e-D}$", ) @@ -12795,6 +12937,22 @@ def plot_plasma_coloumb_logarithms(axis, mfile_data, scan): label=r"$ln \Lambda_{e-T}$", ) + axis.plot( + np.linspace(0, 1, len(plasma_coulomb_log_deuteron_triton_profile)), + plasma_coulomb_log_deuteron_triton_profile, + color="orange", + linestyle="-", + label=r"$ln \Lambda_{D-T}$", + ) + + axis.plot( + np.linspace(0, 1, len(plasma_coulomb_log_electron_alpha_thermal_profile)), + plasma_coulomb_log_electron_alpha_thermal_profile, + color="red", + linestyle="-", + label=r"$ln \Lambda_{e-\alpha,thermal}$", + ) + axis.set_ylabel("Coulomb Logarithm") axis.set_xlabel("$\\rho \\ [r/a]$") axis.grid(True, which="both", linestyle="--", alpha=0.5) @@ -13252,19 +13410,17 @@ def main_plot( plot_debye_length_profile(figs[15].add_subplot(232), m_file, scan) plot_velocity_profile(figs[15].add_subplot(233), m_file, scan) + plot_plasma_coloumb_logarithms(figs[15].add_subplot(231), m_file, scan) - ax_ion = figs[15].add_subplot(414) - ax_electron = figs[15].add_subplot(413, sharex=ax_ion) - - plot_electron_frequency_profile(ax_electron, m_file, scan) + plot_electron_frequency_profile(figs[15].add_subplot(212), m_file, scan) - plot_ion_frequency_profile(ax_ion, m_file, scan) + plot_ion_frequency_profile(figs[16].add_subplot(311), m_file, scan) - plot_plasma_coloumb_logarithms(figs[15].add_subplot(231), m_file, scan) + plot_larmor_radius_profile(figs[16].add_subplot(313), m_file, scan) # Plot poloidal cross-section poloidal_cross_section( - figs[16].add_subplot(121, aspect="equal"), + figs[17].add_subplot(121, aspect="equal"), m_file, scan, demo_ranges, @@ -13274,7 +13430,7 @@ def main_plot( # Plot toroidal cross-section toroidal_cross_section( - figs[16].add_subplot(122, aspect="equal"), + figs[17].add_subplot(122, aspect="equal"), m_file, scan, demo_ranges, @@ -13282,19 +13438,19 @@ def main_plot( ) # Plot color key - ax17 = figs[16].add_subplot(222) + ax17 = figs[17].add_subplot(222) ax17.set_position([0.5, 0.5, 0.5, 0.5]) color_key(ax17, m_file, scan, colour_scheme) - ax18 = figs[17].add_subplot(211) + ax18 = figs[18].add_subplot(211) ax18.set_position([0.1, 0.33, 0.8, 0.6]) plot_radial_build(ax18, m_file, colour_scheme) # Make each axes smaller vertically to leave room for the legend - ax185 = figs[18].add_subplot(211) + ax185 = figs[19].add_subplot(211) ax185.set_position([0.1, 0.61, 0.8, 0.32]) - ax18b = figs[18].add_subplot(212) + ax18b = figs[19].add_subplot(212) ax18b.set_position([0.1, 0.13, 0.8, 0.32]) plot_upper_vertical_build(ax185, m_file, colour_scheme) plot_lower_vertical_build(ax18b, m_file, colour_scheme) @@ -13302,52 +13458,52 @@ def main_plot( # Can only plot WP and turn structure if superconducting coil at the moment if m_file.get("i_tf_sup", scan=scan) == 1: # TF coil with WP - ax19 = figs[19].add_subplot(221, aspect="equal") + ax19 = figs[20].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, scan, figs[19]) + plot_superconducting_tf_wp(ax19, m_file, scan, figs[20]) # TF coil turn structure - ax20 = figs[20].add_subplot(325, aspect="equal") + ax20 = figs[21].add_subplot(325, aspect="equal") ax20.set_position([0.025, 0.5, 0.4, 0.4]) - plot_tf_cable_in_conduit_turn(ax20, figs[20], m_file, scan) - plot_205 = figs[20].add_subplot(223, aspect="equal") + plot_tf_cable_in_conduit_turn(ax20, figs[21], m_file, scan) + plot_205 = figs[21].add_subplot(223, aspect="equal") plot_205.set_position([0.075, 0.1, 0.3, 0.3]) - plot_cable_in_conduit_cable(plot_205, figs[20], m_file, scan) + plot_cable_in_conduit_cable(plot_205, figs[21], m_file, scan) else: - ax19 = figs[19].add_subplot(211, aspect="equal") + ax19 = figs[20].add_subplot(211, aspect="equal") ax19.set_position([0.06, 0.55, 0.675, 0.4]) - plot_resistive_tf_wp(ax19, m_file, scan, figs[19]) - plot_resistive_tf_info(ax19, m_file, scan, figs[19]) + plot_resistive_tf_wp(ax19, m_file, scan, figs[20]) + plot_resistive_tf_info(ax19, m_file, scan, figs[20]) plot_tf_coil_structure( - figs[21].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme + figs[22].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme ) - plot_plasma_outboard_toroidal_ripple_map(figs[22], m_file, scan) + plot_plasma_outboard_toroidal_ripple_map(figs[23], m_file, scan) - plot_tf_stress(figs[23].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) + plot_tf_stress(figs[24].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) - plot_current_profiles_over_time(figs[24].add_subplot(111), m_file, scan) + plot_current_profiles_over_time(figs[25].add_subplot(111), m_file, scan) plot_cs_coil_structure( - figs[25].add_subplot(121, aspect="equal"), figs[25], m_file, scan + figs[26].add_subplot(121, aspect="equal"), figs[26], m_file, scan ) plot_cs_turn_structure( - figs[25].add_subplot(224, aspect="equal"), figs[25], m_file, scan + figs[26].add_subplot(224, aspect="equal"), figs[26], m_file, scan ) plot_first_wall_top_down_cross_section( - figs[26].add_subplot(221, aspect="equal"), m_file, scan + figs[27].add_subplot(221, aspect="equal"), m_file, scan ) - plot_first_wall_poloidal_cross_section(figs[26].add_subplot(122), m_file, scan) - plot_fw_90_deg_pipe_bend(figs[26].add_subplot(337), m_file, scan) + plot_first_wall_poloidal_cross_section(figs[27].add_subplot(122), m_file, scan) + plot_fw_90_deg_pipe_bend(figs[27].add_subplot(337), m_file, scan) - plot_blkt_pipe_bends(figs[27], m_file, scan) - ax_blanket = figs[27].add_subplot(122, aspect="equal") + plot_blkt_pipe_bends(figs[28], m_file, scan) + ax_blanket = figs[28].add_subplot(122, aspect="equal") plot_blanket(ax_blanket, m_file, scan, radial_build, colour_scheme) plot_firstwall(ax_blanket, m_file, scan, radial_build, colour_scheme) ax_blanket.set_xlabel("Radial position [m]") @@ -13390,13 +13546,13 @@ def main_plot( ) plot_main_power_flow( - figs[28].add_subplot(111, aspect="equal"), m_file, scan, figs[28] + figs[29].add_subplot(111, aspect="equal"), m_file, scan, figs[29] ) - ax24 = figs[29].add_subplot(111) + ax24 = figs[30].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, scan, figs[29]) + plot_system_power_profiles_over_time(ax24, m_file, scan, figs[30]) def create_thickness_builds(m_file, scan: int): @@ -13473,7 +13629,7 @@ def main(args=None): # create main plot # Increase range when adding new page - pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(30)] + pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(31)] # run main_plot main_plot( diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 91dd3bef0..d48afc783 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9921,6 +9921,22 @@ def run(self): ) ) + physics_variables.vel_plasma_alpha_thermal_profile = ( + self.calculate_relativistic_particle_speed( + e_kinetic=self.plasma_profile.teprofile.profile_y + * constants.KILOELECTRON_VOLT + * physics_variables.f_temp_plasma_ion_electron, + mass=constants.ALPHA_MASS, + ) + ) + + physics_variables.vel_plasma_alpha_birth = ( + self.calculate_relativistic_particle_speed( + e_kinetic=constants.DT_ALPHA_ENERGY, + mass=constants.ALPHA_MASS, + ) + ) + # ============================ # Plasma frequencies # ============================ @@ -9959,6 +9975,54 @@ def run(self): ) ) + # ============================ + # Upper hybrid frequencies + # ============================ + + physics_variables.freq_plasma_upper_hybrid_profile = self.calculate_upper_hybrid_frequency( + freq_plasma=np.concatenate([ + physics_variables.freq_plasma_electron_profile[::-1], + physics_variables.freq_plasma_electron_profile, + ]), + freq_larmor=physics_variables.freq_plasma_larmor_toroidal_electron_profile, + ) + + # ============================ + # Larmor radii + # ============================ + + # Since isotropic (v⟂)² = 2(v)² for a Maxwellian distribution, + # we can use the total velocity to calculate the Larmor radius for an isotropic profile + physics_variables.radius_plasma_deuteron_toroidal_larmor_isotropic_profile = self.calculate_larmor_radius( + vel_perp=np.sqrt( + 2 + * np.concatenate([ + physics_variables.vel_plasma_deuteron_profile[::-1], + physics_variables.vel_plasma_deuteron_profile, + ]) + ** 2 + ), + freq_larmor=physics_variables.freq_plasma_larmor_toroidal_deuteron_profile + * (2 * np.pi), + ) + + # Since isotropic (v⟂)² = 2(v)² for a Maxwellian distribution, + # we can use the total velocity to calculate the Larmor radius for an isotropic profile + physics_variables.radius_plasma_triton_toroidal_larmor_isotropic_profile = ( + self.calculate_larmor_radius( + vel_perp=np.sqrt( + 2 + * np.concatenate([ + physics_variables.vel_plasma_triton_profile[::-1], + physics_variables.vel_plasma_triton_profile, + ]) + ** 2 + ), + freq_larmor=physics_variables.freq_plasma_larmor_toroidal_triton_profile + * (2 * np.pi), + ) + ) + # ============================ # Coulomb logarithm # ============================ @@ -10056,6 +10120,72 @@ def run(self): for i in range(len(physics_variables.len_plasma_debye_electron_profile)) ]) + physics_variables.plasma_coulomb_log_deuteron_triton_profile = np.array([ + self.calculate_coulomb_log_from_impact( + impact_param_max=physics_variables.len_plasma_debye_electron_profile[i], + impact_param_min=max( + self.calculate_classical_distance_of_closest_approach( + charge1=1, + charge2=1, + m_reduced=self.calculate_reduced_mass( + mass1=constants.TRITON_MASS, + mass2=constants.DEUTERON_MASS, + ), + vel_relative=self.calculate_average_relative_velocity( + velocity_1=physics_variables.vel_plasma_triton_profile[i], + velocity_2=physics_variables.vel_plasma_deuteron_profile[i], + ), + ), + self.calculate_debroglie_wavelength( + mass=self.calculate_reduced_mass( + mass1=constants.TRITON_MASS, + mass2=constants.DEUTERON_MASS, + ), + velocity=self.calculate_average_relative_velocity( + velocity_1=physics_variables.vel_plasma_triton_profile[i], + velocity_2=physics_variables.vel_plasma_deuteron_profile[i], + ), + ), + ), + ) + for i in range(len(physics_variables.len_plasma_debye_electron_profile)) + ]) + + physics_variables.plasma_coulomb_log_electron_alpha_thermal_profile = np.array([ + self.calculate_coulomb_log_from_impact( + impact_param_max=physics_variables.len_plasma_debye_electron_profile[i], + impact_param_min=max( + self.calculate_classical_distance_of_closest_approach( + charge1=1, + charge2=2, + m_reduced=self.calculate_reduced_mass( + mass1=constants.ELECTRON_MASS, + mass2=constants.ALPHA_MASS, + ), + vel_relative=self.calculate_average_relative_velocity( + velocity_1=physics_variables.vel_plasma_electron_profile[i], + velocity_2=physics_variables.vel_plasma_alpha_thermal_profile[ + i + ], + ), + ), + self.calculate_debroglie_wavelength( + mass=self.calculate_reduced_mass( + mass1=constants.ELECTRON_MASS, + mass2=constants.ALPHA_MASS, + ), + velocity=self.calculate_average_relative_velocity( + velocity_1=physics_variables.vel_plasma_electron_profile[i], + velocity_2=physics_variables.vel_plasma_alpha_thermal_profile[ + i + ], + ), + ), + ), + ) + for i in range(len(physics_variables.len_plasma_debye_electron_profile)) + ]) + @staticmethod def calculate_debye_length( temp_plasma_species_kev: float | np.ndarray, @@ -10250,6 +10380,38 @@ def calculate_larmor_frequency( 2 * np.pi * m_particle ) + @staticmethod + def calculate_upper_hybrid_frequency( + freq_plasma: float | np.ndarray, freq_larmor: float | np.ndarray + ) -> float | np.ndarray: + """ + Calculate the upper hybrid frequency for a particle species. + + :param freq_plasma: Plasma frequency of the particle species (Hz). + :type freq_plasma: float | np.ndarray + :param freq_larmor: Larmor frequency of the particle species (Hz). + :type freq_larmor: float | np.ndarray + :returns: Upper hybrid frequency in Hz. + :rtype: float | np.ndarray + """ + return np.sqrt(freq_plasma**2 + freq_larmor**2) + + @staticmethod + def calculate_larmor_radius( + vel_perp: float | np.ndarray, + freq_larmor: float, + ) -> float | np.ndarray: + """ + Calculate the Larmor radius for a particle species. + :param vel_perp: Perpendicular velocity of the particle to the magnetic field (m/s). + :type vel_perp: float | np.ndarray + :param freq_larmor: Larmor frequency of the particle (Hz). + :type freq_larmor: float + :returns: Larmor radius in meters. + :rtype: float | np.ndarray + """ + return vel_perp / (freq_larmor) + @staticmethod def calculate_reduced_mass(mass1: float, mass2: float) -> float: """ @@ -10310,6 +10472,33 @@ def output_detailed_physics(self): physics_variables.len_plasma_debye_electron_profile[i], ) + po.osubhd(self.outfile, "Larmor radii:") + + for i in range( + len( + physics_variables.radius_plasma_deuteron_toroidal_larmor_isotropic_profile + ) + ): + po.ovarre( + self.mfile, + f"Plasma deuteron isotropic Larmor radius at point {i}", + f"(radius_plasma_deuteron_toroidal_larmor_isotropic_profile{i})", + physics_variables.radius_plasma_deuteron_toroidal_larmor_isotropic_profile[ + i + ], + ) + for i in range( + len(physics_variables.radius_plasma_triton_toroidal_larmor_isotropic_profile) + ): + po.ovarre( + self.mfile, + f"Plasma triton isotropic Larmor radius at point {i}", + f"(radius_plasma_triton_toroidal_larmor_isotropic_profile{i})", + physics_variables.radius_plasma_triton_toroidal_larmor_isotropic_profile[ + i + ], + ) + po.osubhd(self.outfile, "Velocities:") for i in range(len(physics_variables.vel_plasma_electron_profile)): @@ -10333,6 +10522,19 @@ def output_detailed_physics(self): f"(vel_plasma_triton_profile{i})", physics_variables.vel_plasma_triton_profile[i], ) + for i in range(len(physics_variables.vel_plasma_alpha_thermal_profile)): + po.ovarre( + self.mfile, + f"Plasma alpha thermal velocity at point {i}", + f"(vel_plasma_alpha_thermal_profile{i})", + physics_variables.vel_plasma_alpha_thermal_profile[i], + ) + po.ovarre( + self.outfile, + "Plasma alpha birth velocity (m/s)", + "(vel_plasma_alpha_birth)", + physics_variables.vel_plasma_alpha_birth, + ) po.osubhd(self.outfile, "Frequencies:") @@ -10371,6 +10573,14 @@ def output_detailed_physics(self): physics_variables.freq_plasma_larmor_toroidal_triton_profile[i], ) + for i in range(len(physics_variables.freq_plasma_upper_hybrid_profile)): + po.ovarre( + self.mfile, + f"Plasma upper hybrid frequency at point {i}", + f"(freq_plasma_upper_hybrid_profile{i})", + physics_variables.freq_plasma_upper_hybrid_profile[i], + ) + po.osubhd(self.outfile, "Coulomb Logarithms:") for i in range( @@ -10402,3 +10612,23 @@ def output_detailed_physics(self): f"(plasma_coulomb_log_electron_triton_profile{i})", physics_variables.plasma_coulomb_log_electron_triton_profile[i], ) + + for i in range( + len(physics_variables.plasma_coulomb_log_deuteron_triton_profile) + ): + po.ovarre( + self.mfile, + f"Deuteron-triton Coulomb log at point {i}", + f"(plasma_coulomb_log_deuteron_triton_profile{i})", + physics_variables.plasma_coulomb_log_deuteron_triton_profile[i], + ) + + for i in range( + len(physics_variables.plasma_coulomb_log_electron_alpha_thermal_profile) + ): + po.ovarre( + self.mfile, + f"Electron-alpha thermal Coulomb log at point {i}", + f"(plasma_coulomb_log_electron_alpha_thermal_profile{i})", + physics_variables.plasma_coulomb_log_electron_alpha_thermal_profile[i], + ) diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 93273518b..b09401ccb 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -3617,7 +3617,7 @@ def test_detailed_physics_run_computes_profiles(): ) # toroidal field profile for larmor frequency calc physics_variables.b_plasma_toroidal_profile = ( - np.ones_like(plasma.neprofile.profile_y) * 5.0 + np.ones(2 * len(plasma.neprofile.profile_y)) * 5.0 ) dp = DetailedPhysics(plasma) @@ -3637,7 +3637,8 @@ def test_detailed_physics_run_computes_profiles(): assert np.all(np.isfinite(physics_variables.freq_plasma_electron_profile)) assert ( - np.shape(physics_variables.freq_plasma_larmor_toroidal_electron_profile)[0] == n + np.shape(physics_variables.freq_plasma_larmor_toroidal_electron_profile)[0] + == n * 2 ) assert np.all( np.isfinite(physics_variables.freq_plasma_larmor_toroidal_electron_profile)