diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index 43cf9bb585f..314e87dafe9 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -95,8 +95,19 @@ int cell_instance_at_level(const GeometryState& p, int level); //! valid geometry coordinate stack. //============================================================================== bool exhaustive_find_cell(GeometryState& p, bool verbose = false); -bool neighbor_list_find_cell( - GeometryState& p, bool verbose = false); // Only usable on surface crossings + +//============================================================================== +//! Locate a particle starting from its current coordinate level. +//============================================================================== +bool neighbor_list_find_cell(GeometryState& p, bool verbose = false); + +//============================================================================== +//! Reconcile the current cell after a direction change near a surface. +//! +//! \param p A particle whose coordinate stack may need to be updated. +//============================================================================== + +void reconcile_cell_after_collision(GeometryState& p); //============================================================================== //! Move a particle into a new lattice tile. diff --git a/src/geometry.cpp b/src/geometry.cpp index bf654f4f9b2..687be57fee9 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -306,6 +306,39 @@ bool neighbor_list_find_cell(GeometryState& p, bool verbose) return found; } +void reconcile_cell_after_collision(GeometryState& p) +{ + // Find the first coordinate level whose current cell is inconsistent with + // the particle's post-collision direction. + int invalid_level = C_NONE; + for (int level = 0; level < p.n_coord(); ++level) { + const auto& coord {p.coord(level)}; + if (coord.cell() == C_NONE || !model::cells[coord.cell()]->contains( + coord.r(), coord.u(), SURFACE_NONE)) { + invalid_level = level; + break; + } + } + + if (invalid_level == C_NONE) + return; + + // Search from the first inconsistent level so that parent cells and + // transformed lower universes are both handled correctly. + if (p.coord(invalid_level).cell() != C_NONE) { + p.n_coord() = invalid_level + 1; + if (neighbor_list_find_cell(p)) + return; + } + + // The current cell may not have a complete neighbor list yet. Fall back to + // the normal exhaustive search used after a failed surface crossing. + p.n_coord() = 1; + if (!exhaustive_find_cell(p)) { + p.mark_as_lost("Could not find particle after a collision near a surface."); + } +} + bool exhaustive_find_cell(GeometryState& p, bool verbose) { int i_universe = p.lowest_coord().universe(); diff --git a/src/particle.cpp b/src/particle.cpp index 3c16bbe1a5a..30ccb903bfc 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -323,6 +323,10 @@ void Particle::event_advance() if (distance == distance_cutoff) { wgt() = 0.0; } + + // Clear surface component if distance is long enough + if (distance > TINY_BIT) + surface() = SURFACE_NONE; } void Particle::event_cross_surface() @@ -402,7 +406,11 @@ void Particle::event_collide() if (!model::active_meshsurf_tallies.empty()) score_meshsurface_tally(*this, model::active_meshsurf_tallies); - // Clear surface component + // Preserve whether the particle is still associated with a recently crossed + // surface so that a direction change during a near-surface collision can be + // reconciled afterward. The surface marker is no longer needed during the + // collision itself. + const bool near_surface = surface() != SURFACE_NONE; surface() = SURFACE_NONE; if (settings::run_CE) { @@ -475,6 +483,9 @@ void Particle::event_collide() #ifdef OPENMC_DAGMC_ENABLED history().reset(); #endif + + if (near_surface && alive()) + reconcile_cell_after_collision(*this); } void Particle::event_revive_from_secondary(const SourceSite& site) diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt index d4e45b5ac14..9f01f0c5416 100644 --- a/tests/cpp_unit_tests/CMakeLists.txt +++ b/tests/cpp_unit_tests/CMakeLists.txt @@ -10,6 +10,7 @@ set(TEST_NAMES test_ray test_region test_tensor + test_geometry # Add additional unit test files here ) diff --git a/tests/cpp_unit_tests/test_geometry.cpp b/tests/cpp_unit_tests/test_geometry.cpp new file mode 100644 index 00000000000..083ea00989d --- /dev/null +++ b/tests/cpp_unit_tests/test_geometry.cpp @@ -0,0 +1,127 @@ +#include + +#include "openmc/cell.h" +#include "openmc/geometry.h" +#include "openmc/particle_data.h" +#include "openmc/surface.h" +#include "openmc/universe.h" + +#include +#include + +#include + +namespace { + +class GeometryFixture { +public: + GeometryFixture() + : root_universe_ {openmc::model::root_universe}, + n_coord_levels_ {openmc::model::n_coord_levels} + { + openmc::model::cells.clear(); + openmc::model::cell_map.clear(); + openmc::model::universes.clear(); + openmc::model::universe_map.clear(); + openmc::model::surfaces.clear(); + openmc::model::surface_map.clear(); + + openmc::model::n_coord_levels = 2; + openmc::model::root_universe = 0; + + pugi::xml_document surface_doc; + auto surface_node = surface_doc.append_child("surface"); + surface_node.append_attribute("id") = 1; + surface_node.append_attribute("type") = "sphere"; + surface_node.append_attribute("coeffs") = "0 0 0 1"; + openmc::model::surfaces.push_back( + std::make_unique(surface_node)); + openmc::model::surface_map[1] = 0; + + openmc::model::cells.push_back(make_cell(1, 0, 1, "")); + openmc::model::cells.push_back(make_cell(2, 1, -1, "-1")); + openmc::model::cells.push_back(make_cell(3, 1, -1, "+1")); + for (int i = 0; i < openmc::model::cells.size(); ++i) + openmc::model::cell_map[openmc::model::cells[i]->id_] = i; + + auto root = std::make_unique(); + root->id_ = 0; + root->cells_ = {0}; + root->n_instances_ = 1; + openmc::model::universes.push_back(std::move(root)); + openmc::model::universe_map[0] = 0; + + auto nested = std::make_unique(); + nested->id_ = 1; + nested->cells_ = {1, 2}; + nested->n_instances_ = 1; + openmc::model::universes.push_back(std::move(nested)); + openmc::model::universe_map[1] = 1; + } + + ~GeometryFixture() + { + openmc::model::cells.clear(); + openmc::model::cell_map.clear(); + openmc::model::universes.clear(); + openmc::model::universe_map.clear(); + openmc::model::surfaces.clear(); + openmc::model::surface_map.clear(); + openmc::model::root_universe = root_universe_; + openmc::model::n_coord_levels = n_coord_levels_; + } + +private: + static std::unique_ptr make_cell( + int id, int universe, int fill, const char* region) + { + pugi::xml_document doc; + auto node = doc.append_child("cell"); + node.append_attribute("id") = id; + node.append_attribute("universe") = universe; + if (fill >= 0) { + const auto fill_value {std::to_string(fill)}; + node.append_child("fill").text() = fill_value.c_str(); + } else { + node.append_child("material").text() = "void"; + } + if (region[0] != '\0') + node.append_child("region").text() = region; + + auto cell = std::make_unique(node); + if (fill < 0) { + cell->type_ = openmc::Fill::MATERIAL; + cell->sqrtkT_.push_back(0.0); + cell->density_mult_.push_back(1.0); + } else { + cell->type_ = openmc::Fill::UNIVERSE; + } + return cell; + } + + int root_universe_; + int n_coord_levels_; +}; + +} // namespace + +TEST_CASE("Reconcile a particle after a collision near a surface") +{ + GeometryFixture fixture; + + openmc::GeometryState p; + p.n_coord() = 2; + p.coord(0).universe() = 0; + p.coord(0).cell() = 0; + p.coord(0).r() = {0.0, 0.0, 0.0}; + p.coord(0).u() = {-1.0, 0.0, 0.0}; + p.coord(1).universe() = 1; + p.coord(1).cell() = 2; + p.coord(1).r() = {1.0 - 1.0e-13, 0.0, 0.0}; + p.coord(1).u() = {-1.0, 0.0, 0.0}; + + openmc::reconcile_cell_after_collision(p); + REQUIRE(p.n_coord() == 2); + REQUIRE(p.coord(0).cell() == 0); + REQUIRE(p.coord(1).cell() == 1); +}