crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge branch 'develop' into reduced-absorption-xs

Browse source
This commit is contained in:
Paul Romano 2022-08-08 12:17:18 -05:00
commit c5a5a3b9a6
48 changed files with 3426 additions and 1084 deletions
Download patch file Download diff file Expand all files Collapse all files

1
.gitignore vendored
View file

@ -3,6 +3,7 @@
*.o
*.log
*.out
*.pkl
# Compiler python objects
*.pyc

84
docs/source/methods/depletion.rst
View file

@ -103,16 +103,17 @@ integrate over the entire timestep.
Our aim here is not to exhaustively describe all integration methods but rather
to give a few examples that elucidate the main considerations one must take into
account when choosing a method. Generally, there is a tradeoff between the
accuracy of the method and its computational expense. The expense is driven
almost entirely by the time to compute a transport solution, i.e., to evaluate
:math:`\mathbf{A}` for a given :math:`\mathbf{n}`. Thus, the cost of a method
scales with the number of :math:`\mathbf{A}` evaluations that are performed per
timestep. On the other hand, methods that require more evaluations generally
achieve higher accuracy. The predictor method only requires one evaluation and
its error converges as :math:`\mathcal{O}(h)`. The CE/CM method requires two
evaluations and is thus twice as expensive as the predictor method, but achieves
an error of :math:`\mathcal{O}(h^2)`. An exhaustive description of time
integration methods and their merits can be found in the `thesis of Colin Josey
accuracy of the method and its computational expense. In the case of
transport-coupled depletion, the expense is driven almost entirely by the time
to compute a transport solution, i.e., to evaluate :math:`\mathbf{A}` for a
given :math:`\mathbf{n}`. Thus, the cost of a method scales with the number of
:math:`\mathbf{A}` evaluations that are performed per timestep. On the other
hand, methods that require more evaluations generally achieve higher accuracy.
The predictor method only requires one evaluation and its error converges as
:math:`\mathcal{O}(h)`. The CE/CM method requires two evaluations and is thus
twice as expensive as the predictor method, but achieves an error of
:math:`\mathcal{O}(h^2)`. An exhaustive description of time integration methods
and their merits can be found in the `thesis of Colin Josey
<http://dspace.mit.edu/handle/1721.1/7582>`_.
OpenMC does not rely on a single time integration method but rather has several
@ -169,12 +170,14 @@ Data Considerations
In principle, solving Eq. :eq:`depletion-matrix` using CRAM is fairly simple:
just construct the burnup matrix at various times and solve a set of sparse
linear systems. However, constructing the burnup matrix itself involves not only
solving the transport equation to estimate transmutation reaction rates but also
a series of choices about what data to include. In OpenMC, the burnup matrix is
constructed based on data inside of a *depletion chain* file, which includes
fundamental data gathered from ENDF incident neutron, decay, and fission product
yield sublibraries. For each nuclide, this file includes:
linear systems. However, constructing the burnup matrix itself involves not
only solving the transport equation to estimate transmutation reaction rates
(in the case of transport-coupled depletion) or to obtain microscopic cross
sections (in the case of transport-independent depletion), but also a series of
choices about what data to include. In OpenMC, the burnup matrix is constructed
based on data inside of a *depletion chain* file, which includes fundamental
data gathered from ENDF incident neutron, decay, and fission product yield
sublibraries. For each nuclide, this file includes:
- What transmutation reactions are possible, their Q values, and their products;
- If a nuclide is not stable, what decay modes are possible, their branching
@ -185,9 +188,12 @@ yield sublibraries. For each nuclide, this file includes:
Transmutation Reactions
-----------------------
OpenMC will setup tallies in a problem based on what transmutation reactions are
available in a depletion chain file, so any arbitrary number of transmutation
reactions can be tracked. The pregenerated chain files that are available on
In transport-coupled depletion, OpenMC will setup tallies in a problem based on
what transmutation reactions are available in a depletion chain file, so any
arbitrary number of transmutation reactions can be tracked. In
transport-independent depletion, OpenMC will calculate reaction rates for every
reaction that is present in both the available cross sections and the depletion
chain file. The pregenerated chain files that are available on
https://openmc.org include the following transmutation reactions: fission, (n,\
:math:`\gamma`\ ), (n,2n), (n,3n), (n,4n), (n,p), and (n,\ :math:`\alpha`\ ).
@ -202,11 +208,12 @@ accurately model the branching of the capture reaction in Am241. This is
complicated by the fact that the branching ratio may depend on the incident
neutron energy causing capture.
OpenMC does not currently allow energy-dependent capture branching ratios.
However, the depletion chain file does allow a transmutation reaction to be
listed multiple times with different branching ratios resulting in different
products. Spectrum-averaged capture branching ratios have been computed in LWR
and SFR spectra and are available at https://openmc.org/depletion-chains.
OpenMC's transport solver does not currently allow energy-dependent capture
branching ratios. However, the depletion chain file does allow a transmutation
reaction to be listed multiple times with different branching ratios resulting
in different products. Spectrum-averaged capture branching ratios have been
computed in LWR and SFR spectra and are available at
https://openmc.org/depletion-chains.
Fission Product Yields
----------------------
@ -217,26 +224,31 @@ energies. It is an open question as to what the best way to handle this energy
dependence is. OpenMC includes three methods for treating the energy dependence
of FPY:
1. Use FPY data corresponding to a specified energy.
1. Use FPY data corresponding to a specified energy. This is used by default in
both transport-coupled and transport-independent depletion.
2. Tally fission rates above and below a specified cutoff energy. Assume that
all fissions below the cutoff energy correspond to thermal FPY data and all
fission above the cutoff energy correspond to fast FPY data.
fission above the cutoff energy correspond to fast FPY data. Only applicable
to transport-coupled depletion.
3. Compute the average energy at which fission events occur and use an effective
FPY by linearly interpolating between FPY provided at neighboring energies.
Only applicable to transport-coupled depletion.
The method can be selected through the ``fission_yield_mode`` argument to the
:class:`openmc.deplete.Operator` constructor.
The method for transport-coupled depletion can be selected through the
``fission_yield_mode`` argument to the :class:`openmc.deplete.CoupledOperator`
constructor.
Power Normalization
-------------------
The reaction rates provided OpenMC are given in units of reactions per source
particle. For depletion, it is necessary to compute an absolute reaction rate in
reactions per second. To do so, the reaction rates are normalized based on a
specified power. A complete description of how this normalization can be
performed is described in :ref:`usersguide_tally_normalization`. Here, we simply
note that the main depletion class, :class:`openmc.deplete.Operator`, allows the
user to choose one of two methods for estimating the heating rate, including:
In transport-coupled depletion, the reaction rates provided OpenMC are given in
units of reactions per source particle. For depletion, it is necessary to
compute an absolute reaction rate in reactions per second. To do so, the
reaction rates are normalized based on a specified power. A complete
description of how this normalization can be performed is described in
:ref:`usersguide_tally_normalization`. Here, we simply note that the main
depletion class, :class:`openmc.deplete.CoupledOperator`, allows the user to
choose one of two methods for estimating the heating rate, including:
1. Using fixed Q values from a depletion chain file (useful for comparisons to
other codes that use fixed Q values), or
@ -244,4 +256,4 @@ user to choose one of two methods for estimating the heating rate, including:
energy-dependent estimate of the true heating rate.
The method for normalization can be chosen through the ``normalization_mode``
argument to the :class:`openmc.deplete.Operator` class.
argument to the :class:`openmc.deplete.CoupledOperator` class.

2
docs/source/pythonapi/data.rst
View file

@ -61,11 +61,13 @@ Core Functions
atomic_mass
atomic_weight
combine_distributions
decay_constant
dose_coefficients
gnd_name
half_life
isotopes
kalbach_slope
linearize
thin
water_density

82
docs/source/pythonapi/deplete.rst
View file

