Limit DAGMC ray fire distance to collision point for BVH pruning

include/openmc/cell.h

@@ -81,8 +81,8 @@ 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;
+  std::pair<double, int32_t> distance(Position r, Direction u,
+    int32_t on_surface, double max_distance = INFINITY) const;
 
   //! Get the BoundingBox for this cell.
   BoundingBox bounding_box(int32_t cell_id) const;
@@ -205,8 +205,9 @@ 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;
+  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.
@@ -421,7 +422,8 @@ class CSGCell : public Cell {
   vector<int32_t> surfaces() const override { return region_.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);
   }

include/openmc/dagmc.h

@@ -74,7 +74,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

@@ -75,7 +75,8 @@ void cross_lattice(
 //! Find the next boundary a particle will intersect.
 //==============================================================================
 
-BoundaryInfo distance_to_boundary(GeometryState& p);
+BoundaryInfo distance_to_boundary(
+  GeometryState& p, double max_distance = INFINITY);
 
 } // namespace openmc
 

src/cell.cpp

@@ -944,7 +944,7 @@ 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
 {
   double min_dist {INFTY};
   int32_t i_surf {std::numeric_limits<int32_t>::max()};

src/dagmc.cpp

@@ -806,8 +806,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
@@ -826,8 +826,10 @@ std::pair<double, int32_t> DAGCell::distance(
     fatal_error("DAGMC call made for particle in a non-DAGMC universe");
 
   // initialize to lost particle conditions
+  // surface idx of 1 and distance of infinity occur when no previous collisions
+  // recorded
   int surf_idx = -1;
-  double dist = INFINITY;
+  double dist = max_distance;
 
   moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
   moab::EntityHandle hit_surf;
@@ -840,7 +842,7 @@ std::pair<double, int32_t> DAGCell::distance(
   if (hit_surf != 0) {
     surf_idx =
       dag_univ->surf_idx_offset_ + dagmc_ptr_->index_by_handle(hit_surf);
-  } else if (!dagmc_ptr_->is_implicit_complement(vol) ||
+  } else if (dist == INFTY && !dagmc_ptr_->is_implicit_complement(vol) ||
              is_root_universe(dag_univ->id_)) {
     // surface boundary conditions are ignored for projection plotting, meaning
     // that the particle may move through the graveyard (bounding) volume and

src/geometry.cpp

@@ -358,7 +358,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;

src/particle.cpp

@@ -267,9 +267,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()) {
@@ -285,9 +282,13 @@ void Particle::event_advance()
   double distance_cutoff =
     (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;
 
+  // Find distance to nearest boundary;
+  // cap bvh traversal distance to collision distance
+  double max_distance = std::min(collision_distance(), distance_cutoff);
+  boundary() = distance_to_boundary(*this, max_distance);
+
   // Select smaller of the three distances
-  double distance =
-    std::min({boundary().distance(), collision_distance(), distance_cutoff});
+  double distance = std::min(boundary().distance(), collision_distance());
 
   // Advance particle in space and time
   this->move_distance(distance);