From 220f12c2a74a797b94817ed7050752afd4c68c2f Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 8 Jul 2019 10:56:40 -0500 Subject: [PATCH 1/6] Add set/get capture branching ratios method for Chain Related to #1237. openmc.deplete.Chain.set_capture_branches takes in a dictionary of parent nuclide names to {daughter: ratio}, e.g. {"Am241": {"Am242": 0.9, "Am242_m1": 0.1}} that will be used to overwrite existing capture reactions with new branching ratios. openmc.deplete.Chain.get_capture_branches returns a similar dictionary of capture reactions with multiple targets and their branching ratios. --- openmc/deplete/chain.py | 137 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 136 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 643235d9f5..eec9d3f507 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -11,7 +11,8 @@ import re from collections import OrderedDict, defaultdict from collections.abc import Mapping -from openmc.checkvalue import check_type +from openmc.checkvalue import check_type, check_less_than +from openmc.data import gnd_name, zam # Try to use lxml if it is available. It preserves the order of attributes and # provides a pretty-printer by default. If not available, @@ -451,3 +452,137 @@ class Chain(object): matrix_dok = sp.dok_matrix((n, n)) dict.update(matrix_dok, matrix) return matrix_dok.tocsr() + + def get_capture_branches(self): + """Return a dictionary with capture branching ratios + + Returns + ------- + capt : + nested dict of parent nuclide keys with capture targets and + branching ratios:: + + {"Am241": {"Am242": 0.91, "Am242_m1": 0.09}} + + See Also + -------- + :meth:`set_capture_branches` + + """ + + capt = {} + for nuclide in self.nuclides: + nuc_capt = {} + for rx in nuclide.reactions: + if rx.type == "(n,gamma)" 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): + """Set the capture branching ratios + + ``branch_ratios`` may be modified in place, only to + insert missing ground state reactions. These will be + inserted only if: + + 1) There is no branch directly to a ground state + target, and + 2) The sum of all ratios on this branch does not + equal 1. + + 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 + parents will be modified. Corresponding values are + dictionaries of ``{target: branching_ratio}`` + + See Also + -------- + :meth:`get_capture_branches` + """ + + # Store some useful information through the validation stage + + sums = {} + capt_ix_map = {} + grounds = {} + + missing = set() + no_capture = set() + + # Check for validity before manipulation + + check_type("branch_ratios", branch_ratios, dict, str) + + for parent, sub in branch_ratios.items(): + if parent not in self: + # TODO How to handle missing branching ratios + missing.add(parent) + continue + + # Make sure this nuclide has capture reactions + + indexes = [] + for ix, rx in enumerate(self[parent].reactions): + if rx.type == "(n,gamma)": + indexes.append(ix) + if "_m" not in rx.target: + grounds[parent] = rx.target + + if len(indexes) == 0: + no_capture.add(parent) + continue + + capt_ix_map[parent] = indexes + + check_type(parent, sub, dict, str) + check_type(parent + " ratios", list(sub.values()), list, float) + this_sum = sum(sub.values()) + check_less_than(parent + " ratios", this_sum, 1.0, True) + sums[parent] = this_sum + + if len(missing) > 0: + print("The following nuclides were not found in {}: {}".format( + self.__class__.__name__, ", ".join(sorted(missing)))) + + if len(no_capture) > 0: + print("The following nuclides did not have capture reactions: " + "{}".format(", ".join(sorted(no_capture)))) + + # Insert new ReactionTuples with updated branch ratios + + for parent_name, capt_index in capt_ix_map.items(): + + parent = self[parent_name] + new_ratios = branch_ratios[parent_name] + capt_index = capt_ix_map[parent_name] + + # Assume Q value is independent of target state + capt_Q = parent.reactions[capt_index[0]].Q + + # Remove existing capture reactions + + for ix in reversed(capt_index): + parent.reactions.pop(ix) + + all_meta = False + + for tgt, br in new_ratios.items(): + all_meta |= ("_m" in tgt) + parent.reactions.append(ReactionTuple( + "(n,gamma)", tgt, capt_Q, br)) + + if all_meta and sums[parent_name] != 1.0: + ground_br = 1.0 - sums[parent_name] + ground_tgt = grounds.get(parent_name, None) + if ground_tgt is None: + pz, pa, pm = zam(parent_name) + 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)) From 5440edd637c7b1cf22c8eb04401e2f2b29c642d2 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 8 Jul 2019 13:00:10 -0500 Subject: [PATCH 2/6] Check that all capture parents and products exist before adding to Chain Before modifying the existing Chain, ensure that all desired products exist. Otherwise, the Chain will be unable to build a depletion matrix in Chain.form_matrix Added an optional argument, strict, that controls the error/print control. If strict, then an KeyError will be raised at the first parent or product that does not exist in the chain. Otherwise, messages will be printed at the end. --- openmc/deplete/chain.py | 49 ++++++++++++++++++++++++++++++++++------- 1 file changed, 41 insertions(+), 8 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index eec9d3f507..0267977d2c 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -480,7 +480,7 @@ class Chain(object): capt[nuclide.name] = nuc_capt return capt - def set_capture_branches(self, branch_ratios): + def set_capture_branches(self, branch_ratios, strict=True): """Set the capture branching ratios ``branch_ratios`` may be modified in place, only to @@ -500,6 +500,13 @@ class Chain(object): ``"Am241"``. The capture 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. See Also -------- @@ -512,7 +519,8 @@ class Chain(object): capt_ix_map = {} grounds = {} - missing = set() + missing_parents = set() + missing_products = {} no_capture = set() # Check for validity before manipulation @@ -521,8 +529,24 @@ class Chain(object): for parent, sub in branch_ratios.items(): if parent not in self: - # TODO How to handle missing branching ratios - missing.add(parent) + if strict: + raise KeyError(parent) + missing_parents.add(parent) + continue + + # Make sure all products are present in the chain + + prod_flag = False + + for product in sub: + if product not in self.nuclide_dict: + if strict: + raise KeyError(product) + missing_products[parent] = product + prod_flag = True + break + + if prod_flag: continue # Make sure this nuclide has capture reactions @@ -535,25 +559,34 @@ class Chain(object): grounds[parent] = rx.target 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) continue capt_ix_map[parent] = indexes - check_type(parent, sub, dict, str) - check_type(parent + " ratios", list(sub.values()), list, float) this_sum = sum(sub.values()) check_less_than(parent + " ratios", this_sum, 1.0, True) sums[parent] = this_sum - if len(missing) > 0: + if len(missing_parents) > 0: print("The following nuclides were not found in {}: {}".format( - self.__class__.__name__, ", ".join(sorted(missing)))) + self.__class__.__name__, ", ".join(sorted(missing_parents)))) if len(no_capture) > 0: print("The following nuclides did not have capture reactions: " "{}".format(", ".join(sorted(no_capture)))) + if len(missing_products) > 0: + tail = ("{} -> {}".format(k, v) + for k, v in sorted(missing_products.items())) + print("The following products were not found in the {} and " + "parents were unmodified: \n{}".format( + self.__class__.__name__, ", ".join(tail))) + # Insert new ReactionTuples with updated branch ratios for parent_name, capt_index in capt_ix_map.items(): From 5b7253bca0f00903c94f24dd4e58fa590e1a0208 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 8 Jul 2019 13:54:05 -0500 Subject: [PATCH 3/6] Add tests for getting, setting capture branching ratios Work with the test chain with isotopes A, B, C to make minor modifications to a depletion chain. Work with the "reference chain" at tests/chain_simple.xml to check inference of ground state, non-construction of reactions that dont' exist. --- tests/unit_tests/test_deplete_chain.py | 67 ++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index dd6817a8e4..2458af7c51 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -243,3 +243,70 @@ def test_set_fiss_q(): for rx in chain_nuc.reactions: if rx.type == 'fission': assert rx.Q == 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 + + # 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 + + # write, re-read + new_chain.export_to_xml("chain_mod.xml") + assert Chain.from_xml("chain_mod.xml").get_capture_branches() == 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 + + # Ensure capture reactions are removed + rem_br = {"A": {"C": 1.0}} + new_chain.set_capture_branches(rem_br) + # A is not in returned dict because there is no branch + assert "A" not in new_chain.get_capture_branches() + + +def test_capture_branch_infer_ground(): + """Ensure the ground state is infered if not given""" + # Make up a metastable capture transition: + infer_br = {"Xe135": {"Xe136_m1": 0.5}} + set_br = {"Xe135": {"Xe136": 0.5, "Xe136_m1": 0.5}} + + chain_file = Path(__file__).parents[1] / "chain_simple.xml" + chain = Chain.from_xml(chain_file) + + # Create nuclide to be added into the chain + xe136m = nuclide.Nuclide() + xe136m.name = "Xe136_m1" + + chain.nuclides.append(xe136m) + chain.nuclide_dict[xe136m.name] = len(chain.nuclides) - 1 + + chain.set_capture_branches(infer_br) + + assert chain.get_capture_branches() == set_br + + +def test_capture_branch_no_rxn(): + """Ensure capture reactions that don't exist aren't created""" + u4br = {"U234": {"U235": 0.5, "U235_m1": 0.5}} + + chain_file = Path(__file__).parents[1] / "chain_simple.xml" + chain = Chain.from_xml(chain_file) + + u5m = nuclide.Nuclide() + u5m.name = "U235_m1" + + 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) From a7b6737069859ba0f1ec8ac578a842453eea3fea Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 8 Jul 2019 14:01:57 -0500 Subject: [PATCH 4/6] Test failure modes for Chain.set_capture_branches Remove one validation check, that the passed item is a dict of strings. The check is covered when every key is inspected to ensure all parents exist in the chain --- openmc/deplete/chain.py | 2 -- tests/unit_tests/test_deplete_chain.py | 19 +++++++++++++++++++ 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 0267977d2c..0a52b78d7f 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -525,8 +525,6 @@ class Chain(object): # Check for validity before manipulation - check_type("branch_ratios", branch_ratios, dict, str) - for parent, sub in branch_ratios.items(): if parent not in self: if strict: diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index 2458af7c51..33b13b596d 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -310,3 +310,22 @@ def test_capture_branch_no_rxn(): phrase = "U234 does not have capture reactions" with pytest.raises(AttributeError, match=phrase): chain.set_capture_branches(u4br) + + +def test_capture_branch_failures(simple_chain): + """Test failure modes for setting capture branch ratios""" + + # Parent isotope not present + br = {"X": {"A": 0.6, "B": 0.7}} + with pytest.raises(KeyError, match="X"): + simple_chain.set_capture_branches(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) + + # 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) From a8dc9dd562fe466b3b7d06216c9532cbc75efa52 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 15 Jul 2019 09:10:13 -0500 Subject: [PATCH 5/6] Use warnings.warn when setting capture branching ratios Remove calls to print for the following cases: 1) Parent nuclides not found in chain 2) Nuclides with requested branching ratios did not have capture reactions 3) Product nuclides not found in chain --- openmc/deplete/chain.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 0a52b78d7f..37c53b8532 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -10,6 +10,7 @@ import math import re 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.data import gnd_name, zam @@ -537,7 +538,7 @@ class Chain(object): prod_flag = False for product in sub: - if product not in self.nuclide_dict: + if product not in self: if strict: raise KeyError(product) missing_products[parent] = product @@ -571,19 +572,19 @@ class Chain(object): sums[parent] = this_sum if len(missing_parents) > 0: - print("The following nuclides were not found in {}: {}".format( - self.__class__.__name__, ", ".join(sorted(missing_parents)))) + warn("The following nuclides were not found in {}: {}".format( + self.__class__.__name__, ", ".join(sorted(missing_parents)))) if len(no_capture) > 0: - print("The following nuclides did not have capture reactions: " - "{}".format(", ".join(sorted(no_capture)))) + warn("The following nuclides did not have capture reactions: " + "{}".format(", ".join(sorted(no_capture)))) if len(missing_products) > 0: tail = ("{} -> {}".format(k, v) for k, v in sorted(missing_products.items())) - print("The following products were not found in the {} and " - "parents were unmodified: \n{}".format( - self.__class__.__name__, ", ".join(tail))) + warn("The following products were not found in the {} and " + "parents were unmodified: \n{}".format( + self.__class__.__name__, ", ".join(tail))) # Insert new ReactionTuples with updated branch ratios From ec2b1816100d97c5838ab54167f123ebafc08dd1 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 15 Jul 2019 09:11:19 -0500 Subject: [PATCH 6/6] Fix all_meta logic in Chain.set_capture_branches --- openmc/deplete/chain.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 37c53b8532..03c7afdd64 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -602,16 +602,16 @@ class Chain(object): for ix in reversed(capt_index): parent.reactions.pop(ix) - all_meta = False + all_meta = True for tgt, br in new_ratios.items(): - all_meta |= ("_m" in tgt) + all_meta = all_meta and ("_m" in tgt) parent.reactions.append(ReactionTuple( "(n,gamma)", tgt, capt_Q, br)) if all_meta and sums[parent_name] != 1.0: ground_br = 1.0 - sums[parent_name] - ground_tgt = grounds.get(parent_name, None) + ground_tgt = grounds.get(parent_name) if ground_tgt is None: pz, pa, pm = zam(parent_name) ground_tgt = gnd_name(pz, pa + 1, 0)