@ -15,16 +15,18 @@ are:
1) A transport operator
2) A time-integration scheme
The former is responsible for executing a transport code, like OpenMC,
and retaining important information required for depletion. The most common examples
are reaction rates and power normalization data. The latter is responsible for
projecting reaction rates and compositions forward in calendar time across
some step size :math:`\Delta t`, and obtaining new compositions given a power
or power density. The :class:`Operator` is provided to handle communicating with
OpenMC. Several classes are provided that implement different time-integration
algorithms for depletion calculations, which are described in detail in Colin
Josey's thesis, `Development and analysis of high order neutron
transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
The former is responsible for calculating and retaining important information
required for depletion. The most common examples are reaction rates and power
normalization data. The latter is responsible for projecting reaction rates and
compositions forward in calendar time across some step size :math:`\Delta t`,
and obtaining new compositions given a power or power density. The
:class:`CoupledOperator` class is provided to obtain reaction rates via tallies
through OpenMC's transport solver, and the :class:`IndependentOperator` class is
provided to obtain reaction rates from cross-section data. Several classes are
provided that implement different time-integration algorithms for depletion
calculations, which are described in detail in Colin Josey's thesis,
`Development and analysis of high order neutron transport-depletion coupling
algorithms <http://hdl.handle.net/1721.1/113721>`_.
.. autosummary::
:toctree: generated
@ -40,18 +42,20 @@ transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
SICELIIntegrator
SILEQIIntegrator
Each of these classes expects a "transport operator" to be passed. An operator
specific to OpenMC is available using the following class:
Each of these classes expects a "transport operator" to be passed. OpenMC
provides The following classes implementing transpor operators:
.. autosummary::
:toctree: generated
:nosignatures:
:template: mycallable.rst
Operator
CoupledOperator
IndependentOperator
The :class:`Operator` must also have some knowledge of how nuclides transmute
and decay. This is handled by the :class:`Chain`.
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`.
Minimal Example
---------------
@ -64,11 +68,12 @@ A minimal example for performing depletion would be:
>>> import openmc.deplete
>>> geometry = openmc.Geometry.from_xml()
>>> settings = openmc.Settings.from_xml()
>>> model = openmc.model.Model(geometry, settings)
# Representation of a depletion chain
>>> chain_file = "chain_casl.xml"
>>> operator = openmc.deplete.Operator(
... geometry, settings, chain_file)
>>> operator = openmc.deplete.CoupledOperator(
... model, chain_file)
# Set up 5 time steps of one day each
>>> dt = [24 * 60 * 60] * 5
@ -131,6 +136,7 @@ data, such as number densities and reaction rates for each material.
:template: myclass.rst
AtomNumber
MicroXS
OperatorResult
ReactionRates
Results
@ -172,7 +178,7 @@ with :func:`cram.CRAM48` being the default.
:class:`multiprocessing.pool.Pool` class. If set to ``None`` (default), the
number returned by :func:`os.cpu_count` is used.
The following classes are used to help the :class:`openmc.deplete.Operator`
The following classes are used to help the :class:`openmc.deplete.CoupledOperator`
compute quantities like effective fission yields, reaction rates, and
total system energy.
@ -189,14 +195,45 @@ total system energy.
helpers.FissionYieldCutoffHelper
helpers.FluxCollapseHelper
The :class:`openmc.deplete.IndependentOperator` uses inner classes subclassed
from those listed above to perform similar calculations.
Intermediate Classes
--------------------
Specific implementations of abstract base classes may utilize some of
the same methods and data structures. These methods and data are stored
in intermediate classes.
Methods common to tally-based implementation of :class:`FissionYieldHelper`
are stored in :class:`helpers.TalliedFissionYieldHelper`
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
helpers.TalliedFissionYieldHelper
Methods common to OpenMC-specific implementations of :class:`TransportOperator`
are stored in :class:`openmc_operator.OpenMCOperator`
.. autosummary::
:toctree: generated
:nosignatures:
:template: mycallable.rst
openmc_operator.OpenMCOperator
Abstract Base Classes
---------------------
A good starting point for extending capabilities in :mod:`openmc.deplete` is
to examine the following abstract base classes. Custom classes can
inherit from :class:`abc.TransportOperator` to implement alternative
schemes for collecting reaction rates and other data from a transport code
prior to depleting materials
schemes for collecting reaction rates and other data prior to depleting
materials
.. autosummary::
:toctree: generated
@ -206,7 +243,9 @@ prior to depleting materials
abc.TransportOperator
The following classes are abstract classes used to pass information from
OpenMC simulations back on to the :class:`abc.TransportOperator`
transport simulations (in the case of transport-coupled depletion) or to
simply calculate these quantities directly (in the case of
transport-independent depletion) back on to the :class:`abc.TransportOperator`
.. autosummary::
:toctree: generated
@ -216,7 +255,6 @@ OpenMC simulations back on to the :class:`abc.TransportOperator`
abc.NormalizationHelper
abc.FissionYieldHelper
abc.ReactionRateHelper
abc.TalliedFissionYieldHelper
Custom integrators or depletion solvers can be developed by subclassing from
the following abstract base classes:

287
docs/source/usersguide/depletion.rst
View file

@ -4,50 +4,30 @@
Depletion and Transmutation
===========================
OpenMC supports coupled depletion, or burnup, calculations through the
:mod:`openmc.deplete` Python module. OpenMC solves the transport equation to
obtain transmutation reaction rates, and then the reaction rates are used to
solve a set of transmutation equations that determine the evolution of nuclide
densities within a material. The nuclide densities predicted as some future time
are then used to determine updated reaction rates, and the process is repeated
for as many timesteps as are requested.
OpenMC supports transport-coupled and transport-independent depletion, or
burnup, calculations through the :mod:`openmc.deplete` Python module. OpenMC
uses transmutation reaction rates to solve a set of transmutation equations
that determine the evolution of nuclide densities within a material. The
nuclide densities predicted at some future time are then used to determine
updated reaction rates, and the process is repeated for as many timesteps as
are requested.
The depletion module is designed such that the flux/reaction rate solution (the
The depletion module is designed such that the reaction rate solution (the
transport "operator") is completely isolated from the solution of the
transmutation equations and the method used for advancing time. At present, the
:mod:`openmc.deplete` module offers a single transport operator,
:class:`openmc.deplete.Operator` (which uses the OpenMC transport solver), but
in principle additional operator classes based on other transport codes could be
implemented and no changes to the depletion solver itself would be needed. The
operator class requires a :class:`openmc.Geometry` instance and a
:class:`openmc.Settings` instance::
geom = openmc.Geometry()
settings = openmc.Settings()
...
op = openmc.deplete.Operator(geom, settings)
Any material that contains a fissionable nuclide is depleted by default, but
this can behavior can be changed with the :attr:`Material.depletable` attribute.
.. important:: The volume must be specified for each material that is depleted by
setting the :attr:`Material.volume` attribute. This is necessary
in order to calculate the proper normalization of tally results
based on the source rate.
transmutation equations and the method used for advancing time.
:mod:`openmc.deplete` supports multiple time-integration methods for determining
material compositions over time. Each method appears as a different class.
For example, :class:`openmc.deplete.CECMIntegrator` runs a depletion calculation
using the CE/CM algorithm (deplete over a timestep using the middle-of-step
reaction rates). An instance of :class:`openmc.deplete.Operator` is passed to
one of these functions along with the timesteps and power level::
reaction rates). An instance of :class:`~openmc.deplete.abc.TransportOperator`
is passed to one of these functions along with the timesteps and power level::
power = 1200.0e6 # watts
timesteps = [10.0, 10.0, 10.0] # days
openmc.deplete.CECMIntegrator(op, timesteps, power, timestep_units='d').integrate()
The coupled transport-depletion problem is executed, and once it is done a
The depletion problem is executed, and once it is done a
``depletion_results.h5`` file is written. The results can be analyzed using the
:class:`openmc.deplete.Results` class. This class has methods that allow for
easy retrieval of k-effective, nuclide concentrations, and reaction rates over
@ -56,11 +36,41 @@ time::
results = openmc.deplete.Results("depletion_results.h5")
time, keff = results.get_keff()
Note that the coupling between the transport solver and the transmutation solver
happens in-memory rather than by reading/writing files on disk.
Note that the coupling between the reaction rate solver and the transmutation
solver happens in-memory rather than by reading/writing files on disk. OpenMC
has two categories of transport operators for obtaining transmutation reaction
rates.
.. _coupled-depletion:
Transport-coupled depletion
===========================
This category of operator solves the transport equation to obtain transmutation
reaction rates. At present, the :mod:`openmc.deplete` module offers a single
transport-coupled operator, :class:`openmc.deplete.CoupledOperator` (which uses
the OpenMC transport solver), but in principle additional transport-coupled
operator classes based on other transport codes could be implemented and no
changes to the depletion solver itself would be needed. The
:class:`openmc.deplete.CoupledOperator` class requires a :class:`~openmc.Model`
instance containing material, geometry, and settings information::
model = openmc.Model()
...
op = openmc.deplete.CoupledOperator(model)
Any material that contains a fissionable nuclide is depleted by default, but
this can behavior can be changed with the :attr:`Material.depletable` attribute.
.. important::
The volume must be specified for each material that is depleted by setting
the :attr:`Material.volume` attribute. This is necessary in order to
calculate the proper normalization of tally results based on the source rate.
Fixed-Source Transmutation
==========================
--------------------------
When the ``power`` or ``power_density`` argument is used for one of the
Integrator classes, it is assumed that OpenMC is running in k-eigenvalue mode,
@ -77,11 +87,11 @@ using the :attr:`Material.depletable` attribute::
mat = openmc.Material()
mat.depletable = True
When constructing the :class:`~openmc.deplete.Operator`, you should indicate
that normalization of tally results will be done based on the source rate rather
than a power or power density::
When constructing the :class:`~openmc.deplete.CoupledOperator`, you should
indicate that normalization of tally results will be done based on the source
rate rather than a power or power density::
op = openmc.deplete.Operator(geometry, settings, normalization_mode='source-rate')
op = openmc.deplete.CoupledOperator(model, normalization_mode='source-rate')
Finally, when creating a depletion integrator, use the ``source_rates`` argument::
@ -92,19 +102,22 @@ timestep in the calculation. A zero source rate for a given timestep will result
in a decay-only step, where all reaction rates are zero.
Caveats
=======
-------
.. _energy-deposition:
Energy Deposition
-----------------
~~~~~~~~~~~~~~~~~
The default energy deposition mode, ``"fission-q"``, instructs the
:class:`openmc.deplete.Operator` to normalize reaction rates using the product
of fission reaction rates and fission Q values taken from the depletion chain.
This approach does not consider indirect contributions to energy deposition,
such as neutron heating and energy from secondary photons. In doing this, the
energy deposited during a transport calculation will be lower than expected.
This causes the reaction rates to be over-adjusted to hit the user-specific
power, or power density, leading to an over-depletion of burnable materials.
:class:`~openmc.deplete.CoupledOperator` to normalize reaction rates using the
product of fission reaction rates and fission Q values taken from the depletion
chain. This approach does not consider indirect contributions to energy
deposition, such as neutron heating and energy from secondary photons. In doing
this, the energy deposited during a transport calculation will be lower than
expected. This causes the reaction rates to be over-adjusted to hit the
user-specific power, or power density, leading to an over-depletion of burnable
materials.
There are some remedies. First, the fission Q values can be directly set in a
variety of ways. This requires knowing what the total fission energy release
@ -113,29 +126,32 @@ should be, including indirect components. Some examples are provided below::
# use a dictionary of fission_q values
fission_q = {"U235": 202e+6} # energy in eV
# create a Model object
model = openmc.Model(geometry, settings)
# create a modified chain and write it to a new file
chain = openmc.deplete.Chain.from_xml("chain.xml", fission_q)
chain.export_to_xml("chain_mod_q.xml")
op = openmc.deplete.Operator(geometry, setting, "chain_mod_q.xml")
op = openmc.deplete.CoupledOperator(model, "chain_mod_q.xml")
# alternatively, pass the modified fission Q directly to the operator
op = openmc.deplete.Operator(geometry, setting, "chain.xml",
op = openmc.deplete.CoupledOperator(model, "chain.xml",
fission_q=fission_q)
A more complete way to model the energy deposition is to use the modified
heating reactions described in :ref:`methods_heating`. These values can be used
heating reactions described in :ref:`methods_heating`. These values can be used
to normalize reaction rates instead of using the fission reaction rates with::
op = openmc.deplete.Operator(geometry, settings, "chain.xml",
op = openmc.deplete.CoupledOperator(model, "chain.xml",
normalization_mode="energy-deposition")
These modified heating libraries can be generated by running the latest version
of :meth:`openmc.data.IncidentNeutron.from_njoy`, and will eventually be bundled
of :meth:`openmc.data.IncidentNeutron.from_njoy()`, and will eventually be bundled
into the distributed libraries.
Local Spectra and Repeated Materials
------------------------------------
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
It is not uncommon to explicitly create a single burnable material across many
locations. From a pure transport perspective, there is nothing wrong with
@ -160,7 +176,7 @@ the next transport step.
This can be countered by instructing the operator to treat repeated instances
of the same material as a unique material definition with::
op = openmc.deplete.Operator(geometry, settings, chain_file,
op = openmc.deplete.CoupledOperator(model, chain_file,
diff_burnable_mats=True)
For our example problem, this would deplete fuel on the outer region of the
@ -177,3 +193,164 @@ across all material instances.
This will increase the total memory usage and run time due to an increased
number of tallies and material definitions.
Transport-independent depletion
===============================
.. warning::
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 pre-calculated one-group microscopic cross
sections to obtain transmutation reaction rates. OpenMC provides the
:class:`~openmc.deplete.IndependentOperator` for this method of calculation.
While the one-group microscopic cross sections can be calculated using a
transport solver, :class:`~openmc.deplete.IndependentOperator` is not directly
coupled to any transport solver. The
:class:`~openmc.deplete.IndependentOperator` class requires a
:class:`openmc.Materials` object, a :class:`~openmc.deplete.MicroXS` object,
and a path to a depletion chain file::
# load in the microscopic cross sections
materials = openmc.Materials()
...
micro_xs = openmc.deplete.MicroXS.from_csv(micro_xs_path)
op = openmc.deplete.IndependentOperator(materials, micro_xs, chain_file)
.. 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.
An alternate constructor,
:meth:`~openmc.deplete.IndependentOperator.from_nuclides`, accepts a volume and
dictionary of nuclide concentrations in place of the :class:`openmc.Materials`
object::
nuclides = {'U234': 8.92e18,
'U235': 9.98e20,
'U238': 2.22e22,
'U236': 4.57e18,
'O16': 4.64e22,
'O17': 1.76e19}
volume = 0.5
op = openmc.deplete.IndependentOperator.from_nuclides(volume,
nuclides,
micro_xs,
chain_file,
nuc_units='atom/cm3')
A user can then define an integrator class as they would for a coupled
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 a nuclide that has a reaction
associated with it in the depletion chain is present in the `nuclides`
parameter but not the cross section data, that reaction will not be
simulated.
Generating Microscopic Cross Sections
-------------------------------------
Users can generate the one-group microscopic cross sections needed by
:class:`~openmc.deplete.IndependentOperator` using the
:class:`~openmc.deplete.MicroXS` class::
import openmc
model = openmc.Model.from_xml()
micro_xs = openmc.deplete.MicroXS.from_model(model,
model.materials[0],
chain_file)
The :meth:`~openmc.deplete.MicroXS.from_model()` method will produce a
:class:`~openmc.deplete.MicroXS` object with microscopic cross section data in
units of barns, which is what :class:`~openmc.deplete.IndependentOperator`
expects the units to be. The :class:`~openmc.deplete.MicroXS` class also
includes functions to read in cross section data directly from a ``.csv`` file
or from data arrays::
micro_xs = MicroXS.from_csv(micro_xs_path)
nuclides = ['U234', 'U235', 'U238']
reactions = ['fission', '(n,gamma)']
data = np.array([[0.1, 0.2],
[0.3, 0.4],
[0.01, 0.5]])
micro_xs = MicroXS.from_array(nuclides, reactions, data)
.. 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
:class:`~openmc.deplete.IndependentOperator` object.
Caveats
-------
Reaction Rate Normalization
~~~~~~~~~~~~~~~~~~~~~~~~~~~
The :class:`~openmc.deplete.IndependentOperator` class supports two methods for
normalizing reaction rates:
.. important::
Make sure you set the correct parameter in the :class:`openmc.abc.Integrator`
class. Use the ``source_rates`` parameter when
``normalization_mode == source-rate``, and use ``power`` or ``power_density``
when ``normalization_mode == fission-q``.
1. ``source-rate`` normalization, which assumes the ``source_rate`` provided by
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 general equation for the flux is
.. math::
\phi = \frac{P}{V \cdot \sum_i (Q_i \cdot \Sigma^f_i \cdot \rho_i)}
where :math:`\sum_i` is the sum over all nuclides :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.
.. warning::
The accuracy of results when using ``fission-q`` is entirely dependent on
your depletion chain. Make sure it has sufficient data to resolve the
dynamics of your particular scenario.
Multiple Materials
~~~~~~~~~~~~~~~~~~
Running a depletion simulation with multiple materials using the
``source-rate`` normalization method treats each material as completely
separate with respect to reaction rates. This can be useful for running many
different cases of a particular scenario. However, running a depletion
simulation with multiple materials using the ``fission-q`` normalization method
treats each material as part of the same "reactor" due to how ``fission-q``
normalization accumulates energy values from each material to a single value.
This behavior may change in the future.
Time integration
~~~~~~~~~~~~~~~~
The one-group microscopic cross sections passed to
:class:`openmc.deplete.IndependentOperator` are fixed values for the entire
depletion simulation. This implicit assumption may produce inaccurate results
for certain scenarios.

75
openmc/data/decay.py
View file

@ -9,7 +9,9 @@ from uncertainties import ufloat, UFloat
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from openmc.stats import Discrete, Tabular, combine_distributions
from .data import ATOMIC_SYMBOL, ATOMIC_NUMBER
from .function import INTERPOLATION_SCHEME
from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record
@ -314,6 +316,12 @@ class Decay(EqualityMixin):
'excited_state', 'mass', 'stable', 'spin', and 'parity'.
spectra : dict
Resulting radiation spectra for each radiation type.
sources : dict
Radioactive decay source distributions represented as a dictionary
mapping particle types (e.g., 'photon') to instances of
:class:`openmc.stats.Univariate`.
.. versionadded:: 0.13.1
"""
def __init__(self, ev_or_filename):
@ -329,6 +337,7 @@ class Decay(EqualityMixin):
self.modes = []
self.spectra = {}
self.average_energies = {}
self._sources = None
# Get head record
items = get_head_record(file_obj)
@ -495,3 +504,69 @@ class Decay(EqualityMixin):
"""
return cls(ev_or_filename)
@property
def sources(self):
"""Radioactive decay source distributions"""
# If property has been computed already, return it
# TODO: Replace with functools.cached_property when support is Python 3.9+
if self._sources is not None:
return self._sources
sources = {}
name = self.nuclide['name']
decay_constant = self.decay_constant.n
for particle, spectra in self.spectra.items():
# Set particle type based on 'particle' above
particle_type = {
'gamma': 'photon',
'beta-': 'electron',
'ec/beta+': 'positron',
'alpha': 'alpha',
'n': 'neutron',
'sf': 'fragment',
'p': 'proton',
'e-': 'electron',
'xray': 'photon',
'anti-neutrino': 'anti-neutrino',
'neutrino': 'neutrino',
}[particle]
if particle_type not in sources:
sources[particle_type] = []
# Create distribution for discrete
if spectra['continuous_flag'] in ('discrete', 'both'):
energies = []
intensities = []
for discrete_data in spectra['discrete']:
energies.append(discrete_data['energy'].n)
intensities.append(discrete_data['intensity'].n)
energies = np.array(energies)
intensity = spectra['discrete_normalization'].n
rates = decay_constant * intensity * np.array(intensities)
dist_discrete = Discrete(energies, rates)
sources[particle_type].append(dist_discrete)
# Create distribution for continuous
if spectra['continuous_flag'] in ('continuous', 'both'):
f = spectra['continuous']['probability']
if len(f.interpolation) > 1:
raise NotImplementedError("Multiple interpolation regions: {name}, {particle}")
interpolation = INTERPOLATION_SCHEME[f.interpolation[0]]
if interpolation not in ('histogram', 'linear-linear'):
raise NotImplementedError("Continuous spectra with {interpolation} interpolation ({name}, {particle}) not supported")
intensity = spectra['continuous_normalization'].n
rates = decay_constant * intensity * f.y
dist_continuous = Tabular(f.x, rates, interpolation)
sources[particle_type].append(dist_continuous)
# Combine discrete distributions
merged_sources = {}
for particle_type, dist_list in sources.items():
merged_sources[particle_type] = combine_distributions(
dist_list, [1.0]*len(dist_list))
self._sources = merged_sources
return self._sources

300
openmc/data/kalbach_mann.py
View file

@ -5,6 +5,7 @@ from warnings import warn
import numpy as np
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from openmc.stats import Tabular, Univariate, Discrete, Mixture
from .function import Tabulated1D, INTERPOLATION_SCHEME
from .angle_energy import AngleEnergy
@ -12,6 +13,240 @@ from .data import EV_PER_MEV
from .endf import get_list_record, get_tab2_record
class _AtomicRepresentation(EqualityMixin):
"""Atomic representation of an isotope or a particle.
Parameters
----------
z : int
Number of protons (atomic number)
a : int
Number of nucleons (mass number)
Raises
------
ValueError
When the number of protons (z) declared is higher than the number
of nucleons (a)
Attributes
----------
z : int
Number of protons (atomic number)
a : int
Number of nucleons (mass number)
n : int
Number of neutrons
za : int
ZA identifier, 1000*Z + A, where Z is the atomic number and A the mass
number
"""
def __init__(self, z, a):
# Sanity checks on values
cv.check_type('z', z, Integral)
cv.check_greater_than('z', z, 0, equality=True)
cv.check_type('a', a, Integral)
cv.check_greater_than('a', a, 0, equality=True)
if z > a:
raise ValueError(f"Number of protons ({z}) must be less than or "
f"equal to number of nucleons ({a}).")
self._z = z
self._a = a
def __add__(self, other):
"""Add two _AtomicRepresentations"""
z = self.z + other.z
a = self.a + other.a
return _AtomicRepresentation(z=z, a=a)
def __sub__(self, other):
"""Substract two _AtomicRepresentations"""
z = self.z - other.z
a = self.a - other.a
return _AtomicRepresentation(z=z, a=a)
@property
def a(self):
return self._a
@property
def z(self):
return self._z
@property
def n(self):
return self.a - self.z
@property
def za(self):
return self.z * 1000 + self.a
@classmethod
def from_za(cls, za):
"""Instantiate an _AtomicRepresentation from a ZA identifier.
Parameters
----------
za : int
ZA identifier, 1000*Z + A, where Z is the atomic number and A the
mass number
Returns
-------
_AtomicRepresentation
Atomic representation of the isotope/particle
"""
z, a = divmod(za, 1000)
return cls(z, a)
def _separation_energy(compound, nucleus, particle):
"""Calculates the separation energy as defined in ENDF-6 manual
BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1
and LANG=2. This function can be used for the incident or emitted
particle of the following reaction: A + a -> C -> B + b
Parameters
----------
compound : _AtomicRepresentation
Atomic representation of the compound (C)
nucleus : _AtomicRepresentation
Atomic representation of the nucleus (A or B)
particle : _AtomicRepresentation
Atomic representation of the particle (a or b)
Returns
-------
separation_energy : float
Separation energy in MeV
"""
# Determine A, Z, and N for compound and nucleus
A_c = compound.a
Z_c = compound.z
N_c = compound.n
A_a = nucleus.a
Z_a = nucleus.z
N_a = nucleus.n
# Determine breakup energy of incident particle (ENDF-6 Formats Manual,
# Appendix H, Table 3) in MeV
za_to_breaking_energy = {
1: 0.0,
1001: 0.0,
1002: 2.224566,
1003: 8.481798,
2003: 7.718043,
2004: 28.29566
}
I_a = za_to_breaking_energy[particle.za]
# Eq. 4 in in doi:10.1103/PhysRevC.37.2350 or ENDF-6 Formats Manual section
# 6.2.3.2
return (
15.68 * (A_c - A_a) -
28.07 * ((N_c - Z_c)**2 / A_c - (N_a - Z_a)**2 / A_a) -
18.56 * (A_c**(2./3.) - A_a**(2./3.)) +
33.22 * ((N_c - Z_c)**2 / A_c**(4./3.) - (N_a - Z_a)**2 / A_a**(4./3.)) -
0.717 * (Z_c**2 / A_c**(1./3.) - Z_a**2 / A_a**(1./3.)) +
1.211 * (Z_c**2 / A_c - Z_a**2 / A_a) -
I_a
)
def kalbach_slope(energy_projectile, energy_emitted, za_projectile,
za_emitted, za_target):
"""Returns Kalbach-Mann slope from calculations.
The associated reaction is defined as:
A + a -> C -> B + b
Where:
- A is the targeted nucleus,
- a is the projectile,
- C is the compound,
- B is the residual nucleus,
- b is the emitted particle.
The Kalbach-Mann slope calculation is done as defined in ENDF-6 manual
BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1 and
LANG=2. One exception to this, is that the entrance and emission channel
energies are not calculated with the AWR number, but approximated with
the number of mass instead.
Parameters
----------
energy_projectile : float
Energy of the projectile in the laboratory system in eV
energy_emitted : float
Energy of the emitted particle in the center of mass system in eV
za_projectile : int
ZA identifier of the projectile
za_emitted : int
ZA identifier of the emitted particle
za_target : int
ZA identifier of the targeted nucleus
Raises
------
NotImplementedError
When the projectile is not a neutron
Returns
-------
slope : float
Kalbach-Mann slope given with the same format as ACE file.
"""
# TODO: develop for photons as projectile
# TODO: test for other particles than neutron
if za_projectile != 1:
raise NotImplementedError(
"Developed and tested for neutron projectile only."
)
# Special handling of elemental carbon
if za_emitted == 6000:
za_emitted = 6012
if za_target == 6000:
za_target = 6012
projectile = _AtomicRepresentation.from_za(za_projectile)
emitted = _AtomicRepresentation.from_za(za_emitted)
target = _AtomicRepresentation.from_za(za_target)
compound = projectile + target
residual = compound - emitted
# Calculate entrance and emission channel energy in MeV, defined in section
# 6.2.3.2 in the ENDF-6 Formats Manual
epsilon_a = energy_projectile * target.a / (target.a + projectile.a) / EV_PER_MEV
epsilon_b = energy_emitted * (residual.a + emitted.a) \
/ (residual.a * EV_PER_MEV)
# Calculate separation energies using Eq. 4 in doi:10.1103/PhysRevC.37.2350
# or ENDF-6 Formats Manual section 6.2.3.2
s_a = _separation_energy(compound, target, projectile)
s_b = _separation_energy(compound, residual, emitted)
# See Eq. 10 in doi:10.1103/PhysRevC.37.2350 or section 6.2.3.2 in the
# ENDF-6 Formats Manual
za_to_M = {1: 1.0, 1001: 1.0, 1002: 1.0, 2004: 0.0}
za_to_m = {1: 0.5, 1001: 1.0, 1002: 1.0, 1003: 1.0, 2003: 1.0, 2004: 2.0}
M = za_to_M[projectile.za]
m = za_to_m[emitted.za]
e_a = epsilon_a + s_a
e_b = epsilon_b + s_b
r_1 = min(e_a, 130.)
r_3 = min(e_a, 41.)
x_1 = r_1 * e_b / e_a
x_3 = r_3 * e_b / e_a
return 0.04 * x_1 + 1.8e-6 * x_1**3 + 6.7e-7 * M * m * x_3**4
class KalbachMann(AngleEnergy):
"""Kalbach-Mann distribution
@ -319,7 +554,7 @@ class KalbachMann(AngleEnergy):
n_energy_out = int(ace.xss[idx + 1])
data = ace.xss[idx + 2:idx + 2 + 5*n_energy_out].copy()
data.shape = (5, n_energy_out)
data[0,:] *= EV_PER_MEV
data[0, :] *= EV_PER_MEV
# Create continuous distribution
eout_continuous = Tabular(data[0][n_discrete_lines:],
@ -352,13 +587,28 @@ class KalbachMann(AngleEnergy):
return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a)
@classmethod
def from_endf(cls, file_obj):
"""Generate Kalbach-Mann distribution from an ENDF evaluation
def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass):
"""Generate Kalbach-Mann distribution from an ENDF evaluation.
If the projectile is a neutron, the slope is calculated when it is
not given explicitly.
Parameters
----------
file_obj : file-like object
ENDF file positioned at the start of the Kalbach-Mann distribution
za_emitted : int
ZA identifier of the emitted particle
za_target : int
ZA identifier of the target
projectile_mass : float
Mass of the projectile
Warns
-----
UserWarning
If the mass of the projectile is not equal to 1 (other than
a neutron), the slope is not calculated and set to 0 if missing.
Returns
-------
@ -374,6 +624,7 @@ class KalbachMann(AngleEnergy):
energy_out = []
precompound = []
slope = []
calculated_slope = []
for i in range(ne):
items, values = get_list_record(file_obj)
energy[i] = items[1]
@ -385,19 +636,46 @@ class KalbachMann(AngleEnergy):
values.shape = (n_energy_out, n_angle + 2)
# Outgoing energy distribution at the i-th incoming energy
eout_i = values[:,0]
eout_p_i = values[:,1]
eout_i = values[:, 0]
eout_p_i = values[:, 1]
energy_out_i = Tabular(eout_i, eout_p_i, INTERPOLATION_SCHEME[lep])
energy_out.append(energy_out_i)
# Precompound and slope factors for Kalbach-Mann
r_i = values[:,2]
# Precompound factors for Kalbach-Mann
r_i = values[:, 2]
# Slope factors for Kalbach-Mann
if n_angle == 2:
a_i = values[:,3]
a_i = values[:, 3]
calculated_slope.append(False)
else:
a_i = np.zeros_like(r_i)
# Check if the projectile is not a neutron
if not np.isclose(projectile_mass, 1.0, atol=1.0e-12, rtol=0.):
warn(
"Kalbach-Mann slope calculation is only available with "
"neutrons as projectile. Slope coefficients are set to 0."
)
a_i = np.zeros_like(r_i)
calculated_slope.append(False)
else:
# TODO: retrieve ZA of the projectile
za_projectile = 1
a_i = [kalbach_slope(energy_projectile=energy[i],
energy_emitted=e,
za_projectile=za_projectile,
za_emitted=za_emitted,
za_target=za_target)
for e in eout_i]
calculated_slope.append(True)
precompound.append(Tabulated1D(eout_i, r_i))
slope.append(Tabulated1D(eout_i, a_i))
return cls(tab2.breakpoints, tab2.interpolation, energy,
energy_out, precompound, slope)
km_distribution = cls(tab2.breakpoints, tab2.interpolation, energy,
energy_out, precompound, slope)
# List of bool to indicate slope calculation by OpenMC
km_distribution._calculated_slope = calculated_slope
return km_distribution

8
openmc/data/njoy.py
View file

@ -127,7 +127,7 @@ acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%%
1 0 1 .{ext} /
'{library}: {zsymam} at {temperature}'/
{mat} {temperature}
1 1/
1 1 {ismooth}/
/
"""
@ -248,7 +248,8 @@ def make_pendf(filename, pendf='pendf', error=0.001, stdout=False):
def make_ace(filename, temperatures=None, acer=True, xsdir=None,
output_dir=None, pendf=False, error=0.001, broadr=True,
heatr=True, gaspr=True, purr=True, evaluation=None, **kwargs):
heatr=True, gaspr=True, purr=True, evaluation=None,
smoothing=True, **kwargs):
"""Generate incident neutron ACE file from an ENDF file
File names can be passed to
@ -298,6 +299,8 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
evaluation : openmc.data.endf.Evaluation, optional
If the ENDF file contains multiple material evaluations, this argument
indicates which evaluation should be used.
smoothing : bool, optional
If the smoothing option (ACER card 6) is on (True) or off (False).
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.run`
@ -380,6 +383,7 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
# acer
if acer:
ismooth = int(smoothing)
nacer_in = nlast
for i, temperature in enumerate(temperatures):
# Extend input with an ACER run for each temperature

29
openmc/data/reaction.py
View file

@ -80,6 +80,14 @@ def _get_products(ev, mt):
mt : int
The MT value of the reaction to get products for
Raises
------
IOError
When the Kalbach-Mann systematics is used, but the product
is not defined in the 'center-of-mass' system. The breakup logic
is not implemented which can lead to this error being raised while
the definition of the product is correct.
Returns
-------
products : list of openmc.data.Product
@ -141,7 +149,26 @@ def _get_products(ev, mt):
if lang == 1:
p.distribution = [CorrelatedAngleEnergy.from_endf(file_obj)]
elif lang == 2:
p.distribution = [KalbachMann.from_endf(file_obj)]
# Products need to be described in the center-of-mass system
product_center_of_mass = False
if reference_frame == 'center-of-mass':
product_center_of_mass = True
elif reference_frame == 'light-heavy':
product_center_of_mass = (awr <= 4.0)
# TODO: 'breakup' logic not implemented
if product_center_of_mass is False:
raise IOError(
"Kalbach-Mann representation must be defined in the "
"'center-of-mass' system"
)
zat = ev.target["atomic_number"] * 1000 + ev.target["mass_number"]
projectile_mass = ev.projectile["mass"]
p.distribution = [KalbachMann.from_endf(file_obj,
za,
zat,
projectile_mass)]
elif law == 2:
# Discrete two-body scattering

5
openmc/deplete/__init__.py
View file

@ -7,7 +7,10 @@ A depletion front-end tool.
from .nuclide import *
from .chain import *
from .operator import *
from .openmc_operator import *
from .coupled_operator import *
from .independent_operator import *
from .microxs import *
from .reaction_rates import *
from .atom_number import *
from .stepresult import *

73
openmc/deplete/abc.py
View file

