Initial implementation of time filter

CMakeLists.txt

@@ -350,6 +350,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_sph_harm.cpp
   src/tallies/filter_sptl_legendre.cpp
   src/tallies/filter_surface.cpp
+  src/tallies/filter_time.cpp
   src/tallies/filter_universe.cpp
   src/tallies/filter_zernike.cpp
   src/tallies/tally.cpp

docs/source/pythonapi/base.rst

@@ -136,6 +136,7 @@ Constructing Tallies
    openmc.LegendreFilter
    openmc.SpatialLegendreFilter
    openmc.SphericalHarmonicsFilter
+   openmc.TimeFilter
    openmc.ZernikeFilter
    openmc.ZernikeRadialFilter
    openmc.ParticleFilter

include/openmc/constants.h

@@ -90,7 +90,7 @@ constexpr double FINE_STRUCTURE {
 constexpr double PLANCK_C {
   1.2398419839593942e4}; // Planck's constant times c in eV-Angstroms
 constexpr double AMU {1.66053906660e-27};      // 1 amu in kg
-constexpr double C_LIGHT {2.99792458e8};       // speed of light in m/s
+constexpr double C_LIGHT {2.99792458e10};      // speed of light in cm/s
 constexpr double N_AVOGADRO {0.602214076};     // Avogadro's number in 10^24/mol
 constexpr double K_BOLTZMANN {8.617333262e-5}; // Boltzmann constant in eV/K
 

include/openmc/particle.h

@@ -35,6 +35,8 @@ public:
 
   Particle() = default;
 
+  double speed() const;
+
   //! create a secondary particle
   //
   //! stores the current phase space attributes of the particle in the

include/openmc/particle_data.h

@@ -227,9 +227,11 @@ private:
   int g_last_;    //!< pre-collision energy group (MG only)
 
   // Other physical data
-  double wgt_ {1.0};  //!< particle weight
-  double mu_;         //!< angle of scatter
-  bool alive_ {true}; //!< is particle alive?
+  double wgt_ {1.0};       //!< particle weight
+  double mu_;              //!< angle of scatter
+  double time_ {0.0};      //!< time in [s]
+  double time_last_ {0.0}; //!< previous time in [s]
+  bool alive_ {true};      //!< is particle alive?
 
   // Other physical data
   Position r_last_current_; //!< coordinates of the last collision or
@@ -349,6 +351,10 @@ public:
   double& wgt() { return wgt_; }
   double& mu() { return mu_; }
   const double& mu() const { return mu_; }
+  double& time() { return time_; }
+  const double& time() const { return time_; }
+  double& time_last() { return time_last_; }
+  const double& time_last() const { return time_last_; }
   bool& alive() { return alive_; }
 
   Position& r_last_current() { return r_last_current_; }

include/openmc/tallies/filter_time.h

@@ -0,0 +1,50 @@
+#ifndef OPENMC_TALLIES_FILTER_TIME_H
+#define OPENMC_TALLIES_FILTER_TIME_H
+
+#include <gsl/gsl-lite.hpp>
+
+#include "openmc/tallies/filter.h"
+#include "openmc/vector.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins the incident particle time.
+//==============================================================================
+
+class TimeFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~TimeFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type() const override { return "time"; }
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator,
+    FilterMatch& match) const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //----------------------------------------------------------------------------
+  // Accessors
+
+  const vector<double>& bins() const { return bins_; }
+  void set_bins(gsl::span<const double> bins);
+
+protected:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  vector<double> bins_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_ENERGY_H

openmc/filter.py

