From e358f6d90a656318b7d298b90d11a07746e43949 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 24 Nov 2018 09:29:24 -0500 Subject: [PATCH 1/8] Switch to cell-based neighbor lists --- include/openmc/cell.h | 3 ++ include/openmc/geometry.h | 9 ++--- src/geometry.F90 | 8 ++-- src/geometry.cpp | 81 ++++++++++++++++++++++++++------------- src/tracking.F90 | 1 - 5 files changed, 66 insertions(+), 36 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index f2c8f0f71..fc8dd3789 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -103,6 +103,9 @@ public: std::vector rpn_; bool simple_; //!< Does the region contain only intersections? + //! \brief Neighboring cells in the same universe. + std::vector neighbors; + Position translation_ {0, 0, 0}; //!< Translation vector for filled universe //! \brief Rotational tranfsormation of the filled universe. diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index 2f0853e0c..a2da2c5f2 100644 --- a/include/openmc/geometry.h +++ b/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. diff --git a/src/geometry.F90 b/src/geometry.F90 index 72caab359..9c8c4df4a 100644 --- a/src/geometry.F90 +++ b/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 @@ -95,9 +95,9 @@ contains integer, optional, intent(in) :: search_surf if (present(search_surf)) then - found = find_cell_c(p, search_surf) + found = find_cell_c(p, logical(.true., kind=C_BOOL)) else - found = find_cell_c(p, 0) + found = find_cell_c(p, logical(.false., kind=C_BOOL)) end if end subroutine find_cell diff --git a/src/geometry.cpp b/src/geometry.cpp index fe3d52c68..a7d694231 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -32,7 +32,8 @@ std::vector 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,29 +62,13 @@ 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* 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. +bool +find_cell_inner(Particle* p, const std::vector* search_cells) +{ + // If a set of cells to search was not specified, search all cells in this + // universe. + if (!search_cells) { + int i_universe = p->coord[p->n_coord-1].universe; search_cells = &model::universes[i_universe]->cells_; } @@ -94,6 +79,7 @@ find_cell(Particle* p, int search_surf) { i_cell = (*search_cells)[i]; // 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}; @@ -205,7 +191,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 +238,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 +247,49 @@ 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. + std::vector* search_cells = &c.neighbors; + bool found = find_cell_inner(p, search_cells); + 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]) { diff --git a/src/tracking.F90 b/src/tracking.F90 index e322f31f2..cb83a368e 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -195,7 +195,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 From 36d85dde4f5854642c74aced1f59f35130ea4191 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 24 Nov 2018 10:29:25 -0500 Subject: [PATCH 2/8] Use neighbor lists for reflective & periodic BCs --- src/geometry.F90 | 8 ++++---- src/tracking.F90 | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 9c8c4df4a..91ed6bbdc 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -89,13 +89,13 @@ 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, logical(.true., kind=C_BOOL)) + if (present(use_neighbor_lists)) then + found = find_cell_c(p, logical(use_neighbor_lists, kind=C_BOOL)) else found = find_cell_c(p, logical(.false., kind=C_BOOL)) end if diff --git a/src/tracking.F90 b/src/tracking.F90 index cb83a368e..283d69b7c 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -403,7 +403,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 +458,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 +495,7 @@ contains end if #endif - call find_cell(p, found, p % surface) + call find_cell(p, found, .true.) if (found) return ! ========================================================================== From a8cb86c6b5aca279f3631ca03aec70050b36a932 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 14 Dec 2018 12:06:02 -0500 Subject: [PATCH 3/8] Implement linked-list/vector hybrid NeighborList --- include/openmc/cell.h | 3 +- include/openmc/geometry_aux.h | 9 ++ include/openmc/neighbor_list.h | 165 ++++++++++++++++++++++++++++++ include/openmc/openmp_interface.h | 50 +++++++++ src/geometry.cpp | 76 ++++++++------ src/geometry_aux.cpp | 8 ++ src/simulation.cpp | 4 + 7 files changed, 285 insertions(+), 30 deletions(-) create mode 100644 include/openmc/neighbor_list.h create mode 100644 include/openmc/openmp_interface.h diff --git a/include/openmc/cell.h b/include/openmc/cell.h index fc8dd3789..59b5e00a2 100644 --- a/include/openmc/cell.h +++ b/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 @@ -104,7 +105,7 @@ public: bool simple_; //!< Does the region contain only intersections? //! \brief Neighboring cells in the same universe. - std::vector neighbors; + NeighborList neighbors; Position translation_ {0, 0, 0}; //!< Translation vector for filled universe diff --git a/include/openmc/geometry_aux.h b/include/openmc/geometry_aux.h index da55c3d9b..4e634ebca 100644 --- a/include/openmc/geometry_aux.h +++ b/include/openmc/geometry_aux.h @@ -88,6 +88,15 @@ distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset); extern "C" int maximum_levels(int32_t univ); +//============================================================================== +//! Perform any geometry operations that should be done between generations. +//! +//! Currently, this just includes converting cell neighbor linked lists into +//! vectors, but similar optimizations could be added here in the future. +//============================================================================== + +void geometry_finalize_generation(); + //============================================================================== //! Deallocates global vectors and maps for cells, universes, and lattices. //============================================================================== diff --git a/include/openmc/neighbor_list.h b/include/openmc/neighbor_list.h new file mode 100644 index 000000000..4b937a69f --- /dev/null +++ b/include/openmc/neighbor_list.h @@ -0,0 +1,165 @@ +#ifndef OPENMC_NEIGHBOR_LIST_H +#define OPENMC_NEIGHBOR_LIST_H + +#include +#include +#include +#include + +#include "openmc/openmp_interface.h" + +namespace openmc { + +// Forward declare the neighbor list iterator type. +template +class NeighborListIter; + +//============================================================================== +//! An efficient, threadsafe, dynamic container for listing neighboring cells. +// +//! This container is a minor improvement upon a linked list. The linked list +//! allows for threadsafe dynamic growth; any number of threads can safely +//! read from the list without locks or reference counting. Write access must +//! be protected with a lock, but write events are fairly rare for neighbor +//! lists. This container also provides the option to flush the linked list +//! into a vector at some convenient time (like between generations where the +//! neighbor lists are not being read) for contiguous data. The resulting +//! performance benefit is small but consistent. +//============================================================================== + +class NeighborList +{ +public: + using value_type = int32_t; + using prefix_iter = std::vector::iterator; + using suffix_iter = std::forward_list::iterator; + using iterator = NeighborListIter; + + void push(int new_elem) + { + // Try to acquire the lock. + std::unique_lock lock(mutex_, std::try_to_lock); + if (lock) { + // Make sure this element isn't already in the suffix. Don't check the + // prefix because we don't expect this function to be called if the "new" + // element is already there. + if (std::find(suffix_.cbegin(), suffix_.cend(), new_elem) + == suffix_.cend()) { + suffix_.push_front(new_elem); + } + } + } + + void make_consecutive() + { + while (!suffix_.empty()) { + prefix_.push_back(suffix_.front()); + suffix_.pop_front(); + } + } + + iterator begin(); + + iterator end(); + +private: + std::vector prefix_; + std::forward_list suffix_; + ThreadMutex mutex_; + + friend class NeighborListIter; +}; + +//============================================================================== + +template +class NeighborListIter +{ +public: + NeighborListIter(NeighborList* nl, T_prefix_iter it) + { + base_ = nl; + if (it != base_->prefix_.end()) { + in_prefix_ = true; + prefix_iter_ = it; + } else { + in_prefix_ = false; + suffix_iter_ = base_->suffix_.begin(); + } + } + + NeighborListIter(NeighborList* nl, T_suffix_iter it) + { + in_prefix_ = false; + suffix_iter_ = it; + } + + T_value operator*() + { + if (in_prefix_) { + return *prefix_iter_; + } else { + return *suffix_iter_; + } + } + + bool operator==(const NeighborListIter& other) + { + if (in_prefix_ != other.in_prefix_) return false; + + if (in_prefix_) { + return prefix_iter_ == other.prefix_iter_; + } else { + return suffix_iter_ == other.suffix_iter_; + } + } + + bool operator!=(const NeighborListIter& other) + {return !(*this == other);} + + NeighborListIter& operator++() + { + if (in_prefix_) { + // We are in the prefix so increment the prefix iterator. + ++prefix_iter_; + + // If we've reached the end of the prefix, switch to the suffix iterator. + if (prefix_iter_ == base_->prefix_.end()) { + in_prefix_ = false; + suffix_iter_ = base_->suffix_.begin(); + } + + } else { + // We are in the suffix so increment the suffix iterator. + ++suffix_iter_; + } + return *this; + } + +private: + NeighborList* base_; + + union { + T_prefix_iter prefix_iter_; + T_suffix_iter suffix_iter_; + }; + + bool in_prefix_; +}; + +//============================================================================== + +inline NeighborList::iterator +NeighborList::begin() +{ + return iterator(this, prefix_.begin()); +} + +inline NeighborList::iterator +NeighborList::end() +{ + return iterator(this, suffix_.end()); +} + +} // namespace openmc +#endif // OPENMC_NEIGHBOR_LIST_H diff --git a/include/openmc/openmp_interface.h b/include/openmc/openmp_interface.h new file mode 100644 index 000000000..563500a69 --- /dev/null +++ b/include/openmc/openmp_interface.h @@ -0,0 +1,50 @@ +#ifndef OPENMC_OPENMP_INTERFACE_H +#define OPENMC_OPENMP_INTERFACE_H + +#ifdef _OPENMP +#include +#endif + +namespace openmc { + +class ThreadMutex +{ + +#ifdef _OPENMP +//============================================================================== +// Implementation for when OpenMP is actually defined. + +public: + ThreadMutex() + {omp_init_lock(&mutex_);} + + ~ThreadMutex() + {omp_destroy_lock(&mutex_);} + + void lock() + {omp_set_lock(&mutex_);} + + bool try_lock() + {return omp_test_lock(&mutex_);} + + void unlock() + {omp_unset_lock(&mutex_);} + +private: + omp_lock_t mutex_; + +#else +//============================================================================== +// Empty implementation for when OpenMP is not defined. + +public: + ThreadMutex() {} + ~ThreadMutex() = default; + void lock() {} + bool try_lock() {return true;} + void unlock() {} +#endif +}; + +} // namespace openmc +#endif // OPENMC_OPENMP_INTERFACE_H diff --git a/src/geometry.cpp b/src/geometry.cpp index a7d694231..baf04f705 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -63,39 +63,58 @@ check_cell_overlap(Particle* p) //============================================================================== bool -find_cell_inner(Particle* p, const std::vector* search_cells) +find_cell_inner(Particle* p, NeighborList* neighbor_list) { - // If a set of cells to search was not specified, search all cells in this - // universe. - if (!search_cells) { - int i_universe = p->coord[p->n_coord-1].universe; - search_cells = &model::universes[i_universe]->cells_; - } - - // Find which cell of this universe the particle is in. + // 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->begin(); it != neighbor_list->end(); ++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; + // 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.begin(); it != cells.end(); 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) { @@ -271,15 +290,14 @@ find_cell(Particle* p, bool use_neighbor_lists) // Search for the particle in that cell's neighbor list. Return if we // found the particle. - std::vector* search_cells = &c.neighbors; - bool found = find_cell_inner(p, search_cells); + 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); + if (found) c.neighbors.push(p->coord[coord_lvl].cell); return found; } else { diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index c326744eb..d411411eb 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -432,6 +432,14 @@ maximum_levels(int32_t univ) //============================================================================== +void +geometry_finalize_generation() +{ + for (Cell* c : model::cells) c->neighbors.make_consecutive(); +} + +//============================================================================== + void free_memory_geometry_c() { diff --git a/src/simulation.cpp b/src/simulation.cpp index 50bf12521..15e2f7cef 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -4,6 +4,7 @@ #include "openmc/container_util.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" +#include "openmc/geometry_aux.h" #include "openmc/message_passing.h" #include "openmc/output.h" #include "openmc/particle.h" @@ -449,6 +450,9 @@ void finalize_generation() // For fixed-source mode, we need to sample the external source fill_source_bank_fixedsource(); } + + // Clean up neighbor lists + geometry_finalize_generation(); } void initialize_history(Particle* p, int64_t index_source) From 066f698da72f46e975ede789394705311f446593 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 15 Dec 2018 15:34:00 -0500 Subject: [PATCH 4/8] Add comments to neighbor list code --- include/openmc/cell.h | 12 +++-- include/openmc/neighbor_list.h | 26 +++++++++-- include/openmc/openmp_interface.h | 77 +++++++++++++++++++++---------- 3 files changed, 81 insertions(+), 34 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 59b5e00a2..3a1011248 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -155,10 +155,11 @@ public: virtual ~Cell() {} }; +//============================================================================== + class CSGCell : public Cell { public: - CSGCell(); explicit CSGCell(pugi::xml_node cell_node); @@ -171,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 { @@ -185,11 +186,12 @@ public: moab::DagMC* dagmc_ptr_; DAGCell(); - std::pair 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 + distance(Position r, Direction u, int32_t on_surface) const; + void to_hdf5(hid_t group_id) const; }; #endif diff --git a/include/openmc/neighbor_list.h b/include/openmc/neighbor_list.h index 4b937a69f..0aeae0d6a 100644 --- a/include/openmc/neighbor_list.h +++ b/include/openmc/neighbor_list.h @@ -35,14 +35,20 @@ public: using suffix_iter = std::forward_list::iterator; using iterator = NeighborListIter; + //! 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 is assumed + //! that returning the transport calculation and possibly re-adding the + //! element later is faster than waiting on the lock to be released. void push(int new_elem) { // Try to acquire the lock. std::unique_lock lock(mutex_, std::try_to_lock); if (lock) { - // Make sure this element isn't already in the suffix. Don't check the - // prefix because we don't expect this function to be called if the "new" - // element is already there. + // It is possible another thread already added this element to the suffix + // while this thread was searching for a cell so make sure the given + // element isn't a duplicate before adding it. if (std::find(suffix_.cbegin(), suffix_.cend(), new_elem) == suffix_.cend()) { suffix_.push_front(new_elem); @@ -50,6 +56,12 @@ public: } } + //! Move data from the linked-list suffix to the consecutive vector prefix. + // + //! The consecutive data slightly improves runtime (likely due to cache + //! locality). Note that this function is not guaranteed threadsafe---the + //! caller is responsible for making sure only one thread at a time calls this + //! function. void make_consecutive() { while (!suffix_.empty()) { @@ -76,8 +88,11 @@ template class NeighborListIter { public: + // Construct from a prefix iterator. NeighborListIter(NeighborList* nl, T_prefix_iter it) { + // If we were given an iterator to the end of the prefix, immediately switch + // over to suffix mode. base_ = nl; if (it != base_->prefix_.end()) { in_prefix_ = true; @@ -88,6 +103,7 @@ public: } } + // Construct from a suffix iterator. NeighborListIter(NeighborList* nl, T_suffix_iter it) { in_prefix_ = false; @@ -139,11 +155,13 @@ public: private: NeighborList* base_; + // This type essentially wraps the implementation for two different external + // iterators. A union is used to contain one iterator or the other, and the + // in_prefix_ flag indicates which version of that union is valid. union { T_prefix_iter prefix_iter_; T_suffix_iter suffix_iter_; }; - bool in_prefix_; }; diff --git a/include/openmc/openmp_interface.h b/include/openmc/openmp_interface.h index 563500a69..ec2ef9946 100644 --- a/include/openmc/openmp_interface.h +++ b/include/openmc/openmp_interface.h @@ -7,43 +7,70 @@ namespace openmc { +//============================================================================== +//! An object used to prevent concurrent access to a piece of data. +// +//! This type meets the C++ "Lockable" requirements. +//============================================================================== + class ThreadMutex { - -#ifdef _OPENMP -//============================================================================== -// Implementation for when OpenMP is actually defined. - public: ThreadMutex() - {omp_init_lock(&mutex_);} + { + #ifdef _OPEMP + omp_init_lock(&mutex_); + #endif + } ~ThreadMutex() - {omp_destroy_lock(&mutex_);} + { + #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. + ThreadMutex(const ThreadMutex&) = delete; + ThreadMutex& operator= (const ThreadMutex&) = delete; + + //! Lock the mutex. + // + //! This function blocks execution until the lock succeeds. void lock() - {omp_set_lock(&mutex_);} + { + #ifdef _OPEMP + omp_set_lock(&mutex_); + #endif + } - bool try_lock() - {return omp_test_lock(&mutex_);} + //! 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 + } - void unlock() - {omp_unset_lock(&mutex_);} + //! Unlock the mutex. + void unlock() noexcept + { + #ifdef _OPEMP + omp_unset_lock(&mutex_); + #endif + } private: - omp_lock_t mutex_; - -#else -//============================================================================== -// Empty implementation for when OpenMP is not defined. - -public: - ThreadMutex() {} - ~ThreadMutex() = default; - void lock() {} - bool try_lock() {return true;} - void unlock() {} -#endif + #ifdef _OPEMP + omp_lock_t mutex_; + #endif }; } // namespace openmc From 35238259c0d476e71ba51d96e34d7bf280ffd43a Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 18 Dec 2018 15:02:07 -0500 Subject: [PATCH 5/8] Simplify neighbor list implementation --- include/openmc/cell.h | 2 +- include/openmc/geometry_aux.h | 9 -- include/openmc/neighbor_list.h | 159 ++++----------------------------- src/geometry.cpp | 10 +-- src/geometry_aux.cpp | 8 -- src/simulation.cpp | 4 - 6 files changed, 24 insertions(+), 168 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 3a1011248..c90cb7f69 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -105,7 +105,7 @@ public: bool simple_; //!< Does the region contain only intersections? //! \brief Neighboring cells in the same universe. - NeighborList neighbors; + NeighborList neighbors_; Position translation_ {0, 0, 0}; //!< Translation vector for filled universe diff --git a/include/openmc/geometry_aux.h b/include/openmc/geometry_aux.h index 4e634ebca..da55c3d9b 100644 --- a/include/openmc/geometry_aux.h +++ b/include/openmc/geometry_aux.h @@ -88,15 +88,6 @@ distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset); extern "C" int maximum_levels(int32_t univ); -//============================================================================== -//! Perform any geometry operations that should be done between generations. -//! -//! Currently, this just includes converting cell neighbor linked lists into -//! vectors, but similar optimizations could be added here in the future. -//============================================================================== - -void geometry_finalize_generation(); - //============================================================================== //! Deallocates global vectors and maps for cells, universes, and lattices. //============================================================================== diff --git a/include/openmc/neighbor_list.h b/include/openmc/neighbor_list.h index 0aeae0d6a..f22e873e1 100644 --- a/include/openmc/neighbor_list.h +++ b/include/openmc/neighbor_list.h @@ -10,174 +10,51 @@ namespace openmc { -// Forward declare the neighbor list iterator type. -template -class NeighborListIter; - //============================================================================== -//! An efficient, threadsafe, dynamic container for listing neighboring cells. +//! A threadsafe, dynamic container for listing neighboring cells. // -//! This container is a minor improvement upon a linked list. The linked list -//! allows for threadsafe dynamic growth; any number of threads can safely -//! read from the list without locks or reference counting. Write access must -//! be protected with a lock, but write events are fairly rare for neighbor -//! lists. This container also provides the option to flush the linked list -//! into a vector at some convenient time (like between generations where the -//! neighbor lists are not being read) for contiguous data. The resulting -//! performance benefit is small but consistent. +//! 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 prefix_iter = std::vector::iterator; - using suffix_iter = std::forward_list::iterator; - using iterator = NeighborListIter; + using const_iterator = std::forward_list::const_iterator; - //! Attempt to add an element. + // 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 is assumed - //! that returning the transport calculation and possibly re-adding the - //! element later is faster than waiting on the lock to be released. + // 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(int new_elem) { // Try to acquire the lock. std::unique_lock lock(mutex_, std::try_to_lock); if (lock) { - // It is possible another thread already added this element to the suffix + // 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(suffix_.cbegin(), suffix_.cend(), new_elem) - == suffix_.cend()) { - suffix_.push_front(new_elem); + if (std::find(list_.cbegin(), list_.cend(), new_elem) == list_.cend()) { + list_.push_front(new_elem); } } } - //! Move data from the linked-list suffix to the consecutive vector prefix. - // - //! The consecutive data slightly improves runtime (likely due to cache - //! locality). Note that this function is not guaranteed threadsafe---the - //! caller is responsible for making sure only one thread at a time calls this - //! function. - void make_consecutive() - { - while (!suffix_.empty()) { - prefix_.push_back(suffix_.front()); - suffix_.pop_front(); - } - } + const_iterator cbegin() const + {return list_.cbegin();} - iterator begin(); + const_iterator cend() const + {return list_.cend();} - iterator end(); private: - std::vector prefix_; - std::forward_list suffix_; + std::forward_list list_; ThreadMutex mutex_; - - friend class NeighborListIter; }; -//============================================================================== - -template -class NeighborListIter -{ -public: - // Construct from a prefix iterator. - NeighborListIter(NeighborList* nl, T_prefix_iter it) - { - // If we were given an iterator to the end of the prefix, immediately switch - // over to suffix mode. - base_ = nl; - if (it != base_->prefix_.end()) { - in_prefix_ = true; - prefix_iter_ = it; - } else { - in_prefix_ = false; - suffix_iter_ = base_->suffix_.begin(); - } - } - - // Construct from a suffix iterator. - NeighborListIter(NeighborList* nl, T_suffix_iter it) - { - in_prefix_ = false; - suffix_iter_ = it; - } - - T_value operator*() - { - if (in_prefix_) { - return *prefix_iter_; - } else { - return *suffix_iter_; - } - } - - bool operator==(const NeighborListIter& other) - { - if (in_prefix_ != other.in_prefix_) return false; - - if (in_prefix_) { - return prefix_iter_ == other.prefix_iter_; - } else { - return suffix_iter_ == other.suffix_iter_; - } - } - - bool operator!=(const NeighborListIter& other) - {return !(*this == other);} - - NeighborListIter& operator++() - { - if (in_prefix_) { - // We are in the prefix so increment the prefix iterator. - ++prefix_iter_; - - // If we've reached the end of the prefix, switch to the suffix iterator. - if (prefix_iter_ == base_->prefix_.end()) { - in_prefix_ = false; - suffix_iter_ = base_->suffix_.begin(); - } - - } else { - // We are in the suffix so increment the suffix iterator. - ++suffix_iter_; - } - return *this; - } - -private: - NeighborList* base_; - - // This type essentially wraps the implementation for two different external - // iterators. A union is used to contain one iterator or the other, and the - // in_prefix_ flag indicates which version of that union is valid. - union { - T_prefix_iter prefix_iter_; - T_suffix_iter suffix_iter_; - }; - bool in_prefix_; -}; - -//============================================================================== - -inline NeighborList::iterator -NeighborList::begin() -{ - return iterator(this, prefix_.begin()); -} - -inline NeighborList::iterator -NeighborList::end() -{ - return iterator(this, suffix_.end()); -} - } // namespace openmc #endif // OPENMC_NEIGHBOR_LIST_H diff --git a/src/geometry.cpp b/src/geometry.cpp index baf04f705..a517d8988 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -63,14 +63,14 @@ check_cell_overlap(Particle* p) //============================================================================== bool -find_cell_inner(Particle* p, NeighborList* neighbor_list) +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; if (neighbor_list) { - for (auto it = neighbor_list->begin(); it != neighbor_list->end(); ++it) { + for (auto it = neighbor_list->cbegin(); it != neighbor_list->cend(); ++it) { i_cell = *it; // Make sure the search cell is in the same universe. @@ -91,7 +91,7 @@ find_cell_inner(Particle* p, NeighborList* neighbor_list) } else { int i_universe = p->coord[p->n_coord-1].universe; const auto& cells {model::universes[i_universe]->cells_}; - for (auto it = cells.begin(); it != cells.end(); it++) { + for (auto it = cells.cbegin(); it != cells.cend(); it++) { i_cell = *it; // Make sure the search cell is in the same universe. @@ -290,14 +290,14 @@ find_cell(Particle* p, bool use_neighbor_lists) // Search for the particle in that cell's neighbor list. Return if we // found the particle. - bool found = find_cell_inner(p, &c.neighbors); + 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(p->coord[coord_lvl].cell); + if (found) c.neighbors_.push(p->coord[coord_lvl].cell); return found; } else { diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index d411411eb..c326744eb 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -432,14 +432,6 @@ maximum_levels(int32_t univ) //============================================================================== -void -geometry_finalize_generation() -{ - for (Cell* c : model::cells) c->neighbors.make_consecutive(); -} - -//============================================================================== - void free_memory_geometry_c() { diff --git a/src/simulation.cpp b/src/simulation.cpp index 15e2f7cef..50bf12521 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -4,7 +4,6 @@ #include "openmc/container_util.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" -#include "openmc/geometry_aux.h" #include "openmc/message_passing.h" #include "openmc/output.h" #include "openmc/particle.h" @@ -450,9 +449,6 @@ void finalize_generation() // For fixed-source mode, we need to sample the external source fill_source_bank_fixedsource(); } - - // Clean up neighbor lists - geometry_finalize_generation(); } void initialize_history(Particle* p, int64_t index_source) From e6466e9867cdb5dbd810d6e48ffd1c73fa5042d5 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 18 Dec 2018 15:06:22 -0500 Subject: [PATCH 6/8] Remove surface-based neighbor lists --- include/openmc/geometry_aux.h | 6 ------ include/openmc/surface.h | 3 --- src/geometry.F90 | 5 +---- src/geometry_aux.cpp | 27 --------------------------- src/input_xml.F90 | 10 +--------- 5 files changed, 2 insertions(+), 49 deletions(-) diff --git a/include/openmc/geometry_aux.h b/include/openmc/geometry_aux.h index da55c3d9b..e456f700c 100644 --- a/include/openmc/geometry_aux.h +++ b/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. //============================================================================== diff --git a/include/openmc/surface.h b/include/openmc/surface.h index d0ada79ea..3df791ba3 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -66,9 +66,6 @@ public: int bc_; //!< Boundary condition std::string name_; //!< User-defined name - std::vector neighbor_pos_; //!< List of cells on positive side - std::vector neighbor_neg_; //!< List of cells on negative side - explicit Surface(pugi::xml_node surf_node); Surface(); diff --git a/src/geometry.F90 b/src/geometry.F90 index 91ed6bbdc..83bf2325f 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -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) & @@ -97,7 +94,7 @@ contains if (present(use_neighbor_lists)) then found = find_cell_c(p, logical(use_neighbor_lists, kind=C_BOOL)) else - found = find_cell_c(p, logical(.false., kind=C_BOOL)) + found = find_cell_c(p, .false._C_BOOL) end if end subroutine find_cell diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index c326744eb..eca9b6c13 100644 --- a/src/geometry_aux.cpp +++ b/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() { diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 1da1bf9df..7bba72d15 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -9,7 +9,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 @@ -132,7 +131,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 @@ -182,12 +181,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() @@ -2297,7 +2290,6 @@ contains logical :: file_exists ! Does multipole library exist? character(7) :: readable ! Is multipole library readable? character(MAX_FILE_LEN) :: filename ! Path to multipole xs library - character(kind=C_CHAR), pointer :: string(:) integer(HID_T) :: file_id integer(HID_T) :: group_id From f5059cfacae0f89227bdfbef6051aa0e3ef69554 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 23 Dec 2018 22:47:02 -0500 Subject: [PATCH 7/8] Append rather than prepend to neighbor lists --- include/openmc/neighbor_list.h | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/include/openmc/neighbor_list.h b/include/openmc/neighbor_list.h index f22e873e1..abf2fab86 100644 --- a/include/openmc/neighbor_list.h +++ b/include/openmc/neighbor_list.h @@ -39,7 +39,15 @@ public: // 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()) { - list_.push_front(new_elem); + // 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); + } } } } From b178939cfe93081b5317b8014f18bf15a7cb8e83 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 29 Dec 2018 12:54:29 -0500 Subject: [PATCH 8/8] Address #1140 comments --- include/openmc/neighbor_list.h | 6 +++--- include/openmc/openmp_interface.h | 10 +++++----- src/geometry.cpp | 2 +- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/include/openmc/neighbor_list.h b/include/openmc/neighbor_list.h index abf2fab86..f6142101d 100644 --- a/include/openmc/neighbor_list.h +++ b/include/openmc/neighbor_list.h @@ -30,10 +30,10 @@ public: // 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(int new_elem) + void push_back(int new_elem) { // Try to acquire the lock. - std::unique_lock lock(mutex_, std::try_to_lock); + std::unique_lock 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 @@ -61,7 +61,7 @@ public: private: std::forward_list list_; - ThreadMutex mutex_; + OpenMPMutex mutex_; }; } // namespace openmc diff --git a/include/openmc/openmp_interface.h b/include/openmc/openmp_interface.h index ec2ef9946..076b19dc0 100644 --- a/include/openmc/openmp_interface.h +++ b/include/openmc/openmp_interface.h @@ -13,17 +13,17 @@ namespace openmc { //! This type meets the C++ "Lockable" requirements. //============================================================================== -class ThreadMutex +class OpenMPMutex { public: - ThreadMutex() + OpenMPMutex() { #ifdef _OPEMP omp_init_lock(&mutex_); #endif } - ~ThreadMutex() + ~OpenMPMutex() { #ifdef _OPEMP omp_destroy_lock(&mutex_); @@ -33,8 +33,8 @@ public: // 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. - ThreadMutex(const ThreadMutex&) = delete; - ThreadMutex& operator= (const ThreadMutex&) = delete; + OpenMPMutex(const OpenMPMutex&) = delete; + OpenMPMutex& operator= (const OpenMPMutex&) = delete; //! Lock the mutex. // diff --git a/src/geometry.cpp b/src/geometry.cpp index a517d8988..2d3aca8c2 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -297,7 +297,7 @@ find_cell(Particle* p, bool use_neighbor_lists) // 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(p->coord[coord_lvl].cell); + if (found) c.neighbors_.push_back(p->coord[coord_lvl].cell); return found; } else {