From 7d5fd4e2cf4e158d0eb76e138050a5d8b64179ae Mon Sep 17 00:00:00 2001 From: Jingang Liang Date: Mon, 28 Jul 2025 16:37:13 +0800 Subject: [PATCH 1/9] Initial Implementation of virtual lattice, with contributions by Li Ruihan --- include/openmc/cell.h | 10 ++ include/openmc/surface.h | 27 ++++ openmc/cell.py | 14 ++ openmc/model/triso.py | 19 ++- src/cell.cpp | 270 +++++++++++++++++++++++++++++++++++++-- src/geometry.cpp | 2 +- src/particle.cpp | 100 ++++++++++++++- src/surface.cpp | 128 +++++++++++++++++++ 8 files changed, 550 insertions(+), 20 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index afc9a33c596..e6fdbe36a96 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -60,6 +60,7 @@ class Universe { public: int32_t id_; //!< Unique ID vector cells_; //!< Cells within this universe + int filled_with_triso_base_ = -1; //! \brief Write universe information to an HDF5 group. //! \param group_id An HDF5 group id. @@ -212,6 +213,15 @@ class Cell { int32_t fill_; //!< Universe # filling this cell int32_t n_instances_ {0}; //!< Number of instances of this cell GeometryType geom_type_; //!< Geometric representation type (CSG, DAGMC) + bool virtual_lattice_; //!< If the cell is the base of a virtual triso lattice + bool triso_particle_; + + + // Specification of the viryual lattice + vector vl_lower_left_; + vector vl_pitch_; + vector vl_shape_; + vector> vl_triso_distribution_; //! \brief Index corresponding to this cell in distribcell arrays int distribcell_index_ {C_NONE}; diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 84a2ed1134d..785e1a7ec19 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -88,6 +88,9 @@ class Surface { std::shared_ptr bc_ {nullptr}; //!< Boundary condition GeometryType geom_type_; //!< Geometry type indicator (CSG or DAGMC) bool surf_source_ {false}; //!< Activate source banking for the surface? + int triso_base_index_; + int triso_particle_index_ = -1; + bool is_triso_surface_ = false; explicit Surface(pugi::xml_node surf_node); Surface(); @@ -126,6 +129,11 @@ class Surface { //! exactly on the surface. virtual double distance(Position r, Direction u, bool coincident) const = 0; + virtual bool triso_in_mesh(vector mesh_center, vector lattice_pitch) const { return {}; }; + virtual void connect_to_triso_base(int triso_index, std::string key) {}; + virtual vector get_center() const { return {}; }; + virtual double get_radius() const { return {}; }; + //! Compute the local outward normal direction of the surface. //! \param r A 3D Cartesian coordinate. //! \return Normal direction @@ -164,6 +172,7 @@ class SurfaceXPlane : public CSGSurface { double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + bool triso_in_mesh(vector mesh_center, vector lattice_pitch) const; BoundingBox bounding_box(bool pos_side) const; double x0_; @@ -182,6 +191,7 @@ class SurfaceYPlane : public CSGSurface { double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + bool triso_in_mesh(vector mesh_center, vector lattice_pitch) const; BoundingBox bounding_box(bool pos_side) const; double y0_; @@ -200,6 +210,7 @@ class SurfaceZPlane : public CSGSurface { double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + bool triso_in_mesh(vector mesh_center, vector lattice_pitch) const; BoundingBox bounding_box(bool pos_side) const; double z0_; @@ -218,6 +229,7 @@ class SurfacePlane : public CSGSurface { double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + bool triso_in_mesh(vector mesh_center, vector lattice_pitch) const; double A_, B_, C_, D_; }; @@ -237,6 +249,7 @@ class SurfaceXCylinder : public CSGSurface { Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; BoundingBox bounding_box(bool pos_side) const; + bool triso_in_mesh(vector mesh_center, vector lattice_pitch) const; double y0_, z0_, radius_; }; @@ -256,6 +269,7 @@ class SurfaceYCylinder : public CSGSurface { Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; BoundingBox bounding_box(bool pos_side) const; + bool triso_in_mesh(vector mesh_center, vector lattice_pitch) const; double x0_, z0_, radius_; }; @@ -275,6 +289,7 @@ class SurfaceZCylinder : public CSGSurface { Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; BoundingBox bounding_box(bool pos_side) const; + bool triso_in_mesh(vector mesh_center, vector lattice_pitch) const; double x0_, y0_, radius_; }; @@ -291,11 +306,16 @@ class SurfaceSphere : public CSGSurface { explicit SurfaceSphere(pugi::xml_node surf_node); double evaluate(Position r) const; double distance(Position r, Direction u, bool coincident) const; + bool triso_in_mesh(vector mesh_center, vector lattice_pitch) const; + vector get_center() const; + double get_radius() const; + void connect_to_triso_base(int triso_index, std::string key); Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; BoundingBox bounding_box(bool pos_side) const; double x0_, y0_, z0_, radius_; + //int triso_base_index_ = -1; }; //============================================================================== @@ -312,6 +332,7 @@ class SurfaceXCone : public CSGSurface { double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + bool triso_in_mesh(vector mesh_center, vector lattice_pitch) const; double x0_, y0_, z0_, radius_sq_; }; @@ -330,6 +351,7 @@ class SurfaceYCone : public CSGSurface { double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + bool triso_in_mesh(vector mesh_center, vector lattice_pitch) const; double x0_, y0_, z0_, radius_sq_; }; @@ -348,6 +370,7 @@ class SurfaceZCone : public CSGSurface { double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + bool triso_in_mesh(vector mesh_center, vector lattice_pitch) const; double x0_, y0_, z0_, radius_sq_; }; @@ -366,6 +389,7 @@ class SurfaceQuadric : public CSGSurface { double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + bool triso_in_mesh(vector mesh_center, vector lattice_pitch) 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_; @@ -384,6 +408,7 @@ class SurfaceXTorus : public CSGSurface { double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + bool triso_in_mesh(vector mesh_center, vector lattice_pitch) const; double x0_, y0_, z0_, A_, B_, C_; }; @@ -401,6 +426,7 @@ class SurfaceYTorus : public CSGSurface { double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + bool triso_in_mesh(vector mesh_center, vector lattice_pitch) const; double x0_, y0_, z0_, A_, B_, C_; }; @@ -417,6 +443,7 @@ class SurfaceZTorus : public CSGSurface { double evaluate(Position r) const; double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; + bool triso_in_mesh(vector mesh_center, vector lattice_pitch) const; void to_hdf5_inner(hid_t group_id) const; double x0_, y0_, z0_, A_, B_, C_; diff --git a/openmc/cell.py b/openmc/cell.py index feba18da54f..c2b81a8aa22 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -112,6 +112,12 @@ def __init__(self, cell_id=None, name='', fill=None, region=None): self._num_instances = None self._volume = None self._atoms = None + self._triso_particle = False + self.virtual_lattice = False + self.lower_left = None + self.pitch = None + self.shape = None + def __contains__(self, point): if self.region is None: @@ -572,6 +578,14 @@ def create_xml_subelement(self, xml_element, memo=None): """ element = ET.Element("cell") element.set("id", str(self.id)) + if self._triso_particle: + element.set("triso_particle", 'true') + if self.virtual_lattice: + element.set("virtual_lattice", str(self.virtual_lattice)) + element.set("lower_left", ' '.join(map(str, self.lower_left))) + element.set("pitch", ' '.join(map(str, self.pitch))) + element.set("shape", ' '.join(map(str, self.shape))) + if len(self._name) > 0: element.set("name", str(self.name)) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index f486482aab6..8250e77a2a5 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -55,6 +55,7 @@ class TRISO(openmc.Cell): def __init__(self, outer_radius, fill, center=(0., 0., 0.)): self._surface = openmc.Sphere(r=outer_radius) + self._triso_particle = False super().__init__(fill=fill, region=-self._surface) self.center = np.asarray(center) @@ -800,7 +801,7 @@ def repel_spheres(self, p, q, d, d_new): q[:] = (q - c)*ll[0]/r + c -def create_triso_lattice(trisos, lower_left, pitch, shape, background): +def create_triso_lattice(trisos, lower_left, pitch, shape, background, virtual=False): """Create a lattice containing TRISO particles for optimized tracking. Parameters @@ -824,9 +825,16 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background): """ + if virtual: + real_pitch = copy.deepcopy(pitch) + real_shape = copy.deepcopy(shape) + pitch = [real_pitch[i]*real_shape[i] for i in range(len(real_pitch))] + shape = [1 for i in range(len(real_shape))] + lattice = openmc.RectLattice() lattice.lower_left = lower_left lattice.pitch = pitch + indices = list(np.broadcast(*np.ogrid[:shape[2], :shape[1], :shape[0]])) triso_locations = {idx: [] for idx in indices} @@ -843,6 +851,8 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background): y0=t._surface.y0, z0=t._surface.z0) t_copy.region = -t_copy._surface + if virtual: + t_copy._triso_particle = True triso_locations[idx].append(t_copy) else: warnings.warn('TRISO particle is partially or completely ' @@ -857,6 +867,13 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background): else: background_cell = openmc.Cell(fill=background) + if virtual: + background_cell.virtual_lattice = True + background_cell.pitch = real_pitch + background_cell.shape = real_shape + background_cell.lower_left = [-pitch[i]/2 for i in range(len(pitch))] + + u = openmc.Universe() u.add_cell(background_cell) for t in triso_list: diff --git a/src/cell.cpp b/src/cell.cpp index 2b6555d128d..95f143428b7 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -109,6 +109,62 @@ vector tokenize(const std::string region_spec) return tokens; } +vector> \ + generate_triso_distribution(\ + vector lattice_shape,\ + vector lattice_pitch,\ + vector lattice_lower_left,\ + vector cell_rpn, int id) +{ + vector> triso_distribution(lattice_shape[0]*lattice_shape[1]*lattice_shape[2]); + vector mesh_center(3); + vector mesh_ind(3); + /* + for (int i=0; i= OP_UNION) continue; + if (model::surfaces[abs(token) - 1]->triso_in_mesh(mesh_center, lattice_pitch)) { + triso_distribution[i+j*lattice_shape[0]+k*lattice_shape[0]*lattice_shape[1]].push_back(token); + model::surfaces[abs(token) - 1]->connect_to_triso_base(id, "base"); + } + } + } + } + }*/ + for (int32_t token : cell_rpn) { + if (token >= OP_UNION) continue; + vector triso_center=model::surfaces[abs(token) - 1]->get_center(); + for (int i=0; i<3; i++) { + mesh_ind[i]=floor((triso_center[i]-lattice_lower_left[i])/lattice_pitch[i]); + } + for (int i=mesh_ind[0]-1; i<=mesh_ind[0]+1; i++) { + for (int j=mesh_ind[1]-1; j<=mesh_ind[1]+1; j++) { + for (int k=mesh_ind[2]-1; k<=mesh_ind[2]+1; k++) { + if (i < 0 || i >= lattice_shape[0] ||\ + j < 0 || j >= lattice_shape[1] ||\ + k < 0 || k >= lattice_shape[2]) continue; + mesh_center[0]=(i+0.5)*lattice_pitch[0]+lattice_lower_left[0]; + mesh_center[1]=(j+0.5)*lattice_pitch[1]+lattice_lower_left[1]; + mesh_center[2]=(k+0.5)*lattice_pitch[2]+lattice_lower_left[2]; + if (model::surfaces[abs(token) - 1]->triso_in_mesh(mesh_center, lattice_pitch)) { + triso_distribution[i+j*lattice_shape[0]+k*lattice_shape[0]*lattice_shape[1]].push_back(token); + model::surfaces[abs(token) - 1]->connect_to_triso_base(id, "base"); + } + } + } + } + } + + + return triso_distribution; +} + + //============================================================================== //! Convert infix region specification to Reverse Polish Notation (RPN) //! @@ -211,6 +267,37 @@ void Universe::to_hdf5(hid_t universes_group) const bool Universe::find_cell(Particle& 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]=floor((r.x-c.vl_lower_left_[0])/c.vl_pitch_[0]); + lat_ind[1]=floor((r.y-c.vl_lower_left_[1])/c.vl_pitch_[1]); + lat_ind[2]=floor((r.z-c.vl_lower_left_[2])/c.vl_pitch_[2]); + 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; + } + } const auto& cells { !partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())}; @@ -481,6 +568,29 @@ CSGCell::CSGCell(pugi::xml_node cell_node) universe_ = 0; } + // Check if the cell is the base of a virtual triso lattice + bool virtual_lattice_present = check_for_node(cell_node, "virtual_lattice"); + if (virtual_lattice_present) { + virtual_lattice_ = get_node_value_bool(cell_node, "virtual_lattice"); + if (virtual_lattice_) { + if (check_for_node(cell_node, "lower_left") && check_for_node(cell_node, "pitch") && check_for_node(cell_node, "shape")) { + vl_lower_left_ = get_node_array(cell_node, "lower_left"); + vl_pitch_ = get_node_array(cell_node, "pitch"); + vl_shape_ = get_node_array(cell_node, "shape"); + } else { + fatal_error(fmt::format("Lower_left, pitch and shape of the virtual lattice must be specified for cell {}", id_)); + } + } + } else { + virtual_lattice_ = false; + } + + if (check_for_node(cell_node, "triso_particle")) { + triso_particle_ = get_node_value_bool(cell_node, "triso_particle"); + } else { + triso_particle_ = false; + } + // Make sure that either material or fill was specified, but not both. bool fill_present = check_for_node(cell_node, "fill"); bool material_present = check_for_node(cell_node, "material"); @@ -603,6 +713,25 @@ CSGCell::CSGCell(pugi::xml_node cell_node) } rpn_.shrink_to_fit(); + if (virtual_lattice_) { + vl_triso_distribution_ = generate_triso_distribution(vl_shape_, vl_pitch_, vl_lower_left_, rpn_, id_); + write_message("successfully generated virtual lattice", 5); + for (int tri_sur_id : vl_triso_distribution_[floor(vl_shape_[0]/2)+floor(vl_shape_[1]/2)*vl_shape_[0]+\ + floor(vl_shape_[2]/2)*vl_shape_[0]*vl_shape_[1]]) + { + std::cout<id_<connect_to_triso_base(id_, "particle"); + } + } + // Read the translation vector. if (check_for_node(cell_node, "translation")) { if (fill_ == C_NONE) { @@ -644,26 +773,136 @@ std::pair CSGCell::distance( { 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 tol_dis = p.collision_distance(); + 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 r_end = {r.x+tol_dis*norm_u[0], r.y+tol_dis*norm_u[1], r.z+tol_dis*norm_u[2]}; + double temp_pos_x, temp_pos_y, temp_pos_z; + vector> passed_lattice={{floor((r.x-vl_lower_left_[0])/vl_pitch_[0]),\ + floor((r.y-vl_lower_left_[1])/vl_pitch_[1]),\ + floor((r.z-vl_lower_left_[2])/vl_pitch_[2]), 0}}; + int index_start; + int index_end; + for (int i = 0; i < 3; i++){ + if (passed_lattice[0][i] == vl_shape_[i]) { + passed_lattice[0][i] = vl_shape_[i]-1; + } + } - for (int32_t token : rpn_) { - // Ignore this token if it corresponds to an operator rather than a region. - if (token >= OP_UNION) - continue; + if (u.x > 0) { + index_start = ceil((r.x-vl_lower_left_[0])/vl_pitch_[0]); + index_end = floor((r_end[0]-vl_lower_left_[0])/vl_pitch_[0]); + if (index_start <= index_end && index_start < vl_shape_[0]) { + for (int i = index_start; i <= index_end; i++) { + if (i >= vl_shape_[0]) break; + temp_pos_x = i*vl_pitch_[0]+vl_lower_left_[0]; + temp_pos_y = (temp_pos_x-r.x)*norm_u[1]/norm_u[0]+r.y; + temp_pos_z = (temp_pos_x-r.x)*norm_u[2]/norm_u[0]+r.z; + passed_lattice.push_back({i, floor((temp_pos_y-vl_lower_left_[1])/vl_pitch_[1]),\ + floor((temp_pos_z-vl_lower_left_[2])/vl_pitch_[2]),\ + (temp_pos_x-r.x)/norm_u[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]); + } + + 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; + } - // Calculate the distance to this surface. - // Note the off-by-one indexing - bool coincident {std::abs(token) == std::abs(on_surface)}; - double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)}; + if (norm_u[mes_bou_crossed] > 0) { + lat_ind[mes_bou_crossed] += 1; + } else { + lat_ind[mes_bou_crossed] += -1; + } - // Check if this distance is the new minimum. - if (d < min_dist) { - if (min_dist - d >= FP_PRECISION * min_dist) { - min_dist = d; - i_surf = -token; + dis_to_bou[mes_bou_crossed] += dis_to_bou_max[mes_bou_crossed]; + + if (tol_dis > max_dis) { + break; } } - } + + } else { + for (int32_t token : rpn_) { + // Ignore this token if it corresponds to an operator rather than a region. + if (token >= OP_UNION) + continue; + + // Calculate the distance to this surface. + // Note the off-by-one indexing + bool coincident {std::abs(token) == std::abs(on_surface)}; + double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)}; + + // Check if this distance is the new minimum. + if (d < min_dist) { + if (min_dist - d >= FP_PRECISION * min_dist) { + min_dist = d; + i_surf = -token; + } + } + } + } + return {min_dist, i_surf}; } @@ -1065,6 +1304,9 @@ void read_cells(pugi::xml_node node) } else { model::universes[it->second]->cells_.push_back(i); } + if (model::cells[i]->virtual_lattice_) { + model::universes[it->second]->filled_with_triso_base_=model::cells[i]->id_; + } } model::universes.shrink_to_fit(); diff --git a/src/geometry.cpp b/src/geometry.cpp index acce3fc78eb..6225bb76306 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -148,7 +148,7 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list) const auto& univ {model::universes[i_universe]}; found = univ->find_cell(p); } - + if (!found) { return found; } diff --git a/src/particle.cpp b/src/particle.cpp index b3dd72d9a94..065b2251610 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -184,7 +184,6 @@ void Particle::event_calculate_xs() void Particle::event_advance() { // Find the distance to the nearest boundary - boundary() = distance_to_boundary(*this); // Sample a distance to collision if (type() == ParticleType::electron || type() == ParticleType::positron) { @@ -195,6 +194,8 @@ void Particle::event_advance() collision_distance() = -std::log(prn(current_seed())) / macro_xs().total; } + boundary() = distance_to_boundary(*this); + // Select smaller of the two distances double distance = std::min(boundary().distance, collision_distance()); @@ -396,7 +397,10 @@ void Particle::cross_surface() if (settings::verbosity >= 10 || trace()) { write_message(1, " Crossing surface {}", surf->id_); } - + /* + if (surf->id_ >= 11675 && surf->id_ <=23340) { + write_message(1, " {}", surf->id_); + }*/ if (surf->surf_source_ && simulation::current_batch == settings::n_batches) { SourceSite site; site.r = r(); @@ -449,8 +453,96 @@ void Particle::cross_surface() } #endif - if (neighbor_list_find_cell(*this)) - return; + + 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; + } + } + } else { + if (neighbor_list_find_cell(*this)) + return; + } + // ========================================================================== // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS diff --git a/src/surface.cpp b/src/surface.cpp index dbfa0bf3dc2..d19387fa463 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -313,10 +313,20 @@ BoundingBox SurfaceXPlane::bounding_box(bool pos_side) const } } +bool SurfaceXPlane::triso_in_mesh(vector mesh_center, vector lattice_pitch) const +{ + return false; +} + //============================================================================== // SurfaceYPlane implementation //============================================================================== +bool SurfaceYPlane::triso_in_mesh(vector mesh_center, vector lattice_pitch) const +{ + return false; +} + SurfaceYPlane::SurfaceYPlane(pugi::xml_node surf_node) : CSGSurface(surf_node) { read_coeffs(surf_node, id_, y0_); @@ -357,6 +367,11 @@ BoundingBox SurfaceYPlane::bounding_box(bool pos_side) const // SurfaceZPlane implementation //============================================================================== +bool SurfaceZPlane::triso_in_mesh(vector mesh_center, vector lattice_pitch) const +{ + return false; +} + SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node) : CSGSurface(surf_node) { read_coeffs(surf_node, id_, z0_); @@ -396,6 +411,10 @@ BoundingBox SurfaceZPlane::bounding_box(bool pos_side) const //============================================================================== // SurfacePlane implementation //============================================================================== +bool SurfacePlane::triso_in_mesh(vector mesh_center, vector lattice_pitch) const +{ + return false; +} SurfacePlane::SurfacePlane(pugi::xml_node surf_node) : CSGSurface(surf_node) { @@ -515,6 +534,11 @@ Direction axis_aligned_cylinder_normal( // SurfaceXCylinder implementation //============================================================================== +bool SurfaceXCylinder::triso_in_mesh(vector mesh_center, vector lattice_pitch) const +{ + return false; +} + SurfaceXCylinder::SurfaceXCylinder(pugi::xml_node surf_node) : CSGSurface(surf_node) { @@ -558,6 +582,11 @@ BoundingBox SurfaceXCylinder::bounding_box(bool pos_side) const // SurfaceYCylinder implementation //============================================================================== +bool SurfaceYCylinder::triso_in_mesh(vector mesh_center, vector lattice_pitch) const +{ + return false; +} + SurfaceYCylinder::SurfaceYCylinder(pugi::xml_node surf_node) : CSGSurface(surf_node) { @@ -602,6 +631,11 @@ BoundingBox SurfaceYCylinder::bounding_box(bool pos_side) const // SurfaceZCylinder implementation //============================================================================== +bool SurfaceZCylinder::triso_in_mesh(vector mesh_center, vector lattice_pitch) const +{ + return false; +} + SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node) : CSGSurface(surf_node) { @@ -720,6 +754,66 @@ 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; + } +} + +vector SurfaceSphere::get_center() const { + return {x0_,y0_,z0_}; +} + +double SurfaceSphere::get_radius() const { + return radius_; +} + +bool SurfaceSphere::triso_in_mesh(vector mesh_center, vector 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_=y_min && y0_<=y_max) { + dis_y=0; + } else if (y0_=z_min && z0_<=z_max) { + dis_z=0; + } else if (z0_ mesh_center, vector lattice_pitch) const +{ + return false; +} + SurfaceXCone::SurfaceXCone(pugi::xml_node surf_node) : CSGSurface(surf_node) { read_coeffs(surf_node, id_, x0_, y0_, z0_, radius_sq_); @@ -844,6 +943,11 @@ void SurfaceXCone::to_hdf5_inner(hid_t group_id) const // SurfaceYCone implementation //============================================================================== +bool SurfaceYCone::triso_in_mesh(vector mesh_center, vector lattice_pitch) const +{ + return false; +} + SurfaceYCone::SurfaceYCone(pugi::xml_node surf_node) : CSGSurface(surf_node) { read_coeffs(surf_node, id_, x0_, y0_, z0_, radius_sq_); @@ -876,6 +980,11 @@ void SurfaceYCone::to_hdf5_inner(hid_t group_id) const // SurfaceZCone implementation //============================================================================== +bool SurfaceZCone::triso_in_mesh(vector mesh_center, vector lattice_pitch) const +{ + return false; +} + SurfaceZCone::SurfaceZCone(pugi::xml_node surf_node) : CSGSurface(surf_node) { read_coeffs(surf_node, id_, x0_, y0_, z0_, radius_sq_); @@ -908,6 +1017,11 @@ void SurfaceZCone::to_hdf5_inner(hid_t group_id) const // SurfaceQuadric implementation //============================================================================== +bool SurfaceQuadric::triso_in_mesh(vector mesh_center, vector lattice_pitch) const +{ + return false; +} + SurfaceQuadric::SurfaceQuadric(pugi::xml_node surf_node) : CSGSurface(surf_node) { read_coeffs(surf_node, id_, A_, B_, C_, D_, E_, F_, G_, H_, J_, K_); @@ -1058,6 +1172,11 @@ double torus_distance(double x1, double x2, double x3, double u1, double u2, // SurfaceXTorus implementation //============================================================================== +bool SurfaceXTorus::triso_in_mesh(vector mesh_center, vector lattice_pitch) const +{ + return false; +} + SurfaceXTorus::SurfaceXTorus(pugi::xml_node surf_node) : CSGSurface(surf_node) { read_coeffs(surf_node, id_, x0_, y0_, z0_, A_, B_, C_); @@ -1111,6 +1230,11 @@ Direction SurfaceXTorus::normal(Position r) const // SurfaceYTorus implementation //============================================================================== +bool SurfaceYTorus::triso_in_mesh(vector mesh_center, vector lattice_pitch) const +{ + return false; +} + SurfaceYTorus::SurfaceYTorus(pugi::xml_node surf_node) : CSGSurface(surf_node) { read_coeffs(surf_node, id_, x0_, y0_, z0_, A_, B_, C_); @@ -1163,6 +1287,10 @@ Direction SurfaceYTorus::normal(Position r) const //============================================================================== // SurfaceZTorus implementation //============================================================================== +bool SurfaceZTorus::triso_in_mesh(vector mesh_center, vector lattice_pitch) const +{ + return false; +} SurfaceZTorus::SurfaceZTorus(pugi::xml_node surf_node) : CSGSurface(surf_node) { From 13bb321d227dea8fca3bd6485321aeb59dba4829 Mon Sep 17 00:00:00 2001 From: Jingang Liang Date: Mon, 28 Jul 2025 16:48:34 +0800 Subject: [PATCH 2/9] fix bug for virtual lattice --- src/cell.cpp | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index 95f143428b7..790f5115b62 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -271,9 +271,11 @@ bool Universe::find_cell(Particle& 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]=floor((r.x-c.vl_lower_left_[0])/c.vl_pitch_[0]); - lat_ind[1]=floor((r.y-c.vl_lower_left_[1])/c.vl_pitch_[1]); - lat_ind[2]=floor((r.z-c.vl_lower_left_[2])/c.vl_pitch_[2]); + 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(); @@ -715,13 +717,6 @@ CSGCell::CSGCell(pugi::xml_node cell_node) if (virtual_lattice_) { vl_triso_distribution_ = generate_triso_distribution(vl_shape_, vl_pitch_, vl_lower_left_, rpn_, id_); - write_message("successfully generated virtual lattice", 5); - for (int tri_sur_id : vl_triso_distribution_[floor(vl_shape_[0]/2)+floor(vl_shape_[1]/2)*vl_shape_[0]+\ - floor(vl_shape_[2]/2)*vl_shape_[0]*vl_shape_[1]]) - { - std::cout<id_< CSGCell::distance( 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; + } } dis_to_bou = {INFTY,INFTY,INFTY}; From 50c3f76f551842ede98203d910fdd020e36ee91e Mon Sep 17 00:00:00 2001 From: Jingang Liang Date: Sun, 3 Aug 2025 09:22:46 +0800 Subject: [PATCH 3/9] clean up --- include/openmc/cell.h | 4 +-- openmc/cell.py | 2 -- openmc/model/triso.py | 4 +-- src/cell.cpp | 67 ++++--------------------------------------- src/geometry.cpp | 2 +- src/particle.cpp | 9 ++---- 6 files changed, 11 insertions(+), 77 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index e6fdbe36a96..075e907682a 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -60,7 +60,7 @@ class Universe { public: int32_t id_; //!< Unique ID vector cells_; //!< Cells within this universe - int filled_with_triso_base_ = -1; + int filled_with_triso_base_ = -1; //!< ID of cell filled with virtual lattice //! \brief Write universe information to an HDF5 group. //! \param group_id An HDF5 group id. @@ -217,7 +217,7 @@ class Cell { bool triso_particle_; - // Specification of the viryual lattice + //! \brief Specification of the virtual lattice vector vl_lower_left_; vector vl_pitch_; vector vl_shape_; diff --git a/openmc/cell.py b/openmc/cell.py index c2b81a8aa22..52ff271418d 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -118,7 +118,6 @@ def __init__(self, cell_id=None, name='', fill=None, region=None): self.pitch = None self.shape = None - def __contains__(self, point): if self.region is None: return True @@ -586,7 +585,6 @@ def create_xml_subelement(self, xml_element, memo=None): element.set("pitch", ' '.join(map(str, self.pitch))) element.set("shape", ' '.join(map(str, self.shape))) - if len(self._name) > 0: element.set("name", str(self.name)) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 8250e77a2a5..76437d10613 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -830,11 +830,10 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background, virtual=F real_shape = copy.deepcopy(shape) pitch = [real_pitch[i]*real_shape[i] for i in range(len(real_pitch))] shape = [1 for i in range(len(real_shape))] - + lattice = openmc.RectLattice() lattice.lower_left = lower_left lattice.pitch = pitch - indices = list(np.broadcast(*np.ogrid[:shape[2], :shape[1], :shape[0]])) triso_locations = {idx: [] for idx in indices} @@ -872,7 +871,6 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background, virtual=F background_cell.pitch = real_pitch background_cell.shape = real_shape background_cell.lower_left = [-pitch[i]/2 for i in range(len(pitch))] - u = openmc.Universe() u.add_cell(background_cell) diff --git a/src/cell.cpp b/src/cell.cpp index 790f5115b62..b59047ffd41 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -119,23 +119,7 @@ vector> \ vector> triso_distribution(lattice_shape[0]*lattice_shape[1]*lattice_shape[2]); vector mesh_center(3); vector mesh_ind(3); - /* - for (int i=0; i= OP_UNION) continue; - if (model::surfaces[abs(token) - 1]->triso_in_mesh(mesh_center, lattice_pitch)) { - triso_distribution[i+j*lattice_shape[0]+k*lattice_shape[0]*lattice_shape[1]].push_back(token); - model::surfaces[abs(token) - 1]->connect_to_triso_base(id, "base"); - } - } - } - } - }*/ + for (int32_t token : cell_rpn) { if (token >= OP_UNION) continue; vector triso_center=model::surfaces[abs(token) - 1]->get_center(); @@ -160,11 +144,9 @@ vector> \ } } - return triso_distribution; } - //============================================================================== //! Convert infix region specification to Reverse Polish Notation (RPN) //! @@ -586,7 +568,7 @@ CSGCell::CSGCell(pugi::xml_node cell_node) } else { virtual_lattice_ = false; } - + if (check_for_node(cell_node, "triso_particle")) { triso_particle_ = get_node_value_bool(cell_node, "triso_particle"); } else { @@ -771,40 +753,6 @@ std::pair CSGCell::distance( double min_dis_vl; int32_t i_surf_vl; if (virtual_lattice_) { - /* - double tol_dis = p.collision_distance(); - 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 r_end = {r.x+tol_dis*norm_u[0], r.y+tol_dis*norm_u[1], r.z+tol_dis*norm_u[2]}; - double temp_pos_x, temp_pos_y, temp_pos_z; - vector> passed_lattice={{floor((r.x-vl_lower_left_[0])/vl_pitch_[0]),\ - floor((r.y-vl_lower_left_[1])/vl_pitch_[1]),\ - floor((r.z-vl_lower_left_[2])/vl_pitch_[2]), 0}}; - int index_start; - int index_end; - for (int i = 0; i < 3; i++){ - if (passed_lattice[0][i] == vl_shape_[i]) { - passed_lattice[0][i] = vl_shape_[i]-1; - } - } - - if (u.x > 0) { - index_start = ceil((r.x-vl_lower_left_[0])/vl_pitch_[0]); - index_end = floor((r_end[0]-vl_lower_left_[0])/vl_pitch_[0]); - if (index_start <= index_end && index_start < vl_shape_[0]) { - for (int i = index_start; i <= index_end; i++) { - if (i >= vl_shape_[0]) break; - temp_pos_x = i*vl_pitch_[0]+vl_lower_left_[0]; - temp_pos_y = (temp_pos_x-r.x)*norm_u[1]/norm_u[0]+r.y; - temp_pos_z = (temp_pos_x-r.x)*norm_u[2]/norm_u[0]+r.z; - passed_lattice.push_back({i, floor((temp_pos_y-vl_lower_left_[1])/vl_pitch_[1]),\ - floor((temp_pos_z-vl_lower_left_[2])/vl_pitch_[2]),\ - (temp_pos_x-r.x)/norm_u[0]}) - } - } - } - */ - double max_dis = p->collision_distance(); double tol_dis = 0; vector dis_to_bou(3), dis_to_bou_max(3); @@ -835,14 +783,10 @@ std::pair CSGCell::distance( } 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)}; @@ -853,7 +797,6 @@ std::pair CSGCell::distance( } } } - int mes_bou_crossed=0; if (dis_to_bou[1] < dis_to_bou[0]) { @@ -876,13 +819,13 @@ std::pair CSGCell::distance( } dis_to_bou[mes_bou_crossed] += dis_to_bou_max[mes_bou_crossed]; - + if (tol_dis > max_dis) { break; } } - + } else { for (int32_t token : rpn_) { // Ignore this token if it corresponds to an operator rather than a region. @@ -903,7 +846,7 @@ std::pair CSGCell::distance( } } } - + return {min_dist, i_surf}; } diff --git a/src/geometry.cpp b/src/geometry.cpp index 6225bb76306..acce3fc78eb 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -148,7 +148,7 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list) const auto& univ {model::universes[i_universe]}; found = univ->find_cell(p); } - + if (!found) { return found; } diff --git a/src/particle.cpp b/src/particle.cpp index 065b2251610..98f0b60aa58 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -397,10 +397,7 @@ void Particle::cross_surface() if (settings::verbosity >= 10 || trace()) { write_message(1, " Crossing surface {}", surf->id_); } - /* - if (surf->id_ >= 11675 && surf->id_ <=23340) { - write_message(1, " {}", surf->id_); - }*/ + if (surf->surf_source_ && simulation::current_batch == settings::n_batches) { SourceSite site; site.r = r(); @@ -453,7 +450,6 @@ void Particle::cross_surface() } #endif - if (surf->is_triso_surface_) { if (surface() > 0){ for (int i = n_coord(); i < model::n_coord_levels; i++) { @@ -477,7 +473,7 @@ void Particle::cross_surface() 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; @@ -542,7 +538,6 @@ void Particle::cross_surface() if (neighbor_list_find_cell(*this)) return; } - // ========================================================================== // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS From 4c1c604b248f25e36fecffe9c8f89037596c12d2 Mon Sep 17 00:00:00 2001 From: Jingang Liang Date: Sun, 3 Aug 2025 16:47:46 +0800 Subject: [PATCH 4/9] add a regression test --- .../triso_virtual_lattice/__init__.py | 0 .../triso_virtual_lattice/inputs_true.dat | 278 ++++++++++++++++++ .../triso_virtual_lattice/results_true.dat | 2 + .../triso_virtual_lattice/test.py | 97 ++++++ 4 files changed, 377 insertions(+) create mode 100644 tests/regression_tests/triso_virtual_lattice/__init__.py create mode 100644 tests/regression_tests/triso_virtual_lattice/inputs_true.dat create mode 100644 tests/regression_tests/triso_virtual_lattice/results_true.dat create mode 100644 tests/regression_tests/triso_virtual_lattice/test.py diff --git a/tests/regression_tests/triso_virtual_lattice/__init__.py b/tests/regression_tests/triso_virtual_lattice/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/triso_virtual_lattice/inputs_true.dat b/tests/regression_tests/triso_virtual_lattice/inputs_true.dat new file mode 100644 index 00000000000..aac4fde5d5f --- /dev/null +++ b/tests/regression_tests/triso_virtual_lattice/inputs_true.dat @@ -0,0 +1,278 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1.0 1.0 1.0 + 12 + 1 1 1 + -0.5 -0.5 -0.5 + +11 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 4 + 0 + + + 0.0 0.0 0.0 + + + diff --git a/tests/regression_tests/triso_virtual_lattice/results_true.dat b/tests/regression_tests/triso_virtual_lattice/results_true.dat new file mode 100644 index 00000000000..3e9ca126f86 --- /dev/null +++ b/tests/regression_tests/triso_virtual_lattice/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.719897E+00 3.608153E-02 diff --git a/tests/regression_tests/triso_virtual_lattice/test.py b/tests/regression_tests/triso_virtual_lattice/test.py new file mode 100644 index 00000000000..481866a9445 --- /dev/null +++ b/tests/regression_tests/triso_virtual_lattice/test.py @@ -0,0 +1,97 @@ +import random +from math import sqrt + +import numpy as np +import openmc +import openmc.model + +from tests.testing_harness import PyAPITestHarness + + +class TRISOVirtualLatticeTestHarness(PyAPITestHarness): + def _build_inputs(self): + # Define TRISO matrials + fuel = openmc.Material() + fuel.set_density('g/cm3', 10.5) + fuel.add_nuclide('U235', 0.14154) + fuel.add_nuclide('U238', 0.85846) + fuel.add_nuclide('C0', 0.5) + fuel.add_nuclide('O16', 1.5) + + porous_carbon = openmc.Material() + porous_carbon.set_density('g/cm3', 1.0) + porous_carbon.add_nuclide('C0', 1.0) + porous_carbon.add_s_alpha_beta('c_Graphite') + + ipyc = openmc.Material() + ipyc.set_density('g/cm3', 1.90) + ipyc.add_nuclide('C0', 1.0) + ipyc.add_s_alpha_beta('c_Graphite') + + sic = openmc.Material() + sic.set_density('g/cm3', 3.20) + sic.add_nuclide('C0', 1.0) + sic.add_element('Si', 1.0) + + opyc = openmc.Material() + opyc.set_density('g/cm3', 1.87) + opyc.add_nuclide('C0', 1.0) + opyc.add_s_alpha_beta('c_Graphite') + + graphite = openmc.Material() + graphite.set_density('g/cm3', 1.1995) + graphite.add_nuclide('C0', 1.0) + graphite.add_s_alpha_beta('c_Graphite') + + # Create TRISO particles + spheres = [openmc.Sphere(r=r*1e-4) + for r in [212.5, 312.5, 347.5, 382.5]] + c1 = openmc.Cell(fill=fuel, region=-spheres[0]) + c2 = openmc.Cell(fill=porous_carbon, region=+spheres[0] & -spheres[1]) + c3 = openmc.Cell(fill=ipyc, region=+spheres[1] & -spheres[2]) + c4 = openmc.Cell(fill=sic, region=+spheres[2] & -spheres[3]) + c5 = openmc.Cell(fill=opyc, region=+spheres[3]) + inner_univ = openmc.Universe(cells=[c1, c2, c3, c4, c5]) + + # Define box to contain lattice and to pack TRISO particles in + min_x = openmc.XPlane(-0.5, boundary_type='reflective') + max_x = openmc.XPlane(0.5, boundary_type='reflective') + min_y = openmc.YPlane(-0.5, boundary_type='reflective') + max_y = openmc.YPlane(0.5, boundary_type='reflective') + min_z = openmc.ZPlane(-0.5, boundary_type='reflective') + max_z = openmc.ZPlane(0.5, boundary_type='reflective') + box_region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z + box = openmc.Cell(region=box_region) + + outer_radius = 422.5*1e-4 + centers = openmc.model.pack_spheres(radius=outer_radius, + region=box_region, num_spheres=100) + trisos = [openmc.model.TRISO(outer_radius, inner_univ, c) + for c in centers] + + # Create lattice + ll, ur = box.region.bounding_box + shape = (3, 3, 3) + pitch = (ur - ll) / shape + lattice = openmc.model.create_triso_lattice( + trisos, ll, pitch, shape, graphite, virtual=True) + box.fill = lattice + + root = openmc.Universe(0, cells=[box]) + geom = openmc.Geometry(root) + geom.export_to_xml() + + settings = openmc.Settings() + settings.batches = 4 + settings.inactive = 0 + settings.particles = 100 + settings.source = openmc.Source(space=openmc.stats.Point()) + settings.export_to_xml() + + mats = openmc.Materials([fuel, porous_carbon, ipyc, sic, opyc, graphite]) + mats.export_to_xml() + + +def test_triso_virtual_lattice(): + harness = TRISOVirtualLatticeTestHarness('statepoint.4.h5') + harness.main() From d1ee645b224537d9ff96e6ab91c5e6831ee56053 Mon Sep 17 00:00:00 2001 From: skywalker_cn <2324714651@qq.com> Date: Thu, 4 Sep 2025 10:46:15 +0000 Subject: [PATCH 5/9] Refactor cell and particle data classes to adapt the virtual lattice function --- include/openmc/cell.h | 24 +++--- include/openmc/particle_data.h | 55 ++++++++++---- include/openmc/surface.h | 3 + include/openmc/universe.h | 1 + src/cell.cpp | 9 ++- src/particle.cpp | 27 ++++--- src/surface.cpp | 130 ++++++++++++++++++++++----------- src/universe.cpp | 56 ++++++++++++++ 8 files changed, 215 insertions(+), 90 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 7f40cd263b5..30089481684 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -99,13 +99,13 @@ class Region { //! 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 generate_postfix(int32_t cell_id) const; + private: //---------------------------------------------------------------------------- // Private Methods - //! Get a vector of the region expression in postfix notation - vector 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 @@ class Cell { //---------------------------------------------------------------------------- // 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 vl_lower_left_; vector vl_pitch_; @@ -380,10 +379,7 @@ class CSGCell : public Cell { vector surfaces() const override { return region_.surfaces(); } std::pair 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 { diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index ec393845f81..2c00e4e3c2f 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -372,22 +372,52 @@ class GeometryState { // 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 @@ class GeometryState { 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 @@ class ParticleData : public GeometryState { bool trace_ {false}; - double collision_distance_; - int n_event_ {0}; int n_split_ {0}; @@ -695,9 +725,6 @@ class ParticleData : public GeometryState { // 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_; } diff --git a/include/openmc/surface.h b/include/openmc/surface.h index fb6b441002a..839bdf43836 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -274,6 +274,9 @@ class SurfaceSphere : public Surface { bool triso_in_mesh( vector mesh_center, vector lattice_pitch) const override; BoundingBox bounding_box(bool pos_side) const override; + vector 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; diff --git a/include/openmc/universe.h b/include/openmc/universe.h index b7450224f4f..2afc4b97ac7 100644 --- a/include/openmc/universe.h +++ b/include/openmc/universe.h @@ -29,6 +29,7 @@ class Universe { public: int32_t id_; //!< Unique ID vector 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. diff --git a/src/cell.cpp b/src/cell.cpp index 9bb400b24f0..a5d68c8c563 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -422,18 +422,19 @@ CSGCell::CSGCell(pugi::xml_node cell_node) // Morgans law Region region(region_spec, id_); region_ = region; + vector 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 CSGCell::distance( } } else { - region_.distance(r, u, on_surface); + return region_.distance(r, u, on_surface); } return {min_dist, i_surf}; diff --git a/src/particle.cpp b/src/particle.cpp index 6d052f2b429..dd813013a07 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -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()) { diff --git a/src/surface.cpp b/src/surface.cpp index e06ed59570b..419859efaeb 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -257,7 +257,8 @@ BoundingBox SurfaceXPlane::bounding_box(bool pos_side) const } } -bool SurfaceXPlane::triso_in_mesh(vector mesh_center, vector lattice_pitch) const +bool SurfaceXPlane::triso_in_mesh( + vector mesh_center, vector lattice_pitch) const { return false; } @@ -265,7 +266,11 @@ bool SurfaceXPlane::triso_in_mesh(vector mesh_center, vector lat //============================================================================== // SurfaceYPlane implementation //============================================================================== - +bool SurfaceYPlane::triso_in_mesh( + vector mesh_center, vector 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 mesh_center, vector 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 mesh_center, vector lattice_pitch) const +bool SurfacePlane::triso_in_mesh( + vector mesh_center, vector lattice_pitch) const { return false; } @@ -468,7 +478,8 @@ Direction axis_aligned_cylinder_normal( // SurfaceXCylinder implementation //============================================================================== -bool SurfaceXCylinder::triso_in_mesh(vector mesh_center, vector lattice_pitch) const +bool SurfaceXCylinder::triso_in_mesh( + vector mesh_center, vector 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 mesh_center, vector lattice_pitch) const +bool SurfaceYCylinder::triso_in_mesh( + vector mesh_center, vector 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 mesh_center, vector lattice_pitch) const +bool SurfaceZCylinder::triso_in_mesh( + vector mesh_center, vector 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 SurfaceSphere::get_center() const { - return {x0_,y0_,z0_}; +vector 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 mesh_center, vector lattice_pitch) const +bool SurfaceSphere::triso_in_mesh( + vector mesh_center, vector 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 && 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 && 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 && 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 mesh_center, vector 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 mesh_center, vector 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 mesh_center, vector 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 mesh_center, vector 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 mesh_center, vector 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 mesh_center, vector 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 mesh_center, vector lattice_pitch) const +bool SurfaceZTorus::triso_in_mesh( + vector mesh_center, vector lattice_pitch) const { return false; } diff --git a/src/universe.cpp b/src/universe.cpp index 78ddd54d4ff..5fa2892b27d 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -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 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; + } + } const auto& cells { !partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())}; From 3f69b394b6d0beacd21c715f4cfc6c05690a14cc Mon Sep 17 00:00:00 2001 From: skywalker_cn <2324714651@qq.com> Date: Fri, 5 Sep 2025 04:56:15 +0000 Subject: [PATCH 6/9] Add TRISO particle and virtual lattice attributes to geometry docs Updates the documentation for the geometry format to include new attributes for TRISO particle and virtual lattice handling. These attributes specify whether a cell contains a TRISO particle, whether a virtual lattice is used, and the lattice's shape if applicable. --- docs/source/io_formats/geometry.rst | 24 ++++++++++++++++++++++++ openmc/model/triso.py | 5 +++++ 2 files changed, 29 insertions(+) diff --git a/docs/source/io_formats/geometry.rst b/docs/source/io_formats/geometry.rst index 6d0a37a24fa..2947502196a 100644 --- a/docs/source/io_formats/geometry.rst +++ b/docs/source/io_formats/geometry.rst @@ -219,6 +219,30 @@ Each ```` element can have the following attributes or sub-elements: *Default*: None + :triso_particle: + If the cell is filled with a TRISO particle, use this element to mark it. + + .. note:: Only cells with spherical area can be marked. + + *Default*: false + + :virtual_lattice: + If the cell is filled with a matrix containing the TRISO particle, use this + element to mark it. This can accelerate the search speed of neutrons in the + region containing a large number of TRISO particles. + + *Default*: false + + :shape: + If the virtual_lattice is True. This element specifies the shape of the + lattice. + + .. note:: The shape of the lattice must be specified if the virtual_lattice + is True. Related methods can be referred to Liang, J., Li, R., Liu, Z., + 2024. Virtual lattice method for efficient Monte Carlo transport simulation + of dispersion nuclear fuels. Computer Physics Communications 295, 108985. + https://doi.org/10.1016/j.cpc.2023.108985 + --------------------- ```` Element diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 0a54880bfa8..ac70044491d 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -819,6 +819,11 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background, virtual=F background : openmc.Material A background material that is used anywhere within the lattice but outside a TRISO particle + virtual : bool + If True, create a virtual lattice where each cell is repeated + according to the pitch and shape. This is useful for creating a + lattice with a very large number of elements. + Default is False. Returns ------- 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 7/9] 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 | 145 ++++++++++++++++++++------------------ src/geometry.cpp | 92 ++++++++++++++++++++++++ src/particle.cpp | 92 ++---------------------- src/universe.cpp | 110 +++++++++++++++-------------- 7 files changed, 237 insertions(+), 211 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 30089481684..d2365c84c44 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -381,6 +381,9 @@ class CSGCell : public Cell { 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 107cc7d1f3e..d41a98e3018 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 2afc4b97ac7..1c9fcceabfe 100644 --- a/include/openmc/universe.h +++ b/include/openmc/universe.h @@ -38,6 +38,8 @@ class Universe { 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 a5d68c8c563..f7d077fb516 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; + } + if (lat_ind[i] == -1 && norm_u[i] > 0) { + lat_ind[i] = 0; } + } - 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]; - } + 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; + 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; - } + 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; - } + 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]; + tol_dis = dis_to_bou[mes_bou_crossed]; - if (min_dist < tol_dis) { - break; - } + 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; - } + 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]; + dis_to_bou[mes_bou_crossed] += dis_to_bou_max[mes_bou_crossed]; - if (tol_dis > max_dis) { - break; - } + if (tol_dis > max_dis) { + break; } - - } else { - return region_.distance(r, u, on_surface); } - return {min_dist, i_surf}; } diff --git a/src/geometry.cpp b/src/geometry.cpp index df485089163..f32b63ed7a5 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 dd813013a07..67cb144935b 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 5fa2892b27d..054f57d0602 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 { From 03f6b5b263bb984528862b8a32da1316a6b070d8 Mon Sep 17 00:00:00 2001 From: skywalker_cn <2324714651@qq.com> Date: Thu, 11 Sep 2025 07:44:43 +0000 Subject: [PATCH 8/9] Update the triso_virtual_lattice test --- .../triso_virtual_lattice/inputs_true.dat | 554 +++++++++--------- .../triso_virtual_lattice/results_true.dat | 2 +- .../triso_virtual_lattice/test.py | 18 +- 3 files changed, 287 insertions(+), 287 deletions(-) diff --git a/tests/regression_tests/triso_virtual_lattice/inputs_true.dat b/tests/regression_tests/triso_virtual_lattice/inputs_true.dat index aac4fde5d5f..cdf1574633f 100644 --- a/tests/regression_tests/triso_virtual_lattice/inputs_true.dat +++ b/tests/regression_tests/triso_virtual_lattice/inputs_true.dat @@ -1,278 +1,278 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1.0 1.0 1.0 - 12 - 1 1 1 - -0.5 -0.5 -0.5 - -11 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - eigenvalue - 100 - 4 - 0 - - - 0.0 0.0 0.0 - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1.0 1.0 1.0 + 4 + 1 1 1 + -0.5 -0.5 -0.5 + +3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 4 + 0 + + + 0.0 0.0 0.0 + + + + diff --git a/tests/regression_tests/triso_virtual_lattice/results_true.dat b/tests/regression_tests/triso_virtual_lattice/results_true.dat index 3e9ca126f86..b5766755996 100644 --- a/tests/regression_tests/triso_virtual_lattice/results_true.dat +++ b/tests/regression_tests/triso_virtual_lattice/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.719897E+00 3.608153E-02 +1.514485E+00 2.297724E-02 diff --git a/tests/regression_tests/triso_virtual_lattice/test.py b/tests/regression_tests/triso_virtual_lattice/test.py index 481866a9445..dbde24bfcc3 100644 --- a/tests/regression_tests/triso_virtual_lattice/test.py +++ b/tests/regression_tests/triso_virtual_lattice/test.py @@ -9,7 +9,8 @@ class TRISOVirtualLatticeTestHarness(PyAPITestHarness): - def _build_inputs(self): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) # Define TRISO matrials fuel = openmc.Material() fuel.set_density('g/cm3', 10.5) @@ -65,7 +66,7 @@ def _build_inputs(self): outer_radius = 422.5*1e-4 centers = openmc.model.pack_spheres(radius=outer_radius, - region=box_region, num_spheres=100) + region=box_region, num_spheres=100, seed=1) trisos = [openmc.model.TRISO(outer_radius, inner_univ, c) for c in centers] @@ -78,20 +79,19 @@ def _build_inputs(self): box.fill = lattice root = openmc.Universe(0, cells=[box]) - geom = openmc.Geometry(root) - geom.export_to_xml() + self._model.geometry = openmc.Geometry(root) settings = openmc.Settings() settings.batches = 4 settings.inactive = 0 settings.particles = 100 - settings.source = openmc.Source(space=openmc.stats.Point()) - settings.export_to_xml() + settings.source = openmc.IndependentSource(space=openmc.stats.Point()) + self._model.settings = settings - mats = openmc.Materials([fuel, porous_carbon, ipyc, sic, opyc, graphite]) - mats.export_to_xml() + self._model.materials = openmc.Materials([fuel, porous_carbon, ipyc, + sic, opyc, graphite]) def test_triso_virtual_lattice(): - harness = TRISOVirtualLatticeTestHarness('statepoint.4.h5') + harness = TRISOVirtualLatticeTestHarness('statepoint.4.h5', model=openmc.Model()) harness.main() From 21781e21c3118528f7c7ad294dd15135580da1d7 Mon Sep 17 00:00:00 2001 From: skywalker_cn <2324714651@qq.com> Date: Thu, 11 Sep 2025 08:13:58 +0000 Subject: [PATCH 9/9] Fix the format error --- include/openmc/universe.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/openmc/universe.h b/include/openmc/universe.h index 1c9fcceabfe..796f9e6508e 100644 --- a/include/openmc/universe.h +++ b/include/openmc/universe.h @@ -27,10 +27,10 @@ extern vector> universes; class Universe { public: - int32_t id_; //!< Unique ID - vector cells_; //!< Cells within this universe + int32_t id_; //!< Unique ID + vector 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 + int32_t n_instances_; //!< Number of instances of this universe //! \brief Write universe information to an HDF5 group. //! \param group_id An HDF5 group id.