Fixes two related bugs affecting depletion restarts

openmc/deplete/abc.py

@@ -882,7 +882,23 @@ def _restore_keff_search_control(self, res: StepResult):
 
     def _get_bos_data(self, step_index, source_rate, bos_conc):
         """Get beginning-of-step concentrations, rates, and control state."""
-        if step_index > 0 or self.operator.prev_res is None:
+        prev_res = self.operator.prev_res
+
+        # Restart data can only be reused for the very first step of this
+        # Integrator, and only if the last StepResult in prev_res actually
+        # came from a real transport/eigenvalue solve. When a previous
+        # integrate() call used final_step=False, its last saved StepResult
+        # is a zero-filled placeholder (see the "final simulation" block in
+        # integrate() below) and must NOT be treated as valid BOS data --
+        # doing so silently zeroes out reaction rates (and therefore
+        # depletion) for the continued run.
+        use_restart = (
+            step_index == 0
+            and prev_res is not None
+            and prev_res[-1].evaluated
+        )
+
+        if not use_restart:
             if self._keff_search_control is not None and source_rate != 0.0:
                 keff_search_root = self._keff_search_control.run(bos_conc)
             else:
@@ -893,7 +909,7 @@ def _get_bos_data(self, step_index, source_rate, bos_conc):
             bos_conc, res = self._get_bos_data_from_restart(
                 source_rate, bos_conc)
             if self._keff_search_control is not None and source_rate != 0.0:
-                keff_search_root = self._restore_keff_search_control(self.operator.prev_res[-1])
+                keff_search_root = self._restore_keff_search_control(prev_res[-1])
             else:
                 keff_search_root = None
 
@@ -953,7 +969,8 @@ def integrate(
                     proc_time,
                     write_rates=write_rates,
                     keff_search_root=keff_search_root,
-                    path=path
+                    path=path,
+                    evaluated=True,
                 )
 
                 # Update for next step
@@ -981,9 +998,21 @@ def integrate(
                 proc_time,
                 write_rates=write_rates,
                 keff_search_root=keff_search_root,
-                path=path
+                path=path,
+                evaluated=final_step,
             )
-            self.operator.write_bos_data(len(self) + self._i_res)
+            # Only document beginning-of-step data (which triggers a
+            # statepoint write) when a real transport/eigenvalue solve was
+            # actually performed above. When final_step=False, res_final is
+            # a zero-source-rate placeholder -- the operator was just
+            # reset(), never run() -- so writing a statepoint here would
+            # capture leftover/inconsistent internal simulation state
+            # (observed as a nonzero k-effective with infinite std. dev.).
+            # The real statepoint for this step gets written on a
+            # subsequent continue_timesteps=True restart, via
+            # _get_bos_data_from_operator's call to write_bos_data.
+            if final_step:
+                self.operator.write_bos_data(len(self) + self._i_res)
 
         self.operator.finalize()
 

openmc/deplete/stepresult.py

@@ -17,7 +17,7 @@
 
 from .reaction_rates import ReactionRates
 
-VERSION_RESULTS = (1, 3)
+VERSION_RESULTS = (1, 4)
 
 
 __all__ = ["StepResult"]
@@ -60,6 +60,11 @@ class StepResult:
         materials and processes
     keff_search_root : float
         The root returned by the keff search control.
+    evaluated : bool
+        Whether this step's ``k`` and ``rates`` came from a real
+        transport/eigenvalue solve (``True``) or are a zero-filled
+        placeholder written when ``Integrator.integrate(final_step=False)``
+        skips the final evaluation (``False``).
 
     """
     def __init__(self):
@@ -69,6 +74,7 @@ def __init__(self):
         self.rates = None
         self.volume = None
         self.proc_time = None
+        self.evaluated = None
 
         self.index_mat = None
         self.index_nuc = None
@@ -192,7 +198,8 @@ def distribute(self, local_materials, ranges):
 
         # Direct transfer
         direct_attrs = ("time", "k", "source_rate", "index_nuc",
-                        "mat_to_hdf5_ind", "mat_to_name", "proc_time")
+                        "mat_to_hdf5_ind", "mat_to_name", "proc_time",
+                        "evaluated")
         for attr in direct_attrs:
             setattr(new, attr, getattr(self, attr))
         # Get applicable slice of data
@@ -375,6 +382,10 @@ def _write_hdf5_metadata(self, handle, write_rates):
             "keff_search_root", (1,), maxshape=(None,),
             dtype="float64")
 
+        handle.create_dataset(
+            "evaluated", (1,), maxshape=(None,),
+            dtype="bool")
+
     def _to_hdf5(self, handle, index, parallel=False, write_rates: bool = False):
         """Converts results object into an hdf5 object.
 
