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/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/docs/source/io_formats/source.rst b/docs/source/io_formats/source.rst index 86ae8a1e5..4ce2229ed 100644 --- a/docs/source/io_formats/source.rst +++ b/docs/source/io_formats/source.rst @@ -20,7 +20,7 @@ following the same format. - **source_bank** (Compound type) -- Source bank information for each particle. The compound type has fields ``r``, ``u``, ``E``, - ``wgt``, ``delayed_group``, ``surf_id`` and ``particle``, - which represent the position, direction, energy, weight, + ``time``, ``wgt``, ``delayed_group``, ``surf_id`` and ``particle``, + which represent the position, direction, energy, time, weight, delayed group, surface ID, and particle type (0=neutron, 1=photon, 2=electron, 3=positron), respectively. diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 46725fa8b..591dc2777 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -52,11 +52,11 @@ The current version of the statepoint file format is 17.0. sum-of-squares for each global tally. - **source_bank** (Compound type) -- Source bank information for each particle. The compound type has fields ``r``, ``u``, ``E``, - ``wgt``, ``delayed_group``, ``surf_id``, and ``particle``, which - represent the position, direction, energy, weight, delayed group, - surface ID, and particle type (0=neutron, 1=photon, 2=electron, - 3=positron), respectively. - Only present when `run_mode` is 'eigenvalue'. + ``time``, ``wgt``, ``delayed_group``, ``surf_id``, and + ``particle``, which represent the position, direction, energy, + time, weight, delayed group, surface ID, and particle type + (0=neutron, 1=photon, 2=electron, 3=positron), respectively. Only + present when `run_mode` is 'eigenvalue'. **/tallies/** 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..06fccc092 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -44,9 +44,10 @@ struct SourceSite { Position r; Direction u; double E; - double wgt; - int delayed_group; - int surf_id; + double time {0.0}; + double wgt {1.0}; + int delayed_group {0}; + int surf_id {0}; ParticleType particle; int64_t parent_id; int64_t progeny_id; @@ -227,9 +228,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 +352,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/openmc/lib/core.py b/openmc/lib/core.py index e277d527a..e585444fb 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -16,6 +16,7 @@ class _SourceSite(Structure): _fields_ = [('r', c_double*3), ('u', c_double*3), ('E', c_double), + ('time', c_double), ('wgt', c_double), ('delayed_group', c_int), ('surf_id', c_int), diff --git a/openmc/source.py b/openmc/source.py index 415aacd96..8f25bcab0 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -253,6 +253,8 @@ class SourceParticle: Directional cosines E : float Energy of particle in [eV] + time : float + Time of particle in [s] wgt : float Weight of the particle delayed_group : int @@ -263,11 +265,12 @@ class SourceParticle: Type of the particle """ - def __init__(self, r=(0., 0., 0.), u=(0., 0., 1.), E=1.0e6, wgt=1.0, + def __init__(self, r=(0., 0., 0.), u=(0., 0., 1.), E=1.0e6, time=0.0, wgt=1.0, delayed_group=0, surf_id=0, particle=ParticleType.NEUTRON): self.r = tuple(r) self.u = tuple(u) self.E = float(E) + self.time = float(time) self.wgt = float(wgt) self.delayed_group = delayed_group self.surf_id = surf_id @@ -282,7 +285,7 @@ class SourceParticle: Source particle attributes """ - return (self.r, self.u, self.E, self.wgt, + return (self.r, self.u, self.E, self.time, self.wgt, self.delayed_group, self.surf_id, self.particle.value) @@ -309,6 +312,7 @@ def write_source_file(source_particles, filename, **kwargs): ('r', pos_dtype), ('u', pos_dtype), ('E', '= 0; --i) { + MPI_Get_address(&b.time, &disp[3]); + MPI_Get_address(&b.wgt, &disp[4]); + MPI_Get_address(&b.delayed_group, &disp[5]); + MPI_Get_address(&b.surf_id, &disp[6]); + MPI_Get_address(&b.particle, &disp[7]); + MPI_Get_address(&b.parent_id, &disp[8]); + MPI_Get_address(&b.progeny_id, &disp[9]); + for (int i = 9; i >= 0; --i) { disp[i] -= disp[0]; } - int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1}; - MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT, - MPI_INT, MPI_INT, MPI_LONG, MPI_LONG}; - MPI_Type_create_struct(9, blocks, disp, types, &mpi::source_site); + int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1}; + MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, + MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG}; + MPI_Type_create_struct(10, blocks, disp, types, &mpi::source_site); MPI_Type_commit(&mpi::source_site); } #endif // OPENMC_MPI diff --git a/src/particle.cpp b/src/particle.cpp index fa35f634a..b3dd72d9a 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 + 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); +} + void Particle::create_secondary( double wgt, Direction u, double E, ParticleType type) { @@ -53,6 +77,7 @@ void Particle::create_secondary( bank.r = r(); bank.u = u; bank.E = settings::run_CE ? E : g(); + bank.time = time(); n_bank_second() += 1; } @@ -87,6 +112,8 @@ void Particle::from_source(const SourceSite* src) E() = data::mg.energy_bin_avg_[g()]; } E_last() = E(); + time() = src->time; + time_last() = src->time; } void Particle::event_calculate_xs() @@ -99,6 +126,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 +198,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()) { @@ -373,6 +402,7 @@ void Particle::cross_surface() site.r = r(); site.u = u(); site.E = E(); + site.time = time(); site.wgt = wgt(); site.delayed_group = delayed_group(); site.surf_id = surf->id_; @@ -655,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) * @@ -667,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/physics.cpp b/src/physics.cpp index 3a52396e1..36161f3ea 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -194,6 +194,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) SourceSite site; site.r = p.r(); site.particle = ParticleType::neutron; + site.time = p.time(); site.wgt = 1. / weight; site.parent_id = p.id(); site.progeny_id = p.n_progeny()++; diff --git a/src/source.cpp b/src/source.cpp index 9c897dcf3..be4270b51 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -147,9 +147,6 @@ SourceSite IndependentSource::sample(uint64_t* seed) const { SourceSite site; - // Set weight to one by default - site.wgt = 1.0; - // Repeat sampling source location until a good site has been found bool found = false; int n_reject = 0; @@ -227,11 +224,6 @@ SourceSite IndependentSource::sample(uint64_t* seed) const break; } - // Set delayed group - site.delayed_group = 0; - // Set surface ID - site.surf_id = 0; - return site; } 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..3f78d1817 --- /dev/null +++ b/src/tallies/filter_time.cpp @@ -0,0 +1,101 @@ +#include "openmc/tallies/filter_time.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 { + +//============================================================================== +// 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()); + + // Ensure time bins are sorted and don't have duplicates + if (std::adjacent_find(bins.cbegin(), bins.cend(), std::greater_equal<>()) != + bins.end()) { + 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; +} + +void TimeFilter::get_all_bins( + const Particle& p, TallyEstimator estimator, FilterMatch& match) const +{ + // Get the start/end time of the particle for this track + 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()) + return; + + if (estimator == TallyEstimator::TRACKLENGTH) { + // ------------------------------------------------------------------------- + // For tracklength estimator, 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); + + // 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 + 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; + } + } else { + // ------------------------------------------------------------------------- + // For collision estimator or surface tallies, find a match based on the + // exact time of the particle + 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); + } +} + +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.cpp b/src/tallies/tally.cpp index 8e8e3126f..2d673d7bc 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -485,6 +485,7 @@ void Tally::set_scores(const vector& scores) fatal_error("Cannot tally mesh surface currents in the same tally as " "normal surface currents"); type_ = TallyType::SURFACE; + estimator_ = TallyEstimator::ANALOG; } else if (meshsurface_present) { type_ = TallyType::MESH_SURFACE; } else { 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: diff --git a/tests/regression_tests/mgxs_library_condense/results_true.dat b/tests/regression_tests/mgxs_library_condense/results_true.dat index 349d9debe..508984e54 100644 --- a/tests/regression_tests/mgxs_library_condense/results_true.dat +++ b/tests/regression_tests/mgxs_library_condense/results_true.dat @@ -162,10 +162,10 @@ 3 2 2 1 1 total 1.0 0.150417 mesh 1 group in nuclide mean std. dev. x y z -0 1 1 1 1 total 4.800872e-07 2.910230e-08 -2 1 2 1 1 total 5.007997e-07 3.613803e-08 -1 2 1 1 1 total 4.872635e-07 4.324235e-08 -3 2 2 1 1 total 4.890451e-07 2.833017e-08 +0 1 1 1 1 total 4.800874e-07 2.910227e-08 +2 1 2 1 1 total 5.008000e-07 3.613802e-08 +1 2 1 1 1 total 4.872638e-07 4.324238e-08 +3 2 2 1 1 total 4.890453e-07 2.833017e-08 mesh 1 group in nuclide mean std. dev. x y z 0 1 1 1 1 total 0.024194 0.001614 diff --git a/tests/regression_tests/mgxs_library_distribcell/results_true.dat b/tests/regression_tests/mgxs_library_distribcell/results_true.dat index 9fff8479d..e6e4b99e4 100644 --- a/tests/regression_tests/mgxs_library_distribcell/results_true.dat +++ b/tests/regression_tests/mgxs_library_distribcell/results_true.dat @@ -49,7 +49,7 @@ sum(distribcell) group out nuclide mean std. dev. 0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, ...),) 1 total 1.0 0.090571 sum(distribcell) group in nuclide mean std. dev. -0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, ...),) 1 total 4.896403e-07 2.047455e-08 +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, ...),) 1 total 4.896406e-07 2.047455e-08 sum(distribcell) group in nuclide mean std. dev. 0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, ...),) 1 total 0.088451 0.003512 sum(distribcell) group in group out nuclide mean std. dev. diff --git a/tests/regression_tests/mgxs_library_hdf5/results_true.dat b/tests/regression_tests/mgxs_library_hdf5/results_true.dat index f49ae97ab..6da73ec9b 100644 --- a/tests/regression_tests/mgxs_library_hdf5/results_true.dat +++ b/tests/regression_tests/mgxs_library_hdf5/results_true.dat @@ -94,8 +94,8 @@ domain=1 type=chi-prompt [1.00000000e+00 0.00000000e+00] [8.07424611e-02 0.00000000e+00] domain=1 type=inverse-velocity -[5.93306434e-08 2.99151552e-06] -[3.31163681e-09 2.75470068e-07] +[5.93309775e-08 2.99151502e-06] +[3.31163250e-09 2.75469881e-07] domain=1 type=prompt-nu-fission [6.07803967e-03 1.32325631e-01] [1.12229103e-04 1.39042100e-02] diff --git a/tests/regression_tests/mgxs_library_mesh/results_true.dat b/tests/regression_tests/mgxs_library_mesh/results_true.dat index 0acc2f407..a8775984f 100644 --- a/tests/regression_tests/mgxs_library_mesh/results_true.dat +++ b/tests/regression_tests/mgxs_library_mesh/results_true.dat @@ -162,10 +162,10 @@ 3 2 2 1 1 total 1.0 0.052287 mesh 1 group in nuclide mean std. dev. x y z -0 1 1 1 1 total 8.523138e-10 3.481466e-11 -2 1 2 1 1 total 9.044552e-10 2.851430e-11 -1 2 1 1 1 total 8.643697e-10 4.046801e-11 -3 2 2 1 1 total 8.667214e-10 1.753703e-11 +0 1 1 1 1 total 8.529639e-10 3.484230e-11 +2 1 2 1 1 total 9.050800e-10 2.850752e-11 +1 2 1 1 1 total 8.650057e-10 4.049297e-11 +3 2 2 1 1 total 8.673543e-10 1.753754e-11 mesh 1 group in nuclide mean std. dev. x y z 0 1 1 1 1 total 0.030356 0.001982 diff --git a/tests/regression_tests/mgxs_library_no_nuclides/results_true.dat b/tests/regression_tests/mgxs_library_no_nuclides/results_true.dat index 5afe3b24c..32ae7f4cf 100644 --- a/tests/regression_tests/mgxs_library_no_nuclides/results_true.dat +++ b/tests/regression_tests/mgxs_library_no_nuclides/results_true.dat @@ -138,8 +138,8 @@ chi-prompt 0 1 2 total 0.0 0.000000 inverse-velocity material group in nuclide mean std. dev. -1 1 1 total 6.148626e-08 3.647428e-09 -0 1 2 total 2.844242e-06 3.098231e-07 +1 1 1 total 6.148661e-08 3.647428e-09 +0 1 2 total 2.844242e-06 3.098230e-07 prompt-nu-fission material group in nuclide mean std. dev. 1 1 1 total 0.019154 0.000623 @@ -440,8 +440,8 @@ chi-prompt 0 2 2 total 0.0 0.0 inverse-velocity material group in nuclide mean std. dev. -1 2 1 total 6.397171e-08 3.117086e-09 -0 2 2 total 2.875630e-06 2.751730e-07 +1 2 1 total 6.397205e-08 3.117085e-09 +0 2 2 total 2.875630e-06 2.751728e-07 prompt-nu-fission material group in nuclide mean std. dev. 1 2 1 total 0.0 0.0 @@ -742,8 +742,8 @@ chi-prompt 0 3 2 total 0.0 0.0 inverse-velocity material group in nuclide mean std. dev. -1 3 1 total 6.447396e-08 2.866174e-09 -0 3 2 total 3.006811e-06 1.808751e-07 +1 3 1 total 6.447429e-08 2.866173e-09 +0 3 2 total 3.006810e-06 1.808750e-07 prompt-nu-fission material group in nuclide mean std. dev. 1 3 1 total 0.0 0.0 diff --git a/tests/regression_tests/mgxs_library_nuclides/results_true.dat b/tests/regression_tests/mgxs_library_nuclides/results_true.dat index fc005b133..c0084547b 100644 --- a/tests/regression_tests/mgxs_library_nuclides/results_true.dat +++ b/tests/regression_tests/mgxs_library_nuclides/results_true.dat @@ -1 +1 @@ -b8706c9586aeeb5ee558829c28772f7f6e18a0d3fe4e30a2059b6e74564d1ebbbd0ab5473965c4e6156dff8b8178ce76c74d8db1d63f032c3d4da5f7eb9eae2c \ No newline at end of file +d36f8abb1131212063470d622dbedeae31602da71006389421bd5dba712c94fe96acbc8ded833cb237687fb1071c1a6c3e0ec67b18ce7cf0f7720451b86d993a \ No newline at end of file diff --git a/tests/regression_tests/surface_source/surface_source_true.h5 b/tests/regression_tests/surface_source/surface_source_true.h5 index 22edb898b..2ea48b2ac 100644 Binary files a/tests/regression_tests/surface_source/surface_source_true.h5 and b/tests/regression_tests/surface_source/surface_source_true.h5 differ diff --git a/tests/regression_tests/tallies/results_true.dat b/tests/regression_tests/tallies/results_true.dat index f6403191d..c73ff33a4 100644 --- a/tests/regression_tests/tallies/results_true.dat +++ b/tests/regression_tests/tallies/results_true.dat @@ -1 +1 @@ -22cd0d5d2764367e2a3ca2b715d0863ac7afd33cbd4f0b2707722e6267d95d3bc07026c9a7fdbc8260c69c40a6617b6e1c831de5f9c4a7cb89c4efa6b74bea76 \ No newline at end of file +de09325b517aad9f58940e5fcd53b005c57ea008a15699769ab91aba723de141014101e55be2621e34c369beefb6db8c6cf847e241ee26e73e394c0f155138b0 \ No newline at end of file diff --git a/tests/unit_tests/test_time_filter.py b/tests/unit_tests/test_time_filter.py new file mode 100644 index 000000000..5b09c6d5f --- /dev/null +++ b/tests/unit_tests/test_time_filter.py @@ -0,0 +1,151 @@ +from math import sqrt +from random import uniform + +import numpy as np +import openmc +import pytest + + +def test_time_filter_basics(): + f = openmc.TimeFilter([1.0, 2.0, 5.0, 9.0]) + np.testing.assert_allclose(f.bins[0], [1.0, 2.0]) + np.testing.assert_allclose(f.bins[1], [2.0, 5.0]) + np.testing.assert_allclose(f.bins[-1], [5.0, 9.0]) + assert len(f.bins) == 3 + + # Make sure __repr__ works + repr(f) + + # to_xml_element() + elem = f.to_xml_element() + assert elem.tag == 'filter' + assert elem.attrib['type'] == 'time' + + +def time(particle, distance, E): + """Return the time it takes a particle at a given energy to travel a certain + distance""" + if particle == 'neutron': + mass = 939.56542052e6 # eV/c² + elif particle == 'photon': + mass = 0.0 + + # Calculate speed via v = c * sqrt(1 - γ^-2) + inv_gamma = mass / (E + mass) + velocity = 2.99792458e10 * sqrt(1 - inv_gamma * inv_gamma) # cm/s + return distance / velocity + + +@pytest.fixture(params=['neutron', 'photon']) +def model(request): + # Select random sphere radius, source position, and source energy + r = uniform(0., 2.) + x = uniform(2., 10.) + E = uniform(0., 20.0e6) + + # Create model + model = openmc.Model() + mat = openmc.Material() + mat.add_nuclide('Zr90', 1.0) + mat.set_density('g/cm3', 1.0) + model.materials.append(mat) + inner_sphere = openmc.Sphere(r=r) + outer_sphere = openmc.Sphere(r=10.0, boundary_type='vacuum') + center_cell = openmc.Cell(fill=mat, region=-inner_sphere) + outer_void = openmc.Cell(region=+inner_sphere & -outer_sphere) + model.geometry = openmc.Geometry([center_cell, outer_void]) + model.settings = openmc.Settings() + model.settings.run_mode = 'fixed source' + model.settings.particles = 1000 + model.settings.batches = 20 + particle = request.param + model.settings.source = openmc.Source( + space=openmc.stats.Point((x, 0., 0.)), + angle=openmc.stats.Monodirectional([-1., 0., 0.]), + energy=openmc.stats.Discrete([E], [1.0]), + particle=particle + ) + + # Calculate time it will take neutrons to reach sphere + t0 = time(particle, x - r, E) + + # Create tally with time filter + tally = openmc.Tally() + tally.filters = [openmc.TimeFilter([0.0, t0, 2*t0])] + tally.scores = ['total'] + model.tallies.append(tally) + return model + + +@pytest.fixture(params=['neutron', 'photon']) +def model_surf(request): + # Select random distance and source energy + x = uniform(50., 100.) + E = uniform(0., 20.0e6) + + # Create model + model = openmc.Model() + mat = openmc.Material() + mat.add_nuclide('Zr90', 1.0) + mat.set_density('g/cm3', 1.0) + model.materials.append(mat) + left = openmc.XPlane(-1., boundary_type='vacuum') + black_surface = openmc.XPlane(x, boundary_type='vacuum') + right = openmc.XPlane(x + 1) + void_cell = openmc.Cell(region=+left & -black_surface) + black_cell = openmc.Cell(region=+black_surface & -right) + model.geometry = openmc.Geometry([void_cell, black_cell]) + model.settings = openmc.Settings() + model.settings.run_mode = 'fixed source' + model.settings.particles = 1000 + model.settings.batches = 20 + particle = request.param + model.settings.source = openmc.Source( + space=openmc.stats.Point((0., 0., 0.)), + angle=openmc.stats.Monodirectional([1., 0., 0.]), + energy=openmc.stats.Discrete([E], [1.0]), + particle=particle + ) + + # Calculate time it will take neutrons to reach purely-absorbing surface + t0 = time(particle, x, E) + + # Create tally with surface and time filters + tally = openmc.Tally() + tally.filters = [ + openmc.SurfaceFilter([black_surface]), + openmc.TimeFilter([0.0, t0*0.999, t0*1.001, 100.0]) + ] + tally.scores = ['current'] + model.tallies.append(tally) + return model + + +def test_time_filter_volume(model, run_in_tmpdir): + sp_filename = model.run() + with openmc.StatePoint(sp_filename) as sp: + t = sp.tallies[model.tallies[0].id] + values = t.mean.ravel() + + # Before t0, the reaction rate should be zero + assert values[0] == pytest.approx(0.0) + + # After t0, the reaction rate should be positive + assert values[1] > 0.0 + + +def test_time_filter_surface(model_surf, run_in_tmpdir): + sp_filename = model_surf.run() + with openmc.StatePoint(sp_filename) as sp: + t = sp.tallies[model_surf.tallies[0].id] + values = t.mean.ravel() + print(values) + + # Before t0-ε, the current should be zero + assert values[0] == 0.0 + + # Between t0-ε and t0+ε, the current should be one + assert values[1] == 1.0 + + # After t0+ε, the current should be zero + assert values[2] == 0.0