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

Merge branch 'develop' into bug-restart-mpi

Browse source
This commit is contained in:
Andrew Johnson 2019-08-07 13:50:42 -05:00
commit becda0c2c2
No known key found for this signature in database
GPG key ID: 253418E91B7F6FEB
34 changed files with 1445 additions and 1402 deletions
Download patch file Download diff file Expand all files Collapse all files

9
docs/source/_templates/myintegrator.rst Normal file
View file

@ -0,0 +1,9 @@
{{ fullname }}
{{ underline }}
.. currentmodule:: {{ module }}
.. autoclass:: {{ objname }}
:members:
:inherited-members:
:special-members: __call__, __len__, __iter__

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

@ -12,18 +12,18 @@ Josey's thesis, `Development and analysis of high order neutron
transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
:toctree: generated
:nosignatures:
:template: myintegrator.rst
integrator.predictor
integrator.cecm
integrator.celi
integrator.leqi
integrator.cf4
integrator.epc_rk4
integrator.si_celi
integrator.si_leqi
PredictorIntegrator
CECMIntegrator
CELIIntegrator
CF4Integrator
EPCRK4Integrator
LEQIIntegrator
SICELIIntegrator
SILEQIIntegrator
Each of these functions expects a "transport operator" to be passed. An operator
specific to OpenMC is available using the following class:
@ -95,7 +95,18 @@ The following classes are abstract classes that can be used to extend the
EnergyHelper
TransportOperator
Each of the integrator functions also relies on a number of "helper" functions
Custom integrators can be developed by subclassing from the following abstract
base classes:
.. autosummary::
:toctree: generated
:nosignatures:
:template: myintegrator.rst
Integrator
SIIntegrator
Each of the integrator classes also relies on a number of "helper" functions
as follows:
.. autosummary::
@ -103,5 +114,5 @@ as follows:
:nosignatures:
:template: myfunction.rst
integrator.CRAM16
integrator.CRAM48
cram.CRAM16
cram.CRAM48

2
docs/source/usersguide/geometry.rst
View file

@ -25,7 +25,7 @@ surface is a locus of zeros of a function of Cartesian coordinates
Defining a surface alone is not sufficient to specify a volume -- in order to
define an actual volume, one must reference the *half-space* of a surface. A
surface half-space is the region whose points satisfy a positive of negative
surface half-space is the region whose points satisfy a positive or negative
inequality of the surface equation. For example, for a sphere of radius one
centered at the origin, the surface equation is :math:`f(x,y,z) = x^2 + y^2 +
z^2 - 1 = 0`. Thus, we say that the negative half-space of the sphere, is

2
openmc/deplete/__init__.py
View file

@ -38,4 +38,4 @@ from .reaction_rates import *
from .abc import *
from .results import *
from .results_list import *
from .integrator import *
from .integrators import *

79
openmc/deplete/_matrix_funcs.py Normal file
View file

@ -0,0 +1,79 @@
"""Functions to form the special matrix for depletion"""
def celi_f1(chain, rates):
return (5 / 12 * chain.form_matrix(rates[0])
+ 1 / 12 * chain.form_matrix(rates[1]))
def celi_f2(chain, rates):
return (1 / 12 * chain.form_matrix(rates[0])
+ 5 / 12 * chain.form_matrix(rates[1]))
def cf4_f1(chain, rates):
return 1 / 2 * chain.form_matrix(rates)
def cf4_f2(chain, rates):
return -1 / 2 * chain.form_matrix(rates[0]) + chain.form_matrix(rates[1])
def cf4_f3(chain, rates):
return (1 / 4 * chain.form_matrix(rates[0])
+ 1 / 6 * chain.form_matrix(rates[1])
+ 1 / 6 * chain.form_matrix(rates[2])
- 1 / 12 * chain.form_matrix(rates[3]))
def cf4_f4(chain, rates):
return (-1 / 12 * chain.form_matrix(rates[0])
+ 1 / 6 * chain.form_matrix(rates[1])
+ 1 / 6 * chain.form_matrix(rates[2])
+ 1 / 4 * chain.form_matrix(rates[3]))
def rk4_f1(chain, rates):
return 1 / 2 * chain.form_matrix(rates)
def rk4_f4(chain, rates):
return (1 / 6 * chain.form_matrix(rates[0])
+ 1 / 3 * chain.form_matrix(rates[1])
+ 1 / 3 * chain.form_matrix(rates[2])
+ 1 / 6 * chain.form_matrix(rates[3]))
def leqi_f1(chain, inputs):
f1 = chain.form_matrix(inputs[0])
f2 = chain.form_matrix(inputs[1])
dt_l, dt = inputs[2], inputs[3]
return -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2
def leqi_f2(chain, inputs):
f1 = chain.form_matrix(inputs[0])
f2 = chain.form_matrix(inputs[1])
dt_l, dt = inputs[2], inputs[3]
return -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2
def leqi_f3(chain, inputs):
f1 = chain.form_matrix(inputs[0])
f2 = chain.form_matrix(inputs[1])
f3 = chain.form_matrix(inputs[2])
dt_l, dt = inputs[3], inputs[4]
return (-dt ** 2 / (12 * dt_l * (dt + dt_l)) * f1
+ (dt ** 2 + 6 * dt * dt_l + 5 * dt_l ** 2)
/ (12 * dt_l * (dt + dt_l)) * f2 + dt_l / (12 * (dt + dt_l)) * f3)
def leqi_f4(chain, inputs):
f1 = chain.form_matrix(inputs[0])
f2 = chain.form_matrix(inputs[1])
f3 = chain.form_matrix(inputs[2])
dt_l, dt = inputs[3], inputs[4]
return (-dt ** 2 / (12 * dt_l * (dt + dt_l)) * f1
+ (dt ** 2 + 2 * dt * dt_l + dt_l ** 2)
/ (12 * dt_l * (dt + dt_l)) * f2
+ (4 * dt * dt_l + 5 * dt_l ** 2) / (12 * dt_l * (dt + dt_l)) * f3)

327
openmc/deplete/abc.py
View file

