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

Separating the universe classes into their own source and header file

Browse source
This commit is contained in:
Patrick Shriwise 2022-02-15 15:36:55 -06:00
parent 09c524a3e7
commit 4699732faa
5 changed files with 302 additions and 268 deletions
Download patch file Download diff file Expand all files Collapse all files

1
CMakeLists.txt
View file

@ -351,6 +351,7 @@ list(APPEND libopenmc_SOURCES
src/timer.cpp
src/thermal.cpp
src/track_output.cpp
src/universe.cpp
src/urr.cpp
src/volume_calc.cpp
src/weight_windows.cpp

62
include/openmc/cell.h
View file

@ -17,6 +17,7 @@
#include "openmc/neighbor_list.h"
#include "openmc/position.h"
#include "openmc/surface.h"
#include "openmc/universe.h"
#include "openmc/vector.h"
namespace openmc {
@ -48,36 +49,8 @@ namespace model {
extern std::unordered_map<int32_t, int32_t> cell_map;
extern vector<unique_ptr<Cell>> cells;
extern std::unordered_map<int32_t, int32_t> universe_map;
extern vector<unique_ptr<Universe>> universes;
} // namespace model
//==============================================================================
//! A geometry primitive that fills all space and contains cells.
//==============================================================================
class Universe {
public:
int32_t id_; //!< Unique ID
vector<int32_t> cells_; //!< Cells within this universe
//! \brief Write universe information to an HDF5 group.
//! \param group_id An HDF5 group id.
virtual void to_hdf5(hid_t group_id) const;
virtual bool find_cell(Particle& p) const;
BoundingBox bounding_box() const;
const GeometryType& geom_type() const { return geom_type_; }
GeometryType& geom_type() { return geom_type_; }
unique_ptr<UniversePartitioner> partitioner_;
private:
GeometryType geom_type_ = GeometryType::CSG;
};
//==============================================================================
//==============================================================================
@ -296,35 +269,6 @@ protected:
vector<int32_t>::iterator start, const vector<int32_t>& rpn);
};
//==============================================================================
//! 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:
explicit UniversePartitioner(const Universe& univ);
//! Return the list of cells that could contain the given coordinates.
const vector<int32_t>& get_cells(Position r, Direction u) const;
private:
//! A sorted vector of indices to surfaces that partition the universe
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]`.
vector<vector<int32_t>> partitions_;
};
//==============================================================================
//! Define an instance of a particular cell
//==============================================================================
@ -356,9 +300,5 @@ struct CellInstanceHash {
void read_cells(pugi::xml_node node);
#ifdef DAGMC
class DAGUniverse;
#endif
} // namespace openmc
#endif // OPENMC_CELL_H

79
include/openmc/universe.h Normal file
View file

@ -0,0 +1,79 @@
#ifndef OPENMC_UNIVERSE_H
#define OPENMC_UNIVERSE_H
#include "openmc/cell.h"
namespace openmc {
#ifdef DAGMC
class DAGUniverse;
#endif
class Universe;
class UniversePartitioner;
namespace model {
extern std::unordered_map<int32_t, int32_t> universe_map;
extern vector<unique_ptr<Universe>> universes;
} // namespace model
//==============================================================================
//! A geometry primitive that fills all space and contains cells.
//==============================================================================
class Universe {
public:
int32_t id_; //!< Unique ID
vector<int32_t> cells_; //!< Cells within this universe
//! \brief Write universe information to an HDF5 group.
//! \param group_id An HDF5 group id.
virtual void to_hdf5(hid_t group_id) const;
virtual bool find_cell(Particle& p) const;
BoundingBox bounding_box() const;
const GeometryType& geom_type() const { return geom_type_; }
GeometryType& geom_type() { return geom_type_; }
unique_ptr<UniversePartitioner> partitioner_;
private:
GeometryType geom_type_ = GeometryType::CSG;
};
//==============================================================================
//! 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:
explicit UniversePartitioner(const Universe& univ);
//! Return the list of cells that could contain the given coordinates.
const vector<int32_t>& get_cells(Position r, Direction u) const;
private:
//! A sorted vector of indices to surfaces that partition the universe
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]`.
vector<vector<int32_t>> partitions_;
};
} // namespace openmc
#endif // OPENMC_UNIVERSE_H

