Reevaluate cell in a close collision after a surface crossing

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.

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();

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)

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
 )
 

tests/cpp_unit_tests/test_geometry.cpp

@@ -0,0 +1,127 @@
+#include <catch2/catch_test_macros.hpp>
+
+#include "openmc/cell.h"
+#include "openmc/geometry.h"
+#include "openmc/particle_data.h"
+#include "openmc/surface.h"
+#include "openmc/universe.h"
+
+#include <memory>
+#include <string>
+
+#include <pugixml.hpp>
+
+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<openmc::SurfaceSphere>(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<openmc::Universe>();
+    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<openmc::Universe>();
+    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<openmc::CSGCell> 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<openmc::CSGCell>(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);
+}