Return energy in J/s/source neutron for EnergyHelper

Addressing comments in review for #1278

- Documentation cleanup
- Better naming convention regarding ReactionRateHelper
  results cache
- The power in Operator tally unpacking and normalizing is
  no longer converted to eV/s, since the EnergyHelper.energy
  property is now returned in J/s/source neutron
This commit is contained in:
Andrew Johnson 2019-07-23 08:29:37 -05:00
parent 808a41b4e3
commit 9f41dce1d2
No known key found for this signature in database
GPG key ID: 253418E91B7F6FEB
4 changed files with 55 additions and 41 deletions

View file

@ -75,8 +75,8 @@ data, such as number densities and reaction rates for each material.
:template: myclass.rst
AtomNumber
ChainFissHelper
DirectRxnRateHelper
ChainFissionHelper
DirectReactionRateHelper
OperatorResult
ReactionRates
Results
@ -92,7 +92,7 @@ The following classes are abstract classes that can be used to extend the
:template: myclass.rst
ReactionRateHelper
FissionEnergyHelper
EnergyHelper
TransportOperator
Each of the integrator functions also relies on a number of "helper" functions

View file

@ -11,9 +11,9 @@ from abc import ABC, abstractmethod
from xml.etree import ElementTree as ET
from warnings import warn
from numpy import zeros, nonzero
from numpy import nonzero, empty
from openmc.data import DataLibrary
from openmc.data import DataLibrary, JOULE_PER_EV
from openmc.checkvalue import check_type
from .chain import Chain
@ -188,7 +188,7 @@ class ReactionRateHelper(ABC):
@abstractmethod
def generate_tallies(self, materials, scores):
"""Use the capi to build tallies needed for reaction rates"""
"""Use the C API to build tallies needed for reaction rates"""
@property
def nuclides(self):
@ -201,13 +201,15 @@ class ReactionRateHelper(ABC):
self._nuclides = nuclides
self._rate_tally.nuclides = nuclides
def _reset_results_cache(self, nnucs, nreact):
"""Cache for results for a given material"""
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 = zeros((nnucs, nreact))
else:
self._results_cache.fill(0.0)
self._results_cache = empty((nnucs, nreact))
return self._results_cache
@abstractmethod
@ -216,6 +218,8 @@ class ReactionRateHelper(ABC):
Parameters
----------
mat_id : int
Unique ID for the requested material
nuc_index : list of str
Ordering of desired nuclides
react_index : list of str
@ -230,17 +234,13 @@ class ReactionRateHelper(ABC):
Parameters
----------
energy : float
Energy produced in this region [W]
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`
2D array ``[n_nuclides, n_rxns]`` of reaction rates normalized by
the number of nuclides
results : numpy.ndarray
Array of reactions rates of shape ``(n_nuclides, n_rxns)``
normalized by the number of nuclides
"""
mask = nonzero(number)
@ -266,7 +266,7 @@ class EnergyHelper(ABC):
All nuclides with desired reaction rates. Ordered to be
consistent with :class:`openmc.deplete.Operator`
energy : float
Total energy [eV/s] produced in a transport simulation.
Total energy [J/s/source neutron] produced in a transport simulation.
Updated in the material iteration with :meth:`update`.
"""
@ -276,7 +276,7 @@ class EnergyHelper(ABC):
@property
def energy(self):
return self._energy
return self._energy * JOULE_PER_EV
def reset(self):
"""Reset energy produced prior to unpacking tallies"""

View file

@ -14,7 +14,14 @@ from .abc import ReactionRateHelper, EnergyHelper
class DirectReactionRateHelper(ReactionRateHelper):
"""Class that generates tallies for one-group rates"""
"""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
@ -41,7 +48,7 @@ class DirectReactionRateHelper(ReactionRateHelper):
Parameters
----------
mat_id : int
Unique id for the requested material
Unique ID for the requested material
nuc_index : iterable of int
Index for each nuclide in :attr:`nuclides` in the
desired reaction rate matrix
@ -50,11 +57,12 @@ class DirectReactionRateHelper(ReactionRateHelper):
Returns
-------
rates : :class:`numpy.ndarray`
2D matrix ``(len(nuc_index), len(react_index))`` with the
rates : numpy.ndarray
Array with shape ``(n_nuclides, n_rxns)`` with the
reaction rates in this material
"""
results = self._reset_results_cache(len(nuc_index), len(react_index))
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)):
@ -63,13 +71,23 @@ class DirectReactionRateHelper(ReactionRateHelper):
return results
# ------------------------------------
# Helpers for obtaining fission energy
# ------------------------------------
# ----------------------------
# Helpers for obtaining energy
# ----------------------------
class ChainFissionHelper(EnergyHelper):
"""Fission Q-values are pulled from chain"""
"""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__()
@ -78,8 +96,8 @@ class ChainFissionHelper(EnergyHelper):
def prepare(self, chain_nucs, rate_index, _materials):
"""Populate the fission Q value vector from a chain.
Paramters
---------
Parameters
----------
chain_nucs : iterable of :class:`openmc.deplete.Nuclide`
Nuclides used in this depletion chain. Do not need
to be ordered

View file

@ -19,7 +19,6 @@ 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
@ -138,7 +137,11 @@ class Operator(TransportOperator):
openmc.reset_auto_ids()
self.burnable_mats, volume, nuclides = self._get_burnable_mats()
self.local_mats = _distribute(self.burnable_mats)
self._mat_index_map = {}
# 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()
@ -158,7 +161,6 @@ class Operator(TransportOperator):
self._rate_helper = DirectReactionRateHelper()
self._energy_helper = ChainFissionHelper()
def __call__(self, vec, power, print_out=True):
"""Runs a simulation.
@ -387,10 +389,6 @@ class Operator(TransportOperator):
self._energy_helper.prepare(
self.chain.nuclides, self.reaction_rates.index_nuc, materials)
# Generate map from local materials => material index
self._mat_index_map = {
lm: self.burnable_mats.index(lm) for lm in self.local_mats}
# Return number density vector
return list(self.number.get_mat_slice(np.s_[:]))
@ -564,11 +562,9 @@ class Operator(TransportOperator):
rates[i] = self._rate_helper.divide_by_adens(number)
# Reduce energy produced from all processes
# J / s / source neutron
energy = comm.allreduce(self._energy_helper.energy)
# Determine power in eV/s
power /= JOULE_PER_EV
# Scale reaction rates to obtain units of reactions/sec
rates *= power / energy