diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index ce1e9c99f..bd3724da9 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1314,15 +1314,33 @@ vel_plasma_electron_profile: list[float] = None """Profile of electron thermal velocity in plasma (m/s)""" +vel_plasma_deuteron_profile: list[float] = None +"""Profile of deuteron thermal velocity in plasma (m/s)""" + +vel_plasma_triton_profile: list[float] = None +"""Profile of triton thermal velocity in plasma (m/s)""" + plasma_coulomb_log_electron_electron_profile: list[float] = None """Profile of electron-electron Coulomb logarithm in plasma""" +plasma_coulomb_log_electron_deuteron_profile: list[float] = None +"""Profile of electron-deuteron Coulomb logarithm in plasma""" + +plasma_coulomb_log_electron_triton_profile: list[float] = None +"""Profile of electron-triton Coulomb logarithm in plasma""" + freq_plasma_electron_profile: list[float] = None """Electron 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)""" +freq_plasma_larmor_toroidal_deuteron_profile: list[float] = None +"""Profile of deuteron Larmor frequency in plasma due to toroidal magnetic field (Hz)""" + +freq_plasma_larmor_toroidal_triton_profile: list[float] = None +"""Profile of triton Larmor frequency in plasma due to toroidal magnetic field (Hz)""" + def init_physics_module(): """Initialise the physics module""" @@ -1643,9 +1661,15 @@ def init_physics_variables(): len_plasma_debye_electron_profile, \ len_plasma_debye_electron_vol_avg, \ vel_plasma_electron_profile, \ + vel_plasma_deuteron_profile, \ + vel_plasma_triton_profile, \ plasma_coulomb_log_electron_electron_profile, \ + plasma_coulomb_log_electron_deuteron_profile, \ + plasma_coulomb_log_electron_triton_profile, \ freq_plasma_electron_profile, \ - freq_plasma_larmor_toroidal_electron_profile + freq_plasma_larmor_toroidal_electron_profile, \ + freq_plasma_larmor_toroidal_deuteron_profile, \ + freq_plasma_larmor_toroidal_triton_profile m_beam_amu = 0.0 m_fuel_amu = 0.0 @@ -1906,6 +1930,12 @@ def init_physics_variables(): len_plasma_debye_electron_profile = [] len_plasma_debye_electron_vol_avg = 0.0 vel_plasma_electron_profile = [] + vel_plasma_deuteron_profile = [] + vel_plasma_triton_profile = [] plasma_coulomb_log_electron_electron_profile = [] + plasma_coulomb_log_electron_deuteron_profile = [] + plasma_coulomb_log_electron_triton_profile = [] freq_plasma_electron_profile = [] freq_plasma_larmor_toroidal_electron_profile = [] + freq_plasma_larmor_toroidal_deuteron_profile = [] + freq_plasma_larmor_toroidal_triton_profile = [] diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 26d2a1e98..67ab13a49 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -12389,7 +12389,7 @@ def plot_ebw_ecrh_coupling_graph(axis: plt.Axes, mfile: mf.MFile, scan: int): axis.minorticks_on() -def plot_debye_length_profile(axis, mfile_data, scan): +def plot_debye_length_profile(axis: plt.Axes, mfile_data: mf.MFile, scan: int): """Plot the Debye length profile on the given axis.""" len_plasma_debye_electron_profile = [ mfile_data.data[f"len_plasma_debye_electron_profile{i}"].get_scan(scan) @@ -12424,6 +12424,14 @@ def plot_velocity_profile(axis, mfile_data, scan): mfile_data.data[f"vel_plasma_electron_profile{i}"].get_scan(scan) for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) ] + vel_plasma_deuteron_profile = [ + mfile_data.data[f"vel_plasma_deuteron_profile{i}"].get_scan(scan) + for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) + ] + vel_plasma_triton_profile = [ + 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))) + ] axis.plot( np.linspace(0, 1, len(vel_plasma_electron_profile)), @@ -12432,7 +12440,22 @@ def plot_velocity_profile(axis, mfile_data, scan): linestyle="-", label=r"$v_{e}$", ) + axis.plot( + np.linspace(0, 1, len(vel_plasma_deuteron_profile)), + vel_plasma_deuteron_profile, + color="red", + linestyle="-", + label=r"$v_{D}$", + ) + axis.plot( + np.linspace(0, 1, len(vel_plasma_triton_profile)), + vel_plasma_triton_profile, + color="green", + linestyle="-", + label=r"$v_{T}$", + ) + axis.set_yscale("log") axis.set_ylabel("Velocity [m/s]") axis.set_xlabel("$\\rho \\ [r/a]$") axis.grid(True, which="both", linestyle="--", alpha=0.5) @@ -12441,7 +12464,7 @@ def plot_velocity_profile(axis, mfile_data, scan): axis.legend() -def plot_frequency_profile(axis, mfile_data, scan): +def plot_electron_frequency_profile(axis, mfile_data, scan): """Plot the electron thermal frequency profile on the given axis.""" freq_plasma_electron_profile = [ mfile_data.data[f"freq_plasma_electron_profile{i}"].get_scan(scan) @@ -12472,6 +12495,44 @@ def plot_frequency_profile(axis, mfile_data, scan): axis.set_xlim(-1.025, 1.025) axis.set_ylabel("Frequency [GHz]") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + + axis.legend() + + +def plot_ion_frequency_profile(axis, mfile_data, scan): + freq_plasma_larmor_toroidal_deuteron_profile = [ + mfile_data.data[f"freq_plasma_larmor_toroidal_deuteron_profile{i}"].get_scan( + scan + ) + for i in range( + 2 * int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + freq_plasma_larmor_toroidal_triton_profile = [ + mfile_data.data[f"freq_plasma_larmor_toroidal_triton_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_deuteron_profile)), + np.array(freq_plasma_larmor_toroidal_deuteron_profile) / 1e6, + color="red", + linestyle="-", + label=r"$f_{Larmor,toroidal,D}$", + ) + axis.plot( + np.linspace(-1, 1, len(freq_plasma_larmor_toroidal_triton_profile)), + np.array(freq_plasma_larmor_toroidal_triton_profile) / 1e6, + color="green", + linestyle="-", + label=r"$f_{Larmor,toroidal,T}$", + ) + + axis.set_ylabel("Frequency [MHz]") axis.set_xlabel("$\\rho \\ [r/a]$") axis.grid(True, which="both", linestyle="--", alpha=0.5) axis.minorticks_on() @@ -12487,6 +12548,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_electron_deuteron_profile = [ + mfile_data.data[f"plasma_coulomb_log_electron_deuteron_profile{i}"].get_scan( + scan + ) + for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) + ] + + plasma_coulomb_log_electron_triton_profile = [ + mfile_data.data[f"plasma_coulomb_log_electron_triton_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, @@ -12495,6 +12568,22 @@ def plot_plasma_coloumb_logarithms(axis, mfile_data, scan): label=r"$ln \Lambda_{e-e}$", ) + axis.plot( + np.linspace(0, 1, len(plasma_coulomb_log_electron_deuteron_profile)), + plasma_coulomb_log_electron_deuteron_profile, + color="red", + linestyle="-", + label=r"$ln \Lambda_{e-D}$", + ) + + axis.plot( + np.linspace(0, 1, len(plasma_coulomb_log_electron_triton_profile)), + plasma_coulomb_log_electron_triton_profile, + color="green", + linestyle="-", + label=r"$ln \Lambda_{e-T}$", + ) + axis.set_ylabel("Coulomb Logarithm") axis.set_xlabel("$\\rho \\ [r/a]$") axis.grid(True, which="both", linestyle="--", alpha=0.5) @@ -12924,7 +13013,14 @@ 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_frequency_profile(figs[15].add_subplot(212), 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_ion_frequency_profile(ax_ion, m_file, scan) + plot_plasma_coloumb_logarithms(figs[15].add_subplot(231), m_file, scan) # Plot poloidal cross-section diff --git a/process/physics.py b/process/physics.py index 94d7eeefb..9d803888b 100644 --- a/process/physics.py +++ b/process/physics.py @@ -9280,6 +9280,24 @@ def run(self): ) ) + physics_variables.vel_plasma_deuteron_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.DEUTERON_MASS, + ) + ) + + physics_variables.vel_plasma_triton_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.TRITON_MASS, + ) + ) + # ============================ # Plasma frequencies # ============================ @@ -9302,6 +9320,22 @@ def run(self): ) ) + physics_variables.freq_plasma_larmor_toroidal_deuteron_profile = ( + self.calculate_larmor_frequency( + b_field=physics_variables.b_plasma_toroidal_profile, + m_particle=constants.DEUTERON_MASS, + z_particle=1.0, + ) + ) + + physics_variables.freq_plasma_larmor_toroidal_triton_profile = ( + self.calculate_larmor_frequency( + b_field=physics_variables.b_plasma_toroidal_profile, + m_particle=constants.TRITON_MASS, + z_particle=1.0, + ) + ) + # ============================ # Coulomb logarithm # ============================ @@ -9313,12 +9347,86 @@ def run(self): self.calculate_classical_distance_of_closest_approach( charge1=1, charge2=1, - e_kinetic=self.plasma_profile.teprofile.profile_y[i] - * constants.KILOELECTRON_VOLT, + m_reduced=self.calculate_reduced_mass( + mass1=constants.ELECTRON_MASS, + mass2=constants.ELECTRON_MASS, + ), + vel_relative=self.calculate_average_relative_velocity( + velocity_1=physics_variables.vel_plasma_electron_profile[i], + velocity_2=physics_variables.vel_plasma_electron_profile[i], + ), + ), + self.calculate_debroglie_wavelength( + mass=self.calculate_reduced_mass( + mass1=constants.ELECTRON_MASS, + mass2=constants.ELECTRON_MASS, + ), + velocity=self.calculate_average_relative_velocity( + velocity_1=physics_variables.vel_plasma_electron_profile[i], + velocity_2=physics_variables.vel_plasma_electron_profile[i], + ), + ), + ), + ) + for i in range(len(physics_variables.len_plasma_debye_electron_profile)) + ]) + + physics_variables.plasma_coulomb_log_electron_deuteron_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.DEUTERON_MASS, + mass2=constants.ELECTRON_MASS, + ), + vel_relative=self.calculate_average_relative_velocity( + velocity_1=physics_variables.vel_plasma_deuteron_profile[i], + velocity_2=physics_variables.vel_plasma_electron_profile[i], + ), + ), + self.calculate_debroglie_wavelength( + mass=self.calculate_reduced_mass( + mass1=constants.DEUTERON_MASS, + mass2=constants.ELECTRON_MASS, + ), + velocity=self.calculate_average_relative_velocity( + velocity_1=physics_variables.vel_plasma_deuteron_profile[i], + velocity_2=physics_variables.vel_plasma_electron_profile[i], + ), + ), + ), + ) + for i in range(len(physics_variables.len_plasma_debye_electron_profile)) + ]) + + physics_variables.plasma_coulomb_log_electron_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.ELECTRON_MASS, + ), + vel_relative=self.calculate_average_relative_velocity( + velocity_1=physics_variables.vel_plasma_triton_profile[i], + velocity_2=physics_variables.vel_plasma_electron_profile[i], + ), ), self.calculate_debroglie_wavelength( - mass=constants.ELECTRON_MASS, - velocity=physics_variables.vel_plasma_electron_profile[i], + mass=self.calculate_reduced_mass( + mass1=constants.TRITON_MASS, + mass2=constants.ELECTRON_MASS, + ), + velocity=self.calculate_average_relative_velocity( + velocity_1=physics_variables.vel_plasma_triton_profile[i], + velocity_2=physics_variables.vel_plasma_electron_profile[i], + ), ), ), ) @@ -9394,7 +9502,8 @@ def calculate_coulomb_log_from_impact( def calculate_classical_distance_of_closest_approach( charge1: float, charge2: float, - e_kinetic: float | np.ndarray, + m_reduced: float, + vel_relative: float | np.ndarray, ) -> float | np.ndarray: """ Calculate the classical distance of closest approach for two charged particles. @@ -9403,14 +9512,16 @@ def calculate_classical_distance_of_closest_approach( :type charge1: float :param charge2: Charge of particle 2 in units of elementary charge. :type charge2: float - :param e_kinetic: Kinetic energy of the particles in Joules. - :type e_kinetic: float | np.ndarray + :param m_reduced: Reduced mass of the two-particle system in kg. + :type m_reduced: float + :param vel_relative: Relative velocity of the two particles in m/s. + :type vel_relative: float | np.ndarray :returns: Distance of closest approach in meters. :rtype: float | np.ndarray """ return (charge1 * charge2 * constants.ELECTRON_CHARGE**2) / ( - 4 * np.pi * constants.EPSILON0 * e_kinetic + 2 * np.pi * constants.EPSILON0 * m_reduced * vel_relative**2 ) @staticmethod @@ -9472,6 +9583,34 @@ def calculate_larmor_frequency( 2 * np.pi * m_particle ) + @staticmethod + def calculate_reduced_mass(mass1: float, mass2: float) -> float: + """ + Calculate the reduced mass of two particles. + :param mass1: Mass of particle 1 (kg). + :type mass1: float + :param mass2: Mass of particle 2 (kg). + :type mass2: float + :returns: Reduced mass (kg). + :rtype: float + """ + return (mass1 * mass2) / (mass1 + mass2) + + @staticmethod + def calculate_average_relative_velocity( + velocity_1: float | np.ndarray, velocity_2: float | np.ndarray + ) -> float | np.ndarray: + """ + Calculate the average relative velocity between two particles. + :param velocity_1: Velocity of particle 1 (m/s). + :type velocity_1: float | np.ndarray + :param velocity_2: Velocity of particle 2 (m/s). + :type velocity_2: float | np.ndarray + :returns: Average relative velocity (m/s). + :rtype: float | np.ndarray + """ + return np.sqrt(velocity_1**2 + velocity_2**2) + def output_detailed_physics(self): """Outputs detailed physics variables to file.""" @@ -9503,6 +9642,20 @@ def output_detailed_physics(self): f"(vel_plasma_electron_profile{i})", physics_variables.vel_plasma_electron_profile[i], ) + for i in range(len(physics_variables.vel_plasma_deuteron_profile)): + po.ovarre( + self.mfile, + f"Plasma deuteron thermal velocity at point {i}", + f"(vel_plasma_deuteron_profile{i})", + physics_variables.vel_plasma_deuteron_profile[i], + ) + for i in range(len(physics_variables.vel_plasma_triton_profile)): + po.ovarre( + self.mfile, + f"Plasma triton thermal velocity at point {i}", + f"(vel_plasma_triton_profile{i})", + physics_variables.vel_plasma_triton_profile[i], + ) po.osubhd(self.outfile, "Frequencies:") @@ -9518,10 +9671,28 @@ def output_detailed_physics(self): ): po.ovarre( self.mfile, - f"Plasma electron Larmor frequency at point {i}", + f"Plasma electron toroidal Larmor frequency at point {i}", f"(freq_plasma_larmor_toroidal_electron_profile{i})", physics_variables.freq_plasma_larmor_toroidal_electron_profile[i], ) + for i in range( + len(physics_variables.freq_plasma_larmor_toroidal_deuteron_profile) + ): + po.ovarre( + self.mfile, + f"Plasma deuteron toroidal Larmor frequency at point {i}", + f"(freq_plasma_larmor_toroidal_deuteron_profile{i})", + physics_variables.freq_plasma_larmor_toroidal_deuteron_profile[i], + ) + for i in range( + len(physics_variables.freq_plasma_larmor_toroidal_triton_profile) + ): + po.ovarre( + self.mfile, + f"Plasma triton toroidal Larmor frequency at point {i}", + f"(freq_plasma_larmor_toroidal_triton_profile{i})", + physics_variables.freq_plasma_larmor_toroidal_triton_profile[i], + ) po.osubhd(self.outfile, "Coulomb Logarithms:") @@ -9534,3 +9705,23 @@ def output_detailed_physics(self): f"(plasma_coulomb_log_electron_electron_profile{i})", physics_variables.plasma_coulomb_log_electron_electron_profile[i], ) + + for i in range( + len(physics_variables.plasma_coulomb_log_electron_deuteron_profile) + ): + po.ovarre( + self.mfile, + f"Electron-deuteron Coulomb log at point {i}", + f"(plasma_coulomb_log_electron_deuteron_profile{i})", + physics_variables.plasma_coulomb_log_electron_deuteron_profile[i], + ) + + for i in range( + len(physics_variables.plasma_coulomb_log_electron_triton_profile) + ): + po.ovarre( + self.mfile, + f"Electron-triton Coulomb log at point {i}", + f"(plasma_coulomb_log_electron_triton_profile{i})", + physics_variables.plasma_coulomb_log_electron_triton_profile[i], + )