Allow MOAB k-d tree to be configured (#2976)

docs/source/io_formats/statepoint.rst

@@ -96,6 +96,8 @@ The current version of the statepoint file format is 18.1.
               - **library** (*char[]*) -- Mesh library used to represent the
                                           mesh ("moab" or "libmesh").
               - **length_multiplier** (*double*) Scaling factor applied to the mesh.
+              - **options** (*char[]*) -- Special options that control spatial
+                                          search data structures used.
               - **volumes** (*double[]*) -- Volume of each mesh cell.
               - **vertices** (*double[]*) -- x, y, z values of the mesh vertices.
               - **connectivity** (*int[]*) -- Connectivity array for the mesh

docs/source/io_formats/tallies.rst

@@ -364,6 +364,10 @@ attributes/sub-elements:
     The mesh library used to represent an unstructured mesh. This can be either
     "moab" or "libmesh". (For unstructured mesh only.)
 
+  :options:
+    Special options that control spatial search data structures used. (For
+    unstructured mesh using MOAB only)
+
   :filename:
     The name of the mesh file to be loaded at runtime. (For unstructured mesh
     only.)

include/openmc/mesh.h

@@ -659,11 +659,6 @@ class UnstructuredMesh : public Mesh {
   //! Set the length multiplier to apply to each point in the mesh
   void set_length_multiplier(const double length_multiplier);
 
-  // Data members
-  double length_multiplier_ {
-    1.0}; //!< Constant multiplication factor to apply to mesh coordinates
-  bool specified_length_multiplier_ {false};
-
   //! Sample barycentric coordinates given a seed and the vertex positions and
   //! return the sampled position
   //
@@ -672,6 +667,11 @@ class UnstructuredMesh : public Mesh {
   //! \return Sampled position within the tetrahedron
   Position sample_tet(std::array<Position, 4> coords, uint64_t* seed) const;
 
+  // Data members
+  double length_multiplier_ {
+    -1.0};              //!< Multiplicative factor applied to mesh coordinates
+  std::string options_; //!< Options for search data structures
+
 private:
   //! Setup method for the mesh. Builds data structures,
   //! sets up element mapping, creates bounding boxes, etc.

openmc/filter.py

@@ -1,3 +1,4 @@
+from __future__ import annotations
 from abc import ABCMeta
 from collections.abc import Iterable
 import hashlib
@@ -804,8 +805,8 @@ class MeshFilter(Filter):
     id : int
         Unique identifier for the filter
     translation : Iterable of float
-        This array specifies a vector that is used to translate (shift)
-        the mesh for this filter
+        This array specifies a vector that is used to translate (shift) the mesh
+        for this filter
     bins : list of tuple
         A list of mesh indices for each filter bin, e.g. [(1, 1, 1), (2, 1, 1),
         ...]
@@ -846,7 +847,6 @@ def from_hdf5(cls, group, **kwargs):
         mesh_obj = kwargs['meshes'][mesh_id]
         filter_id = int(group.name.split('/')[-1].lstrip('filter '))
 
-
         out = cls(mesh_obj, filter_id=filter_id)
 
         translation = group.get('translation')
@@ -972,7 +972,7 @@ def to_xml_element(self):
         return element
 
     @classmethod
-    def from_xml_element(cls, elem, **kwargs):
+    def from_xml_element(cls, elem: ET.Element, **kwargs) -> MeshFilter:
         mesh_id = int(get_text(elem, 'bins'))
         mesh_obj = kwargs['meshes'][mesh_id]
         filter_id = int(elem.get('id'))

openmc/mesh.py

@@ -1952,6 +1952,11 @@ class UnstructuredMesh(MeshBase):
         Name of the mesh
     length_multiplier: float
         Constant multiplier to apply to mesh coordinates
+    options : str, optional
+        Special options that control spatial search data structures used. This
+        is currently only used to set `parameters
+        <https://tinyurl.com/kdtree-params>`_ for MOAB's AdaptiveKDTree. If
+        None, OpenMC internally uses a default of "MAX_DEPTH=20;PLANE_SET=2;".
 
     Attributes
     ----------
@@ -1965,6 +1970,11 @@ class UnstructuredMesh(MeshBase):
         Multiplicative factor to apply to mesh coordinates
     library : {'moab', 'libmesh'}
         Mesh library used for the unstructured mesh tally
+    options : str
+        Special options that control spatial search data structures used. This
+        is currently only used to set `parameters
+        <https://tinyurl.com/kdtree-params>`_ for MOAB's AdaptiveKDTree. If
+        None, OpenMC internally uses a default of "MAX_DEPTH=20;PLANE_SET=2;".
     output : bool
         Indicates whether or not automatic tally output should be generated for
         this mesh
@@ -1998,7 +2008,8 @@ class UnstructuredMesh(MeshBase):
     _LINEAR_HEX = 1
 
     def __init__(self, filename: PathLike, library: str, mesh_id: Optional[int] = None,
-                 name: str = '', length_multiplier: float = 1.0):
+                 name: str = '', length_multiplier: float = 1.0,
+                 options: Optional[str] = None):
         super().__init__(mesh_id, name)
         self.filename = filename
         self._volumes = None
@@ -2008,6 +2019,7 @@ def __init__(self, filename: PathLike, library: str, mesh_id: Optional[int] = No
         self.library = library
         self._output = False
         self.length_multiplier = length_multiplier
+        self.options = options
         self._has_statepoint_data = False
 
     @property
@@ -2028,6 +2040,15 @@ def library(self, lib: str):
         cv.check_value('Unstructured mesh library', lib, ('moab', 'libmesh'))
         self._library = lib
 
+    @property
+    def options(self) -> Optional[str]:
+        return self._options
+
+    @options.setter
+    def options(self, options: Optional[str]):
+        cv.check_type('options', options, (str, type(None)))
+        self._options = options
+
     @property
     @require_statepoint_data
     def size(self):
@@ -2139,6 +2160,8 @@ def __repr__(self):
         if self.length_multiplier != 1.0:
             string += '{: <16}=\t{}\n'.format('\tLength multiplier',
                                               self.length_multiplier)
+        if self.options is not None:
+            string += '{: <16}=\t{}\n'.format('\tOptions', self.options)
         return string
 
     @property
@@ -2294,8 +2317,12 @@ def from_hdf5(cls, group: h5py.Group):
         mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
         filename = group['filename'][()].decode()
         library = group['library'][()].decode()
+        if 'options' in group.attrs:
+            options = group.attrs['options'].decode()
+        else:
+            options = None
 
-        mesh = cls(filename=filename, library=library, mesh_id=mesh_id)
+        mesh = cls(filename=filename, library=library, mesh_id=mesh_id, options=options)
         mesh._has_statepoint_data = True
         vol_data = group['volumes'][()]
         mesh.volumes = np.reshape(vol_data, (vol_data.shape[0],))
@@ -2326,6 +2353,8 @@ def to_xml_element(self):
         element.set("id", str(self._id))
         element.set("type", "unstructured")
         element.set("library", self._library)
+        if self.options is not None:
+            element.set('options', self.options)
         subelement = ET.SubElement(element, "filename")
         subelement.text = str(self.filename)
 
@@ -2352,8 +2381,9 @@ def from_xml_element(cls, elem: ET.Element):
         filename = get_text(elem, 'filename')
         library = get_text(elem, 'library')
         length_multiplier = float(get_text(elem, 'length_multiplier', 1.0))
+        options = elem.get('options')
 
-        return cls(filename, library, mesh_id, '', length_multiplier)
+        return cls(filename, library, mesh_id, '', length_multiplier, options)
 
 
 def _read_meshes(elem):

src/mesh.cpp

@@ -44,6 +44,10 @@
 #include "libmesh/numeric_vector.h"
 #endif
 
+#ifdef DAGMC
+#include "moab/FileOptions.hpp"
+#endif
+
 namespace openmc {
 
 //==============================================================================
@@ -291,7 +295,6 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)
   // check if a length unit multiplier was specified
   if (check_for_node(node, "length_multiplier")) {
     length_multiplier_ = std::stod(get_node_value(node, "length_multiplier"));
-    specified_length_multiplier_ = true;
   }
 
   // get the filename of the unstructured mesh to load
@@ -305,6 +308,10 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)
       "No filename supplied for unstructured mesh with ID: {}", id_));
   }
 
