Merge 73b9727f1e into db673b9acb

include/openmc/cell.h

@@ -81,8 +81,8 @@ public:
   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 @@ public:
   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 @@ public:
   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 @@ public:
   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;
 
@@ -105,13 +106,15 @@ public:
   //! \param[in] auto_geom_ids Whether or not to automatically assign cell and
   //! surface IDs \param[in] auto_mat_ids Whether or not to automatically assign
   //! material IDs
+  //! \param[in] use_collision_distance_cap Whether or not to cap the distance to
+  //! boundaries to the collision distance when finding the next boundary
   explicit DAGUniverse(const std::string& filename, bool auto_geom_ids = false,
-    bool auto_mat_ids = false);
+    bool auto_mat_ids = false, bool use_collision_distance_cap = false);
 
   //! Alternative DAGMC universe constructor for external DAGMC instance
   explicit DAGUniverse(std::shared_ptr<moab::DagMC> external_dagmc_ptr,
     const std::string& filename = "", bool auto_geom_ids = false,
-    bool auto_mat_ids = false);
+    bool auto_mat_ids = false, bool use_collision_distance_cap = false);
 
   //! Initialize the DAGMC accel. data structures, indices, material
   //! assignments, etc.
@@ -190,6 +193,7 @@ public:
   // Accessors
   moab::DagMC* dagmc_ptr() const { return dagmc_instance_.get(); }
   bool has_graveyard() const { return has_graveyard_; }
+  bool use_collision_distance_cap() const { return use_collision_distance_cap_; }
 
 private:
   void set_id();        //!< Deduce the universe id from model::universes
@@ -211,6 +215,9 @@ private:
                              //!< generate new material IDs for the universe
   bool has_graveyard_; //!< Indicates if the DAGMC geometry has a "graveyard"
                        //!< volume
+  bool use_collision_distance_cap_; //!< Indicates whether or not to cap the
+                                  //!< distance to boundaries to the collision
+                                  //!< distance when finding the next boundary
 };
 
 //==============================================================================

include/openmc/geometry.h

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

@@ -80,6 +80,11 @@ DAGUniverse::DAGUniverse(pugi::xml_node node)
     adjust_material_ids_ = get_node_value_bool(node, "auto_mat_ids");
   }
 
+  use_collision_distance_cap_ = false;
+  if (check_for_node(node, "use_collision_distance_cap")) {
+    use_collision_distance_cap_ = get_node_value_bool(node, "use_collision_distance_cap");
+  }
+
   // Get material assignment overrides from nested DAGMC cell elements.
   if (node.child("cell")) {
     for (pugi::xml_node cell_node : node.children("cell")) {
@@ -154,18 +159,18 @@ DAGUniverse::DAGUniverse(pugi::xml_node node)
 }
 
 DAGUniverse::DAGUniverse(
-  const std::string& filename, bool auto_geom_ids, bool auto_mat_ids)
+  const std::string& filename, bool auto_geom_ids, bool auto_mat_ids, bool use_collision_distance_cap)
   : filename_(filename), adjust_geometry_ids_(auto_geom_ids),
-    adjust_material_ids_(auto_mat_ids)
+    adjust_material_ids_(auto_mat_ids), use_collision_distance_cap_(use_collision_distance_cap)
 {
   set_id();
   initialize();
 }
 
 DAGUniverse::DAGUniverse(std::shared_ptr<moab::DagMC> dagmc_ptr,
-  const std::string& filename, bool auto_geom_ids, bool auto_mat_ids)
+  const std::string& filename, bool auto_geom_ids, bool auto_mat_ids, bool use_collision_distance_cap)
   : dagmc_instance_(dagmc_ptr), filename_(filename),
-    adjust_geometry_ids_(auto_geom_ids), adjust_material_ids_(auto_mat_ids)
+    adjust_geometry_ids_(auto_geom_ids), adjust_material_ids_(auto_mat_ids), use_collision_distance_cap_(use_collision_distance_cap)
 {
   MaterialOverrides material_overrides;
   TemperatureOverrides temperature_overrides;
@@ -806,8 +811,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 +831,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 = INFTY;
 
   moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
   moab::EntityHandle hit_surf;
@@ -836,12 +843,12 @@ std::pair<double, int32_t> DAGCell::distance(
   double pnt[3] = {r.x, r.y, r.z};
   double dir[3] = {u.x, u.y, u.z};
   MB_CHK_ERR_CONT(
-    dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &p->history()));
+    dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &p->history(), max_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) ||
-             is_root_universe(dag_univ->id_)) {
+  } else if (max_distance == 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
     // into the implicit complement on the other side where no intersection will

src/geometry.cpp

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

@@ -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,19 @@ 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 = INFTY;
+  const auto& univ = model::universes[lowest_coord().universe()];
+  if (auto* dag_univ = dynamic_cast<DAGUniverse*>(univ.get())) {
+    if (dag_univ->use_collision_distance_cap()) {
+      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);