Prepare point query data structures on meshes when applying Weight Wi…

…ndows (#3157)

Co-authored-by: Paul Romano <[email protected]>

include/openmc/mesh.h

@@ -83,8 +83,8 @@ class Mesh {
   virtual ~Mesh() = default;
 
   // Methods
-  //! Perform any preparation needed to support use in mesh filters
-  virtual void prepare_for_tallies() {};
+  //! Perform any preparation needed to support point location within the mesh
+  virtual void prepare_for_point_location() {};
 
   //! Update a position to the local coordinates of the mesh
   virtual void local_coords(Position& r) const {};
@@ -737,7 +737,7 @@ class MOABMesh : public UnstructuredMesh {
   // Overridden Methods
 
   //! Perform any preparation needed to support use in mesh filters
-  void prepare_for_tallies() override;
+  void prepare_for_point_location() override;
 
   Position sample_element(int32_t bin, uint64_t* seed) const override;
 

openmc/weight_windows.py

@@ -328,8 +328,9 @@ def to_xml_element(self) -> ET.Element:
         subelement = ET.SubElement(element, 'particle_type')
         subelement.text = self.particle_type
 
-        subelement = ET.SubElement(element, 'energy_bounds')
-        subelement.text = ' '.join(str(e) for e in self.energy_bounds)
+        if self.energy_bounds is not None:
+            subelement = ET.SubElement(element, 'energy_bounds')
+            subelement.text = ' '.join(str(e) for e in self.energy_bounds)
 
         subelement = ET.SubElement(element, 'lower_ww_bounds')
         subelement.text = ' '.join(str(b) for b in self.lower_ww_bounds.ravel('F'))

src/mesh.cpp

@@ -2315,7 +2315,7 @@ void MOABMesh::initialize()
   this->determine_bounds();
 }
 
-void MOABMesh::prepare_for_tallies()
+void MOABMesh::prepare_for_point_location()
 {
   // if the KDTree has already been constructed, do nothing
   if (kdtree_)
@@ -2365,7 +2365,8 @@ 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);
+  write_message(
+    7, "Building adaptive k-d tree for tet mesh with ID {}...", id_);
   kdtree_ = make_unique<moab::AdaptiveKDTree>(mbi_.get());
 
   // Determine what options to use

src/tallies/filter_mesh.cpp

@@ -80,7 +80,7 @@ void MeshFilter::set_mesh(int32_t mesh)
   // perform any additional perparation for mesh tallies here
   mesh_ = mesh;
   n_bins_ = model::meshes[mesh_]->n_bins();
-  model::meshes[mesh_]->prepare_for_tallies();
+  model::meshes[mesh_]->prepare_for_point_location();
 }
 
 void MeshFilter::set_translation(const Position& translation)

src/weight_windows.cpp

@@ -147,8 +147,8 @@ WeightWindows::WeightWindows(int32_t id)
 WeightWindows::WeightWindows(pugi::xml_node node)
 {
   // Make sure required elements are present
-  const vector<std::string> required_elems {"id", "particle_type",
-    "energy_bounds", "lower_ww_bounds", "upper_ww_bounds"};
+  const vector<std::string> required_elems {
+    "id", "particle_type", "lower_ww_bounds", "upper_ww_bounds"};
   for (const auto& elem : required_elems) {
     if (!check_for_node(node, elem.c_str())) {
       fatal_error(fmt::format("Must specify <{}> for weight windows.", elem));
@@ -165,7 +165,7 @@ WeightWindows::WeightWindows(pugi::xml_node node)
 
   // Determine associated mesh
   int32_t mesh_id = std::stoi(get_node_value(node, "mesh"));
-  mesh_idx_ = model::mesh_map.at(mesh_id);
+  set_mesh(model::mesh_map.at(mesh_id));
 
   // energy bounds
   if (check_for_node(node, "energy_bounds"))
@@ -340,6 +340,7 @@ void WeightWindows::set_mesh(int32_t mesh_idx)
     fatal_error(fmt::format("Could not find a mesh for index {}", mesh_idx));
 
   mesh_idx_ = mesh_idx;
+  model::meshes[mesh_idx_]->prepare_for_point_location();
   allocate_ww_bounds();
 }
 

tests/unit_tests/weightwindows/test.py

@@ -296,4 +296,40 @@ def test_ww_attrs_capi(run_in_tmpdir, model):
     assert wws.id == 2
     assert wws.particle == openmc.ParticleType.PHOTON
 
-    openmc.lib.finalize()
+    openmc.lib.finalize()
+
+
+@pytest.mark.parametrize('library', ('libmesh', 'moab'))
+def test_unstructured_mesh_applied_wws(request, run_in_tmpdir, library):
+    """
+    Ensure that weight windows on unstructured mesh work when
+    they aren't part of a tally or weight window generator
+    """
+
+    if library == 'libmesh' and not openmc.lib._libmesh_enabled():
+        pytest.skip('LibMesh not enabled in this build.')
+    if library == 'moab' and not openmc.lib._dagmc_enabled():
+        pytest.skip('DAGMC (and MOAB) mesh not enabled in this build.')
+
+    water = openmc.Material(name='water')
+    water.add_nuclide('H1', 2.0)
+    water.add_nuclide('O16', 1.0)
+    water.set_density('g/cc', 1.0)
+    box = openmc.model.RectangularParallelepiped(*(3*[-10, 10]), boundary_type='vacuum')
+    cell = openmc.Cell(region=-box, fill=water)
+
+    geometry = openmc.Geometry([cell])
+    mesh_file = str(request.fspath.dirpath() / 'test_mesh_tets.exo')
+    mesh = openmc.UnstructuredMesh(mesh_file, library)
+
+    dummy_wws = np.ones((12_000,))
+
+    wws = openmc.WeightWindows(mesh, dummy_wws, upper_bound_ratio=5.0)
+
+    model = openmc.Model(geometry)
+    model.settings.weight_windows = wws
+    model.settings.weight_windows_on = True
+    model.settings.run_mode = 'fixed source'
+    model.settings.particles = 100
+    model.settings.batches = 2
+    model.run()

tests/unit_tests/weightwindows/test_mesh_tets.exo

@@ -0,0 +1 @@
+../../regression_tests/unstructured_mesh/test_mesh_tets.e