diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 452361e73..99c32044d 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -16,8 +16,6 @@ transport-depletion coupling algorithms `_. :nosignatures: :template: myfunction.rst - integrator.celi - integrator.leqi integrator.cf4 integrator.epc_rk4 integrator.si_celi @@ -30,6 +28,8 @@ transport-depletion coupling algorithms `_. 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: diff --git a/openmc/deplete/integrator/celi.py b/openmc/deplete/integrator/celi.py index 2e80c8b73..9346347a3 100644 --- a/openmc/deplete/integrator/celi.py +++ b/openmc/deplete/integrator/celi.py @@ -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]) diff --git a/openmc/deplete/integrator/leqi.py b/openmc/deplete/integrator/leqi.py index 19b1025e5..4b1c072b5 100644 --- a/openmc/deplete/integrator/leqi.py +++ b/openmc/deplete/integrator/leqi.py @@ -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 `_. + + "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 `_. - - "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)) diff --git a/tests/unit_tests/test_deplete_celi.py b/tests/unit_tests/test_deplete_celi.py index 47d3319fd..570ebe151 100644 --- a/tests/unit_tests/test_deplete_celi.py +++ b/tests/unit_tests/test_deplete_celi.py @@ -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") diff --git a/tests/unit_tests/test_deplete_leqi.py b/tests/unit_tests/test_deplete_leqi.py index a26329cf4..0ab4350f6 100644 --- a/tests/unit_tests/test_deplete_leqi.py +++ b/tests/unit_tests/test_deplete_leqi.py @@ -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") diff --git a/tests/unit_tests/test_deplete_restart.py b/tests/unit_tests/test_deplete_restart.py index cf7eeb0ed..a914c8649 100644 --- a/tests/unit_tests/test_deplete_restart.py +++ b/tests/unit_tests/test_deplete_restart.py @@ -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")