From c1aba0cc81304b2bf5f0e781885016cfcb5bd831 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 16 Dec 2025 14:09:29 -0600 Subject: [PATCH 01/10] Compute bounding boxes during mesh material volume calculation --- include/openmc/capi.h | 2 + include/openmc/mesh.h | 78 +++++- openmc/lib/mesh.py | 32 ++- openmc/mesh.py | 116 +++++++- src/mesh.cpp | 496 +++++++++++++++++++++++++++++++++- tests/unit_tests/test_mesh.py | 69 +++++ 6 files changed, 767 insertions(+), 26 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index d8041ef4149..5b4731fdda5 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -117,6 +117,8 @@ int openmc_mesh_get_n_elements(int32_t index, size_t* n); int openmc_mesh_get_volumes(int32_t index, double* volumes); int openmc_mesh_material_volumes(int32_t index, int nx, int ny, int nz, int max_mats, int32_t* materials, double* volumes); +int openmc_mesh_material_volumes_bbox(int32_t index, int nx, int ny, int nz, + int max_mats, int32_t* materials, double* volumes, double* bboxes); int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh); int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh); int openmc_new_filter(const char* type, int32_t* index); diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 7ceb623dae1..72cc840b0f7 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -87,8 +87,12 @@ namespace detail { class MaterialVolumes { public: + MaterialVolumes(int32_t* mats, double* vols, double* bboxes, int table_size) + : materials_(mats), volumes_(vols), bboxes_(bboxes), table_size_(table_size) + {} + MaterialVolumes(int32_t* mats, double* vols, int table_size) - : materials_(mats), volumes_(vols), table_size_(table_size) + : MaterialVolumes(mats, vols, nullptr, table_size) {} //! Add volume for a given material in a mesh element @@ -99,6 +103,17 @@ class MaterialVolumes { void add_volume(int index_elem, int index_material, double volume); void add_volume_unsafe(int index_elem, int index_material, double volume); + //! Add volume and union a bounding box for a given material in a mesh element + // + //! \param[in] index_elem Index of the mesh element + //! \param[in] index_material Index of the material within the model + //! \param[in] volume Volume to add + //! \param[in] bbox Bounding box to union into the result + void add_volume_bbox( + int index_elem, int index_material, double volume, const BoundingBox& bbox); + void add_volume_bbox_unsafe( + int index_elem, int index_material, double volume, const BoundingBox& bbox); + // Accessors int32_t& materials(int i, int j) { return materials_[i * table_size_ + j]; } const int32_t& materials(int i, int j) const @@ -112,11 +127,23 @@ class MaterialVolumes { return volumes_[i * table_size_ + j]; } + double& bboxes(int i, int j, int k) + { + return bboxes_[(i * table_size_ + j) * 6 + k]; + } + const double& bboxes(int i, int j, int k) const + { + return bboxes_[(i * table_size_ + j) * 6 + k]; + } + + bool has_bboxes() const { return bboxes_ != nullptr; } + bool table_full() const { return table_full_; } private: int32_t* materials_; //!< material index (bins, table_size) double* volumes_; //!< volume in [cm^3] (bins, table_size) + double* bboxes_; //!< bounding boxes (bins, table_size, 6) int table_size_; //!< Size of hash table for each mesh element bool table_full_ {false}; //!< Whether the hash table is full }; @@ -164,6 +191,34 @@ class Mesh { virtual void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const = 0; + //! Determine which bins were crossed by a particle with segment start offsets + // + //! The default implementation assumes that the returned segments are + //! contiguous (i.e., no gaps along the track) and constructs `start` + //! from a cumulative sum of `lengths`. Meshes whose `bins_crossed` omits + //! gaps (e.g., cylindrical/spherical or unstructured meshes) should + //! override this method to provide correct `start` values. + // + //! \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 + //! \param[out] start Fractional start offset for each segment + virtual void bins_crossed_with_start(Position r0, Position r1, + const Direction& u, vector& bins, vector& lengths, + vector& start) const + { + this->bins_crossed(r0, r1, u, bins, lengths); + start.clear(); + start.reserve(lengths.size()); + double cumulative {0.0}; + for (double frac : lengths) { + start.push_back(cumulative); + cumulative += frac; + } + } + //! Determine which surface bins were crossed by a particle // //! \param[in] r0 Previous position of the particle @@ -239,6 +294,19 @@ class Mesh { void material_volumes(int nx, int ny, int nz, int max_materials, int32_t* materials, double* volumes) const; + //! Determine volume and bounding boxes of materials within each mesh element + // + //! \param[in] nx Number of samples in x direction + //! \param[in] ny Number of samples in y direction + //! \param[in] nz Number of samples in z direction + //! \param[in] max_materials Maximum number of materials in a single mesh + //! element + //! \param[inout] materials Array storing material indices + //! \param[inout] volumes Array storing volumes + //! \param[inout] bboxes Array storing bounding boxes (n_elems, table_size, 6) + void material_volumes(int nx, int ny, int nz, int max_materials, + int32_t* materials, double* volumes, double* bboxes) const; + //! Determine bounding box of mesh // //! \return Bounding box of mesh @@ -297,6 +365,10 @@ class StructuredMesh : public Mesh { void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const override; + void bins_crossed_with_start(Position r0, Position r1, const Direction& u, + vector& bins, vector& lengths, + vector& start) const override; + void surface_bins_crossed(Position r0, Position r1, const Direction& u, vector& bins) const override; @@ -803,6 +875,10 @@ class MOABMesh : public UnstructuredMesh { void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const override; + void bins_crossed_with_start(Position r0, Position r1, const Direction& u, + vector& bins, vector& lengths, + vector& start) const override; + int get_bin(Position r) const override; int n_bins() const override; diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py index 3a720f98a31..5442620bd87 100644 --- a/openmc/lib/mesh.py +++ b/openmc/lib/mesh.py @@ -24,6 +24,7 @@ arr_2d_int32 = np.ctypeslib.ndpointer(dtype=np.int32, ndim=2, flags='CONTIGUOUS') arr_2d_double = np.ctypeslib.ndpointer(dtype=np.double, ndim=2, flags='CONTIGUOUS') +arr_3d_double = np.ctypeslib.ndpointer(dtype=np.double, ndim=3, flags='CONTIGUOUS') # Mesh functions _dll.openmc_extend_meshes.argtypes = [c_int32, c_char_p, POINTER(c_int32), @@ -50,6 +51,11 @@ c_int32, c_int, c_int, c_int, c_int, arr_2d_int32, arr_2d_double] _dll.openmc_mesh_material_volumes.restype = c_int _dll.openmc_mesh_material_volumes.errcheck = _error_handler +_dll.openmc_mesh_material_volumes_bbox.argtypes = [ + c_int32, c_int, c_int, c_int, c_int, arr_2d_int32, arr_2d_double, + arr_3d_double] +_dll.openmc_mesh_material_volumes_bbox.restype = c_int +_dll.openmc_mesh_material_volumes_bbox.errcheck = _error_handler _dll.openmc_mesh_get_plot_bins.argtypes = [ c_int32, _Position, _Position, c_int, POINTER(c_int), POINTER(c_int32) ] @@ -188,6 +194,7 @@ def material_volumes( n_samples: int | tuple[int, int, int] = 10_000, max_materials: int = 4, output: bool = True, + bounding_boxes: bool = False, ) -> MeshMaterialVolumes: """Determine volume of materials in each mesh element. @@ -213,6 +220,11 @@ def material_volumes( Estimated maximum number of materials in any given mesh element. output : bool, optional Whether or not to show output. + bounding_boxes : bool, optional + Whether or not to compute an axis-aligned bounding box for each + (mesh element, material) combination. When enabled, the bounding + box encloses the ray-estimator prisms used for the volume + estimation. Returns ------- @@ -243,23 +255,37 @@ def material_volumes( table_size = slot_factor*max_materials materials = np.full((n, table_size), EMPTY_SLOT, dtype=np.int32) volumes = np.zeros((n, table_size), dtype=np.float64) + bboxes = None + if bounding_boxes: + bboxes = np.empty((n, table_size, 6), dtype=np.float64) + bboxes[..., 0:3] = np.inf + bboxes[..., 3:6] = -np.inf # Run material volume calculation while True: try: with quiet_dll(output): - _dll.openmc_mesh_material_volumes( - self._index, nx, ny, nz, table_size, materials, volumes) + if bounding_boxes: + _dll.openmc_mesh_material_volumes_bbox( + self._index, nx, ny, nz, table_size, materials, + volumes, bboxes) + else: + _dll.openmc_mesh_material_volumes( + self._index, nx, ny, nz, table_size, materials, volumes) except AllocationError: # Increase size of result array and try again table_size *= 2 materials = np.full((n, table_size), EMPTY_SLOT, dtype=np.int32) volumes = np.zeros((n, table_size), dtype=np.float64) + if bounding_boxes: + bboxes = np.empty((n, table_size, 6), dtype=np.float64) + bboxes[..., 0:3] = np.inf + bboxes[..., 3:6] = -np.inf else: # If no error, break out of loop break - return MeshMaterialVolumes(materials, volumes) + return MeshMaterialVolumes(materials, volumes, bboxes) def get_plot_bins( self, diff --git a/openmc/mesh.py b/openmc/mesh.py index ce5218b5e97..6524c1823db 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -17,6 +17,7 @@ import openmc.checkvalue as cv from openmc.checkvalue import PathLike from openmc.utility_funcs import change_directory +from .bounding_box import BoundingBox from ._xml import get_elem_list, get_text from .mixin import IDManagerMixin from .surface import _BOUNDARY_TYPES @@ -40,6 +41,11 @@ class MeshMaterialVolumes(Mapping): Array of shape (elements, max_materials) storing material IDs volumes : numpy.ndarray Array of shape (elements, max_materials) storing material volumes + bboxes : numpy.ndarray, optional + Array of shape (elements, max_materials, 6) storing axis-aligned + bounding boxes for each (element, material) combination with ordering + (xmin, ymin, zmin, xmax, ymax, zmax). Bounding boxes enclose the + ray-estimator prisms used to compute volumes. See Also -------- @@ -64,9 +70,30 @@ class MeshMaterialVolumes(Mapping): [(2, 31.87963824195591), (1, 6.129949130817542)] """ - def __init__(self, materials: np.ndarray, volumes: np.ndarray): + def __init__( + self, + materials: np.ndarray, + volumes: np.ndarray, + bboxes: np.ndarray | None = None + ): self._materials = materials self._volumes = volumes + self._bboxes = bboxes + + if self._bboxes is not None: + if self._bboxes.shape[:2] != self._materials.shape: + raise ValueError( + 'bboxes must have shape (elements, max_materials, 6) ' + 'matching materials/volumes.' + ) + if self._bboxes.shape[2] != 6: + raise ValueError( + 'bboxes must have shape (elements, max_materials, 6).' + ) + + @property + def has_bounding_boxes(self) -> bool: + return self._bboxes is not None @property def num_elements(self) -> int: @@ -92,7 +119,11 @@ def __getitem__(self, material_id: int) -> np.ndarray: volumes[indices] = self._volumes[indices, i] return volumes - def by_element(self, index_elem: int) -> list[tuple[int | None, float]]: + def by_element( + self, + index_elem: int, + include_bboxes: bool = False + ) -> list[tuple[int | None, float] | tuple[int | None, float, BoundingBox | None]]: """Get a list of volumes for each material within a specific element. Parameters @@ -102,15 +133,66 @@ def by_element(self, index_elem: int) -> list[tuple[int | None, float]]: Returns ------- - list of tuple of (material ID, volume) + list of tuple + If ``include_bboxes`` is False (default), returns tuples of + (material ID, volume). If ``include_bboxes`` is True, returns + tuples of (material ID, volume, bounding box). """ table_size = self._volumes.shape[1] - return [ - (m if m > -1 else None, self._volumes[index_elem, i]) - for i in range(table_size) - if (m := self._materials[index_elem, i]) != -2 - ] + if include_bboxes and self._bboxes is None: + raise ValueError('Bounding boxes were not computed for this object.') + + results = [] + for i in range(table_size): + m = self._materials[index_elem, i] + if m == -2: + continue + mat_id = m if m > -1 else None + vol = self._volumes[index_elem, i] + + if include_bboxes: + bbox = None + vals = self._bboxes[index_elem, i] + if (vals[0] <= vals[3]) and (vals[1] <= vals[4]) and (vals[2] <= vals[5]): + bbox = BoundingBox(vals[0:3], vals[3:6]) + results.append((mat_id, vol, bbox)) + else: + results.append((mat_id, vol)) + + return results + + def bounding_box( + self, + material_id: int | None, + index_elem: int + ) -> BoundingBox | None: + """Get the bounding box for a (material, element) pair. + + Parameters + ---------- + material_id : int or None + Material ID. Use None for void. + index_elem : int + Mesh element index + + Returns + ------- + openmc.BoundingBox or None + Bounding box if present; otherwise None. + + """ + if self._bboxes is None: + raise ValueError('Bounding boxes were not computed for this object.') + + mat = -1 if material_id is None else material_id + for i in range(self._materials.shape[1]): + if self._materials[index_elem, i] == mat: + vals = self._bboxes[index_elem, i] + if (vals[0] <= vals[3]) and (vals[1] <= vals[4]) and (vals[2] <= vals[5]): + return BoundingBox(vals[0:3], vals[3:6]) + return None + return None def save(self, filename: PathLike): """Save material volumes to a .npz file. @@ -120,8 +202,10 @@ def save(self, filename: PathLike): filename : path-like Filename where data will be saved """ - np.savez_compressed( - filename, materials=self._materials, volumes=self._volumes) + kwargs = {'materials': self._materials, 'volumes': self._volumes} + if self._bboxes is not None: + kwargs['bboxes'] = self._bboxes + np.savez_compressed(filename, **kwargs) @classmethod def from_npz(cls, filename: PathLike) -> MeshMaterialVolumes: @@ -134,7 +218,8 @@ def from_npz(cls, filename: PathLike) -> MeshMaterialVolumes: """ filedata = np.load(filename) - return cls(filedata['materials'], filedata['volumes']) + bboxes = filedata['bboxes'] if 'bboxes' in filedata.files else None + return cls(filedata['materials'], filedata['volumes'], bboxes) class MeshBase(IDManagerMixin, ABC): @@ -366,6 +451,7 @@ def material_volumes( model: openmc.Model, n_samples: int | tuple[int, int, int] = 10_000, max_materials: int = 4, + bounding_boxes: bool = False, **kwargs ) -> MeshMaterialVolumes: """Determine volume of materials in each mesh element. @@ -388,6 +474,11 @@ def material_volumes( the x, y, and z dimensions. max_materials : int, optional Estimated maximum number of materials in any given mesh element. + bounding_boxes : bool, optional + Whether to compute an axis-aligned bounding box for each + (mesh element, material) combination. When enabled, the bounding + box encloses the ray-estimator prisms used for the volume + estimation. **kwargs : dict Keyword arguments passed to :func:`openmc.lib.init` @@ -419,7 +510,8 @@ def material_volumes( # Compute material volumes volumes = mesh.material_volumes( - n_samples, max_materials, output=kwargs['output']) + n_samples, max_materials, output=kwargs['output'], + bounding_boxes=bounding_boxes) # Restore original tallies model.tallies = original_tallies diff --git a/src/mesh.cpp b/src/mesh.cpp index 9a7a8e7517b..558d273a171 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1,6 +1,8 @@ #include "openmc/mesh.h" #include // for copy, equal, min, min_element #include +#include // for uint64_t +#include // for memcpy #define _USE_MATH_DEFINES // to make M_PI declared in Intel and MSVC compilers #include // for ceil #include // for size_t @@ -140,6 +142,96 @@ inline bool atomic_cas_int32(int32_t* ptr, int32_t& expected, int32_t desired) #endif } +inline uint64_t double_to_uint64(double value) +{ + uint64_t bits; + std::memcpy(&bits, &value, sizeof(value)); + return bits; +} + +inline double uint64_to_double(uint64_t bits) +{ + double value; + std::memcpy(&value, &bits, sizeof(value)); + return value; +} + +inline void atomic_min_double(double* ptr, double value) +{ +#if defined(__GNUC__) || defined(__clang__) + using may_alias_uint64_t [[gnu::may_alias]] = uint64_t; + auto* bits_ptr = reinterpret_cast(ptr); + uint64_t current_bits = __atomic_load_n(bits_ptr, __ATOMIC_SEQ_CST); + double current = uint64_to_double(current_bits); + while (value < current) { + uint64_t desired_bits = double_to_uint64(value); + uint64_t expected_bits = current_bits; + if (__atomic_compare_exchange_n(bits_ptr, &expected_bits, desired_bits, + false, __ATOMIC_SEQ_CST, __ATOMIC_SEQ_CST)) { + return; + } + current_bits = expected_bits; + current = uint64_to_double(current_bits); + } + +#elif defined(_MSC_VER) + auto* bits_ptr = reinterpret_cast(ptr); + long long current_bits = *bits_ptr; + double current = uint64_to_double(current_bits); + while (value < current) { + long long desired_bits = double_to_uint64(value); + long long old_bits = + _InterlockedCompareExchange64(bits_ptr, desired_bits, current_bits); + if (old_bits == current_bits) { + return; + } + current_bits = old_bits; + current = uint64_to_double(current_bits); + } + +#else +#error "No compare-and-swap implementation available for this compiler." +#endif +} + +inline void atomic_max_double(double* ptr, double value) +{ +#if defined(__GNUC__) || defined(__clang__) + using may_alias_uint64_t [[gnu::may_alias]] = uint64_t; + auto* bits_ptr = reinterpret_cast(ptr); + uint64_t current_bits = __atomic_load_n(bits_ptr, __ATOMIC_SEQ_CST); + double current = uint64_to_double(current_bits); + while (value > current) { + uint64_t desired_bits = double_to_uint64(value); + uint64_t expected_bits = current_bits; + if (__atomic_compare_exchange_n(bits_ptr, &expected_bits, desired_bits, + false, __ATOMIC_SEQ_CST, __ATOMIC_SEQ_CST)) { + return; + } + current_bits = expected_bits; + current = uint64_to_double(current_bits); + } + +#elif defined(_MSC_VER) + auto* bits_ptr = reinterpret_cast(ptr); + long long current_bits = *bits_ptr; + double current = uint64_to_double(current_bits); + while (value > current) { + long long desired_bits = double_to_uint64(value); + long long old_bits = + _InterlockedCompareExchange64(bits_ptr, desired_bits, current_bits); + if (old_bits == current_bits) { + return; + } + current_bits = old_bits; + current = uint64_to_double(current_bits); + } + +#else +#error "No compare-and-swap implementation available for this compiler." +#endif +} + namespace detail { //============================================================================== @@ -225,6 +317,123 @@ void MaterialVolumes::add_volume_unsafe( table_full_ = true; } +void MaterialVolumes::add_volume_bbox( + int index_elem, int index_material, double volume, const BoundingBox& bbox) +{ + // Loop for linear probing + for (int attempt = 0; attempt < table_size_; ++attempt) { + // Determine slot to check, making sure it is positive + int slot = (index_material + attempt) % table_size_; + if (slot < 0) + slot += table_size_; + int32_t* slot_ptr = &this->materials(index_elem, slot); + + // Non-atomic read of current material + int32_t current_val = *slot_ptr; + + // Found the desired material; accumulate volume and bbox + if (current_val == index_material) { +#pragma omp atomic + this->volumes(index_elem, slot) += volume; + if (bboxes_) { + atomic_min_double(&this->bboxes(index_elem, slot, 0), bbox.min.x); + atomic_min_double(&this->bboxes(index_elem, slot, 1), bbox.min.y); + atomic_min_double(&this->bboxes(index_elem, slot, 2), bbox.min.z); + atomic_max_double(&this->bboxes(index_elem, slot, 3), bbox.max.x); + atomic_max_double(&this->bboxes(index_elem, slot, 4), bbox.max.y); + atomic_max_double(&this->bboxes(index_elem, slot, 5), bbox.max.z); + } + return; + } + + // Slot appears to be empty; attempt to claim + if (current_val == EMPTY) { + // Attempt compare-and-swap from EMPTY to index_material + int32_t expected_val = EMPTY; + bool claimed_slot = + atomic_cas_int32(slot_ptr, expected_val, index_material); + + // If we claimed the slot or another thread claimed it but the same + // material was inserted, proceed to accumulate + if (claimed_slot || (expected_val == index_material)) { +#pragma omp atomic + this->volumes(index_elem, slot) += volume; + if (bboxes_) { + atomic_min_double(&this->bboxes(index_elem, slot, 0), bbox.min.x); + atomic_min_double(&this->bboxes(index_elem, slot, 1), bbox.min.y); + atomic_min_double(&this->bboxes(index_elem, slot, 2), bbox.min.z); + atomic_max_double(&this->bboxes(index_elem, slot, 3), bbox.max.x); + atomic_max_double(&this->bboxes(index_elem, slot, 4), bbox.max.y); + atomic_max_double(&this->bboxes(index_elem, slot, 5), bbox.max.z); + } + return; + } + } + } + + // If table is full, set a flag that can be checked later + table_full_ = true; +} + +void MaterialVolumes::add_volume_bbox_unsafe( + int index_elem, int index_material, double volume, const BoundingBox& bbox) +{ + // Linear probe + for (int attempt = 0; attempt < table_size_; ++attempt) { + // Determine slot to check, making sure it is positive + int slot = (index_material + attempt) % table_size_; + if (slot < 0) + slot += table_size_; + + // Read current material + int32_t current_val = this->materials(index_elem, slot); + + // Found the desired material; accumulate volume and bbox + if (current_val == index_material) { + this->volumes(index_elem, slot) += volume; + if (bboxes_) { + this->bboxes(index_elem, slot, 0) = + std::min(this->bboxes(index_elem, slot, 0), bbox.min.x); + this->bboxes(index_elem, slot, 1) = + std::min(this->bboxes(index_elem, slot, 1), bbox.min.y); + this->bboxes(index_elem, slot, 2) = + std::min(this->bboxes(index_elem, slot, 2), bbox.min.z); + this->bboxes(index_elem, slot, 3) = + std::max(this->bboxes(index_elem, slot, 3), bbox.max.x); + this->bboxes(index_elem, slot, 4) = + std::max(this->bboxes(index_elem, slot, 4), bbox.max.y); + this->bboxes(index_elem, slot, 5) = + std::max(this->bboxes(index_elem, slot, 5), bbox.max.z); + } + return; + } + + // Claim empty slot + if (current_val == EMPTY) { + this->materials(index_elem, slot) = index_material; + this->volumes(index_elem, slot) += volume; + if (bboxes_) { + this->bboxes(index_elem, slot, 0) = + std::min(this->bboxes(index_elem, slot, 0), bbox.min.x); + this->bboxes(index_elem, slot, 1) = + std::min(this->bboxes(index_elem, slot, 1), bbox.min.y); + this->bboxes(index_elem, slot, 2) = + std::min(this->bboxes(index_elem, slot, 2), bbox.min.z); + this->bboxes(index_elem, slot, 3) = + std::max(this->bboxes(index_elem, slot, 3), bbox.max.x); + this->bboxes(index_elem, slot, 4) = + std::max(this->bboxes(index_elem, slot, 4), bbox.max.y); + this->bboxes(index_elem, slot, 5) = + std::max(this->bboxes(index_elem, slot, 5), bbox.max.z); + } + return; + } + } + + // If table is full, set a flag that can be checked later + table_full_ = true; +} + } // namespace detail //============================================================================== @@ -333,6 +542,12 @@ vector Mesh::volumes() const void Mesh::material_volumes(int nx, int ny, int nz, int table_size, int32_t* materials, double* volumes) const +{ + this->material_volumes(nx, ny, nz, table_size, materials, volumes, nullptr); +} + +void Mesh::material_volumes(int nx, int ny, int nz, int table_size, + int32_t* materials, double* volumes, double* bboxes) const { if (mpi::master) { header("MESH MATERIAL VOLUMES CALCULATION", 7); @@ -351,7 +566,8 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, timer.start(); // Create object for keeping track of materials/volumes - detail::MaterialVolumes result(materials, volumes, table_size); + detail::MaterialVolumes result(materials, volumes, bboxes, table_size); + bool compute_bboxes = bboxes != nullptr; // Determine bounding box auto bbox = this->bounding_box(); @@ -370,8 +586,9 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, #pragma omp parallel { // Preallocate vector for mesh indices and length fractions and particle - std::vector bins; - std::vector length_fractions; + vector bins; + vector length_fractions; + vector start_fractions; Particle p; SourceSite site; @@ -446,7 +663,13 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, // Determine what mesh elements were crossed by particle bins.clear(); length_fractions.clear(); - this->bins_crossed(r0, p.r(), p.u(), bins, length_fractions); + if (compute_bboxes) { + start_fractions.clear(); + this->bins_crossed_with_start( + r0, p.r(), p.u(), bins, length_fractions, start_fractions); + } else { + this->bins_crossed(r0, p.r(), p.u(), bins, length_fractions); + } // Add volumes to any mesh elements that were crossed int i_material = p.material(); @@ -456,9 +679,32 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, for (int i_bin = 0; i_bin < bins.size(); i_bin++) { int mesh_index = bins[i_bin]; double length = distance * length_fractions[i_bin]; - - // Add volume to result - result.add_volume(mesh_index, i_material, length * d1 * d2); + double volume = length * d1 * d2; + + if (compute_bboxes) { + double axis_start = + r0[axis] + distance * start_fractions[i_bin]; + double axis_end = axis_start + length; + + Position contrib_min = site.r; + Position contrib_max = site.r; + + contrib_min[ax1] = site.r[ax1] - 0.5 * d1; + contrib_max[ax1] = site.r[ax1] + 0.5 * d1; + contrib_min[ax2] = site.r[ax2] - 0.5 * d2; + contrib_max[ax2] = site.r[ax2] + 0.5 * d2; + contrib_min[axis] = std::min(axis_start, axis_end); + contrib_max[axis] = std::max(axis_start, axis_end); + + BoundingBox contrib_bbox {contrib_min, contrib_max}; + contrib_bbox &= bbox; + + result.add_volume_bbox( + mesh_index, i_material, volume, contrib_bbox); + } else { + // Add volume to result + result.add_volume(mesh_index, i_material, volume); + } } if (distance == max_distance) @@ -505,10 +751,15 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, // Combine results from multiple MPI processes if (mpi::n_procs > 1) { int total = this->n_bins() * table_size; + int total_bbox = total * 6; if (mpi::master) { // Allocate temporary buffer for receiving data - std::vector mats(total); - std::vector vols(total); + vector mats(total); + vector vols(total); + vector recv_bboxes; + if (compute_bboxes) { + recv_bboxes.resize(total_bbox); + } for (int i = 1; i < mpi::n_procs; ++i) { // Receive material indices and volumes from process i @@ -516,6 +767,10 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, MPI_STATUS_IGNORE); MPI_Recv(vols.data(), total, MPI_DOUBLE, i, i, mpi::intracomm, MPI_STATUS_IGNORE); + if (compute_bboxes) { + MPI_Recv(recv_bboxes.data(), total_bbox, MPI_DOUBLE, i, i, + mpi::intracomm, MPI_STATUS_IGNORE); + } // Combine with existing results; we can call thread unsafe version of // add_volume because each thread is operating on a different element @@ -524,7 +779,18 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, for (int k = 0; k < table_size; ++k) { int index = index_elem * table_size + k; if (mats[index] != EMPTY) { - result.add_volume_unsafe(index_elem, mats[index], vols[index]); + if (compute_bboxes) { + int bbox_index = index * 6; + BoundingBox slot_bbox { + {recv_bboxes[bbox_index + 0], recv_bboxes[bbox_index + 1], + recv_bboxes[bbox_index + 2]}, + {recv_bboxes[bbox_index + 3], recv_bboxes[bbox_index + 4], + recv_bboxes[bbox_index + 5]}}; + result.add_volume_bbox_unsafe( + index_elem, mats[index], vols[index], slot_bbox); + } else { + result.add_volume_unsafe(index_elem, mats[index], vols[index]); + } } } } @@ -533,6 +799,9 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, // Send material indices and volumes to process 0 MPI_Send(materials, total, MPI_INT32_T, 0, mpi::rank, mpi::intracomm); MPI_Send(volumes, total, MPI_DOUBLE, 0, mpi::rank, mpi::intracomm); + if (compute_bboxes) { + MPI_Send(bboxes, total_bbox, MPI_DOUBLE, 0, mpi::rank, mpi::intracomm); + } } } @@ -1146,6 +1415,116 @@ void StructuredMesh::bins_crossed(Position r0, Position r1, const Direction& u, raytrace_mesh(r0, r1, u, TrackAggregator(this, bins, lengths)); } +void StructuredMesh::bins_crossed_with_start(Position r0, Position r1, + const Direction& u, vector& bins, vector& lengths, + vector& start) const +{ + bins.clear(); + lengths.clear(); + start.clear(); + + // Compute the length of the entire track. + double total_distance = (r1 - r0).norm(); + if (total_distance == 0.0 && settings::solver_type != SolverType::RANDOM_RAY) + return; + + // keep a copy of the original global position to pass to get_indices, + // which performs its own transformation to local coordinates + Position global_r = r0; + Position local_r = local_coords(r0); + + const int n = n_dimension_; + + // Flag if position is inside the mesh + bool in_mesh; + + // Position is r = r0 + u * traveled_distance, start at r0 + double traveled_distance {0.0}; + + // Calculate index of current cell. Offset the position a tiny bit in + // direction of flight + MeshIndex ijk = get_indices(global_r + TINY_BIT * u, in_mesh); + + // If track is very short, assume that it is completely inside one cell. + // Only the current cell will score. + if (total_distance < 2 * TINY_BIT) { + if (in_mesh) { + bins.push_back(get_bin_from_indices(ijk)); + lengths.push_back(1.0); + start.push_back(0.0); + } + return; + } + + // Calculate initial distances to next surfaces in all three dimensions + std::array distances; + for (int k = 0; k < n; ++k) { + distances[k] = distance_to_grid_boundary(ijk, k, local_r, u, 0.0); + } + + // Loop until r = r1 is eventually reached + while (true) { + if (in_mesh) { + // find surface with minimal distance to current position + const auto k = std::min_element(distances.begin(), distances.end()) - + distances.begin(); + + // Segment bounds along track + double seg_end = std::min(distances[k].distance, total_distance); + double seg_len = seg_end - traveled_distance; + if (seg_len > 0.0) { + bins.push_back(get_bin_from_indices(ijk)); + start.push_back(traveled_distance / total_distance); + lengths.push_back(seg_len / total_distance); + } + + // update position and leave, if we have reached end position + traveled_distance = distances[k].distance; + if (traveled_distance >= total_distance) + return; + + // Update cell and calculate distance to next surface in k-direction. + // The two other directions are still valid! + ijk[k] = distances[k].next_index; + distances[k] = + distance_to_grid_boundary(ijk, k, local_r, u, traveled_distance); + + // Check if we have left the interior of the mesh + in_mesh = ((ijk[k] >= 1) && (ijk[k] <= shape_[k])); + + } else { // not inside mesh + // For all directions outside the mesh, find the distance that we need + // to travel to reach the next surface. Use the largest distance, as + // only this will cross all outer surfaces. + int k_max {-1}; + for (int k = 0; k < n; ++k) { + if ((ijk[k] < 1 || ijk[k] > shape_[k]) && + (distances[k].distance > traveled_distance)) { + traveled_distance = distances[k].distance; + k_max = k; + } + } + + // Assure some distance is traveled + if (k_max == -1) { + traveled_distance += TINY_BIT; + } + + // If r1 is not inside the mesh, exit here + if (traveled_distance >= total_distance) + return; + + // Calculate the new cell index and update all distances to next + // surfaces. + ijk = get_indices(global_r + (traveled_distance + TINY_BIT) * u, in_mesh); + for (int k = 0; k < n; ++k) { + distances[k] = + distance_to_grid_boundary(ijk, k, local_r, u, traveled_distance); + } + } + } +} + void StructuredMesh::surface_bins_crossed( Position r0, Position r1, const Direction& u, vector& bins) const { @@ -2448,6 +2827,27 @@ extern "C" int openmc_mesh_material_volumes(int32_t index, int nx, int ny, return 0; } +extern "C" int openmc_mesh_material_volumes_bbox(int32_t index, int nx, int ny, + int nz, int table_size, int32_t* materials, double* volumes, double* bboxes) +{ + if (int err = check_mesh(index)) + return err; + + try { + model::meshes[index]->material_volumes( + nx, ny, nz, table_size, materials, volumes, bboxes); + } catch (const std::exception& e) { + set_errmsg(e.what()); + if (starts_with(e.what(), "Mesh")) { + return OPENMC_E_GEOMETRY; + } else { + return OPENMC_E_ALLOCATE; + } + } + + return 0; +} + extern "C" int openmc_mesh_get_plot_bins(int32_t index, Position origin, Position width, int basis, int* pixels, int32_t* data) { @@ -2977,6 +3377,82 @@ void MOABMesh::bins_crossed(Position r0, Position r1, const Direction& u, } }; +void MOABMesh::bins_crossed_with_start(Position r0, Position r1, + const Direction& u, vector& bins, vector& lengths, + vector& start) const +{ + moab::CartVect start_vec(r0.x, r0.y, r0.z); + moab::CartVect end_vec(r1.x, r1.y, r1.z); + moab::CartVect dir(u.x, u.y, u.z); + dir.normalize(); + + double track_len = (end_vec - start_vec).length(); + if (track_len == 0.0) + return; + + start_vec -= TINY_BIT * dir; + end_vec += TINY_BIT * dir; + + vector hits; + intersect_track(start_vec, dir, track_len, hits); + + bins.clear(); + lengths.clear(); + start.clear(); + + // if there are no intersections the track may lie entirely within a single + // tet. If this is the case, apply entire score to that tet and return. + if (hits.size() == 0) { + Position midpoint = r0 + u * (track_len * 0.5); + int bin = this->get_bin(midpoint); + if (bin != -1) { + bins.push_back(bin); + start.push_back(0.0); + lengths.push_back(1.0); + } + return; + } + + // for each segment in the set of tracks, try to look up a tet at the + // midpoint of the segment + Position current = r0; + double last_dist = 0.0; + for (const auto& hit : hits) { + double segment_start = last_dist; + double segment_length = hit - last_dist; + last_dist = hit; + + Position midpoint = current + u * (segment_length * 0.5); + int bin = this->get_bin(midpoint); + + // determine the start point for this segment + current = r0 + u * hit; + + if (bin == -1) { + continue; + } + + bins.push_back(bin); + start.push_back(segment_start / track_len); + lengths.push_back(segment_length / track_len); + } + + // tally remaining portion of track after last hit if 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) { + double segment_start = hits.back(); + double segment_length = track_len - hits.back(); + Position segment_start_pos = r0 + u * segment_start; + Position midpoint = segment_start_pos + u * (segment_length * 0.5); + int bin = this->get_bin(midpoint); + if (bin != -1) { + bins.push_back(bin); + start.push_back(segment_start / track_len); + lengths.push_back(segment_length / track_len); + } + } +} + moab::EntityHandle MOABMesh::get_tet(const Position& r) const { moab::CartVect pos(r.x, r.y, r.z); diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 9aca8b59656..19ec1940b30 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -691,6 +691,54 @@ def test_mesh_material_volumes_serialize(): assert new_volumes.by_element(3) == [(2, 1.0)] +def test_mesh_material_volumes_serialize_with_bboxes(): + materials = np.array([ + [1, -1, -2], + [-1, -2, -2], + [2, 1, -2], + [2, -2, -2] + ]) + volumes = np.array([ + [0.5, 0.5, 0.0], + [1.0, 0.0, 0.0], + [0.5, 0.5, 0.0], + [1.0, 0.0, 0.0] + ]) + + # (xmin, ymin, zmin, xmax, ymax, zmax) + bboxes = np.empty((4, 3, 6)) + bboxes[..., 0:3] = np.inf + bboxes[..., 3:6] = -np.inf + bboxes[0, 0] = [-1.0, -2.0, -3.0, 1.0, 2.0, 3.0] # material 1 + bboxes[0, 1] = [-5.0, -6.0, -7.0, 5.0, 6.0, 7.0] # void + bboxes[1, 0] = [0.0, 0.0, 0.0, 10.0, 1.0, 2.0] # void + bboxes[2, 0] = [-1.0, -1.0, -1.0, 0.0, 0.0, 0.0] # material 2 + bboxes[2, 1] = [0.0, 0.0, 0.0, 1.0, 1.0, 1.0] # material 1 + bboxes[3, 0] = [-2.0, -2.0, -2.0, 2.0, 2.0, 2.0] # material 2 + + mmv = openmc.MeshMaterialVolumes(materials, volumes, bboxes) + with TemporaryDirectory() as tmpdir: + path = f'{tmpdir}/volumes_bboxes.npz' + mmv.save(path) + loaded = openmc.MeshMaterialVolumes.from_npz(path) + + assert loaded.has_bounding_boxes + bb = loaded.bounding_box(1, 0) + assert isinstance(bb, openmc.BoundingBox) + np.testing.assert_array_equal(bb.lower_left, (-1.0, -2.0, -3.0)) + np.testing.assert_array_equal(bb.upper_right, (1.0, 2.0, 3.0)) + + bb_void = loaded.bounding_box(None, 0) + assert isinstance(bb_void, openmc.BoundingBox) + np.testing.assert_array_equal(bb_void.lower_left, (-5.0, -6.0, -7.0)) + np.testing.assert_array_equal(bb_void.upper_right, (5.0, 6.0, 7.0)) + + first = loaded.by_element(0, include_bboxes=True)[0][2] + assert isinstance(first, openmc.BoundingBox) + np.testing.assert_array_equal(first.lower_left, (-1.0, -2.0, -3.0)) + np.testing.assert_array_equal(first.upper_right, (1.0, 2.0, 3.0)) + + def test_mesh_material_volumes_boundary_conditions(sphere_model): """Test the material volumes method using a regular mesh that overlaps with a vacuum boundary condition.""" @@ -718,6 +766,27 @@ def test_mesh_material_volumes_boundary_conditions(sphere_model): assert evaluated[1] == pytest.approx(expected[1], rel=1e-2) +def test_mesh_material_volumes_bounding_boxes(sphere_model): + # Choose r_max such that the spherical mesh bounding box is fully contained + # in the geometry (the sphere_model has an outer sphere boundary at r=25). + mesh = openmc.SphericalMesh([0.0, 14.0], origin=(0.0, 0.0, 0.0)) + + mmv = mesh.material_volumes(sphere_model, (0, 50, 50), bounding_boxes=True) + assert mmv.has_bounding_boxes + + mesh_bb = mesh.bounding_box + found = False + for _, vol, bb in mmv.by_element(0, include_bboxes=True): + if vol > 0.0: + assert bb is not None + assert np.all(bb.lower_left >= mesh_bb.lower_left - 1e-12) + assert np.all(bb.upper_right <= mesh_bb.upper_right + 1e-12) + found = True + break + + assert found + + def test_raytrace_mesh_infinite_loop(run_in_tmpdir): # Create a model with one large spherical cell sphere = openmc.Sphere(r=100, boundary_type='vacuum') From de48189188e452b81789a269312953ad9005468c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 27 Dec 2025 16:11:14 -0600 Subject: [PATCH 02/10] Utilize element-material bounding boxes in R2SManager --- openmc/deplete/r2s.py | 45 +++++++++++++++++++------------------------ 1 file changed, 20 insertions(+), 25 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 97277cbda84..57bbe437ffb 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -269,6 +269,7 @@ def step1_neutron_transport( # Compute material volume fractions on the mesh if mat_vol_kwargs is None: mat_vol_kwargs = {} + mat_vol_kwargs.setdefault('bounding_boxes', True) self.results['mesh_material_volumes'] = mmv = comm.bcast( self.domains.material_volumes(self.neutron_model, **mat_vol_kwargs)) @@ -539,14 +540,15 @@ def step3_photon_transport( def get_decay_photon_source_mesh( self, time_index: int = -1 - ) -> list[openmc.MeshSource]: + ) -> list[openmc.IndependentSource]: """Create decay photon source for a mesh-based calculation. - This function creates N :class:`MeshSource` objects where N is the - maximum number of unique materials that appears in a single mesh - element. For each mesh element-material combination, and - IndependentSource instance is created with a spatial constraint limited - the sampled decay photons to the correct region. + For each mesh element-material combination, an + :class:`~openmc.IndependentSource` is created with a + :class:`~openmc.stats.Box` spatial distribution based on the bounding + box of the material within the mesh element. A material constraint is + also applied so that sampled source sites are limited to the correct + region. When the photon transport model is different from the neutron model, the photon MeshMaterialVolumes is used to determine whether an (element, @@ -559,19 +561,15 @@ def get_decay_photon_source_mesh( Returns ------- - list of openmc.MeshSource - A list of MeshSource objects, each containing IndependentSource - instances for the decay photons in the corresponding mesh element. + list of openmc.IndependentSource + A list of IndependentSource objects for the decay photons, one for + each mesh element-material combination with non-zero source strength. """ mat_dict = self.neutron_model._get_all_materials() - # Some MeshSource objects will have empty positions; create a "null source" - # that is used for this case - null_source = openmc.IndependentSource(particle='photon', strength=0.0) - - # List to hold sources for each MeshSource (length = N) - source_lists = [] + # List to hold all sources + sources = [] # Index in the overall list of activated materials index_mat = 0 @@ -594,7 +592,7 @@ def get_decay_photon_source_mesh( if mat_id is not None } - for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)): + for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True): # Skip void volume if mat_id is None: continue @@ -604,30 +602,27 @@ def get_decay_photon_source_mesh( index_mat += 1 continue - # Check whether a new MeshSource object is needed - if j >= len(source_lists): - source_lists.append([null_source]*n_elements) - # Get activated material composition original_mat = materials[index_mat] activated_mat = results[time_index].get_material(str(original_mat.id)) - # Create decay photon source source + # Create decay photon source energy = activated_mat.get_decay_photon_energy() if energy is not None: strength = energy.integral() - source_lists[j][index_elem] = openmc.IndependentSource( + space = openmc.stats.Box(*bbox) + sources.append(openmc.IndependentSource( + space=space, energy=energy, particle='photon', strength=strength, constraints={'domains': [mat_dict[mat_id]]} - ) + )) # Increment index of activated material index_mat += 1 - # Return list of mesh sources - return [openmc.MeshSource(self.domains, sources) for sources in source_lists] + return sources def load_results(self, path: PathLike): """Load results from a previous R2S calculation. From 6e9f39647e43da512a28ac261d0c2a9f65f9ba3c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 31 Dec 2025 13:42:43 -0500 Subject: [PATCH 03/10] Keep one version of openmc_mesh_material_volumes --- include/openmc/capi.h | 2 -- openmc/lib/mesh.py | 24 +++++++++--------------- src/mesh.cpp | 21 --------------------- 3 files changed, 9 insertions(+), 38 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 5b4731fdda5..749dc98a09f 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -116,8 +116,6 @@ int openmc_mesh_set_id(int32_t index, int32_t id); int openmc_mesh_get_n_elements(int32_t index, size_t* n); int openmc_mesh_get_volumes(int32_t index, double* volumes); int openmc_mesh_material_volumes(int32_t index, int nx, int ny, int nz, - int max_mats, int32_t* materials, double* volumes); -int openmc_mesh_material_volumes_bbox(int32_t index, int nx, int ny, int nz, int max_mats, int32_t* materials, double* volumes, double* bboxes); int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh); int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh); diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py index 5442620bd87..19e6f74d7ad 100644 --- a/openmc/lib/mesh.py +++ b/openmc/lib/mesh.py @@ -1,5 +1,5 @@ from collections.abc import Mapping, Sequence -from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER, +from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER, c_void_p, create_string_buffer, c_size_t) from math import sqrt import sys @@ -24,7 +24,6 @@ arr_2d_int32 = np.ctypeslib.ndpointer(dtype=np.int32, ndim=2, flags='CONTIGUOUS') arr_2d_double = np.ctypeslib.ndpointer(dtype=np.double, ndim=2, flags='CONTIGUOUS') -arr_3d_double = np.ctypeslib.ndpointer(dtype=np.double, ndim=3, flags='CONTIGUOUS') # Mesh functions _dll.openmc_extend_meshes.argtypes = [c_int32, c_char_p, POINTER(c_int32), @@ -48,14 +47,10 @@ _dll.openmc_mesh_bounding_box.restype = c_int _dll.openmc_mesh_bounding_box.errcheck = _error_handler _dll.openmc_mesh_material_volumes.argtypes = [ - c_int32, c_int, c_int, c_int, c_int, arr_2d_int32, arr_2d_double] + c_int32, c_int, c_int, c_int, c_int, arr_2d_int32, arr_2d_double, + c_void_p] _dll.openmc_mesh_material_volumes.restype = c_int _dll.openmc_mesh_material_volumes.errcheck = _error_handler -_dll.openmc_mesh_material_volumes_bbox.argtypes = [ - c_int32, c_int, c_int, c_int, c_int, arr_2d_int32, arr_2d_double, - arr_3d_double] -_dll.openmc_mesh_material_volumes_bbox.restype = c_int -_dll.openmc_mesh_material_volumes_bbox.errcheck = _error_handler _dll.openmc_mesh_get_plot_bins.argtypes = [ c_int32, _Position, _Position, c_int, POINTER(c_int), POINTER(c_int32) ] @@ -264,14 +259,13 @@ def material_volumes( # Run material volume calculation while True: try: + bboxes_ptr = None + if bboxes is not None: + bboxes_ptr = bboxes.ctypes.data_as(POINTER(c_double)) with quiet_dll(output): - if bounding_boxes: - _dll.openmc_mesh_material_volumes_bbox( - self._index, nx, ny, nz, table_size, materials, - volumes, bboxes) - else: - _dll.openmc_mesh_material_volumes( - self._index, nx, ny, nz, table_size, materials, volumes) + _dll.openmc_mesh_material_volumes( + self._index, nx, ny, nz, table_size, materials, + volumes, bboxes_ptr) except AllocationError: # Increase size of result array and try again table_size *= 2 diff --git a/src/mesh.cpp b/src/mesh.cpp index 558d273a171..3d12720a722 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2807,27 +2807,6 @@ extern "C" int openmc_mesh_bounding_box(int32_t index, double* ll, double* ur) } extern "C" int openmc_mesh_material_volumes(int32_t index, int nx, int ny, - int nz, int table_size, int32_t* materials, double* volumes) -{ - if (int err = check_mesh(index)) - return err; - - try { - model::meshes[index]->material_volumes( - nx, ny, nz, table_size, materials, volumes); - } catch (const std::exception& e) { - set_errmsg(e.what()); - if (starts_with(e.what(), "Mesh")) { - return OPENMC_E_GEOMETRY; - } else { - return OPENMC_E_ALLOCATE; - } - } - - return 0; -} - -extern "C" int openmc_mesh_material_volumes_bbox(int32_t index, int nx, int ny, int nz, int table_size, int32_t* materials, double* volumes, double* bboxes) { if (int err = check_mesh(index)) From 2c563f1b265ac79486d78246168bf21949eeb2fa Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 4 Jan 2026 15:38:03 -0600 Subject: [PATCH 04/10] Combine add_volume and add_volume_bbox --- include/openmc/mesh.h | 18 ++---- src/mesh.cpp | 143 +++++++++++------------------------------- 2 files changed, 40 insertions(+), 121 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 72cc840b0f7..2fdad26921e 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -100,19 +100,11 @@ class MaterialVolumes { //! \param[in] index_elem Index of the mesh element //! \param[in] index_material Index of the material within the model //! \param[in] volume Volume to add - void add_volume(int index_elem, int index_material, double volume); - void add_volume_unsafe(int index_elem, int index_material, double volume); - - //! Add volume and union a bounding box for a given material in a mesh element - // - //! \param[in] index_elem Index of the mesh element - //! \param[in] index_material Index of the material within the model - //! \param[in] volume Volume to add - //! \param[in] bbox Bounding box to union into the result - void add_volume_bbox( - int index_elem, int index_material, double volume, const BoundingBox& bbox); - void add_volume_bbox_unsafe( - int index_elem, int index_material, double volume, const BoundingBox& bbox); + //! \param[in] bbox Bounding box to union into the result (optional) + void add_volume(int index_elem, int index_material, double volume, + const BoundingBox* bbox = nullptr); + void add_volume_unsafe(int index_elem, int index_material, double volume, + const BoundingBox* bbox = nullptr); // Accessors int32_t& materials(int i, int j) { return materials_[i * table_size_ + j]; } diff --git a/src/mesh.cpp b/src/mesh.cpp index 3d12720a722..4144ada4833 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -239,7 +239,7 @@ namespace detail { //============================================================================== void MaterialVolumes::add_volume( - int index_elem, int index_material, double volume) + int index_elem, int index_material, double volume, const BoundingBox* bbox) { // This method handles adding elements to the materials hash table, // implementing open addressing with linear probing. Consistency across @@ -247,79 +247,6 @@ void MaterialVolumes::add_volume( // Ideally, we would use #pragma omp atomic compare, but it was introduced in // OpenMP 5.1 and is not widely supported yet. - // Loop for linear probing - for (int attempt = 0; attempt < table_size_; ++attempt) { - // Determine slot to check, making sure it is positive - int slot = (index_material + attempt) % table_size_; - if (slot < 0) - slot += table_size_; - int32_t* slot_ptr = &this->materials(index_elem, slot); - - // Non-atomic read of current material - int32_t current_val = *slot_ptr; - - // Found the desired material; accumulate volume - if (current_val == index_material) { -#pragma omp atomic - this->volumes(index_elem, slot) += volume; - return; - } - - // Slot appears to be empty; attempt to claim - if (current_val == EMPTY) { - // Attempt compare-and-swap from EMPTY to index_material - int32_t expected_val = EMPTY; - bool claimed_slot = - atomic_cas_int32(slot_ptr, expected_val, index_material); - - // If we claimed the slot or another thread claimed it but the same - // material was inserted, proceed to accumulate - if (claimed_slot || (expected_val == index_material)) { -#pragma omp atomic - this->volumes(index_elem, slot) += volume; - return; - } - } - } - - // If table is full, set a flag that can be checked later - table_full_ = true; -} - -void MaterialVolumes::add_volume_unsafe( - int index_elem, int index_material, double volume) -{ - // Linear probe - for (int attempt = 0; attempt < table_size_; ++attempt) { - // Determine slot to check, making sure it is positive - int slot = (index_material + attempt) % table_size_; - if (slot < 0) - slot += table_size_; - - // Read current material - int32_t current_val = this->materials(index_elem, slot); - - // Found the desired material; accumulate volume - if (current_val == index_material) { - this->volumes(index_elem, slot) += volume; - return; - } - - // Claim empty slot - if (current_val == EMPTY) { - this->materials(index_elem, slot) = index_material; - this->volumes(index_elem, slot) += volume; - return; - } - } - - // If table is full, set a flag that can be checked later - table_full_ = true; -} - -void MaterialVolumes::add_volume_bbox( - int index_elem, int index_material, double volume, const BoundingBox& bbox) -{ // Loop for linear probing for (int attempt = 0; attempt < table_size_; ++attempt) { // Determine slot to check, making sure it is positive @@ -335,13 +262,13 @@ void MaterialVolumes::add_volume_bbox( if (current_val == index_material) { #pragma omp atomic this->volumes(index_elem, slot) += volume; - if (bboxes_) { - atomic_min_double(&this->bboxes(index_elem, slot, 0), bbox.min.x); - atomic_min_double(&this->bboxes(index_elem, slot, 1), bbox.min.y); - atomic_min_double(&this->bboxes(index_elem, slot, 2), bbox.min.z); - atomic_max_double(&this->bboxes(index_elem, slot, 3), bbox.max.x); - atomic_max_double(&this->bboxes(index_elem, slot, 4), bbox.max.y); - atomic_max_double(&this->bboxes(index_elem, slot, 5), bbox.max.z); + if (bbox) { + atomic_min_double(&this->bboxes(index_elem, slot, 0), bbox->min.x); + atomic_min_double(&this->bboxes(index_elem, slot, 1), bbox->min.y); + atomic_min_double(&this->bboxes(index_elem, slot, 2), bbox->min.z); + atomic_max_double(&this->bboxes(index_elem, slot, 3), bbox->max.x); + atomic_max_double(&this->bboxes(index_elem, slot, 4), bbox->max.y); + atomic_max_double(&this->bboxes(index_elem, slot, 5), bbox->max.z); } return; } @@ -358,13 +285,13 @@ void MaterialVolumes::add_volume_bbox( if (claimed_slot || (expected_val == index_material)) { #pragma omp atomic this->volumes(index_elem, slot) += volume; - if (bboxes_) { - atomic_min_double(&this->bboxes(index_elem, slot, 0), bbox.min.x); - atomic_min_double(&this->bboxes(index_elem, slot, 1), bbox.min.y); - atomic_min_double(&this->bboxes(index_elem, slot, 2), bbox.min.z); - atomic_max_double(&this->bboxes(index_elem, slot, 3), bbox.max.x); - atomic_max_double(&this->bboxes(index_elem, slot, 4), bbox.max.y); - atomic_max_double(&this->bboxes(index_elem, slot, 5), bbox.max.z); + if (bbox) { + atomic_min_double(&this->bboxes(index_elem, slot, 0), bbox->min.x); + atomic_min_double(&this->bboxes(index_elem, slot, 1), bbox->min.y); + atomic_min_double(&this->bboxes(index_elem, slot, 2), bbox->min.z); + atomic_max_double(&this->bboxes(index_elem, slot, 3), bbox->max.x); + atomic_max_double(&this->bboxes(index_elem, slot, 4), bbox->max.y); + atomic_max_double(&this->bboxes(index_elem, slot, 5), bbox->max.z); } return; } @@ -375,8 +302,8 @@ void MaterialVolumes::add_volume_bbox( table_full_ = true; } -void MaterialVolumes::add_volume_bbox_unsafe( - int index_elem, int index_material, double volume, const BoundingBox& bbox) +void MaterialVolumes::add_volume_unsafe( + int index_elem, int index_material, double volume, const BoundingBox* bbox) { // Linear probe for (int attempt = 0; attempt < table_size_; ++attempt) { @@ -391,19 +318,19 @@ void MaterialVolumes::add_volume_bbox_unsafe( // Found the desired material; accumulate volume and bbox if (current_val == index_material) { this->volumes(index_elem, slot) += volume; - if (bboxes_) { + if (bbox) { this->bboxes(index_elem, slot, 0) = - std::min(this->bboxes(index_elem, slot, 0), bbox.min.x); + std::min(this->bboxes(index_elem, slot, 0), bbox->min.x); this->bboxes(index_elem, slot, 1) = - std::min(this->bboxes(index_elem, slot, 1), bbox.min.y); + std::min(this->bboxes(index_elem, slot, 1), bbox->min.y); this->bboxes(index_elem, slot, 2) = - std::min(this->bboxes(index_elem, slot, 2), bbox.min.z); + std::min(this->bboxes(index_elem, slot, 2), bbox->min.z); this->bboxes(index_elem, slot, 3) = - std::max(this->bboxes(index_elem, slot, 3), bbox.max.x); + std::max(this->bboxes(index_elem, slot, 3), bbox->max.x); this->bboxes(index_elem, slot, 4) = - std::max(this->bboxes(index_elem, slot, 4), bbox.max.y); + std::max(this->bboxes(index_elem, slot, 4), bbox->max.y); this->bboxes(index_elem, slot, 5) = - std::max(this->bboxes(index_elem, slot, 5), bbox.max.z); + std::max(this->bboxes(index_elem, slot, 5), bbox->max.z); } return; } @@ -412,19 +339,19 @@ void MaterialVolumes::add_volume_bbox_unsafe( if (current_val == EMPTY) { this->materials(index_elem, slot) = index_material; this->volumes(index_elem, slot) += volume; - if (bboxes_) { + if (bbox) { this->bboxes(index_elem, slot, 0) = - std::min(this->bboxes(index_elem, slot, 0), bbox.min.x); + std::min(this->bboxes(index_elem, slot, 0), bbox->min.x); this->bboxes(index_elem, slot, 1) = - std::min(this->bboxes(index_elem, slot, 1), bbox.min.y); + std::min(this->bboxes(index_elem, slot, 1), bbox->min.y); this->bboxes(index_elem, slot, 2) = - std::min(this->bboxes(index_elem, slot, 2), bbox.min.z); + std::min(this->bboxes(index_elem, slot, 2), bbox->min.z); this->bboxes(index_elem, slot, 3) = - std::max(this->bboxes(index_elem, slot, 3), bbox.max.x); + std::max(this->bboxes(index_elem, slot, 3), bbox->max.x); this->bboxes(index_elem, slot, 4) = - std::max(this->bboxes(index_elem, slot, 4), bbox.max.y); + std::max(this->bboxes(index_elem, slot, 4), bbox->max.y); this->bboxes(index_elem, slot, 5) = - std::max(this->bboxes(index_elem, slot, 5), bbox.max.z); + std::max(this->bboxes(index_elem, slot, 5), bbox->max.z); } return; } @@ -699,8 +626,8 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, BoundingBox contrib_bbox {contrib_min, contrib_max}; contrib_bbox &= bbox; - result.add_volume_bbox( - mesh_index, i_material, volume, contrib_bbox); + result.add_volume( + mesh_index, i_material, volume, &contrib_bbox); } else { // Add volume to result result.add_volume(mesh_index, i_material, volume); @@ -786,8 +713,8 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, recv_bboxes[bbox_index + 2]}, {recv_bboxes[bbox_index + 3], recv_bboxes[bbox_index + 4], recv_bboxes[bbox_index + 5]}}; - result.add_volume_bbox_unsafe( - index_elem, mats[index], vols[index], slot_bbox); + result.add_volume_unsafe( + index_elem, mats[index], vols[index], &slot_bbox); } else { result.add_volume_unsafe(index_elem, mats[index], vols[index]); } From 956d3be2c2b2f8bb490cae9020d73be457adc855 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 4 Jan 2026 15:52:18 -0600 Subject: [PATCH 05/10] Remove bins_crossed_with_start --- include/openmc/mesh.h | 36 -------- src/mesh.cpp | 200 +----------------------------------------- 2 files changed, 4 insertions(+), 232 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 2fdad26921e..852d301367e 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -183,34 +183,6 @@ class Mesh { virtual void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const = 0; - //! Determine which bins were crossed by a particle with segment start offsets - // - //! The default implementation assumes that the returned segments are - //! contiguous (i.e., no gaps along the track) and constructs `start` - //! from a cumulative sum of `lengths`. Meshes whose `bins_crossed` omits - //! gaps (e.g., cylindrical/spherical or unstructured meshes) should - //! override this method to provide correct `start` values. - // - //! \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 - //! \param[out] start Fractional start offset for each segment - virtual void bins_crossed_with_start(Position r0, Position r1, - const Direction& u, vector& bins, vector& lengths, - vector& start) const - { - this->bins_crossed(r0, r1, u, bins, lengths); - start.clear(); - start.reserve(lengths.size()); - double cumulative {0.0}; - for (double frac : lengths) { - start.push_back(cumulative); - cumulative += frac; - } - } - //! Determine which surface bins were crossed by a particle // //! \param[in] r0 Previous position of the particle @@ -357,10 +329,6 @@ class StructuredMesh : public Mesh { void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const override; - void bins_crossed_with_start(Position r0, Position r1, const Direction& u, - vector& bins, vector& lengths, - vector& start) const override; - void surface_bins_crossed(Position r0, Position r1, const Direction& u, vector& bins) const override; @@ -867,10 +835,6 @@ class MOABMesh : public UnstructuredMesh { void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const override; - void bins_crossed_with_start(Position r0, Position r1, const Direction& u, - vector& bins, vector& lengths, - vector& start) const override; - int get_bin(Position r) const override; int n_bins() const override; diff --git a/src/mesh.cpp b/src/mesh.cpp index 4144ada4833..7389a640573 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -515,7 +515,6 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, // Preallocate vector for mesh indices and length fractions and particle vector bins; vector length_fractions; - vector start_fractions; Particle p; SourceSite site; @@ -590,28 +589,23 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, // Determine what mesh elements were crossed by particle bins.clear(); length_fractions.clear(); - if (compute_bboxes) { - start_fractions.clear(); - this->bins_crossed_with_start( - r0, p.r(), p.u(), bins, length_fractions, start_fractions); - } else { - this->bins_crossed(r0, p.r(), p.u(), bins, length_fractions); - } + this->bins_crossed(r0, p.r(), p.u(), bins, length_fractions); // Add volumes to any mesh elements that were crossed int i_material = p.material(); if (i_material != C_NONE) { i_material = model::materials[i_material]->id(); } + double cumulative_frac = 0.0; for (int i_bin = 0; i_bin < bins.size(); i_bin++) { int mesh_index = bins[i_bin]; double length = distance * length_fractions[i_bin]; double volume = length * d1 * d2; if (compute_bboxes) { - double axis_start = - r0[axis] + distance * start_fractions[i_bin]; + double axis_start = r0[axis] + distance * cumulative_frac; double axis_end = axis_start + length; + cumulative_frac += length_fractions[i_bin]; Position contrib_min = site.r; Position contrib_max = site.r; @@ -1342,116 +1336,6 @@ void StructuredMesh::bins_crossed(Position r0, Position r1, const Direction& u, raytrace_mesh(r0, r1, u, TrackAggregator(this, bins, lengths)); } -void StructuredMesh::bins_crossed_with_start(Position r0, Position r1, - const Direction& u, vector& bins, vector& lengths, - vector& start) const -{ - bins.clear(); - lengths.clear(); - start.clear(); - - // Compute the length of the entire track. - double total_distance = (r1 - r0).norm(); - if (total_distance == 0.0 && settings::solver_type != SolverType::RANDOM_RAY) - return; - - // keep a copy of the original global position to pass to get_indices, - // which performs its own transformation to local coordinates - Position global_r = r0; - Position local_r = local_coords(r0); - - const int n = n_dimension_; - - // Flag if position is inside the mesh - bool in_mesh; - - // Position is r = r0 + u * traveled_distance, start at r0 - double traveled_distance {0.0}; - - // Calculate index of current cell. Offset the position a tiny bit in - // direction of flight - MeshIndex ijk = get_indices(global_r + TINY_BIT * u, in_mesh); - - // If track is very short, assume that it is completely inside one cell. - // Only the current cell will score. - if (total_distance < 2 * TINY_BIT) { - if (in_mesh) { - bins.push_back(get_bin_from_indices(ijk)); - lengths.push_back(1.0); - start.push_back(0.0); - } - return; - } - - // Calculate initial distances to next surfaces in all three dimensions - std::array distances; - for (int k = 0; k < n; ++k) { - distances[k] = distance_to_grid_boundary(ijk, k, local_r, u, 0.0); - } - - // Loop until r = r1 is eventually reached - while (true) { - if (in_mesh) { - // find surface with minimal distance to current position - const auto k = std::min_element(distances.begin(), distances.end()) - - distances.begin(); - - // Segment bounds along track - double seg_end = std::min(distances[k].distance, total_distance); - double seg_len = seg_end - traveled_distance; - if (seg_len > 0.0) { - bins.push_back(get_bin_from_indices(ijk)); - start.push_back(traveled_distance / total_distance); - lengths.push_back(seg_len / total_distance); - } - - // update position and leave, if we have reached end position - traveled_distance = distances[k].distance; - if (traveled_distance >= total_distance) - return; - - // Update cell and calculate distance to next surface in k-direction. - // The two other directions are still valid! - ijk[k] = distances[k].next_index; - distances[k] = - distance_to_grid_boundary(ijk, k, local_r, u, traveled_distance); - - // Check if we have left the interior of the mesh - in_mesh = ((ijk[k] >= 1) && (ijk[k] <= shape_[k])); - - } else { // not inside mesh - // For all directions outside the mesh, find the distance that we need - // to travel to reach the next surface. Use the largest distance, as - // only this will cross all outer surfaces. - int k_max {-1}; - for (int k = 0; k < n; ++k) { - if ((ijk[k] < 1 || ijk[k] > shape_[k]) && - (distances[k].distance > traveled_distance)) { - traveled_distance = distances[k].distance; - k_max = k; - } - } - - // Assure some distance is traveled - if (k_max == -1) { - traveled_distance += TINY_BIT; - } - - // If r1 is not inside the mesh, exit here - if (traveled_distance >= total_distance) - return; - - // Calculate the new cell index and update all distances to next - // surfaces. - ijk = get_indices(global_r + (traveled_distance + TINY_BIT) * u, in_mesh); - for (int k = 0; k < n; ++k) { - distances[k] = - distance_to_grid_boundary(ijk, k, local_r, u, traveled_distance); - } - } - } -} - void StructuredMesh::surface_bins_crossed( Position r0, Position r1, const Direction& u, vector& bins) const { @@ -3283,82 +3167,6 @@ void MOABMesh::bins_crossed(Position r0, Position r1, const Direction& u, } }; -void MOABMesh::bins_crossed_with_start(Position r0, Position r1, - const Direction& u, vector& bins, vector& lengths, - vector& start) const -{ - moab::CartVect start_vec(r0.x, r0.y, r0.z); - moab::CartVect end_vec(r1.x, r1.y, r1.z); - moab::CartVect dir(u.x, u.y, u.z); - dir.normalize(); - - double track_len = (end_vec - start_vec).length(); - if (track_len == 0.0) - return; - - start_vec -= TINY_BIT * dir; - end_vec += TINY_BIT * dir; - - vector hits; - intersect_track(start_vec, dir, track_len, hits); - - bins.clear(); - lengths.clear(); - start.clear(); - - // if there are no intersections the track may lie entirely within a single - // tet. If this is the case, apply entire score to that tet and return. - if (hits.size() == 0) { - Position midpoint = r0 + u * (track_len * 0.5); - int bin = this->get_bin(midpoint); - if (bin != -1) { - bins.push_back(bin); - start.push_back(0.0); - lengths.push_back(1.0); - } - return; - } - - // for each segment in the set of tracks, try to look up a tet at the - // midpoint of the segment - Position current = r0; - double last_dist = 0.0; - for (const auto& hit : hits) { - double segment_start = last_dist; - double segment_length = hit - last_dist; - last_dist = hit; - - Position midpoint = current + u * (segment_length * 0.5); - int bin = this->get_bin(midpoint); - - // determine the start point for this segment - current = r0 + u * hit; - - if (bin == -1) { - continue; - } - - bins.push_back(bin); - start.push_back(segment_start / track_len); - lengths.push_back(segment_length / track_len); - } - - // tally remaining portion of track after last hit if 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) { - double segment_start = hits.back(); - double segment_length = track_len - hits.back(); - Position segment_start_pos = r0 + u * segment_start; - Position midpoint = segment_start_pos + u * (segment_length * 0.5); - int bin = this->get_bin(midpoint); - if (bin != -1) { - bins.push_back(bin); - start.push_back(segment_start / track_len); - lengths.push_back(segment_length / track_len); - } - } -} - moab::EntityHandle MOABMesh::get_tet(const Position& r) const { moab::CartVect pos(r.x, r.y, r.z); From 97b6f089e5fd3c1fab507ea4116273a1c719220d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 11 Jan 2026 14:01:30 -0600 Subject: [PATCH 06/10] Remove bounding_box method on MeshMaterialVolumes --- openmc/mesh.py | 32 -------------------------------- tests/unit_tests/test_mesh.py | 15 +++++---------- 2 files changed, 5 insertions(+), 42 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 6524c1823db..e685cf46452 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -162,38 +162,6 @@ def by_element( return results - def bounding_box( - self, - material_id: int | None, - index_elem: int - ) -> BoundingBox | None: - """Get the bounding box for a (material, element) pair. - - Parameters - ---------- - material_id : int or None - Material ID. Use None for void. - index_elem : int - Mesh element index - - Returns - ------- - openmc.BoundingBox or None - Bounding box if present; otherwise None. - - """ - if self._bboxes is None: - raise ValueError('Bounding boxes were not computed for this object.') - - mat = -1 if material_id is None else material_id - for i in range(self._materials.shape[1]): - if self._materials[index_elem, i] == mat: - vals = self._bboxes[index_elem, i] - if (vals[0] <= vals[3]) and (vals[1] <= vals[4]) and (vals[2] <= vals[5]): - return BoundingBox(vals[0:3], vals[3:6]) - return None - return None - def save(self, filename: PathLike): """Save material volumes to a .npz file. diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 19ec1940b30..7cf15a382d4 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -723,21 +723,16 @@ def test_mesh_material_volumes_serialize_with_bboxes(): loaded = openmc.MeshMaterialVolumes.from_npz(path) assert loaded.has_bounding_boxes - bb = loaded.bounding_box(1, 0) - assert isinstance(bb, openmc.BoundingBox) - np.testing.assert_array_equal(bb.lower_left, (-1.0, -2.0, -3.0)) - np.testing.assert_array_equal(bb.upper_right, (1.0, 2.0, 3.0)) - - bb_void = loaded.bounding_box(None, 0) - assert isinstance(bb_void, openmc.BoundingBox) - np.testing.assert_array_equal(bb_void.lower_left, (-5.0, -6.0, -7.0)) - np.testing.assert_array_equal(bb_void.upper_right, (5.0, 6.0, 7.0)) - first = loaded.by_element(0, include_bboxes=True)[0][2] assert isinstance(first, openmc.BoundingBox) np.testing.assert_array_equal(first.lower_left, (-1.0, -2.0, -3.0)) np.testing.assert_array_equal(first.upper_right, (1.0, 2.0, 3.0)) + second = loaded.by_element(0, include_bboxes=True)[1][2] + assert isinstance(second, openmc.BoundingBox) + np.testing.assert_array_equal(second.lower_left, (-5.0, -6.0, -7.0)) + np.testing.assert_array_equal(second.upper_right, (5.0, 6.0, 7.0)) + def test_mesh_material_volumes_boundary_conditions(sphere_model): """Test the material volumes method using a regular mesh From 198c2b5bd005709869ed87b39dc36688cfd1c4c8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 11 Jan 2026 14:05:32 -0600 Subject: [PATCH 07/10] Simplification in MeshMaterialVolumes.by_element --- openmc/mesh.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index e685cf46452..296a9c1f53b 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -152,10 +152,8 @@ def by_element( vol = self._volumes[index_elem, i] if include_bboxes: - bbox = None vals = self._bboxes[index_elem, i] - if (vals[0] <= vals[3]) and (vals[1] <= vals[4]) and (vals[2] <= vals[5]): - bbox = BoundingBox(vals[0:3], vals[3:6]) + bbox = BoundingBox(vals[0:3], vals[3:6]) results.append((mat_id, vol, bbox)) else: results.append((mat_id, vol)) From adaab31ff242930ab865c5222a145697dad66d55 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Jan 2026 22:54:04 -0600 Subject: [PATCH 08/10] Combine implementation of atomic_min/max_double --- src/mesh.cpp | 44 ++++++++------------------------------------ 1 file changed, 8 insertions(+), 36 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 7389a640573..c61cdcfce2b 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -156,14 +156,14 @@ inline double uint64_to_double(uint64_t bits) return value; } -inline void atomic_min_double(double* ptr, double value) +inline void atomic_update_double(double* ptr, double value, bool is_min) { #if defined(__GNUC__) || defined(__clang__) using may_alias_uint64_t [[gnu::may_alias]] = uint64_t; auto* bits_ptr = reinterpret_cast(ptr); uint64_t current_bits = __atomic_load_n(bits_ptr, __ATOMIC_SEQ_CST); double current = uint64_to_double(current_bits); - while (value < current) { + while (is_min ? (value < current) : (value > current)) { uint64_t desired_bits = double_to_uint64(value); uint64_t expected_bits = current_bits; if (__atomic_compare_exchange_n(bits_ptr, &expected_bits, desired_bits, @@ -178,7 +178,7 @@ inline void atomic_min_double(double* ptr, double value) auto* bits_ptr = reinterpret_cast(ptr); long long current_bits = *bits_ptr; double current = uint64_to_double(current_bits); - while (value < current) { + while (is_min ? (value < current) : (value > current)) { long long desired_bits = double_to_uint64(value); long long old_bits = _InterlockedCompareExchange64(bits_ptr, desired_bits, current_bits); @@ -196,40 +196,12 @@ inline void atomic_min_double(double* ptr, double value) inline void atomic_max_double(double* ptr, double value) { -#if defined(__GNUC__) || defined(__clang__) - using may_alias_uint64_t [[gnu::may_alias]] = uint64_t; - auto* bits_ptr = reinterpret_cast(ptr); - uint64_t current_bits = __atomic_load_n(bits_ptr, __ATOMIC_SEQ_CST); - double current = uint64_to_double(current_bits); - while (value > current) { - uint64_t desired_bits = double_to_uint64(value); - uint64_t expected_bits = current_bits; - if (__atomic_compare_exchange_n(bits_ptr, &expected_bits, desired_bits, - false, __ATOMIC_SEQ_CST, __ATOMIC_SEQ_CST)) { - return; - } - current_bits = expected_bits; - current = uint64_to_double(current_bits); - } - -#elif defined(_MSC_VER) - auto* bits_ptr = reinterpret_cast(ptr); - long long current_bits = *bits_ptr; - double current = uint64_to_double(current_bits); - while (value > current) { - long long desired_bits = double_to_uint64(value); - long long old_bits = - _InterlockedCompareExchange64(bits_ptr, desired_bits, current_bits); - if (old_bits == current_bits) { - return; - } - current_bits = old_bits; - current = uint64_to_double(current_bits); - } + atomic_update_double(ptr, value, false); +} -#else -#error "No compare-and-swap implementation available for this compiler." -#endif +inline void atomic_min_double(double* ptr, double value) +{ + atomic_update_double(ptr, value, true); } namespace detail { From cf40820f3612eb31458bd18ee7ce85da4d09486d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Jan 2026 22:59:18 -0600 Subject: [PATCH 09/10] Simplify bit casting --- src/mesh.cpp | 29 ++++++++++++----------------- 1 file changed, 12 insertions(+), 17 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index c61cdcfce2b..b84a45a6a87 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -142,18 +142,13 @@ inline bool atomic_cas_int32(int32_t* ptr, int32_t& expected, int32_t desired) #endif } -inline uint64_t double_to_uint64(double value) +// Helper function equivalent to std::bit_cast in C++20 +template +inline To bit_cast_value(const From& value) { - uint64_t bits; - std::memcpy(&bits, &value, sizeof(value)); - return bits; -} - -inline double uint64_to_double(uint64_t bits) -{ - double value; - std::memcpy(&value, &bits, sizeof(value)); - return value; + To out; + std::memcpy(&out, &value, sizeof(To)); + return out; } inline void atomic_update_double(double* ptr, double value, bool is_min) @@ -162,31 +157,31 @@ inline void atomic_update_double(double* ptr, double value, bool is_min) using may_alias_uint64_t [[gnu::may_alias]] = uint64_t; auto* bits_ptr = reinterpret_cast(ptr); uint64_t current_bits = __atomic_load_n(bits_ptr, __ATOMIC_SEQ_CST); - double current = uint64_to_double(current_bits); + double current = bit_cast_value(current_bits); while (is_min ? (value < current) : (value > current)) { - uint64_t desired_bits = double_to_uint64(value); + uint64_t desired_bits = bit_cast_value(value); uint64_t expected_bits = current_bits; if (__atomic_compare_exchange_n(bits_ptr, &expected_bits, desired_bits, false, __ATOMIC_SEQ_CST, __ATOMIC_SEQ_CST)) { return; } current_bits = expected_bits; - current = uint64_to_double(current_bits); + current = bit_cast_value(current_bits); } #elif defined(_MSC_VER) auto* bits_ptr = reinterpret_cast(ptr); long long current_bits = *bits_ptr; - double current = uint64_to_double(current_bits); + double current = bit_cast_value(current_bits); while (is_min ? (value < current) : (value > current)) { - long long desired_bits = double_to_uint64(value); + long long desired_bits = bit_cast_value(value); long long old_bits = _InterlockedCompareExchange64(bits_ptr, desired_bits, current_bits); if (old_bits == current_bits) { return; } current_bits = old_bits; - current = uint64_to_double(current_bits); + current = bit_cast_value(current_bits); } #else From dc515205c0a58a1910c017b394bd76619b8ffd15 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 15 Jan 2026 00:09:50 -0600 Subject: [PATCH 10/10] Update MMV bounding box test --- tests/unit_tests/test_mesh.py | 60 +++++++++++++++++++++++++---------- 1 file changed, 44 insertions(+), 16 deletions(-) diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 7cf15a382d4..c5855a7b051 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -1,6 +1,8 @@ from math import pi from tempfile import TemporaryDirectory from pathlib import Path +import itertools +import random import h5py import numpy as np @@ -761,25 +763,51 @@ def test_mesh_material_volumes_boundary_conditions(sphere_model): assert evaluated[1] == pytest.approx(expected[1], rel=1e-2) -def test_mesh_material_volumes_bounding_boxes(sphere_model): - # Choose r_max such that the spherical mesh bounding box is fully contained - # in the geometry (the sphere_model has an outer sphere boundary at r=25). - mesh = openmc.SphericalMesh([0.0, 14.0], origin=(0.0, 0.0, 0.0)) +def test_mesh_material_volumes_bounding_boxes(): + # Create a model with 8 spherical cells at known locations with random radii + box = openmc.model.RectangularParallelepiped( + -10, 10, -10, 10, -10, 10, boundary_type='vacuum') - mmv = mesh.material_volumes(sphere_model, (0, 50, 50), bounding_boxes=True) + mat = openmc.Material() + mat.add_nuclide('H1', 1.0) + + sph_cells = [] + for x, y, z in itertools.product((-5., 5.), repeat=3): + mat_i = mat.clone() + sph = openmc.Sphere(x, y, z, r=random.uniform(0.5, 1.5)) + sph_cells.append(openmc.Cell(region=-sph, fill=mat_i)) + background = openmc.Cell(region=-box & openmc.Intersection([~c.region for c in sph_cells])) + + model = openmc.Model() + model.geometry = openmc.Geometry(sph_cells + [background]) + model.settings.particles = 1000 + model.settings.batches = 10 + + # Create a one-element mesh that encompasses the entire geometry + mesh = openmc.RegularMesh() + mesh.lower_left = (-10., -10., -10.) + mesh.upper_right = (10., 10., 10.) + mesh.dimension = (1, 1, 1) + + # Run material volume calculation with bounding boxes + n_samples = (400, 400, 400) + mmv = mesh.material_volumes(model, n_samples, max_materials=10, bounding_boxes=True) assert mmv.has_bounding_boxes - mesh_bb = mesh.bounding_box - found = False - for _, vol, bb in mmv.by_element(0, include_bboxes=True): - if vol > 0.0: - assert bb is not None - assert np.all(bb.lower_left >= mesh_bb.lower_left - 1e-12) - assert np.all(bb.upper_right <= mesh_bb.upper_right + 1e-12) - found = True - break - - assert found + # Create a mapping of material ID to bounding box + bbox_by_mat = { + mat_id: bbox + for mat_id, vol, bbox in mmv.by_element(0, include_bboxes=True) + if mat_id is not None and vol > 0.0 + } + + # Match the mesh ray spacing used for the bounding box estimator. + tol = 0.5 * mesh.bounding_box.width[0] / n_samples[0] + for cell in sph_cells: + bbox = bbox_by_mat[cell.fill.id] + cell_bbox = cell.bounding_box + np.testing.assert_allclose(bbox.lower_left, cell_bbox.lower_left, atol=tol) + np.testing.assert_allclose(bbox.upper_right, cell_bbox.upper_right, atol=tol) def test_raytrace_mesh_infinite_loop(run_in_tmpdir):