From d28546bcd671d7905d59ff3fa66cbb8ae36a54c5 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 26 Jul 2019 12:05:06 -0500 Subject: [PATCH] Add SI_Integrator, SI_CELI_Integrator classes; remove si_celi The si_celi depletion function has been removed. The SI-CELI method can be implemented with >>> from openmc.deplete import SI_CELI_Integrator >>> SI_CELI_Integrator(op, dt, power, stages).integrate() The stages parameter is optional, and defaults to 10 to be consistent with the previous default for si_celi. The SI_CELI_Integrator inherits from the new abstract base class SI_Integrator. There are a few differences in how the SI-based methods perform the integration. - The initial, t=0.0, i=0, no restart transport solution is scaled by the number of stages. - The SI methods also do not perform a new transport solution at each successive iteration [t>0, i>0]. Instead, the reaction rates are pulled from the last stage of the previous step. - There is no final transport solution once all the depletion steps have been taken. The final reaction rates are saved as those from the last stage of the last depletion step. The si_celi function has been removed from documentation and testing, and replaced with the SI_CELI_Integrator. --- docs/source/pythonapi/deplete.rst | 2 +- openmc/deplete/integrator/abc.py | 84 +++++++++++++++ openmc/deplete/integrator/si_celi.py | 131 +++++++++-------------- tests/unit_tests/test_deplete_restart.py | 6 +- tests/unit_tests/test_deplete_si_celi.py | 6 +- 5 files changed, 143 insertions(+), 86 deletions(-) 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")