207
src/cell.cpp
View file

@ -34,8 +34,6 @@ namespace model {
std::unordered_map<int32_t, int32_t> cell_map;
vector<unique_ptr<Cell>> cells;
std::unordered_map<int32_t, int32_t> universe_map;
vector<unique_ptr<Universe>> universes;
} // namespace model
//==============================================================================
@ -186,66 +184,6 @@ vector<int32_t> generate_rpn(int32_t cell_id, vector<int32_t> infix)
return rpn;
}
//==============================================================================
// Universe implementation
//==============================================================================
void Universe::to_hdf5(hid_t universes_group) const
{
// Create a group for this universe.
auto group = create_group(universes_group, fmt::format("universe {}", id_));
// Write the geometry representation type.
write_string(group, "geom_type", "csg", false);
// Write the contained cells.
if (cells_.size() > 0) {
vector<int32_t> cell_ids;
for (auto i_cell : cells_)
cell_ids.push_back(model::cells[i_cell]->id_);
write_dataset(group, "cells", cell_ids);
}
close_group(group);
}
bool Universe::find_cell(Particle& p) const
{
const auto& cells {
!partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())};
for (auto it = cells.begin(); it != cells.end(); it++) {
int32_t i_cell = *it;
int32_t i_univ = p.coord(p.n_coord() - 1).universe;
if (model::cells[i_cell]->universe_ != i_univ)
continue;
// Check if this cell contains the particle;
Position r {p.r_local()};
Direction u {p.u_local()};
auto surf = p.surface();
if (model::cells[i_cell]->contains(r, u, surf)) {
p.coord(p.n_coord() - 1).cell = i_cell;
return true;
}
}
return false;
}
BoundingBox Universe::bounding_box() const
{
BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY};
if (cells_.size() == 0) {
return {};
} else {
for (const auto& cell : cells_) {
auto& c = model::cells[cell];
bbox |= c->bounding_box();
}
}
return bbox;
}
//==============================================================================
// Cell implementation
//==============================================================================
@ -876,151 +814,6 @@ bool CSGCell::contains_complex(
}
}
//==============================================================================
// UniversePartitioner implementation
//==============================================================================
UniversePartitioner::UniversePartitioner(const Universe& univ)
{
// Define an ordered set of surface indices that point to z-planes. Use a
// functor to to order the set by the z0_ values of the corresponding planes.
struct compare_surfs {
bool operator()(const int32_t& i_surf, const int32_t& j_surf) const
{
const auto* surf = model::surfaces[i_surf].get();
const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf);
double zi = zplane->z0_;
surf = model::surfaces[j_surf].get();
zplane = dynamic_cast<const SurfaceZPlane*>(surf);
double zj = zplane->z0_;
return zi < zj;
}
};
std::set<int32_t, compare_surfs> surf_set;
// Find all of the z-planes in this universe. A set is used here for the
// O(log(n)) insertions that will ensure entries are not repeated.
for (auto i_cell : univ.cells_) {
for (auto token : model::cells[i_cell]->rpn_) {
if (token < OP_UNION) {
auto i_surf = std::abs(token) - 1;
const auto* surf = model::surfaces[i_surf].get();
if (const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf))
surf_set.insert(i_surf);
}
}
}
// Populate the surfs_ vector from the ordered set.
surfs_.insert(surfs_.begin(), surf_set.begin(), surf_set.end());
// Populate the partition lists.
partitions_.resize(surfs_.size() + 1);
for (auto i_cell : univ.cells_) {
// It is difficult to determine the bounds of a complex cell, so add complex
// cells to all partitions.
if (!model::cells[i_cell]->simple_) {
for (auto& p : partitions_)
p.push_back(i_cell);
continue;
}
// Find the tokens for bounding z-planes.
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 (lower_token == 0 || zplane->z0_ < min_z) {
lower_token = token;
min_z = zplane->z0_;
}
if (upper_token == 0 || zplane->z0_ > max_z) {
upper_token = token;
max_z = zplane->z0_;
}
}
}
}
// If there are no bounding z-planes, add this cell to all partitions.
if (lower_token == 0) {
for (auto& p : partitions_)
p.push_back(i_cell);
continue;
}
// 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;
}
}
}
// 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; ++i) {
partitions_[i].push_back(i_cell);
}
}
}
const 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];
}
}
}
}
//==============================================================================
// Non-method functions
//==============================================================================

