From 51b06c473cb1923cc1ae3ed7d9e176e0fd62b3b2 Mon Sep 17 00:00:00 2001 From: GuySten Date: Sun, 3 May 2026 11:05:07 +0300 Subject: [PATCH 1/2] recalc cell in a collision after a surface crossing --- src/particle.cpp | 13 +++++++++++++ tests/unit_tests/test_geometry.py | 27 +++++++++++++++++++++++++++ 2 files changed, 40 insertions(+) diff --git a/src/particle.cpp b/src/particle.cpp index 46df63cd13a..2977d359819 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -296,6 +296,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() @@ -374,6 +378,8 @@ void Particle::event_collide() if (!model::active_meshsurf_tallies.empty()) score_meshsurface_tally(*this, model::active_meshsurf_tallies); + auto last_surface = std::abs(surface()); + // Clear surface component surface() = SURFACE_NONE; @@ -383,6 +389,13 @@ void Particle::event_collide() collision_mg(*this); } + if (last_surface != SURFACE_NONE) { + const auto& surf {*model::surfaces[last_surface - 1].get()}; + Direction normal = surf.normal(r()); + if (normal.dot(u()) * normal.dot(u_last()) < 0.0) + neighbor_list_find_cell(*this); + } + // Collision track feature to recording particle interaction if (settings::collision_track) { collision_track_record(*this); diff --git a/tests/unit_tests/test_geometry.py b/tests/unit_tests/test_geometry.py index 6cc577c820c..9269f157e3b 100644 --- a/tests/unit_tests/test_geometry.py +++ b/tests/unit_tests/test_geometry.py @@ -403,3 +403,30 @@ def test_redundant_surfaces(): geom = openmc.Geometry([c3]) redundant_surfs = geom.remove_redundant_surfaces() assert len(redundant_surfs) == 0 + + +def test_sphere_overlap(run_in_tmpdir): + model = openmc.model.Model() + shell_material = openmc.Material() + shell_material.add_element("Au", 1.0, percent_type="ao") + shell_material.set_density("g/cm3", 1930) + # revolver density + model.materials = openmc.Materials([shell_material]) + surface_inner_shell = openmc.Sphere(r=0.0035) + surface_outer_shell = openmc.Sphere(r=0.0095) + sphere_surface_detector_1 = openmc.Sphere(r=20000, boundary_type="vacuum") + fuel_region = -surface_inner_shell + shell_region = +surface_inner_shell & -surface_outer_shell + void_region_1 = +surface_outer_shell & -sphere_surface_detector_1 + fuel_cell = openmc.Cell(region=fuel_region) + shell_cell = openmc.Cell(region=shell_region, fill=shell_material) + void_cell_1 = openmc.Cell(region=void_region_1) + model.geometry = openmc.Geometry([fuel_cell, shell_cell, void_cell_1]) + source = openmc.Source() + source.angle = openmc.stats.Isotropic() + model.settings = openmc.Settings() + model.settings.batches = 100 + model.settings.particles = 80000 + model.settings.run_mode = "fixed source" + model.settings.source = source + model.run(geometry_debug=True) From 4dfb153ff97a2289fe94a97e928b49de1f41811d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 3 Aug 2026 10:33:37 -0500 Subject: [PATCH 2/2] Add reconcile_cell_after_collision function --- include/openmc/geometry.h | 15 ++- src/geometry.cpp | 33 +++++++ src/particle.cpp | 18 ++-- tests/cpp_unit_tests/CMakeLists.txt | 1 + tests/cpp_unit_tests/test_geometry.cpp | 127 +++++++++++++++++++++++++ tests/unit_tests/test_geometry.py | 27 ------ 6 files changed, 182 insertions(+), 39 deletions(-) create mode 100644 tests/cpp_unit_tests/test_geometry.cpp diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index 107cc7d1f3e..8f261675810 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -61,8 +61,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 ddb61385f18..5bb34dca648 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -281,6 +281,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 2977d359819..5eedaf53726 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -378,9 +378,11 @@ void Particle::event_collide() if (!model::active_meshsurf_tallies.empty()) score_meshsurface_tally(*this, model::active_meshsurf_tallies); - auto last_surface = std::abs(surface()); - - // 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) { @@ -389,13 +391,6 @@ void Particle::event_collide() collision_mg(*this); } - if (last_surface != SURFACE_NONE) { - const auto& surf {*model::surfaces[last_surface - 1].get()}; - Direction normal = surf.normal(r()); - if (normal.dot(u()) * normal.dot(u_last()) < 0.0) - neighbor_list_find_cell(*this); - } - // Collision track feature to recording particle interaction if (settings::collision_track) { collision_track_record(*this); @@ -460,6 +455,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() diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt index ce7e539ea57..57e39186241 100644 --- a/tests/cpp_unit_tests/CMakeLists.txt +++ b/tests/cpp_unit_tests/CMakeLists.txt @@ -8,6 +8,7 @@ set(TEST_NAMES test_mesh 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); +} diff --git a/tests/unit_tests/test_geometry.py b/tests/unit_tests/test_geometry.py index 9269f157e3b..6cc577c820c 100644 --- a/tests/unit_tests/test_geometry.py +++ b/tests/unit_tests/test_geometry.py @@ -403,30 +403,3 @@ def test_redundant_surfaces(): geom = openmc.Geometry([c3]) redundant_surfs = geom.remove_redundant_surfaces() assert len(redundant_surfs) == 0 - - -def test_sphere_overlap(run_in_tmpdir): - model = openmc.model.Model() - shell_material = openmc.Material() - shell_material.add_element("Au", 1.0, percent_type="ao") - shell_material.set_density("g/cm3", 1930) - # revolver density - model.materials = openmc.Materials([shell_material]) - surface_inner_shell = openmc.Sphere(r=0.0035) - surface_outer_shell = openmc.Sphere(r=0.0095) - sphere_surface_detector_1 = openmc.Sphere(r=20000, boundary_type="vacuum") - fuel_region = -surface_inner_shell - shell_region = +surface_inner_shell & -surface_outer_shell - void_region_1 = +surface_outer_shell & -sphere_surface_detector_1 - fuel_cell = openmc.Cell(region=fuel_region) - shell_cell = openmc.Cell(region=shell_region, fill=shell_material) - void_cell_1 = openmc.Cell(region=void_region_1) - model.geometry = openmc.Geometry([fuel_cell, shell_cell, void_cell_1]) - source = openmc.Source() - source.angle = openmc.stats.Isotropic() - model.settings = openmc.Settings() - model.settings.batches = 100 - model.settings.particles = 80000 - model.settings.run_mode = "fixed source" - model.settings.source = source - model.run(geometry_debug=True)