diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 64b62f3e5b..2ca9fc430f 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -74,6 +74,8 @@ extern "C" { int openmc_material_filter_set_bins(int32_t index, size_t n, const int32_t* bins); int openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh); int openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh); + int openmc_mesh_filter_get_translation(int32_t index, double translation[3]); + int openmc_mesh_filter_set_translation(int32_t index, double translation[3]); int openmc_mesh_get_id(int32_t index, int32_t* id); int openmc_mesh_set_id(int32_t index, int32_t id); int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh); diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index e9519ca0a1..2af2e87df2 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -71,18 +71,28 @@ public: //! Determine which bins were crossed by a particle // - //! \param[in] p Particle to check + //! \param[in] r0 Previous position of the particle + //! \param[in] r1 Current position of the particle + //! \param[in] u Particle direction //! \param[out] bins Bins that were crossed //! \param[out] lengths Fraction of tracklength in each bin - virtual void bins_crossed(const Particle& p, std::vector& bins, + virtual void bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins, std::vector& lengths) const = 0; //! Determine which surface bins were crossed by a particle // - //! \param[in] p Particle to check + //! \param[in] r0 Previous position of the particle + //! \param[in] r1 Current position of the particle + //! \param[in] u Particle direction //! \param[out] bins Surface bins that were crossed virtual void - surface_bins_crossed(const Particle& p, std::vector& bins) const = 0; + surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const = 0; //! Get bin at a given position in space // @@ -137,7 +147,10 @@ public: int n_surface_bins() const override; - void bins_crossed(const Particle& p, std::vector& bins, + void bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins, std::vector& lengths) const override; //! Count number of bank sites in each mesh bin / energy bin @@ -219,9 +232,11 @@ public: RegularMesh(pugi::xml_node node); // Overriden methods - - void surface_bins_crossed(const Particle& p, std::vector& bins) - const override; + void + surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const override; int get_index_in_direction(double r, int i) const override; @@ -259,9 +274,11 @@ public: RectilinearMesh(pugi::xml_node node); // Overriden methods - - void surface_bins_crossed(const Particle& p, std::vector& bins) - const override; + void + surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const override; int get_index_in_direction(double r, int i) const override; @@ -324,8 +341,9 @@ public: std::string bin_label(int bin) const override; - void surface_bins_crossed(const Particle& p, - std::vector& bins) const override; + void surface_bins_crossed(Position r0, + Position r1, + std::vector& bins) const; void to_hdf5(hid_t group) const override; @@ -348,10 +366,18 @@ public: MOABMesh(pugi::xml_node); MOABMesh(const std::string& filename); - void bins_crossed(const Particle& p, + void bins_crossed(Position r0, + Position r1, + const Direction& u, std::vector& bins, std::vector& lengths) const override; + void + surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const override; + int get_bin(Position r) const; int n_bins() const override; @@ -496,10 +522,18 @@ public: LibMesh(const std::string& filename); // Methods - void bins_crossed(const Particle& p, + void bins_crossed(Position r0, + Position r1, + const Direction& u, std::vector& bins, std::vector& lengths) const override; + void + surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const override; + int get_bin(Position r) const override; int n_bins() const override; diff --git a/include/openmc/tallies/filter_mesh.h b/include/openmc/tallies/filter_mesh.h index fa64877e36..8a712f4cb0 100644 --- a/include/openmc/tallies/filter_mesh.h +++ b/include/openmc/tallies/filter_mesh.h @@ -4,6 +4,7 @@ #include #include "openmc/tallies/filter.h" +#include "openmc/position.h" namespace openmc { @@ -42,11 +43,22 @@ public: virtual void set_mesh(int32_t mesh); + virtual void set_translation(const Position& translation); + + virtual void set_translation(const double translation[3]); + + virtual const Position& translation() const {return translation_;} + + virtual bool translated() const {return translated_;} + + protected: //---------------------------------------------------------------------------- // Data members int32_t mesh_; + bool translated_ {false}; + Position translation_ {0.0, 0.0, 0.0}; }; } // namespace openmc diff --git a/openmc/filter.py b/openmc/filter.py index 2d8faf3784..8da3cf78f6 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -716,6 +716,9 @@ class MeshFilter(Filter): The mesh object that events will be tallied onto 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), ...] @@ -727,6 +730,7 @@ class MeshFilter(Filter): def __init__(self, mesh, filter_id=None): self.mesh = mesh self.id = filter_id + self._translation = None def __hash__(self): string = type(self).__name__ + '\n' @@ -737,6 +741,7 @@ class MeshFilter(Filter): string = type(self).__name__ + '\n' string += '{: <16}=\t{}\n'.format('\tMesh ID', self.mesh.id) string += '{: <16}=\t{}\n'.format('\tID', self.id) + string += '{: <16}=\t{}\n'.format('\tTranslation', self.translation) return string @classmethod @@ -754,8 +759,13 @@ class MeshFilter(Filter): mesh_obj = kwargs['meshes'][mesh_id] filter_id = int(group.name.split('/')[-1].lstrip('filter ')) + out = cls(mesh_obj, filter_id=filter_id) + translation = group.get('translation') + if translation: + out.translation = translation[()] + return out @property @@ -774,6 +784,16 @@ class MeshFilter(Filter): else: self.bins = list(mesh.indices) + @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): # Mesh filters cannot have more than one bin return False @@ -854,6 +874,8 @@ class MeshFilter(Filter): """ element = super().to_xml_element() element[0].text = str(self.mesh.id) + if self.translation is not None: + element.set('translation', ' '.join(map(str, self.translation))) return element @@ -873,6 +895,9 @@ class MeshSurfaceFilter(MeshFilter): The mesh ID mesh : openmc.MeshBase The mesh object that events will be tallied onto + translation : Iterable of float + This array specifies a vector that is used to translate (shift) + the mesh for this filter id : int Unique identifier for the filter bins : list of tuple @@ -971,9 +996,9 @@ class CollisionFilter(Filter): Parameters ---------- bins : Iterable of int - A list or iterable of the number of collisions, as integer values. - The events whose post-scattering collision number equals one of - the provided values will be counted. + A list or iterable of the number of collisions, as integer values. + The events whose post-scattering collision number equals one of + the provided values will be counted. filter_id : int Unique identifier for the filter diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index 497d4592b9..5d2231cc2e 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -84,6 +84,12 @@ _dll.openmc_meshsurface_filter_get_mesh.errcheck = _error_handler _dll.openmc_meshsurface_filter_set_mesh.argtypes = [c_int32, c_int32] _dll.openmc_meshsurface_filter_set_mesh.restype = c_int _dll.openmc_meshsurface_filter_set_mesh.errcheck = _error_handler +_dll.openmc_mesh_filter_get_translation.argtypes = [c_int32, POINTER(c_double*3)] +_dll.openmc_mesh_filter_get_translation.restype = c_int +_dll.openmc_mesh_filter_get_translation.errcheck = _error_handler +_dll.openmc_mesh_filter_set_translation.argtypes = [c_int32, POINTER(c_double*3)] +_dll.openmc_mesh_filter_set_translation.restype = c_int +_dll.openmc_mesh_filter_set_translation.errcheck = _error_handler _dll.openmc_new_filter.argtypes = [c_char_p, POINTER(c_int32)] _dll.openmc_new_filter.restype = c_int _dll.openmc_new_filter.errcheck = _error_handler @@ -325,6 +331,16 @@ class MeshFilter(Filter): def mesh(self, mesh): _dll.openmc_mesh_filter_set_mesh(self._index, mesh._index) + @property + def translation(self): + translation = (c_double*3)() + _dll.openmc_mesh_filter_get_translation(self._index, translation) + return tuple(translation) + + @translation.setter + def translation(self, translation): + _dll.openmc_mesh_filter_set_translation(self._index, (c_double*3)(*translation)) + class MeshSurfaceFilter(Filter): filter_type = 'meshsurface' @@ -344,6 +360,16 @@ class MeshSurfaceFilter(Filter): def mesh(self, mesh): _dll.openmc_meshsurface_filter_set_mesh(self._index, mesh._index) + @property + def translation(self): + translation = (c_double*3)() + _dll.openmc_mesh_filter_get_translation(self._index, translation) + return tuple(translation) + + @translation.setter + def translation(self, translation): + _dll.openmc_mesh_filter_set_translation(self._index, (c_double*3)(*translation)) + class MuFilter(Filter): filter_type = 'mu' diff --git a/openmc/mesh.py b/openmc/mesh.py index 459b48a694..16e25cddc4 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -471,6 +471,8 @@ class RectilinearMesh(MeshBase): Unique identifier for the mesh name : str Name of the mesh + dimension : Iterable of int + The number of mesh cells in each direction. n_dimension : int Number of mesh dimensions (always 3 for a RectilinearMesh). x_grid : Iterable of float @@ -492,6 +494,12 @@ class RectilinearMesh(MeshBase): self._y_grid = None self._z_grid = None + @property + def dimension(self): + return (len(self.x_grid) - 1, + len(self.y_grid) - 1, + len(self.z_grid) - 1) + @property def n_dimension(self): return 3 diff --git a/src/mesh.cpp b/src/mesh.cpp index 16ce630ae0..9e974477cb 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -190,7 +190,8 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) { } void -UnstructuredMesh::surface_bins_crossed(const Particle& p, +UnstructuredMesh::surface_bins_crossed(Position r0, + Position r1, std::vector& bins) const { fatal_error("Unstructured mesh surface tallies are not implemented."); } @@ -581,26 +582,24 @@ bool StructuredMesh::intersects_3d(Position& r0, Position r1, int* ijk) const return min_dist < INFTY; } -void StructuredMesh::bins_crossed(const Particle& p, std::vector& bins, - std::vector& lengths) const +void StructuredMesh::bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins, + std::vector& lengths) const { // ======================================================================== // Determine where the track intersects the mesh and if it intersects at all. - // Copy the starting and ending coordinates of the particle. - Position last_r {p.r_last_}; - Position r {p.r()}; - Direction u {p.u()}; - // Compute the length of the entire track. - double total_distance = (r - last_r).norm(); + double total_distance = (r1 - r0).norm(); if (total_distance == 0.0) return; // While determining if this track intersects the mesh, offset the starting // and ending coords by a bit. This avoid finite-precision errors that can // occur when the mesh surfaces coincide with lattice or geometric surfaces. - Position r0 = last_r + TINY_BIT*u; - Position r1 = r - TINY_BIT*u; + Position last_r = r0 + TINY_BIT*u; + Position r = r1 - TINY_BIT*u; // Determine the mesh indices for the starting and ending coords. Here, we // use arrays for ijk0 and ijk1 instead of std::vector because we obtain a @@ -611,21 +610,21 @@ void StructuredMesh::bins_crossed(const Particle& p, std::vector& bins, int n = n_dimension_; int ijk0[3], ijk1[3]; bool start_in_mesh; - get_indices(r0, ijk0, &start_in_mesh); + get_indices(last_r, ijk0, &start_in_mesh); bool end_in_mesh; - get_indices(r1, ijk1, &end_in_mesh); + get_indices(r, ijk1, &end_in_mesh); // Reset coordinates and check for a mesh intersection if necessary. if (start_in_mesh) { // The initial coords lie in the mesh, use those coords for tallying. - r0 = last_r; + last_r = r0; } else { // The initial coords do not lie in the mesh. Check to see if the particle // eventually intersects the mesh and compute the relevant coords and // indices. - if (!intersects(r0, r1, ijk0)) return; + if (!intersects(last_r, r, ijk0)) return; } - r1 = r; + r = r1; // The TINY_BIT offsets above mean that the preceding logic cannot always find // the correct ijk0 and ijk1 indices. For tracks shorter than 2*TINY_BIT, just @@ -644,7 +643,7 @@ void StructuredMesh::bins_crossed(const Particle& p, std::vector& bins, if (std::equal(ijk0, ijk0 + n, ijk1)) { // The track ends in this cell. Use the particle end location rather // than the mesh surface and stop iterating. - double distance = (r1 - r0).norm(); + double distance = (last_r - r).norm(); bins.push_back(get_bin_from_indices(ijk0)); lengths.push_back(distance / total_distance); break; @@ -657,10 +656,10 @@ void StructuredMesh::bins_crossed(const Particle& p, std::vector& bins, d[k] = INFTY; } else if (u[k] > 0) { double xyz_cross = positive_grid_boundary(ijk0, k); - d[k] = (xyz_cross - r0[k]) / u[k]; + d[k] = (xyz_cross - last_r[k]) / u[k]; } else { double xyz_cross = negative_grid_boundary(ijk0, k); - d[k] = (xyz_cross - r0[k]) / u[k]; + d[k] = (xyz_cross - last_r[k]) / u[k]; } } @@ -671,7 +670,7 @@ void StructuredMesh::bins_crossed(const Particle& p, std::vector& bins, lengths.push_back(distance / total_distance); // Translate to the oncoming mesh surface. - r0 += distance * u; + last_r += distance * u; // Increment the indices into the next mesh cell. if (u[j] > 0.0) { @@ -797,17 +796,12 @@ double RegularMesh::negative_grid_boundary(int* ijk, int i) const return lower_left_[i] + (ijk[i] - 1) * width_[i]; } -void RegularMesh::surface_bins_crossed(const Particle& p, - std::vector& bins) const +void +RegularMesh::surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const { - // ======================================================================== - // Determine if the track intersects the tally mesh. - - // Copy the starting and ending coordinates of the particle. - Position r0 {p.r_last_current_}; - Position r1 {p.r()}; - Direction u {p.u()}; - // Determine indices for starting and ending location. int n = n_dimension_; std::vector ijk0(n), ijk1(n); @@ -1107,17 +1101,11 @@ int RectilinearMesh::set_grid() return 0; } -void RectilinearMesh::surface_bins_crossed(const Particle& p, +void RectilinearMesh::surface_bins_crossed(Position r0, + Position r1, + const Direction& u, std::vector& bins) const { - // ======================================================================== - // Determine if the track intersects the tally mesh. - - // Copy the starting and ending coordinates of the particle. - Position r0 {p.r_last_current_}; - Position r1 {p.r()}; - Direction u {p.u()}; - // Determine indices for starting and ending location. int ijk0[3], ijk1[3]; bool start_in_mesh; @@ -1444,7 +1432,7 @@ openmc_get_mesh_index(int32_t id, int32_t* index) return 0; } -// Return the ID of a mesh +//! Return the ID of a mesh extern "C" int openmc_mesh_get_id(int32_t index, int32_t* id) { @@ -1701,26 +1689,25 @@ MOABMesh::intersect_track(const moab::CartVect& start, } void -MOABMesh::bins_crossed(const Particle& p, +MOABMesh::bins_crossed(Position r0, + Position r1, + const Direction& u, std::vector& bins, std::vector& lengths) const { - Position last_r{p.r_last_}; - Position r{p.r()}; - Direction u{p.u()}; - u /= u.norm(); - moab::CartVect r0(last_r.x, last_r.y, last_r.z); - moab::CartVect r1(r.x, r.y, r.z); + moab::CartVect start(r0.x, r0.y, r0.z); + moab::CartVect end(r1.x, r1.y, r1.z); moab::CartVect dir(u.x, u.y, u.z); + dir.normalize(); - double track_len = (r1 - r0).length(); + double track_len = (end - start).length(); if (track_len == 0.0) return; - r0 -= TINY_BIT * dir; - r1 += TINY_BIT * dir; + start -= TINY_BIT * dir; + end += TINY_BIT * dir; std::vector hits; - intersect_track(r0, dir, track_len, hits); + intersect_track(start, dir, track_len, hits); bins.clear(); lengths.clear(); @@ -1729,7 +1716,7 @@ MOABMesh::bins_crossed(const Particle& p, // within a single tet. If this is the case, apply entire // score to that tet and return. if (hits.size() == 0) { - Position midpoint = last_r + u * (track_len * 0.5); + Position midpoint = r0 + u * (track_len * 0.5); int bin = this->get_bin(midpoint); if (bin != -1) { bins.push_back(bin); @@ -1740,7 +1727,7 @@ MOABMesh::bins_crossed(const Particle& p, // for each segment in the set of tracks, try to look up a tet // at the midpoint of the segment - Position current = last_r; + Position current = r0; double last_dist = 0.0; for (const auto& hit : hits) { // get the segment length @@ -1752,7 +1739,7 @@ MOABMesh::bins_crossed(const Particle& p, int bin = this->get_bin(midpoint); // determine the start point for this segment - current = last_r + u * hit; + current = r0 + u * hit; if (bin == -1) { continue; @@ -1767,7 +1754,7 @@ MOABMesh::bins_crossed(const Particle& p, // the last segment of the track is in the mesh but doesn't // reach the other side of the tet if (hits.back() < track_len) { - Position segment_start = last_r + u * hits.back(); + Position segment_start = r0 + u * hits.back(); double segment_length = track_len - hits.back(); Position midpoint = segment_start + u * (segment_length * 0.5); int bin = this->get_bin(midpoint); @@ -1833,6 +1820,16 @@ double MOABMesh::tet_volume(moab::EntityHandle tet) const return 1.0 / 6.0 * (((p[1] - p[0]) * (p[2] - p[0])) % (p[3] - p[0])); } +void MOABMesh::surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const +{ + + // TODO: Implement triangle crossings here + throw std::runtime_error{"Unstructured mesh surface tallies are not implemented."}; +} + int MOABMesh::get_bin(Position r) const { moab::EntityHandle tet = get_tet(r); @@ -2212,8 +2209,18 @@ int LibMesh::n_bins() const return m_->n_elem(); } -int LibMesh::n_surface_bins() const +void +LibMesh::surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const { + // TODO: Implement triangle crossings here + throw std::runtime_error{"Unstructured mesh surface tallies are not implemented."}; +} + +int +LibMesh::n_surface_bins() const { int n_bins = 0; for (int i = 0; i < this->n_bins(); i++) { const libMesh::Elem& e = get_element_from_bin(i); @@ -2300,7 +2307,9 @@ void LibMesh::write(const std::string& filename) const } void -LibMesh::bins_crossed(const Particle& p, +LibMesh::bins_crossed(Position r0, + Position r1, + const Direction& u, std::vector& bins, std::vector& lengths) const { diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp index 7897173868..1595de70c9 100644 --- a/src/tallies/filter_mesh.cpp +++ b/src/tallies/filter_mesh.cpp @@ -28,20 +28,36 @@ MeshFilter::from_xml(pugi::xml_node node) 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")); + } + } void MeshFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { + + Position last_r = p.r_last_; + Position r = p.r(); + Position u = p.u(); + + // apply translation if present + if (translated_) { + last_r -= translation(); + r -= translation(); + } + if (estimator != TallyEstimator::TRACKLENGTH) { - auto bin = model::meshes[mesh_]->get_bin(p.r()); + auto bin = model::meshes[mesh_]->get_bin(r); if (bin >= 0) { match.bins_.push_back(bin); match.weights_.push_back(1.0); } } else { - model::meshes[mesh_]->bins_crossed(p, match.bins_, match.weights_); + model::meshes[mesh_]->bins_crossed(last_r, r, u, match.bins_, match.weights_); } } @@ -50,13 +66,17 @@ MeshFilter::to_statepoint(hid_t filter_group) const { Filter::to_statepoint(filter_group); write_dataset(filter_group, "bins", model::meshes[mesh_]->id_); + if (translated_) { + write_dataset(filter_group, "translation", translation_); + } } std::string MeshFilter::text_label(int bin) const { auto& mesh = *model::meshes.at(mesh_); - return mesh.bin_label(bin); + std::string label = mesh.bin_label(bin); + return label; } void @@ -66,6 +86,17 @@ MeshFilter::set_mesh(int32_t mesh) n_bins_ = model::meshes[mesh_]->n_bins(); } +void MeshFilter::set_translation(const Position& translation) +{ + translated_ = true; + translation_ = translation; +} + +void MeshFilter::set_translation(const double translation[3]) +{ + this->set_translation({translation[0], translation[1], translation[2]}); +} + //============================================================================== // C-API functions //============================================================================== @@ -123,4 +154,47 @@ openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh) return 0; } +extern "C" int +openmc_mesh_filter_get_translation(int32_t index, double translation[3]) +{ + // Make sure this is a valid index to an allocated filter + if (int err = verify_filter(index)) return err; + + // Check the filter type + const auto& filter = model::tally_filters[index]; + if (filter->type() != "mesh" && filter->type() != "meshsurface") { + set_errmsg("Tried to get a translation from a non-mesh-based filter."); + return OPENMC_E_INVALID_TYPE; + } + + // Get translation from the mesh filter and set value + auto mesh_filter = dynamic_cast(filter.get()); + const auto& t = mesh_filter->translation(); + for (int i = 0; i < 3; i++) { translation[i] = t[i]; } + + return 0; +} + +extern "C" int +openmc_mesh_filter_set_translation(int32_t index, double translation[3]) +{ + // Make sure this is a valid index to an allocated filter + if (int err = verify_filter(index)) return err; + + const auto& filter = model::tally_filters[index]; + // Check the filter type + if (filter->type() != "mesh" && filter->type() != "meshsurface") { + set_errmsg("Tried to set mesh on a non-mesh-based filter."); + return OPENMC_E_INVALID_TYPE; + } + + // Get a pointer to the filter and downcast + auto mesh_filter = dynamic_cast(filter.get()); + + // Set the translation + mesh_filter->set_translation(translation); + + return 0; +} + } // namespace openmc diff --git a/src/tallies/filter_meshsurface.cpp b/src/tallies/filter_meshsurface.cpp index 997fef7293..af595ad52a 100644 --- a/src/tallies/filter_meshsurface.cpp +++ b/src/tallies/filter_meshsurface.cpp @@ -11,7 +11,15 @@ void MeshSurfaceFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - model::meshes[mesh_]->surface_bins_crossed(p, match.bins_); + Position r0 = p.r_last_current_; + Position r1 = p.r(); + if (translated_) { + r0 -= translation(); + r1 -= translation(); + } + + Direction u = p.u(); + model::meshes[mesh_]->surface_bins_crossed(r0, r1, u, match.bins_); for (auto b : match.bins_) match.weights_.push_back(1.0); } diff --git a/tests/regression_tests/filter_translations/__init__.py b/tests/regression_tests/filter_translations/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/filter_translations/inputs_true.dat b/tests/regression_tests/filter_translations/inputs_true.dat new file mode 100644 index 0000000000..d38bec39f1 --- /dev/null +++ b/tests/regression_tests/filter_translations/inputs_true.dat @@ -0,0 +1,84 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 5 + 0 + + + + + 3 4 5 + -9 -9 -9 + 9 9 9 + + + -9.0 0.0 9.0 + -9.0 -3.0 3.0 9.0 + -9.0 -4.5 0.0 4.5 9.0 + + + 3 4 5 + -19 -4 -9 + -1 14 9 + + + -19.0 -10.0 -1.0 + -4.0 2.0 8.0 14.0 + -9.0 -4.5 0.0 4.5 9.0 + + + 1 + + + 2 + + + 3 + + + 4 + + + 1 + total + + + 2 + total + + + 3 + total + + + 4 + total + + diff --git a/tests/regression_tests/filter_translations/results_true.dat b/tests/regression_tests/filter_translations/results_true.dat new file mode 100644 index 0000000000..e9813b183e --- /dev/null +++ b/tests/regression_tests/filter_translations/results_true.dat @@ -0,0 +1,342 @@ +k-combined: +7.656044E-01 5.168921E-02 +tally 1: +4.195134E-02 +3.639654E-04 +7.345058E-02 +1.082591E-03 +5.325577E-02 +5.772001E-04 +1.390109E-01 +3.954675E-03 +3.829062E-01 +3.001419E-02 +1.408217E-01 +4.113283E-03 +1.412139E-01 +4.218445E-03 +3.567122E-01 +2.663446E-02 +1.393013E-01 +3.978946E-03 +5.117767E-02 +5.394412E-04 +8.111129E-02 +1.362452E-03 +5.594852E-02 +6.333549E-04 +6.298061E-02 +8.073054E-04 +9.082908E-02 +1.668351E-03 +6.783427E-02 +9.250122E-04 +1.999222E-01 +8.028088E-03 +5.909662E-01 +7.061368E-02 +1.896895E-01 +7.242074E-03 +1.877543E-01 +7.145844E-03 +6.189100E-01 +7.667032E-02 +1.847604E-01 +6.905108E-03 +6.017074E-02 +7.311068E-04 +1.109867E-01 +2.474797E-03 +7.322493E-02 +1.112742E-03 +6.859899E-02 +9.576740E-04 +1.062208E-01 +2.360724E-03 +6.762553E-02 +9.247454E-04 +2.367276E-01 +1.130443E-02 +9.674870E-01 +2.606562E-01 +2.113618E-01 +8.988266E-03 +2.079097E-01 +8.765092E-03 +9.535566E-01 +2.216359E-01 +2.131094E-01 +9.143533E-03 +6.644946E-02 +9.024189E-04 +1.064224E-01 +2.284439E-03 +8.143154E-02 +1.362849E-03 +6.209959E-02 +7.875797E-04 +9.702493E-02 +1.902685E-03 +6.416842E-02 +8.369934E-04 +1.727520E-01 +6.027045E-03 +5.966710E-01 +7.197863E-02 +1.889340E-01 +7.213441E-03 +1.914452E-01 +7.420265E-03 +5.726219E-01 +6.723631E-02 +1.646287E-01 +5.489910E-03 +5.663491E-02 +6.494568E-04 +9.923269E-02 +2.053123E-03 +6.809958E-02 +9.386140E-04 +4.850339E-02 +4.856214E-04 +8.169408E-02 +1.354206E-03 +5.474566E-02 +6.128204E-04 +1.208954E-01 +3.074917E-03 +3.393101E-01 +2.398955E-02 +1.358915E-01 +3.758521E-03 +1.347810E-01 +3.718991E-03 +3.404279E-01 +2.380238E-02 +1.315829E-01 +3.585161E-03 +4.754269E-02 +4.593490E-04 +7.591439E-02 +1.170776E-03 +5.121227E-02 +5.346435E-04 +tally 2: +2.194240E-01 +9.811952E-03 +2.718491E-01 +1.501408E-02 +5.823038E-01 +6.967280E-02 +5.955012E-01 +7.246929E-02 +2.496521E-01 +1.268521E-02 +2.623766E-01 +1.397567E-02 +3.475667E-01 +2.453188E-02 +3.290723E-01 +2.187833E-02 +1.123843E+00 +2.718539E-01 +1.060797E+00 +2.384769E-01 +3.373095E-01 +2.280788E-02 +3.563621E-01 +2.594087E-02 +3.442434E-01 +2.400708E-02 +3.319540E-01 +2.206590E-02 +1.051427E+00 +2.459968E-01 +1.065112E+00 +2.465202E-01 +3.182627E-01 +2.049613E-02 +3.724575E-01 +2.824300E-02 +2.274581E-01 +1.052009E-02 +2.558629E-01 +1.318940E-02 +5.434276E-01 +6.081708E-02 +5.623197E-01 +6.408372E-02 +2.320052E-01 +1.096050E-02 +2.380175E-01 +1.168927E-02 +tally 3: +4.195134E-02 +3.639654E-04 +7.345058E-02 +1.082591E-03 +5.325577E-02 +5.772001E-04 +1.390109E-01 +3.954675E-03 +3.829062E-01 +3.001419E-02 +1.408217E-01 +4.113283E-03 +1.412139E-01 +4.218445E-03 +3.567122E-01 +2.663446E-02 +1.393013E-01 +3.978946E-03 +5.117767E-02 +5.394412E-04 +8.111129E-02 +1.362452E-03 +5.594852E-02 +6.333549E-04 +6.298061E-02 +8.073054E-04 +9.082908E-02 +1.668351E-03 +6.783427E-02 +9.250122E-04 +1.999222E-01 +8.028088E-03 +5.909662E-01 +7.061368E-02 +1.896895E-01 +7.242074E-03 +1.877543E-01 +7.145844E-03 +6.189100E-01 +7.667032E-02 +1.847604E-01 +6.905108E-03 +6.017074E-02 +7.311068E-04 +1.109867E-01 +2.474797E-03 +7.322493E-02 +1.112742E-03 +6.859899E-02 +9.576740E-04 +1.062208E-01 +2.360724E-03 +6.762553E-02 +9.247454E-04 +2.367276E-01 +1.130443E-02 +9.674870E-01 +2.606562E-01 +2.113618E-01 +8.988266E-03 +2.079097E-01 +8.765092E-03 +9.535566E-01 +2.216359E-01 +2.131094E-01 +9.143533E-03 +6.644946E-02 +9.024189E-04 +1.064224E-01 +2.284439E-03 +8.143154E-02 +1.362849E-03 +6.209959E-02 +7.875797E-04 +9.702493E-02 +1.902685E-03 +6.416842E-02 +8.369934E-04 +1.727520E-01 +6.027045E-03 +5.966710E-01 +7.197863E-02 +1.889340E-01 +7.213441E-03 +1.914452E-01 +7.420265E-03 +5.726219E-01 +6.723631E-02 +1.646287E-01 +5.489910E-03 +5.663491E-02 +6.494568E-04 +9.923269E-02 +2.053123E-03 +6.809958E-02 +9.386140E-04 +4.850339E-02 +4.856214E-04 +8.169408E-02 +1.354206E-03 +5.474566E-02 +6.128204E-04 +1.208954E-01 +3.074917E-03 +3.393101E-01 +2.398955E-02 +1.358915E-01 +3.758521E-03 +1.347810E-01 +3.718991E-03 +3.404279E-01 +2.380238E-02 +1.315829E-01 +3.585161E-03 +4.754269E-02 +4.593490E-04 +7.591439E-02 +1.170776E-03 +5.121227E-02 +5.346435E-04 +tally 4: +2.194240E-01 +9.811952E-03 +2.718491E-01 +1.501408E-02 +5.823038E-01 +6.967280E-02 +5.955012E-01 +7.246929E-02 +2.496521E-01 +1.268521E-02 +2.623766E-01 +1.397567E-02 +3.475667E-01 +2.453188E-02 +3.290723E-01 +2.187833E-02 +1.123843E+00 +2.718539E-01 +1.060797E+00 +2.384769E-01 +3.373095E-01 +2.280788E-02 +3.563621E-01 +2.594087E-02 +3.442434E-01 +2.400708E-02 +3.319540E-01 +2.206590E-02 +1.051427E+00 +2.459968E-01 +1.065112E+00 +2.465202E-01 +3.182627E-01 +2.049613E-02 +3.724575E-01 +2.824300E-02 +2.274581E-01 +1.052009E-02 +2.558629E-01 +1.318940E-02 +5.434276E-01 +6.081708E-02 +5.623197E-01 +6.408372E-02 +2.320052E-01 +1.096050E-02 +2.380175E-01 +1.168927E-02 diff --git a/tests/regression_tests/filter_translations/test.py b/tests/regression_tests/filter_translations/test.py new file mode 100644 index 0000000000..d61667b014 --- /dev/null +++ b/tests/regression_tests/filter_translations/test.py @@ -0,0 +1,90 @@ +import numpy as np + +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def model(): + + model = openmc.model.Model() + + fuel = openmc.Material() + fuel.set_density('g/cm3', 10.0) + fuel.add_nuclide('U235', 1.0) + zr = openmc.Material() + zr.set_density('g/cm3', 1.0) + zr.add_nuclide('Zr90', 1.0) + model.materials.extend([fuel, zr]) + + box1 = openmc.model.rectangular_prism(10.0, 10.0) + box2 = openmc.model.rectangular_prism(20.0, 20.0, boundary_type='reflective') + top = openmc.ZPlane(z0=10.0, boundary_type='vacuum') + bottom = openmc.ZPlane(z0=-10.0, boundary_type='vacuum') + cell1 = openmc.Cell(fill=fuel, region=box1 & +bottom & -top) + cell2 = openmc.Cell(fill=zr, region=~box1 & box2 & +bottom & -top) + model.geometry = openmc.Geometry([cell1, cell2]) + + model.settings.batches = 5 + model.settings.inactive = 0 + model.settings.particles = 1000 + + translation = np.array((10, -5, 0)) + + llc = np.array([-9, -9, -9]) + urc = np.array([9, 9, 9]) + + mesh_dims = (3, 4, 5) + + filters = [] + + # un-translated meshes + reg_mesh = openmc.RegularMesh() + reg_mesh.dimension = mesh_dims + reg_mesh.lower_left = llc + reg_mesh.upper_right = urc + + filters.append(openmc.MeshFilter(reg_mesh)) + + recti_mesh = openmc.RectilinearMesh() + recti_mesh.x_grid = np.linspace(llc[0], urc[0], mesh_dims[0]) + recti_mesh.y_grid = np.linspace(llc[1], urc[1], mesh_dims[1]) + recti_mesh.z_grid = np.linspace(llc[2], urc[2], mesh_dims[2]) + + filters.append(openmc.MeshFilter(recti_mesh)) + + llc = np.array(llc - translation) + urc = np.array(urc - translation) + + # translated meshes + translated_reg_mesh = openmc.RegularMesh() + translated_reg_mesh.dimension = mesh_dims + translated_reg_mesh.lower_left = llc + translated_reg_mesh.upper_right = urc + + filters.append(openmc.MeshFilter(translated_reg_mesh)) + filters[-1].translation = translation + + translated_recti_mesh = openmc.RectilinearMesh() + translated_recti_mesh.x_grid = np.linspace(llc[0], urc[0], mesh_dims[0]) + translated_recti_mesh.y_grid = np.linspace(llc[1], urc[1], mesh_dims[1]) + translated_recti_mesh.z_grid = np.linspace(llc[2], urc[2], mesh_dims[2]) + + filters.append(openmc.MeshFilter(translated_recti_mesh)) + filters[-1].translation = translation + + # Create tallies + for f in filters: + tally = openmc.Tally() + tally.filters = [f] + tally.scores = ['total'] + model.tallies.append(tally) + + return model + + +def test_filter_mesh_translations(model): + harness = PyAPITestHarness('statepoint.5.h5', model) + harness.main() diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index 4a6d692555..6d8f567bc8 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -470,11 +470,17 @@ def test_regular_mesh(lib_init): assert isinstance(mesh, openmc.lib.RegularMesh) assert mesh_id == mesh.id + translation = (1.0, 2.0, 3.0) + mf = openmc.lib.MeshFilter(mesh) assert mf.mesh == mesh + mf.translation = translation + assert mf.translation == translation msf = openmc.lib.MeshSurfaceFilter(mesh) assert msf.mesh == mesh + msf.translation = translation + assert msf.translation == translation def test_rectilinear_mesh(lib_init):