Remove celi, leqi functions for CELIIntegrator, LEQIIntegrator

The depletion functions openmc.deplete.celi and
openmc.deplete.leqi hvae been removed in favor of a class-based
approach. The following syntax will replicate the behavior
of the integration schemes:

>>> from openmc.deplete import CELIIntegrator
>>> celi = CELIIntegrator(operator, time, power)
>>> celi.integrate()

The expression can be reduced to a single line:

>>> LEQIIntegrator(operator, time, power).integrate()
This commit is contained in:
Andrew Johnson 2019-07-17 16:46:30 -05:00
parent 1d791a16c9
commit 769f903dbc
No known key found for this signature in database
GPG key ID: 253418E91B7F6FEB
6 changed files with 160 additions and 276 deletions

View file

@ -16,8 +16,6 @@ transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
:nosignatures:
:template: myfunction.rst
integrator.celi
integrator.leqi
integrator.cf4
integrator.epc_rk4
integrator.si_celi
@ -30,6 +28,8 @@ transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
integrator.PredictorIntegrator
integrator.CECMIntegrator
integrator.CELIIntegrator
integrator.LEQIIntegrator
Each of these functions expects a "transport operator" to be passed. An operator
specific to OpenMC is available using the following class:

View file

@ -1,25 +1,10 @@
"""The CE/LI CFQ4 integrator."""
import copy
from collections.abc import Iterable
from .cram import timed_deplete
from ..results import Results
from .abc import Integrator
# Functions to form the special matrix for depletion
def _celi_f1(chain, rates):
return 5/12 * chain.form_matrix(rates[0]) + \
1/12 * chain.form_matrix(rates[1])
def _celi_f2(chain, rates):
return 1/12 * chain.form_matrix(rates[0]) + \
5/12 * chain.form_matrix(rates[1])
def celi(operator, timesteps, power=None, power_density=None,
print_out=True):
class CELIIntegrator(Integrator):
r"""Deplete using the CE/LI CFQ4 algorithm.
Implements the CE/LI Predictor-Corrector algorithm using the `fourth order
@ -37,131 +22,57 @@ def celi(operator, timesteps, power=None, power_density=None,
y_{n+1} &= \text{expm}(\frac{h}{12} A_0 + \frac{5h}{12} A1)
\text{expm}(\frac{5h}{12} A_0 + \frac{h}{12} A1) y_n
\end{aligned}
Parameters
----------
operator : openmc.deplete.TransportOperator
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not cumulative.
power : float or iterable of float, optional
Power of the reactor in [W]. A single value indicates that the power is
constant over all timesteps. An iterable indicates potentially different
power levels for each timestep. For a 2D problem, the power can be given
in [W/cm] as long as the "volume" assigned to a depletion material is
actually an area in [cm^2]. Either `power` or `power_density` must be
specified.
power_density : float or iterable of float, optional
Power density of the reactor in [W/gHM]. It is multiplied by initial
heavy metal inventory to get total power if `power` is not speficied.
print_out : bool, optional
Whether or not to print out time.
"""
if power is None:
if power_density is None:
raise ValueError(
"Neither power nor power density was specified.")
if not isinstance(power_density, Iterable):
power = power_density*operator.heavy_metal
else:
power = [i*operator.heavy_metal for i in power_density]
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
# Generate initial conditions
with operator as vec:
# Initialize time and starting index
if operator.prev_res is None:
t = 0.0
i_res = 0
else:
t = operator.prev_res[-1].time[-1]
i_res = len(operator.prev_res)
for i, (dt, p) in enumerate(zip(timesteps, power)):
vec, t, _ = celi_inner(operator, vec, p, i, i_res, t, dt,
print_out)
# Perform one last simulation
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], power[-1])]
# Create results, write to disk
Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps))
def celi_inner(operator, vec, p, i, i_res, t, dt, print_out):
""" The inner loop of CE/LI CFQ4.
Parameters
----------
operator : Operator
The operator object to simulate on.
x : list of nuclide vector
Nuclide vector, beginning of time.
p : float
Power of the reactor in [W]
i : int
Current iteration number.
i_res : int
Starting index, for restart calculation.
t : float
Time at start of step.
dt : float
Time step.
print_out : bool
Whether or not to print out time.
Returns
-------
list of numpy.array
Nuclide vector, end of time.
float
Next time
OperatorResult
Operator result from beginning of step.
"""
chain = operator.chain
def __call__(self, bos_conc, rates, dt, power, _i=-1):
"""Perform the integration across one time step
# Get beginning-of-timestep concentrations and reaction rates
# Avoid doing first transport run if already done in previous
# calculation
if i > 0 or operator.prev_res is None:
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], p)]
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 [W]
_i : int
Current iteration count. Not used
else:
# Get initial concentration
x = [operator.prev_res[-1].data[0]]
Returns
-------
proc_time : float
Time spent in CRAM routines for all materials
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)
# Get rates
op_results = [operator.prev_res[-1]]
op_results[0].rates = op_results[0].rates[0]
# deplete using two matrix exponeitials
list_rates = list(zip(rates, res_ce.rates))
# Set first stage value of keff
op_results[0].k = op_results[0].k[0]
time_le1, conc_inter = timed_deplete(
self.chain, bos_conc, list_rates, dt, matrix_func=_celi_f1)
# Scale reaction rates by ratio of powers
power_res = operator.prev_res[-1].power
ratio_power = p / power_res
op_results[0].rates *= ratio_power[0]
time_le2, conc_end = timed_deplete(
self.chain, conc_inter, list_rates, dt, matrix_func=_celi_f2)
# Deplete to end
proc_time, x_new = timed_deplete(chain, x[0], op_results[0].rates, dt, print_out)
x.append(x_new)
op_results.append(operator(x[1], p))
return proc_time + time_le1 + time_le1, [conc_ce, conc_end], [res_ce]
# Deplete with two matrix exponentials
rates = list(zip(op_results[0].rates, op_results[1].rates))
time_1, x_end = timed_deplete(chain, x[0], rates, dt, print_out,
matrix_func=_celi_f1)
time_2, x_end = timed_deplete(chain, x_end, rates, dt, print_out,
matrix_func=_celi_f2)
# Create results, write to disk
Results.save(operator, x, op_results, [t, t + dt], p, i_res + i, proc_time + time_1 + time_2)
# 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])
# return updated time and vectors
return x_end, t + dt, op_results[0]
def _celi_f2(chain, rates):
return 1/12 * chain.form_matrix(rates[0]) + \
5/12 * chain.form_matrix(rates[1])

View file

@ -1,12 +1,109 @@
"""The LE/QI CFQ4 integrator."""
import copy
from collections.abc import Iterable
from itertools import repeat
from .celi import celi_inner
from .abc import Integrator
from .celi import CELIIntegrator
from .cram import timed_deplete
from ..results import Results
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.
"""
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 [W]
i : int
Current depletion step index
Returns
-------
proc_time : float
Time spent in CRAM routines for all materials
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._ires <= 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]
prevdt = self.timesteps[i] - prev_res.time[0]
self._prev_rates = prev_res.rates[0]
else:
prevdt = 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(prevdt), 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(prevdt), 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
@ -42,129 +139,3 @@ def _leqi_f4(chain, inputs):
return -dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + \
(dt**2 + 2*dt*dt_l + dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + \
(4 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f3
def leqi(operator, timesteps, power=None, power_density=None, print_out=True):
r"""Deplete using the LE/QI CFQ4 algorithm.
Implements the LE/QI Predictor-Corrector algorithm using the `fourth order
commutator-free integrator <https://doi.org/10.1137/05063042>`_.
"LE/QI" stands for linear extrapolation on predictor and quadratic
interpolation on corrector. This algorithm is mathematically defined as:
.. math::
\begin{aligned}
y' &= A(y, t) y(t) \\
A_{last} &= A(y_{n-1}, t_n - h_1) \\
A_0 &= A(y_n, t_n) \\
F_1 &= \frac{-h_2^2}{12h_1} A_{last} + \frac{h_2(6h_1+h_2)}{12h_1} A_0 \\
F_2 &= \frac{-5h_2^2}{12h_1} A_{last} + \frac{h_2(6h_1+5h_2)}{12h_1} A_0 \\
y_p &= \text{expm}(F_2) \text{expm}(F_1) y_n \\
A_1 &= A(y_p, t_n + h_2) \\
F_3 &= \frac{-h_2^3}{12 h_1 (h_1 + h_2)} A_{last} +
\frac{h_2 (5 h_1^2 + 6 h_2 h_1 + h_2^2)}{12 h_1 (h_1 + h_2)} A_0 +
\frac{h_2 h_1)}{12 (h_1 + h_2)} A_1 \\
F_4 &= \frac{-h_2^3}{12 h_1 (h_1 + h_2)} A_{last} +
\frac{h_2 (h_1^2 + 2 h_2 h_1 + h_2^2)}{12 h_1 (h_1 + h_2)} A_0 +
\frac{h_2 (5 h_1^2 + 4 h_2 h_1)}{12 h_1 (h_1 + h_2)} A_1 \\
y_{n+1} &= \text{expm}(F_4) \text{expm}(F_3) y_n
\end{aligned}
It is initialized using the CE/LI algorithm.
Parameters
----------
operator : openmc.deplete.TransportOperator
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not cumulative.
power : float or iterable of float, optional
Power of the reactor in [W]. A single value indicates that the power is
constant over all timesteps. An iterable indicates potentially different
power levels for each timestep. For a 2D problem, the power can be given
in [W/cm] as long as the "volume" assigned to a depletion material is
actually an area in [cm^2]. Either `power` or `power_density` must be
specified.
power_density : float or iterable of float, optional
Power density of the reactor in [W/gHM]. It is multiplied by initial
heavy metal inventory to get total power if `power` is not speficied.
print_out : bool, optional
Whether or not to print out time.
"""
if power is None:
if power_density is None:
raise ValueError(
"Neither power nor power density was specified.")
if not isinstance(power_density, Iterable):
power = power_density*operator.heavy_metal
else:
power = [i*operator.heavy_metal for i in power_density]
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
# Generate initial conditions
with operator as vec:
# Initialize time and starting index
if operator.prev_res is None:
t = 0.0
i_res = 0
else:
t = operator.prev_res[-1].time[-1]
i_res = len(operator.prev_res)
chain = operator.chain
for i, (dt, p) in enumerate(zip(timesteps, power)):
# LE/QI needs the last step results to start
# Perform CE/LI CFQ4 or restore results for the first step
if i == 0:
if i_res <= 1:
dt_l = dt
x_new, t, op_res_last = celi_inner(operator, vec, p, i,
i_res, t, dt, print_out)
continue
else:
dt_l = t - operator.prev_res[-2].time[0]
op_res_last = operator.prev_res[-2]
op_res_last.rates = op_res_last.rates[0]
x_new = operator.prev_res[-1].data[0]
# Perform remaining LE/QI
x = [copy.deepcopy(x_new)]
op_results = [operator(x[0], p)]
inputs = list(zip(op_res_last.rates, op_results[0].rates,
repeat(dt_l), repeat(dt)))
time_1, x_new = timed_deplete(
chain, x[0], inputs, dt, print_out, matrix_func=_leqi_f1)
time_2, x_new = timed_deplete(
chain, x_new, inputs, dt, print_out, matrix_func=_leqi_f2)
x.append(x_new)
op_results.append(operator(x[1], p))
inputs = list(zip(op_res_last.rates, op_results[0].rates,
op_results[1].rates, repeat(dt_l), repeat(dt)))
time_3, x_new = timed_deplete(
chain, x[0], inputs, dt, print_out, matrix_func=_leqi_f3)
time_4, x_new = timed_deplete(
chain, x_new, inputs, dt, print_out, matrix_func=_leqi_f4)
# Create results, write to disk
Results.save(
operator, x, op_results, [t, t+dt], p, i_res+i,
time_1 + time_2 + time_3 + time_4)
# update results
op_res_last = copy.deepcopy(op_results[0])
t += dt
dt_l = dt
# Perform one last simulation
x = [copy.deepcopy(x_new)]
op_results = [operator(x[0], power[-1])]
# Create results, write to disk
Results.save(
operator, x, op_results, [t, t], p, i_res + len(timesteps))

View file

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

View file

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

View file

@ -6,7 +6,9 @@ problem described in dummy_geometry.py.
from pytest import approx
import openmc.deplete
from openmc.deplete import CECMIntegrator, PredictorIntegrator
from openmc.deplete import (
CECMIntegrator, PredictorIntegrator, CELIIntegrator, LEQIIntegrator,
)
from tests import dummy_operator
@ -261,7 +263,7 @@ def test_restart_celi(run_in_tmpdir):
# Perform simulation
dt = [0.75]
power = 1.0
openmc.deplete.celi(op, dt, power, print_out=False)
CELIIntegrator(op, dt, power).integrate()
# Load the files
prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
@ -271,7 +273,7 @@ def test_restart_celi(run_in_tmpdir):
op.output_dir = output_dir
# Perform restarts simulation
openmc.deplete.celi(op, dt, power, print_out=False)
CELIIntegrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
@ -300,7 +302,7 @@ def test_restart_leqi(run_in_tmpdir):
# Perform simulation
dt = [0.75]
power = 1.0
openmc.deplete.leqi(op, dt, power, print_out=False)
LEQIIntegrator(op, dt, power).integrate()
# Load the files
prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
@ -310,7 +312,7 @@ def test_restart_leqi(run_in_tmpdir):
op.output_dir = output_dir
# Perform restarts simulation
openmc.deplete.leqi(op, dt, power, print_out=False)
LEQIIntegrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")