diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index c1d7a36bd3d..aabf8562311 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -760,7 +760,8 @@ def _get_bos_data_from_restart(self, step_index, source_rate, bos_conc):
         k = ufloat(res.k[0, 0], res.k[0, 1])
 
         # Scale reaction rates by ratio of source rates
-        rates *= source_rate / res.source_rate
+        if res.source_rate != 0.0:
+            rates *= source_rate / res.source_rate
         return bos_conc, OperatorResult(k, rates)
 
     def _get_start_data(self):
@@ -808,8 +809,12 @@ def integrate(self, final_step=True, output=True):
                 # Solve transport equation (or obtain result from restart)
                 if i > 0 or self.operator.prev_res is None:
                     # Update geometry/material according to batchwise definition
-                    if self.batchwise and source_rate != 0.0:
-                        n, root = self._get_bos_from_batchwise(i, n)
+                    if self.batchwise:
+                        if source_rate != 0.0:
+                            n, root = self._get_bos_from_batchwise(i, n)
+                        else:
+                            # Store root at previous timestep
+                            root = self.batchwise._get_cell_attrib()
                     else:
                         root = None
                     n, res = self._get_bos_data_from_operator(i, source_rate, n)
@@ -817,7 +822,15 @@ def integrate(self, final_step=True, output=True):
                     n, res = self._get_bos_data_from_restart(i, source_rate, n)
                     if self.batchwise:
                         root = self.operator.prev_res[-1].batchwise
-                        self.batchwise.update_from_restart(n, root)
+                        #TODO: this is just temporary (import math)
+                        import math
+                        if math.isnan(root):
+                            prev_res_ts = -2
+                            while (math.isnan(root)):
+                                root = self.operator.prev_res[prev_res_ts].batchwise
+                                prev_res_ts -= 1
+
+                        self.batchwise.update_from_restart(i, n, root)
                     else:
                         root = None
 
diff --git a/openmc/deplete/batchwise.py b/openmc/deplete/batchwise.py
index dd0ef39d360..ba34120a051 100644
--- a/openmc/deplete/batchwise.py
+++ b/openmc/deplete/batchwise.py
@@ -180,7 +180,7 @@ def _model_builder(self, param):
         """
 
     @abstractmethod
-    def update_from_restart(self, x, root):
+    def update_from_restart(self, step_index, x, root):
         """
         """
 
@@ -812,9 +812,9 @@ def search_for_keff(self, x, step_index):
 
         return x, root
 
-    def update_from_restart(self, x, root):
+    def update_from_restart(self, step_index, x, root):
         self._set_cell_attrib(root)
-        super()._update_materials(x)
+        super()._update_materials(x, step_index)
 
 class BatchwiseCellGeometrical(BatchwiseCell):
     """
@@ -1222,8 +1222,8 @@ def _calculate_volumes(self, res):
                     volumes[mat_id] = number_i.get_mat_volume(mat_id)
         return volumes
 
-    def update_from_restart(self, x, root):
-        super()._update_materials(x)
+    def update_from_restart(self, step_index, x, root):
+        super()._update_materials(x, step_index)
 
 class BatchwiseMaterialRefuel(BatchwiseMaterial):
     """
@@ -1938,7 +1938,7 @@ def set_density_function(self, mats, density_func, oxidation_states):
     def add_material(self, mat, value, mat_vector, timestep, units):
         self.bw_geom.add_material(mat, value, mat_vector, timestep, units)
 
-    def update_from_restart(self, x, root):
+    def update_from_restart(self, step_index, x, root):
         """
         This is for a restart. TODO update abc class
         Parameters
@@ -1946,7 +1946,7 @@ def update_from_restart(self, x, root):
         x : list of numpy.ndarray
             Total atom concentrations
         """
-        self.bw_geom.update_from_restart(x, root)
+        self.bw_geom.update_from_restart(step_index, x, root)
 
     def search_for_keff(self, x, step_index):
         """
@@ -2046,7 +2046,7 @@ def set_density_function(self, mats, density_func, oxidation_states):
     def add_material(self, mat, value, mat_vector, timestep, units):
         self.bw_geom.add_material(mat, value, mat_vector, timestep, units)
 
-    def update_from_restart(self, x, root):
+    def update_from_restart(self, step_index, x, root):
         """
         This is for a restart. TODO update abc class
         Parameters
@@ -2054,7 +2054,7 @@ def update_from_restart(self, x, root):
         x : list of numpy.ndarray
             Total atom concentrations
         """
-        self.bw_geom.update_from_restart(x, root)
+        self.bw_geom.update_from_restart(step_index, x, root)
 
     def search_for_keff(self, x, step_index):
         """