Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use multiply_density in depletion reaction rate helper classes #2559

Merged
merged 4 commits into from
Jun 15, 2023
+43 −39
Merged

Use multiply_density in depletion reaction rate helper classes #2559

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
71f3a29
Work in progress using multiply_density = False for direct reaction r…
paulromano Jun 6, 2023
7735a33
Make sure FluxCollapseHelper returns reaction rates not multiply by a…
paulromano Jun 7, 2023
e9306a8
Only multiply by volume in _IndependentRateHelper.get_material_rates
paulromano Jun 7, 2023
dac50d4
Update openmc/deplete/helpers.py
paulromano Jun 13, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 1 addition & 2 deletions openmc/deplete/abc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
-------
Expand Down
35 changes: 15 additions & 20 deletions openmc/deplete/helpers.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
"ConstantFissionYieldHelper", "FissionYieldCutoffHelper",
"AveragedFissionYieldHelper", "FluxCollapseHelper")


class TalliedFissionYieldHelper(FissionYieldHelper):
"""Abstract class for computing fission yields with tallies

Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -304,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
Expand Down Expand Up @@ -375,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
Expand All @@ -396,8 +392,7 @@ 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 * density
self._results_cache[i_nuc, i_rx] = rate_per_nuc

return self._results_cache

Expand Down
22 changes: 14 additions & 8 deletions openmc/deplete/independent_operator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

"""
Expand All @@ -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):
Expand All @@ -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.
paulromano marked this conversation as resolved.
Show resolved Hide resolved
self._results_cache[i_nuc, i_react] = \
self._op.cross_sections[rx][nuc] * volume_b_cm

return self._results_cache

Expand Down
22 changes: 13 additions & 9 deletions openmc/deplete/openmc_operator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
-------
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
volume_b_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 / volume_b_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 / volume_b_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
Expand Down