Implement Virtual Lattice Method for Efficient Simulation of Dispersed Fuels

include/openmc/cell.h

@@ -60,6 +60,7 @@ class Universe {
 public:
   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
 
   //! \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<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

@@ -88,6 +88,9 @@ class Surface {
   std::shared_ptr<BoundaryCondition> 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<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
@@ -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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> lattice_pitch) const;
+  vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> 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<double> mesh_center, vector<double> lattice_pitch) const;
   void to_hdf5_inner(hid_t group_id) const;
 
   double x0_, y0_, z0_, A_, B_, C_;

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

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: