🏋️ Calculate plasma fuelling masses

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,20 +4237,55 @@ 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)",
             "(m_fuel_amu)",
             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)",
             "(m_beam_amu)",
             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

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

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)