Resolve merge conflicts for PR #3516 - Implement Virtual Lattice Method for Efficient Simulation of Dispersed Fuels

docs/source/io_formats/geometry.rst

@@ -219,6 +219,30 @@ Each ``<cell>`` 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
+
 
 ---------------------
 ``<lattice>`` Element

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<int32_t> generate_postfix(int32_t cell_id) const;
+
 private:
   //----------------------------------------------------------------------------
   // Private Methods
 
-  //! Get a vector of the region expression in postfix notation
-  vector<int32_t> 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<double> vl_lower_left_;
+  vector<double> vl_pitch_;
+  vector<int32_t> vl_shape_;
+  vector<vector<int32_t>> 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<int32_t> surfaces() const override { return region_.surfaces(); }
 
   std::pair<double, int32_t> 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<double, int32_t> 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
   {

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.
 //==============================================================================

include/openmc/particle_data.h

@@ -1,3 +1,3 @@
 #ifndef OPENMC_PARTICLE_DATA_H
 #define OPENMC_PARTICLE_DATA_H
 
@@ -372,22 +372,52 @@
   // 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 @@
   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 @@
 
   bool trace_ {false};
 
-  double collision_distance_;
-
   int n_event_ {0};
 
   int n_split_ {0};
@@ -695,9 +725,6 @@
   // 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_; }
 

include/openmc/surface.h

@@ -39,6 +39,9 @@ class Surface {
   std::string name_;                 //!< User-defined name
   unique_ptr<BoundaryCondition> 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<double> mesh_center, vector<double> lattice_pitch) const
+  {
+    return {};
+  };
+  virtual void connect_to_triso_base(int triso_index, std::string key) {};
+  virtual vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> lattice_pitch) const override;
   BoundingBox bounding_box(bool pos_side) const override;
+  vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> lattice_pitch) const override;
 
   double x0_, y0_, z0_, A_, B_, C_;
 };

include/openmc/universe.h

@@ -27,16 +27,19 @@ extern vector<unique_ptr<Universe>> universes;
 
 class Universe {
 public:
-  int32_t id_;            //!< Unique ID
-  vector<int32_t> cells_; //!< Cells within this universe
-  int32_t n_instances_;   //!< Number of instances of this universe
+  int32_t id_;                      //!< Unique ID
+  vector<int32_t> 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.
   virtual void to_hdf5(hid_t group_id) const;
 
   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