diff --git a/include/openmc/neighbor_list.h b/include/openmc/neighbor_list.h index f781424ff59..5b53be78d19 100644 --- a/include/openmc/neighbor_list.h +++ b/include/openmc/neighbor_list.h @@ -2,8 +2,9 @@ #define OPENMC_NEIGHBOR_LIST_H #include +#include #include -#include +#include #include #include "openmc/openmp_interface.h" @@ -13,51 +14,128 @@ 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. +//! This container is a reduced interface for a growable array 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. +// +//! Entries are published through an atomic length: a writer fills the next slot +//! and then releases the new length, and a reader acquires the length before +//! reading any slot. That pairing is what makes the lock-free reads well +//! defined. The first INLINE_CAPACITY entries live in the object itself, so a +//! typical cell never allocates and a search reads the entries from the same +//! cache line as the length. Beyond that the entries are copied into a heap +//! block, which doubles in capacity as needed; a replaced block is retained so +//! that a reader still traversing it stays valid, and since the capacity +//! doubles, the retained blocks hold fewer entries than the live one. //============================================================================== class NeighborList { public: using value_type = int32_t; - using const_iterator = std::forward_list::const_iterator; + + //! Number of entries held in the object itself, before any allocation + static constexpr int INLINE_CAPACITY {8}; + + //! A snapshot of the entries that were published when it was taken + class View { + public: + View(const value_type* data, int size) : data_ {data}, size_ {size} {} + const value_type* begin() const { return data_; } + const value_type* end() const { return data_ + size_; } + int size() const { return size_; } + + private: + const value_type* data_; + int size_; + }; + + NeighborList() = default; + NeighborList(const NeighborList&) = delete; + NeighborList& operator=(const NeighborList&) = delete; + + ~NeighborList() + { + void* block = block_.load(std::memory_order_relaxed); + while (block) { + void* previous = *static_cast(block); + std::free(block); + block = previous; + } + } + + //! Take a snapshot of the entries for lock-free traversal + // + //! The length is read first: a reader that observes a length also observes + //! the storage that held that many entries, since a spill to a larger block + //! is published before the length that counts the entry which caused it. + View view() const + { + int size = length_.load(std::memory_order_acquire); + void* block = block_.load(std::memory_order_acquire); + return View {block ? entries(block) : inline_, size}; + } // 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) + // found that returning to the transport calculation and possibly re-adding + // the element later is slightly faster than waiting on the lock. + void push_back(value_type 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); - } - } - } - } + if (!lock) + return; + + // Only one thread writes at a time and taking the lock orders this thread + // against previous writers, so the entries can be read without ordering. + int size = length_.load(std::memory_order_relaxed); + void* block = block_.load(std::memory_order_relaxed); + value_type* data = block ? entries(block) : inline_; + + // 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(data, data + size, new_elem) != data + size) + return; + + if (size == capacity_) { + int capacity = 2 * capacity_; + void* new_block = + std::malloc(sizeof(void*) + capacity * sizeof(value_type)); + if (!new_block) + return; - const_iterator cbegin() const { return list_.cbegin(); } + // Chain the old block so that every block is freed on destruction, and + // copy the entries a reader may already be traversing. + *static_cast(new_block) = block; + value_type* new_data = entries(new_block); + std::copy(data, data + size, new_data); - const_iterator cend() const { return list_.cend(); } + block_.store(new_block, std::memory_order_release); + capacity_ = capacity; + data = new_data; + } + + // Write the element before publishing it, so a reader that sees the new + // length is guaranteed to see the element too. + data[size] = new_elem; + length_.store(size + 1, std::memory_order_release); + } private: - std::forward_list list_; + //! Entries follow the pointer to the previous block at the start of a block + static value_type* entries(void* block) + { + return reinterpret_cast(static_cast(block) + 1); + } + + value_type inline_[INLINE_CAPACITY]; //!< entries before any allocation + std::atomic block_ {nullptr}; //!< spilled entries, chained to older + std::atomic length_ {0}; //!< entries published in that storage + int capacity_ {INLINE_CAPACITY}; //!< slots in that storage; writers only OpenMPMutex mutex_; }; diff --git a/src/geometry.cpp b/src/geometry.cpp index ecacf0bffbf..da328884fce 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -132,7 +132,9 @@ bool find_cell_inner( bool found = false; int32_t i_cell = C_NONE; if (neighbor_list) { - for (auto it = neighbor_list->cbegin(); it != neighbor_list->cend(); ++it) { + // Take one snapshot of the list rather than re-reading it each iteration + auto neighbors = neighbor_list->view(); + for (auto it = neighbors.begin(); it != neighbors.end(); ++it) { i_cell = *it; // Make sure the search cell is in the same universe.