Allow decay only nuclides to be specified in Material for a depletion simulation

openmc/deplete/coupled_operator.py

@@ -10,6 +10,7 @@
 
 import copy
 from warnings import warn
+from typing import List, Dict, Tuple
 
 import numpy as np
 from uncertainties import ufloat
@@ -416,7 +417,81 @@ def _generate_materials_xml(self):
         for mat in self.materials:
             mat._nuclides.sort(key=lambda x: nuclides.index(x[0]))
 
+        # Remove decay-only nuclides
+        #mat_ids = [mat.id for mat in self.materials]
+        material_decay_nuclides = self._remove_decay_nuclides(self.burnable_mats)
+
         self.materials.export_to_xml()
+        if bool(material_decay_nuclides):
+            self._add_decay_nuclides(self.burnable_mats, material_decay_nuclides)
+
+    def _remove_decay_nuclides(
+        self, material_ids: List[str]
+    ) -> Dict[str, List[Tuple[str, float, str]]]:
+        """ Remove any decay-only nuclides (nuclides in the chain file without
+        cross section data) from the specified materials,
+
+        .. versionadded:: 0.13.4
+
+        Parameters
+        ----------
+        material_ids : list of str
+            Material IDs as strings
+
+        Returns
+        -------
+        material_decay_nuclides : dict of str to list of 3-tuples
+            Dictionary mapping material ID's as strings to a list containing
+            decay-only nudclides as 3-tuples of (nuclide, density percent,
+            density percent type)
+
+        """
+
+        material_decay_nuclides = {}
+        for mat in self.materials:
+            mat_id = str(mat.id)
+            if mat_id in material_ids:
+                decay_nuclides = []
+
+                # Remove decay-only nuclides
+                for nuc_tuple in mat._nuclides:
+                    if nuc_tuple.name in self._decay_nucs:
+                        decay_nuclides += [nuc_tuple]
+
+                if len(decay_nuclides) != 0:
+                    for nuc_tuple in decay_nuclides:
+                        mat._nuclides.remove(nuc_tuple)
+                    material_decay_nuclides[mat_id] = decay_nuclides
+
+        return material_decay_nuclides
+
+    def _add_decay_nuclides(
+            self,
+            material_ids: List[str],
+            material_decay_nuclides: Dict[str, List[Tuple[str, float, str]]]):
+        """Add decay-only nuclides to the specified materials
+
+        .. versionadded:: 0.13.4
+
+        Parameters
+        ----------
+        material_ids : list of str
+            Material IDs as strings
+
+        material_decay_nuclides : dict of str to list of 3-tuples
+            Dictionary mapping material ID's as strings to a list containing
+            decay-only nudclides as 3-tuples of (nuclide, density percent,
+            density percent type)
+
+        """
+        # Add decay-only nuclides back in
+        nuclides = list(self.number.nuclides)
+        for mat in self.materials:
+            mat_id = str(mat.id)
+            if mat_id in material_ids:
+                mat._nuclides += material_decay_nuclides[mat_id]
+                # Sort nuclides according to order in AtomNumber object
+                mat._nuclides.sort(key=lambda x: nuclides.index(x[0]))
 
     def __call__(self, vec, source_rate):
         """Runs a simulation.

openmc/deplete/openmc_operator.py

@@ -6,6 +6,7 @@
 """
 
 from abc import abstractmethod
+from warnings import warn
 
 import numpy as np
 
@@ -132,6 +133,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)
 
@@ -141,13 +154,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,
@@ -187,7 +193,13 @@ def _get_burnable_mats(self):
         # 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:

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>

tests/regression_tests/deplete_decay_only/test.py

@@ -0,0 +1,109 @@
+""" 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]]),
+                                 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)