From 24b741f2b3ce081f4c51fa15b130c50a363a8fc6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 10 Sep 2021 13:57:20 -0500 Subject: [PATCH] Initial implementation of time filter --- CMakeLists.txt | 1 + docs/source/pythonapi/base.rst | 1 + include/openmc/constants.h | 2 +- include/openmc/particle.h | 2 + include/openmc/particle_data.h | 12 +++- include/openmc/tallies/filter_time.h | 50 +++++++++++++++++ openmc/filter.py | 46 ++++++++++++++- src/particle.cpp | 28 +++++++++- src/tallies/filter.cpp | 3 + src/tallies/filter_time.cpp | 84 ++++++++++++++++++++++++++++ src/tallies/tally_scoring.cpp | 2 +- 11 files changed, 224 insertions(+), 7 deletions(-) create mode 100644 include/openmc/tallies/filter_time.h create mode 100644 src/tallies/filter_time.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index fb5da6f2a..9516971b9 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index 496201ad7..2dea91bb0 100644 --- a/docs/source/pythonapi/base.rst +++ b/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 diff --git a/include/openmc/constants.h b/include/openmc/constants.h index b01238611..2596a0334 100644 --- a/include/openmc/constants.h +++ b/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 diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 952c8c4e7..da2f6a61f 100644 --- a/include/openmc/particle.h +++ b/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 diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 67dabe3bf..733ed3b81 100644 --- a/include/openmc/particle_data.h +++ b/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_; } diff --git a/include/openmc/tallies/filter_time.h b/include/openmc/tallies/filter_time.h new file mode 100644 index 000000000..c66481c57 --- /dev/null +++ b/include/openmc/tallies/filter_time.h @@ -0,0 +1,50 @@ +#ifndef OPENMC_TALLIES_FILTER_TIME_H +#define OPENMC_TALLIES_FILTER_TIME_H + +#include + +#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& 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 95be1774d..699fa1aef 100644 --- a/openmc/filter.py +++ b/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 diff --git a/src/particle.cpp b/src/particle.cpp index fa35f634a..501792e31 100644 --- a/src/particle.cpp +++ b/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()) { diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index e3934d7ea..a1e5c709f 100644 --- a/src/tallies/filter.cpp +++ b/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(id); } else if (type == "sphericalharmonics") { return Filter::create(id); + } else if (type == "time") { + return Filter::create(id); } else if (type == "universe") { return Filter::create(id); } else if (type == "zernike") { diff --git a/src/tallies/filter_time.cpp b/src/tallies/filter_time.cpp new file mode 100644 index 000000000..496c8d2aa --- /dev/null +++ b/src/tallies/filter_time.cpp @@ -0,0 +1,84 @@ +#include "openmc/tallies/filter_time.h" + +#include // for min, max + +#include + +#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(node, "bins"); + this->set_bins(bins); +} + +void TimeFilter::set_bins(gsl::span 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 diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index a0fc79b5b..35d80d0e8 100644 --- a/src/tallies/tally_scoring.cpp +++ b/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: