diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 7b24326d3..532e06655 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -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 `_. .. autosummary:: :toctree: generated diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py index 926e27257..42024863a 100644 --- a/openmc/deplete/integrator/cecm.py +++ b/openmc/deplete/integrator/cecm.py @@ -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 - `_. This algorithm is mathematically - defined as: + `_. + + "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) diff --git a/openmc/deplete/integrator/celi.py b/openmc/deplete/integrator/celi.py index 61edaeb72..0ba86d585 100644 --- a/openmc/deplete/integrator/celi.py +++ b/openmc/deplete/integrator/celi.py @@ -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 `_. - 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: diff --git a/openmc/deplete/integrator/cf4.py b/openmc/deplete/integrator/cf4.py index 79cb3b3d9..e93b22d5c 100644 --- a/openmc/deplete/integrator/cf4.py +++ b/openmc/deplete/integrator/cf4.py @@ -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 + `_. 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, diff --git a/openmc/deplete/integrator/epc_rk4.py b/openmc/deplete/integrator/epc_rk4.py index 98e8ac58b..3789f3698 100644 --- a/openmc/deplete/integrator/epc_rk4.py +++ b/openmc/deplete/integrator/epc_rk4.py @@ -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, diff --git a/openmc/deplete/integrator/leqi.py b/openmc/deplete/integrator/leqi.py index c247bbe63..9d221e877 100644 --- a/openmc/deplete/integrator/leqi.py +++ b/openmc/deplete/integrator/leqi.py @@ -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 `_. - 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: diff --git a/openmc/deplete/integrator/si_celi.py b/openmc/deplete/integrator/si_celi.py index 96aa7e2ba..4f34bace6 100644 --- a/openmc/deplete/integrator/si_celi.py +++ b/openmc/deplete/integrator/si_celi.py @@ -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 `_. - 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 + `_. 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. diff --git a/openmc/deplete/integrator/si_leqi.py b/openmc/deplete/integrator/si_leqi.py index 71c603ed3..05002880a 100644 --- a/openmc/deplete/integrator/si_leqi.py +++ b/openmc/deplete/integrator/si_leqi.py @@ -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 `_. - 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 + `_. 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: diff --git a/openmc/model/model.py b/openmc/model/model.py index dfa3decd3..b6a89792d 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -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.,