@ -1,7 +1,6 @@
"""function module.
"""abc module.
This module contains the Operator class, which is then passed to an integrator
to run a full depletion simulation.
This module contains Abstract Base Classes for implementing operator, integrator, depletion system solver, and operator helper classes
"""
from abc import ABC, abstractmethod
@ -29,8 +28,8 @@ from .pool import deplete
__all__ = [
"OperatorResult", "TransportOperator", "ReactionRateHelper",
"NormalizationHelper", "FissionYieldHelper",
"OperatorResult", "TransportOperator",
"ReactionRateHelper", "NormalizationHelper", "FissionYieldHelper",
"Integrator", "SIIntegrator", "DepSystemSolver", "add_params"]
@ -83,7 +82,8 @@ class TransportOperator(ABC):
operator that takes a vector of material compositions and returns an
eigenvalue and reaction rates. This abstract class sets the requirements
for such a transport operator. Users should instantiate
:class:`openmc.deplete.Operator` rather than this class.
:class:`openmc.deplete.CoupledOperator` or
:class:`openmc.deplete.IndependentOperator` rather than this class.
Parameters
----------
@ -221,9 +221,11 @@ class ReactionRateHelper(ABC):
Parameters
----------
n_nucs : int
Number of burnable nuclides tracked by :class:`openmc.deplete.Operator`
Number of burnable nuclides tracked by
:class:`openmc.deplete.abc.TransportOperator`
n_react : int
Number of reactions tracked by :class:`openmc.deplete.Operator`
Number of reactions tracked by
:class:`openmc.deplete.abc.TransportOperator`
Attributes
----------
@ -292,9 +294,9 @@ class NormalizationHelper(ABC):
"""Abstract class for obtaining normalization factor on tallies
This helper class determines how reaction rates calculated by an instance of
:class:`openmc.deplete.Operator` should be normalized for the purpose of
constructing a burnup matrix. Based on the method chosen, the power or
source rate provided by the user, and reaction rates from a
:class:`openmc.deplete.abc.TransportOperator` should be normalized for the
purpose of constructing a burnup matrix. Based on the method chosen, the
power or source rate provided by the user, and reaction rates from a
:class:`ReactionRateHelper`, this class will scale reaction rates to the
correct values.
@ -302,7 +304,7 @@ class NormalizationHelper(ABC):
----------
nuclides : list of str
All nuclides with desired reaction rates. Ordered to be
consistent with :class:`openmc.deplete.Operator`
consistent with :class:`openmc.deplete.abc.TransportOperator`
"""
@ -316,9 +318,9 @@ class NormalizationHelper(ABC):
def prepare(self, chain_nucs, rate_index):
"""Perform work needed to obtain energy produced
This method is called prior to the transport simulations
in :meth:`openmc.deplete.Operator.initial_condition`. Only used for
energy-based normalization.
This method is called prior to calculating the reaction rates
in :meth:`openmc.deplete.abc.TransportOperator.initial_condition`. Only
used for energy-based normalization.
Parameters
----------
@ -329,7 +331,7 @@ class NormalizationHelper(ABC):
`fission_rates` for :meth:`update`.
"""
def update(self, fission_rates):
def update(self, fission_rates, mat_index=None):
"""Update the normalization based on fission rates (only used for
energy-based normalization)
@ -339,6 +341,8 @@ class NormalizationHelper(ABC):
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
Material index
"""
@property
@ -429,7 +433,7 @@ class FissionYieldHelper(ABC):
def unpack():
"""Unpack tally data prior to compute fission yields.
Called after a :meth:`openmc.deplete.Operator.__call__`
Called after a :meth:`openmc.deplete.abc.TransportOperator.__call__`
routine during the normalization of reaction rates.
Not necessary for all subclasses to implement, unless tallies
@ -450,7 +454,7 @@ class FissionYieldHelper(ABC):
mat_indexes : iterable of int
Indices of tallied materials that will have their fission
yields computed by this helper. Necessary as the
:class:`openmc.deplete.Operator` that uses this helper
:class:`openmc.deplete.CoupledOperator` that uses this helper
may only burn a subset of all materials when running
in parallel mode.
"""
@ -462,14 +466,15 @@ class FissionYieldHelper(ABC):
----------
nuclides : iterable of str
Nuclides with non-zero densities from the
:class:`openmc.deplete.Operator`
:class:`openmc.deplete.abc.TransportOperator`
Returns
-------
nuclides : list of str
Union of nuclides that the :class:`openmc.deplete.Operator`
says have non-zero densities at this stage and those that
have yield data. Sorted by nuclide name
Union of nuclides that the
:class:`openmc.deplete.abc.TransportOperator` says have non-zero
densities at this stage and those that have yield data. Sorted by
nuclide name
"""
return sorted(self._chain_set & set(nuclides))
@ -483,7 +488,7 @@ class FissionYieldHelper(ABC):
Parameters
----------
operator : openmc.deplete.TransportOperator
operator : openmc.deplete.abc.TransportOperator
Operator with a depletion chain
kwargs: optional
Additional keyword arguments to be used in constuction
@ -504,7 +509,7 @@ class Integrator(ABC):
_params = r"""
Parameters
----------
operator : openmc.deplete.TransportOperator
operator : openmc.deplete.abc.TransportOperator
Operator to perform transport simulations
timesteps : iterable of float or iterable of tuple
Array of timesteps. Note that values are not cumulative. The units are
@ -522,9 +527,10 @@ class Integrator(ABC):
power_density : float or iterable of float, optional
Power density of the reactor in [W/gHM]. It is multiplied by
initial heavy metal inventory to get total power if ``power``
is not speficied.
is not specified.
source_rates : float or iterable of float, optional
Source rate in [neutron/sec] for each interval in :attr:`timesteps`
Source rate in [neutron/sec] or neutron flux in [neutron/s-cm^2] for
each interval in :attr:`timesteps`
.. versionadded:: 0.12.1
timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
@ -546,7 +552,7 @@ class Integrator(ABC):
Attributes
----------
operator : openmc.deplete.TransportOperator
operator : openmc.deplete.abc.TransportOperator
Operator to perform transport simulations
chain : openmc.deplete.Chain
Depletion chain
@ -841,8 +847,8 @@ class SIIntegrator(Integrator):
_params = r"""
Parameters
----------
operator : openmc.deplete.TransportOperator
The operator object to simulate on.
operator : openmc.deplete.abc.TransportOperator
Operator to perform transport simulations
timesteps : iterable of float or iterable of tuple
Array of timesteps. Note that values are not cumulative. The units are
specified by the `timestep_units` argument when `timesteps` is an
@ -859,9 +865,10 @@ class SIIntegrator(Integrator):
power_density : float or iterable of float, optional
Power density of the reactor in [W/gHM]. It is multiplied by
initial heavy metal inventory to get total power if ``power``
is not speficied.
is not specified.
source_rates : float or iterable of float, optional
Source rate in [neutron/sec] for each interval in :attr:`timesteps`
Source rate in [neutron/sec] or neutron flux in [neutron/s-cm^2] for
each interval in :attr:`timesteps`
.. versionadded:: 0.12.1
timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
@ -886,7 +893,7 @@ class SIIntegrator(Integrator):
Attributes
----------
operator : openmc.deplete.TransportOperator
operator : openmc.deplete.abc.TransportOperator
Operator to perform transport simulations
chain : openmc.deplete.Chain
Depletion chain
@ -1006,7 +1013,7 @@ class DepSystemSolver(ABC):
Parameters
----------
A : scipy.sparse.csr_matrix
Sparse transmutation matrix ``A[j, i]`` desribing rates at
Sparse transmutation matrix ``A[j, i]`` describing rates at
which isotope ``i`` transmutes to isotope ``j``
n0 : numpy.ndarray
Initial compositions, typically given in number of atoms in some

17
openmc/deplete/atom_number.py
View file

@ -126,6 +126,23 @@ class AtomNumber:
return [nuc for nuc, ind in self.index_nuc.items()
if ind < self.n_nuc_burn]
def get_mat_volume(self, mat):
"""Return material volume
Parameters
----------
mat : str, int, openmc.Material, or slice
Material index.
Returns
-------
float
Material volume in [cm^3]
"""
mat = self._get_mat_index(mat)
return self.volume[mat]
def get_atom_density(self, mat, nuc):
"""Return atom density of given material and nuclide

3
openmc/deplete/chain.py
View file

@ -264,7 +264,8 @@ class Chain:
requires a list of ENDF incident neutron, decay, and neutron fission product
yield sublibrary files. The depletion chain used during a depletion
simulation is indicated by either an argument to
:class:`openmc.deplete.Operator` or through the
:class:`openmc.deplete.CoupledOperator` or
:class:`openmc.deplete.IndependentOperator`, or through the
``depletion_chain`` item in the :envvar:`OPENMC_CROSS_SECTIONS`
environment variable.

546
openmc/deplete/coupled_operator.py Normal file
View file

@ -0,0 +1,546 @@
"""Transport-coupled transport operator for depletion.
This module implements a transport operator coupled to OpenMC's transport solver
so that it can be used by depletion integrators. The implementation makes use of
the Python bindings to OpenMC's C API so that reading tally results and updating
material number densities is all done in-memory instead of through the
filesystem.
"""
import copy
import os
from warnings import warn
import numpy as np
from uncertainties import ufloat
import openmc
from openmc.checkvalue import check_value
from openmc.data import DataLibrary
from openmc.exceptions import DataError
import openmc.lib
from openmc.mpi import comm
from .abc import OperatorResult
from .chain import _find_chain_file
from .openmc_operator import OpenMCOperator, _distribute
from .results import Results
from .helpers import (
DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper,
FissionYieldCutoffHelper, AveragedFissionYieldHelper, EnergyScoreHelper,
SourceRateHelper, FluxCollapseHelper)
__all__ = ["CoupledOperator", "Operator", "OperatorResult"]
def _find_cross_sections(model):
"""Determine cross sections to use for depletion
Parameters
----------
model : openmc.model.Model
Reactor model
"""
if model.materials and model.materials.cross_sections is not None:
# Prefer info from Model class if available
return model.materials.cross_sections
# otherwise fallback to environment variable
cross_sections = os.environ.get("OPENMC_CROSS_SECTIONS")
if cross_sections is None:
raise DataError(
"Cross sections were not specified in Model.materials and "
"the OPENMC_CROSS_SECTIONS environment variable is not set."
)
return cross_sections
def _get_nuclides_with_data(cross_sections):
"""Loads cross_sections.xml file to find nuclides with neutron data
Parameters
----------
cross_sections : str
Path to cross_sections.xml file
Returns
-------
nuclides : set of str
Set of nuclide names that have cross secton data
"""
nuclides = set()
data_lib = DataLibrary.from_xml(cross_sections)
for library in data_lib.libraries:
if library['type'] != 'neutron':
continue
for name in library['materials']:
if name not in nuclides:
nuclides.add(name)
return nuclides
class CoupledOperator(OpenMCOperator):
"""Transport-coupled transport operator.
Instances of this class can be used to perform transport-coupled depletion
using OpenMC's transport solver. Normally, a user needn't call methods of
this class directly. Instead, an instance of this class is passed to an
integrator class, such as :class:`openmc.deplete.CECMIntegrator`.
.. versionchanged:: 0.13.0
The geometry and settings parameters have been replaced with a
model parameter that takes a :class:`~openmc.model.Model` object
Parameters
----------
model : openmc.model.Model
OpenMC model object
chain_file : str, optional
Path to the depletion chain XML file. Defaults to the file
listed under ``depletion_chain`` in
:envvar:`OPENMC_CROSS_SECTIONS` environment variable.
prev_results : Results, optional
Results from a previous depletion calculation. If this argument is
specified, the depletion calculation will start from the latest state
in the previous results.
diff_burnable_mats : bool, optional
Whether to differentiate burnable materials with multiple instances.
Volumes are divided equally from the original material volume.
normalization_mode : {"energy-deposition", "fission-q", "source-rate"}
Indicate how tally results should be normalized. ``"energy-deposition"``
computes the total energy deposited in the system and uses the ratio of
the power to the energy produced as a normalization factor.
``"fission-q"`` uses the fission Q values from the depletion chain to
compute the total energy deposited. ``"source-rate"`` normalizes
tallies based on the source rate (for fixed source calculations).
fission_q : dict, optional
Dictionary of nuclides and their fission Q values [eV]. If not given,
values will be pulled from the ``chain_file``. Only applicable
if ``"normalization_mode" == "fission-q"``
dilute_initial : float, optional
Initial atom density [atoms/cm^3] to add for nuclides that are zero
in initial condition to ensure they exist in the decay chain.
Only done for nuclides with reaction rates.
fission_yield_mode : {"constant", "cutoff", "average"}
Key indicating what fission product yield scheme to use. The
key determines what fission energy helper is used:
* "constant": :class:`~openmc.deplete.helpers.ConstantFissionYieldHelper`
* "cutoff": :class:`~openmc.deplete.helpers.FissionYieldCutoffHelper`
* "average": :class:`~openmc.deplete.helpers.AveragedFissionYieldHelper`
The documentation on these classes describe their methodology
and differences. Default: ``"constant"``
fission_yield_opts : dict of str to option, optional
Optional arguments to pass to the helper determined by
``fission_yield_mode``. Will be passed directly on to the
helper. Passing a value of None will use the defaults for
the associated helper.
reaction_rate_mode : {"direct", "flux"}, optional
Indicate how one-group reaction rates should be calculated. The "direct"
method tallies transmutation reaction rates directly. The "flux" method
tallies a multigroup flux spectrum and then collapses one-group reaction
rates after a transport solve (with an option to tally some reaction
rates directly).
.. versionadded:: 0.12.1
reaction_rate_opts : dict, optional
Keyword arguments that are passed to the reaction rate helper class.
When ``reaction_rate_mode`` is set to "flux", energy group boundaries
can be set using the "energies" key. See the
:class:`~openmc.deplete.helpers.FluxCollapseHelper` class for all
options.
.. versionadded:: 0.12.1
reduce_chain : bool, optional
If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the
depletion chain up to ``reduce_chain_level``.
.. versionadded:: 0.12
reduce_chain_level : int, optional
Depth of the search when reducing the depletion chain. Only used
if ``reduce_chain`` evaluates to true. The default value of
``None`` implies no limit on the depth.
.. versionadded:: 0.12
Attributes
----------
model : openmc.model.Model
OpenMC model object
geometry : openmc.Geometry
OpenMC geometry object
settings : openmc.Settings
OpenMC settings object
dilute_initial : float
Initial atom density [atoms/cm^3] to add for nuclides that
are zero in initial condition to ensure they exist in the decay
chain. Only done for nuclides with reaction rates.
output_dir : pathlib.Path
Path to output directory to save results.
round_number : bool
Whether or not to round output to OpenMC to 8 digits.
Useful in testing, as OpenMC is incredibly sensitive to exact values.
number : openmc.deplete.AtomNumber
Total number of atoms in simulation.
nuclides_with_data : set of str
A set listing all unique nuclides available from cross_sections.xml.
chain : openmc.deplete.Chain
The depletion chain information necessary to form matrices and tallies.
reaction_rates : openmc.deplete.ReactionRates
Reaction rates from the last operator step.
burnable_mats : list of str
All burnable material IDs
heavy_metal : float
Initial heavy metal inventory [g]
local_mats : list of str
All burnable material IDs being managed by a single process
prev_res : Results or None
Results from a previous depletion calculation. ``None`` if no
results are to be used.
cleanup_when_done : bool
Whether to finalize and clear the shared library memory when the
depletion operation is complete. Defaults to clearing the library.
"""
_fission_helpers = {
"average": AveragedFissionYieldHelper,
"constant": ConstantFissionYieldHelper,
"cutoff": FissionYieldCutoffHelper,
}
def __init__(self, model, chain_file=None, prev_results=None,
diff_burnable_mats=False, normalization_mode="fission-q",
fission_q=None, dilute_initial=1.0e3,
fission_yield_mode="constant", fission_yield_opts=None,
reaction_rate_mode="direct", reaction_rate_opts=None,
reduce_chain=False, reduce_chain_level=None):
# check for old call to constructor
if isinstance(model, openmc.Geometry):
msg = "As of version 0.13.0 openmc.deplete.CoupledOperator " \
"requires an openmc.Model object rather than the " \
"openmc.Geometry and openmc.Settings parameters. Please use " \
"the geometry and settings objects passed here to create a " \
" model with which to generate the transport Operator."
raise TypeError(msg)
# Determine cross sections / depletion chain
cross_sections = _find_cross_sections(model)
if chain_file is None:
chain_file = _find_chain_file(cross_sections)
check_value('fission yield mode', fission_yield_mode,
self._fission_helpers.keys())
check_value('normalization mode', normalization_mode,
('energy-deposition', 'fission-q', 'source-rate'))
if normalization_mode != "fission-q":
if fission_q is not None:
warn("Fission Q dictionary will not be used")
fission_q = None
self.model = model
self.settings = model.settings
self.geometry = model.geometry
# determine set of materials in the model
if not model.materials:
model.materials = openmc.Materials(
model.geometry.get_all_materials().values()
)
self.cleanup_when_done = True
if reaction_rate_opts is None:
reaction_rate_opts = {}
if fission_yield_opts is None:
fission_yield_opts = {}
helper_kwargs = {
'reaction_rate_mode': reaction_rate_mode,
'normalization_mode': normalization_mode,
'fission_yield_mode': fission_yield_mode,
'reaction_rate_opts': reaction_rate_opts,
'fission_yield_opts': fission_yield_opts
}
super().__init__(
model.materials,
cross_sections,
chain_file,
prev_results,
diff_burnable_mats,
fission_q,
dilute_initial,
helper_kwargs,
reduce_chain,
reduce_chain_level)
def _differentiate_burnable_mats(self):
"""Assign distribmats for each burnable material"""
# Count the number of instances for each cell and material
self.geometry.determine_paths(instances_only=True)
# Extract all burnable materials which have multiple instances
distribmats = set(
[mat for mat in self.materials
if mat.depletable and mat.num_instances > 1])
for mat in distribmats:
if mat.volume is None:
raise RuntimeError("Volume not specified for depletable "
"material with ID={}.".format(mat.id))
mat.volume /= mat.num_instances
if distribmats:
# Assign distribmats to cells
for cell in self.geometry.get_all_material_cells().values():
if cell.fill in distribmats:
mat = cell.fill
cell.fill = [mat.clone()
for i in range(cell.num_instances)]
self.materials = openmc.Materials(
self.model.geometry.get_all_materials().values()
)
def _load_previous_results(self):
"""Load results from a previous depletion simulation"""
# Reload volumes into geometry
self.prev_res[-1].transfer_volumes(self.model)
# Store previous results in operator
# Distribute reaction rates according to those tracked
# on this process
if comm.size != 1:
prev_results = self.prev_res
self.prev_res = Results()
mat_indexes = _distribute(range(len(self.burnable_mats)))
for res_obj in prev_results:
new_res = res_obj.distribute(self.local_mats, mat_indexes)
self.prev_res.append(new_res)
def _get_nuclides_with_data(self, cross_sections):
"""Loads cross_sections.xml file to find nuclides with neutron data
Parameters
----------
cross_sections : str
Path to cross_sections.xml file
Returns
-------
nuclides : set of str
Set of nuclide names that have cross secton data
"""
return _get_nuclides_with_data(cross_sections)
def _get_helper_classes(self, helper_kwargs):
"""Create the ``_rate_helper``, ``_normalization_helper``, and
``_yield_helper`` objects.
Parameters
----------
helper_kwargs : dict
Keyword arguments for helper classes
"""
reaction_rate_mode = helper_kwargs['reaction_rate_mode']
normalization_mode = helper_kwargs['normalization_mode']
fission_yield_mode = helper_kwargs['fission_yield_mode']
reaction_rate_opts = helper_kwargs['reaction_rate_opts']
fission_yield_opts = helper_kwargs['fission_yield_opts']
# Get classes to assist working with tallies
if reaction_rate_mode == "direct":
self._rate_helper = DirectReactionRateHelper(
self.reaction_rates.n_nuc, self.reaction_rates.n_react)
elif reaction_rate_mode == "flux":
# Ensure energy group boundaries were specified
if 'energies' not in reaction_rate_opts:
raise ValueError(
"Energy group boundaries must be specified in the "
"reaction_rate_opts argument when reaction_rate_mode is"
"set to 'flux'.")
self._rate_helper = FluxCollapseHelper(
self.reaction_rates.n_nuc,
self.reaction_rates.n_react,
**reaction_rate_opts
)
else:
raise ValueError("Invalid reaction rate mode.")
if normalization_mode == "fission-q":
self._normalization_helper = ChainFissionHelper()
elif normalization_mode == "energy-deposition":
score = "heating" if self.settings.photon_transport else "heating-local"
self._normalization_helper = EnergyScoreHelper(score)
else:
self._normalization_helper = SourceRateHelper()
# Select and create fission yield helper
fission_helper = self._fission_helpers[fission_yield_mode]
self._yield_helper = fission_helper.from_operator(
self, **fission_yield_opts)
def initial_condition(self):
"""Performs final setup and returns initial condition.
Returns
-------
list of numpy.ndarray
Total density for initial conditions.
"""
# Create XML files
if comm.rank == 0:
self.geometry.export_to_xml()
self.settings.export_to_xml()
self._generate_materials_xml()
# Initialize OpenMC library
comm.barrier()
if not openmc.lib.is_initialized:
openmc.lib.init(intracomm=comm)
# Generate tallies in memory
materials = [openmc.lib.materials[int(i)] for i in self.burnable_mats]
return super().initial_condition(materials)
def _generate_materials_xml(self):
"""Creates materials.xml from self.number.
Due to uncertainty with how MPI interacts with OpenMC API, this
constructs the XML manually. The long term goal is to do this
through direct memory writing.
"""
# Sort nuclides according to order in AtomNumber object
nuclides = list(self.number.nuclides)
for mat in self.materials:
mat._nuclides.sort(key=lambda x: nuclides.index(x[0]))
self.materials.export_to_xml()
def __call__(self, vec, source_rate):
"""Runs a simulation.
Simulation will abort under the following circumstances:
1) No energy is computed using OpenMC tallies.
Parameters
----------
vec : list of numpy.ndarray
Total atoms to be used in function.
source_rate : float
Power in [W] or source rate in [neutron/sec]
Returns
-------
openmc.deplete.OperatorResult
Eigenvalue and reaction rates resulting from transport operator
"""
# Reset results in OpenMC
openmc.lib.reset()
self._update_materials_and_nuclides(vec)
# If the source rate is zero, return zero reaction rates without running
# a transport solve
if source_rate == 0.0:
rates = self.reaction_rates.copy()
rates.fill(0.0)
return OperatorResult(ufloat(0.0, 0.0), rates)
# Run OpenMC
openmc.lib.run()
openmc.lib.reset_timers()
# Extract results
rates = self._calculate_reaction_rates(source_rate)
# Get k and uncertainty
keff = ufloat(*openmc.lib.keff())
op_result = OperatorResult(keff, rates)
return copy.deepcopy(op_result)
def _update_materials(self):
"""Updates material compositions in OpenMC on all processes."""
for rank in range(comm.size):
number_i = comm.bcast(self.number, root=rank)
for mat in number_i.materials:
nuclides = []
densities = []
for nuc in number_i.nuclides:
if nuc in self.nuclides_with_data:
val = 1.0e-24 * number_i.get_atom_density(mat, nuc)
# If nuclide is zero, do not add to the problem.
if val > 0.0:
if self.round_number:
val_magnitude = np.floor(np.log10(val))
val_scaled = val / 10**val_magnitude
val_round = round(val_scaled, 8)
val = val_round * 10**val_magnitude
nuclides.append(nuc)
densities.append(val)
else:
# Only output warnings if values are significantly
# negative. CRAM does not guarantee positive
# values.
if val < -1.0e-21:
print(f'WARNING: nuclide {nuc} in material'
f'{mat} is negative (density = {val}'
' atom/b-cm)')
number_i[mat, nuc] = 0.0
# Update densities on C API side
mat_internal = openmc.lib.materials[int(mat)]
mat_internal.set_densities(nuclides, densities)
# TODO Update densities on the Python side, otherwise the
# summary.h5 file contains densities at the first time step
@staticmethod
def write_bos_data(step):
"""Write a state-point file with beginning of step data
Parameters
----------
step : int
Current depletion step including restarts
"""
openmc.lib.statepoint_write(
"openmc_simulation_n{}.h5".format(step),
write_source=False)
def finalize(self):
"""Finalize a depletion simulation and release resources."""
if self.cleanup_when_done:
openmc.lib.finalize()
# Retain deprecated name for the time being
def Operator(*args, **kwargs):
# warn of name change
warn(
"The Operator(...) class has been renamed and will "
"be removed in a future version of OpenMC. Use "
"CoupledOperator(...) instead.",
FutureWarning
)
return CoupledOperator(*args, **kwargs)

40
openmc/deplete/helpers.py
View file

@ -1,5 +1,5 @@
"""
Class for normalizing fission energy deposition
Classes for collecting and calculating quantities for reaction rate operators
"""
import bisect
from abc import abstractmethod
@ -68,7 +68,7 @@ class TalliedFissionYieldHelper(FissionYieldHelper):
mat_indexes : iterable of int
Indices of tallied materials that will have their fission
yields computed by this helper. Necessary as the
:class:`openmc.deplete.Operator` that uses this helper
:class:`openmc.deplete.CoupledOperator` that uses this helper
may only burn a subset of all materials when running
in parallel mode.
"""
@ -133,9 +133,11 @@ class DirectReactionRateHelper(ReactionRateHelper):
Parameters
----------
n_nucs : int
Number of burnable nuclides tracked by :class:`openmc.deplete.Operator`
Number of burnable nuclides tracked by
:class:`openmc.deplete.CoupledOperator`
n_react : int
Number of reactions tracked by :class:`openmc.deplete.Operator`
Number of reactions tracked by an instance of
:class:`openmc.deplete.CoupledOperator`
Attributes
----------
@ -217,9 +219,10 @@ class FluxCollapseHelper(ReactionRateHelper):
Parameters
----------
n_nucs : int
Number of burnable nuclides tracked by :class:`openmc.deplete.Operator`
Number of burnable nuclides tracked by
:class:`openmc.deplete.CoupledOperator`
n_react : int
Number of reactions tracked by :class:`openmc.deplete.Operator`
Number of reactions tracked by :class:`openmc.deplete.CoupledOperator`
energies : iterable of float
Energy group boundaries for flux spectrum in [eV]
reactions : iterable of str
@ -385,7 +388,7 @@ class ChainFissionHelper(EnergyNormalizationHelper):
----------
nuclides : list of str
All nuclides with desired reaction rates. Ordered to be
consistent with :class:`openmc.deplete.Operator`
consistent with :class:`openmc.deplete.CoupledOperator`
energy : float
Total energy [J/s/source neutron] produced in a transport simulation.
Updated in the material iteration with :meth:`update`.
@ -423,7 +426,7 @@ class ChainFissionHelper(EnergyNormalizationHelper):
self._fission_q_vector = fission_qs
def update(self, fission_rates):
def update(self, fission_rates, mat_index=None):
"""Update energy produced with fission rates in a material
Parameters
@ -432,6 +435,8 @@ class ChainFissionHelper(EnergyNormalizationHelper):
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
Unused
"""
self._energy += dot(fission_rates, self._fission_q_vector)
@ -552,7 +557,7 @@ class ConstantFissionYieldHelper(FissionYieldHelper):
Parameters
----------
operator : openmc.deplete.TransportOperator
operator : openmc.deplete.abc.TransportOperator
operator with a depletion chain
kwargs:
Additional keyword arguments to be used in construction
@ -610,7 +615,6 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper):
Default: 0.0253 [eV]
fast_energy : float, optional
Energy of yield data corresponding to fast yields.
Default: 500 [kev]
Attributes
----------
@ -632,10 +636,10 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper):
Array of fission rate fractions with shape
``(n_mats, 2, n_nucs)``. ``results[:, 0]``
corresponds to the fraction of all fissions
that occured below ``cutoff``. The number
that occurred below ``cutoff``. The number
of materials in the first axis corresponds
to the number of materials burned by the
:class:`openmc.deplete.Operator`
:class:`openmc.deplete.CoupledOperator`
"""
def __init__(self, chain_nuclides, n_bmats, cutoff=112.0,
@ -692,7 +696,7 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper):
Parameters
----------
operator : openmc.deplete.Operator
operator : openmc.deplete.CoupledOperator
Operator with a chain and burnable materials
kwargs:
Additional keyword arguments to be used in construction
@ -718,7 +722,7 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper):
mat_indexes : iterable of int
Indices of tallied materials that will have their fission
yields computed by this helper. Necessary as the
:class:`openmc.deplete.Operator` that uses this helper
:class:`openmc.deplete.CoupledOperator` that uses this helper
may only burn a subset of all materials when running
in parallel mode.
"""
@ -788,7 +792,7 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper):
class AveragedFissionYieldHelper(TalliedFissionYieldHelper):
r"""Class that computes fission yields based on average fission energy
Computes average energy at which fission events occured with
Computes average energy at which fission events occurred with
.. math::
@ -843,7 +847,7 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper):
mat_indexes : iterable of int
Indices of tallied materials that will have their fission
yields computed by this helper. Necessary as the
:class:`openmc.deplete.Operator` that uses this helper
:class:`openmc.deplete.CoupledOperator` that uses this helper
may only burn a subset of all materials when running
in parallel mode.
"""
@ -906,7 +910,7 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper):
Use the computed average energy of fission
events to determine fission yields. If average
energy is between two sets of yields, linearly
interpolate bewteen the two.
interpolate between the two.
Otherwise take the closet set of yields.
Parameters
@ -956,7 +960,7 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper):
Parameters
----------
operator : openmc.deplete.TransportOperator
operator : openmc.deplete.CoupledOperator
Operator with a depletion chain
kwargs :
Additional keyword arguments to be used in construction

