Merge pull request #1140 from smharper/neighbor_lists_pr

Switch from surface-based to cell-based neighbor lists

include/openmc/cell.h

@@ -11,6 +11,7 @@
 #include "pugixml.hpp"
 
 #include "openmc/constants.h"
+#include "openmc/neighbor_list.h"
 #include "openmc/position.h"
 
 #ifdef DAGMC
@@ -103,6 +104,9 @@ public:
   std::vector<std::int32_t> rpn_;
   bool simple_;  //!< Does the region contain only intersections?
 
+  //! \brief Neighboring cells in the same universe.
+  NeighborList neighbors_;
+
   Position translation_ {0, 0, 0}; //!< Translation vector for filled universe
 
   //! \brief Rotational tranfsormation of the filled universe.
@@ -151,10 +155,11 @@ public:
   virtual ~Cell() {}
 };
 
+//==============================================================================
+
 class CSGCell : public Cell
 {
 public:
-
   CSGCell();
 
   explicit CSGCell(pugi::xml_node cell_node);
@@ -167,13 +172,13 @@ public:
 
   void to_hdf5(hid_t group_id) const;
 
-
-
 protected:
   bool contains_simple(Position r, Direction u, int32_t on_surface) const;
   bool contains_complex(Position r, Direction u, int32_t on_surface) const;
 };
 
+//==============================================================================
+
 #ifdef DAGMC
 class DAGCell : public Cell
 {
@@ -181,11 +186,12 @@ public:
   moab::DagMC* dagmc_ptr_;
   DAGCell();
 
-  std::pair<double, int32_t> distance(Position r, Direction u, int32_t on_surface) const;
   bool contains(Position r, Direction u, int32_t on_surface) const;
 
-  void to_hdf5(hid_t group_id) const;
+  std::pair<double, int32_t>
+  distance(Position r, Direction u, int32_t on_surface) const;
 
+  void to_hdf5(hid_t group_id) const;
 };
 #endif
 

include/openmc/geometry.h

@@ -33,16 +33,15 @@ check_cell_overlap(Particle* p);
 //!
 //! \param p A particle to be located.  This function will populate the
 //!   geometry-dependent data fields of the particle.
-//! \param search_surf A surface that the particle is expected to be on.  This
-//!   value should be the signed, 1-based index of a surface.  If positive, the
-//!   cells on the positive half-space of the surface will be searched.  If
-//!   negative, the negative half-space will be searched.
+//! \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.
 //==============================================================================
 
 extern "C" bool
-find_cell(Particle* p, int search_surf);
+find_cell(Particle* p, bool use_neighbor_lists);
 
 //==============================================================================
 //! Move a particle into a new lattice tile.

include/openmc/geometry_aux.h

@@ -32,12 +32,6 @@ extern "C" void assign_temperatures();
 
 extern "C" int32_t find_root_universe();
 
-//!=============================================================================
-//! Build a list of neighboring cells to each surface to speed up tracking.
-//!=============================================================================
-
-extern "C" void neighbor_lists();
-
 //==============================================================================
 //! Populate all data structures needed for distribcells.
 //==============================================================================

include/openmc/neighbor_list.h

@@ -0,0 +1,68 @@
+#ifndef OPENMC_NEIGHBOR_LIST_H
+#define OPENMC_NEIGHBOR_LIST_H
+
+#include <algorithm>
+#include <cstdint>
+#include <forward_list>
+#include <mutex>
+
+#include "openmc/openmp_interface.h"
+
+namespace openmc {
+
+//==============================================================================
+//! A threadsafe, dynamic container for listing neighboring cells.
+//
+//! This container is a reduced interface for a linked list with an added OpenMP
+//! lock for write operations.  It allows for threadsafe dynamic growth; any
+//! number of threads can safely read data without locks or reference counting.
+//==============================================================================
+
+class NeighborList
+{
+public:
+  using value_type = int32_t;
+  using const_iterator = std::forward_list<value_type>::const_iterator;
+
+  // Attempt to add an element.
+  //
+  // If the relevant OpenMP lock is currently owned by another thread, this
+  // function will return without actually modifying the data.  It has been
+  // found that returning the transport calculation and possibly re-adding the
+  // element later is slightly faster than waiting on the lock to be released.
+  void push_back(int new_elem)
+  {
+    // Try to acquire the lock.
+    std::unique_lock<OpenMPMutex> lock(mutex_, std::try_to_lock);
+    if (lock) {
+      // It is possible another thread already added this element to the list
+      // while this thread was searching for a cell so make sure the given
+      // element isn't a duplicate before adding it.
+      if (std::find(list_.cbegin(), list_.cend(), new_elem) == list_.cend()) {
+        // Find the end of the list and add the the new element there.
+        if (!list_.empty()) {
+          auto it1 = list_.cbegin();
+          auto it2 = ++list_.cbegin();
+          while (it2 != list_.cend()) it1 = it2++;
+          list_.insert_after(it1, new_elem);
+        } else {
+          list_.push_front(new_elem);
+        }
+      }
+    }
+  }
+
+  const_iterator cbegin() const
+  {return list_.cbegin();}
+
+  const_iterator cend() const
+  {return list_.cend();}
+
+
+private:
+  std::forward_list<value_type> list_;
+  OpenMPMutex mutex_;
+};
+
+} // namespace openmc
+#endif // OPENMC_NEIGHBOR_LIST_H

