From fb0a5cd9f3f29abfcc344f6c62607ccc76b58c7a Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 3 Oct 2026 08:53:49 +0000 Subject: [PATCH] Skip cells whose bounding boxes do not contain the point in cell searches A bounding box is computed for each CSG cell at initialization, enlarged slightly for roundoff, and the searches over the cells of a universe and over neighbor lists skip cells whose boxes do not contain the point before evaluating whether the cell contains it. Cells are checked in the same order as before, so the same cell is found, while expensive containment checks are replaced by a box test for cells that cannot contain the point. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014RN7JroEkrbkLKVHbMDap9 --- docs/source/methods/geometry.rst | 12 ++- include/openmc/cell.h | 13 +++ src/geometry.cpp | 7 +- src/geometry_aux.cpp | 21 +++++ src/universe.cpp | 5 +- tests/unit_tests/test_cell_search.py | 128 +++++++++++++++++++++++++++ 6 files changed, 181 insertions(+), 5 deletions(-) create mode 100644 tests/unit_tests/test_cell_search.py diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst index 05cda4b6423..bfd48d6ed63 100644 --- a/docs/source/methods/geometry.rst +++ b/docs/source/methods/geometry.rst @@ -537,6 +537,14 @@ lattice is determined, and then whatever universe fills that lattice position is recursively searched. The search ends once a cell containing a normal material is found that contains the specified point. +Checking whether a point is inside a cell requires evaluating the senses of the +point with respect to the surfaces of the cell. To avoid this for cells that +cannot contain the point, a bounding box is computed for each cell at +initialization, enlarged slightly to allow for roundoff in positions on the +boundary of the cell. During a search, a cell is only checked if its bounding +box contains the point. Since cells are still checked in the same order, the +same cell is found as without the bounding boxes. + .. _cell-contains: ---------------------- @@ -728,7 +736,9 @@ cell-based neighbor lists in OpenMC grow dynamically as particles are transported through the geometry and cross surfaces. Special care must be taken to ensure that these dynamic neighbor lists are populated in a threadsafe manner. Full details of the implementation in OpenMC can be found in a paper by -`Harper et al `_. +`Harper et al `_. When a neighbor +list is searched, cells whose bounding boxes do not contain the particle's +position are skipped, as described in :ref:`find-cell`. .. _reflection: diff --git a/include/openmc/cell.h b/include/openmc/cell.h index c2a09d774b5..05eb160c82c 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -400,6 +400,19 @@ class Cell { //! \brief Neighboring cells in the same universe. NeighborList neighbors_; + //! \brief Bounding box of the cell, enlarged slightly to allow for roundoff + //! in positions on its boundary, used to skip the cell when searching for + //! the cell containing a point outside of it + BoundingBox search_box_; + + //! \brief Whether a point may be in the cell based on its bounding box + bool may_contain(Position r) const + { + return r.x >= search_box_.min.x && r.x <= search_box_.max.x && + r.y >= search_box_.min.y && r.y <= search_box_.max.y && + r.z >= search_box_.min.z && r.z <= search_box_.max.z; + } + Position translation_ {0, 0, 0}; //!< Translation vector for filled universe //! \brief Rotational tranfsormation of the filled universe. diff --git a/src/geometry.cpp b/src/geometry.cpp index ecacf0bffbf..e17d51b1c62 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -137,14 +137,17 @@ bool find_cell_inner( // Make sure the search cell is in the same universe. int i_universe = p.lowest_coord().universe(); - if (model::cells[i_cell]->universe_ != i_universe) + const auto& c = *model::cells[i_cell]; + if (c.universe_ != i_universe) continue; // Check if this cell contains the particle. Position r {p.r_local()}; Direction u {p.u_local()}; auto surf = p.surface(); - if (model::cells[i_cell]->contains(r, u, surf)) { + if (!c.may_contain(r)) + continue; + if (c.contains(r, u, surf)) { p.lowest_coord().cell() = i_cell; found = true; break; diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index a740740c1e6..8b803632251 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -141,6 +141,26 @@ void adjust_indices() } } +//============================================================================== +//! Set the bounding boxes used to skip CSG cells that cannot contain a point +//! when searching for the cell containing it. + +void set_cell_search_boxes() +{ + for (auto& c : model::cells) { + if (c->geom_type() != GeometryType::CSG) + continue; + BoundingBox b = c->bounding_box(); + for (int i = 0; i < 3; ++i) { + double pad = + 1e-9 * std::max({1.0, std::abs(b.min[i]), std::abs(b.max[i])}); + b.min[i] -= pad; + b.max[i] += pad; + } + c->search_box_ = b; + } +} + //============================================================================== //! Partition some universes with many z-planes for faster find_cell searches. @@ -276,6 +296,7 @@ void finalize_geometry() adjust_indices(); count_universe_instances(); partition_universes(); + set_cell_search_boxes(); // Assign temperatures to cells that don't have temperatures already assigned assign_temperatures(); diff --git a/src/universe.cpp b/src/universe.cpp index 5d14308fe29..ff87b84533c 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -48,10 +48,11 @@ bool Universe::find_cell(GeometryState& p) const int32_t i_univ = p.lowest_coord().universe(); for (auto i_cell : cells) { - if (model::cells[i_cell]->universe_ != i_univ) + const auto& c = *model::cells[i_cell]; + if (c.universe_ != i_univ || !c.may_contain(r)) continue; // Check if this cell contains the particle - if (model::cells[i_cell]->contains(r, u, surf)) { + if (c.contains(r, u, surf)) { p.lowest_coord().cell() = i_cell; return true; } diff --git a/tests/unit_tests/test_cell_search.py b/tests/unit_tests/test_cell_search.py new file mode 100644 index 00000000000..5caf31d56ba --- /dev/null +++ b/tests/unit_tests/test_cell_search.py @@ -0,0 +1,128 @@ +import numpy as np +import pytest +import openmc +import openmc.lib + + +@pytest.fixture +def many_cells_model(): + """Spheres, finite rods and their overlaps in a box, with a background + cell outside of all of them, so that most cells are skipped based on their + bounding boxes when searching for the cell containing a point.""" + rng = np.random.default_rng(1) + mat = openmc.Material() + mat.add_nuclide('H1', 1.0) + mat.set_density('g/cm3', 0.1) + + box = openmc.model.RectangularParallelepiped( + -10, 10, -10, 10, -10, 10, boundary_type='vacuum') + objects = [] + for _ in range(40): + center = rng.uniform(-8, 8, 3) + objects.append(-openmc.Sphere(*center, r=rng.uniform(0.3, 1.0))) + for _ in range(20): + x0, y0 = rng.uniform(-8, 8, 2) + z0 = rng.uniform(-8, 6) + objects.append(-openmc.ZCylinder(x0, y0, r=rng.uniform(0.2, 0.8)) & + +openmc.ZPlane(z0) & -openmc.ZPlane(z0 + 2.0)) + + # Objects may overlap, so each cell is the part of an object outside of + # the preceding objects + cells = [] + for i, obj in enumerate(objects): + region = openmc.Intersection([obj]) + for other in objects[:i]: + region &= ~other + cells.append(openmc.Cell(fill=mat, region=region)) + background = openmc.Intersection([-box]) + for obj in objects: + background &= ~obj + cells.append(openmc.Cell(fill=mat, region=background)) + + model = openmc.Model() + model.geometry = openmc.Geometry(cells) + model.materials = openmc.Materials([mat]) + model.settings.run_mode = 'fixed source' + model.settings.particles = 1000 + model.settings.batches = 2 + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box((-9, -9, -9), (9, 9, 9))) + return model + + +def test_find_cell_many_cells(run_in_tmpdir, many_cells_model): + model = many_cells_model + model.export_to_model_xml() + cells = model.geometry.root_universe.cells + points = np.random.default_rng(2).uniform(-9.9, 9.9, size=(2000, 3)) + openmc.lib.init() + try: + for p in points: + cell, _ = openmc.lib.find_cell(p) + expected = [c.id for c in cells.values() if tuple(p) in c.region] + assert [cell.id] == expected + finally: + openmc.lib.finalize() + + +def test_transport_many_cells(run_in_tmpdir, many_cells_model): + model = many_cells_model + model.settings.max_lost_particles = 1 + model.run() + + +def random_region(rng, surfaces, depth): + """Random region expression with unions, intersections and complements""" + if depth == 0 or rng.random() < 0.25: + s = surfaces[rng.integers(len(surfaces))] + return -s if rng.random() < 0.5 else +s + n = rng.integers(2, 4) + terms = [random_region(rng, surfaces, depth - 1) for _ in range(n)] + region = openmc.Union(terms) if rng.random() < 0.5 else openmc.Intersection(terms) + return ~region if rng.random() < 0.3 else region + + +@pytest.mark.parametrize('seed', range(5)) +def test_find_cell_random_regions(run_in_tmpdir, seed): + """Cells found for random points agree with Python's evaluation of random + region expressions, and the points lie within the bounding boxes of the + cells, which the search relies on.""" + rng = np.random.default_rng(seed) + surfaces = [ + openmc.XPlane(-1.0), openmc.XPlane(2.0), openmc.YPlane(0.5), + openmc.ZPlane(1.5), openmc.Sphere(x0=1.0, r=3.0), + openmc.ZCylinder(y0=-1.0, r=2.0), openmc.XCylinder(z0=1.0, r=2.5), + openmc.Plane(a=1.0, b=1.0, c=0.0, d=0.5), + ] + outer = openmc.Sphere(r=8.0, boundary_type='vacuum') + + # Split space into cells: each cell is a random region outside of the + # preceding ones + regions = [random_region(rng, surfaces, 3) for _ in range(6)] + cells = [] + remaining = openmc.Intersection([-outer]) + for region in regions: + cells.append(openmc.Cell(region=remaining & region)) + remaining = remaining & ~region + cells.append(openmc.Cell(region=remaining)) + + model = openmc.Model() + model.geometry = openmc.Geometry(cells) + model.settings.run_mode = 'fixed source' + model.settings.particles = 1 + model.settings.batches = 1 + model.export_to_model_xml() + + points = rng.uniform(-5.6, 5.6, size=(500, 3)) + openmc.lib.init() + try: + for p in points: + expected = [c.id for c in cells if tuple(p) in c.region] + if len(expected) != 1: + continue + cell, _ = openmc.lib.find_cell(p) + assert cell.id == expected[0] + lower_left, upper_right = cell.bounding_box + assert np.all(lower_left <= p) and np.all(p <= upper_right) + finally: + openmc.lib.finalize()