diff --git a/include/openmc/cell.h b/include/openmc/cell.h index f2c8f0f71..c90cb7f69 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 @@ -103,6 +104,9 @@ public: std::vector 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 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/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/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/neighbor_list.h b/include/openmc/neighbor_list.h new file mode 100644 index 000000000..f6142101d --- /dev/null +++ b/include/openmc/neighbor_list.h @@ -0,0 +1,68 @@ +#ifndef OPENMC_NEIGHBOR_LIST_H +#define OPENMC_NEIGHBOR_LIST_H + +#include +#include +#include +#include + +#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::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 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 list_; + OpenMPMutex mutex_; +}; + +} // 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..076b19dc0 --- /dev/null +++ b/include/openmc/openmp_interface.h @@ -0,0 +1,77 @@ +#ifndef OPENMC_OPENMP_INTERFACE_H +#define OPENMC_OPENMP_INTERFACE_H + +#ifdef _OPENMP +#include +#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 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 72caab359..83bf2325f 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 @@ -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 diff --git a/src/geometry.cpp b/src/geometry.cpp index fe3d52c68..2d3aca8c2 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,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* 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]) { 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 190ce52ab..d723aae99 100644 --- a/src/input_xml.F90 +++ b/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() diff --git a/src/tracking.F90 b/src/tracking.F90 index d345d7cc0..b94b35dba 100644 --- a/src/tracking.F90 +++ b/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 ! ==========================================================================