Refactor cell and particle data classes to adapt the virtual lattice function

This commit is contained in:
skywalker_cn 2025-09-04 10:46:15 +00:00
parent 309841ead1
commit d1ee645b22
8 changed files with 215 additions and 90 deletions

View file

@ -99,13 +99,13 @@ public:
//! Get Boolean of if the cell is simple or not
bool is_simple() const { return simple_; }
//! Get a vector of the region expression in postfix notation
vector<int32_t> generate_postfix(int32_t cell_id) const;
private:
//----------------------------------------------------------------------------
// Private Methods
//! Get a vector of the region expression in postfix notation
vector<int32_t> generate_postfix(int32_t cell_id) const;
//! Determine if a particle is inside the cell for a simple cell (only
//! intersection operators)
bool contains_simple(Position r, Direction u, int32_t on_surface) const;
@ -314,15 +314,14 @@ public:
//----------------------------------------------------------------------------
// Data members
int32_t id_; //!< Unique ID
std::string name_; //!< User-defined name
Fill type_; //!< Material, universe, or lattice
int32_t universe_; //!< Universe # this cell is in
int32_t fill_; //!< Universe # filling this cell
bool virtual_lattice_; //!< If the cell is the base of a virtual triso lattice
int32_t id_; //!< Unique ID
std::string name_; //!< User-defined name
Fill type_; //!< Material, universe, or lattice
int32_t universe_; //!< Universe # this cell is in
int32_t fill_; //!< Universe # filling this cell
bool virtual_lattice_; //!< If the cell is the base of a virtual triso lattice
bool triso_particle_;
//! \brief Specification of the virtual lattice
vector<double> vl_lower_left_;
vector<double> vl_pitch_;
@ -380,10 +379,7 @@ public:
vector<int32_t> surfaces() const override { return region_.surfaces(); }
std::pair<double, int32_t> distance(Position r, Direction u,
int32_t on_surface, GeometryState* p) const override
{
return region_.distance(r, u, on_surface);
}
int32_t on_surface, GeometryState* p) const override;
bool contains(Position r, Direction u, int32_t on_surface) const override
{

View file

@ -372,22 +372,52 @@ public:
// Boundary information
BoundaryInfo& boundary() { return boundary_; }
// Distance to the next collision
double& collision_distance() { return collision_distance_; }
#ifdef OPENMC_DAGMC_ENABLED
// DagMC state variables
moab::DagMC::RayHistory& history() { return history_; }
Direction& last_dir() { return last_dir_; }
moab::DagMC::RayHistory& history()
{
return history_;
}
Direction& last_dir()
{
return last_dir_;
}
#endif
// material of current and last cell
int& material() { return material_; }
const int& material() const { return material_; }
int& material_last() { return material_last_; }
const int& material_last() const { return material_last_; }
int& material()
{
return material_;
}
const int& material() const
{
return material_;
}
int& material_last()
{
return material_last_;
}
const int& material_last() const
{
return material_last_;
}
// temperature of current and last cell
double& sqrtkT() { return sqrtkT_; }
const double& sqrtkT() const { return sqrtkT_; }
double& sqrtkT_last() { return sqrtkT_last_; }
double& sqrtkT()
{
return sqrtkT_;
}
const double& sqrtkT() const
{
return sqrtkT_;
}
double& sqrtkT_last()
{
return sqrtkT_last_;
}
private:
int64_t id_ {-1}; //!< Unique ID
@ -417,6 +447,8 @@ private:
double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV
double sqrtkT_last_ {0.0}; //!< last temperature
double collision_distance_ {INFTY};
#ifdef OPENMC_DAGMC_ENABLED
moab::DagMC::RayHistory history_;
Direction last_dir_;
@ -528,8 +560,6 @@ private:
bool trace_ {false};
double collision_distance_;
int n_event_ {0};
int n_split_ {0};
@ -695,9 +725,6 @@ public:
// Shows debug info
bool& trace() { return trace_; }
// Distance to the next collision
double& collision_distance() { return collision_distance_; }
// Number of events particle has undergone
int& n_event() { return n_event_; }

View file

@ -274,6 +274,9 @@ public:
bool triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const override;
BoundingBox bounding_box(bool pos_side) const override;
vector<double> get_center() const override;
double get_radius() const override;
void connect_to_triso_base(int triso_index, std::string key) override;
double x0_, y0_, z0_, radius_;
// int triso_base_index_ = -1;

View file

@ -29,6 +29,7 @@ class Universe {
public:
int32_t id_; //!< Unique ID
vector<int32_t> cells_; //!< Cells within this universe
int filled_with_triso_base_ = -1; //!< ID of cell filled with virtual lattice
int32_t n_instances_; //!< Number of instances of this universe
//! \brief Write universe information to an HDF5 group.

View file

@ -422,18 +422,19 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
// Morgans law
Region region(region_spec, id_);
region_ = region;
vector<int32_t> rpn = region_.generate_postfix(id_);
if (virtual_lattice_) {
vl_triso_distribution_ = generate_triso_distribution(
vl_shape_, vl_pitch_, vl_lower_left_, rpn_, id_);
vl_shape_, vl_pitch_, vl_lower_left_, rpn, id_);
}
if (triso_particle_) {
if (rpn_.size() != 1) {
if (rpn.size() != 1) {
fatal_error(
fmt::format("Wrong surface definition of triso particle cell {}", id_));
} else {
model::surfaces[abs(rpn_[0]) - 1]->connect_to_triso_base(id_, "particle");
model::surfaces[abs(rpn[0]) - 1]->connect_to_triso_base(id_, "particle");
}
}
@ -591,7 +592,7 @@ std::pair<double, int32_t> CSGCell::distance(
}
} else {
region_.distance(r, u, on_surface);
return region_.distance(r, u, on_surface);
}
return {min_dist, i_surf};

View file

@ -225,8 +225,6 @@ void Particle::event_calculate_xs()
void Particle::event_advance()
{
// Find the distance to the nearest boundary
// Sample a distance to collision
if (type() == ParticleType::electron || type() == ParticleType::positron) {
collision_distance() = 0.0;
@ -236,6 +234,7 @@ void Particle::event_advance()
collision_distance() = -std::log(prn(current_seed())) / macro_xs().total;
}
// Find the distance to the nearest boundary
boundary() = distance_to_boundary(*this);
double speed = this->speed();
@ -575,14 +574,14 @@ void Particle::cross_surface(const Surface& surf)
return;
}
#endif
int i_surface = std::abs(surface());
bool verbose = settings::verbosity >= 10 || trace();
if (surf->is_triso_surface_) {
if (surf.is_triso_surface_) {
if (surface() > 0) {
for (int i = n_coord(); i < model::n_coord_levels; i++) {
coord(i).reset();
}
coord(n_coord() - 1).cell =
coord(n_coord() - 1).cell() =
model::cell_map[model::surfaces[i_surface - 1]->triso_base_index_];
} else if (surface() < 0) {
for (int i = n_coord(); i < model::n_coord_levels; i++) {
@ -592,20 +591,20 @@ void Particle::cross_surface(const Surface& surf)
fatal_error(fmt::format("Particle cell of surface {} is not defined",
model::surfaces[i_surface - 1]->id_));
}
coord(n_coord() - 1).cell =
coord(n_coord() - 1).cell() =
model::cell_map[model::surfaces[i_surface - 1]->triso_particle_index_];
}
// find material
bool found = true;
int i_cell = coord(n_coord() - 1).cell;
int i_cell = coord(n_coord() - 1).cell();
for (;; ++n_coord()) {
if (i_cell == C_NONE) {
int i_universe = coord(n_coord() - 1).universe;
int i_universe = coord(n_coord() - 1).universe();
const auto& univ {model::universes[i_universe]};
if (univ->filled_with_triso_base_ != -1) {
coord(n_coord() - 1).cell =
coord(n_coord() - 1).cell() =
model::cell_map[univ->filled_with_triso_base_];
found = true;
} else {
@ -616,7 +615,7 @@ void Particle::cross_surface(const Surface& surf)
}
}
i_cell = coord(n_coord() - 1).cell;
i_cell = coord(n_coord() - 1).cell();
Cell& c {*model::cells[i_cell]};
if (c.type_ == Fill::MATERIAL) {
@ -650,14 +649,14 @@ void Particle::cross_surface(const Surface& surf)
// Set the lower coordinate level universe.
auto& coor {coord(n_coord())};
coor.universe = c.fill_;
coor.universe() = c.fill_;
// Set the position and direction.
coor.r = r_local();
coor.u = u_local();
coor.r() = r_local();
coor.u() = u_local();
// Apply translation.
coor.r -= c.translation_;
coor.r() -= c.translation_;
// Apply rotation.
if (!c.rotation_.empty()) {

View file

@ -257,7 +257,8 @@ BoundingBox SurfaceXPlane::bounding_box(bool pos_side) const
}
}
bool SurfaceXPlane::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
bool SurfaceXPlane::triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const
{
return false;
}
@ -265,7 +266,11 @@ bool SurfaceXPlane::triso_in_mesh(vector<double> mesh_center, vector<double> lat
//==============================================================================
// SurfaceYPlane implementation
//==============================================================================
bool SurfaceYPlane::triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const
{
return false;
}
SurfaceYPlane::SurfaceYPlane(pugi::xml_node surf_node) : Surface(surf_node)
{
read_coeffs(surf_node, id_, {&y0_});
@ -305,7 +310,11 @@ BoundingBox SurfaceYPlane::bounding_box(bool pos_side) const
//==============================================================================
// SurfaceZPlane implementation
//==============================================================================
bool SurfaceZPlane::triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const
{
return false;
}
SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node) : Surface(surf_node)
{
read_coeffs(surf_node, id_, {&z0_});
@ -345,7 +354,8 @@ BoundingBox SurfaceZPlane::bounding_box(bool pos_side) const
//==============================================================================
// SurfacePlane implementation
//==============================================================================
bool SurfacePlane::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
bool SurfacePlane::triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const
{
return false;
}
@ -468,7 +478,8 @@ Direction axis_aligned_cylinder_normal(
// SurfaceXCylinder implementation
//==============================================================================
bool SurfaceXCylinder::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
bool SurfaceXCylinder::triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const
{
return false;
}
@ -516,7 +527,8 @@ BoundingBox SurfaceXCylinder::bounding_box(bool pos_side) const
// SurfaceYCylinder implementation
//==============================================================================
bool SurfaceYCylinder::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
bool SurfaceYCylinder::triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const
{
return false;
}
@ -565,7 +577,8 @@ BoundingBox SurfaceYCylinder::bounding_box(bool pos_side) const
// SurfaceZCylinder implementation
//==============================================================================
bool SurfaceZCylinder::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
bool SurfaceZCylinder::triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const
{
return false;
}
@ -690,58 +703,61 @@ BoundingBox SurfaceSphere::bounding_box(bool pos_side) const
void SurfaceSphere::connect_to_triso_base(int triso_index, std::string key)
{
if (key=="base") {
triso_base_index_=triso_index;
is_triso_surface_=true;
} else if (key=="particle") {
triso_particle_index_=triso_index;
if (key == "base") {
triso_base_index_ = triso_index;
is_triso_surface_ = true;
} else if (key == "particle") {
triso_particle_index_ = triso_index;
}
}
vector<double> SurfaceSphere::get_center() const {
return {x0_,y0_,z0_};
vector<double> SurfaceSphere::get_center() const
{
return {x0_, y0_, z0_};
}
double SurfaceSphere::get_radius() const {
double SurfaceSphere::get_radius() const
{
return radius_;
}
bool SurfaceSphere::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
bool SurfaceSphere::triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const
{
double dis_x;
double dis_y;
double dis_z;
double x_min=mesh_center[0]-lattice_pitch[0]/2;
double x_max=mesh_center[0]+lattice_pitch[0]/2;
double y_min=mesh_center[1]-lattice_pitch[1]/2;
double y_max=mesh_center[1]+lattice_pitch[1]/2;
double z_min=mesh_center[2]-lattice_pitch[2]/2;
double z_max=mesh_center[2]+lattice_pitch[2]/2;
if (x0_>=x_min && x0_<=x_max) {
dis_x=0;
} else if (x0_<x_min) {
dis_x=pow(x_min-x0_, 2);
double x_min = mesh_center[0] - lattice_pitch[0] / 2;
double x_max = mesh_center[0] + lattice_pitch[0] / 2;
double y_min = mesh_center[1] - lattice_pitch[1] / 2;
double y_max = mesh_center[1] + lattice_pitch[1] / 2;
double z_min = mesh_center[2] - lattice_pitch[2] / 2;
double z_max = mesh_center[2] + lattice_pitch[2] / 2;
if (x0_ >= x_min && x0_ <= x_max) {
dis_x = 0;
} else if (x0_ < x_min) {
dis_x = pow(x_min - x0_, 2);
} else {
dis_x=pow(x_max-x0_, 2);
dis_x = pow(x_max - x0_, 2);
}
if (y0_>=y_min && y0_<=y_max) {
dis_y=0;
} else if (y0_<y_min) {
dis_y=pow(y_min-y0_, 2);
if (y0_ >= y_min && y0_ <= y_max) {
dis_y = 0;
} else if (y0_ < y_min) {
dis_y = pow(y_min - y0_, 2);
} else {
dis_y=pow(y_max-y0_, 2);
dis_y = pow(y_max - y0_, 2);
}
if (z0_>=z_min && z0_<=z_max) {
dis_z=0;
} else if (z0_<z_min) {
dis_z=pow(z_min-z0_, 2);
if (z0_ >= z_min && z0_ <= z_max) {
dis_z = 0;
} else if (z0_ < z_min) {
dis_z = pow(z_min - z0_, 2);
} else {
dis_z=pow(z_max-z0_, 2);
dis_z = pow(z_max - z0_, 2);
}
if (sqrt(dis_x+dis_y+dis_z) < radius_) {
if (sqrt(dis_x + dis_y + dis_z) < radius_) {
return true;
} else {
return false;
@ -839,7 +855,11 @@ Direction axis_aligned_cone_normal(
//==============================================================================
// SurfaceXCone implementation
//==============================================================================
bool SurfaceXCone::triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const
{
return false;
}
SurfaceXCone::SurfaceXCone(pugi::xml_node surf_node) : Surface(surf_node)
{
read_coeffs(surf_node, id_, {&x0_, &y0_, &z0_, &radius_sq_});
@ -871,7 +891,11 @@ void SurfaceXCone::to_hdf5_inner(hid_t group_id) const
//==============================================================================
// SurfaceYCone implementation
//==============================================================================
bool SurfaceYCone::triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const
{
return false;
}
SurfaceYCone::SurfaceYCone(pugi::xml_node surf_node) : Surface(surf_node)
{
read_coeffs(surf_node, id_, {&x0_, &y0_, &z0_, &radius_sq_});
@ -903,6 +927,11 @@ void SurfaceYCone::to_hdf5_inner(hid_t group_id) const
//==============================================================================
// SurfaceZCone implementation
//==============================================================================
bool SurfaceZCone::triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const
{
return false;
}
SurfaceZCone::SurfaceZCone(pugi::xml_node surf_node) : Surface(surf_node)
{
@ -935,7 +964,11 @@ void SurfaceZCone::to_hdf5_inner(hid_t group_id) const
//==============================================================================
// SurfaceQuadric implementation
//==============================================================================
bool SurfaceQuadric::triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const
{
return false;
}
SurfaceQuadric::SurfaceQuadric(pugi::xml_node surf_node) : Surface(surf_node)
{
read_coeffs(
@ -1095,7 +1128,11 @@ double torus_distance(double x1, double x2, double x3, double u1, double u2,
//==============================================================================
// SurfaceXTorus implementation
//==============================================================================
bool SurfaceXTorus::triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const
{
return false;
}
SurfaceXTorus::SurfaceXTorus(pugi::xml_node surf_node) : Surface(surf_node)
{
read_coeffs(surf_node, id_, {&x0_, &y0_, &z0_, &A_, &B_, &C_});
@ -1148,7 +1185,11 @@ Direction SurfaceXTorus::normal(Position r) const
//==============================================================================
// SurfaceYTorus implementation
//==============================================================================
bool SurfaceYTorus::triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const
{
return false;
}
SurfaceYTorus::SurfaceYTorus(pugi::xml_node surf_node) : Surface(surf_node)
{
read_coeffs(surf_node, id_, {&x0_, &y0_, &z0_, &A_, &B_, &C_});
@ -1201,7 +1242,8 @@ Direction SurfaceYTorus::normal(Position r) const
//==============================================================================
// SurfaceZTorus implementation
//==============================================================================
bool SurfaceZTorus::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
bool SurfaceZTorus::triso_in_mesh(
vector<double> mesh_center, vector<double> lattice_pitch) const
{
return false;
}

View file

@ -39,6 +39,62 @@ void Universe::to_hdf5(hid_t universes_group) const
bool Universe::find_cell(GeometryState& p) const
{
if (filled_with_triso_base_ != -1) {
Cell& c {*model::cells[model::cell_map[filled_with_triso_base_]]};
vector<int> lat_ind(3);
Position r {p.r_local()};
lat_ind[0] = std::max<int>(
std::min<int>(floor((r.x - c.vl_lower_left_[0]) / c.vl_pitch_[0]),
c.vl_shape_[0] - 1),
0);
lat_ind[1] = std::max<int>(
std::min<int>(floor((r.y - c.vl_lower_left_[1]) / c.vl_pitch_[1]),
c.vl_shape_[1] - 1),
0);
lat_ind[2] = std::max<int>(
std::min<int>(floor((r.z - c.vl_lower_left_[2]) / c.vl_pitch_[2]),
c.vl_shape_[2] - 1),
0);
int32_t i_univ = p.coord(p.n_coord() - 1).universe();
for (int token :
c.vl_triso_distribution_[lat_ind[0] + lat_ind[1] * c.vl_shape_[0] +
lat_ind[2] * c.vl_shape_[0] * c.vl_shape_[1]]) {
vector<double> triso_center =
model::surfaces[abs(token) - 1]->get_center();
double triso_radius = model::surfaces[abs(token) - 1]->get_radius();
if (model::cells[model::cell_map[model::surfaces[abs(token) - 1]
->triso_base_index_]]
->universe_ != i_univ)
continue;
if (abs(token) == abs(p.surface())) {
if (p.surface() < 0) {
p.coord(p.n_coord() - 1).cell() =
model::cell_map[model::surfaces[abs(token) - 1]
->triso_particle_index_];
return true;
} else {
p.coord(p.n_coord() - 1).cell() =
model::cell_map[filled_with_triso_base_];
return true;
}
}
if (pow(r.x - triso_center[0], 2) + pow(r.y - triso_center[1], 2) +
pow(r.z - triso_center[2], 2) <
pow(triso_radius, 2)) {
p.coord(p.n_coord() - 1).cell() =
model::cell_map[model::surfaces[abs(token) - 1]
->triso_particle_index_];
return true;
}
}
if (model::cells[model::cell_map[filled_with_triso_base_]]->universe_ ==
i_univ) {
p.coord(p.n_coord() - 1).cell() =
model::cell_map[filled_with_triso_base_];
return true;
}
}
const auto& cells {
!partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())};