paulromano · Mar 6, 2025 · Jan 22, 2025 · Jan 28, 2025 · Jan 28, 2025 · Jan 28, 2025
@@ -1,6 +1,7 @@
 """abc module.
 
-This module contains Abstract Base Classes for implementing operator, integrator, depletion system solver, and operator helper classes
+This module contains Abstract Base Classes for implementing operator,
+integrator, depletion system solver, and operator helper classes
 """
 
 from __future__ import annotations
@@ -24,7 +25,8 @@
 from openmc import Material
 from .stepresult import StepResult
 from .chain import Chain
-from .results import Results
+from .results import Results, _SECONDS_PER_MINUTE, _SECONDS_PER_HOUR, \
+    _SECONDS_PER_DAY, _SECONDS_PER_JULIAN_YEAR
 from .pool import deplete
 from .reaction_rates import ReactionRates
 from .transfer_rates import TransferRates
@@ -36,12 +38,6 @@
     "Integrator", "SIIntegrator", "DepSystemSolver", "add_params"]
 
 
-_SECONDS_PER_MINUTE = 60
-_SECONDS_PER_HOUR = 60*60
-_SECONDS_PER_DAY = 24*60*60
-_SECONDS_PER_JULIAN_YEAR = 365.25*24*60*60
-
-
 def _normalize_timesteps(
         timesteps: Sequence[float] | Sequence[tuple[float, str]],
         source_rates: float | Sequence[float],
@@ -572,6 +568,18 @@ class Integrator(ABC):
         :attr:`solver`.
 
         .. versionadded:: 0.12
+    continue_timesteps : bool, optional
+        Whether or not to treat the current solve as a continuation of a
+        previous simulation. Defaults to `False`. When `False`, the depletion
+        steps provided are appended to any previous steps. If `True`, the
+        timesteps provided to the `Integrator` must exacly match any that
+        exist in the `prev_results` passed to the `Operator`. The `power`,
+        `power_density`, or `source_rates` must match as well. The
+        method of specifying `power`, `power_density`, or
+        `source_rates` should be the same as the initial run.
+
+        .. versionadded:: 0.15.1
+
     Attributes
     ----------
     operator : openmc.deplete.abc.TransportOperator
@@ -613,7 +621,8 @@ def __init__(
             power_density: Optional[Union[float, Sequence[float]]] = None,
             source_rates: Optional[Union[float, Sequence[float]]] = None,
             timestep_units: str = 's',
-            solver: str = "cram48"
+            solver: str = "cram48",
+            continue_timesteps: bool = False,
         ):
         # Check number of stages previously used
         if operator.prev_res is not None:
@@ -625,6 +634,8 @@ def __init__(
                     "this uses {}".format(
                         self.__class__.__name__, res.data.shape[0],
                         self._num_stages))
+        elif continue_timesteps:
+            raise ValueError("Continuation run requires passing prev_results.")
         self.operator = operator
         self.chain = operator.chain
 
@@ -642,6 +653,35 @@ def __init__(
         # Normalize timesteps and source rates
         seconds, source_rates = _normalize_timesteps(
             timesteps, source_rates, timestep_units, operator)
+
+        if continue_timesteps:
+            # Get timesteps and source rates from previous results
+            prev_times = operator.prev_res.get_times(timestep_units)
+            prev_source_rates = operator.prev_res.get_source_rates()
+            prev_timesteps = np.diff(prev_times)
+
+            # Make sure parameters from the previous results are consistent with
+            # those passed to operator
+            num_prev = len(prev_timesteps)
+            if not np.array_equal(prev_timesteps, timesteps[:num_prev]):
+                raise ValueError(
+                    "You are attempting to continue a run in which the previous timesteps "
+                    "do not have the same initial timesteps as those provided to the "
+                    "Integrator. Please make sure you are using the correct timesteps."
+                )
+            if not np.array_equal(prev_source_rates, source_rates[:num_prev]):
+                raise ValueError(
+                    "You are attempting to continue a run in which the previous results "
+                    "do not have the same initial source rates, powers, or power densities "
+                    "as those provided to the Integrator. Please make sure you are using "
+                    "the correct powers, power densities, or source rates and previous "
+                    "results file."
+                )
+
+            # Run with only the new time steps and source rates provided
+            seconds = seconds[num_prev:]
+            source_rates = source_rates[num_prev:]
+
         self.timesteps = np.asarray(seconds)
         self.source_rates = np.asarray(source_rates)
 
@@ -932,6 +972,18 @@ class SIIntegrator(Integrator):
         :attr:`solver`.
 
         .. versionadded:: 0.12
+    continue_timesteps : bool, optional
+        Whether or not to treat the current solve as a continuation of a
+        previous simulation. Defaults to `False`. If `False`, all time
+        steps and source rates will be run in an append fashion and will run
+        after whatever time steps exist, if any. If `True`, the timesteps
+        provided to the `Integrator` must match exactly those that exist
+        in the `prev_results` passed to the `Opereator`. The `power`,
+        `power_density`, or `source_rates` must match as well. The
+        method of specifying `power`, `power_density`, or
+        `source_rates` should be the same as the initial run.
+
+        .. versionadded:: 0.15.1
 
     Attributes
     ----------
@@ -973,13 +1025,14 @@ def __init__(
             source_rates: Optional[Sequence[float]] = None,
             timestep_units: str = 's',
             n_steps: int = 10,
-            solver: str = "cram48"
+            solver: str = "cram48",
+            continue_timesteps: bool = False,
         ):
         check_type("n_steps", n_steps, Integral)
         check_greater_than("n_steps", n_steps, 0)
         super().__init__(
             operator, timesteps, power, power_density, source_rates,
-            timestep_units=timestep_units, solver=solver)
+            timestep_units=timestep_units, solver=solver, continue_timesteps=continue_timesteps)
         self.n_steps = n_steps
 
     def _get_bos_data_from_operator(self, step_index, step_power, n_bos):

@@ -17,6 +17,11 @@
 
 __all__ = ["Results", "ResultsList"]
 
+_SECONDS_PER_MINUTE = 60
+_SECONDS_PER_HOUR = 60*60
+_SECONDS_PER_DAY = 24*60*60
+_SECONDS_PER_JULIAN_YEAR = 365.25*24*60*60  # 365.25 due to the leap year
+
 
 def _get_time_as(seconds: float, units: str) -> float:
     """Converts the time in seconds to time in different units
@@ -31,13 +36,13 @@ def _get_time_as(seconds: float, units: str) -> float:
 
     """
     if units == "a":
-        return seconds / (60 * 60 * 24 * 365.25)  # 365.25 due to the leap year
+        return seconds / _SECONDS_PER_JULIAN_YEAR
     if units == "d":
-        return seconds / (60 * 60 * 24)
+        return seconds / _SECONDS_PER_DAY
     elif units == "h":
-        return seconds / (60 * 60)
+        return seconds / _SECONDS_PER_HOUR
     elif units == "min":
-        return seconds / 60
+        return seconds / _SECONDS_PER_MINUTE
     else:
         return seconds
 
@@ -71,7 +76,6 @@ def __init__(self, filename='depletion_results.h5'):
                     data.append(StepResult.from_hdf5(fh, i))
         super().__init__(data)
 
-
     @classmethod
     def from_hdf5(cls, filename: PathLike):
         """Load in depletion results from a previous file
@@ -460,6 +464,26 @@ def get_times(self, time_units: str = "d") -> np.ndarray:
 
         return _get_time_as(times, time_units)
 
+    def get_source_rates(self) -> np.ndarray:
+        """
+        .. versionadded:: 0.15.1
+
+        Returns
+        -------
+        numpy.ndarray
+            1-D vector of source rates at each point in the depletion simulation
+            with the units originally defined by the user.
+
+        """
+        # Results duplicate the final source rate at the final simulation time
+        source_rates = np.fromiter(
+            (r.source_rate for r in self),
+            dtype=self[0].source_rate.dtype,
+            count=len(self)-1,
+        )
+
+        return source_rates
+
     def get_step_where(
         self, time, time_units: str = "d", atol: float = 1e-6, rtol: float = 1e-3
     ) -> int:

@@ -0,0 +1,90 @@
+"""Unit tests for openmc.deplete continue run capability.
+
+These tests run in two steps: first a normal run and then a continue run using the previous results
+"""
+
+import pytest
+
+import openmc.deplete
+
+from tests import dummy_operator
+
+
+@pytest.mark.parametrize("scheme", dummy_operator.SCHEMES)
+def test_continue(run_in_tmpdir, scheme):
+    """Test to ensure that a properly defined continue run works"""
+    # set up the problem
+    bundle = dummy_operator.SCHEMES[scheme]
+    operator = dummy_operator.DummyOperator()
+
+    # initial depletion
+    bundle.solver(operator, [1.0, 2.0], [1.0, 2.0]).integrate()
+
+    # set up continue run
+    prev_res = openmc.deplete.Results(operator.output_dir / "depletion_results.h5")
+    operator = dummy_operator.DummyOperator(prev_res)
+
+    # if continue run happens, test passes
+    bundle.solver(operator, [1.0, 2.0, 3.0, 4.0], [1.0, 2.0, 3.0, 4.0], continue_timesteps=True).integrate()
+
+
+@pytest.mark.parametrize("scheme", dummy_operator.SCHEMES)
+def test_continue_for_null_previous(run_in_tmpdir, scheme):
+    """Test to ensure that a continue run works even if there are no previous results"""
+    # set up the problem
+    bundle = dummy_operator.SCHEMES[scheme]
+    operator = dummy_operator.DummyOperator()
+
+    # initial depletion
+    bundle.solver(operator, [1.0, 2.0], [1.0, 2.0], continue_timesteps=True).integrate()
+
+
+@pytest.mark.parametrize("scheme", dummy_operator.SCHEMES)
+def test_mismatched_initial_times(run_in_tmpdir, scheme):
+    """Test to ensure that a continue run with different initial steps is properly caught"""
+    # set up the problem
+    bundle = dummy_operator.SCHEMES[scheme]
+    operator = dummy_operator.DummyOperator()
+
+    # perform initial steps
+    bundle.solver(operator, [0.75, 0.75], [1.0, 1.0]).integrate()
+
+    # restart
+    prev_res = openmc.deplete.Results(operator.output_dir / "depletion_results.h5")
+    operator = dummy_operator.DummyOperator(prev_res)
+
+    with pytest.raises(
+        ValueError,
+        match="You are attempting to continue a run in which the previous timesteps "
+        "do not have the same initial timesteps as those provided to the "
+        "Integrator. Please make sure you are using the correct timesteps.",
+    ):
+        bundle.solver(
+            operator, [0.75, 0.5, 0.75], [1.0, 1.0, 1.0], continue_timesteps=True
+        ).integrate()
+
+
+@pytest.mark.parametrize("scheme", dummy_operator.SCHEMES)
+def test_mismatched_initial_source_rates(run_in_tmpdir, scheme):
+    """Test to ensure that a continue run with different initial steps is properly caught"""
+    # set up the problem
+    bundle = dummy_operator.SCHEMES[scheme]
+    operator = dummy_operator.DummyOperator()
+
+    # perform initial steps
+    bundle.solver(operator, [0.75, 0.75], [1.0, 1.0]).integrate()
+
+    # restart
+    prev_res = openmc.deplete.Results(operator.output_dir / "depletion_results.h5")
+    operator = dummy_operator.DummyOperator(prev_res)
+
+    with pytest.raises(
+        ValueError,
+        match="You are attempting to continue a run in which the previous results "
+        "do not have the same initial source rates, powers, or power densities "
+        "as those provided to the Integrator. Please make sure you are using "
+        "the correct powers, power densities, or source rates and previous results file.",
+    ):
+        bundle.solver(
+            operator, [0.75, 0.75, 0.75], [1.0, 2.0, 1.0], continue_timesteps=True
+        ).integrate()