@ -5,18 +5,22 @@ to run a full depletion simulation.
"""
from collections import namedtuple
from collections.abc import Iterable
import os
from pathlib import Path
from abc import ABC, abstractmethod
from xml.etree import ElementTree as ET
from copy import deepcopy
from warnings import warn
from numbers import Real
from numbers import Real, Integral
from numpy import nonzero, empty
from uncertainties import ufloat
from openmc.data import DataLibrary, JOULE_PER_EV
from openmc.checkvalue import check_type, check_greater_than
from .results import Results
from .chain import Chain
from .results_list import ResultsList
OperatorResult = namedtuple('OperatorResult', ['k', 'rates'])
OperatorResult.__doc__ = """\
@ -24,8 +28,8 @@ Result of applying transport operator
Parameters
----------
k : float
Resulting eigenvalue
k : uncertainties.ufloat
Resulting eigenvalue and standard deviation
rates : openmc.deplete.ReactionRates
Resulting reaction rates
@ -43,8 +47,8 @@ class TransportOperator(ABC):
Each depletion integrator is written to work with a generic transport
operator that takes a vector of material compositions and returns an
eigenvalue and reaction rates. This abstract class sets the requirements for
such a transport operator. Users should instantiate
eigenvalue and reaction rates. This abstract class sets the requirements
for such a transport operator. Users should instantiate
:class:`openmc.deplete.Operator` rather than this class.
Parameters
@ -61,6 +65,8 @@ class TransportOperator(ABC):
in initial condition to ensure they exist in the decay chain.
Only done for nuclides with reaction rates.
Defaults to 1.0e3.
prev_results : ResultsList, optional
Results from a previous depletion calculation.
Attributes
----------
@ -68,8 +74,12 @@ class TransportOperator(ABC):
Initial atom density [atoms/cm^3] to add for nuclides that are zero
in initial condition to ensure they exist in the decay chain.
Only done for nuclides with reaction rates.
prev_res : ResultsList or None
Results from a previous depletion calculation. ``None`` if no
results are to be used.
"""
def __init__(self, chain_file=None, fission_q=None, dilute_initial=1.0e3):
def __init__(self, chain_file=None, fission_q=None, dilute_initial=1.0e3,
prev_results=None):
self.dilute_initial = dilute_initial
self.output_dir = '.'
@ -93,6 +103,11 @@ class TransportOperator(ABC):
"of adding depletion_chain to OPENMC_CROSS_SECTIONS",
FutureWarning)
self.chain = Chain.from_xml(chain_file, fission_q)
if prev_results is None:
self.prev_res = None
else:
check_type("previous results", prev_results, ResultsList)
self.prev_results = prev_results
@property
def dilute_initial(self):
@ -106,15 +121,15 @@ class TransportOperator(ABC):
self._dilute_initial = value
@abstractmethod
def __call__(self, vec, print_out=True):
def __call__(self, vec, power):
"""Runs a simulation.
Parameters
----------
vec : list of numpy.ndarray
Total atoms to be used in function.
print_out : bool, optional
Whether or not to print out time.
power : float
Power of the reactor in [W]
Returns
-------
@ -122,7 +137,6 @@ class TransportOperator(ABC):
Eigenvalue and reaction rates resulting from transport operator
"""
pass
def __enter__(self):
# Save current directory and move to specific output directory
@ -157,8 +171,6 @@ class TransportOperator(ABC):
Total density for initial conditions.
"""
pass
@abstractmethod
def get_results_info(self):
"""Returns volume list, cell lists, and nuc lists.
@ -170,16 +182,28 @@ class TransportOperator(ABC):
nuc_list : list of str
A list of all nuclide names. Used for sorting the simulation.
burn_list : list of int
A list of all cell IDs to be burned. Used for sorting the simulation.
A list of all cell IDs to be burned. Used for sorting the
simulation.
full_burn_list : list of int
All burnable materials in the geometry.
"""
pass
def finalize(self):
pass
@abstractmethod
def write_bos_data(self, step):
"""Document beginning of step data for a given step
Called at the beginning of a depletion step and at
the final point in the simulation.
Parameters
----------
step : int
Current depletion step including restarts
"""
class ReactionRateHelper(ABC):
"""Abstract class for generating reaction rates for operators
@ -246,7 +270,8 @@ class ReactionRateHelper(ABC):
Parameters
----------
number : iterable of float
Number density [atoms/b-cm] of each nuclide tracked in the calculation.
Number density [atoms/b-cm] of each nuclide tracked in the
calculation.
Returns
-------
@ -337,3 +362,271 @@ class EnergyHelper(ABC):
def nuclides(self, nuclides):
check_type("nuclides", nuclides, list, str)
self._nuclides = nuclides
class Integrator(ABC):
"""Abstract class for solving the time-integration for depletion
Parameters
----------
operator : openmc.deplete.TransportOperator
Operator to perform transport simulations
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not
cumulative.
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.
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`
"""
def __init__(self, operator, timesteps, power=None, power_density=None):
# Check number of stages previously used
if operator.prev_res is not None:
res = operator.prev_res[-1]
if res.data.shape[0] != self._num_stages:
raise ValueError(
"{} incompatible with previous restart calculation. "
"Previous scheme used {} intermediate solutions, while "
"this uses {}".format(
self.__class__.__name__, res.data.shape[0],
self._num_stages))
self.operator = operator
self.chain = operator.chain
if not isinstance(timesteps, Iterable):
self.timesteps = [timesteps]
else:
self.timesteps = timesteps
if power is None:
if power_density is None:
raise ValueError("Either power or power density must be set")
if not isinstance(power_density, Iterable):
power = 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(self.timesteps)
elif len(power) != len(self.timesteps):
raise ValueError(
"Number of time steps != number of powers. {} vs {}".format(
len(self.timesteps), len(power)))
self.power = power
@abstractmethod
def __call__(self, conc, rates, dt, power, i):
"""Perform the integration across one time step
Parameters
----------
conc : numpy.ndarray
Initial concentrations for all nuclides in [atom]
rates : openmc.deplete.ReactionRates
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system in [W]
i : int
Current depletion step index
Returns
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
conc_list : list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulations
"""
@property
@abstractmethod
def _num_stages(self):
"""Number of intermediate transport solutions
Needed to ensure schemes are consistent with restarts
"""
def __iter__(self):
"""Return pairs of time steps in [s] and powers in [W]"""
return zip(self.timesteps, self.power)
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):
"""Get beginning of step concentrations, reaction rates from Operator
"""
x = deepcopy(bos_conc)
res = self.operator(x, step_power)
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):
"""Get beginning of step concentrations, reaction rates from restart"""
res = self.operator.prev_res[-1]
# Depletion methods expect list of arrays
bos_conc = list(res.data[0])
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]
return bos_conc, OperatorResult(k, rates)
def _get_start_data(self):
if self.operator.prev_res is None:
return 0.0, 0
return (self.operator.prev_res[-1].time[-1],
len(self.operator.prev_res) - 1)
def integrate(self):
"""Perform the entire depletion process across all steps"""
with self.operator as conc:
t, self._i_res = self._get_start_data()
for i, (dt, p) in enumerate(self):
if i > 0 or self.operator.prev_res is None:
conc, res = self._get_bos_data_from_operator(i, p, 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)
# Insert BOS concentration, transport results
conc_list.insert(0, conc)
res_list.insert(0, res)
# Remove actual EOS concentration for next step
conc = conc_list.pop()
Results.save(self.operator, conc_list, res_list, [t, t + dt],
p, self._i_res + i, proc_time)
t += dt
# Final simulation
res_list = [self.operator(conc, p)]
Results.save(self.operator, [conc], res_list, [t, t],
p, self._i_res + len(self), proc_time)
self.operator.write_bos_data(len(self) + self._i_res)
class SIIntegrator(Integrator):
"""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
Parameters
----------
operator : openmc.deplete.TransportOperator
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not
cumulative.
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.
n_steps : int, optional
Number of stochastic iterations per depletion interval.
Must be greater than zero. Default : 10
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
"""
def __init__(self, operator, timesteps, power=None, power_density=None,
n_steps=10):
check_type("n_steps", n_steps, Integral)
check_greater_than("n_steps", n_steps, 0)
super().__init__(operator, timesteps, power, power_density)
self.n_steps = n_steps
def _get_bos_data_from_operator(self, step_index, step_power, bos_conc):
reset_particles = False
if step_index == 0 and hasattr(self.operator, "settings"):
reset_particles = True
self.operator.settings.particles *= self.n_stages
inherited = super()._get_bos_data_from_operator(
step_index, step_power, bos_conc)
if reset_particles:
self.operator.settings.particles //= self.n_stages
return inherited
def integrate(self):
"""Perform the entire depletion process across all steps"""
with self.operator as conc:
t, self._i_res = self._get_start_data()
for i, (dt, p) in enumerate(self):
if i == 0:
if self.operator.prev_res is None:
conc, res = self._get_bos_data_from_operator(i, p, conc)
else:
conc, res = self._get_bos_data_from_restart(i, p, conc)
else:
# Pull rates, k from previous iteration w/o
# re-running transport
res = res_list[-1] # defined in previous i iteration
proc_time, conc_list, res_list = self(conc, res.rates, dt, p, i)
# Insert BOS concentration, transport results
conc_list.insert(0, conc)
res_list.insert(0, res)
# Remove actual EOS concentration for next step
conc = conc_list.pop()
Results.save(self.operator, conc_list, res_list, [t, t + dt],
p, self._i_res + i, proc_time)
t += dt
# No final simulation for SIE, use last iteration results
Results.save(self.operator, [conc], [res_list[-1]], [t, t],
p, self._i_res + len(self), proc_time)
self.operator.write_bos_data(self._i_res + len(self))

21
openmc/deplete/integrator/cram.py → openmc/deplete/cram.py
View file

@ -11,10 +11,12 @@ import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as sla
from .. import comm
from . import comm
__all__ = ["deplete", "timed_deplete", "CRAM16", "CRAM48"]
def deplete(chain, x, rates, dt, print_out=True, matrix_func=None):
def deplete(chain, x, rates, dt, matrix_func=None):
"""Deplete materials using given reaction rates for a specified time
Parameters
@ -27,29 +29,22 @@ def deplete(chain, x, rates, dt, print_out=True, matrix_func=None):
Reaction rates (from transport operator)
dt : float
Time in [s] to deplete for
print_out : bool, optional
Whether to show elapsed time
maxtrix_func : function, optional
Function to form the depletion matrix
maxtrix_func : Callable, optional
Function of two variables: ``chain`` and ``rates``.
Expected to return the depletion matrix required by
:func:`CRAM48`.
Returns
-------
x_result : list of numpy.ndarray
Updated atom number vectors for each material
"""
t_start = time.time()
# Use multiprocessing pool to distribute work
with Pool() as pool:
iters = zip(repeat(chain), x, rates, repeat(dt), repeat(matrix_func))
x_result = list(pool.starmap(_cram_wrapper, iters))
t_end = time.time()
if comm.rank == 0:
if print_out:
print("Time to matexp: ", t_end - t_start)
return x_result

16
openmc/deplete/integrator/__init__.py
View file

@ -1,16 +0,0 @@
"""
Integrator
===========
The integrator subcomponents.
"""
from .cf4 import *
from .cecm import *
from .celi import *
from .cram import *
from .epc_rk4 import *
from .leqi import *
from .predictor import *
from .si_celi import *
from .si_leqi import *

123
openmc/deplete/integrator/cecm.py
View file

@ -1,123 +0,0 @@
"""The CE/CM integrator."""
import copy
from collections.abc import Iterable
from .cram import timed_deplete
from ..results import Results
def cecm(operator, timesteps, power=None, power_density=None, print_out=True):
r"""Deplete using the CE/CM algorithm.
Implements the second order `CE/CM predictor-corrector algorithm
<https://doi.org/10.13182/NSE14-92>`_.
"CE/CM" stands for constant extrapolation on predictor and constant
midpoint on corrector. This algorithm is mathematically defined as:
.. math::
\begin{aligned}
y' &= A(y, t) y(t) \\
A_p &= A(y_n, t_n) \\
y_m &= \text{expm}(A_p h/2) y_n \\
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
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not cumulative.
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.
print_out : bool, optional
Whether or not to print out time.
"""
if power is None:
if power_density is None:
raise ValueError(
"Neither power nor power density was specified.")
if not isinstance(power_density, Iterable):
power = power_density*operator.heavy_metal
else:
power = [i*operator.heavy_metal for i in power_density]
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
# Generate initial conditions
with operator as vec:
# Initialize time and starting index
if operator.prev_res is None:
t = 0.0
i_res = 0
else:
t = operator.prev_res[-1].time[-1]
i_res = len(operator.prev_res)
chain = operator.chain
for i, (dt, p) in enumerate(zip(timesteps, power)):
# Get beginning-of-timestep concentrations and reaction rates
# Avoid doing first transport run if already done in previous
# calculation
if i > 0 or operator.prev_res is None:
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], p)]
else:
# Get initial concentration
x = [operator.prev_res[-1].data[0]]
# Get rates
op_results = [operator.prev_res[-1]]
op_results[0].rates = op_results[0].rates[0]
# Set first stage value of keff
op_results[0].k = op_results[0].k[0]
# Scale reaction rates by ratio of powers
power_res = operator.prev_res[-1].power
ratio_power = p / power_res
op_results[0].rates *= ratio_power[0]
# Deplete for first half of timestep
proc_time, x_middle = timed_deplete(
chain, x[0], op_results[0].rates, dt/2, print_out)
# Get middle-of-timestep reaction rates
x.append(x_middle)
op_results.append(operator(x_middle, p))
# Deplete for full timestep using beginning-of-step materials
# and middle-of-timestep reaction rates
pt_end, x_end = timed_deplete(
chain, x[0], op_results[1].rates, dt, print_out)
# Create results, write to disk
Results.save(
operator, x, op_results, [t, t + dt], p, i_res + i,
proc_time + pt_end)
# Advance time, update vector
t += dt
vec = copy.deepcopy(x_end)
# Perform one last simulation
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], power[-1])]
# Create results, write to disk
Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps))

167
openmc/deplete/integrator/celi.py
View file

@ -1,167 +0,0 @@
"""The CE/LI CFQ4 integrator."""
import copy
from collections.abc import Iterable
from .cram import timed_deplete
from ..results import Results
# Functions to form the special matrix for depletion
def _celi_f1(chain, rates):
return 5/12 * chain.form_matrix(rates[0]) + \
1/12 * chain.form_matrix(rates[1])
def _celi_f2(chain, rates):
return 1/12 * chain.form_matrix(rates[0]) + \
5/12 * chain.form_matrix(rates[1])
def celi(operator, timesteps, power=None, power_density=None,
print_out=True):
r"""Deplete using the CE/LI CFQ4 algorithm.
Implements the CE/LI Predictor-Corrector algorithm using the `fourth order
commutator-free integrator <https://doi.org/10.1137/05063042>`_.
"CE/LI" stands for constant extrapolation on predictor and linear
interpolation on corrector. This algorithm is mathematically defined as:
.. math::
\begin{aligned}
y' &= A(y, t) y(t) \\
A_0 &= A(y_n, t_n) \\
y_p &= \text{expm}(h A_0) y_n \\
A_1 &= A(y_p, t_n + h) \\
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
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not cumulative.
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.
print_out : bool, optional
Whether or not to print out time.
"""
if power is None:
if power_density is None:
raise ValueError(
"Neither power nor power density was specified.")
if not isinstance(power_density, Iterable):
power = power_density*operator.heavy_metal
else:
power = [i*operator.heavy_metal for i in power_density]
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
# Generate initial conditions
with operator as vec:
# Initialize time and starting index
if operator.prev_res is None:
t = 0.0
i_res = 0
else:
t = operator.prev_res[-1].time[-1]
i_res = len(operator.prev_res)
for i, (dt, p) in enumerate(zip(timesteps, power)):
vec, t, _ = celi_inner(operator, vec, p, i, i_res, t, dt,
print_out)
# Perform one last simulation
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], power[-1])]
# Create results, write to disk
Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps))
def celi_inner(operator, vec, p, i, i_res, t, dt, print_out):
""" The inner loop of CE/LI CFQ4.
Parameters
----------
operator : Operator
The operator object to simulate on.
x : list of nuclide vector
Nuclide vector, beginning of time.
p : float
Power of the reactor in [W]
i : int
Current iteration number.
i_res : int
Starting index, for restart calculation.
t : float
Time at start of step.
dt : float
Time step.
print_out : bool
Whether or not to print out time.
Returns
-------
list of numpy.array
Nuclide vector, end of time.
float
Next time
OperatorResult
Operator result from beginning of step.
"""
chain = operator.chain
# Get beginning-of-timestep concentrations and reaction rates
# Avoid doing first transport run if already done in previous
# calculation
if i > 0 or operator.prev_res is None:
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], p)]
else:
# Get initial concentration
x = [operator.prev_res[-1].data[0]]
# Get rates
op_results = [operator.prev_res[-1]]
op_results[0].rates = op_results[0].rates[0]
# Set first stage value of keff
op_results[0].k = op_results[0].k[0]
# Scale reaction rates by ratio of powers
power_res = operator.prev_res[-1].power
ratio_power = p / power_res
op_results[0].rates *= ratio_power[0]
# Deplete to end
proc_time, x_new = timed_deplete(chain, x[0], op_results[0].rates, dt, print_out)
x.append(x_new)
op_results.append(operator(x[1], p))
# Deplete with two matrix exponentials
rates = list(zip(op_results[0].rates, op_results[1].rates))
time_1, x_end = timed_deplete(chain, x[0], rates, dt, print_out,
matrix_func=_celi_f1)
time_2, x_end = timed_deplete(chain, x_end, rates, dt, print_out,
matrix_func=_celi_f2)
# Create results, write to disk
Results.save(operator, x, op_results, [t, t + dt], p, i_res + i, proc_time + time_1 + time_2)
# return updated time and vectors
return x_end, t + dt, op_results[0]

164
openmc/deplete/integrator/cf4.py
View file

@ -1,164 +0,0 @@
"""The CF4 integrator."""
import copy
from collections.abc import Iterable
from .cram import timed_deplete
from ..results import Results
# Functions to form the special matrix for depletion
def _cf4_f1(chain, rates):
return 1/2 * chain.form_matrix(rates)
def _cf4_f2(chain, rates):
return -1/2 * chain.form_matrix(rates[0]) + \
chain.form_matrix(rates[1])
def _cf4_f3(chain, rates):
return 1/4 * chain.form_matrix(rates[0]) + \
1/6 * chain.form_matrix(rates[1]) + \
1/6 * chain.form_matrix(rates[2]) + \
-1/12 * chain.form_matrix(rates[3])
def _cf4_f4(chain, rates):
return -1/12 * chain.form_matrix(rates[0]) + \
1/6 * chain.form_matrix(rates[1]) + \
1/6 * chain.form_matrix(rates[2]) + \
1/4 * chain.form_matrix(rates[3])
def cf4(operator, timesteps, power=None, power_density=None, print_out=True):
r"""Deplete using the CF4 algorithm.
Implements the fourth order `commutator-free Lie algorithm
<https://doi.org/10.1016/S0167-739X(02)00161-9>`_.
This algorithm is mathematically defined as:
.. math::
\begin{aligned}
F_1 &= h A(y_0) \\
y_1 &= \text{expm}(1/2 F_1) y_0 \\
F_2 &= h A(y_1) \\
y_2 &= \text{expm}(1/2 F_2) y_0 \\
F_3 &= h A(y_2) \\
y_3 &= \text{expm}(-1/2 F_1 + F_3) y_1 \\
F_4 &= h A(y_3) \\
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
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not cumulative.
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.
print_out : bool, optional
Whether or not to print out time.
"""
if power is None:
if power_density is None:
raise ValueError(
"Neither power nor power density was specified.")
if not isinstance(power_density, Iterable):
power = power_density*operator.heavy_metal
else:
power = [i*operator.heavy_metal for i in power_density]
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
# Generate initial conditions
with operator as vec:
# Initialize time and starting index
if operator.prev_res is None:
t = 0.0
i_res = 0
else:
t = operator.prev_res[-1].time[-1]
i_res = len(operator.prev_res)
chain = operator.chain
for i, (dt, p) in enumerate(zip(timesteps, power)):
# Get beginning-of-timestep concentrations and reaction rates
# Avoid doing first transport run if already done in previous
# calculation
if i > 0 or operator.prev_res is None:
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], p)]
else:
# Get initial concentration
x = [operator.prev_res[-1].data[0]]
# Get rates
op_results = [operator.prev_res[-1]]
op_results[0].rates = op_results[0].rates[0]
# Set first stage value of keff
op_results[0].k = op_results[0].k[0]
# Scale reaction rates by ratio of powers
power_res = operator.prev_res[-1].power
ratio_power = p / power_res
op_results[0].rates *= ratio_power[0]
# Step 1: deplete with matrix 1/2*A(y0)
time_1, x_new = timed_deplete(
chain, x[0], op_results[0].rates, dt, print_out,
matrix_func=_cf4_f1)
x.append(x_new)
op_results.append(operator(x_new, p))
# Step 2: deplete with matrix 1/2*A(y1)
time_2, x_new = timed_deplete(
chain, x[0], op_results[1].rates, dt, print_out,
matrix_func=_cf4_f1)
x.append(x_new)
op_results.append(operator(x_new, p))
# Step 3: deplete with matrix -1/2*A(y0)+A(y2)
rates = list(zip(op_results[0].rates, op_results[2].rates))
time_3, x_new = timed_deplete(
chain, x[1], rates, dt, print_out, matrix_func=_cf4_f2)
x.append(x_new)
op_results.append(operator(x_new, p))
# Step 4: deplete with two matrix exponentials
rates = list(zip(op_results[0].rates, op_results[1].rates,
op_results[2].rates, op_results[3].rates))
time_4, x_end = timed_deplete(
chain, x[0], rates, dt, print_out, matrix_func=_cf4_f3)
time_5, x_end = timed_deplete(
chain, x_end, rates, dt, print_out, matrix_func=_cf4_f4)
# Create results, write to disk
Results.save(
operator, x, op_results, [t, t + dt], p, i_res + i,
time_1 + time_2 + time_3 + time_4 + time_5)
# Advance time, update vector
t += dt
vec = copy.deepcopy(x_end)
# Perform one last simulation
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], power[-1])]
# Create results, write to disk
Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps))

148
openmc/deplete/integrator/epc_rk4.py
View file

@ -1,148 +0,0 @@
"""The EPC-RK4 integrator."""
import copy
from collections.abc import Iterable
from .cram import timed_deplete
from ..results import Results
# Functions to form the special matrix for depletion
def _rk4_f1(chain, rates):
return 1/2 * chain.form_matrix(rates)
def _rk4_f4(chain, rates):
return 1/6 * chain.form_matrix(rates[0]) + \
1/3 * chain.form_matrix(rates[1]) + \
1/3 * chain.form_matrix(rates[2]) + \
1/6 * chain.form_matrix(rates[3])
def epc_rk4(operator, timesteps, power=None, power_density=None, print_out=True):
r"""Deplete using the EPC-RK4 algorithm.
Implements an extended predictor-corrector algorithm with traditional
Runge-Kutta 4 method.
This algorithm is mathematically defined as:
.. math::
\begin{aligned}
F_1 &= h A(y_0) \\
y_1 &= \text{expm}(1/2 F_1) y_0 \\
F_2 &= h A(y_1) \\
y_2 &= \text{expm}(1/2 F_2) y_0 \\
F_3 &= h A(y_2) \\
y_3 &= \text{expm}(F_3) y_0 \\
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
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not cumulative.
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.
print_out : bool, optional
Whether or not to print out time.
"""
if power is None:
if power_density is None:
raise ValueError(
"Neither power nor power density was specified.")
if not isinstance(power_density, Iterable):
power = power_density*operator.heavy_metal
else:
power = [i*operator.heavy_metal for i in power_density]
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
# Generate initial conditions
with operator as vec:
# Initialize time and starting index
if operator.prev_res is None:
t = 0.0
i_res = 0
else:
t = operator.prev_res[-1].time[-1]
i_res = len(operator.prev_res)
chain = operator.chain
for i, (dt, p) in enumerate(zip(timesteps, power)):
# Get beginning-of-timestep concentrations and reaction rates
# Avoid doing first transport run if already done in previous
# calculation
if i > 0 or operator.prev_res is None:
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], p)]
else:
# Get initial concentration
x = [operator.prev_res[-1].data[0]]
# Get rates
op_results = [operator.prev_res[-1]]
op_results[0].rates = op_results[0].rates[0]
# Set first stage value of keff
op_results[0].k = op_results[0].k[0]
# Scale reaction rates by ratio of powers
power_res = operator.prev_res[-1].power
ratio_power = p / power_res
op_results[0].rates *= ratio_power[0]
# Step 1: deplete with matrix 1/2*A(y0)
time_1, x_new = timed_deplete(
chain, x[0], op_results[0].rates, dt, print_out,
matrix_func=_rk4_f1)
x.append(x_new)
op_results.append(operator(x[1], p))
# Step 2: deplete with matrix 1/2*A(y1)
time_2, x_new = timed_deplete(
chain, x[0], op_results[1].rates, dt, print_out,
matrix_func=_rk4_f1)
x.append(x_new)
op_results.append(operator(x[2], p))
# Step 3: deplete with matrix A(y2)
time_3, x_new = timed_deplete(
chain, x[0], op_results[2].rates, dt, print_out)
x.append(x_new)
op_results.append(operator(x[3], p))
# Step 4: deplete with matrix 1/6*A(y0)+1/3*A(y1)+1/3*A(y2)+1/6*A(y3)
rates = list(zip(op_results[0].rates, op_results[1].rates,
op_results[2].rates, op_results[3].rates))
time_4, x_end = timed_deplete(
chain, x[0], rates, dt, print_out, matrix_func=_rk4_f4)
# Create results, write to disk
Results.save(
operator, x, op_results, [t, t + dt], p, i_res + i,
time_1 + time_2 + time_3 + time_4)
# Advance time, update vector
t += dt
vec = copy.deepcopy(x_end)
# Perform one last simulation
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], power[-1])]
# Create results, write to disk
Results.save(
operator, x, op_results, [t, t], p, i_res + len(timesteps))

170
openmc/deplete/integrator/leqi.py
View file

@ -1,170 +0,0 @@
"""The LE/QI CFQ4 integrator."""
import copy
from collections.abc import Iterable
from itertools import repeat
from .celi import celi_inner
from .cram import timed_deplete
from ..results import Results
# Functions to form the special matrix for depletion
def _leqi_f1(chain, inputs):
f1 = chain.form_matrix(inputs[0])
f2 = chain.form_matrix(inputs[1])
dt_l, dt = inputs[2], inputs[3]
return -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2
def _leqi_f2(chain, inputs):
f1 = chain.form_matrix(inputs[0])
f2 = chain.form_matrix(inputs[1])
dt_l, dt = inputs[2], inputs[3]
return -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2
def _leqi_f3(chain, inputs):
f1 = chain.form_matrix(inputs[0])
f2 = chain.form_matrix(inputs[1])
f3 = chain.form_matrix(inputs[2])
dt_l, dt = inputs[3], inputs[4]
return -dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + \
(dt**2 + 6*dt*dt_l + 5*dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + \
dt_l / (12 * (dt + dt_l)) * f3
def _leqi_f4(chain, inputs):
f1 = chain.form_matrix(inputs[0])
f2 = chain.form_matrix(inputs[1])
f3 = chain.form_matrix(inputs[2])
dt_l, dt = inputs[3], inputs[4]
return -dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + \
(dt**2 + 2*dt*dt_l + dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + \
(4 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f3
def leqi(operator, timesteps, power=None, power_density=None, print_out=True):
r"""Deplete using the LE/QI CFQ4 algorithm.
Implements the LE/QI Predictor-Corrector algorithm using the `fourth order
commutator-free integrator <https://doi.org/10.1137/05063042>`_.
"LE/QI" stands for linear extrapolation on predictor and quadratic
interpolation on corrector. This algorithm is mathematically defined as:
.. math::
\begin{aligned}
y' &= A(y, t) y(t) \\
A_{last} &= A(y_{n-1}, t_n - h_1) \\
A_0 &= A(y_n, t_n) \\
F_1 &= \frac{-h_2^2}{12h_1} A_{last} + \frac{h_2(6h_1+h_2)}{12h_1} A_0 \\
F_2 &= \frac{-5h_2^2}{12h_1} A_{last} + \frac{h_2(6h_1+5h_2)}{12h_1} A_0 \\
y_p &= \text{expm}(F_2) \text{expm}(F_1) y_n \\
A_1 &= A(y_p, t_n + h_2) \\
F_3 &= \frac{-h_2^3}{12 h_1 (h_1 + h_2)} A_{last} +
\frac{h_2 (5 h_1^2 + 6 h_2 h_1 + h_2^2)}{12 h_1 (h_1 + h_2)} A_0 +
\frac{h_2 h_1)}{12 (h_1 + h_2)} A_1 \\
F_4 &= \frac{-h_2^3}{12 h_1 (h_1 + h_2)} A_{last} +
\frac{h_2 (h_1^2 + 2 h_2 h_1 + h_2^2)}{12 h_1 (h_1 + h_2)} A_0 +
\frac{h_2 (5 h_1^2 + 4 h_2 h_1)}{12 h_1 (h_1 + h_2)} A_1 \\
y_{n+1} &= \text{expm}(F_4) \text{expm}(F_3) y_n
\end{aligned}
It is initialized using the CE/LI algorithm.
Parameters
----------
operator : openmc.deplete.TransportOperator
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not cumulative.
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.
print_out : bool, optional
Whether or not to print out time.
"""
if power is None:
if power_density is None:
raise ValueError(
"Neither power nor power density was specified.")
if not isinstance(power_density, Iterable):
power = power_density*operator.heavy_metal
else:
power = [i*operator.heavy_metal for i in power_density]
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
# Generate initial conditions
with operator as vec:
# Initialize time and starting index
if operator.prev_res is None:
t = 0.0
i_res = 0
else:
t = operator.prev_res[-1].time[-1]
i_res = len(operator.prev_res)
chain = operator.chain
for i, (dt, p) in enumerate(zip(timesteps, power)):
# LE/QI needs the last step results to start
# Perform CE/LI CFQ4 or restore results for the first step
if i == 0:
if i_res <= 1:
dt_l = dt
x_new, t, op_res_last = celi_inner(operator, vec, p, i,
i_res, t, dt, print_out)
continue
else:
dt_l = t - operator.prev_res[-2].time[0]
op_res_last = operator.prev_res[-2]
op_res_last.rates = op_res_last.rates[0]
x_new = operator.prev_res[-1].data[0]
# Perform remaining LE/QI
x = [copy.deepcopy(x_new)]
op_results = [operator(x[0], p)]
inputs = list(zip(op_res_last.rates, op_results[0].rates,
repeat(dt_l), repeat(dt)))
time_1, x_new = timed_deplete(
chain, x[0], inputs, dt, print_out, matrix_func=_leqi_f1)
time_2, x_new = timed_deplete(
chain, x_new, inputs, dt, print_out, matrix_func=_leqi_f2)
x.append(x_new)
op_results.append(operator(x[1], p))
inputs = list(zip(op_res_last.rates, op_results[0].rates,
op_results[1].rates, repeat(dt_l), repeat(dt)))
time_3, x_new = timed_deplete(
chain, x[0], inputs, dt, print_out, matrix_func=_leqi_f3)
time_4, x_new = timed_deplete(
chain, x_new, inputs, dt, print_out, matrix_func=_leqi_f4)
# Create results, write to disk
Results.save(
operator, x, op_results, [t, t+dt], p, i_res+i,
time_1 + time_2 + time_3 + time_4)
# update results
op_res_last = copy.deepcopy(op_results[0])
t += dt
dt_l = dt
# Perform one last simulation
x = [copy.deepcopy(x_new)]
op_results = [operator(x[0], power[-1])]
# Create results, write to disk
Results.save(
operator, x, op_results, [t, t], p, i_res + len(timesteps))

106
openmc/deplete/integrator/predictor.py
View file

@ -1,106 +0,0 @@
"""First-order predictor algorithm."""
import copy
from collections.abc import Iterable
from .cram import timed_deplete
from ..results import Results
def predictor(operator, timesteps, power=None, power_density=None,
print_out=True):
r"""Deplete using a first-order predictor algorithm.
Implements the first-order predictor algorithm. This algorithm is
mathematically defined as:
.. math::
\begin{aligned}
y' &= A(y, t) y(t) \\
A_p &= A(y_n, t_n) \\
y_{n+1} &= \text{expm}(A_p h) y_n
\end{aligned}
Parameters
----------
operator : openmc.deplete.TransportOperator
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not cumulative.
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.
print_out : bool, optional
Whether or not to print out time.
"""
if power is None:
if power_density is None:
raise ValueError(
"Neither power nor power density was specified.")
if not isinstance(power_density, Iterable):
power = power_density*operator.heavy_metal
else:
power = [i*operator.heavy_metal for i in power_density]
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
proc_time = None
# Generate initial conditions
with operator as vec:
# Initialize time and starting index
if operator.prev_res is None:
t = 0.0
i_res = 0
else:
t = operator.prev_res[-1].time[-1]
i_res = len(operator.prev_res) - 1
chain = operator.chain
for i, (dt, p) in enumerate(zip(timesteps, power)):
# Get beginning-of-timestep concentrations and reaction rates
# Avoid doing first transport run if already done in previous
# calculation
if i > 0 or operator.prev_res is None:
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], p)]
# Create results, write to disk
Results.save(operator, x, op_results, [t, t + dt], p, i_res + i, proc_time)
else:
# Get initial concentration
x = [operator.prev_res[-1].data[0]]
# Get rates
op_results = [operator.prev_res[-1]]
op_results[0].rates = op_results[0].rates[0]
# Scale reaction rates by ratio of powers
power_res = operator.prev_res[-1].power
ratio_power = p / power_res
op_results[0].rates *= ratio_power[0]
# Deplete for full timestep
proc_time, x_end = timed_deplete(
chain, x[0], op_results[0].rates, dt, print_out)
# Advance time, update vector
t += dt
vec = copy.deepcopy(x_end)
# Perform one last simulation
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], power[-1])]
# Create results, write to disk
Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps), proc_time)

166
openmc/deplete/integrator/si_celi.py
View file

@ -1,166 +0,0 @@
"""The SI-CE/LI CFQ4 integrator."""
import copy
from collections.abc import Iterable
from .cram import timed_deplete
from ..results import Results
from ..abc import OperatorResult
from .celi import _celi_f1, _celi_f2
def si_celi(operator, timesteps, power=None, power_density=None,
print_out=True, m=10):
r"""Deplete using the SI-CE/LI CFQ4 algorithm.
Implements the Stochastic Implicit CE/LI Predictor-Corrector algorithm using
the `fourth order commutator-free integrator <https://doi.org/10.1137/05063042>`_.
Detailed algorithm can be found in Section 3.2 in `Colin Josey's thesis
<http://hdl.handle.net/1721.1/113721>`_.
Parameters
----------
operator : openmc.deplete.TransportOperator
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not cumulative.
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.
print_out : bool, optional
Whether or not to print out time.
m : int, optional
Number of stages.
"""
if power is None:
if power_density is None:
raise ValueError(
"Neither power nor power density was specified.")
if not isinstance(power_density, Iterable):
power = power_density*operator.heavy_metal
else:
power = [i*operator.heavy_metal for i in power_density]
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
# Generate initial conditions
with operator as vec:
# Initialize time and starting index
if operator.prev_res is None:
t = 0.0
i_res = 0
else:
t = operator.prev_res[-1].time[-1]
i_res = len(operator.prev_res)
# Get the concentrations and reaction rates for the first
# beginning-of-timestep (BOS). Compute with m (stage number) times as
# many neutrons as later simulations for statistics reasons if no
# previous calculation results present
if operator.prev_res is None:
x = [copy.deepcopy(vec)]
if hasattr(operator, "settings"):
operator.settings.particles *= m
op_results = [operator(x[0], power[0])]
if hasattr(operator, "settings"):
operator.settings.particles //= m
else:
# Get initial concentration
x = [operator.prev_res[-1].data[0]]
# Get rates
op_results = [operator.prev_res[-1]]
op_results[0].rates = op_results[0].rates[0]
# Set first stage value of keff
op_results[0].k = op_results[0].k[0]
# Scale reaction rates by ratio of powers
power_res = operator.prev_res[-1].power
ratio_power = power[0] / power_res
op_results[0].rates *= ratio_power[0]
for i, (dt, p) in enumerate(zip(timesteps, power)):
x, t, op_results = si_celi_inner(operator, x, op_results, p,
i, i_res, t, dt, print_out, m)
# Create results for last point, write to disk
Results.save(
operator, x, op_results, [t, t], p, i_res + len(timesteps))
def si_celi_inner(operator, x, op_results, p, i, i_res, t, dt, print_out, m=10):
""" The inner loop of SI-CE/LI CFQ4.
Parameters
----------
operator : Operator
The operator object to simulate on.
x : list of nuclide vector
Nuclide vector, beginning of time.
op_results : list of OperatorResult
Operator result at BOS.
p : float
Power of the reactor in [W]
i : int
Current iteration number.
i_res : int
Starting index, for restart calculation.
t : float
Time at start of step.
dt : float
Time step.
print_out : bool
Whether or not to print out time.
m : int, optional
Number of stages.
Returns
-------
list of nuclide vector (numpy.array)
Nuclide vector, end of time.
float
Next time
list of OperatorResult
Operator result at end of time.
"""
chain = operator.chain
# Deplete to end
proc_time, x_new = timed_deplete(
chain, x[0], op_results[0].rates, dt, print_out)
x.append(x_new)
for j in range(m + 1):
op_res = operator(x_new, p)
if j <= 1:
op_res_bar = copy.deepcopy(op_res)
else:
rates = 1/j * op_res.rates + (1 - 1/j) * op_res_bar.rates
k = 1/j * op_res.k + (1 - 1/j) * op_res_bar.k
op_res_bar = OperatorResult(k, rates)
rates = list(zip(op_results[0].rates, op_res_bar.rates))
time_1, x_new = timed_deplete(
chain, x[0], rates, dt, print_out, matrix_func=_celi_f1)
time_2, x_new = timed_deplete(
chain, x_new, rates, dt, print_out, matrix_func=_celi_f2)
proc_time += time_1 + time_2
# Create results, write to disk
op_results.append(op_res_bar)
Results.save(operator, x, op_results, [t, t+dt], p, i_res+i, proc_time)
# return updated time and vectors
return [x_new], t + dt, [op_res_bar]

157
openmc/deplete/integrator/si_leqi.py
View file

@ -1,157 +0,0 @@
"""The SI-LE/QI CFQ4 integrator."""
import copy
from collections.abc import Iterable
from itertools import repeat
from .si_celi import si_celi_inner
from .leqi import _leqi_f1, _leqi_f2, _leqi_f3, _leqi_f4
from .cram import timed_deplete
from ..results import Results
from ..abc import OperatorResult
def si_leqi(operator, timesteps, power=None, power_density=None,
print_out=True, m=10):
r"""Deplete using the SI-LE/QI CFQ4 algorithm.
Implements the Stochastic Implicit LE/QI Predictor-Corrector algorithm using
the `fourth order commutator-free integrator <https://doi.org/10.1137/05063042>`_.
Detailed algorithm can be found in Section 3.2 in `Colin Josey's thesis
<http://hdl.handle.net/1721.1/113721>`_.
Parameters
----------
operator : openmc.deplete.TransportOperator
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not cumulative.
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.
print_out : bool, optional
Whether or not to print out time.
m : int, optional
Number of stages.
"""
if power is None:
if power_density is None:
raise ValueError(
"Neither power nor power density was specified.")
if not isinstance(power_density, Iterable):
power = power_density*operator.heavy_metal
else:
power = [i*operator.heavy_metal for i in power_density]
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
# Generate initial conditions
with operator as vec:
# Initialize time and starting index
if operator.prev_res is None:
t = 0.0
i_res = 0
else:
t = operator.prev_res[-1].time[-1]
i_res = len(operator.prev_res)
# Get the concentrations and reaction rates for the first
# beginning-of-timestep (BOS). Compute with m (stage number) times as
# many neutrons as later simulations for statistics reasons if no
# previous calculation results present
if operator.prev_res is None:
x = [copy.deepcopy(vec)]
if hasattr(operator, "settings"):
operator.settings.particles *= m
op_results = [operator(x[0], power[0])]
if hasattr(operator, "settings"):
operator.settings.particles //= m
else:
# Get initial concentration
x = [operator.prev_res[-1].data[0]]
# Get rates
op_results = [operator.prev_res[-1]]
op_results[0].rates = op_results[0].rates[0]
# Set first stage value of keff
op_results[0].k = op_results[0].k[0]
# Scale reaction rates by ratio of powers
power_res = operator.prev_res[-1].power
ratio_power = power[0] / power_res
op_results[0].rates *= ratio_power[0]
chain = operator.chain
for i, (dt, p) in enumerate(zip(timesteps, power)):
# LE/QI needs the last step results to start
# Perform SI-CE/LI CFQ4 or restore results for the first step
if i == 0:
dt_l = dt
if i_res <= 1:
op_res_last = copy.deepcopy(op_results[0])
x, t, op_results = si_celi_inner(operator, x, op_results, p,
i, i_res, t, dt, print_out)
continue
else:
dt_l = t - operator.prev_res[-2].time[0]
op_res_last = operator.prev_res[-2]
op_res_last.rates = op_res_last.rates[0]
x = [operator.prev_res[-1].data[0]]
# Perform remaining LE/QI
inputs = list(zip(op_res_last.rates, op_results[0].rates,
repeat(dt_l), repeat(dt)))
proc_time, x_new = timed_deplete(
chain, x[0], inputs, dt, print_out, matrix_func=_leqi_f1)
time_1, x_new = timed_deplete(
chain, x_new, inputs, dt, print_out, matrix_func=_leqi_f2)
x.append(x_new)
proc_time += time_1
# Loop on inner
for j in range(m + 1):
op_res = operator(x_new, p)
if j <= 1:
op_res_bar = copy.deepcopy(op_res)
else:
rates = 1/j * op_res.rates + (1 - 1/j) * op_res_bar.rates
k = 1/j * op_res.k + (1 - 1/j) * op_res_bar.k
op_res_bar = OperatorResult(k, rates)
inputs = list(zip(op_res_last.rates, op_results[0].rates,
op_res_bar.rates, repeat(dt_l), repeat(dt)))
time_1, x_new = timed_deplete(
chain, x[0], inputs, dt, print_out, matrix_func=_leqi_f3)
time_2, x_new = timed_deplete(
chain, x_new, inputs, dt, print_out, matrix_func=_leqi_f4)
proc_time += time_1 + time_2
# Create results, write to disk
op_results.append(op_res_bar)
Results.save(
operator, x, op_results, [t, t+dt], p, i_res+i, proc_time)
# update results
x = [x_new]
op_res_last = copy.deepcopy(op_results[0])
op_results = [op_res_bar]
t += dt
dt_l = dt
# Create results for last point, write to disk
Results.save(
operator, x, op_results, [t, t], p, i_res+len(timesteps))

836
openmc/deplete/integrators.py Normal file
View file

@ -0,0 +1,836 @@
import copy
from itertools import repeat
from .abc import Integrator, SIIntegrator, OperatorResult
from .cram import timed_deplete
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
)
__all__ = [
"PredictorIntegrator", "CECMIntegrator", "CF4Integrator",
"CELIIntegrator", "EPCRK4Integrator", "LEQIIntegrator",
"SICELIIntegrator", "SILEQIIntegrator"]
class PredictorIntegrator(Integrator):
r"""Deplete using a first-order predictor algorithm.
Implements the first-order predictor algorithm. This algorithm is
mathematically defined as:
.. math::
\begin{aligned}
y' &= A(y, t) y(t) \\
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
Array of timesteps in units of [s]. Note that values are not
cumulative.
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.
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`
"""
_num_stages = 1
def __call__(self, conc, rates, dt, power, _i=None):
"""Perform the integration across one time step
Parameters
----------
conc : numpy.ndarray
Initial concentrations for all nuclides in [atom]
rates : openmc.deplete.ReactionRates
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system in [W]
_i : int or None
Iteration index. Not used
Returns
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
conc_list : list of numpy.ndarray
Concentrations at end of interval
op_results : empty list
Kept for consistency with API. No intermediate calls to
operator with predictor
"""
proc_time, conc_end = timed_deplete(self.chain, conc, rates, dt)
return proc_time, [conc_end], []
class CECMIntegrator(Integrator):
r"""Deplete using the CE/CM algorithm.
Implements the second order `CE/CM predictor-corrector algorithm
<https://doi.org/10.13182/NSE14-92>`_.
"CE/CM" stands for constant extrapolation on predictor and constant
midpoint on corrector. This algorithm is mathematically defined as:
.. math::
\begin{aligned}
y' &= A(y, t) y(t) \\
A_p &= A(y_n, t_n) \\
y_m &= \text{expm}(A_p h/2) y_n \\
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
Array of timesteps in units of [s]. Note that values are not
cumulative.
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.
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`
"""
_num_stages = 2
def __call__(self, conc, rates, dt, power, _i=None):
"""Integrate using CE/CM
Parameters
----------
conc : numpy.ndarray
Initial concentrations for all nuclides in [atom]
rates : openmc.deplete.ReactionRates
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system [W]
_i : int, optional
Current iteration count. Not used
Returns
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
conc_list : list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from transport simulations
"""
# deplete across first half of inteval
time0, x_middle = timed_deplete(self.chain, conc, rates, dt / 2)
res_middle = self.operator(x_middle, power)
# deplete across entire interval with BOS concentrations,
# MOS reaction rates
time1, x_end = timed_deplete(self.chain, conc, res_middle.rates, dt)
return time0 + time1, [x_middle, x_end], [res_middle]
class CF4Integrator(Integrator):
r"""Deplete using the CF4 algorithm.
Implements the fourth order `commutator-free Lie algorithm
<https://doi.org/10.1016/S0167-739X(02)00161-9>`_.
This algorithm is mathematically defined as:
.. math::
\begin{aligned}
F_1 &= h A(y_0) \\
y_1 &= \text{expm}(1/2 F_1) y_0 \\
F_2 &= h A(y_1) \\
y_2 &= \text{expm}(1/2 F_2) y_0 \\
F_3 &= h A(y_2) \\
y_3 &= \text{expm}(-1/2 F_1 + F_3) y_1 \\
F_4 &= h A(y_3) \\
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
Array of timesteps in units of [s]. Note that values are not
cumulative.
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.
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`
"""
_num_stages = 4
def __call__(self, bos_conc, bos_rates, dt, power, _i=None):
"""Perform the integration across one time step
Parameters
----------
bos_conc : numpy.ndarray
Initial concentrations for all nuclides in [atom]
bos_rates : openmc.deplete.ReactionRates
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system in [W]
_i : int, optional
Current depletion step index. Not used
Returns
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
conc_list : list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulations
"""
# Step 1: deplete with matrix 1/2*A(y0)
time1, conc_eos1 = timed_deplete(
self.chain, bos_conc, bos_rates, dt, matrix_func=cf4_f1)
res1 = self.operator(conc_eos1, power)
# Step 2: deplete with matrix 1/2*A(y1)
time2, conc_eos2 = timed_deplete(
self.chain, bos_conc, res1.rates, dt, matrix_func=cf4_f1)
res2 = self.operator(conc_eos2, power)
# Step 3: deplete with matrix -1/2*A(y0)+A(y2)
list_rates = list(zip(bos_rates, res2.rates))
time3, conc_eos3 = timed_deplete(
self.chain, conc_eos1, list_rates, dt, matrix_func=cf4_f2)
res3 = self.operator(conc_eos3, power)
# Step 4: deplete with two matrix exponentials
list_rates = list(zip(bos_rates, res1.rates, res2.rates, res3.rates))
time4, conc_inter = timed_deplete(
self.chain, bos_conc, list_rates, dt, matrix_func=cf4_f3)
time5, conc_eos5 = timed_deplete(
self.chain, conc_inter, list_rates, dt, matrix_func=cf4_f4)
return (time1 + time2 + time3 + time4 + time5,
[conc_eos1, conc_eos2, conc_eos3, conc_eos5],
[res1, res2, res3])
class CELIIntegrator(Integrator):
r"""Deplete using the CE/LI CFQ4 algorithm.
Implements the CE/LI Predictor-Corrector algorithm using the `fourth order
commutator-free integrator <https://doi.org/10.1137/05063042>`_.
"CE/LI" stands for constant extrapolation on predictor and linear
interpolation on corrector. This algorithm is mathematically defined as:
.. math::
\begin{aligned}
y' &= A(y, t) y(t) \\
A_0 &= A(y_n, t_n) \\
y_p &= \text{expm}(h A_0) y_n \\
A_1 &= A(y_p, t_n + h) \\
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
Array of timesteps in units of [s]. Note that values are not
cumulative.
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.
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`
"""
_num_stages = 2
def __call__(self, bos_conc, rates, dt, power, _i=None):
"""Perform the integration across one time step
Parameters
----------
bos_conc : numpy.ndarray
Initial concentrations for all nuclides in [atom]
rates : openmc.deplete.ReactionRates
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system in [W]
_i : int, optional
Current iteration count. Not used
Returns
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
conc_list : list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulation
"""
# deplete to end using BOS rates
proc_time, conc_ce = timed_deplete(self.chain, bos_conc, rates, dt)
res_ce = self.operator(conc_ce, power)
# deplete using two matrix exponentials
list_rates = list(zip(rates, res_ce.rates))
time_le1, conc_inter = timed_deplete(
self.chain, bos_conc, list_rates, dt, matrix_func=celi_f1)
time_le2, conc_end = timed_deplete(
self.chain, conc_inter, list_rates, dt, matrix_func=celi_f2)
return proc_time + time_le1 + time_le1, [conc_ce, conc_end], [res_ce]
class EPCRK4Integrator(Integrator):
r"""Deplete using the EPC-RK4 algorithm.
Implements an extended predictor-corrector algorithm with traditional
Runge-Kutta 4 method. This algorithm is mathematically defined as:
.. math::
\begin{aligned}
F_1 &= h A(y_0) \\
y_1 &= \text{expm}(1/2 F_1) y_0 \\
F_2 &= h A(y_1) \\
y_2 &= \text{expm}(1/2 F_2) y_0 \\
F_3 &= h A(y_2) \\
y_3 &= \text{expm}(F_3) y_0 \\
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
Array of timesteps in units of [s]. Note that values are not
cumulative.
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.
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`
"""
_num_stages = 4
def __call__(self, conc, rates, dt, power, _i=None):
"""Perform the integration across one time step
Parameters
----------
conc : numpy.ndarray
Initial concentrations for all nuclides in [atom]
rates : openmc.deplete.ReactionRates
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system in [W]
_i : int, optional
Current depletion step index, unused.
Returns
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
conc_list : list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulations
"""
# Step 1: deplete with matrix A(y0) / 2
time1, conc1 = timed_deplete(
self.chain, conc, rates, dt, matrix_func=rk4_f1)
res1 = self.operator(conc1, power)
# Step 2: deplete with matrix A(y1) / 2
time2, conc2 = timed_deplete(
self.chain, conc, res1.rates, dt, matrix_func=rk4_f1)
res2 = self.operator(conc2, power)
# Step 3: deplete with matrix A(y2)
time3, conc3 = timed_deplete(
self.chain, conc, res2.rates, dt)
res3 = self.operator(conc3, power)
# Step 4: deplete with matrix built from weighted rates
list_rates = list(zip(rates, res1.rates, res2.rates, res3.rates))
time4, conc4 = timed_deplete(
self.chain, conc, list_rates, dt, matrix_func=rk4_f4)
return (time1 + time2 + time3 + time4, [conc1, conc2, conc3, conc4],
[res1, res2, res3])
class LEQIIntegrator(Integrator):
r"""Deplete using the LE/QI CFQ4 algorithm.
Implements the LE/QI Predictor-Corrector algorithm using the `fourth order
commutator-free integrator <https://doi.org/10.1137/05063042>`_.
"LE/QI" stands for linear extrapolation on predictor and quadratic
interpolation on corrector. This algorithm is mathematically defined as:
.. math::
\begin{aligned}
y' &= A(y, t) y(t) \\
A_{last} &= A(y_{n-1}, t_n - h_1) \\
A_0 &= A(y_n, t_n) \\
F_1 &= \frac{-h_2^2}{12h_1} A_{last} + \frac{h_2(6h_1+h_2)}{12h_1} A_0 \\
F_2 &= \frac{-5h_2^2}{12h_1} A_{last} + \frac{h_2(6h_1+5h_2)}{12h_1} A_0 \\
y_p &= \text{expm}(F_2) \text{expm}(F_1) y_n \\
A_1 &= A(y_p, t_n + h_2) \\
F_3 &= \frac{-h_2^3}{12 h_1 (h_1 + h_2)} A_{last} +
\frac{h_2 (5 h_1^2 + 6 h_2 h_1 + h_2^2)}{12 h_1 (h_1 + h_2)} A_0 +
\frac{h_2 h_1)}{12 (h_1 + h_2)} A_1 \\
F_4 &= \frac{-h_2^3}{12 h_1 (h_1 + h_2)} A_{last} +
\frac{h_2 (h_1^2 + 2 h_2 h_1 + h_2^2)}{12 h_1 (h_1 + h_2)} A_0 +
\frac{h_2 (5 h_1^2 + 4 h_2 h_1)}{12 h_1 (h_1 + h_2)} A_1 \\
y_{n+1} &= \text{expm}(F_4) \text{expm}(F_3) y_n
\end{aligned}
It is initialized using the CE/LI algorithm.
Parameters
----------
operator : openmc.deplete.TransportOperator
Operator to perform transport simulations
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not
cumulative.
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.
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`
"""
_num_stages = 2
def __call__(self, bos_conc, bos_rates, dt, power, i):
"""Perform the integration across one time step
Parameters
----------
conc : numpy.ndarray
Initial concentrations for all nuclides in [atom]
rates : openmc.deplete.ReactionRates
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system in [W]
i : int
Current depletion step index
Returns
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
conc_list : list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulation
"""
if i == 0:
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)
prev_res = self.operator.prev_res[-2]
prev_dt = self.timesteps[i] - prev_res.time[0]
self._prev_rates = prev_res.rates[0]
else:
prev_dt = self.timesteps[i - 1]
# Remaining LE/QI
bos_res = self.operator(bos_conc, power)
le_inputs = list(zip(
self._prev_rates, bos_res.rates, repeat(prev_dt), repeat(dt)))
time1, conc_inter = timed_deplete(
self.chain, bos_conc, le_inputs, dt, matrix_func=leqi_f1)
time2, conc_eos0 = timed_deplete(
self.chain, conc_inter, le_inputs, dt, matrix_func=leqi_f2)
res_inter = self.operator(conc_eos0, power)
qi_inputs = list(zip(
self._prev_rates, bos_res.rates, res_inter.rates,
repeat(prev_dt), repeat(dt)))
time3, conc_inter = timed_deplete(
self.chain, bos_conc, qi_inputs, dt, matrix_func=leqi_f3)
time4, conc_eos1 = timed_deplete(
self.chain, conc_inter, qi_inputs, dt, matrix_func=leqi_f4)
# store updated rates
self._prev_rates = copy.deepcopy(bos_res.rates)
return (
time1 + time2 + time3 + time4, [conc_eos0, conc_eos1],
[bos_res, res_inter])
class SICELIIntegrator(SIIntegrator):
r"""Deplete using the SI-CE/LI CFQ4 algorithm.
Implements the stochastic implicit CE/LI predictor-corrector algorithm
using the `fourth order commutator-free integrator
<https://doi.org/10.1137/05063042>`_.
Detailed algorithm can be found in section 3.2 in `Colin Josey's thesis
<http://hdl.handle.net/1721.1/113721>`_.
Parameters
----------
operator : openmc.deplete.TransportOperator
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not
cumulative.
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.
n_steps : int, optional
Number of stochastic iterations per depletion interval.
Must be greater than zero. Default : 10
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
"""
_num_stages = 2
def __call__(self, bos_conc, bos_rates, dt, power, _i=None):
"""Perform the integration across one time step
Parameters
----------
bos_conc : numpy.ndarray
Initial bos_concentrations for all nuclides in [atom]
bos_rates : openmc.deplete.ReactionRates
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system in [W]
_i : int, optional
Current depletion step index. Not used
Returns
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
bos_conc_list : list of numpy.ndarray
Concentrations at each of the intermediate points with
the final bos_concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulations
"""
proc_time, eos_conc = timed_deplete(
self.chain, bos_conc, bos_rates, dt)
inter_conc = copy.deepcopy(eos_conc)
# Begin iteration
for j in range(self.n_steps + 1):
inter_res = self.operator(inter_conc, power)
if j <= 1:
res_bar = copy.deepcopy(inter_res)
else:
rates = 1/j * inter_res.rates + (1 - 1 / j) * res_bar.rates
k = 1/j * inter_res.k + (1 - 1 / j) * res_bar.k
res_bar = OperatorResult(k, rates)
list_rates = list(zip(bos_rates, res_bar.rates))
time1, inter_conc = timed_deplete(
self.chain, bos_conc, list_rates, dt, matrix_func=celi_f1)
time2, inter_conc = timed_deplete(
self.chain, inter_conc, list_rates, dt, matrix_func=celi_f2)
proc_time += time1 + time2
# end iteration
return proc_time, [eos_conc, inter_conc], [res_bar]
class SILEQIIntegrator(SIIntegrator):
r"""Deplete using the SI-LE/QI CFQ4 algorithm.
Implements the Stochastic Implicit LE/QI Predictor-Corrector algorithm
using the `fourth order commutator-free integrator
<https://doi.org/10.1137/05063042>`_.
Detailed algorithm can be found in Section 3.2 in `Colin Josey's thesis
<http://hdl.handle.net/1721.1/113721>`_.
Parameters
----------
operator : openmc.deplete.TransportOperator
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not
cumulative.
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.
n_steps : int, optional
Number of stochastic iterations per depletion interval.
Must be greater than zero. Default : 10
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
"""
_num_stages = 2
def __call__(self, bos_conc, bos_rates, dt, power, i):
"""Perform the integration across one time step
Parameters
----------
bos_conc : list of numpy.ndarray
Initial concentrations for all nuclides in [atom] for
all depletable materials
bos_rates : list of openmc.deplete.ReactionRates
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]
i : int
Current depletion step index
Returns
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
conc_list : list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulation
"""
if i == 0:
if self._i_res < 1:
self._prev_rates = bos_rates
# Perform CELI for initial steps
return SICELIIntegrator.__call__(
self, bos_conc, bos_rates, dt, power, i)
prev_res = self.operator.prev_res[-2]
prev_dt = self.timesteps[i] - prev_res.time[0]
self._prev_rates = prev_res.rates[0]
else:
prev_dt = self.timesteps[i - 1]
# Perform remaining LE/QI
inputs = list(zip(self._prev_rates, bos_rates,
repeat(prev_dt), repeat(dt)))
proc_time, inter_conc = timed_deplete(
self.chain, bos_conc, inputs, dt, matrix_func=leqi_f1)
time1, eos_conc = timed_deplete(
self.chain, inter_conc, inputs, dt, matrix_func=leqi_f2)
proc_time += time1
inter_conc = copy.deepcopy(eos_conc)
for j in range(self.n_steps + 1):
inter_res = self.operator(inter_conc, power)
if j <= 1:
res_bar = copy.deepcopy(inter_res)
else:
rates = 1 / j * inter_res.rates + (1 - 1 / j) * res_bar.rates
k = 1 / j * inter_res.k + (1 - 1 / j) * res_bar.k
res_bar = OperatorResult(k, rates)
inputs = list(zip(self._prev_rates, bos_rates, res_bar.rates,
repeat(prev_dt), repeat(dt)))
time1, inter_conc = timed_deplete(
self.chain, bos_conc, inputs, dt, matrix_func=leqi_f3)
time2, inter_conc = timed_deplete(
self.chain, inter_conc, inputs, dt, matrix_func=leqi_f4)
proc_time += time1 + time2
return proc_time, [eos_conc, inter_conc], [res_bar]

38
openmc/deplete/operator.py
View file

@ -16,6 +16,7 @@ import xml.etree.ElementTree as ET
import h5py
import numpy as np
from uncertainties import ufloat
import openmc
import openmc.capi
@ -56,7 +57,7 @@ class Operator(TransportOperator):
Instances of this class can be used to perform depletion using OpenMC as the
transport operator. Normally, a user needn't call methods of this class
directly. Instead, an instance of this class is passed to an integrator
function, such as :func:`openmc.deplete.integrator.cecm`.
class, such as :class:`openmc.deplete.CECMIntegrator`.
Parameters
----------
@ -113,15 +114,16 @@ class Operator(TransportOperator):
Initial heavy metal inventory
local_mats : list of str
All burnable material IDs being managed by a single process
prev_res : ResultsList
Results from a previous depletion calculation
prev_res : ResultsList or None
Results from a previous depletion calculation. ``None`` if no
results are to be used.
diff_burnable_mats : bool
Whether to differentiate burnable materials with multiple instances
"""
def __init__(self, geometry, settings, chain_file=None, prev_results=None,
diff_burnable_mats=False, fission_q=None,
dilute_initial=1.0e3):
super().__init__(chain_file, fission_q, dilute_initial)
super().__init__(chain_file, fission_q, dilute_initial, prev_results)
self.round_number = False
self.prev_res = None
self.settings = settings
@ -141,7 +143,7 @@ class Operator(TransportOperator):
self._mat_index_map = {
lm: self.burnable_mats.index(lm) for lm in self.local_mats}
if prev_results is not None:
if self.prev_res is not None:
# Reload volumes into geometry
prev_results[-1].transfer_volumes(geometry)
@ -176,7 +178,7 @@ class Operator(TransportOperator):
self.reaction_rates.n_nuc, self.reaction_rates.n_react)
self._energy_helper = ChainFissionHelper()
def __call__(self, vec, power, print_out=True):
def __call__(self, vec, power):
"""Runs a simulation.
Parameters
@ -185,8 +187,6 @@ class Operator(TransportOperator):
Total atoms to be used in function.
power : float
Power of the reactor in [W]
print_out : bool, optional
Whether or not to print out time.
Returns
-------
@ -216,15 +216,21 @@ class Operator(TransportOperator):
# Extract results
op_result = self._unpack_tallies_and_normalize(power)
if comm.rank == 0:
time_unpack = time.time()
if print_out:
print("Time to openmc: ", time_openmc - time_start)
print("Time to unpack: ", time_unpack - time_openmc)
return copy.deepcopy(op_result)
@staticmethod
def write_bos_data(step):
"""Write a state-point file with beginning of step data
Parameters
----------
step : int
Current depletion step including restarts
"""
openmc.capi.statepoint_write(
"openmc_simulation_n{}.h5".format(step),
write_source=False)
def _differentiate_burnable_mats(self):
"""Assign distribmats for each burnable material
@ -534,7 +540,7 @@ class Operator(TransportOperator):
rates.fill(0.0)
# Get k and uncertainty
k_combined = openmc.capi.keff()
k_combined = ufloat(*openmc.capi.keff())
# Extract tally bins
nuclides = self._rate_helper.nuclides

