From 71f3a29010f7b2e5993341a14b64d46f57b12dd0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 6 Jun 2023 10:49:39 -0500 Subject: [PATCH 1/4] Work in progress using multiply_density = False for direct reaction rate helper --- openmc/deplete/abc.py | 3 +-- openmc/deplete/helpers.py | 25 +++++++++++++------------ openmc/deplete/openmc_operator.py | 22 +++++++++++++--------- 3 files changed, 27 insertions(+), 23 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index d4cc1ba4d7f..802b3833f16 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -274,8 +274,7 @@ def divide_by_atoms(self, number): Parameters ---------- number : iterable of float - Number density [atoms/b-cm] of each nuclide tracked in the - calculation. + Number of each nuclide in [atom] tracked in the calculation. Returns ------- diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index ec0f6565025..dd189fca7c5 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -26,6 +26,7 @@ "ConstantFissionYieldHelper", "FissionYieldCutoffHelper", "AveragedFissionYieldHelper", "FluxCollapseHelper") + class TalliedFissionYieldHelper(FissionYieldHelper): """Abstract class for computing fission yields with tallies @@ -159,22 +160,23 @@ def nuclides(self, nuclides): def generate_tallies(self, materials, scores): """Produce one-group reaction rate tally - Uses the :mod:`openmc.lib` to generate a tally - of relevant reactions across all burnable materials. + Uses the :mod:`openmc.lib` to generate a tally of relevant reactions + across all burnable materials. Parameters ---------- - materials : iterable of :class:`openmc.Material` - Burnable materials in the problem. Used to - construct a :class:`openmc.MaterialFilter` + materials : iterable of :class:`openmc.lib.Material` + Burnable materials in the problem. Used to construct a + :class:`openmc.lib.MaterialFilter` scores : iterable of str - Reaction identifiers, e.g. ``"(n, fission)"``, - ``"(n, gamma)"``, needed for the reaction rate tally. + Reaction identifiers, e.g. ``"(n, fission)"``, ``"(n, gamma)"``, + needed for the reaction rate tally. """ self._rate_tally = Tally() self._rate_tally.writable = False self._rate_tally.scores = scores self._rate_tally.filters = [MaterialFilter(materials)] + self._rate_tally.multiply_density = False self._rate_tally_means_cache = None @property @@ -194,7 +196,7 @@ def reset_tally_means(self): """ self._rate_tally_means_cache = None - def get_material_rates(self, mat_id, nuc_index, react_index): + def get_material_rates(self, mat_id, nuc_index, rx_index): """Return an array of reaction rates for a material Parameters @@ -204,7 +206,7 @@ def get_material_rates(self, mat_id, nuc_index, react_index): nuc_index : iterable of int Index for each nuclide in :attr:`nuclides` in the desired reaction rate matrix - react_index : iterable of int + rx_index : iterable of int Index for each reaction scored in the tally Returns @@ -215,9 +217,8 @@ def get_material_rates(self, mat_id, nuc_index, react_index): """ self._results_cache.fill(0.0) full_tally_res = self.rate_tally_means[mat_id] - for i_tally, (i_nuc, i_react) in enumerate( - product(nuc_index, react_index)): - self._results_cache[i_nuc, i_react] = full_tally_res[i_tally] + for i_tally, (i_nuc, i_rx) in enumerate(product(nuc_index, rx_index)): + self._results_cache[i_nuc, i_rx] = full_tally_res[i_tally] return self._results_cache diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py index 8a13bff7fa0..2c6af805bc4 100644 --- a/openmc/deplete/openmc_operator.py +++ b/openmc/deplete/openmc_operator.py @@ -371,8 +371,8 @@ def initial_condition(self, materials): Parameters ---------- - materials : list of str - list of material IDs + materials : list of openmc.lib.Material + list of materials Returns ------- @@ -503,7 +503,7 @@ def _calculate_reaction_rates(self, source_rate): # Form fast map nuc_ind = [rates.index_nuc[nuc] for nuc in rxn_nuclides] - react_ind = [rates.index_rx[react] for react in self.chain.reactions] + rx_ind = [rates.index_rx[react] for react in self.chain.reactions] # Keep track of energy produced from all reactions in eV per source # particle @@ -534,21 +534,25 @@ def _calculate_reaction_rates(self, source_rate): for nuc, i_nuc_results in zip(rxn_nuclides, nuc_ind): number[i_nuc_results] = self.number[mat, nuc] + # Get microscopic reaction rates in [(reactions/src)*b-cm/atom]. 2D + # array with shape (nuclides, reactions). tally_rates = self._rate_helper.get_material_rates( - mat_index, nuc_ind, react_ind) + mat_index, nuc_ind, rx_ind) # Compute fission yields for this material fission_yields.append(self._yield_helper.weighted_yields(i)) # Accumulate energy from fission + barn_cm = 1e24 * self.number.get_mat_volume(mat) if fission_ind is not None: - self._normalization_helper.update( - tally_rates[:, fission_ind]) + atom_per_bcm = number / barn_cm + fission_rates = tally_rates[:, fission_ind] * atom_per_bcm + self._normalization_helper.update(fission_rates) - # Divide by total number of atoms and store - rates[i] = self._rate_helper.divide_by_atoms(number) + # Divide by [b-cm] to get [(reactions/src)/atom] + rates[i] = tally_rates / barn_cm - # Scale reaction rates to obtain units of (reactions/sec)/atom + # Scale reaction rates to obtain units of [(reactions/sec)/atom] rates *= self._normalization_helper.factor(source_rate) # Store new fission yields on the chain From 7735a33f35b543f54521e5755993513b550688f9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 7 Jun 2023 06:48:14 -0500 Subject: [PATCH 2/4] Make sure FluxCollapseHelper returns reaction rates not multiply by atom/b-cm --- openmc/deplete/helpers.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index dd189fca7c5..a224bfd6e6e 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -305,6 +305,7 @@ def generate_tallies(self, materials, scores): self._rate_tally.writable = False self._rate_tally.scores = self._reactions_direct self._rate_tally.filters = [MaterialFilter(materials)] + self._rate_tally.multiply_density = False self._rate_tally_means_cache = None if self._nuclides_direct is not None: # check if any direct tally nuclides are requested that are not @@ -376,13 +377,7 @@ def get_material_rates(self, mat_index, nuc_index, react_index): mat = self._materials[mat_index] - # Build nucname: density mapping to enable O(1) lookup in loop below - densities = dict(zip(mat.nuclides, mat.densities)) - for name, i_nuc in zip(self.nuclides, nuc_index): - # Determine density of nuclide - density = densities[name] - for mt, score, i_rx in zip(self._mts, self._scores, react_index): if score in self._reactions_direct and name in nuclides_direct: # Determine index in rx_rates @@ -398,7 +393,7 @@ def get_material_rates(self, mat_index, nuc_index, react_index): mt, mat.temperature, self._energies, flux) # Multiply by density to get absolute reaction rate - self._results_cache[i_nuc, i_rx] = rate_per_nuc * density + self._results_cache[i_nuc, i_rx] = rate_per_nuc return self._results_cache From e9306a852db8c57e154a8d48b24a6d4f4c945b53 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 7 Jun 2023 10:37:40 -0500 Subject: [PATCH 3/4] Only multiply by volume in _IndependentRateHelper.get_material_rates --- openmc/deplete/independent_operator.py | 22 ++++++++++++++-------- openmc/deplete/openmc_operator.py | 6 +++--- 2 files changed, 17 insertions(+), 11 deletions(-) diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index 0c950b31b6f..6caef263521 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -295,7 +295,7 @@ class _IndependentRateHelper(ReactionRateHelper): ---------- nuc_ind_map : dict of int to str Dictionary mapping the nuclide index to nuclide name - rxn_ind_map : dict of int to str + rx_ind_map : dict of int to str Dictionary mapping reaction index to reaction name """ @@ -305,7 +305,7 @@ def __init__(self, op): super().__init__(rates.n_nuc, rates.n_react) self.nuc_ind_map = {ind: nuc for nuc, ind in rates.index_nuc.items()} - self.rxn_ind_map = {ind: rxn for rxn, ind in rates.index_rx.items()} + self.rx_ind_map = {ind: rxn for rxn, ind in rates.index_rx.items()} self._op = op def generate_tallies(self, materials, scores): @@ -330,14 +330,20 @@ def get_material_rates(self, mat_id, nuc_index, react_index): """ self._results_cache.fill(0.0) + # Get volume in units of [b-cm] + volume_b_cm = 1e24 * self._op.number.get_mat_volume(mat_id) + for i_nuc, i_react in product(nuc_index, react_index): nuc = self.nuc_ind_map[i_nuc] - rxn = self.rxn_ind_map[i_react] - - # Sigma^j_i * V = sigma^j_i * n_i * V = sigma^j_i * N_i - self._results_cache[i_nuc,i_react] = \ - self._op.cross_sections[rxn][nuc] * \ - self._op.number[mat_id, nuc] + rx = self.rx_ind_map[i_react] + + # OK, this is kind of weird, but we multiply by volume in [b-cm] + # only because OpenMCOperator._calculate_reaction_rates has to + # divide it out later. It might make more sense to account for + # the source rate (flux) here rather than in the normalization + # helper. + self._results_cache[i_nuc, i_react] = \ + self._op.cross_sections[rx][nuc] * volume_b_cm return self._results_cache diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py index 2c6af805bc4..adf1c2eb7b3 100644 --- a/openmc/deplete/openmc_operator.py +++ b/openmc/deplete/openmc_operator.py @@ -543,14 +543,14 @@ def _calculate_reaction_rates(self, source_rate): fission_yields.append(self._yield_helper.weighted_yields(i)) # Accumulate energy from fission - barn_cm = 1e24 * self.number.get_mat_volume(mat) + volume_b_cm = 1e24 * self.number.get_mat_volume(mat) if fission_ind is not None: - atom_per_bcm = number / barn_cm + atom_per_bcm = number / volume_b_cm fission_rates = tally_rates[:, fission_ind] * atom_per_bcm self._normalization_helper.update(fission_rates) # Divide by [b-cm] to get [(reactions/src)/atom] - rates[i] = tally_rates / barn_cm + rates[i] = tally_rates / volume_b_cm # Scale reaction rates to obtain units of [(reactions/sec)/atom] rates *= self._normalization_helper.factor(source_rate) From dac50d485a03b316b9836d4e946dfd395b8b2dc5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 12 Jun 2023 23:05:46 -0500 Subject: [PATCH 4/4] Update openmc/deplete/helpers.py Co-authored-by: Andrew Johnson --- openmc/deplete/helpers.py | 1 - 1 file changed, 1 deletion(-) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index a224bfd6e6e..9a1fc6a3fa4 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -392,7 +392,6 @@ def get_material_rates(self, mat_index, nuc_index, react_index): rate_per_nuc = nuc.collapse_rate( mt, mat.temperature, self._energies, flux) - # Multiply by density to get absolute reaction rate self._results_cache[i_nuc, i_rx] = rate_per_nuc return self._results_cache