Skip cells whose bounding boxes do not contain the point in cell searches

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 <https://doi.org/10.1080/00295639.2020.1719765>`_.
+`Harper et al <https://doi.org/10.1080/00295639.2020.1719765>`_. 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:
 

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.

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;

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();

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;
     }

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()