From 51cf475723339f619bd119c6da635385ef2d26ae Mon Sep 17 00:00:00 2001
From: Olek <45364492+yardasol@users.noreply.github.com>
Date: Mon, 26 Jun 2023 15:11:41 +0200
Subject: [PATCH] Add `get_mass()` function to `openmc.deplete.Results` (#2565)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
---
 openmc/deplete/results.py                    | 54 ++++++++++++++++++++
 tests/unit_tests/test_deplete_resultslist.py | 29 +++++++++++
 2 files changed, 83 insertions(+)

diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py
index 9cc3b439a60..f250ab36397 100644
--- a/openmc/deplete/results.py
+++ b/openmc/deplete/results.py
@@ -10,6 +10,7 @@
 
 from .stepresult import StepResult, VERSION_RESULTS
 import openmc.checkvalue as cv
+from openmc.data import atomic_mass, AVOGADRO
 from openmc.data.library import DataLibrary
 from openmc.material import Material, Materials
 from openmc.exceptions import DataError
@@ -157,6 +158,59 @@ def get_atoms(
 
         return times, concentrations
 
+    def get_mass(self,
+        mat: typing.Union[Material, str],
+        nuc: str,
+        mass_units: str = "g",
+        time_units: str = "s"
+    ) -> Tuple[np.ndarray, np.ndarray]:
+        """Get mass of nuclides over time from a single material
+
+        .. versionadded:: 0.13.4
+
+        Parameters
+        ----------
+        mat : openmc.Material, str
+            Material object or material id to evaluate
+        nuc : str
+            Nuclide name to evaluate
+        mass_units : {"g", "g/cm3", "kg"}, optional
+            Units for the returned mass.
+        time_units : {"s", "min", "h", "d", "a"}, optional
+            Units for the returned time array. Default is ``"s"`` to
+            return the value in seconds. Other options are minutes ``"min"``,
+            hours ``"h"``, days ``"d"``, and Julian years ``"a"``.
+
+        Returns
+        -------
+        times : numpy.ndarray
+            Array of times in units of ``time_units``
+        mass : numpy.ndarray
+            Mass of specified nuclide in units of ``mass_units``
+
+        """
+        cv.check_value("mass_units", mass_units, {"g", "g/cm3", "kg"})
+
+        if isinstance(mat, Material):
+            mat_id = str(mat.id)
+        elif isinstance(mat, str):
+            mat_id = mat
+        else:
+            raise TypeError('mat should be of type openmc.Material or str')
+
+        times, atoms = self.get_atoms(mat, nuc, time_units=time_units)
+
+        mass = atoms * atomic_mass(nuc) / AVOGADRO
+
+        # Unit conversions
+        if mass_units == "g/cm3":
+            # Divide by volume to get density
+            mass /= self[0].volume[mat_id]
+        elif mass_units == "kg":
+            mass *= 1e3
+
+        return times, mass
+
     def get_reaction_rate(
         self,
         mat: typing.Union[Material, str],
diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py
index 6e4a4210bdd..457a6fea294 100644
--- a/tests/unit_tests/test_deplete_resultslist.py
+++ b/tests/unit_tests/test_deplete_resultslist.py
@@ -43,6 +43,35 @@ def test_get_atoms(res):
     assert t_hour == pytest.approx(t_ref / (60 * 60))
 
 
+def test_get_mass(res):
+    """Tests evaluating single nuclide concentration."""
+    t, n = res.get_mass("1", "Xe135")
+
+    t_ref = np.array([0.0, 1296000.0, 2592000.0, 3888000.0])
+    n_ref = np.array(
+        [6.67473282e+08, 3.88942731e+14, 3.73091215e+14, 3.26987387e+14])
+
+    # Get g
+    n_ref *= openmc.data.atomic_mass('Xe135') / openmc.data.AVOGADRO
+
+    np.testing.assert_allclose(t, t_ref)
+    np.testing.assert_allclose(n, n_ref)
+
+    # Check alternate units
+    volume = res[0].volume["1"]
+    t_days, n_cm3 = res.get_mass("1", "Xe135", mass_units="g/cm3", time_units="d")
+
+    assert t_days == pytest.approx(t_ref / (60 * 60 * 24))
+    assert n_cm3 == pytest.approx(n_ref / volume)
+
+    t_min, n_bcm = res.get_mass("1", "Xe135", mass_units="kg", time_units="min")
+    assert n_bcm == pytest.approx(n_ref * 1e3)
+    assert t_min == pytest.approx(t_ref / 60)
+
+    t_hour, _n = res.get_mass("1", "Xe135", time_units="h")
+    assert t_hour == pytest.approx(t_ref / (60 * 60))
+
+
 def test_get_reaction_rate(res):
     """Tests evaluating reaction rate."""
     t, r = res.get_reaction_rate("1", "Xe135", "(n,gamma)")