+  if (check_for_node(node, "options")) {
+    options_ = get_node_value(node, "options");
+  }
+
   // check if mesh tally data should be written with
   // statepoint files
   if (check_for_node(node, "output")) {
@@ -367,8 +374,11 @@ void UnstructuredMesh::to_hdf5(hid_t group) const
   write_dataset(mesh_group, "type", mesh_type);
   write_dataset(mesh_group, "filename", filename_);
   write_dataset(mesh_group, "library", this->library());
+  if (!options_.empty()) {
+    write_attribute(mesh_group, "options", options_);
+  }
 
-  if (specified_length_multiplier_)
+  if (length_multiplier_ > 0.0)
     write_dataset(mesh_group, "length_multiplier", length_multiplier_);
 
   // write vertex coordinates
@@ -428,9 +438,6 @@ void UnstructuredMesh::to_hdf5(hid_t group) const
 void UnstructuredMesh::set_length_multiplier(double length_multiplier)
 {
   length_multiplier_ = length_multiplier;
-
-  if (length_multiplier_ != 1.0)
-    specified_length_multiplier_ = true;
 }
 
 ElementType UnstructuredMesh::element_type(int bin) const
@@ -2231,7 +2238,7 @@ void MOABMesh::initialize()
     fatal_error("Failed to add tetrahedra to an entity set.");
   }
 
