Merge pull request #2158 from yardasol/independent-operator-hotfix

`IndependentOperator` hotfix

docs/source/pythonapi/deplete.rst

@@ -43,7 +43,7 @@ algorithms <http://hdl.handle.net/1721.1/113721>`_.
     SILEQIIntegrator
 
 Each of these classes expects a "transport operator" to be passed. OpenMC
-provides The following classes implementing transpor operators:
+provides the following transport operator classes:
 
 .. autosummary::
    :toctree: generated
@@ -55,7 +55,7 @@ provides The following classes implementing transpor operators:
 
 The :class:`CoupledOperator` and :class:`IndependentOperator` classes must also
 have some knowledge of how nuclides transmute and decay. This is handled by the 
-:class:`Chain`.
+:class:`Chain` class.
 
 Minimal Example
 ---------------

docs/source/usersguide/depletion.rst

@@ -21,7 +21,8 @@ material compositions over time. Each method appears as a different class.
 For example, :class:`openmc.deplete.CECMIntegrator` runs a depletion calculation
 using the CE/CM algorithm (deplete over a timestep using the middle-of-step
 reaction rates). An instance of :class:`~openmc.deplete.abc.TransportOperator`
-is passed to one of these functions along with the timesteps and power level::
+is passed to one of these Integrator classes along with the timesteps and power
+level::
 
     power = 1200.0e6  # watts
     timesteps = [10.0, 10.0, 10.0]  # days
@@ -318,7 +319,7 @@ normalizing reaction rates:
 
    .. math::
 
-      \phi = \frac{P}{V \cdot \sum_i (Q_i \cdot \Sigma^f_i \cdot \rho_i)}
+      \phi = \frac{P}{V \cdot \sum_i (Q_i \cdot \sigma^f_i \cdot \n_i)}
 
    where :math:`\sum_i` is the sum over all nuclides :math:`i`. This equation
    makes the same assumptions and issues as discussed in

openmc/deplete/abc.py

@@ -331,7 +331,7 @@ class NormalizationHelper(ABC):
             `fission_rates` for :meth:`update`.
         """
 
-    def update(self, fission_rates, mat_index=None):
+    def update(self, fission_rates):
         """Update the normalization based on fission rates (only used for
         energy-based normalization)
 
@@ -341,8 +341,6 @@ class NormalizationHelper(ABC):
             fission reaction rate for each isotope in the specified
             material. Should be ordered corresponding to initial
             ``rate_index`` used in :meth:`prepare`
-        mat_index : int
-            Material index
         """
 
     @property

openmc/deplete/helpers.py

@@ -426,7 +426,7 @@ class ChainFissionHelper(EnergyNormalizationHelper):
 
         self._fission_q_vector = fission_qs
 
-    def update(self, fission_rates, mat_index=None):
+    def update(self, fission_rates):
         """Update energy produced with fission rates in a material
 
         Parameters
@@ -435,8 +435,6 @@ class ChainFissionHelper(EnergyNormalizationHelper):
             fission reaction rate for each isotope in the specified
             material. Should be ordered corresponding to initial
             ``rate_index`` used in :meth:`prepare`
