mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 13:45:36 -04:00
Merge pull request #1278 from drewejohnson/feat-dep-tally-op
Tally and Q value helper for depletion Operators
This commit is contained in:
Paul Romano
GitHub
commit
b4b063482e
4 changed files with 355 additions and 66 deletions
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
139
openmc/deplete/helpers.py
Normal file
139
openmc/deplete/helpers.py
Normal file
|
|
@ -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)
|
||||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue