Add get_decay_heat() function to deplete.Results (#2625)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

openmc/deplete/results.py

@@ -99,7 +99,7 @@ class Results(list):
         self,
         mat: typing.Union[Material, str],
         units: str = "Bq/cm3",
-        by_nuclide: bool = False, 
+        by_nuclide: bool = False,
         volume: Optional[float] = None
     ) -> Tuple[np.ndarray, typing.Union[np.ndarray, List[dict]]]:
         """Get activity of material over time.
@@ -122,11 +122,11 @@ class Results(list):
         -------
         times : numpy.ndarray
             Array of times in [s]
-        activities : numpy.ndarray
+        activities : numpy.ndarray or List[dict]
             Array of total activities if by_nuclide = False (default)
-            or array of dictionaries of activities by nuclide if
+            or list of dictionaries of activities by nuclide if
             by_nuclide = True.
-            
+
         """
         if isinstance(mat, Material):
             mat_id = str(mat.id)
@@ -134,7 +134,7 @@ class Results(list):
             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)
@@ -145,7 +145,7 @@ class Results(list):
         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(
@@ -211,6 +211,60 @@ class Results(list):
 
         return times, concentrations
 
+    def get_decay_heat(
+            self,
+            mat: typing.Union[Material, str],
+            units: str = "W",
+            by_nuclide: bool = False,
+            volume: Optional[float] = None
+    ) -> Tuple[np.ndarray, typing.Union[np.ndarray, List[dict]]]:
+        """Get decay heat of material over time.
+
+        Parameters
+        ----------
+        mat : openmc.Material, str
+            Material object or material id to evaluate.
+        units : {'W', 'W/g', 'W/cm3'}
+            Specifies the units of decay heat to return. Options include total
+            heat [W], specific [W/g] or volumetric heat [W/cm3].
+        by_nuclide : bool
+            Specifies if the decay heat 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]
+        decay_heat : numpy.ndarray or List[dict]
+            Array of total decay heat values if by_nuclide = False (default)
+            or list of dictionaries of decay heat values 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:
+            decay_heat = [None] * len(self)
+        else:
+            decay_heat = np.empty_like(self, dtype=float)
+
+        # Evaluate decay heat for each depletion time
+        for i, result in enumerate(self):
+            times[i] = result.time[0]
+            decay_heat[i] = result.get_material(mat_id).get_decay_heat(
+                units, by_nuclide, volume)
+
+        return times, decay_heat
+
     def get_mass(self,
         mat: typing.Union[Material, str],
         nuc: str,

tests/unit_tests/test_deplete_resultslist.py

@@ -30,7 +30,7 @@ def test_get_activity(res):
     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)
@@ -64,6 +64,31 @@ def test_get_atoms(res):
     assert t_hour == pytest.approx(t_ref / (60 * 60))
 
 
+def test_get_decay_heat(res):
+    """Tests evaluating decay heat."""
+    # Set chain file for testing
+    openmc.config['chain_file'] = Path(__file__).parents[1] / 'chain_simple.xml'
+
+    t_ref = np.array([0.0, 1296000.0, 2592000.0, 3888000.0])
+    dh_ref = np.array(
+        [1.27933813e-09, 5.85347232e-03, 7.38773010e-03, 5.79954067e-03])
+
+    t, dh = res.get_decay_heat("1")
+
+    np.testing.assert_allclose(t, t_ref)
+    np.testing.assert_allclose(dh, dh_ref)
+
+    # Check by nuclide
+    dh_xe135_ref = np.array(
+        [1.27933813e-09, 7.45481920e-04, 7.15099509e-04, 6.26732849e-04])
+    t_nuc, dh_nuc = res.get_decay_heat("1", by_nuclide=True)
+
+    dh_nuc_xe135 = np.array([dh_nuc_i["Xe135"] for dh_nuc_i in dh_nuc])
+
+    np.testing.assert_allclose(t_nuc, t_ref)
+    np.testing.assert_allclose(dh_nuc_xe135, dh_xe135_ref)
+
+
 def test_get_mass(res):
     """Tests evaluating single nuclide concentration."""
     t, n = res.get_mass("1", "Xe135")