Make sure actinides missing FPY data are assigned something else

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

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