From d1878a205e8020bc72d06c0f5f64503fdd808929 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Fri, 20 Mar 2020 18:07:17 -0400
Subject: [PATCH 01/10] Replace "Nothing" depletion targets with None

This change is more declarative and mitigates any changes
that there are errors if a target is Nothing, None, or their lower
cased varieties. This is helpful as the depletion chain gets
reduced down and we want to keep total removal rates for
isotopes that are present, but their products don't exist
in the new reduced chain.

Chain.form_matrix has been updated to handle targets that
are None for non-fission reactions and for decay reactions.
Documentation for ReactionTuple and DecayTuple have been
updated to indicate that the target attribute could be None.
---
 openmc/deplete/chain.py   | 23 +++++++++++------------
 openmc/deplete/nuclide.py | 22 ++++++++++++++++------
 2 files changed, 27 insertions(+), 18 deletions(-)

diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py
index 814de12fc4..8c31914acb 100644
--- a/openmc/deplete/chain.py
+++ b/openmc/deplete/chain.py
@@ -500,7 +500,7 @@ class Chain(object):
                 # Gain
                 for _, target, branching_ratio in nuc.decay_modes:
                     # Allow for total annihilation for debug purposes
-                    if target != 'Nothing':
+                    if target is not None:
                         branch_val = branching_ratio * decay_constant
 
                         if branch_val != 0.0:
@@ -525,17 +525,16 @@ class Chain(object):
                             matrix[i, i] -= path_rate
 
                     # Gain term; allow for total annihilation for debug purposes
-                    if target != 'Nothing':
-                        if r_type != 'fission':
-                            if path_rate != 0.0:
-                                k = self.nuclide_dict[target]
-                                matrix[k, i] += path_rate * br
-                        else:
-                            for product, y in fission_yields[nuc.name].items():
-                                yield_val = y * path_rate
-                                if yield_val != 0.0:
-                                    k = self.nuclide_dict[product]
-                                    matrix[k, i] += yield_val
+                    if r_type != 'fission':
+                        if target is not None and path_rate != 0.0:
+                            k = self.nuclide_dict[target]
+                            matrix[k, i] += path_rate * br
+                    else:
+                        for product, y in fission_yields[nuc.name].items():
+                            yield_val = y * path_rate
+                            if yield_val != 0.0:
+                                k = self.nuclide_dict[product]
+                                matrix[k, i] += yield_val
 
                 # Clear set of reactions
                 reactions.clear()
diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py
index 0aae5adb72..7235745aa2 100644
--- a/openmc/deplete/nuclide.py
+++ b/openmc/deplete/nuclide.py
@@ -30,8 +30,10 @@ Parameters
 ----------
 type : str
     Type of the decay mode, e.g., 'beta-'
-target : str
-    Nuclide resulting from decay
+target : str or None
+    Nuclide resulting from decay. A value of ``None`` implies the
+    target does not exist in the currently configured depletion
+    chain
 branching_ratio : float
     Branching ratio of the decay mode
 
@@ -53,8 +55,11 @@ Parameters
 ----------
 type : str
     Type of the reaction, e.g., 'fission'
-target : str
-    nuclide resulting from reaction
+target : str or None
+    Nuclide resulting from reaction. A value of ``None``
+    implies either no single target, e.g. from fission,
+    or that the target nuclide is not considered
+    in the current depletion chain
 Q : float
     Q value of the reaction in [eV]
 branching_ratio : float
@@ -179,6 +184,8 @@ class Nuclide(object):
         for decay_elem in element.iter('decay'):
             d_type = decay_elem.get('type')
             target = decay_elem.get('target')
+            if target is not None and target.lower() == "nothing":
+                target = None
             branching_ratio = float(decay_elem.get('branching_ratio'))
             nuc.decay_modes.append(DecayTuple(d_type, target, branching_ratio))
 
@@ -192,6 +199,8 @@ class Nuclide(object):
             # just set null values
             if r_type != 'fission':
                 target = reaction_elem.get('target')
+                if target is not None and target.lower() == "nothing":
+                    target = None
             else:
                 target = None
                 if fission_q is not None:
@@ -226,7 +235,7 @@ class Nuclide(object):
             for mode, daughter, br in self.decay_modes:
                 mode_elem = ET.SubElement(elem, 'decay')
                 mode_elem.set('type', mode)