452
openmc/deplete/independent_operator.py Normal file
View file

@ -0,0 +1,452 @@
"""Transport-independent transport operator for depletion.
This module implements a transport operator that runs independently of any
transport solver by using user-provided one-group cross sections.
"""
import copy
from collections import OrderedDict
from warnings import warn
from itertools import product
import numpy as np
from uncertainties import ufloat
import openmc
from openmc.checkvalue import check_type
from openmc.mpi import comm
from .abc import ReactionRateHelper, OperatorResult
from .openmc_operator import OpenMCOperator, _distribute
from .microxs import MicroXS
from .results import Results
from .helpers import ChainFissionHelper, ConstantFissionYieldHelper, SourceRateHelper
class IndependentOperator(OpenMCOperator):
"""Transport-independent transport operator that uses one-group cross
sections to calculate reaction rates.
Instances of this class can be used to perform depletion using one-group
cross sections and constant flux or constant power. Normally, a user needn't
call methods of this class directly. Instead, an instance of this class is
passed to an integrator class, such as
:class:`openmc.deplete.CECMIntegrator`.
Parameters
----------
materials : openmc.Materials
Materials to deplete.
micro_xs : MicroXS
One-group microscopic cross sections in [b] .
chain_file : str
Path to the depletion chain XML file.
keff : 2-tuple of float, optional
keff eigenvalue and uncertainty from transport calculation.
Default is None.
prev_results : Results, optional
Results from a previous depletion calculation.
normalization_mode : {"fission-q", "source-rate"}
Indicate how reaction rates should be calculated.
``"fission-q"`` uses the fission Q values from the depletion chain to
compute the flux based on the power. ``"source-rate"`` uses a the
source rate (assumed to be neutron flux) to calculate the
reaction rates.
fission_q : dict, optional
Dictionary of nuclides and their fission Q values [eV]. If not given,
values will be pulled from the ``chain_file``. Only applicable
if ``"normalization_mode" == "fission-q"``.
reduce_chain : bool, optional
If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the
depletion chain up to ``reduce_chain_level``.
reduce_chain_level : int, optional
Depth of the search when reducing the depletion chain. Only used
if ``reduce_chain`` evaluates to true. The default value of
``None`` implies no limit on the depth.
fission_yield_opts : dict of str to option, optional
Optional arguments to pass to the
:class:`openmc.deplete.helpers.FissionYieldHelper` object. Will be
passed directly on to the helper. Passing a value of None will use
the defaults for the associated helper.
Attributes
----------
materials : openmc.Materials
All materials present in the model
cross_sections : MicroXS
Object containing one-group cross-sections in [cm^2].
dilute_initial : float
Initial atom density [atoms/cm^3] to add for nuclides that
are zero in initial condition to ensure they exist in the decay
chain. Only done for nuclides with reaction rates.
output_dir : pathlib.Path
Path to output directory to save results.
round_number : bool
Whether or not to round output to OpenMC to 8 digits.
Useful in testing, as OpenMC is incredibly sensitive to exact values.
number : openmc.deplete.AtomNumber
Total number of atoms in simulation.
nuclides_with_data : set of str
A set listing all unique nuclides available from cross_sections.xml.
chain : openmc.deplete.Chain
The depletion chain information necessary to form matrices and tallies.
reaction_rates : openmc.deplete.ReactionRates
Reaction rates from the last operator step.
burnable_mats : list of str
All burnable material IDs
heavy_metal : float
Initial heavy metal inventory [g]
local_mats : list of str
All burnable material IDs being managed by a single process
prev_res : Results or None
Results from a previous depletion calculation. ``None`` if no
results are to be used.
"""
def __init__(self,
materials,
micro_xs,
chain_file,
keff=None,
normalization_mode='source-rate',
fission_q=None,
prev_results=None,
reduce_chain=False,
reduce_chain_level=None,
fission_yield_opts=None):
# Validate micro-xs parameters
check_type('materials', materials, openmc.Materials)
check_type('micro_xs', micro_xs, MicroXS)
if keff is not None:
check_type('keff', keff, tuple, float)
keff = ufloat(*keff)
self._keff = keff
if fission_yield_opts is None:
fission_yield_opts = {}
helper_kwargs = {'normalization_mode': normalization_mode,
'fission_yield_opts': fission_yield_opts}
cross_sections = micro_xs * 1e-24
super().__init__(
materials,
cross_sections,
chain_file,
prev_results,
fission_q=fission_q,
helper_kwargs=helper_kwargs,
reduce_chain=reduce_chain,
reduce_chain_level=reduce_chain_level)
@classmethod
def from_nuclides(cls, volume, nuclides,
micro_xs,
chain_file,
nuc_units='atom/b-cm',
keff=None,
normalization_mode='source-rate',
fission_q=None,
prev_results=None,
reduce_chain=False,
reduce_chain_level=None,
fission_yield_opts=None):
"""
Alternate constructor from a dictionary of nuclide concentrations
volume : float
Volume of the material being depleted in [cm^3]
nuclides : dict of str to float
Dictionary with nuclide names as keys and nuclide concentrations as
values.
micro_xs : MicroXS
One-group microscopic cross sections.
chain_file : str
Path to the depletion chain XML file.
nuc_units : {'atom/cm3', 'atom/b-cm'}
Units for nuclide concentration.
keff : 2-tuple of float, optional
keff eigenvalue and uncertainty from transport calculation.
Default is None.
normalization_mode : {"fission-q", "source-rate"}
Indicate how reaction rates should be calculated.
``"fission-q"`` uses the fission Q values from the depletion
chain to compute the flux based on the power. ``"source-rate"`` uses
the source rate (assumed to be neutron flux) to calculate the
reaction rates.
fission_q : dict, optional
Dictionary of nuclides and their fission Q values [eV]. If not
given, values will be pulled from the ``chain_file``. Only
applicable if ``"normalization_mode" == "fission-q"``.
prev_results : Results, optional
Results from a previous depletion calculation.
reduce_chain : bool, optional
If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the
depletion chain up to ``reduce_chain_level``. Default is False.
reduce_chain_level : int, optional
Depth of the search when reducing the depletion chain. Only used
if ``reduce_chain`` evaluates to true. The default value of
``None`` implies no limit on the depth.
fission_yield_opts : dict of str to option, optional
Optional arguments to pass to the
:class:`openmc.deplete.helpers.FissionYieldHelper` class. Will be
passed directly on to the helper. Passing a value of None will use
the defaults for the associated helper.
"""
check_type('nuclides', nuclides, dict, str)
materials = cls._consolidate_nuclides_to_material(nuclides, nuc_units, volume)
return cls(materials,
micro_xs,
chain_file,
keff=keff,
normalization_mode=normalization_mode,
fission_q=fission_q,
prev_results=prev_results,
reduce_chain=reduce_chain,
reduce_chain_level=reduce_chain_level,
fission_yield_opts=fission_yield_opts)
@staticmethod
def _consolidate_nuclides_to_material(nuclides, nuc_units, volume):
"""Puts nuclide list into an openmc.Materials object.
"""
openmc.reset_auto_ids()
mat = openmc.Material()
if nuc_units == 'atom/b-cm':
for nuc, conc in nuclides.items():
mat.add_nuclide(nuc, conc)
elif nuc_units == 'atom/cm3':
for nuc, conc in nuclides.items():
mat.add_nuclide(nuc, conc * 1e-24) # convert to at/b-cm
else:
raise ValueError(f"Unit '{nuc_units}' is invalid.")
mat.volume = volume
mat.depletable = True
return openmc.Materials([mat])
def _load_previous_results(self):
"""Load results from a previous depletion simulation"""
# Reload volumes into geometry
model = openmc.Model(materials=self.materials)
self.prev_res[-1].transfer_volumes(model)
self.materials = model.materials
# Store previous results in operator
# Distribute reaction rates according to those tracked
# on this process
if comm.size != 1:
prev_results = self.prev_res
self.prev_res = Results()
mat_indexes = _distribute(range(len(self.burnable_mats)))
for res_obj in prev_results:
new_res = res_obj.distribute(self.local_mats, mat_indexes)
self.prev_res.append(new_res)
def _get_nuclides_with_data(self, cross_sections):
"""Finds nuclides with cross section data"""
return set(cross_sections.index)
class _IndependentNormalizationHelper(ChainFissionHelper):
"""Class for calculating one-group flux based on a power.
flux = Power / X, where X = volume * sum_i(Q_i * fission_micro_xs_i * density_i)
Parameters
----------
op : openmc.deplete.IndependentOperator
Reference to the object encapsulating _IndependentNormalizationHelper.
We pass this so we don't have to duplicate :attr:`IndependentOperator.number`.
"""
def __init__(self, op):
rates = op.reaction_rates
self.nuc_ind_map = {ind: nuc for nuc, ind in rates.index_nuc.items()}
self._op = op
super().__init__()
def update(self, fission_rates, mat_index=None):
"""Update 'energy' produced with fission rates in a material. What
this actually calculates is the quantity X.
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
Material index
"""
volume = self._op.number.get_mat_volume(mat_index)
densities = np.empty_like(fission_rates)
for i_nuc, nuc in self.nuc_ind_map.items():
densities[i_nuc] = self._op.number.get_atom_density(mat_index, nuc)
fission_rates *= volume * densities
super().update(fission_rates)
class _IndependentRateHelper(ReactionRateHelper):
"""Class for generating one-group reaction rates with flux and
one-group cross sections.
This class does not generate tallies, and instead stores cross sections
for each nuclide and transmutation reaction relevant for a depletion
calculation. The reaction rate is calculated by multiplying the flux by the
cross sections.
Parameters
----------
op : openmc.deplete.IndependentOperator
Reference to the object encapsulate _IndependentRateHelper.
We pass this so we don't have to duplicate the :attr:`IndependentOperator.number` object.
Attributes
----------
nuc_ind_map : dict of int to str
Dictionary mapping the nuclide index to nuclide name
rxn_ind_map : dict of int to str
Dictionary mapping reaction index to reaction name
"""
def __init__(self, op):
rates = op.reaction_rates
super().__init__(rates.n_nuc, rates.n_react)
self.nuc_ind_map = {ind: nuc for nuc, ind in rates.index_nuc.items()}
self.rxn_ind_map = {ind: rxn for rxn, ind in rates.index_rx.items()}
self._op = op
def generate_tallies(self, materials, scores):
"""Unused in this case"""
pass
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
"""
self._results_cache.fill(0.0)
volume = self._op.number.get_mat_volume(mat_id)
for i_nuc, i_react in product(nuc_index, react_index):
nuc = self.nuc_ind_map[i_nuc]
rxn = self.rxn_ind_map[i_react]
density = self._op.number.get_atom_density(mat_id, nuc)
self._results_cache[i_nuc,
i_react] = self._op.cross_sections[rxn][nuc] * density * volume
return self._results_cache
def _get_helper_classes(self, helper_kwargs):
"""Get helper classes for calculating reaction rates and fission yields
Parameters
----------
helper_kwargs : dict
Keyword arguments for helper classes
"""
normalization_mode = helper_kwargs['normalization_mode']
fission_yield_opts = helper_kwargs['fission_yield_opts']
self._rate_helper = self._IndependentRateHelper(self)
if normalization_mode == "fission-q":
self._normalization_helper = self._IndependentNormalizationHelper(self)
else:
self._normalization_helper = SourceRateHelper()
# Select and create fission yield helper
fission_helper = ConstantFissionYieldHelper
self._yield_helper = fission_helper.from_operator(
self, **fission_yield_opts)
def initial_condition(self):
"""Performs final setup and returns initial condition.
Returns
-------
list of numpy.ndarray
Total density for initial conditions.
"""
# Return number density vector
return super().initial_condition(self.materials)
def __call__(self, vec, source_rate):
"""Obtain the reaction rates
Parameters
----------
vec : list of numpy.ndarray
Total atoms to be used in function.
source_rate : float
Power in [W] or flux in [neutron/cm^2-s]
Returns
-------
openmc.deplete.OperatorResult
Eigenvalue and reaction rates resulting from transport operator
"""
self._update_materials_and_nuclides(vec)
rates = self._calculate_reaction_rates(source_rate)
keff = self._keff
op_result = OperatorResult(keff, rates)
return copy.deepcopy(op_result)
def _update_materials(self):
"""Updates material compositions in OpenMC on all processes."""
for rank in range(comm.size):
number_i = comm.bcast(self.number, root=rank)
for mat in number_i.materials:
nuclides = []
densities = []
for nuc in number_i.nuclides:
if nuc in self.nuclides_with_data:
val = 1.0e-24 * number_i.get_atom_density(mat, nuc)
# If nuclide is zero, do not add to the problem.
if val > 0.0:
if self.round_number:
val_magnitude = np.floor(np.log10(val))
val_scaled = val / 10**val_magnitude
val_round = round(val_scaled, 8)
val = val_round * 10**val_magnitude
nuclides.append(nuc)
densities.append(val)
else:
# Only output warnings if values are significantly
# negative. CRAM does not guarantee positive
# values.
if val < -1.0e-21:
print(f'WARNING: nuclide {nuc} in material'
f'{mat} is negative (density = {val}'
' atom/b-cm)')
number_i[mat, nuc] = 0.0

2
openmc/deplete/integrators.py
View file

@ -103,7 +103,7 @@ class CECMIntegrator(Integrator):
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from transport simulations
"""
# deplete across first half of inteval
# deplete across first half of interval
time0, x_middle = self._timed_deplete(conc, rates, dt / 2)
res_middle = self.operator(x_middle, source_rate)

230
openmc/deplete/microxs.py Normal file
View file