include/openmc/openmp_interface.h

@@ -0,0 +1,77 @@
+#ifndef OPENMC_OPENMP_INTERFACE_H
+#define OPENMC_OPENMP_INTERFACE_H
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+namespace openmc {
+
+//==============================================================================
+//! An object used to prevent concurrent access to a piece of data.
+//
+//! This type meets the C++ "Lockable" requirements.
+//==============================================================================
+
+class OpenMPMutex
+{
+public:
+  OpenMPMutex()
+  {
+    #ifdef _OPEMP
+      omp_init_lock(&mutex_);
+    #endif
+  }
+
+  ~OpenMPMutex()
+  {
+    #ifdef _OPEMP
+      omp_destroy_lock(&mutex_);
+    #endif
+  }
+
+  // Mutexes cannot be copied.  We need to explicitly delete the copy
+  // constructor and copy assignment operator to ensure the compiler doesn't
+  // "help" us by implicitly trying to copy the underlying mutexes.
+  OpenMPMutex(const OpenMPMutex&) = delete;
+  OpenMPMutex& operator= (const OpenMPMutex&) = delete;
+
+  //! Lock the mutex.
+  //
+  //! This function blocks execution until the lock succeeds.
+  void lock()
+  {
+    #ifdef _OPEMP
+      omp_set_lock(&mutex_);
+    #endif
+  }
+
+  //! Try to lock the mutex and indicate success.
+  //
+  //! This function does not block.  It returns immediately and gives false if
+  //! the lock is unavailable.
+  bool try_lock() noexcept
+  {
+    #ifdef _OPEMP
+      return omp_test_lock(&mutex_);
+    #else
+      return true;
+    #endif
+  }
+
+  //! Unlock the mutex.
+  void unlock() noexcept
+  {
+    #ifdef _OPEMP
+      omp_unset_lock(&mutex_);
+    #endif
+  }
+
+private:
+  #ifdef _OPEMP
+    omp_lock_t mutex_;
+  #endif
+};
+
+} // namespace openmc
+#endif // OPENMC_OPENMP_INTERFACE_H

include/openmc/surface.h

@@ -66,9 +66,6 @@ public:
   int bc_;                    //!< Boundary condition
   std::string name_;          //!< User-defined name
 
-  std::vector<int> neighbor_pos_; //!< List of cells on positive side
-  std::vector<int> neighbor_neg_; //!< List of cells on negative side
-
   explicit Surface(pugi::xml_node surf_node);
   Surface();
 

src/geometry.F90

@@ -28,11 +28,11 @@ module geometry
       type(Particle), intent(in) :: p
     end subroutine check_cell_overlap
 
-    function find_cell_c(p, search_surf) &
+    function find_cell_c(p, use_neighbor_lists) &
          bind(C, name="find_cell") result(found)
       import Particle, C_INT, C_BOOL
       type(Particle),  intent(inout)        :: p
-      integer(C_INT),  intent(in), value    :: search_surf
+      logical(C_BOOL), intent(in), value    :: use_neighbor_lists
       logical(C_BOOL)                       :: found
     end function find_cell_c
 
@@ -53,9 +53,6 @@ module geometry
       integer(C_INT), intent(out)   :: next_level
     end subroutine distance_to_boundary
 
-    subroutine neighbor_lists() bind(C)
-    end subroutine neighbor_lists
-
 #ifdef DAGMC
 
     function next_cell_c(current_cell, surface_crossed) &
@@ -89,15 +86,15 @@ contains
 ! as it's within the geometry
 !===============================================================================
 
-  subroutine find_cell(p, found, search_surf)
+  subroutine find_cell(p, found, use_neighbor_lists)
     type(Particle),    intent(inout) :: p
     logical,           intent(inout) :: found
