Merge pull request #2301 from paulromano/decay-products

Produce light particles from decay

openmc/deplete/chain.py

@@ -630,15 +630,27 @@ class Chain:
 
             # Gain from radioactive decay
             if nuc.n_decay_modes != 0:
-                for _, target, branching_ratio in nuc.decay_modes:
-                    # Allow for total annihilation for debug purposes
-                    if target is not None:
-                        branch_val = branching_ratio * decay_constant
+                for decay_type, target, branching_ratio in nuc.decay_modes:
+                    branch_val = branching_ratio * decay_constant
 
-                        if branch_val != 0.0:
+                    # Allow for total annihilation for debug purposes
+                    if branch_val != 0.0:
+                        if target is not None:
                             k = self.nuclide_dict[target]
                             matrix[k, i] += branch_val
 
+                        # Produce alphas and protons from decay
+                        if 'alpha' in decay_type:
+                            k = self.nuclide_dict.get('He4')
+                            if k is not None:
+                                count = decay_type.count('alpha')
+                                matrix[k, i] += count * branch_val
+                        elif 'p' in decay_type:
+                            k = self.nuclide_dict.get('H1')
+                            if k is not None:
+                                count = decay_type.count('p')
+                                matrix[k, i] += count * branch_val
+
             if nuc.name in rates.index_nuc:
                 # Extract all reactions for this nuclide in this cell
                 nuc_ind = rates.index_nuc[nuc.name]

openmc/deplete/independent_operator.py

@@ -31,6 +31,9 @@ class IndependentOperator(OpenMCOperator):
     passed to an integrator class, such as
     :class:`openmc.deplete.CECMIntegrator`.
 
+    Note that passing an empty :class:`~openmc.deplete.MicroXS` instance to the
+    ``micro_xs`` argument allows a decay-only calculation to be run.
+
     .. versionadded:: 0.13.1
 
     Parameters
@@ -38,9 +41,11 @@ class IndependentOperator(OpenMCOperator):
     materials : openmc.Materials
         Materials to deplete.
     micro_xs : MicroXS
-        One-group microscopic cross sections in [b] .
+        One-group microscopic cross sections in [b]. If the
+        :class:`~openmc.deplete.MicroXS` object is empty, a decay-only calculation will
+        be run.
     chain_file : str
-        Path to the depletion chain XML file.  Defaults to
+        Path to the depletion chain XML file. Defaults to
         ``openmc.config['chain_file']``.
     keff : 2-tuple of float, optional
        keff eigenvalue and uncertainty from transport calculation.
@@ -150,7 +155,7 @@ class IndependentOperator(OpenMCOperator):
     @classmethod
     def from_nuclides(cls, volume, nuclides,
                       micro_xs,
-                      chain_file,
+                      chain_file=None,
                       nuc_units='atom/b-cm',
                       keff=None,
                       normalization_mode='fission-q',
@@ -169,10 +174,13 @@ class IndependentOperator(OpenMCOperator):
             Dictionary with nuclide names as keys and nuclide concentrations as
             values.
         micro_xs : MicroXS
-            One-group microscopic cross sections.
-        chain_file : str
-            Path to the depletion chain XML file.
-        nuc_units : {'atom/cm3', 'atom/b-cm'}
+            One-group microscopic cross sections in [b]. If the
+            :class:`~openmc.deplete.MicroXS` object is empty, a decay-only calculation
+            will be run.
+        chain_file : str, optional
+            Path to the depletion chain XML file. Defaults to
+            ``openmc.config['chain_file']``.
+        nuc_units : {'atom/cm3', 'atom/b-cm'}, optional
             Units for nuclide concentration.
         keff : 2-tuple of float, optional
            keff eigenvalue and uncertainty from transport calculation.

openmc/deplete/stepresult.py

@@ -318,10 +318,11 @@ class StepResult:
                               chunks=(1, 1, n_mats, n_nuc_number),
                               dtype='float64')
 