@ -0,0 +1,230 @@
"""MicroXS module
A pandas.DataFrame storing microscopic cross section data with
nuclides names as row indices and reaction names as column indices.
"""
import tempfile
from pathlib import Path
from copy import deepcopy
from pandas import DataFrame, read_csv, concat
import numpy as np
from openmc.checkvalue import check_type, check_value, check_iterable_type
from openmc.mgxs import EnergyGroups, ArbitraryXS, FissionXS
from openmc.data import DataLibrary
from openmc import Tallies, StatePoint, Materials, Material
from .chain import Chain, REACTIONS
from .coupled_operator import _find_cross_sections, _get_nuclides_with_data
_valid_rxns = list(REACTIONS)
_valid_rxns.append('fission')
class MicroXS(DataFrame):
"""Stores microscopic cross section data for use in
independent depletion.
"""
@classmethod
def from_model(cls,
model,
reaction_domain,
chain_file,
dilute_initial=1.0e3,
energy_bounds=(0, 20e6)):
"""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.
reaction_domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh
Domain in which to tally reaction rates.
chain_file : str
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.
dilute_initial : float
Initial atom density [atoms/cm^3] to add for nuclides that
are zero in initial condition to ensure they exist in the cross
section data. Only done for nuclides with reaction rates.
reactions : list of str, optional
Reaction names to tally
energy_bound : 2-tuple of float, optional
Bounds for the energy group.
Returns
-------
MicroXS
Cross section data in [b]
"""
groups = EnergyGroups(energy_bounds)
# Set up the reaction tallies
original_tallies = model.tallies
original_materials = deepcopy(model.materials)
tallies = Tallies()
xs = {}
reactions, diluted_materials = cls._add_dilute_nuclides(chain_file,
model,
dilute_initial)
model.materials = diluted_materials
for rxn in reactions:
if rxn == 'fission':
xs[rxn] = FissionXS(domain=reaction_domain, groups=groups, by_nuclide=True)
else:
xs[rxn] = ArbitraryXS(rxn, domain=reaction_domain, groups=groups, by_nuclide=True)
tallies += xs[rxn].tallies.values()
model.tallies = tallies
# create temporary run
with tempfile.TemporaryDirectory() as temp_dir:
statepoint_path = model.run(cwd=temp_dir)
with StatePoint(statepoint_path) as sp:
for rxn in xs:
xs[rxn].load_from_statepoint(sp)
# Build the DataFrame
series = {}
for rx in xs:
df = xs[rx].get_pandas_dataframe(xs_type='micro')
series[rx] = df.set_index('nuclide')['mean']
# Revert to the original tallies and materials
model.tallies = original_tallies
model.materials = original_materials
return cls(series)
@classmethod
def _add_dilute_nuclides(cls, chain_file, model, dilute_initial):
"""
Add nuclides not present in burnable materials that have neutron data
and are present in the depletion chain to those materials. This allows
us to tally those specific nuclides for reactions to create one-group
cross sections.
Parameters
----------
chain_file : str
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.
model : openmc.Model
Model object
dilute_initial : float
Initial atom density [atoms/cm^3] to add for nuclides that
are zero in initial condition to ensure they exist in the cross
section data. Only done for nuclides with reaction rates.
Returns
-------
reactions : list of str
List of reaction names
diluted_materials : openmc.Materials
:class:`openmc.Materials` object with nuclides added to burnable
materials.
"""
chain = Chain.from_xml(chain_file)
reactions = chain.reactions
cross_sections = _find_cross_sections(model)
nuclides_with_data = _get_nuclides_with_data(cross_sections)
burnable_nucs = [nuc.name for nuc in chain.nuclides
if nuc.name in nuclides_with_data]
diluted_materials = Materials()
for material in model.materials:
if material.depletable:
nuc_densities = material.get_nuclide_atom_densities()
dilute_density = 1.e-24 * dilute_initial
material.set_density('sum')
for nuc, density in nuc_densities.items():
material.remove_nuclide(nuc)
material.add_nuclide(nuc, density)
for burn_nuc in burnable_nucs:
if burn_nuc not in nuc_densities:
material.add_nuclide(burn_nuc,
dilute_density)
diluted_materials.append(material)
return reactions, diluted_materials
@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
"""
# Validate inputs
if data.shape != (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}')
cls._validate_micro_xs_inputs(
nuclides, reactions, data)
micro_xs = cls(index=nuclides, columns=reactions, data=data)
return micro_xs
@classmethod
def from_csv(cls, csv_file, **kwargs):
"""
Load a ``MicroXS`` object from a ``.csv`` file.
Parameters
----------
csv_file : str
Relative path to csv-file containing microscopic cross section
data. Cross section values are assumed to be in [b]
**kwargs : dict
Keyword arguments to pass to :func:`pandas.read_csv()`.
Returns
-------
MicroXS
"""
if 'float_precision' not in kwargs:
kwargs['float_precision'] = 'round_trip'
micro_xs = cls(read_csv(csv_file, index_col=0, **kwargs))
cls._validate_micro_xs_inputs(list(micro_xs.index),
list(micro_xs.columns),
micro_xs.to_numpy())
return micro_xs
@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)

571
openmc/deplete/openmc_operator.py Normal file
View file

@ -0,0 +1,571 @@
"""OpenMC transport operator
This module implements functions shared by both OpenMC transport-coupled and
transport-independent transport operators.
"""
from abc import abstractmethod
from collections import OrderedDict
import numpy as np
import openmc
from openmc.mpi import comm
from .abc import TransportOperator, OperatorResult
from .atom_number import AtomNumber
from .reaction_rates import ReactionRates
__all__ = ["OpenMCOperator", "OperatorResult"]
def _distribute(items):
"""Distribute items across MPI communicator
Parameters
----------
items : list
List of items of distribute
Returns
-------
list
Items assigned to process that called
"""
min_size, extra = divmod(len(items), comm.size)
j = 0
for i in range(comm.size):
chunk_size = min_size + int(i < extra)
if comm.rank == i:
return items[j:j + chunk_size]
j += chunk_size
class OpenMCOperator(TransportOperator):
"""Abstract class holding OpenMC-specific functions for running
depletion calculations.
Specific classes for running transport-coupled or transport-independent
depletion calculations are implemented as subclasses of OpenMCOperator.
Parameters
----------
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.
chain_file : str, optional
Path to the depletion chain XML file. Defaults to the file
listed under ``depletion_chain`` in
:envvar:`OPENMC_CROSS_SECTIONS` environment variable.
prev_results : Results, optional
Results from a previous depletion calculation. If this argument is
specified, the depletion calculation will start from the latest state
in the previous results.
diff_burnable_mats : bool, optional
Whether to differentiate burnable materials with multiple instances.
Volumes are divided equally from the original material volume.
fission_q : dict, optional
Dictionary of nuclides and their fission Q values [eV].
dilute_initial : float, optional
Initial atom density [atoms/cm^3] to add for nuclides that are zero
in initial condition to ensure they exist in the decay chain.
Only done for nuclides with reaction rates.
helper_kwargs : dict
Keyword arguments for helper classes
reduce_chain : bool, optional
If True, use :meth:`openmc.deplete.Chain.reduce()` to reduce the
depletion chain up to ``reduce_chain_level``.
reduce_chain_level : int, optional
Depth of the search when reducing the depletion chain. Only used
if ``reduce_chain`` evaluates to true. The default value of
``None`` implies no limit on the depth.
Attributes
----------
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.
dilute_initial : float
Initial atom density [atoms/cm^3] to add for nuclides that
are zero in initial condition to ensure they exist in the decay
chain. Only done for nuclides with reaction rates.
output_dir : pathlib.Path
Path to output directory to save results.
round_number : bool
Whether or not to round output to OpenMC to 8 digits.
Useful in testing, as OpenMC is incredibly sensitive to exact values.
number : openmc.deplete.AtomNumber
Total number of atoms in simulation.
nuclides_with_data : set of str
A set listing all unique nuclides available from cross_sections.xml.
chain : openmc.deplete.Chain
The depletion chain information necessary to form matrices and tallies.
reaction_rates : openmc.deplete.ReactionRates
Reaction rates from the last operator step.
burnable_mats : list of str
All burnable material IDs
heavy_metal : float
Initial heavy metal inventory [g]
local_mats : list of str
All burnable material IDs being managed by a single process
prev_res : Results or None
Results from a previous depletion calculation. ``None`` if no
results are to be used.
"""
def __init__(
self,
materials=None,
cross_sections=None,
chain_file=None,
prev_results=None,
diff_burnable_mats=False,
fission_q=None,
dilute_initial=0.0,
helper_kwargs=None,
reduce_chain=False,
reduce_chain_level=None):
super().__init__(chain_file, fission_q, dilute_initial, prev_results)
self.round_number = False
self.materials = materials
self.cross_sections = cross_sections
# Reduce the chain to only those nuclides present
if reduce_chain:
init_nuclides = set()
for material in self.materials:
if not material.depletable:
continue
for name, _dens_percent, _dens_type in material.nuclides:
init_nuclides.add(name)
self.chain = self.chain.reduce(init_nuclides, reduce_chain_level)
if diff_burnable_mats:
self._differentiate_burnable_mats()
# Determine which nuclides have cross section data
# This nuclides variables contains every nuclides
# for which there is an entry in the micro_xs parameter
openmc.reset_auto_ids()
self.burnable_mats, volumes, all_nuclides = self._get_burnable_mats()
self.local_mats = _distribute(self.burnable_mats)
self._mat_index_map = {
lm: self.burnable_mats.index(lm) for lm in self.local_mats}
if self.prev_res is not None:
self._load_previous_results()
self.nuclides_with_data = self._get_nuclides_with_data(
self.cross_sections)
# Select nuclides with data that are also in the chain
self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides
if nuc.name in self.nuclides_with_data]
# Extract number densities from the geometry / previous depletion run
self._extract_number(self.local_mats,
volumes,
all_nuclides,
self.prev_res)
# Create reaction rates array
self.reaction_rates = ReactionRates(
self.local_mats, self._burnable_nucs, self.chain.reactions)
self._get_helper_classes(helper_kwargs)
def _differentiate_burnable_mats(self):
"""Assign distribmats for each burnable material"""
pass
def _get_burnable_mats(self):
"""Determine depletable materials, volumes, and nuclides
Returns
-------
burnable_mats : list of str
List of burnable material IDs
volume : OrderedDict of str to float
Volume of each material in [cm^3]
nuclides : list of str
Nuclides in order of how they'll appear in the simulation.
"""
burnable_mats = set()
model_nuclides = set()
volume = OrderedDict()
self.heavy_metal = 0.0
# Iterate once through the geometry to get dictionaries
for mat in self.materials:
for nuclide in mat.get_nuclides():
model_nuclides.add(nuclide)
if mat.depletable:
burnable_mats.add(str(mat.id))
if mat.volume is None:
raise RuntimeError("Volume not specified for depletable "
"material with ID={}.".format(mat.id))
volume[str(mat.id)] = mat.volume
self.heavy_metal += mat.fissionable_mass
# Make sure there are burnable materials
if not burnable_mats:
raise RuntimeError(
"No depletable materials were found in the model.")
# Sort the sets
burnable_mats = sorted(burnable_mats, key=int)
model_nuclides = sorted(model_nuclides)
# Construct a global nuclide dictionary, burned first
nuclides = list(self.chain.nuclide_dict)
for nuc in model_nuclides:
if nuc not in nuclides:
nuclides.append(nuc)
return burnable_mats, volume, nuclides
def _load_previous_results(self):
"""Load results from a previous depletion simulation"""
pass
@abstractmethod
def _get_nuclides_with_data(self, cross_sections):
"""Find nuclides with cross section data
Parameters
----------
cross_sections : str or pandas.DataFrame
Path to continuous energy cross section library, or object
containing one-group cross-sections.
Returns
-------
nuclides : set of str
Set of nuclide names that have cross secton data
"""
def _extract_number(self, local_mats, volume, all_nuclides, prev_res=None):
"""Construct AtomNumber using geometry
Parameters
----------
local_mats : list of str
Material IDs to be managed by this process
volume : OrderedDict of str to float
Volumes for the above materials in [cm^3]
all_nuclides : list of str
Nuclides to be used in the simulation.
prev_res : Results, optional
Results from a previous depletion calculation
"""
self.number = AtomNumber(local_mats, all_nuclides, volume, len(self.chain))
if self.dilute_initial != 0.0:
for nuc in self._burnable_nucs:
self.number.set_atom_density(
np.s_[:], nuc, self.dilute_initial)
# Now extract and store the number densities
# From the geometry if no previous depletion results
if prev_res is None:
for mat in self.materials:
if str(mat.id) in local_mats:
self._set_number_from_mat(mat)
# Else from previous depletion results
else:
for mat in self.materials:
if str(mat.id) in local_mats:
self._set_number_from_results(mat, prev_res)
def _set_number_from_mat(self, mat):
"""Extracts material and number densities from openmc.Material
Parameters
----------
mat : openmc.Material
The material to read from
"""
mat_id = str(mat.id)
for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items():
atom_per_cc = atom_per_bcm * 1.0e24
self.number.set_atom_density(mat_id, nuclide, atom_per_cc)
def _set_number_from_results(self, mat, prev_res):
"""Extracts material nuclides and number densities.
If the nuclide concentration's evolution is tracked, the densities come
from depletion results. Else, densities are extracted from the geometry
in the summary.
Parameters
----------
mat : openmc.Material
The material to read from
prev_res : Results
Results from a previous depletion calculation
"""
mat_id = str(mat.id)
# Get nuclide lists from geometry and depletion results
depl_nuc = prev_res[-1].nuc_to_ind
geom_nuc_densities = mat.get_nuclide_atom_densities()
# Merge lists of nuclides, with the same order for every calculation
geom_nuc_densities.update(depl_nuc)
for nuclide, atom_per_bcm in geom_nuc_densities.items():
if nuclide in depl_nuc:
concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1]
volume = prev_res[-1].volume[mat_id]
atom_per_cc = concentration / volume
else:
atom_per_cc = atom_per_bcm * 1.0e24
self.number.set_atom_density(mat_id, nuclide, atom_per_cc)
@abstractmethod
def _get_helper_classes(self, helper_kwargs):
"""Create the ``_rate_helper``, ``_normalization_helper``, and
``_yield_helper`` objects.
Parameters
----------
helper_kwargs : dict
Keyword arguments for helper classes
"""
def initial_condition(self, materials):
"""Performs final setup and returns initial condition.
Parameters
----------
materials : list of str
list of material IDs
Returns
-------
list of numpy.ndarray
Total density for initial conditions.
"""
self._rate_helper.generate_tallies(materials, self.chain.reactions)
self._normalization_helper.prepare(
self.chain.nuclides, self.reaction_rates.index_nuc)
# Tell fission yield helper what materials this process is
# responsible for
self._yield_helper.generate_tallies(
materials, tuple(sorted(self._mat_index_map.values())))
# Return number density vector
return list(self.number.get_mat_slice(np.s_[:]))
def _update_materials_and_nuclides(self, vec):
"""Update the number density, material compositions, and nuclide
lists in helper objects
Parameters
----------
vec : list of numpy.ndarray
Total atoms.
"""
# Update the number densities regardless of the source rate
self.number.set_density(vec)
self._update_materials()
# Prevent OpenMC from complaining about re-creating tallies
openmc.reset_auto_ids()
# Update tally nuclides data in preparation for transport solve
nuclides = self._get_reaction_nuclides()
self._rate_helper.nuclides = nuclides
self._normalization_helper.nuclides = nuclides
self._yield_helper.update_tally_nuclides(nuclides)
@abstractmethod
def _update_materials(self):
"""Updates material compositions in OpenMC on all processes."""
def write_bos_data(self, step):
"""Document beginning of step data for a given step
Called at the beginning of a depletion step and at
the final point in the simulation.
Parameters
----------
step : int
Current depletion step including restarts
"""
# Since we aren't running a transport simulation, we simply pass
pass
def _get_reaction_nuclides(self):
"""Determine nuclides that should be tallied for reaction rates.
This method returns a list of all nuclides that have cross section data
and are listed in the depletion chain. Technically, we should count
nuclides that may not appear in the depletion chain because we still
need to get the fission reaction rate for these nuclides in order to
normalize power, but that is left as a future exercise.
Returns
-------
list of str
Nuclides with reaction rates
"""
nuc_set = set()
# Create the set of all nuclides in the decay chain in materials marked
# for burning in which the number density is greater than zero.
for nuc in self.number.nuclides:
if nuc in self.nuclides_with_data:
if np.sum(self.number[:, nuc]) > 0.0:
nuc_set.add(nuc)
# Communicate which nuclides have nonzeros to rank 0
if comm.rank == 0:
for i in range(1, comm.size):
nuc_newset = comm.recv(source=i, tag=i)
nuc_set |= nuc_newset
else:
comm.send(nuc_set, dest=0, tag=comm.rank)
if comm.rank == 0:
# Sort nuclides in the same order as self.number
nuc_list = [nuc for nuc in self.number.nuclides
if nuc in nuc_set]
else:
nuc_list = None
# Store list of nuclides on each process
nuc_list = comm.bcast(nuc_list)
return [nuc for nuc in nuc_list if nuc in self.chain]
def _calculate_reaction_rates(self, source_rate):
"""Unpack tallies from OpenMC and return an operator result
This method uses OpenMC's C API bindings to determine the k-effective
value and reaction rates from the simulation. The reaction rates are
normalized by a helper class depending on the method being used.
Parameters
----------
source_rate : float
Power in [W] or source rate in [neutron/sec]
Returns
-------
rates : openmc.deplete.ReactionRates
Reaction rates for nuclides
"""
rates = self.reaction_rates
rates.fill(0.0)
# Extract reaction nuclides
rxn_nuclides = self._rate_helper.nuclides
# Form fast map
nuc_ind = [rates.index_nuc[nuc] for nuc in rxn_nuclides]
react_ind = [rates.index_rx[react] for react in self.chain.reactions]
# Keep track of energy produced from all reactions in eV per source
# particle
self._normalization_helper.reset()
self._yield_helper.unpack()
# Store fission yield dictionaries
fission_yields = []
# Create arrays to store fission Q values, reaction rates, and nuclide
# numbers, zeroed out in material iteration
number = np.empty(rates.n_nuc)
fission_ind = rates.index_rx.get("fission")
# Extract results
for i, mat in enumerate(self.local_mats):
# Get tally index
mat_index = self._mat_index_map[mat]
# Zero out reaction rates and nuclide numbers
number.fill(0.0)
# Get new number densities
for nuc, i_nuc_results in zip(rxn_nuclides, nuc_ind):
number[i_nuc_results] = self.number[mat, nuc]
tally_rates = self._rate_helper.get_material_rates(
mat_index, nuc_ind, react_ind)
# Compute fission yields for this material
fission_yields.append(self._yield_helper.weighted_yields(i))
# Accumulate energy from fission
if fission_ind is not None:
self._normalization_helper.update(
tally_rates[:, fission_ind],
mat_index=mat_index)
# Divide by total number and store
rates[i] = self._rate_helper.divide_by_adens(number)
# Scale reaction rates to obtain units of reactions/sec
rates *= self._normalization_helper.factor(source_rate)
# Store new fission yields on the chain
self.chain.fission_yields = fission_yields
return rates
def get_results_info(self):
"""Returns volume list, material lists, and nuc lists.
Returns
-------
volume : dict of str float
Volumes corresponding to materials in full_burn_dict
nuc_list : list of str
A list of all nuclide names. Used for sorting the simulation.
burn_list : list of int
A list of all material IDs to be burned. Used for sorting the simulation.
full_burn_list : list
List of all burnable material IDs
"""
nuc_list = self.number.burnable_nuclides
burn_list = self.local_mats
volume = {}
for i, mat in enumerate(burn_list):
volume[mat] = self.number.volume[i]
# Combine volume dictionaries across processes
volume_list = comm.allgather(volume)
volume = {k: v for d in volume_list for k, v in d.items()}
return volume, nuc_list, burn_list, self.burnable_mats

822
openmc/deplete/operator.py
View file

@ -1,822 +0,0 @@
"""OpenMC transport operator
This module implements a transport operator for OpenMC so that it can be used by
depletion integrators. The implementation makes use of the Python bindings to
OpenMC's C API so that reading tally results and updating material number
densities is all done in-memory instead of through the filesystem.
"""
import copy
from collections import OrderedDict
import os
from warnings import warn
import numpy as np
from uncertainties import ufloat
import openmc
from openmc.checkvalue import check_value
from openmc.data import DataLibrary
from openmc.exceptions import DataError
import openmc.lib
from openmc.mpi import comm
from .abc import TransportOperator, OperatorResult
from .atom_number import AtomNumber
from .chain import _find_chain_file
from .reaction_rates import ReactionRates
from .results import Results
from .helpers import (
DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper,
FissionYieldCutoffHelper, AveragedFissionYieldHelper, EnergyScoreHelper,
SourceRateHelper, FluxCollapseHelper)
__all__ = ["Operator", "OperatorResult"]
def _distribute(items):
"""Distribute items across MPI communicator
Parameters
----------
items : list
List of items of distribute
Returns
-------
list
Items assigned to process that called
"""
min_size, extra = divmod(len(items), comm.size)
j = 0
for i in range(comm.size):
chunk_size = min_size + int(i < extra)
if comm.rank == i:
return items[j:j + chunk_size]
j += chunk_size
def _find_cross_sections(model):
"""Determine cross sections to use for depletion"""
if model.materials and model.materials.cross_sections is not None:
# Prefer info from Model class if available
return model.materials.cross_sections
# otherwise fallback to environment variable
cross_sections = os.environ.get("OPENMC_CROSS_SECTIONS")
if cross_sections is None:
raise DataError(
"Cross sections were not specified in Model.materials and "
"the OPENMC_CROSS_SECTIONS environment variable is not set."
)
return cross_sections
class Operator(TransportOperator):
"""OpenMC transport operator for depletion.
Instances of this class can be used to perform depletion using OpenMC as the
transport operator. Normally, a user needn't call methods of this class
directly. Instead, an instance of this class is passed to an integrator
class, such as :class:`openmc.deplete.CECMIntegrator`.
.. versionchanged:: 0.13.0
The geometry and settings parameters have been replaced with a
model parameter that takes a :class:`~openmc.model.Model` object
Parameters
----------
model : openmc.model.Model
OpenMC model object
chain_file : str, optional
Path to the depletion chain XML file. Defaults to the file
listed under ``depletion_chain`` in
:envvar:`OPENMC_CROSS_SECTIONS` environment variable.
prev_results : Results, optional
Results from a previous depletion calculation. If this argument is
specified, the depletion calculation will start from the latest state
in the previous results.
diff_burnable_mats : bool, optional
Whether to differentiate burnable materials with multiple instances.
Volumes are divided equally from the original material volume.
Default: False.
normalization_mode : {"energy-deposition", "fission-q", "source-rate"}
Indicate how tally results should be normalized. ``"energy-deposition"``
computes the total energy deposited in the system and uses the ratio of
the power to the energy produced as a normalization factor.
``"fission-q"`` uses the fission Q values from the depletion chain to
compute the total energy deposited. ``"source-rate"`` normalizes
tallies based on the source rate (for fixed source calculations).
fission_q : dict, optional
Dictionary of nuclides and their fission Q values [eV]. If not given,
values will be pulled from the ``chain_file``. Only applicable
if ``"normalization_mode" == "fission-q"``
dilute_initial : float, optional
Initial atom density [atoms/cm^3] to add for nuclides that are zero
in initial condition to ensure they exist in the decay chain.
Only done for nuclides with reaction rates.
Defaults to 1.0e3.
fission_yield_mode : {"constant", "cutoff", "average"}
Key indicating what fission product yield scheme to use. The
key determines what fission energy helper is used:
* "constant": :class:`~openmc.deplete.helpers.ConstantFissionYieldHelper`
* "cutoff": :class:`~openmc.deplete.helpers.FissionYieldCutoffHelper`
* "average": :class:`~openmc.deplete.helpers.AveragedFissionYieldHelper`
The documentation on these classes describe their methodology
and differences. Default: ``"constant"``
fission_yield_opts : dict of str to option, optional
Optional arguments to pass to the helper determined by
``fission_yield_mode``. Will be passed directly on to the
helper. Passing a value of None will use the defaults for
the associated helper.
reaction_rate_mode : {"direct", "flux"}, optional
Indicate how one-group reaction rates should be calculated. The "direct"
method tallies transmutation reaction rates directly. The "flux" method
tallies a multigroup flux spectrum and then collapses one-group reaction
rates after a transport solve (with an option to tally some reaction
rates directly).
.. versionadded:: 0.12.1
reaction_rate_opts : dict, optional
Keyword arguments that are passed to the reaction rate helper class.
When ``reaction_rate_mode`` is set to "flux", energy group boundaries
can be set using the "energies" key. See the
:class:`~openmc.deplete.helpers.FluxCollapseHelper` class for all
options.
.. versionadded:: 0.12.1
reduce_chain : bool, optional
If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the
depletion chain up to ``reduce_chain_level``. Default is False.
.. versionadded:: 0.12
reduce_chain_level : int, optional
Depth of the search when reducing the depletion chain. Only used
if ``reduce_chain`` evaluates to true. The default value of
``None`` implies no limit on the depth.
.. versionadded:: 0.12
Attributes
----------
model : openmc.model.Model
OpenMC model object
geometry : openmc.Geometry
OpenMC geometry object
settings : openmc.Settings
OpenMC settings object
dilute_initial : float
Initial atom density [atoms/cm^3] to add for nuclides that
are zero in initial condition to ensure they exist in the decay
chain. Only done for nuclides with reaction rates.
output_dir : pathlib.Path
Path to output directory to save results.
round_number : bool
Whether or not to round output to OpenMC to 8 digits.
Useful in testing, as OpenMC is incredibly sensitive to exact values.
number : openmc.deplete.AtomNumber
Total number of atoms in simulation.
nuclides_with_data : set of str
A set listing all unique nuclides available from cross_sections.xml.
chain : openmc.deplete.Chain
The depletion chain information necessary to form matrices and tallies.
reaction_rates : openmc.deplete.ReactionRates
Reaction rates from the last operator step.
burnable_mats : list of str
All burnable material IDs
heavy_metal : float
Initial heavy metal inventory [g]
local_mats : list of str
All burnable material IDs being managed by a single process
prev_res : Results or None
Results from a previous depletion calculation. ``None`` if no
results are to be used.
cleanup_when_done : bool
Whether to finalize and clear the shared library memory when the
depletion operation is complete. Defaults to clearing the library.
"""
_fission_helpers = {
"average": AveragedFissionYieldHelper,
"constant": ConstantFissionYieldHelper,
"cutoff": FissionYieldCutoffHelper,
}
def __init__(self, model, chain_file=None, prev_results=None,
diff_burnable_mats=False, normalization_mode="fission-q",
fission_q=None, dilute_initial=1.0e3,
fission_yield_mode="constant", fission_yield_opts=None,
reaction_rate_mode="direct", reaction_rate_opts=None,
reduce_chain=False, reduce_chain_level=None):
# check for old call to constructor
if isinstance(model, openmc.Geometry):
msg = "As of version 0.13.0 openmc.deplete.Operator requires an " \
"openmc.Model object rather than the openmc.Geometry and " \
"openmc.Settings parameters. Please use the geometry and " \
"settings objects passed here to create a model with which " \
"to generate the depletion Operator."
raise TypeError(msg)
# Determine cross sections / depletion chain
cross_sections = _find_cross_sections(model)
if chain_file is None:
chain_file = _find_chain_file(cross_sections)
check_value('fission yield mode', fission_yield_mode,
self._fission_helpers.keys())
check_value('normalization mode', normalization_mode,
('energy-deposition', 'fission-q', 'source-rate'))
if normalization_mode != "fission-q":
if fission_q is not None:
warn("Fission Q dictionary will not be used")
fission_q = None
super().__init__(chain_file, fission_q, dilute_initial, prev_results)
self.round_number = False
self.model = model
self.settings = model.settings
self.geometry = model.geometry
# determine set of materials in the model
if not model.materials:
model.materials = openmc.Materials(
model.geometry.get_all_materials().values()
)
self.materials = model.materials
self.cleanup_when_done = True
# Reduce the chain before we create more materials
if reduce_chain:
all_isotopes = set()
for material in self.materials:
if not material.depletable:
continue
for name, _dens_percent, _dens_type in material.nuclides:
all_isotopes.add(name)
self.chain = self.chain.reduce(all_isotopes, reduce_chain_level)
# Differentiate burnable materials with multiple instances
if diff_burnable_mats:
self._differentiate_burnable_mats()
self.materials = openmc.Materials(
model.geometry.get_all_materials().values()
)
# Clear out OpenMC, create task lists, distribute
openmc.reset_auto_ids()
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}
if self.prev_res is not None:
# Reload volumes into geometry
prev_results[-1].transfer_volumes(self.model)
# Store previous results in operator
# Distribute reaction rates according to those tracked
# on this process
if comm.size == 1:
self.prev_res = prev_results
else:
self.prev_res = Results()
mat_indexes = _distribute(range(len(self.burnable_mats)))
for res_obj in prev_results:
new_res = res_obj.distribute(self.local_mats, mat_indexes)
self.prev_res.append(new_res)
# Determine which nuclides have incident neutron data
self.nuclides_with_data = self._get_nuclides_with_data(cross_sections)
# Select nuclides with data that are also in the chain
self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides
if nuc.name in self.nuclides_with_data]
# Extract number densities from the geometry / previous depletion run
self._extract_number(self.local_mats, volume, nuclides, self.prev_res)
# Create reaction rates array
self.reaction_rates = ReactionRates(
self.local_mats, self._burnable_nucs, self.chain.reactions)
# Get classes to assist working with tallies
if reaction_rate_mode == "direct":
self._rate_helper = DirectReactionRateHelper(
self.reaction_rates.n_nuc, self.reaction_rates.n_react)
elif reaction_rate_mode == "flux":
if reaction_rate_opts is None:
reaction_rate_opts = {}
# Ensure energy group boundaries were specified
if 'energies' not in reaction_rate_opts:
raise ValueError(
"Energy group boundaries must be specified in the "
"reaction_rate_opts argument when reaction_rate_mode is"
"set to 'flux'.")
self._rate_helper = FluxCollapseHelper(
self.reaction_rates.n_nuc,
self.reaction_rates.n_react,
**reaction_rate_opts
)
else:
raise ValueError("Invalid reaction rate mode.")
if normalization_mode == "fission-q":
self._normalization_helper = ChainFissionHelper()
elif normalization_mode == "energy-deposition":
score = "heating" if self.settings.photon_transport else "heating-local"
self._normalization_helper = EnergyScoreHelper(score)
else:
self._normalization_helper = SourceRateHelper()
# Select and create fission yield helper
fission_helper = self._fission_helpers[fission_yield_mode]
fission_yield_opts = (
{} if fission_yield_opts is None else fission_yield_opts)
self._yield_helper = fission_helper.from_operator(
self, **fission_yield_opts)
def __call__(self, vec, source_rate):
"""Runs a simulation.
Simulation will abort under the following circumstances:
1) No energy is computed using OpenMC tallies.
Parameters
----------
vec : list of numpy.ndarray
Total atoms to be used in function.
source_rate : float
Power in [W] or source rate in [neutron/sec]
Returns
-------
openmc.deplete.OperatorResult
Eigenvalue and reaction rates resulting from transport operator
"""
# Reset results in OpenMC
openmc.lib.reset()
# Update the number densities regardless of the source rate
self.number.set_density(vec)
self._update_materials()
# If the source rate is zero, return zero reaction rates without running
# a transport solve
if source_rate == 0.0:
rates = self.reaction_rates.copy()
rates.fill(0.0)
return OperatorResult(ufloat(0.0, 0.0), rates)
# Prevent OpenMC from complaining about re-creating tallies
openmc.reset_auto_ids()
# Update tally nuclides data in preparation for transport solve
nuclides = self._get_tally_nuclides()
self._rate_helper.nuclides = nuclides
self._normalization_helper.nuclides = nuclides
self._yield_helper.update_tally_nuclides(nuclides)
# Run OpenMC
openmc.lib.run()
openmc.lib.reset_timers()
# Extract results
op_result = self._unpack_tallies_and_normalize(source_rate)
return copy.deepcopy(op_result)
@staticmethod
def write_bos_data(step):
"""Write a state-point file with beginning of step data
Parameters
----------
step : int
Current depletion step including restarts
"""
openmc.lib.statepoint_write(
"openmc_simulation_n{}.h5".format(step),
write_source=False)
def _differentiate_burnable_mats(self):
"""Assign distribmats for each burnable material
"""
# Count the number of instances for each cell and material
self.geometry.determine_paths(instances_only=True)
# Extract all burnable materials which have multiple instances
distribmats = set(
[mat for mat in self.materials
if mat.depletable and mat.num_instances > 1])
for mat in distribmats:
if mat.volume is None:
raise RuntimeError("Volume not specified for depletable "
"material with ID={}.".format(mat.id))
mat.volume /= mat.num_instances
if distribmats:
# Assign distribmats to cells
for cell in self.geometry.get_all_material_cells().values():
if cell.fill in distribmats:
mat = cell.fill
cell.fill = [mat.clone()
for i in range(cell.num_instances)]
def _get_burnable_mats(self):
"""Determine depletable materials, volumes, and nuclides
Returns
-------
burnable_mats : list of str
List of burnable material IDs
volume : OrderedDict of str to float
Volume of each material in [cm^3]
nuclides : list of str
Nuclides in order of how they'll appear in the simulation.
"""
burnable_mats = set()
model_nuclides = set()
volume = OrderedDict()
self.heavy_metal = 0.0
# Iterate once through the geometry to get dictionaries
for mat in self.materials:
for nuclide in mat.get_nuclides():
model_nuclides.add(nuclide)
if mat.depletable:
burnable_mats.add(str(mat.id))
if mat.volume is None:
raise RuntimeError("Volume not specified for depletable "
"material with ID={}.".format(mat.id))
volume[str(mat.id)] = mat.volume
self.heavy_metal += mat.fissionable_mass
# Make sure there are burnable materials
if not burnable_mats:
raise RuntimeError(
"No depletable materials were found in the model.")
# Sort the sets
burnable_mats = sorted(burnable_mats, key=int)
model_nuclides = sorted(model_nuclides)
# Construct a global nuclide dictionary, burned first
nuclides = list(self.chain.nuclide_dict)
for nuc in model_nuclides:
if nuc not in nuclides:
nuclides.append(nuc)
return burnable_mats, volume, nuclides
def _extract_number(self, local_mats, volume, nuclides, prev_res=None):
"""Construct AtomNumber using geometry
Parameters
----------
local_mats : list of str
Material IDs to be managed by this process
volume : OrderedDict of str to float
Volumes for the above materials in [cm^3]
nuclides : list of str
Nuclides to be used in the simulation.
prev_res : Results, optional
Results from a previous depletion calculation
"""
self.number = AtomNumber(local_mats, nuclides, volume, len(self.chain))
if self.dilute_initial != 0.0:
for nuc in self._burnable_nucs:
self.number.set_atom_density(np.s_[:], nuc, self.dilute_initial)
# Now extract and store the number densities
# From the geometry if no previous depletion results
if prev_res is None:
for mat in self.materials:
if str(mat.id) in local_mats:
self._set_number_from_mat(mat)
# Else from previous depletion results
else:
for mat in self.materials:
if str(mat.id) in local_mats:
self._set_number_from_results(mat, prev_res)
def _set_number_from_mat(self, mat):
"""Extracts material and number densities from openmc.Material
Parameters
----------
mat : openmc.Material
The material to read from
"""
mat_id = str(mat.id)
for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items():
atom_per_cc = atom_per_bcm * 1.0e24
self.number.set_atom_density(mat_id, nuclide, atom_per_cc)
def _set_number_from_results(self, mat, prev_res):
"""Extracts material nuclides and number densities.
If the nuclide concentration's evolution is tracked, the densities come
from depletion results. Else, densities are extracted from the geometry
in the summary.
Parameters
----------
mat : openmc.Material
The material to read from
prev_res : Results
Results from a previous depletion calculation
"""
mat_id = str(mat.id)
# Get nuclide lists from geometry and depletion results
depl_nuc = prev_res[-1].nuc_to_ind
geom_nuc_densities = mat.get_nuclide_atom_densities()
# Merge lists of nuclides, with the same order for every calculation
geom_nuc_densities.update(depl_nuc)
for nuclide, atom_per_bcm in geom_nuc_densities.items():
if nuclide in depl_nuc:
concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1]
volume = prev_res[-1].volume[mat_id]
atom_per_cc = concentration / volume
else:
atom_per_cc = atom_per_bcm * 1.0e24
self.number.set_atom_density(mat_id, nuclide, atom_per_cc)
def initial_condition(self):
"""Performs final setup and returns initial condition.
Returns
-------
list of numpy.ndarray
Total density for initial conditions.
"""
# Create XML files
if comm.rank == 0:
self.geometry.export_to_xml()
self.settings.export_to_xml()
self._generate_materials_xml()
# Initialize OpenMC library
comm.barrier()
if not openmc.lib.is_initialized:
openmc.lib.init(intracomm=comm)
# Generate tallies in memory
materials = [openmc.lib.materials[int(i)]
for i in self.burnable_mats]
self._rate_helper.generate_tallies(materials, self.chain.reactions)
self._normalization_helper.prepare(
self.chain.nuclides, self.reaction_rates.index_nuc)
# Tell fission yield helper what materials this process is
# responsible for
self._yield_helper.generate_tallies(
materials, tuple(sorted(self._mat_index_map.values())))
# Return number density vector
return list(self.number.get_mat_slice(np.s_[:]))
def finalize(self):
"""Finalize a depletion simulation and release resources."""
if self.cleanup_when_done:
openmc.lib.finalize()
def _update_materials(self):
"""Updates material compositions in OpenMC on all processes."""
for rank in range(comm.size):
number_i = comm.bcast(self.number, root=rank)
for mat in number_i.materials:
nuclides = []
densities = []
for nuc in number_i.nuclides:
if nuc in self.nuclides_with_data:
val = 1.0e-24 * number_i.get_atom_density(mat, nuc)
# If nuclide is zero, do not add to the problem.
if val > 0.0:
if self.round_number:
val_magnitude = np.floor(np.log10(val))
val_scaled = val / 10**val_magnitude
val_round = round(val_scaled, 8)
val = val_round * 10**val_magnitude
nuclides.append(nuc)
densities.append(val)
else:
# Only output warnings if values are significantly
# negative. CRAM does not guarantee positive values.
if val < -1.0e-21:
print("WARNING: nuclide ", nuc, " in material ", mat,
" is negative (density = ", val, " at/barn-cm)")
number_i[mat, nuc] = 0.0
# Update densities on C API side
mat_internal = openmc.lib.materials[int(mat)]
mat_internal.set_densities(nuclides, densities)
#TODO Update densities on the Python side, otherwise the
# summary.h5 file contains densities at the first time step
def _generate_materials_xml(self):
"""Creates materials.xml from self.number.
Due to uncertainty with how MPI interacts with OpenMC API, this
constructs the XML manually. The long term goal is to do this
through direct memory writing.
"""
# Sort nuclides according to order in AtomNumber object
nuclides = list(self.number.nuclides)
for mat in self.materials:
mat._nuclides.sort(key=lambda x: nuclides.index(x[0]))
self.materials.export_to_xml()
def _get_tally_nuclides(self):
"""Determine nuclides that should be tallied for reaction rates.
This method returns a list of all nuclides that have neutron data and
are listed in the depletion chain. Technically, we should tally nuclides
that may not appear in the depletion chain because we still need to get
the fission reaction rate for these nuclides in order to normalize
power, but that is left as a future exercise.
Returns
-------
list of str
Tally nuclides
"""
nuc_set = set()
# Create the set of all nuclides in the decay chain in materials marked
# for burning in which the number density is greater than zero.
for nuc in self.number.nuclides:
if nuc in self.nuclides_with_data:
if np.sum(self.number[:, nuc]) > 0.0:
nuc_set.add(nuc)
# Communicate which nuclides have nonzeros to rank 0
if comm.rank == 0:
for i in range(1, comm.size):
nuc_newset = comm.recv(source=i, tag=i)
nuc_set |= nuc_newset
else:
comm.send(nuc_set, dest=0, tag=comm.rank)
if comm.rank == 0:
# Sort nuclides in the same order as self.number
nuc_list = [nuc for nuc in self.number.nuclides
if nuc in nuc_set]
else:
nuc_list = None
# Store list of tally nuclides on each process
nuc_list = comm.bcast(nuc_list)
return [nuc for nuc in nuc_list if nuc in self.chain]
def _unpack_tallies_and_normalize(self, source_rate):
"""Unpack tallies from OpenMC and return an operator result
This method uses OpenMC's C API bindings to determine the k-effective
value and reaction rates from the simulation. The reaction rates are
normalized by a helper class depending on the method being used.
Parameters
----------
source_rate : float
Power in [W] or source rate in [neutron/sec]
Returns
-------
openmc.deplete.OperatorResult
Eigenvalue and reaction rates resulting from transport operator
"""
rates = self.reaction_rates
rates.fill(0.0)
# Get k and uncertainty
keff = ufloat(*openmc.lib.keff())
# Extract tally bins
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]
# Keep track of energy produced from all reactions in eV per source
# particle
self._normalization_helper.reset()
self._yield_helper.unpack()
# Store fission yield dictionaries
fission_yields = []
# Create arrays to store fission Q values, reaction rates, and nuclide
# numbers, zeroed out in material iteration
number = np.empty(rates.n_nuc)
fission_ind = rates.index_rx.get("fission")
# Extract results
for i, mat in enumerate(self.local_mats):
# Get tally index
mat_index = self._mat_index_map[mat]
# Zero out reaction rates and nuclide numbers
number.fill(0.0)
# Get new number densities
for nuc, i_nuc_results in zip(nuclides, nuc_ind):
number[i_nuc_results] = self.number[mat, nuc]
tally_rates = self._rate_helper.get_material_rates(
mat_index, nuc_ind, react_ind)
# Compute fission yields for this material
fission_yields.append(self._yield_helper.weighted_yields(i))
# Accumulate energy from fission
if fission_ind is not None:
self._normalization_helper.update(tally_rates[:, fission_ind])
# Divide by total number and store
rates[i] = self._rate_helper.divide_by_adens(number)
# Scale reaction rates to obtain units of reactions/sec
rates *= self._normalization_helper.factor(source_rate)
# Store new fission yields on the chain
self.chain.fission_yields = fission_yields
return OperatorResult(keff, rates)
def _get_nuclides_with_data(self, cross_sections):
"""Loads cross_sections.xml file to find nuclides with neutron data"""
nuclides = set()
data_lib = DataLibrary.from_xml(cross_sections)
for library in data_lib.libraries:
if library['type'] != 'neutron':
continue
for name in library['materials']:
if name not in nuclides:
nuclides.add(name)
return nuclides
def get_results_info(self):
"""Returns volume list, material lists, and nuc lists.
Returns
-------
volume : dict of str float
Volumes corresponding to materials in full_burn_dict
nuc_list : list of str
A list of all nuclide names. Used for sorting the simulation.
burn_list : list of int
A list of all material IDs to be burned. Used for sorting the simulation.
full_burn_list : list
List of all burnable material IDs
"""
nuc_list = self.number.burnable_nuclides
burn_list = self.local_mats
volume = {}
for i, mat in enumerate(burn_list):
volume[mat] = self.number.volume[i]
# Combine volume dictionaries across processes
volume_list = comm.allgather(volume)
volume = {k: v for d in volume_list for k, v in d.items()}
return volume, nuc_list, burn_list, self.burnable_mats

50
openmc/deplete/results.py
View file

@ -16,6 +16,19 @@ __all__ = ["Results", "ResultsList"]
def _get_time_as(seconds, units):
"""Converts the time in seconds to time in different units
Parameters
----------
seconds : float
The time to convert expressed in seconds
units : {"s", "min", "h", "d", "a"}
The units to convert time into. Available options are seconds ``"s"``,
minutes ``"min"``, hours ``"h"`` days ``"d"``, Julian years ``"a"``
"""
if units == "a":
return seconds / (60 * 60 * 24 * 365.25) # 365.25 due to the leap year
if units == "d":
return seconds / (60 * 60 * 24)
elif units == "h":
@ -81,9 +94,9 @@ class Results(list):
.. note::
Initial values for some isotopes that do not appear in
initial concentrations may be non-zero, depending on the
value of the :attr:`openmc.deplete.Operator.dilute_initial`
attribute. The :class:`openmc.deplete.Operator` class adds isotopes
according to this setting, which can be set to zero.
value of the :attr:`openmc.deplete.CoupledOperator.dilute_initial`
attribute. The :class:`openmc.deplete.CoupledOperator` class adds
isotopes according to this setting, which can be set to zero.
Parameters
----------
@ -95,9 +108,10 @@ class Results(list):
Units for the returned concentration. Default is ``"atoms"``
.. versionadded:: 0.12
time_units : {"s", "min", "h", "d"}, optional
time_units : {"s", "min", "h", "d", "a"}, optional
Units for the returned time array. Default is ``"s"`` to
return the value in seconds.
return the value in seconds. Other options are minutes ``"min"``,
hours ``"h"``, days ``"d"``, and Julian years ``"a"``.
.. versionadded:: 0.12
@ -109,7 +123,7 @@ class Results(list):
Concentration of specified nuclide in units of ``nuc_units``
"""
cv.check_value("time_units", time_units, {"s", "d", "min", "h"})
cv.check_value("time_units", time_units, {"s", "d", "min", "h", "a"})
cv.check_value("nuc_units", nuc_units,
{"atoms", "atom/b-cm", "atom/cm3"})
@ -145,8 +159,8 @@ class Results(list):
Initial values for some isotopes that do not appear in
initial concentrations may be non-zero, depending on the
value of :class:`openmc.deplete.Operator` ``dilute_initial``
The :class:`openmc.deplete.Operator` adds isotopes according
value of :class:`openmc.deplete.CoupledOperator` ``dilute_initial``
The :class:`openmc.deplete.CoupledOperator` adds isotopes according
to this setting, which can be set to zero.
Parameters
@ -190,8 +204,10 @@ class Results(list):
Parameters
----------
time_units : {"s", "d", "h", "min"}, optional
Desired units for the times array
time_units : {"s", "d", "min", "h", "a"}, optional
Desired units for the times array. Options are seconds ``"s"``,
minutes ``"min"``, hours ``"h"``, days ``"d"``, and Julian years
``"a"``.
Returns
-------
@ -203,7 +219,7 @@ class Results(list):
1 contains the associated uncertainty
"""
cv.check_value("time_units", time_units, {"s", "d", "min", "h"})
cv.check_value("time_units", time_units, {"s", "d", "min", "h", "a"})
times = np.empty_like(self, dtype=float)
eigenvalues = np.empty((len(self), 2), dtype=float)
@ -257,9 +273,10 @@ class Results(list):
Parameters
----------
time_units : {"s", "d", "h", "min"}, optional
time_units : {"s", "d", "min", "h", "a"}, optional
Return the vector in these units. Default is to
convert to days
convert to days ``"d"``. Other options are seconds ``"s"``, minutes
``"min"``, hours ``"h"``, days ``"d"``, and Julian years ``"a"``.
Returns
-------
@ -267,7 +284,7 @@ class Results(list):
1-D vector of time points
"""
cv.check_value("time_units", time_units, {"s", "d", "min", "h"})
cv.check_value("time_units", time_units, {"s", "d", "min", "h", "a"})
times = np.fromiter(
(r.time[0] for r in self),
@ -296,8 +313,9 @@ class Results(list):
----------
time : float
Desired point in time
time_units : {"s", "d", "min", "h"}, optional
Units on ``time``. Default: days
time_units : {"s", "d", "min", "h", "a"}, optional
Units on ``time``. Default: days ``"d"``. Other options are seconds
``"s"``, minutes ``"min"``, hours ``"h"`` and Julian years ``"a"``.
atol : float, optional
Absolute tolerance (in ``time_units``) if ``time`` is not
found.

10
openmc/deplete/stepresult.py
View file

@ -462,7 +462,7 @@ class StepResult:
Parameters
----------
op : openmc.deplete.TransportOperator
op : openmc.deplete.abc.TransportOperator
The operator used to generate these results.
x : list of list of numpy.array
The prior x vectors. Indexed [i][cell] using the above equation.
@ -495,7 +495,13 @@ class StepResult:
for mat_i in range(n_mat):
results[i, mat_i, :] = x[i][mat_i]
results.k = [(r.k.nominal_value, r.k.std_dev) for r in op_results]
ks = []
for r in op_results:
if isinstance(r.k, type(None)):
ks += [(None, None)]
else:
ks += [(r.k.nominal_value, r.k.std_dev)]
results.k = ks
results.rates = [r.rates for r in op_results]
results.time = t
results.source_rate = source_rate

4
openmc/material.py
View file

@ -861,7 +861,7 @@ class Material(IDManagerMixin):
elif self.density_units == 'atom/b-cm':
density = self.density
elif self.density_units == 'atom/cm3' or self.density_units == 'atom/cc':
density = 1.E-24 * self.density
density = 1.e-24 * self.density
# For ease of processing split out nuc, nuc_density,
# and nuc_density_type into separate arrays
@ -897,7 +897,7 @@ class Material(IDManagerMixin):
# Convert the mass density to an atom density
if not density_in_atom:
density = -density / self.average_molar_mass * 1.E-24 \
density = -density / self.average_molar_mass * 1.e-24 \
* openmc.data.AVOGADRO
nuc_densities = density * nuc_densities

4
openmc/model/surface_composite.py
View file

@ -451,10 +451,10 @@ class RectangularParallelepiped(CompositeSurface):
self.zmax = openmc.ZPlane(z0=zmax, **kwargs)
def __neg__(self):
return +self.xmin & -self.xmax & +self.ymin & -self.ymax & +self.zmin & -self.zmax
return -self.xmax & +self.xmin & -self.ymax & +self.ymin & -self.zmax & +self.zmin
def __pos__(self):
return -self.xmin | +self.xmax | -self.ymin | +self.ymax | -self.zmin | +self.zmax
return +self.xmax | -self.xmin | +self.ymax | -self.ymin | +self.zmax | -self.zmin
class XConeOneSided(CompositeSurface):

135
openmc/stats/univariate.py
View file

@ -1,6 +1,7 @@
from abc import ABC, abstractmethod
from collections import defaultdict
from collections.abc import Iterable
from copy import deepcopy
from numbers import Real
from xml.etree import ElementTree as ET
@ -95,9 +96,9 @@ class Discrete(Univariate):
Attributes
----------
x : Iterable of float
x : numpy.ndarray
Values of the random variable
p : Iterable of float
p : numpy.ndarray
Discrete probability for each value
"""
@ -122,7 +123,7 @@ class Discrete(Univariate):
if isinstance(x, Real):
x = [x]
cv.check_type('discrete values', x, Iterable, Real)
self._x = x
self._x = np.array(x, dtype=float)
@p.setter
def p(self, p):
@ -131,14 +132,15 @@ class Discrete(Univariate):
cv.check_type('discrete probabilities', p, Iterable, Real)
for pk in p:
cv.check_greater_than('discrete probability', pk, 0.0, True)
self._p = p
self._p = np.array(p, dtype=float)
def cdf(self):
return np.insert(np.cumsum(self.p), 0, 0.0)
def sample(self, n_samples=1, seed=None):
np.random.seed(seed)
return np.random.choice(self.x, n_samples, p=self.p)
p = self.p / self.p.sum()
return np.random.choice(self.x, n_samples, p=p)
def normalize(self):
"""Normalize the probabilities stored on the distribution"""
@ -220,6 +222,18 @@ class Discrete(Univariate):
p_arr = np.array([p_merged[x] for x in x_arr])
return cls(x_arr, p_arr)
def integral(self):
"""Return integral of distribution
.. versionadded:: 0.13.1
Returns
-------
float
Integral of discrete distribution
"""
return np.sum(self.p)
class Uniform(Univariate):
"""Distribution with constant probability over a finite interval [a,b]
@ -825,9 +839,9 @@ class Tabular(Univariate):
Attributes
----------
x : Iterable of float
x : numpy.ndarray
Tabulated values of the random variable
p : Iterable of float
p : numpy.ndarray
Tabulated probabilities
interpolation : {'histogram', 'linear-linear', 'linear-log', 'log-linear', 'log-log'}, optional
Indicate whether the density function is constant between tabulated
@ -860,7 +874,7 @@ class Tabular(Univariate):
@x.setter
def x(self, x):
cv.check_type('tabulated values', x, Iterable, Real)
self._x = x
self._x = np.array(x, dtype=float)
@p.setter
def p(self, p):
@ -868,7 +882,7 @@ class Tabular(Univariate):
if not self._ignore_negative:
for pk in p:
cv.check_greater_than('tabulated probability', pk, 0.0, True)
self._p = p
self._p = np.array(p, dtype=float)
@interpolation.setter
def interpolation(self, interpolation):
@ -881,8 +895,8 @@ class Tabular(Univariate):
'distributions using histogram or '
'linear-linear interpolation')
c = np.zeros_like(self.x)
x = np.asarray(self.x)
p = np.asarray(self.p)
x = self.x
p = self.p
if self.interpolation == 'histogram':
c[1:] = p[:-1] * np.diff(x)
@ -922,7 +936,7 @@ class Tabular(Univariate):
def normalize(self):
"""Normalize the probabilities stored on the distribution"""
self.p = np.asarray(self.p) / self.cdf().max()
self.p /= self.cdf().max()
def sample(self, n_samples=1, seed=None):
if not self.interpolation in ('histogram', 'linear-linear'):
@ -931,10 +945,11 @@ class Tabular(Univariate):
'linear-linear interpolation')
np.random.seed(seed)
xi = np.random.rand(n_samples)
cdf = self.cdf()
cdf /= cdf.max()
# always use normalized probabilities when sampling
cdf = self.cdf()
p = self.p / cdf.max()
cdf /= cdf.max()
# get CDF bins that are above the
# sampled values
@ -949,7 +964,7 @@ class Tabular(Univariate):
# the random number is less than the next cdf
# entry
x_i = self.x[cdf_idx]
p_i = self.p[cdf_idx]
p_i = p[cdf_idx]
if self.interpolation == 'histogram':
# mask where probability is greater than zero
@ -967,7 +982,7 @@ class Tabular(Univariate):
# get variable and probability values for the
# next entry
x_i1 = self.x[cdf_idx + 1]
p_i1 = self.p[cdf_idx + 1]
p_i1 = p[cdf_idx + 1]
# compute slope between entries
m = (p_i1 - p_i) / (x_i1 - x_i)
# set values for zero slope
@ -1027,6 +1042,24 @@ class Tabular(Univariate):
p = params[len(params)//2:]
return cls(x, p, interpolation)
def integral(self):
"""Return integral of distribution
.. versionadded: 0.13.1
Returns
-------
float
Integral of tabular distrbution
"""
if self.interpolation == 'histogram':
return np.sum(np.diff(self.x) * self.p[:-1])
elif self.interpolation == 'linear-linear':
return np.trapz(self.p, self.x)
else:
raise NotImplementedError(
f'integral() not supported for {self.inteprolation} interpolation')
class Legendre(Univariate):
r"""Probability density given by a Legendre polynomial expansion
@ -1135,9 +1168,10 @@ class Mixture(Univariate):
def sample(self, n_samples=1, seed=None):
np.random.seed(seed)
idx = np.random.choice(self.distribution, n_samples, p=self.probability)
idx = np.random.choice(range(len(self.distribution)),
n_samples, p=self.probability)
out = np.zeros_like(idx)
out = np.empty_like(idx, dtype=float)
for i in np.unique(idx):
n_dist_samples = np.count_nonzero(idx == i)
samples = self.distribution[i].sample(n_dist_samples)
@ -1199,3 +1233,68 @@ class Mixture(Univariate):
distribution.append(Univariate.from_xml_element(pair.find("dist")))
return cls(probability, distribution)
def integral(self):
"""Return integral of the distribution
.. versionadded:: 0.13.1
Returns
-------
float
Integral of the distribution
"""
return sum([
p*dist.integral()
for p, dist in zip(self.probability, self.distribution)
])
def combine_distributions(dists, probs):
"""Combine distributions with specified probabilities
This function can be used to combine multiple instances of
:class:`~openmc.stats.Discrete` and `~openmc.stats.Tabular`. Multiple
discrete distributions are merged into a single distribution and the
remainder of the distributions are put into a :class:`~openmc.stats.Mixture`
distribution.
.. versionadded:: 0.13.1
Parameters
----------
dists : iterable of openmc.stats.Univariate
Distributions to combine
probs : iterable of float
Probability (or intensity) of each distribution
"""
# Get copy of distribution list so as not to modify the argument
dist_list = deepcopy(dists)
# Get list of discrete/continuous distribution indices
discrete_index = [i for i, d in enumerate(dist_list) if isinstance(d, Discrete)]
cont_index = [i for i, d in enumerate(dist_list) if isinstance(d, Tabular)]
# Apply probabilites to continuous distributions
for i in cont_index:
dist = dist_list[i]
dist.p *= probs[i]
if discrete_index:
# Create combined discrete distribution
dist_discrete = [dist_list[i] for i in discrete_index]
discrete_probs = [probs[i] for i in discrete_index]
combined_dist = Discrete.merge(dist_discrete, discrete_probs)
# Replace multiple discrete distributions with merged
for idx in reversed(discrete_index):
dist_list.pop(idx)
dist_list.append(combined_dist)
# Combine discrete and continuous if present
if len(dist_list) > 1:
probs = [d.integral() for d in dist_list]
dist_list[:] = [Mixture(probs, dist_list.copy())]
return dist_list[0]

9
openmc/universe.py
View file

@ -334,6 +334,13 @@ class Universe(UniverseBase):
if seed is not None:
model.settings.seed = seed
# Determine whether any materials contains macroscopic data and if
# so, set energy mode accordingly
for mat in self.get_all_materials().values():
if mat._macroscopic is not None:
model.settings.energy_mode = 'multi-group'
break
# Create plot object matching passed arguments
plot = openmc.Plot()
plot.origin = origin
@ -785,7 +792,7 @@ class DAGMCUniverse(UniverseBase):
Universe instance
"""
bounding_cell = openmc.Cell(fill=self, cell_id =bounding_cell_id, region=self.bounding_region(**kwargs))
bounding_cell = openmc.Cell(fill=self, cell_id=bounding_cell_id, region=self.bounding_region(**kwargs))
return openmc.Universe(cells=[bounding_cell])
@classmethod

13
tests/micro_xs_simple.csv Normal file
View file

@ -0,0 +1,13 @@
nuclide,"(n,gamma)",fission
U234,22.231989822002465,0.49620744663749855
U235,10.479008971197121,48.41787337164604
U238,0.8673334105437324,0.10467880588762352
U236,8.65171044607122,0.31948392400019293
O16,7.497851000107524e-05,0.0
O17,0.0004079227797153371,0.0
I135,6.842395323713927,0.0
Xe135,227463.86426990604,0.0
Xe136,0.02317896034753588,0.0
Cs135,2.1721665580713623,0.0
Gd157,12786.09939237018,0.0
Gd156,3.4006085445846983,0.0
1 nuclide (n,gamma) fission
2 U234 22.231989822002465 0.49620744663749855
3 U235 10.479008971197121 48.41787337164604
4 U238 0.8673334105437324 0.10467880588762352
5 U236 8.65171044607122 0.31948392400019293
6 O16 7.497851000107524e-05 0.0
7 O17 0.0004079227797153371 0.0
8 I135 6.842395323713927 0.0
9 Xe135 227463.86426990604 0.0
10 Xe136 0.02317896034753588 0.0
11 Cs135 2.1721665580713623 0.0
12 Gd157 12786.09939237018 0.0
13 Gd156 3.4006085445846983 0.0

4
tests/regression_tests/cpp_driver/test.py
View file

@ -48,8 +48,8 @@ def cpp_driver(request):
finally:
# Remove local build directory when test is complete
shutil.rmtree('build')
os.remove('CMakeLists.txt')
shutil.rmtree(request.node.path.parent / 'build')
os.remove(request.node.path.parent / 'CMakeLists.txt')
@pytest.fixture

0
tests/regression_tests/deplete_no_transport/__init__.py Normal file
View file

106
tests/regression_tests/deplete_no_transport/test.py Normal file
View file

@ -0,0 +1,106 @@
""" Transport-free depletion test suite """
from pathlib import Path
import numpy as np
import pytest
import openmc
import openmc.deplete
from openmc.deplete import IndependentOperator, MicroXS
@pytest.fixture(scope="module")
def fuel():
fuel = openmc.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.volume = np.pi * 0.42 ** 2
return fuel
@pytest.mark.parametrize("multiproc, from_nuclides, normalization_mode, power, flux", [
(True, True,'source-rate', None, 1164719970082145.0),
(False, True, 'source-rate', None, 1164719970082145.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, 'fission-q', 174, None),
(False, False, 'fission-q', 174, None)])
def test_no_transport_from_nuclides(run_in_tmpdir, fuel, multiproc, from_nuclides, normalization_mode, power, flux):
"""Transport free system test suite.
Runs an OpenMC transport-free depletion calculation and verifies
that the outputs match a reference file.
"""
# Create operator
micro_xs_file = Path(__file__).parents[2] / 'micro_xs_simple.csv'
micro_xs = MicroXS.from_csv(micro_xs_file)
chain_file = Path(__file__).parents[2] / 'chain_simple.xml'
if from_nuclides:
nuclides = {}
for nuc, dens in fuel.get_nuclide_atom_densities().items():
nuclides[nuc] = dens
op = IndependentOperator.from_nuclides(
fuel.volume, nuclides, micro_xs, chain_file, normalization_mode=normalization_mode)
else:
op = IndependentOperator(openmc.Materials([fuel]), micro_xs, chain_file, normalization_mode=normalization_mode)
# Power and timesteps
dt = [30] # single step
# Perform simulation using the predictor algorithm
openmc.deplete.pool.USE_MULTIPROCESSING = multiproc
openmc.deplete.PredictorIntegrator(
op, dt, power=power, source_rates=flux, timestep_units='d').integrate()
# Get path to test and reference results
path_test = op.output_dir / 'depletion_results.h5'
if flux is not None:
ref_path = 'test_reference_source_rate.h5'
else:
ref_path = 'test_reference_fission_q.h5'
path_reference = Path(__file__).with_name(ref_path)
# Load the reference/test results
res_test = openmc.deplete.Results(path_test)
res_ref = openmc.deplete.Results(path_reference)
# Assert same mats
for mat in res_ref[0].mat_to_ind:
assert mat in res_test[0].mat_to_ind, \
"Material {} not in new results.".format(mat)
for nuc in res_ref[0].nuc_to_ind:
assert nuc in res_test[0].nuc_to_ind, \
"Nuclide {} not in new results.".format(nuc)
for mat in res_test[0].mat_to_ind:
assert mat in res_ref[0].mat_to_ind, \
"Material {} not in old results.".format(mat)
for nuc in res_test[0].nuc_to_ind:
assert nuc in res_ref[0].nuc_to_ind, \
"Nuclide {} not in old results.".format(nuc)
tol = 1.0e-6
for mat in res_test[0].mat_to_ind:
for nuc in res_test[0].nuc_to_ind:
_, y_test = res_test.get_atoms(mat, nuc)
_, y_old = res_ref.get_atoms(mat, nuc)
# Test each point
correct = True
for i, ref in enumerate(y_old):
if ref != y_test[i]:
if ref != 0.0:
correct = np.abs(y_test[i] - ref) / ref <= tol
else:
correct = False
assert correct, "Discrepancy in mat {} and nuc {}\n{}\n{}".format(
mat, nuc, y_old, y_test)

BIN
tests/regression_tests/deplete_no_transport/test_reference_fission_q.h5 Normal file
View file

Binary file not shown.

BIN
tests/regression_tests/deplete_no_transport/test_reference_source_rate.h5 Normal file
View file

Binary file not shown.

0
tests/regression_tests/deplete_with_transport/__init__.py Normal file
View file

0
tests/regression_tests/deplete/example_geometry.py → tests/regression_tests/deplete_with_transport/example_geometry.py
View file

0
tests/regression_tests/deplete/last_step_reference_materials.xml → tests/regression_tests/deplete_with_transport/last_step_reference_materials.xml
View file

2
tests/regression_tests/deplete/test.py → tests/regression_tests/deplete_with_transport/test.py
View file

@ -12,7 +12,7 @@ from openmc.data import JOULE_PER_EV
import openmc.deplete
from tests.regression_tests import config
from example_geometry import generate_problem
from .example_geometry import generate_problem
@pytest.fixture(scope="module")

0
tests/regression_tests/deplete/test_reference.h5 → tests/regression_tests/deplete_with_transport/test_reference.h5
View file

0
tests/regression_tests/microxs/__init.py__ Normal file
View file

52
tests/regression_tests/microxs/test.py Normal file
View file

@ -0,0 +1,52 @@
"""Test one-group cross section generation"""
from pathlib import Path
import numpy as np
import pytest
import openmc
from openmc.deplete import MicroXS
CHAIN_FILE = Path(__file__).parents[2] / "chain_simple.xml"
@pytest.fixture(scope="module")
def model():
fuel = openmc.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)
clad = openmc.Material(name="clad")
clad.add_element("Zr", 1)
clad.set_density("g/cc", 6)
water = openmc.Material(name="water")
water.add_element("O", 1)
water.add_element("H", 2)
water.set_density("g/cc", 1.0)
water.add_s_alpha_beta("c_H_in_H2O")
radii = [0.42, 0.45]
fuel.volume = np.pi * radii[0] ** 2
materials = openmc.Materials([fuel, clad, water])
pin_surfaces = [openmc.ZCylinder(r=r) for r in radii]
pin_univ = openmc.model.pin(pin_surfaces, materials)
bound_box = openmc.rectangular_prism(1.24, 1.24, boundary_type="reflective")
root_cell = openmc.Cell(fill=pin_univ, region=bound_box)
geometry = openmc.Geometry([root_cell])
settings = openmc.Settings()
settings.particles = 1000
settings.inactive = 10
settings.batches = 50
return openmc.Model(geometry, materials, settings)
def test_from_model(model):
ref_xs = MicroXS.from_csv('test_reference.csv')
test_xs = MicroXS.from_model(model, model.materials[0], CHAIN_FILE)
np.testing.assert_allclose(test_xs, ref_xs, rtol=1e-11)

