diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 355e23ff3..2cd0bd3fa 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -196,40 +196,69 @@ public: }; #endif +//============================================================================== +//! Speeds up geometry searches by grouping cells in a search tree. +// +//! Currently this object only works with universes that are divided up by a +//! bunch of z-planes. It could be generalized to other planes, cylinders, +//! and spheres. //============================================================================== class UniversePartitioner { public: - UniversePartitioner(const Universe& univ); + explicit UniversePartitioner(const Universe& univ); + //! Return the list of cells that could contain the given coordinates. const std::vector& get_cells(Position r, Direction u) const {return get_cells(r, u, 0, surfs_.size()-1, (surfs_.size()-1)/2);} - const std::vector& get_cells(Position r, Direction u, - int left, int right, int middle) const +private: + //! Perform a binary search for the given coordinates. + // + //! \param left The lower search bound of the `surfs_` vector + //! \param right The upper search bound of the `surfs_` vector + //! \param middle The index of the `Surface` to test against + const std::vector& + get_cells(Position r, Direction u, int left, int right, int middle) const { + // Check the sense of the coordinates for the current surface. auto i_surf = surfs_[middle]; const auto& surf = *model::surfaces[i_surf]; if (surf.sense(r, u)) { + // The coordinates lie in the positive halfspace. Recurse if there are + // more surfaces to check. Otherwise, return the cells on the positive + // side of this surface. int right_leaf = right - (right - middle) / 2; if (right_leaf != middle) { return get_cells(r, u, middle+1, right, right_leaf); } else { - return cells_[middle+1]; + return partitions_[middle+1]; } } else { + // The coordinates lie in the negative halfspace. Recurse if there are + // more surfaces to check. Otherwise, return the cells on the negative + // side of this surface. int left_leaf = left + (middle - left) / 2; if (left_leaf != middle) { return get_cells(r, u, left, middle-1, left_leaf); } else { - return cells_[middle]; + return partitions_[middle]; } } } - std::vector> cells_; + //! A sorted vector of indices to surfaces that partition the universe std::vector surfs_; + + //! Vectors listing the indices of the cells that lie within each partition + // + //! There are n+1 partitions with n surfaces. `partitions_.front()` gives the + //! cells that lie on the negative side of `surfs_.front()`. + //! `partitions_.back()` gives the cells that lie on the positive side of + //! `surfs_.back()`. Otherwise, `partitions_[i]` gives cells sandwiched + //! between `surfs_[i-1]` and `surfs_[i]`. + std::vector> partitions_; }; //============================================================================== diff --git a/include/openmc/geometry_aux.h b/include/openmc/geometry_aux.h index 8236af2f7..ffcdcf885 100644 --- a/include/openmc/geometry_aux.h +++ b/include/openmc/geometry_aux.h @@ -10,11 +10,6 @@ namespace openmc { -extern int search_numerator1; -extern int search_denominator1; -extern int search_numerator2; -extern int search_denominator2; - void read_geometry_xml(); //============================================================================== diff --git a/include/openmc/surface.h b/include/openmc/surface.h index c862e14c9..4db397a59 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -180,7 +180,6 @@ public: class SurfaceXPlane : public PeriodicSurface { - double x0_; public: explicit SurfaceXPlane(pugi::xml_node surf_node); double evaluate(Position r) const; @@ -190,6 +189,8 @@ public: bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const; BoundingBox bounding_box() const; + + double x0_; }; //============================================================================== @@ -200,7 +201,6 @@ public: class SurfaceYPlane : public PeriodicSurface { - double y0_; public: explicit SurfaceYPlane(pugi::xml_node surf_node); double evaluate(Position r) const; @@ -210,6 +210,8 @@ public: bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const; BoundingBox bounding_box() const; + + double y0_; }; //============================================================================== @@ -221,7 +223,6 @@ public: class SurfaceZPlane : public PeriodicSurface { public: - double z0_; explicit SurfaceZPlane(pugi::xml_node surf_node); double evaluate(Position r) const; double distance(Position r, Direction u, bool coincident) const; @@ -230,6 +231,8 @@ public: bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const; BoundingBox bounding_box() const; + + double z0_; }; //============================================================================== @@ -240,7 +243,6 @@ public: class SurfacePlane : public PeriodicSurface { - double A_, B_, C_, D_; public: explicit SurfacePlane(pugi::xml_node surf_node); double evaluate(Position r) const; @@ -250,6 +252,8 @@ public: bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const; BoundingBox bounding_box() const; + + double A_, B_, C_, D_; }; //============================================================================== @@ -261,13 +265,14 @@ public: class SurfaceXCylinder : public CSGSurface { - double y0_, z0_, radius_; public: explicit SurfaceXCylinder(pugi::xml_node surf_node); double evaluate(Position r) const; double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + + double y0_, z0_, radius_; }; //============================================================================== @@ -279,13 +284,14 @@ public: class SurfaceYCylinder : public CSGSurface { - double x0_, z0_, radius_; public: explicit SurfaceYCylinder(pugi::xml_node surf_node); double evaluate(Position r) const; double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + + double x0_, z0_, radius_; }; //============================================================================== @@ -297,13 +303,14 @@ public: class SurfaceZCylinder : public CSGSurface { - double x0_, y0_, radius_; public: explicit SurfaceZCylinder(pugi::xml_node surf_node); double evaluate(Position r) const; double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + + double x0_, y0_, radius_; }; //============================================================================== @@ -315,13 +322,14 @@ public: class SurfaceSphere : public CSGSurface { - double x0_, y0_, z0_, radius_; public: explicit SurfaceSphere(pugi::xml_node surf_node); double evaluate(Position r) const; double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + + double x0_, y0_, z0_, radius_; }; //============================================================================== @@ -333,13 +341,14 @@ public: class SurfaceXCone : public CSGSurface { - double x0_, y0_, z0_, radius_sq_; public: explicit SurfaceXCone(pugi::xml_node surf_node); double evaluate(Position r) const; double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + + double x0_, y0_, z0_, radius_sq_; }; //============================================================================== @@ -351,13 +360,14 @@ public: class SurfaceYCone : public CSGSurface { - double x0_, y0_, z0_, radius_sq_; public: explicit SurfaceYCone(pugi::xml_node surf_node); double evaluate(Position r) const; double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + + double x0_, y0_, z0_, radius_sq_; }; //============================================================================== @@ -369,13 +379,14 @@ public: class SurfaceZCone : public CSGSurface { - double x0_, y0_, z0_, radius_sq_; public: explicit SurfaceZCone(pugi::xml_node surf_node); double evaluate(Position r) const; double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + + double x0_, y0_, z0_, radius_sq_; }; //============================================================================== @@ -386,14 +397,15 @@ public: class SurfaceQuadric : public CSGSurface { - // Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0 - double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_; public: explicit SurfaceQuadric(pugi::xml_node surf_node); double evaluate(Position r) const; double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + + // Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0 + double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_; }; //============================================================================== diff --git a/src/cell.cpp b/src/cell.cpp index db1190813..5e803c4cd 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -636,10 +636,13 @@ void DAGCell::to_hdf5(hid_t group_id) const { return; } #endif //============================================================================== +// UniversePartitioner implementation //============================================================================== UniversePartitioner::UniversePartitioner(const Universe& univ) { + // Find all of the z-planes in this universe. Add them to the surfs_ member. + // Use surf_set for O(1) searches that ensure surfs_ entries are not repeated. std::unordered_set surf_set; for (auto i_cell : univ.cells_) { for (auto token : model::cells[i_cell]->rpn_) { @@ -656,6 +659,8 @@ UniversePartitioner::UniversePartitioner(const Universe& univ) } } + // Define a functor for comparing z-planes by their location, and use it to + // sort the surfs_ member. struct { bool operator()(int32_t i_surf, int32_t j_surf) const { @@ -670,9 +675,10 @@ UniversePartitioner::UniversePartitioner(const Universe& univ) } compare_surfs; std::sort(surfs_.begin(), surfs_.end(), compare_surfs); - cells_.resize(surfs_.size() + 1); - + // Populate the partition lists. + partitions_.resize(surfs_.size() + 1); for (auto i_cell : univ.cells_) { + // Find the tokens for bounding z-planes. int32_t min_token = 0, max_token = 0; double min_z, max_z; for (auto token : model::cells[i_cell]->rpn_) { @@ -691,25 +697,28 @@ UniversePartitioner::UniversePartitioner(const Universe& univ) } } + // If there are no bounding z-planes, add this cell to all partitions. if (min_token == 0) { - min_token = -(surfs_.front() + 1); - max_token = surfs_.back() + 1; + for (auto& p : partitions_) p.push_back(i_cell); + continue; } + // Iterate over partitions, and add this cell to each appropriate one. + // Since surfs_ is sorted, we know this cell belongs to every partition + // between min_token and max_token. bool in_partition = min_token < 0; - for (auto i = 0; i < surfs_.size(); ++i) { if (in_partition) { - cells_[i].push_back(i_cell); - - if (max_token == -(surfs_[i] + 1)) + partitions_[i].push_back(i_cell); + if (max_token == -(surfs_[i] + 1)) { in_partition = false; + break; + } } else { in_partition = min_token == surfs_[i] + 1; } } - if (in_partition) - cells_.back().push_back(i_cell); + if (in_partition) partitions_.back().push_back(i_cell); } } diff --git a/src/geometry.cpp b/src/geometry.cpp index 5e6ea34ca..bbd7b9d0e 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -11,8 +11,6 @@ #include "openmc/simulation.h" #include "openmc/surface.h" -#include "openmc/geometry_aux.h" - namespace openmc { @@ -71,7 +69,6 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) bool found = false; int32_t i_cell; if (neighbor_list) { - ++search_numerator1; for (auto it = neighbor_list->cbegin(); it != neighbor_list->cend(); ++it) { i_cell = *it; @@ -83,7 +80,6 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) Position r {p->r_local()}; Direction u {p->u_local()}; auto surf = p->surface_; - ++search_denominator1; if (model::cells[i_cell]->contains(r, u, surf)) { p->coord_[p->n_coord_-1].cell = i_cell; found = true; @@ -92,13 +88,11 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) } } else { - ++search_numerator2; int i_universe = p->coord_[p->n_coord_-1].universe; - //const auto& cells {model::universes[i_universe]->cells_}; const auto& univ {*model::universes[i_universe]}; const auto& cells { - !univ.partitioner_ ? - model::universes[i_universe]->cells_ + !univ.partitioner_ + ? model::universes[i_universe]->cells_ : univ.partitioner_->get_cells(p->r_local(), p->u_local()) }; for (auto it = cells.cbegin(); it != cells.cend(); it++) { @@ -112,7 +106,6 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) Position r {p->r_local()}; Direction u {p->u_local()}; auto surf = p->surface_; - ++search_denominator2; if (model::cells[i_cell]->contains(r, u, surf)) { p->coord_[p->n_coord_-1].cell = i_cell; found = true; diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 1108bca29..38d3a8a0e 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -23,11 +23,6 @@ namespace openmc { -int search_denominator1 {0}; -int search_numerator1 {0}; -int search_denominator2 {0}; -int search_numerator2 {0}; - void read_geometry_xml() { #ifdef DAGMC @@ -125,6 +120,41 @@ adjust_indices() } } +//============================================================================== +//! Partition some universes with many z-planes for faster find_cell searches. + +void +partition_universes() +{ + // Iterate over universes with more than 10 cells. (Fewer than 10 is likely + // not worth partitioning.) + for (const auto& univ : model::universes) { + if (univ->cells_.size() > 10) { + // Collect the set of surfaces in this universe. + std::unordered_set surf_inds; + for (auto i_cell : univ->cells_) { + for (auto token : model::cells[i_cell]->rpn_) { + if (token < OP_UNION) surf_inds.insert(std::abs(token) - 1); + } + } + + // Partition the universe if there are more than 5 z-planes. (Fewer than + // five is likely no worth it.) + int n_zplanes = 0; + for (auto i_surf : surf_inds) { + if (dynamic_cast(model::surfaces[i_surf].get())) { + ++n_zplanes; + if (n_zplanes > 5) { + univ->partitioner_ = std::make_unique( + UniversePartitioner(*univ)); + break; + } + } + } + } + } +} + //============================================================================== void @@ -215,6 +245,7 @@ void finalize_geometry(std::vector>& nuc_temps, // Perform some final operations to set up the geometry adjust_indices(); count_cell_instances(model::root_universe); + partition_universes(); // Assign temperatures to cells that don't have temperatures already assigned assign_temperatures(); @@ -224,28 +255,6 @@ void finalize_geometry(std::vector>& nuc_temps, // Determine number of nested coordinate levels in the geometry model::n_coord_levels = maximum_levels(model::root_universe); - - for (const auto& univ : model::universes) { - if (univ->cells_.size() > 10) { - std::unordered_set surf_inds; - for (auto i_cell : univ->cells_) { - for (auto token : model::cells[i_cell]->rpn_) { - if (token < OP_UNION) surf_inds.insert(std::abs(token) - 1); - } - } - int n_zplanes = 0; - for (auto i_surf : surf_inds) { - if (dynamic_cast(model::surfaces[i_surf].get())) { - ++n_zplanes; - if (n_zplanes > 5) { - univ->partitioner_ = std::make_unique( - UniversePartitioner(*univ)); - break; - } - } - } - } - } } //============================================================================== @@ -546,8 +555,6 @@ maximum_levels(int32_t univ) void free_memory_geometry() { - std::cout << "SEARCH EFFICIENCY = " << static_cast(search_numerator1) / search_denominator1 << "\n"; - std::cout << "SEARCH EFFICIENCY = " << static_cast(search_numerator2) / search_denominator2 << "\n"; model::cells.clear(); model::cell_map.clear();