From 288c0f25319d00679ba3dae3fdc92e9d95e656f2 Mon Sep 17 00:00:00 2001 From: skywalker_cn <2324714651@qq.com> Date: Tue, 9 Sep 2025 13:58:55 +0000 Subject: [PATCH] Adds support for virtual lattice geometry handling Introduces functionality to calculate distances and locate cells within virtual lattice geometries. Refactors existing geometry handling logic to reduce code duplication and improve maintainability. Key updates include: - Addition of `distance_in_virtual_lattice` and `find_cell_in_virtual_lattice` methods for specialized handling. - Integration of virtual lattice checks in standard cell distance computations. - Refactoring of repetitive code in surface crossing handlers for improved clarity. Enhances support for TRISO particle geometries within virtual lattices, ensuring accurate material and cell identification during particle transport. --- include/openmc/cell.h | 3 + include/openmc/geometry.h | 4 + include/openmc/universe.h | 2 + src/cell.cpp | 167 ++++++++++++++++++++------------------ src/geometry.cpp | 92 +++++++++++++++++++++ src/particle.cpp | 92 +-------------------- src/universe.cpp | 110 +++++++++++++------------ 7 files changed, 248 insertions(+), 222 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 3008948168..d2365c84c4 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -381,6 +381,9 @@ public: std::pair distance(Position r, Direction u, int32_t on_surface, GeometryState* p) const override; + std::pair distance_in_virtual_lattice( + Position r, Direction u, int32_t on_surface, GeometryState* p) const; + bool contains(Position r, Direction u, int32_t on_surface) const override { return region_.contains(r, u, on_surface); diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index 107cc7d1f3..d41a98e301 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -64,6 +64,10 @@ bool exhaustive_find_cell(GeometryState& p, bool verbose = false); bool neighbor_list_find_cell( GeometryState& p, bool verbose = false); // Only usable on surface crossings +bool find_cell_in_virtual_lattice(GeometryState& p, + bool verbose = + false); // Only usable on triso surface crossings in virtual lattice + //============================================================================== //! Move a particle into a new lattice tile. //============================================================================== diff --git a/include/openmc/universe.h b/include/openmc/universe.h index 2afc4b97ac..1c9fcceabf 100644 --- a/include/openmc/universe.h +++ b/include/openmc/universe.h @@ -38,6 +38,8 @@ public: virtual bool find_cell(GeometryState& p) const; + virtual bool find_cell_in_virtual_lattice(GeometryState& p) const; + BoundingBox bounding_box() const; /* By default, universes are CSG universes. The DAGMC diff --git a/src/cell.cpp b/src/cell.cpp index a5d68c8c56..f7d077fb51 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -505,96 +505,101 @@ vector::iterator CSGCell::find_left_parenthesis( } std::pair CSGCell::distance( Position r, Direction u, int32_t on_surface, GeometryState* p) const +{ + if (virtual_lattice_) { + return distance_in_virtual_lattice(r, u, on_surface, p); + } else { + return region_.distance(r, u, on_surface); + } +} + +std::pair CSGCell::distance_in_virtual_lattice( + Position r, Direction u, int32_t on_surface, GeometryState* p) const { double min_dist {INFTY}; int32_t i_surf {std::numeric_limits::max()}; double min_dis_vl; int32_t i_surf_vl; - if (virtual_lattice_) { - double max_dis = p->collision_distance(); - double tol_dis = 0; - vector dis_to_bou(3), dis_to_bou_max(3); - double u_value = sqrt(pow(u.x, 2) + pow(u.y, 2) + - pow(u.z, 2)); // don't know if u has been normalized - vector norm_u = {u.x / u_value, u.y / u_value, u.z / u_value}; - vector lat_ind(3); - vector temp_pos = {r.x, r.y, r.z}; - int loop_time; - for (int i = 0; i < 3; i++) { - lat_ind[i] = floor((temp_pos[i] - vl_lower_left_[i]) / vl_pitch_[i]); - if (lat_ind[i] == vl_shape_[i] && norm_u[i] < 0) { - lat_ind[i] = vl_shape_[i] - 1; - } - if (lat_ind[i] == -1 && norm_u[i] > 0) { - lat_ind[i] = 0; - } + + double max_dis = p->collision_distance(); + double tol_dis = 0; + vector dis_to_bou(3), dis_to_bou_max(3); + double u_value = sqrt(pow(u.x, 2) + pow(u.y, 2) + + pow(u.z, 2)); // don't know if u has been normalized + vector norm_u = {u.x / u_value, u.y / u_value, u.z / u_value}; + vector lat_ind(3); + vector temp_pos = {r.x, r.y, r.z}; + int loop_time; + for (int i = 0; i < 3; i++) { + lat_ind[i] = floor((temp_pos[i] - vl_lower_left_[i]) / vl_pitch_[i]); + if (lat_ind[i] == vl_shape_[i] && norm_u[i] < 0) { + lat_ind[i] = vl_shape_[i] - 1; } - - dis_to_bou = {INFTY, INFTY, INFTY}; - for (int i = 0; i < 3; i++) { - if (norm_u[i] > 0) { - dis_to_bou[i] = std::abs( - ((lat_ind[i] + 1) * vl_pitch_[i] + vl_lower_left_[i] - temp_pos[i]) / - norm_u[i]); - dis_to_bou_max[i] = vl_pitch_[i] / norm_u[i]; - } else if (norm_u[i] < 0) { - dis_to_bou[i] = std::abs( - (lat_ind[i] * vl_pitch_[i] + vl_lower_left_[i] - temp_pos[i]) / - norm_u[i]); - dis_to_bou_max[i] = -vl_pitch_[i] / norm_u[i]; - } + if (lat_ind[i] == -1 && norm_u[i] > 0) { + lat_ind[i] = 0; } - - while (true) { - if (lat_ind[0] < 0 || lat_ind[0] >= vl_shape_[0] || lat_ind[1] < 0 || - lat_ind[1] >= vl_shape_[1] || lat_ind[2] < 0 || - lat_ind[2] >= vl_shape_[2]) - break; - - for (int token : - vl_triso_distribution_[lat_ind[0] + lat_ind[1] * vl_shape_[0] + - lat_ind[2] * vl_shape_[0] * vl_shape_[1]]) { - bool coincident {std::abs(token) == std::abs(on_surface)}; - double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)}; - if (d < min_dist) { - if (min_dist - d >= FP_PRECISION * min_dist) { - min_dist = d; - i_surf = -token; - } - } - } - - int mes_bou_crossed = 0; - if (dis_to_bou[1] < dis_to_bou[0]) { - mes_bou_crossed = 1; - } - if (dis_to_bou[2] < dis_to_bou[mes_bou_crossed]) { - mes_bou_crossed = 2; - } - - tol_dis = dis_to_bou[mes_bou_crossed]; - - if (min_dist < tol_dis) { - break; - } - - if (norm_u[mes_bou_crossed] > 0) { - lat_ind[mes_bou_crossed] += 1; - } else { - lat_ind[mes_bou_crossed] += -1; - } - - dis_to_bou[mes_bou_crossed] += dis_to_bou_max[mes_bou_crossed]; - - if (tol_dis > max_dis) { - break; - } - } - - } else { - return region_.distance(r, u, on_surface); } + dis_to_bou = {INFTY, INFTY, INFTY}; + for (int i = 0; i < 3; i++) { + if (norm_u[i] > 0) { + dis_to_bou[i] = std::abs( + ((lat_ind[i] + 1) * vl_pitch_[i] + vl_lower_left_[i] - temp_pos[i]) / + norm_u[i]); + dis_to_bou_max[i] = vl_pitch_[i] / norm_u[i]; + } else if (norm_u[i] < 0) { + dis_to_bou[i] = + std::abs((lat_ind[i] * vl_pitch_[i] + vl_lower_left_[i] - temp_pos[i]) / + norm_u[i]); + dis_to_bou_max[i] = -vl_pitch_[i] / norm_u[i]; + } + } + + while (true) { + if (lat_ind[0] < 0 || lat_ind[0] >= vl_shape_[0] || lat_ind[1] < 0 || + lat_ind[1] >= vl_shape_[1] || lat_ind[2] < 0 || + lat_ind[2] >= vl_shape_[2]) + break; + + for (int token : + vl_triso_distribution_[lat_ind[0] + lat_ind[1] * vl_shape_[0] + + lat_ind[2] * vl_shape_[0] * vl_shape_[1]]) { + bool coincident {std::abs(token) == std::abs(on_surface)}; + double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)}; + if (d < min_dist) { + if (min_dist - d >= FP_PRECISION * min_dist) { + min_dist = d; + i_surf = -token; + } + } + } + + int mes_bou_crossed = 0; + if (dis_to_bou[1] < dis_to_bou[0]) { + mes_bou_crossed = 1; + } + if (dis_to_bou[2] < dis_to_bou[mes_bou_crossed]) { + mes_bou_crossed = 2; + } + + tol_dis = dis_to_bou[mes_bou_crossed]; + + if (min_dist < tol_dis) { + break; + } + + if (norm_u[mes_bou_crossed] > 0) { + lat_ind[mes_bou_crossed] += 1; + } else { + lat_ind[mes_bou_crossed] += -1; + } + + dis_to_bou[mes_bou_crossed] += dis_to_bou_max[mes_bou_crossed]; + + if (tol_dis > max_dis) { + break; + } + } return {min_dist, i_surf}; } diff --git a/src/geometry.cpp b/src/geometry.cpp index df48508916..f32b63ed7a 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -294,6 +294,98 @@ bool exhaustive_find_cell(GeometryState& p, bool verbose) return find_cell_inner(p, nullptr, verbose); } +bool find_cell_in_virtual_lattice(GeometryState& p, bool verbose) +{ + int i_surface = std::abs(p.surface()); + if (p.surface() > 0) { + for (int i = p.n_coord(); i < model::n_coord_levels; i++) { + p.coord(i).reset(); + } + p.coord(p.n_coord() - 1).cell() = + model::cell_map[model::surfaces[i_surface - 1]->triso_base_index_]; + } else if (p.surface() < 0) { + for (int i = p.n_coord(); i < model::n_coord_levels; i++) { + p.coord(i).reset(); + } + if (model::surfaces[i_surface - 1]->triso_particle_index_ == -1) { + fatal_error(fmt::format("Particle cell of surface {} is not defined", + model::surfaces[i_surface - 1]->id_)); + } + p.lowest_coord().cell() = + model::cell_map[model::surfaces[i_surface - 1]->triso_particle_index_]; + } + + // find material + bool found = true; + int i_cell = p.lowest_coord().cell(); + for (;; ++p.n_coord()) { + if (i_cell == C_NONE) { + int i_universe = p.lowest_coord().universe(); + const auto& univ {model::universes[i_universe]}; + + if (univ->filled_with_triso_base_ != -1) { + p.lowest_coord().cell() = + model::cell_map[univ->filled_with_triso_base_]; + found = true; + } else { + found = univ->find_cell(p); + } + if (!found) { + return found; + } + } + + i_cell = p.lowest_coord().cell(); + + Cell& c {*model::cells[i_cell]}; + if (c.type_ == Fill::MATERIAL) { + // Found a material cell which means this is the lowest coord level. + + p.cell_instance() = 0; + // Find the distribcell instance number. + if (c.distribcell_index_ >= 0) { + p.cell_instance() = cell_instance_at_level(p, p.n_coord() - 1); + } + + // Set the material and temperature. + p.material_last() = p.material(); + if (c.material_.size() > 1) { + p.material() = c.material_[p.cell_instance()]; + } else { + p.material() = c.material_[0]; + } + p.sqrtkT_last() = p.sqrtkT(); + if (c.sqrtkT_.size() > 1) { + p.sqrtkT() = c.sqrtkT_[p.cell_instance()]; + } else { + p.sqrtkT() = c.sqrtkT_[0]; + } + return found; + + } else if (c.type_ == Fill::UNIVERSE) { + //======================================================================== + //! Found a lower universe, update this coord level then search the + //! next. + + // Set the lower coordinate level universe. + auto& coor {p.coord(p.n_coord())}; + coor.universe() = c.fill_; + + // Set the position and direction. + coor.r() = p.r_local(); + coor.u() = p.u_local(); + + // Apply translation. + coor.r() -= c.translation_; + + // Apply rotation. + if (!c.rotation_.empty()) { + coor.rotate(c.rotation_); + } + i_cell = C_NONE; + } + } +} //============================================================================== void cross_lattice(GeometryState& p, const BoundaryInfo& boundary, bool verbose) diff --git a/src/particle.cpp b/src/particle.cpp index dd813013a0..67cb144935 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -577,97 +577,13 @@ void Particle::cross_surface(const Surface& surf) int i_surface = std::abs(surface()); bool verbose = settings::verbosity >= 10 || trace(); 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() = - 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++) { - coord(i).reset(); - } - if (model::surfaces[i_surface - 1]->triso_particle_index_ == -1) { - fatal_error(fmt::format("Particle cell of surface {} is not defined", - model::surfaces[i_surface - 1]->id_)); - } - 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(); - for (;; ++n_coord()) { - if (i_cell == C_NONE) { - 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() = - model::cell_map[univ->filled_with_triso_base_]; - found = true; - } else { - found = univ->find_cell(*this); - } - if (!found) { - break; - } - } - - i_cell = coord(n_coord() - 1).cell(); - - Cell& c {*model::cells[i_cell]}; - if (c.type_ == Fill::MATERIAL) { - // Found a material cell which means this is the lowest coord level. - - cell_instance() = 0; - // Find the distribcell instance number. - if (c.distribcell_index_ >= 0) { - cell_instance() = cell_instance_at_level(*this, n_coord() - 1); - } - - // Set the material and temperature. - material_last() = material(); - if (c.material_.size() > 1) { - material() = c.material_[cell_instance()]; - } else { - material() = c.material_[0]; - } - sqrtkT_last() = sqrtkT(); - if (c.sqrtkT_.size() > 1) { - sqrtkT() = c.sqrtkT_[cell_instance()]; - } else { - sqrtkT() = c.sqrtkT_[0]; - } - return; - - } else if (c.type_ == Fill::UNIVERSE) { - //======================================================================== - //! Found a lower universe, update this coord level then search the - //! next. - - // Set the lower coordinate level universe. - auto& coor {coord(n_coord())}; - coor.universe() = c.fill_; - - // Set the position and direction. - coor.r() = r_local(); - coor.u() = u_local(); - - // Apply translation. - coor.r() -= c.translation_; - - // Apply rotation. - if (!c.rotation_.empty()) { - coor.rotate(c.rotation_); - } - i_cell = C_NONE; - } + if (find_cell_in_virtual_lattice(*this, verbose)) { + return; } } else { - if (neighbor_list_find_cell(*this, verbose)) + if (neighbor_list_find_cell(*this, verbose)) { return; + } } // ========================================================================== diff --git a/src/universe.cpp b/src/universe.cpp index 5fa2892b27..054f57d060 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -40,59 +40,9 @@ 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 lat_ind(3); - Position r {p.r_local()}; - lat_ind[0] = std::max( - std::min(floor((r.x - c.vl_lower_left_[0]) / c.vl_pitch_[0]), - c.vl_shape_[0] - 1), - 0); - lat_ind[1] = std::max( - std::min(floor((r.y - c.vl_lower_left_[1]) / c.vl_pitch_[1]), - c.vl_shape_[1] - 1), - 0); - lat_ind[2] = std::max( - std::min(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 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; + bool found = find_cell_in_virtual_lattice(p); + if (found) { + return found; } } const auto& cells { @@ -114,6 +64,60 @@ bool Universe::find_cell(GeometryState& p) const } return false; } +bool Universe::find_cell_in_virtual_lattice(GeometryState& p) const +{ + Cell& c {*model::cells[model::cell_map[filled_with_triso_base_]]}; + vector lat_ind(3); + Position r {p.r_local()}; + lat_ind[0] = std::max( + std::min( + floor((r.x - c.vl_lower_left_[0]) / c.vl_pitch_[0]), c.vl_shape_[0] - 1), + 0); + lat_ind[1] = std::max( + std::min( + floor((r.y - c.vl_lower_left_[1]) / c.vl_pitch_[1]), c.vl_shape_[1] - 1), + 0); + lat_ind[2] = std::max( + std::min( + floor((r.z - c.vl_lower_left_[2]) / c.vl_pitch_[2]), c.vl_shape_[2] - 1), + 0); + + int32_t i_univ = p.lowest_coord().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 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.lowest_coord().cell() = + model::cell_map[model::surfaces[abs(token) - 1] + ->triso_particle_index_]; + return true; + } else { + p.lowest_coord().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.lowest_coord().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.lowest_coord().cell() = model::cell_map[filled_with_triso_base_]; + return true; + } + return false; +} BoundingBox Universe::bounding_box() const {