17
openmc/deplete/results.py
View file

@ -21,8 +21,8 @@ class Results(object):
Attributes
----------
k : list of float
Eigenvalue for each substep.
k : list of (float, float)
Eigenvalue and uncertainty for each substep.
time : list of float
Time at beginning, end of step, in seconds.
power : float
@ -456,16 +456,7 @@ class Results(object):
# Get indexing terms
vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info()
# For a restart calculation, limit number of stages saved to meet the
# format of the hdf5 file
stages = len(x)
offset = 0
if op.prev_res is not None and op.prev_res[0].n_stages < stages:
offset = stages - op.prev_res[0].n_stages
stages = min(stages, op.prev_res[0].n_stages)
warn("Number of restart integrator stages saved limited by initial"
" depletion integrator choice to {}"
.format(op.prev_res[0].n_stages))
# Create results
results = Results()
@ -475,9 +466,9 @@ class Results(object):
for i in range(stages):
for mat_i in range(n_mat):
results[i, mat_i, :] = x[offset + i][mat_i][:]
results[i, mat_i, :] = x[i][mat_i]
results.k = [r.k for r in op_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

7
tests/dummy_operator.py
View file

@ -1,5 +1,7 @@
import numpy as np
import scipy.sparse as sp
from uncertainties import ufloat
from openmc.deplete.reaction_rates import ReactionRates
from openmc.deplete.abc import TransportOperator, OperatorResult
@ -48,7 +50,7 @@ class DummyOperator(TransportOperator):
reaction_rates[0, 1, 0] = vec[0][1]
# Create a fake rates object
return OperatorResult(0.0, reaction_rates)
return OperatorResult(ufloat(0.0, 0.0), reaction_rates)
@property
def chain(self):
@ -114,6 +116,9 @@ class DummyOperator(TransportOperator):
"""Maps cell name to index in global geometry."""
return self.local_mats
@staticmethod
def write_bos_data(_step):
"""Empty method but avoids calls to C API"""
@property
def reaction_rates(self):

19
tests/regression_tests/deplete/test.py
View file

@ -39,7 +39,7 @@ def test_full(run_in_tmpdir):
space = openmc.stats.Box(lower_left, upper_right)
settings.source = openmc.Source(space=space)
settings.seed = 1
settings.verbosity = 3
settings.verbosity = 1
# Create operator
chain_file = Path(__file__).parents[2] / 'chain_simple.xml'
@ -54,7 +54,7 @@ def test_full(run_in_tmpdir):
power = 2.337e15*4*JOULE_PER_EV*1e6 # MeV/second cm from CASMO
# Perform simulation using the predictor algorithm
openmc.deplete.integrator.predictor(op, dt, power)
openmc.deplete.PredictorIntegrator(op, dt, power).integrate()
# Get path to test and reference results
path_test = op.output_dir / 'depletion_results.h5'
@ -101,3 +101,18 @@ def test_full(run_in_tmpdir):
assert correct, "Discrepancy in mat {} and nuc {}\n{}\n{}".format(
mat, nuc, y_old, y_test)
# Compare statepoint files with depletion results
t_test, k_test = res_test.get_eigenvalue()
t_ref, k_ref = res_ref.get_eigenvalue()
k_state = np.empty_like(k_ref)
# Get statepoint files for all BOS points and EOL
for n in range(N + 1):
statepoint = openmc.StatePoint("openmc_simulation_n{}.h5".format(n))
k_n = statepoint.k_combined
k_state[n] = [k_n.nominal_value, k_n.std_dev]
# Look for exact match pulling from statepoint and depletion_results
assert np.all(k_state == k_test)
assert np.allclose(k_test, k_ref)

7
tests/unit_tests/test_deplete_cecm.py
View file

@ -4,7 +4,7 @@ These tests integrate a simple test problem described in dummy_geometry.py.
"""
from pytest import approx
import openmc.deplete
from openmc.deplete import CECMIntegrator, ResultsList
from tests import dummy_operator
@ -18,10 +18,11 @@ def test_cecm(run_in_tmpdir):
# Perform simulation using the MCNPX/MCNP6 algorithm
dt = [0.75, 0.75]
power = 1.0
openmc.deplete.cecm(op, dt, power, print_out=False)
integrator = CECMIntegrator(op, dt, power)
integrator.integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
res = ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
_, y1 = res.get_atoms("1", "1")
_, y2 = res.get_atoms("1", "2")

6
tests/unit_tests/test_deplete_celi.py
View file

@ -4,7 +4,7 @@ These tests integrate a simple test problem described in dummy_geometry.py.
"""
from pytest import approx
import openmc.deplete
from openmc.deplete import ResultsList, CELIIntegrator
from tests import dummy_operator
@ -18,10 +18,10 @@ def test_celi(run_in_tmpdir):
# Perform simulation using the celi algorithm
dt = [0.75, 0.75]
power = 1.0
openmc.deplete.celi(op, dt, power, print_out=False)
CELIIntegrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
res = ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
_, y1 = res.get_atoms("1", "1")
_, y2 = res.get_atoms("1", "2")

6
tests/unit_tests/test_deplete_cf4.py
View file

@ -4,7 +4,7 @@ These tests integrate a simple test problem described in dummy_geometry.py.
"""
from pytest import approx
import openmc.deplete
from openmc.deplete import CF4Integrator, ResultsList
from tests import dummy_operator
@ -18,10 +18,10 @@ def test_cf4(run_in_tmpdir):
# Perform simulation using the cf4 algorithm
dt = [0.75, 0.75]
power = 1.0
openmc.deplete.cf4(op, dt, power, print_out=False)
CF4Integrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
res = ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
_, y1 = res.get_atoms("1", "1")
_, y2 = res.get_atoms("1", "2")

2
tests/unit_tests/test_deplete_cram.py
View file

@ -6,7 +6,7 @@ Compares a few Mathematica matrix exponentials to CRAM16/CRAM48.
from pytest import approx
import numpy as np
import scipy.sparse as sp
from openmc.deplete.integrator import CRAM16, CRAM48
from openmc.deplete.cram import CRAM16, CRAM48
def test_CRAM16():

6
tests/unit_tests/test_deplete_epc_rk4.py
View file

@ -4,7 +4,7 @@ These tests integrate a simple test problem described in dummy_geometry.py.
"""
from pytest import approx
import openmc.deplete
from openmc.deplete import EPCRK4Integrator, ResultsList
from tests import dummy_operator
@ -18,10 +18,10 @@ def test_epc_rk4(run_in_tmpdir):
# Perform simulation using the epc_rk4 algorithm
dt = [0.75, 0.75]
power = 1.0
openmc.deplete.epc_rk4(op, dt, power, print_out=False)
EPCRK4Integrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
res = ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
_, y1 = res.get_atoms("1", "1")
_, y2 = res.get_atoms("1", "2")

43
tests/unit_tests/test_deplete_integrator.py
View file

@ -11,8 +11,12 @@ import os
from unittest.mock import MagicMock
import numpy as np
from openmc.deplete import (ReactionRates, Results, ResultsList, comm,
OperatorResult)
from uncertainties import ufloat
import pytest
from openmc.deplete import (
ReactionRates, Results, ResultsList, comm, OperatorResult,
PredictorIntegrator, SICELIIntegrator)
def test_results_save(run_in_tmpdir):
@ -74,8 +78,10 @@ def test_results_save(run_in_tmpdir):
t1 = [0.0, 1.0]
t2 = [1.0, 2.0]
op_result1 = [OperatorResult(k, rates) for k, rates in zip(eigvl1, rate1)]
op_result2 = [OperatorResult(k, rates) for k, rates in zip(eigvl2, rate2)]
op_result1 = [OperatorResult(ufloat(*k), rates)
for k, rates in zip(eigvl1, rate1)]
op_result2 = [OperatorResult(ufloat(*k), rates)
for k, rates in zip(eigvl2, rate2)]
Results.save(op, x1, op_result1, t1, 0, 0)
Results.save(op, x2, op_result2, t2, 0, 1)
@ -95,3 +101,32 @@ def test_results_save(run_in_tmpdir):
np.testing.assert_array_equal(res[1].k, eigvl2)
np.testing.assert_array_equal(res[1].time, t2)
@pytest.mark.parametrize("timesteps", (1, [1]))
def test_bad_integrator_inputs(timesteps):
"""Test failure modes for Integrator inputs"""
op = MagicMock()
op.prev_res = None
op.chain = None
op.heavy_metal = 1.0
# No power nor power density given
with pytest.raises(ValueError, match="Either power or power density"):
PredictorIntegrator(op, timesteps)
# Length of power != length time
with pytest.raises(ValueError, match="number of powers"):
PredictorIntegrator(op, timesteps, power=[1, 2])
# Length of power density != length time
with pytest.raises(ValueError, match="number of powers"):
PredictorIntegrator(op, timesteps, power_density=[1, 2])
# SI integrator with bad steps
with pytest.raises(TypeError, match="n_steps"):
SICELIIntegrator(op, timesteps, [1], n_steps=2.5)
with pytest.raises(ValueError, match="n_steps"):
SICELIIntegrator(op, timesteps, [1], n_steps=0)

6
tests/unit_tests/test_deplete_leqi.py
View file

@ -4,7 +4,7 @@ These tests integrate a simple test problem described in dummy_geometry.py.
"""
from pytest import approx
import openmc.deplete
from openmc.deplete import LEQIIntegrator, ResultsList
from tests import dummy_operator
@ -18,10 +18,10 @@ def test_leqi(run_in_tmpdir):
# Perform simulation using the leqi algorithm
dt = [0.75, 0.75]
power = 1.0
openmc.deplete.leqi(op, dt, power, print_out=False)
LEQIIntegrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
res = ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
_, y1 = res.get_atoms("1", "1")
_, y2 = res.get_atoms("1", "2")

14
tests/unit_tests/test_deplete_operator.py
View file

@ -37,14 +37,20 @@ def bare_xs(run_in_tmpdir):
class BareDepleteOperator(TransportOperator):
"""Very basic class for testing the initialization."""
# declare abstract methods so object can be created
def __call__(self, *args, **kwargs):
@staticmethod
def __call__(*args, **kwargs):
pass
def initial_condition(self):
@staticmethod
def initial_condition():
pass
def get_results_info(self):
@staticmethod
def get_results_info():
pass
@staticmethod
def write_bos_data():
pass

7
tests/unit_tests/test_deplete_predictor.py
View file

@ -4,7 +4,7 @@ These tests integrate a simple test problem described in dummy_geometry.py.
"""
from pytest import approx
import openmc.deplete
from openmc.deplete import PredictorIntegrator, ResultsList
from tests import dummy_operator
@ -18,10 +18,11 @@ def test_predictor(run_in_tmpdir):
# Perform simulation using the predictor algorithm
dt = [0.75, 0.75]
power = 1.0
openmc.deplete.predictor(op, dt, power, print_out=False)
PredictorIntegrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
res = ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
_, y1 = res.get_atoms("1", "1")
_, y2 = res.get_atoms("1", "2")

120
tests/unit_tests/test_deplete_restart.py
View file

@ -4,8 +4,12 @@ These tests run in two steps, a first run then a restart run, a simple test
problem described in dummy_geometry.py.
"""
from pytest import approx
from pytest import approx, raises
import openmc.deplete
from openmc.deplete import (
CECMIntegrator, PredictorIntegrator, CELIIntegrator, LEQIIntegrator,
EPCRK4Integrator, CF4Integrator, SICELIIntegrator, SILEQIIntegrator
)
from tests import dummy_operator
@ -20,7 +24,7 @@ def test_restart_predictor(run_in_tmpdir):
# Perform simulation using the predictor algorithm
dt = [0.75]
power = 1.0
openmc.deplete.predictor(op, dt, power, print_out=False)
PredictorIntegrator(op, dt, power).integrate()
# Load the files
prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -30,7 +34,7 @@ def test_restart_predictor(run_in_tmpdir):
op.output_dir = output_dir
# Perform restarts simulation using the predictor algorithm
openmc.deplete.predictor(op, dt, power, print_out=False)
PredictorIntegrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -59,7 +63,8 @@ def test_restart_cecm(run_in_tmpdir):
# Perform simulation using the MCNPX/MCNP6 algorithm
dt = [0.75]
power = 1.0
openmc.deplete.cecm(op, dt, power, print_out=False)
cecm = CECMIntegrator(op, dt, power)
cecm.integrate()
# Load the files
prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -69,7 +74,8 @@ def test_restart_cecm(run_in_tmpdir):
op.output_dir = output_dir
# Perform restarts simulation using the MCNPX/MCNP6 algorithm
openmc.deplete.cecm(op, dt, power, print_out=False)
cecm_restart = CECMIntegrator(op, dt, power)
cecm_restart.integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -84,13 +90,12 @@ def test_restart_cecm(run_in_tmpdir):
assert y1[1] == approx(s1[0])
assert y2[1] == approx(s1[1])
assert y1[3] == approx(s2[0])
assert y2[3] == approx(s2[1])
assert y1[2] == approx(s2[0])
assert y2[2] == approx(s2[1])
def test_restart_predictor_cecm(run_in_tmpdir):
"""Integral regression test of integrator algorithm using predictor
for the first run then CE/CM for the restart run."""
"""Test to ensure that schemes with different stages are not compatible"""
op = dummy_operator.DummyOperator()
output_dir = "test_restart_predictor_cecm"
@ -99,7 +104,7 @@ def test_restart_predictor_cecm(run_in_tmpdir):
# Perform simulation using the predictor algorithm
dt = [0.75]
power = 1.0
openmc.deplete.predictor(op, dt, power, print_out=False)
PredictorIntegrator(op, dt, power).integrate()
# Load the files
prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -108,24 +113,9 @@ def test_restart_predictor_cecm(run_in_tmpdir):
op = dummy_operator.DummyOperator(prev_res)
op.output_dir = output_dir
# Perform restarts simulation using the MCNPX/MCNP6 algorithm
openmc.deplete.cecm(op, dt, power, print_out=False)
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
_, y1 = res.get_atoms("1", "1")
_, y2 = res.get_atoms("1", "2")
# Test solution
s1 = [2.46847546272295, 0.986431226850467]
s2 = [3.09106948392, 0.607102912398]
assert y1[1] == approx(s1[0])
assert y2[1] == approx(s1[1])
assert y1[2] == approx(s2[0])
assert y2[2] == approx(s2[1])
# check ValueError is raised, indicating previous and current stages
with raises(ValueError, match="incompatible.* 1.*2"):
CECMIntegrator(op, dt, power)
def test_restart_cecm_predictor(run_in_tmpdir):
@ -139,7 +129,8 @@ def test_restart_cecm_predictor(run_in_tmpdir):
# Perform simulation using the MCNPX/MCNP6 algorithm
dt = [0.75]
power = 1.0
openmc.deplete.cecm(op, dt, power, print_out=False)
cecm = CECMIntegrator(op, dt, power)
cecm.integrate()
# Load the files
prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -148,25 +139,9 @@ def test_restart_cecm_predictor(run_in_tmpdir):
op = dummy_operator.DummyOperator(prev_res)
op.output_dir = output_dir
# Perform restarts simulation using the predictor algorithm
openmc.deplete.predictor(op, dt, power, print_out=False)
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
_, y1 = res.get_atoms("1", "1")
_, y2 = res.get_atoms("1", "2")
# Test solution
s1 = [1.86872629872102, 1.395525772416039]
s2 = [3.32776806576, 2.391425905]
assert y1[1] == approx(s1[0])
assert y2[1] == approx(s1[1])
assert y1[2] == approx(s2[0])
assert y2[2] == approx(s2[1])
# check ValueError is raised, indicating previous and current stages
with raises(ValueError, match="incompatible.* 2.*1"):
PredictorIntegrator(op, dt, power)
def test_restart_cf4(run_in_tmpdir):
"""Integral regression test of integrator algorithm using CF4."""
@ -178,7 +153,7 @@ def test_restart_cf4(run_in_tmpdir):
# Perform simulation
dt = [0.75]
power = 1.0
openmc.deplete.cf4(op, dt, power, print_out=False)
CF4Integrator(op, dt, power).integrate()
# Load the files
prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -188,7 +163,7 @@ def test_restart_cf4(run_in_tmpdir):
op.output_dir = output_dir
# Perform restarts simulation
openmc.deplete.cf4(op, dt, power, print_out=False)
CF4Integrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -203,8 +178,8 @@ def test_restart_cf4(run_in_tmpdir):
assert y1[1] == approx(s1[0])
assert y2[1] == approx(s1[1])
assert y1[3] == approx(s2[0])
assert y2[3] == approx(s2[1])
assert y1[2] == approx(s2[0])
assert y2[2] == approx(s2[1])
def test_restart_epc_rk4(run_in_tmpdir):
@ -217,7 +192,7 @@ def test_restart_epc_rk4(run_in_tmpdir):
# Perform simulation
dt = [0.75]
power = 1.0
openmc.deplete.epc_rk4(op, dt, power, print_out=False)
EPCRK4Integrator(op, dt, power).integrate()
# Load the files
prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -227,7 +202,7 @@ def test_restart_epc_rk4(run_in_tmpdir):
op.output_dir = output_dir
# Perform restarts simulation
openmc.deplete.epc_rk4(op, dt, power, print_out=False)
EPCRK4Integrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -242,8 +217,8 @@ def test_restart_epc_rk4(run_in_tmpdir):
assert y1[1] == approx(s1[0])
assert y2[1] == approx(s1[1])
assert y1[3] == approx(s2[0])
assert y2[3] == approx(s2[1])
assert y1[2] == approx(s2[0])
assert y2[2] == approx(s2[1])
def test_restart_celi(run_in_tmpdir):
@ -256,7 +231,7 @@ def test_restart_celi(run_in_tmpdir):
# Perform simulation
dt = [0.75]
power = 1.0
openmc.deplete.celi(op, dt, power, print_out=False)
CELIIntegrator(op, dt, power).integrate()
# Load the files
prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -266,7 +241,7 @@ def test_restart_celi(run_in_tmpdir):
op.output_dir = output_dir
# Perform restarts simulation
openmc.deplete.celi(op, dt, power, print_out=False)
CELIIntegrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -281,8 +256,8 @@ def test_restart_celi(run_in_tmpdir):
assert y1[1] == approx(s1[0])
assert y2[1] == approx(s1[1])
assert y1[3] == approx(s2[0])
assert y2[3] == approx(s2[1])
assert y1[2] == approx(s2[0])
assert y2[2] == approx(s2[1])
def test_restart_leqi(run_in_tmpdir):
@ -295,7 +270,7 @@ def test_restart_leqi(run_in_tmpdir):
# Perform simulation
dt = [0.75]
power = 1.0
openmc.deplete.leqi(op, dt, power, print_out=False)
LEQIIntegrator(op, dt, power).integrate()
# Load the files
prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -305,7 +280,7 @@ def test_restart_leqi(run_in_tmpdir):
op.output_dir = output_dir
# Perform restarts simulation
openmc.deplete.leqi(op, dt, power, print_out=False)
LEQIIntegrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -320,8 +295,8 @@ def test_restart_leqi(run_in_tmpdir):
assert y1[1] == approx(s1[0])
assert y2[1] == approx(s1[1])
assert y1[3] == approx(s2[0])
assert y2[3] == approx(s2[1])
assert y1[2] == approx(s2[0])
assert y2[2] == approx(s2[1])
def test_restart_si_celi(run_in_tmpdir):
"""Integral regression test of integrator algorithm using SI-CELI."""
@ -333,7 +308,7 @@ def test_restart_si_celi(run_in_tmpdir):
# Perform simulation
dt = [0.75]
power = 1.0
openmc.deplete.si_celi(op, dt, power, print_out=False)
SICELIIntegrator(op, dt, power).integrate()
# Load the files
prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -343,7 +318,7 @@ def test_restart_si_celi(run_in_tmpdir):
op.output_dir = output_dir
# Perform restarts simulation
openmc.deplete.si_celi(op, dt, power, print_out=False)
SICELIIntegrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -358,8 +333,8 @@ def test_restart_si_celi(run_in_tmpdir):
assert y1[1] == approx(s1[0])
assert y2[1] == approx(s1[1])
assert y1[3] == approx(s2[0])
assert y2[3] == approx(s2[1])
assert y1[2] == approx(s2[0])
assert y2[2] == approx(s2[1])
def test_restart_si_leqi(run_in_tmpdir):
@ -372,7 +347,8 @@ def test_restart_si_leqi(run_in_tmpdir):
# Perform simulation
dt = [0.75]
power = 1.0
openmc.deplete.si_leqi(op, dt, power, print_out=False)
nstages = 10
SILEQIIntegrator(op, dt, power, nstages).integrate()
# Load the files
prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -382,7 +358,7 @@ def test_restart_si_leqi(run_in_tmpdir):
op.output_dir = output_dir
# Perform restarts simulation
openmc.deplete.si_leqi(op, dt, power, print_out=False)
SILEQIIntegrator(op, dt, power, nstages).integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
@ -397,5 +373,5 @@ def test_restart_si_leqi(run_in_tmpdir):
assert y1[1] == approx(s1[0])
assert y2[1] == approx(s1[1])
assert y1[3] == approx(s2[0])
assert y2[3] == approx(s2[1])
assert y1[2] == approx(s2[0])
assert y2[2] == approx(s2[1])

6
tests/unit_tests/test_deplete_si_celi.py
View file

@ -4,7 +4,7 @@ These tests integrate a simple test problem described in dummy_geometry.py.
"""
from pytest import approx
import openmc.deplete
from openmc.deplete import SICELIIntegrator, ResultsList
from tests import dummy_operator
@ -18,10 +18,10 @@ def test_si_celi(run_in_tmpdir):
# Perform simulation using the si_celi algorithm
dt = [0.75, 0.75]
power = 1.0
openmc.deplete.si_celi(op, dt, power, print_out=False)
SICELIIntegrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
res = ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
_, y1 = res.get_atoms("1", "1")
_, y2 = res.get_atoms("1", "2")

6
tests/unit_tests/test_deplete_si_leqi.py
View file

@ -4,7 +4,7 @@ These tests integrate a simple test problem described in dummy_geometry.py.
"""
from pytest import approx
import openmc.deplete
from openmc.deplete import SILEQIIntegrator, ResultsList
from tests import dummy_operator
@ -18,10 +18,10 @@ def test_si_leqi(run_in_tmpdir):
# Perform simulation using the si_leqi algorithm
dt = [0.75, 0.75]
power = 1.0
openmc.deplete.si_leqi(op, dt, power, print_out=False)
SILEQIIntegrator(op, dt, power, 10).integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
res = ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
_, y1 = res.get_atoms("1", "1")
_, y2 = res.get_atoms("1", "2")

5
tests/unit_tests/test_transfer_volumes.py
View file

@ -4,7 +4,8 @@
"""
from pytest import approx
import openmc.deplete
import openmc
from openmc.deplete import PredictorIntegrator, ResultsList
from tests import dummy_operator
@ -18,7 +19,7 @@ def test_transfer_volumes(run_in_tmpdir):
# Perform simulation using the predictor algorithm
dt = [0.75]
power = 1.0
openmc.deplete.predictor(op, dt, power, print_out=False)
PredictorIntegrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")