-  if (specified_length_multiplier_) {
+  if (length_multiplier_ > 0.0) {
     // get the connectivity of all tets
     moab::Range adj;
     rval = mbi_->get_adjacencies(ehs_, 0, true, adj, moab::Interface::UNION);
@@ -2291,6 +2298,7 @@ void MOABMesh::build_kdtree(const moab::Range& all_tets)
 {
   moab::Range all_tris;
   int adj_dim = 2;
+  write_message("Getting tet adjacencies...", 7);
   moab::ErrorCode rval = mbi_->get_adjacencies(
     all_tets, adj_dim, true, all_tris, moab::Interface::UNION);
   if (rval != moab::MB_SUCCESS) {
@@ -2309,10 +2317,20 @@ void MOABMesh::build_kdtree(const moab::Range& all_tets)
   all_tets_and_tris.merge(all_tris);
 
   // create a kd-tree instance
+  write_message("Building adaptive k-d tree for tet mesh...", 7);
   kdtree_ = make_unique<moab::AdaptiveKDTree>(mbi_.get());
 
-  // build the tree
-  rval = kdtree_->build_tree(all_tets_and_tris, &kdtree_root_);
+  // Determine what options to use
+  std::ostringstream options_stream;
+  if (options_.empty()) {
+    options_stream << "MAX_DEPTH=20;PLANE_SET=2;";
+  } else {
+    options_stream << options_;
+  }
+  moab::FileOptions file_opts(options_stream.str().c_str());
+
+  // Build the k-d tree
+  rval = kdtree_->build_tree(all_tets_and_tris, &kdtree_root_, &file_opts);
   if (rval != moab::MB_SUCCESS) {
     fatal_error("Failed to construct KDTree for the "
                 "unstructured mesh file: " +
@@ -2899,7 +2917,7 @@ void LibMesh::initialize()
   // assuming that unstructured meshes used in OpenMC are 3D
   n_dimension_ = 3;
 
-  if (specified_length_multiplier_) {
+  if (length_multiplier_ > 0.0) {
     libMesh::MeshTools::Modification::scale(*m_, length_multiplier_);
   }
   // if OpenMC is managing the libMesh::MeshBase instance, prepare the mesh.

tests/regression_tests/unstructured_mesh/inputs_true10.dat

@@ -67,7 +67,7 @@
       <lower_left>-10.0 -10.0 -10.0</lower_left>
       <upper_right>10.0 10.0 10.0</upper_right>
     </mesh>
-    <mesh id="2" library="moab" type="unstructured">
+    <mesh id="2" library="moab" options="MAX_DEPTH=15;PLANE_SET=2" type="unstructured">
       <filename>test_mesh_tets_w_holes.e</filename>
     </mesh>
     <filter id="1" type="mesh">

tests/regression_tests/unstructured_mesh/inputs_true11.dat

@@ -67,7 +67,7 @@
       <lower_left>-10.0 -10.0 -10.0</lower_left>
       <upper_right>10.0 10.0 10.0</upper_right>
     </mesh>
-    <mesh id="2" library="moab" type="unstructured">
+    <mesh id="2" library="moab" options="MAX_DEPTH=15;PLANE_SET=2" type="unstructured">
       <filename>test_mesh_tets.e</filename>
     </mesh>
     <filter id="1" type="mesh">

tests/regression_tests/unstructured_mesh/inputs_true12.dat

@@ -67,7 +67,7 @@
       <lower_left>-10.0 -10.0 -10.0</lower_left>
       <upper_right>10.0 10.0 10.0</upper_right>
     </mesh>
-    <mesh id="2" library="moab" type="unstructured">
+    <mesh id="2" library="moab" options="MAX_DEPTH=15;PLANE_SET=2" type="unstructured">
       <filename>test_mesh_tets_w_holes.e</filename>
     </mesh>
     <filter id="1" type="mesh">

tests/regression_tests/unstructured_mesh/inputs_true13.dat

@@ -67,7 +67,7 @@
       <lower_left>-10.0 -10.0 -10.0</lower_left>
       <upper_right>10.0 10.0 10.0</upper_right>
     </mesh>
-    <mesh id="2" library="moab" type="unstructured">
+    <mesh id="2" library="moab" options="MAX_DEPTH=15;PLANE_SET=2" type="unstructured">
       <filename>test_mesh_tets.e</filename>
     </mesh>
     <filter id="1" type="mesh">

tests/regression_tests/unstructured_mesh/inputs_true14.dat

@@ -67,7 +67,7 @@
       <lower_left>-10.0 -10.0 -10.0</lower_left>
       <upper_right>10.0 10.0 10.0</upper_right>
     </mesh>
-    <mesh id="2" library="moab" type="unstructured">
+    <mesh id="2" library="moab" options="MAX_DEPTH=15;PLANE_SET=2" type="unstructured">
       <filename>test_mesh_tets_w_holes.e</filename>
     </mesh>
     <filter id="1" type="mesh">

tests/regression_tests/unstructured_mesh/inputs_true15.dat

@@ -67,7 +67,7 @@
       <lower_left>-10.0 -10.0 -10.0</lower_left>
       <upper_right>10.0 10.0 10.0</upper_right>
     </mesh>
-    <mesh id="2" library="moab" type="unstructured">
+    <mesh id="2" library="moab" options="MAX_DEPTH=15;PLANE_SET=2" type="unstructured">
       <filename>test_mesh_tets.e</filename>
     </mesh>
     <filter id="1" type="mesh">

tests/regression_tests/unstructured_mesh/inputs_true8.dat

@@ -67,7 +67,7 @@
       <lower_left>-10.0 -10.0 -10.0</lower_left>
       <upper_right>10.0 10.0 10.0</upper_right>
     </mesh>
-    <mesh id="2" library="moab" type="unstructured">
+    <mesh id="2" library="moab" options="MAX_DEPTH=15;PLANE_SET=2" type="unstructured">
       <filename>test_mesh_tets_w_holes.e</filename>
     </mesh>
     <filter id="1" type="mesh">

tests/regression_tests/unstructured_mesh/inputs_true9.dat

@@ -67,7 +67,7 @@
       <lower_left>-10.0 -10.0 -10.0</lower_left>
       <upper_right>10.0 10.0 10.0</upper_right>
     </mesh>
-    <mesh id="2" library="moab" type="unstructured">
+    <mesh id="2" library="moab" options="MAX_DEPTH=15;PLANE_SET=2" type="unstructured">
       <filename>test_mesh_tets.e</filename>
     </mesh>
     <filter id="1" type="mesh">

tests/regression_tests/unstructured_mesh/test.py

@@ -277,6 +277,8 @@ def test_unstructured_mesh_tets(model, test_opts):
 
     # add analagous unstructured mesh tally
     uscd_mesh = openmc.UnstructuredMesh(mesh_filename, test_opts['library'])
+    if test_opts['library'] == 'moab':
+        uscd_mesh.options = 'MAX_DEPTH=15;PLANE_SET=2'
     uscd_filter = openmc.MeshFilter(mesh=uscd_mesh)
 
     # create tallies