diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index ffbffbc734..3976fe7c91 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -24,6 +24,7 @@ from . import comm from .abc import TransportOperator, OperatorResult from .atom_number import AtomNumber from .reaction_rates import ReactionRates +from .tally_helpers import ChainFissTallyHelper def _distribute(items): @@ -152,6 +153,10 @@ class Operator(TransportOperator): self.reaction_rates = ReactionRates( self.local_mats, self._burnable_nucs, self.chain.reactions) + # Get class to assist working with tallies + self._tally_helper = ChainFissTallyHelper(len(self.local_mats)) + + def __call__(self, vec, power, print_out=True): """Runs a simulation. @@ -180,7 +185,7 @@ class Operator(TransportOperator): # Update material compositions and tally nuclides self._update_materials() - self._tally.nuclides = self._get_tally_nuclides() + self._tally_helper.reaction_tally.nuclides = self._get_tally_nuclides() # Run OpenMC openmc.capi.reset() @@ -373,7 +378,9 @@ 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._tally_helper.generate_tallies(materials, self.chain.reactions) # Return number density vector return list(self.number.get_mat_slice(np.s_[:])) @@ -482,27 +489,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 @@ -529,7 +515,7 @@ class Operator(TransportOperator): # Extract tally bins materials = self.burnable_mats - nuclides = self._tally.nuclides + nuclides = self._tally_helper.reaction_tally.nuclides # Form fast map nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides] @@ -544,19 +530,14 @@ class Operator(TransportOperator): # 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) 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 + self._tally_helper.set_fission_q(self.chain.nuclides, rates.index_nuc) + + tally_results = self._tally_helper.reaction_tally.results # Extract results for i, mat in enumerate(self.local_mats): @@ -564,7 +545,7 @@ class Operator(TransportOperator): slab = materials.index(mat) # Get material results hyperslab - results = self._tally.results[slab, :, 1] + results = tally_results[slab, :, 1] # Zero out reaction rates and nuclide numbers rates_expanded[:] = 0.0 @@ -579,7 +560,7 @@ class Operator(TransportOperator): j += 1 # Accumulate energy from fission - energy += np.dot(rates_expanded[:, fission_ind], fission_Q) + energy += self._tally_helper.get_fiss_energy(rates_expanded[:, fission_ind], i) # Divide by total number and store for i_nuc_results in nuc_ind: diff --git a/openmc/deplete/tally_helpers.py b/openmc/deplete/tally_helpers.py new file mode 100644 index 0000000000..47c78ba8e0 --- /dev/null +++ b/openmc/deplete/tally_helpers.py @@ -0,0 +1,50 @@ +""" +Class for normalizing fission energy deposition +""" + +from numpy import dot, zeros + +from openmc.capi import Tally, MaterialFilter + + +class ChainFissTallyHelper(object): + """Fission Q-values are pulled from chain""" + + def __init__(self, n_materials): + self._fiss_q = None + self.n_materials = n_materials + self._rx_tally = None + + def set_fission_q(self, chain_nucs, rate_index): + if (self._fiss_q is not None + and self._fiss_q.shape == (len(rate_index), )): + return + + fq = zeros(len(rate_index)) + + for nuclide in chain_nucs: + if nuclide.name in rate_index: + for rx in nuclide.reactions: + if rx.type == "fission": + fq[rate_index[nuclide.name]] = rx.Q + break + + self._fiss_q = fq + + @property + def reaction_tally(self): + if self._rx_tally is None: + raise AttributeError( + "Reaction tally for {} not set.".format( + self.__class__.__name__ + ) + ) + return self._rx_tally + + def generate_tallies(self, materials, scores): + self._rx_tally = Tally() + self._rx_tally.scores = scores + self._rx_tally.filters = [MaterialFilter(materials)] + + def get_fiss_energy(self, fiss_rates, _mat_index): + return dot(fiss_rates, self._fiss_q)