Stop searching complex cells for a boundary beyond the collision site

include/openmc/cell.h

@@ -81,8 +81,10 @@ class Region {
   bool contains(Position r, Direction u, int32_t on_surface) const;
 
   //! Find the oncoming boundary of this cell.
-  std::pair<double, int32_t> distance(
-    Position r, Direction u, int32_t on_surface) const;
+  //! \param max_distance Distance beyond which the boundary is not needed. If
+  //!   the boundary is not nearer, the distance returned may be INFTY.
+  std::pair<double, int32_t> distance(Position r, Direction u,
+    int32_t on_surface, double max_distance = INFTY) const;
 
   //! Get the BoundingBox for this cell.
   BoundingBox bounding_box(int32_t cell_id) const;
@@ -126,7 +128,7 @@ class Region {
 
   //! Find the oncoming boundary of this cell for a complex cell.
   std::pair<double, int32_t> distance_complex(
-    Position r, Direction u, int32_t on_surface) const;
+    Position r, Direction u, int32_t on_surface, double max_distance) const;
 
   //! BoundingBox if the particle is in a simple cell.
   BoundingBox bounding_box_simple() const;
@@ -216,8 +218,11 @@ class Cell {
   virtual bool contains(Position r, Direction u, int32_t on_surface) const = 0;
 
   //! Find the oncoming boundary of this cell.
-  virtual std::pair<double, int32_t> distance(
-    Position r, Direction u, int32_t on_surface, GeometryState* p) const = 0;
+  //! \param max_distance Distance beyond which the boundary is not needed. If
+  //!   the boundary is not nearer, the distance returned may be INFTY.
+  virtual std::pair<double, int32_t> distance(Position r, Direction u,
+    int32_t on_surface, GeometryState* p,
+    double max_distance = INFTY) const = 0;
 
   //! Write all information needed to reconstruct the cell to an HDF5 group.
   //! \param group_id An HDF5 group id.
@@ -437,9 +442,10 @@ class CSGCell : public Cell {
   int n_surfaces() const override { return region_.n_surfaces(); }
 
   std::pair<double, int32_t> distance(Position r, Direction u,
-    int32_t on_surface, GeometryState* p) const override
+    int32_t on_surface, GeometryState* p,
+    double max_distance = INFTY) const override
   {
-    return region_.distance(r, u, on_surface);
+    return region_.distance(r, u, on_surface, max_distance);
   }
 
   bool contains(Position r, Direction u, int32_t on_surface) const override

include/openmc/dagmc.h

@@ -69,7 +69,8 @@ class DAGCell : public Cell {
   bool contains(Position r, Direction u, int32_t on_surface) const override;
 
   std::pair<double, int32_t> distance(Position r, Direction u,
-    int32_t on_surface, GeometryState* p) const override;
+    int32_t on_surface, GeometryState* p,
+    double max_distance = INFTY) const override;
 
   BoundingBox bounding_box() const override;
 

include/openmc/geometry.h

@@ -120,7 +120,10 @@ void cross_lattice(
 //! Find the next boundary a particle will intersect.
 //==============================================================================
 
-BoundaryInfo distance_to_boundary(GeometryState& p);
+//! \param max_distance Distance beyond which a boundary is not needed. The
+//!   distance to a farther boundary may be returned as INFTY.
+BoundaryInfo distance_to_boundary(
+  GeometryState& p, double max_distance = INFTY);
 
 } // namespace openmc
 

src/cell.cpp

@@ -948,12 +948,12 @@ std::string Region::str() const
 //==============================================================================
 
 std::pair<double, int32_t> Region::distance(
-  Position r, Direction u, int32_t on_surface) const
+  Position r, Direction u, int32_t on_surface, double max_distance) const
 {
   if (simple_) {
     return distance_to_nearest_surface(r, u, on_surface, false);
   } else {
-    return distance_complex(r, u, on_surface);
+    return distance_complex(r, u, on_surface, max_distance);
   }
 }
 
@@ -997,7 +997,7 @@ std::pair<double, int32_t> Region::distance_to_nearest_surface(Position r,
 //==============================================================================
 
 std::pair<double, int32_t> Region::distance_complex(
-  Position r, Direction u, int32_t on_surface) const
+  Position r, Direction u, int32_t on_surface, double max_distance) const
 {
   const bool in_region = contains_complex(r, u, on_surface);
   double total_distance {0.0};
@@ -1015,6 +1015,13 @@ std::pair<double, int32_t> Region::distance_complex(
     // the wrong side of a curved surface.
     r += distance * u;
     total_distance += distance;
+
+    // Stop searching once the boundary is known to be beyond the distance of
+    // interest
+    if (total_distance >= max_distance) {
+      return {INFTY, std::numeric_limits<int32_t>::max()};
+    }
+
     i_surf = std::abs(i_surf);
     const auto& surf {*model::surfaces[i_surf - 1]};
     if (u.dot(surf.normal(r)) <= 0.0) {

src/dagmc.cpp

@@ -821,8 +821,8 @@ void DAGUniverse::override_assign_material(std::unique_ptr<DAGCell>& c,
 DAGCell::DAGCell(std::shared_ptr<moab::DagMC> dag_ptr, int32_t dag_idx)
   : Cell {}, dagmc_ptr_(dag_ptr), dag_index_(dag_idx) {};
 
-std::pair<double, int32_t> DAGCell::distance(
-  Position r, Direction u, int32_t on_surface, GeometryState* p) const
+std::pair<double, int32_t> DAGCell::distance(Position r, Direction u,
+  int32_t on_surface, GeometryState* p, double /*max_distance*/) const
 {
   // if we've changed direction or we're not on a surface,
   // reset the history and update last direction

src/geometry.cpp

@@ -426,7 +426,7 @@ void cross_lattice(GeometryState& p, const BoundaryInfo& boundary, bool verbose)
 
 //==============================================================================
 
-BoundaryInfo distance_to_boundary(GeometryState& p)
+BoundaryInfo distance_to_boundary(GeometryState& p, double max_distance)
 {
   BoundaryInfo info;
   double d_lat = INFINITY;
@@ -442,7 +442,7 @@ BoundaryInfo distance_to_boundary(GeometryState& p)
     Cell& c {*model::cells[coord.cell()]};
 
     // Find the oncoming surface in this cell and the distance to it.
-    auto surface_distance = c.distance(r, u, p.surface(), &p);
+    auto surface_distance = c.distance(r, u, p.surface(), &p, max_distance);
     d_surf = surface_distance.first;
     level_surf_cross = surface_distance.second;
 

src/particle.cpp

@@ -276,9 +276,6 @@ void Particle::event_calculate_xs()
 
 void Particle::event_advance()
 {
-  // Find the distance to the nearest boundary
-  boundary() = distance_to_boundary(*this);
-
   // Sample a distance to collision
   if (type() == ParticleType::electron() ||
       type() == ParticleType::positron()) {
@@ -294,6 +291,18 @@ void Particle::event_advance()
   double distance_cutoff =
     (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;
 
+  // Find the distance to the nearest boundary. A boundary beyond the
+  // collision site is not needed, which lets the search through complex cells
+  // stop early, and no boundary is needed at all for a particle that collides
+  // where it is, such as an electron or positron whose energy is deposited
+  // locally. The search uses no random numbers, so results do not depend on
+  // the order in which the two distances are found.
+  if (collision_distance() > 0.0) {
+    boundary() = distance_to_boundary(*this, collision_distance());
+  } else {
+    boundary().reset();
+  }
+
   // Select smaller of the three distances
   double distance =
     std::min({boundary().distance(), collision_distance(), distance_cutoff});

tests/cpp_unit_tests/test_region.cpp

@@ -183,6 +183,23 @@ TEST_CASE("Find boundary after virtual surface crossings")
     REQUIRE(surface == 2);
   }
 
+  SECTION("Boundary beyond a maximum distance")
+  {
+    // A boundary nearer than the maximum distance is found as before
+    auto [distance, surface] =
+      region.distance({1.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, 0, 5.5);
+    REQUIRE(distance == Catch::Approx(5.0));
+    REQUIRE(surface == 2);
+
+    // The search stops once it passes the maximum distance, whether before or
+    // after the virtual crossings
+    for (double max_distance : {2.0, 4.5}) {
+      auto [far_distance, far_surface] =
+        region.distance({1.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, 0, max_distance);
+      REQUIRE(far_distance == openmc::INFTY);
+    }
+  }
+
   SECTION("Starting outside the region")
   {
     // Along -x from x=7, entering the sphere at x=6 is the first boundary.