221
src/universe.cpp Normal file
View file

@ -0,0 +1,221 @@
#include "openmc/universe.h"
#include <set>
#include "openmc/hdf5_interface.h"
namespace openmc {
namespace model {
std::unordered_map<int32_t, int32_t> universe_map;
vector<unique_ptr<Universe>> universes;
} // namespace model
//==============================================================================
// Universe implementation
//==============================================================================
void Universe::to_hdf5(hid_t universes_group) const
{
// Create a group for this universe.
auto group = create_group(universes_group, fmt::format("universe {}", id_));
// Write the geometry representation type.
write_string(group, "geom_type", "csg", false);
// Write the contained cells.
if (cells_.size() > 0) {
vector<int32_t> cell_ids;
for (auto i_cell : cells_)
cell_ids.push_back(model::cells[i_cell]->id_);
write_dataset(group, "cells", cell_ids);
}
close_group(group);
}
bool Universe::find_cell(Particle& p) const
{
const auto& cells {
!partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())};
for (auto it = cells.begin(); it != cells.end(); it++) {
int32_t i_cell = *it;
int32_t i_univ = p.coord(p.n_coord() - 1).universe;
if (model::cells[i_cell]->universe_ != i_univ)
continue;
// Check if this cell contains the particle;
Position r {p.r_local()};
Direction u {p.u_local()};
auto surf = p.surface();
if (model::cells[i_cell]->contains(r, u, surf)) {
p.coord(p.n_coord() - 1).cell = i_cell;
return true;
}
}
return false;
}
BoundingBox Universe::bounding_box() const
{
BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY};
if (cells_.size() == 0) {
return {};
} else {
for (const auto& cell : cells_) {
auto& c = model::cells[cell];
bbox |= c->bounding_box();
}
}
return bbox;
}
//==============================================================================
// UniversePartitioner implementation
//==============================================================================
UniversePartitioner::UniversePartitioner(const Universe& univ)
{
// Define an ordered set of surface indices that point to z-planes. Use a
// functor to to order the set by the z0_ values of the corresponding planes.
struct compare_surfs {
bool operator()(const int32_t& i_surf, const int32_t& j_surf) const
{
const auto* surf = model::surfaces[i_surf].get();
const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf);
double zi = zplane->z0_;
surf = model::surfaces[j_surf].get();
zplane = dynamic_cast<const SurfaceZPlane*>(surf);
double zj = zplane->z0_;
return zi < zj;
}
};
std::set<int32_t, compare_surfs> surf_set;
// Find all of the z-planes in this universe. A set is used here for the
// O(log(n)) insertions that will ensure entries are not repeated.
for (auto i_cell : univ.cells_) {
for (auto token : model::cells[i_cell]->rpn_) {
if (token < OP_UNION) {
auto i_surf = std::abs(token) - 1;
const auto* surf = model::surfaces[i_surf].get();
if (const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf))
surf_set.insert(i_surf);
}
}
}
// Populate the surfs_ vector from the ordered set.
surfs_.insert(surfs_.begin(), surf_set.begin(), surf_set.end());
// Populate the partition lists.
partitions_.resize(surfs_.size() + 1);
for (auto i_cell : univ.cells_) {
// It is difficult to determine the bounds of a complex cell, so add complex
// cells to all partitions.
if (!model::cells[i_cell]->simple_) {
for (auto& p : partitions_)
p.push_back(i_cell);
continue;
}
// Find the tokens for bounding z-planes.
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 (lower_token == 0 || zplane->z0_ < min_z) {
lower_token = token;
min_z = zplane->z0_;
}
if (upper_token == 0 || zplane->z0_ > max_z) {
upper_token = token;
max_z = zplane->z0_;
}
}
}
}
// If there are no bounding z-planes, add this cell to all partitions.
if (lower_token == 0) {
for (auto& p : partitions_)
p.push_back(i_cell);
continue;
}
// 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;
}
}
}
// 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; ++i) {
partitions_[i].push_back(i_cell);
}
}
}
const 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];
}
}
}
}
}