diff --git a/process/physics.py b/process/physics.py index 41bfdb23f5..1909f655a1 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1548,6 +1548,22 @@ def physics(self): # Issue #261 Remove old radiation model (imprad_model=0) self.plasma_composition() + ( + physics_variables.m_plasma_fuel_ions, + physics_variables.m_plasma_ions_total, + physics_variables.m_plasma_alpha, + physics_variables.m_plasma_electron, + physics_variables.m_plasma, + ) = self.calculate_plasma_masses( + physics_variables.m_fuel_amu, + physics_variables.m_ions_total_amu, + physics_variables.nd_ions_total, + physics_variables.nd_fuel_ions, + physics_variables.nd_alphas, + physics_variables.vol_plasma, + physics_variables.dene, + ) + # Define coulomb logarithm # (collisions: ion-electron, electron-electron) physics_variables.dlamee = ( @@ -4221,6 +4237,13 @@ def outplas(self): physics_variables.m_ions_total_amu, "OP ", ) + po.ovarre( + self.outfile, + "Total mass of all ions in plasma (kg)", + "(m_plasma_ions_total)", + physics_variables.m_plasma_ions_total, + "OP ", + ) po.ovarre( self.outfile, "Average mass of all fuel ions (amu)", @@ -4228,6 +4251,13 @@ def outplas(self): physics_variables.m_fuel_amu, "OP ", ) + po.ovarre( + self.outfile, + "Total mass of all fuel ions in plasma (kg)", + "(m_plasma_fuel_ions)", + physics_variables.m_plasma_fuel_ions, + "OP ", + ) po.ovarre( self.outfile, "Average mass of all beam ions (amu)", @@ -4235,6 +4265,27 @@ def outplas(self): physics_variables.m_beam_amu, "OP ", ) + po.ovarre( + self.outfile, + "Total mass of all alpha particles in plasma (kg)", + "(m_plasma_alpha)", + physics_variables.m_plasma_alpha, + "OP ", + ) + po.ovarre( + self.outfile, + "Total mass of all electrons in plasma (kg)", + "(m_plasma_electron)", + physics_variables.m_plasma_electron, + "OP ", + ) + po.ovarre( + self.outfile, + "Total mass of the plasma (kg)", + "(m_plasma)", + physics_variables.m_plasma, + "OP ", + ) po.oblnkl(self.outfile) po.ovarrf( @@ -7620,6 +7671,61 @@ def calculate_confinement_time( p_plasma_loss_mw, ) + @staticmethod + def calculate_plasma_masses( + m_fuel_amu: float, + m_ions_total_amu: float, + nd_ions_total: float, + nd_fuel_ions: float, + nd_alphas: float, + vol_plasma: float, + dene: float, + ) -> tuple[float, float, float, float, float]: + """ + Calculate the plasma masses. + + :param m_fuel_amu: Average mass of fuel (amu). + :type m_fuel_amu: float + :param m_ions_total_amu: Average mass of all ions (amu). + :type m_ions_total_amu: float + :param nd_ions_total: Total ion density (/m3). + :type nd_ions_total: float + :param nd_fuel_ions: Fuel ion density (/m3). + :type nd_fuel_ions: float + :param nd_alphas: Alpha ash density (/m3). + :type nd_alphas: float + :param vol_plasma: Plasma volume (m3). + :type vol_plasma: float + :param dene: Volume averaged electron density (/m3). + :type dene: float + + :returns: A tuple containing: + :rtype: tuple[float, float, float, float, float] + """ + + # Calculate mass of fuel ions + m_plasma_fuel_ions = (m_fuel_amu * constants.atomic_mass_unit) * ( + nd_fuel_ions * vol_plasma + ) + + m_plasma_ions_total = (m_ions_total_amu * constants.atomic_mass_unit) * ( + nd_ions_total * vol_plasma + ) + + m_plasma_alpha = (nd_alphas * vol_plasma) * constants.alpha_mass + + m_plasma_electron = constants.electron_mass * (dene * vol_plasma) + + m_plasma = m_plasma_electron + m_plasma_ions_total + + return ( + m_plasma_fuel_ions, + m_plasma_ions_total, + m_plasma_alpha, + m_plasma_electron, + m_plasma, + ) + def calculate_poloidal_beta(btot, bp, beta): """Calculates total poloidal beta diff --git a/source/fortran/physics_variables.f90 b/source/fortran/physics_variables.f90 index a4edcbbc96..1b9b3f43fe 100644 --- a/source/fortran/physics_variables.f90 +++ b/source/fortran/physics_variables.f90 @@ -26,6 +26,21 @@ module physics_variables real(dp) :: m_ions_total_amu !! average mass of all ions (amu) + real(dp) :: m_plasma_fuel_ions + !! Mass of the plasma fuel ions (kg) + + real(dp) :: m_plasma_ions_total + !! Mass of all ions in plasma (kg) + + real(dp) :: m_plasma_alpha + !! Mass of the alpha particles in the plasma (kg) + + real(dp) :: m_plasma_electron + !! Mass of the electrons in the plasma (kg) + + real(dp) :: m_plasma + !! Total mass of the plasma (kg) + real(dp) :: alphaj !! current profile index (calculated from q_0 and q if `iprofile=1`) @@ -903,6 +918,11 @@ subroutine init_physics_variables m_beam_amu = 0.0D0 m_fuel_amu = 0.0D0 m_ions_total_amu = 0.0D0 + m_plasma_fuel_ions = 0.0D0 + m_plasma_ions_total = 0.0D0 + m_plasma_alpha = 0.0D0 + m_plasma_electron = 0.0D0 + m_plasma = 0.0D0 alphaj = 1.0D0 alphan = 0.25D0 alphap = 0.0D0 diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 63c3a46403..22f6f2afcd 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -3055,3 +3055,36 @@ def test_calculate_confinement_time(confinementtimeparam, monkeypatch, physics): assert t_ion_energy_confinement == pytest.approx( confinementtimeparam.expected_t_ion_energy_confinement ) + + +def test_calculate_plasma_masses(): + """Test calculate_plasma_masses()""" + m_fuel_amu = 2.5 + m_ions_total_amu = 3.0 + nd_ions_total = 1.0e20 + nd_fuel_ions = 0.8e20 + nd_alphas = 0.1e20 + vol_plasma = 100.0 + dene = 1.0e20 + + ( + m_plasma_fuel_ions, + m_plasma_ions_total, + m_plasma_alpha, + m_plasma_electron, + m_plasma, + ) = Physics.calculate_plasma_masses( + m_fuel_amu=m_fuel_amu, + m_ions_total_amu=m_ions_total_amu, + nd_ions_total=nd_ions_total, + nd_fuel_ions=nd_fuel_ions, + nd_alphas=nd_alphas, + vol_plasma=vol_plasma, + dene=dene, + ) + + assert m_plasma_fuel_ions == pytest.approx(3.32107813784e-05, abs=1e-30) + assert m_plasma_ions_total == pytest.approx(4.9816172067599995e-05, abs=1e-30) + assert m_plasma_alpha == pytest.approx(6.644657345e-06, abs=1e-30) + assert m_plasma_electron == pytest.approx(9.1093837139e-09, abs=1e-34) + assert m_plasma == pytest.approx(4.982528145131389e-05, abs=1e-30)