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

Address #1210 comments

Browse source
This commit is contained in:
Sterling Harper 2019-04-03 19:50:12 -04:00
parent 5d5b82f0ae
commit d0a0d34b89
3 changed files with 76 additions and 62 deletions
Download patch file Download diff file Expand all files Collapse all files

37
include/openmc/cell.h
View file

@ -210,44 +210,9 @@ public:
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) 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 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 partitions_[middle];
}
}
}
//! A sorted vector of indices to surfaces that partition the universe
std::vector<int32_t> surfs_;

96
src/cell.cpp
View file

@ -659,10 +659,10 @@ UniversePartitioner::UniversePartitioner(const Universe& univ)
}
}
// Define a functor for comparing z-planes by their location, and use it to
// Define a lambda 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
std::sort(surfs_.begin(), surfs_.end(),
[](int32_t i_surf, int32_t j_surf) -> bool
{
const auto* surf = model::surfaces[i_surf].get();
const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf);
@ -672,25 +672,24 @@ UniversePartitioner::UniversePartitioner(const Universe& univ)
double zj = zplane->z0_;
return zi < zj;
}
} compare_surfs;
std::sort(surfs_.begin(), surfs_.end(), compare_surfs);
);
// 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;
int32_t lower_token = 0, upper_token = 0;
double min_z, max_z;
for (auto token : model::cells[i_cell]->rpn_) {
if (token < OP_UNION) {
const auto* surf = model::surfaces[std::abs(token) - 1].get();
if (const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf)) {
if (min_token == 0 || zplane->z0_ < min_z) {
min_token = token;
if (lower_token == 0 || zplane->z0_ < min_z) {
lower_token = token;
min_z = zplane->z0_;
}
if (max_token == 0 || zplane->z0_ > max_z) {
max_token = token;
if (upper_token == 0 || zplane->z0_ > max_z) {
upper_token = token;
max_z = zplane->z0_;
}
}
@ -698,27 +697,78 @@ UniversePartitioner::UniversePartitioner(const Universe& univ)
}
// If there are no bounding z-planes, add this cell to all partitions.
if (min_token == 0) {
if (lower_token == 0) {
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) {
partitions_[i].push_back(i_cell);
if (max_token == -(surfs_[i] + 1)) {
in_partition = false;
// Find the first partition this cell lies in. If the lower_token indicates
// a negative halfspace, then the cell is unbounded in the lower direction
// and it lies in the first partition onward. Otherwise, it is bounded by
// the positive halfspace given by the lower_token.
int first_partition = 0;
if (lower_token > 0) {
for (int i = 0; i < surfs_.size(); ++i) {
if (lower_token == surfs_[i] + 1) {
first_partition = i + 1;
break;
}
} else {
in_partition = min_token == surfs_[i] + 1;
}
}
if (in_partition) partitions_.back().push_back(i_cell);
// Find the last partition this cell lies in. The logic is analogous to the
// logic for first_partition.
int last_partition = surfs_.size();
if (upper_token < 0) {
for (int i = first_partition; i < surfs_.size(); ++i) {
if (upper_token == -(surfs_[i] + 1)) {
last_partition = i;
break;
}
}
}
// Add the cell to all relevant partitions.
for (int i = first_partition; i < last_partition + 1; ++i) {
partitions_[i].push_back(i_cell);
}
}
}
const std::vector<int32_t>&
UniversePartitioner::get_cells(Position r, Direction u) const
{
// Perform a binary search for the partition containing the given coordinates.
int left = 0;
int middle = (surfs_.size() - 1) / 2;
int right = surfs_.size() - 1;
while (true) {
// Check the sense of the coordinates for the current surface.
const auto& surf = *model::surfaces[surfs_[middle]];
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) {
left = middle + 1;
middle = right_leaf;
} else {
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) {
right = middle-1;
middle = left_leaf;
} else {
return partitions_[middle];
}
}
}
}

5
src/geometry_aux.cpp
View file

@ -139,14 +139,13 @@ partition_universes()
}
// Partition the universe if there are more than 5 z-planes. (Fewer than
// five is likely no worth it.)
// 5 is likely not 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));
univ->partitioner_ = std::make_unique<UniversePartitioner>(*univ);
break;
}
}