-        mat_index : int
-            Unused
         """
         self._energy += dot(fission_rates, self._fission_q_vector)
 

openmc/deplete/independent_operator.py

@@ -55,6 +55,10 @@ 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``.
@@ -108,8 +112,9 @@ class IndependentOperator(OpenMCOperator):
                  micro_xs,
                  chain_file,
                  keff=None,
-                 normalization_mode='source-rate',
+                 normalization_mode='fission-q',
                  fission_q=None,
+                 dilute_initial=1.0e3,
                  prev_results=None,
                  reduce_chain=False,
                  reduce_chain_level=None,
@@ -135,6 +140,7 @@ 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)
@@ -145,8 +151,9 @@ class IndependentOperator(OpenMCOperator):
                       chain_file,
                       nuc_units='atom/b-cm',
                       keff=None,
-                      normalization_mode='source-rate',
+                      normalization_mode='fission-q',
                       fission_q=None,
+                      dilute_initial=1.0e3,
                       prev_results=None,
                       reduce_chain=False,
                       reduce_chain_level=None,
@@ -178,6 +185,10 @@ 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
@@ -202,6 +213,7 @@ 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,
@@ -253,47 +265,6 @@ class IndependentOperator(OpenMCOperator):
         """Finds nuclides with cross section data"""
         return set(cross_sections.index)
 
-    class _IndependentNormalizationHelper(ChainFissionHelper):
-        """Class for calculating one-group flux based on a power.
-
-        flux = Power / X, where X = volume * sum_i(Q_i * fission_micro_xs_i * density_i)
-
-        Parameters
-        ----------
-        op : openmc.deplete.IndependentOperator
-            Reference to the object encapsulating _IndependentNormalizationHelper.
-            We pass this so we don't have to duplicate :attr:`IndependentOperator.number`.
-
-        """
-
-        def __init__(self, op):
-            rates = op.reaction_rates
-            self.nuc_ind_map = {ind: nuc for nuc, ind in rates.index_nuc.items()}
-            self._op = op
-            super().__init__()
-
-        def update(self, fission_rates, mat_index=None):
-            """Update 'energy' produced with fission rates in a material. What
-            this actually calculates is the quantity X.
-
-            Parameters
-            ----------
-            fission_rates : numpy.ndarray
-                fission reaction rate for each isotope in the specified
-                material. Should be ordered corresponding to initial
-                ``rate_index`` used in :meth:`prepare`
-            mat_index : int
-                Material index
-
-            """
-            volume = self._op.number.get_mat_volume(mat_index)
-            densities = np.empty_like(fission_rates)
-            for i_nuc, nuc in self.nuc_ind_map.items():
-                densities[i_nuc] = self._op.number.get_atom_density(mat_index, nuc)
-            fission_rates *= volume * densities
-
-            super().update(fission_rates)
-
     class _IndependentRateHelper(ReactionRateHelper):
         """Class for generating one-group reaction rates with flux and
         one-group cross sections.
@@ -350,8 +321,10 @@ class IndependentOperator(OpenMCOperator):
                 nuc = self.nuc_ind_map[i_nuc]
                 rxn = self.rxn_ind_map[i_react]
                 density = self._op.number.get_atom_density(mat_id, nuc)
-                self._results_cache[i_nuc,
-                                    i_react] = self._op.cross_sections[rxn][nuc] * density * volume
+
+                # Sigma^j_i * V = sigma^j_i * rho * V
+                self._results_cache[i_nuc,i_react] = \
+                    self._op.cross_sections[rxn][nuc] * density * volume
 
             return self._results_cache
 
@@ -370,7 +343,7 @@ class IndependentOperator(OpenMCOperator):
 
         self._rate_helper = self._IndependentRateHelper(self)
         if normalization_mode == "fission-q":
-            self._normalization_helper = self._IndependentNormalizationHelper(self)
+            self._normalization_helper = ChainFissionHelper()
         else:
             self._normalization_helper = SourceRateHelper()
 

openmc/deplete/microxs.py

@@ -1,7 +1,7 @@
 """MicroXS module
 
 A pandas.DataFrame storing microscopic cross section data with
-nuclides names as row indices and reaction names as column indices.
+nuclide names as row indices and reaction names as column indices.
 """
 
 import tempfile
@@ -26,7 +26,7 @@ _valid_rxns.append('fission')
 
 class MicroXS(DataFrame):
     """Stores microscopic cross section data for use in
-    independent depletion.
+    transport-independent depletion.
     """
 
     @classmethod
@@ -35,7 +35,8 @@ class MicroXS(DataFrame):
                    reaction_domain,
                    chain_file,
                    dilute_initial=1.0e3,
