diff --git a/include/openmc/cell.h b/include/openmc/cell.h index afc9a33c596..075e907682a 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; //!< ID of cell filled with virtual lattice //! \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_; + + + //! \brief Specification of the virtual 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..52ff271418d 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -112,6 +112,11 @@ 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 +577,13 @@ 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..76437d10613 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,6 +825,12 @@ 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 @@ -843,6 +850,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 +866,12 @@ 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..b59047ffd41 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -109,6 +109,44 @@ 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 (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 +249,39 @@ 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]=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())}; @@ -481,6 +552,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 +697,18 @@ 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_); + } + + if (triso_particle_) { + 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"); + } + } + // Read the translation vector. if (check_for_node(cell_node, "translation")) { if (fill_ == C_NONE) { @@ -644,22 +750,99 @@ 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 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; + } + } - for (int32_t token : rpn_) { - // Ignore this token if it corresponds to an operator rather than a region. - if (token >= OP_UNION) - continue; + 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]; - // 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 (min_dist < tol_dis) { + break; + } - // 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; + if (norm_u[mes_bou_crossed] > 0) { + lat_ind[mes_bou_crossed] += 1; + } else { + lat_ind[mes_bou_crossed] += -1; + } + + dis_to_bou[mes_bou_crossed] += dis_to_bou_max[mes_bou_crossed]; + + if (tol_dis > max_dis) { + break; + } + } + + + } else { + 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; + } } } } @@ -1065,6 +1248,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/particle.cpp b/src/particle.cpp index b3dd72d9a94..98f0b60aa58 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()); @@ -449,8 +450,94 @@ 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) { 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()