From 3ccd57841f03ae8a613b4bef68778ee325fa49ac Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 3 Oct 2026 01:00:18 +0000 Subject: [PATCH 1/2] Stop searching complex cells for a boundary beyond the collision site The distance to collision is now sampled before the distance to the nearest boundary and passed to the boundary search as a maximum distance. The search through a complex cell, which moves from one surface crossing to the next until one changes whether the ray is in the cell, stops once it passes the collision site, since a farther boundary is not used. The boundary search uses no random numbers, so results are unchanged. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014RN7JroEkrbkLKVHbMDap9 --- include/openmc/cell.h | 20 +++++++++++++------- include/openmc/dagmc.h | 3 ++- include/openmc/geometry.h | 5 ++++- src/cell.cpp | 13 ++++++++++--- src/dagmc.cpp | 4 ++-- src/geometry.cpp | 4 ++-- src/particle.cpp | 9 ++++++--- tests/cpp_unit_tests/test_region.cpp | 17 +++++++++++++++++ 8 files changed, 56 insertions(+), 19 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index c2a09d774b5..e59f3f65d8a 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -81,8 +81,10 @@ class Region { bool contains(Position r, Direction u, int32_t on_surface) const; //! Find the oncoming boundary of this cell. - std::pair distance( - Position r, Direction u, int32_t on_surface) const; + //! \param max_distance Distance beyond which the boundary is not needed. If + //! the boundary is not nearer, the distance returned may be INFTY. + std::pair distance(Position r, Direction u, + int32_t on_surface, double max_distance = INFTY) const; //! Get the BoundingBox for this cell. BoundingBox bounding_box(int32_t cell_id) const; @@ -126,7 +128,7 @@ class Region { //! Find the oncoming boundary of this cell for a complex cell. std::pair distance_complex( - Position r, Direction u, int32_t on_surface) const; + Position r, Direction u, int32_t on_surface, double max_distance) const; //! BoundingBox if the particle is in a simple cell. BoundingBox bounding_box_simple() const; @@ -216,8 +218,11 @@ class Cell { virtual bool contains(Position r, Direction u, int32_t on_surface) const = 0; //! Find the oncoming boundary of this cell. - virtual std::pair distance( - Position r, Direction u, int32_t on_surface, GeometryState* p) const = 0; + //! \param max_distance Distance beyond which the boundary is not needed. If + //! the boundary is not nearer, the distance returned may be INFTY. + virtual std::pair distance(Position r, Direction u, + int32_t on_surface, GeometryState* p, + double max_distance = INFTY) const = 0; //! Write all information needed to reconstruct the cell to an HDF5 group. //! \param group_id An HDF5 group id. @@ -437,9 +442,10 @@ class CSGCell : public Cell { int n_surfaces() const override { return region_.n_surfaces(); } std::pair distance(Position r, Direction u, - int32_t on_surface, GeometryState* p) const override + int32_t on_surface, GeometryState* p, + double max_distance = INFTY) const override { - return region_.distance(r, u, on_surface); + return region_.distance(r, u, on_surface, max_distance); } bool contains(Position r, Direction u, int32_t on_surface) const override diff --git a/include/openmc/dagmc.h b/include/openmc/dagmc.h index 1031f5fbcb9..41a5afbee00 100644 --- a/include/openmc/dagmc.h +++ b/include/openmc/dagmc.h @@ -69,7 +69,8 @@ class DAGCell : public Cell { bool contains(Position r, Direction u, int32_t on_surface) const override; std::pair distance(Position r, Direction u, - int32_t on_surface, GeometryState* p) const override; + int32_t on_surface, GeometryState* p, + double max_distance = INFTY) const override; BoundingBox bounding_box() const override; diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index e8504d48261..b9dfb6c4a98 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -120,7 +120,10 @@ void cross_lattice( //! Find the next boundary a particle will intersect. //============================================================================== -BoundaryInfo distance_to_boundary(GeometryState& p); +//! \param max_distance Distance beyond which a boundary is not needed. The +//! distance to a farther boundary may be returned as INFTY. +BoundaryInfo distance_to_boundary( + GeometryState& p, double max_distance = INFTY); } // namespace openmc diff --git a/src/cell.cpp b/src/cell.cpp index 8be8fb8d48f..66c22a3af79 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -948,12 +948,12 @@ std::string Region::str() const //============================================================================== std::pair Region::distance( - Position r, Direction u, int32_t on_surface) const + Position r, Direction u, int32_t on_surface, double max_distance) const { if (simple_) { return distance_to_nearest_surface(r, u, on_surface, false); } else { - return distance_complex(r, u, on_surface); + return distance_complex(r, u, on_surface, max_distance); } } @@ -997,7 +997,7 @@ std::pair Region::distance_to_nearest_surface(Position r, //============================================================================== std::pair Region::distance_complex( - Position r, Direction u, int32_t on_surface) const + Position r, Direction u, int32_t on_surface, double max_distance) const { const bool in_region = contains_complex(r, u, on_surface); double total_distance {0.0}; @@ -1015,6 +1015,13 @@ std::pair Region::distance_complex( // the wrong side of a curved surface. r += distance * u; total_distance += distance; + + // Stop searching once the boundary is known to be beyond the distance of + // interest + if (total_distance >= max_distance) { + return {INFTY, std::numeric_limits::max()}; + } + i_surf = std::abs(i_surf); const auto& surf {*model::surfaces[i_surf - 1]}; if (u.dot(surf.normal(r)) <= 0.0) { diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 9fce4df4aef..751d8390845 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -821,8 +821,8 @@ void DAGUniverse::override_assign_material(std::unique_ptr& c, DAGCell::DAGCell(std::shared_ptr dag_ptr, int32_t dag_idx) : Cell {}, dagmc_ptr_(dag_ptr), dag_index_(dag_idx) {}; -std::pair DAGCell::distance( - Position r, Direction u, int32_t on_surface, GeometryState* p) const +std::pair DAGCell::distance(Position r, Direction u, + int32_t on_surface, GeometryState* p, double /*max_distance*/) const { // if we've changed direction or we're not on a surface, // reset the history and update last direction diff --git a/src/geometry.cpp b/src/geometry.cpp index ecacf0bffbf..55fdb08964e 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -426,7 +426,7 @@ void cross_lattice(GeometryState& p, const BoundaryInfo& boundary, bool verbose) //============================================================================== -BoundaryInfo distance_to_boundary(GeometryState& p) +BoundaryInfo distance_to_boundary(GeometryState& p, double max_distance) { BoundaryInfo info; double d_lat = INFINITY; @@ -442,7 +442,7 @@ BoundaryInfo distance_to_boundary(GeometryState& p) Cell& c {*model::cells[coord.cell()]}; // Find the oncoming surface in this cell and the distance to it. - auto surface_distance = c.distance(r, u, p.surface(), &p); + auto surface_distance = c.distance(r, u, p.surface(), &p, max_distance); d_surf = surface_distance.first; level_surf_cross = surface_distance.second; diff --git a/src/particle.cpp b/src/particle.cpp index 998b71883a8..ed4a6d03e49 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -276,9 +276,6 @@ void Particle::event_calculate_xs() void Particle::event_advance() { - // Find the distance to the nearest boundary - boundary() = distance_to_boundary(*this); - // Sample a distance to collision if (type() == ParticleType::electron() || type() == ParticleType::positron()) { @@ -294,6 +291,12 @@ void Particle::event_advance() double distance_cutoff = (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY; + // Find the distance to the nearest boundary. A boundary beyond the + // collision site is not needed, which lets the search through complex cells + // stop early. The search uses no random numbers, so results do not depend + // on the order in which the two distances are found. + boundary() = distance_to_boundary(*this, collision_distance()); + // Select smaller of the three distances double distance = std::min({boundary().distance(), collision_distance(), distance_cutoff}); diff --git a/tests/cpp_unit_tests/test_region.cpp b/tests/cpp_unit_tests/test_region.cpp index ba784c842c7..5d9499446da 100644 --- a/tests/cpp_unit_tests/test_region.cpp +++ b/tests/cpp_unit_tests/test_region.cpp @@ -183,6 +183,23 @@ TEST_CASE("Find boundary after virtual surface crossings") REQUIRE(surface == 2); } + SECTION("Boundary beyond a maximum distance") + { + // A boundary nearer than the maximum distance is found as before + auto [distance, surface] = + region.distance({1.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, 0, 5.5); + REQUIRE(distance == Catch::Approx(5.0)); + REQUIRE(surface == 2); + + // The search stops once it passes the maximum distance, whether before or + // after the virtual crossings + for (double max_distance : {2.0, 4.5}) { + auto [far_distance, far_surface] = + region.distance({1.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, 0, max_distance); + REQUIRE(far_distance == openmc::INFTY); + } + } + SECTION("Starting outside the region") { // Along -x from x=7, entering the sphere at x=6 is the first boundary. From 1ac28b4608527704f42c9268a7849a3a1d932d0c Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 3 Oct 2026 01:14:47 +0000 Subject: [PATCH 2/2] Skip the boundary search for particles that collide in place Electrons and positrons whose energy is deposited locally have a distance to collision of zero, so the distance to the nearest boundary is never used. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014RN7JroEkrbkLKVHbMDap9 --- src/particle.cpp | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index ed4a6d03e49..343d4a0d7b4 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -293,9 +293,15 @@ void Particle::event_advance() // Find the distance to the nearest boundary. A boundary beyond the // collision site is not needed, which lets the search through complex cells - // stop early. The search uses no random numbers, so results do not depend - // on the order in which the two distances are found. - boundary() = distance_to_boundary(*this, collision_distance()); + // stop early, and no boundary is needed at all for a particle that collides + // where it is, such as an electron or positron whose energy is deposited + // locally. The search uses no random numbers, so results do not depend on + // the order in which the two distances are found. + if (collision_distance() > 0.0) { + boundary() = distance_to_boundary(*this, collision_distance()); + } else { + boundary().reset(); + } // Select smaller of the three distances double distance =