-        handle.create_dataset("reaction rates", (1, n_stages, n_mats, n_nuc_rxn, n_rxn),
-                              maxshape=(None, n_stages, n_mats, n_nuc_rxn, n_rxn),
-                              chunks=(1, 1, n_mats, n_nuc_rxn, n_rxn),
-                              dtype='float64')
+        if n_nuc_rxn > 0 and n_rxn > 0:
+            handle.create_dataset("reaction rates", (1, n_stages, n_mats, n_nuc_rxn, n_rxn),
+                                maxshape=(None, n_stages, n_mats, n_nuc_rxn, n_rxn),
+                                chunks=(1, 1, n_mats, n_nuc_rxn, n_rxn),
+                                dtype='float64')
 
         handle.create_dataset("eigenvalues", (1, n_stages, 2),
                               maxshape=(None, n_stages, 2), dtype='float64')
@@ -358,7 +359,9 @@ class StepResult:
 
         # Grab handles
         number_dset = handle["/number"]
-        rxn_dset = handle["/reaction rates"]
+        has_reactions = ("reaction rates" in handle)
+        if has_reactions:
+            rxn_dset = handle["/reaction rates"]
         eigenvalues_dset = handle["/eigenvalues"]
         time_dset = handle["/time"]
         source_rate_dset = handle["/source_rate"]
@@ -375,9 +378,10 @@ class StepResult:
             number_shape[0] = new_shape
             number_dset.resize(number_shape)
 
-            rxn_shape = list(rxn_dset.shape)
-            rxn_shape[0] = new_shape
-            rxn_dset.resize(rxn_shape)
+            if has_reactions:
+                rxn_shape = list(rxn_dset.shape)
+                rxn_shape[0] = new_shape
+                rxn_dset.resize(rxn_shape)
 
             eigenvalues_shape = list(eigenvalues_dset.shape)
             eigenvalues_shape[0] = new_shape
@@ -407,7 +411,8 @@ class StepResult:
         high = max(inds)
         for i in range(n_stages):
             number_dset[index, i, low:high+1] = self.data[i]
-            rxn_dset[index, i, low:high+1] = self.rates[i]
+            if has_reactions:
+                rxn_dset[index, i, low:high+1] = self.rates[i]
             if comm.rank == 0:
                 eigenvalues_dset[index, i] = self.k[i]
         if comm.rank == 0:
@@ -483,7 +488,8 @@ class StepResult:
         for i in range(results.n_stages):
             rate = ReactionRates(results.index_mat, rxn_nuc_to_ind, rxn_to_ind, True)
 
-            rate[:] = handle["/reaction rates"][step, i, :, :, :]
+            if "reaction rates" in handle:
+                rate[:] = handle["/reaction rates"][step, i, :, :, :]
             results.rates.append(rate)
 
         return results

tests/unit_tests/test_deplete_decay_products.py

@@ -0,0 +1,43 @@
+import openmc.deplete
+import pytest
+
+
+def test_deplete_decay_products(run_in_tmpdir):
+    # Create chain file with H1, He4, and Li5
+    with open('test_chain.xml', 'w') as chain_file:
+        chain_file.write("""
+<depletion_chain>
+  <nuclide name="H1" reactions="0">
+  </nuclide>
+  <nuclide name="He4" reactions="0"/>
+  <nuclide name="Li5" half_life="3.06868e-22" decay_modes="1" decay_energy="1965000.0" reactions="0">
+    <decay type="alpha" target="H1" branching_ratio="1.0"/>
+  </nuclide>
+</depletion_chain>
+        """)
+
+    # Create depletion operator with no reactions
+    micro_xs = openmc.deplete.MicroXS()
+    op = openmc.deplete.IndependentOperator.from_nuclides(
+        volume=1.0,
+        nuclides={'Li5': 1.0},
+        micro_xs=micro_xs,
+        chain_file='test_chain.xml',
+        normalization_mode='source-rate'
+    )
+
+    # Create time-integrator and integrate
+    integrator = openmc.deplete.PredictorIntegrator(
+        op, timesteps=[1.0], source_rates=[0.0], timestep_units='d'
+    )
+    integrator.integrate(final_step=False)
+
+    # Get concentration of H1 and He4
+    results = openmc.deplete.Results('depletion_results.h5')
+    _, h1 = results.get_atoms("1", "H1")
+    _, he4 = results.get_atoms("1", "He4")
+
+    # Since we started with 1e24 atoms of Li5, we should have 1e24 atoms of both
+    # H1 and He4
+    assert h1[1] == pytest.approx(1e24)
+    assert he4[1] == pytest.approx(1e24)