13
tests/regression_tests/microxs/test_reference.csv Normal file
View file

@ -0,0 +1,13 @@
nuclide,"(n,gamma)",fission
U234,22.231989822002465,0.49620744663749855
U235,10.479008971197121,48.41787337164604
U238,0.8673334105437324,0.10467880588762352
U236,8.65171044607122,0.31948392400019293
O16,7.497851000107524e-05,0.0
O17,0.0004079227797153371,0.0
I135,6.842395323713927,0.0
Xe135,227463.86426990604,0.0
Xe136,0.02317896034753588,0.0
Cs135,2.1721665580713623,0.0
Gd157,12786.09939237018,0.0
Gd156,3.4006085445846983,0.0
1 nuclide (n,gamma) fission
2 U234 22.231989822002465 0.49620744663749855
3 U235 10.479008971197121 48.41787337164604
4 U238 0.8673334105437324 0.10467880588762352
5 U236 8.65171044607122 0.31948392400019293
6 O16 7.497851000107524e-05 0.0
7 O17 0.0004079227797153371 0.0
8 I135 6.842395323713927 0.0
9 Xe135 227463.86426990604 0.0
10 Xe136 0.02317896034753588 0.0
11 Cs135 2.1721665580713623 0.0
12 Gd157 12786.09939237018 0.0
13 Gd156 3.4006085445846983 0.0

