From e37ff271c6326fde875619d8ff7cc3bc01c67995 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 20 Aug 2019 16:56:15 -0500 Subject: [PATCH 1/4] Generalize Chain.set_capture_branches -> set_branch_ratios Generalized method that allows the modification of branching ratios for an arbitrary reaction. A similar replacement has been made with the companion get_capture_branches -> get_branch_ratios. Both of the new methods accept an additional reaction argument, which defaults to "(n,gamma)", indicating what reaction to modify/retrieve. Also added a tolerance argument, defaults to 1e-5, for checking the sum of user-supplied branching ratios. This tolerance is used as 1 - tolerance < sum_br < 1 + tolerance. Additional logic is included if the user has not specified the ground state, which can be inferred. For this case, the lower end is allowed to be less than 1 - tolerance, as the ground state will have a branch ratio such the sum of ratios is unity. Products are not checked for consistency, that is left up to the user. tests/unit_test/test_deplete_chain.py has been modified with the new notation, and also with a variety of calls with and without the optional reaction argument. A minor test is added that modifies (n,alpha) and (n,2n) reactions just to ensure that arbitrary reactions can be used. Minor changes: * Internal variables changed to reflect the generality of the method * Added documentation on what exceptions are raised in set_branch_ratios Related: openmc issue #1237 --- openmc/deplete/chain.py | 138 +++++++++++++++++-------- tests/unit_tests/test_deplete_chain.py | 70 ++++++++++--- 2 files changed, 146 insertions(+), 62 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 1f16d9caf..3b1870b97 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -454,70 +454,96 @@ class Chain(object): dict.update(matrix_dok, matrix) return matrix_dok.tocsr() - def get_capture_branches(self): - """Return a dictionary with capture branching ratios + def get_branch_ratios(self, reaction="(n,gamma)"): + """Return a dictionary with reaction branching ratios + + Parameters + ---------- + reaction : str, optional + Reaction name like ``"(n,gamma)"`` [default], or + ``"(n,alpha)"``. Returns ------- - capt : - nested dict of parent nuclide keys with capture targets and - branching ratios:: + branches : dict + nested dict of parent nuclide keys with reaction targets and + branching ratios. Consider the capture, ``"(n,gamma)"``, + reaction for Am241:: {"Am241": {"Am242": 0.91, "Am242_m1": 0.09}} See Also -------- - :meth:`set_capture_branches` - + :meth:`set_branch_ratios` """ capt = {} for nuclide in self.nuclides: nuc_capt = {} for rx in nuclide.reactions: - if rx.type == "(n,gamma)" and rx.branching_ratio != 1.0: + if rx.type == reaction and rx.branching_ratio != 1.0: nuc_capt[rx.target] = rx.branching_ratio if len(nuc_capt) > 0: capt[nuclide.name] = nuc_capt return capt - def set_capture_branches(self, branch_ratios, strict=True): - """Set the capture branching ratios - - To provide a buffer around floating point precisions, - the sum of all branching ratios from a single parent - cannot be greater than 1.00001. + def set_branch_ratios(self, branch_ratios, reaction="(n,gamma)", + strict=True, tolerance=1e-5): + """Set the branching ratios for a given reactions Parameters ---------- branch_ratios : dict of {str: {str: float}} Capture branching ratios to be inserted. First layer keys are names of parent nuclides, e.g. - ``"Am241"``. The capture branching ratios for these + ``"Am241"``. The branching ratios for these parents will be modified. Corresponding values are dictionaries of ``{target: branching_ratio}`` - strict : bool - If this evalutes to ``True``, then all parents and - products must exist in the :class:`Chain`. A - :class:`KeyError` will be raised at the first - nuclide that does not exist. Otherwise, print - a warning message for missing parents and/or - products. + reaction : str, optional + Reaction name like ``"(n,gamma)"`` [default], or + ``"(n, alpha)"``. + strict : bool, optional + Error control. If this evalutes to ``True``, then errors will + be raised if inconsistencies are found. Otherwise, warnings + will be raised for most issues. + tolerance : float, optional + Tolerance on the sum of all branching ratios for a + single parent. Will be checked with:: + + 1 - tol < sum_br < 1 + tol + + Raises + ------ + IndexError + If no isotopes were found on the chain that have the requested + reaction + KeyError + If ``strict`` evaluates to ``False`` and a parent isotope in + ``branch_ratios`` does not exist on the chain + AttributeError + If ``strict`` evaluates to ``False`` and a parent isotope in + ``branch_ratios`` does not have the requested reaction + ValueError + If ``strict`` evalutes to ``False`` and the sum of one parents + branch ratios is outside 1 +/- ``tolerance`` See Also -------- - :meth:`get_capture_branches` + :meth:`get_branch_ratios` """ # Store some useful information through the validation stage sums = {} - capt_ix_map = {} + rxn_ix_map = {} grounds = {} + tolerance = abs(tolerance) + missing_parents = set() missing_products = {} - no_capture = set() + missing_reaction = set() + bad_sums = {} # Check for validity before manipulation @@ -543,11 +569,11 @@ class Chain(object): if prod_flag: continue - # Make sure this nuclide has capture reactions + # Make sure this nuclide has the reaction indexes = [] for ix, rx in enumerate(self[parent].reactions): - if rx.type == "(n,gamma)": + if rx.type == reaction: indexes.append(ix) if "_m" not in rx.target: grounds[parent] = rx.target @@ -555,24 +581,39 @@ class Chain(object): if len(indexes) == 0: if strict: raise AttributeError( - "Nuclide {} does not have capture reactions in " - "this {}".format(parent, self.__class__.__name__)) - no_capture.add(parent) + "Nuclide {} does not have {} reactions".format( + parent, reaction)) + missing_reaction.add(parent) continue - capt_ix_map[parent] = indexes - this_sum = sum(sub.values()) - check_less_than(parent + " ratios", this_sum, 1.00001) - sums[parent] = this_sum + # sum of branching ratios can be lower than 1 if no ground + # target is given, but never greater + if (this_sum >= 1 + tolerance or (grounds[parent] in sub + and this_sum <= 1 - tolerance)): + if strict: + msg = ("Sum of {} branching ratios for {} " + "({:7.3f}) outside tolerance of 1 +/- " + "{:5.3e}".format( + reaction, parent, this_sum, tolerance)) + raise ValueError(msg) + bad_sums[parent] = this_sum + else: + rxn_ix_map[parent] = indexes + sums[parent] = this_sum + + if len(rxn_ix_map) == 0: + raise IndexError( + "No {} reactions found in this {}".format( + reaction, self.__class__.__name__)) if len(missing_parents) > 0: warn("The following nuclides were not found in {}: {}".format( self.__class__.__name__, ", ".join(sorted(missing_parents)))) - if len(no_capture) > 0: - warn("The following nuclides did not have capture reactions: " - "{}".format(", ".join(sorted(no_capture)))) + if len(missing_reaction) > 0: + warn("The following nuclides did not have {} reactions: " + "{}".format(reaction, ", ".join(sorted(missing_reaction)))) if len(missing_products) > 0: tail = ("{} -> {}".format(k, v) @@ -581,28 +622,35 @@ class Chain(object): "parents were unmodified: \n{}".format( self.__class__.__name__, ", ".join(tail))) + if len(bad_sums) > 0: + tail = ("{}: {:5.3f}".format(k, s) + for k, s in sorted(bad_sums.items())) + warn("The following parent nuclides were given {} branch ratios " + "with a sum outside tolerance of 1 +/- {:5.3e}:\n{}".format( + reaction, tolerance, "\n".join(tail))) + # Insert new ReactionTuples with updated branch ratios - for parent_name, capt_index in capt_ix_map.items(): + for parent_name, rxn_index in rxn_ix_map.items(): parent = self[parent_name] new_ratios = branch_ratios[parent_name] - capt_index = capt_ix_map[parent_name] + rxn_index = rxn_ix_map[parent_name] # Assume Q value is independent of target state - capt_Q = parent.reactions[capt_index[0]].Q + rxn_Q = parent.reactions[rxn_index[0]].Q - # Remove existing capture reactions + # Remove existing reactions - for ix in reversed(capt_index): + for ix in reversed(rxn_index): parent.reactions.pop(ix) all_meta = True for tgt, br in new_ratios.items(): - all_meta = all_meta and ("_m" in tgt) + all_meta = all_meta and ("_m" in tgt) parent.reactions.append(ReactionTuple( - "(n,gamma)", tgt, capt_Q, br)) + reaction, tgt, rxn_Q, br)) if all_meta and sums[parent_name] != 1.0: ground_br = 1.0 - sums[parent_name] @@ -612,4 +660,4 @@ class Chain(object): ground_tgt = gnd_name(pz, pa + 1, 0) new_ratios[ground_tgt] = ground_br parent.reactions.append(ReactionTuple( - "(n,gamma)", ground_tgt, capt_Q, ground_br)) + reaction, ground_tgt, rxn_Q, ground_br)) diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index 33b13b596..686c6d508 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -248,29 +248,29 @@ def test_set_fiss_q(): def test_get_set_chain_br(simple_chain): """Test minor modifications to capture branch ratios""" expected = {"C": {"A": 0.7, "B": 0.3}} - assert simple_chain.get_capture_branches() == expected + assert simple_chain.get_branch_ratios() == expected # safely modify new_chain = Chain.from_xml("chain_test.xml") new_br = {"C": {"A": 0.5, "B": 0.5}, "A": {"C": 0.99, "B": 0.01}} - new_chain.set_capture_branches(new_br) - assert new_chain.get_capture_branches() == new_br + new_chain.set_branch_ratios(new_br) + assert new_chain.get_branch_ratios() == new_br # write, re-read new_chain.export_to_xml("chain_mod.xml") - assert Chain.from_xml("chain_mod.xml").get_capture_branches() == new_br + assert Chain.from_xml("chain_mod.xml").get_branch_ratios() == new_br # Test non-strict [warn, not error] setting bad_br = {"B": {"X": 0.6, "A": 0.4}, "X": {"A": 0.5, "C": 0.5}} bad_br.update(new_br) - new_chain.set_capture_branches(bad_br, strict=False) - assert new_chain.get_capture_branches() == new_br + new_chain.set_branch_ratios(bad_br, strict=False) + assert new_chain.get_branch_ratios() == new_br # Ensure capture reactions are removed rem_br = {"A": {"C": 1.0}} - new_chain.set_capture_branches(rem_br) + new_chain.set_branch_ratios(rem_br) # A is not in returned dict because there is no branch - assert "A" not in new_chain.get_capture_branches() + assert "A" not in new_chain.get_branch_ratios() def test_capture_branch_infer_ground(): @@ -289,9 +289,9 @@ def test_capture_branch_infer_ground(): chain.nuclides.append(xe136m) chain.nuclide_dict[xe136m.name] = len(chain.nuclides) - 1 - chain.set_capture_branches(infer_br) + chain.set_branch_ratios(infer_br, "(n,gamma)") - assert chain.get_capture_branches() == set_br + assert chain.get_branch_ratios("(n,gamma)") == set_br def test_capture_branch_no_rxn(): @@ -307,9 +307,8 @@ def test_capture_branch_no_rxn(): chain.nuclides.append(u5m) chain.nuclide_dict[u5m.name] = len(chain.nuclides) - 1 - phrase = "U234 does not have capture reactions" - with pytest.raises(AttributeError, match=phrase): - chain.set_capture_branches(u4br) + with pytest.raises(AttributeError, match="U234"): + chain.set_branch_ratios(u4br) def test_capture_branch_failures(simple_chain): @@ -318,14 +317,51 @@ def test_capture_branch_failures(simple_chain): # Parent isotope not present br = {"X": {"A": 0.6, "B": 0.7}} with pytest.raises(KeyError, match="X"): - simple_chain.set_capture_branches(br) + simple_chain.set_branch_ratios(br) # Product isotope not present br = {"C": {"X": 0.4, "A": 0.2, "B": 0.4}} with pytest.raises(KeyError, match="X"): - simple_chain.set_capture_branches(br) + simple_chain.set_branch_ratios(br) # Sum of ratios > 1.0 br = {"C": {"A": 1.0, "B": 1.0}} - with pytest.raises(ValueError, match="C ratios"): - simple_chain.set_capture_branches(br) + with pytest.raises(ValueError, match=r"Sum of \(n,gamma\).*for C"): + simple_chain.set_branch_ratios(br, "(n,gamma)") + + +@pytest.mark.parametrize("reaction", ("(n,alpha)", "(n,2n)")) +def test_set_general_branches(reaction): + """Test setting of non-capture branching ratios""" + # Build a mock chain + chain = Chain() + + parent = nuclide.Nuclide() + parent.name = "A" + + ground_tgt = nuclide.Nuclide() + ground_tgt.name = "B" + + meta_tgt = nuclide.Nuclide() + meta_tgt.name = "B_m1" + + for ix, nuc in enumerate((parent, ground_tgt, meta_tgt)): + chain.nuclides.append(nuc) + chain.nuclide_dict[nuc.name] = ix + + # add reactions to parent + parent.reactions.append(nuclide.ReactionTuple( + reaction, ground_tgt.name, 1.0, 0.6)) + parent.reactions.append(nuclide.ReactionTuple( + reaction, meta_tgt.name, 1.0, 0.4)) + + expected_ref = {"A": {"B": 0.6, "B_m1": 0.4}} + + assert chain.get_branch_ratios(reaction) == expected_ref + + # alter and check again + + altered = {"A": {"B": 0.5, "B_m1": 0.5}} + + chain.set_branch_ratios(altered, reaction) + assert chain.get_branch_ratios(reaction) == altered From 34add135be85bba0e4a3c1659f235198875ee18c Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 20 Aug 2019 17:19:53 -0500 Subject: [PATCH 2/4] Leave He4 alone in Chain.set_branch_ratios In an (n, alpha) reaction, the branching ratio for helium should always be one, as it is produced for every reaction. It is the branching ratio of the other targets that must be modified. He4 is treated as a secondary particle for this case, and its reaction is not modified nor removed from the parent. --- openmc/deplete/chain.py | 5 ++++- tests/unit_tests/test_deplete_chain.py | 30 ++++++++++++++++++-------- 2 files changed, 25 insertions(+), 10 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 3b1870b97..b557a6111 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -545,6 +545,9 @@ class Chain(object): missing_reaction = set() bad_sums = {} + # Secondary products, like alpha particles, should not be modified + secondary = "He4" if reaction == "(n,a)" else None + # Check for validity before manipulation for parent, sub in branch_ratios.items(): @@ -573,7 +576,7 @@ class Chain(object): indexes = [] for ix, rx in enumerate(self[parent].reactions): - if rx.type == reaction: + if rx.type == reaction and rx.target != secondary: indexes.append(ix) if "_m" not in rx.target: grounds[parent] = rx.target diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index 686c6d508..3cf726ed1 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -330,38 +330,50 @@ def test_capture_branch_failures(simple_chain): simple_chain.set_branch_ratios(br, "(n,gamma)") -@pytest.mark.parametrize("reaction", ("(n,alpha)", "(n,2n)")) -def test_set_general_branches(reaction): - """Test setting of non-capture branching ratios""" + +def test_set_alpha_branches(): + """Test setting of alpha reaction branching ratios""" # Build a mock chain chain = Chain() parent = nuclide.Nuclide() parent.name = "A" + he4 = nuclide.Nuclide() + he4.name = "He4" + ground_tgt = nuclide.Nuclide() ground_tgt.name = "B" meta_tgt = nuclide.Nuclide() meta_tgt.name = "B_m1" - for ix, nuc in enumerate((parent, ground_tgt, meta_tgt)): + for ix, nuc in enumerate((parent, ground_tgt, meta_tgt, he4)): chain.nuclides.append(nuc) chain.nuclide_dict[nuc.name] = ix # add reactions to parent parent.reactions.append(nuclide.ReactionTuple( - reaction, ground_tgt.name, 1.0, 0.6)) + "(n,a)", ground_tgt.name, 1.0, 0.6)) parent.reactions.append(nuclide.ReactionTuple( - reaction, meta_tgt.name, 1.0, 0.4)) + "(n,a)", meta_tgt.name, 1.0, 0.4)) + parent.reactions.append(nuclide.ReactionTuple( + "(n,a)", he4.name, 1.0, 1.0)) expected_ref = {"A": {"B": 0.6, "B_m1": 0.4}} - assert chain.get_branch_ratios(reaction) == expected_ref + assert chain.get_branch_ratios("(n,a)") == expected_ref # alter and check again altered = {"A": {"B": 0.5, "B_m1": 0.5}} - chain.set_branch_ratios(altered, reaction) - assert chain.get_branch_ratios(reaction) == altered + chain.set_branch_ratios(altered, "(n,a)") + assert chain.get_branch_ratios("(n,a)") == altered + + # make sure that alpha particle still produced + for r in parent.reactions: + if r.target == he4.name: + break + else: + raise ValueError("Helium has been removed and should not have been") From 9287db4a2636bc725050b8f49b427e493c165f0b Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 22 Aug 2019 09:41:39 -0500 Subject: [PATCH 3/4] Support more secondary nuclides in Chain.set_branch_ratios Reactions are pulled from reaction scores table on https://docs.openmc.org/en/latest/usersguide/tallies.html#scores for reactions that produce hydrogen and helium isotopes --- openmc/deplete/chain.py | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index b557a6111..3cb3dbf01 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -12,7 +12,7 @@ from collections import OrderedDict, defaultdict from collections.abc import Mapping from warnings import warn -from openmc.checkvalue import check_type, check_less_than +from openmc.checkvalue import check_type from openmc.data import gnd_name, zam # Try to use lxml if it is available. It preserves the order of attributes and @@ -103,6 +103,19 @@ def replace_missing(product, decay_data): return product +_secondary_particles = { + "p": ["H1"], "d": ["H2"], "t": ["H3"], "3He": ["He3"], "a": ["He4"], + "2nd": ["H2"], "na": ["He4"], "n3a": ["He4"] * 3, "2na": ["He4"], + "3na": ["He4"], "np": ["H1"], "n2a": ["He4"] * 2, + "2n2a": ["He4"] * 2, "nd": ["H2"], "nt": ["H3"], + "nHe-3": ["He3"], "nd2a": ["H2", "He4"], "nt2a": ["H3", "He4", "He4"], + "2np": ["H1"], "3np": ["H1"], "n2p": ["H1"] * 2, + "2a": ["He4"] * 2, "3a": ["He4"] * 3, "2p": ["H1"] * 2, + "pa": ["H1", "He4"], "t2a": ["H3", "He4", "He4"], + "d2a": ["H2", "He4", "He4"], "pd": ["H1", "H2"], "pt": ["H1", "H3"], + "da": ["H2", "He4"]} + + class Chain(object): """Full representation of a depletion chain. @@ -122,7 +135,6 @@ class Chain(object): Reactions that are tracked in the depletion chain nuclide_dict : OrderedDict of str to int Maps a nuclide name to an index in nuclides. - """ def __init__(self): @@ -546,7 +558,8 @@ class Chain(object): bad_sums = {} # Secondary products, like alpha particles, should not be modified - secondary = "He4" if reaction == "(n,a)" else None + secondary = _secondary_particles.get( + reaction[reaction.index(",") + 1:-1], []) # Check for validity before manipulation @@ -576,7 +589,7 @@ class Chain(object): indexes = [] for ix, rx in enumerate(self[parent].reactions): - if rx.type == reaction and rx.target != secondary: + if rx.type == reaction and rx.target not in secondary: indexes.append(ix) if "_m" not in rx.target: grounds[parent] = rx.target From c8e205aff87af51ecc83aafcb8406d091cec9bcb Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 23 Aug 2019 11:43:19 -0500 Subject: [PATCH 4/4] Update dictionary of secondary particles for Chain --- openmc/deplete/chain.py | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 04fdf658e..b3e2dae29 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -104,17 +104,17 @@ def replace_missing(product, decay_data): return product -_secondary_particles = { - "p": ["H1"], "d": ["H2"], "t": ["H3"], "3He": ["He3"], "a": ["He4"], - "2nd": ["H2"], "na": ["He4"], "n3a": ["He4"] * 3, "2na": ["He4"], - "3na": ["He4"], "np": ["H1"], "n2a": ["He4"] * 2, - "2n2a": ["He4"] * 2, "nd": ["H2"], "nt": ["H3"], - "nHe-3": ["He3"], "nd2a": ["H2", "He4"], "nt2a": ["H3", "He4", "He4"], - "2np": ["H1"], "3np": ["H1"], "n2p": ["H1"] * 2, - "2a": ["He4"] * 2, "3a": ["He4"] * 3, "2p": ["H1"] * 2, - "pa": ["H1", "He4"], "t2a": ["H3", "He4", "He4"], - "d2a": ["H2", "He4", "He4"], "pd": ["H1", "H2"], "pt": ["H1", "H3"], - "da": ["H2", "He4"]} +_SECONDARY_PARTICLES = { + "(n,p)": ["H1"], "(n,d)": ["H2"], "(n,t)": ["H3"], "(n,3He)": ["He3"], + "(n,a)": ["He4"], "(n,2nd)": ["H2"], "(n,na)": ["He4"], "(n,3na)": ["He4"], + "(n,n3a)": ["He4"] * 3, "(n,2na)": ["He4"], "(n,np)": ["H1"], + "(n,n2a)": ["He4"] * 2, "(n,2n2a)": ["He4"] * 2, "(n,nd)": ["H2"], + "(n,nt)": ["H3"], "(n,nHe-3)": ["He3"], "(n,nd2a)": ["H2", "He4"], + "(n,nt2a)": ["H3", "He4", "He4"], "(n,2np)": ["H1"], "(n,3np)": ["H1"], + "(n,n2p)": ["H1"] * 2, "(n,2a)": ["He4"] * 2, "(n,3a)": ["He4"] * 3, + "(n,2p)": ["H1"] * 2, "(n,pa)": ["H1", "He4"], + "(n,t2a)": ["H3", "He4", "He4"], "(n,d2a)": ["H2", "He4", "He4"], + "(n,pd)": ["H1", "H2"], "(n,pt)": ["H1", "H3"], "(n,da)": ["H2", "He4"]} class Chain(object): @@ -559,8 +559,7 @@ class Chain(object): bad_sums = {} # Secondary products, like alpha particles, should not be modified - secondary = _secondary_particles.get( - reaction[reaction.index(",") + 1:-1], []) + secondary = _SECONDARY_PARTICLES.get(reaction, []) # Check for validity before manipulation