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

Merge pull request #1560 from drewejohnson/pick-cram

Browse source 
Allow the depletion solver / CRAM order to be configured
This commit is contained in:
Paul Romano 2020-05-15 10:48:30 -05:00 • committed by GitHub
commit 2813eaddda
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
7 changed files with 489 additions and 143 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -143,8 +143,7 @@ with :func:`cram.CRAM48` being the default.
cram.CRAM16
cram.CRAM48
cram.deplete
cram.timed_deplete
pool.deplete
The following classes are used to help the :class:`openmc.deplete.Operator`
compute quantities like effective fission yields, reaction rates, and

145
openmc/deplete/abc.py
View file

@ -4,15 +4,16 @@ This module contains the Operator class, which is then passed to an integrator
to run a full depletion simulation.
"""
from collections import namedtuple
from collections import defaultdict
from collections.abc import Iterable
from collections import namedtuple, defaultdict
from collections.abc import Iterable, Callable
import os
from pathlib import Path
from abc import ABC, abstractmethod
from copy import deepcopy
from warnings import warn
from numbers import Real, Integral
from inspect import signature
import time
from numpy import nonzero, empty, asarray
from uncertainties import ufloat
@ -23,6 +24,7 @@ from openmc.checkvalue import check_type, check_greater_than
from .results import Results
from .chain import Chain
from .results_list import ResultsList
from .pool import deplete
__all__ = [
@ -595,7 +597,7 @@ class TalliedFissionYieldHelper(FissionYieldHelper):
class Integrator(ABC):
"""Abstract class for solving the time-integration for depletion
r"""Abstract class for solving the time-integration for depletion
Parameters
----------
@ -623,6 +625,17 @@ class Integrator(ABC):
seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
that the values are given in burnup (MW-d of energy deposited per
kilogram of initial heavy metal).
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
----------
@ -634,10 +647,27 @@ class Integrator(ABC):
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
def __init__(self, operator, timesteps, power=None, power_density=None,
timestep_units='s'):
timestep_units='s', solver="cram48"):
# Check number of stages previously used
if operator.prev_res is not None:
res = operator.prev_res[-1]
@ -677,24 +707,24 @@ class Integrator(ABC):
# Determine number of seconds for each timestep
seconds = []
for time, unit, watts in zip(times, units, power):
for timestep, unit, watts in zip(times, units, power):
# Make sure values passed make sense
check_type('timestep', time, Real)
check_greater_than('timestep', time, 0.0, False)
check_type('timestep', timestep, Real)
check_greater_than('timestep', timestep, 0.0, False)
check_type('timestep units', unit, str)
check_type('power', watts, Real)
check_greater_than('power', watts, 0.0, True)
if unit in ('s', 'sec'):
seconds.append(time)
seconds.append(timestep)
elif unit in ('min', 'minute'):
seconds.append(time*_SECONDS_PER_MINUTE)
seconds.append(timestep*_SECONDS_PER_MINUTE)
elif unit in ('h', 'hr', 'hour'):
seconds.append(time*_SECONDS_PER_HOUR)
seconds.append(timestep*_SECONDS_PER_HOUR)
elif unit in ('d', 'day'):
seconds.append(time*_SECONDS_PER_DAY)
seconds.append(timestep*_SECONDS_PER_DAY)
elif unit.lower() == 'mwd/kg':
watt_days_per_kg = 1e6*time
watt_days_per_kg = 1e6*timestep
kilograms = 1e-3*operator.heavy_metal
days = watt_days_per_kg * kilograms / watts
seconds.append(days*_SECONDS_PER_DAY)
@ -704,6 +734,59 @@ class Integrator(ABC):
self.timesteps = asarray(seconds)
self.power = asarray(power)
if isinstance(solver, str):
# Delay importing of cram module, which requires this file
if solver == "cram48":
from .cram import CRAM48
self._solver = CRAM48
elif solver == "cram16":
from .cram import CRAM16
self._solver = CRAM16
else:
raise ValueError(
"Solver {} not understood. Expected 'cram48' or "
"'cram16'".format(solver))
else:
self.solver = solver
@property
def solver(self):
return self._solver
@solver.setter
def solver(self, func):
if not isinstance(func, Callable):
raise TypeError(
"Solver must be callable, not {}".format(type(func)))
try:
sig = signature(func)
except ValueError:
# Guard against callables that aren't introspectable, e.g.
# fortran functions wrapped by F2PY
warn("Could not determine arguments to {}. Proceeding "
"anyways".format(func))
self._solver = func
return
# Inspect arguments
if len(sig.parameters) != 3:
raise ValueError("Function {} does not support three arguments: "
"{!s}".format(func, sig))
for ix, param in enumerate(sig.parameters.values()):
if param.kind in {param.KEYWORD_ONLY, param.VAR_KEYWORD}:
raise ValueError(
"Keyword arguments like {} at position {} are not "
"allowed".format(ix, param))
self._solver = func
def _timed_deplete(self, concs, rates, dt, matrix_func=None):
start = time.time()
results = deplete(
self._solver, self.chain, concs, rates, dt, matrix_func)
return time.time() - start, results
@abstractmethod
def __call__(self, conc, rates, dt, power, i):
"""Perform the integration across one time step
@ -807,7 +890,7 @@ class Integrator(ABC):
class SIIntegrator(Integrator):
"""Abstract class for the Stochastic Implicit Euler integrators
r"""Abstract class for the Stochastic Implicit Euler integrators
Does not provide a ``__call__`` method, but scales and resets
the number of particles used in initial transport calculation
@ -841,6 +924,17 @@ class SIIntegrator(Integrator):
n_steps : int, optional
Number of stochastic iterations per depletion interval.
Must be greater than zero. Default : 10
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
----------
@ -854,12 +948,31 @@ class SIIntegrator(Integrator):
Power of the reactor in [W] for each interval in :attr:`timesteps`
n_steps : int
Number of stochastic iterations per depletion interval
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
def __init__(self, operator, timesteps, power=None, power_density=None,
timestep_units='s', n_steps=10):
timestep_units='s', n_steps=10, solver="cram48"):
check_type("n_steps", n_steps, Integral)
check_greater_than("n_steps", n_steps, 0)
super().__init__(operator, timesteps, power, power_density, timestep_units)
super().__init__(
operator, timesteps, power, power_density, timestep_units,
solver=solver)
self.n_steps = n_steps
def _get_bos_data_from_operator(self, step_index, step_power, bos_conc):