2
tests/regression_tests/torus/inputs_true.dat
View file

@ -3,7 +3,7 @@
<cell id="1" material="1" region="-2" universe="1" />
<cell id="2" material="1" region="-3" universe="1" />
<cell id="3" material="1" region="-1" universe="1" />
<cell id="4" material="2" region="4 -5 6 -7 8 -9 2 3 1" universe="1" />
<cell id="4" material="2" region="-5 4 -7 6 -9 8 2 3 1" universe="1" />
<surface coeffs="0.0 0.0 0.0 3 1.5 1" id="1" type="z-torus" />
<surface coeffs="6 0.0 0.0 3 1.5 1" id="2" type="x-torus" />
<surface coeffs="6 0.0 0.0 6 1 0.75" id="3" type="y-torus" />

35
tests/unit_tests/test_data_decay.py
View file

@ -23,6 +23,14 @@ def nb90():
return openmc.data.Decay.from_endf(filename)
@pytest.fixture(scope='module')
def ba137m():
"""Ba137_m1 decay data."""
endf_data = os.environ['OPENMC_ENDF_DATA']
filename = os.path.join(endf_data, 'decay', 'dec-056_Ba_137m1.endf')
return openmc.data.Decay.from_endf(filename)
@pytest.fixture(scope='module')
def u235_yields():
"""U235 fission product yield data."""
@ -48,6 +56,7 @@ def test_nb90_halflife(nb90):
ufloat_close(nb90.decay_constant, log(2.)/nb90.half_life)
ufloat_close(nb90.decay_energy, ufloat(2265527.5, 25159.400474401213))
def test_nb90_nuclide(nb90):
assert nb90.nuclide['atomic_number'] == 41
assert nb90.nuclide['mass_number'] == 90
@ -91,3 +100,29 @@ def test_fpy(u235_yields):
assert len(u235_yields.independent) == 3
thermal = u235_yields.independent[0]
ufloat_close(thermal['I135'], ufloat(0.0292737, 0.000819663))
def test_sources(ba137m, nb90):
# Running .sources twice should give same objects
sources = ba137m.sources
sources2 = ba137m.sources
for key in sources:
assert sources[key] is sources2[key]
# Each source should be a univariate distribution
for dist in sources.values():
assert isinstance(dist, openmc.stats.Univariate)
# Check for presence of 662 keV gamma ray in decay of Ba137m
gamma_source = ba137m.sources['photon']
assert isinstance(gamma_source, openmc.stats.Discrete)
b = np.isclose(gamma_source.x, 661657.)
assert np.count_nonzero(b) == 1
# Check value of decay/s/atom
idx = np.flatnonzero(b)[0]
assert gamma_source.p[idx] == pytest.approx(0.004069614)
# Nb90 decays by β+ and should emit positrons, electrons, and photons
sources = nb90.sources
assert len(set(sources.keys()) ^ {'positron', 'electron', 'photon'}) == 0

