diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 60628a61b..aa1288411 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -11,13 +11,6 @@ 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 - :nosignatures: - :template: myfunction.rst - - integrator.si_leqi - .. autosummary:: :toctree: generated :nosignatures: @@ -30,6 +23,7 @@ transport-depletion coupling algorithms `_. integrator.EPC_RK4_Integrator integrator.LEQIIntegrator integrator.SI_CELI_Integrator + integrator.SI_LEQI_Integrator 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/si_celi.py b/openmc/deplete/integrator/si_celi.py index 4d6ce2f08..dde3cfb81 100644 --- a/openmc/deplete/integrator/si_celi.py +++ b/openmc/deplete/integrator/si_celi.py @@ -71,71 +71,3 @@ class SI_CELI_Integrator(SI_Integrator): # end iteration return proc_time, [eos_conc, inter_conc], [res_bar] - - -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. - - Parameters - ---------- - operator : Operator - The operator object to simulate on. - x : list of nuclide vector - Nuclide vector, beginning of time. - op_results : list of OperatorResult - Operator result at BOS. - 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. - m : int, optional - Number of stages. - - Returns - ------- - list of nuclide vector (numpy.array) - Nuclide vector, end of time. - float - Next time - list of OperatorResult - Operator result at end of time. - """ - - chain = operator.chain - - # Deplete to end - proc_time, x_new = timed_deplete( - chain, x[0], op_results[0].rates, dt, print_out) - x.append(x_new) - - for j in range(m + 1): - op_res = operator(x_new, p) - - if j <= 1: - op_res_bar = copy.deepcopy(op_res) - else: - rates = 1/j * op_res.rates + (1 - 1/j) * op_res_bar.rates - k = 1/j * op_res.k + (1 - 1/j) * op_res_bar.k - op_res_bar = OperatorResult(k, rates) - - rates = list(zip(op_results[0].rates, op_res_bar.rates)) - time_1, x_new = timed_deplete( - chain, x[0], rates, dt, print_out, matrix_func=_celi_f1) - time_2, x_new = timed_deplete( - chain, x_new, rates, dt, print_out, matrix_func=_celi_f2) - proc_time += time_1 + time_2 - - # Create results, write to disk - op_results.append(op_res_bar) - Results.save(operator, x, op_results, [t, t+dt], p, i_res+i, proc_time) - - # return updated time and vectors - return [x_new], t + dt, [op_res_bar] diff --git a/openmc/deplete/integrator/si_leqi.py b/openmc/deplete/integrator/si_leqi.py index 86b184038..bd42104ea 100644 --- a/openmc/deplete/integrator/si_leqi.py +++ b/openmc/deplete/integrator/si_leqi.py @@ -6,15 +6,15 @@ from itertools import repeat from uncertainties import ufloat -from .si_celi import si_celi_inner +from .abc import SI_Integrator +from .si_celi import SI_CELI_Integrator from .leqi import _leqi_f1, _leqi_f2, _leqi_f3, _leqi_f4 from .cram import timed_deplete from ..results import Results from ..abc import OperatorResult -def si_leqi(operator, timesteps, power=None, power_density=None, - print_out=True, m=10): +class SI_LEQI_Integrator(SI_Integrator): r"""Deplete using the SI-LE/QI CFQ4 algorithm. Implements the Stochastic Implicit LE/QI Predictor-Corrector algorithm using @@ -22,138 +22,75 @@ def si_leqi(operator, timesteps, power=None, power_density=None, Detailed algorithm can be found in Section 3.2 in `Colin Josey's thesis `_. - - 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. - m : int, optional - Number of stages. """ - 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 + + def __call__(self, bos_conc, bos_rates, dt, power, i): + """Perform the integration across one time step + + Parameters + ---------- + bos_conc : list of numpy.ndarray + Initial concentrations for all nuclides in [atom] for + all depletable materials + bos_rates : list of openmc.deplete.ReactionRates + Reaction rates from operator for all depletable materials + 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: + self._prev_rates = bos_rates + # Perform CELI for initial steps + return SI_CELI_Integrator.__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: - power = [i*operator.heavy_metal for i in power_density] + prevdt = self.timesteps[i - 1] - if not isinstance(power, Iterable): - power = [power]*len(timesteps) + # Perform remaining LE/QI + inputs = list(zip(self._prev_rates, bos_rates, + repeat(prevdt), 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) - # 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) + proc_time += time1 + inter_conc = copy.deepcopy(eos_conc) - # Get the concentrations and reaction rates for the first - # beginning-of-timestep (BOS). Compute with m (stage number) times as - # many neutrons as later simulations for statistics reasons if no - # previous calculation results present - if operator.prev_res is None: - x = [copy.deepcopy(vec)] - if hasattr(operator, "settings"): - operator.settings.particles *= m - op_results = [operator(x[0], power[0])] - if hasattr(operator, "settings"): - operator.settings.particles //= m - else: - # Get initial concentration - x = [operator.prev_res[-1].data[0]] + for j in range(self.n_stages + 1): + inter_res = self.operator(inter_conc, power) - # Get rates - op_results = [operator.prev_res[-1]] - op_results[0].rates = op_results[0].rates[0] + if j <= 1: + res_bar = copy.deepcopy(inter_res) + else: + rates = 1 / j * inter_res.rates + (1 - 1 / j) * res_bar.rates + k = 1 / j * inter_res.k + (1 - 1 / j) * res_bar.k + res_bar = OperatorResult(k, rates) - # Set first stage value of keff - op_results[0].k = ufloat(*op_results[0].k[0]) + inputs = list(zip(self._prev_rates, bos_rates, res_bar.rates, + repeat(prevdt), 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) + proc_time += time1 + time2 - # Scale reaction rates by ratio of powers - power_res = operator.prev_res[-1].power - ratio_power = power[0] / power_res - op_results[0].rates *= ratio_power[0] - - chain = operator.chain - - for i, (dt, p) in enumerate(zip(timesteps, power)): - # LE/QI needs the last step results to start - # Perform SI-CE/LI CFQ4 or restore results for the first step - if i == 0: - dt_l = dt - if i_res <= 1: - op_res_last = copy.deepcopy(op_results[0]) - x, t, op_results = si_celi_inner(operator, x, op_results, 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 = [operator.prev_res[-1].data[0]] - - # Perform remaining LE/QI - inputs = list(zip(op_res_last.rates, op_results[0].rates, - repeat(dt_l), repeat(dt))) - proc_time, x_new = timed_deplete( - chain, x[0], inputs, dt, print_out, matrix_func=_leqi_f1) - time_1, x_new = timed_deplete( - chain, x_new, inputs, dt, print_out, matrix_func=_leqi_f2) - x.append(x_new) - - proc_time += time_1 - - # Loop on inner - for j in range(m + 1): - op_res = operator(x_new, p) - - if j <= 1: - op_res_bar = copy.deepcopy(op_res) - else: - rates = 1/j * op_res.rates + (1 - 1/j) * op_res_bar.rates - k = 1/j * op_res.k + (1 - 1/j) * op_res_bar.k - op_res_bar = OperatorResult(k, rates) - - inputs = list(zip(op_res_last.rates, op_results[0].rates, - op_res_bar.rates, repeat(dt_l), repeat(dt))) - time_1, x_new = timed_deplete( - chain, x[0], inputs, dt, print_out, matrix_func=_leqi_f3) - time_2, x_new = timed_deplete( - chain, x_new, inputs, dt, print_out, matrix_func=_leqi_f4) - - proc_time += time_1 + time_2 - - # Create results, write to disk - op_results.append(op_res_bar) - Results.save( - operator, x, op_results, [t, t+dt], p, i_res+i, proc_time) - - # update results - x = [x_new] - op_res_last = copy.deepcopy(op_results[0]) - op_results = [op_res_bar] - t += dt - dt_l = dt - - # Create results for last point, write to disk - Results.save( - operator, x, op_results, [t, t], p, i_res+len(timesteps)) + return proc_time, [eos_conc, inter_conc], [res_bar] diff --git a/tests/unit_tests/test_deplete_restart.py b/tests/unit_tests/test_deplete_restart.py index 01ee56a82..78a910d41 100644 --- a/tests/unit_tests/test_deplete_restart.py +++ b/tests/unit_tests/test_deplete_restart.py @@ -8,7 +8,7 @@ from pytest import approx import openmc.deplete from openmc.deplete import ( CECMIntegrator, PredictorIntegrator, CELIIntegrator, LEQIIntegrator, - EPC_RK4_Integrator, CF4Integrator, SI_CELI_Integrator + EPC_RK4_Integrator, CF4Integrator, SI_CELI_Integrator, SI_LEQI_Integrator ) from tests import dummy_operator @@ -380,7 +380,8 @@ def test_restart_si_leqi(run_in_tmpdir): # Perform simulation dt = [0.75] power = 1.0 - openmc.deplete.si_leqi(op, dt, power, print_out=False) + nstages = 10 + SI_LEQI_Integrator(op, dt, power, nstages).integrate() # Load the files prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") @@ -390,7 +391,7 @@ def test_restart_si_leqi(run_in_tmpdir): op.output_dir = output_dir # Perform restarts simulation - openmc.deplete.si_leqi(op, dt, power, print_out=False) + SI_LEQI_Integrator(op, dt, power, nstages).integrate() # Load the files res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") diff --git a/tests/unit_tests/test_deplete_si_leqi.py b/tests/unit_tests/test_deplete_si_leqi.py index 1396065c0..c26c06d03 100644 --- a/tests/unit_tests/test_deplete_si_leqi.py +++ b/tests/unit_tests/test_deplete_si_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 SI_LEQI_Integrator, ResultsList from tests import dummy_operator @@ -18,10 +18,10 @@ def test_si_leqi(run_in_tmpdir): # Perform simulation using the si_leqi algorithm dt = [0.75, 0.75] power = 1.0 - openmc.deplete.si_leqi(op, dt, power, print_out=False) + SI_LEQI_Integrator(op, dt, power, 10).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")