diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index 107cc7d1f3e..466673b6636 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -3,6 +3,8 @@ #include #include +#include +#include #include "openmc/array.h" #include "openmc/constants.h" @@ -13,6 +15,19 @@ namespace openmc { class BoundaryInfo; class GeometryState; +// OverlapKey to store cell and universe data of a single overlap +struct OverlapKey { + int universe_id; + int cell1_id; + int cell2_id; + + bool operator==(const OverlapKey& other) const + { + return universe_id == other.universe_id && cell1_id == other.cell1_id && + cell2_id == other.cell2_id; + } +}; + //============================================================================== // Global variables //============================================================================== @@ -24,6 +39,10 @@ extern "C" int n_coord_levels; //!< Number of CSG coordinate levels extern vector overlap_check_count; +// Overlap data structures get cleared every slice_data run +extern vector overlap_keys; +extern std::unordered_map overlap_key_index; + } // namespace model //============================================================================== @@ -38,8 +57,7 @@ inline bool coincident(double d1, double d2) //============================================================================== //! Check for overlapping cells at a particle's position. //============================================================================== - -bool check_cell_overlap(GeometryState& p, bool error = true); +size_t check_cell_overlap(GeometryState& p, bool error = true); //============================================================================== //! Get the cell instance for a particle at the specified universe level @@ -79,4 +97,16 @@ BoundaryInfo distance_to_boundary(GeometryState& p); } // namespace openmc +// Hash specialization for use in std::unordered_map +template<> +struct std::hash { + size_t operator()(const openmc::OverlapKey& k) const noexcept + { + size_t h = std::hash {}(k.universe_id); + h ^= std::hash {}(k.cell1_id) + 0x9e3779b9 + (h << 6) + (h >> 2); + h ^= std::hash {}(k.cell2_id) + 0x9e3779b9 + (h << 6) + (h >> 2); + return h; + } +}; + #endif // OPENMC_GEOMETRY_H diff --git a/include/openmc/plot.h b/include/openmc/plot.h index f97d313847f..de56fe19b6d 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -5,6 +5,7 @@ #include #include #include +#include #include "openmc/tensor.h" #include "pugixml.hpp" @@ -155,7 +156,7 @@ struct IdData { // Methods void set_value(size_t y, size_t x, const Particle& p, int level, Filter* filter = nullptr, FilterMatch* match = nullptr); - void set_overlap(size_t y, size_t x); + void set_overlap(size_t y, size_t x, size_t overlap_idx); // Members tensor::Tensor data_; //!< 2D array of cell & material ids @@ -168,7 +169,7 @@ struct PropertyData { // Methods void set_value(size_t y, size_t x, const Particle& p, int level, Filter* filter = nullptr, FilterMatch* match = nullptr); - void set_overlap(size_t y, size_t x); + void set_overlap(size_t y, size_t x, size_t overlap_idx); // Members tensor::Tensor data_; //!< 2D array of temperature & density data @@ -181,7 +182,7 @@ struct RasterData { // Methods void set_value(size_t y, size_t x, const Particle& p, int level, Filter* filter = nullptr, FilterMatch* match = nullptr); - void set_overlap(size_t y, size_t x); + void set_overlap(size_t y, size_t x, size_t overlap_idx); // Members tensor::Tensor @@ -278,8 +279,11 @@ T SlicePlotBase::get_map(int32_t filter_index) const if (found_cell) { data.set_value(y, x, p, j, filter, &match); } - if (show_overlaps_ && check_cell_overlap(p, false)) { - data.set_overlap(y, x); + if (show_overlaps_) { + auto overlap_idx = check_cell_overlap(p, false); + if (overlap_idx != SIZE_MAX) { + data.set_overlap(y, x, overlap_idx); + } } } // inner for } diff --git a/openmc/lib/plot.py b/openmc/lib/plot.py index 44d6ac273d3..bff33602073 100644 --- a/openmc/lib/plot.py +++ b/openmc/lib/plot.py @@ -1,6 +1,7 @@ from collections.abc import Mapping from ctypes import (c_bool, c_int, c_size_t, c_int32, - c_double, c_uint8, Structure, POINTER) + c_double, c_uint8, Structure, POINTER, + byref) from weakref import WeakValueDictionary from ..exceptions import AllocationError, InvalidIDError @@ -255,6 +256,33 @@ def property_map(plot): return prop_data +# Python wrappings for overlap functions + +def slice_data_overlap_count(): + count = c_size_t() + _dll.openmc_slice_data_overlap_count(byref(count)) + return count.value + +def slice_data_overlap_info(): + n = slice_data_overlap_count() + overlap_info = np.empty(n * 3, dtype=np.int32) + + if n > 0: + _dll.openmc_slice_data_overlap_info( + n, + overlap_info.ctypes.data_as(POINTER(c_int32)), + ) + + return overlap_info, n + +_dll.openmc_slice_data_overlap_count.argtypes = [POINTER(c_size_t)] +_dll.openmc_slice_data_overlap_count.restype = c_int +_dll.openmc_slice_data_overlap_count.errcheck = _error_handler + +_dll.openmc_slice_data_overlap_info.argtypes = [c_size_t, POINTER(c_int32)] +_dll.openmc_slice_data_overlap_info.restype = c_int +_dll.openmc_slice_data_overlap_info.errcheck = _error_handler + _dll.openmc_get_plot_index.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_get_plot_index.restype = c_int _dll.openmc_get_plot_index.errcheck = _error_handler diff --git a/src/geometry.cpp b/src/geometry.cpp index ddb61385f18..cdd09f12e02 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -26,16 +26,22 @@ int n_coord_levels; vector overlap_check_count; +vector overlap_keys; +std::unordered_map overlap_key_index; + } // namespace model //============================================================================== // Non-member functions //============================================================================== -bool check_cell_overlap(GeometryState& p, bool error) +size_t check_cell_overlap(GeometryState& p, bool error) { int n_coord = p.n_coord(); + // If no overlap found, return a nonphysical index + size_t overlap_index = SIZE_MAX; + // Loop through each coordinate level for (int j = 0; j < n_coord; j++) { Universe& univ = *model::universes[p.coord(j).universe()]; @@ -44,21 +50,39 @@ bool check_cell_overlap(GeometryState& p, bool error) for (auto index_cell : univ.cells_) { Cell& c = *model::cells[index_cell]; if (c.contains(p.coord(j).r(), p.coord(j).u(), p.surface())) { +#pragma omp atomic + ++model::overlap_check_count[index_cell]; + if (index_cell != p.coord(j).cell()) { if (error) { fatal_error( fmt::format("Overlapping cells detected: {}, {} on universe {}", c.id_, model::cells[p.coord(j).cell()]->id_, univ.id_)); } - return true; + + // With no fatal error (plotter is calling), now adds overlaps and + // calls them; ensures order does not matter when making overlap key + int cell_a = model::cells[index_cell]->id_; + int cell_b = model::cells[p.coord(j).cell()]->id_; + int a = std::min(cell_a, cell_b); + int b = std::max(cell_a, cell_b); + OverlapKey key {univ.id_, a, b}; + + auto it = model::overlap_key_index.find(key); + if (it != model::overlap_key_index.end()) { + overlap_index = it->second; // already exists, reuse index + } else { + size_t idx = model::overlap_keys.size(); + model::overlap_keys.push_back(key); + model::overlap_key_index[key] = idx; + overlap_index = idx; + } + break; } -#pragma omp atomic - ++model::overlap_check_count[index_cell]; } } } - - return false; + return overlap_index; } //============================================================================== diff --git a/src/plot.cpp b/src/plot.cpp index b9a1136afdb..f2fc79a960e 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -73,7 +73,7 @@ void IdData::set_value(size_t y, size_t x, const Particle& p, int level, } } -void IdData::set_overlap(size_t y, size_t x) +void IdData::set_overlap(size_t y, size_t x, size_t /*overlap_idx*/) { for (size_t k = 0; k < data_.shape(2); ++k) data_(y, x, k) = OVERLAP; @@ -94,7 +94,7 @@ void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level, } } -void PropertyData::set_overlap(size_t y, size_t x) +void PropertyData::set_overlap(size_t y, size_t x, size_t /*overlap_idx*/) { data_(y, x) = OVERLAP; } @@ -154,12 +154,12 @@ void RasterData::set_value(size_t y, size_t x, const Particle& p, int level, } } -void RasterData::set_overlap(size_t y, size_t x) +void RasterData::set_overlap(size_t y, size_t x, size_t overlap_idx) { // Set cell, instance, and material to OVERLAP, but preserve filter bin id_data_(y, x, 0) = OVERLAP; id_data_(y, x, 1) = OVERLAP; - id_data_(y, x, 2) = OVERLAP; + id_data_(y, x, 2) = OVERLAP - overlap_idx - 1; // Note: id_data_(y, x, 3) is NOT overwritten - preserves filter bin for tally // plotting @@ -177,6 +177,8 @@ std::unordered_map plot_map; vector> plots; uint64_t plotter_seed = 1; +std::unique_ptr last_slice_data; + } // namespace model //============================================================================== @@ -1980,6 +1982,10 @@ extern "C" int openmc_slice_data(const double origin[3], const double u_span[3], model::overlap_check_count.resize(model::cells.size()); } + if (color_overlaps) { + settings::check_overlaps = true; + } + try { // Create a temporary SlicePlotBase object to reuse get_map logic SlicePlotBase plot_params; @@ -1991,10 +1997,12 @@ extern "C" int openmc_slice_data(const double origin[3], const double u_span[3], plot_params.show_overlaps_ = color_overlaps; plot_params.slice_level_ = level; + // Clear overlap data structures on new slice call + model::overlap_keys.clear(); + model::overlap_key_index.clear(); + // Use get_map to generate data auto data = plot_params.get_map(filter_index); - - // Copy geometry data std::copy(data.id_data_.begin(), data.id_data_.end(), geom_data); // Copy property data if requested @@ -2002,6 +2010,7 @@ extern "C" int openmc_slice_data(const double origin[3], const double u_span[3], std::copy( data.property_data_.begin(), data.property_data_.end(), property_data); } + } catch (const std::exception& e) { set_errmsg(e.what()); return OPENMC_E_UNASSIGNED; @@ -2010,6 +2019,32 @@ extern "C" int openmc_slice_data(const double origin[3], const double u_span[3], return 0; } +// Gets the number of overlaps that we need data for +extern "C" int openmc_slice_data_overlap_count(size_t* count) +{ + if (!count) { + set_errmsg("Null pointer passed for overlap count."); + return OPENMC_E_INVALID_ARGUMENT; + } + *count = model::overlap_keys.size(); + + return 0; +} + +// Plotter pre-allocates array size based on what is returned with +// overlap_count; populates an array of size 3*count +extern "C" int openmc_slice_data_overlap_info( + size_t count, int32_t* overlap_info) +{ + for (size_t i = 0; i < count; ++i) { + overlap_info[i * 3] = model::overlap_keys[i].universe_id; + overlap_info[i * 3 + 1] = model::overlap_keys[i].cell1_id; + overlap_info[i * 3 + 2] = model::overlap_keys[i].cell2_id; + } + + return 0; +} + extern "C" int openmc_get_plot_index(int32_t id, int32_t* index) { auto it = model::plot_map.find(id);