crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #1628 from paulromano/ccfe-spp-2

Browse source 
Initial implementation of fixed-source activation calculations
This commit is contained in:
Patrick Shriwise 2020-09-10 14:19:56 -05:00 • committed by GitHub
commit 364422fea8
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
15 changed files with 613 additions and 895 deletions
Download patch file Download diff file Expand all files Collapse all files

9
docs/source/io_formats/depletion_results.rst
View file

@ -4,7 +4,7 @@
Depletion Results File Format
=============================
The current version of the depletion results file format is 1.0.
The current version of the depletion results file format is 1.1.
**/**
@ -14,7 +14,7 @@ The current version of the depletion results file format is 1.0.
:Datasets: - **eigenvalues** (*double[][][2]*) -- k-eigenvalues at each
time/stage. This array has shape (number of timesteps, number of
stages, value). The last axis contains the eigenvalue and the
stages, value). The last axis contains the eigenvalue and the
associated uncertainty
- **number** (*double[][][][]*) -- Total number of atoms. This array
has shape (number of timesteps, number of stages, number of
@ -25,7 +25,10 @@ The current version of the depletion results file format is 1.0.
nuclides, number of reactions).
- **time** (*double[][2]*) -- Time in [s] at beginning/end of each
step.
- **depletion time** (*double[]*) -- Average process time in [s]
- **source_rate** (*double[][]*) -- Power in [W] or source rate in
[neutron/sec]. This array has shape (number of timesteps, number
of stages).
- **depletion time** (*double[]*) -- Average process time in [s]
spent depleting a material across all burnable materials and,
if applicable, MPI processes.

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

@ -243,5 +243,5 @@ user to choose one of two methods for estimating the heating rate, including:
2. Using the ``heating`` or ``heating-local`` scores to obtain an nuclide- and
energy-dependent estimate of the true heating rate.
The method for normalization can be chosen through the ``energy_mode`` argument
to the :class:`openmc.deplete.Operator` class.
The method for normalization can be chosen through the ``normalization_mode``
argument to the :class:`openmc.deplete.Operator` class.

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

@ -207,7 +207,7 @@ OpenMC simulations back on to the :class:`abc.TransportOperator`
:nosignatures:
:template: myclass.rst
abc.EnergyHelper
abc.NormalizationHelper
abc.FissionYieldHelper
abc.ReactionRateHelper
abc.TalliedFissionYieldHelper

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

@ -1,8 +1,8 @@
.. _usersguide_depletion:
=========
Depletion
=========
===========================
Depletion and Transmutation
===========================
OpenMC supports coupled depletion, or burnup, calculations through the
:mod:`openmc.deplete` Python module. OpenMC solves the transport equation to
@ -28,6 +28,8 @@ operator class requires a :class:`openmc.Geometry` instance and a
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.
: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
@ -51,6 +53,38 @@ time::
Note that the coupling between the transport solver and the transmutation solver
happens in-memory rather than by reading/writing files on disk.
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,
and normalization of tally results is performed based on energy deposition. It
is also possible to run a fixed-source simulation and perform normalization
based on a known source rate. First, as with all fixed-source calculations, we
need to set the run mode::
settings.run_mode = 'fixed source'
Additionally, all materials that you wish to deplete need to be marked as such
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::
op = openmc.deplete.Operator(geometry, settings, normalization_mode='source-rate')
Finally, when creating a depletion integrator, use the ``source_rates`` argument::
integrator = openmc.deplete.PredictorIntegrator(op, timesteps, sources_rates=...)
As with the ``power`` argument, you can provide a different source rate for each
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
=======
@ -61,14 +95,14 @@ 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.
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 should
be, including indirect components. Some examples are provided below::
variety of ways. This requires knowing what the total fission energy release
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
@ -83,39 +117,39 @@ be, including indirect components. Some examples are provided below::
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 to normalize
reaction rates instead of using the fission reaction rates with::
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
to normalize reaction rates instead of using the fission reaction rates with::
op = openmc.deplete.Operator(geometry, settings, "chain.xml",
energy_mode="energy-deposition")
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 into
the distributed libraries.
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 creating a single
3.5 wt.% enriched fuel ``fuel_3``, and placing that fuel in every fuel pin in an assembly
or even full core problem. This certainly expedites the model making process, but can pose
issues with depletion.
Under this setup, :mod:`openmc.deplete` will deplete a single ``fuel_3`` material using
a single set of reaction rates, and produce a single new composition for the next time
step. This can be problematic if the same ``fuel_3`` is used in very different regions
of the problem.
It is not uncommon to explicitly create a single burnable material across many
locations. From a pure transport perspective, there is nothing wrong with
creating a single 3.5 wt.% enriched fuel ``fuel_3``, and placing that fuel in
every fuel pin in an assembly or even full core problem. This certainly
expedites the model making process, but can pose issues with depletion. Under
this setup, :mod:`openmc.deplete` will deplete a single ``fuel_3`` material
using a single set of reaction rates, and produce a single new composition for
the next time step. This can be problematic if the same ``fuel_3`` is used in
very different regions of the problem.
As an example, consider a full-scale power reactor core with vacuum boundary
conditions, and with fuel pins solely composed of the same ``fuel_3`` material.
The fuel pins towards the center of the problem will surely experience a more intense
neutron flux and greater reaction rates than those towards the edge of the domain.
This indicates that the fuel in the center should be at a more depleted state than
periphery pins, at least for the fist depletion step.
However, without any other instructions, OpenMC will deplete ``fuel_3`` as a single
material, and all of the fuel pins will have an identical composition at the next
transport step.
The fuel pins towards the center of the problem will surely experience a more
intense neutron flux and greater reaction rates than those towards the edge of
the domain. This indicates that the fuel in the center should be at a more
depleted state than periphery pins, at least for the fist depletion step.
However, without any other instructions, OpenMC will deplete ``fuel_3`` as a
single material, and all of the fuel pins will have an identical composition at
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::
@ -123,12 +157,13 @@ of the same material as a unique material definition with::
op = openmc.deplete.Operator(geometry, settings, chain_file,
diff_burnable_mats=True)
For our example problem, this would deplete fuel on the outer region of the problem
with different reaction rates than those in the center. Materials will be depleted
corresponding to their local neutron spectra, and have unique compositions at each
transport step. The volume of the original ``fuel_3`` material must represent
the volume of **all** the ``fuel_3`` in the problem. When creating the unique
materials, this volume will be equally distributed across all material instances.
For our example problem, this would deplete fuel on the outer region of the
problem with different reaction rates than those in the center. Materials will
be depleted corresponding to their local neutron spectra, and have unique
compositions at each transport step. The volume of the original ``fuel_3``
material must represent the volume of **all** the ``fuel_3`` in the problem.
When creating the unique materials, this volume will be equally distributed
across all material instances.
.. note::

194
openmc/deplete/abc.py
View file

@ -4,23 +4,25 @@ This module contains the Operator class, which is then passed to an integrator
to run a full depletion simulation.
"""
from abc import ABC, abstractmethod
from collections import namedtuple, defaultdict
from collections.abc import Iterable, Callable
from copy import deepcopy
from inspect import signature
from numbers import Real, Integral
import os
from pathlib import Path
from abc import ABC, abstractmethod
from copy import deepcopy
from warnings import warn
from numbers import Real, Integral
from inspect import signature
import sys
import time
from warnings import warn
from numpy import nonzero, empty, asarray
from uncertainties import ufloat
from openmc.data import DataLibrary, JOULE_PER_EV
from openmc.data import DataLibrary
from openmc.lib import MaterialFilter, Tally
from openmc.checkvalue import check_type, check_greater_than
from . import comm
from .results import Results
from .chain import Chain
from .results_list import ResultsList
@ -29,8 +31,8 @@ from .pool import deplete
__all__ = [
"OperatorResult", "TransportOperator", "ReactionRateHelper",
"EnergyHelper", "FissionYieldHelper", "TalliedFissionYieldHelper",
"Integrator", "SIIntegrator", "DepSystemSolver"]
"NormalizationHelper", "FissionYieldHelper", "TalliedFissionYieldHelper",
"Integrator", "SIIntegrator", "DepSystemSolver", "add_params"]
_SECONDS_PER_MINUTE = 60
@ -136,15 +138,15 @@ class TransportOperator(ABC):
self._dilute_initial = value
@abstractmethod
def __call__(self, vec, power):
def __call__(self, vec, source_rate):
"""Runs a simulation.
Parameters
----------
vec : list of numpy.ndarray
Total atoms to be used in function.
power : float
Power of the reactor in [W]
source_rate : float
Power in [W] or source rate in [neutron/sec]
Returns
-------
@ -302,14 +304,14 @@ class ReactionRateHelper(ABC):
return results
class EnergyHelper(ABC):
"""Abstract class for obtaining energy produced
class NormalizationHelper(ABC):
"""Abstract class for obtaining normalization factor on tallies
The ultimate goal of this helper is to provide instances of
:class:`openmc.deplete.Operator` with the total energy produced
in a transport simulation. This information, provided with the
power requested by the user and reaction rates from a
:class:`ReactionRateHelper` will scale reaction rates to the
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:`ReactionRateHelper`, this class will scale reaction rates to the
correct values.
Attributes
@ -317,29 +319,22 @@ class EnergyHelper(ABC):
nuclides : list of str
All nuclides with desired reaction rates. Ordered to be
consistent with :class:`openmc.deplete.Operator`
energy : float
Total energy [J/s/source neutron] produced in a transport simulation.
Updated in the material iteration with :meth:`update`.
"""
def __init__(self):
self._nuclides = None
self._energy = 0.0
@property
def energy(self):
return self._energy * JOULE_PER_EV
def reset(self):
"""Reset energy produced prior to unpacking tallies"""
self._energy = 0.0
"""Reset state for normalization"""
@abstractmethod
def prepare(self, chain_nucs, rate_index, materials):
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`.
in :meth:`openmc.deplete.Operator.initial_condition`. Only used for
energy-based normalization.
Parameters
----------
@ -348,14 +343,11 @@ class EnergyHelper(ABC):
rate_index : dict of str to int
Mapping from nuclide name to index in the
`fission_rates` for :meth:`update`.
materials : list of str
All materials tracked on the operator helped by this
object. Should correspond to
:attr:`openmc.deplete.Operator.burnable_materials`
"""
def update(self, fission_rates, mat_index):
"""Update the energy produced
def update(self, fission_rates):
"""Update the normalization based on fission rates (only used for
energy-based normalization)
Parameters
----------
@ -363,9 +355,6 @@ class EnergyHelper(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
Index for the specific material in the list of all burnable
materials.
"""
@property
@ -378,6 +367,22 @@ class EnergyHelper(ABC):
check_type("nuclides", nuclides, list, str)
self._nuclides = nuclides
@abstractmethod
def factor(self, source_rate):
"""Return normalization factor
Parameters
----------
source_rate : float
Power in [W] or source rate in [neutron/sec]
Returns
-------
float
Normalization factor for tallies
"""
class FissionYieldHelper(ABC):
"""Abstract class for processing energy dependent fission yields
@ -596,9 +601,17 @@ class TalliedFissionYieldHelper(FissionYieldHelper):
"""
def add_params(cls):
cls.__doc__ += cls._params
return cls
@add_params
class Integrator(ABC):
r"""Abstract class for solving the time-integration for depletion
"""
_params = r"""
Parameters
----------
operator : openmc.deplete.TransportOperator
@ -614,12 +627,16 @@ class Integrator(ABC):
indicates potentially different power levels for each timestep.
For a 2D problem, the power can be given in [W/cm] as long
as the "volume" assigned to a depletion material is actually
an area in [cm^2]. Either ``power`` or ``power_density`` must be
specified.
an area in [cm^2]. Either ``power``, ``power_density``, or
``source_rates`` must be specified.
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.
source_rates : float or iterable of float, optional
Source rate in [neutron/sec] for each interval in :attr:`timesteps`
.. versionadded:: 0.12.1
timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
Units for values specified in the `timesteps` argument. 's' means
seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
@ -645,8 +662,9 @@ class Integrator(ABC):
Depletion chain
timesteps : iterable of float
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
source_rates : iterable of float
Source rate in [W] or [neutron/sec] for each interval in
:attr:`timesteps`
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
@ -667,7 +685,7 @@ class Integrator(ABC):
"""
def __init__(self, operator, timesteps, power=None, power_density=None,
timestep_units='s', solver="cram48"):
source_rates=None, timestep_units='s', solver="cram48"):
# Check number of stages previously used
if operator.prev_res is not None:
res = operator.prev_res[-1]
@ -681,22 +699,25 @@ class Integrator(ABC):
self.operator = operator
self.chain = operator.chain
# Determine power and normalize units to W
if power is None:
if power_density is None:
raise ValueError("Either power or power density must be set")
# Determine source rate and normalize units to W in using power
if power is not None:
source_rates = power
elif power_density is not None:
if not isinstance(power_density, Iterable):
power = power_density * operator.heavy_metal
source_rates = power_density * operator.heavy_metal
else:
power = [p*operator.heavy_metal for p in power_density]
if not isinstance(power, Iterable):
# Ensure that power is single value if that is the case
power = [power] * len(timesteps)
source_rates = [p*operator.heavy_metal for p in power_density]
elif source_rates is None:
raise ValueError("Either power, power_density, or source_rates must be set")
if len(power) != len(timesteps):
if not isinstance(source_rates, Iterable):
# Ensure that rate is single value if that is the case
source_rates = [source_rates] * len(timesteps)
if len(source_rates) != len(timesteps):
raise ValueError(
"Number of time steps ({}) != number of powers ({})".format(
len(timesteps), len(power)))
len(timesteps), len(source_rates)))
# Get list of times / units
if isinstance(timesteps[0], Iterable):
@ -707,13 +728,13 @@ class Integrator(ABC):
# Determine number of seconds for each timestep
seconds = []
for timestep, unit, watts in zip(times, units, power):
for timestep, unit, rate in zip(times, units, source_rates):
# Make sure values passed make sense
check_type('timestep', timestep, Real)
check_greater_than('timestep', timestep, 0.0, False)
check_type('timestep units', unit, str)
check_type('power', watts, Real)
check_greater_than('power', watts, 0.0, True)
check_type('source rate', rate, Real)
check_greater_than('source rate', rate, 0.0, True)
if unit in ('s', 'sec'):
seconds.append(timestep)
@ -726,13 +747,13 @@ class Integrator(ABC):
elif unit.lower() == 'mwd/kg':
watt_days_per_kg = 1e6*timestep
kilograms = 1e-3*operator.heavy_metal
days = watt_days_per_kg * kilograms / watts
days = watt_days_per_kg * kilograms / rate
seconds.append(days*_SECONDS_PER_DAY)
else:
raise ValueError("Invalid timestep unit '{}'".format(unit))
self.timesteps = asarray(seconds)
self.power = asarray(power)
self.source_rates = asarray(source_rates)
if isinstance(solver, str):
# Delay importing of cram module, which requires this file
@ -788,7 +809,7 @@ class Integrator(ABC):
return time.time() - start, results
@abstractmethod
def __call__(self, conc, rates, dt, power, i):
def __call__(self, conc, rates, dt, source_rate, i):
"""Perform the integration across one time step
Parameters
@ -799,8 +820,8 @@ class Integrator(ABC):
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system in [W]
source_rate : float
Power in [W] or source rate in [neutron/sec]
i : int
Current depletion step index
@ -825,22 +846,22 @@ class Integrator(ABC):
"""
def __iter__(self):
"""Return pairs of time steps in [s] and powers in [W]"""
return zip(self.timesteps, self.power)
"""Return pair of time step in [s] and source rate in [W] or [neutron/sec]"""
return zip(self.timesteps, self.source_rates)
def __len__(self):
"""Return integer number of depletion intervals"""
return len(self.timesteps)
def _get_bos_data_from_operator(self, step_index, step_power, bos_conc):
def _get_bos_data_from_operator(self, step_index, source_rate, bos_conc):
"""Get beginning of step concentrations, reaction rates from Operator
"""
x = deepcopy(bos_conc)
res = self.operator(x, step_power)
res = self.operator(x, source_rate)
self.operator.write_bos_data(step_index + self._i_res)
return x, res
def _get_bos_data_from_restart(self, step_index, step_power, bos_conc):
def _get_bos_data_from_restart(self, step_index, source_rate, bos_conc):
"""Get beginning of step concentrations, reaction rates from restart"""
res = self.operator.prev_res[-1]
# Depletion methods expect list of arrays
@ -848,8 +869,8 @@ class Integrator(ABC):
rates = res.rates[0]
k = ufloat(res.k[0, 0], res.k[0, 1])
# Scale rates by ratio of powers
rates *= step_power / res.power[0]
# Scale reaction rates by ratio of source rates
rates *= source_rate / res.source_rate[0]
return bos_conc, OperatorResult(k, rates)
def _get_start_data(self):
@ -873,15 +894,15 @@ class Integrator(ABC):
with self.operator as conc:
t, self._i_res = self._get_start_data()
for i, (dt, p) in enumerate(self):
for i, (dt, source_rate) in enumerate(self):
# Solve transport equation (or obtain result from restart)
if i > 0 or self.operator.prev_res is None:
conc, res = self._get_bos_data_from_operator(i, p, conc)
conc, res = self._get_bos_data_from_operator(i, source_rate, conc)
else:
conc, res = self._get_bos_data_from_restart(i, p, conc)
conc, res = self._get_bos_data_from_restart(i, source_rate, conc)
# Solve Bateman equations over time interval
proc_time, conc_list, res_list = self(conc, res.rates, dt, p, i)
proc_time, conc_list, res_list = self(conc, res.rates, dt, source_rate, i)
# Insert BOS concentration, transport results
conc_list.insert(0, conc)
@ -891,7 +912,7 @@ class Integrator(ABC):
conc = conc_list.pop()
Results.save(self.operator, conc_list, res_list, [t, t + dt],
p, self._i_res + i, proc_time)
source_rate, self._i_res + i, proc_time)
t += dt
@ -899,18 +920,21 @@ class Integrator(ABC):
# source rate is passed to the transport operator (which knows to
# just return zero reaction rates without actually doing a transport
# solve)
res_list = [self.operator(conc, p if final_step else 0.0)]
res_list = [self.operator(conc, source_rate if final_step else 0.0)]
Results.save(self.operator, [conc], res_list, [t, t],
p, self._i_res + len(self), proc_time)
source_rate, self._i_res + len(self), proc_time)
self.operator.write_bos_data(len(self) + self._i_res)
@add_params
class SIIntegrator(Integrator):
r"""Abstract class for the Stochastic Implicit Euler integrators
Does not provide a ``__call__`` method, but scales and resets
the number of particles used in initial transport calculation
"""
_params = r"""
Parameters
----------
operator : openmc.deplete.TransportOperator
@ -926,12 +950,16 @@ class SIIntegrator(Integrator):
indicates potentially different power levels for each timestep.
For a 2D problem, the power can be given in [W/cm] as long
as the "volume" assigned to a depletion material is actually
an area in [cm^2]. Either ``power`` or ``power_density`` must be
specified.
an area in [cm^2]. Either ``power``, ``power_density``, or
``source_rates`` must be specified.
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.
source_rates : float or iterable of float, optional
Source rate in [neutron/sec] for each interval in :attr:`timesteps`
.. versionadded:: 0.12.1
timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
Units for values specified in the `timesteps` argument. 's' means
seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
@ -982,13 +1010,15 @@ class SIIntegrator(Integrator):
.. versionadded:: 0.12
"""
def __init__(self, operator, timesteps, power=None, power_density=None,
timestep_units='s', n_steps=10, solver="cram48"):
source_rates=None, timestep_units='s', n_steps=10,
solver="cram48"):
check_type("n_steps", n_steps, Integral)
check_greater_than("n_steps", n_steps, 0)
super().__init__(
operator, timesteps, power, power_density, timestep_units,
solver=solver)
operator, timesteps, power, power_density, source_rates,
timestep_units=timestep_units, solver=solver)
self.n_steps = n_steps
def _get_bos_data_from_operator(self, step_index, step_power, bos_conc):

137
openmc/deplete/chain.py
View file

@ -284,6 +284,23 @@ class Chain:
"""Number of nuclides in chain."""
return len(self.nuclides)
def add_nuclide(self, nuclide):
"""Add a nuclide to the depletion chain
Parameters
----------
nuclide : openmc.deplete.Nuclide
Nuclide to add
"""
self.nuclide_dict[nuclide.name] = len(self.nuclides)
self.nuclides.append(nuclide)
# Check for reaction paths
for rx in nuclide.reactions:
if rx.type not in self.reactions:
self.reactions.append(rx.type)
@classmethod
def from_endf(cls, decay_files, fpy_files, neutron_files,
reactions=('(n,2n)', '(n,3n)', '(n,4n)', '(n,gamma)', '(n,p)', '(n,a)'),
@ -376,9 +393,6 @@ class Chain:
nuclide = Nuclide(parent)
chain.nuclides.append(nuclide)
chain.nuclide_dict[parent] = idx
if not data.nuclide['stable'] and data.half_life.nominal_value != 0.0:
nuclide.half_life = data.half_life.nominal_value
nuclide.decay_energy = sum(E.nominal_value for E in
@ -401,7 +415,7 @@ class Chain:
for m in data.modes[:-1])
# Append decay mode
nuclide.decay_modes.append(DecayTuple(type_, target, br))
nuclide.add_decay_mode(type_, target, br)
fissionable = False
if parent in reactions:
@ -414,9 +428,6 @@ class Chain:
Z = data.nuclide['atomic_number'] + delta_Z
daughter = '{}{}'.format(openmc.data.ATOMIC_SYMBOL[Z], A)
if name not in chain.reactions:
chain.reactions.append(name)
if daughter not in decay_data:
daughter = replace_missing(daughter, decay_data)
if daughter is None:
@ -430,49 +441,47 @@ class Chain:
else:
q_value = 0.0
nuclide.reactions.append(ReactionTuple(
name, daughter, q_value, 1.0))
nuclide.add_reaction(name, daughter, q_value, 1.0)
if any(mt in reactions_available for mt in openmc.data.FISSION_MTS):
q_value = reactions[parent][18]
nuclide.reactions.append(
ReactionTuple('fission', None, q_value, 1.0))
nuclide.add_reaction('fission', None, q_value, 1.0)
fissionable = True
if 'fission' not in chain.reactions:
chain.reactions.append('fission')
if fissionable:
if parent in fpy_data:
fpy = fpy_data[parent]
if fpy.energies is not None:
yield_energies = fpy.energies
else:
yield_energies = [0.0]
yield_data = {}
for E, yield_table in zip(yield_energies, fpy.independent):
yield_replace = 0.0
yields = defaultdict(float)
for product, y in yield_table.items():
# Handle fission products that have no decay data
if product not in decay_data:
daughter = replace_missing(product, decay_data)
product = daughter
yield_replace += y.nominal_value
yields[product] += y.nominal_value
if yield_replace > 0.0:
missing_fp.append((parent, E, yield_replace))
yield_data[E] = yields
nuclide.yield_data = FissionYieldDistribution(yield_data)
else:
nuclide._fpy = replace_missing_fpy(parent, fpy_data, decay_data)
missing_fpy.append((parent, nuclide._fpy))
continue
if fpy.energies is not None:
yield_energies = fpy.energies
else:
yield_energies = [0.0]
yield_data = {}
for E, yield_table in zip(yield_energies, fpy.independent):
yield_replace = 0.0
yields = defaultdict(float)
for product, y in yield_table.items():
# Handle fission products that have no decay data
if product not in decay_data:
daughter = replace_missing(product, decay_data)
product = daughter
yield_replace += y.nominal_value
yields[product] += y.nominal_value
if yield_replace > 0.0:
missing_fp.append((parent, E, yield_replace))
yield_data[E] = yields
nuclide.yield_data = FissionYieldDistribution(yield_data)
# Add nuclide to chain
chain.add_nuclide(nuclide)
# Replace missing FPY data
for nuclide in chain.nuclides:
@ -532,14 +541,7 @@ class Chain:
this_q = fission_q.get(nuclide_elem.get("name"))
nuc = Nuclide.from_xml(nuclide_elem, root, this_q)
chain.nuclide_dict[nuc.name] = i
# Check for reaction paths
for rx in nuc.reactions:
if rx.type not in chain.reactions:
chain.reactions.append(rx.type)
chain.nuclides.append(nuc)
chain.add_nuclide(nuc)
return chain
@ -876,8 +878,7 @@ class Chain:
all_meta = True
for target, br in new_ratios.items():
all_meta = all_meta and ("_m" in target)
parent.reactions.append(ReactionTuple(
reaction, target, rxn_Q, br))
parent.add_reaction(reaction, target, rxn_Q, br)
# If branching ratios don't add to unity, add reaction to ground
# with remainder of branching ratio
@ -888,8 +889,7 @@ class Chain:
pz, pa, pm = zam(parent_name)
ground_target = gnd_name(pz, pa + 1, 0)
new_ratios[ground_target] = ground_br
parent.reactions.append(ReactionTuple(
reaction, ground_target, rxn_Q, ground_br))
parent.add_reaction(reaction, ground_target, rxn_Q, ground_br)
@property
def fission_yields(self):
@ -1013,9 +1013,7 @@ class Chain:
# Avoid re-sorting for fission yields
name_sort = sorted(all_isotopes)
nuclides = []
nuclide_dict = {}
reactions = set()
new_chain = type(self)()
for idx, iso in enumerate(sorted(all_isotopes, key=openmc.data.zam)):
previous = self[iso]
@ -1025,44 +1023,29 @@ class Chain:
if hasattr(previous, '_fpy'):
new_nuclide._fpy = previous._fpy
new_decay = []
for mode in previous.decay_modes:
if mode.target in all_isotopes:
new_decay.append(mode)
new_nuclide.add_decay_mode(*mode)
else:
new_decay.append(DecayTuple(
mode.type, None, mode.branching_ratio))
new_nuclide.decay_modes = new_decay
new_nuclide.add_decay_mode(mode.type, None, mode.branching_ratio)
new_reactions = []
for rxn in previous.reactions:
if rxn.target in all_isotopes:
new_reactions.append(rxn)
reactions.add(rxn.type)
elif rxn.type == "fission":
for rx in previous.reactions:
if rx.target in all_isotopes:
new_nuclide.add_reaction(*rx)
elif rx.type == "fission":
new_yields = new_nuclide.yield_data = (
previous.yield_data.restrict_products(name_sort))
if new_yields is not None:
new_reactions.append(rxn)
reactions.add("fission")
new_nuclide.add_reaction(*rx)
# Maintain total destruction rates but set no target
else:
new_reactions.append(ReactionTuple(
rxn.type, None, rxn.Q, rxn.branching_ratio))
reactions.add(rxn.type)
new_nuclide.add_reaction(rx.type, None, rx.Q, rx.branching_ratio)
new_nuclide.reactions = new_reactions
nuclides.append(new_nuclide)
nuclide_dict[iso] = idx
new_chain = type(self)()
new_chain.nuclides = nuclides
new_chain.nuclide_dict = nuclide_dict
new_chain.add_nuclide(new_nuclide)
# Doesn't appear that the ordering matters for the reactions,
# just the contents
new_chain.reactions = sorted(reactions)
new_chain.reactions = sorted(new_chain.reactions)
return new_chain

74
openmc/deplete/helpers.py
View file

@ -11,15 +11,16 @@ from numpy import dot, zeros, newaxis
from . import comm
from openmc.checkvalue import check_type, check_greater_than
from openmc.data import JOULE_PER_EV
from openmc.lib import (
Tally, MaterialFilter, EnergyFilter, EnergyFunctionFilter)
from .abc import (
ReactionRateHelper, EnergyHelper, FissionYieldHelper,
ReactionRateHelper, NormalizationHelper, FissionYieldHelper,
TalliedFissionYieldHelper)
__all__ = (
"DirectReactionRateHelper", "ChainFissionHelper",
"ConstantFissionYieldHelper", "FissionYieldCutoffHelper",
"DirectReactionRateHelper", "ChainFissionHelper", "EnergyScoreHelper"
"SourceRateHelper", "ConstantFissionYieldHelper", "FissionYieldCutoffHelper",
"AveragedFissionYieldHelper")
# -------------------------------------
@ -83,7 +84,7 @@ class DirectReactionRateHelper(ReactionRateHelper):
reaction rates in this material
"""
self._results_cache.fill(0.0)
full_tally_res = self._rate_tally.results[mat_id, :, 1]
full_tally_res = self._rate_tally.mean[mat_id]
for i_tally, (i_nuc, i_react) in enumerate(
product(nuc_index, react_index)):
self._results_cache[i_nuc, i_react] = full_tally_res[i_tally]
@ -91,13 +92,39 @@ class DirectReactionRateHelper(ReactionRateHelper):
return self._results_cache
# ----------------------------
# Helpers for obtaining energy
# ----------------------------
# ------------------------------------------
# Helpers for obtaining normalization factor
# ------------------------------------------
class ChainFissionHelper(EnergyHelper):
"""Computes energy using fission Q values from depletion chain
class EnergyNormalizationHelper(NormalizationHelper):
"""Compute energy-based normalization."""
def reset(self):
"""Reset energy produced prior to unpacking tallies"""
self._energy = 0.0
def factor(self, source_rate):
# Reduce energy produced from all processes
# J / source neutron
energy = comm.allreduce(self._energy) * JOULE_PER_EV
# Guard against divide by zero
if energy == 0:
if comm.rank == 0:
sys.stderr.flush()
print("No energy reported from OpenMC tallies. Do your HDF5 "
"files have heating data?\n", file=sys.stderr, flush=True)
comm.barrier()
comm.Abort(1)
# Return normalization factor for scaling reaction rates. In this case,
# the source rate is the power in [W], so [W] / [J/src] = [src/s]
return source_rate / energy
class ChainFissionHelper(EnergyNormalizationHelper):
"""Computes normalization using fission Q values from depletion chain
Attributes
----------
@ -113,7 +140,7 @@ class ChainFissionHelper(EnergyHelper):
super().__init__()
self._fission_q_vector = None
def prepare(self, chain_nucs, rate_index, _materials):
def prepare(self, chain_nucs, rate_index):
"""Populate the fission Q value vector from a chain.
Parameters
@ -125,8 +152,6 @@ class ChainFissionHelper(EnergyHelper):
Dictionary mapping names of nuclides, e.g. ``"U235"``,
to a corresponding index in the desired fission Q
vector.
_materials : list of str
Unused. Materials to be tracked for this helper.
"""
if (self._fission_q_vector is not None
and self._fission_q_vector.shape == (len(rate_index),)):
@ -143,7 +168,7 @@ class ChainFissionHelper(EnergyHelper):
self._fission_q_vector = fission_qs
def update(self, fission_rates, _mat_index):
def update(self, fission_rates):
"""Update energy produced with fission rates in a material
Parameters
@ -152,14 +177,11 @@ class ChainFissionHelper(EnergyHelper):
fission reaction rate for each isotope in the specified
material. Should be ordered corresponding to initial
``rate_index`` used in :meth:`prepare`
_mat_index : int
index for the material requested. Unused, as identical
isotopes in all materials have the same Q value.
"""
self._energy += dot(fission_rates, self._fission_q_vector)
class EnergyScoreHelper(EnergyHelper):
class EnergyScoreHelper(EnergyNormalizationHelper):
"""Class responsible for obtaining system energy via a tally score
Parameters
@ -172,7 +194,7 @@ class EnergyScoreHelper(EnergyHelper):
----------
nuclides : list of str
List of nuclides with reaction rates. Not needed, but provided
for a consistent API across other :class:`EnergyHelper`
for a consistent API across other :class:`NormalizationHelper`
energy : float
System energy [eV] computed from the tally. Will be zero for
all MPI processes that are not the "master" process to avoid
@ -208,7 +230,15 @@ class EnergyScoreHelper(EnergyHelper):
"""
super().reset()
if comm.rank == 0:
self._energy = self._tally.results[0, 0, 1]
self._energy = self._tally.mean[0, 0]
class SourceRateHelper(NormalizationHelper):
def prepare(self, *args, **kwargs):
pass
def factor(self, source_rate):
return source_rate
# ------------------------------------
# Helper for collapsing fission yields
@ -447,7 +477,7 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper):
if not self._tally_nucs or self._local_indexes.size == 0:
self.results = None
return
fission_rates = self._fission_rate_tally.results[..., 1].reshape(
fission_rates = self._fission_rate_tally.mean.reshape(
self.n_bmats, 2, len(self._tally_nucs))
self.results = fission_rates[self._local_indexes]
total_fission = self.results.sum(axis=1)
@ -609,9 +639,9 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper):
self.results = None
return
fission_results = (
self._fission_rate_tally.results[self._local_indexes, :, 1])
self._fission_rate_tally.mean[self._local_indexes])
self.results = (
self._weighted_tally.results[self._local_indexes, :, 1]).copy()
self._weighted_tally.mean[self._local_indexes]).copy()
nz_mat, nz_nuc = fission_results.nonzero()
self.results[nz_mat, nz_nuc] /= fission_results[nz_mat, nz_nuc]

636
openmc/deplete/integrators.py
View file

@ -1,7 +1,7 @@
import copy
from itertools import repeat
from .abc import Integrator, SIIntegrator, OperatorResult
from .abc import Integrator, SIIntegrator, OperatorResult, add_params
from ._matrix_funcs import (
cf4_f1, cf4_f2, cf4_f3, cf4_f4, celi_f1, celi_f2,
leqi_f1, leqi_f2, leqi_f3, leqi_f4, rk4_f1, rk4_f4
@ -13,6 +13,7 @@ __all__ = [
"SICELIIntegrator", "SILEQIIntegrator"]
@add_params
class PredictorIntegrator(Integrator):
r"""Deplete using a first-order predictor algorithm.
@ -25,78 +26,10 @@ class PredictorIntegrator(Integrator):
A_p &= A(y_n, t_n) \\
y_{n+1} &= \text{expm}(A_p h) y_n
\end{aligned}
Parameters
----------
operator : openmc.deplete.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
iterable of float. Alternatively, units can be specified for each step
by passing an iterable of (value, unit) tuples.
power : float or iterable of float, optional
Power of the reactor in [W]. A single value indicates that
the power is constant over all timesteps. An iterable
indicates potentially different power levels for each timestep.
For a 2D problem, the power can be given in [W/cm] as long
as the "volume" assigned to a depletion material is actually
an area in [cm^2]. Either ``power`` or ``power_density`` must be
specified.
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.
timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
Units for values specified in the `timesteps` argument. 's' means
seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
that the values are given in burnup (MW-d of energy deposited per
kilogram of initial heavy metal).
.. versionadded:: 0.12
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
----------
operator : openmc.deplete.TransportOperator
Operator to perform transport simulations
chain : openmc.deplete.Chain
Depletion chain
timesteps : iterable of float
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
_num_stages = 1
def __call__(self, conc, rates, dt, power, _i=None):
def __call__(self, conc, rates, dt, source_rate, _i=None):
"""Perform the integration across one time step
Parameters
@ -107,8 +40,8 @@ class PredictorIntegrator(Integrator):
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system in [W]
source_rate : float
Power in [W] or source rate in [neutron/sec]
_i : int or None
Iteration index. Not used
@ -127,6 +60,7 @@ class PredictorIntegrator(Integrator):
return proc_time, [conc_end], []
@add_params
class CECMIntegrator(Integrator):
r"""Deplete using the CE/CM algorithm.
@ -144,78 +78,10 @@ class CECMIntegrator(Integrator):
A_c &= A(y_m, t_n + h/2) \\
y_{n+1} &= \text{expm}(A_c h) y_n
\end{aligned}
Parameters
----------
operator : openmc.deplete.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
iterable of float. Alternatively, units can be specified for each step
by passing an iterable of (value, unit) tuples.
power : float or iterable of float, optional
Power of the reactor in [W]. A single value indicates that
the power is constant over all timesteps. An iterable
indicates potentially different power levels for each timestep.
For a 2D problem, the power can be given in [W/cm] as long
as the "volume" assigned to a depletion material is actually
an area in [cm^2]. Either ``power`` or ``power_density`` must be
specified.
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.
timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
Units for values specified in the `timesteps` argument. 's' means
seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
that the values are given in burnup (MW-d of energy deposited per
kilogram of initial heavy metal).
.. versionadded:: 0.12
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
----------
operator : openmc.deplete.TransportOperator
Operator to perform transport simulations
chain : openmc.deplete.Chain
Depletion chain
timesteps : iterable of float
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
_num_stages = 2
def __call__(self, conc, rates, dt, power, _i=None):
def __call__(self, conc, rates, dt, source_rate, _i=None):
"""Integrate using CE/CM
Parameters
@ -226,8 +92,8 @@ class CECMIntegrator(Integrator):
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system [W]
source_rate : float
Power in [W] or source rate in [neutron/sec]
_i : int, optional
Current iteration count. Not used
@ -243,7 +109,7 @@ class CECMIntegrator(Integrator):
"""
# deplete across first half of inteval
time0, x_middle = self._timed_deplete(conc, rates, dt / 2)
res_middle = self.operator(x_middle, power)
res_middle = self.operator(x_middle, source_rate)
# deplete across entire interval with BOS concentrations,
# MOS reaction rates
@ -252,6 +118,7 @@ class CECMIntegrator(Integrator):
return time0 + time1, [x_middle, x_end], [res_middle]
@add_params
class CF4Integrator(Integrator):
r"""Deplete using the CF4 algorithm.
@ -271,78 +138,10 @@ class CF4Integrator(Integrator):
y_4 &= \text{expm}( 1/4 F_1 + 1/6 F_2 + 1/6 F_3 - 1/12 F_4)
\text{expm}(-1/12 F_1 + 1/6 F_2 + 1/6 F_3 + 1/4 F_4) y_0
\end{aligned}
Parameters
----------
operator : openmc.deplete.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
iterable of float. Alternatively, units can be specified for each step
by passing an iterable of (value, unit) tuples.
power : float or iterable of float, optional
Power of the reactor in [W]. A single value indicates that
the power is constant over all timesteps. An iterable
indicates potentially different power levels for each timestep.
For a 2D problem, the power can be given in [W/cm] as long
as the "volume" assigned to a depletion material is actually
an area in [cm^2]. Either ``power`` or ``power_density`` must be
specified.
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.
timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
Units for values specified in the `timesteps` argument. 's' means
seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
that the values are given in burnup (MW-d of energy deposited per
kilogram of initial heavy metal).
.. versionadded:: 0.12
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
----------
operator : openmc.deplete.TransportOperator
Operator to perform transport simulations
chain : openmc.deplete.Chain
Depletion chain
timesteps : iterable of float
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
_num_stages = 4
def __call__(self, bos_conc, bos_rates, dt, power, _i=None):
def __call__(self, bos_conc, bos_rates, dt, source_rate, _i=None):
"""Perform the integration across one time step
Parameters
@ -353,8 +152,8 @@ class CF4Integrator(Integrator):
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system in [W]
source_rate : float
Power in [W] or source rate in [neutron/sec]
_i : int, optional
Current depletion step index. Not used
@ -372,18 +171,18 @@ class CF4Integrator(Integrator):
# Step 1: deplete with matrix 1/2*A(y0)
time1, conc_eos1 = self._timed_deplete(
bos_conc, bos_rates, dt, matrix_func=cf4_f1)
res1 = self.operator(conc_eos1, power)
res1 = self.operator(conc_eos1, source_rate)
# Step 2: deplete with matrix 1/2*A(y1)
time2, conc_eos2 = self._timed_deplete(
bos_conc, res1.rates, dt, matrix_func=cf4_f1)
res2 = self.operator(conc_eos2, power)
res2 = self.operator(conc_eos2, source_rate)
# Step 3: deplete with matrix -1/2*A(y0)+A(y2)
list_rates = list(zip(bos_rates, res2.rates))
time3, conc_eos3 = self._timed_deplete(
conc_eos1, list_rates, dt, matrix_func=cf4_f2)
res3 = self.operator(conc_eos3, power)
res3 = self.operator(conc_eos3, source_rate)
# Step 4: deplete with two matrix exponentials
list_rates = list(zip(bos_rates, res1.rates, res2.rates, res3.rates))
@ -397,6 +196,7 @@ class CF4Integrator(Integrator):
[res1, res2, res3])
@add_params
class CELIIntegrator(Integrator):
r"""Deplete using the CE/LI CFQ4 algorithm.
@ -415,78 +215,10 @@ class CELIIntegrator(Integrator):
y_{n+1} &= \text{expm}(\frac{h}{12} A_0 + \frac{5h}{12} A1)
\text{expm}(\frac{5h}{12} A_0 + \frac{h}{12} A1) y_n
\end{aligned}
Parameters
----------
operator : openmc.deplete.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
iterable of float. Alternatively, units can be specified for each step
by passing an iterable of (value, unit) tuples.
power : float or iterable of float, optional
Power of the reactor in [W]. A single value indicates that
the power is constant over all timesteps. An iterable
indicates potentially different power levels for each timestep.
For a 2D problem, the power can be given in [W/cm] as long
as the "volume" assigned to a depletion material is actually
an area in [cm^2]. Either ``power`` or ``power_density`` must be
specified.
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.
timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
Units for values specified in the `timesteps` argument. 's' means
seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
that the values are given in burnup (MW-d of energy deposited per
kilogram of initial heavy metal).
.. versionadded:: 0.12
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
----------
operator : openmc.deplete.TransportOperator
Operator to perform transport simulations
chain : openmc.deplete.Chain
Depletion chain
timesteps : iterable of float
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
_num_stages = 2
def __call__(self, bos_conc, rates, dt, power, _i=None):
def __call__(self, bos_conc, rates, dt, source_rate, _i=None):
"""Perform the integration across one time step
Parameters
@ -497,8 +229,8 @@ class CELIIntegrator(Integrator):
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system in [W]
source_rate : float
Power in [W] or source rate in [neutron/sec]
_i : int, optional
Current iteration count. Not used
@ -515,7 +247,7 @@ class CELIIntegrator(Integrator):
"""
# deplete to end using BOS rates
proc_time, conc_ce = self._timed_deplete(bos_conc, rates, dt)
res_ce = self.operator(conc_ce, power)
res_ce = self.operator(conc_ce, source_rate)
# deplete using two matrix exponentials
list_rates = list(zip(rates, res_ce.rates))
@ -529,6 +261,7 @@ class CELIIntegrator(Integrator):
return proc_time + time_le1 + time_le1, [conc_ce, conc_end], [res_ce]
@add_params
class EPCRK4Integrator(Integrator):
r"""Deplete using the EPC-RK4 algorithm.
@ -546,73 +279,10 @@ class EPCRK4Integrator(Integrator):
F_4 &= h A(y_3) \\
y_4 &= \text{expm}(1/6 F_1 + 1/3 F_2 + 1/3 F_3 + 1/6 F_4) y_0
\end{aligned}
Parameters
----------
operator : openmc.deplete.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
iterable of float. Alternatively, units can be specified for each step
by passing an iterable of (value, unit) tuples.
power : float or iterable of float, optional
Power of the reactor in [W]. A single value indicates that
the power is constant over all timesteps. An iterable
indicates potentially different power levels for each timestep.
For a 2D problem, the power can be given in [W/cm] as long
as the "volume" assigned to a depletion material is actually
an area in [cm^2]. Either ``power`` or ``power_density`` must be
specified.
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.
timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
Units for values specified in the `timesteps` argument. 's' means
seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
that the values are given in burnup (MW-d of energy deposited per
kilogram of initial heavy metal).
.. versionadded:: 0.12
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
----------
operator : openmc.deplete.TransportOperator
Operator to perform transport simulations
chain : openmc.deplete.Chain
Depletion chain
timesteps : iterable of float
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
"""
_num_stages = 4
def __call__(self, conc, rates, dt, power, _i=None):
def __call__(self, conc, rates, dt, source_rate, _i=None):
"""Perform the integration across one time step
Parameters
@ -623,8 +293,8 @@ class EPCRK4Integrator(Integrator):
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system in [W]
source_rate : float
Power in [W] or source rate in [neutron/sec]
_i : int, optional
Current depletion step index, unused.
@ -643,16 +313,16 @@ class EPCRK4Integrator(Integrator):
# Step 1: deplete with matrix A(y0) / 2
time1, conc1 = self._timed_deplete(
conc, rates, dt, matrix_func=rk4_f1)
res1 = self.operator(conc1, power)
res1 = self.operator(conc1, source_rate)
# Step 2: deplete with matrix A(y1) / 2
time2, conc2 = self._timed_deplete(
conc, res1.rates, dt, matrix_func=rk4_f1)
res2 = self.operator(conc2, power)
res2 = self.operator(conc2, source_rate)
# Step 3: deplete with matrix A(y2)
time3, conc3 = self._timed_deplete(conc, res2.rates, dt)
res3 = self.operator(conc3, power)
res3 = self.operator(conc3, source_rate)
# Step 4: deplete with matrix built from weighted rates
list_rates = list(zip(rates, res1.rates, res2.rates, res3.rates))
@ -663,6 +333,7 @@ class EPCRK4Integrator(Integrator):
[res1, res2, res3])
@add_params
class LEQIIntegrator(Integrator):
r"""Deplete using the LE/QI CFQ4 algorithm.
@ -691,78 +362,10 @@ class LEQIIntegrator(Integrator):
\end{aligned}
It is initialized using the CE/LI algorithm.
Parameters
----------
operator : openmc.deplete.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
iterable of float. Alternatively, units can be specified for each step
by passing an iterable of (value, unit) tuples.
power : float or iterable of float, optional
Power of the reactor in [W]. A single value indicates that
the power is constant over all timesteps. An iterable
indicates potentially different power levels for each timestep.
For a 2D problem, the power can be given in [W/cm] as long
as the "volume" assigned to a depletion material is actually
an area in [cm^2]. Either ``power`` or ``power_density`` must be
specified.
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.
timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
Units for values specified in the `timesteps` argument. 's' means
seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
that the values are given in burnup (MW-d of energy deposited per
kilogram of initial heavy metal).
.. versionadded:: 0.12
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
----------
operator : openmc.deplete.TransportOperator
Operator to perform transport simulations
chain : openmc.deplete.Chain
Depletion chain
timesteps : iterable of float
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
_num_stages = 2
def __call__(self, bos_conc, bos_rates, dt, power, i):
def __call__(self, bos_conc, bos_rates, dt, source_rate, i):
"""Perform the integration across one time step
Parameters
@ -773,8 +376,8 @@ class LEQIIntegrator(Integrator):
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system in [W]
source_rate : float
Power in [W] or source rate in [neutron/sec]
i : int
Current depletion step index
@ -793,7 +396,7 @@ class LEQIIntegrator(Integrator):
if self._i_res < 1: # need at least previous transport solution
self._prev_rates = bos_rates
return CELIIntegrator.__call__(
self, bos_conc, bos_rates, dt, power, i)
self, bos_conc, bos_rates, dt, source_rate, i)
prev_res = self.operator.prev_res[-2]
prev_dt = self.timesteps[i] - prev_res.time[0]
self._prev_rates = prev_res.rates[0]
@ -801,7 +404,7 @@ class LEQIIntegrator(Integrator):
prev_dt = self.timesteps[i - 1]
# Remaining LE/QI
bos_res = self.operator(bos_conc, power)
bos_res = self.operator(bos_conc, source_rate)
le_inputs = list(zip(
self._prev_rates, bos_res.rates, repeat(prev_dt), repeat(dt)))
@ -811,7 +414,7 @@ class LEQIIntegrator(Integrator):
time2, conc_eos0 = self._timed_deplete(
conc_inter, le_inputs, dt, matrix_func=leqi_f2)
res_inter = self.operator(conc_eos0, power)
res_inter = self.operator(conc_eos0, source_rate)
qi_inputs = list(zip(
self._prev_rates, bos_res.rates, res_inter.rates,
@ -830,6 +433,7 @@ class LEQIIntegrator(Integrator):
[bos_res, res_inter])
@add_params
class SICELIIntegrator(SIIntegrator):
r"""Deplete using the SI-CE/LI CFQ4 algorithm.
@ -839,83 +443,10 @@ class SICELIIntegrator(SIIntegrator):
Detailed algorithm can be found in section 3.2 in `Colin Josey's thesis
<http://hdl.handle.net/1721.1/113721>`_.
Parameters
----------
operator : openmc.deplete.TransportOperator
The operator object to simulate on.
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
iterable of float. Alternatively, units can be specified for each step
by passing an iterable of (value, unit) tuples.
power : float or iterable of float, optional
Power of the reactor in [W]. A single value indicates that
the power is constant over all timesteps. An iterable
indicates potentially different power levels for each timestep.
For a 2D problem, the power can be given in [W/cm] as long
as the "volume" assigned to a depletion material is actually
an area in [cm^2]. Either ``power`` or ``power_density`` must be
specified.
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.
timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
Units for values specified in the `timesteps` argument. 's' means
seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
that the values are given in burnup (MW-d of energy deposited per
kilogram of initial heavy metal).
.. versionadded:: 0.12
n_steps : int, optional
Number of stochastic iterations per depletion interval.
Must be greater than zero. Default : 10
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
----------
operator : openmc.deplete.TransportOperator
Operator to perform transport simulations
chain : openmc.deplete.Chain
Depletion chain
timesteps : iterable of float
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
n_steps : int
Number of stochastic iterations per depletion interval
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
_num_stages = 2
def __call__(self, bos_conc, bos_rates, dt, power, _i=None):
def __call__(self, bos_conc, bos_rates, dt, source_rate, _i=None):
"""Perform the integration across one time step
Parameters
@ -926,8 +457,8 @@ class SICELIIntegrator(SIIntegrator):
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system in [W]
source_rate : float
Power in [W] or source rate in [neutron/sec]
_i : int, optional
Current depletion step index. Not used
@ -947,7 +478,7 @@ class SICELIIntegrator(SIIntegrator):
# Begin iteration
for j in range(self.n_steps + 1):
inter_res = self.operator(inter_conc, power)
inter_res = self.operator(inter_conc, source_rate)
if j <= 1:
res_bar = copy.deepcopy(inter_res)
@ -967,6 +498,7 @@ class SICELIIntegrator(SIIntegrator):
return proc_time, [eos_conc, inter_conc], [res_bar]
@add_params
class SILEQIIntegrator(SIIntegrator):
r"""Deplete using the SI-LE/QI CFQ4 algorithm.
@ -976,84 +508,10 @@ class SILEQIIntegrator(SIIntegrator):
Detailed algorithm can be found in Section 3.2 in `Colin Josey's thesis
<http://hdl.handle.net/1721.1/113721>`_.
Parameters
----------
operator : openmc.deplete.TransportOperator
The operator object to simulate on.
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
iterable of float. Alternatively, units can be specified for each step
by passing an iterable of (value, unit) tuples.
power : float or iterable of float, optional
Power of the reactor in [W]. A single value indicates that
the power is constant over all timesteps. An iterable
indicates potentially different power levels for each timestep.
For a 2D problem, the power can be given in [W/cm] as long
as the "volume" assigned to a depletion material is actually
an area in [cm^2]. Either ``power`` or ``power_density`` must be
specified.
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.
timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
Units for values specified in the `timesteps` argument. 's' means
seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
that the values are given in burnup (MW-d of energy deposited per
kilogram of initial heavy metal).
.. versionadded:: 0.12
n_steps : int, optional
Number of stochastic iterations per depletion interval.
Must be greater than zero. Default : 10
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
----------
operator : openmc.deplete.TransportOperator
Operator to perform transport simulations
chain : openmc.deplete.Chain
Depletion chain
timesteps : iterable of float
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
n_steps : int
Number of stochastic iterations per depletion interval
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
_num_stages = 2
def __call__(self, bos_conc, bos_rates, dt, power, i):
def __call__(self, bos_conc, bos_rates, dt, source_rate, i):
"""Perform the integration across one time step
Parameters
@ -1065,8 +523,8 @@ class SILEQIIntegrator(SIIntegrator):
Reaction rates from operator for all depletable materials
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system in [W]
source_rate : float
Power in [W] or source rate in [neutron/sec]
i : int
Current depletion step index
@ -1086,7 +544,7 @@ class SILEQIIntegrator(SIIntegrator):
self._prev_rates = bos_rates
# Perform CELI for initial steps
return SICELIIntegrator.__call__(
self, bos_conc, bos_rates, dt, power, i)
self, bos_conc, bos_rates, dt, source_rate, i)
prev_res = self.operator.prev_res[-2]
prev_dt = self.timesteps[i] - prev_res.time[0]
self._prev_rates = prev_res.rates[0]
@ -1105,7 +563,7 @@ class SILEQIIntegrator(SIIntegrator):
inter_conc = copy.deepcopy(eos_conc)
for j in range(self.n_steps + 1):
inter_res = self.operator(inter_conc, power)
inter_res = self.operator(inter_conc, source_rate)
if j <= 1:
res_bar = copy.deepcopy(inter_res)

41
openmc/deplete/nuclide.py
View file

@ -154,6 +154,47 @@ class Nuclide:
return None
return self.yield_data.energies
def add_decay_mode(self, type, target, branching_ratio):
"""Add decay mode to the nuclide
Parameters
----------
type : str
Type of the decay mode, e.g., 'beta-'
target : str or None
Nuclide resulting from decay. A value of ``None`` implies the
target does not exist in the currently configured depletion
chain
branching_ratio : float
Branching ratio of the decay mode
"""
self.decay_modes.append(
DecayTuple(type, target, branching_ratio)
)
def add_reaction(self, type, target, Q, branching_ratio):
"""Add transmutation reaction to the nuclide
Parameters
----------
type : str
Type of the reaction, e.g., 'fission'
target : str or None
Nuclide resulting from reaction. A value of ``None``
implies either no single target, e.g. from fission,
or that the target nuclide is not considered
in the current depletion chain
Q : float
Q value of the reaction in [eV]
branching_ratio : float
Branching ratio of the reaction
"""
self.reactions.append(
ReactionTuple(type, target, Q, branching_ratio)
)
@classmethod
def from_xml(cls, element, root=None, fission_q=None):
"""Read nuclide from an XML element.

98
openmc/deplete/operator.py
View file

@ -7,7 +7,6 @@ densities is all done in-memory instead of through the filesystem.
"""
import sys
import copy
from collections import OrderedDict
import os
@ -18,6 +17,7 @@ import numpy as np
from uncertainties import ufloat
import openmc
from openmc.checkvalue import check_value
import openmc.lib
from . import comm
from .abc import TransportOperator, OperatorResult
@ -26,7 +26,8 @@ from .reaction_rates import ReactionRates
from .results_list import ResultsList
from .helpers import (
DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper,
FissionYieldCutoffHelper, AveragedFissionYieldHelper, EnergyScoreHelper)
FissionYieldCutoffHelper, AveragedFissionYieldHelper, EnergyScoreHelper,
SourceRateHelper)
__all__ = ["Operator", "OperatorResult"]
@ -81,15 +82,17 @@ class Operator(TransportOperator):
Whether to differentiate burnable materials with multiple instances.
Volumes are divided equally from the original material volume.
Default: False.
energy_mode : {"energy-deposition", "fission-q"}
Indicator for computing system energy. ``"energy-deposition"`` will
compute with a single energy deposition tally, taking fission energy
release data and heating into consideration. ``"fission-q"`` will
use the fission Q values from the depletion chain
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 ``"energy_mode" == "fission-q"``
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.
@ -164,23 +167,18 @@ class Operator(TransportOperator):
}
def __init__(self, geometry, settings, chain_file=None, prev_results=None,
diff_burnable_mats=False, energy_mode="fission-q",
diff_burnable_mats=False, normalization_mode="fission-q",
fission_q=None, dilute_initial=1.0e3,
fission_yield_mode="constant", fission_yield_opts=None,
reduce_chain=False, reduce_chain_level=None):
if fission_yield_mode not in self._fission_helpers:
raise KeyError(
"fission_yield_mode must be one of {}, not {}".format(
", ".join(self._fission_helpers), fission_yield_mode))
if energy_mode == "energy-deposition":
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 not used if energy deposition "
"is used")
warn("Fission Q dictionary will not be used")
fission_q = None
elif energy_mode != "fission-q":
raise ValueError(
"energy_mode {} not supported. Must be energy-deposition "
"or fission-q".format(energy_mode))
super().__init__(chain_file, fission_q, dilute_initial, prev_results)
self.round_number = False
self.settings = settings
@ -243,11 +241,14 @@ class Operator(TransportOperator):
# Get classes to assist working with tallies
self._rate_helper = DirectReactionRateHelper(
self.reaction_rates.n_nuc, self.reaction_rates.n_react)
if energy_mode == "fission-q":
self._energy_helper = ChainFissionHelper()
else:
if normalization_mode == "fission-q":
self._normalization_helper = ChainFissionHelper()
elif normalization_mode == "energy-deposition":
score = "heating" if settings.photon_transport else "heating-local"
self._energy_helper = EnergyScoreHelper(score)
self._normalization_helper = EnergyScoreHelper(score)
else:
self._normalization_helper = SourceRateHelper()
# Select and create fission yield helper
fission_helper = self._fission_helpers[fission_yield_mode]
@ -256,7 +257,7 @@ class Operator(TransportOperator):
self._yield_helper = fission_helper.from_operator(
self, **fission_yield_opts)
def __call__(self, vec, power):
def __call__(self, vec, source_rate):
"""Runs a simulation.
Simulation will abort under the following circumstances:
@ -267,8 +268,8 @@ class Operator(TransportOperator):
----------
vec : list of numpy.ndarray
Total atoms to be used in function.
power : float
Power of the reactor in [W]
source_rate : float
Power in [W] or source rate in [neutron/sec]
Returns
-------
@ -281,7 +282,7 @@ class Operator(TransportOperator):
# If the source rate is zero, return zero reaction rates without running
# a transport solve
if power == 0.0:
if source_rate == 0.0:
rates = self.reaction_rates.copy()
rates.fill(0.0)
return OperatorResult(ufloat(0.0, 0.0), rates)
@ -296,7 +297,7 @@ class Operator(TransportOperator):
self._update_materials()
nuclides = self._get_tally_nuclides()
self._rate_helper.nuclides = nuclides
self._energy_helper.nuclides = nuclides
self._normalization_helper.nuclides = nuclides
self._yield_helper.update_tally_nuclides(nuclides)
# Run OpenMC
@ -304,7 +305,7 @@ class Operator(TransportOperator):
openmc.lib.reset_timers()
# Extract results
op_result = self._unpack_tallies_and_normalize(power)
op_result = self._unpack_tallies_and_normalize(source_rate)
return copy.deepcopy(op_result)
@ -502,8 +503,8 @@ class Operator(TransportOperator):
materials = [openmc.lib.materials[int(i)]
for i in self.burnable_mats]
self._rate_helper.generate_tallies(materials, self.chain.reactions)
self._energy_helper.prepare(
self.chain.nuclides, self.reaction_rates.index_nuc, materials)
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(
@ -616,18 +617,17 @@ class Operator(TransportOperator):
nuc_list = comm.bcast(nuc_list)
return [nuc for nuc in nuc_list if nuc in self.chain]
def _unpack_tallies_and_normalize(self, power):
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 the user-specified power, summing the product of the
fission reaction rate times the fission Q value for each material.
normalized by a helper class depending on the method being used.
Parameters
----------
power : float
Power of the reactor in [W]
source_rate : float
Power in [W] or source rate in [neutron/sec]
Returns
-------
@ -648,11 +648,9 @@ class Operator(TransportOperator):
nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides]
react_ind = [rates.index_rx[react] for react in self.chain.reactions]
# Compute fission power
# Keep track of energy produced from all reactions in eV per source
# particle
self._energy_helper.reset()
self._normalization_helper.reset()
self._yield_helper.unpack()
# Store fission yield dictionaries
@ -662,7 +660,7 @@ class Operator(TransportOperator):
# numbers, zeroed out in material iteration
number = np.empty(rates.n_nuc)
fission_ind = rates.index_rx["fission"]
fission_ind = rates.index_rx.get("fission")
# Extract results
for i, mat in enumerate(self.local_mats):
@ -683,26 +681,14 @@ class Operator(TransportOperator):
fission_yields.append(self._yield_helper.weighted_yields(i))
# Accumulate energy from fission
self._energy_helper.update(tally_rates[:, fission_ind], mat_index)
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)
# Reduce energy produced from all processes
# J / s / source neutron
energy = comm.allreduce(self._energy_helper.energy)
# Guard against divide by zero
if energy == 0:
if comm.rank == 0:
sys.stderr.flush()
print(" No energy reported from OpenMC tallies. Do your HDF5 "
"files have heating data?\n", file=sys.stderr, flush=True)
comm.barrier()
comm.Abort(1)
# Scale reaction rates to obtain units of reactions/sec
rates *= power / energy
rates *= self._normalization_helper.factor(source_rate)
# Store new fission yields on the chain
self.chain.fission_yields = fission_yields

38
openmc/deplete/results.py
View file

@ -12,7 +12,7 @@ import numpy as np
from . import comm, MPI
from .reaction_rates import ReactionRates
VERSION_RESULTS = (1, 0)
VERSION_RESULTS = (1, 1)
__all__ = ["Results"]
@ -27,8 +27,8 @@ class Results:
Eigenvalue and uncertainty for each substep.
time : list of float
Time at beginning, end of step, in seconds.
power : float
Power during time step, in Watts
source_rate : float
Source rate during timestep in [W] or [neutron/sec]
n_mat : int
Number of mats.
n_nuc : int
@ -57,7 +57,7 @@ class Results:
def __init__(self):
self.k = None
self.time = None
self.power = None
self.source_rate = None
self.rates = None
self.volume = None
self.proc_time = None
@ -177,7 +177,7 @@ class Results:
new.mat_to_ind = {mat: idx for (idx, mat) in enumerate(local_materials)}
# Direct transfer
direct_attrs = ("time", "k", "power", "nuc_to_ind",
direct_attrs = ("time", "k", "source_rate", "nuc_to_ind",
"mat_to_hdf5_ind", "proc_time")
for attr in direct_attrs:
setattr(new, attr, getattr(self, attr))
@ -287,7 +287,7 @@ class Results:
handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64')
handle.create_dataset("power", (1, n_stages), maxshape=(None, n_stages),
handle.create_dataset("source_rate", (1, n_stages), maxshape=(None, n_stages),
dtype='float64')
handle.create_dataset(
@ -320,7 +320,7 @@ class Results:
rxn_dset = handle["/reaction rates"]
eigenvalues_dset = handle["/eigenvalues"]
time_dset = handle["/time"]
power_dset = handle["/power"]
source_rate_dset = handle["/source_rate"]
proc_time_dset = handle["/depletion time"]
# Get number of results stored
@ -346,9 +346,9 @@ class Results:
time_shape[0] = new_shape
time_dset.resize(time_shape)
power_shape = list(power_dset.shape)
power_shape[0] = new_shape
power_dset.resize(power_shape)
source_rate_shape = list(source_rate_dset.shape)
source_rate_shape[0] = new_shape
source_rate_dset.resize(source_rate_shape)
proc_shape = list(proc_time_dset.shape)
proc_shape[0] = new_shape
@ -371,7 +371,7 @@ class Results:
eigenvalues_dset[index, i] = self.k[i]
if comm.rank == 0:
time_dset[index] = self.time
power_dset[index] = self.power
source_rate_dset[index] = self.source_rate
if self.proc_time is not None:
proc_time_dset[index] = (
self.proc_time / (comm.size * self.n_hdf5_mats)
@ -394,12 +394,16 @@ class Results:
number_dset = handle["/number"]
eigenvalues_dset = handle["/eigenvalues"]
time_dset = handle["/time"]
power_dset = handle["/power"]
if "source_rate" in handle:
source_rate_dset = handle["/source_rate"]
else:
# Older versions used "power" instead of "source_rate"
source_rate_dset = handle["/power"]
results.data = number_dset[step, :, :, :]
results.k = eigenvalues_dset[step, :]
results.time = time_dset[step, :]
results.power = power_dset[step, :]
results.source_rate = source_rate_dset[step, :]
if "depletion time" in handle:
proc_time_dset = handle["/depletion time"]
@ -444,7 +448,7 @@ class Results:
return results
@staticmethod
def save(op, x, op_results, t, power, step_ind, proc_time=None):
def save(op, x, op_results, t, source_rate, step_ind, proc_time=None):
"""Creates and writes depletion results to disk
Parameters
@ -457,8 +461,8 @@ class Results:
Results of applying transport operator
t : list of float
Time indices.
power : float
Power during time step
source_rate : float
Source rate during time step in [W] or [neutron/sec]
step_ind : int
Step index.
proc_time : float or None
@ -485,7 +489,7 @@ class Results:
results.k = [(r.k.nominal_value, r.k.std_dev) for r in op_results]
results.rates = [r.rates for r in op_results]
results.time = t
results.power = power
results.source_rate = source_rate
results.proc_time = proc_time
if results.proc_time is not None:
results.proc_time = comm.reduce(proc_time, op=MPI.SUM)

BIN
tests/regression_tests/deplete/test_reference.h5
View file

Binary file not shown.

149
tests/unit_tests/test_deplete_activation.py Normal file
View file

@ -0,0 +1,149 @@
from math import pi, log
from random import uniform, normalvariate
import openmc.deplete
import openmc
import pytest
@pytest.fixture
def model():
"""Sphere of single nuclide"""
model = openmc.model.Model()
w = openmc.Material(name='tungsten')
w.add_nuclide('W186', 1.0)
w.set_density('g/cm3', 19.3)
w.depletable = True
r = uniform(1.0, 10.0)
w.volume = 4/3 * pi * r**3
surf = openmc.Sphere(r=r, boundary_type='vacuum')
cell = openmc.Cell(fill=w, region=-surf)
model.geometry = openmc.Geometry([cell])
model.settings.batches = 10
model.settings.particles = 1000
model.settings.source = openmc.Source(
space=openmc.stats.Point(),
energy=openmc.stats.Discrete([1.0e6], [1.0])
)
model.settings.run_mode = 'fixed source'
rx_tally = openmc.Tally()
rx_tally.scores = ['(n,gamma)']
model.tallies.append(rx_tally)
return model
def test_activation(run_in_tmpdir, model):
# Determine (n.gamma) reaction rate using initial run
sp = model.run()
with openmc.StatePoint(sp) as sp:
tally = sp.tallies[1]
capture_rate = tally.mean.flat[0]
# Create one-nuclide depletion chain
chain = openmc.deplete.Chain()
w186 = openmc.deplete.Nuclide('W186')
w186.add_reaction('(n,gamma)', None, 0.0, 1.0)
chain.add_nuclide(w186)
chain.export_to_xml('test_chain.xml')
# Create transport operator
op = openmc.deplete.Operator(
model.geometry, model.settings, 'test_chain.xml',
normalization_mode="source-rate"
)
# To determine the source rate necessary to reduce W186 density in half, we
# start with the single-nuclide transmutation equation:
#
# dn/dt = -f * sigma * phi * n
# n(t) = n0 * exp(-f * sigma * phi * t)
#
# where f is the source rate. The capture rate, r, is sigma * phi * n0,
# meaning that:
#
# n(t) = n0 * exp(-f * r * t / n0)
#
# To reduce the density by half, we would need:
#
# n(t)/n0 = exp(-f * r * t / n0) = 1/2
# f = n0 / (r * t) ln(2)
#
# So we need to know the initial number of atoms (n0), the capture rate (r),
# and choose an irradiation time (t)
w = model.geometry.get_materials_by_name('tungsten')[0]
atom_densities = w.get_nuclide_atom_densities()
atom_per_cc = 1e24 * atom_densities['W186'][1] # Density in atom/cm^3
n0 = atom_per_cc * w.volume # Absolute number of atoms
# Pick a random irradiation time and then determine necessary source rate to
# reduce material by half
t = uniform(1.0, 5.0) * 86400
source_rates = [n0/(capture_rate*t) * log(2.0)]
# Now activate the material
integrator = openmc.deplete.PredictorIntegrator(
op, [t], source_rates=source_rates
)
integrator.integrate()
# Get resulting number of atoms
results = openmc.deplete.ResultsList.from_hdf5('depletion_results.h5')
_, atoms = results.get_atoms(str(w.id), "W186")
assert atoms[0] == pytest.approx(n0)
assert atoms[1] / atoms[0] == pytest.approx(0.5, rel=1e-3)
def test_decay(run_in_tmpdir):
"""Test decay-only timesteps where no transport solve is performed"""
# Create a model with a single nuclide, Sr89
mat = openmc.Material()
mat.add_nuclide('Sr89', 1.0)
mat.set_density('g/cm3', 1.0)
mat.depletable = True
r = 5.0
mat.volume = 4/3 * pi * r**3
surf = openmc.Sphere(r=r, boundary_type='vacuum')
cell = openmc.Cell(fill=mat, region=-surf)
geometry = openmc.Geometry([cell])
settings = openmc.Settings()
settings.batches = 10
settings.particles = 1000
settings.run_mode = 'fixed source'
# Create depletion chain with only Sr89 and sample its half-life. Note that
# currently at least one reaction has to exist in the depletion chain
chain = openmc.deplete.Chain()
sr89 = openmc.deplete.Nuclide('Sr89')
sr89.half_life = normalvariate(4365792.0, 6048.0)
sr89.add_decay_mode('beta-', None, 1.0)
sr89.add_reaction('(n,gamma)', None, 0.0, 1.0)
chain.add_nuclide(sr89)
chain.export_to_xml('test_chain.xml')
# Create transport operator
op = openmc.deplete.Operator(
geometry, settings, 'test_chain.xml', normalization_mode="source-rate"
)
# Deplete with two decay steps
integrator = openmc.deplete.PredictorIntegrator(
op, [sr89.half_life, 2*sr89.half_life], source_rates=[0.0, 0.0]
)
integrator.integrate()
# Get resulting number of atoms
results = openmc.deplete.ResultsList.from_hdf5('depletion_results.h5')
_, atoms = results.get_atoms(str(mat.id), "Sr89")
# Ensure density goes down by a factor of 2 after each half-life
assert atoms[1] / atoms[0] == pytest.approx(0.5)
assert atoms[2] / atoms[1] == pytest.approx(0.25)

15
tests/unit_tests/test_deplete_fission_yields.py
View file

@ -54,8 +54,7 @@ def materials(tmpdir_factory):
def proxy_tally_data(tally, fill=None):
"""Construct an empty matrix built from a C tally
The shape of tally.results will be
``(n_bins, n_nuc * n_scores, 3)``
The shape of tally.mean will be ``(n_bins, n_nuc * n_scores)``
"""
n_nucs = max(len(tally.nuclides), 1)
n_scores = max(len(tally.scores), 1)
@ -67,7 +66,7 @@ def proxy_tally_data(tally, fill=None):
if isinstance(tfilter, lib.EnergyFilter):
this_bins -= 1
n_bins *= max(this_bins, 1)
data = np.empty((n_bins, n_nucs * n_scores, 3))
data = np.empty((n_bins, n_nucs * n_scores))
if fill is not None:
data.fill(fill)
return data
@ -203,9 +202,9 @@ def test_cutoff_helper(materials, nuclide_bundle, therm_frac):
tally_data = proxy_tally_data(fission_tally)
helper._fission_rate_tally = Mock()
helper_flux = 1e6
tally_data[0, :, 1] = therm_frac * helper_flux
tally_data[1, :, 1] = (1 - therm_frac) * helper_flux
helper._fission_rate_tally.results = tally_data
tally_data[0] = therm_frac * helper_flux
tally_data[1] = (1 - therm_frac) * helper_flux
helper._fission_rate_tally.mean = tally_data
helper.unpack()
# expected results of shape (n_mats, 2, n_tnucs)
@ -261,8 +260,8 @@ def test_averaged_helper(materials, nuclide_bundle, avg_energy):
helper._fission_rate_tally = Mock()
helper._weighted_tally = Mock()
helper._fission_rate_tally.results = fission_results
helper._weighted_tally.results = weighted_results
helper._fission_rate_tally.mean = fission_results
helper._weighted_tally.mean = weighted_results
helper.unpack()
expected_results = np.ones((1, len(tallied_nucs))) * avg_energy

2
tests/unit_tests/test_deplete_integrator.py
View file

@ -130,7 +130,7 @@ def test_bad_integrator_inputs():
timesteps = [1]
# No power nor power density given
with pytest.raises(ValueError, match="Either power or power density"):
with pytest.raises(ValueError, match="Either power"):
PredictorIntegrator(op, timesteps)
# Length of power != length time