diff --git a/docs/source/io_formats/particle_restart.rst b/docs/source/io_formats/particle_restart.rst index 6ebf7e7ba..53f86735f 100644 --- a/docs/source/io_formats/particle_restart.rst +++ b/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]. diff --git a/src/particle.cpp b/src/particle.cpp index 7397c318f..b3dd72d9a 100644 --- a/src/particle.cpp +++ b/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 diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index f7159ce44..32d187dd8 100644 --- a/src/particle_restart.cpp +++ b/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); diff --git a/src/tallies/filter_time.cpp b/src/tallies/filter_time.cpp index 225ef3ad3..6d3b9ab78 100644 --- a/src/tallies/filter_time.cpp +++ b/src/tallies/filter_time.cpp @@ -1,6 +1,7 @@ #include "openmc/tallies/filter_time.h" -#include // for min, max +#include // for min, max, copy +#include // for back_inserter #include @@ -25,14 +26,13 @@ void TimeFilter::set_bins(gsl::span 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())) { + 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;