-                mode_elem.set('target', daughter)
+                mode_elem.set('target', daughter or "Nothing")
                 mode_elem.set('branching_ratio', str(br))
 
         elem.set('reactions', str(len(self.reactions)))
@@ -234,7 +243,7 @@ class Nuclide(object):
             rx_elem = ET.SubElement(elem, 'reaction')
             rx_elem.set('type', rx)
             rx_elem.set('Q', str(Q))
-            if rx != 'fission':
+            if rx != 'fission' or daughter is not None:
                 rx_elem.set('target', daughter)
             if br != 1.0:
                 rx_elem.set('branching_ratio', str(br))
@@ -461,6 +470,7 @@ class FissionYieldDistribution(Mapping):
             data_elem.text = " ".join(map(str, yield_obj.yields))
 
 
+
 class FissionYield(Mapping):
     """Mapping for fission yields of a parent at a specific energy
 

From 58b644cd3b7dfc2872895b73f7795a9e9d96561b Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Wed, 18 Mar 2020 20:08:03 -0400
Subject: [PATCH 02/10] Minor fission yield touch ups

---
 openmc/deplete/nuclide.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py
index 7235745aa2..960ba6eff5 100644
--- a/openmc/deplete/nuclide.py
+++ b/openmc/deplete/nuclide.py
@@ -406,9 +406,7 @@ class FissionYieldDistribution(Mapping):
         for g_index, energy in enumerate(energies):
             prod_map = fission_yields[energy]
             for prod_ix, product in enumerate(ordered_prod):
-                yield_val = prod_map.get(product)
-                yield_matrix[g_index, prod_ix] = (
-                    0.0 if yield_val is None else yield_val)
+                yield_matrix[g_index, prod_ix] = prod_map.get(product, 0.0)
         self.energies = tuple(energies)
         self.products = tuple(ordered_prod)
         self.yield_matrix = yield_matrix
@@ -444,7 +442,7 @@ class FissionYieldDistribution(Mapping):
         FissionYieldDistribution
         """
         all_yields = {}
-        for elem_index, yield_elem in enumerate(element.iter("fission_yields")):
+        for yield_elem in element.iter("fission_yields"):
             energy = float(yield_elem.get("energy"))
             products = yield_elem.find("products").text.split()
             yields = map(float, yield_elem.find("data").text.split())

From 2c2d9d7f312e8b7dfcea10de2ca9aee8e8361f3c Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Wed, 18 Mar 2020 20:09:19 -0400
Subject: [PATCH 03/10] Provide FissionYieldDistribution.restrict_products

---
 openmc/deplete/nuclide.py                | 32 +++++++++++++++++++++++-
 tests/unit_tests/test_deplete_nuclide.py | 15 +++++++++++
 2 files changed, 46 insertions(+), 1 deletion(-)

diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py
index 960ba6eff5..7ba72b6d56 100644
--- a/openmc/deplete/nuclide.py
+++ b/openmc/deplete/nuclide.py
@@ -13,7 +13,7 @@ try:
 except ImportError:
     import xml.etree.ElementTree as ET
 
-from numpy import empty
+from numpy import empty, searchsorted
 
 from openmc.checkvalue import check_type
 
@@ -467,6 +467,36 @@ class FissionYieldDistribution(Mapping):
             data_elem = ET.SubElement(yield_element, "data")
             data_elem.text = " ".join(map(str, yield_obj.yields))
 
+    def restrict_products(self, possible_products):
+        """Return a new distribution with select products
+
+        Parameters
+        ----------
+        possible_products : iterable of str
+            Candidate pool of fission products. Existing products
+            not contained here will not exist in the new instance
+
+        Returns
+        -------
+        FissionYieldDistribution or None
+            A value of None indicates no values in
+            ``possible_products`` exist in :attr:`products`
+
+        """
+
+        overlap = set(self.products).intersection(possible_products)
+        if not overlap:
+            return None
+
+        products = sorted(overlap)
+        indices = searchsorted(self.products, products)
+
+        # coerce back to dictionary to pass back to __init__
+        new_yields = {}
+        for ene, yields in zip(self.energies, self.yield_matrix.copy()):
+            new_yields[ene] = dict(zip(products, yields[indices]))
+
+        return type(self)(new_yields)
 
 
 class FissionYield(Mapping):
diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py
index b8a258df00..92b004845a 100644
--- a/tests/unit_tests/test_deplete_nuclide.py
+++ b/tests/unit_tests/test_deplete_nuclide.py
@@ -176,6 +176,21 @@ def test_fission_yield_distribution():
     with pytest.raises(TypeError):
         orig_yields *= similar
 
+    # Test restriction of fission products
+    strict_restrict = yield_dist.restrict_products(["Xe135", "Sm149"])
+    with_extras = yield_dist.restrict_products(
+        ["Xe135", "Sm149", "H1", "U235"])
+
+    assert strict_restrict.products == ("Sm149", "Xe135", )
+    assert strict_restrict.energies == yield_dist.energies
+    assert with_extras.products == ("Sm149", "Xe135", )
+    assert with_extras.energies == yield_dist.energies
+    for ene, new_yields in strict_restrict.items():
+        for product in strict_restrict.products:
+            assert new_yields[product] == yield_dist[ene][product]
+            assert with_extras[ene][product] == yield_dist[ene][product]
+
+    assert yield_dist.restrict_products(["U235"]) is None
 
 def test_validate():
 

From 789926d107ce447886034a4e6f01bed0587a8091 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Wed, 18 Mar 2020 20:10:49 -0400
Subject: [PATCH 04/10] Provide Chain.reduce for following paths

---
 openmc/deplete/chain.py | 142 +++++++++++++++++++++++++++++++++++++++-
 1 file changed, 141 insertions(+), 1 deletion(-)

diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py
index 8c31914acb..d820fe520b 100644
--- a/openmc/deplete/chain.py
+++ b/openmc/deplete/chain.py
@@ -10,7 +10,7 @@ import math
 import re
 from collections import OrderedDict, defaultdict
 from collections.abc import Mapping, Iterable
-from numbers import Real
+from numbers import Real, Integral
 from warnings import warn
 
 from openmc.checkvalue import check_type, check_greater_than
@@ -820,3 +820,143 @@ class Chain(object):
                 return stat
             valid = valid and stat
         return valid
+
+    def reduce(self, initial_isotopes, level=None):
+        """Reduce the size of the chain by follwing transmutation paths
+
+        Parameters
+        ----------
+        initial_isotopes : iterable of str
+            Start the search based on the contents of these isotopes
+        level : int, optional
+            Depth of transmuation path to follow. Must be greater than
+            or equal to zero. A value of zero returns a chain with
+            ``initial_isotopes``. The default value of None implies
+            that all isotopes that appear in the transmutation paths
+            of the initial isotopes and their progeny should be
+            explored
+
+        Returns
+        -------
+        Chain
+
+        """
+        check_type("initial_isotopes", initial_isotopes, Iterable, str)
+        if level is None:
+            level = math.inf
+        else:
+            check_type("level", level, Integral)
+            check_greater_than("level", level, 0, equality=True)
+
+        all_isotopes = self._follow(set(initial_isotopes), level)
+
+        # Avoid re-sorting for fission yields
+        name_sort = sorted(all_isotopes)
+
+        nuclides = []
+        nuclide_dict = {}
+        reactions = set()
+
+        for idx, iso in enumerate(sorted(all_isotopes, key=openmc.data.zam)):
+            previous = self[iso]
+            new_nuclide = Nuclide(previous.name)
+            new_nuclide.half_life = previous.half_life
+            new_nuclide.decay_energy = new_nuclide.decay_energy
+
+            new_decay = []
+            for mode in previous.decay_modes:
+                if mode.target in all_isotopes:
+                    new_decay.append(mode)
+                else:
+                    new_decay.append(DecayTuple(
+                        mode.type, None, mode.branching_ratio))
+            new_nuclide.decay_modes = new_decay
+
+            new_reactions = []
+            for rxn in previous.reactions:
+                if rxn.target in all_isotopes:
+                    new_reactions.append(rxn)
+                    reactions.add(rxn.type)
+                elif rxn.type == "fission":
+                    new_yields = new_nuclide.yield_data = (
+                        previous.yield_data.restrict_products(name_sort))
+                    if new_yields is not None:
+                        new_reactions.append(rxn)
+                        reactions.add("fission")
+                # Maintain total destruction rates but set no target
+                else:
+                    new_reactions.append(ReactionTuple(
+                        rxn.type, None, rxn.Q, rxn.branching_ratio))
+                    reactions.add(rxn.type)
+
+            new_nuclide.reactions = new_reactions
+
+            nuclides.append(new_nuclide)
+            nuclide_dict[iso] = idx
+
+        new_chain = type(self)()
+        new_chain.nuclides = nuclides
+        new_chain.nuclide_dict = nuclide_dict
+
+        # Doesn't appear that the ordering matters for the reactions,
+        # just the contents
+        new_chain.reactions = sorted(reactions)
+
+        return new_chain
+
+    def _follow(self, isotopes, level):
+        """Return all isotopes present up to depth level"""
+        found = set(isotopes)
+        remaining = set(self.nuclide_dict)
+        if not found.issubset(remaining):
+            raise IndexError(
+                "The following isotopes were not found in the chain: "
+                "{}".format(", ".join(found - remaining)))
+
+        if level == 0:
+            return found
+
+        remaining.difference_update(found)
+
+        depth = 0
+        next_iso = set()
+
+        while depth < level and remaining:
+            # Exhaust all isotopes at this level
+            while isotopes:
+                iso = isotopes.pop()
+                found.add(iso)
+                nuclide = self[iso]
+
+                # Follow all transmutation paths for this nuclide
+                for rxn in nuclide.reactions + nuclide.decay_modes:
+                    if rxn.type == "fission" or rxn.target is None:
+                        continue
+                    # Skip if we've already come across this isotope
+                    elif (rxn.target in next_iso
+                          or rxn.target in found or rxn.target in isotopes):
+                        continue
+                    next_iso.add(rxn.target)
+
+                if nuclide.yield_data is not None:
+                    for product in nuclide.yield_data.products:
+                        if (product in next_iso
+                                or product in found or product in isotopes):
+                            continue
+                        next_iso.add(product)
+
+            if not next_iso:
+                # No additional isotopes to process, nor to update the
+                # current set of discovered isotopes
+                return found
+
+            # Prepare for next dig
+            depth += 1
+            isotopes.update(next_iso)
+            remaining.difference_update(next_iso)
+            next_iso.clear()
+
+        # Process isotope that would have started next depth
+        found.update(isotopes)
+
+        return found

From b01a3e011be15f096f265fec7c19620861d08985 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Wed, 18 Mar 2020 20:11:47 -0400
Subject: [PATCH 05/10] Add unit tests for Chain.reduce

---
 tests/unit_tests/test_deplete_chain.py | 96 +++++++++++++++++++++++++-
 1 file changed, 94 insertions(+), 2 deletions(-)

diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py
index 08c60d3ed1..b7a1e6c0f2 100644
--- a/tests/unit_tests/test_deplete_chain.py
+++ b/tests/unit_tests/test_deplete_chain.py
@@ -6,7 +6,6 @@ from pathlib import Path
 from itertools import product
 
 import numpy as np
-from openmc.data import zam, ATOMIC_SYMBOL
 from openmc.deplete import comm, Chain, reaction_rates, nuclide, cram
 import pytest
 
@@ -405,7 +404,8 @@ def test_fission_yield_attribute(simple_chain):
     dummy_conc = [[1, 2]] * (len(empty_chain.fission_yields) + 1)
     with pytest.raises(
             ValueError, match="fission yield.*not equal.*compositions"):
-       cram.deplete(empty_chain, dummy_conc, None, 0.5)
+        cram.deplete(empty_chain, dummy_conc, None, 0.5)
+
 
 def test_validate(simple_chain):
     """Test the validate method"""
@@ -457,3 +457,95 @@ def test_validate_inputs():
 
     with pytest.raises(ValueError, match="tolerance"):
         c.validate(tolerance=-1)
