Refactor MicroXS class and usage in IndependentOperator (#2595)

Co-authored-by: Jonathan Shimwell <drshimwell@gmail.com>
Co-authored-by: Olek <45364492+yardasol@users.noreply.github.com>
This commit is contained in:
Paul Romano 2023-08-31 09:55:21 -05:00 committed by GitHub
parent c93feed91f
commit d21b653d32
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
14 changed files with 447 additions and 326 deletions

View file

@ -57,6 +57,16 @@ The :class:`CoupledOperator` and :class:`IndependentOperator` classes must also
have some knowledge of how nuclides transmute and decay. This is handled by the
:class:`Chain` class.
The :class:`IndependentOperator` class requires a set of fluxes and microscopic
cross sections. The following function can be used to generate this information:
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
get_microxs_and_flux
Minimal Example
---------------

View file

@ -197,43 +197,54 @@ across all material instances.
Transport-independent depletion
===============================
.. warning::
This category of operator uses multigroup microscopic cross sections along with
multigroup flux spectra to obtain transmutation reaction rates. The cross
sections are pre-calculated, so there is no need for direct coupling between a
transport-independent operator and a transport solver. The :mod:`openmc.deplete`
module offers a single transport-independent operator,
:class:`~openmc.deplete.IndependentOperator`, and only one operator is needed
since, in theory, any transport code could calculate the multigroup microscopic
cross sections. The :class:`~openmc.deplete.IndependentOperator` class has two
constructors. The default constructor requires a :class:`openmc.Materials`
instance, a list of multigroup flux arrays, and a list of
:class:`~openmc.deplete.MicroXS` instances containing multigroup microscopic
cross sections in units of barns. This might look like the following::
This feature is still under heavy development and has yet to be rigorously
verified. API changes and feature additions are possible and likely in
the near future.
This category of operator uses one-group microscopic cross sections to obtain
transmutation reaction rates. The cross sections are pre-calculated, so there is
no need for direct coupling between a transport-independent operator and a
transport solver. The :mod:`openmc.deplete` module offers a single
transport-independent operator, :class:`~openmc.deplete.IndependentOperator`,
and only one operator is needed since, in theory, any transport code could
calcuate the one-group microscopic cross sections.
The :class:`~openmc.deplete.IndependentOperator` class has two constructors.
The default constructor requires a :class:`openmc.Materials` instance, a
:class:`~openmc.deplete.MicroXS` instance containing one-group microscoic cross
sections in units of barns, and a path to a depletion chain file::
materials = openmc.Materials()
materials = openmc.Materials([m1, m2, m3])
...
# load in the microscopic cross sections
micro_xs = openmc.deplete.MicroXS.from_csv(micro_xs_path)
# Assign fluxes (generated from any code)
flux_m1 = numpy.array([...])
flux_m2 = numpy.array([...])
flux_m3 = numpy.array([...])
fluxes = [flux_m1, flux_m2, flux_m3]
op = openmc.deplete.IndependentOperator(materials, micro_xs, chain_file)
# Assign microscopic cross sections
micro_m1 = openmc.deplete.MicroXS.from_csv('xs_m1.csv')
micro_m2 = openmc.deplete.MicroXS.from_csv('xs_m2.csv')
micro_m3 = openmc.deplete.MicroXS.from_csv('xs_m3.csv')
micros = [micro_m1, micro_m2, micro_m3]
# Create operator
op = openmc.deplete.IndependentOperator(materials, fluxes, micros)
For more details on the :class:`~openmc.deplete.MicroXS` class, including how to
use OpenMC's transport solver to generate microscopic cross sections and fluxes
for use with :class:`~openmc.deplete.IndependentOperator`, see :ref:`micros`.
.. note::
The same statements from :ref:`coupled-depletion` about which
materials are depleted and the requirement for depletable materials to have
a specified volume also apply here.
The same statements from :ref:`coupled-depletion` about which materials are
depleted and the requirement for depletable materials to have a specified
volume also apply here.
An alternate constructor,
:meth:`~openmc.deplete.IndependentOperator.from_nuclides`, accepts a volume and
dictionary of nuclide concentrations in place of the :class:`openmc.Materials`
instance::
instance. Note that while the normal constructor allows multiple materials to be
depleted with a single operator, the
:meth:`~openmc.deplete.IndependentOperator.from_nuclides` classmethod only works
for a single material::
nuclides = {'U234': 8.92e18,
'U235': 9.98e20,
@ -244,6 +255,7 @@ instance::
volume = 0.5
op = openmc.deplete.IndependentOperator.from_nuclides(volume,
nuclides,
flux,
micro_xs,
chain_file,
nuc_units='atom/cm3')
@ -253,18 +265,21 @@ transport-depletion calculation and follow the same steps from there.
.. note::
Ideally, one-group cross section data should be available for every
reaction in the depletion chain. If cross section data is not present for
a nuclide in the depletion chain with at least one reaction, that reaction
will not be simulated.
Ideally, multigroup cross section data should be available for every reaction
in the depletion chain. If cross section data is not present for a nuclide in
the depletion chain with at least one reaction, that reaction will not be
simulated.
.. _micros:
Loading and Generating Microscopic Cross Sections
-------------------------------------------------
As mentioned earlier, any transport code could be used to calculate one-group
microscopic cross sections. The :mod:`openmc.deplete` module provides the
:class:`~openmc.deplete.MicroXS` class, which contains methods to read in
pre-calculated cross sections from a ``.csv`` file or from data arrays::
As mentioned above, any transport code could be used to calculate multigroup
microscopic cross sections and fluxes. The :mod:`openmc.deplete` module provides
the :class:`~openmc.deplete.MicroXS` class, which can either be instantiated
from pre-calculated cross sections in a ``.csv`` file or from data arrays
directly::
micro_xs = MicroXS.from_csv(micro_xs_path)
@ -273,37 +288,31 @@ pre-calculated cross sections from a ``.csv`` file or from data arrays::
data = np.array([[0.1, 0.2],
[0.3, 0.4],
[0.01, 0.5]])
micro_xs = MicroXS.from_array(nuclides, reactions, data)
micro_xs = MicroXS(data, nuclides, reactions)
.. important::
Both :meth:`~openmc.deplete.MicroXS.from_csv()` and
:meth:`~openmc.deplete.MicroXS.from_array()` assume the cross section values
provided are in barns by defualt, but have no way of verifying this. Make
sure your cross sections are in the correct units before passing to a
The cross section values are assumed to be in units of barns. Make sure your
cross sections are in the correct units before passing to a
:class:`~openmc.deplete.IndependentOperator` object.
The :class:`~openmc.deplete.MicroXS` class also contains a method to generate one-group microscopic cross sections using OpenMC's transport solver. The
:meth:`~openmc.deplete.MicroXS.from_model()` method will produce a
:class:`~openmc.deplete.MicroXS` instance with microscopic cross section data in
units of barns::
Additionally, a convenience function,
:func:`~openmc.deplete.get_microxs_and_flux`, can provide the needed fluxes and
cross sections using OpenMC's transport solver::
import openmc
model = openmc.Model()
...
model = openmc.Model.from_xml()
fluxes, micros = openmc.deplete.get_microxs_and_flux(model, materials)
micro_xs = openmc.deplete.MicroXS.from_model(model,
model.materials[0],
chain_file)
If you are running :meth:`~openmc.deplete.MicroXS.from_model()` on a cluster
If you are running :func:`~openmc.deplete.get_microxs_and_flux` on a cluster
where temporary files are created on a local filesystem that is not shared
across nodes, you'll need to set an environment variable pointing to a local
directoy so that each MPI process knows where to store output files used to
calculate the microscopic cross sections. In order of priority, they are
:envvar:`TMPDIR`. :envvar:`TEMP`, and :envvar:`TMP`. Users interested in
further details can read the documentation for the `tempfile <https://docs.python.org/3/library/tempfile.html#tempfile.gettempdir>`_ module.
:envvar:`TMPDIR`. :envvar:`TEMP`, and :envvar:`TMP`. Users interested in further
details can read the documentation for the `tempfile
<https://docs.python.org/3/library/tempfile.html#tempfile.gettempdir>`_ module.
Caveats
-------
@ -325,24 +334,30 @@ normalizing reaction rates:
the time integrator is a flux, and obtains the reaction rates by multiplying
the cross sections by the ``source-rate``.
2. ``fission-q`` normalization, which uses the ``power`` or ``power_density``
provided by the time integrator to obtain reaction rates by computing a value
for the flux based on this power. The equation we use for this calculation is
provided by the time integrator to obtain normalized reaction rates by
computing a normalization factor as the ratio of the user-specified power to
the "observed" power based on fission reaction rates. The equation for the
normalization factor is
.. math::
:label: fission-q
\phi = \frac{P}{\sum\limits_i (Q_i \sigma^f_i N_i)}
f = \frac{P}{\sum\limits_m \sum\limits_i \left(Q_i N_{i,m} \sum\limits_g
\sigma^f_{i,g,m} \phi_{g,m} \right)}
where :math:`P` is the power, :math:`Q_i` is the fission Q value for nuclide
:math:`i`, :math:`\sigma_i^f` is the microscopic fission cross section for
nuclide :math:`i`, and :math:`N_i` is the number of atoms of nuclide
:math:`i`. This equation makes the same assumptions and issues as discussed
in :ref:`energy-deposition`. Unfortunately, the proposed solution in that
section does not apply here since we are decoupled from transport code.
However, there is a method to converge to a more accurate value for flux by
using substeps during time integration. `This paper
<https://doi.org/10.1016/j.anucene.2016.05.031>`_ provides a good discussion
of this method.
:math:`i`, :math:`\sigma_{i,g,m}^f` is the microscopic fission cross section
for nuclide :math:`i` in energy group :math:`g` for material :math:`m`,
:math:`\phi_{g,m}` is the neutron flux in group :math:`g` for material
:math:`m`, and :math:`N_{i,m}` is the number of atoms of nuclide :math:`i`
for material :math:`m`. Reaction rates are then multiplied by :math:`f` so
that the total fission power matches :math:`P`. This equation makes the same
assumptions and issues as discussed in :ref:`energy-deposition`.
Unfortunately, the proposed solution in that section does not apply here
since we are decoupled from transport code. However, there is a method to
converge to a more accurate value for flux by using substeps during time
integration. `This paper <https://doi.org/10.1016/j.anucene.2016.05.031>`_
provides a good discussion of this method.
.. warning::
@ -359,14 +374,14 @@ useful for running many different cases of a particular scenario. A
transport-independent depletion simulation using ``fission-q`` normalization
will sum the fission energy values across all materials into :math:`Q_i` in
Equation :math:numref:`fission-q`, and Equation :math:numref:`fission-q`
provides the flux we use to calculate the reaction rates in each material.
provides the normalization factor applied to reaction rates in each material.
This can be useful for running a scenario with multiple depletable materials
that are part of the same reactor. This behavior may change in the future.
Time integration
~~~~~~~~~~~~~~~~
The values of the one-group microscopic cross sections passed to
The values of the microscopic cross sections passed to
:class:`openmc.deplete.IndependentOperator` are fixed for the entire depletion
simulation. This implicit assumption may produce inaccurate results for certain
scenarios.

View file

@ -196,13 +196,13 @@ class DirectReactionRateHelper(ReactionRateHelper):
"""
self._rate_tally_means_cache = None
def get_material_rates(self, mat_id, nuc_index, rx_index):
def get_material_rates(self, mat_index, nuc_index, rx_index):
"""Return an array of reaction rates for a material
Parameters
----------
mat_id : int
Unique ID for the requested material
mat_index : int
Index for the material
nuc_index : iterable of int
Index for each nuclide in :attr:`nuclides` in the
desired reaction rate matrix
@ -216,7 +216,7 @@ class DirectReactionRateHelper(ReactionRateHelper):
reaction rates in this material
"""
self._results_cache.fill(0.0)
full_tally_res = self.rate_tally_means[mat_id]
full_tally_res = self.rate_tally_means[mat_index]
for i_tally, (i_nuc, i_rx) in enumerate(product(nuc_index, rx_index)):
self._results_cache[i_nuc, i_rx] = full_tally_res[i_tally]

View file

@ -5,8 +5,10 @@ transport solver by using user-provided one-group cross sections.
"""
from __future__ import annotations
from collections.abc import Iterable
import copy
from itertools import product
from typing import List, Set
import numpy as np
from uncertainties import ufloat
@ -37,14 +39,20 @@ class IndependentOperator(OpenMCOperator):
.. versionadded:: 0.13.1
.. versionchanged:: 0.13.4
Arguments updated to include list of fluxes and microscopic cross
sections.
Parameters
----------
materials : openmc.Materials
Materials to deplete.
micro_xs : MicroXS
One-group microscopic cross sections in [b]. If the
:class:`~openmc.deplete.MicroXS` object is empty, a decay-only calculation will
be run.
fluxes : list of numpy.ndarray
Flux in each group in [n-cm/src] for each domain
micros : list of MicroXS
Cross sections in [b] for each domain. If the
:class:`~openmc.deplete.MicroXS` object is empty, a decay-only
calculation will be run.
chain_file : str
Path to the depletion chain XML file. Defaults to
``openmc.config['chain_file']``.
@ -109,7 +117,8 @@ class IndependentOperator(OpenMCOperator):
def __init__(self,
materials,
micro_xs,
fluxes,
micros,
chain_file=None,
keff=None,
normalization_mode='fission-q',
@ -120,7 +129,7 @@ class IndependentOperator(OpenMCOperator):
fission_yield_opts=None):
# Validate micro-xs parameters
check_type('materials', materials, openmc.Materials)
check_type('micro_xs', micro_xs, MicroXS)
check_type('micros', micros, Iterable, MicroXS)
if keff is not None:
check_type('keff', keff, tuple, float)
keff = ufloat(*keff)
@ -132,10 +141,15 @@ class IndependentOperator(OpenMCOperator):
helper_kwargs = {'normalization_mode': normalization_mode,
'fission_yield_opts': fission_yield_opts}
cross_sections = micro_xs * 1e-24
# Sort fluxes and micros in same order that materials get sorted
index_sort = np.argsort([mat.id for mat in materials])
fluxes = [fluxes[i] for i in index_sort]
micros = [micros[i] for i in index_sort]
self.fluxes = fluxes
super().__init__(
materials,
cross_sections,
micros,
chain_file,
prev_results,
fission_q=fission_q,
@ -145,6 +159,7 @@ class IndependentOperator(OpenMCOperator):
@classmethod
def from_nuclides(cls, volume, nuclides,
flux,
micro_xs,
chain_file=None,
nuc_units='atom/b-cm',
@ -163,10 +178,11 @@ class IndependentOperator(OpenMCOperator):
nuclides : dict of str to float
Dictionary with nuclide names as keys and nuclide concentrations as
values.
flux : numpy.ndarray
Flux in each group in [n-cm/src]
micro_xs : MicroXS
One-group microscopic cross sections in [b]. If the
:class:`~openmc.deplete.MicroXS` object is empty, a decay-only calculation
will be run.
Cross sections in [b]. If the :class:`~openmc.deplete.MicroXS`
object is empty, a decay-only calculation will be run.
chain_file : str, optional
Path to the depletion chain XML file. Defaults to
``openmc.config['chain_file']``.
@ -203,8 +219,11 @@ class IndependentOperator(OpenMCOperator):
"""
check_type('nuclides', nuclides, dict, str)
materials = cls._consolidate_nuclides_to_material(nuclides, nuc_units, volume)
fluxes = [flux]
micros = [micro_xs]
return cls(materials,
micro_xs,
fluxes,
micros,
chain_file,
keff=keff,
normalization_mode=normalization_mode,
@ -256,9 +275,9 @@ class IndependentOperator(OpenMCOperator):
self.prev_res.append(new_res)
def _get_nuclides_with_data(self, cross_sections):
def _get_nuclides_with_data(self, cross_sections: List[MicroXS]) -> Set[str]:
"""Finds nuclides with cross section data"""
return set(cross_sections.index)
return set(cross_sections[0].nuclides)
class _IndependentRateHelper(ReactionRateHelper):
"""Class for generating one-group reaction rates with flux and
@ -285,7 +304,7 @@ class IndependentOperator(OpenMCOperator):
"""
def __init__(self, op):
def __init__(self, op: IndependentOperator):
rates = op.reaction_rates
super().__init__(rates.n_nuc, rates.n_react)
@ -301,13 +320,13 @@ class IndependentOperator(OpenMCOperator):
"""Unused in this case"""
pass
def get_material_rates(self, mat_id, nuc_index, react_index):
def get_material_rates(self, mat_index, nuc_index, react_index):
"""Return 2D array of [nuclide, reaction] reaction rates
Parameters
----------
mat_id : int
Unique ID for the requested material
mat_index : int
Index for the material
nuc_index : list of str
Ordering of desired nuclides
react_index : list of str
@ -315,20 +334,18 @@ class IndependentOperator(OpenMCOperator):
"""
self._results_cache.fill(0.0)
# Get volume in units of [b-cm]
volume_b_cm = 1e24 * self._op.number.get_mat_volume(mat_id)
# Get flux and microscopic cross sections from operator
flux = self._op.fluxes[mat_index]
xs = self._op.cross_sections[mat_index]
for i_nuc, i_react in product(nuc_index, react_index):
for i_nuc in nuc_index:
nuc = self.nuc_ind_map[i_nuc]
rx = self.rx_ind_map[i_react]
for i_rx in react_index:
rx = self.rx_ind_map[i_rx]
# OK, this is kind of weird, but we multiply by volume in [b-cm]
# only because OpenMCOperator._calculate_reaction_rates has to
# divide it out later. It might make more sense to account for
# the source rate (flux) here rather than in the normalization
# helper.
self._results_cache[i_nuc, i_react] = \
self._op.cross_sections[rx][nuc] * volume_b_cm
# Determine reaction rate by multiplying xs in [b] by flux
# in [n-cm/src] to give [(reactions/src)*b-cm/atom]
self._results_cache[i_nuc, i_rx] = (xs[nuc, rx] * flux).sum()
return self._results_cache

View file

@ -1,15 +1,17 @@
"""MicroXS module
A pandas.DataFrame storing microscopic cross section data with
nuclide names as row indices and reaction names as column indices.
A class for storing microscopic cross section data that can be used with the
IndependentOperator class for depletion.
"""
from __future__ import annotations
import tempfile
from typing import List, Tuple, Iterable, Optional, Union
from pandas import DataFrame, read_csv, Series
import pandas as pd
import numpy as np
from openmc.checkvalue import check_type, check_value, check_iterable_type
from openmc.checkvalue import check_type, check_value, check_iterable_type, PathLike
from openmc.exceptions import DataError
from openmc import StatePoint
import openmc
@ -20,168 +22,178 @@ _valid_rxns = list(REACTIONS)
_valid_rxns.append('fission')
class MicroXS(DataFrame):
def get_microxs_and_flux(
model: openmc.Model,
domains,
nuclides: Optional[Iterable[str]] = None,
reactions: Optional[Iterable[str]] = None,
energies: Optional[Union[Iterable[float], str]] = None,
chain_file: Optional[PathLike] = None,
run_kwargs=None
) -> Tuple[List[np.ndarray], List[MicroXS]]:
"""Generate a microscopic cross sections and flux from a Model
.. versionadded:: 0.13.4
Parameters
----------
model : openmc.Model
OpenMC model object. Must contain geometry, materials, and settings.
domains : list of openmc.Material or openmc.Cell or openmc.Universe
Domains in which to tally reaction rates.
nuclides : list of str
Nuclides to get cross sections for. If not specified, all burnable
nuclides from the depletion chain file are used.
reactions : list of str
Reactions to get cross sections for. If not specified, all neutron
reactions listed in the depletion chain file are used.
energies : iterable of float or str
Energy group boundaries in [eV] or the name of the group structure
chain_file : str, optional
Path to the depletion chain XML file that will be used in depletion
simulation. Used to determine cross sections for materials not
present in the inital composition. Defaults to
``openmc.config['chain_file']``.
run_kwargs : dict, optional
Keyword arguments passed to :meth:`openmc.Model.run`
Returns
-------
list of numpy.ndarray
Flux in each group in [n-cm/src] for each domain
list of MicroXS
Cross section data in [b] for each domain
"""
# Save any original tallies on the model
original_tallies = model.tallies
# Determine what reactions and nuclides are available in chain
if chain_file is None:
chain_file = openmc.config.get('chain_file')
if chain_file is None:
raise DataError(
"No depletion chain specified and could not find depletion "
"chain in openmc.config['chain_file']"
)
chain = Chain.from_xml(chain_file)
if reactions is None:
reactions = chain.reactions
if not nuclides:
cross_sections = _find_cross_sections(model)
nuclides_with_data = _get_nuclides_with_data(cross_sections)
nuclides = [nuc.name for nuc in chain.nuclides
if nuc.name in nuclides_with_data]
# Set up the reaction rate and flux tallies
if energies is None:
energy_filter = openmc.EnergyFilter([0.0, 100.0e6])
elif isinstance(energies, str):
energy_filter = openmc.EnergyFilter.from_group_structure(energies)
else:
energy_filter = openmc.EnergyFilter(energies)
if isinstance(domains[0], openmc.Material):
domain_filter = openmc.MaterialFilter(domains)
elif isinstance(domains[0], openmc.Cell):
domain_filter = openmc.CellFilter(domains)
elif isinstance(domains[0], openmc.Universe):
domain_filter = openmc.UniverseFilter(domains)
else:
raise ValueError(f"Unsupported domain type: {type(domains[0])}")
rr_tally = openmc.Tally(name='MicroXS RR')
rr_tally.filters = [domain_filter, energy_filter]
rr_tally.nuclides = nuclides
rr_tally.multiply_density = False
rr_tally.scores = reactions
flux_tally = openmc.Tally(name='MicroXS flux')
flux_tally.filters = [domain_filter, energy_filter]
flux_tally.scores = ['flux']
model.tallies = openmc.Tallies([rr_tally, flux_tally])
# create temporary run
with tempfile.TemporaryDirectory() as temp_dir:
if run_kwargs is None:
run_kwargs = {}
else:
run_kwargs = dict(run_kwargs)
run_kwargs.setdefault('cwd', temp_dir)
statepoint_path = model.run(**run_kwargs)
with StatePoint(statepoint_path) as sp:
rr_tally = sp.tallies[rr_tally.id]
rr_tally._read_results()
flux_tally = sp.tallies[flux_tally.id]
flux_tally._read_results()
# Get reaction rates and flux values
reaction_rates = rr_tally.get_reshaped_data() # (domains, groups, nuclides, reactions)
flux = flux_tally.get_reshaped_data() # (domains, groups, 1, 1)
# Make energy groups last dimension
reaction_rates = np.moveaxis(reaction_rates, 1, -1) # (domains, nuclides, reactions, groups)
flux = np.moveaxis(flux, 1, -1) # (domains, 1, 1, groups)
# Divide RR by flux to get microscopic cross sections
xs = np.empty_like(reaction_rates) # (domains, nuclides, reactions, groups)
d, _, _, g = np.nonzero(flux)
xs[d, ..., g] = reaction_rates[d, ..., g] / flux[d, :, :, g]
# Reset tallies
model.tallies = original_tallies
# Create lists where each item corresponds to one domain
fluxes = list(flux.squeeze((1, 2)))
micros = [MicroXS(xs_i, nuclides, reactions) for xs_i in xs]
return fluxes, micros
class MicroXS:
"""Microscopic cross section data for use in transport-independent depletion.
.. versionadded:: 0.13.1
.. versionchanged:: 0.13.4
Class was heavily refactored and no longer subclasses pandas.DataFrame.
Parameters
----------
data : numpy.ndarray of floats
3D array containing microscopic cross section values for each
nuclide, reaction, and energy group. Cross section values are assumed to
be in [b], and indexed by [nuclide, reaction, energy group]
nuclides : list of str
List of nuclide symbols for that have data for at least one
reaction.
reactions : list of str
List of reactions. All reactions must match those in
:data:`openmc.deplete.chain.REACTIONS`
"""
@classmethod
def from_model(cls,
model,
domain,
nuclides=None,
reactions=None,
chain_file=None,
energy_bounds=(0, 20e6),
run_kwargs=None):
"""Generate a one-group cross-section dataframe using OpenMC.
Note that the ``openmc`` executable must be compiled.
Parameters
----------
model : openmc.Model
OpenMC model object. Must contain geometry, materials, and settings.
domain : openmc.Material or openmc.Cell or openmc.Universe
Domain in which to tally reaction rates.
nuclides : list of str
Nuclides to get cross sections for. If not specified, all burnable
nuclides from the depletion chain file are used.
reactions : list of str
Reactions to get cross sections for. If not specified, all neutron
reactions listed in the depletion chain file are used.
chain_file : str, optional
Path to the depletion chain XML file that will be used in depletion
simulation. Used to determine cross sections for materials not
present in the inital composition. Defaults to
``openmc.config['chain_file']``.
energy_bound : 2-tuple of float, optional
Bounds for the energy group.
run_kwargs : dict, optional
Keyword arguments passed to :meth:`openmc.model.Model.run`
Returns
-------
MicroXS
Cross section data in [b]
"""
# Save any original tallies on the model
original_tallies = model.tallies
# Determine what reactions and nuclides are available in chain
if chain_file is None:
chain_file = openmc.config.get('chain_file')
if chain_file is None:
raise DataError(
"No depletion chain specified and could not find depletion "
"chain in openmc.config['chain_file']"
)
chain = Chain.from_xml(chain_file)
if reactions is None:
reactions = chain.reactions
if not nuclides:
cross_sections = _find_cross_sections(model)
nuclides_with_data = _get_nuclides_with_data(cross_sections)
nuclides = [nuc.name for nuc in chain.nuclides
if nuc.name in nuclides_with_data]
# Set up the reaction rate and flux tallies
energy_filter = openmc.EnergyFilter(energy_bounds)
if isinstance(domain, openmc.Material):
domain_filter = openmc.MaterialFilter([domain])
elif isinstance(domain, openmc.Cell):
domain_filter = openmc.CellFilter([domain])
elif isinstance(domain, openmc.Universe):
domain_filter = openmc.UniverseFilter([domain])
else:
raise ValueError(f"Unsupported domain type: {type(domain)}")
rr_tally = openmc.Tally(name='MicroXS RR')
rr_tally.filters = [domain_filter, energy_filter]
rr_tally.nuclides = nuclides
rr_tally.multiply_density = False
rr_tally.scores = reactions
flux_tally = openmc.Tally(name='MicroXS flux')
flux_tally.filters = [domain_filter, energy_filter]
flux_tally.scores = ['flux']
tallies = openmc.Tallies([rr_tally, flux_tally])
model.tallies = tallies
# create temporary run
with tempfile.TemporaryDirectory() as temp_dir:
if run_kwargs is None:
run_kwargs = {}
run_kwargs.setdefault('cwd', temp_dir)
statepoint_path = model.run(**run_kwargs)
with StatePoint(statepoint_path) as sp:
rr_tally = sp.tallies[rr_tally.id]
rr_tally._read_results()
flux_tally = sp.tallies[flux_tally.id]
flux_tally._read_results()
# Get reaction rates and flux values
reaction_rates = rr_tally.mean.sum(axis=0) # (nuclides, reactions)
flux = flux_tally.mean[0, 0, 0]
# Divide RR by flux to get microscopic cross sections
xs = reaction_rates / flux
# Build Series objects
series = {}
for i, rx in enumerate(reactions):
series[rx] = Series(xs[..., i], index=rr_tally.nuclides)
# Revert to the original tallies and materials
model.tallies = original_tallies
return cls(series).rename_axis('nuclide')
@classmethod
def from_array(cls, nuclides, reactions, data):
"""
Creates a ``MicroXS`` object from arrays.
Parameters
----------
nuclides : list of str
List of nuclide symbols for that have data for at least one
reaction.
reactions : list of str
List of reactions. All reactions must match those in
:data:`openmc.deplete.chain.REACTIONS`
data : ndarray of floats
Array containing one-group microscopic cross section values for
each nuclide and reaction. Cross section values are assumed to be
in [b].
Returns
-------
MicroXS
"""
def __init__(self, data: np.ndarray, nuclides: List[str], reactions: List[str]):
# Validate inputs
if data.shape != (len(nuclides), len(reactions)):
if data.shape[:2] != (len(nuclides), len(reactions)):
raise ValueError(
f'Nuclides list of length {len(nuclides)} and '
f'reactions array of length {len(reactions)} do not '
f'match dimensions of data array of shape {data.shape}')
check_iterable_type('nuclides', nuclides, str)
check_iterable_type('reactions', reactions, str)
check_type('data', data, np.ndarray, expected_iter_type=float)
for reaction in reactions:
check_value('reactions', reaction, _valid_rxns)
cls._validate_micro_xs_inputs(
nuclides, reactions, data)
micro_xs = cls(index=nuclides, columns=reactions, data=data)
self.data = data
self.nuclides = nuclides
self.reactions = reactions
return micro_xs
# TODO: Add a classmethod for generating MicroXS directly from cross section
# data using a known flux spectrum
@classmethod
def from_csv(cls, csv_file, **kwargs):
"""
Load a ``MicroXS`` object from a ``.csv`` file.
"""Load data from a comma-separated values (csv) file.
Parameters
----------
@ -199,17 +211,37 @@ class MicroXS(DataFrame):
if 'float_precision' not in kwargs:
kwargs['float_precision'] = 'round_trip'
micro_xs = cls(read_csv(csv_file, index_col=0, **kwargs))
df = pd.read_csv(csv_file, **kwargs)
df.set_index(['nuclides', 'reactions', 'groups'], inplace=True)
nuclides = list(df.index.unique(level='nuclides'))
reactions = list(df.index.unique(level='reactions'))
groups = list(df.index.unique(level='groups'))
shape = (len(nuclides), len(reactions), len(groups))
data = df.values.reshape(shape)
return cls(data, nuclides, reactions)
cls._validate_micro_xs_inputs(list(micro_xs.index),
list(micro_xs.columns),
micro_xs.to_numpy())
return micro_xs
def __getitem__(self, index):
nuc, rx = index
i_nuc = self.nuclides.index(nuc)
i_rx = self.reactions.index(rx)
return self.data[i_nuc, i_rx]
def to_csv(self, *args, **kwargs):
"""Write data to a comma-separated values (csv) file
Parameters
----------
*args
Positional arguments passed to :meth:`pandas.DataFrame.to_csv`
**kwargs
Keyword arguments passed to :meth:`pandas.DataFrame.to_csv`
"""
groups = self.data.shape[2]
multi_index = pd.MultiIndex.from_product(
[self.nuclides, self.reactions, range(1, groups + 1)],
names=['nuclides', 'reactions', 'groups']
)
df = pd.DataFrame({'xs': self.data.flatten()}, index=multi_index)
df.to_csv(*args, **kwargs)
@staticmethod
def _validate_micro_xs_inputs(nuclides, reactions, data):
check_iterable_type('nuclides', nuclides, str)
check_iterable_type('reactions', reactions, str)
check_type('data', data, np.ndarray, expected_iter_type=float)
for reaction in reactions:
check_value('reactions', reaction, _valid_rxns)

View file

@ -6,6 +6,7 @@ transport-independent transport operators.
"""
from abc import abstractmethod
from typing import List, Tuple, Dict
import numpy as np
@ -31,9 +32,9 @@ class OpenMCOperator(TransportOperator):
----------
materials : openmc.Materials
List of all materials in the model
cross_sections : str or pandas.DataFrame
Path to continuous energy cross section library, or object containing
one-group cross-sections.
cross_sections : str or list of MicroXS
Path to continuous energy cross section library, or list of objects
containing cross sections.
chain_file : str, optional
Path to the depletion chain XML file. Defaults to
openmc.config['chain_file'].
@ -60,9 +61,9 @@ class OpenMCOperator(TransportOperator):
----------
materials : openmc.Materials
All materials present in the model
cross_sections : str or MicroXS
Path to continuous energy cross section library, or object
containing one-group cross-sections.
cross_sections : str or list of MicroXS
Path to continuous energy cross section library, or list of objects
containing cross sections.
output_dir : pathlib.Path
Path to output directory to save results.
round_number : bool
@ -164,7 +165,7 @@ class OpenMCOperator(TransportOperator):
"""Assign distribmats for each burnable material"""
pass
def _get_burnable_mats(self):
def _get_burnable_mats(self) -> Tuple[List[str], Dict[str, float], List[str]]:
"""Determine depletable materials, volumes, and nuclides
Returns

View file

@ -2,6 +2,8 @@
An ndarray to store reaction rates with string, integer, or slice indexing.
"""
from typing import Dict
import numpy as np
@ -51,6 +53,10 @@ class ReactionRates(np.ndarray):
# the __array_finalize__ method (discussed here:
# https://docs.scipy.org/doc/numpy/user/basics.subclassing.html)
index_mat: Dict[str, int]
index_nuc: Dict[str, int]
index_rx: Dict[str, int]
def __new__(cls, local_mats, nuclides, reactions, from_results=False):
# Create appropriately-sized zeroed-out ndarray
shape = (len(local_mats), len(nuclides), len(reactions))

View file

@ -1,13 +1,25 @@
nuclide,"(n,gamma)",fission
U234,22.231989822002454,0.4962074466374984
U235,10.479008971197121,48.41787337164606
U238,0.8673334105437321,0.1046788058876236
U236,8.651710446071224,0.31948392400019293
O16,7.497851000107522e-05,0.0
O17,0.0004079227797153372,0.0
I135,6.842395323713929,0.0
Xe135,227463.8642699061,0.0
Xe136,0.023178960347535887,0.0
Cs135,2.1721665580713623,0.0
Gd157,12786.099392370175,0.0
Gd156,3.4006085445846983,0.0
nuclides,reactions,groups,xs
U234,"(n,gamma)",1,22.23198982200245
U234,fission,1,0.4962074466374984
U235,"(n,gamma)",1,10.47900897119712
U235,fission,1,48.41787337164606
U238,"(n,gamma)",1,0.8673334105437321
U238,fission,1,0.1046788058876236
U236,"(n,gamma)",1,8.651710446071224
U236,fission,1,0.3194839240001929
O16,"(n,gamma)",1,7.497851000107522e-05
O16,fission,1,0.0
O17,"(n,gamma)",1,0.0004079227797153
O17,fission,1,0.0
I135,"(n,gamma)",1,6.842395323713929
I135,fission,1,0.0
Xe135,"(n,gamma)",1,227463.8642699061
Xe135,fission,1,0.0
Xe136,"(n,gamma)",1,0.0231789603475358
Xe136,fission,1,0.0
Cs135,"(n,gamma)",1,2.1721665580713623
Cs135,fission,1,0.0
Gd157,"(n,gamma)",1,12786.099392370175
Gd157,fission,1,0.0
Gd156,"(n,gamma)",1,3.4006085445846983
Gd156,fission,1,0.0

1 nuclide nuclides (n,gamma) reactions fission groups xs
2 U234 U234 22.231989822002454 (n,gamma) 0.4962074466374984 1 22.23198982200245
3 U235 U234 10.479008971197121 fission 48.41787337164606 1 0.4962074466374984
4 U238 U235 0.8673334105437321 (n,gamma) 0.1046788058876236 1 10.47900897119712
5 U236 U235 8.651710446071224 fission 0.31948392400019293 1 48.41787337164606
6 O16 U238 7.497851000107522e-05 (n,gamma) 0.0 1 0.8673334105437321
7 O17 U238 0.0004079227797153372 fission 0.0 1 0.1046788058876236
8 I135 U236 6.842395323713929 (n,gamma) 0.0 1 8.651710446071224
9 Xe135 U236 227463.8642699061 fission 0.0 1 0.3194839240001929
10 Xe136 O16 0.023178960347535887 (n,gamma) 0.0 1 7.497851000107522e-05
11 Cs135 O16 2.1721665580713623 fission 0.0 1 0.0
12 Gd157 O17 12786.099392370175 (n,gamma) 0.0 1 0.0004079227797153
13 Gd156 O17 3.4006085445846983 fission 0.0 1 0.0
14 I135 (n,gamma) 1 6.842395323713929
15 I135 fission 1 0.0
16 Xe135 (n,gamma) 1 227463.8642699061
17 Xe135 fission 1 0.0
18 Xe136 (n,gamma) 1 0.0231789603475358
19 Xe136 fission 1 0.0
20 Cs135 (n,gamma) 1 2.1721665580713623
21 Cs135 fission 1 0.0
22 Gd157 (n,gamma) 1 12786.099392370175
23 Gd157 fission 1 0.0
24 Gd156 (n,gamma) 1 3.4006085445846983
25 Gd156 fission 1 0.0

View file

@ -36,13 +36,16 @@ def chain_file():
return Path(__file__).parents[2] / 'chain_simple.xml'
@pytest.mark.parametrize("multiproc, from_nuclides, normalization_mode, power, flux", [
(True, True, 'source-rate', None, 1164719970082145.0),
(False, True, 'source-rate', None, 1164719970082145.0),
neutron_per_cm2_sec = 1164719970082145.0
@pytest.mark.parametrize("multiproc, from_nuclides, normalization_mode, power, source_rate", [
(True, True, 'source-rate', None, 1.0),
(False, True, 'source-rate', None, 1.0),
(True, True, 'fission-q', 174, None),
(False, True, 'fission-q', 174, None),
(True, False, 'source-rate', None, 1164719970082145.0),
(False, False, 'source-rate', None, 1164719970082145.0),
(True, False, 'source-rate', None, 1.0),
(False, False, 'source-rate', None, 1.0),
(True, False, 'fission-q', 174, None),
(False, False, 'fission-q', 174, None)])
def test_against_self(run_in_tmpdir,
@ -53,7 +56,7 @@ def test_against_self(run_in_tmpdir,
from_nuclides,
normalization_mode,
power,
flux):
source_rate):
"""Transport free system test suite.
Runs an OpenMC transport-free depletion calculation and verifies
@ -61,8 +64,10 @@ def test_against_self(run_in_tmpdir,
"""
# Create operator
flux = neutron_per_cm2_sec * fuel.volume
op = _create_operator(from_nuclides,
fuel,
flux,
micro_xs,
chain_file,
normalization_mode)
@ -75,12 +80,12 @@ def test_against_self(run_in_tmpdir,
openmc.deplete.PredictorIntegrator(op,
dt,
power=power,
source_rates=flux,
source_rates=source_rate,
timestep_units='s').integrate()
# Get path to test and reference results
path_test = op.output_dir / 'depletion_results.h5'
if flux is not None:
if power is None:
ref_path = 'test_reference_source_rate.h5'
else:
ref_path = 'test_reference_fission_q.h5'
@ -99,8 +104,8 @@ def test_against_self(run_in_tmpdir,
_assert_same_mats(res_test, res_ref)
tol = 1.0e-14
assert_atoms_equal(res_test, res_ref, tol)
assert_reaction_rates_equal(res_test, res_ref, tol)
assert_atoms_equal(res_ref, res_test, tol)
assert_reaction_rates_equal(res_ref, res_test, tol)
@pytest.mark.parametrize("multiproc, dt, time_units, time_type, atom_tol, rx_tol ", [
@ -123,7 +128,8 @@ def test_against_coupled(run_in_tmpdir,
atom_tol,
rx_tol):
# Create operator
op = _create_operator(False, fuel, micro_xs, chain_file, 'fission-q')
flux = neutron_per_cm2_sec * fuel.volume
op = _create_operator(False, fuel, flux, micro_xs, chain_file, 'fission-q')
# Power and timesteps
dt = [dt] # single step
@ -151,12 +157,13 @@ def test_against_coupled(run_in_tmpdir,
# Assert same mats
_assert_same_mats(res_test, res_ref)
assert_atoms_equal(res_test, res_ref, atom_tol)
assert_reaction_rates_equal(res_test, res_ref, rx_tol)
assert_atoms_equal(res_ref, res_test, atom_tol)
assert_reaction_rates_equal(res_ref, res_test, rx_tol)
def _create_operator(from_nuclides,
fuel,
flux,
micro_xs,
chain_file,
normalization_mode):
@ -167,13 +174,15 @@ def _create_operator(from_nuclides,
op = IndependentOperator.from_nuclides(fuel.volume,
nuclides,
flux,
micro_xs,
chain_file,
normalization_mode=normalization_mode)
else:
op = IndependentOperator(openmc.Materials([fuel]),
micro_xs,
[flux],
[micro_xs],
chain_file,
normalization_mode=normalization_mode)

View file

@ -4,7 +4,7 @@ from pathlib import Path
import numpy as np
import pytest
import openmc
from openmc.deplete import MicroXS
from openmc.deplete import MicroXS, get_microxs_and_flux
from tests.regression_tests import config
@ -47,13 +47,13 @@ def model():
def test_from_model(model):
fuel = model.materials[0]
domains = model.materials[:1]
nuclides = ['U234', 'U235', 'U238', 'U236', 'O16', 'O17', 'I135', 'Xe135',
'Xe136', 'Cs135', 'Gd157', 'Gd156']
test_xs = MicroXS.from_model(model, fuel, nuclides, chain_file=CHAIN_FILE)
_, test_xs = get_microxs_and_flux(model, domains, nuclides, chain_file=CHAIN_FILE)
if config['update']:
test_xs.to_csv('test_reference.csv')
test_xs[0].to_csv('test_reference.csv')
ref_xs = MicroXS.from_csv('test_reference.csv')
np.testing.assert_allclose(test_xs, ref_xs, rtol=1e-11)
np.testing.assert_allclose(test_xs[0].data, ref_xs.data, rtol=1e-11)

View file

@ -1,13 +1,25 @@
nuclide,"(n,gamma)",fission
U234,21.418670317831197,0.5014588470882195
U235,10.343944102483244,47.46718472611891
U238,0.8741166723597251,0.10829568455139126
U236,9.083486784689326,0.3325287927011428
O16,7.548646353912453e-05,0.0
O17,0.0004018486221310307,0.0
I135,6.6912565089429235,0.0
Xe135,223998.64185667288,0.0
Xe136,0.022934362666193576,0.0
Cs135,2.28453952223533,0.0
Gd157,12582.079620036275,0.0
Gd156,2.9421127515332417,0.0
nuclides,reactions,groups,xs
U234,"(n,gamma)",1,21.418670317831076
U234,fission,1,0.5014588470882162
U235,"(n,gamma)",1,10.343944102483215
U235,fission,1,47.46718472611895
U238,"(n,gamma)",1,0.8741166723597229
U238,fission,1,0.10829568455139067
U236,"(n,gamma)",1,9.08348678468935
U236,fission,1,0.3325287927011424
O16,"(n,gamma)",1,7.548646353912426e-05
O16,fission,1,0.0
O17,"(n,gamma)",1,0.00040184862213103105
O17,fission,1,0.0
I135,"(n,gamma)",1,6.691256508942912
I135,fission,1,0.0
Xe135,"(n,gamma)",1,223998.6418566729
Xe135,fission,1,0.0
Xe136,"(n,gamma)",1,0.022934362666193503
Xe136,fission,1,0.0
Cs135,"(n,gamma)",1,2.2845395222353204
Cs135,fission,1,0.0
Gd157,"(n,gamma)",1,12582.07962003624
Gd157,fission,1,0.0
Gd156,"(n,gamma)",1,2.942112751533234
Gd156,fission,1,0.0

1 nuclide nuclides (n,gamma) reactions fission groups xs
2 U234 U234 21.418670317831197 (n,gamma) 0.5014588470882195 1 21.418670317831076
3 U235 U234 10.343944102483244 fission 47.46718472611891 1 0.5014588470882162
4 U238 U235 0.8741166723597251 (n,gamma) 0.10829568455139126 1 10.343944102483215
5 U236 U235 9.083486784689326 fission 0.3325287927011428 1 47.46718472611895
6 O16 U238 7.548646353912453e-05 (n,gamma) 0.0 1 0.8741166723597229
7 O17 U238 0.0004018486221310307 fission 0.0 1 0.10829568455139067
8 I135 U236 6.6912565089429235 (n,gamma) 0.0 1 9.08348678468935
9 Xe135 U236 223998.64185667288 fission 0.0 1 0.3325287927011424
10 Xe136 O16 0.022934362666193576 (n,gamma) 0.0 1 7.548646353912426e-05
11 Cs135 O16 2.28453952223533 fission 0.0 1 0.0
12 Gd157 O17 12582.079620036275 (n,gamma) 0.0 1 0.00040184862213103105
13 Gd156 O17 2.9421127515332417 fission 0.0 1 0.0
14 I135 (n,gamma) 1 6.691256508942912
15 I135 fission 1 0.0
16 Xe135 (n,gamma) 1 223998.6418566729
17 Xe135 fission 1 0.0
18 Xe136 (n,gamma) 1 0.022934362666193503
19 Xe136 fission 1 0.0
20 Cs135 (n,gamma) 1 2.2845395222353204
21 Cs135 fission 1 0.0
22 Gd157 (n,gamma) 1 12582.07962003624
23 Gd157 fission 1 0.0
24 Gd156 (n,gamma) 1 2.942112751533234
25 Gd156 fission 1 0.0

View file

@ -1,4 +1,5 @@
import openmc.deplete
import numpy as np
import pytest
@ -16,11 +17,14 @@ def test_deplete_decay_products(run_in_tmpdir):
</depletion_chain>
""")
# Create MicroXS object with no cross sections
micro_xs = openmc.deplete.MicroXS(np.empty((0, 0)), [], [])
# Create depletion operator with no reactions
micro_xs = openmc.deplete.MicroXS()
op = openmc.deplete.IndependentOperator.from_nuclides(
volume=1.0,
nuclides={'Li5': 1.0},
flux=0.0,
micro_xs=micro_xs,
chain_file='test_chain.xml',
normalization_mode='source-rate'

View file

@ -4,14 +4,13 @@
from pathlib import Path
import pytest
from openmc.deplete import IndependentOperator, MicroXS
from openmc import Material, Materials
import numpy as np
CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml"
ONE_GROUP_XS = Path(__file__).parents[1] / "micro_xs_simple.csv"
def test_operator_init():
"""The test uses a temporary dummy chain. This file will be removed
at the end of the test, and only contains a depletion_chain node."""
@ -22,15 +21,18 @@ def test_operator_init():
'U236': 4.5724195495061115e+18,
'O16': 4.639065406771322e+22,
'O17': 1.7588724018066158e+19}
flux = 1.0
micro_xs = MicroXS.from_csv(ONE_GROUP_XS)
IndependentOperator.from_nuclides(
volume, nuclides, micro_xs, CHAIN_PATH, nuc_units='atom/cm3')
volume, nuclides, flux, micro_xs, CHAIN_PATH, nuc_units='atom/cm3')
fuel = Material(name="uo2")
fuel.add_element("U", 1, percent_type="ao", enrichment=4.25)
fuel.add_element("O", 2)
fuel.set_density("g/cc", 10.4)
fuel.depletable=True
fuel.depletable = True
fuel.volume = 1
materials = Materials([fuel])
IndependentOperator(materials, micro_xs, CHAIN_PATH)
fluxes = [1.0]
micros = [micro_xs]
IndependentOperator(materials, fluxes, micros, CHAIN_PATH)

View file

@ -43,17 +43,18 @@ def test_from_array():
[0., 0.],
[0., 0.1],
[0., 0.1]])
data.shape = (12, 2, 1)
MicroXS.from_array(nuclides, reactions, data)
MicroXS(data, nuclides, reactions)
with pytest.raises(ValueError, match=r'Nuclides list of length \d* and '
r'reactions array of length \d* do not '
r'match dimensions of data array of shape \(\d*\,d*\)'):
MicroXS.from_array(nuclides, reactions, data[:, 0])
r'match dimensions of data array of shape \(\d*\, \d*\)'):
MicroXS(data[:, 0], nuclides, reactions)
def test_csv():
ref_xs = MicroXS.from_csv(ONE_GROUP_XS)
ref_xs.to_csv('temp_xs.csv')
temp_xs = MicroXS.from_csv('temp_xs.csv')
assert np.all(ref_xs == temp_xs)
assert np.all(ref_xs.data == temp_xs.data)
remove('temp_xs.csv')