crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Cleanup UniversePartitioner code

Browse source
This commit is contained in:
Sterling Harper 2019-03-31 12:29:18 -04:00
parent 4b76c28a3b
commit 5d5b82f0ae
6 changed files with 117 additions and 72 deletions
Download patch file Download diff file Expand all files Collapse all files

41
include/openmc/cell.h
View file

@ -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<int32_t>& get_cells(Position r, Direction u) const
{return get_cells(r, u, 0, surfs_.size()-1, (surfs_.size()-1)/2);}
const std::vector<int32_t>& 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<int32_t>&
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<std::vector<int32_t>> cells_;
//! A sorted vector of indices to surfaces that partition the universe
std::vector<int32_t> 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<std::vector<int32_t>> partitions_;
};
//==============================================================================

5
include/openmc/geometry_aux.h
View file

@ -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();
//==============================================================================

38
include/openmc/surface.h
View file

@ -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_;
};
//==============================================================================

29
src/cell.cpp
View file

@ -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<int32_t> 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);
}
}

11
src/geometry.cpp
View file

@ -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;

65
src/geometry_aux.cpp
View file

@ -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<int32_t> 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<const SurfaceZPlane*>(model::surfaces[i_surf].get())) {
++n_zplanes;
if (n_zplanes > 5) {
univ->partitioner_ = std::make_unique<UniversePartitioner>(
UniversePartitioner(*univ));
break;
}
}
}
}
}
}
//==============================================================================
void
@ -215,6 +245,7 @@ void finalize_geometry(std::vector<std::vector<double>>& 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<std::vector<double>>& 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<int32_t> 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<const SurfaceZPlane*>(model::surfaces[i_surf].get())) {
++n_zplanes;
if (n_zplanes > 5) {
univ->partitioner_ = std::make_unique<UniversePartitioner>(
UniversePartitioner(*univ));
break;
}
}
}
}
}
}
//==============================================================================
@ -546,8 +555,6 @@ maximum_levels(int32_t univ)
void
free_memory_geometry()
{
std::cout << "SEARCH EFFICIENCY = " << static_cast<double>(search_numerator1) / search_denominator1 << "\n";
std::cout << "SEARCH EFFICIENCY = " << static_cast<double>(search_numerator2) / search_denominator2 << "\n";
model::cells.clear();
model::cell_map.clear();