@@ -408,6 +419,7 @@ def _to_hdf5(self, handle, index, parallel=False, write_rates: bool = False):
         source_rate_dset = handle["/source_rate"]
         proc_time_dset = handle["/depletion time"]
         keff_search_root_dset = handle["/keff_search_root"]
+        evaluated_dset = handle["/evaluated"]
 
         # Get number of results stored
         number_shape = list(number_dset.shape)
@@ -445,6 +457,10 @@ def _to_hdf5(self, handle, index, parallel=False, write_rates: bool = False):
             keff_search_root_shape[0] = new_shape
             keff_search_root_dset.resize(keff_search_root_shape)
 
+            evaluated_shape = list(evaluated_dset.shape)
+            evaluated_shape[0] = new_shape
+            evaluated_dset.resize(evaluated_shape)
+
         # If nothing to write, just return
         if len(self.index_mat) == 0:
             return
@@ -465,6 +481,9 @@ def _to_hdf5(self, handle, index, parallel=False, write_rates: bool = False):
                     self.proc_time / (comm.size * self.n_hdf5_mats)
                 )
             keff_search_root_dset[index] = self.keff_search_root
+            evaluated_dset[index] = (
+                True if self.evaluated is None else self.evaluated
+            )
 
     @classmethod
     def from_hdf5(cls, handle, step):
@@ -517,6 +536,18 @@ def from_hdf5(cls, handle, step):
             keff_search_root_dset = handle["/keff_search_root"]
             results.keff_search_root = keff_search_root_dset[step]
 
+        if "evaluated" in handle:
+            evaluated_dset = handle["/evaluated"]
+            results.evaluated = bool(evaluated_dset[step])
+        else:
+            # Older results files (pre VERSION_RESULTS (1, 4)) did not track
+            # whether a step's k/rates came from a real solve or were a
+            # final_step=False placeholder. Assume valid so that restarts
+            # against pre-existing files keep behaving as before -- this
+            # only matters for continuation runs, and files from before the
+            # continue_timesteps feature existed never hit this code path.
+            results.evaluated = True
+
         if results.proc_time is None:
             results.proc_time = np.array([np.nan])
 
@@ -572,7 +603,8 @@ def save(
         proc_time=None,
         write_rates: bool = False,
         keff_search_root=None,
-        path: PathLike = "depletion_results.h5"
+        path: PathLike = "depletion_results.h5",
+        evaluated: bool = True,
     ):
         """Creates and writes depletion results to disk
 
@@ -602,6 +634,11 @@ def save(
             Path to file to write. Defaults to 'depletion_results.h5'.
 
             .. versionadded:: 0.14.0
+        evaluated : bool, optional
+            Whether ``op_results`` came from a real transport/eigenvalue
+            solve. Should be set to ``False`` only for the zero-filled
+            placeholder written when ``final_step=False``. Defaults to
+            ``True``.
         """
         # Get indexing terms
         vol_dict, nuc_list, burn_list, full_burn_list, name_list = op.get_results_info()
@@ -626,6 +663,7 @@ def save(
         if results.proc_time is not None:
             results.proc_time = comm.reduce(proc_time, op=MPI.SUM)
         results.keff_search_root = keff_search_root
+        results.evaluated = evaluated
 
         if not Path(path).is_file():
             Path(path).parent.mkdir(parents=True, exist_ok=True)