Allow decay only nuclides to be specified in Material for a depleti…

…on simulation (#2616) Co-authored-by: Jonathan Shimwell <drshimwell@gmail.com> Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
openmc-dev · Sep 14, 2023 · bbd7756 · bbd7756
1 parent 0cf29f1
commit bbd7756
Show file tree

Hide file tree

Showing 7 changed files with 233 additions and 20 deletions.
diff --git a/openmc/deplete/coupled_operator.py b/openmc/deplete/coupled_operator.py
@@ -416,7 +416,7 @@ def _generate_materials_xml(self):
         for mat in self.materials:
             mat._nuclides.sort(key=lambda x: nuclides.index(x[0]))
 
-        self.materials.export_to_xml()
+        self.materials.export_to_xml(nuclides_to_ignore=self._decay_nucs)
 
     def __call__(self, vec, source_rate):
         """Runs a simulation.

diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py
@@ -6,6 +6,7 @@
 """
 
 from abc import abstractmethod
+from warnings import warn
 from typing import List, Tuple, Dict
 
 import numpy as np
@@ -133,6 +134,18 @@ def __init__(
         # This nuclides variables contains every nuclides
         # for which there is an entry in the micro_xs parameter
         openmc.reset_auto_ids()
+
+        self.nuclides_with_data = self._get_nuclides_with_data(
+            self.cross_sections)
+
+        # Select nuclides with data that are also in the chain
+        self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides
+                               if nuc.name in self.nuclides_with_data]
+
+        # Select nuclides without data that are also in the chain
+        self._decay_nucs = [nuc.name for nuc in self.chain.nuclides
+                            if nuc.name not in self.nuclides_with_data]
+
         self.burnable_mats, volumes, all_nuclides = self._get_burnable_mats()
         self.local_mats = _distribute(self.burnable_mats)
 
@@ -142,13 +155,6 @@ def __init__(
         if self.prev_res is not None:
             self._load_previous_results()
 
-        self.nuclides_with_data = self._get_nuclides_with_data(
-            self.cross_sections)
-
-        # Select nuclides with data that are also in the chain
-        self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides
-                               if nuc.name in self.nuclides_with_data]
-
         # Extract number densities from the geometry / previous depletion run
         self._extract_number(self.local_mats,
                              volumes,
@@ -171,7 +177,7 @@ def _get_burnable_mats(self) -> Tuple[List[str], Dict[str, float], List[str]]:
         Returns
         -------
         burnable_mats : list of str
-            List of burnable material IDs
+            list of burnable material IDs
         volume : dict of str to float
             Volume of each material in [cm^3]
         nuclides : list of str
@@ -188,7 +194,13 @@ def _get_burnable_mats(self) -> Tuple[List[str], Dict[str, float], List[str]]:
         # Iterate once through the geometry to get dictionaries
         for mat in self.materials:
             for nuclide in mat.get_nuclides():
-                model_nuclides.add(nuclide)
+                if nuclide in self.nuclides_with_data or self._decay_nucs:
+                    model_nuclides.add(nuclide)
+                else:
+                    msg = (f"Nuclilde {nuclide} in material {mat.id} is not "
+                           "present in the depletion chain and has no cross "
+                           "section data.")
+                    raise warn(msg)
             if mat.depletable:
                 burnable_mats.add(str(mat.id))
                 if mat.volume is None:

diff --git a/openmc/material.py b/openmc/material.py
@@ -1250,7 +1250,7 @@ def clone(self, memo: Optional[dict] = None) -> Material:
 
         return memo[self]
 
-    def _get_nuclide_xml(self, nuclide: str) -> ET.Element:
+    def _get_nuclide_xml(self, nuclide: NuclideTuple) -> ET.Element:
         xml_element = ET.Element("nuclide")
         xml_element.set("name", nuclide.name)
 
@@ -1267,15 +1267,28 @@ def _get_macroscopic_xml(self, macroscopic: str) -> ET.Element:
 
         return xml_element
 
-    def _get_nuclides_xml(self, nuclides: typing.Iterable[str]) -> List[ET.Element]:
+    def _get_nuclides_xml(
+            self, nuclides: typing.Iterable[NuclideTuple],
+            nuclides_to_ignore: Optional[typing.Iterable[str]] = None)-> List[ET.Element]:
         xml_elements = []
-        for nuclide in nuclides:
-            xml_elements.append(self._get_nuclide_xml(nuclide))
+
+        # Remove any nuclides to ignore from the XML export
+        if nuclides_to_ignore:
+            nuclides = [nuclide for nuclide in nuclides if nuclide.name not in nuclides_to_ignore]
+
+        xml_elements = [self._get_nuclide_xml(nuclide) for nuclide in nuclides]
+
         return xml_elements
 
-    def to_xml_element(self) -> ET.Element:
+    def to_xml_element(
+            self, nuclides_to_ignore: Optional[typing.Iterable[str]] = None) -> ET.Element:
         """Return XML representation of the material
 
+        Parameters
+        ----------
+        nuclides_to_ignore : list of str
+            Nuclides to ignore when exporting to XML.
+
         Returns
         -------
         element : lxml.etree._Element
@@ -1320,7 +1333,8 @@ def to_xml_element(self) -> ET.Element:
 
         if self._macroscopic is None:
             # Create nuclide XML subelements
-            subelements = self._get_nuclides_xml(self._nuclides)
+            subelements = self._get_nuclides_xml(self._nuclides,
+                                                 nuclides_to_ignore=nuclides_to_ignore)
             for subelement in subelements:
                 element.append(subelement)
         else:
@@ -1576,7 +1590,8 @@ def make_isotropic_in_lab(self):
         for material in self:
             material.make_isotropic_in_lab()
 
-    def _write_xml(self, file, header=True, level=0, spaces_per_level=2, trailing_indent=True):
+    def _write_xml(self, file, header=True, level=0, spaces_per_level=2,
+                   trailing_indent=True, nuclides_to_ignore=None):
         """Writes XML content of the materials to an open file handle.
 
         Parameters
@@ -1591,6 +1606,8 @@ def _write_xml(self, file, header=True, level=0, spaces_per_level=2, trailing_in
             Number of spaces per indentation
         trailing_indentation : bool
             Whether or not to write a trailing indentation for the materials element
+        nuclides_to_ignore : list of str
+            Nuclides to ignore when exporting to XML.
 
         """
         indentation = level*spaces_per_level*' '
@@ -1611,7 +1628,7 @@ def _write_xml(self, file, header=True, level=0, spaces_per_level=2, trailing_in
 
         # Write the <material> elements.
         for material in sorted(self, key=lambda x: x.id):
-            element = material.to_xml_element()
+            element = material.to_xml_element(nuclides_to_ignore=nuclides_to_ignore)
             clean_indentation(element, level=level+1)
             element.tail = element.tail.strip(' ')
             file.write((level+1)*spaces_per_level*' ')
@@ -1626,13 +1643,16 @@ def _write_xml(self, file, header=True, level=0, spaces_per_level=2, trailing_in
         if trailing_indent:
             file.write(indentation)
 
-    def export_to_xml(self, path: PathLike = 'materials.xml'):
+    def export_to_xml(self, path: PathLike = 'materials.xml',
+                      nuclides_to_ignore: Optional[typing.Iterable[str]] = None):
         """Export material collection to an XML file.
 
         Parameters
         ----------
         path : str
             Path to file to write. Defaults to 'materials.xml'.
+        nuclides_to_ignore : list of str
+            Nuclides to ignore when exporting to XML.
 
         """
         # Check if path is a directory
@@ -1645,7 +1665,7 @@ def export_to_xml(self, path: PathLike = 'materials.xml'):
         # one go.
         with open(str(p), 'w', encoding='utf-8',
                   errors='xmlcharrefreplace') as fh:
-            self._write_xml(fh)
+            self._write_xml(fh, nuclides_to_ignore=nuclides_to_ignore)
 
     @classmethod
     def from_xml_element(cls, elem) -> Material:

diff --git a/tests/chain_simple_decay.xml b/tests/chain_simple_decay.xml
@@ -0,0 +1,66 @@
+<?xml version="1.0"?>
+<depletion_chain>
+  <nuclide name="I135" decay_modes="1" reactions="1" half_life="2.36520E+04" decay_energy="1916827.5">
+    <decay type="beta" target="Xe135" branching_ratio="1.0" />
+    <reaction type="(n,gamma)" Q="0.0" target="Xe136" /> <!-- Not precisely true, but whatever -->
+    <source type="discrete" particle="photon">
+      <parameters>3696.125 4095.822 4477.27 5097.122 29452.1 29781.3 33566.5 33629.4 33865.1 33878.5 34395.3 34408.0 34486.2 34488.2 112780.0 113150.0 162650.0 165740.0 184490.0 197190.0 220502.0 229720.0 247500.0 254740.0 264260.0 288451.0 290270.0 304910.0 305830.0 326000.0 333600.0 342520.0 361850.0 403030.0 414830.0 417633.0 429930.0 433741.0 451630.0 530800.0 546557.0 575970.0 588280.0 616900.0 649850.0 656090.0 679220.0 684600.0 690130.0 707920.0 785480.0 795500.0 797710.0 807200.0 836804.0 960290.0 961430.0 971960.0 972620.0 995090.0 1038760.0 1096860.0 1101580.0 1124000.0 1131511.0 1151510.0 1159900.0 1169040.0 1225600.0 1240470.0 1254800.0 1260409.0 1315770.0 1334800.0 1343660.0 1367890.0 1441800.0 1448350.0 1457560.0 1502790.0 1521990.0 1543700.0 1566410.0 1678027.0 1706459.0 1791196.0 1830690.0 1845300.0 1927300.0 1948490.0 2045880.0 2112400.0 2151500.0 2189400.0 2255457.0 2408650.0 2466070.0 2477100.0 9.714352819815078e-10 7.460941651551526e-09 6.047882056201745e-09 8.510389107205747e-10 3.7979729633684727e-08 7.033747061198817e-08 6.602815946851458e-09 1.2800909198245071e-08 6.513060244589454e-11 8.894667887850525e-11 1.3843783302275766e-09 2.700005498135763e-09 1.2141115256573815e-11 1.654893875618862e-11 3.7007705885806645e-09 2.0186021392258173e-09 2.859686363903241e-09 9.16781804898392e-09 6.896890642354875e-09 9.588360161322631e-09 5.130613770532286e-07 7.06510748729036e-08 8.410842246774237e-09 6.7286737974193905e-09 5.3829390379355124e-08 9.083709626516178e-07 8.915492781580693e-08 9.251926471451662e-09 2.783988783682273e-08 6.728673797419391e-10 1.093409492080651e-08 2.5232526740322717e-10 5.467047460403255e-08 6.812782219887132e-08 8.83138435911295e-08 1.0345335963532313e-06 8.915492781580693e-08 1.623292553627428e-07 9.251926471451663e-08 9.251926471451662e-09 2.0942997194467853e-06 3.784879011048407e-08 1.513951604419363e-08 1.093409492080651e-08 1.337323917237104e-07 2.186818984161302e-08 1.5980600268871052e-08 6.7286737974193905e-09 3.784879011048407e-08 1.9344937167580748e-07 4.457746390790346e-08 6.7286737974193905e-09 5.0465053480645434e-08 1.3457347594838781e-08 1.9597262434983976e-06 1.0093010696129087e-08 4.289529545854862e-08 2.607361096500014e-07 3.53255374364518e-07 4.541854813258089e-08 2.329803302356464e-06 2.607361096500014e-08 4.7100716581935733e-07 1.059766123093554e-06 6.6193328482113254e-06 4.205421123387119e-10 3.027903208838726e-08 2.565306885266143e-07 1.2616263370161358e-08 2.649415307733885e-07 3.3643368987096953e-09 8.410842246774237e-06 1.934493716758075e-08 9.251926471451662e-09 2.2709274066290446e-08 1.7830985563161385e-07 5.046505348064544e-09 9.251926471451663e-08 2.5400743585258203e-06 3.154065842540339e-07 1.093409492080651e-08 7.569758022096815e-09 3.784879011048407e-07 2.8008104681758214e-06 1.2027504412887161e-06 2.26251656438227e-06 1.6989901338483962e-07 1.6821684493548476e-09 8.663167514177466e-08 1.8503852942903323e-08 2.5568960430193683e-07 2.0186021392258175e-08 6.560456952483906e-09 3.784879011048407e-09 1.7999202408096872e-07 2.800810468175822e-07 2.1027105616935597e-08 4.205421123387119e-10</parameters>
+    </source>
+  </nuclide>
+  <nuclide name="Xe135" reactions="1" decay_energy="567890.1">
+    <reaction type="(n,gamma)" Q="0.0" target="Xe136" />
+    <source type="tabular" interpolation="histogram" particle="photon">
+      <parameters>0.0 10000.0 20000.0 50000.0 100000.0 200000.0 300000.0 400000.0 600000.0 800000.0 1000000.0 1220000.0 1.1612176249914943e-11 0.0 3.1976524142990123e-11 0.0 6.418466676129901e-13 1.894551923065874e-10 5.099949011192198e-13 3.170644340540729e-13 2.756798470867507e-12 6.67627508515641e-14 3.711746246444946e-15 0.0</parameters>
+    </source>
+  </nuclide>
+  <nuclide name="Xe135_m1" half_life="917.4" decay_modes="1" decay_energy="526729.1900000001" reactions="0">
+    <decay type="IT" target="Xe135" branching_ratio="1"/>
+  </nuclide>
+  <nuclide name="Xe136" decay_modes="0" reactions="0" />
+  <nuclide name="Cs135" decay_modes="0" reactions="0" />
+  <nuclide name="Cs135_m1" half_life="3180.0" decay_modes="1" decay_energy="1633299.89" reactions="0">
+    <decay type="IT" target="Cs135" branching_ratio="1.0"/>
+  </nuclide>
+  <nuclide name="Gd157" decay_modes="0" reactions="1"  >
+    <reaction type="(n,gamma)" Q="0.0" target="Nothing" />
+  </nuclide>
+  <nuclide name="Gd156" decay_modes="0" reactions="1">
+    <reaction type="(n,gamma)" Q="0.0" target="Gd157" />
+  </nuclide>
+  <nuclide name="U234" decay_modes="0" reactions="1">
+    <reaction type="fission" Q="191840000."/>
+    <neutron_fission_yields>
+      <energies>2.53000e-02</energies>
+      <fission_yields energy="2.53000e-02">
+        <products>Gd157 Gd156 I135 Xe135 Xe136 Cs135</products>
+        <data>1.093250e-04 2.087260e-04 2.780820e-02 6.759540e-03 2.392300e-02 4.356330e-05</data>
+      </fission_yields>
+    </neutron_fission_yields>
+  </nuclide>
+  <nuclide name="U235" decay_modes="0" reactions="1">
+    <reaction type="fission" Q="193410000."/>
+    <source type="discrete" particle="photon">
+      <parameters>3065.349 12960.11 13197.49 16125.43 19185.05 19590.0 31600.0 34700.0 41400.0 41960.0 51220.0 54100.0 54250.0 64350.0 72700.0 75020.0 76198.0 90330.0 93795.0 96090.0 105278.0 106074.0 106608.0 106771.0 108948.0 109154.0 109160.0 109395.0 109433.0 115450.0 120350.0 136550.0 140760.0 142400.0 143760.0 150930.0 163330.0 173300.0 182610.0 185715.0 194940.0 198900.0 202110.0 205311.0 215280.0 221380.0 228780.0 233500.0 240870.0 246840.0 266450.0 275129.0 275430.0 281420.0 282920.0 289560.0 291650.0 301700.0 317100.0 343500.0 345900.0 356030.0 387820.0 410290.0 428710.0 448400.0 1.4211820389290126e-18 3.695705148646686e-18 4.752472005970374e-19 4.0825252294602915e-18 8.998127221991115e-19 1.9035285506819052e-21 5.305446177665699e-21 1.1547147563154756e-20 9.362552078233586e-21 1.86835843588509e-20 1.0610892355331398e-20 2.736290732002608e-22 4.776811520523089e-21 3.977552295559136e-21 3.432935762018982e-20 1.872510415646717e-20 2.4966805541956234e-21 1.0265454754703584e-18 1.7575878976520235e-18 2.839974130397521e-20 2.1057468800839102e-19 4.1342252420362975e-19 7.098565272539688e-21 8.20050332279016e-21 5.276915360632628e-20 1.0602918581811435e-19 4.806110066826575e-19 2.0521431485853304e-21 2.375669880669523e-21 9.362552078233586e-21 8.267152210184412e-21 3.745020831293435e-21 6.865871524037964e-20 1.5604253463722645e-21 3.420452359248004e-18 2.4966805541956232e-20 1.5853921519142206e-18 1.8725104156467174e-21 1.0610892355331397e-19 1.7851265962498703e-17 1.966135936429053e-19 1.3107572909527022e-20 3.370518748164091e-19 1.5635461970650089e-18 9.050467008959134e-21 3.745020831293434e-20 2.18459548492117e-21 9.050467008959134e-21 2.3406380195583967e-20 1.6540508671546e-20 1.8725104156467174e-21 1.622842360227155e-20 2.18459548492117e-21 1.8725104156467174e-21 1.8725104156467174e-21 2.18459548492117e-21 1.2483402770978116e-20 1.5604253463722645e-21 3.120850692744529e-22 9.362552078233587e-22 1.2483402770978116e-20 1.5604253463722645e-21 1.2483402770978116e-20 9.362552078233587e-22 3.120850692744529e-22 3.120850692744529e-22</parameters>
+    </source>
+    <neutron_fission_yields>
+      <energies>2.53000e-02</energies>
+      <fission_yields energy="2.53000e-02">
+        <products>Gd157 Gd156 I135 Xe135 Xe136 Cs135</products>
+        <data>6.142710e-5 1.483250e-04 0.0292737 0.002566345 0.0219242 4.9097e-6</data>
+      </fission_yields>
+    </neutron_fission_yields>
+  </nuclide>
+  <nuclide name="U238" decay_modes="0" reactions="1">
+    <reaction type="fission" Q="197790000."/>
+    <source type="discrete" particle="photon">
+      <parameters>3061.32 12959.8 13440.07 16150.05 19148.83 49550.0 90330.0 93795.0 105278.0 106074.0 106608.0 106771.0 108948.0 109154.0 109395.0 109433.0 113500.0 5.3130501476983077e-20 1.4936931167066328e-19 1.943509007980312e-20 1.8224649880365157e-19 4.0452400597373603e-20 3.146222282132249e-21 3.300026341588747e-23 5.444173877278485e-23 6.3486292118621375e-24 1.2400176384733938e-23 2.124537722121886e-25 2.4542156071495764e-25 1.5792541400719878e-24 3.1731991718126067e-24 6.141568457919309e-26 7.109808533300877e-26 5.014291762148272e-22</parameters>
+    </source>
+    <neutron_fission_yields>
+      <energies>2.53000e-02</energies>
+      <fission_yields energy="2.53000e-02">
+        <products>Gd157 Gd156 I135 Xe135 Xe136 Cs135</products>
+        <data>4.141120e-04 7.605360e-04 0.0135457 0.00026864 0.0024432 3.7100E-07</data>
+      </fission_yields>
+    </neutron_fission_yields>
+  </nuclide>
+</depletion_chain>
diff --git a/tests/regression_tests/deplete_decay_only/__init__.py b/tests/regression_tests/deplete_decay_only/__init__.py
diff --git a/tests/regression_tests/deplete_decay_only/test.py b/tests/regression_tests/deplete_decay_only/test.py
@@ -0,0 +1,103 @@
+""" Transport-free depletion test suite """
+
+from pathlib import Path
+import shutil
+
+import numpy as np
+import pytest
+import openmc
+import openmc.deplete
+from openmc.deplete import CoupledOperator, IndependentOperator, MicroXS
+
+
+@pytest.fixture(scope="module")
+def model():
+    fuel = openmc.Material(name="uo2")
+    fuel.add_element("U", 1, percent_type="ao", enrichment=4.25)
+    fuel.add_element("O", 2)
+    fuel.add_nuclide("Xe135_m1", 1)
+    fuel.add_nuclide("Cs135_m1", 1)
+    fuel.set_density("g/cc", 10.4)
+
+    clad = openmc.Material(name="clad")
+    clad.add_element("Zr", 1)
+    clad.set_density("g/cc", 6)
+
+    water = openmc.Material(name="water")
+    water.add_element("O", 1)
+    water.add_element("H", 2)
+    water.set_density("g/cc", 1.0)
+    water.add_s_alpha_beta("c_H_in_H2O")
+
+    radii = [0.42, 0.45]
+    fuel.volume = np.pi * radii[0] ** 2
+
+    materials = openmc.Materials([fuel, clad, water])
+
+    pin_surfaces = [openmc.ZCylinder(r=r) for r in radii]
+    pin_univ = openmc.model.pin(pin_surfaces, materials)
+    bound_box = openmc.rectangular_prism(1.24, 1.24, boundary_type="reflective")
+    root_cell = openmc.Cell(fill=pin_univ, region=bound_box)
+    geometry = openmc.Geometry([root_cell])
+
+    settings = openmc.Settings()
+    settings.particles = 1000
+    settings.inactive = 5
+    settings.batches = 10
+
+    return openmc.Model(geometry, materials, settings)
+
+@pytest.fixture(scope="module")
+def micro_xs():
+    micro_xs_file = Path(__file__).parents[2] / 'micro_xs_simple.csv'
+    return MicroXS.from_csv(micro_xs_file)
+
+
+@pytest.fixture(scope="module")
+def chain_file():
+    return Path(__file__).parents[2] / 'chain_simple_decay.xml'
+
+
+@pytest.mark.parametrize("operator_type", ["coupled", "independent"])
+def test_decay_only(run_in_tmpdir, operator_type, model, micro_xs, chain_file):
+    """Transport free system test suite.
+
+    """
+    # Create operator
+    if operator_type == "coupled":
+        op = CoupledOperator(model, chain_file=chain_file)
+    else:
+        op = IndependentOperator(openmc.Materials([model.materials[0]]),
+                                 [1e15],
+                                 [micro_xs],
+                                 chain_file)
+
+    # Power and timesteps
+    dt = [917.4, 2262.6]  # one Xe135_m1 half life and one Cs135_m1 half life
+
+    # Perform simulation using the predictor algorithm
+    openmc.deplete.PredictorIntegrator(op,
+                                       dt,
+                                       power=0.0,
+                                       timestep_units='s').integrate()
+
+    # Get path to test and reference results
+    path_test = op.output_dir / 'depletion_results.h5'
+
+    # Load the reference/test results
+    res_test = openmc.deplete.Results(path_test)
+
+    _, xe135m1_atoms = res_test.get_atoms('1', 'Xe135_m1')
+    _, xe135_atoms = res_test.get_atoms('1', 'Xe135')
+    _, cs135m1_atoms = res_test.get_atoms('1', 'Cs135_m1')
+    _, cs135_atoms = res_test.get_atoms('1', 'Cs135')
+
+    tol = 1.0e-14
+    assert xe135m1_atoms[0] == pytest.approx(xe135m1_atoms[1] * 2, rel=tol)
+
+    # WARNING: this is generally not true as Xe135_m1 has two
+    # decay modes, and Xe135 will also decay, but we've modified the depletion chain so
+    # that Xe135_m1 only decays to Xe135, and that Xe135 has has no decay modes
+    assert xe135_atoms[1] == pytest.approx(xe135m1_atoms[1], rel=tol)
+    assert cs135m1_atoms[0] == pytest.approx(cs135m1_atoms[2] * 2, rel=tol)
+    assert cs135_atoms[2] == pytest.approx(cs135m1_atoms[2], rel=tol)