From e7d8574d395f8dc0232d8131915d47fa70e66bbd Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sun, 12 Jul 2026 00:25:08 -0500
Subject: [PATCH 1/7] Split R2S photon source and transport steps

---
 docs/source/usersguide/decay_sources.rst |   6 +-
 openmc/deplete/r2s.py                    | 173 +++++++++++++----------
 tests/unit_tests/test_r2s.py             |  10 ++
 3 files changed, 114 insertions(+), 75 deletions(-)

diff --git a/docs/source/usersguide/decay_sources.rst b/docs/source/usersguide/decay_sources.rst
index 21981fdaa82..914fd278f67 100644
--- a/docs/source/usersguide/decay_sources.rst
+++ b/docs/source/usersguide/decay_sources.rst
@@ -148,12 +148,13 @@ a directory that is named "r2s_<timestamp>/". The ``output_dir`` argument to the
 :meth:`~openmc.deplete.R2SManager.run` method enables you to override the
 default output directory name if desired.
 
-The :meth:`~openmc.deplete.R2SManager.run` method actually runs three
+The :meth:`~openmc.deplete.R2SManager.run` method actually runs four
 lower-level methods under the hood::
 
     r2s.step1_neutron_transport(...)
     r2s.step2_activation(...)
-    r2s.step3_photon_transport(...)
+    r2s.step3_photon_source(...)
+    r2s.step4_photon_transport(...)
 
 For users looking for more control over the calculation, these lower-level
 methods can be used in lieu of the :meth:`openmc.deplete.R2SManager.run` method.
@@ -255,4 +256,3 @@ relevant tallies. This can be done with the aid of the
 
     # Apply time correction factors
     tally = d1s.apply_time_correction(dose_tally, factors, time_index)
-
diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py
index d41e230eedc..6cee45b1a9f 100644
--- a/openmc/deplete/r2s.py
+++ b/openmc/deplete/r2s.py
@@ -256,9 +256,12 @@ def run(
                 timesteps, source_rates, timestep_units,
                 output_dir / 'activation', operator_kwargs=operator_kwargs
             )