@@ -23,7 +23,7 @@ _FILTER_TYPES = (
     'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
     'energyfunction', 'cellfrom', 'legendre', 'spatiallegendre',
     'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance',
-    'collision'
+    'collision', 'time'
 )
 
 _CURRENT_NAMES = (
@@ -1308,6 +1308,50 @@ class EnergyoutFilter(EnergyFilter):
 
     """
 
+class TimeFilter(RealFilter):
+    """Bins tally events based on the particle's time.
+
+    Parameters
+    ----------
+    values : iterable of float
+        A list of values for which each successive pair constitutes a range of
+        time in [s] for a single bin
+    filter_id : int
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    values : numpy.ndarray
+        An array of values for which each successive pair constitutes a range of
+        time in [s] for a single bin
+    id : int
+        Unique identifier for the filter
+    bins : numpy.ndarray
+        An array of shape (N, 2) where each row is a pair of time in [s]
+        for a single filter bin
+    num_bins : int
+        The number of filter bins
+
+    """
+    units = 's'
+
+    def get_bin_index(self, filter_bin):
+        # Use lower energy bound to find index for RealFilters
+        deltas = np.abs(self.bins[:, 1] - filter_bin[1]) / filter_bin[1]
+        min_delta = np.min(deltas)
+        if min_delta < 1e-3:
+            return deltas.argmin()
+        else:
+            msg = ('Unable to get the bin index for Filter since '
+                   f'"{filter_bin}" is not one of the bins')
+            raise ValueError(msg)
+
+    def check_bins(self, bins):
+        super().check_bins(bins)
+        for v0, v1 in bins:
+            cv.check_greater_than('filter value', v0, 0., equality=True)
+            cv.check_greater_than('filter value', v1, 0., equality=True)
+
 
 def _path_to_levels(path):
     """Convert distribcell path to list of levels

src/particle.cpp

@@ -36,6 +36,30 @@
 
 namespace openmc {
 
+double Particle::speed() const
+{
+  // Determine mass in eV/c^2
+  double mass;
+  switch (this->type()) {
+  case ParticleType::neutron:
+    mass = MASS_NEUTRON_EV;
+    break;
+  case ParticleType::photon:
+    mass = 0.0;
+    break;
+  case ParticleType::electron:
+  case ParticleType::positron:
+    mass = MASS_ELECTRON_EV;
+    break;
+  }
+
+  // Calculate inverse of Lorentz factor
+  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);
+}
+
 void Particle::create_secondary(
   double wgt, Direction u, double E, ParticleType type)
 {
@@ -99,6 +123,7 @@ void Particle::event_calculate_xs()
   E_last() = E();
   u_last() = u();
   r_last() = r();
+  time_last() = time();
 
   // Reset event variables
   event() = TallyEvent::KILL;
@@ -170,10 +195,11 @@ void Particle::event_advance()
   // Select smaller of the two distances
   double distance = std::min(boundary().distance, collision_distance());
 
-  // Advance particle
+  // Advance particle in space and time
   for (int j = 0; j < n_coord(); ++j) {
     coord(j).r += distance * coord(j).u;
   }
+  this->time() += distance / this->speed();
 
   // Score track-length tallies
   if (!model::active_tracklength_tallies.empty()) {

src/tallies/filter.cpp

@@ -29,6 +29,7 @@
 #include "openmc/tallies/filter_sph_harm.h"
 #include "openmc/tallies/filter_sptl_legendre.h"
 #include "openmc/tallies/filter_surface.h"
+#include "openmc/tallies/filter_time.h"
 #include "openmc/tallies/filter_universe.h"
 #include "openmc/tallies/filter_zernike.h"
 #include "openmc/xml_interface.h"
@@ -151,6 +152,8 @@ Filter* Filter::create(const std::string& type, int32_t id)
     return Filter::create<SpatialLegendreFilter>(id);
   } else if (type == "sphericalharmonics") {
     return Filter::create<SphericalHarmonicsFilter>(id);
+  } else if (type == "time") {
+    return Filter::create<TimeFilter>(id);
   } else if (type == "universe") {
     return Filter::create<UniverseFilter>(id);
   } else if (type == "zernike") {

src/tallies/filter_time.cpp

@@ -0,0 +1,84 @@
+#include "openmc/tallies/filter_time.h"
+
+#include <algorithm> // for min, max
+
+#include <fmt/core.h>
+
+#include "openmc/search.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+//==============================================================================
+// TimeFilter implementation
+//==============================================================================
+
+void TimeFilter::from_xml(pugi::xml_node node)
+{
+  auto bins = get_node_array<double>(node, "bins");
+  this->set_bins(bins);
+}
+
+void TimeFilter::set_bins(gsl::span<const double> bins)
+{
+  // Clear existing 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]);
+  }
+
+  n_bins_ = bins_.size() - 1;
+}
+
+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();
+
+  // If time interval is entirely out of time bin range, exit
+  if (t_end < bins_.front() || t_start >= bins_.back())
+    return;
+
+  // Determine first bin containing a portion of time interval
+  int i_bin = 0;
+  while (bins_[i_bin + 1] < t_start) {
+    ++i_bin;
+  }
+
+  // Find matching bins
+  double dt_total = t_end - t_start;
+  for (; i_bin < bins_.size() - 1; ++i_bin) {
+    double t_left = std::max(t_start, bins_[i_bin]);
+    double t_right = std::min(t_end, bins_[i_bin + 1]);
+
+    // Add match with weight equal to the fraction of the time interval within
+    // the current time bin
+    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;
+  }
+}
+
+void TimeFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "bins", bins_);
+}
+
+std::string TimeFilter::text_label(int bin) const
+{
+  return fmt::format("Time [{}, {})", bins_[bin], bins_[bin + 1]);
+}
+
+} // namespace openmc

src/tallies/tally_scoring.cpp

@@ -628,7 +628,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
         score = flux;
       }
       // Score inverse velocity in units of s/cm.
-      score /= std::sqrt(2. * E / MASS_NEUTRON_EV) * C_LIGHT * 100.;
+      score /= p.speed();
       break;
 
     case SCORE_SCATTER: