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

address review comments from @paulromano

Browse source
This commit is contained in:
liangjg 2019-01-31 12:01:55 -05:00
parent 2f3a1d5361
commit 760baa43d6
9 changed files with 106 additions and 144 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -7,7 +7,9 @@
.. module:: openmc.deplete
Several functions are provided that implement different time-integration
algorithms for depletion calculations.
algorithms for depletion calculations, which are described in detail in Colin
Josey's thesis, `Development and analysis of high order neutron
transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
.. autosummary::
:toctree: generated

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

@ -11,8 +11,10 @@ def cecm(operator, timesteps, power=None, power_density=None, print_out=True):
r"""Deplete using the CE/CM algorithm.
Implements the second order `CE/CM predictor-corrector algorithm
<https://doi.org/10.13182/NSE14-92>`_. This algorithm is mathematically
defined as:
<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::
y' &= A(y, t) y(t)

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

@ -5,35 +5,38 @@ from collections.abc import Iterable
from .cram import deplete
from ..results import Results
from ..abc import OperatorResult
# 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 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])
return 1/12 * chain.form_matrix(rates[0]) + \
5/12 * chain.form_matrix(rates[1])
def celi(operator, timesteps, power=None, power_density=None,
print_out=True):
r"""Deplete using the CE/LI CFQ4 algorithm.
Implements the CE/LI Predictor-Corrector algorithm using the [fourth order
commutator-free integrator]_.
Implements the CE/LI Predictor-Corrector algorithm using the `fourth order
commutator-free integrator <https://doi.org/10.1137/05063042>`_.
The CE/LI algorithm is mathematically defined as:
"CE/LI" stands for constant extrapolation on predictor and linear
interpolation on corrector. This algorithm is mathematically defined as:
.. math:
y' = A(y, t) y(t)
A_p = A(y_n, t_n)
y_p = expm(A_p h) y_n
A_c = A(y_p, t_n)
A(t) = t/dt * A_c + (dt - t)/dt * A_p
.. math::
y' &= A(y, t) y(t)
Here, A(t) is integrated using the fourth order algorithm CFQ4.
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
Parameters
----------
@ -53,13 +56,6 @@ def celi(operator, timesteps, power=None, power_density=None,
heavy metal inventory to get total power if `power` is not speficied.
print_out : bool, optional
Whether or not to print out time.
References
----------
.. [fourth order commutator-free integrator]
Thalhammer, Mechthild. "A fourth-order commutator-free exponential
integrator for nonautonomous differential equations." SIAM journal on
numerical analysis 44.2 (2006): 851-864.
"""
if power is None:
if power_density is None:

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

@ -9,32 +9,33 @@ from ..results import Results
# Functions to form the special matrix for depletion
def _cf4_f1(chain, rates):
return 1/2 * chain.form_matrix(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])
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])
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])
return -1/12 * chain.form_matrix(rates[0]) + \
1/6 * chain.form_matrix(rates[1]) + \
1/6 * chain.form_matrix(rates[2]) + \
1/4 * chain.form_matrix(rates[3])
def cf4(operator, timesteps, power=None, power_density=None, print_out=True):
r"""Deplete using the CF4 algorithm.
Implements the fourth order [commutator-free Lie algorithm]_.
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::
F_1 = h A(y_0)
F_1 &= h A(y_0)
y_1 &= \text{expm}(1/2 F_1) y_0
@ -69,13 +70,6 @@ def cf4(operator, timesteps, power=None, power_density=None, print_out=True):
heavy metal inventory to get total power if `power` is not speficied.
print_out : bool, optional
Whether or not to print out time.
References
----------
.. [commutator-free Lie algorithm]
Celledoni, E., Marthinsen, A. and Owren, B., 2003. Commutator-free Lie
group methods. Future Generation Computer Systems, 19(3), pp.341-352.
"""
if power is None:
if power_density is None:
@ -129,20 +123,20 @@ def cf4(operator, timesteps, power=None, power_density=None, print_out=True):
x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out,
matrix_func=_cf4_f1)
x.append(x_new)
op_results.append(operator(x[1], p))
op_results.append(operator(x_new, p))
# Step 2: deplete with matrix 1/2*A(y1)
x_new = deplete(chain, x[0], op_results[1].rates, dt, print_out,
matrix_func=_cf4_f1)
x.append(x_new)
op_results.append(operator(x[2], p))
op_results.append(operator(x_new, p))
# Step 3: deplete with matrix -1/2*A(y0)+A(y2)
rates = list(zip(op_results[0].rates, op_results[2].rates))
x_new = deplete(chain, x[1], rates, dt, print_out,
matrix_func=_cf4_f2)
x.append(x_new)
op_results.append(operator(x[3], p))
op_results.append(operator(x_new, p))
# Step 4: deplete with two matrix exponentials
rates = list(zip(op_results[0].rates, op_results[1].rates,

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

@ -9,13 +9,13 @@ from ..results import Results
# Functions to form the special matrix for depletion
def _rk4_f1(chain, rates):
return 1/2 * chain.form_matrix(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])
return 1/6 * chain.form_matrix(rates[0]) + \
1/3 * chain.form_matrix(rates[1]) + \
1/3 * chain.form_matrix(rates[2]) + \
1/6 * chain.form_matrix(rates[3])
def epc_rk4(operator, timesteps, power=None, power_density=None, print_out=True):
r"""Deplete using the EPC-RK4 algorithm.
@ -25,7 +25,7 @@ def epc_rk4(operator, timesteps, power=None, power_density=None, print_out=True)
This algorithm is mathematically defined as:
.. math::
F_1 = h A(y_0)
F_1 &= h A(y_0)
y_1 &= \text{expm}(1/2 F_1) y_0
@ -85,12 +85,6 @@ def epc_rk4(operator, timesteps, power=None, power_density=None, print_out=True)
chain = operator.chain
# Initialize starting index for saving results
if operator.prev_res is None:
i_res = 0
else:
i_res = len(operator.prev_res)
for i, (dt, p) in enumerate(zip(timesteps, power)):
# Get beginning-of-timestep concentrations and reaction rates
# Avoid doing first transport run if already done in previous
@ -132,7 +126,7 @@ def epc_rk4(operator, timesteps, power=None, power_density=None, print_out=True)
x.append(x_new)
op_results.append(operator(x[3], p))
# Step 4: deplete with two matrix exponentials
# Step 4: deplete with matrix 1/6*A(y0)+1/3*A(y1)+1/3*A(y2)+1/6*A(y3)
rates = list(zip(op_results[0].rates, op_results[1].rates,
op_results[2].rates, op_results[3].rates))
x_end = deplete(chain, x[0], rates, dt, print_out,

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

@ -7,58 +7,72 @@ from itertools import repeat
from .celi import celi_inner
from .cram import deplete
from ..results import Results
from ..abc import OperatorResult
# 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
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
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)
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)
f1 = chain.form_matrix(inputs[0])
f2 = chain.form_matrix(inputs[1])
f3 = chain.form_matrix(inputs[2])
dt_l, dt = inputs[3], inputs[4]
return -dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + \
(dt**2 + 2*dt*dt_l + dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + \
(4 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f3
def leqi(operator, timesteps, power=None, power_density=None, print_out=True):
r"""Deplete using the LE/QI CFQ4 algorithm.
Implements the LE/QI Predictor-Corrector algorithm using the [fourth order
commutator-free integrator]_.
Implements the LE/QI Predictor-Corrector algorithm using the `fourth order
commutator-free integrator <https://doi.org/10.1137/05063042>`_.
The LE/QI algorithm is mathematically defined as:
"LE/QI" stands for linear extrapolation on predictor and quadratic
interpolation on corrector. This algorithm is mathematically defined as:
.. math:
y' = A(y, t) y(t)
A_m1 = A(y_n-1, t_n-1)
A_0 = A(y_n, t_n)
A_l(t) linear extrapolation of A_m1, A_0
Integrate to t_n+1 to get y_p
A_c = A(y_p, y_n+1)
A_q(t) quadratic interpolation of A_m1, A_0, A_c
.. math::
y' &= A(y, t) y(t)
Here, A(t) is integrated using the fourth order algorithm CFQ4.
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
It is initialized using the CE/LI algorithm.
@ -80,13 +94,6 @@ def leqi(operator, timesteps, power=None, power_density=None, print_out=True):
heavy metal inventory to get total power if `power` is not speficied.
print_out : bool, optional
Whether or not to print out time.
References
----------
.. [fourth order commutator-free integrator]
Thalhammer, Mechthild. "A fourth-order commutator-free exponential
integrator for nonautonomous differential equations." SIAM journal on
numerical analysis 44.2 (2006): 851-864.
"""
if power is None:
if power_density is None:

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

@ -14,18 +14,10 @@ def si_celi(operator, timesteps, power=None, power_density=None,
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]_.
the `fourth order commutator-free integrator <https://doi.org/10.1137/05063042>`_.
The CE/LI algorithm is mathematically defined as:
.. math:
y' = A(y, t) y(t)
A_p = A(y_n, t_n)
y_p = expm(A_p h) y_n
A_c = A(y_p, t_n)
A(t) = t/dt * A_c + (dt - t)/dt * A_p
Here, A(t) is integrated using the fourth order algorithm CFQ4.
Detailed algorithm can be found in Section 3.2 in `Colin Josey's thesis
<http://hdl.handle.net/1721.1/113721>`_.
Parameters
----------
@ -47,13 +39,6 @@ def si_celi(operator, timesteps, power=None, power_density=None,
Whether or not to print out time.
m : int, optional
Number of stages.
References
----------
.. [fourth order commutator-free integrator]
Thalhammer, Mechthild. "A fourth-order commutator-free exponential
integrator for nonautonomous differential equations." SIAM journal on
numerical analysis 44.2 (2006): 851-864.
"""
if power is None:
if power_density is None:
@ -111,6 +96,7 @@ def si_celi(operator, timesteps, power=None, power_density=None,
# Create results for last point, write to disk
Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps))
def si_celi_inner(operator, x, op_results, p, i, i_res, t, dt, print_out, m=10):
""" The inner loop of SI-CE/LI CFQ4.

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

@ -16,22 +16,10 @@ def si_leqi(operator, timesteps, power=None, power_density=None,
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]_.
the `fourth order commutator-free integrator <https://doi.org/10.1137/05063042>`_.
The LE/QI algorithm is mathematically defined as:
.. math:
y' = A(y, t) y(t)
A_m1 = A(y_n-1, t_n-1)
A_0 = A(y_n, t_n)
A_l(t) linear extrapolation of A_m1, A_0
Integrate to t_n+1 to get y_p
A_c = A(y_p, y_n+1)
A_q(t) quadratic interpolation of A_m1, A_0, A_c
Here, A(t) is integrated using the fourth order algorithm CFQ4.
It is initialized using the CE/LI algorithm.
Detailed algorithm can be found in Section 3.2 in `Colin Josey's thesis
<http://hdl.handle.net/1721.1/113721>`_.
Parameters
----------
@ -53,13 +41,6 @@ def si_leqi(operator, timesteps, power=None, power_density=None,
Whether or not to print out time.
m : int, optional
Number of stages.
References
----------
.. [fourth order commutator-free integrator]
Thalhammer, Mechthild. "A fourth-order commutator-free exponential
integrator for nonautonomous differential equations." SIAM journal on
numerical analysis 44.2 (2006): 851-864.
"""
if power is None:
if power_density is None:

2
openmc/model/model.py
View file

@ -134,7 +134,7 @@ class Model(object):
chain_file : str, optional
Path to the depletion chain XML file. Defaults to the
:envvar:`OPENMC_DEPLETE_CHAIN` environment variable if it exists.
method : string
method : str
Integration method used for depletion (e.g., 'cecm', 'predictor')
**kwargs
Keyword arguments passed to integration function (e.g.,