-    integer, optional, intent(in)    :: search_surf
+    logical, optional, intent(in)    :: use_neighbor_lists
 
-    if (present(search_surf)) then
-      found = find_cell_c(p, search_surf)
+    if (present(use_neighbor_lists)) then
+      found = find_cell_c(p, logical(use_neighbor_lists, kind=C_BOOL))
     else
-      found = find_cell_c(p, 0)
+      found = find_cell_c(p, .false._C_BOOL)
     end if
 
   end subroutine find_cell

src/geometry.cpp

@@ -32,7 +32,8 @@ std::vector<int64_t> overlap_check_count;
 //==============================================================================
 
 extern "C" bool
-check_cell_overlap(Particle* p) {
+check_cell_overlap(Particle* p)
+{
   int n_coord = p->n_coord;
 
   // Loop through each coordinate level
@@ -61,55 +62,59 @@ check_cell_overlap(Particle* p) {
 
 //==============================================================================
 
-extern "C" bool
-find_cell(Particle* p, int search_surf) {
-  for (int i = p->n_coord; i < MAX_COORD; i++) {
-    p->coord[i].reset();
-  }
-
-  // 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[p->n_coord-1].universe = model::root_universe;
-    i_universe = model::root_universe;
-  }
-
-  // If a surface was indicated, only search cells from the neighbor list of
-  // that surface.  The surface index is signed, and the sign signifies whether
-  // the positive or negative side of the surface should be searched.
-  const std::vector<int>* search_cells;
-  if (search_surf > 0) {
-    search_cells = &model::surfaces[search_surf-1]->neighbor_pos_;
-  } else if (search_surf < 0) {
-    search_cells = &model::surfaces[-search_surf-1]->neighbor_neg_;
-  } else {
-    // No surface was indicated, search all cells in the universe.
-    search_cells = &model::universes[i_universe]->cells_;
-  }
-
-  // Find which cell of this universe the particle is in.
+bool
+find_cell_inner(Particle* p, const NeighborList* neighbor_list)
+{
+  // Find which cell of this universe the particle is in.  Use the neighbor list
+  // to shorten the search if one was provided.
   bool found = false;
   int32_t i_cell;
-  for (int i = 0; i < search_cells->size(); i++) {
-    i_cell = (*search_cells)[i];
+  if (neighbor_list) {
+    for (auto it = neighbor_list->cbegin(); it != neighbor_list->cend(); ++it) {
+      i_cell = *it;
 
-    // Make sure the search cell is in the same universe.
-    if (model::cells[i_cell]->universe_ != i_universe) continue;
+      // Make sure the search cell is in the same universe.
+      int i_universe = p->coord[p->n_coord-1].universe;
+      if (model::cells[i_cell]->universe_ != i_universe) continue;
 
-    Position r {p->coord[p->n_coord-1].xyz};
-    Direction u {p->coord[p->n_coord-1].uvw};
-    int32_t surf = p->surface;
-    if (model::cells[i_cell]->contains(r, u, surf)) {
-      p->coord[p->n_coord-1].cell = i_cell;
-
-      if (settings::verbosity >= 10 || simulation::trace) {
-        std::stringstream msg;
-        msg << "    Entering cell " << model::cells[i_cell]->id_;
-        write_message(msg, 1);
+      // Check if this cell contains the particle.
+      Position r {p->coord[p->n_coord-1].xyz};
+      Direction u {p->coord[p->n_coord-1].uvw};
+      auto surf = p->surface;
+      if (model::cells[i_cell]->contains(r, u, surf)) {
+        p->coord[p->n_coord-1].cell = i_cell;
+        found = true;
+        break;
       }
-      found = true;
-      break;
     }
+
+  } else {
+    int i_universe = p->coord[p->n_coord-1].universe;
+    const auto& cells {model::universes[i_universe]->cells_};
+    for (auto it = cells.cbegin(); it != cells.cend(); it++) {
+      i_cell = *it;
+
+      // Make sure the search cell is in the same universe.
+      int i_universe = p->coord[p->n_coord-1].universe;
+      if (model::cells[i_cell]->universe_ != i_universe) continue;
+
+      // Check if this cell contains the particle.
+      Position r {p->coord[p->n_coord-1].xyz};
+      Direction u {p->coord[p->n_coord-1].uvw};
+      auto surf = p->surface;
+      if (model::cells[i_cell]->contains(r, u, surf)) {
+        p->coord[p->n_coord-1].cell = i_cell;
+        found = true;
+        break;
+      }
+    }
+  }
+
+  // Announce the cell that the particle is entering.
+  if (found && (settings::verbosity >= 10 || simulation::trace)) {
+    std::stringstream msg;
+    msg << "    Entering cell " << model::cells[i_cell]->id_;
+    write_message(msg, 1);
   }
 
   if (found) {
@@ -205,7 +210,7 @@ find_cell(Particle* p, int search_surf) {
 
       // Update the coordinate level and recurse.
       ++p->n_coord;
-      return find_cell(p, 0);
+      return find_cell_inner(p, nullptr);
 
     } else if (c.type_ == FILL_LATTICE) {
       //========================================================================
@@ -252,7 +257,7 @@ find_cell(Particle* p, int search_surf) {
 
       // Update the coordinate level and recurse.
       ++p->n_coord;
-      return find_cell(p, 0);
+      return find_cell_inner(p, nullptr);
     }
   }
 
@@ -261,6 +266,48 @@ find_cell(Particle* p, int search_surf) {
 
 //==============================================================================
 
+extern "C" bool
+find_cell(Particle* p, bool use_neighbor_lists)
+{
+  // 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 < MAX_COORD; 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);
+  }
+}
+
+//==============================================================================
+
 extern "C" void
 cross_lattice(Particle* p, int lattice_translation[3])
 {

src/geometry_aux.cpp

@@ -155,33 +155,6 @@ find_root_universe()
 
 //==============================================================================
 
-void
-neighbor_lists()
-{
-  write_message("Building neighboring cells lists for each surface...", 6);
-
-  for (int i = 0; i < model::cells.size(); i++) {
-    for (auto token : model::cells[i]->region_) {
-      // Skip operator tokens.
-      if (std::abs(token) >= OP_UNION) continue;
-
-      // This token is a surface index.  Add the cell to the surface's list.
-      if (token > 0) {
-        model::surfaces[std::abs(token)-1]->neighbor_pos_.push_back(i);
-      } else {
-        model::surfaces[std::abs(token)-1]->neighbor_neg_.push_back(i);
-      }
-    }
-  }
-
-  for (Surface* surf : model::surfaces) {
-    surf->neighbor_pos_.shrink_to_fit();
-    surf->neighbor_neg_.shrink_to_fit();
-  }
-}
-
-//==============================================================================
-
 void
 prepare_distribcell()
 {

src/input_xml.F90

@@ -7,7 +7,6 @@ module input_xml
   use dict_header,      only: DictIntInt, DictCharInt, DictEntryCI
   use endf,             only: reaction_name
   use error,            only: fatal_error, warning, write_message, openmc_err_msg
-  use geometry,         only: neighbor_lists
   use geometry_header
 #ifdef DAGMC
   use dagmc_header
@@ -129,7 +128,7 @@ contains
     call read_materials_xml()
     call read_geometry_xml()
 
-    ! Set up neighbor lists, convert user IDs -> indices, assign temperatures
+    ! Convert user IDs -> indices, assign temperatures
     call finalize_geometry(nuc_temps, sab_temps)
 
     if (run_mode /= MODE_PLOTTING) then
@@ -177,12 +176,6 @@ contains
     call adjust_indices()
     call count_cell_instances(root_universe)
 
-    ! After reading input and basic geometry setup is complete, build lists of
-    ! neighboring cells for efficient tracking
-    if (.not. dagmc) then
-      call neighbor_lists()
-    end if
-
     ! Assign temperatures to cells that don't have temperatures already assigned
     call assign_temperatures()
 

src/tracking.F90

@@ -194,7 +194,6 @@ contains
         end do
         p % last_n_coord = p % n_coord
 
-        p % coord(p % n_coord) % cell = C_NONE
         if (any(lattice_translation /= 0)) then
           ! Particle crosses lattice boundary
           p % surface = ERROR_INT
@@ -403,7 +402,7 @@ contains
       ! the lower universes.
 
       p % n_coord = 1
-      call find_cell(p, found)
+      call find_cell(p, found, .true.)
       if (.not. found) then
         call particle_mark_as_lost(p, "Couldn't find particle after reflecting&
              & from surface " // trim(to_str(surf % id())) // ".")
@@ -458,7 +457,7 @@ contains
 
       ! Figure out what cell particle is in now
       p % n_coord = 1
-      call find_cell(p, found)
+      call find_cell(p, found, .true.)
       if (.not. found) then
         call particle_mark_as_lost(p, "Couldn't find particle after hitting &
              &periodic boundary on surface " // trim(to_str(surf % id())) &
@@ -495,7 +494,7 @@ contains
     end if
 #endif
 
-    call find_cell(p, found, p % surface)
+    call find_cell(p, found, .true.)
     if (found) return
 
     ! ==========================================================================