Remove initial dilute nuclides in depletion (#2568)

docs/source/io_formats/depletion_results.rst

@@ -47,10 +47,3 @@ The current version of the depletion results file format is 1.1.
 **/reactions/<name>/**
 
 :Attributes: - **index** (*int*) -- Index user in results for this reaction
-
-.. note::
-
-    The reaction rates for some isotopes not originally present may
-    be non-zero, but should be negligible compared to other atoms.
-    This can be controlled by changing the
-    :class:`openmc.deplete.Operator` ``dilute_initial`` attribute.

openmc/deplete/abc.py

@@ -93,20 +93,11 @@ class TransportOperator(ABC):
     fission_q : dict, optional
         Dictionary of nuclides and their fission Q values [eV]. If not given,
         values will be pulled from the ``chain_file``.
-    dilute_initial : float, optional
-        Initial atom density [atoms/cm^3] to add for nuclides that are zero
-        in initial condition to ensure they exist in the decay chain.
-        Only done for nuclides with reaction rates.
-        Defaults to 1.0e3.
     prev_results : Results, optional
         Results from a previous depletion calculation.
 
     Attributes
     ----------
-    dilute_initial : float
-        Initial atom density [atoms/cm^3] to add for nuclides that are zero
-        in initial condition to ensure they exist in the decay chain.
-        Only done for nuclides with reaction rates.
     output_dir : pathlib.Path
         Path to output directory to save results.
     prev_res : Results or None
@@ -116,9 +107,7 @@ class TransportOperator(ABC):
         The depletion chain information necessary to form matrices and tallies.
 
     """
-    def __init__(self, chain_file, fission_q=None, dilute_initial=1.0e3,
-                 prev_results=None):
-        self.dilute_initial = dilute_initial
+    def __init__(self, chain_file, fission_q=None, prev_results=None):
         self.output_dir = '.'
 
         # Read depletion chain
@@ -129,17 +118,6 @@ class TransportOperator(ABC):
             check_type("previous results", prev_results, Results)
             self.prev_res = prev_results
 
-    @property
-    def dilute_initial(self):
-        """Initial atom density for nuclides with zero initial concentration"""
-        return self._dilute_initial
-
-    @dilute_initial.setter
-    def dilute_initial(self, value):
-        check_type("dilute_initial", value, Real)
-        check_greater_than("dilute_initial", value, 0.0, equality=True)
-        self._dilute_initial = value
-
     @abstractmethod
     def __call__(self, vec, source_rate):
         """Runs a simulation.

openmc/deplete/coupled_operator.py

@@ -122,10 +122,6 @@ class CoupledOperator(OpenMCOperator):
         Dictionary of nuclides and their fission Q values [eV]. If not given,
         values will be pulled from the ``chain_file``. Only applicable
         if ``"normalization_mode" == "fission-q"``
-    dilute_initial : float, optional
-        Initial atom density [atoms/cm^3] to add for nuclides that are zero
-        in initial condition to ensure they exist in the decay chain.
-        Only done for nuclides with reaction rates.
     fission_yield_mode : {"constant", "cutoff", "average"}
         Key indicating what fission product yield scheme to use. The
         key determines what fission energy helper is used:
@@ -177,10 +173,6 @@ class CoupledOperator(OpenMCOperator):
         OpenMC geometry object
     settings : openmc.Settings
         OpenMC settings object
-    dilute_initial : float
-        Initial atom density [atoms/cm^3] to add for nuclides that
-        are zero in initial condition to ensure they exist in the decay
-        chain. Only done for nuclides with reaction rates.
     output_dir : pathlib.Path
         Path to output directory to save results.
     round_number : bool
@@ -215,7 +207,7 @@ class CoupledOperator(OpenMCOperator):
 
     def __init__(self, model, chain_file=None, prev_results=None,
                  diff_burnable_mats=False, normalization_mode="fission-q",
-                 fission_q=None, dilute_initial=1.0e3,
+                 fission_q=None,
                  fission_yield_mode="constant", fission_yield_opts=None,
                  reaction_rate_mode="direct", reaction_rate_opts=None,
                  reduce_chain=False, reduce_chain_level=None):
@@ -271,7 +263,6 @@ class CoupledOperator(OpenMCOperator):
             prev_results,
             diff_burnable_mats,
             fission_q,
-            dilute_initial,
             helper_kwargs,
             reduce_chain,
             reduce_chain_level)

openmc/deplete/independent_operator.py

@@ -63,10 +63,6 @@ class IndependentOperator(OpenMCOperator):
         Dictionary of nuclides and their fission Q values [eV]. If not given,
         values will be pulled from the ``chain_file``. Only applicable
         if ``"normalization_mode" == "fission-q"``.
-    dilute_initial : float, optional
-        Initial atom density [atoms/cm^3] to add for nuclides that are zero
-        in initial condition to ensure they exist in the decay chain.
-        Only done for nuclides with reaction rates.
     reduce_chain : bool, optional
         If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the
         depletion chain up to ``reduce_chain_level``.
@@ -86,10 +82,6 @@ class IndependentOperator(OpenMCOperator):
         All materials present in the model
     cross_sections : MicroXS
         Object containing one-group cross-sections in [cm^2].
-    dilute_initial : float
-        Initial atom density [atoms/cm^3] to add for nuclides that
-        are zero in initial condition to ensure they exist in the decay
-        chain. Only done for nuclides with reaction rates.
     output_dir : pathlib.Path
         Path to output directory to save results.
     round_number : bool
@@ -122,7 +114,6 @@ class IndependentOperator(OpenMCOperator):
                  keff=None,
                  normalization_mode='fission-q',
                  fission_q=None,
-                 dilute_initial=1.0e3,
                  prev_results=None,
                  reduce_chain=False,
                  reduce_chain_level=None,
@@ -148,7 +139,6 @@ class IndependentOperator(OpenMCOperator):
             chain_file,
             prev_results,
             fission_q=fission_q,
-            dilute_initial=dilute_initial,
             helper_kwargs=helper_kwargs,
             reduce_chain=reduce_chain,
             reduce_chain_level=reduce_chain_level)
@@ -161,7 +151,6 @@ class IndependentOperator(OpenMCOperator):
                       keff=None,
                       normalization_mode='fission-q',
                       fission_q=None,
-                      dilute_initial=1.0e3,
                       prev_results=None,
                       reduce_chain=False,
                       reduce_chain_level=None,
@@ -196,10 +185,6 @@ class IndependentOperator(OpenMCOperator):
             Dictionary of nuclides and their fission Q values [eV]. If not
             given, values will be pulled from the ``chain_file``. Only
             applicable if ``"normalization_mode" == "fission-q"``.
-        dilute_initial : float
-            Initial atom density [atoms/cm^3] to add for nuclides that
-            are zero in initial condition to ensure they exist in the decay
-            chain. Only done for nuclides with reaction rates.
         prev_results : Results, optional
             Results from a previous depletion calculation.
         reduce_chain : bool, optional
@@ -224,7 +209,6 @@ class IndependentOperator(OpenMCOperator):
                    keff=keff,
                    normalization_mode=normalization_mode,
                    fission_q=fission_q,
-                   dilute_initial=dilute_initial,
                    prev_results=prev_results,
                    reduce_chain=reduce_chain,
                    reduce_chain_level=reduce_chain_level,

openmc/deplete/openmc_operator.py

@@ -47,10 +47,6 @@ class OpenMCOperator(TransportOperator):
         Volumes are divided equally from the original material volume.
     fission_q : dict, optional
         Dictionary of nuclides and their fission Q values [eV].
-    dilute_initial : float, optional
-        Initial atom density [atoms/cm^3] to add for nuclides that are zero
-        in initial condition to ensure they exist in the decay chain.
-        Only done for nuclides with reaction rates.
     helper_kwargs : dict
         Keyword arguments for helper classes
     reduce_chain : bool, optional
@@ -68,10 +64,6 @@ class OpenMCOperator(TransportOperator):
     cross_sections : str or MicroXS
             Path to continuous energy cross section library, or object
             containing one-group cross-sections.
-    dilute_initial : float
-        Initial atom density [atoms/cm^3] to add for nuclides that
-        are zero in initial condition to ensure they exist in the decay
-        chain. Only done for nuclides with reaction rates.
     output_dir : pathlib.Path
         Path to output directory to save results.
     round_number : bool
@@ -105,7 +97,6 @@ class OpenMCOperator(TransportOperator):
             prev_results=None,
             diff_burnable_mats=False,
             fission_q=None,
-            dilute_initial=0.0,
             helper_kwargs=None,
             reduce_chain=False,
             reduce_chain_level=None):
@@ -119,7 +110,7 @@ class OpenMCOperator(TransportOperator):
                     "chain in openmc.config['chain_file']"
                 )
 
-        super().__init__(chain_file, fission_q, dilute_initial, prev_results)
+        super().__init__(chain_file, fission_q, prev_results)
         self.round_number = False
         self.materials = materials
         self.cross_sections = cross_sections
@@ -265,11 +256,6 @@ class OpenMCOperator(TransportOperator):
         """
         self.number = AtomNumber(local_mats, all_nuclides, volume, len(self.chain))
 
-        if self.dilute_initial != 0.0:
-            for nuc in self._burnable_nucs:
-                self.number.set_atom_density(
-                    np.s_[:], nuc, self.dilute_initial)
-
         # Now extract and store the number densities
         # From the geometry if no previous depletion results
         if prev_res is None:
@@ -433,8 +419,7 @@ class OpenMCOperator(TransportOperator):
         # for burning in which the number density is greater than zero.
         for nuc in self.number.nuclides:
             if nuc in self.nuclides_with_data:
-                if np.sum(self.number[:, nuc]) > 0.0:
-                    nuc_set.add(nuc)
+                nuc_set.add(nuc)
 
         # Communicate which nuclides have nonzeros to rank 0
         if comm.rank == 0:

openmc/deplete/results.py

@@ -103,13 +103,6 @@ class Results(list):
     ) -> Tuple[np.ndarray, np.ndarray]:
         """Get number of nuclides over time from a single material
 
-        .. note::
-            Initial values for some isotopes that do not appear in
-            initial concentrations may be non-zero, depending on the
-            value of the :attr:`openmc.deplete.CoupledOperator.dilute_initial`
-            attribute. The :class:`openmc.deplete.CoupledOperator` class adds
-            isotopes according to this setting, which can be set to zero.
-
         Parameters
         ----------
         mat : openmc.Material, str
@@ -172,14 +165,6 @@ class Results(list):
     ) -> Tuple[np.ndarray, np.ndarray]:
         """Get reaction rate in a single material/nuclide over time
 
-        .. note::
-
-            Initial values for some isotopes that do not appear in
-            initial concentrations may be non-zero, depending on the
-            value of :class:`openmc.deplete.CoupledOperator` ``dilute_initial``
-            The :class:`openmc.deplete.CoupledOperator` adds isotopes according
-            to this setting, which can be set to zero.
-
         Parameters
         ----------
         mat : openmc.Material, str

src/tallies/tally.cpp

@@ -1258,10 +1258,13 @@ extern "C" int openmc_tally_set_nuclides(
     } else {
       auto search = data::nuclide_map.find(word);
       if (search == data::nuclide_map.end()) {
-        set_errmsg("Nuclide \"" + word + "\" has not been loaded yet");
-        return OPENMC_E_DATA;
+        int err = openmc_load_nuclide(word.c_str(), nullptr, 0);
+        if (err < 0) {
+          set_errmsg(openmc_err_msg);
+          return OPENMC_E_DATA;
+        }
       }
-      nucs.push_back(search->second);
+      nucs.push_back(data::nuclide_map.at(word));
     }
   }
 

tests/regression_tests/__init__.py

@@ -1,3 +1,5 @@
+import pytest
+
 # Test configuration options for regression tests
 config = {
     'event' : False,
@@ -8,3 +10,25 @@ config = {
     'update': False,
     'build_inputs': False
 }
+
+
+def assert_atoms_equal(res_ref, res_test, tol=1e-5):
+    for mat in res_test[0].index_mat:
+        for nuc in res_test[0].index_nuc:
+            _, y_test = res_test.get_atoms(mat, nuc)
+            _, y_ref = res_ref.get_atoms(mat, nuc)
+            assert y_test == pytest.approx(y_ref, rel=tol), \
+                f'Atoms not equal for material {mat}, nuclide {nuc}\n' \
+                f'y_ref={y_ref}\ny_test={y_test}'
+
+
+def assert_reaction_rates_equal(res_ref, res_test, tol=1e-5):
+    for reactions in res_test[0].rates:
+        for mat in reactions.index_mat:
+            for nuc in reactions.index_nuc:
+                for rx in reactions.index_rx:
+                    y_test = res_test.get_reaction_rate(mat, nuc, rx)[1]
+                    y_ref = res_ref.get_reaction_rate(mat, nuc, rx)[1]
+                    assert y_test == pytest.approx(y_ref, rel=tol), \
+                        f'Reaction rate not equal for material {mat}, nuclide '\
+                        f'{nuc}, {rx}\ny_ref={y_ref}\ny_test={y_test}'

tests/regression_tests/deplete_no_transport/test.py

@@ -1,6 +1,7 @@
 """ Transport-free depletion test suite """
 
 from pathlib import Path
+import shutil
 
 import numpy as np
 import pytest
@@ -8,6 +9,9 @@ import openmc
 import openmc.deplete
 from openmc.deplete import IndependentOperator, MicroXS
 
+from tests.regression_tests import config, assert_atoms_equal, \
+    assert_reaction_rates_equal
+
 
 @pytest.fixture(scope="module")
 def fuel():
@@ -82,6 +86,11 @@ def test_against_self(run_in_tmpdir,
         ref_path = 'test_reference_fission_q.h5'
     path_reference = Path(__file__).with_name(ref_path)
 
+    # If updating results, do so and return
+    if config['update']:
+        shutil.copyfile(str(path_test), str(path_reference))
+        return
+
     # Load the reference/test results
     res_test = openmc.deplete.Results(path_test)
     res_ref = openmc.deplete.Results(path_reference)
@@ -90,8 +99,8 @@ def test_against_self(run_in_tmpdir,
     _assert_same_mats(res_test, res_ref)
 
     tol = 1.0e-14
-    _assert_atoms_equal(res_test, res_ref, tol)
-    _assert_reaction_rates_equal(res_test, res_ref, tol)
+    assert_atoms_equal(res_test, res_ref, tol)
+    assert_reaction_rates_equal(res_test, res_ref, tol)
 
 
 @pytest.mark.parametrize("multiproc, dt, time_units, time_type, atom_tol, rx_tol ", [
@@ -130,6 +139,11 @@ def test_against_coupled(run_in_tmpdir,
     ref_path = f'test_reference_coupled_{time_type}.h5'
     path_reference = Path(__file__).with_name(ref_path)
 
+    # If updating results, do so and return
+    if config['update']:
+        shutil.copyfile(str(path_test), str(path_reference))
+        return
+
     # Load the reference/test results
     res_test = openmc.deplete.Results(path_test)
     res_ref = openmc.deplete.Results(path_reference)
@@ -137,8 +151,8 @@ def test_against_coupled(run_in_tmpdir,
     # Assert same mats
     _assert_same_mats(res_test, res_ref)
 
-    _assert_atoms_equal(res_test, res_ref, atom_tol)
-    _assert_reaction_rates_equal(res_test, res_ref, rx_tol)
+    assert_atoms_equal(res_test, res_ref, atom_tol)
+    assert_reaction_rates_equal(res_test, res_ref, rx_tol)
 
 
 def _create_operator(from_nuclides,
@@ -181,45 +195,3 @@ def _assert_same_mats(res_ref, res_test):
         assert nuc in res_ref[0].index_nuc, \
             f"Nuclide {nuc} not in old results."
 
-
-def _assert_atoms_equal(res_ref, res_test, tol):
-    for mat in res_test[0].index_mat:
-        for nuc in res_test[0].index_nuc:
-            _, y_test = res_test.get_atoms(mat, nuc)
-            _, y_old = res_ref.get_atoms(mat, nuc)
-
-            # Test each point
-            correct = True
-            for i, ref in enumerate(y_old):
-                if ref != y_test[i]:
-                    if ref != 0.0:
-                        correct = np.abs(y_test[i] - ref) / ref <= tol
-                    else:
-                        correct = False
-
-            assert correct, "Discrepancy in mat {} and nuc {}\n{}\n{}".format(
-                mat, nuc, y_old, y_test)
-
-
-def _assert_reaction_rates_equal(res_ref, res_test, tol):
-    for reactions in res_test[0].rates:
-        for mat in reactions.index_mat:
-            for nuc in reactions.index_nuc:
-                for rx in reactions.index_rx:
-                    y_test = res_test.get_reaction_rate(mat, nuc, rx)[1] / \
-                        res_test.get_atoms(mat, nuc)[1]
-                    y_old = res_ref.get_reaction_rate(mat, nuc, rx)[1] / \
-                        res_ref.get_atoms(mat, nuc)[1]
-
-                    # Test each point
-                    correct = True
-                    for i, ref in enumerate(y_old):
-                        if ref != y_test[i]:
-                            if ref != 0.0:
-                                correct = np.abs(y_test[i] - ref) / ref <= tol
-                            else:
-                                if y_test[i] != 0.0:
-                                    correct = False
-
-                    assert correct, "Discrepancy in mat {}, nuc {}, and rx {}\n{}\n{}".format(
-                        mat, nuc, rx, y_old, y_test)

tests/regression_tests/deplete_with_transfer_rates/test.py

@@ -1,6 +1,7 @@
 """ TransferRates depletion test suite """
 
 from pathlib import Path
+import shutil
 
 import numpy as np
 import pytest
@@ -8,6 +9,10 @@ import openmc
 import openmc.deplete
 from openmc.deplete import CoupledOperator
 
+from tests.regression_tests import config, assert_reaction_rates_equal, \
+    assert_atoms_equal
+
+
 @pytest.fixture
 def model():
     f = openmc.Material(name="f")
@@ -34,9 +39,9 @@ def model():
     geometry = openmc.Geometry([cell_f, cell_w])
 
     settings = openmc.Settings()
-    settings.particles = 100
+    settings.particles = 500
     settings.inactive = 0
-    settings.batches = 10
+    settings.batches = 2
 
     return openmc.Model(geometry, materials, settings)
 
@@ -66,15 +71,14 @@ def test_transfer_rates(run_in_tmpdir, model, rate, dest_mat, power, ref_result)
     path_test = op.output_dir / 'depletion_results.h5'
     path_reference = Path(__file__).with_name(f'ref_{ref_result}.h5')
 
+    # If updating results, do so and return
+    if config['update']:
+        shutil.copyfile(str(path_test), str(path_reference))
+        return
+
     # Load the reference/test results
     res_ref = openmc.deplete.Results(path_reference)
     res_test = openmc.deplete.Results(path_test)
 
-    for mat in res_test[0].rates[0].index_mat:
-        for nuc in res_test[0].rates[0].index_nuc:
-            for rx in res_test[0].rates[0].index_rx:
-                y_test = res_test.get_reaction_rate(mat, nuc, rx)[1] / \
-                    res_test.get_atoms(mat, nuc)[1]
-                y_ref = res_ref.get_reaction_rate(mat, nuc, rx)[1] / \
-                    res_ref.get_atoms(mat, nuc)[1]
-                assert y_test == pytest.approx(y_ref, abs=1e-5)
+    assert_atoms_equal(res_ref, res_test)
+    assert_reaction_rates_equal(res_ref, res_test)

tests/regression_tests/deplete_with_transport/test.py

@@ -12,7 +12,7 @@ import openmc
 from openmc.data import JOULE_PER_EV
 import openmc.deplete
 
-from tests.regression_tests import config
+from tests.regression_tests import config, assert_atoms_equal
 from .example_geometry import generate_problem
 
 
@@ -96,23 +96,7 @@ def test_full(run_in_tmpdir, problem, multiproc):
         assert nuc in res_ref[0].index_nuc, \
             "Nuclide {} not in old results.".format(nuc)
 
-    tol = 1.0e-6
-    for mat in res_test[0].index_mat:
-        for nuc in res_test[0].index_nuc:
-            _, y_test = res_test.get_atoms(mat, nuc)
-            _, y_old = res_ref.get_atoms(mat, nuc)
-
-            # Test each point
-            correct = True
-            for i, ref in enumerate(y_old):
-                if ref != y_test[i]:
-                    if ref != 0.0:
-                        correct = np.abs(y_test[i] - ref) / ref <= tol
-                    else:
-                        correct = False
-
-            assert correct, "Discrepancy in mat {} and nuc {}\n{}\n{}".format(
-                mat, nuc, y_old, y_test)
+    assert_atoms_equal(res_ref, res_test, tol=1e-6)
 
     # Compare statepoint files with depletion results