-                   energy_bounds=(0, 20e6)):
+                   energy_bounds=(0, 20e6),
+                   run_kwargs=None):
         """Generate a one-group cross-section dataframe using
         OpenMC. Note that the ``openmc`` executable must be compiled.
 
@@ -57,6 +58,8 @@ class MicroXS(DataFrame):
             Reaction names to tally
         energy_bound : 2-tuple of float, optional
             Bounds for the energy group.
+        run_kwargs : dict, optional
+            Keyword arguments for :meth:`openmc.model.Model.run()`
 
         Returns
         -------
@@ -89,7 +92,10 @@ class MicroXS(DataFrame):
 
         # create temporary run
         with tempfile.TemporaryDirectory() as temp_dir:
-            statepoint_path = model.run(cwd=temp_dir)
+            if run_kwargs is None:
+                run_kwargs = {}
+            run_kwargs.setdefault('cwd', temp_dir)
+            statepoint_path = model.run(**run_kwargs)
 
             with StatePoint(statepoint_path) as sp:
                 for rx in xs:
@@ -146,7 +152,7 @@ class MicroXS(DataFrame):
         for material in model.materials:
             if material.depletable:
                 nuc_densities = material.get_nuclide_atom_densities()
-                dilute_density = 1.e-24 * dilute_initial
+                dilute_density = 1.0e-24 * dilute_initial
                 material.set_density('sum')
                 for nuc, density in nuc_densities.items():
                     material.remove_nuclide(nuc)

openmc/deplete/openmc_operator.py

@@ -528,8 +528,7 @@ class OpenMCOperator(TransportOperator):
             # Accumulate energy from fission
             if fission_ind is not None:
                 self._normalization_helper.update(
-                    tally_rates[:, fission_ind],
-                    mat_index=mat_index)
+                    tally_rates[:, fission_ind])
 
             # Divide by total number and store
             rates[i] = self._rate_helper.divide_by_adens(number)

openmc/deplete/results.py

@@ -417,7 +417,16 @@ class Results(list):
             mat_id = str(mat.id)
             if mat_id in result.mat_to_ind:
                 mat.volume = result.volume[mat_id]
+
+                # Change density of all nuclides in material to atom/b-cm
+                atoms_per_barn_cm = mat.get_nuclide_atom_densities()
+                for nuc, value in atoms_per_barn_cm.items():
+                    mat.remove_nuclide(nuc)
+                    mat.add_nuclide(nuc, value)
                 mat.set_density('sum')
+
+                # For nuclides in chain that have cross sections, replace
+                # density in original material with new density from results
                 for nuc in result.nuc_to_ind:
                     if nuc not in available_cross_sections:
                         continue

tests/micro_xs_simple.csv

@@ -1,13 +1,13 @@
 nuclide,"(n,gamma)",fission
-U234,22.231989822002465,0.49620744663749855
-U235,10.479008971197121,48.41787337164604
-U238,0.8673334105437324,0.10467880588762352
-U236,8.65171044607122,0.31948392400019293
-O16,7.497851000107524e-05,0.0
-O17,0.0004079227797153371,0.0
-I135,6.842395323713927,0.0
-Xe135,227463.86426990604,0.0
-Xe136,0.02317896034753588,0.0
+U234,22.231989822002454,0.4962074466374984
+U235,10.479008971197121,48.41787337164606
+U238,0.8673334105437321,0.1046788058876236
+U236,8.651710446071224,0.31948392400019293
+O16,7.497851000107522e-05,0.0
+O17,0.0004079227797153372,0.0
+I135,6.842395323713929,0.0
+Xe135,227463.8642699061,0.0
+Xe136,0.023178960347535887,0.0
 Cs135,2.1721665580713623,0.0
-Gd157,12786.09939237018,0.0
+Gd157,12786.099392370175,0.0
 Gd156,3.4006085445846983,0.0

tests/regression_tests/deplete_no_transport/test.py

@@ -20,6 +20,15 @@ def fuel():
 
     return fuel
 
+@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.xml'
+
 @pytest.mark.parametrize("multiproc, from_nuclides, normalization_mode, power, flux", [
     (True, True,'source-rate', None, 1164719970082145.0),
     (False, True, 'source-rate', None, 1164719970082145.0),
@@ -29,7 +38,15 @@ def fuel():
     (False, False, 'source-rate', None, 1164719970082145.0),
     (True, False, 'fission-q', 174, None),
     (False, False, 'fission-q', 174, None)])
-def test_no_transport_from_nuclides(run_in_tmpdir, fuel, multiproc, from_nuclides, normalization_mode, power, flux):
+def test_against_self(run_in_tmpdir,
+                      fuel,
+                      micro_xs,
+                      chain_file,
+                      multiproc,
+                      from_nuclides,
+                      normalization_mode,
+                      power,
+                      flux):
     """Transport free system test suite.
 
     Runs an OpenMC transport-free depletion calculation and verifies
