diff --git a/docs/source/_templates/myintegrator.rst b/docs/source/_templates/myintegrator.rst new file mode 100644 index 0000000000..803dfde067 --- /dev/null +++ b/docs/source/_templates/myintegrator.rst @@ -0,0 +1,9 @@ +{{ fullname }} +{{ underline }} + +.. currentmodule:: {{ module }} + +.. autoclass:: {{ objname }} + :members: + :inherited-members: + :special-members: __call__, __len__, __iter__ diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 0dcd47ef54..a72504ac92 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -12,18 +12,18 @@ Josey's thesis, `Development and analysis of high order neutron transport-depletion coupling algorithms `_. .. autosummary:: - :toctree: generated - :nosignatures: - :template: myfunction.rst + :toctree: generated + :nosignatures: + :template: myintegrator.rst - integrator.predictor - integrator.cecm - integrator.celi - integrator.leqi - integrator.cf4 - integrator.epc_rk4 - integrator.si_celi - integrator.si_leqi + PredictorIntegrator + CECMIntegrator + CELIIntegrator + CF4Integrator + EPCRK4Integrator + LEQIIntegrator + SICELIIntegrator + SILEQIIntegrator Each of these functions expects a "transport operator" to be passed. An operator specific to OpenMC is available using the following class: @@ -95,7 +95,18 @@ The following classes are abstract classes that can be used to extend the EnergyHelper TransportOperator -Each of the integrator functions also relies on a number of "helper" functions +Custom integrators can be developed by subclassing from the following abstract +base classes: + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myintegrator.rst + + Integrator + SIIntegrator + +Each of the integrator classes also relies on a number of "helper" functions as follows: .. autosummary:: @@ -103,5 +114,5 @@ as follows: :nosignatures: :template: myfunction.rst - integrator.CRAM16 - integrator.CRAM48 + cram.CRAM16 + cram.CRAM48 diff --git a/docs/source/usersguide/geometry.rst b/docs/source/usersguide/geometry.rst index 80e60bd34c..51daf7b5d1 100644 --- a/docs/source/usersguide/geometry.rst +++ b/docs/source/usersguide/geometry.rst @@ -25,7 +25,7 @@ surface is a locus of zeros of a function of Cartesian coordinates Defining a surface alone is not sufficient to specify a volume -- in order to define an actual volume, one must reference the *half-space* of a surface. A -surface half-space is the region whose points satisfy a positive of negative +surface half-space is the region whose points satisfy a positive or negative inequality of the surface equation. For example, for a sphere of radius one centered at the origin, the surface equation is :math:`f(x,y,z) = x^2 + y^2 + z^2 - 1 = 0`. Thus, we say that the negative half-space of the sphere, is diff --git a/openmc/deplete/__init__.py b/openmc/deplete/__init__.py index b5eb2dec2e..8c19615ece 100644 --- a/openmc/deplete/__init__.py +++ b/openmc/deplete/__init__.py @@ -38,4 +38,4 @@ from .reaction_rates import * from .abc import * from .results import * from .results_list import * -from .integrator import * +from .integrators import * diff --git a/openmc/deplete/_matrix_funcs.py b/openmc/deplete/_matrix_funcs.py new file mode 100644 index 0000000000..1cacbd2c06 --- /dev/null +++ b/openmc/deplete/_matrix_funcs.py @@ -0,0 +1,79 @@ +"""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 cf4_f1(chain, 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]) + + +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])) + + +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])) + + +def rk4_f1(chain, 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])) + + +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 + + +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 + + +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) + + +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) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index a33d794441..0f1203aeb8 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -5,18 +5,22 @@ to run a full depletion simulation. """ from collections import namedtuple +from collections.abc import Iterable import os from pathlib import Path from abc import ABC, abstractmethod -from xml.etree import ElementTree as ET +from copy import deepcopy from warnings import warn -from numbers import Real +from numbers import Real, Integral from numpy import nonzero, empty +from uncertainties import ufloat from openmc.data import DataLibrary, JOULE_PER_EV from openmc.checkvalue import check_type, check_greater_than +from .results import Results from .chain import Chain +from .results_list import ResultsList OperatorResult = namedtuple('OperatorResult', ['k', 'rates']) OperatorResult.__doc__ = """\ @@ -24,8 +28,8 @@ Result of applying transport operator Parameters ---------- -k : float - Resulting eigenvalue +k : uncertainties.ufloat + Resulting eigenvalue and standard deviation rates : openmc.deplete.ReactionRates Resulting reaction rates @@ -43,8 +47,8 @@ class TransportOperator(ABC): Each depletion integrator is written to work with a generic transport operator that takes a vector of material compositions and returns an - eigenvalue and reaction rates. This abstract class sets the requirements for - such a transport operator. Users should instantiate + eigenvalue and reaction rates. This abstract class sets the requirements + for such a transport operator. Users should instantiate :class:`openmc.deplete.Operator` rather than this class. Parameters @@ -61,6 +65,8 @@ class TransportOperator(ABC): in initial condition to ensure they exist in the decay chain. Only done for nuclides with reaction rates. Defaults to 1.0e3. + prev_results : ResultsList, optional + Results from a previous depletion calculation. Attributes ---------- @@ -68,8 +74,12 @@ class TransportOperator(ABC): Initial atom density [atoms/cm^3] to add for nuclides that are zero in initial condition to ensure they exist in the decay chain. Only done for nuclides with reaction rates. + prev_res : ResultsList or None + Results from a previous depletion calculation. ``None`` if no + results are to be used. """ - def __init__(self, chain_file=None, fission_q=None, dilute_initial=1.0e3): + def __init__(self, chain_file=None, fission_q=None, dilute_initial=1.0e3, + prev_results=None): self.dilute_initial = dilute_initial self.output_dir = '.' @@ -93,6 +103,11 @@ class TransportOperator(ABC): "of adding depletion_chain to OPENMC_CROSS_SECTIONS", FutureWarning) self.chain = Chain.from_xml(chain_file, fission_q) + if prev_results is None: + self.prev_res = None + else: + check_type("previous results", prev_results, ResultsList) + self.prev_results = prev_results @property def dilute_initial(self): @@ -106,15 +121,15 @@ class TransportOperator(ABC): self._dilute_initial = value @abstractmethod - def __call__(self, vec, print_out=True): + def __call__(self, vec, power): """Runs a simulation. Parameters ---------- vec : list of numpy.ndarray Total atoms to be used in function. - print_out : bool, optional - Whether or not to print out time. + power : float + Power of the reactor in [W] Returns ------- @@ -122,7 +137,6 @@ class TransportOperator(ABC): Eigenvalue and reaction rates resulting from transport operator """ - pass def __enter__(self): # Save current directory and move to specific output directory @@ -157,8 +171,6 @@ class TransportOperator(ABC): Total density for initial conditions. """ - pass - @abstractmethod def get_results_info(self): """Returns volume list, cell lists, and nuc lists. @@ -170,16 +182,28 @@ class TransportOperator(ABC): nuc_list : list of str A list of all nuclide names. Used for sorting the simulation. burn_list : list of int - A list of all cell IDs to be burned. Used for sorting the simulation. + A list of all cell IDs to be burned. Used for sorting the + simulation. full_burn_list : list of int All burnable materials in the geometry. """ - pass - def finalize(self): pass + @abstractmethod + def write_bos_data(self, step): + """Document beginning of step data for a given step + + Called at the beginning of a depletion step and at + the final point in the simulation. + + Parameters + ---------- + step : int + Current depletion step including restarts + """ + class ReactionRateHelper(ABC): """Abstract class for generating reaction rates for operators @@ -246,7 +270,8 @@ class ReactionRateHelper(ABC): Parameters ---------- number : iterable of float - Number density [atoms/b-cm] of each nuclide tracked in the calculation. + Number density [atoms/b-cm] of each nuclide tracked in the + calculation. Returns ------- @@ -337,3 +362,271 @@ class EnergyHelper(ABC): def nuclides(self, nuclides): check_type("nuclides", nuclides, list, str) self._nuclides = nuclides + + +class Integrator(ABC): + """Abstract class for solving the time-integration for depletion + + Parameters + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + 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. + + Attributes + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + chain : openmc.deplete.Chain + Depletion chain + timesteps : iterable of float + Size of each depletion interval in [s] + power : iterable of float + Power of the reactor in [W] for each interval in :attr:`timesteps` + """ + + def __init__(self, operator, timesteps, power=None, power_density=None): + # Check number of stages previously used + if operator.prev_res is not None: + res = operator.prev_res[-1] + if res.data.shape[0] != self._num_stages: + raise ValueError( + "{} incompatible with previous restart calculation. " + "Previous scheme used {} intermediate solutions, while " + "this uses {}".format( + self.__class__.__name__, res.data.shape[0], + self._num_stages)) + self.operator = operator + self.chain = operator.chain + if not isinstance(timesteps, Iterable): + self.timesteps = [timesteps] + else: + self.timesteps = timesteps + if power is None: + if power_density is None: + raise ValueError("Either power or power density must be set") + if not isinstance(power_density, Iterable): + power = power_density * operator.heavy_metal + else: + power = [p * operator.heavy_metal for p in power_density] + + if not isinstance(power, Iterable): + # Ensure that power is single value if that is the case + power = [power] * len(self.timesteps) + elif len(power) != len(self.timesteps): + raise ValueError( + "Number of time steps != number of powers. {} vs {}".format( + len(self.timesteps), len(power))) + + self.power = power + + @abstractmethod + def __call__(self, conc, 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 in [W] + i : int + Current depletion step index + + Returns + ------- + proc_time : float + Time spent in CRAM routines for all materials in [s] + 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 + simulations + """ + + @property + @abstractmethod + def _num_stages(self): + """Number of intermediate transport solutions + + Needed to ensure schemes are consistent with restarts + """ + + def __iter__(self): + """Return pairs of time steps in [s] and powers in [W]""" + return zip(self.timesteps, self.power) + + def __len__(self): + """Return integer number of depletion intervals""" + return len(self.timesteps) + + def _get_bos_data_from_operator(self, step_index, step_power, bos_conc): + """Get beginning of step concentrations, reaction rates from Operator + """ + x = deepcopy(bos_conc) + res = self.operator(x, step_power) + self.operator.write_bos_data(step_index + self._i_res) + return x, res + + def _get_bos_data_from_restart(self, step_index, step_power, bos_conc): + """Get beginning of step concentrations, reaction rates from restart""" + res = self.operator.prev_res[-1] + # Depletion methods expect list of arrays + bos_conc = list(res.data[0]) + rates = res.rates[0] + k = ufloat(res.k[0, 0], res.k[0, 1]) + + # Scale rates by ratio of powers + rates *= step_power / res.power[0] + return bos_conc, OperatorResult(k, rates) + + def _get_start_data(self): + if self.operator.prev_res is None: + return 0.0, 0 + return (self.operator.prev_res[-1].time[-1], + len(self.operator.prev_res) - 1) + + def integrate(self): + """Perform the entire depletion process across all steps""" + with self.operator as conc: + t, self._i_res = self._get_start_data() + + for i, (dt, p) in enumerate(self): + if i > 0 or self.operator.prev_res is None: + conc, res = self._get_bos_data_from_operator(i, p, conc) + else: + conc, res = self._get_bos_data_from_restart(i, p, conc) + 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() + + Results.save(self.operator, conc_list, res_list, [t, t + dt], + p, self._i_res + i, proc_time) + + t += dt + + # Final simulation + res_list = [self.operator(conc, p)] + Results.save(self.operator, [conc], res_list, [t, t], + p, self._i_res + len(self), proc_time) + self.operator.write_bos_data(len(self) + self._i_res) + + +class SIIntegrator(Integrator): + """Abstract class 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 + + 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_steps : int, optional + Number of stochastic iterations per depletion interval. + Must be greater than zero. Default : 10 + + Attributes + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + chain : openmc.deplete.Chain + Depletion chain + timesteps : iterable of float + Size of each depletion interval in [s] + power : iterable of float + Power of the reactor in [W] for each interval in :attr:`timesteps` + n_steps : int + Number of stochastic iterations per depletion interval + """ + def __init__(self, operator, timesteps, power=None, power_density=None, + n_steps=10): + check_type("n_steps", n_steps, Integral) + check_greater_than("n_steps", n_steps, 0) + super().__init__(operator, timesteps, power, power_density) + self.n_steps = n_steps + + def _get_bos_data_from_operator(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_operator( + 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._i_res = 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_operator(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() + + Results.save(self.operator, conc_list, res_list, [t, t + dt], + p, self._i_res + i, proc_time) + + t += dt + + # No final simulation for SIE, use last iteration results + Results.save(self.operator, [conc], [res_list[-1]], [t, t], + p, self._i_res + len(self), proc_time) + self.operator.write_bos_data(self._i_res + len(self)) diff --git a/openmc/deplete/integrator/cram.py b/openmc/deplete/cram.py similarity index 95% rename from openmc/deplete/integrator/cram.py rename to openmc/deplete/cram.py index f5da6d696b..41a4e4c710 100644 --- a/openmc/deplete/integrator/cram.py +++ b/openmc/deplete/cram.py @@ -11,10 +11,12 @@ import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as sla -from .. import comm +from . import comm + +__all__ = ["deplete", "timed_deplete", "CRAM16", "CRAM48"] -def deplete(chain, x, rates, dt, print_out=True, matrix_func=None): +def deplete(chain, x, rates, dt, matrix_func=None): """Deplete materials using given reaction rates for a specified time Parameters @@ -27,29 +29,22 @@ def deplete(chain, x, rates, dt, print_out=True, matrix_func=None): Reaction rates (from transport operator) dt : float Time in [s] to deplete for - print_out : bool, optional - Whether to show elapsed time - maxtrix_func : function, optional - Function to form the depletion matrix + maxtrix_func : Callable, optional + Function of two variables: ``chain`` and ``rates``. + Expected to return the depletion matrix required by + :func:`CRAM48`. Returns ------- x_result : list of numpy.ndarray Updated atom number vectors for each material - """ - t_start = time.time() # Use multiprocessing pool to distribute work with Pool() as pool: iters = zip(repeat(chain), x, rates, repeat(dt), repeat(matrix_func)) x_result = list(pool.starmap(_cram_wrapper, iters)) - t_end = time.time() - if comm.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) - return x_result diff --git a/openmc/deplete/integrator/__init__.py b/openmc/deplete/integrator/__init__.py deleted file mode 100644 index db906b3c13..0000000000 --- a/openmc/deplete/integrator/__init__.py +++ /dev/null @@ -1,16 +0,0 @@ -""" -Integrator -=========== - -The integrator subcomponents. -""" - -from .cf4 import * -from .cecm import * -from .celi import * -from .cram import * -from .epc_rk4 import * -from .leqi import * -from .predictor import * -from .si_celi import * -from .si_leqi import * diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py deleted file mode 100644 index 758302b027..0000000000 --- a/openmc/deplete/integrator/cecm.py +++ /dev/null @@ -1,123 +0,0 @@ -"""The CE/CM integrator.""" - -import copy -from collections.abc import Iterable - -from .cram import timed_deplete -from ..results import Results - - -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 - `_. - - "CE/CM" stands for constant extrapolation on predictor and constant - midpoint on corrector. This algorithm is mathematically defined as: - - .. math:: - \begin{aligned} - y' &= A(y, t) y(t) \\ - A_p &= A(y_n, t_n) \\ - y_m &= \text{expm}(A_p h/2) y_n \\ - A_c &= A(y_m, t_n + h/2) \\ - y_{n+1} &= \text{expm}(A_c h) 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) - - chain = operator.chain - - 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 - # calculation - if i > 0 or operator.prev_res is None: - x = [copy.deepcopy(vec)] - op_results = [operator(x[0], p)] - - else: - # Get initial concentration - x = [operator.prev_res[-1].data[0]] - - # Get rates - op_results = [operator.prev_res[-1]] - op_results[0].rates = op_results[0].rates[0] - - # Set first stage value of keff - op_results[0].k = op_results[0].k[0] - - # 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] - - # Deplete for first half of timestep - proc_time, x_middle = timed_deplete( - chain, x[0], op_results[0].rates, dt/2, print_out) - - # Get middle-of-timestep reaction rates - x.append(x_middle) - op_results.append(operator(x_middle, p)) - - # Deplete for full timestep using beginning-of-step materials - # and middle-of-timestep reaction rates - pt_end, x_end = timed_deplete( - chain, x[0], op_results[1].rates, dt, print_out) - - # Create results, write to disk - Results.save( - operator, x, op_results, [t, t + dt], p, i_res + i, - proc_time + pt_end) - - # Advance time, update vector - t += dt - vec = copy.deepcopy(x_end) - - # 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)) diff --git a/openmc/deplete/integrator/celi.py b/openmc/deplete/integrator/celi.py deleted file mode 100644 index 2e80c8b739..0000000000 --- a/openmc/deplete/integrator/celi.py +++ /dev/null @@ -1,167 +0,0 @@ -"""The CE/LI CFQ4 integrator.""" - -import copy -from collections.abc import Iterable - -from .cram import timed_deplete -from ..results import Results - - -# 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): - r"""Deplete using the CE/LI CFQ4 algorithm. - - Implements the CE/LI Predictor-Corrector algorithm using the `fourth order - commutator-free integrator `_. - - "CE/LI" stands for constant extrapolation on predictor and linear - interpolation on corrector. This algorithm is mathematically defined as: - - .. math:: - \begin{aligned} - y' &= A(y, t) y(t) \\ - 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 - \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 - - # 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)] - - else: - # Get initial concentration - x = [operator.prev_res[-1].data[0]] - - # Get rates - op_results = [operator.prev_res[-1]] - op_results[0].rates = op_results[0].rates[0] - - # Set first stage value of keff - op_results[0].k = op_results[0].k[0] - - # 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] - - # 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)) - - # 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) - - # return updated time and vectors - return x_end, t + dt, op_results[0] diff --git a/openmc/deplete/integrator/cf4.py b/openmc/deplete/integrator/cf4.py deleted file mode 100644 index 8b7f3a05a0..0000000000 --- a/openmc/deplete/integrator/cf4.py +++ /dev/null @@ -1,164 +0,0 @@ -"""The CF4 integrator.""" - -import copy -from collections.abc import Iterable - -from .cram import timed_deplete -from ..results import Results - - -# Functions to form the special matrix for depletion -def _cf4_f1(chain, 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]) - - -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]) - - -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]) - - -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 - `_. - This algorithm is mathematically defined as: - - .. math:: - \begin{aligned} - F_1 &= h A(y_0) \\ - y_1 &= \text{expm}(1/2 F_1) y_0 \\ - F_2 &= h A(y_1) \\ - y_2 &= \text{expm}(1/2 F_2) y_0 \\ - F_3 &= h A(y_2) \\ - y_3 &= \text{expm}(-1/2 F_1 + F_3) y_1 \\ - F_4 &= h A(y_3) \\ - y_4 &= \text{expm}( 1/4 F_1 + 1/6 F_2 + 1/6 F_3 - 1/12 F_4) - \text{expm}(-1/12 F_1 + 1/6 F_2 + 1/6 F_3 + 1/4 F_4) y_0 - \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) - - chain = operator.chain - - 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 - # calculation - if i > 0 or operator.prev_res is None: - x = [copy.deepcopy(vec)] - op_results = [operator(x[0], p)] - - else: - # Get initial concentration - x = [operator.prev_res[-1].data[0]] - - # Get rates - op_results = [operator.prev_res[-1]] - op_results[0].rates = op_results[0].rates[0] - - # Set first stage value of keff - op_results[0].k = op_results[0].k[0] - - # 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] - - # Step 1: deplete with matrix 1/2*A(y0) - time_1, x_new = timed_deplete( - chain, x[0], op_results[0].rates, dt, print_out, - matrix_func=_cf4_f1) - x.append(x_new) - op_results.append(operator(x_new, p)) - - # Step 2: deplete with matrix 1/2*A(y1) - time_2, x_new = timed_deplete( - chain, x[0], op_results[1].rates, dt, print_out, - matrix_func=_cf4_f1) - x.append(x_new) - 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)) - time_3, x_new = timed_deplete( - chain, x[1], rates, dt, print_out, matrix_func=_cf4_f2) - x.append(x_new) - 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, - op_results[2].rates, op_results[3].rates)) - time_4, x_end = timed_deplete( - chain, x[0], rates, dt, print_out, matrix_func=_cf4_f3) - time_5, x_end = timed_deplete( - chain, x_end, rates, dt, print_out, matrix_func=_cf4_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 + time_5) - - # Advance time, update vector - t += dt - vec = copy.deepcopy(x_end) - - # 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)) diff --git a/openmc/deplete/integrator/epc_rk4.py b/openmc/deplete/integrator/epc_rk4.py deleted file mode 100644 index 373c323e42..0000000000 --- a/openmc/deplete/integrator/epc_rk4.py +++ /dev/null @@ -1,148 +0,0 @@ -"""The EPC-RK4 integrator.""" - -import copy -from collections.abc import Iterable - -from .cram import timed_deplete -from ..results import Results - - -# Functions to form the special matrix for depletion -def _rk4_f1(chain, 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]) - -def epc_rk4(operator, timesteps, power=None, power_density=None, print_out=True): - r"""Deplete using the EPC-RK4 algorithm. - - Implements an extended predictor-corrector algorithm with traditional - Runge-Kutta 4 method. - This algorithm is mathematically defined as: - - .. math:: - \begin{aligned} - F_1 &= h A(y_0) \\ - y_1 &= \text{expm}(1/2 F_1) y_0 \\ - F_2 &= h A(y_1) \\ - y_2 &= \text{expm}(1/2 F_2) y_0 \\ - F_3 &= h A(y_2) \\ - y_3 &= \text{expm}(F_3) y_0 \\ - F_4 &= h A(y_3) \\ - y_4 &= \text{expm}(1/6 F_1 + 1/3 F_2 + 1/3 F_3 + 1/6 F_4) y_0 - \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) - - chain = operator.chain - - 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 - # calculation - if i > 0 or operator.prev_res is None: - x = [copy.deepcopy(vec)] - op_results = [operator(x[0], p)] - - else: - # Get initial concentration - x = [operator.prev_res[-1].data[0]] - - # Get rates - op_results = [operator.prev_res[-1]] - op_results[0].rates = op_results[0].rates[0] - - # Set first stage value of keff - op_results[0].k = op_results[0].k[0] - - # 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] - - # Step 1: deplete with matrix 1/2*A(y0) - time_1, x_new = timed_deplete( - chain, x[0], op_results[0].rates, dt, print_out, - matrix_func=_rk4_f1) - x.append(x_new) - op_results.append(operator(x[1], p)) - - # Step 2: deplete with matrix 1/2*A(y1) - time_2, x_new = timed_deplete( - chain, x[0], op_results[1].rates, dt, print_out, - matrix_func=_rk4_f1) - x.append(x_new) - op_results.append(operator(x[2], p)) - - # Step 3: deplete with matrix A(y2) - time_3, x_new = timed_deplete( - chain, x[0], op_results[2].rates, dt, print_out) - x.append(x_new) - op_results.append(operator(x[3], p)) - - # 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)) - time_4, x_end = timed_deplete( - chain, x[0], rates, dt, print_out, matrix_func=_rk4_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) - - # Advance time, update vector - t += dt - vec = copy.deepcopy(x_end) - - # 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)) diff --git a/openmc/deplete/integrator/leqi.py b/openmc/deplete/integrator/leqi.py deleted file mode 100644 index 19b1025e5c..0000000000 --- a/openmc/deplete/integrator/leqi.py +++ /dev/null @@ -1,170 +0,0 @@ -"""The LE/QI CFQ4 integrator.""" - -import copy -from collections.abc import Iterable -from itertools import repeat - -from .celi import celi_inner -from .cram import timed_deplete -from ..results import Results - - -# 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 - - -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 - - -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 - - -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 - - -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/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py deleted file mode 100644 index b30484c695..0000000000 --- a/openmc/deplete/integrator/predictor.py +++ /dev/null @@ -1,106 +0,0 @@ -"""First-order predictor algorithm.""" - -import copy -from collections.abc import Iterable - -from .cram import timed_deplete -from ..results import Results - - -def predictor(operator, timesteps, power=None, power_density=None, - print_out=True): - r"""Deplete using a first-order predictor algorithm. - - Implements the first-order predictor algorithm. This algorithm is - mathematically defined as: - - .. math:: - \begin{aligned} - y' &= A(y, t) y(t) \\ - A_p &= A(y_n, t_n) \\ - y_{n+1} &= \text{expm}(A_p h) 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) - - proc_time = None - - # 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) - 1 - - chain = operator.chain - - 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 - # calculation - if i > 0 or operator.prev_res is None: - x = [copy.deepcopy(vec)] - op_results = [operator(x[0], p)] - - # Create results, write to disk - Results.save(operator, x, op_results, [t, t + dt], p, i_res + i, proc_time) - else: - # Get initial concentration - x = [operator.prev_res[-1].data[0]] - - # Get rates - op_results = [operator.prev_res[-1]] - op_results[0].rates = op_results[0].rates[0] - - # 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] - - # Deplete for full timestep - proc_time, x_end = timed_deplete( - chain, x[0], op_results[0].rates, dt, print_out) - - # Advance time, update vector - t += dt - vec = copy.deepcopy(x_end) - - # 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), proc_time) diff --git a/openmc/deplete/integrator/si_celi.py b/openmc/deplete/integrator/si_celi.py deleted file mode 100644 index 26087c1f2c..0000000000 --- a/openmc/deplete/integrator/si_celi.py +++ /dev/null @@ -1,166 +0,0 @@ -"""The SI-CE/LI CFQ4 integrator.""" - -import copy -from collections.abc import Iterable - -from .cram import timed_deplete -from ..results import Results -from ..abc import OperatorResult -from .celi import _celi_f1, _celi_f2 - - -def si_celi(operator, timesteps, power=None, power_density=None, - print_out=True, m=10): - 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 `_. - - 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] - - 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) - - # 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]] - - # Get rates - op_results = [operator.prev_res[-1]] - op_results[0].rates = op_results[0].rates[0] - - # Set first stage value of keff - op_results[0].k = op_results[0].k[0] - - # 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)) - - -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 deleted file mode 100644 index 7d01782bae..0000000000 --- a/openmc/deplete/integrator/si_leqi.py +++ /dev/null @@ -1,157 +0,0 @@ -"""The SI-LE/QI CFQ4 integrator.""" - -import copy -from collections.abc import Iterable -from itertools import repeat - -from .si_celi import si_celi_inner -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): - 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 `_. - - 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] - - 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) - - # 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]] - - # Get rates - op_results = [operator.prev_res[-1]] - op_results[0].rates = op_results[0].rates[0] - - # Set first stage value of keff - op_results[0].k = op_results[0].k[0] - - # 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)) diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py new file mode 100644 index 0000000000..67106aa3e1 --- /dev/null +++ b/openmc/deplete/integrators.py @@ -0,0 +1,836 @@ +import copy +from itertools import repeat + +from .abc import Integrator, SIIntegrator, OperatorResult +from .cram import timed_deplete +from ._matrix_funcs import ( + cf4_f1, cf4_f2, cf4_f3, cf4_f4, celi_f1, celi_f2, + leqi_f1, leqi_f2, leqi_f3, leqi_f4, rk4_f1, rk4_f4 +) + +__all__ = [ + "PredictorIntegrator", "CECMIntegrator", "CF4Integrator", + "CELIIntegrator", "EPCRK4Integrator", "LEQIIntegrator", + "SICELIIntegrator", "SILEQIIntegrator"] + + +class PredictorIntegrator(Integrator): + r"""Deplete using a first-order predictor algorithm. + + Implements the first-order predictor algorithm. This algorithm is + mathematically defined as: + + .. math:: + \begin{aligned} + y' &= A(y, t) y(t) \\ + A_p &= A(y_n, t_n) \\ + y_{n+1} &= \text{expm}(A_p h) y_n + \end{aligned} + + Parameters + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + 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. + + Attributes + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + chain : openmc.deplete.Chain + Depletion chain + timesteps : iterable of float + Size of each depletion interval in [s] + power : iterable of float + Power of the reactor in [W] for each interval in :attr:`timesteps` + """ + _num_stages = 1 + + def __call__(self, conc, rates, dt, power, _i=None): + """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 in [W] + _i : int or None + Iteration index. Not used + + Returns + ------- + proc_time : float + Time spent in CRAM routines for all materials in [s] + conc_list : list of numpy.ndarray + Concentrations at end of interval + op_results : empty list + Kept for consistency with API. No intermediate calls to + operator with predictor + + """ + proc_time, conc_end = timed_deplete(self.chain, conc, rates, dt) + return proc_time, [conc_end], [] + + +class CECMIntegrator(Integrator): + r"""Deplete using the CE/CM algorithm. + + Implements the second order `CE/CM predictor-corrector algorithm + `_. + + "CE/CM" stands for constant extrapolation on predictor and constant + midpoint on corrector. This algorithm is mathematically defined as: + + .. math:: + \begin{aligned} + y' &= A(y, t) y(t) \\ + A_p &= A(y_n, t_n) \\ + y_m &= \text{expm}(A_p h/2) y_n \\ + A_c &= A(y_m, t_n + h/2) \\ + y_{n+1} &= \text{expm}(A_c h) y_n + \end{aligned} + + Parameters + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + 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. + + Attributes + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + chain : openmc.deplete.Chain + Depletion chain + timesteps : iterable of float + Size of each depletion interval in [s] + power : iterable of float + Power of the reactor in [W] for each interval in :attr:`timesteps` + """ + _num_stages = 2 + + def __call__(self, conc, rates, dt, power, _i=None): + """Integrate using CE/CM + + 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, optional + Current iteration count. Not used + + Returns + ------- + proc_time : float + Time spent in CRAM routines for all materials in [s] + 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 transport simulations + """ + # deplete across first half of inteval + time0, x_middle = timed_deplete(self.chain, conc, rates, dt / 2) + res_middle = self.operator(x_middle, power) + + # deplete across entire interval with BOS concentrations, + # MOS reaction rates + time1, x_end = timed_deplete(self.chain, conc, res_middle.rates, dt) + + return time0 + time1, [x_middle, x_end], [res_middle] + + +class CF4Integrator(Integrator): + r"""Deplete using the CF4 algorithm. + + Implements the fourth order `commutator-free Lie algorithm + `_. + This algorithm is mathematically defined as: + + .. math:: + \begin{aligned} + F_1 &= h A(y_0) \\ + y_1 &= \text{expm}(1/2 F_1) y_0 \\ + F_2 &= h A(y_1) \\ + y_2 &= \text{expm}(1/2 F_2) y_0 \\ + F_3 &= h A(y_2) \\ + y_3 &= \text{expm}(-1/2 F_1 + F_3) y_1 \\ + F_4 &= h A(y_3) \\ + y_4 &= \text{expm}( 1/4 F_1 + 1/6 F_2 + 1/6 F_3 - 1/12 F_4) + \text{expm}(-1/12 F_1 + 1/6 F_2 + 1/6 F_3 + 1/4 F_4) y_0 + \end{aligned} + + Parameters + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + 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. + + Attributes + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + chain : openmc.deplete.Chain + Depletion chain + timesteps : iterable of float + Size of each depletion interval in [s] + power : iterable of float + Power of the reactor in [W] for each interval in :attr:`timesteps` + """ + _num_stages = 4 + + def __call__(self, bos_conc, bos_rates, dt, power, _i=None): + """Perform the integration across one time step + + Parameters + ---------- + bos_conc : numpy.ndarray + Initial 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 in [W] + _i : int, optional + Current depletion step index. Not used + + Returns + ------- + proc_time : float + Time spent in CRAM routines for all materials in [s] + 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 + simulations + """ + # Step 1: deplete with matrix 1/2*A(y0) + time1, conc_eos1 = timed_deplete( + self.chain, bos_conc, bos_rates, dt, matrix_func=cf4_f1) + res1 = self.operator(conc_eos1, power) + + # Step 2: deplete with matrix 1/2*A(y1) + time2, conc_eos2 = timed_deplete( + self.chain, bos_conc, res1.rates, dt, matrix_func=cf4_f1) + res2 = self.operator(conc_eos2, power) + + # Step 3: deplete with matrix -1/2*A(y0)+A(y2) + list_rates = list(zip(bos_rates, res2.rates)) + time3, conc_eos3 = timed_deplete( + self.chain, conc_eos1, list_rates, dt, matrix_func=cf4_f2) + res3 = self.operator(conc_eos3, power) + + # Step 4: deplete with two matrix exponentials + list_rates = list(zip(bos_rates, res1.rates, res2.rates, res3.rates)) + time4, conc_inter = timed_deplete( + self.chain, bos_conc, list_rates, dt, matrix_func=cf4_f3) + time5, conc_eos5 = timed_deplete( + self.chain, conc_inter, list_rates, dt, matrix_func=cf4_f4) + + return (time1 + time2 + time3 + time4 + time5, + [conc_eos1, conc_eos2, conc_eos3, conc_eos5], + [res1, res2, res3]) + + +class CELIIntegrator(Integrator): + r"""Deplete using the CE/LI CFQ4 algorithm. + + Implements the CE/LI Predictor-Corrector algorithm using the `fourth order + commutator-free integrator `_. + + "CE/LI" stands for constant extrapolation on predictor and linear + interpolation on corrector. This algorithm is mathematically defined as: + + .. math:: + \begin{aligned} + y' &= A(y, t) y(t) \\ + 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 + \end{aligned} + + Parameters + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + 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. + + Attributes + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + chain : openmc.deplete.Chain + Depletion chain + timesteps : iterable of float + Size of each depletion interval in [s] + power : iterable of float + Power of the reactor in [W] for each interval in :attr:`timesteps` + """ + _num_stages = 2 + + def __call__(self, bos_conc, rates, dt, power, _i=None): + """Perform the integration across one time step + + 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 in [W] + _i : int, optional + Current iteration count. Not used + + Returns + ------- + proc_time : float + Time spent in CRAM routines for all materials in [s] + 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) + + # deplete using two matrix exponentials + list_rates = list(zip(rates, res_ce.rates)) + + time_le1, conc_inter = timed_deplete( + self.chain, bos_conc, list_rates, dt, matrix_func=celi_f1) + + time_le2, conc_end = timed_deplete( + self.chain, conc_inter, list_rates, dt, matrix_func=celi_f2) + + return proc_time + time_le1 + time_le1, [conc_ce, conc_end], [res_ce] + + +class EPCRK4Integrator(Integrator): + r"""Deplete using the EPC-RK4 algorithm. + + Implements an extended predictor-corrector algorithm with traditional + Runge-Kutta 4 method. This algorithm is mathematically defined as: + + .. math:: + \begin{aligned} + F_1 &= h A(y_0) \\ + y_1 &= \text{expm}(1/2 F_1) y_0 \\ + F_2 &= h A(y_1) \\ + y_2 &= \text{expm}(1/2 F_2) y_0 \\ + F_3 &= h A(y_2) \\ + y_3 &= \text{expm}(F_3) y_0 \\ + F_4 &= h A(y_3) \\ + y_4 &= \text{expm}(1/6 F_1 + 1/3 F_2 + 1/3 F_3 + 1/6 F_4) y_0 + \end{aligned} + + Parameters + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + 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. + + Attributes + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + chain : openmc.deplete.Chain + Depletion chain + timesteps : iterable of float + Size of each depletion interval in [s] + power : iterable of float + Power of the reactor in [W] for each interval in :attr:`timesteps` + """ + _num_stages = 4 + + def __call__(self, conc, rates, dt, power, _i=None): + """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 in [W] + _i : int, optional + Current depletion step index, unused. + + Returns + ------- + proc_time : float + Time spent in CRAM routines for all materials in [s] + 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 + simulations + """ + + # Step 1: deplete with matrix A(y0) / 2 + time1, conc1 = timed_deplete( + self.chain, conc, rates, dt, matrix_func=rk4_f1) + res1 = self.operator(conc1, power) + + # Step 2: deplete with matrix A(y1) / 2 + time2, conc2 = timed_deplete( + self.chain, conc, res1.rates, dt, matrix_func=rk4_f1) + res2 = self.operator(conc2, power) + + # Step 3: deplete with matrix A(y2) + time3, conc3 = timed_deplete( + self.chain, conc, res2.rates, dt) + res3 = self.operator(conc3, power) + + # Step 4: deplete with matrix built from weighted rates + list_rates = list(zip(rates, res1.rates, res2.rates, res3.rates)) + time4, conc4 = timed_deplete( + self.chain, conc, list_rates, dt, matrix_func=rk4_f4) + + return (time1 + time2 + time3 + time4, [conc1, conc2, conc3, conc4], + [res1, res2, res3]) + + +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. + + Parameters + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + 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. + + Attributes + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + chain : openmc.deplete.Chain + Depletion chain + timesteps : iterable of float + Size of each depletion interval in [s] + power : iterable of float + Power of the reactor in [W] for each interval in :attr:`timesteps` + """ + _num_stages = 2 + + 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 in [W] + i : int + Current depletion step index + + Returns + ------- + proc_time : float + Time spent in CRAM routines for all materials in [s] + 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._i_res < 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] + prev_dt = self.timesteps[i] - prev_res.time[0] + self._prev_rates = prev_res.rates[0] + else: + prev_dt = 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(prev_dt), 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(prev_dt), 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]) + + +class SICELIIntegrator(SIIntegrator): + 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 + `_. + + 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. + n_steps : int, optional + Number of stochastic iterations per depletion interval. + Must be greater than zero. Default : 10 + + Attributes + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + chain : openmc.deplete.Chain + Depletion chain + timesteps : iterable of float + Size of each depletion interval in [s] + power : iterable of float + Power of the reactor in [W] for each interval in :attr:`timesteps` + n_steps : int + Number of stochastic iterations per depletion interval + """ + _num_stages = 2 + + def __call__(self, bos_conc, bos_rates, dt, power, _i=None): + """Perform the integration across one time step + + 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 in [W] + _i : int, optional + Current depletion step index. Not used + + Returns + ------- + proc_time : float + Time spent in CRAM routines for all materials in [s] + 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) + + # Begin iteration + for j in range(self.n_steps + 1): + inter_res = self.operator(inter_conc, power) + + 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) + + 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 + + # end iteration + return proc_time, [eos_conc, inter_conc], [res_bar] + + +class SILEQIIntegrator(SIIntegrator): + 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 + `_. + + 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. + n_steps : int, optional + Number of stochastic iterations per depletion interval. + Must be greater than zero. Default : 10 + + Attributes + ---------- + operator : openmc.deplete.TransportOperator + Operator to perform transport simulations + chain : openmc.deplete.Chain + Depletion chain + timesteps : iterable of float + Size of each depletion interval in [s] + power : iterable of float + Power of the reactor in [W] for each interval in :attr:`timesteps` + n_steps : int + Number of stochastic iterations per depletion interval + """ + _num_stages = 2 + + 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 in [W] + i : int + Current depletion step index + + Returns + ------- + proc_time : float + Time spent in CRAM routines for all materials in [s] + 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._i_res < 1: + self._prev_rates = bos_rates + # Perform CELI for initial steps + return SICELIIntegrator.__call__( + self, bos_conc, bos_rates, dt, power, i) + prev_res = self.operator.prev_res[-2] + prev_dt = self.timesteps[i] - prev_res.time[0] + self._prev_rates = prev_res.rates[0] + else: + prev_dt = self.timesteps[i - 1] + + # Perform remaining LE/QI + inputs = list(zip(self._prev_rates, bos_rates, + repeat(prev_dt), 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) + + proc_time += time1 + inter_conc = copy.deepcopy(eos_conc) + + for j in range(self.n_steps + 1): + inter_res = self.operator(inter_conc, power) + + 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) + + inputs = list(zip(self._prev_rates, bos_rates, res_bar.rates, + repeat(prev_dt), 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 + + return proc_time, [eos_conc, inter_conc], [res_bar] diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 88f73b9958..f8288841f7 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -16,6 +16,7 @@ import xml.etree.ElementTree as ET import h5py import numpy as np +from uncertainties import ufloat import openmc import openmc.capi @@ -56,7 +57,7 @@ class Operator(TransportOperator): Instances of this class can be used to perform depletion using OpenMC as the transport operator. Normally, a user needn't call methods of this class directly. Instead, an instance of this class is passed to an integrator - function, such as :func:`openmc.deplete.integrator.cecm`. + class, such as :class:`openmc.deplete.CECMIntegrator`. Parameters ---------- @@ -113,15 +114,16 @@ class Operator(TransportOperator): Initial heavy metal inventory local_mats : list of str All burnable material IDs being managed by a single process - prev_res : ResultsList - Results from a previous depletion calculation + prev_res : ResultsList or None + Results from a previous depletion calculation. ``None`` if no + results are to be used. diff_burnable_mats : bool Whether to differentiate burnable materials with multiple instances """ def __init__(self, geometry, settings, chain_file=None, prev_results=None, diff_burnable_mats=False, fission_q=None, dilute_initial=1.0e3): - super().__init__(chain_file, fission_q, dilute_initial) + super().__init__(chain_file, fission_q, dilute_initial, prev_results) self.round_number = False self.prev_res = None self.settings = settings @@ -141,7 +143,7 @@ class Operator(TransportOperator): self._mat_index_map = { lm: self.burnable_mats.index(lm) for lm in self.local_mats} - if prev_results is not None: + if self.prev_res is not None: # Reload volumes into geometry prev_results[-1].transfer_volumes(geometry) @@ -176,7 +178,7 @@ class Operator(TransportOperator): self.reaction_rates.n_nuc, self.reaction_rates.n_react) self._energy_helper = ChainFissionHelper() - def __call__(self, vec, power, print_out=True): + def __call__(self, vec, power): """Runs a simulation. Parameters @@ -185,8 +187,6 @@ class Operator(TransportOperator): Total atoms to be used in function. power : float Power of the reactor in [W] - print_out : bool, optional - Whether or not to print out time. Returns ------- @@ -216,15 +216,21 @@ class Operator(TransportOperator): # Extract results op_result = self._unpack_tallies_and_normalize(power) - if comm.rank == 0: - time_unpack = time.time() - - if print_out: - print("Time to openmc: ", time_openmc - time_start) - print("Time to unpack: ", time_unpack - time_openmc) - return copy.deepcopy(op_result) + @staticmethod + def write_bos_data(step): + """Write a state-point file with beginning of step data + + Parameters + ---------- + step : int + Current depletion step including restarts + """ + openmc.capi.statepoint_write( + "openmc_simulation_n{}.h5".format(step), + write_source=False) + def _differentiate_burnable_mats(self): """Assign distribmats for each burnable material @@ -534,7 +540,7 @@ class Operator(TransportOperator): rates.fill(0.0) # Get k and uncertainty - k_combined = openmc.capi.keff() + k_combined = ufloat(*openmc.capi.keff()) # Extract tally bins nuclides = self._rate_helper.nuclides diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index 19911ba7f6..6593d5fc5f 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -21,8 +21,8 @@ class Results(object): Attributes ---------- - k : list of float - Eigenvalue for each substep. + k : list of (float, float) + Eigenvalue and uncertainty for each substep. time : list of float Time at beginning, end of step, in seconds. power : float @@ -456,16 +456,7 @@ class Results(object): # Get indexing terms vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info() - # For a restart calculation, limit number of stages saved to meet the - # format of the hdf5 file stages = len(x) - offset = 0 - if op.prev_res is not None and op.prev_res[0].n_stages < stages: - offset = stages - op.prev_res[0].n_stages - stages = min(stages, op.prev_res[0].n_stages) - warn("Number of restart integrator stages saved limited by initial" - " depletion integrator choice to {}" - .format(op.prev_res[0].n_stages)) # Create results results = Results() @@ -475,9 +466,9 @@ class Results(object): for i in range(stages): for mat_i in range(n_mat): - results[i, mat_i, :] = x[offset + i][mat_i][:] + results[i, mat_i, :] = x[i][mat_i] - results.k = [r.k for r in op_results] + results.k = [(r.k.nominal_value, r.k.std_dev) for r in op_results] results.rates = [r.rates for r in op_results] results.time = t results.power = power diff --git a/tests/dummy_operator.py b/tests/dummy_operator.py index c66070ad40..c23abef631 100644 --- a/tests/dummy_operator.py +++ b/tests/dummy_operator.py @@ -1,5 +1,7 @@ import numpy as np import scipy.sparse as sp +from uncertainties import ufloat + from openmc.deplete.reaction_rates import ReactionRates from openmc.deplete.abc import TransportOperator, OperatorResult @@ -48,7 +50,7 @@ class DummyOperator(TransportOperator): reaction_rates[0, 1, 0] = vec[0][1] # Create a fake rates object - return OperatorResult(0.0, reaction_rates) + return OperatorResult(ufloat(0.0, 0.0), reaction_rates) @property def chain(self): @@ -114,6 +116,9 @@ class DummyOperator(TransportOperator): """Maps cell name to index in global geometry.""" return self.local_mats + @staticmethod + def write_bos_data(_step): + """Empty method but avoids calls to C API""" @property def reaction_rates(self): diff --git a/tests/regression_tests/deplete/test.py b/tests/regression_tests/deplete/test.py index d81bf607b9..c7a9c162e0 100644 --- a/tests/regression_tests/deplete/test.py +++ b/tests/regression_tests/deplete/test.py @@ -39,7 +39,7 @@ def test_full(run_in_tmpdir): space = openmc.stats.Box(lower_left, upper_right) settings.source = openmc.Source(space=space) settings.seed = 1 - settings.verbosity = 3 + settings.verbosity = 1 # Create operator chain_file = Path(__file__).parents[2] / 'chain_simple.xml' @@ -54,7 +54,7 @@ def test_full(run_in_tmpdir): power = 2.337e15*4*JOULE_PER_EV*1e6 # MeV/second cm from CASMO # Perform simulation using the predictor algorithm - openmc.deplete.integrator.predictor(op, dt, power) + openmc.deplete.PredictorIntegrator(op, dt, power).integrate() # Get path to test and reference results path_test = op.output_dir / 'depletion_results.h5' @@ -101,3 +101,18 @@ def test_full(run_in_tmpdir): assert correct, "Discrepancy in mat {} and nuc {}\n{}\n{}".format( mat, nuc, y_old, y_test) + + # Compare statepoint files with depletion results + + t_test, k_test = res_test.get_eigenvalue() + t_ref, k_ref = res_ref.get_eigenvalue() + k_state = np.empty_like(k_ref) + + # Get statepoint files for all BOS points and EOL + for n in range(N + 1): + statepoint = openmc.StatePoint("openmc_simulation_n{}.h5".format(n)) + k_n = statepoint.k_combined + k_state[n] = [k_n.nominal_value, k_n.std_dev] + # Look for exact match pulling from statepoint and depletion_results + assert np.all(k_state == k_test) + assert np.allclose(k_test, k_ref) diff --git a/tests/unit_tests/test_deplete_cecm.py b/tests/unit_tests/test_deplete_cecm.py index f979e0b4d9..65dd61643b 100644 --- a/tests/unit_tests/test_deplete_cecm.py +++ b/tests/unit_tests/test_deplete_cecm.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 CECMIntegrator, ResultsList from tests import dummy_operator @@ -18,10 +18,11 @@ def test_cecm(run_in_tmpdir): # Perform simulation using the MCNPX/MCNP6 algorithm dt = [0.75, 0.75] power = 1.0 - openmc.deplete.cecm(op, dt, power, print_out=False) + integrator = CECMIntegrator(op, dt, power) + integrator.integrate() # Load the files - res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") + res = ResultsList.from_hdf5(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_celi.py b/tests/unit_tests/test_deplete_celi.py index ffa5c304c3..01ec60d372 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.from_hdf5(op.output_dir / "depletion_results.h5") + res = ResultsList.from_hdf5(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_cf4.py b/tests/unit_tests/test_deplete_cf4.py index aee373b406..817dfb5e17 100644 --- a/tests/unit_tests/test_deplete_cf4.py +++ b/tests/unit_tests/test_deplete_cf4.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 CF4Integrator, ResultsList from tests import dummy_operator @@ -18,10 +18,10 @@ def test_cf4(run_in_tmpdir): # Perform simulation using the cf4 algorithm dt = [0.75, 0.75] power = 1.0 - openmc.deplete.cf4(op, dt, power, print_out=False) + CF4Integrator(op, dt, power).integrate() # Load the files - res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") + res = ResultsList.from_hdf5(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_cram.py b/tests/unit_tests/test_deplete_cram.py index 21b93d17e3..8987fbd7aa 100644 --- a/tests/unit_tests/test_deplete_cram.py +++ b/tests/unit_tests/test_deplete_cram.py @@ -6,7 +6,7 @@ Compares a few Mathematica matrix exponentials to CRAM16/CRAM48. from pytest import approx import numpy as np import scipy.sparse as sp -from openmc.deplete.integrator import CRAM16, CRAM48 +from openmc.deplete.cram import CRAM16, CRAM48 def test_CRAM16(): diff --git a/tests/unit_tests/test_deplete_epc_rk4.py b/tests/unit_tests/test_deplete_epc_rk4.py index 1b0af0365a..0c656d7dbe 100644 --- a/tests/unit_tests/test_deplete_epc_rk4.py +++ b/tests/unit_tests/test_deplete_epc_rk4.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 EPCRK4Integrator, ResultsList from tests import dummy_operator @@ -18,10 +18,10 @@ def test_epc_rk4(run_in_tmpdir): # Perform simulation using the epc_rk4 algorithm dt = [0.75, 0.75] power = 1.0 - openmc.deplete.epc_rk4(op, dt, power, print_out=False) + EPCRK4Integrator(op, dt, power).integrate() # Load the files - res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") + res = ResultsList.from_hdf5(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_integrator.py b/tests/unit_tests/test_deplete_integrator.py index b8429392ff..fd4e2ba3ff 100644 --- a/tests/unit_tests/test_deplete_integrator.py +++ b/tests/unit_tests/test_deplete_integrator.py @@ -11,8 +11,12 @@ import os from unittest.mock import MagicMock import numpy as np -from openmc.deplete import (ReactionRates, Results, ResultsList, comm, - OperatorResult) +from uncertainties import ufloat +import pytest + +from openmc.deplete import ( + ReactionRates, Results, ResultsList, comm, OperatorResult, + PredictorIntegrator, SICELIIntegrator) def test_results_save(run_in_tmpdir): @@ -74,8 +78,10 @@ def test_results_save(run_in_tmpdir): t1 = [0.0, 1.0] t2 = [1.0, 2.0] - op_result1 = [OperatorResult(k, rates) for k, rates in zip(eigvl1, rate1)] - op_result2 = [OperatorResult(k, rates) for k, rates in zip(eigvl2, rate2)] + op_result1 = [OperatorResult(ufloat(*k), rates) + for k, rates in zip(eigvl1, rate1)] + op_result2 = [OperatorResult(ufloat(*k), rates) + for k, rates in zip(eigvl2, rate2)] Results.save(op, x1, op_result1, t1, 0, 0) Results.save(op, x2, op_result2, t2, 0, 1) @@ -95,3 +101,32 @@ def test_results_save(run_in_tmpdir): np.testing.assert_array_equal(res[1].k, eigvl2) np.testing.assert_array_equal(res[1].time, t2) + + +@pytest.mark.parametrize("timesteps", (1, [1])) +def test_bad_integrator_inputs(timesteps): + """Test failure modes for Integrator inputs""" + + op = MagicMock() + op.prev_res = None + op.chain = None + op.heavy_metal = 1.0 + + # No power nor power density given + with pytest.raises(ValueError, match="Either power or power density"): + PredictorIntegrator(op, timesteps) + + # Length of power != length time + with pytest.raises(ValueError, match="number of powers"): + PredictorIntegrator(op, timesteps, power=[1, 2]) + + # Length of power density != length time + with pytest.raises(ValueError, match="number of powers"): + PredictorIntegrator(op, timesteps, power_density=[1, 2]) + + # SI integrator with bad steps + with pytest.raises(TypeError, match="n_steps"): + SICELIIntegrator(op, timesteps, [1], n_steps=2.5) + + with pytest.raises(ValueError, match="n_steps"): + SICELIIntegrator(op, timesteps, [1], n_steps=0) diff --git a/tests/unit_tests/test_deplete_leqi.py b/tests/unit_tests/test_deplete_leqi.py index ac477a36ec..113cb6e306 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.from_hdf5(op.output_dir / "depletion_results.h5") + res = ResultsList.from_hdf5(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_operator.py b/tests/unit_tests/test_deplete_operator.py index 45c2aada4e..3574d293dc 100644 --- a/tests/unit_tests/test_deplete_operator.py +++ b/tests/unit_tests/test_deplete_operator.py @@ -37,14 +37,20 @@ def bare_xs(run_in_tmpdir): class BareDepleteOperator(TransportOperator): """Very basic class for testing the initialization.""" - # declare abstract methods so object can be created - def __call__(self, *args, **kwargs): + @staticmethod + def __call__(*args, **kwargs): pass - def initial_condition(self): + @staticmethod + def initial_condition(): pass - def get_results_info(self): + @staticmethod + def get_results_info(): + pass + + @staticmethod + def write_bos_data(): pass diff --git a/tests/unit_tests/test_deplete_predictor.py b/tests/unit_tests/test_deplete_predictor.py index 092ebbc35e..38bee86c39 100644 --- a/tests/unit_tests/test_deplete_predictor.py +++ b/tests/unit_tests/test_deplete_predictor.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 PredictorIntegrator, ResultsList from tests import dummy_operator @@ -18,10 +18,11 @@ def test_predictor(run_in_tmpdir): # Perform simulation using the predictor algorithm dt = [0.75, 0.75] power = 1.0 - openmc.deplete.predictor(op, dt, power, print_out=False) + PredictorIntegrator(op, dt, power).integrate() + # Load the files - res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") + res = ResultsList.from_hdf5(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 4c0c68a2e3..c55f7fd0f4 100644 --- a/tests/unit_tests/test_deplete_restart.py +++ b/tests/unit_tests/test_deplete_restart.py @@ -4,8 +4,12 @@ These tests run in two steps, a first run then a restart run, a simple test problem described in dummy_geometry.py. """ -from pytest import approx +from pytest import approx, raises import openmc.deplete +from openmc.deplete import ( + CECMIntegrator, PredictorIntegrator, CELIIntegrator, LEQIIntegrator, + EPCRK4Integrator, CF4Integrator, SICELIIntegrator, SILEQIIntegrator +) from tests import dummy_operator @@ -20,7 +24,7 @@ def test_restart_predictor(run_in_tmpdir): # Perform simulation using the predictor algorithm dt = [0.75] power = 1.0 - openmc.deplete.predictor(op, dt, power, print_out=False) + PredictorIntegrator(op, dt, power).integrate() # Load the files prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") @@ -30,7 +34,7 @@ def test_restart_predictor(run_in_tmpdir): op.output_dir = output_dir # Perform restarts simulation using the predictor algorithm - openmc.deplete.predictor(op, dt, power, print_out=False) + PredictorIntegrator(op, dt, power).integrate() # Load the files res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") @@ -59,7 +63,8 @@ def test_restart_cecm(run_in_tmpdir): # Perform simulation using the MCNPX/MCNP6 algorithm dt = [0.75] power = 1.0 - openmc.deplete.cecm(op, dt, power, print_out=False) + cecm = CECMIntegrator(op, dt, power) + cecm.integrate() # Load the files prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") @@ -69,7 +74,8 @@ def test_restart_cecm(run_in_tmpdir): op.output_dir = output_dir # Perform restarts simulation using the MCNPX/MCNP6 algorithm - openmc.deplete.cecm(op, dt, power, print_out=False) + cecm_restart = CECMIntegrator(op, dt, power) + cecm_restart.integrate() # Load the files res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") @@ -84,13 +90,12 @@ def test_restart_cecm(run_in_tmpdir): assert y1[1] == approx(s1[0]) assert y2[1] == approx(s1[1]) - assert y1[3] == approx(s2[0]) - assert y2[3] == approx(s2[1]) + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) def test_restart_predictor_cecm(run_in_tmpdir): - """Integral regression test of integrator algorithm using predictor - for the first run then CE/CM for the restart run.""" + """Test to ensure that schemes with different stages are not compatible""" op = dummy_operator.DummyOperator() output_dir = "test_restart_predictor_cecm" @@ -99,7 +104,7 @@ def test_restart_predictor_cecm(run_in_tmpdir): # Perform simulation using the predictor algorithm dt = [0.75] power = 1.0 - openmc.deplete.predictor(op, dt, power, print_out=False) + PredictorIntegrator(op, dt, power).integrate() # Load the files prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") @@ -108,24 +113,9 @@ def test_restart_predictor_cecm(run_in_tmpdir): op = dummy_operator.DummyOperator(prev_res) op.output_dir = output_dir - # Perform restarts simulation using the MCNPX/MCNP6 algorithm - openmc.deplete.cecm(op, dt, power, print_out=False) - - # Load the files - res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") - - _, y1 = res.get_atoms("1", "1") - _, y2 = res.get_atoms("1", "2") - - # Test solution - s1 = [2.46847546272295, 0.986431226850467] - s2 = [3.09106948392, 0.607102912398] - - assert y1[1] == approx(s1[0]) - assert y2[1] == approx(s1[1]) - - assert y1[2] == approx(s2[0]) - assert y2[2] == approx(s2[1]) + # check ValueError is raised, indicating previous and current stages + with raises(ValueError, match="incompatible.* 1.*2"): + CECMIntegrator(op, dt, power) def test_restart_cecm_predictor(run_in_tmpdir): @@ -139,7 +129,8 @@ def test_restart_cecm_predictor(run_in_tmpdir): # Perform simulation using the MCNPX/MCNP6 algorithm dt = [0.75] power = 1.0 - openmc.deplete.cecm(op, dt, power, print_out=False) + cecm = CECMIntegrator(op, dt, power) + cecm.integrate() # Load the files prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") @@ -148,25 +139,9 @@ def test_restart_cecm_predictor(run_in_tmpdir): op = dummy_operator.DummyOperator(prev_res) op.output_dir = output_dir - # Perform restarts simulation using the predictor algorithm - openmc.deplete.predictor(op, dt, power, print_out=False) - - # Load the files - res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") - - _, y1 = res.get_atoms("1", "1") - _, y2 = res.get_atoms("1", "2") - - # Test solution - s1 = [1.86872629872102, 1.395525772416039] - s2 = [3.32776806576, 2.391425905] - - assert y1[1] == approx(s1[0]) - assert y2[1] == approx(s1[1]) - - assert y1[2] == approx(s2[0]) - assert y2[2] == approx(s2[1]) - + # check ValueError is raised, indicating previous and current stages + with raises(ValueError, match="incompatible.* 2.*1"): + PredictorIntegrator(op, dt, power) def test_restart_cf4(run_in_tmpdir): """Integral regression test of integrator algorithm using CF4.""" @@ -178,7 +153,7 @@ def test_restart_cf4(run_in_tmpdir): # Perform simulation dt = [0.75] power = 1.0 - openmc.deplete.cf4(op, dt, power, print_out=False) + CF4Integrator(op, dt, power).integrate() # Load the files prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") @@ -188,7 +163,7 @@ def test_restart_cf4(run_in_tmpdir): op.output_dir = output_dir # Perform restarts simulation - openmc.deplete.cf4(op, dt, power, print_out=False) + CF4Integrator(op, dt, power).integrate() # Load the files res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") @@ -203,8 +178,8 @@ def test_restart_cf4(run_in_tmpdir): assert y1[1] == approx(s1[0]) assert y2[1] == approx(s1[1]) - assert y1[3] == approx(s2[0]) - assert y2[3] == approx(s2[1]) + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) def test_restart_epc_rk4(run_in_tmpdir): @@ -217,7 +192,7 @@ def test_restart_epc_rk4(run_in_tmpdir): # Perform simulation dt = [0.75] power = 1.0 - openmc.deplete.epc_rk4(op, dt, power, print_out=False) + EPCRK4Integrator(op, dt, power).integrate() # Load the files prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") @@ -227,7 +202,7 @@ def test_restart_epc_rk4(run_in_tmpdir): op.output_dir = output_dir # Perform restarts simulation - openmc.deplete.epc_rk4(op, dt, power, print_out=False) + EPCRK4Integrator(op, dt, power).integrate() # Load the files res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") @@ -242,8 +217,8 @@ def test_restart_epc_rk4(run_in_tmpdir): assert y1[1] == approx(s1[0]) assert y2[1] == approx(s1[1]) - assert y1[3] == approx(s2[0]) - assert y2[3] == approx(s2[1]) + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) def test_restart_celi(run_in_tmpdir): @@ -256,7 +231,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.from_hdf5(op.output_dir / "depletion_results.h5") @@ -266,7 +241,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.from_hdf5(op.output_dir / "depletion_results.h5") @@ -281,8 +256,8 @@ def test_restart_celi(run_in_tmpdir): assert y1[1] == approx(s1[0]) assert y2[1] == approx(s1[1]) - assert y1[3] == approx(s2[0]) - assert y2[3] == approx(s2[1]) + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) def test_restart_leqi(run_in_tmpdir): @@ -295,7 +270,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.from_hdf5(op.output_dir / "depletion_results.h5") @@ -305,7 +280,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.from_hdf5(op.output_dir / "depletion_results.h5") @@ -320,8 +295,8 @@ def test_restart_leqi(run_in_tmpdir): assert y1[1] == approx(s1[0]) assert y2[1] == approx(s1[1]) - assert y1[3] == approx(s2[0]) - assert y2[3] == approx(s2[1]) + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) def test_restart_si_celi(run_in_tmpdir): """Integral regression test of integrator algorithm using SI-CELI.""" @@ -333,7 +308,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) + SICELIIntegrator(op, dt, power).integrate() # Load the files prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") @@ -343,7 +318,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) + SICELIIntegrator(op, dt, power).integrate() # Load the files res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") @@ -358,8 +333,8 @@ def test_restart_si_celi(run_in_tmpdir): assert y1[1] == approx(s1[0]) assert y2[1] == approx(s1[1]) - assert y1[3] == approx(s2[0]) - assert y2[3] == approx(s2[1]) + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) def test_restart_si_leqi(run_in_tmpdir): @@ -372,7 +347,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 + SILEQIIntegrator(op, dt, power, nstages).integrate() # Load the files prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") @@ -382,7 +358,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) + SILEQIIntegrator(op, dt, power, nstages).integrate() # Load the files res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") @@ -397,5 +373,5 @@ def test_restart_si_leqi(run_in_tmpdir): assert y1[1] == approx(s1[0]) assert y2[1] == approx(s1[1]) - assert y1[3] == approx(s2[0]) - assert y2[3] == approx(s2[1]) + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) diff --git a/tests/unit_tests/test_deplete_si_celi.py b/tests/unit_tests/test_deplete_si_celi.py index e4ea4f6acf..5cc0d4e6c4 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 SICELIIntegrator, 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) + SICELIIntegrator(op, dt, power).integrate() # Load the files - res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") + res = ResultsList.from_hdf5(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_si_leqi.py b/tests/unit_tests/test_deplete_si_leqi.py index cef64fbe88..8e3c3a4648 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 SILEQIIntegrator, 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) + SILEQIIntegrator(op, dt, power, 10).integrate() # Load the files - res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") + res = ResultsList.from_hdf5(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_transfer_volumes.py b/tests/unit_tests/test_transfer_volumes.py index 38520ea1e3..2e735c911e 100644 --- a/tests/unit_tests/test_transfer_volumes.py +++ b/tests/unit_tests/test_transfer_volumes.py @@ -4,7 +4,8 @@ """ from pytest import approx -import openmc.deplete +import openmc +from openmc.deplete import PredictorIntegrator, ResultsList from tests import dummy_operator @@ -18,7 +19,7 @@ def test_transfer_volumes(run_in_tmpdir): # Perform simulation using the predictor algorithm dt = [0.75] power = 1.0 - openmc.deplete.predictor(op, dt, power, print_out=False) + PredictorIntegrator(op, dt, power).integrate() # Load the files res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")