+
+
+@pytest.fixture
+def gnd_simple_chain():
+    chainfile = Path(__file__).parents[1] / "chain_simple.xml"
+    return Chain.from_xml(chainfile)
+
+
+def test_reduce(gnd_simple_chain):
+    ref_U5 = gnd_simple_chain["U235"]
+    ref_iodine = gnd_simple_chain["I135"]
+    ref_U5_yields = ref_U5.yield_data
+
+    no_depth = gnd_simple_chain.reduce(["U235", "I135"], 0)
+    # We should get a chain just containing U235 and I135
+    assert len(no_depth) == 2
+    assert set(no_depth.reactions) == set(gnd_simple_chain.reactions)
+
+    u5_round0 = no_depth["U235"]
+    assert u5_round0.n_decay_modes == ref_U5.n_decay_modes
+    for newmode, refmode in zip(u5_round0.decay_modes, ref_U5.decay_modes):
+        assert newmode.target is None
+        assert newmode.type == refmode.type, newmode
+        assert newmode.branching_ratio == refmode.branching_ratio, newmode
+
+    assert u5_round0.n_reaction_paths == ref_U5.n_reaction_paths
+    for newrxn, refrxn in zip(u5_round0.reactions, ref_U5.reactions):
+        assert newrxn.target is None, newrxn
+        assert newrxn.type == refrxn.type, newrxn
+        assert newrxn.Q == refrxn.Q, newrxn
+        assert newrxn.branching_ratio == refrxn.branching_ratio, newrxn
+
+    assert u5_round0.yield_data is not None
+    assert u5_round0.yield_data.products == ("I135", )
+    assert u5_round0.yield_data.yield_matrix == (
+            ref_U5_yields.yield_matrix[:, ref_U5_yields.products.index("I135")]
+    )
+
+    bareI5 = no_depth["I135"]
+    assert bareI5.n_decay_modes == ref_iodine.n_decay_modes
+    for newmode, refmode in zip(bareI5.decay_modes, ref_iodine.decay_modes):
+        assert newmode.target is None
+        assert newmode.type == refmode.type, newmode
+        assert newmode.branching_ratio == refmode.branching_ratio, newmode
+
+    assert bareI5.n_reaction_paths == ref_iodine.n_reaction_paths
+    for newrxn, refrxn in zip(bareI5.reactions, ref_iodine.reactions):
+        assert newrxn.target is None, newrxn
+        assert newrxn.type == refrxn.type, newrxn
+        assert newrxn.Q == refrxn.Q, newrxn
+        assert newrxn.branching_ratio == refrxn.branching_ratio, newrxn
+
+    follow_u5 = gnd_simple_chain.reduce(["U235"], 1)
+    u5_round1 = follow_u5["U235"]
+    assert u5_round1.decay_modes == ref_U5.decay_modes
+    assert u5_round1.reactions == ref_U5.reactions
+    assert u5_round1.yield_data is not None
+    assert (
+        u5_round1.yield_data.yield_matrix == ref_U5_yields.yield_matrix
+    ).all()
+
+    # Per the chain_simple.xml
+    # I135 -> Xe135 -> Cs135
+    # I135 -> Xe136
+    # No limit on depth
+    iodine_chain = gnd_simple_chain.reduce(["I135"])
+    truncated_iodine = gnd_simple_chain.reduce(["I135"], 1)
+    assert len(iodine_chain) == 4
+    assert len(truncated_iodine) == 3
+    assert set(iodine_chain.nuclide_dict) == {
+        "I135", "Xe135", "Xe136", "Cs135"}
+    assert set(truncated_iodine.nuclide_dict) == {"I135", "Xe135", "Xe136"}
+    assert iodine_chain.reactions == ["(n,gamma)"]
+    assert iodine_chain["I135"].decay_modes == ref_iodine.decay_modes
+    assert iodine_chain["I135"].reactions == ref_iodine.reactions
+    for mode in truncated_iodine["Xe135"].decay_modes:
+        assert mode.target is None
+
+    # Test that no FissionYieldDistribution is made if there are no
+    # fission products
+    u5_noyields = gnd_simple_chain.reduce(["U235"], 0)["U235"]
+    assert u5_noyields.yield_data is None
+
+    # Check early termination if the eventual full chain
+    # is specified by using the iodine isotopes
+    new_iodine = gnd_simple_chain.reduce(set(iodine_chain.nuclide_dict))
+    assert set(iodine_chain.nuclide_dict) == set(new_iodine.nuclide_dict)
+
+    # Failure if some requested isotopes not in chain
+
+    with pytest.raises(IndexError, match=".*not found.*Xx999"):
+        gnd_simple_chain.reduce(["U235", "Xx999"])

From d67d0cd0f9deb79d5627fa79bc985854d89f5d03 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Fri, 20 Mar 2020 18:11:54 -0400
Subject: [PATCH 06/10] Allow Operator to reduce depletion chain from geometry