@@ -37,28 +54,22 @@ def test_no_transport_from_nuclides(run_in_tmpdir, fuel, multiproc, from_nuclide
 
     """
     # Create operator
-    micro_xs_file = Path(__file__).parents[2] / 'micro_xs_simple.csv'
-    micro_xs = MicroXS.from_csv(micro_xs_file)
-    chain_file = Path(__file__).parents[2] / 'chain_simple.xml'
-
-    if from_nuclides:
-        nuclides = {}
-        for nuc, dens in fuel.get_nuclide_atom_densities().items():
-            nuclides[nuc] = dens
-
-        op = IndependentOperator.from_nuclides(
-            fuel.volume, nuclides, micro_xs, chain_file, normalization_mode=normalization_mode)
-
-    else:
-        op = IndependentOperator(openmc.Materials([fuel]), micro_xs, chain_file, normalization_mode=normalization_mode)
+    op = _create_operator(from_nuclides,
+                          fuel,
+                          micro_xs,
+                          chain_file,
+                          normalization_mode)
 
     # Power and timesteps
-    dt = [30]  # single step
+    dt = [360]  # single step
 
     # Perform simulation using the predictor algorithm
     openmc.deplete.pool.USE_MULTIPROCESSING = multiproc
-    openmc.deplete.PredictorIntegrator(
-        op, dt, power=power, source_rates=flux, timestep_units='d').integrate()
+    openmc.deplete.PredictorIntegrator(op,
+                                       dt,
+                                       power=power,
+                                       source_rates=flux,
+                                       timestep_units='s').integrate()
 
     # Get path to test and reference results
     path_test = op.output_dir / 'depletion_results.h5'
@@ -73,9 +84,86 @@ def test_no_transport_from_nuclides(run_in_tmpdir, fuel, multiproc, from_nuclide
     res_ref = openmc.deplete.Results(path_reference)
 
     # Assert same mats
+    _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)
+
+@pytest.mark.parametrize("multiproc, dt, time_units, time_type, atom_tol, rx_tol ", [
+    (True, 360, 's', 'minutes', 2.0e-3, 3.0e-2),
+    (False, 360, 's', 'minutes', 2.0e-3, 3.0e-2),
+    (True, 4, 'h', 'hours', 2.0e-3, 6.0e-2),
+    (False,4, 'h', 'hours', 2.0e-3, 6.0e-2),
+    (True, 5, 'd', 'days', 2.0e-3, 5.0e-2),
+    (False,5, 'd', 'days', 2.0e-3, 5.0e-2),
+    (True, 100, 'd', 'months', 4.0e-3, 9.0e-2),
+    (False, 100, 'd', 'months', 4.0e-3, 9.0e-2)])
+def test_against_coupled(run_in_tmpdir,
+                         fuel,
+                         micro_xs,
+                         chain_file,
+                         multiproc,
+                         dt,
+                         time_units,
+                         time_type,
+                         atom_tol,
+                         rx_tol):
+    # Create operator
+    op = _create_operator(False, fuel, micro_xs, chain_file, 'fission-q')
+
+    # Power and timesteps
+    dt = [dt]  # single step
+
+    # Perform simulation using the predictor algorithm
+    openmc.deplete.pool.USE_MULTIPROCESSING = multiproc
+    openmc.deplete.PredictorIntegrator(
+        op, dt, power=174, timestep_units=time_units).integrate()
+
+    # Get path to test and reference results
+    path_test = op.output_dir / 'depletion_results.h5'
+
+    ref_path = f'test_reference_coupled_{time_type}.h5'
+    path_reference = Path(__file__).with_name(ref_path)
+
+    # Load the reference/test results
+    res_test = openmc.deplete.Results(path_test)
+    res_ref = openmc.deplete.Results(path_reference)
+
+    # 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)
+
+def _create_operator(from_nuclides,
+                     fuel,
+                     micro_xs,
+                     chain_file,
+                     normalization_mode):
+    if from_nuclides:
+        nuclides = {}
+        for nuc, dens in fuel.get_nuclide_atom_densities().items():
+            nuclides[nuc] = dens
+
+        op = IndependentOperator.from_nuclides(fuel.volume,
+                                               nuclides,
+                                               micro_xs,
+                                               chain_file,
+                                               normalization_mode=normalization_mode)
+
+    else:
+        op = IndependentOperator(openmc.Materials([fuel]),
+                                 micro_xs,
+                                 chain_file,
+                                 normalization_mode=normalization_mode)
+
+    return op
+
+def _assert_same_mats(res_ref, res_test):
     for mat in res_ref[0].mat_to_ind:
-        assert mat in res_test[0].mat_to_ind, \
-            "Material {} not in new results.".format(mat)
+            assert mat in res_test[0].mat_to_ind, \
+                "Material {} not in new results.".format(mat)
     for nuc in res_ref[0].nuc_to_ind:
         assert nuc in res_test[0].nuc_to_ind, \
             "Nuclide {} not in new results.".format(nuc)
@@ -87,7 +175,7 @@ def test_no_transport_from_nuclides(run_in_tmpdir, fuel, multiproc, from_nuclide
         assert nuc in res_ref[0].nuc_to_ind, \
             "Nuclide {} not in old results.".format(nuc)
 
-    tol = 1.0e-6
+def _assert_atoms_equal(res_ref, res_test, tol):
     for mat in res_test[0].mat_to_ind:
         for nuc in res_test[0].nuc_to_ind:
             _, y_test = res_test.get_atoms(mat, nuc)
@@ -104,3 +192,26 @@ def test_no_transport_from_nuclides(run_in_tmpdir, fuel, multiproc, from_nuclide
 
             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)