From e7d1a13f33bfb2a991ddf9d3618c2db78dae8885 Mon Sep 17 00:00:00 2001
From: kevinm387 <65828579+kevinm387@users.noreply.github.com>
Date: Tue, 25 Jul 2023 14:11:21 -0700
Subject: [PATCH] Add get_activity() function to deplete.Results class (#2617)

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

diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py
index f250ab3639..bd1f884d58 100644
--- a/openmc/deplete/results.py
+++ b/openmc/deplete/results.py
@@ -2,7 +2,7 @@ import numbers
 import bisect
 import math
 import typing  # required to prevent typing.Union namespace overwriting Union
-from typing import Iterable, Optional, Tuple
+from typing import Iterable, Optional, Tuple, List
 from warnings import warn
 
 import h5py
@@ -95,6 +95,59 @@ class Results(list):
         )
         return cls(filename)
 
+    def get_activity(
+        self,
+        mat: typing.Union[Material, str],
+        units: str = "Bq/cm3",
+        by_nuclide: bool = False, 
+        volume: Optional[float] = None
+    ) -> Tuple[np.ndarray, typing.Union[np.ndarray, List[dict]]]:
+        """Get activity of material over time.
+
+        Parameters
+        ----------
+        mat : openmc.Material, str
+            Material object or material id to evaluate
+        units : {'Bq', 'Bq/g', 'Bq/cm3'}
+            Specifies the type of activity to return, options include total
+            activity [Bq], specific [Bq/g] or volumetric activity [Bq/cm3].
+        by_nuclide : bool
+            Specifies if the activity should be returned for the material as a
+            whole or per nuclide. Default is False.
+        volume : float, optional
+            Volume of the material. If not passed, defaults to using the
+            :attr:`Material.volume` attribute.
+
+        Returns
+        -------
+        times : numpy.ndarray
+            Array of times in [s]
+        activities : numpy.ndarray
+            Array of total activities if by_nuclide = False (default)
+            or array of dictionaries of activities by nuclide if
+            by_nuclide = True.
+            
+        """
+        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 = np.empty_like(self, dtype=float)
+        if by_nuclide:
+            activities = [None] * len(self)
+        else:
+            activities = np.empty_like(self, dtype=float)
+
+        # Evaluate activity for each depletion time
+        for i, result in enumerate(self):
+            times[i] = result.time[0]
+            activities[i] = result.get_material(mat_id).get_activity(units, by_nuclide, volume)
+        
+        return times, activities
+
     def get_atoms(
         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 457a6fea29..f5716eb8a7 100644
--- a/tests/unit_tests/test_deplete_resultslist.py
+++ b/tests/unit_tests/test_deplete_resultslist.py
@@ -15,6 +15,27 @@ def res():
                 / 'test_reference.h5')
     return openmc.deplete.Results(filename)
 
+def test_get_activity(res):
+    """Tests evaluating activity"""
+    t, a = res.get_activity("1")
+
+    t_ref = np.array([0.0, 1296000.0, 2592000.0, 3888000.0])
+    a_ref = np.array(
+        [1.25167956e+06, 3.71938527e+11, 4.43264300e+11, 3.55547176e+11])
+
+    np.testing.assert_allclose(t, t_ref)
+    np.testing.assert_allclose(a, a_ref)
+
+    # Check by_nuclide
+    a_xe135_ref = np.array(
+        [2.106574218e+05, 1.227519888e+11, 1.177491828e+11, 1.031986176e+11])
+    t_nuc, a_nuc = res.get_activity("1", by_nuclide=True)
+    
+    a_xe135 = np.array([a_nuc_i["Xe135"] for a_nuc_i in a_nuc])
+
+    np.testing.assert_allclose(t_nuc, t_ref)
+    np.testing.assert_allclose(a_xe135, a_xe135_ref)
+
 
 def test_get_atoms(res):
     """Tests evaluating single nuclide concentration."""