diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 532e06655..70d48f931 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -75,10 +75,24 @@ data, such as number densities and reaction rates for each material. :template: myclass.rst AtomNumber + ChainFissionHelper + DirectReactionRateHelper OperatorResult ReactionRates Results ResultsList + + +The following classes are abstract classes that can be used to extend the +:mod:`openmc.deplete` capabilities: + +.. autosummary:: + :toctree:generated + :nosignatures: + :template: myclass.rst + + ReactionRateHelper + EnergyHelper TransportOperator Each of the integrator functions also relies on a number of "helper" functions diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 1be3e27ff..c5c086854 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -7,11 +7,14 @@ to run a full depletion simulation. from collections import namedtuple import os from pathlib import Path -from abc import ABCMeta, abstractmethod +from abc import ABC, abstractmethod from xml.etree import ElementTree as ET from warnings import warn -from openmc.data import DataLibrary +from numpy import nonzero, empty + +from openmc.data import DataLibrary, JOULE_PER_EV +from openmc.checkvalue import check_type from .chain import Chain OperatorResult = namedtuple('OperatorResult', ['k', 'rates']) @@ -34,7 +37,7 @@ except AttributeError: pass -class TransportOperator(metaclass=ABCMeta): +class TransportOperator(ABC): """Abstract class defining a transport operator Each depletion integrator is written to work with a generic transport @@ -160,3 +163,167 @@ class TransportOperator(metaclass=ABCMeta): def finalize(self): pass + + +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. + + Reaction rates are passed back to the operator for be used in + an :class:`openmc.deplete.OperatorResult` instance + + Attributes + ---------- + nuclides : list of str + All nuclides with desired reaction rates. Ordered to be + consistent with :class:`openmc.deplete.Operator` + """ + + def __init__(self): + self._nuclides = None + self._rate_tally = None + self._results_cache = None + + @abstractmethod + def generate_tallies(self, materials, scores): + """Use the C API to build tallies needed for reaction rates""" + + @property + def nuclides(self): + """List of nuclides with requested reaction rates""" + return self._nuclides + + @nuclides.setter + def nuclides(self, nuclides): + check_type("nuclides", nuclides, list, str) + self._nuclides = nuclides + self._rate_tally.nuclides = nuclides + + def _get_results_cache(self, nnucs, nreact): + """Cache for results for a given material + + Creates an empty array of shape ``(nnucs, nreact)`` + if the shape does not match the current cache. + """ + if (self._results_cache is None + or self._results_cache.shape != (nnucs, nreact)): + self._results_cache = empty((nnucs, nreact)) + return self._results_cache + + @abstractmethod + def get_material_rates(self, mat_id, nuc_index, react_index): + """Return 2D array of [nuclide, reaction] reaction rates + + Parameters + ---------- + mat_id : int + Unique ID for the requested material + nuc_index : list of str + Ordering of desired nuclides + react_index : list of str + Ordering of reactions + """ + + def divide_by_adens(self, number): + """Normalize reaction rates by number of nuclides + + Acts on the current material examined by + :meth:`get_material_rates` + + Parameters + ---------- + number : iterable of float + Number density [atoms/b-cm] of each nuclide tracked in the calculation. + Ordered identically to :attr:`nuclides` + + Returns + ------- + results : numpy.ndarray + Array of reactions rates of shape ``(n_nuclides, n_rxns)`` + normalized by the number of nuclides + """ + + mask = nonzero(number) + results = self._results_cache + for col in range(results.shape[1]): + results[mask, col] /= number[mask] + return results + + +class EnergyHelper(ABC): + """Abstract class for obtaining energy produced + + The ultimate goal of this helper is to provide instances of + :class:`openmc.deplete.Operator` with the total energy produced + in a transport simulation. This information, provided with the + power requested by the user and reaction rates from a + :class:`ReactionRateHelper` will scale reaction rates to the + correct values. + + Attributes + ---------- + nuclides : list of str + All nuclides with desired reaction rates. Ordered to be + consistent with :class:`openmc.deplete.Operator` + energy : float + Total energy [J/s/source neutron] produced in a transport simulation. + Updated in the material iteration with :meth:`update`. + """ + + def __init__(self): + self._nuclides = None + self._energy = 0.0 + + @property + def energy(self): + return self._energy * JOULE_PER_EV + + def reset(self): + """Reset energy produced prior to unpacking tallies""" + self._energy = 0.0 + + @abstractmethod + def prepare(self, chain_nucs, rate_index, materials): + """Perform work needed to obtain energy produced + + This method is called prior to the transport simulations + in :meth:`openmc.deplete.Operator.initial_condition`. + + Parameters + ---------- + chain_nucs : list of str + All nuclides to be tracked in this problem + rate_index : dict of str to int + Mapping from nuclide name to index in the + `fission_rates` for :meth:`update`. + materials : list of str + All materials tracked on the operator helped by this + object. Should correspond to + :attr:`openmc.deplete.Operator.burnable_materials` + """ + + def update(self, fission_rates, mat_index): + """Update the energy produced + + 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 + Index for the specific material in the list of all burnable + materials. + """ + + @property + def nuclides(self): + """List of nuclides with requested reaction rates""" + return self._nuclides + + @nuclides.setter + def nuclides(self, nuclides): + check_type("nuclides", nuclides, list, str) + self._nuclides = nuclides diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py new file mode 100644 index 000000000..edf0feb06 --- /dev/null +++ b/openmc/deplete/helpers.py @@ -0,0 +1,139 @@ +""" +Class for normalizing fission energy deposition +""" +from itertools import product + +from numpy import dot, zeros + +from openmc.capi import Tally, MaterialFilter +from .abc import ReactionRateHelper, EnergyHelper + +# ------------------------------------- +# Helpers for generating reaction rates +# ------------------------------------- + + +class DirectReactionRateHelper(ReactionRateHelper): + """Class that generates tallies for one-group rates + + Attributes + ---------- + nuclides : list of str + All nuclides with desired reaction rates. Ordered to be + consistent with :class:`openmc.deplete.Operator` + """ + + def generate_tallies(self, materials, scores): + """Produce one-group reaction rate tally + + Uses the :mod:`openmc.capi` 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` + scores : iterable of str + Reaction identifiers, e.g. ``"(n, fission)"``, + ``"(n, gamma)"``, needed for the reaction rate tally. + """ + self._rate_tally = Tally() + self._rate_tally.scores = scores + self._rate_tally.filters = [MaterialFilter(materials)] + + def get_material_rates(self, mat_id, nuc_index, react_index): + """Return an array of reaction rates for a material + + Parameters + ---------- + mat_id : int + Unique ID for the requested 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 + """ + results = self._get_results_cache(len(nuc_index), len(react_index)) + results.fill(0.0) + full_tally_res = self._rate_tally.results[mat_id, :, 1] + for i_tally, (i_nuc, i_react) in enumerate( + product(nuc_index, react_index)): + results[i_nuc, i_react] = full_tally_res[i_tally] + + return results + + +# ---------------------------- +# Helpers for obtaining energy +# ---------------------------- + + +class ChainFissionHelper(EnergyHelper): + """Computes energy using fission Q values from depletion chain + + Attributes + ---------- + nuclides : list of str + All nuclides with desired reaction rates. Ordered to be + consistent with :class:`openmc.deplete.Operator` + energy : float + Total energy [J/s/source neutron] produced in a transport simulation. + Updated in the material iteration with :meth:`update`. + """ + + def __init__(self): + super().__init__() + self._fission_q_vector = None + + def prepare(self, chain_nucs, rate_index, _materials): + """Populate the fission Q value vector from a chain. + + Parameters + ---------- + chain_nucs : iterable of :class:`openmc.deplete.Nuclide` + Nuclides used in this depletion chain. Do not need + to be ordered + rate_index : dict of str to int + Dictionary mapping names of nuclides, e.g. ``"U235"``, + to a corresponding index in the desired fission Q + vector. + _materials : list of str + Unused. Materials to be tracked for this helper. + """ + if (self._fission_q_vector is not None + and self._fission_q_vector.shape == (len(rate_index),)): + return + + fission_qs = zeros(len(rate_index)) + + for nuclide in chain_nucs: + if nuclide.name in rate_index: + for rx in nuclide.reactions: + if rx.type == "fission": + fission_qs[rate_index[nuclide.name]] = rx.Q + break + + self._fission_q_vector = fission_qs + + def update(self, fission_rates, _mat_index): + """Update energy produced with fission rates in a material + + 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 + index for the material requested. Unused, as identical + isotopes in all materials have the same Q value. + """ + self._energy += dot(fission_rates, self._fission_q_vector) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 21c96b047..bb3c96e4f 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -19,11 +19,11 @@ import numpy as np import openmc import openmc.capi -from openmc.data import JOULE_PER_EV from . import comm from .abc import TransportOperator, OperatorResult from .atom_number import AtomNumber from .reaction_rates import ReactionRates +from .helpers import DirectReactionRateHelper, ChainFissionHelper def _distribute(items): @@ -138,6 +138,11 @@ class Operator(TransportOperator): self.burnable_mats, volume, nuclides = self._get_burnable_mats() self.local_mats = _distribute(self.burnable_mats) + # Generate map from local materials => material index + self._mat_index_map = { + lm: self.burnable_mats.index(lm) for lm in self.local_mats} + + # Determine which nuclides have incident neutron data self.nuclides_with_data = self._get_nuclides_with_data() @@ -152,6 +157,10 @@ class Operator(TransportOperator): self.reaction_rates = ReactionRates( self.local_mats, self._burnable_nucs, self.chain.reactions) + # Get classes to assist working with tallies + self._rate_helper = DirectReactionRateHelper() + self._energy_helper = ChainFissionHelper() + def __call__(self, vec, power, print_out=True): """Runs a simulation. @@ -180,7 +189,8 @@ class Operator(TransportOperator): # Update material compositions and tally nuclides self._update_materials() - self._tally.nuclides = self._get_tally_nuclides() + self._rate_helper.nuclides = self._get_tally_nuclides() + self._energy_helper.nuclides = self._rate_helper.nuclides # Run OpenMC openmc.capi.reset() @@ -373,7 +383,11 @@ class Operator(TransportOperator): openmc.capi.init(intracomm=comm) # Generate tallies in memory - self._generate_tallies() + materials = [openmc.capi.materials[int(i)] + for i in self.burnable_mats] + self._rate_helper.generate_tallies(materials, self.chain.reactions) + self._energy_helper.prepare( + self.chain.nuclides, self.reaction_rates.index_nuc, materials) # Return number density vector return list(self.number.get_mat_slice(np.s_[:])) @@ -482,27 +496,6 @@ class Operator(TransportOperator): nuc_list = comm.bcast(nuc_list) return [nuc for nuc in nuc_list if nuc in self.chain] - def _generate_tallies(self): - """Generates depletion tallies. - - Using information from the depletion chain as well as the nuclides - currently in the problem, this function automatically generates a - tally.xml for the simulation. - - """ - # Create tallies for depleting regions - materials = [openmc.capi.materials[int(i)] - for i in self.burnable_mats] - mat_filter = openmc.capi.MaterialFilter(materials) - - # Set up a tally that has a material filter covering each depletable - # material and scores corresponding to all reactions that cause - # transmutation. The nuclides for the tally are set later when eval() is - # called. - self._tally = openmc.capi.Tally() - self._tally.scores = self.chain.reactions - self._tally.filters = [mat_filter] - def _unpack_tallies_and_normalize(self, power): """Unpack tallies from OpenMC and return an operator result @@ -523,78 +516,54 @@ class Operator(TransportOperator): """ rates = self.reaction_rates - rates[:, :, :] = 0.0 + rates.fill(0.0) # Get k and uncertainty k_combined = openmc.capi.keff() # Extract tally bins - materials = self.burnable_mats - nuclides = self._tally.nuclides + nuclides = self._rate_helper.nuclides # Form fast map nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides] react_ind = [rates.index_rx[react] for react in self.chain.reactions] # Compute fission power - # TODO : improve this calculation # Keep track of energy produced from all reactions in eV per source # particle - energy = 0.0 + self._energy_helper.reset() # Create arrays to store fission Q values, reaction rates, and nuclide - # numbers - fission_Q = np.zeros(rates.n_nuc) - rates_expanded = np.zeros((rates.n_nuc, rates.n_react)) - number = np.zeros(rates.n_nuc) + # numbers, zeroed out in material iteration + number = np.empty(rates.n_nuc) fission_ind = rates.index_rx["fission"] - for nuclide in self.chain.nuclides: - if nuclide.name in rates.index_nuc: - for rx in nuclide.reactions: - if rx.type == 'fission': - ind = rates.index_nuc[nuclide.name] - fission_Q[ind] = rx.Q - break - # Extract results for i, mat in enumerate(self.local_mats): # Get tally index - slab = materials.index(mat) - - # Get material results hyperslab - results = self._tally.results[slab, :, 1] + mat_index = self._mat_index_map[mat] # Zero out reaction rates and nuclide numbers - rates_expanded[:] = 0.0 - number[:] = 0.0 + number.fill(0.0) - # Expand into our memory layout - j = 0 + # Get new number densities for nuc, i_nuc_results in zip(nuclides, nuc_ind): number[i_nuc_results] = self.number[mat, nuc] - for react in react_ind: - rates_expanded[i_nuc_results, react] = results[j] - j += 1 + + tally_rates = self._rate_helper.get_material_rates( + mat_index, nuc_ind, react_ind) # Accumulate energy from fission - energy += np.dot(rates_expanded[:, fission_ind], fission_Q) + self._energy_helper.update(tally_rates[:, fission_ind], mat_index) # Divide by total number and store - for i_nuc_results in nuc_ind: - if number[i_nuc_results] != 0.0: - for react in react_ind: - rates_expanded[i_nuc_results, react] /= number[i_nuc_results] - - rates[i, :, :] = rates_expanded + rates[i] = self._rate_helper.divide_by_adens(number) # Reduce energy produced from all processes - energy = comm.allreduce(energy) - - # Determine power in eV/s - power /= JOULE_PER_EV + # J / s / source neutron + energy = comm.allreduce(self._energy_helper.energy) # Scale reaction rates to obtain units of reactions/sec rates *= power / energy