diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 6ba481504a..2405a67a27 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -47,6 +47,7 @@ _REACTIONS = [ __all__ = ["Chain"] + def replace_missing(product, decay_data): """Replace missing product with suitable decay daughter. @@ -109,6 +110,54 @@ def replace_missing(product, decay_data): return product +def replace_missing_fpy(actinide, fpy_data, decay_data): + """Replace missing fission product yields + + Parameters + ---------- + actinide : str + Name of actinide missing FPY data + fpy_data : dict + Dictionary of FPY data + decay_data : dict + Dictionary of decay data + + Returns + ------- + str + Actinide that can be used as replacement for FPY purposes + + """ + + # Check if metastable state has data (e.g., Am242m) + Z, A, m = zam(actinide) + if m == 0: + metastable = gnd_name(Z, A, 1) + if metastable in fpy_data: + return metastable + + # Try increasing Z, holding N constant + isotone = actinide + while isotone in decay_data: + Z += 1 + A += 1 + isotone = gnd_name(Z, A, 0) + if isotone in fpy_data: + return isotone + + # Try decreasing Z, holding N constant + isotone = actinide + while isotone in decay_data: + Z -= 1 + A -= 1 + isotone = gnd_name(Z, A, 0) + if isotone in fpy_data: + return isotone + + # If all else fails, use U235 yields + return 'U235' + + _SECONDARY_PARTICLES = { '(n,p)': ['H1'], '(n,d)': ['H2'], @@ -290,6 +339,7 @@ class Chain: # Append decay mode nuclide.decay_modes.append(DecayTuple(type_, target, br)) + fissionable = False if parent in reactions: reactions_available = set(reactions[parent].keys()) for name, mts, changes in _REACTIONS: @@ -317,18 +367,21 @@ class Chain: name, daughter, q_value, 1.0)) if any(mt in reactions_available for mt in openmc.data.FISSION_MTS): - if parent in fpy_data: - q_value = reactions[parent][18] - nuclide.reactions.append( - ReactionTuple('fission', 0, q_value, 1.0)) + q_value = reactions[parent][18] + nuclide.reactions.append( + ReactionTuple('fission', None, q_value, 1.0)) + fissionable = True - if 'fission' not in chain.reactions: - chain.reactions.append('fission') - else: - missing_fpy.append(parent) + if 'fission' not in chain.reactions: + chain.reactions.append('fission') - if parent in fpy_data: - fpy = fpy_data[parent] + if fissionable: + if parent in fpy_data: + fpy = fpy_data[parent] + else: + nuclide._fpy = replace_missing_fpy(parent, fpy_data, decay_data) + missing_fpy.append((parent, nuclide._fpy)) + continue if fpy.energies is not None: yield_energies = fpy.energies @@ -354,6 +407,11 @@ class Chain: nuclide.yield_data = FissionYieldDistribution(yield_data) + # Replace missing FPY data + for nuclide in chain.nuclides: + if hasattr(nuclide, '_fpy'): + nuclide.yield_data = chain[nuclide._fpy].yield_data + # Display warnings if missing_daughter: print('The following decay modes have daughters with no decay data:') @@ -369,8 +427,8 @@ class Chain: if missing_fpy: print('The following fissionable nuclides have no fission product yields:') - for parent in missing_fpy: - print(' ' + parent) + for parent, replacement in missing_fpy: + print(' {}, replaced with {}'.format(parent, replacement)) print('') if missing_fp: @@ -406,7 +464,7 @@ class Chain: for i, nuclide_elem in enumerate(root.findall('nuclide')): this_q = fission_q.get(nuclide_elem.get("name")) - nuc = Nuclide.from_xml(nuclide_elem, this_q) + nuc = Nuclide.from_xml(nuclide_elem, root, this_q) chain.nuclide_dict[nuc.name] = i # Check for reaction paths diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 7fb032aaa0..0b089e0e2c 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -155,13 +155,15 @@ class Nuclide: return self.yield_data.energies @classmethod - def from_xml(cls, element, fission_q=None): + def from_xml(cls, element, root, fission_q=None): """Read nuclide from an XML element. Parameters ---------- element : xml.etree.ElementTree.Element - XML element to write nuclide data to + XML element to read nuclide data from + root : xml.etree.ElementTree.Element + Root XML element for chain file fission_q : None or float User-supplied fission Q value [eV]. Will be read from the element if not given @@ -212,6 +214,14 @@ class Nuclide: fpy_elem = element.find('neutron_fission_yields') if fpy_elem is not None: + # Check for use of FPY from other nuclide + parent = fpy_elem.get('parent') + if parent is not None: + fpy_elem = root.find( + './/nuclide[@name="{}"]/neutron_fission_yields'.format(parent) + ) + nuc._fpy = parent + nuc.yield_data = FissionYieldDistribution.from_xml_element(fpy_elem) return nuc @@ -250,9 +260,14 @@ class Nuclide: if self.yield_data: fpy_elem = ET.SubElement(elem, 'neutron_fission_yields') - energy_elem = ET.SubElement(fpy_elem, 'energies') - energy_elem.text = ' '.join(str(E) for E in self.yield_energies) - self.yield_data.to_xml_element(fpy_elem) + + if hasattr(self, '_fpy'): + # Check for link to other nuclide data + fpy_elem.set('parent', self._fpy) + else: + energy_elem = ET.SubElement(fpy_elem, 'energies') + energy_elem.text = ' '.join(str(E) for E in self.yield_energies) + self.yield_data.to_xml_element(fpy_elem) return elem