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/include/openmc/cell.h b/include/openmc/cell.h index e9fcfe3911b..d2365c84c44 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,12 +314,19 @@ 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 - + 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_; + vector vl_shape_; + vector> vl_triso_distribution_; //! \brief Index corresponding to this cell in distribcell arrays int distribcell_index_ {C_NONE}; @@ -372,10 +379,10 @@ 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; + + 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 { 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/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 549e8f6abac..839bdf43836 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -39,6 +39,9 @@ class Surface { std::string name_; //!< User-defined name unique_ptr bc_; //!< Boundary condition 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(); @@ -78,6 +81,15 @@ 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 @@ -116,6 +128,8 @@ class SurfaceXPlane : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + bool triso_in_mesh( + vector mesh_center, vector lattice_pitch) const override; BoundingBox bounding_box(bool pos_side) const override; double x0_; @@ -134,6 +148,8 @@ class SurfaceYPlane : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + bool triso_in_mesh( + vector mesh_center, vector lattice_pitch) const override; BoundingBox bounding_box(bool pos_side) const override; double y0_; @@ -152,6 +168,8 @@ class SurfaceZPlane : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + bool triso_in_mesh( + vector mesh_center, vector lattice_pitch) const override; BoundingBox bounding_box(bool pos_side) const override; double z0_; @@ -169,6 +187,8 @@ class SurfacePlane : public Surface { double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; + bool triso_in_mesh( + vector mesh_center, vector lattice_pitch) const override; void to_hdf5_inner(hid_t group_id) const override; double A_, B_, C_, D_; @@ -188,6 +208,8 @@ class SurfaceXCylinder : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + bool triso_in_mesh( + vector mesh_center, vector lattice_pitch) const override; BoundingBox bounding_box(bool pos_side) const override; double y0_, z0_, radius_; @@ -207,6 +229,8 @@ class SurfaceYCylinder : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + bool triso_in_mesh( + vector mesh_center, vector lattice_pitch) const override; BoundingBox bounding_box(bool pos_side) const override; double x0_, z0_, radius_; @@ -226,6 +250,8 @@ class SurfaceZCylinder : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + bool triso_in_mesh( + vector mesh_center, vector lattice_pitch) const override; BoundingBox bounding_box(bool pos_side) const override; double x0_, y0_, radius_; @@ -245,9 +271,15 @@ class SurfaceSphere : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + 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; }; //============================================================================== @@ -264,6 +296,8 @@ class SurfaceXCone : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + bool triso_in_mesh( + vector mesh_center, vector lattice_pitch) const override; double x0_, y0_, z0_, radius_sq_; }; @@ -282,6 +316,8 @@ class SurfaceYCone : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + bool triso_in_mesh( + vector mesh_center, vector lattice_pitch) const override; double x0_, y0_, z0_, radius_sq_; }; @@ -300,6 +336,8 @@ class SurfaceZCone : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + bool triso_in_mesh( + vector mesh_center, vector lattice_pitch) const override; double x0_, y0_, z0_, radius_sq_; }; @@ -318,6 +356,8 @@ class SurfaceQuadric : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + bool triso_in_mesh( + vector mesh_center, vector lattice_pitch) const override; // 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_; @@ -336,6 +376,8 @@ class SurfaceXTorus : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + bool triso_in_mesh( + vector mesh_center, vector lattice_pitch) const override; double x0_, y0_, z0_, A_, B_, C_; }; @@ -353,6 +395,8 @@ class SurfaceYTorus : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + bool triso_in_mesh( + vector mesh_center, vector lattice_pitch) const override; double x0_, y0_, z0_, A_, B_, C_; }; @@ -370,6 +414,8 @@ class SurfaceZTorus : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + bool triso_in_mesh( + vector mesh_center, vector lattice_pitch) const override; double x0_, y0_, z0_, A_, B_, C_; }; diff --git a/include/openmc/universe.h b/include/openmc/universe.h index b7450224f4f..796f9e6508e 100644 --- a/include/openmc/universe.h +++ b/include/openmc/universe.h @@ -27,9 +27,10 @@ extern vector> universes; class Universe { public: - int32_t id_; //!< Unique ID - vector cells_; //!< Cells within this universe - int32_t n_instances_; //!< Number of instances of 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 //! \brief Write universe information to an HDF5 group. //! \param group_id An HDF5 group id. @@ -37,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/openmc/cell.py b/openmc/cell.py index 88aaebc8f27..7d450eb7feb 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -114,6 +114,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: @@ -593,6 +598,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 f344b8edd42..ac70044491d 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) @@ -802,7 +803,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 @@ -818,6 +819,11 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background): 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 ------- @@ -826,6 +832,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 @@ -845,6 +857,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 ' @@ -859,6 +873,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 d838bfbb4d1..f7d077fb516 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -36,6 +36,47 @@ vector> cells; } // namespace model +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; +} + //============================================================================== // Cell implementation //============================================================================== @@ -272,6 +313,33 @@ 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"); @@ -354,6 +422,21 @@ 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_); + } + + 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")) { @@ -420,6 +503,105 @@ vector::iterator CSGCell::find_left_parenthesis( } return it; } +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; + + 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]; + } + } + + while (true) { + if (lat_ind[0] < 0 || lat_ind[0] >= vl_shape_[0] || lat_ind[1] < 0 || + lat_ind[1] >= vl_shape_[1] || lat_ind[2] < 0 || + lat_ind[2] >= vl_shape_[2]) + break; + + for (int token : + vl_triso_distribution_[lat_ind[0] + lat_ind[1] * vl_shape_[0] + + lat_ind[2] * vl_shape_[0] * vl_shape_[1]]) { + bool coincident {std::abs(token) == std::abs(on_surface)}; + double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)}; + if (d < min_dist) { + if (min_dist - d >= FP_PRECISION * min_dist) { + min_dist = d; + i_surf = -token; + } + } + } + + int mes_bou_crossed = 0; + if (dis_to_bou[1] < dis_to_bou[0]) { + mes_bou_crossed = 1; + } + if (dis_to_bou[2] < dis_to_bou[mes_bou_crossed]) { + mes_bou_crossed = 2; + } + + tol_dis = dis_to_bou[mes_bou_crossed]; + + if (min_dist < tol_dis) { + break; + } + + if (norm_u[mes_bou_crossed] > 0) { + lat_ind[mes_bou_crossed] += 1; + } else { + lat_ind[mes_bou_crossed] += -1; + } + + dis_to_bou[mes_bou_crossed] += dis_to_bou_max[mes_bou_crossed]; + + if (tol_dis > max_dis) { + break; + } + } + return {min_dist, i_surf}; +} //============================================================================== // Region implementation @@ -1042,6 +1224,10 @@ void populate_universes() model::universes[it->second]->cells_.push_back(index_cell); } + if (model::cells[index_cell]->virtual_lattice_) { + model::universes[it->second]->filled_with_triso_base_ = + model::cells[index_cell]->id_; + } } // Add DAGUniverse implicit complement cells last 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 fb41d82e950..67cb144935b 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -225,9 +225,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) { collision_distance() = 0.0; @@ -237,6 +234,9 @@ 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(); double time_cutoff = settings::time_cutoff[static_cast(type())]; double distance_cutoff = @@ -574,10 +574,16 @@ void Particle::cross_surface(const Surface& surf) return; } #endif - + int i_surface = std::abs(surface()); bool verbose = settings::verbosity >= 10 || trace(); - if (neighbor_list_find_cell(*this, verbose)) { - return; + if (surf.is_triso_surface_) { + if (find_cell_in_virtual_lattice(*this, verbose)) { + return; + } + } else { + if (neighbor_list_find_cell(*this, verbose)) { + return; + } } // ========================================================================== diff --git a/src/surface.cpp b/src/surface.cpp index ea19514a311..419859efaeb 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -257,10 +257,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) : Surface(surf_node) { read_coeffs(surf_node, id_, {&y0_}); @@ -300,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_}); @@ -340,6 +354,11 @@ 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) : Surface(surf_node) { @@ -459,6 +478,12 @@ 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) : Surface(surf_node) { @@ -502,6 +527,12 @@ 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) : Surface(surf_node) { @@ -546,6 +577,12 @@ 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) : Surface(surf_node) { @@ -664,6 +701,69 @@ 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_ < x_min) { + dis_x = pow(x_min - x0_, 2); + } else { + dis_x = pow(x_max - x0_, 2); + } + + if (y0_ >= y_min && y0_ <= y_max) { + dis_y = 0; + } else if (y0_ < y_min) { + dis_y = pow(y_min - y0_, 2); + } else { + dis_y = pow(y_max - y0_, 2); + } + + if (z0_ >= z_min && z0_ <= z_max) { + dis_z = 0; + } else if (z0_ < z_min) { + dis_z = pow(z_min - z0_, 2); + } else { + dis_z = pow(z_max - z0_, 2); + } + + if (sqrt(dis_x + dis_y + dis_z) < radius_) { + return true; + } else { + return false; + } +} + //============================================================================== // Generic functions for x-, y-, and z-, cones //============================================================================== @@ -755,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_}); @@ -787,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_}); @@ -819,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) { @@ -851,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( @@ -1011,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_}); @@ -1064,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_}); @@ -1117,6 +1242,11 @@ 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) : Surface(surf_node) { diff --git a/src/universe.cpp b/src/universe.cpp index 78ddd54d4ff..054f57d0602 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -39,6 +39,12 @@ void Universe::to_hdf5(hid_t universes_group) const bool Universe::find_cell(GeometryState& p) const { + if (filled_with_triso_base_ != -1) { + bool found = find_cell_in_virtual_lattice(p); + if (found) { + return found; + } + } const auto& cells { !partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())}; @@ -58,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 { 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..cdf1574633f --- /dev/null +++ b/tests/regression_tests/triso_virtual_lattice/inputs_true.dat @@ -0,0 +1,278 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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 new file mode 100644 index 00000000000..b5766755996 --- /dev/null +++ b/tests/regression_tests/triso_virtual_lattice/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +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 new file mode 100644 index 00000000000..dbde24bfcc3 --- /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 __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + # 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, seed=1) + 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]) + self._model.geometry = openmc.Geometry(root) + + settings = openmc.Settings() + settings.batches = 4 + settings.inactive = 0 + settings.particles = 100 + settings.source = openmc.IndependentSource(space=openmc.stats.Point()) + self._model.settings = settings + + self._model.materials = openmc.Materials([fuel, porous_carbon, ipyc, + sic, opyc, graphite]) + + +def test_triso_virtual_lattice(): + harness = TRISOVirtualLatticeTestHarness('statepoint.4.h5', model=openmc.Model()) + harness.main()