crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-25 12:35:29 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Move openmc.deplete.integrator into openmc.deplete

Browse source 
Much of the previous API is intact, with the major change
being CRAM functions are imported from openmc.deplete.cram
in the test_deplete_cram.

The integrator abstract classes are placed in
openmc/deplete/abc.py with all concrete classes going in to
openmc/deplete/integrators.py

Documentation updated accordingly
This commit is contained in:
Andrew Johnson 2019-08-07 10:51:19 -05:00
parent 5adc3b5816
commit 0bc7800a92
No known key found for this signature in database
GPG key ID: 253418E91B7F6FEB
17 changed files with 1210 additions and 1267 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -16,14 +16,14 @@ transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
:nosignatures:
:template: myintegrator.rst
integrator.PredictorIntegrator
integrator.CECMIntegrator
integrator.CELIIntegrator
integrator.CF4Integrator
integrator.EPCRK4Integrator
integrator.LEQIIntegrator
integrator.SICELIIntegrator
integrator.SILEQIIntegrator
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:
@ -114,5 +114,5 @@ as follows:
:nosignatures:
:template: myfunction.rst
integrator.CRAM16
integrator.CRAM48
cram.CRAM16
cram.CRAM48

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)

287
openmc/deplete/abc.py
View file

@ -5,17 +5,20 @@ 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
@ -44,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
@ -104,7 +107,7 @@ class TransportOperator(ABC):
self.prev_res = None
else:
check_type("previous results", prev_results, ResultsList)
self.prev_results = prev_res
self.prev_results = prev_results
@property
def dilute_initial(self):
@ -179,7 +182,8 @@ 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.
"""
@ -266,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
-------
@ -357,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))

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

@ -11,7 +11,9 @@ 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, matrix_func=None):

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

@ -1,17 +0,0 @@
"""
Integrator
===========
The integrator subcomponents.
"""
from .abc import Integrator
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 *

277
openmc/deplete/integrator/abc.py
View file

@ -1,277 +0,0 @@
from copy import deepcopy
from abc import ABC, abstractmethod
from collections.abc import Iterable
from numbers import Integral
from uncertainties import ufloat
from openmc.checkvalue import check_type, check_greater_than
from openmc.capi import statepoint_write
from openmc.deplete import Results, OperatorResult
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))

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

@ -1,92 +0,0 @@
"""The CE/CM integrator."""
from .abc import Integrator
from .cram import timed_deplete
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=-1):
"""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]

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

@ -1,110 +0,0 @@
"""The CE/LI CFQ4 integrator."""
from .cram import timed_deplete
from .abc import Integrator
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=-1):
"""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
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]
# 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])

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

@ -1,140 +0,0 @@
"""The CF4 integrator."""
import copy
from collections.abc import Iterable
from .cram import timed_deplete
from .abc import Integrator
from ..results import Results
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):
"""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
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
"""
# 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])
# 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])

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

@ -1,119 +0,0 @@
"""The EPC-RK4 integrator."""
from .cram import timed_deplete
from .abc import Integrator
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):
"""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, 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])
# 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])

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

@ -1,173 +0,0 @@
"""The LE/QI CFQ4 integrator."""
import copy
from itertools import repeat
from .abc import Integrator
from .celi import CELIIntegrator
from .cram import timed_deplete
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])
# 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

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

@ -1,81 +0,0 @@
"""First-order predictor algorithm."""
from .abc import Integrator
from .cram import timed_deplete
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], []

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

@ -1,109 +0,0 @@
"""The SI-CE/LI CFQ4 integrator."""
import copy
from .abc import SIIntegrator
from .cram import timed_deplete
from .celi import _celi_f1, _celi_f2
from ..abc import OperatorResult
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):
"""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
Current depletion step index. Unused
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]

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

@ -1,129 +0,0 @@
"""The SI-LE/QI CFQ4 integrator."""
import copy
from itertools import repeat
from .abc import SIIntegrator
from .si_celi import SICELIIntegrator
from .leqi import _leqi_f1, _leqi_f2, _leqi_f3, _leqi_f4
from .cram import timed_deplete
from ..abc import OperatorResult
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]

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=-1):
"""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):
"""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
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
"""
# 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=-1):
"""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
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):
"""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, 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):
"""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
Current depletion step index. Unused
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]

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():