From 76fc6fb3e756c7da237e433f22289594b2b01055 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 4 Aug 2020 14:35:12 -0500 Subject: [PATCH] Implement FluxCollapseHelper that is used by Operator --- openmc/deplete/abc.py | 14 ++--- openmc/deplete/helpers.py | 116 ++++++++++++++++++++++++++++++++++++- openmc/deplete/operator.py | 35 ++++++++++- 3 files changed, 150 insertions(+), 15 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 4df2c09904..2615caaa18 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -225,11 +225,11 @@ class TransportOperator(ABC): class ReactionRateHelper(ABC): """Abstract class for generating reaction rates for operators - Responsible for generating reaction rate tallies for burnable - materials, given nuclides and scores from the operator. + Responsible for generating reaction rate tallies for burnable materials, + given nuclides and scores from the operator. - Reaction rates are passed back to the operator for be used in - an :class:`openmc.deplete.OperatorResult` instance + Reaction rates are passed back to the operator for be used in an + :class:`openmc.deplete.OperatorResult` instance. Parameters ---------- @@ -246,7 +246,6 @@ class ReactionRateHelper(ABC): def __init__(self, n_nucs, n_react): self._nuclides = None - self._rate_tally = None self._results_cache = empty((n_nucs, n_react)) @abstractmethod @@ -262,7 +261,6 @@ class ReactionRateHelper(ABC): def nuclides(self, nuclides): check_type("nuclides", nuclides, list, str) self._nuclides = nuclides - self._rate_tally.nuclides = nuclides @abstractmethod def get_material_rates(self, mat_id, nuc_index, react_index): @@ -281,8 +279,7 @@ class ReactionRateHelper(ABC): def divide_by_adens(self, number): """Normalize reaction rates by number of nuclides - Acts on the current material examined by - :meth:`get_material_rates` + Acts on the current material examined by :meth:`get_material_rates` Parameters ---------- @@ -559,7 +556,6 @@ class TalliedFissionYieldHelper(FissionYieldHelper): self._fission_rate_tally = Tally() self._fission_rate_tally.writable = False self._fission_rate_tally.scores = ['fission'] - self._fission_rate_tally.filters = [MaterialFilter(materials)] def update_tally_nuclides(self, nuclides): diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index edcf6a5b3f..8dca85afb0 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -7,13 +7,14 @@ from numbers import Real import bisect from collections import defaultdict -from numpy import dot, zeros, newaxis +from numpy import dot, zeros, newaxis, asarray from . import comm from openmc.checkvalue import check_type, check_greater_than -from openmc.data import JOULE_PER_EV +from openmc.data import JOULE_PER_EV, REACTION_NAME from openmc.lib import ( Tally, MaterialFilter, EnergyFilter, EnergyFunctionFilter) +import openmc.lib from .abc import ( ReactionRateHelper, NormalizationHelper, FissionYieldHelper, TalliedFissionYieldHelper) @@ -21,7 +22,7 @@ from .abc import ( __all__ = ( "DirectReactionRateHelper", "ChainFissionHelper", "EnergyScoreHelper" "SourceRateHelper", "ConstantFissionYieldHelper", "FissionYieldCutoffHelper", - "AveragedFissionYieldHelper") + "AveragedFissionYieldHelper", "FluxCollapseHelper") # ------------------------------------- # Helpers for generating reaction rates @@ -43,6 +44,14 @@ class DirectReactionRateHelper(ReactionRateHelper): nuclides : list of str All nuclides with desired reaction rates. """ + def __init__(self, n_nuc, n_react): + super().__init__(n_nuc, n_react) + self._rate_tally = None + + @ReactionRateHelper.nuclides.setter + def nuclides(self, nuclides): + ReactionRateHelper.nuclides.fset(self, nuclides) + self._rate_tally.nuclides = nuclides def generate_tallies(self, materials, scores): """Produce one-group reaction rate tally @@ -92,6 +101,107 @@ class DirectReactionRateHelper(ReactionRateHelper): return self._results_cache +class FluxCollapseHelper(ReactionRateHelper): + """Class that generates tallies for one-group rates + + .. versionadded:: 0.12.1 + + Parameters + ---------- + n_nucs : int + Number of burnable nuclides tracked by :class:`openmc.deplete.Operator` + n_react : int + Number of reactions tracked by :class:`openmc.deplete.Operator` + energies : iterable of float + Energy group boundaries for flux spectrum in [eV] + + Attributes + ---------- + nuclides : list of str + All nuclides with desired reaction rates. + """ + def __init__(self, n_nucs, n_reacts, energies): + super().__init__(n_nucs, n_reacts) + self._energies = asarray(energies) + + def generate_tallies(self, materials, scores): + """Produce multigroup flux spectrum tally + + Uses the :mod:`openmc.lib` module to generate a multigroup flux tally + for each burnable material. + + Parameters + ---------- + materials : iterable of :class:`openmc.Material` + Burnable materials in the problem. Used to construct a + :class:`openmc.MaterialFilter` + scores : iterable of str + Reaction identifiers, e.g. ``"(n, fission)"``, ``"(n, gamma)"``, + needed for the reaction rate tally. + """ + self._materials = materials + + # Convert reactions to MT values (needed when collapsing) + mt_values = {v: k for k, v in REACTION_NAME.items()} + mt_values['fission'] = 18 + self._mts = [mt_values[x] for x in scores] + + # Create flux tally with material and energy filters + self._flux_tally = Tally() + self._flux_tally.writable = False + self._flux_tally.filters = [ + MaterialFilter(materials), + EnergyFilter(self._energies) + ] + self._flux_tally.scores = ['flux'] + + def get_material_rates(self, mat_index, nuc_index, react_index): + """Return an array of reaction rates for a material + + Parameters + ---------- + mat_index : int + Index for material + nuc_index : iterable of int + Index for each nuclide in :attr:`nuclides` in the + desired reaction rate matrix + react_index : iterable of int + Index for each reaction scored in the tally + + Returns + ------- + rates : numpy.ndarray + Array with shape ``(n_nuclides, n_rxns)`` with the reaction rates in + this material + + """ + self._results_cache.fill(0.0) + + # Get flux for specified material + shape = (len(self._materials), len(self._energies) - 1) + mean_value = self._flux_tally.mean.reshape(shape) + flux = mean_value[mat_index] + + mat = self._materials[mat_index] + for name, i_nuc in zip(self.nuclides, nuc_index): + # Determine density of nuclide + # TODO: O(N) search for nuclide is not ideal + for nucname, density in zip(mat.nuclides, mat.densities): + if nucname == name: + break + else: + raise ValueError("Couldn't find {} in material {}".format(name, mat_id)) + + for mt, i_react in zip(self._mts, react_index): + # Use flux to collapse reaction rate (per N) + nuc = openmc.lib.nuclides[name] + rate_per_nuc = nuc.collapse_rate(mt, self._energies, flux) + + # Multiply by density to get absolute reaction rate + self._results_cache[i_nuc, i_react] = rate_per_nuc * density + + return self._results_cache + # ------------------------------------------ # Helpers for obtaining normalization factor # ------------------------------------------ diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index c1869a7302..f1a2912e9a 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -27,7 +27,7 @@ from .results_list import ResultsList from .helpers import ( DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper, FissionYieldCutoffHelper, AveragedFissionYieldHelper, EnergyScoreHelper, - SourceRateHelper) + SourceRateHelper, FluxCollapseHelper) __all__ = ["Operator", "OperatorResult"] @@ -113,6 +113,18 @@ class Operator(TransportOperator): ``fission_yield_mode``. Will be passed directly on to the helper. Passing a value of None will use the defaults for the associated helper. + reaction_rate_mode : {"direct", "flux"} + Indicate how one-group reaction rates should be calculated. The "direct" + method tallies transmutation reaction rates directly. The "flux" method + tallies a multigroup flux spectrum and then collapses one-group reaction + rates after a transport solve. + + .. versionadded:: 0.12.1 + reaction_rate_energies : iterable of float + Energy group boundaries that are to be used for calculating a multigroup + flux spectrum when the "flux" based ``reaction_rate_mode`` is being used. + + .. versionadded:: 0.12.1 reduce_chain : bool, optional If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the depletion chain up to ``reduce_chain_level``. Default is False. @@ -170,6 +182,7 @@ class Operator(TransportOperator): diff_burnable_mats=False, normalization_mode="fission-q", fission_q=None, dilute_initial=1.0e3, fission_yield_mode="constant", fission_yield_opts=None, + reaction_rate_mode="direct", reaction_rate_energies=None, reduce_chain=False, reduce_chain_level=None): check_value('fission yield mode', fission_yield_mode, self._fission_helpers.keys()) @@ -239,8 +252,24 @@ class Operator(TransportOperator): self.local_mats, self._burnable_nucs, self.chain.reactions) # Get classes to assist working with tallies - self._rate_helper = DirectReactionRateHelper( - self.reaction_rates.n_nuc, self.reaction_rates.n_react) + if reaction_rate_mode == "direct": + self._rate_helper = DirectReactionRateHelper( + self.reaction_rates.n_nuc, self.reaction_rates.n_react) + elif reaction_rate_mode == "flux": + # Ensure energy group boundaries were specified + if reaction_rate_energies is None: + raise ValueError( + "Energy group boundaries must be specified in the " + "reaction_rate_energies argument when reaction_rate_mode is" + "set to 'flux'.") + + self._rate_helper = FluxCollapseHelper( + self.reaction_rates.n_nuc, + self.reaction_rates.n_react, + reaction_rate_energies + ) + else: + raise ValueError("Invalid reaction rate mode.") if normalization_mode == "fission-q": self._normalization_helper = ChainFissionHelper()