From 7d5fd4e2cf4e158d0eb76e138050a5d8b64179ae Mon Sep 17 00:00:00 2001 From: Jingang Liang Date: Mon, 28 Jul 2025 16:37:13 +0800 Subject: [PATCH 1/4] 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/4] 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/4] 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/4] 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()