diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index feb2f6dd6..60628a61b 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -16,7 +16,6 @@ transport-depletion coupling algorithms `_. :nosignatures: :template: myfunction.rst - integrator.si_celi integrator.si_leqi .. autosummary:: @@ -30,6 +29,7 @@ transport-depletion coupling algorithms `_. integrator.CELIIntegrator integrator.EPC_RK4_Integrator integrator.LEQIIntegrator + integrator.SI_CELI_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/abc.py b/openmc/deplete/integrator/abc.py index 07151d13b..422311964 100644 --- a/openmc/deplete/integrator/abc.py +++ b/openmc/deplete/integrator/abc.py @@ -149,3 +149,87 @@ class Integrator(ABC): Results.save( self.operator, conc_list, results_list, time_list, power, index, proc_time) + + +class SI_Integrator(Integrator): + """Abstract for the Stochastic Implicit Euler integrators + + Does not provide a ``__call__`` method, but scales and resets + the number of particles used in initial transport calculation + """ + def __init__(self, operator, timesteps, power=None, power_density=None, + n_steps=10): + """ + 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. + n_stages : int, optional + Number of stages. Default : 10 + """ + + def __init__(self, operator, timesteps, power=None, power_density=None, + n_stages=10): + super().__init__(operator, timesteps, power, power_density) + self.n_stages = n_stages + + def _get_bos_data_from_openmc(self, step_index, step_power, bos_conc): + reset_particles = False + if step_index == 0 and hasattr(self.operator, "settings"): + reset_particles = True + self.operator.settings.particles *= self.n_stages + inherited = super()._get_bos_data_from_openmc( + step_index, step_power, bos_conc) + if reset_particles: + self.operator.settings.particles //= self.n_stages + return inherited + + def integrate(self): + """Perform the entire depletion process across all steps""" + with self.operator as conc: + t, self._ires = self._get_start_data() + + for i, (dt, p) in enumerate(self): + if i == 0: + if self.operator.prev_res is None: + conc, res = self._get_bos_data_from_openmc(i, p, conc) + else: + conc, res = self._get_bos_data_from_restart(i, p, conc) + else: + # Pull rates, k from previous iteration w/o + # re-running transport + res = res_list[-1] # defined in previous i iteration + + proc_time, conc_list, res_list = self(conc, res.rates, dt, p, i) + + # Insert BOS concentration, transport results + conc_list.insert(0, conc) + res_list.insert(0, res) + + # Remove actual EOS concentration for next step + conc = conc_list.pop() + + self._save_results( + conc_list, res_list, [t, t + dt], p, self._ires + i, + proc_time) + + t += dt + + # No final simulation for SIE, use last iteration results + self._save_results( + [conc], [res_list[-1]], [t, t], p, self._ires + len(self)) diff --git a/openmc/deplete/integrator/si_celi.py b/openmc/deplete/integrator/si_celi.py index 872863ef2..4d6ce2f08 100644 --- a/openmc/deplete/integrator/si_celi.py +++ b/openmc/deplete/integrator/si_celi.py @@ -1,103 +1,76 @@ """The SI-CE/LI CFQ4 integrator.""" import copy -from collections.abc import Iterable from uncertainties import ufloat +from .abc import SI_Integrator from .cram import timed_deplete -from ..results import Results -from ..abc import OperatorResult from .celi import _celi_f1, _celi_f2 +from ..abc import OperatorResult +from ..results import Results -def si_celi(operator, timesteps, power=None, power_density=None, - print_out=True, m=10): - r"""Deplete using the SI-CE/LI CFQ4 algorithm. +class SI_CELI_Integrator(SI_Integrator): + r"""Deplete using the si-ce/li cfq4 algorithm. - Implements the Stochastic Implicit CE/LI Predictor-Corrector algorithm using + Implements the stochastic implicit ce/li predictor-corrector algorithm using the `fourth order commutator-free integrator `_. - Detailed algorithm can be found in Section 3.2 in `Colin Josey's thesis + 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 - else: - power = [i*operator.heavy_metal for i in power_density] + def __call__(self, bos_conc, bos_rates, dt, power, _i): + """Perform the integration across one time step - if not isinstance(power, Iterable): - power = [power]*len(timesteps) + Parameters + ---------- + bos_conc : numpy.ndarray + Initial bos_concentrations for all nuclides in [atom] + bos_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. Unused - # 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) + Returns + ------- + proc_time : float + Time spent in CRAM routines for all materials + bos_conc_list : list of numpy.ndarray + Concentrations at each of the intermediate points with + the final bos_concentration as the last element + op_results : list of openmc.deplete.OperatorResult + Eigenvalue and reaction rates from intermediate transport + simulations + """ + proc_time, eos_conc = timed_deplete( + self.chain, bos_conc, bos_rates, dt) + 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]] + # Begin iteration + 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]) + list_rates = list(zip(bos_rates, res_bar.rates)) + time1, inter_conc = timed_deplete( + self.chain, bos_conc, list_rates, dt, matrix_func=_celi_f1) + time2, inter_conc = timed_deplete( + self.chain, inter_conc, list_rates, dt, matrix_func=_celi_f2) + 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] - - for i, (dt, p) in enumerate(zip(timesteps, power)): - x, t, op_results = si_celi_inner(operator, x, op_results, p, - i, i_res, t, dt, print_out, m) - - # Create results for last point, write to disk - Results.save( - operator, x, op_results, [t, t], p, i_res + len(timesteps)) + # 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): diff --git a/tests/unit_tests/test_deplete_restart.py b/tests/unit_tests/test_deplete_restart.py index 617575720..01ee56a82 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 + EPC_RK4_Integrator, CF4Integrator, SI_CELI_Integrator ) from tests import dummy_operator @@ -341,7 +341,7 @@ def test_restart_si_celi(run_in_tmpdir): # Perform simulation dt = [0.75] power = 1.0 - openmc.deplete.si_celi(op, dt, power, print_out=False) + SI_CELI_Integrator(op, dt, power).integrate() # Load the files prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") @@ -351,7 +351,7 @@ def test_restart_si_celi(run_in_tmpdir): op.output_dir = output_dir # Perform restarts simulation - openmc.deplete.si_celi(op, dt, power, print_out=False) + SI_CELI_Integrator(op, dt, power).integrate() # Load the files res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") diff --git a/tests/unit_tests/test_deplete_si_celi.py b/tests/unit_tests/test_deplete_si_celi.py index 71679439c..954c9d890 100644 --- a/tests/unit_tests/test_deplete_si_celi.py +++ b/tests/unit_tests/test_deplete_si_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 SI_CELI_Integrator, ResultsList from tests import dummy_operator @@ -18,10 +18,10 @@ def test_si_celi(run_in_tmpdir): # Perform simulation using the si_celi algorithm dt = [0.75, 0.75] power = 1.0 - openmc.deplete.si_celi(op, dt, power, print_out=False) + SI_CELI_Integrator(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")