D filter

CMakeLists.txt

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

include/openmc/tallies/filter.h

@@ -24,6 +24,7 @@ enum class FilterType {
   CELL,
   CELL_INSTANCE,
   COLLISION,
+  D,
   DELAYED_GROUP,
   DISTRIBCELL,
   ENERGY_FUNCTION,

include/openmc/tallies/filter_d.h

@@ -0,0 +1,51 @@
+#ifndef OPENMC_TALLIES_FILTER_D_H
+#define OPENMC_TALLIES_FILTER_D_H
+
+#include <gsl/gsl-lite.hpp>
+
+#include "openmc/tallies/filter.h"
+#include "openmc/vector.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins the incident particle time.
+//==============================================================================
+
+class DFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~DFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "d"; }
+  FilterType type() const override { return FilterType::D; }
+
+  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

@@ -25,7 +25,7 @@
     'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
     'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre',
     'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance',
-    'collision', 'time'
+    'collision', 'time' , "d"
 )
 
 _CURRENT_NAMES = (
@@ -1569,6 +1569,52 @@ def _path_to_levels(path):
     return path_items
 
 
+class DFilter(RealFilter):
+    """Bins tally events based on the particle's time multiplied by velocity.
+
+    Parameters
+    ----------
+    values : iterable of float
+        A list of values for which each successive pair constitutes a range of
+        time × velocity in [cm] for a single bin.
+    filter_id : int, optional
+        Unique identifier for the filter.
+
+    Attributes
+    ----------
+    values : numpy.ndarray
+        An array of values for which each successive pair constitutes a range of
+        time × velocity in [m] 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 × velocity in [m]
+        for a single filter bin.
+    num_bins : int
+        The number of filter bins.
+
+    """
+    units = 'cm'  # Modified time units (since time * velocity = distance)
+
+    def get_bin_index(self, filter_bin):
+        """Finds the corresponding bin index for a given time × velocity."""
+        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):
+        """Check that bins are valid and strictly increasing."""
+        super().check_bins(bins)
+        for v0, v1 in bins:
+            cv.check_greater_than('filter value', v0, 0., equality=False)
+            cv.check_greater_than('filter value', v1, 0., equality=False)
+
+
 class DistribcellFilter(Filter):
     """Bins tally event locations on instances of repeated cells.
 

src/tallies/filter.cpp

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

src/tallies/filter_d.cpp

@@ -0,0 +1,70 @@
+#include "openmc/tallies/filter_d.h"
+
+#include <algorithm>  // for min, max, copy, adjacent_find
+#include <functional> // for greater_equal
+#include <iterator>   // for back_inserter
+
+#include <fmt/core.h>
+
+#include "openmc/search.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+//==============================================================================
+// DFilter implementation
+//==============================================================================
+
+// Load filter parameters from an XML input file
+void DFilter::from_xml(pugi::xml_node node) {
+  auto bins = get_node_array<double>(node, "bins");
+  this->set_bins(bins);
+}
+
+// Set bin values for the filter
+void DFilter::set_bins(gsl::span<const double> bins) {
+  // Clear existing bins
+  bins_.clear();
+  bins_.reserve(bins.size());
+
+  // Ensure bins are sorted and unique
+  if (std::adjacent_find(bins.cbegin(), bins.cend(), std::greater_equal<>()) != bins.end()) {
+    throw std::runtime_error {"DFilter bins must be monotonically increasing."};
+  }
+
+  // Copy bins
+  std::copy(bins.cbegin(), bins.cend(), std::back_inserter(bins_));
+  n_bins_ = bins_.size() - 1;
+}
+
+// Determine which bin the particle belongs to
+void DFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const {
+  // Get particle properties
+  double particle_time = p.time();
+  double particle_velocity = p.speed(); // Assuming OpenMC provides a `speed()` function
+
+  // Compute the modified time
+  double modified_time = particle_time * particle_velocity;
+  // double modified_time = 25.12;
+
+  // If modified time is outside bin range, exit
+  if (modified_time < bins_.front() || modified_time >= bins_.back()) return;
+
+  // Find the appropriate bin
+  auto i_bin = lower_bound_index(bins_.begin(), bins_.end(), modified_time);
+  match.bins_.push_back(i_bin);
+  match.weights_.push_back(1.0);
+}
+
+// Write filter state to a statepoint file
+void DFilter::to_statepoint(hid_t filter_group) const {
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "bins", bins_);
+}
+
+// Provide a text label for the filter (used in output files)
+std::string DFilter::text_label(int bin) const {
+  return fmt::format("D Filter [{}, {})", bins_[bin], bins_[bin + 1]);
+}
+
+} // namespace openmc