diff --git a/CMakeLists.txt b/CMakeLists.txt index 94caec757ab..84c64777fe3 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index 210ab284ba9..c558f53cc90 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -24,6 +24,7 @@ enum class FilterType { CELL, CELL_INSTANCE, COLLISION, + D, DELAYED_GROUP, DISTRIBCELL, ENERGY_FUNCTION, diff --git a/include/openmc/tallies/filter_d.h b/include/openmc/tallies/filter_d.h new file mode 100644 index 00000000000..d10e7719d23 --- /dev/null +++ b/include/openmc/tallies/filter_d.h @@ -0,0 +1,51 @@ +#ifndef OPENMC_TALLIES_FILTER_D_H +#define OPENMC_TALLIES_FILTER_D_H + +#include + +#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& bins() const { return bins_; } + void set_bins(gsl::span bins); + +protected: + //---------------------------------------------------------------------------- + // Data members + + vector bins_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_ENERGY_H diff --git a/openmc/filter.py b/openmc/filter.py index b5926af5ad9..334e5c8aaed 100644 --- a/openmc/filter.py +++ b/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. diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index 074212db44d..a691a37566a 100644 --- a/src/tallies/filter.cpp +++ b/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(id); } else if (type == "time") { return Filter::create(id); + } else if (type == "d") { + return Filter::create(id); } else if (type == "universe") { return Filter::create(id); } else if (type == "zernike") { diff --git a/src/tallies/filter_d.cpp b/src/tallies/filter_d.cpp new file mode 100644 index 00000000000..be5e9ea40b0 --- /dev/null +++ b/src/tallies/filter_d.cpp @@ -0,0 +1,70 @@ +#include "openmc/tallies/filter_d.h" + +#include // for min, max, copy, adjacent_find +#include // for greater_equal +#include // for back_inserter + +#include + +#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(node, "bins"); + this->set_bins(bins); +} + +// Set bin values for the filter +void DFilter::set_bins(gsl::span 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