Address comments on #1934 from @ojschumann

docs/source/io_formats/particle_restart.rst

@@ -32,3 +32,4 @@ The current version of the particle restart file format is 2.0.
              multi-group mode.
            - **xyz** (*double[3]*) -- Position of the particle.
            - **uvw** (*double[3]*) -- Direction of the particle.
+           - **time** (*double*) -- Time of the particle in [s].

src/particle.cpp

@@ -54,7 +54,7 @@ double Particle::speed() const
   }
 
   // Calculate inverse of Lorentz factor
-  double inv_gamma = mass / (this->E() + mass);
+  const double inv_gamma = mass / (this->E() + mass);
 
   // Calculate speed via v = c * sqrt(1 - γ^-2)
   return C_LIGHT * std::sqrt(1 - inv_gamma * inv_gamma);
@@ -685,6 +685,7 @@ void Particle::write_restart() const
       write_dataset(file_id, "energy", simulation::source_bank[i - 1].E);
       write_dataset(file_id, "xyz", simulation::source_bank[i - 1].r);
       write_dataset(file_id, "uvw", simulation::source_bank[i - 1].u);
+      write_dataset(file_id, "time", simulation::source_bank[i - 1].time);
     } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
       // re-sample using rng random number seed used to generate source particle
       int64_t id = (simulation::total_gen + overall_generation() - 1) *
@@ -697,6 +698,7 @@ void Particle::write_restart() const
       write_dataset(file_id, "energy", site.E);
       write_dataset(file_id, "xyz", site.r);
       write_dataset(file_id, "uvw", site.u);
+      write_dataset(file_id, "time", site.time);
     }
 
     // Close file

src/particle_restart.cpp

@@ -50,6 +50,7 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode)
   read_dataset(file_id, "energy", p.E());
   read_dataset(file_id, "xyz", p.r());
   read_dataset(file_id, "uvw", p.u());
+  read_dataset(file_id, "time", p.time());
 
   // Set energy group and average energy in multi-group mode
   if (!settings::run_CE) {
@@ -64,6 +65,7 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode)
   p.u_last() = p.u();
   p.E_last() = p.E();
   p.g_last() = p.g();
+  p.time_last() = p.time();
 
   // Close hdf5 file
   file_close(file_id);

src/tallies/filter_time.cpp

@@ -1,6 +1,7 @@
 #include "openmc/tallies/filter_time.h"
 
-#include <algorithm> // for min, max
+#include <algorithm> // for min, max, copy
+#include <iterator>  // for back_inserter
 
 #include <fmt/core.h>
 
@@ -25,14 +26,13 @@ void TimeFilter::set_bins(gsl::span<const double> bins)
   bins_.clear();
   bins_.reserve(bins.size());
 
-  // Copy bins, ensuring they are valid
-  for (gsl::index i = 0; i < bins.size(); ++i) {
-    if (i > 0 && bins[i] <= bins[i - 1]) {
-      throw std::runtime_error {"Time bins must be monotonically increasing."};
-    }
-    bins_.push_back(bins[i]);
+  // Ensure time bins are sorted
+  if (!std::is_sorted(bins.cbegin(), bins.cend(), std::less_equal<double>())) {
+    throw std::runtime_error {"Time bins must be monotonically increasing."};
   }
 
+  // Copy bins
+  std::copy(bins.cbegin(), bins.cend(), std::back_inserter(bins_));
   n_bins_ = bins_.size() - 1;
 }
 
@@ -40,8 +40,8 @@ void TimeFilter::get_all_bins(
   const Particle& p, TallyEstimator estimator, FilterMatch& match) const
 {
   // Get the start/end time of the particle for this track
-  auto t_start = p.time_last();
-  auto t_end = p.time();
+  const auto t_start = p.time_last();
+  const auto t_end = p.time();
 
   // If time interval is entirely out of time bin range, exit
   if (t_end < bins_.front() || t_start >= bins_.back())
@@ -51,7 +51,7 @@ void TimeFilter::get_all_bins(
     // -------------------------------------------------------------------------
     // For surface tallies, find a match based on the exact time the particle
     // crosses the surface
-    auto i_bin = lower_bound_index(bins_.begin(), bins_.end(), t_end);
+    const auto i_bin = lower_bound_index(bins_.begin(), bins_.end(), t_end);
     match.bins_.push_back(i_bin);
     match.weights_.push_back(1.0);
 
@@ -60,6 +60,10 @@ void TimeFilter::get_all_bins(
     // For volume tallies, we have to check the start/end time of the current
     // track and find where it overlaps with time bins and score accordingly
 
+    // Skip if time interval is zero
+    if (t_start == t_end)
+      return;
+
     // Determine first bin containing a portion of time interval
     auto i_bin = lower_bound_index(bins_.begin(), bins_.end(), t_start);
 
@@ -71,11 +75,9 @@ void TimeFilter::get_all_bins(
 
       // Add match with weight equal to the fraction of the time interval within
       // the current time bin
-      if (dt_total > 0.0) {
-        double fraction = (t_right - t_left) / dt_total;
-        match.bins_.push_back(i_bin);
-        match.weights_.push_back(fraction);
-      }
+      const double fraction = (t_right - t_left) / dt_total;
+      match.bins_.push_back(i_bin);
+      match.weights_.push_back(fraction);
 
       if (t_end < bins_[i_bin + 1])
         break;