Merge branch 'develop' into adding_outlines_to_plot

include/openmc/geometry_aux.h

@@ -126,6 +126,14 @@ std::string distribcell_path(
 
 int maximum_levels(int32_t univ);
 
+//==============================================================================
+//! Check whether or not a universe is the root universe using its ID.
+//! \param univ_id The ID of the universe to check.
+//! \return Whether or not it is the root universe.
+//==============================================================================
+
+bool is_root_universe(int32_t univ_id);
+
 //==============================================================================
 //! Deallocates global vectors and maps for cells, universes, and lattices.
 //==============================================================================

openmc/lattice.py

@@ -34,16 +34,16 @@ class Lattice(IDManagerMixin, ABC):
         Name of the lattice
     pitch : Iterable of float
         Pitch of the lattice in each direction in cm
-    outer : openmc.Universe
+    outer : openmc.UniverseBase
         A universe to fill all space outside the lattice
-    universes : Iterable of Iterable of openmc.Universe
+    universes : Iterable of Iterable of openmc.UniverseBase
         A two-or three-dimensional list/array of universes filling each element
         of the lattice
 
     """
 
     next_id = 1
-    used_ids = openmc.Universe.used_ids
+    used_ids = openmc.UniverseBase.used_ids
 
     def __init__(self, lattice_id=None, name=''):
         # Initialize Lattice class attributes
@@ -79,7 +79,7 @@ class Lattice(IDManagerMixin, ABC):
 
     @outer.setter
     def outer(self, outer):
-        cv.check_type('outer universe', outer, openmc.Universe)
+        cv.check_type('outer universe', outer, openmc.UniverseBase)
         self._outer = outer
 
     @staticmethod
@@ -92,7 +92,7 @@ class Lattice(IDManagerMixin, ABC):
             Group in HDF5 file
         universes : dict
             Dictionary mapping universe IDs to instances of
-            :class:`openmc.Universe`.
+            :class:`openmc.UniverseBase`.
 
         Returns
         -------
@@ -115,20 +115,20 @@ class Lattice(IDManagerMixin, ABC):
         -------
         universes : collections.OrderedDict
             Dictionary whose keys are universe IDs and values are
-            :class:`openmc.Universe` instances
+            :class:`openmc.UniverseBase` instances
 
         """
 
         univs = OrderedDict()
         for k in range(len(self._universes)):
             for j in range(len(self._universes[k])):
-                if isinstance(self._universes[k][j], openmc.Universe):
+                if isinstance(self._universes[k][j], openmc.UniverseBase):
                     u = self._universes[k][j]
                     univs[u._id] = u
                 else:
                     for i in range(len(self._universes[k][j])):
                         u = self._universes[k][j][i]
-                        assert isinstance(u, openmc.Universe)
+                        assert isinstance(u, openmc.UniverseBase)
                         univs[u._id] = u
 
         if self.outer is not None:
@@ -246,7 +246,7 @@ class Lattice(IDManagerMixin, ABC):
 
         Returns
         -------
-        openmc.Universe
+        openmc.UniverseBase
             Universe with given indices
 
         """
@@ -370,9 +370,9 @@ class RectLattice(Lattice):
     pitch : Iterable of float
         Pitch of the lattice in the x, y, and (if applicable) z directions in
         cm.
-    outer : openmc.Universe
+    outer : openmc.UniverseBase
         A universe to fill all space outside the lattice
-    universes : Iterable of Iterable of openmc.Universe
+    universes : Iterable of Iterable of openmc.UniverseBase
         A two- or three-dimensional list/array of universes filling each element
         of the lattice. The first dimension corresponds to the z-direction (if
         applicable), the second dimension corresponds to the y-direction, and
@@ -633,11 +633,11 @@ class RectLattice(Lattice):
 
         cv.check_value('strategy', strategy, ('degenerate', 'lns'))
         cv.check_type('universes_to_ignore', universes_to_ignore, Iterable,
-                      openmc.Universe)
+                      openmc.UniverseBase)
         cv.check_type('materials_to_clone', materials_to_clone, Iterable,
                       openmc.Material)
         cv.check_type('lattice_neighbors', lattice_neighbors, Iterable,
-                      openmc.Universe)
+                      openmc.UniverseBase)
         cv.check_value('number of lattice_neighbors', len(lattice_neighbors),
                        (0, 8))
         cv.check_type('key', key, types.FunctionType)
@@ -973,7 +973,7 @@ class RectLattice(Lattice):
             Group in HDF5 file
         universes : dict
             Dictionary mapping universe IDs to instances of
-            :class:`openmc.Universe`.
+            :class:`openmc.UniverseBase`.
 
         Returns
         -------
@@ -999,7 +999,7 @@ class RectLattice(Lattice):
             lattice.outer = universes[outer]
 
         # Build array of Universe pointers for the Lattice
-        uarray = np.empty(universe_ids.shape, dtype=openmc.Universe)
+        uarray = np.empty(universe_ids.shape, dtype=openmc.UniverseBase)
 
         for z in range(universe_ids.shape[0]):
             for y in range(universe_ids.shape[1]):
@@ -1054,9 +1054,9 @@ class HexLattice(Lattice):
         Pitch of the lattice in cm. The first item in the iterable specifies the
         pitch in the radial direction and, if the lattice is 3D, the second item
         in the iterable specifies the pitch in the axial direction.
-    outer : openmc.Universe
+    outer : openmc.UniverseBase
         A universe to fill all space outside the lattice
-    universes : Nested Iterable of openmc.Universe
+    universes : Nested Iterable of openmc.UniverseBase
         A two- or three-dimensional list/array of universes filling each element
         of the lattice. Each sub-list corresponds to one ring of universes and
         should be ordered from outermost ring to innermost ring. The universes
@@ -1173,7 +1173,7 @@ class HexLattice(Lattice):
 
     @property
     def ndim(self):
-        return 2 if isinstance(self.universes[0][0], openmc.Universe) else 3
+        return 2 if isinstance(self.universes[0][0], openmc.UniverseBase) else 3
 
     @center.setter
     def center(self, center):
@@ -1196,7 +1196,7 @@ class HexLattice(Lattice):
 
     @Lattice.universes.setter
     def universes(self, universes):
-        cv.check_iterable_type('lattice universes', universes, openmc.Universe,
+        cv.check_iterable_type('lattice universes', universes, openmc.UniverseBase,
                                min_depth=2, max_depth=3)
         self._universes = universes
 
@@ -2052,7 +2052,7 @@ class HexLattice(Lattice):
             Group in HDF5 file
         universes : dict
             Dictionary mapping universe IDs to instances of
-            :class:`openmc.Universe`.
+            :class:`openmc.UniverseBase`.
 
         Returns
         -------

src/dagmc.cpp

@@ -587,32 +587,42 @@ std::pair<double, int32_t> DAGCell::distance(
   if (!dag_univ)
     fatal_error("DAGMC call made for particle in a non-DAGMC universe");
 
-  moab::ErrorCode rval;
+  // initialize to lost particle conditions
+  int surf_idx = -1;
+  double dist = INFINITY;
+
   moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
   moab::EntityHandle hit_surf;
-  double dist;
+
+  // create the ray
   double pnt[3] = {r.x, r.y, r.z};
   double dir[3] = {u.x, u.y, u.z};
-  rval = dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &p->history());
-  MB_CHK_ERR_CONT(rval);
-  int surf_idx;
+  MB_CHK_ERR_CONT(
+    dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &p->history()));
   if (hit_surf != 0) {
     surf_idx =
       dag_univ->surf_idx_offset_ + dagmc_ptr_->index_by_handle(hit_surf);
-  } else {
-    // indicate that particle is lost
-    surf_idx = -1;
-    dist = INFINITY;
-    if (settings::run_mode == RunMode::PLOTTING) return {dist, surf_idx};
-    if (!dagmc_ptr_->is_implicit_complement(vol) ||
-        model::universe_map[dag_univ->id_] == model::root_universe) {
-      std::string material_id = p->material() == MATERIAL_VOID
-                              ? "-1 (VOID)"
-                              : std::to_string(model::materials[p->material()]->id());
-      p->mark_as_lost(
-        fmt::format("No intersection found with DAGMC cell {}, material {}",
-          id_, material_id));
-    }
+  } else if (!dagmc_ptr_->is_implicit_complement(vol) ||
+             is_root_universe(dag_univ->id_)) {
+    // surface boundary conditions are ignored for projection plotting, meaning
+    // that the particle may move through the graveyard (bounding) volume and
+    // into the implicit complement on the other side where no intersection will
+    // be found. Treating this as a lost particle is problematic when plotting.
+    // Instead, the infinite distance and invalid surface index are returned.
+    if (settings::run_mode == RunMode::PLOTTING) return {INFTY, -1};
+
+    // the particle should be marked as lost immediately if an intersection
+    // isn't found in a volume that is not the implicit complement. In the case
+    // that the DAGMC model is the root universe of the geometry, even a missing
+    // intersection in the implicit complement should trigger this condition.
+    std::string material_id =
+      p->material() == MATERIAL_VOID
+        ? "-1 (VOID)"
+        : std::to_string(model::materials[p->material()]->id());
+    auto lost_particle_msg = fmt::format(
+      "No intersection found with DAGMC cell {}, filled with material {}", id_,
+      material_id);
+    p->mark_as_lost(lost_particle_msg);
   }
 
   return {dist, surf_idx};

src/geometry_aux.cpp

@@ -610,6 +610,11 @@ int maximum_levels(int32_t univ)
   return levels_below;
 }
 
+bool is_root_universe(int32_t univ_id)
+{
+  return model::universe_map[univ_id] == model::root_universe;
+}
+
 //==============================================================================
 
 void free_memory_geometry()

tests/unit_tests/dagmc/test_lost_particles.py

@@ -0,0 +1,81 @@
+import numpy as np
+from pathlib import Path
+
+import openmc
+import openmc.lib
+
+import pytest
+
+pytestmark = pytest.mark.skipif(
+    not openmc.lib._dagmc_enabled(),
+    reason="DAGMC CAD geometry is not enabled.")
+
+
+@pytest.fixture
+def broken_dagmc_model(request):
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    ### MATERIALS ###
+    fuel = openmc.Material(name='no-void fuel')
+    fuel.set_density('g/cc', 10.29769)
+    fuel.add_nuclide('U233', 1.0)
+
+    cladding = openmc.Material(name='clad')
+    cladding.set_density('g/cc', 6.55)
+    cladding.add_nuclide('Zr90', 1.0)
+
+    h1 = openmc.Material(name='water')
+    h1.set_density('g/cc', 0.75)
+    h1.add_nuclide('H1', 1.0)
+
+    model.materials = openmc.Materials([fuel, cladding, h1])
+
+    ### GEOMETRY ###
+    # create the DAGMC universe using a model that has many triangles
+    # removed
+    dagmc_file = Path(request.fspath).parent / "broken_model.h5m"
+    pincell_univ = openmc.DAGMCUniverse(filename=dagmc_file, auto_geom_ids=True)
+
+    # create a 2 x 2 lattice using the DAGMC pincell
+    pitch = np.asarray((24.0, 24.0))
+    lattice = openmc.RectLattice()
+    lattice.pitch = pitch
+    lattice.universes = [[pincell_univ] * 2] * 2
+    lattice.lower_left = -pitch
+
+    # clip the DAGMC geometry at +/- 10 cm w/ CSG planes
+    rpp = openmc.model.RectangularParallelepiped(
+        -pitch[0], pitch[0], -pitch[1], pitch[1], -10.0, 10.0, boundary_type='reflective')
+    bounding_cell = openmc.Cell(fill=lattice, region=-rpp)
+
+    model.geometry = openmc.Geometry(root=[bounding_cell])
+
+    # settings
+    model.settings.particles = 100
+    model.settings.batches = 10
+    model.settings.inactive = 2
+    model.settings.output = {'summary': False}
+
+    model.export_to_xml()
+
+    return model
+
+
+def test_lost_particles(run_in_tmpdir, broken_dagmc_model):
+    broken_dagmc_model.export_to_xml()
+    # ensure that particles will be lost when cell intersections can't be found
+    # due to the removed triangles in this model
+    with pytest.raises(RuntimeError, match='Maximum number of lost particles has been reached.'):
+        openmc.run()
+
+    # run this again, but with the dagmc universe as the root unvierse
+    for univ in broken_dagmc_model.geometry.get_all_universes().values():
+        if isinstance(univ, openmc.DAGMCUniverse):
+            broken_dagmc_model.geometry.root_unvierse = univ
+            break
+
+    broken_dagmc_model.export_to_xml()
+    with pytest.raises(RuntimeError, match='Maximum number of lost particles has been reached.'):
+        openmc.run()
+

tools/ci/gha-install-vectfit.sh

@@ -39,7 +39,10 @@ cd $HOME
 git clone -b $XTENSOR_BLAS_BRANCH $XTENSOR_BLAS_REPO
 cd xtensor-blas && mkdir build && cd build && cmake .. && sudo make install
 
+# Install wheel (remove when vectfit supports installation with build isolation)
+pip install wheel
+
 # Install vectfit
 cd $HOME
 git clone https://github.com/liangjg/vectfit.git
-pip install ./vectfit
+pip install --no-build-isolation ./vectfit