diff --git a/documentation/physics-models/detailed_physics.md b/documentation/physics-models/detailed_physics.md new file mode 100644 index 0000000000..dd01a72ffd --- /dev/null +++ b/documentation/physics-models/detailed_physics.md @@ -0,0 +1,85 @@ +# Detailed Plasma Physics + +It can sometimes be useful to calculate rough values for key plasma paramters that are normally used in higher fidelity codes. The `DetailedPhysics()` class stores functions that are called and the end of the run to show rough values for key plasma behavior parameters. The calculation is done at the end as no other methods currently depend on these values. + +## Detailed Plasma Physics | `DetailedPhysics()` + + +------------------ + +### Debye length | `calculate_debye_length()` + +Calculates the Debye lenght given by: + +$$ +\lambda_{D} = \sqrt{\frac{\epsilon_0 k_B T_e}{n e^2}} +$$ + +------------------- + +### Relativistic particle speed | `calculate_relativistic_particle_speed()` + +$$ +v = c \times \sqrt{\left(1- \frac{1}{\left(1+\frac{E}{mc^2}\right)^2}\right)} +$$ + +------------------ + +### Coulomb Logarithm | `calculate_coulomb_log_from_impact()` + +Calculates the Coulomb logarithm assuming a straight line Landau-Spitzer method + +$$ +\ln \Lambda = \ln{\left(\frac{b_{\text{max}}}{b_{\text{min}}}\right)} +$$ + +The maximum impact parameter is given by the Debye length calculated by [`calculate_debye_length()`](#debye-length--calculate_debye_length) +$$ +b_{\text{max}} = \lambda_{\text{Debye}} +$$ + +The minimum impact paramter is the largest of either the classical distance of closest approach or the Debye length. + +$$ +\begin{split}b_{\text{min}} ≡ +\left\{ + \begin{array}{ll} + λ_{\text{de Broglie}} & \mbox{if } λ_{\text{de Broglie}} ≥ ρ_⟂ \\ + ρ_⟂ & \mbox{if } ρ_⟂ ≥ λ_{\text{de Broglie}} + \end{array} +\right.\end{split} +$$ + +$ρ_⟂$ is the classical distance of closest approach calculated by [`calculate_classical_distance_of_closest_approach()`](#classical-distance-of-closest-approach----calculate_classical_distance_of_closest_approach) + +------------------ + +### Classical distance of closest approach | `calculate_classical_distance_of_closest_approach()` + +$$ +\frac{Z_1Z_2e^2}{4\pi \epsilon_0 E_{\text{kinetic}}} +$$ + +--------------------- + +### DeBroglie Wavelength | `calculate_debroglie_wavelength()` + +$$ +\lambda_{\text{DeBroglie}} = \frac{h}{2\pi m v} +$$ + +---------------------- + +### Plasma Frequency | `calculate_plasma_frequency()` + +$$ +\omega_p = \sqrt{\frac{n_ie^2}{\epsilon_0 m_i}} +$$ + +--------------------- + +### Larmor Frequency | `calculate_larmor_frequency()` + +$$ +f_{\text{Larmor}} = \frac{Z_ieB}{2\pi m_i} +$$ \ No newline at end of file diff --git a/mkdocs.yml b/mkdocs.yml index 22dcfea2dc..ba7cb96aa7 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -70,6 +70,7 @@ nav: - Confinement time: physics-models/plasma_confinement.md - L-H transition: physics-models/plasma_h_mode.md - Plasma Core Power Balance: physics-models/plasma_power_balance.md + - Detailed Plasma Physics: physics-models/detailed_physics.md - Pulsed Plant Operation: physics-models/pulsed-plant.md - Engineering Models: - Machine Build: eng-models/machine-build.md diff --git a/process/constants.py b/process/constants.py index 98025ec20b..100fe4d657 100644 --- a/process/constants.py +++ b/process/constants.py @@ -189,12 +189,15 @@ https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=W """ -SPEED_LIGHT = 299792458 +SPEED_LIGHT = 299792458.0 """Speed of light in vacuum (c) [m/s] Reference: National Institute of Standards and Technology (NIST) https://physics.nist.gov/cgi-bin/cuu/Value?c|search_for=light """ +PLANCK_CONSTANT = 6.62607015e-34 +"""Planck's constant [J.s]""" + D_T_ENERGY = ( (DEUTERON_MASS + TRITON_MASS) - (ALPHA_MASS + NEUTRON_MASS) ) * SPEED_LIGHT**2 diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index e14f132539..e7fc6ab59d 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1320,6 +1320,24 @@ n_charge_plasma_effective_mass_weighted_vol_avg: float = None """Plasma mass-weighted volume averaged plasma effective charge""" +len_plasma_debye_electron_profile: list[float] = None +"""Profile of electron Debye length in plasma (m)""" + +len_plasma_debye_electron_vol_avg: float = None +"""Volume averaged electron Debye length in plasma (m)""" + +vel_plasma_electron_profile: list[float] = None +"""Profile of electron thermal velocity in plasma (m/s)""" + +plasma_coulomb_log_electron_electron_profile: list[float] = None +"""Profile of electron-electron 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)""" + def init_physics_module(): """Initialise the physics module""" @@ -1639,7 +1657,13 @@ def init_physics_variables(): a_plasma_poloidal, \ n_charge_plasma_effective_vol_avg, \ n_charge_plasma_effective_profile, \ - n_charge_plasma_effective_mass_weighted_vol_avg + n_charge_plasma_effective_mass_weighted_vol_avg, \ + len_plasma_debye_electron_profile, \ + len_plasma_debye_electron_vol_avg, \ + vel_plasma_electron_profile, \ + plasma_coulomb_log_electron_electron_profile, \ + freq_plasma_electron_profile, \ + freq_plasma_larmor_toroidal_electron_profile m_beam_amu = 0.0 m_fuel_amu = 0.0 @@ -1900,3 +1924,9 @@ def init_physics_variables(): n_charge_plasma_effective_vol_avg = 0.0 n_charge_plasma_effective_profile = [] n_charge_plasma_effective_mass_weighted_vol_avg = 0.0 + len_plasma_debye_electron_profile = [] + len_plasma_debye_electron_vol_avg = 0.0 + vel_plasma_electron_profile = [] + plasma_coulomb_log_electron_electron_profile = [] + freq_plasma_electron_profile = [] + freq_plasma_larmor_toroidal_electron_profile = [] diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 7d175b1741..b64fd5ea91 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -12269,6 +12269,119 @@ 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): + """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) + for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) + ] + + # Convert to micrometres (1e-6 m) + len_plasma_debye_electron_profile_um = [ + length * 1e6 for length in len_plasma_debye_electron_profile + ] + + axis.plot( + np.linspace(0, 1, len(len_plasma_debye_electron_profile_um)), + len_plasma_debye_electron_profile_um, + color="blue", + linestyle="-", + label=r"$\lambda_{Debye,e}$", + ) + + axis.set_ylabel(r"Debye Length [$\mu$m]") + + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.set_xlim([0, 1.025]) + axis.minorticks_on() + axis.legend() + + +def plot_velocity_profile(axis, mfile_data, scan): + """Plot the electron thermal velocity profile on the given axis.""" + vel_plasma_electron_profile = [ + 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))) + ] + + axis.plot( + np.linspace(0, 1, len(vel_plasma_electron_profile)), + vel_plasma_electron_profile, + color="blue", + linestyle="-", + label=r"$v_{e}$", + ) + + axis.set_ylabel("Velocity [m/s]") + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.set_xlim([0, 1.025]) + axis.minorticks_on() + axis.legend() + + +def plot_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) + for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) + ] + freq_plasma_larmor_toroidal_electron_profile = [ + mfile_data.data[f"freq_plasma_larmor_toroidal_electron_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}$", + ) + 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}$") + # mirrored across the y-axis (drawn at negative rho) + axis.plot(-x, y, color="blue", linestyle="-", label="_nolegend_") + axis.set_xlim(-1.025, 1.025) + + axis.set_ylabel("Frequency [GHz]") + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.minorticks_on() + axis.legend() + + +def plot_plasma_coloumb_logarithms(axis, mfile_data, scan): + """Plot the plasma coloumb logarithms on the given axis.""" + plasma_coulomb_log_electron_electron_profile = [ + mfile_data.data[f"plasma_coulomb_log_electron_electron_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, + color="blue", + linestyle="-", + label=r"$ln \Lambda_{e-e}$", + ) + + axis.set_ylabel("Coulomb Logarithm") + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.minorticks_on() + axis.legend() + + def main_plot( figs: list[Axes], m_file: mf.MFile, @@ -12440,9 +12553,14 @@ def main_plot( plot_density_limit_comparison(figs[13].add_subplot(221), m_file, scan) plot_confinement_time_comparison(figs[13].add_subplot(224), m_file, scan) + plot_debye_length_profile(figs[14].add_subplot(232), m_file, scan) + plot_velocity_profile(figs[14].add_subplot(233), m_file, scan) + plot_frequency_profile(figs[14].add_subplot(212), m_file, scan) + plot_plasma_coloumb_logarithms(figs[14].add_subplot(231), m_file, scan) + # Plot poloidal cross-section poloidal_cross_section( - figs[14].add_subplot(121, aspect="equal"), + figs[15].add_subplot(121, aspect="equal"), m_file, scan, demo_ranges, @@ -12452,7 +12570,7 @@ def main_plot( # Plot toroidal cross-section toroidal_cross_section( - figs[14].add_subplot(122, aspect="equal"), + figs[15].add_subplot(122, aspect="equal"), m_file, scan, demo_ranges, @@ -12460,19 +12578,19 @@ def main_plot( ) # Plot color key - ax17 = figs[14].add_subplot(222) + ax17 = figs[15].add_subplot(222) ax17.set_position([0.5, 0.5, 0.5, 0.5]) color_key(ax17, m_file, scan, colour_scheme) - ax18 = figs[15].add_subplot(211) + ax18 = figs[16].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[16].add_subplot(211) + ax185 = figs[17].add_subplot(211) ax185.set_position([0.1, 0.61, 0.8, 0.32]) - ax18b = figs[16].add_subplot(212) + ax18b = figs[17].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) @@ -12480,48 +12598,51 @@ 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[17].add_subplot(221, aspect="equal") - # Half height, a bit wider, top left - ax19.set_position([0.025, 0.45, 0.5, 0.5]) - plot_superconducting_tf_wp(ax19, m_file, scan, figs[17]) + ax19 = figs[18].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[18]) # TF coil turn structure - ax20 = figs[18].add_subplot(325, aspect="equal") + ax20 = figs[19].add_subplot(325, aspect="equal") ax20.set_position([0.025, 0.5, 0.4, 0.4]) - plot_tf_cable_in_conduit_turn(ax20, figs[18], m_file, scan) - plot_205 = figs[18].add_subplot(223, aspect="equal") + plot_tf_cable_in_conduit_turn(ax20, figs[19], m_file, scan) + plot_205 = figs[19].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[18], m_file, scan) + plot_cable_in_conduit_cable(plot_205, figs[19], m_file, scan) else: - ax19 = figs[17].add_subplot(211, aspect="equal") + ax19 = figs[18].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[17]) - + plot_resistive_tf_wp(ax19, m_file, scan, figs[18]) plot_tf_coil_structure( - figs[19].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme + figs[20].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme ) - plot_plasma_outboard_toroidal_ripple_map(figs[20], m_file, scan) + plot_plasma_outboard_toroidal_ripple_map(figs[21], m_file, scan) - plot_tf_stress(figs[21].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) + plot_tf_stress(figs[22].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) - plot_current_profiles_over_time(figs[22].add_subplot(111), m_file, scan) + plot_current_profiles_over_time(figs[23].add_subplot(111), m_file, scan) plot_cs_coil_structure( - figs[23].add_subplot(121, aspect="equal"), figs[23], m_file, scan + figs[24].add_subplot(121, aspect="equal"), figs[24], m_file, scan ) plot_cs_turn_structure( - figs[23].add_subplot(224, aspect="equal"), figs[23], m_file, scan + figs[24].add_subplot(224, aspect="equal"), figs[24], m_file, scan ) plot_first_wall_top_down_cross_section( - figs[24].add_subplot(221, aspect="equal"), m_file, scan + figs[25].add_subplot(221, aspect="equal"), m_file, scan ) - plot_first_wall_poloidal_cross_section(figs[24].add_subplot(122), m_file, scan) - plot_fw_90_deg_pipe_bend(figs[24].add_subplot(337), m_file, scan) + plot_first_wall_poloidal_cross_section(figs[25].add_subplot(122), m_file, scan) + plot_fw_90_deg_pipe_bend(figs[25].add_subplot(337), m_file, scan) - plot_blkt_pipe_bends(figs[25], m_file, scan) - ax_blanket = figs[25].add_subplot(122, aspect="equal") + plot_blkt_pipe_bends(figs[26], m_file, scan) + ax_blanket = figs[26].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]") @@ -12564,13 +12685,13 @@ def main_plot( ) plot_main_power_flow( - figs[26].add_subplot(111, aspect="equal"), m_file, scan, figs[26] + figs[27].add_subplot(111, aspect="equal"), m_file, scan, figs[27] ) - ax24 = figs[27].add_subplot(111) + ax24 = figs[28].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[27]) + plot_system_power_profiles_over_time(ax24, m_file, scan, figs[28]) def create_thickness_builds(m_file, scan: int): @@ -12647,7 +12768,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(28)] + pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(29)] # run main_plot main_plot( diff --git a/process/main.py b/process/main.py index e16b7f2109..97c074e6b1 100644 --- a/process/main.py +++ b/process/main.py @@ -94,7 +94,7 @@ ) from process.log import logging_model_handler, show_errors from process.pfcoil import PFCoil -from process.physics import Physics +from process.physics import DetailedPhysics, Physics from process.plasma_geometry import PlasmaGeom from process.plasma_profiles import PlasmaProfile from process.power import Power @@ -678,6 +678,9 @@ def __init__(self): self.physics = Physics( plasma_profile=self.plasma_profile, current_drive=self.current_drive ) + self.physics_detailed = DetailedPhysics( + plasma_profile=self.plasma_profile, + ) self.neoclassics = Neoclassics() self.stellarator = Stellarator( availability=self.availability, diff --git a/process/output.py b/process/output.py index 5cd8422f73..1fe44b8281 100644 --- a/process/output.py +++ b/process/output.py @@ -44,6 +44,11 @@ def write(models, _outfile): models.physics.calculate_effective_charge_ionisation_profiles() models.physics.outplas() + # Detailed physics, currently only done at final point as values are not used + # by any other functions + models.physics_detailed.run() + models.physics_detailed.output_detailed_physics() + # TODO what is this? Not in caller.f90? models.current_drive.output_current_drive() diff --git a/process/physics.py b/process/physics.py index 727859b55b..82f3963961 100644 --- a/process/physics.py +++ b/process/physics.py @@ -9046,3 +9046,296 @@ def reinke_tsep(b_plasma_toroidal_on_axis, flh, qstar, rmajor, eps, fgw, kappa, * (eps**0.15 * (1.0 + kappa**2.0) ** 0.34) * (lhat**0.29 * kappa_0 ** (-0.29) * 0.285) ) + + +class DetailedPhysics: + """Class to hold detailed physics models for plasma processing.""" + + def __init__(self, plasma_profile): + self.outfile = constants.NOUT + self.mfile = constants.MFILE + self.plasma_profile = plasma_profile + + def run(self): + # --------------------------- + # Debye length calculation + # --------------------------- + + physics_variables.len_plasma_debye_electron_vol_avg = self.calculate_debye_length( + temp_plasma_species_kev=physics_variables.temp_plasma_electron_vol_avg_kev, + nd_plasma_species=physics_variables.nd_plasma_electrons_vol_avg, + ) + + physics_variables.len_plasma_debye_electron_profile = ( + self.calculate_debye_length( + temp_plasma_species_kev=self.plasma_profile.teprofile.profile_y, + nd_plasma_species=self.plasma_profile.neprofile.profile_y, + ) + ) + + # ============================ + # Particle relativistic speeds + # ============================ + + physics_variables.vel_plasma_electron_profile = ( + self.calculate_relativistic_particle_speed( + e_kinetic=self.plasma_profile.teprofile.profile_y + * constants.KILOELECTRON_VOLT, + mass=constants.ELECTRON_MASS, + ) + ) + + # ============================ + # Plasma frequencies + # ============================ + + physics_variables.freq_plasma_electron_profile = self.calculate_plasma_frequency( + nd_particle=self.plasma_profile.neprofile.profile_y, + m_particle=constants.ELECTRON_MASS, + z_particle=1.0, + ) + + # ============================ + # Larmor frequencies + # ============================ + + physics_variables.freq_plasma_larmor_toroidal_electron_profile = ( + self.calculate_larmor_frequency( + b_field=physics_variables.b_plasma_toroidal_profile, + m_particle=constants.ELECTRON_MASS, + z_particle=1.0, + ) + ) + + # ============================ + # Coulomb logarithm + # ============================ + + physics_variables.plasma_coulomb_log_electron_electron_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, + e_kinetic=self.plasma_profile.teprofile.profile_y[i] + * constants.KILOELECTRON_VOLT, + ), + self.calculate_debroglie_wavelength( + mass=constants.ELECTRON_MASS, + velocity=physics_variables.vel_plasma_electron_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, + nd_plasma_species: float | np.ndarray, + ) -> float | np.ndarray: + """ + Calculate the Debye length for a plasma. + + :param temp_plasma_species_kev: Species temperature in keV. + :type temp_plasma_species_kev: float | np.ndarray + :param nd_plasma_species: Species number density (/m^3). + :type nd_plasma_species: float | np.ndarray + + :returns: Debye length in meters. + :rtype: float | np.ndarray + """ + return ( + (constants.EPSILON0 * temp_plasma_species_kev * constants.KILOELECTRON_VOLT) + / (nd_plasma_species * constants.ELECTRON_CHARGE**2) + ) ** 0.5 + + @staticmethod + def calculate_lorentz_factor(velocity: float | np.ndarray) -> float | np.ndarray: + """ + Calculate the Lorentz factor for a given velocity. + :param velocity: Velocity in m/s. + :type velocity: float | np.ndarray + :returns: Lorentz factor (dimensionless). + :rtype: float | np.ndarray + """ + return 1 / (1 - (velocity / constants.SPEED_LIGHT) ** 2) ** 0.5 + + @staticmethod + def calculate_relativistic_particle_speed( + e_kinetic: float | np.ndarray, mass: float + ) -> float | np.ndarray: + """ + Calculate the speed of a particle given its kinetic energy and mass using relativistic mechanics. + :param e_kinetic: Kinetic energy in Joules. + :type e_kinetic: float | np.ndarray + :param mass: Mass of the particle in kg. + :type mass: float + :returns: Speed of the particle in m/s. + :rtype: float | np.ndarray + """ + return ( + constants.SPEED_LIGHT + * (1 - (1 / ((e_kinetic / (mass * constants.SPEED_LIGHT**2)) + 1) ** 2)) + ** 0.5 + ) + + def calculate_coulomb_log_from_impact( + self, impact_param_max: float, impact_param_min: float + ) -> float: + """ + Calculate the Coulomb logarithm from maximum and minimum impact parameters. + :param impact_param_max: Maximum impact parameter in meters. + :type impact_param_max: float + :param impact_param_min: Minimum impact parameter in meters. + :type impact_param_min: float + :returns: Coulomb logarithm (dimensionless). + :rtype: float + """ + return np.log(impact_param_max / impact_param_min) + + @staticmethod + def calculate_classical_distance_of_closest_approach( + charge1: float, + charge2: float, + e_kinetic: float | np.ndarray, + ) -> float | np.ndarray: + """ + Calculate the classical distance of closest approach for two charged particles. + + :param charge1: Charge of particle 1 in units of elementary charge. + :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 + :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 + ) + + @staticmethod + def calculate_debroglie_wavelength( + mass: float, velocity: float | np.ndarray + ) -> float | np.ndarray: + """ + Calculate the de Broglie wavelength of a particle. + :param mass: Mass of the particle in kg. + :type mass: float + :param velocity: Velocity of the particle in m/s. + :type velocity: float | np.ndarray + :returns: de Broglie wavelength in meters. + :rtype: float | np.ndarray + + :note: Reduced Planck constant (h-bar) is used in the calculation as this is for scattering. + """ + return (constants.PLANCK_CONSTANT / (2 * np.pi)) / (mass * velocity) + + @staticmethod + def calculate_plasma_frequency( + nd_particle: float | np.ndarray, m_particle: float, z_particle: float + ) -> float | np.ndarray: + """ + Calculate the plasma frequency for a particle species. + :param nd_particle: Number density of the particle species (/m^3). + :type nd_particle: float | np.ndarray + :param m_particle: Mass of the particle species (kg). + :type m_particle: float + :param Z_particle: Charge state of the particle species (dimensionless). + :type Z_particle: float + :returns: Plasma frequency in Hz. + :rtype: float | np.ndarray + """ + return ( + ( + (nd_particle * z_particle**2 * constants.ELECTRON_CHARGE**2) + / (m_particle * constants.EPSILON0) + ) + ** 0.5 + ) / (2 * np.pi) + + @staticmethod + def calculate_larmor_frequency( + b_field: float | np.ndarray, m_particle: float, z_particle: float + ) -> float | np.ndarray: + """ + Calculate the Larmor frequency for a particle species. + :param b_field: Magnetic field strength (T). + :type b_field: float | np.ndarray + :param m_particle: Mass of the particle species (kg). + :type m_particle: float + :param Z_particle: Charge state of the particle species (dimensionless). + :type Z_particle: float + :returns: Larmor frequency in Hz. + :rtype: float | np.ndarray + """ + return (z_particle * constants.ELECTRON_CHARGE * b_field) / ( + 2 * np.pi * m_particle + ) + + def output_detailed_physics(self): + """Outputs detailed physics variables to file.""" + + po.oheadr(self.outfile, "Detailed Plasma") + + po.osubhd(self.outfile, "Debye lengths:") + + po.ovarrf( + self.outfile, + "Plasma volume averaged electron Debye length (m)", + "(len_plasma_debye_electron_vol_avg)", + physics_variables.len_plasma_debye_electron_vol_avg, + "OP ", + ) + for i in range(len(physics_variables.len_plasma_debye_electron_profile)): + po.ovarre( + self.mfile, + f"Plasma electron Debye length at point {i}", + f"(len_plasma_debye_electron_profile{i})", + physics_variables.len_plasma_debye_electron_profile[i], + ) + + po.osubhd(self.outfile, "Velocities:") + + for i in range(len(physics_variables.vel_plasma_electron_profile)): + po.ovarre( + self.mfile, + f"Plasma electron thermal velocity at point {i}", + f"(vel_plasma_electron_profile{i})", + physics_variables.vel_plasma_electron_profile[i], + ) + + po.osubhd(self.outfile, "Frequencies:") + + for i in range(len(physics_variables.freq_plasma_electron_profile)): + po.ovarre( + self.mfile, + f"Plasma electron frequency at point {i}", + f"(freq_plasma_electron_profile{i})", + physics_variables.freq_plasma_electron_profile[i], + ) + for i in range( + len(physics_variables.freq_plasma_larmor_toroidal_electron_profile) + ): + po.ovarre( + self.mfile, + f"Plasma electron Larmor frequency at point {i}", + f"(freq_plasma_larmor_toroidal_electron_profile{i})", + physics_variables.freq_plasma_larmor_toroidal_electron_profile[i], + ) + + po.osubhd(self.outfile, "Coulomb Logarithms:") + + for i in range( + len(physics_variables.plasma_coulomb_log_electron_electron_profile) + ): + po.ovarre( + self.mfile, + f"Electron-electron Coulomb log at point {i}", + f"(plasma_coulomb_log_electron_electron_profile{i})", + physics_variables.plasma_coulomb_log_electron_electron_profile[i], + ) diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index a262ba33c1..3caf874982 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -21,6 +21,7 @@ ) from process.impurity_radiation import initialise_imprad from process.physics import ( + DetailedPhysics, Physics, calculate_beta_limit, calculate_current_coefficient_hastie, @@ -3466,3 +3467,168 @@ def test_calculate_internal_inductance_iter_3(): b_plasma_poloidal_vol_avg=1.0, c_plasma=1.5e7, vol_plasma=1000.0, rmajor=6.2 ) assert result == pytest.approx(0.9078959099585583, abs=0.00001) + + +@pytest.mark.parametrize( + "b_field, m_particle, z_particle, expected", + [ + ( + 1.0, + constants.ELECTRON_MASS, + 1.0, + 2.799249e10, + ), # typical electron in 1T field + (2.0, constants.ELECTRON_MASS, 1.0, 5.598498e10), # double field + (1.0, constants.PROTON_MASS, 1.0, 15245186.43761083), # proton in 1T field + (0.5, constants.ELECTRON_MASS, 2.0, 2.799249e10), # half field, double charge + (0.0, constants.ELECTRON_MASS, 1.0, 0.0), # zero field + ], +) +def test_calculate_larmor_frequency(b_field, m_particle, z_particle, expected): + """Test calculate_larmor_frequency for various particles and fields.""" + result = DetailedPhysics.calculate_larmor_frequency( + b_field=b_field, m_particle=m_particle, z_particle=z_particle + ) + assert result == pytest.approx(expected, rel=1e-5) + + +@pytest.mark.parametrize( + "nd_particle,m_particle,z_particle, expected", + ( + (1.0e20, constants.ELECTRON_MASS, 1.0, 89786628157.96086), + (1.0e19, constants.PROTON_MASS, 1.0, 662608904.2919972), + (5.0e19, constants.PROTON_MASS, 2.0, 2963277104.987116), + (0.0, constants.ELECTRON_MASS, 1.0, 0.0), + ), +) +def test_calculate_plasma_frequency(nd_particle, m_particle, z_particle, expected): + """Parametrised tests for DetailedPhysics.calculate_plasma_frequency().""" + + result = DetailedPhysics.calculate_plasma_frequency( + nd_particle, m_particle, z_particle + ) + assert result == pytest.approx(expected, rel=1e-12, abs=1e-12) + + +@pytest.mark.parametrize( + "mass,velocity,expected", + ( + ( + constants.ELECTRON_MASS, + constants.SPEED_LIGHT * 0.1, + 3.861592674352376e-12, + ), + ( + constants.PROTON_MASS, + constants.SPEED_LIGHT * 0.5, + 4.2061782010279145e-16, + ), + ( + constants.PROTON_MASS, + 1e3, + 6.304902508360882e-11, + ), + ), +) +def test_calculate_debroglie_wavelength(mass, velocity, expected): + """Test DetailedPhysics.calculate_debroglie_wavelength with several parameters.""" + result = DetailedPhysics.calculate_debroglie_wavelength(mass, velocity) + assert result == pytest.approx(expected, rel=1e-12) + + +@pytest.mark.parametrize( + "e_kev, mass, expected", + ( + (0.0, constants.ELECTRON_MASS, 0.0), + (1.6e-19, constants.ELECTRON_MASS, 592693.0770572403), + (1.0e-13, constants.ELECTRON_MASS, 267699064.11978555), + (1.0e-10, constants.PROTON_MASS, 239716127.82335472), + ), +) +def test_calculate_relativistic_particle_speed(e_kev, mass, expected): + """Parametrised tests for DetailedPhysics.calculate_relativistic_particle_speed""" + + result = DetailedPhysics.calculate_relativistic_particle_speed( + e_kinetic=e_kev, mass=mass + ) + assert result == pytest.approx(expected, rel=1e-5) + + +@pytest.mark.parametrize( + "velocity, expected_gamma", + ( + (0.0, 1.0), + (0.6 * constants.SPEED_LIGHT, 1.25), + (0.99 * constants.SPEED_LIGHT, 7.088812050083354), + ), +) +def test_calculate_lorentz_factor(velocity, expected_gamma): + """Test DetailedPhysics.calculate_lorentz_factor for several velocities.""" + + result = DetailedPhysics.calculate_lorentz_factor(velocity=velocity) + assert result == pytest.approx(expected_gamma, rel=1e-12) + + +@pytest.mark.parametrize( + "temp_keV, nd, expected", + ( + (1.0, 1e19, 7.433941993525029e-05), + (10.0, 1e20, 7.433941993525029e-05), + (0.1, 1e18, 7.433941993525029e-05), + ), +) +def test_calculate_debye_length_parametrized(temp_keV, nd, expected): + """Parametrized test for DetailedPhysics.calculate_debye_length.""" + result = DetailedPhysics.calculate_debye_length(temp_keV, nd) + assert result == pytest.approx(expected, rel=1e-12) + + +def test_detailed_physics_run_computes_profiles(): + # Minimal plasma profile + plasma = PlasmaProfile() + plasma.teprofile.profile_x = np.array([0.0, 0.5, 1.0]) + plasma.teprofile.profile_y = np.array([1.0, 2.0, 3.0]) # keV + plasma.neprofile.profile_x = plasma.teprofile.profile_x + plasma.neprofile.profile_y = np.array([1.0e19, 2.0e19, 3.0e19]) # m^-3 + + # Set global physics variables required by DetailedPhysics.run + physics_variables.temp_plasma_electron_vol_avg_kev = float( + np.mean(plasma.teprofile.profile_y) + ) + physics_variables.nd_plasma_electrons_vol_avg = float( + np.mean(plasma.neprofile.profile_y) + ) + # toroidal field profile for larmor frequency calc + physics_variables.b_plasma_toroidal_profile = ( + np.ones_like(plasma.neprofile.profile_y) * 5.0 + ) + + dp = DetailedPhysics(plasma) + + # Run should complete without error and populate physics_variables + dp.run() + + n = len(plasma.teprofile.profile_y) + assert physics_variables.len_plasma_debye_electron_vol_avg > 0 + assert hasattr(physics_variables, "len_plasma_debye_electron_profile") + assert np.shape(physics_variables.len_plasma_debye_electron_profile)[0] == n + + assert np.shape(physics_variables.vel_plasma_electron_profile)[0] == n + assert np.all(np.isfinite(physics_variables.vel_plasma_electron_profile)) + + assert np.shape(physics_variables.freq_plasma_electron_profile)[0] == n + assert np.all(np.isfinite(physics_variables.freq_plasma_electron_profile)) + + assert ( + np.shape(physics_variables.freq_plasma_larmor_toroidal_electron_profile)[0] == n + ) + assert np.all( + np.isfinite(physics_variables.freq_plasma_larmor_toroidal_electron_profile) + ) + + assert ( + np.shape(physics_variables.plasma_coulomb_log_electron_electron_profile)[0] == n + ) + assert np.all( + np.isfinite(physics_variables.plasma_coulomb_log_electron_electron_profile) + )