Allow pure decay IndependentOperator (#2966)

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

openmc/deplete/abc.py

@@ -848,7 +848,7 @@ class Integrator(ABC):
                 print(f"[openmc.deplete] t={t} (final operator evaluation)")
             res_list = [self.operator(n, source_rate if final_step else 0.0)]
             StepResult.save(self.operator, [n], res_list, [t, t],
-                         source_rate, self._i_res + len(self), proc_time)
+                         source_rate, self._i_res + len(self), proc_time, path)
             self.operator.write_bos_data(len(self) + self._i_res)
 
         self.operator.finalize()

openmc/deplete/independent_operator.py

@@ -241,7 +241,6 @@ class IndependentOperator(OpenMCOperator):
                    reduce_chain_level=reduce_chain_level,
                    fission_yield_opts=fission_yield_opts)
 
-
     @staticmethod
     def _consolidate_nuclides_to_material(nuclides, nuc_units, volume):
         """Puts nuclide list into an openmc.Materials object.
@@ -270,7 +269,6 @@ class IndependentOperator(OpenMCOperator):
         self.prev_res[-1].transfer_volumes(model)
         self.materials = model.materials
 
-
         # Store previous results in operator
         # Distribute reaction rates according to those tracked
         # on this process
@@ -282,7 +280,6 @@ class IndependentOperator(OpenMCOperator):
                 new_res = res_obj.distribute(self.local_mats, mat_indexes)
                 self.prev_res.append(new_res)
 
-
     def _get_nuclides_with_data(self, cross_sections: List[MicroXS]) -> Set[str]:
         """Finds nuclides with cross section data"""
         return set(cross_sections[0].nuclides)
@@ -412,6 +409,12 @@ class IndependentOperator(OpenMCOperator):
 
         self._update_materials_and_nuclides(vec)
 
+        # If the source rate is zero, return zero reaction rates
+        if source_rate == 0.0:
+            rates = self.reaction_rates.copy()
+            rates.fill(0.0)
+            return OperatorResult(ufloat(0.0, 0.0), rates)
+
         rates = self._calculate_reaction_rates(source_rate)
         keff = self._keff
 

tests/unit_tests/test_deplete_decay.py

@@ -1,3 +1,5 @@
+from pathlib import Path
+
 import openmc.deplete
 import numpy as np
 import pytest
@@ -45,3 +47,37 @@ def test_deplete_decay_products(run_in_tmpdir):
     # H1 and He4
     assert h1[1] == pytest.approx(1e24)
     assert he4[1] == pytest.approx(1e24)
+
+
+def test_deplete_decay_step_fissionable(run_in_tmpdir):
+    """Ensures that power is not computed in zero power cases with
+    fissionable material present. This tests decay calculations without
+    power, although this specific example does not exhibit any decay.
+
+    Proves github issue #2963 is fixed
+    """
+
+    # Set up a pure decay operator
+    micro_xs = openmc.deplete.MicroXS(np.empty((0, 0)), [], [])
+    mat = openmc.Material()
+    mat.name = 'I do not decay.'
+    mat.add_nuclide('U238', 1.0, 'ao')
+    mat.volume = 10.0
+    mat.set_density('g/cc', 1.0)
+    original_atoms = mat.get_nuclide_atoms()['U238']
+
+    mats = openmc.Materials([mat])
+    op = openmc.deplete.IndependentOperator(
+        mats, [1.0], [micro_xs], Path(__file__).parents[1] / "chain_simple.xml")
+
+    # Create time integrator and integrate
+    integrator = openmc.deplete.PredictorIntegrator(
+        op, [1.0], power=[0.0], timestep_units='s'
+    )
+    integrator.integrate()
+
+    # Get concentration of U238. It should be unchanged since this chain has no U238 decay.
+    results = openmc.deplete.Results('depletion_results.h5')
+    _, u238 = results.get_atoms("1", "U238")
+
+    assert u238[1] == pytest.approx(original_atoms)