diff --git a/include/openmc/geometry_aux.h b/include/openmc/geometry_aux.h index cf62debc0..0c870c9fb 100644 --- a/include/openmc/geometry_aux.h +++ b/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. //============================================================================== diff --git a/openmc/lattice.py b/openmc/lattice.py index 464daaf7e..b067b47a3 100644 --- a/openmc/lattice.py +++ b/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 ------- diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 04c8424aa..4d63f0892 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -587,32 +587,42 @@ std::pair 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}; diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index eabdc9f70..3ab5fcc2d 100644 --- a/src/geometry_aux.cpp +++ b/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() diff --git a/tests/unit_tests/dagmc/broken_model.h5m b/tests/unit_tests/dagmc/broken_model.h5m new file mode 100644 index 000000000..4bc6c7d42 Binary files /dev/null and b/tests/unit_tests/dagmc/broken_model.h5m differ diff --git a/tests/unit_tests/dagmc/test_lost_particles.py b/tests/unit_tests/dagmc/test_lost_particles.py new file mode 100644 index 000000000..502bd795e --- /dev/null +++ b/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() + diff --git a/tools/ci/gha-install-vectfit.sh b/tools/ci/gha-install-vectfit.sh index 5884c6fac..8444c3036 100755 --- a/tools/ci/gha-install-vectfit.sh +++ b/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