74
openmc/deplete/cram.py
View file

@ -4,9 +4,6 @@ Implements two different forms of CRAM for use in openmc.deplete.
"""
import numbers
from itertools import repeat
from multiprocessing import Pool
import time
import numpy as np
import scipy.sparse as sp
@ -15,76 +12,7 @@ import scipy.sparse.linalg as sla
from openmc.checkvalue import check_type, check_length
from .abc import DepSystemSolver
__all__ = [
"deplete", "timed_deplete", "CRAM16", "CRAM48",
"Cram16Solver", "Cram48Solver", "IPFCramSolver"]
def deplete(chain, x, rates, dt, matrix_func=None):
"""Deplete materials using given reaction rates for a specified time
Parameters
----------
chain : openmc.deplete.Chain
Depletion chain
x : list of numpy.ndarray
Atom number vectors for each material
rates : openmc.deplete.ReactionRates
Reaction rates (from transport operator)
dt : float
Time in [s] to deplete for
maxtrix_func : Callable, optional
Function to form the depletion matrix after calling
``matrix_func(chain, rates, fission_yields)``, where
``fission_yields = {parent: {product: yield_frac}}``
Expected to return the depletion matrix required by
:func:`CRAM48`.
Returns
-------
x_result : list of numpy.ndarray
Updated atom number vectors for each material
"""
fission_yields = chain.fission_yields
if len(fission_yields) == 1:
fission_yields = repeat(fission_yields[0])
elif len(fission_yields) != len(x):
raise ValueError(
"Number of material fission yield distributions {} is not equal "
"to the number of compositions {}".format(len(fission_yields),
len(x)))
if matrix_func is None:
matrices = map(chain.form_matrix, rates, fission_yields)
else:
matrices = map(matrix_func, repeat(chain), rates, fission_yields)
# Use multiprocessing pool to distribute work
with Pool() as pool:
inputs = zip(matrices, x, repeat(dt))
x_result = list(pool.starmap(CRAM48, inputs))
return x_result
def timed_deplete(*args, **kwargs):
"""Wrapper over :func:`deplete` that also returns process time
All arguments and keyword arguments are passed onto
:func:`deplete` directly.
Returns
-------
proc_time: float
Process time required to return from deplete
results: list of numpy arrays
Output from :func:`deplete` call
"""
start = time.time()
results = deplete(*args, **kwargs)
return time.time() - start, results
__all__ = ["CRAM16", "CRAM48", "Cram16Solver", "Cram48Solver", "IPFCramSolver"]
class IPFCramSolver(DepSystemSolver):

315
openmc/deplete/integrators.py
View file

@ -2,7 +2,6 @@ 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
@ -54,6 +53,17 @@ class PredictorIntegrator(Integrator):
that the values are given in burnup (MW-d of energy deposited per
kilogram of initial heavy metal).
.. versionadded:: 0.12
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
@ -66,6 +76,23 @@ class PredictorIntegrator(Integrator):
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
_num_stages = 1
@ -96,7 +123,7 @@ class PredictorIntegrator(Integrator):
operator with predictor
"""
proc_time, conc_end = timed_deplete(self.chain, conc, rates, dt)
proc_time, conc_end = self._timed_deplete(conc, rates, dt)
return proc_time, [conc_end], []
@ -145,6 +172,17 @@ class CECMIntegrator(Integrator):
that the values are given in burnup (MW-d of energy deposited per
kilogram of initial heavy metal).
.. versionadded:: 0.12
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
@ -157,6 +195,23 @@ class CECMIntegrator(Integrator):
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
_num_stages = 2
@ -187,12 +242,12 @@ class CECMIntegrator(Integrator):
Eigenvalue and reaction rates from transport simulations
"""
# deplete across first half of inteval
time0, x_middle = timed_deplete(self.chain, conc, rates, dt / 2)
time0, x_middle = self._timed_deplete(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)
time1, x_end = self._timed_deplete(conc, res_middle.rates, dt)
return time0 + time1, [x_middle, x_end], [res_middle]
@ -244,6 +299,17 @@ class CF4Integrator(Integrator):
that the values are given in burnup (MW-d of energy deposited per
kilogram of initial heavy metal).
.. versionadded:: 0.12
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
@ -256,6 +322,23 @@ class CF4Integrator(Integrator):
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
_num_stages = 4
@ -287,27 +370,27 @@ class CF4Integrator(Integrator):
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)
time1, conc_eos1 = self._timed_deplete(
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)
time2, conc_eos2 = self._timed_deplete(
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)
time3, conc_eos3 = self._timed_deplete(
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)
time4, conc_inter = self._timed_deplete(
bos_conc, list_rates, dt, matrix_func=cf4_f3)
time5, conc_eos5 = self._timed_deplete(
conc_inter, list_rates, dt, matrix_func=cf4_f4)
return (time1 + time2 + time3 + time4 + time5,
[conc_eos1, conc_eos2, conc_eos3, conc_eos5],
@ -360,6 +443,17 @@ class CELIIntegrator(Integrator):
that the values are given in burnup (MW-d of energy deposited per
kilogram of initial heavy metal).
.. versionadded:: 0.12
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
@ -372,6 +466,23 @@ class CELIIntegrator(Integrator):
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
_num_stages = 2
@ -403,17 +514,17 @@ class CELIIntegrator(Integrator):
simulation
"""
# deplete to end using BOS rates
proc_time, conc_ce = timed_deplete(self.chain, bos_conc, rates, dt)
proc_time, conc_ce = self._timed_deplete(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_le1, conc_inter = self._timed_deplete(
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)
time_le2, conc_end = self._timed_deplete(
conc_inter, list_rates, dt, matrix_func=celi_f2)
return proc_time + time_le1 + time_le1, [conc_ce, conc_end], [res_ce]
@ -463,6 +574,17 @@ class EPCRK4Integrator(Integrator):
that the values are given in burnup (MW-d of energy deposited per
kilogram of initial heavy metal).
.. versionadded:: 0.12
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
@ -475,6 +597,18 @@ class EPCRK4Integrator(Integrator):
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
"""
_num_stages = 4
@ -507,24 +641,23 @@ class EPCRK4Integrator(Integrator):
"""
# Step 1: deplete with matrix A(y0) / 2
time1, conc1 = timed_deplete(
self.chain, conc, rates, dt, matrix_func=rk4_f1)
time1, conc1 = self._timed_deplete(
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)
time2, conc2 = self._timed_deplete(
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)
time3, conc3 = self._timed_deplete(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)
time4, conc4 = self._timed_deplete(
conc, list_rates, dt, matrix_func=rk4_f4)
return (time1 + time2 + time3 + time4, [conc1, conc2, conc3, conc4],
[res1, res2, res3])
@ -586,6 +719,17 @@ class LEQIIntegrator(Integrator):
that the values are given in burnup (MW-d of energy deposited per
kilogram of initial heavy metal).
.. versionadded:: 0.12
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
@ -598,6 +742,23 @@ class LEQIIntegrator(Integrator):
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
_num_stages = 2
@ -645,10 +806,10 @@ class LEQIIntegrator(Integrator):
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)
time1, conc_inter = self._timed_deplete(
bos_conc, le_inputs, dt, matrix_func=leqi_f1)
time2, conc_eos0 = self._timed_deplete(
conc_inter, le_inputs, dt, matrix_func=leqi_f2)
res_inter = self.operator(conc_eos0, power)
@ -656,10 +817,10 @@ class LEQIIntegrator(Integrator):
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)
time3, conc_inter = self._timed_deplete(
bos_conc, qi_inputs, dt, matrix_func=leqi_f3)
time4, conc_eos1 = self._timed_deplete(
conc_inter, qi_inputs, dt, matrix_func=leqi_f4)
# store updated rates
self._prev_rates = copy.deepcopy(bos_res.rates)
@ -710,6 +871,17 @@ class SICELIIntegrator(SIIntegrator):
n_steps : int, optional
Number of stochastic iterations per depletion interval.
Must be greater than zero. Default : 10
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
----------
@ -723,6 +895,23 @@ class SICELIIntegrator(SIIntegrator):
Power of the reactor in [W] for each interval in :attr:`timesteps`
n_steps : int
Number of stochastic iterations per depletion interval
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
_num_stages = 2
@ -753,8 +942,7 @@ class SICELIIntegrator(SIIntegrator):
Eigenvalue and reaction rates from intermediate transport
simulations
"""
proc_time, eos_conc = timed_deplete(
self.chain, bos_conc, bos_rates, dt)
proc_time, eos_conc = self._timed_deplete(bos_conc, bos_rates, dt)
inter_conc = copy.deepcopy(eos_conc)
# Begin iteration
@ -769,10 +957,10 @@ class SICELIIntegrator(SIIntegrator):
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)
time1, inter_conc = self._timed_deplete(
bos_conc, list_rates, dt, matrix_func=celi_f1)
time2, inter_conc = self._timed_deplete(
inter_conc, list_rates, dt, matrix_func=celi_f2)
proc_time += time1 + time2
# end iteration
@ -820,6 +1008,17 @@ class SILEQIIntegrator(SIIntegrator):
n_steps : int, optional
Number of stochastic iterations per depletion interval.
Must be greater than zero. Default : 10
solver : str or callable, optional
If a string, must be the name of the solver responsible for
solving the Bateman equations. Current options are:
* ``cram16`` - 16th order IPF CRAM
* ``cram48`` - 48th order IPF CRAM [default]
If a function or other callable, must adhere to the requirements in
:attr:`solver`.
.. versionadded:: 0.12
Attributes
----------
@ -833,6 +1032,24 @@ class SILEQIIntegrator(SIIntegrator):
Power of the reactor in [W] for each interval in :attr:`timesteps`
n_steps : int
Number of stochastic iterations per depletion interval
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
size :math:`t_i`. Can be configured using the ``solver`` argument.
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
* ``t`` is a float of the time step size in seconds, and
* ``n1`` is a :class:`numpy.ndarray` of compositions at the
next time step. Expected to be of the same shape as ``n0``
.. versionadded:: 0.12
"""
_num_stages = 2
@ -879,10 +1096,10 @@ class SILEQIIntegrator(SIIntegrator):
# 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, inter_conc = self._timed_deplete(
bos_conc, inputs, dt, matrix_func=leqi_f1)
time1, eos_conc = self._timed_deplete(
inter_conc, inputs, dt, matrix_func=leqi_f2)
proc_time += time1
inter_conc = copy.deepcopy(eos_conc)
@ -899,10 +1116,10 @@ class SILEQIIntegrator(SIIntegrator):
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)
time1, inter_conc = self._timed_deplete(
bos_conc, inputs, dt, matrix_func=leqi_f3)
time2, inter_conc = self._timed_deplete(
inter_conc, inputs, dt, matrix_func=leqi_f4)
proc_time += time1 + time2
return proc_time, [eos_conc, inter_conc], [res_bar]

