separate find_cell into two functions

include/openmc/geometry.h

@@ -43,14 +43,11 @@ bool check_cell_overlap(Particle& p, bool error=true);
 //!
 //! \param p A particle to be located.  This function will populate the
 //!   geometry-dependent data fields of the particle.
-//! \param use_neighbor_lists If true, neighbor lists will be used to accelerate
-//!   the geometry search, but this only works if the cell attribute of the
-//!   particle's lowest coordinate level is valid and meaningful.
 //! \return True if the particle's location could be found and ascribed to a
 //!   valid geometry coordinate stack.
 //==============================================================================
-
-bool find_cell(Particle& p, bool use_neighbor_lists);
+bool brute_force_find_cell(Particle& p);
+bool neighbor_list_find_cell(Particle& p); // Only usable on surface crossings
 
 //==============================================================================
 //! Move a particle into a new lattice tile.

include/openmc/plot.h

@@ -182,7 +182,7 @@ T PlotBase::get_map() const {
         p.r()[in_i] = xyz[in_i] + in_pixel * x;
         p.n_coord_ = 1;
         // local variables
-        bool found_cell = find_cell(p, 0);
+        bool found_cell = brute_force_find_cell(p);
         j = p.n_coord_ - 1;
         if (level >= 0) { j = level; }
         if (found_cell) {

src/geometry.cpp

@@ -246,44 +246,37 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list)
 
 //==============================================================================
 
-bool
-find_cell(Particle& p, bool use_neighbor_lists)
+bool neighbor_list_find_cell(Particle& p)
+{
+  // Get the cell this particle was in previously.
+  auto coord_lvl = p.n_coord_ - 1;
+  auto i_cell = p.coord_[coord_lvl].cell;
+  Cell& c {*model::cells[i_cell]};
+
+  // Search for the particle in that cell's neighbor list.  Return if we
+  // found the particle.
+  bool found = find_cell_inner(p, &c.neighbors_);
+  if (found)
+    return found;
+
+  // The particle could not be found in the neighbor list.  Try searching all
+  // cells in this universe, and update the neighbor list if we find a new
+  // neighboring cell.
+  found = find_cell_inner(p, nullptr);
+  if (found)
+    c.neighbors_.push_back(p.coord_[coord_lvl].cell);
+  return found;
+}
+
+bool brute_force_find_cell(Particle& p)
 {
-  // Determine universe (if not yet set, use root universe).
   int i_universe = p.coord_[p.n_coord_-1].universe;
   if (i_universe == C_NONE) {
     p.coord_[0].universe = model::root_universe;
     p.n_coord_ = 1;
     i_universe = model::root_universe;
   }
-
-  // Reset all the deeper coordinate levels.
-  for (int i = p.n_coord_; i < p.coord_.size(); i++) {
-    p.coord_[i].reset();
-  }
-
-  if (use_neighbor_lists) {
-    // Get the cell this particle was in previously.
-    auto coord_lvl = p.n_coord_ - 1;
-    auto i_cell = p.coord_[coord_lvl].cell;
-    Cell& c {*model::cells[i_cell]};
-
-    // Search for the particle in that cell's neighbor list.  Return if we
-    // found the particle.
-    bool found = find_cell_inner(p, &c.neighbors_);
-    if (found) return found;
-
-    // The particle could not be found in the neighbor list.  Try searching all
-    // cells in this universe, and update the neighbor list if we find a new
-    // neighboring cell.
-    found = find_cell_inner(p, nullptr);
-    if (found) c.neighbors_.push_back(p.coord_[coord_lvl].cell);
-    return found;
-
-  } else {
-    // Search all cells in this universe for the particle.
-    return find_cell_inner(p, nullptr);
-  }
+  return find_cell_inner(p, nullptr);
 }
 
 //==============================================================================
@@ -319,7 +312,7 @@ cross_lattice(Particle& p, const BoundaryInfo& boundary)
   if (!lat.are_valid_indices(i_xyz)) {
     // The particle is outside the lattice.  Search for it from the base coords.
     p.n_coord_ = 1;
-    bool found = find_cell(p, 0);
+    bool found = brute_force_find_cell(p);
     if (!found && p.alive_) {
       p.mark_as_lost(fmt::format("Could not locate particle {} after "
         "crossing a lattice boundary", p.id_));
@@ -328,13 +321,13 @@ cross_lattice(Particle& p, const BoundaryInfo& boundary)
   } else {
     // Find cell in next lattice element.
     p.coord_[p.n_coord_-1].universe = lat[i_xyz];
-    bool found = find_cell(p, 0);
+    bool found = brute_force_find_cell(p);
 
     if (!found) {
       // A particle crossing the corner of a lattice tile may not be found.  In
       // this case, search for it from the base coords.
       p.n_coord_ = 1;
-      bool found = find_cell(p, 0);
+      bool found = brute_force_find_cell(p);
       if (!found && p.alive_) {
         p.mark_as_lost(fmt::format("Could not locate particle {} after "
           "crossing a lattice boundary", p.id_));
@@ -450,7 +443,7 @@ openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance)
   p.r() = Position{xyz};
   p.u() = {0.0, 0.0, 1.0};
 
-  if (!find_cell(p, false)) {
+  if (!brute_force_find_cell(p)) {
     set_errmsg(fmt::format("Could not find cell at position {}.", p.r()));
     return OPENMC_E_GEOMETRY;
   }

src/particle.cpp

@@ -156,7 +156,7 @@ Particle::event_calculate_xs()
   // initiate a search for the current cell. This generally happens at the
   // beginning of the history and again for any secondary particles
   if (coord_[n_coord_ - 1].cell == C_NONE) {
-    if (!find_cell(*this, false)) {
+    if (!brute_force_find_cell(*this)) {
       this->mark_as_lost("Could not find the cell containing particle "
         + std::to_string(id_));
       return;
@@ -454,7 +454,8 @@ Particle::cross_surface()
   }
 #endif
 
-  if (find_cell(*this, true)) return;
+  if (neighbor_list_find_cell(*this))
+    return;
 
   // ==========================================================================
   // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
@@ -462,7 +463,7 @@ Particle::cross_surface()
   // Remove lower coordinate levels and assignment of surface
   surface_ = 0;
   n_coord_ = 1;
-  bool found = find_cell(*this, false);
+  bool found = brute_force_find_cell(*this);
 
   if (settings::run_mode != RunMode::PLOTTING && (!found)) {
     // If a cell is still not found, there are two possible causes: 1) there is
@@ -476,7 +477,7 @@ Particle::cross_surface()
     // Couldn't find next cell anywhere! This probably means there is an actual
     // undefined region in the geometry.
 
-    if (!find_cell(*this, false)) {
+    if (!brute_force_find_cell(*this)) {
       this->mark_as_lost("After particle " + std::to_string(id_) +
         " crossed surface " + std::to_string(surf->id_) +
         " it could not be located in any cell and it did not leak.");
@@ -553,7 +554,7 @@ Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
   // (unless we're using a dagmc model, which has exactly one universe)
   if (!settings::dagmc) {
     n_coord_ = 1;
-    if (!find_cell(*this, true)) {
+    if (!neighbor_list_find_cell(*this)) {
       this->mark_as_lost("Couldn't find particle after reflecting from surface "
                          + std::to_string(surf.id_) + ".");
       return;
@@ -601,7 +602,7 @@ Particle::cross_periodic_bc(const Surface& surf, Position new_r,
   // Figure out what cell particle is in now
   n_coord_ = 1;
 
-  if (!find_cell(*this, true)) {
+  if (!neighbor_list_find_cell(*this)) {
     this->mark_as_lost("Couldn't find particle after hitting periodic "
       "boundary on surface " + std::to_string(surf.id_) + ". The normal vector "
       "of one periodic surface may need to be reversed.");

src/volume_calc.cpp

@@ -135,7 +135,8 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
         p.u() = {0.5, 0.5, 0.5};
 
         // If this location is not in the geometry at all, move on to next block
-        if (!find_cell(p, false)) continue;
+        if (!brute_force_find_cell(p))
+          continue;
 
         if (domain_type_ == TallyDomain::MATERIAL) {
           if (p.material_ != MATERIAL_VOID) {