Merge branch 'virtual_lattice' into virtual_lattice_0.15.2

Author: skywalker_cn

include/openmc/cell.h

@@ -319,7 +319,15 @@ class Cell {
   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};
 

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,12 @@ 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;
 
   double x0_, y0_, z0_, radius_;
+  // int triso_base_index_ = -1;
 };
 
 //==============================================================================
@@ -264,6 +293,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 +313,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 +333,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 +353,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 +373,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 +392,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 +411,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_;
 };

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))

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
@@ -826,6 +827,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 +852,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 +868,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: