From bf8405f6cc6e5bc073cd565ffb1232ae9d68ec9e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 13 Jul 2020 16:15:14 -0500 Subject: [PATCH 01/21] Change EnergyHelper -> NormalizationHelper --- docs/source/pythonapi/deplete.rst | 2 +- openmc/deplete/abc.py | 16 +++++------ openmc/deplete/helpers.py | 16 +++++------ openmc/deplete/operator.py | 44 +++++++++++++++---------------- 4 files changed, 38 insertions(+), 40 deletions(-) diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index f74c4524a3..e9cd6d3052 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -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 diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 4924e4f6fb..ac24e565ff 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -29,7 +29,7 @@ from .pool import deplete __all__ = [ "OperatorResult", "TransportOperator", "ReactionRateHelper", - "EnergyHelper", "FissionYieldHelper", "TalliedFissionYieldHelper", + "NormalizationHelper", "FissionYieldHelper", "TalliedFissionYieldHelper", "Integrator", "SIIntegrator", "DepSystemSolver"] @@ -302,14 +302,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 diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 8c17edac8a..2c1afb4da8 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -14,7 +14,7 @@ from openmc.checkvalue import check_type, check_greater_than from openmc.lib import ( Tally, MaterialFilter, EnergyFilter, EnergyFunctionFilter) from .abc import ( - ReactionRateHelper, EnergyHelper, FissionYieldHelper, + ReactionRateHelper, NormalizationHelper, FissionYieldHelper, TalliedFissionYieldHelper) __all__ = ( @@ -91,13 +91,13 @@ 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 ChainFissionHelper(NormalizationHelper): + """Computes normalization using fission Q values from depletion chain Attributes ---------- @@ -159,7 +159,7 @@ class ChainFissionHelper(EnergyHelper): self._energy += dot(fission_rates, self._fission_q_vector) -class EnergyScoreHelper(EnergyHelper): +class EnergyScoreHelper(NormalizationHelper): """Class responsible for obtaining system energy via a tally score Parameters @@ -172,7 +172,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 diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index a11f6d95dc..b9d870d9db 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -18,6 +18,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 @@ -81,7 +82,7 @@ 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"} + normalization_mode : {"energy-deposition", "fission-q", "source-rate"} 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 @@ -89,7 +90,7 @@ class Operator(TransportOperator): 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 +165,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 +239,13 @@ 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 = ... # Select and create fission yield helper fission_helper = self._fission_helpers[fission_yield_mode] @@ -286,7 +284,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 @@ -493,7 +491,7 @@ 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._normalization_helper.prepare( self.chain.nuclides, self.reaction_rates.index_nuc, materials) # Tell fission yield helper what materials this process is # responsible for @@ -643,7 +641,7 @@ class Operator(TransportOperator): # 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 @@ -674,14 +672,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) + self._normalization_helper.update(tally_rates[:, fission_ind], mat_index) # 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) + energy = comm.allreduce(self._normalization_helper.energy) # Guard against divide by zero if energy == 0: From cdf58828af90ddb7f0f1ce41a0e4ede3113fa259 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 13 Jul 2020 16:28:22 -0500 Subject: [PATCH 02/21] Remove unused arguments on NormalizationHelper classes --- openmc/deplete/abc.py | 11 ++--------- openmc/deplete/helpers.py | 9 ++------- openmc/deplete/operator.py | 4 ++-- 3 files changed, 6 insertions(+), 18 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index ac24e565ff..f07174a95b 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -335,7 +335,7 @@ class NormalizationHelper(ABC): self._energy = 0.0 @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 @@ -348,13 +348,9 @@ class NormalizationHelper(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): + def update(self, fission_rates): """Update the energy produced Parameters @@ -363,9 +359,6 @@ 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 - Index for the specific material in the list of all burnable - materials. """ @property diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 2c1afb4da8..204a214256 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -113,7 +113,7 @@ class ChainFissionHelper(NormalizationHelper): 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 +125,6 @@ class ChainFissionHelper(NormalizationHelper): 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 +141,7 @@ class ChainFissionHelper(NormalizationHelper): 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,9 +150,6 @@ class ChainFissionHelper(NormalizationHelper): 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) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index b9d870d9db..05c9cf5e4e 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -492,7 +492,7 @@ class Operator(TransportOperator): 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, materials) + self.chain.nuclides, self.reaction_rates.index_nuc) # Tell fission yield helper what materials this process is # responsible for self._yield_helper.generate_tallies( @@ -672,7 +672,7 @@ class Operator(TransportOperator): fission_yields.append(self._yield_helper.weighted_yields(i)) # Accumulate energy from fission - self._normalization_helper.update(tally_rates[:, fission_ind], mat_index) + self._normalization_helper.update(tally_rates[:, fission_ind]) # Divide by total number and store rates[i] = self._rate_helper.divide_by_adens(number) From 04501e3213cc81a5fb4f68dbc0f7d78af485cad3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 13 Jul 2020 22:53:59 -0500 Subject: [PATCH 03/21] Move calculation of normalization factor into NormalizationHelper --- openmc/deplete/abc.py | 38 +++++++++++++++++++++++++------------- openmc/deplete/operator.py | 16 +--------------- 2 files changed, 26 insertions(+), 28 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index f07174a95b..2fb85a6771 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -4,16 +4,17 @@ 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 @@ -21,6 +22,7 @@ from uncertainties import ufloat from openmc.data import DataLibrary, JOULE_PER_EV 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 @@ -317,18 +319,11 @@ class NormalizationHelper(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""" @@ -371,6 +366,23 @@ class NormalizationHelper(ABC): check_type("nuclides", nuclides, list, str) self._nuclides = nuclides + def factor(self, power): + # Reduce energy produced from all processes + # J / s / 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 + return power / energy + class FissionYieldHelper(ABC): """Abstract class for processing energy dependent fission yields diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 05c9cf5e4e..281360634f 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -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 @@ -677,21 +676,8 @@ class Operator(TransportOperator): # 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._normalization_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(power) # Store new fission yields on the chain self.chain.fission_yields = fission_yields From e2ae3c74012759fc0870aa4fffbd5340a04d5c29 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 15 Jul 2020 08:09:37 -0500 Subject: [PATCH 04/21] Change power to source_rate throughout deplete module --- openmc/deplete/abc.py | 97 ++++++++++-------- openmc/deplete/operator.py | 8 +- openmc/deplete/results.py | 32 +++--- .../deplete/test_reference.h5 | Bin 165504 -> 165832 bytes tests/unit_tests/test_deplete_integrator.py | 2 +- 5 files changed, 74 insertions(+), 65 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 2fb85a6771..d5ad97d673 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -138,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 ------- @@ -366,9 +366,9 @@ class NormalizationHelper(ABC): check_type("nuclides", nuclides, list, str) self._nuclides = nuclides - def factor(self, power): + def factor(self, source_rate): # Reduce energy produced from all processes - # J / s / source neutron + # J / source neutron energy = comm.allreduce(self._energy) * JOULE_PER_EV # Guard against divide by zero @@ -380,8 +380,9 @@ class NormalizationHelper(ABC): comm.barrier() comm.Abort(1) - # Return normalization factor for scaling reaction rates - return power / energy + # 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 FissionYieldHelper(ABC): @@ -619,12 +620,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 : iterable of float + Source rate in [neutrons/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 @@ -650,8 +655,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 [neutrons/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 @@ -672,7 +678,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] @@ -686,22 +692,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): @@ -712,13 +721,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) @@ -731,13 +740,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 @@ -830,22 +839,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 [neutrons/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 @@ -853,8 +862,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): @@ -868,12 +877,12 @@ 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): 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) - proc_time, conc_list, res_list = self(conc, res.rates, dt, p, i) + conc, res = self._get_bos_data_from_restart(i, source_rate, conc) + proc_time, conc_list, res_list = self(conc, res.rates, dt, source_rate, i) # Insert BOS concentration, transport results conc_list.insert(0, conc) @@ -883,14 +892,14 @@ 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 # Final simulation - res_list = [self.operator(conc, p)] + res_list = [self.operator(conc, source_rate)] 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) @@ -976,8 +985,8 @@ class SIIntegrator(Integrator): 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, + timestep_units=timestep_units, solver=solver) self.n_steps = n_steps def _get_bos_data_from_operator(self, step_index, step_power, bos_conc): diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 281360634f..fecf9f6944 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -604,7 +604,7 @@ 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 @@ -614,8 +614,8 @@ class Operator(TransportOperator): Parameters ---------- - power : float - Power of the reactor in [W] + source_rate : float + Power in [W] or source rate in [neutrons/sec] Returns ------- @@ -677,7 +677,7 @@ class Operator(TransportOperator): rates[i] = self._rate_helper.divide_by_adens(number) # Scale reaction rates to obtain units of reactions/sec - rates *= self._normalization_helper.factor(power) + rates *= self._normalization_helper.factor(source_rate) # Store new fission yields on the chain self.chain.fission_yields = fission_yields diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index 11191211cf..c8c4da1d9d 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -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 [neutrons/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,12 @@ class Results: number_dset = handle["/number"] eigenvalues_dset = handle["/eigenvalues"] time_dset = handle["/time"] - power_dset = handle["/power"] + source_rate_dset = handle["/source_rate"] 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 +444,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 +457,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 [neutrons/sec] step_ind : int Step index. proc_time : float or None @@ -485,7 +485,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) diff --git a/tests/regression_tests/deplete/test_reference.h5 b/tests/regression_tests/deplete/test_reference.h5 index 3ee9d4447a2c9243d4b695e43a9345016ef66619..922f60e5d8b81f5823f232e624afc75174ee3586 100644 GIT binary patch delta 133 zcmZqZ>vMEgXo6h@|r{fgWc3=q)J#(}}ZI9T4 b{rp`R8CWKlvj;L3Ob%pMpL~K{185QeML{by delta 104 zcmX@n&(+Y&H9>=^p>3j8%Jg;Cj8dE{*ciYdVY4B}cgD#EY#fs}FmW&zY~IQAm6?%a z@_M#~j2Y8!wlS$s?`UHRnAXmu(fomZ`v-Q$9kV$ZA$l1cCLT1Oe!!JcX#20ljBSqq D#V8_H diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py index 0560cb9f73..3195c59768 100644 --- a/tests/unit_tests/test_deplete_integrator.py +++ b/tests/unit_tests/test_deplete_integrator.py @@ -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 From 8d489cc23bd9efeb7b6917f6756816cff5241539 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 15 Jul 2020 21:28:35 -0500 Subject: [PATCH 05/21] Simplify Integrator docstrings, few more power -> source_rate --- openmc/deplete/abc.py | 20 +- openmc/deplete/integrators.py | 560 +--------------------------------- openmc/deplete/operator.py | 15 +- openmc/deplete/results.py | 4 +- 4 files changed, 33 insertions(+), 566 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index d5ad97d673..5ccb90dec7 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -32,7 +32,7 @@ from .pool import deplete __all__ = [ "OperatorResult", "TransportOperator", "ReactionRateHelper", "NormalizationHelper", "FissionYieldHelper", "TalliedFissionYieldHelper", - "Integrator", "SIIntegrator", "DepSystemSolver"] + "Integrator", "SIIntegrator", "DepSystemSolver", "add_params"] _SECONDS_PER_MINUTE = 60 @@ -602,9 +602,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 @@ -802,7 +810,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 @@ -813,8 +821,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 [neutrons/sec] i : int Current depletion step index @@ -903,12 +911,15 @@ class Integrator(ABC): 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 @@ -980,6 +991,7 @@ class SIIntegrator(Integrator): .. versionadded:: 0.12 """ + def __init__(self, operator, timesteps, power=None, power_density=None, timestep_units='s', n_steps=10, solver="cram48"): check_type("n_steps", n_steps, Integral) diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py index 1fe6369003..10aeb6f88c 100644 --- a/openmc/deplete/integrators.py +++ b/openmc/deplete/integrators.py @@ -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,74 +26,6 @@ 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 @@ -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,74 +78,6 @@ 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 @@ -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,74 +138,6 @@ 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 @@ -397,6 +196,7 @@ class CF4Integrator(Integrator): [res1, res2, res3]) +@add_params class CELIIntegrator(Integrator): r"""Deplete using the CE/LI CFQ4 algorithm. @@ -415,74 +215,6 @@ 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 @@ -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,69 +279,6 @@ 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 @@ -663,6 +333,7 @@ class EPCRK4Integrator(Integrator): [res1, res2, res3]) +@add_params class LEQIIntegrator(Integrator): r"""Deplete using the LE/QI CFQ4 algorithm. @@ -691,74 +362,6 @@ 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 @@ -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,79 +443,6 @@ class SICELIIntegrator(SIIntegrator): Detailed algorithm can be found in section 3.2 in `Colin Josey's thesis `_. - - 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 @@ -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,80 +508,6 @@ class SILEQIIntegrator(SIIntegrator): Detailed algorithm can be found in Section 3.2 in `Colin Josey's thesis `_. - - 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 diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index fecf9f6944..d6acd20394 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -253,7 +253,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: @@ -264,8 +264,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 + Source rate in [W] in [neutron/sec] Returns ------- @@ -292,7 +292,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) @@ -609,13 +609,12 @@ class Operator(TransportOperator): 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 ---------- source_rate : float - Power in [W] or source rate in [neutrons/sec] + Power in [W] or source rate in [neutron/sec] Returns ------- @@ -636,8 +635,6 @@ 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._normalization_helper.reset() diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index c8c4da1d9d..90e7573f21 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -28,7 +28,7 @@ class Results: time : list of float Time at beginning, end of step, in seconds. source_rate : float - Source rate during timestep in [W] or [neutrons/sec] + Source rate during timestep in [W] or [neutron/sec] n_mat : int Number of mats. n_nuc : int @@ -458,7 +458,7 @@ class Results: t : list of float Time indices. source_rate : float - Source rate during time step in [W] or [neutrons/sec] + Source rate during time step in [W] or [neutron/sec] step_ind : int Step index. proc_time : float or None From 14fdcf9d4da56ba02ffeafa11c3b17aa0bea792e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 15 Jul 2020 21:53:06 -0500 Subject: [PATCH 06/21] More renames of power to source_rate --- openmc/deplete/abc.py | 4 +- openmc/deplete/integrators.py | 76 +++++++++++++++++------------------ 2 files changed, 40 insertions(+), 40 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 5ccb90dec7..a30852e0af 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -822,7 +822,7 @@ class Integrator(ABC): dt : float Time in [s] for the entire depletion interval source_rate : float - Power in [W] or source rate in [neutrons/sec] + Power in [W] or source rate in [neutron/sec] i : int Current depletion step index @@ -847,7 +847,7 @@ class Integrator(ABC): """ def __iter__(self): - """Return pair of time step in [s] and source rate in [W] or [neutrons/sec]""" + """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): diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py index 10aeb6f88c..fbf43a808e 100644 --- a/openmc/deplete/integrators.py +++ b/openmc/deplete/integrators.py @@ -29,7 +29,7 @@ class PredictorIntegrator(Integrator): """ _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 @@ -40,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 @@ -81,7 +81,7 @@ class CECMIntegrator(Integrator): """ _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 @@ -92,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 @@ -109,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 @@ -141,7 +141,7 @@ class CF4Integrator(Integrator): """ _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 @@ -152,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 @@ -171,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)) @@ -218,7 +218,7 @@ class CELIIntegrator(Integrator): """ _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 @@ -229,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 @@ -247,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)) @@ -282,7 +282,7 @@ class EPCRK4Integrator(Integrator): """ _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 @@ -293,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. @@ -313,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)) @@ -365,7 +365,7 @@ class LEQIIntegrator(Integrator): """ _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 @@ -376,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 @@ -396,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] @@ -404,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))) @@ -414,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, @@ -446,7 +446,7 @@ class SICELIIntegrator(SIIntegrator): """ _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 @@ -457,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 @@ -478,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) @@ -511,7 +511,7 @@ class SILEQIIntegrator(SIIntegrator): """ _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 @@ -523,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 @@ -544,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] @@ -563,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) From baecead5da7c2dc3507107b5d624fc05e4af7273 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 17 Jul 2020 14:19:49 -0500 Subject: [PATCH 07/21] Start creating SourceRateHelper --- openmc/deplete/helpers.py | 12 ++++++++++-- openmc/deplete/operator.py | 6 ++++-- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 204a214256..b52b32ce80 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -18,8 +18,8 @@ from .abc import ( TalliedFissionYieldHelper) __all__ = ( - "DirectReactionRateHelper", "ChainFissionHelper", - "ConstantFissionYieldHelper", "FissionYieldCutoffHelper", + "DirectReactionRateHelper", "ChainFissionHelper", "EnergyScoreHelper" + "SourceRateHelper", "ConstantFissionYieldHelper", "FissionYieldCutoffHelper", "AveragedFissionYieldHelper") # ------------------------------------- @@ -205,6 +205,14 @@ class EnergyScoreHelper(NormalizationHelper): if comm.rank == 0: self._energy = self._tally.results[0, 0, 1] + +class SourceRateHelper(NormalizationHelper): + def prepare(self, *args, **kwargs): + pass + + def factor(self, source_rate): + return source_rate + # ------------------------------------ # Helper for collapsing fission yields # ------------------------------------ diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index d6acd20394..499738e8c5 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -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"] @@ -238,13 +239,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 normalization_mode == "fission-q": self._normalization_helper = ChainFissionHelper() elif normalization_mode == "energy-deposition": score = "heating" if settings.photon_transport else "heating-local" self._normalization_helper = EnergyScoreHelper(score) else: - self._normalization_helper = ... + self._normalization_helper = SourceRateHelper() # Select and create fission yield helper fission_helper = self._fission_helpers[fission_yield_mode] From 06e598f38e440decb37d27cebcb2491cb2e4e263 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 22 Jul 2020 22:50:05 -0500 Subject: [PATCH 08/21] Use tally mean in helpers rather than direct access to results --- openmc/deplete/helpers.py | 10 +++++----- tests/unit_tests/test_deplete_fission_yields.py | 15 +++++++-------- 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index b52b32ce80..684ac80c98 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -83,7 +83,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] @@ -203,7 +203,7 @@ class EnergyScoreHelper(NormalizationHelper): """ super().reset() if comm.rank == 0: - self._energy = self._tally.results[0, 0, 1] + self._energy = self._tally.mean[0, 0] class SourceRateHelper(NormalizationHelper): @@ -450,7 +450,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) @@ -612,9 +612,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] diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index 854c530f95..8839b42b3d 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -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 From 1d187b3a8ad0daaf664458347b71c765e5ed94c8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Jul 2020 13:40:50 -0500 Subject: [PATCH 09/21] Add Chain.add_nuclide method --- openmc/deplete/chain.py | 37 +++++++++++++++++++++---------------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 8293273d47..0e31705cdb 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -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 @@ -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: @@ -439,8 +450,6 @@ class Chain: ReactionTuple('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: @@ -474,6 +483,9 @@ class Chain: nuclide.yield_data = FissionYieldDistribution(yield_data) + # Add nuclide to chain + chain.add_nuclide(nuclide) + # Replace missing FPY data for nuclide in chain.nuclides: if hasattr(nuclide, '_fpy'): @@ -532,14 +544,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 From d83b9568c9adb3dcdad47136f69c364489884ba8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Jul 2020 14:15:10 -0500 Subject: [PATCH 10/21] Add two methods on Nuclide: add_decay_mode and add_reaction --- openmc/deplete/chain.py | 102 +++++++++++++++----------------------- openmc/deplete/nuclide.py | 41 +++++++++++++++ 2 files changed, 81 insertions(+), 62 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 0e31705cdb..88ccd95ef1 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -415,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: @@ -441,47 +441,44 @@ 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 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) @@ -881,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 @@ -893,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): @@ -1018,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] @@ -1030,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 diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 767163b057..edba9748c7 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -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. From 72cae7ffa3d1fe65c4621f241f11ff1d3cf70537 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Jul 2020 22:42:08 -0500 Subject: [PATCH 11/21] Allow a depletion chain with no fission to function --- openmc/deplete/operator.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 499738e8c5..eb41fe9dfd 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -649,7 +649,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): @@ -670,7 +670,8 @@ class Operator(TransportOperator): fission_yields.append(self._yield_helper.weighted_yields(i)) # Accumulate energy from fission - self._normalization_helper.update(tally_rates[:, fission_ind]) + 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) From 7ff79778a28e0654e8728af86dbdce3257e5fc25 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Jul 2020 22:42:29 -0500 Subject: [PATCH 12/21] Add new activation test --- tests/unit_tests/test_deplete_activation.py | 101 ++++++++++++++++++++ 1 file changed, 101 insertions(+) create mode 100644 tests/unit_tests/test_deplete_activation.py diff --git a/tests/unit_tests/test_deplete_activation.py b/tests/unit_tests/test_deplete_activation.py new file mode 100644 index 0000000000..d30557f0d4 --- /dev/null +++ b/tests/unit_tests/test_deplete_activation.py @@ -0,0 +1,101 @@ +from math import pi, log +from random import uniform + +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.get_values().ravel()[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) From 673052f7ed21049f47ffb3863b0f81b241321943 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 27 Jul 2020 11:36:38 -0500 Subject: [PATCH 13/21] Fix (n,na) score --- src/tallies/tally.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index ecdf926139..79e1b1bfb3 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -137,7 +137,7 @@ score_str_to_int(std::string score_str) if (score_str == "(n,2nd)") return N_2ND; if (score_str == "(n,na)") - return N_2NA; + return N_NA; if (score_str == "(n,n3a)") return N_N3A; if (score_str == "(n,2na)") @@ -233,7 +233,6 @@ score_str_to_int(std::string score_str) || score_str.rfind("flux-y", 0) == 0) fatal_error(score_str + " is no longer an available score"); - // Assume the given string is a reaction MT number. Make sure it's a natural // number then return. int MT; From e266bd244b24470ca3ca3940c6a57d14b5b82c43 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 3 Aug 2020 14:25:58 -0500 Subject: [PATCH 14/21] Move energy-related details to intermediate EnergyNormalizationHelper class --- openmc/deplete/abc.py | 37 ++++++++++++++++++------------------- openmc/deplete/helpers.py | 31 +++++++++++++++++++++++++++++-- 2 files changed, 47 insertions(+), 21 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index a30852e0af..7a586a69c9 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -19,7 +19,7 @@ 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 @@ -326,15 +326,15 @@ class NormalizationHelper(ABC): self._nuclides = None 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): """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 ---------- @@ -346,7 +346,8 @@ class NormalizationHelper(ABC): """ def update(self, fission_rates): - """Update the energy produced + """Update the normalization based on fission rates (only used for + energy-based normalization) Parameters ---------- @@ -366,23 +367,21 @@ class NormalizationHelper(ABC): check_type("nuclides", nuclides, list, str) self._nuclides = nuclides + @abstractmethod def factor(self, source_rate): - # Reduce energy produced from all processes - # J / source neutron - energy = comm.allreduce(self._energy) * JOULE_PER_EV + """Return normalization factor - # 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) + Parameters + ---------- + source_rate : float + Power in [W] or source rate in [neutron/sec] - # 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 + Returns + ------- + float + Normalization factor for tallies + + """ class FissionYieldHelper(ABC): diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 684ac80c98..edcf6a5b3f 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -11,6 +11,7 @@ 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 ( @@ -96,7 +97,33 @@ class DirectReactionRateHelper(ReactionRateHelper): # ------------------------------------------ -class ChainFissionHelper(NormalizationHelper): +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 @@ -154,7 +181,7 @@ class ChainFissionHelper(NormalizationHelper): self._energy += dot(fission_rates, self._fission_q_vector) -class EnergyScoreHelper(NormalizationHelper): +class EnergyScoreHelper(EnergyNormalizationHelper): """Class responsible for obtaining system energy via a tally score Parameters From 32d2aada20c1dca674edfb8d89eb5203fffe04bd Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 3 Aug 2020 15:01:09 -0500 Subject: [PATCH 15/21] Update documentation (energy_mode -> normalization_mode) --- docs/source/methods/depletion.rst | 4 ++-- docs/source/usersguide/depletion.rst | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/source/methods/depletion.rst b/docs/source/methods/depletion.rst index ebd238eaa6..dce1f2503e 100644 --- a/docs/source/methods/depletion.rst +++ b/docs/source/methods/depletion.rst @@ -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. diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst index 2f59d72522..3d5517416d 100644 --- a/docs/source/usersguide/depletion.rst +++ b/docs/source/usersguide/depletion.rst @@ -88,7 +88,7 @@ reactions described in :ref:`methods_heating`. These values can be used to norm 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 From e682aa88e4ae5a04779d6388c078e887536e2817 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 6 Aug 2020 08:12:36 -0500 Subject: [PATCH 16/21] Address @drewejohnson comments on #1628 --- openmc/deplete/abc.py | 15 ++++++++++----- tests/unit_tests/test_deplete_activation.py | 2 +- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 7a586a69c9..af32fe3979 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -633,7 +633,7 @@ class Integrator(ABC): 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 : iterable of float + source_rates : float or iterable of float, optional Source rate in [neutrons/sec] for each interval in :attr:`timesteps` .. versionadded:: 0.12.1 @@ -934,12 +934,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 [neutrons/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 @@ -992,11 +996,12 @@ class SIIntegrator(Integrator): """ 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, + operator, timesteps, power, power_density, source_rates, timestep_units=timestep_units, solver=solver) self.n_steps = n_steps diff --git a/tests/unit_tests/test_deplete_activation.py b/tests/unit_tests/test_deplete_activation.py index d30557f0d4..239ae52502 100644 --- a/tests/unit_tests/test_deplete_activation.py +++ b/tests/unit_tests/test_deplete_activation.py @@ -43,7 +43,7 @@ def test_activation(run_in_tmpdir, model): sp = model.run() with openmc.StatePoint(sp) as sp: tally = sp.tallies[1] - capture_rate = tally.get_values().ravel()[0] + capture_rate = tally.mean.flat[0] # Create one-nuclide depletion chain chain = openmc.deplete.Chain() From c3d3b5a71f9c478a0b49d59aa7a678138732dc1d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 6 Aug 2020 08:23:23 -0500 Subject: [PATCH 17/21] Increment version of depletion results file, and support reading older version --- docs/source/io_formats/depletion_results.rst | 9 ++++++--- openmc/deplete/results.py | 8 ++++++-- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/docs/source/io_formats/depletion_results.rst b/docs/source/io_formats/depletion_results.rst index de28b04777..ad6f858f2c 100644 --- a/docs/source/io_formats/depletion_results.rst +++ b/docs/source/io_formats/depletion_results.rst @@ -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 + [neutrons/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. diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index 90e7573f21..ce8331e825 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -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"] @@ -394,7 +394,11 @@ class Results: number_dset = handle["/number"] eigenvalues_dset = handle["/eigenvalues"] time_dset = handle["/time"] - source_rate_dset = handle["/source_rate"] + 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, :] From 14175b50ed2ba5f64ee5889678ef90392884d0c9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 7 Aug 2020 11:10:35 -0500 Subject: [PATCH 18/21] Update depletion documentation to mention fixed-source mode --- docs/source/usersguide/depletion.rst | 99 ++++++++++++++++++---------- 1 file changed, 63 insertions(+), 36 deletions(-) diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst index 3d5517416d..ad1ba779e6 100644 --- a/docs/source/usersguide/depletion.rst +++ b/docs/source/usersguide/depletion.rst @@ -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 @@ -51,6 +51,32 @@ 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' + +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 +87,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 +109,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", 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 +149,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:: From c3968ff9a98e146a218468c892201b15293d7b78 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 10 Aug 2020 07:29:37 -0500 Subject: [PATCH 19/21] Expand description of normalization_mode in docstring --- docs/source/io_formats/depletion_results.rst | 2 +- openmc/deplete/abc.py | 6 +++--- openmc/deplete/operator.py | 12 +++++++----- 3 files changed, 11 insertions(+), 9 deletions(-) diff --git a/docs/source/io_formats/depletion_results.rst b/docs/source/io_formats/depletion_results.rst index ad6f858f2c..172c17ae63 100644 --- a/docs/source/io_formats/depletion_results.rst +++ b/docs/source/io_formats/depletion_results.rst @@ -26,7 +26,7 @@ The current version of the depletion results file format is 1.1. - **time** (*double[][2]*) -- Time in [s] at beginning/end of each step. - **source_rate** (*double[][]*) -- Power in [W] or source rate in - [neutrons/sec]. This array has shape (number of timesteps, number + [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, diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index af32fe3979..9b9fe182ab 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -634,7 +634,7 @@ class Integrator(ABC): initial heavy metal inventory to get total power if ``power`` is not speficied. source_rates : float or iterable of float, optional - Source rate in [neutrons/sec] for each interval in :attr:`timesteps` + Source rate in [neutron/sec] for each interval in :attr:`timesteps` .. versionadded:: 0.12.1 timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'} @@ -663,7 +663,7 @@ class Integrator(ABC): timesteps : iterable of float Size of each depletion interval in [s] source_rates : iterable of float - Source rate in [W] or [neutrons/sec] for each interval in + Source rate in [W] or [neutron/sec] for each interval in :attr:`timesteps` solver : callable Function that will solve the Bateman equations @@ -941,7 +941,7 @@ class SIIntegrator(Integrator): initial heavy metal inventory to get total power if ``power`` is not speficied. source_rates : float or iterable of float, optional - Source rate in [neutrons/sec] for each interval in :attr:`timesteps` + Source rate in [neutron/sec] for each interval in :attr:`timesteps` .. versionadded:: 0.12.1 timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'} diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index eb41fe9dfd..73974525af 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -83,10 +83,12 @@ class Operator(TransportOperator): Volumes are divided equally from the original material volume. Default: False. normalization_mode : {"energy-deposition", "fission-q", "source-rate"} - 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 + 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 @@ -267,7 +269,7 @@ class Operator(TransportOperator): vec : list of numpy.ndarray Total atoms to be used in function. source_rate : float - Source rate in [W] in [neutron/sec] + Power in [W] or source rate in [neutron/sec] Returns ------- From efbd375fcee23c6095fb5e63f578e27d65944b3d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 10 Aug 2020 08:03:23 -0500 Subject: [PATCH 20/21] Add note about depletable attribute in user's guide documentation --- docs/source/usersguide/depletion.rst | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst index ad1ba779e6..d550222cf4 100644 --- a/docs/source/usersguide/depletion.rst +++ b/docs/source/usersguide/depletion.rst @@ -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 @@ -63,6 +65,12 @@ 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:: From 1d1e7fc5c15a104ed222a696b2afe2e41b52cefb Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 13 Aug 2020 09:24:47 -0500 Subject: [PATCH 21/21] Add depletion decay-only timestep test --- tests/unit_tests/test_deplete_activation.py | 50 ++++++++++++++++++++- 1 file changed, 49 insertions(+), 1 deletion(-) diff --git a/tests/unit_tests/test_deplete_activation.py b/tests/unit_tests/test_deplete_activation.py index 239ae52502..52936e16a3 100644 --- a/tests/unit_tests/test_deplete_activation.py +++ b/tests/unit_tests/test_deplete_activation.py @@ -1,5 +1,5 @@ from math import pi, log -from random import uniform +from random import uniform, normalvariate import openmc.deplete import openmc @@ -99,3 +99,51 @@ def test_activation(run_in_tmpdir, model): 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)