-            self.step3_photon_transport(
+            self.step3_photon_source(
                 photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
-                mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs
+                mat_vol_kwargs=mat_vol_kwargs
+            )
+            self.step4_photon_transport(
+                output_dir / 'photon_transport', run_kwargs=run_kwargs
             )
 
         return output_dir
@@ -438,23 +441,20 @@ def step2_activation(
         # Get depletion results
         self.results['depletion_results'] = Results(output_path)
 
-    def step3_photon_transport(
+    def step3_photon_source(
         self,
         time_indices: Sequence[int] | None = None,
         bounding_boxes: dict[int, openmc.BoundingBox] | None = None,
         output_dir: PathLike = 'photon_transport',
         mat_vol_kwargs: dict | None = None,
-        run_kwargs: dict | None = None,
     ):
-        """Run the photon transport step.
+        """Create decay photon sources.
 
-        This step performs photon transport calculations using decay photon
-        sources created from the activated materials. For each specified time,
-        it creates appropriate photon sources and runs a transport calculation.
-        In mesh-based mode, the sources are created using the mesh material
-        volumes, while in cell-based mode, they are created using bounding boxes
-        for each cell.  This step will populate the 'photon_tallies' key in the
-        results dictionary.
+        This step creates decay photon sources from the activated materials for
+        each specified time. In mesh-based mode, the sources are created using
+        mesh material volumes, while in cell-based mode, they are created using
+        bounding boxes for each cell. This step will populate the
+        'photon_sources' key in the results dictionary.
 
         Parameters
         ----------
@@ -470,13 +470,10 @@ def step3_photon_transport(
             source sampling in cell-based R2S calculations. Required if method
             is 'cell-based'.
         output_dir : PathLike, optional
-            Path to directory where photon transport outputs will be saved.
+            Path to directory where photon source outputs will be saved.
         mat_vol_kwargs : dict, optional
             Additional keyword arguments passed to
             :meth:`openmc.MeshBase.material_volumes`.
-        run_kwargs : dict, optional
-            Additional keyword arguments passed to :meth:`openmc.Model.run`
-            during the photon transport step. By default, output is disabled.
         """
 
         # TODO: Automatically determine bounding box for each cell
@@ -484,13 +481,6 @@ def step3_photon_transport(
             raise ValueError("bounding_boxes must be provided for cell-based "
                              "R2S calculations.")
 
-        # Set default run arguments if not provided
-        if run_kwargs is None:
-            run_kwargs = {}
-        run_kwargs.setdefault('output', False)
-
-        # Write out JSON file with tally IDs that can be used for loading
-        # results
         output_dir = Path(output_dir)
         output_dir.mkdir(parents=True, exist_ok=True)
 
@@ -498,6 +488,8 @@ def step3_photon_transport(
         if time_indices is None:
             n_steps = len(self.results['depletion_results'])
             time_indices = list(range(n_steps))
+        else:
+            time_indices = list(time_indices)
 
         # Check whether the photon model is different
         neutron_univ = self.neutron_model.geometry.root_universe
@@ -523,13 +515,6 @@ def step3_photon_transport(
 
             self.results['mesh_material_volumes_photon'] = photon_mmv_list
 
-        if comm.rank == 0:
-            tally_ids = [tally.id for tally in self.photon_model.tallies]
-            with open(output_dir / 'tally_ids.json', 'w') as f:
-                json.dump(tally_ids, f)
-
-        self.results['photon_tallies'] = {}
-
         # Get dictionary of cells in the photon model
         if different_photon_model:
             photon_cells = self.photon_model.geometry.get_all_cells()
@@ -547,16 +532,57 @@ def step3_photon_transport(
                         continue
                 work_items.append((cell, original_mat, bounding_boxes[cell.id]))
 
-        # Ensure photon transport is enabled in settings
-        self.photon_model.settings.photon_transport = True
+        # Normalize negative time indices before creating source lists.
+        n_steps = len(self.results['depletion_results'])
+        time_indices = [i + n_steps if i < 0 else i for i in time_indices]
+
+        self.results['photon_sources'] = {
+            time_index: self._create_photon_sources(time_index, work_items)
+            for time_index in dict.fromkeys(time_indices)
+        }
+
+    def step4_photon_transport(
+        self,
+        output_dir: PathLike = 'photon_transport',
+        run_kwargs: dict | None = None,
+    ):
+        """Run photon transport using prepared decay photon sources.
 
-        for time_index in time_indices:
-            # Convert time_index (which may be negative) to a normal index
-            if time_index < 0:
-                time_index += len(self.results['depletion_results'])
+        This step runs a photon transport calculation for each source list
+        created by :meth:`step3_photon_source`. It will populate the
+        'photon_tallies' key in the results dictionary.
 
-            # Build decay photon sources and assign to the photon model
-            sources = self._create_photon_sources(time_index, work_items)
+        Parameters
+        ----------
+        output_dir : PathLike, optional
+            Path to directory where photon transport outputs will be saved.
+        run_kwargs : dict, optional
+            Additional keyword arguments passed to :meth:`openmc.Model.run`.
+            By default, output is disabled.
+        """
+        if 'photon_sources' not in self.results:
+            raise RuntimeError(
+                'Photon sources must be created with step3_photon_source '
+                'before running photon transport.')
+
+        if run_kwargs is None:
+            run_kwargs = {}
+        run_kwargs.setdefault('output', False)
+
+        output_dir = Path(output_dir)
+        output_dir.mkdir(parents=True, exist_ok=True)
+
+        if comm.rank == 0:
+            tally_ids = [tally.id for tally in self.photon_model.tallies]
+            with open(output_dir / 'tally_ids.json', 'w') as f:
+                json.dump(tally_ids, f)
+
+        self.results['photon_tallies'] = {}
+
+        # Ensure photon transport is enabled in settings.
+        self.photon_model.settings.photon_transport = True
+
+        for time_index, sources in self.results['photon_sources'].items():
             self.photon_model.settings.source = sources
 
             # Run photon transport calculation
@@ -635,45 +661,48 @@ def _create_photon_sources(self, time_index, work_items):
         list of openmc.IndependentSource
             Photon sources for each activated region.
         """
-        step_result = self.results['depletion_results'][time_index]
-        materials = self.results['activation_materials']
-        mesh_based = self.method == 'mesh-based'
-        if mesh_based:
-            mat_dict = self.neutron_model._get_all_materials()
-
+        mat_dict = (
+            self.neutron_model._get_all_materials()
+            if self.method == 'mesh-based' else None
+        )
         sources = []
         for item in work_items:
-            if mesh_based:
-                index_mat, domain_id, bbox = item
-                original_mat = materials[index_mat]
-                domain = mat_dict[domain_id]
-            else:
-                cell, original_mat, bbox = item
-                domain = cell
-
-            activated_mat = step_result.get_material(str(original_mat.id))
-            nuclides = activated_mat.get_nuclide_atom_densities()
-            if not nuclides:
-                continue
-
-            # Eliminate nuclides with zero density
-            nuclides = {nuclide: density for nuclide, density in nuclides.items()
-                        if density > 0}
-
-            energy = openmc.stats.DecaySpectrum(nuclides, activated_mat.volume)
-            energy.clip(inplace=True)
-            if not energy.nuclides:
-                continue
-
-            sources.append(openmc.IndependentSource(
-                space=openmc.stats.Box(bbox.lower_left, bbox.upper_right),
-                energy=energy,
-                particle='photon',
-                constraints={'domains': [domain]},
-            ))
+            source = self._create_photon_source(time_index, item, mat_dict)
+            if source is not None:
+                sources.append(source)
 
         return sources
 
+    def _create_photon_source(self, time_index, item, mat_dict=None):
+        """Create a decay photon source for one activation region."""
+        step_result = self.results['depletion_results'][time_index]
+        if self.method == 'mesh-based':
+            index_mat, domain_id, bbox = item
+            original_mat = self.results['activation_materials'][index_mat]
+            domain = mat_dict[domain_id]
+        else:
+            domain, original_mat, bbox = item
+
+        activated_mat = step_result.get_material(str(original_mat.id))
+        nuclides = activated_mat.get_nuclide_atom_densities()
+        if not nuclides:
+            return None
+
+        # Eliminate nuclides with zero density.
+        nuclides = {nuclide: density for nuclide, density in nuclides.items()
+                    if density > 0}
+        energy = openmc.stats.DecaySpectrum(nuclides, activated_mat.volume)
+        energy.clip(inplace=True)
+        if not energy.nuclides:
+            return None
+
+        return openmc.IndependentSource(
+            space=openmc.stats.Box(bbox.lower_left, bbox.upper_right),
+            energy=energy,
+            particle='photon',
+            constraints={'domains': [domain]},
+        )
+
     def load_results(self, path: PathLike):
         """Load results from a previous R2S calculation.
 
diff --git a/tests/unit_tests/test_r2s.py b/tests/unit_tests/test_r2s.py
index 266bd49763a..caf3cb8cd76 100644
--- a/tests/unit_tests/test_r2s.py
+++ b/tests/unit_tests/test_r2s.py
@@ -82,6 +82,8 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
     assert len(r2s.results['mesh_material_volumes'][0]) == 2
     assert len(r2s.results['activation_materials']) == 2
     assert len(r2s.results['depletion_results']) == 2
+    assert list(r2s.results['photon_sources']) == [1]
+    assert r2s.results['photon_sources'][1]
 
     # Check activation materials
     amats = r2s.results['activation_materials']
@@ -217,3 +219,11 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path):
     assert len(r2s_loaded.results['micros']) == 2
     assert len(r2s_loaded.results['activation_materials']) == 2
     assert len(r2s_loaded.results['depletion_results']) == 2
+
+
+def test_step4_requires_photon_sources(simple_model_and_mesh):
+    model, (c1, c2), _ = simple_model_and_mesh
+    r2s = R2SManager(model, [c1, c2])
+
+    with pytest.raises(RuntimeError, match='step3_photon_source'):
+        r2s.step4_photon_transport()

From a5b6ea0e0d0385eb1bfb8753985509b5e6343be9 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sun, 12 Jul 2026 00:28:22 -0500
Subject: [PATCH 2/7] Add radionuclide dose breakdown for R2S

---
 docs/source/usersguide/decay_sources.rst | 13 +++++
 openmc/deplete/r2s.py                    | 31 +++++++++++-
 tests/unit_tests/test_r2s.py             | 64 ++++++++++++++++++++++++
 3 files changed, 107 insertions(+), 1 deletion(-)

diff --git a/docs/source/usersguide/decay_sources.rst b/docs/source/usersguide/decay_sources.rst
index 914fd278f67..5a961805f2d 100644
--- a/docs/source/usersguide/decay_sources.rst
+++ b/docs/source/usersguide/decay_sources.rst
@@ -141,6 +141,19 @@ time (after the 5 hour decay), we would run::
     r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes,
             photon_time_indices=[2])
 
+To obtain a dose-rate breakdown by radionuclide, set
+``nuclide_dose_breakdown=True``. This automatically adds a
+:class:`openmc.ParentNuclideFilter` to every photon tally that does not already
+have one. The filter bins are the union of radionuclides contributing to the
+prepared decay photon sources. The resulting bins can be used directly when
+inspecting the tally results::
+
+    r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes,
+            photon_time_indices=[2], nuclide_dose_breakdown=True)
+
+    dose_tally = r2s.results['photon_tallies'][2][0]
+    dose_by_nuclide = dose_tally.get_pandas_dataframe()
+
 After an R2S calculation has been run, the :class:`~openmc.deplete.R2SManager`
 instance will have a ``results`` dictionary that allows you to directly access
 results from each of the steps. It will also write out all the output files into
diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py
index 6cee45b1a9f..90347fafa40 100644
--- a/openmc/deplete/r2s.py
+++ b/openmc/deplete/r2s.py
@@ -156,6 +156,7 @@ def run(
         mat_vol_kwargs: dict | None = None,
         run_kwargs: dict | None = None,
         operator_kwargs: dict | None = None,
+        nuclide_dose_breakdown: bool = False,
     ):
         """Run the R2S calculation.
 
@@ -207,6 +208,10 @@ def run(
         operator_kwargs : dict, optional
             Additional keyword arguments passed to
             :class:`openmc.deplete.IndependentOperator`.
+        nuclide_dose_breakdown : bool, optional
+            Whether to add a :class:`~openmc.ParentNuclideFilter` to photon
+            tallies that do not already contain one. The filter bins are
+            determined automatically from the prepared decay photon sources.
 
         Returns
         -------
@@ -258,7 +263,8 @@ def run(
             )
             self.step3_photon_source(
                 photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
-                mat_vol_kwargs=mat_vol_kwargs
+                mat_vol_kwargs=mat_vol_kwargs,
+                nuclide_dose_breakdown=nuclide_dose_breakdown,
             )
             self.step4_photon_transport(
                 output_dir / 'photon_transport', run_kwargs=run_kwargs
@@ -447,6 +453,7 @@ def step3_photon_source(
         bounding_boxes: dict[int, openmc.BoundingBox] | None = None,
         output_dir: PathLike = 'photon_transport',
         mat_vol_kwargs: dict | None = None,
+        nuclide_dose_breakdown: bool = False,
     ):
         """Create decay photon sources.
 
@@ -474,6 +481,10 @@ def step3_photon_source(
         mat_vol_kwargs : dict, optional
             Additional keyword arguments passed to
             :meth:`openmc.MeshBase.material_volumes`.
+        nuclide_dose_breakdown : bool, optional
+            Whether to add a :class:`~openmc.ParentNuclideFilter` to photon
+            tallies that do not already contain one. The filter bins are
+            determined automatically from the prepared decay photon sources.
         """
 
         # TODO: Automatically determine bounding box for each cell
@@ -540,6 +551,14 @@ def step3_photon_source(
             time_index: self._create_photon_sources(time_index, work_items)
             for time_index in dict.fromkeys(time_indices)
         }
+        if nuclide_dose_breakdown:
+            radionuclides = sorted({
+                nuclide
+                for sources in self.results['photon_sources'].values()
+                for source in sources
+                for nuclide in source.energy.nuclides
+            })
+            self._add_parent_nuclide_filters(radionuclides)
 
     def step4_photon_transport(
         self,
@@ -640,6 +659,16 @@ def _get_mesh_work_items(self):
 
         return work_items
 
+    def _add_parent_nuclide_filters(self, radionuclides):
+        """Add parent nuclide filters to photon tallies when needed."""
+        if not radionuclides:
+            return
+
+        parent_filter = openmc.ParentNuclideFilter(radionuclides)
+        for tally in self.photon_model.tallies:
+            if not tally.contains_filter(openmc.ParentNuclideFilter):
+                tally.filters.append(parent_filter)
+
     def _create_photon_sources(self, time_index, work_items):
         """Create decay photon sources for a set of regions.
 
diff --git a/tests/unit_tests/test_r2s.py b/tests/unit_tests/test_r2s.py
index caf3cb8cd76..0e163626efc 100644
--- a/tests/unit_tests/test_r2s.py
+++ b/tests/unit_tests/test_r2s.py
@@ -169,6 +169,9 @@ def test_r2s_multi_mesh(simple_model_and_mesh, tmp_path):
 
 def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path):
     model, (c1, c2), _ = simple_model_and_mesh
+    tally = openmc.Tally()
+    tally.scores = ['flux']
+    model.tallies = [tally]
 
     # Use cell-based domains
     r2s = R2SManager(model, [c1, c2])
@@ -202,6 +205,8 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path):
     assert len(r2s.results['micros']) == 2
     assert len(r2s.results['activation_materials']) == 2
     assert len(r2s.results['depletion_results']) == 2
+    assert not r2s.photon_model.tallies[0].contains_filter(
+        openmc.ParentNuclideFilter)
 
     # Check activation materials
     amats = r2s.results['activation_materials']
@@ -227,3 +232,62 @@ def test_step4_requires_photon_sources(simple_model_and_mesh):
 
     with pytest.raises(RuntimeError, match='step3_photon_source'):
         r2s.step4_photon_transport()
+
+
+def test_r2s_nuclide_dose_breakdown(simple_model_and_mesh, tmp_path, monkeypatch):
+    model, (c1, c2), _ = simple_model_and_mesh
+    tally = openmc.Tally()
+    tally.scores = ['flux']
+    model.tallies = [tally]
+    r2s = R2SManager(model, [c1, c2])
+    chain = Chain.from_xml(Path(__file__).parents[1] / 'chain_ni.xml')
+    source_counts = {}
+    create_sources = r2s._create_photon_sources
+
+    def count_source_creation(time_index, work_items):
+        source_counts[time_index] = source_counts.get(time_index, 0) + 1
+        return create_sources(time_index, work_items)
+
+    monkeypatch.setattr(r2s, '_create_photon_sources', count_source_creation)
+
+    r2s.run(
+        timesteps=[(1.0, 'd'), (1.0, 'd')],
+        source_rates=[1.0, 0.0],
+        photon_time_indices=[1, 2],
+        nuclide_dose_breakdown=True,
+        output_dir=tmp_path,
+        bounding_boxes={c1.id: c1.bounding_box, c2.id: c2.bounding_box},
+        chain_file=chain,
+    )
+
+    filters = [
+        filter for filter in r2s.photon_model.tallies[0].filters
+        if isinstance(filter, openmc.ParentNuclideFilter)
+    ]
+    assert len(filters) == 1
+    assert list(filters[0].bins) == sorted(filters[0].bins)
+    assert filters[0].bins.size > 0
+    assert source_counts == {1: 1, 2: 1}
+
+    for time_index in (1, 2):
+        result_tally = r2s.results['photon_tallies'][time_index][0]
+        result_filter = next(
+            filter for filter in result_tally.filters
+            if isinstance(filter, openmc.ParentNuclideFilter)
+        )
+        assert list(result_filter.bins) == list(filters[0].bins)
+
+
+def test_r2s_preserves_parent_nuclide_filters(simple_model_and_mesh):
+    model, (c1, c2), _ = simple_model_and_mesh
+    tally = openmc.Tally()
+    tally.filters = [openmc.ParentNuclideFilter(['Co60'])]
+    unfiltered_tally = openmc.Tally()
+    model.tallies = [tally, unfiltered_tally]
+    r2s = R2SManager(model, [c1, c2])
+
+    r2s._add_parent_nuclide_filters(['Co60', 'Ni65'])
+
+    assert len(r2s.photon_model.tallies[0].filters) == 1
+    assert list(r2s.photon_model.tallies[0].filters[0].bins) == ['Co60']
+    assert list(r2s.photon_model.tallies[1].filters[0].bins) == ['Co60', 'Ni65']

From 1abeb91993c5038f03d8ee3f8df80f490e69f55d Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sun, 12 Jul 2026 16:33:01 -0500
Subject: [PATCH 3/7] Handle R2S timesteps at which there is no decay photon
 source

---
 docs/source/usersguide/decay_sources.rst |  5 +-
 openmc/deplete/r2s.py                    | 68 ++++++++++++++---
 tests/unit_tests/test_r2s.py             | 93 +++++++++++++++++++++++-
 3 files changed, 150 insertions(+), 16 deletions(-)

diff --git a/docs/source/usersguide/decay_sources.rst b/docs/source/usersguide/decay_sources.rst
index 5a961805f2d..c7ebf383449 100644
--- a/docs/source/usersguide/decay_sources.rst
+++ b/docs/source/usersguide/decay_sources.rst
@@ -132,8 +132,9 @@ can be run::
     r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes)
 
 If not specified otherwise, a photon transport calculation is run at each time
-in the depletion schedule. That means in the case above, we would see three
-photon transport calculations. To specify specific times at which photon
+in the depletion schedule for which a decay photon source exists. Times without
+a decay photon source, such as the initial state of a model containing only
+stable nuclides, are omitted. To specify particular times at which photon
 transport calculations should be run, pass the ``photon_time_indices`` argument.
 For example, if we wanted to run a photon transport calculation only on the last
 time (after the 5 hour decay), we would run::
diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py
index 90347fafa40..dff11e6c93b 100644
--- a/openmc/deplete/r2s.py
+++ b/openmc/deplete/r2s.py
@@ -4,6 +4,7 @@
 import copy
 from datetime import datetime
 import json
+from numbers import Integral
 from pathlib import Path
 
 import numpy as np
@@ -180,7 +181,7 @@ def run(
             timesteps. For example, if two timesteps are specified, the array of
             times would contain three entries, and [2] would indicate computing
             photon results at the last time. A value of None indicates to run
-            photon transport for each time.
+            photon transport at each time that has a decay photon source.
         output_dir : PathLike, optional
             Path to directory where R2S calculation outputs will be saved. If
             not provided, a timestamped directory 'r2s_YYYY-MM-DDTHH-MM-SS' is
@@ -471,7 +472,7 @@ def step3_photon_source(
             timesteps. For example, if two timesteps are specified, the array of
             times would contain three entries, and [2] would indicate computing
             photon results at the last time. A value of None indicates to run
-            photon transport for each time.
+            photon transport at each time that has a decay photon source.
         bounding_boxes : dict[int, openmc.BoundingBox], optional
             Dictionary mapping cell IDs to bounding boxes used for spatial
             source sampling in cell-based R2S calculations. Required if method
@@ -487,6 +488,10 @@ def step3_photon_source(
             determined automatically from the prepared decay photon sources.
         """
 
+        # Do not retain sources from an earlier successful call if this source
+        # preparation attempt fails.
+        self.results.pop('photon_sources', None)
+
         # TODO: Automatically determine bounding box for each cell
         if bounding_boxes is None and self.method == 'cell-based':
             raise ValueError("bounding_boxes must be provided for cell-based "
@@ -495,12 +500,30 @@ def step3_photon_source(
         output_dir = Path(output_dir)
         output_dir.mkdir(parents=True, exist_ok=True)
 
-        # Get default time indices if not provided
+        # Determine and validate time indices before preparing source data.
+        n_steps = len(self.results['depletion_results'])
+        implicit_time_indices = time_indices is None
         if time_indices is None:
-            n_steps = len(self.results['depletion_results'])
             time_indices = list(range(n_steps))
         else:
             time_indices = list(time_indices)
+            if not time_indices:
+                raise ValueError('time_indices must contain at least one index')
+
+            normalized_indices = []
+            for index in time_indices:
+                if isinstance(index, bool) or not isinstance(index, Integral):
+                    raise TypeError('time_indices must contain only integers')
+                index = int(index)
+                if index < -n_steps or index >= n_steps:
+                    raise IndexError(
+                        f'Photon time index {index} is out of range for '
+                        f'{n_steps} depletion results')
+                normalized_index = index + n_steps if index < 0 else index
+                normalized_indices.append(normalized_index)
+
+            # Remove duplicates while preserving order
+            time_indices = list(dict.fromkeys(normalized_indices))
 
         # Check whether the photon model is different
         neutron_univ = self.neutron_model.geometry.root_universe
@@ -543,14 +566,33 @@ def step3_photon_source(
                         continue
                 work_items.append((cell, original_mat, bounding_boxes[cell.id]))
 
-        # Normalize negative time indices before creating source lists.
-        n_steps = len(self.results['depletion_results'])
-        time_indices = [i + n_steps if i < 0 else i for i in time_indices]
-
-        self.results['photon_sources'] = {
+        # Create decay photon sources for each time index
+        photon_sources = {
             time_index: self._create_photon_sources(time_index, work_items)
-            for time_index in dict.fromkeys(time_indices)
+            for time_index in time_indices
         }
+
+        # Determine if any times have no decay photon sources. If the user
+        # didn't specify any specific time indices, remove those times from the
+        # photon_sources dictionary. If the user did specify time indices, raise
+        # an error if any of those times have no decay photon sources.
+        empty_indices = [
+            time_index for time_index, sources in photon_sources.items()
+            if not sources
+        ]
+        if implicit_time_indices:
+            for time_index in empty_indices:
+                del photon_sources[time_index]
+            if not photon_sources:
+                raise RuntimeError(
+                    'No decay photon sources were found at any depletion time')
+        elif empty_indices:
+            indices = ', '.join(str(index) for index in empty_indices)
+            raise RuntimeError(
+                f'No decay photon source was found for requested time '
+                f'indices: {indices}')
+
+        self.results['photon_sources'] = photon_sources
         if nuclide_dose_breakdown:
             radionuclides = sorted({
                 nuclide
@@ -583,6 +625,10 @@ def step4_photon_transport(
             raise RuntimeError(
                 'Photon sources must be created with step3_photon_source '
                 'before running photon transport.')
+        photon_sources = self.results['photon_sources']
+        if not photon_sources:
+            raise RuntimeError(
+                'No decay photon sources are available for transport')
 
         if run_kwargs is None:
             run_kwargs = {}
@@ -601,7 +647,7 @@ def step4_photon_transport(
         # Ensure photon transport is enabled in settings.
         self.photon_model.settings.photon_transport = True
 
-        for time_index, sources in self.results['photon_sources'].items():
+        for time_index, sources in photon_sources.items():
             self.photon_model.settings.source = sources
 
             # Run photon transport calculation
diff --git a/tests/unit_tests/test_r2s.py b/tests/unit_tests/test_r2s.py
index 0e163626efc..a6e212a5648 100644
--- a/tests/unit_tests/test_r2s.py
+++ b/tests/unit_tests/test_r2s.py
@@ -46,6 +46,16 @@ def simple_model_and_mesh():
     return model, (c1, c2), mesh
 
 
+@pytest.fixture
+def source_stage_manager(simple_model_and_mesh):
+    model, (c1, c2), _ = simple_model_and_mesh
+    r2s = R2SManager(model, [c1, c2])
+    r2s.results['depletion_results'] = [None, None]
+    r2s.results['activation_materials'] = [c1.fill, c2.fill]
+    bounding_boxes = {c1.id: c1.bounding_box, c2.id: c2.bounding_box}
+    return r2s, bounding_boxes
+
+
 def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
     model, (c1, c2), mesh = simple_model_and_mesh
 
@@ -59,7 +69,6 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
     outdir = r2s.run(
         timesteps=[(1.0, 'd')],
         source_rates=[1.0],
-        photon_time_indices=[1],
         output_dir=tmp_path,
         chain_file=chain,
     )
@@ -73,6 +82,7 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
     assert (act / 'depletion_results.h5').exists()
     pt = Path(outdir) / 'photon_transport'
     assert (pt / 'tally_ids.json').exists()
+    assert not (pt / 'time_0').exists()
     assert (pt / 'time_1' / 'statepoint.10.h5').exists()
 
     # Basic results structure checks
@@ -226,12 +236,89 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path):
     assert len(r2s_loaded.results['depletion_results']) == 2
 
 
-def test_step4_requires_photon_sources(simple_model_and_mesh):
+def test_step4_requires_photon_sources(simple_model_and_mesh, tmp_path):
     model, (c1, c2), _ = simple_model_and_mesh
     r2s = R2SManager(model, [c1, c2])
+    output_dir = tmp_path / 'photon'
 
     with pytest.raises(RuntimeError, match='step3_photon_source'):
-        r2s.step4_photon_transport()
+        r2s.step4_photon_transport(output_dir)
+
+    r2s.results['photon_sources'] = {}
+    with pytest.raises(RuntimeError, match='No decay photon sources'):
+        r2s.step4_photon_transport(output_dir)
+
+    assert not output_dir.exists()
+
+
+def test_default_photon_times_skip_empty_sources(
+    source_stage_manager, tmp_path, monkeypatch
+):
+    r2s, bounding_boxes = source_stage_manager
+    source = object()
+    sources_by_time = {0: [], 1: [source]}
+    monkeypatch.setattr(
+        r2s, '_create_photon_sources',
+        lambda time_index, work_items: sources_by_time[time_index])
+
+    r2s.step3_photon_source(
+        bounding_boxes=bounding_boxes, output_dir=tmp_path)
+
+    assert r2s.results['photon_sources'] == {1: [source]}
+
+
+def test_explicit_empty_photon_source_fails(
+    source_stage_manager, tmp_path, monkeypatch
+):
+    r2s, bounding_boxes = source_stage_manager
+    source = object()
+    sources_by_time = {0: [], 1: [source]}
+    monkeypatch.setattr(
+        r2s, '_create_photon_sources',
+        lambda time_index, work_items: sources_by_time[time_index])
+    r2s.results['photon_sources'] = {99: [source]}
+
+    with pytest.raises(RuntimeError, match='requested time indices: 0'):
+        r2s.step3_photon_source(
+            [0, 1], bounding_boxes, output_dir=tmp_path)
+
+    assert 'photon_sources' not in r2s.results
+
+
+def test_default_photon_times_require_a_source(
+    source_stage_manager, tmp_path, monkeypatch
+):
+    r2s, bounding_boxes = source_stage_manager
+    monkeypatch.setattr(
+        r2s, '_create_photon_sources',
+        lambda time_index, work_items: [])
+
+    with pytest.raises(RuntimeError, match='at any depletion time'):
+        r2s.step3_photon_source(
+            bounding_boxes=bounding_boxes, output_dir=tmp_path)
+
+    assert 'photon_sources' not in r2s.results
+
+
+@pytest.mark.parametrize(
+    ('time_indices', 'exception'),
+    [
+        ([], ValueError),
+        ([2], IndexError),
+        ([-3], IndexError),
+        ([1.0], TypeError),
+    ],
+)
+def test_photon_time_index_validation(
+    source_stage_manager, tmp_path, time_indices, exception
+):
+    r2s, bounding_boxes = source_stage_manager
+
+    with pytest.raises(exception):
+        r2s.step3_photon_source(
+            time_indices, bounding_boxes, output_dir=tmp_path)
+
+    assert 'photon_sources' not in r2s.results
 
 
 def test_r2s_nuclide_dose_breakdown(simple_model_and_mesh, tmp_path, monkeypatch):

From 3f7c6406ce1ee81de8414dee6068946c8cffcdbd Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sun, 12 Jul 2026 17:01:55 -0500
Subject: [PATCH 4/7] Rename to by_parent_nuclide. Join photon source loops

---
 docs/source/usersguide/decay_sources.rst | 10 +--
 openmc/deplete/r2s.py                    | 89 ++++++++++++------------
 tests/unit_tests/test_r2s.py             |  4 +-
 3 files changed, 51 insertions(+), 52 deletions(-)

diff --git a/docs/source/usersguide/decay_sources.rst b/docs/source/usersguide/decay_sources.rst
index c7ebf383449..19f1e5d5f94 100644
--- a/docs/source/usersguide/decay_sources.rst
+++ b/docs/source/usersguide/decay_sources.rst
@@ -142,18 +142,18 @@ time (after the 5 hour decay), we would run::
     r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes,
             photon_time_indices=[2])
 
-To obtain a dose-rate breakdown by radionuclide, set
-``nuclide_dose_breakdown=True``. This automatically adds a
+To attribute photon tally results to their parent radionuclides, set
+``by_parent_nuclide=True``. This automatically adds a
 :class:`openmc.ParentNuclideFilter` to every photon tally that does not already
 have one. The filter bins are the union of radionuclides contributing to the
 prepared decay photon sources. The resulting bins can be used directly when
 inspecting the tally results::
 
     r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes,
-            photon_time_indices=[2], nuclide_dose_breakdown=True)
+            photon_time_indices=[2], by_parent_nuclide=True)
 
-    dose_tally = r2s.results['photon_tallies'][2][0]
-    dose_by_nuclide = dose_tally.get_pandas_dataframe()
+    photon_tally = r2s.results['photon_tallies'][2][0]
+    tally_by_parent = photon_tally.get_pandas_dataframe()
 
 After an R2S calculation has been run, the :class:`~openmc.deplete.R2SManager`
 instance will have a ``results`` dictionary that allows you to directly access
diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py
index dff11e6c93b..a07e1202bc6 100644
--- a/openmc/deplete/r2s.py
+++ b/openmc/deplete/r2s.py
@@ -157,7 +157,7 @@ def run(
         mat_vol_kwargs: dict | None = None,
         run_kwargs: dict | None = None,
         operator_kwargs: dict | None = None,
-        nuclide_dose_breakdown: bool = False,
+        by_parent_nuclide: bool = False,
     ):
         """Run the R2S calculation.
 
@@ -209,10 +209,11 @@ def run(
         operator_kwargs : dict, optional
             Additional keyword arguments passed to
             :class:`openmc.deplete.IndependentOperator`.
-        nuclide_dose_breakdown : bool, optional
-            Whether to add a :class:`~openmc.ParentNuclideFilter` to photon
-            tallies that do not already contain one. The filter bins are
-            determined automatically from the prepared decay photon sources.
+        by_parent_nuclide : bool, optional
+            Whether to score photon tallies separately for each parent
+            radionuclide. A :class:`~openmc.ParentNuclideFilter` is added to
+            tallies that do not already contain one, with bins determined from
+            the prepared decay photon sources.
 
         Returns
         -------
@@ -265,7 +266,7 @@ def run(
             self.step3_photon_source(
                 photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
                 mat_vol_kwargs=mat_vol_kwargs,
-                nuclide_dose_breakdown=nuclide_dose_breakdown,
+                by_parent_nuclide=by_parent_nuclide,
             )
             self.step4_photon_transport(
                 output_dir / 'photon_transport', run_kwargs=run_kwargs
@@ -454,7 +455,7 @@ def step3_photon_source(
         bounding_boxes: dict[int, openmc.BoundingBox] | None = None,
         output_dir: PathLike = 'photon_transport',
         mat_vol_kwargs: dict | None = None,
-        nuclide_dose_breakdown: bool = False,
+        by_parent_nuclide: bool = False,
     ):
         """Create decay photon sources.
 
@@ -482,10 +483,11 @@ def step3_photon_source(
         mat_vol_kwargs : dict, optional
             Additional keyword arguments passed to
             :meth:`openmc.MeshBase.material_volumes`.
-        nuclide_dose_breakdown : bool, optional
-            Whether to add a :class:`~openmc.ParentNuclideFilter` to photon
-            tallies that do not already contain one. The filter bins are
-            determined automatically from the prepared decay photon sources.
+        by_parent_nuclide : bool, optional
+            Whether to score photon tallies separately for each parent
+            radionuclide. A :class:`~openmc.ParentNuclideFilter` is added to
+            tallies that do not already contain one, with bins determined from
+            the prepared decay photon sources.
         """
 
         # Do not retain sources from an earlier successful call if this source
@@ -593,7 +595,7 @@ def step3_photon_source(
                 f'indices: {indices}')
 
         self.results['photon_sources'] = photon_sources
-        if nuclide_dose_breakdown:
+        if by_parent_nuclide:
             radionuclides = sorted({
                 nuclide
                 for sources in self.results['photon_sources'].values()
@@ -740,44 +742,41 @@ def _create_photon_sources(self, time_index, work_items):
             self.neutron_model._get_all_materials()
             if self.method == 'mesh-based' else None
         )
+        step_result = self.results['depletion_results'][time_index]
         sources = []
         for item in work_items:
-            source = self._create_photon_source(time_index, item, mat_dict)
-            if source is not None:
-                sources.append(source)
+            if self.method == 'mesh-based':
+                index_mat, domain_id, bbox = item
+                original_mat = self.results['activation_materials'][index_mat]
+                domain = mat_dict[domain_id]
+            else:
+                domain, original_mat, bbox = item
+
+            activated_mat = step_result.get_material(str(original_mat.id))
+            nuclides = activated_mat.get_nuclide_atom_densities()
+            if not nuclides:
+                continue
+
+            # Eliminate nuclides with zero density.
+            nuclides = {
+                nuclide: density for nuclide, density in nuclides.items()
+                if density > 0
+            }
+            energy = openmc.stats.DecaySpectrum(
+                nuclides, activated_mat.volume)
+            energy.clip(inplace=True)
+            if not energy.nuclides:
+                continue
+
+            sources.append(openmc.IndependentSource(
+                space=openmc.stats.Box(bbox.lower_left, bbox.upper_right),
+                energy=energy,
+                particle='photon',
+                constraints={'domains': [domain]},
+            ))
 
         return sources
 
-    def _create_photon_source(self, time_index, item, mat_dict=None):
-        """Create a decay photon source for one activation region."""
-        step_result = self.results['depletion_results'][time_index]
-        if self.method == 'mesh-based':
-            index_mat, domain_id, bbox = item
-            original_mat = self.results['activation_materials'][index_mat]
-            domain = mat_dict[domain_id]
-        else:
-            domain, original_mat, bbox = item
-
-        activated_mat = step_result.get_material(str(original_mat.id))
-        nuclides = activated_mat.get_nuclide_atom_densities()
-        if not nuclides:
-            return None
-
-        # Eliminate nuclides with zero density.
-        nuclides = {nuclide: density for nuclide, density in nuclides.items()
-                    if density > 0}
-        energy = openmc.stats.DecaySpectrum(nuclides, activated_mat.volume)
-        energy.clip(inplace=True)
-        if not energy.nuclides:
-            return None
-
-        return openmc.IndependentSource(
-            space=openmc.stats.Box(bbox.lower_left, bbox.upper_right),
-            energy=energy,
-            particle='photon',
-            constraints={'domains': [domain]},
-        )
-
     def load_results(self, path: PathLike):
         """Load results from a previous R2S calculation.
 
diff --git a/tests/unit_tests/test_r2s.py b/tests/unit_tests/test_r2s.py
index a6e212a5648..6ff9397dd63 100644
--- a/tests/unit_tests/test_r2s.py
+++ b/tests/unit_tests/test_r2s.py
@@ -321,7 +321,7 @@ def test_photon_time_index_validation(
     assert 'photon_sources' not in r2s.results
 
 
-def test_r2s_nuclide_dose_breakdown(simple_model_and_mesh, tmp_path, monkeypatch):
+def test_r2s_by_parent_nuclide(simple_model_and_mesh, tmp_path, monkeypatch):
     model, (c1, c2), _ = simple_model_and_mesh
     tally = openmc.Tally()
     tally.scores = ['flux']
@@ -341,7 +341,7 @@ def count_source_creation(time_index, work_items):
         timesteps=[(1.0, 'd'), (1.0, 'd')],
         source_rates=[1.0, 0.0],
         photon_time_indices=[1, 2],
-        nuclide_dose_breakdown=True,
+        by_parent_nuclide=True,
         output_dir=tmp_path,
         bounding_boxes={c1.id: c1.bounding_box, c2.id: c2.bounding_box},
         chain_file=chain,

From d6c547d8277a8327502f6a2649d52d5e6dcc7f3e Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sun, 12 Jul 2026 17:10:14 -0500
Subject: [PATCH 5/7] Move by_parent_nuclide to step4_photon_transport

---
 openmc/deplete/r2s.py        | 48 +++++++++++++++++-------------------
 tests/unit_tests/test_r2s.py | 15 -----------
 2 files changed, 22 insertions(+), 41 deletions(-)

diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py
index a07e1202bc6..fc7edc5e7e2 100644
--- a/openmc/deplete/r2s.py
+++ b/openmc/deplete/r2s.py
@@ -266,10 +266,10 @@ def run(
             self.step3_photon_source(
                 photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
                 mat_vol_kwargs=mat_vol_kwargs,
-                by_parent_nuclide=by_parent_nuclide,
             )
             self.step4_photon_transport(
-                output_dir / 'photon_transport', run_kwargs=run_kwargs
+                output_dir / 'photon_transport', run_kwargs=run_kwargs,
+                by_parent_nuclide=by_parent_nuclide,
             )
 
         return output_dir
@@ -455,7 +455,6 @@ def step3_photon_source(
         bounding_boxes: dict[int, openmc.BoundingBox] | None = None,
         output_dir: PathLike = 'photon_transport',
         mat_vol_kwargs: dict | None = None,
-        by_parent_nuclide: bool = False,
     ):
         """Create decay photon sources.
 
@@ -483,11 +482,6 @@ def step3_photon_source(
         mat_vol_kwargs : dict, optional
             Additional keyword arguments passed to
             :meth:`openmc.MeshBase.material_volumes`.
-        by_parent_nuclide : bool, optional
-            Whether to score photon tallies separately for each parent
-            radionuclide. A :class:`~openmc.ParentNuclideFilter` is added to
-            tallies that do not already contain one, with bins determined from
-            the prepared decay photon sources.
         """
 
         # Do not retain sources from an earlier successful call if this source
@@ -595,19 +589,12 @@ def step3_photon_source(
                 f'indices: {indices}')
 
         self.results['photon_sources'] = photon_sources
-        if by_parent_nuclide:
-            radionuclides = sorted({
-                nuclide
-                for sources in self.results['photon_sources'].values()
-                for source in sources
-                for nuclide in source.energy.nuclides
-            })
-            self._add_parent_nuclide_filters(radionuclides)
 
     def step4_photon_transport(
         self,
         output_dir: PathLike = 'photon_transport',
         run_kwargs: dict | None = None,
+        by_parent_nuclide: bool = False,
     ):
         """Run photon transport using prepared decay photon sources.
 
@@ -622,6 +609,11 @@ def step4_photon_transport(
         run_kwargs : dict, optional
             Additional keyword arguments passed to :meth:`openmc.Model.run`.
             By default, output is disabled.
+        by_parent_nuclide : bool, optional
+            Whether to score photon tallies separately for each parent
+            radionuclide. A :class:`~openmc.ParentNuclideFilter` is added to
+            tallies that do not already contain one, with bins determined from
+            the prepared decay photon sources.
         """
         if 'photon_sources' not in self.results:
             raise RuntimeError(
@@ -632,6 +624,20 @@ def step4_photon_transport(
             raise RuntimeError(
                 'No decay photon sources are available for transport')
 
+        if by_parent_nuclide:
+            radionuclides = sorted({
+                nuclide
+                for sources in photon_sources.values()
+                for source in sources
+                for nuclide in source.energy.nuclides
+            })
+
+            if radionuclides:
+                parent_filter = openmc.ParentNuclideFilter(radionuclides)
+                for tally in self.photon_model.tallies:
+                    if not tally.contains_filter(openmc.ParentNuclideFilter):
+                        tally.filters.append(parent_filter)
+
         if run_kwargs is None:
             run_kwargs = {}
         run_kwargs.setdefault('output', False)
@@ -707,16 +713,6 @@ def _get_mesh_work_items(self):
 
         return work_items
 
-    def _add_parent_nuclide_filters(self, radionuclides):
-        """Add parent nuclide filters to photon tallies when needed."""
-        if not radionuclides:
-            return
-
-        parent_filter = openmc.ParentNuclideFilter(radionuclides)
-        for tally in self.photon_model.tallies:
-            if not tally.contains_filter(openmc.ParentNuclideFilter):
-                tally.filters.append(parent_filter)
-
     def _create_photon_sources(self, time_index, work_items):
         """Create decay photon sources for a set of regions.
 
diff --git a/tests/unit_tests/test_r2s.py b/tests/unit_tests/test_r2s.py
index 6ff9397dd63..1f4520db548 100644
--- a/tests/unit_tests/test_r2s.py
+++ b/tests/unit_tests/test_r2s.py
@@ -363,18 +363,3 @@ def count_source_creation(time_index, work_items):
             if isinstance(filter, openmc.ParentNuclideFilter)
         )
         assert list(result_filter.bins) == list(filters[0].bins)
-
-
-def test_r2s_preserves_parent_nuclide_filters(simple_model_and_mesh):
-    model, (c1, c2), _ = simple_model_and_mesh
-    tally = openmc.Tally()
-    tally.filters = [openmc.ParentNuclideFilter(['Co60'])]
-    unfiltered_tally = openmc.Tally()
-    model.tallies = [tally, unfiltered_tally]
-    r2s = R2SManager(model, [c1, c2])
-
-    r2s._add_parent_nuclide_filters(['Co60', 'Ni65'])
-
-    assert len(r2s.photon_model.tallies[0].filters) == 1
-    assert list(r2s.photon_model.tallies[0].filters[0].bins) == ['Co60']
-    assert list(r2s.photon_model.tallies[1].filters[0].bins) == ['Co60', 'Ni65']

From 7a58205f4951417d1c16cf1ad61273a53bf1601f Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sun, 12 Jul 2026 17:13:04 -0500
Subject: [PATCH 6/7] Preserve original code in _create_photon_sources

---
 openmc/deplete/r2s.py | 28 ++++++++++++++--------------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py
index fc7edc5e7e2..e04297d973f 100644
--- a/openmc/deplete/r2s.py
+++ b/openmc/deplete/r2s.py
@@ -734,32 +734,32 @@ def _create_photon_sources(self, time_index, work_items):
         list of openmc.IndependentSource
             Photon sources for each activated region.
         """
-        mat_dict = (
-            self.neutron_model._get_all_materials()
-            if self.method == 'mesh-based' else None
-        )
         step_result = self.results['depletion_results'][time_index]
+        materials = self.results['activation_materials']
+        mesh_based = self.method == 'mesh-based'
+        if mesh_based:
+            mat_dict = self.neutron_model._get_all_materials()
+
         sources = []
         for item in work_items:
-            if self.method == 'mesh-based':
+            if mesh_based:
                 index_mat, domain_id, bbox = item
-                original_mat = self.results['activation_materials'][index_mat]
+                original_mat = materials[index_mat]
                 domain = mat_dict[domain_id]
             else:
-                domain, original_mat, bbox = item
+                cell, original_mat, bbox = item
+                domain = cell
 
             activated_mat = step_result.get_material(str(original_mat.id))
             nuclides = activated_mat.get_nuclide_atom_densities()
             if not nuclides:
                 continue
 
-            # Eliminate nuclides with zero density.
-            nuclides = {
-                nuclide: density for nuclide, density in nuclides.items()
-                if density > 0
-            }
-            energy = openmc.stats.DecaySpectrum(
-                nuclides, activated_mat.volume)
+            # Eliminate nuclides with zero density
+            nuclides = {nuclide: density for nuclide, density in nuclides.items()
+                        if density > 0}
+
+            energy = openmc.stats.DecaySpectrum(nuclides, activated_mat.volume)
             energy.clip(inplace=True)
             if not energy.nuclides:
                 continue

From afe5a289ec3404351f9a9fb9c2762cc1d7f2e27e Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sun, 12 Jul 2026 17:36:46 -0500
Subject: [PATCH 7/7] Simplify/speed up tests

---
 tests/unit_tests/test_r2s.py | 62 +++++++-----------------------------
 1 file changed, 11 insertions(+), 51 deletions(-)

diff --git a/tests/unit_tests/test_r2s.py b/tests/unit_tests/test_r2s.py
index 1f4520db548..5e617919d9c 100644
--- a/tests/unit_tests/test_r2s.py
+++ b/tests/unit_tests/test_r2s.py
@@ -33,8 +33,8 @@ def simple_model_and_mesh():
 
     # Simple settings with a point source
     settings = openmc.Settings()
-    settings.batches = 10
-    settings.particles = 1000
+    settings.batches = 2
+    settings.particles = 250
     settings.run_mode = 'fixed source'
     settings.source = openmc.IndependentSource()
     model = openmc.Model(geometry, settings=settings)
@@ -71,6 +71,7 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
         source_rates=[1.0],
         output_dir=tmp_path,
         chain_file=chain,
+        micro_kwargs={'nuclides': ['Ni58'], 'reactions': ['(n,p)']},
     )
 
     # Check directories and files exist
@@ -83,7 +84,7 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
     pt = Path(outdir) / 'photon_transport'
     assert (pt / 'tally_ids.json').exists()
     assert not (pt / 'time_0').exists()
-    assert (pt / 'time_1' / 'statepoint.10.h5').exists()
+    assert (pt / 'time_1' / 'statepoint.2.h5').exists()
 
     # Basic results structure checks
     assert len(r2s.results['fluxes']) == 2
@@ -137,6 +138,7 @@ def test_r2s_multi_mesh(simple_model_and_mesh, tmp_path):
         photon_time_indices=[1],
         output_dir=tmp_path,
         chain_file=chain,
+        micro_kwargs={'nuclides': ['Ni58'], 'reactions': ['(n,p)']},
     )
 
     # Check that per-mesh MMV files were written
@@ -149,7 +151,7 @@ def test_r2s_multi_mesh(simple_model_and_mesh, tmp_path):
     assert (act / 'depletion_results.h5').exists()
     pt = Path(outdir) / 'photon_transport'
     assert (pt / 'tally_ids.json').exists()
-    assert (pt / 'time_1' / 'statepoint.10.h5').exists()
+    assert (pt / 'time_1' / 'statepoint.2.h5').exists()
 
     # Two meshes, each with 1 element containing both materials →
     # 2 element-material combinations per mesh, 4 total
@@ -195,9 +197,11 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path):
         timesteps=[(1.0, 'd')],
         source_rates=[1.0],
         photon_time_indices=[1],
+        by_parent_nuclide=True,
         output_dir=tmp_path,
         bounding_boxes=bounding_boxes,
-        chain_file=chain
+        chain_file=chain,
+        micro_kwargs={'nuclides': ['Ni58'], 'reactions': ['(n,p)']},
     )
 
     # Check directories and files exist
@@ -208,14 +212,14 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path):
     assert (act / 'depletion_results.h5').exists()
     pt = Path(outdir) / 'photon_transport'
     assert (pt / 'tally_ids.json').exists()
-    assert (pt / 'time_1' / 'statepoint.10.h5').exists()
+    assert (pt / 'time_1' / 'statepoint.2.h5').exists()
 
     # Basic results structure checks
     assert len(r2s.results['fluxes']) == 2
     assert len(r2s.results['micros']) == 2
     assert len(r2s.results['activation_materials']) == 2
     assert len(r2s.results['depletion_results']) == 2
-    assert not r2s.photon_model.tallies[0].contains_filter(
+    assert r2s.photon_model.tallies[0].contains_filter(
         openmc.ParentNuclideFilter)
 
     # Check activation materials
@@ -319,47 +323,3 @@ def test_photon_time_index_validation(
             time_indices, bounding_boxes, output_dir=tmp_path)
 
     assert 'photon_sources' not in r2s.results
-
-
-def test_r2s_by_parent_nuclide(simple_model_and_mesh, tmp_path, monkeypatch):
-    model, (c1, c2), _ = simple_model_and_mesh
-    tally = openmc.Tally()
-    tally.scores = ['flux']
-    model.tallies = [tally]
-    r2s = R2SManager(model, [c1, c2])
-    chain = Chain.from_xml(Path(__file__).parents[1] / 'chain_ni.xml')
-    source_counts = {}
-    create_sources = r2s._create_photon_sources
-
-    def count_source_creation(time_index, work_items):
-        source_counts[time_index] = source_counts.get(time_index, 0) + 1
-        return create_sources(time_index, work_items)
-
-    monkeypatch.setattr(r2s, '_create_photon_sources', count_source_creation)
-
-    r2s.run(
-        timesteps=[(1.0, 'd'), (1.0, 'd')],
-        source_rates=[1.0, 0.0],
-        photon_time_indices=[1, 2],
-        by_parent_nuclide=True,
-        output_dir=tmp_path,
-        bounding_boxes={c1.id: c1.bounding_box, c2.id: c2.bounding_box},
-        chain_file=chain,
-    )
-
-    filters = [
-        filter for filter in r2s.photon_model.tallies[0].filters
-        if isinstance(filter, openmc.ParentNuclideFilter)
-    ]
-    assert len(filters) == 1
-    assert list(filters[0].bins) == sorted(filters[0].bins)
-    assert filters[0].bins.size > 0
-    assert source_counts == {1: 1, 2: 1}
-
-    for time_index in (1, 2):
-        result_tally = r2s.results['photon_tallies'][time_index][0]
-        result_filter = next(
-            filter for filter in result_tally.filters
-            if isinstance(filter, openmc.ParentNuclideFilter)
-        )
-        assert list(result_filter.bins) == list(filters[0].bins)