Reverting the sorting parts in filter and mesh python files

Author: Practice2001

openmc/filter.py

@@ -1112,20 +1112,12 @@ def from_volumes(cls, mesh: openmc.MeshBase, volumes: openmc.MeshMaterialVolumes
             A new MeshMaterialFilter instance
 
         """
-        # Build bins sorted by volume ascending (mat ID ascending as tiebreaker)
-        # so the bin ordering is deterministic regardless of the ray-traversal
-        # order stored in the raw _materials array.
-        bins = []
-        for i in range(volumes.num_elements):
-            entries = volumes.by_element(i)
-            # sort by (volume asc, mat_id asc), excluding void (None)
-            entries = sorted(
-                ((mat_id, vol)
-                 for mat_id, vol in entries if mat_id is not None),
-                key=lambda t: (t[1], t[0])
-            )
-            for mat_id, _ in entries:
-                bins.append((i, mat_id))
+        # Get flat arrays of material IDs and element indices
+        mat_ids = volumes._materials[volumes._materials > -1]
+        elems, _ = np.where(volumes._materials > -1)
+
+        # Stack them into a 2D array of (element, material) pairs
+        bins = np.column_stack((elems, mat_ids))
         return cls(mesh, bins)
 
     def __hash__(self):

openmc/mesh.py

@@ -506,33 +506,7 @@ def material_volumes(
 
         # Restore original tallies
         model.tallies = original_tallies
-
-        # Sort each element's materials by material ID (ascending), with
-        # None (-1) after positive IDs and empty slots (-2) last.  This
-        # gives a deterministic ordering that is independent of the ray
-        # traversal order used by the C library.
-        mat_arr = volumes._materials.copy()
-        vol_arr = volumes._volumes.copy()
-        bbox_arr = volumes._bboxes.copy() if volumes._bboxes is not None else None
-        n_elements, table_size = mat_arr.shape
-
-        def _sort_key(m):
-            m = int(m)
-            if m == -2:
-                return (2, 0)   # empty slot – always last
-            if m == -1:
-                return (1, 0)   # vacuum / None – just before empty
-            return (0, m)       # real material – ascending ID
-
-        for i in range(n_elements):
-            indices = sorted(range(table_size),
-                             key=lambda j: _sort_key(mat_arr[i, j]))
-            mat_arr[i] = mat_arr[i, indices]
-            vol_arr[i] = vol_arr[i, indices]
-            if bbox_arr is not None:
-                bbox_arr[i] = bbox_arr[i, indices]
-
-        return MeshMaterialVolumes(mat_arr, vol_arr, bbox_arr)
+        return volumes
 
 
 class StructuredMesh(MeshBase):
@@ -3366,6 +3340,7 @@ def append_dataset(dset, array):
                 )
 
         with h5py.File(filename, "w") as f:
+
             root = f.create_group("VTKHDF")
             vtk_file_format_version = (2, 1)
             root.attrs["Version"] = vtk_file_format_version

tests/unit_tests/conftest.py

@@ -1,20 +1,9 @@
-import os
-
 import openmc
 import pytest
 
 from tests.regression_tests import config
 
 
-@pytest.fixture
-def run_in_tmpdir(tmp_path):
-    """Run test in a temporary directory, restoring cwd afterwards."""
-    orig = os.getcwd()
-    os.chdir(tmp_path)
-    yield tmp_path
-    os.chdir(orig)
-
-
 @pytest.fixture(scope='module')
 def mpi_intracomm():
     if config['mpi']:

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
@@ -642,6 +644,7 @@ def test_mesh_get_homogenized_materials():
 
 @pytest.fixture
 def sphere_model():
+    openmc.reset_auto_ids()
     # Model with three materials separated by planes x=0 and z=0
     mats = []
     for i in range(3):
@@ -768,6 +771,49 @@ def test_mesh_material_volumes_serialize_with_bboxes():
     np.testing.assert_array_equal(second.upper_right, (5.0, 6.0, 7.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
+    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
     that overlaps with a vacuum boundary condition."""
@@ -841,6 +887,53 @@ def test_mesh_material_volumes_bounding_boxes():
         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_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')
+
+    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
+
+    # 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):
     # Create a model with one large spherical cell
     sphere = openmc.Sphere(r=100, boundary_type='vacuum')
@@ -953,10 +1046,6 @@ def test_filter_time_mesh(run_in_tmpdir):
     )
 
 
-# =============================================================================
-# VTKHDF Format Tests for StructuredMeshes
-# =============================================================================
-
 def test_regular_mesh_get_indices_at_coords():
     """Test get_indices_at_coords method for RegularMesh"""
     # Create a 10x10x10 mesh from (0,0,0) to (1,1,1)