A new parameter, reduce_chain, can be used to instruct the
operator to create a smaller depletion chain. The geometry is
used to find all nuclides in burnable materials. The argument
can be a boolean, or a non-negative integer. A value of True
indicates no depth for following the paths, while an integer
can be used to directly specify the depth.
---
 openmc/deplete/operator.py | 20 +++++++++++++++++++-
 1 file changed, 19 insertions(+), 1 deletion(-)

diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index 8eb033783f..2572150012 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -113,6 +113,12 @@ class Operator(TransportOperator):
         ``fission_yield_mode``. Will be passed directly on to the
         helper. Passing a value of None will use the defaults for
         the associated helper.
+    reduce_chain : bool or int, optional
+        If ``True`` or an integer, create a reduced depletion chain by
+        following transmuation paths for isotopes in burnable materials.
+        A value of ``True`` implies to follow all paths to completion,
+        while a non-negative integer indicates the depth. See
+        :meth:`openmc.deplete.Chain.reduce`
 
     Attributes
     ----------
@@ -158,7 +164,8 @@ class Operator(TransportOperator):
     def __init__(self, geometry, settings, chain_file=None, prev_results=None,
                  diff_burnable_mats=False, energy_mode="fission-q",
                  fission_q=None, dilute_initial=1.0e3,
-                 fission_yield_mode="constant", fission_yield_opts=None):
+                 fission_yield_mode="constant", fission_yield_opts=None,
+                 reduce_chain=False):
         if fission_yield_mode not in self._fission_helpers:
             raise KeyError(
                 "fission_yield_mode must be one of {}, not {}".format(
@@ -179,6 +186,17 @@ class Operator(TransportOperator):
         self.geometry = geometry
         self.diff_burnable_mats = diff_burnable_mats
 
+        # Reduce the chain before we create more materials
+        if reduce_chain is not False:
+            all_isotopes = set()
+            for material in geometry.get_all_materials().values():
+                if not material.depletable:
+                    continue
+                for name, _dens_percent, _dens_type in material.nuclides:
+                    all_isotopes.add(name)
+            level = None if reduce_chain is True else reduce_chain
+            self.chain = self.chain.reduce(all_isotopes, level)
+
         # Differentiate burnable materials with multiple instances
         if self.diff_burnable_mats:
             self._differentiate_burnable_mats()

From 3669709061f6dc51591bdb46ee6dbf6170b500d7 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Wed, 25 Mar 2020 06:29:21 -0400
Subject: [PATCH 07/10] Apply suggestions from code review

* Use set operations (|=) rather than in-place methods (update)
* Remove second arguments to some asserts in Chain.reduce testing

Co-Authored-By: Paul Romano <paul.k.romano@gmail.com>
---
 openmc/deplete/chain.py                  |  8 ++++----
 tests/unit_tests/test_deplete_chain.py   | 26 ++++++++++++------------
 tests/unit_tests/test_deplete_nuclide.py |  4 ++--
 3 files changed, 19 insertions(+), 19 deletions(-)

diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py
index d820fe520b..fce5051b0b 100644
--- a/openmc/deplete/chain.py
+++ b/openmc/deplete/chain.py
@@ -906,7 +906,7 @@ class Chain(object):
 
     def _follow(self, isotopes, level):
         """Return all isotopes present up to depth level"""
-        found = set(isotopes)
+        found = isotopes.copy()
         remaining = set(self.nuclide_dict)
         if not found.issubset(remaining):
             raise IndexError(
@@ -916,7 +916,7 @@ class Chain(object):
         if level == 0:
             return found
 
-        remaining.difference_update(found)
+        remaining -= found
 
         depth = 0
         next_iso = set()
@@ -952,8 +952,8 @@ class Chain(object):
 
             # Prepare for next dig
             depth += 1
-            isotopes.update(next_iso)
-            remaining.difference_update(next_iso)
+            isotopes |= next_iso
+            remaining -= next_iso
             next_iso.clear()
 
         # Process isotope that would have started next depth
diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py
index b7a1e6c0f2..ba3be4e37a 100644
--- a/tests/unit_tests/test_deplete_chain.py
+++ b/tests/unit_tests/test_deplete_chain.py
@@ -479,18 +479,18 @@ def test_reduce(gnd_simple_chain):
     assert u5_round0.n_decay_modes == ref_U5.n_decay_modes
     for newmode, refmode in zip(u5_round0.decay_modes, ref_U5.decay_modes):
         assert newmode.target is None
-        assert newmode.type == refmode.type, newmode
-        assert newmode.branching_ratio == refmode.branching_ratio, newmode
+        assert newmode.type == refmode.type
+        assert newmode.branching_ratio == refmode.branching_ratio
 
     assert u5_round0.n_reaction_paths == ref_U5.n_reaction_paths
     for newrxn, refrxn in zip(u5_round0.reactions, ref_U5.reactions):
-        assert newrxn.target is None, newrxn
-        assert newrxn.type == refrxn.type, newrxn
-        assert newrxn.Q == refrxn.Q, newrxn
-        assert newrxn.branching_ratio == refrxn.branching_ratio, newrxn
+        assert newrxn.target is None
+        assert newrxn.type == refrxn.type
+        assert newrxn.Q == refrxn.Q
+        assert newrxn.branching_ratio == refrxn.branching_ratio
 
     assert u5_round0.yield_data is not None
-    assert u5_round0.yield_data.products == ("I135", )
+    assert u5_round0.yield_data.products == ("I135",)
     assert u5_round0.yield_data.yield_matrix == (
             ref_U5_yields.yield_matrix[:, ref_U5_yields.products.index("I135")]
     )
@@ -499,15 +499,15 @@ def test_reduce(gnd_simple_chain):
     assert bareI5.n_decay_modes == ref_iodine.n_decay_modes
     for newmode, refmode in zip(bareI5.decay_modes, ref_iodine.decay_modes):
         assert newmode.target is None
-        assert newmode.type == refmode.type, newmode
-        assert newmode.branching_ratio == refmode.branching_ratio, newmode
+        assert newmode.type == refmode.type
+        assert newmode.branching_ratio == refmode.branching_ratio
 
     assert bareI5.n_reaction_paths == ref_iodine.n_reaction_paths
     for newrxn, refrxn in zip(bareI5.reactions, ref_iodine.reactions):
-        assert newrxn.target is None, newrxn
-        assert newrxn.type == refrxn.type, newrxn
-        assert newrxn.Q == refrxn.Q, newrxn
-        assert newrxn.branching_ratio == refrxn.branching_ratio, newrxn
+        assert newrxn.target is None
+        assert newrxn.type == refrxn.type
+        assert newrxn.Q == refrxn.Q
+        assert newrxn.branching_ratio == refrxn.branching_ratio
 
     follow_u5 = gnd_simple_chain.reduce(["U235"], 1)
     u5_round1 = follow_u5["U235"]
diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py
index 92b004845a..126206cc70 100644
--- a/tests/unit_tests/test_deplete_nuclide.py
+++ b/tests/unit_tests/test_deplete_nuclide.py
@@ -181,9 +181,9 @@ def test_fission_yield_distribution():
     with_extras = yield_dist.restrict_products(
         ["Xe135", "Sm149", "H1", "U235"])
 
-    assert strict_restrict.products == ("Sm149", "Xe135", )
+    assert strict_restrict.products == ("Sm149", "Xe135")
     assert strict_restrict.energies == yield_dist.energies
-    assert with_extras.products == ("Sm149", "Xe135", )
+    assert with_extras.products == ("Sm149", "Xe135")
     assert with_extras.energies == yield_dist.energies
     for ene, new_yields in strict_restrict.items():
         for product in strict_restrict.products:

From c2562322ab7a788bdc09ed1846a2934582756060 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Wed, 25 Mar 2020 06:48:43 -0400
Subject: [PATCH 08/10] Expand Chain.reduce docstring: methodology and return

The methodology is explained through a simple example mostly
borrowed from the chain_simple.xml test chain. A few different
cases are presented, showing the effect of initial isotopes and
search depth.

A description of the returned chain is also included.
---
 openmc/deplete/chain.py | 27 ++++++++++++++++++++++++++-
 1 file changed, 26 insertions(+), 1 deletion(-)

diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py
index fce5051b0b..a23ba2bc1c 100644
--- a/openmc/deplete/chain.py
+++ b/openmc/deplete/chain.py
@@ -822,7 +822,30 @@ class Chain(object):
         return valid
 
     def reduce(self, initial_isotopes, level=None):
-        """Reduce the size of the chain by follwing transmutation paths
+        """Reduce the size of the chain by following transmutation paths
+
+        As an example, consider a simple chain with the following
+        isotopes and transmutation paths::
+
+            U235 (n,gamma) U236
+                 (n,fission) (Xe135, I135, Cs135)
+            I135 (beta decay) Xe135 (beta decay) Cs135
+            Xe135 (n,gamma) Xe136
+
+        Calling ``chain.reduce(["I135"])`` will produce a depletion
+        chain that contains only isotopes that would originate from
+        I135: I135, Xe135, Cs135, and Xe136. U235 and U236 will not
+        be included, but multiple isotopes can be used to start
+        the search.
+
+        The ``level`` value controls the depth of the search.
+        ``chain.reduce(["U235"], level=1)`` would return a chain
+        with all isotopes except Xe136, since it is two transmutations
+        removed from U235 in this case.
+
+        While targets will not be included in the new chain, the
+        total destruction rate and decay rate of included isotopes
+        will be preserved.
 
         Parameters
         ----------
@@ -839,6 +862,8 @@ class Chain(object):
         Returns
         -------
         Chain
+            Depletion chain containing isotopes that would appear
+            after following up to ``level`` reactions and decay paths
 
         """
         check_type("initial_isotopes", initial_isotopes, Iterable, str)

From d3368db5ee030bdedf5e6b039ddda2c126aaeafe Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Wed, 25 Mar 2020 07:10:32 -0400
Subject: [PATCH 09/10] Use reduce_chain and reduce_chain_level with Operator

Control if the chain is to be reduced, and how far down the paths
to follow. Defaults to not reducing the depletion chain. If the
chain is reduced, the default depth is None, implying no limit
---
 openmc/deplete/operator.py | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index 2572150012..5fa6745c8e 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -113,12 +113,13 @@ class Operator(TransportOperator):
         ``fission_yield_mode``. Will be passed directly on to the
         helper. Passing a value of None will use the defaults for
         the associated helper.
-    reduce_chain : bool or int, optional
-        If ``True`` or an integer, create a reduced depletion chain by
-        following transmuation paths for isotopes in burnable materials.
-        A value of ``True`` implies to follow all paths to completion,
-        while a non-negative integer indicates the depth. See
-        :meth:`openmc.deplete.Chain.reduce`
+    reduce_chain : bool, optional
+        If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the
+        depletion chain up to ``reduce_chain_level``. Default is False.
+    reduce_chain_level : int, optional
+        Depth of the search when reducing the depletion chain. Only used
+        if ``reduce_chain`` evaluates to true. The default value of
+        ``None`` implies no limit on the depth.
 
     Attributes
     ----------
@@ -165,7 +166,7 @@ class Operator(TransportOperator):
                  diff_burnable_mats=False, energy_mode="fission-q",
                  fission_q=None, dilute_initial=1.0e3,
                  fission_yield_mode="constant", fission_yield_opts=None,
-                 reduce_chain=False):
+                 reduce_chain=False, reduce_chain_level=None):
         if fission_yield_mode not in self._fission_helpers:
             raise KeyError(
                 "fission_yield_mode must be one of {}, not {}".format(
@@ -187,15 +188,14 @@ class Operator(TransportOperator):
         self.diff_burnable_mats = diff_burnable_mats
 
         # Reduce the chain before we create more materials
-        if reduce_chain is not False:
+        if reduce_chain:
             all_isotopes = set()
             for material in geometry.get_all_materials().values():
                 if not material.depletable:
                     continue
                 for name, _dens_percent, _dens_type in material.nuclides:
                     all_isotopes.add(name)
-            level = None if reduce_chain is True else reduce_chain
-            self.chain = self.chain.reduce(all_isotopes, level)
+            self.chain = self.chain.reduce(all_isotopes, reduce_chain_level)
 
         # Differentiate burnable materials with multiple instances
         if self.diff_burnable_mats:

From f882cf666bb02a2650342bda9842fcf804f91cbc Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Wed, 25 Mar 2020 07:12:40 -0400
Subject: [PATCH 10/10] Remove unused imports from operator.py

time, itertools.chain, and h5py not used
---
 openmc/deplete/operator.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index 5fa6745c8e..ce7e709d9c 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -10,13 +10,10 @@ densities is all done in-memory instead of through the filesystem.
 import sys
 import copy
 from collections import OrderedDict
-from itertools import chain
 import os
-import time
 import xml.etree.ElementTree as ET
 from warnings import warn
 
-import h5py
 import numpy as np
 from uncertainties import ufloat