From e37ff271c6326fde875619d8ff7cc3bc01c67995 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 20 Aug 2019 16:56:15 -0500 Subject: [PATCH] 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