176
tests/unit_tests/test_data_kalbach_mann.py Normal file
View file

@ -0,0 +1,176 @@
"""Test of the Kalbach-Mann slope calculation when data are
retrieved from ENDF files."""
import os
from pathlib import Path
import pytest
import numpy as np
from openmc.data import IncidentNeutron
from openmc.data.kalbach_mann import _separation_energy, _AtomicRepresentation
from openmc.data import kalbach_slope
from openmc.data import KalbachMann
from . import needs_njoy
@pytest.fixture(scope='module')
def neutron():
"""Neutron AtomicRepresentation."""
return _AtomicRepresentation(z=0, a=1)
@pytest.fixture(scope='module')
def triton():
"""Triton AtomicRepresentation."""
return _AtomicRepresentation(z=1, a=3)
@pytest.fixture(scope='module')
def b10():
"""B10 AtomicRepresentation."""
return _AtomicRepresentation(z=5, a=10)
@pytest.fixture(scope='module')
def c12():
"""C12 AtomicRepresentation."""
return _AtomicRepresentation(z=6, a=12)
@pytest.fixture(scope='module')
def c13():
"""C13 AtomicRepresentation."""
return _AtomicRepresentation(z=6, a=13)
@pytest.fixture(scope='module')
def na23():
"""Na23 AtomicRepresentation."""
return _AtomicRepresentation(z=11, a=23)
def test_atomic_representation(neutron, triton, b10, c12, c13, na23):
"""Test the _AtomicRepresentation class."""
# Test instantiation from_za
assert b10 == _AtomicRepresentation.from_za(5010)
# Test addition
assert c13 + b10 == na23
# Test substraction
assert c13 - c12 == neutron
assert c13 - b10 == triton
# Test properties when no information for Kalbach-Mann are given
assert c13.a == 13
assert c13.z == 6
assert c13.n == 7
assert c13.za == 6013
# Test properties when information for Kalbach-Mann are given
assert triton.a == 3
assert triton.z == 1
assert triton.n == 2
assert triton.za == 1003
# Test instantiation errors
with pytest.raises(ValueError):
_AtomicRepresentation(z=5, a=1)
with pytest.raises(ValueError):
_AtomicRepresentation(z=-1, a=1)
with pytest.raises(ValueError):
_AtomicRepresentation(z=5, a=0)
with pytest.raises(ValueError):
_AtomicRepresentation(z=5, a=-2)
with pytest.raises(ValueError):
neutron - triton
def test_separation_energy(triton, b10, c13):
"""Comparison to hand-calculations on a simple example."""
assert _separation_energy(
compound=c13,
nucleus=b10,
particle=triton
) == pytest.approx(18.6880713)
def test_kalbach_slope():
"""Comparison to hand-calculations for n + c12 -> c13 -> triton + b10."""
energy_projectile = 10.2 # [eV]
energy_emitted = 5.4 # [eV]
# Check that NotImplementedError is raised if the projectile is not
# a neutron
with pytest.raises(NotImplementedError):
kalbach_slope(
energy_projectile=energy_projectile,
energy_emitted=energy_emitted,
za_projectile=1000,
za_emitted=1,
za_target=6012
)
assert kalbach_slope(
energy_projectile=energy_projectile,
energy_emitted=energy_emitted,
za_projectile=1,
za_emitted=1003,
za_target=6012
) == pytest.approx(0.8409921475)
@pytest.mark.parametrize(
"hdf5_filename, endf_filename", [
('O16.h5', 'n-008_O_016.endf'),
('Ca46.h5', 'n-020_Ca_046.endf'),
('Hg204.h5', 'n-080_Hg_204.endf')
]
)
def test_comparison_slope_hdf5(hdf5_filename, endf_filename):
"""Test the calculation of the Kalbach-Mann slope done by OpenMC
by comparing it to HDF5 data. The test is based on the first product
of MT=5 (neutron). The isotopes tested have been selected because the
corresponding products in ENDF/B-VII.1 are described using MF=6, LAW=1,
LANG=2 (i.e., Kalbach-Mann systematics) and the slope is not given
explicitly.
If an error occurs during the "validity check", this means that
the nuclear data evaluation has evolved and the distribution might
no longer be described using Kalbach-Mann systematics. Another
isotope needs to be identified and tested.
Warning: This test is valid as long as ENDF files are not directly
used to generate the HDF5 files used in the tests.
"""
# HDF5 data
hdf5_directory = Path(os.environ['OPENMC_CROSS_SECTIONS']).parent
hdf5_data = IncidentNeutron.from_hdf5(hdf5_directory / hdf5_filename)
hdf5_product = hdf5_data[5].products[0]
hdf5_distribution = hdf5_product.distribution[0]
# ENDF data
endf_directory = Path(os.environ['OPENMC_ENDF_DATA'])
endf_path = endf_directory / 'neutrons' / endf_filename
endf_data = IncidentNeutron.from_endf(endf_path)
endf_product = endf_data[5].products[0]
endf_distribution = endf_product.distribution[0]
# Validity check
assert isinstance(endf_distribution, KalbachMann)
assert isinstance(hdf5_distribution, KalbachMann)
assert endf_product.particle == hdf5_product.particle
assert len(endf_distribution.slope) == len(hdf5_distribution.slope)
# Results check
for i, hdf5_slope in enumerate(hdf5_distribution.slope):
assert endf_distribution._calculated_slope[i]
np.testing.assert_array_almost_equal(
endf_distribution.slope[i].y,
hdf5_slope.y,
decimal=5
)

55
tests/unit_tests/test_deplete_coupled_operator.py Normal file
View file

@ -0,0 +1,55 @@
"""Basic unit tests for openmc.deplete.CoupledOperator instantiation
"""
from pathlib import Path
import pytest
from openmc.deplete import CoupledOperator
import openmc
import numpy as np
CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml"
@pytest.fixture(scope="module")
def model():
fuel = openmc.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)
clad = openmc.Material(name="clad")
clad.add_element("Zr", 1)
clad.set_density("g/cc", 6)
water = openmc.Material(name="water")
water.add_element("O", 1)
water.add_element("H", 2)
water.set_density("g/cc", 1.0)
water.add_s_alpha_beta("c_H_in_H2O")
radii = [0.42, 0.45]
fuel.volume = np.pi * radii[0] ** 2
clad.volume = np.pi * (radii[1]**2 - radii[0]**2)
water.volume = 1.24**2 - (np.pi * radii[1]**2)
materials = openmc.Materials([fuel, clad, water])
pin_surfaces = [openmc.ZCylinder(r=r) for r in radii]
pin_univ = openmc.model.pin(pin_surfaces, materials)
bound_box = openmc.rectangular_prism(1.24, 1.24, boundary_type="reflective")
root_cell = openmc.Cell(fill=pin_univ, region=bound_box)
geometry = openmc.Geometry([root_cell])
settings = openmc.Settings()
settings.particles = 1000
settings.inactive = 10
settings.batches = 50
return openmc.Model(geometry, materials, settings)
def test_operator_init(model):
"""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."""
CoupledOperator(model, CHAIN_PATH)

36
tests/unit_tests/test_deplete_independent_operator.py Normal file
View file

@ -0,0 +1,36 @@
"""Basic unit tests for openmc.deplete.IndependentOperator instantiation
"""
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."""
volume = 1
nuclides = {'U234': 8.922411359424315e+18,
'U235': 9.98240191860822e+20,
'U238': 2.2192386373095893e+22,
'U236': 4.5724195495061115e+18,
'O16': 4.639065406771322e+22,
'O17': 1.7588724018066158e+19}
micro_xs = MicroXS.from_csv(ONE_GROUP_XS)
IndependentOperator.from_nuclides(
volume, nuclides, 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.volume = 1
materials = Materials([fuel])
IndependentOperator(materials, micro_xs, CHAIN_PATH)

59
tests/unit_tests/test_deplete_microxs.py Normal file
View file

@ -0,0 +1,59 @@
"""Basic unit tests for openmc.deplete.IndependentOperator instantiation
Modifies and resets environment variable OPENMC_CROSS_SECTIONS
to a custom file with new depletion_chain node
"""
from os import remove
from pathlib import Path
import pytest
from openmc.deplete import MicroXS
import numpy as np
ONE_GROUP_XS = Path(__file__).parents[1] / "micro_xs_simple.csv"
def test_from_array():
nuclides = [
'U234',
'U235',
'U238',
'U236',
'O16',
'O17',
'I135',
'Xe135',
'Xe136',
'Cs135',
'Gd157',
'Gd156']
reactions = ['fission', '(n,gamma)']
# These values are placeholders and are not at all
# physically meaningful.
data = np.array([[0.1, 0.],
[0.1, 0.],
[0.9, 0.],
[0.4, 0.],
[0., 0.],
[0., 0.],
[0., 0.1],
[0., 0.9],
[0., 0.],
[0., 0.],
[0., 0.1],
[0., 0.1]])
MicroXS.from_array(nuclides, reactions, data)
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])
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)
remove('temp_xs.csv')

8
tests/unit_tests/test_deplete_resultslist.py
View file

@ -11,7 +11,7 @@ import openmc.deplete
@pytest.fixture
def res():
"""Load the reference results"""
filename = (Path(__file__).parents[1] / 'regression_tests' / 'deplete'
filename = (Path(__file__).parents[1] / 'regression_tests' / 'deplete_with_transport'
/ 'test_reference.h5')
return openmc.deplete.Results(filename)
@ -70,14 +70,16 @@ def test_get_keff(res):
np.testing.assert_allclose(k[:, 1], u_ref)
@pytest.mark.parametrize("unit", ("s", "d", "min", "h"))
@pytest.mark.parametrize("unit", ("s", "d", "min", "h", "a"))
def test_get_steps(unit):
# Make a Results full of near-empty Result instances
# Just fill out a time schedule
results = openmc.deplete.Results()
# Time in units of unit
times = np.linspace(0, 100, num=5)
if unit == "d":
if unit == "a":
conversion_to_seconds = 60 * 60 * 24 * 365.25
elif unit == "d":
conversion_to_seconds = 60 * 60 * 24
elif unit == "h":
conversion_to_seconds = 60 * 60

108
tests/unit_tests/test_stats.py
View file

@ -6,6 +6,12 @@ import openmc
import openmc.stats
def assert_sample_mean(samples, expected_mean):
std_dev = samples.std() / np.sqrt(samples.size)
assert np.abs(expected_mean - samples.mean()) < 3*std_dev
def test_discrete():
x = [0.0, 1.0, 10.0]
p = [0.3, 0.2, 0.5]
@ -13,8 +19,8 @@ def test_discrete():
elem = d.to_xml_element('distribution')
d = openmc.stats.Discrete.from_xml_element(elem)
assert d.x == x
assert d.p == p
np.testing.assert_array_equal(d.x, x)
np.testing.assert_array_equal(d.p, p)
assert len(d) == len(x)
d = openmc.stats.Univariate.from_xml_element(elem)
@ -33,13 +39,11 @@ def test_discrete():
d3 = openmc.stats.Discrete(vals, probs)
# sample discrete distribution
# sample discrete distribution and check that the mean of the samples is
# within 3 std. dev. of the expected mean
n_samples = 1_000_000
samples = d3.sample(n_samples, seed=100)
# check that the mean of the samples is within 3 std. dev.
# of the expected mean
std_dev = samples.std() / np.sqrt(n_samples)
assert np.abs(exp_mean - samples.mean()) < 3*std_dev
assert_sample_mean(samples, exp_mean)
def test_merge_discrete():
@ -76,17 +80,16 @@ def test_uniform():
assert len(d) == 2
t = d.to_tabular()
assert t.x == [a, b]
assert t.p == [1/(b-a), 1/(b-a)]
np.testing.assert_array_equal(t.x, [a, b])
np.testing.assert_array_equal(t.p, [1/(b-a), 1/(b-a)])
assert t.interpolation == 'histogram'
# Sample distribution and check that the mean of the samples is within 3
# std. dev. of the expected mean
exp_mean = 0.5 * (a + b)
n_samples = 1_000_000
samples = d.sample(n_samples, seed=100)
# check that the mean of the samples is within 3 std. dev.
# of the expected mean
std_dev = samples.std() / np.sqrt(n_samples)
assert np.abs(exp_mean - samples.mean()) < 3*std_dev
assert_sample_mean(samples, exp_mean)
def test_powerlaw():
@ -102,13 +105,11 @@ def test_powerlaw():
exp_mean = 100.0 * (n+1) / (n+2)
# sample power law distribution
# sample power law distribution and check that the mean of the samples is
# within 3 std. dev. of the expected mean
n_samples = 1_000_000
samples = d.sample(n_samples, seed=100)
# check that the mean of the samples is within 3 std. dev.
# of the expected mean
std_dev = samples.std() / np.sqrt(n_samples)
assert np.abs(exp_mean - samples.mean()) < 3*std_dev
assert_sample_mean(samples, exp_mean)
def test_maxwell():
@ -122,21 +123,15 @@ def test_maxwell():
exp_mean = 3/2 * theta
# sample maxwell distribution
# sample maxwell distribution and check that the mean of the samples is
# within 3 std. dev. of the expected mean
n_samples = 1_000_000
samples = d.sample(n_samples, seed=100)
# check that the mean of the samples is within 3 std. dev.
# of the expected mean
std_dev = samples.std() / np.sqrt(n_samples)
assert np.abs(exp_mean - samples.mean()) < 3*std_dev
assert_sample_mean(samples, exp_mean)
# A second sample with a different seed
samples_2 = d.sample(n_samples, seed=200)
# check that the mean of the samples is within 3 std. dev.
# of the expected mean
std_dev = samples_2.std() / np.sqrt(n_samples)
assert np.abs(exp_mean - samples_2.mean()) < 3*std_dev
assert_sample_mean(samples_2, exp_mean)
assert samples_2.mean() != samples.mean()
@ -156,13 +151,11 @@ def test_watt():
# https://doi.org/10.1016/j.physletb.2003.09.048
exp_mean = 3/2 * a + a**2 * b / 4
# sample Watt distribution
# sample Watt distribution and check that the mean of the samples is within
# 3 std. dev. of the expected mean
n_samples = 1_000_000
samples = d.sample(n_samples, seed=100)
# check that the mean of the samples is within 3 std. dev.
# of the expected mean
std_dev = samples.std() / np.sqrt(n_samples)
assert np.abs(exp_mean - samples.mean()) < 3*std_dev
assert_sample_mean(samples, exp_mean)
def test_tabular():
@ -228,6 +221,11 @@ def test_mixture():
assert mix.distribution == [d1, d2]
assert len(mix) == 4
# Sample and make sure sample mean is close to expected mean
n_samples = 1_000_000
samples = mix.sample(n_samples)
assert_sample_mean(samples, (2.5 + 5.0)/2)
elem = mix.to_xml_element('distribution')
d = openmc.stats.Mixture.from_xml_element(elem)
@ -396,3 +394,47 @@ def test_muir():
assert within_2_sigma / n_samples >= 0.95
within_3_sigma = np.count_nonzero(samples < 3*d.std_dev)
assert within_3_sigma / n_samples >= 0.99
def test_combine_distributions():
# Combine two discrete (same data as in test_merge_discrete)
x1 = [0.0, 1.0, 10.0]
p1 = [0.3, 0.2, 0.5]
d1 = openmc.stats.Discrete(x1, p1)
x2 = [0.5, 1.0, 5.0]
p2 = [0.4, 0.5, 0.1]
d2 = openmc.stats.Discrete(x2, p2)
# Merged distribution should have x values sorted and probabilities
# appropriately combined. Duplicate x values should appear once.
merged = openmc.stats.combine_distributions([d1, d2], [0.6, 0.4])
assert isinstance(merged, openmc.stats.Discrete)
assert merged.x == pytest.approx([0.0, 0.5, 1.0, 5.0, 10.0])
assert merged.p == pytest.approx(
[0.6*0.3, 0.4*0.4, 0.6*0.2 + 0.4*0.5, 0.4*0.1, 0.6*0.5])
# Probabilities add up but are not normalized
d1 = openmc.stats.Discrete([3.0], [1.0])
triple = openmc.stats.combine_distributions([d1, d1, d1], [1.0, 2.0, 3.0])
assert triple.x == pytest.approx([3.0])
assert triple.p == pytest.approx([6.0])
# Combine discrete and tabular
t1 = openmc.stats.Tabular(x2, p2)
mixed = openmc.stats.combine_distributions([d1, t1], [0.5, 0.5])
assert isinstance(mixed, openmc.stats.Mixture)
assert len(mixed.distribution) == 2
assert len(mixed.probability) == 2
# Combine 1 discrete and 2 tabular -- the tabular distributions should
# combine to produce a uniform distribution with mean 0.5. The combined
# distribution should have a mean of 0.25.
t1 = openmc.stats.Tabular([0., 1.], [2.0, 0.0])
t2 = openmc.stats.Tabular([0., 1.], [0.0, 2.0])
d1 = openmc.stats.Discrete([0.0], [1.0])
combined = openmc.stats.combine_distributions([t1, t2, d1], [0.25, 0.25, 0.5])
# Sample the combined distribution and make sure the sample mean is within
# uncertainty of the expected value
samples = combined.sample(1000)
assert_sample_mean(samples, 0.25)