From 24ed0180afecd7c6271bd0c321bc58055025e8ef Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 2 Aug 2024 10:50:08 +0700 Subject: [PATCH 01/17] add TimedMeshFilter Python API --- openmc/filter.py | 155 ++++++++++++++++++++++++++++++++++++++++++- openmc/lib/filter.py | 2 + 2 files changed, 155 insertions(+), 2 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index 0f884cfcd16..5e065d351b8 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -25,7 +25,7 @@ 'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre', 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance', - 'collision', 'time' + 'collision', 'time', 'timedmesh' ) _CURRENT_NAMES = ( @@ -765,7 +765,7 @@ def __eq__(self, other): @Filter.bins.setter def bins(self, bins): cv.check_type('bins', bins, Sequence, str) - bins = np.atleast_1d(bins) + bins = np.atleast_1d(bins) for edge in bins: cv.check_value('filter bin', edge, _PARTICLES) self._bins = bins @@ -1119,6 +1119,157 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs): return pd.concat([df, pd.DataFrame(filter_dict)]) +class TimedMeshFilter(Filter): + """Bins tally event locations by mesh elements and time grids. + + Parameters + ---------- + mesh : openmc.MeshBase + The mesh object that events will be tallied onto + time_grid : iterable of float + A list of values for which each successive pair constitutes a range of + time in [s] for a single time bin + filter_id : int + Unique identifier for the filter + + Attributes + ---------- + mesh : openmc.MeshBase + The mesh object that events will be tallied onto + time_grid : numpy.ndarray + An array of values for which each successive pair constitutes a range of + time in [s] for a single time bin + id : int + Unique identifier for the filter + translation : Iterable of float + This array specifies a vector that is used to translate (shift) the mesh + for this filter + bins : list of tuple + A list of mesh indices for each filter bin, e.g. [(1, 1, 1), (2, 1, 1), + ...] + num_bins : Integral + The number of filter bins + + """ + + def __init__(self, mesh, time_grid, filter_id=None): + self._time_grid = np.array([0.0, np.inf]) # Temporary + self.mesh = mesh + self.time_grid = time_grid + self.id = filter_id + self._translation = None + + def __hash__(self): + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tMesh ID', self.mesh.id) + string += '{: <16}=\t{}\n'.format('\tTime grid', self.time_grid) + return hash(string) + + def __repr__(self): + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tMesh ID', self.mesh.id) + string += '{: <16}=\t{}\n'.format('\tTime grid', self.time_grid) + string += '{: <16}=\t{}\n'.format('\tID', self.id) + string += '{: <16}=\t{}\n'.format('\tTranslation', self.translation) + return string + + @classmethod + def from_hdf5(cls, group, **kwargs): + if group['type'][()].decode() != cls.short_name.lower(): + raise ValueError("Expected HDF5 data for filter type '" + + cls.short_name.lower() + "' but got '" + + group['type'][()].decode() + " instead") + + if 'meshes' not in kwargs: + raise ValueError(cls.__name__ + " requires a 'meshes' keyword " + "argument.") + + mesh_id = group['mesh_bins'][()] + mesh_obj = kwargs['meshes'][mesh_id] + filter_id = int(group.name.split('/')[-1].lstrip('filter ')) + + time_grid = group['time_bins'][()] + + out = cls(mesh_obj, time_grid, filter_id=filter_id) + + translation = group.get('translation') + if translation: + out.translation = translation[()] + + return out + + @property + def mesh(self): + return self._mesh + + @property + def time_grid(self): + return self._time_grid + + @mesh.setter + def mesh(self, mesh): + cv.check_type('filter mesh', mesh, openmc.MeshBase) + self._mesh = mesh + self._reset_bins() + + @time_grid.setter + def time_grid(self, time_grid): + self._time_grid = time_grid + self._reset_bins() + + def _reset_bins(self): + mesh = self._mesh + time_grid = self._time_grid + + # Mesh bins + if isinstance(mesh, openmc.UnstructuredMesh): + if mesh.has_statepoint_data: + mesh_bins = list(range(len(mesh.volumes))) + else: + mesh_bins = [] + else: + mesh_bins = list(mesh.indices) + + # Reset bins + self.bins = [] + for i in range(len(self._time_grid) - 1): + for j in range(len(mesh_bins)): + new_bin = () + new_bin += ([float(time_grid[i]), float(time_grid[i+1])],) + new_bin += (mesh_bins[j],) + self.bins.append(new_bin) + + @property + def shape(self): + return (self.time_grid - 1,) + self.mesh.dimension + + @property + def translation(self): + return self._translation + + @translation.setter + def translation(self, t): + cv.check_type('mesh filter translation', t, Iterable, Real) + cv.check_length('mesh filter translation', t, 3) + self._translation = np.asarray(t) + + def can_merge(self, other): + return False + + def to_xml_element(self): + element = super().to_xml_element() + element[0].text = None + + time_bins = ET.SubElement(element[0], 'time_bins') + mesh_bins = ET.SubElement(element[0], 'mesh_bins') + + time_bins.text = ' '.join(str(x) for x in self.time_grid) + mesh_bins.text = str(self.mesh.id) + if self.translation is not None: + element.set('translation', ' '.join(map(str, self.translation))) + return element + + class CollisionFilter(Filter): """Bins tally events based on the number of collisions. diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index 340c2fa3448..37d759fbb66 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -647,6 +647,8 @@ class ZernikeRadialFilter(ZernikeFilter): 'sphericalharmonics': SphericalHarmonicsFilter, 'spatiallegendre': SpatialLegendreFilter, 'surface': SurfaceFilter, + 'time': TimeFilter, + 'timedmesh': TimedMeshFilter, 'universe': UniverseFilter, 'zernike': ZernikeFilter, 'zernikeradial': ZernikeRadialFilter From 75aa7989f659c9817c5deb1cd2bb5b3432d23c29 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 2 Aug 2024 10:56:17 +0700 Subject: [PATCH 02/17] autopep8 --- openmc/filter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/filter.py b/openmc/filter.py index 5e065d351b8..499d2ab612d 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1153,7 +1153,7 @@ class TimedMeshFilter(Filter): """ def __init__(self, mesh, time_grid, filter_id=None): - self._time_grid = np.array([0.0, np.inf]) # Temporary + self._time_grid = np.array([0.0, np.inf]) # Temporary self.mesh = mesh self.time_grid = time_grid self.id = filter_id From 319d2a111889d51e0518238be73a51a7ac8458f9 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 2 Aug 2024 11:15:52 +0700 Subject: [PATCH 03/17] add timed_mesh_filter.cpp --- CMakeLists.txt | 1 + include/openmc/tallies/filter.h | 1 + include/openmc/tallies/filter_timed_mesh.h | 74 +++++++ src/tallies/filter.cpp | 3 + src/tallies/filter_timed_mesh.cpp | 230 +++++++++++++++++++++ src/tallies/tally.cpp | 9 + 6 files changed, 318 insertions(+) create mode 100644 include/openmc/tallies/filter_timed_mesh.h create mode 100644 src/tallies/filter_timed_mesh.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 0f4cc1b527b..1a780c9f254 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -422,6 +422,7 @@ list(APPEND libopenmc_SOURCES src/tallies/filter_sptl_legendre.cpp src/tallies/filter_surface.cpp src/tallies/filter_time.cpp + src/tallies/filter_timed_mesh.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 1166c0eee53..de8c62d302b 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -42,6 +42,7 @@ enum class FilterType { SPATIAL_LEGENDRE, SURFACE, TIME, + TIMED_MESH, UNIVERSE, ZERNIKE, ZERNIKE_RADIAL diff --git a/include/openmc/tallies/filter_timed_mesh.h b/include/openmc/tallies/filter_timed_mesh.h new file mode 100644 index 00000000000..a73ed9a759f --- /dev/null +++ b/include/openmc/tallies/filter_timed_mesh.h @@ -0,0 +1,74 @@ +#ifndef OPENMC_TALLIES_FILTER_TIMED_MESH_H +#define OPENMC_TALLIES_FILTER_TIMED_MESH_H + +#include + +#include "openmc/constants.h" +#include "openmc/position.h" +#include "openmc/tallies/filter.h" + +namespace openmc { + +//============================================================================== +//! Indexes the time-location of particle events to a time gridn and a regular +// mesh. For tracklength tallies, it will produce multiple valid bins and the +// bin weight will correspond to the fraction of the track length that lies in +//! that bin. +//============================================================================== + +class TimedMeshFilter : public Filter { +public: + //---------------------------------------------------------------------------- + // Constructors, destructors + + ~TimedMeshFilter() = default; + + //---------------------------------------------------------------------------- + // Methods + + std::string type_str() const override { return "timedmesh"; } + FilterType type() const override { return FilterType::TIMED_MESH; } + + 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 + + int32_t mesh() const { return mesh_; } + + void set_mesh(int32_t mesh); + + void set_translation(const Position& translation); + + void set_translation(const double translation[3]); + + const Position& translation() const { return translation_; } + + bool translated() const { return translated_; } + + const vector& time_grid() const { return time_grid_; } + + void set_time_grid(gsl::span time_grid); + + void reset_bins(); + +protected: + //---------------------------------------------------------------------------- + // Data members + + int32_t mesh_; //!< Index of the mesh + int mesh_n_bins_; + bool translated_ {false}; //!< Whether or not the filter is translated + Position translation_ {0.0, 0.0, 0.0}; //!< Filter translation + vector time_grid_ {0.0, INFTY}; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_TIMED_MESH_H diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index ff7a3416b90..f1496ddb74b 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -32,6 +32,7 @@ #include "openmc/tallies/filter_sptl_legendre.h" #include "openmc/tallies/filter_surface.h" #include "openmc/tallies/filter_time.h" +#include "openmc/tallies/filter_timed_mesh.h" #include "openmc/tallies/filter_universe.h" #include "openmc/tallies/filter_zernike.h" #include "openmc/xml_interface.h" @@ -145,6 +146,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 == "timedmesh") { + return Filter::create(id); } else if (type == "universe") { return Filter::create(id); } else if (type == "zernike") { diff --git a/src/tallies/filter_timed_mesh.cpp b/src/tallies/filter_timed_mesh.cpp new file mode 100644 index 00000000000..38a6462c805 --- /dev/null +++ b/src/tallies/filter_timed_mesh.cpp @@ -0,0 +1,230 @@ +#include "openmc/tallies/filter_timed_mesh.h" + +#include +#include + +#include "openmc/capi.h" +#include "openmc/constants.h" +#include "openmc/error.h" +#include "openmc/mesh.h" +#include "openmc/search.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void TimedMeshFilter::from_xml(pugi::xml_node node) +{ + pugi::xml_node bin_node = node.child("bins"); + + //---------------------------------------------------------------------------- + // Mesh + + auto mesh_bins_ = get_node_array(bin_node, "mesh_bins"); + if (mesh_bins_.size() != 1) { + fatal_error( + "Only one mesh can be specified per " + type_str() + " mesh filter."); + } + + auto id = mesh_bins_[0]; + auto search = model::mesh_map.find(id); + if (search != model::mesh_map.end()) { + set_mesh(search->second); + } else { + fatal_error( + fmt::format("Could not find mesh {} specified on tally filter.", id)); + } + + if (check_for_node(node, "translation")) { + set_translation(get_node_array(node, "translation")); + } + + //---------------------------------------------------------------------------- + // Time bins + + auto time_grid = get_node_array(bin_node, "time_bins"); + this->set_time_grid(time_grid); +} + +void TimedMeshFilter::set_mesh(int32_t mesh) +{ + // perform any additional perparation for mesh tallies here + mesh_ = mesh; + model::meshes[mesh_]->prepare_for_tallies(); + + reset_bins(); +} + +void TimedMeshFilter::set_time_grid(gsl::span time_grid) +{ + // Clear existing bins + time_grid_.clear(); + time_grid_.reserve(time_grid.size()); + + // Ensure time grid is sorted and don't have duplicates + if (std::adjacent_find(time_grid.cbegin(), time_grid.cend(), std::greater_equal<>()) != + time_grid.end()) { + throw std::runtime_error {"Time grid must be monotonically increasing."}; + } + + // Copy grid + std::copy(time_grid.cbegin(), time_grid.cend(), std::back_inserter(time_grid_)); + + reset_bins(); +} + +void TimedMeshFilter::reset_bins() +{ + mesh_n_bins_ = model::meshes[mesh_]->n_bins(); + n_bins_ = (time_grid_.size() - 1) * mesh_n_bins_; +} + +void TimedMeshFilter::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 < time_grid_.front() || t_start >= time_grid_.back()) + return; + + // Get the start/end positions, direction, and speed + Position last_r = p.r_last(); + Position r = p.r(); + Position u = p.u(); + const auto speed = p.speed(); + + // apply translation if present + if (translated_) { + last_r -= translation(); + r -= translation(); + } + + if (estimator != TallyEstimator::TRACKLENGTH) { + // ------------------------------------------------------------------------- + // Non-tracklength estimators + // Find a match based on the exact time-position of the particle + + // Proceed only if inside time grid + if (t_end < time_grid_.back()) + return; + + auto time_bin = lower_bound_index(time_grid_.begin(), time_grid_.end(), t_end); + auto mesh_bin = model::meshes[mesh_]->get_bin(r); + + // Inside the mesh? + if (mesh_bin < 0) + return; + + auto bin = time_bin * mesh_n_bins_ + mesh_bin; + match.bins_.push_back(bin); + match.weights_.push_back(1.0); + + } else { + // ------------------------------------------------------------------------- + // For tracklength estimator, we have to check the start/end time-position + // of the current track and find where it overlaps with time-mesh bins and + // score accordingly. + + // Determine first time bin containing a portion of time interval + auto i_time_bin = lower_bound_index(time_grid_.begin(), time_grid_.end(), t_start); + + //std::cout<bins_crossed(last_r, r, u, match.bins_, match.weights_); + // Offset the bin location accordingly + match.bins_.back() += i_time_bin * mesh_n_bins_; + + /* + std::cout<<"end = start\n"; + double sum {0.0}; + for (int i = 0; i < match.bins_.size(); i++) { + std::cout<<" "<bins_crossed( + r_start, r_end, u, match.bins_, match.weights_); + const auto n_match = match.bins_.size() - n_match_old; + + // Update the newly-added bins and weights + const auto offset = i_time_bin * mesh_n_bins_; + for (int i = 0; i < n_match; i++) { + match.bins_[match.bins_.size() - i - 1] += offset; + match.weights_[match.weights_.size() - i - 1] *= fraction; + } + + if (t_end < time_grid_[i_time_bin + 1]) + break; + } + + /* + double sum {0.0}; + for (int i = 0; i < match.bins_.size(); i++) { + std::cout<<" "<id_); + if (translated_) { + write_dataset(filter_group, "translation", translation_); + } +} + +std::string TimedMeshFilter::text_label(int bin) const +{ + int bin_time = bin / mesh_n_bins_; + int bin_mesh = bin % mesh_n_bins_; + + auto& mesh = *model::meshes.at(mesh_); + + std::string label_time = fmt::format("Time [{}, {}) : ", time_grid_[bin], time_grid_[bin + 1]); + std::string label_mesh = mesh.bin_label(bin_mesh); + + return label_time + label_mesh; +} + +void TimedMeshFilter::set_translation(const Position& translation) +{ + translated_ = true; + translation_ = translation; +} + +void TimedMeshFilter::set_translation(const double translation[3]) +{ + this->set_translation({translation[0], translation[1], translation[2]}); +} + +} // namespace openmc diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 674987b8f13..710df0f2049 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -31,6 +31,7 @@ #include "openmc/tallies/filter_particle.h" #include "openmc/tallies/filter_sph_harm.h" #include "openmc/tallies/filter_surface.h" +#include "openmc/tallies/filter_timed_mesh.h" #include "openmc/xml_interface.h" #include "xtensor/xadapt.hpp" @@ -328,6 +329,14 @@ Tally::Tally(pugi::xml_node node) "an unstructured LibMesh tally."); } } + auto df2 = dynamic_cast(model::tally_filters[i].get()); + if (df2) { + auto lm = dynamic_cast(model::meshes[df2->mesh()].get()); + if (lm && estimator_ == TallyEstimator::TRACKLENGTH) { + fatal_error("A tracklength estimator cannot be used with " + "an unstructured LibMesh tally."); + } + } } #endif } From 46beb5491a003e8355bf3cb0cff88660baeb5dba Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 2 Aug 2024 11:29:04 +0700 Subject: [PATCH 04/17] add timed_mesh_filter.cpp (previously didn't get added...) --- include/openmc/tallies/filter_timed_mesh.h | 10 ++++---- src/tallies/filter_timed_mesh.cpp | 30 +++++++++++++--------- 2 files changed, 23 insertions(+), 17 deletions(-) diff --git a/include/openmc/tallies/filter_timed_mesh.h b/include/openmc/tallies/filter_timed_mesh.h index a73ed9a759f..bd6483cda9f 100644 --- a/include/openmc/tallies/filter_timed_mesh.h +++ b/include/openmc/tallies/filter_timed_mesh.h @@ -10,7 +10,7 @@ namespace openmc { //============================================================================== -//! Indexes the time-location of particle events to a time gridn and a regular +//! Indexes the time-location of particle events to a time gridn and a regular // mesh. For tracklength tallies, it will produce multiple valid bins and the // bin weight will correspond to the fraction of the track length that lies in //! that bin. @@ -52,18 +52,18 @@ class TimedMeshFilter : public Filter { const Position& translation() const { return translation_; } bool translated() const { return translated_; } - + const vector& time_grid() const { return time_grid_; } - + void set_time_grid(gsl::span time_grid); - + void reset_bins(); protected: //---------------------------------------------------------------------------- // Data members - int32_t mesh_; //!< Index of the mesh + int32_t mesh_; //!< Index of the mesh int mesh_n_bins_; bool translated_ {false}; //!< Whether or not the filter is translated Position translation_ {0.0, 0.0, 0.0}; //!< Filter translation diff --git a/src/tallies/filter_timed_mesh.cpp b/src/tallies/filter_timed_mesh.cpp index 38a6462c805..0a0e84b6037 100644 --- a/src/tallies/filter_timed_mesh.cpp +++ b/src/tallies/filter_timed_mesh.cpp @@ -37,7 +37,7 @@ void TimedMeshFilter::from_xml(pugi::xml_node node) if (check_for_node(node, "translation")) { set_translation(get_node_array(node, "translation")); } - + //---------------------------------------------------------------------------- // Time bins @@ -61,13 +61,14 @@ void TimedMeshFilter::set_time_grid(gsl::span time_grid) time_grid_.reserve(time_grid.size()); // Ensure time grid is sorted and don't have duplicates - if (std::adjacent_find(time_grid.cbegin(), time_grid.cend(), std::greater_equal<>()) != - time_grid.end()) { + if (std::adjacent_find(time_grid.cbegin(), time_grid.cend(), + std::greater_equal<>()) != time_grid.end()) { throw std::runtime_error {"Time grid must be monotonically increasing."}; } // Copy grid - std::copy(time_grid.cbegin(), time_grid.cend(), std::back_inserter(time_grid_)); + std::copy( + time_grid.cbegin(), time_grid.cend(), std::back_inserter(time_grid_)); reset_bins(); } @@ -94,7 +95,7 @@ void TimedMeshFilter::get_all_bins( Position r = p.r(); Position u = p.u(); const auto speed = p.speed(); - + // apply translation if present if (translated_) { last_r -= translation(); @@ -110,7 +111,8 @@ void TimedMeshFilter::get_all_bins( if (t_end < time_grid_.back()) return; - auto time_bin = lower_bound_index(time_grid_.begin(), time_grid_.end(), t_end); + auto time_bin = + lower_bound_index(time_grid_.begin(), time_grid_.end(), t_end); auto mesh_bin = model::meshes[mesh_]->get_bin(r); // Inside the mesh? @@ -123,18 +125,21 @@ void TimedMeshFilter::get_all_bins( } else { // ------------------------------------------------------------------------- - // For tracklength estimator, we have to check the start/end time-position + // For tracklength estimator, we have to check the start/end time-position // of the current track and find where it overlaps with time-mesh bins and // score accordingly. // Determine first time bin containing a portion of time interval - auto i_time_bin = lower_bound_index(time_grid_.begin(), time_grid_.end(), t_start); + auto i_time_bin = + lower_bound_index(time_grid_.begin(), time_grid_.end(), t_start); - //std::cout<bins_crossed(last_r, r, u, match.bins_, match.weights_); + model::meshes[mesh_]->bins_crossed( + last_r, r, u, match.bins_, match.weights_); // Offset the bin location accordingly match.bins_.back() += i_time_bin * mesh_n_bins_; @@ -181,7 +186,7 @@ void TimedMeshFilter::get_all_bins( if (t_end < time_grid_[i_time_bin + 1]) break; } - + /* double sum {0.0}; for (int i = 0; i < match.bins_.size(); i++) { @@ -210,7 +215,8 @@ std::string TimedMeshFilter::text_label(int bin) const auto& mesh = *model::meshes.at(mesh_); - std::string label_time = fmt::format("Time [{}, {}) : ", time_grid_[bin], time_grid_[bin + 1]); + std::string label_time = + fmt::format("Time [{}, {}) : ", time_grid_[bin], time_grid_[bin + 1]); std::string label_mesh = mesh.bin_label(bin_mesh); return label_time + label_mesh; From 19d4594f147b0d818ba38477916df375b233d817 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 2 Aug 2024 11:29:41 +0700 Subject: [PATCH 05/17] fix MG timing for fission secondaries --- src/physics_mg.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 361cf5affcf..37690ab1787 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -127,6 +127,7 @@ void create_fission_sites(Particle& p) 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()++; From ddde280e8147b7f469b7dd90e6094138f34df2c1 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 2 Aug 2024 11:34:22 +0700 Subject: [PATCH 06/17] fix minor bug --- src/physics_mg.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 37690ab1787..8fc877a427c 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -127,7 +127,7 @@ void create_fission_sites(Particle& p) SourceSite site; site.r = p.r(); site.particle = ParticleType::neutron; - site.time = p.time; + site.time = p.time(); site.wgt = 1. / weight; site.parent_id = p.id(); site.progeny_id = p.n_progeny()++; From cc7c75b5373ad45f0ab58c4218a6cccffdd4fc3b Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 2 Aug 2024 11:43:35 +0700 Subject: [PATCH 07/17] add MG speed getter --- src/particle.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/particle.cpp b/src/particle.cpp index 7d30e26bdf8..f81143e7b88 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -44,6 +44,15 @@ namespace openmc { double Particle::speed() const { + // Multigroup speed? + if (!settings::run_CE) { + auto& macro_xs = data::mg.macro_xs_[this->material()]; + int macro_t = this->mg_xs_cache().t; + int macro_a = macro_xs.get_angle_index(this->u()); + return 1.0 / macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, this->g(), nullptr, + nullptr, nullptr, macro_t, macro_a); + } + // Determine mass in eV/c^2 double mass; switch (this->type()) { From 9280e10ac2731d39bd25e89f50dc9ad2aa7cda3f Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 2 Aug 2024 11:48:42 +0700 Subject: [PATCH 08/17] fix bug in tracklength tallying at time-cutoff crossing --- src/particle.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/particle.cpp b/src/particle.cpp index f81143e7b88..d7fda63381a 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -253,6 +253,9 @@ void Particle::event_advance() double push_back_distance = speed() * dt; this->move_distance(-push_back_distance); hit_time_boundary = true; + + // Reduce the distance traveled for tallying + distance -= push_back_distance; } // Score track-length tallies From 7cb54697da5d037060f2de51d954c1a96c8e515b Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Sun, 4 Aug 2024 06:18:00 +0700 Subject: [PATCH 09/17] remove time filters from src/lib/filter.py --- openmc/lib/filter.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index 37d759fbb66..340c2fa3448 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -647,8 +647,6 @@ class ZernikeRadialFilter(ZernikeFilter): 'sphericalharmonics': SphericalHarmonicsFilter, 'spatiallegendre': SpatialLegendreFilter, 'surface': SurfaceFilter, - 'time': TimeFilter, - 'timedmesh': TimedMeshFilter, 'universe': UniverseFilter, 'zernike': ZernikeFilter, 'zernikeradial': ZernikeRadialFilter From a8f7521c9df6a7daba8bc008ad0dd7f06fcdf68f Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Sun, 22 Sep 2024 17:33:41 +0700 Subject: [PATCH 10/17] also generate mesh xml element for timed mesh filter --- openmc/tallies.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index f39ffa0789f..6a6f3d3680c 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3174,7 +3174,7 @@ def _create_mesh_subelements(self, root_element, memo=None): already_written = memo if memo else set() for tally in self: for f in tally.filters: - if isinstance(f, openmc.MeshFilter): + if isinstance(f, openmc.MeshFilter, openmc.TimedMeshFilter): if f.mesh.id in already_written: continue if len(f.mesh.name) > 0: From 4c86e6ef43591928d90be3734a4c0541743a5ba2 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Sat, 5 Oct 2024 20:25:40 +0700 Subject: [PATCH 11/17] fix a mistake in instance checking --- openmc/tallies.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 6a6f3d3680c..b69499969a3 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3174,7 +3174,8 @@ def _create_mesh_subelements(self, root_element, memo=None): already_written = memo if memo else set() for tally in self: for f in tally.filters: - if isinstance(f, openmc.MeshFilter, openmc.TimedMeshFilter): + if isinstance(f, openmc.MeshFilter) or\ + isinstance(f, openmc.TimedMeshFilter): if f.mesh.id in already_written: continue if len(f.mesh.name) > 0: From 3b30ffad5081f186ae99442c7dbf987f6fbc7b3e Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Sun, 10 Nov 2024 22:45:56 -0800 Subject: [PATCH 12/17] update timed mesh filter --- src/tallies/filter_timed_mesh.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tallies/filter_timed_mesh.cpp b/src/tallies/filter_timed_mesh.cpp index 0a0e84b6037..16ad96e3bd3 100644 --- a/src/tallies/filter_timed_mesh.cpp +++ b/src/tallies/filter_timed_mesh.cpp @@ -49,7 +49,7 @@ void TimedMeshFilter::set_mesh(int32_t mesh) { // perform any additional perparation for mesh tallies here mesh_ = mesh; - model::meshes[mesh_]->prepare_for_tallies(); + model::meshes[mesh_]->prepare_for_point_location(); reset_bins(); } From 36321b9a26123fde6fbf48468b14fc858c7ca0b9 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Mon, 2 Jun 2025 14:05:32 -0700 Subject: [PATCH 13/17] add a unit test --- docs/source/pythonapi/base.rst | 1 + openmc/filter.py | 2 +- openmc/tallies.py | 3 +-- tests/unit_tests/test_filters.py | 24 ++++++++++++++++++++++++ 4 files changed, 27 insertions(+), 3 deletions(-) diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index 1d8a93d47f1..e40965cf9e7 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -143,6 +143,7 @@ Constructing Tallies openmc.SpatialLegendreFilter openmc.SphericalHarmonicsFilter openmc.TimeFilter + openmc.TimedMeshFilter openmc.WeightFilter openmc.ZernikeFilter openmc.ZernikeRadialFilter diff --git a/openmc/filter.py b/openmc/filter.py index 5a70fa8b567..6e48e3ea344 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1268,7 +1268,7 @@ def _reset_bins(self): @property def shape(self): - return (self.time_grid - 1,) + self.mesh.dimension + return (len(self.time_grid) - 1,) + self.mesh.dimension @property def translation(self): diff --git a/openmc/tallies.py b/openmc/tallies.py index 6fb4f09321f..9d9cdc74566 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3285,8 +3285,7 @@ def _create_mesh_subelements(self, root_element, memo=None): already_written = memo if memo else set() for tally in self: for f in tally.filters: - if isinstance(f, openmc.MeshFilter) or\ - isinstance(f, openmc.TimedMeshFilter): + if type(f) in {openmc.MeshFilter, openmc.TimedMeshFilter}: if f.mesh.id in already_written: continue if len(f.mesh.name) > 0: diff --git a/tests/unit_tests/test_filters.py b/tests/unit_tests/test_filters.py index 60a60e3f56e..1b3cd291ae4 100644 --- a/tests/unit_tests/test_filters.py +++ b/tests/unit_tests/test_filters.py @@ -333,3 +333,27 @@ def test_weight(): new_f = openmc.Filter.from_xml_element(elem) assert new_f.id == f.id assert np.allclose(new_f.bins, f.bins) + +def test_timed_mesh(): + mesh = openmc.RegularMesh() + mesh.lower_left = (-1., -1., -1.) + mesh.upper_right = (1., 1., 1.) + mesh.dimension = (2, 4, 1) + + time_grid = np.linspace(0.0, 20.0, 41) + + f = openmc.TimedMeshFilter(mesh, time_grid) + + # Attributes + assert f.mesh == mesh + assert np.allclose(f.time_grid, time_grid) + assert f.shape == (40, 2, 4, 1) + + # to_xml_element() + elem = f.to_xml_element() + assert elem.tag == 'filter' + assert elem.attrib['type'] == 'timedmesh' + + # Test hash and str + hash(f) + str(f) From 5ec152fa9655e8203746cd4b0e4452c4b7aa2e2f Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Mon, 2 Jun 2025 15:39:18 -0700 Subject: [PATCH 14/17] fix gsl::span --- include/openmc/tallies/filter_timed_mesh.h | 3 ++- src/tallies/filter_timed_mesh.cpp | 3 +-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/include/openmc/tallies/filter_timed_mesh.h b/include/openmc/tallies/filter_timed_mesh.h index bd6483cda9f..82dec157814 100644 --- a/include/openmc/tallies/filter_timed_mesh.h +++ b/include/openmc/tallies/filter_timed_mesh.h @@ -5,6 +5,7 @@ #include "openmc/constants.h" #include "openmc/position.h" +#include "openmc/span.h" #include "openmc/tallies/filter.h" namespace openmc { @@ -55,7 +56,7 @@ class TimedMeshFilter : public Filter { const vector& time_grid() const { return time_grid_; } - void set_time_grid(gsl::span time_grid); + void set_time_grid(span time_grid); void reset_bins(); diff --git a/src/tallies/filter_timed_mesh.cpp b/src/tallies/filter_timed_mesh.cpp index 16ad96e3bd3..c506d639142 100644 --- a/src/tallies/filter_timed_mesh.cpp +++ b/src/tallies/filter_timed_mesh.cpp @@ -1,7 +1,6 @@ #include "openmc/tallies/filter_timed_mesh.h" #include -#include #include "openmc/capi.h" #include "openmc/constants.h" @@ -54,7 +53,7 @@ void TimedMeshFilter::set_mesh(int32_t mesh) reset_bins(); } -void TimedMeshFilter::set_time_grid(gsl::span time_grid) +void TimedMeshFilter::set_time_grid(span time_grid) { // Clear existing bins time_grid_.clear(); From 0d7ba71860b242ceb1e2391006d809dc62b48a77 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Mon, 2 Jun 2025 15:39:31 -0700 Subject: [PATCH 15/17] add a unit test --- tests/unit_tests/test_filter_timedmesh.py | 123 ++++++++++++++++++++++ 1 file changed, 123 insertions(+) create mode 100644 tests/unit_tests/test_filter_timedmesh.py diff --git a/tests/unit_tests/test_filter_timedmesh.py b/tests/unit_tests/test_filter_timedmesh.py new file mode 100644 index 00000000000..6a2249d75ef --- /dev/null +++ b/tests/unit_tests/test_filter_timedmesh.py @@ -0,0 +1,123 @@ +import math + +import numpy as np +import pytest +from uncertainties import unumpy + +import openmc + + +def test_filter_timed_Mesh(run_in_tmpdir): + """ + Test that TimeFilter+MeshFilter with collision estimator agree + with TimedMeshFilter with tracklength estimator. + """ + + # =========================================================================== + # Set Material + # =========================================================================== + + mat = openmc.Material() + mat.add_nuclide('Fe56', 1.0) + mat.set_density('g/cm3', 7.8) + + # =========================================================================== + # Set geometry + # =========================================================================== + + # Instantiate ZCylinder surfaces + surf_Z1 = openmc.XPlane(surface_id=1, x0=-1e10, boundary_type="reflective") + surf_Z2 = openmc.XPlane(surface_id=2, x0=1e10, boundary_type="reflective") + + # Instantiate Cells + cell_F = openmc.Cell(cell_id=1, name="F") + + # Use surface half-spaces to define regions + cell_F.region = +surf_Z1 & -surf_Z2 + + # Register Materials with Cells + cell_F.fill = mat + + # Instantiate Universes + root = openmc.Universe(universe_id=0, name="root universe", cells=[cell_F]) + + # Instantiate a Geometry, register the root Universe, and export to XML + geometry = openmc.Geometry(root) + + # =========================================================================== + # Settings + # =========================================================================== + + # Instantiate a Settings object, set all runtime parameters, and export to XML + settings_file = openmc.Settings() + settings_file.run_mode = "fixed source" + settings_file.particles = 100000 + settings_file.batches = 10 + settings_file.output = {"tallies": False} + settings_file.cutoff = {"time_neutron": 1E-7} + + # Create an initial uniform spatial source distribution over fissionable zones + delta_dist = openmc.stats.Point() + isotropic = openmc.stats.Isotropic() + settings_file.source = openmc.IndependentSource(space=delta_dist, angle=isotropic) + + # =========================================================================== + # Set tallies + # =========================================================================== + + # Create a mesh filter that can be used in a tally + mesh = openmc.RegularMesh() + mesh.dimension = (21, 1, 1) + mesh.lower_left = (-20.5, -1e10, -1e10) + mesh.upper_right = (20.5, 1e10, 1e10) + time_grid = np.linspace(0.0, 1E-7, 21) + + mesh_filter = openmc.MeshFilter(mesh) + time_filter = openmc.TimeFilter(time_grid) + timed_mesh_filter = openmc.TimedMeshFilter(mesh, time_grid) + + # Now use the mesh filter in a tally and indicate what scores are desired + tally1 = openmc.Tally(name="collision") + tally1.estimator = "collision" + tally1.filters = [time_filter, mesh_filter] + tally1.scores = ["flux"] + + tally2 = openmc.Tally(name="tracklength") + tally2.estimator = "tracklength" + tally2.filters = [timed_mesh_filter] + tally2.scores = ["flux"] + + # Instantiate a tallies collection + tallies = openmc.Tallies([tally1, tally2]) + + # =========================================================================== + # Set the model + # =========================================================================== + + model = openmc.Model() + model.geometry = geometry + model.settings = settings_file + model.tallies = tallies + + # =========================================================================== + # Run and post-process + # =========================================================================== + + sp_filename = model.run() + + # Get radial flux distribution + with openmc.StatePoint(sp_filename) as sp: + flux_collision = sp.tallies[tally1.id].mean.ravel() + flux_collision_unc = sp.tallies[tally1.id].std_dev.ravel() + flux_tracklength = sp.tallies[tally2.id].mean.ravel() + flux_tracklength_unc = sp.tallies[tally2.id].std_dev.ravel() + + # Construct arrays with uncertainties + collision = unumpy.uarray(flux_collision, flux_collision_unc) + tracklength = unumpy.uarray(flux_tracklength, flux_tracklength_unc) + delta = collision - tracklength + + # Check that difference is within uncertainty + diff = unumpy.nominal_values(delta) + std_dev = unumpy.std_devs(delta) + assert np.all(diff < 3*std_dev) From fcd84604d7f2d2c35656d02d76d0671deaad64a0 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 12 Jun 2025 09:41:57 +0700 Subject: [PATCH 16/17] remove debugging prints --- src/tallies/filter_timed_mesh.cpp | 43 +++++++------------------------ 1 file changed, 10 insertions(+), 33 deletions(-) diff --git a/src/tallies/filter_timed_mesh.cpp b/src/tallies/filter_timed_mesh.cpp index c506d639142..36d1410b3e2 100644 --- a/src/tallies/filter_timed_mesh.cpp +++ b/src/tallies/filter_timed_mesh.cpp @@ -82,8 +82,8 @@ void TimedMeshFilter::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(); + const double t_start = p.time_last(); + const double t_end = p.time(); // If time interval is entirely out of time bin range, exit if (t_end < time_grid_.front() || t_start >= time_grid_.back()) @@ -93,7 +93,7 @@ void TimedMeshFilter::get_all_bins( Position last_r = p.r_last(); Position r = p.r(); Position u = p.u(); - const auto speed = p.speed(); + const double speed = p.speed(); // apply translation if present if (translated_) { @@ -110,15 +110,15 @@ void TimedMeshFilter::get_all_bins( if (t_end < time_grid_.back()) return; - auto time_bin = + int time_bin = lower_bound_index(time_grid_.begin(), time_grid_.end(), t_end); - auto mesh_bin = model::meshes[mesh_]->get_bin(r); + int mesh_bin = model::meshes[mesh_]->get_bin(r); // Inside the mesh? if (mesh_bin < 0) return; - auto bin = time_bin * mesh_n_bins_ + mesh_bin; + int bin = time_bin * mesh_n_bins_ + mesh_bin; match.bins_.push_back(bin); match.weights_.push_back(1.0); @@ -129,29 +129,15 @@ void TimedMeshFilter::get_all_bins( // score accordingly. // Determine first time bin containing a portion of time interval - auto i_time_bin = + int i_time_bin = lower_bound_index(time_grid_.begin(), time_grid_.end(), t_start); - // std::cout<bins_crossed( last_r, r, u, match.bins_, match.weights_); // Offset the bin location accordingly match.bins_.back() += i_time_bin * mesh_n_bins_; - - /* - std::cout<<"end = start\n"; - double sum {0.0}; - for (int i = 0; i < match.bins_.size(); i++) { - std::cout<<" "<bins_crossed( r_start, r_end, u, match.bins_, match.weights_); - const auto n_match = match.bins_.size() - n_match_old; + const int n_match = match.bins_.size() - n_match_old; // Update the newly-added bins and weights - const auto offset = i_time_bin * mesh_n_bins_; + const int offset = i_time_bin * mesh_n_bins_; for (int i = 0; i < n_match; i++) { match.bins_[match.bins_.size() - i - 1] += offset; match.weights_[match.weights_.size() - i - 1] *= fraction; @@ -185,15 +171,6 @@ void TimedMeshFilter::get_all_bins( if (t_end < time_grid_[i_time_bin + 1]) break; } - - /* - double sum {0.0}; - for (int i = 0; i < match.bins_.size(); i++) { - std::cout<<" "< Date: Thu, 12 Jun 2025 13:41:40 +0700 Subject: [PATCH 17/17] skip timed-mesh filtering if zero time interval --- src/tallies/filter_timed_mesh.cpp | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/tallies/filter_timed_mesh.cpp b/src/tallies/filter_timed_mesh.cpp index 36d1410b3e2..980b604fd08 100644 --- a/src/tallies/filter_timed_mesh.cpp +++ b/src/tallies/filter_timed_mesh.cpp @@ -132,12 +132,8 @@ void TimedMeshFilter::get_all_bins( int i_time_bin = lower_bound_index(time_grid_.begin(), time_grid_.end(), t_start); - // If time interval is zero, add a match corresponding to the starting time + // If time interval is zero, tracklength is zero, no score if (t_end == t_start) { - model::meshes[mesh_]->bins_crossed( - last_r, r, u, match.bins_, match.weights_); - // Offset the bin location accordingly - match.bins_.back() += i_time_bin * mesh_n_bins_; return; }