58
openmc/deplete/pool.py Normal file
View file

@ -0,0 +1,58 @@
"""Dedicated module containing depletion function
Provided to avoid some circular imports
"""
from itertools import repeat
from multiprocessing import Pool
def deplete(func, chain, x, rates, dt, matrix_func=None):
"""Deplete materials using given reaction rates for a specified time
Parameters
----------
func : callable
Function to use to get new compositions. Expected to have the
signature ``func(A, n0, t) -> n1``
chain : openmc.deplete.Chain
Depletion chain
x : list of numpy.ndarray
Atom number vectors for each material
rates : openmc.deplete.ReactionRates
Reaction rates (from transport operator)
dt : float
Time in [s] to deplete for
maxtrix_func : callable, optional
Function to form the depletion matrix after calling
``matrix_func(chain, rates, fission_yields)``, where
``fission_yields = {parent: {product: yield_frac}}``
Expected to return the depletion matrix required by
``func``
Returns
-------
x_result : list of numpy.ndarray
Updated atom number vectors for each material
"""
fission_yields = chain.fission_yields
if len(fission_yields) == 1:
fission_yields = repeat(fission_yields[0])
elif len(fission_yields) != len(x):
raise ValueError(
"Number of material fission yield distributions {} is not equal "
"to the number of compositions {}".format(
len(fission_yields), len(x)))
if matrix_func is None:
matrices = map(chain.form_matrix, rates, fission_yields)
else:
matrices = map(matrix_func, repeat(chain), rates, fission_yields)
# Use multiprocessing pool to distribute work
with Pool() as pool:
inputs = zip(matrices, x, repeat(dt))
x_result = list(pool.starmap(func, inputs))
return x_result

4
tests/unit_tests/test_deplete_chain.py
View file

@ -7,7 +7,7 @@ import os
from pathlib import Path
import numpy as np
from openmc.deplete import comm, Chain, reaction_rates, nuclide, cram
from openmc.deplete import comm, Chain, reaction_rates, nuclide, cram, pool
import pytest
from tests import cdtemp
@ -414,7 +414,7 @@ def test_fission_yield_attribute(simple_chain):
dummy_conc = [[1, 2]] * (len(empty_chain.fission_yields) + 1)
with pytest.raises(
ValueError, match="fission yield.*not equal.*compositions"):
cram.deplete(empty_chain, dummy_conc, None, 0.5)
pool.deplete(cram.CRAM48, empty_chain, dummy_conc, None, 0.5)
def test_validate(simple_chain):

33
tests/unit_tests/test_deplete_integrator.py
View file

@ -17,7 +17,8 @@ import pytest
from openmc.deplete import (
ReactionRates, Results, ResultsList, comm, OperatorResult,
PredictorIntegrator, CECMIntegrator, CF4Integrator, CELIIntegrator,
EPCRK4Integrator, LEQIIntegrator, SICELIIntegrator, SILEQIIntegrator)
EPCRK4Integrator, LEQIIntegrator, SICELIIntegrator, SILEQIIntegrator,
cram)
from tests import dummy_operator
@ -147,6 +148,23 @@ def test_bad_integrator_inputs():
with pytest.raises(ValueError, match="n_steps"):
SICELIIntegrator(op, timesteps, [1], n_steps=0)
with pytest.raises(ValueError, match="Solver failure"):
PredictorIntegrator(op, timesteps, power=1, solver="failure")
with pytest.raises(TypeError, match=".*callable.*NoneType"):
PredictorIntegrator(op, timesteps, power=1, solver=None)
with pytest.raises(ValueError, match=".*arguments"):
PredictorIntegrator(op, timesteps, power=1, solver=mock_bad_solver_nargs)
def mock_good_solver(A, n, t):
pass
def mock_bad_solver_nargs(A, n):
pass
@pytest.mark.parametrize("scheme", dummy_operator.SCHEMES)
def test_integrator(run_in_tmpdir, scheme):
@ -174,6 +192,19 @@ def test_integrator(run_in_tmpdir, scheme):
assert dep_time.shape == (2, )
assert all(dep_time > 0)
integrator = bundle.solver(operator, [0.75], 1, solver=cram.CRAM48)
assert integrator.solver is cram.CRAM48
integrator = bundle.solver(operator, [0.75], 1, solver="cram16")
assert integrator.solver is cram.CRAM16
integrator.solver = mock_good_solver
assert integrator.solver is mock_good_solver
lfunc = lambda A, n, t: mock_good_solver(A, n, t)
integrator.solver = lfunc
assert integrator.solver is lfunc
@pytest.mark.parametrize("integrator", INTEGRATORS)
def test_timesteps(integrator):