From 728d4191e9824553f69e1fcb4d9bc4f707eaa390 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 20 Jul 2026 23:07:09 -0500
Subject: [PATCH] Fix handling of surface source IDs

---
 src/particle.cpp                              | 17 ++++++--
 tests/unit_tests/test_surface_source_write.py | 40 +++++++++++++++++++
 2 files changed, 54 insertions(+), 3 deletions(-)

diff --git a/src/particle.cpp b/src/particle.cpp
index 10e8f213dd..3dfcd233ef 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -182,10 +182,21 @@ void Particle::from_source(const SourceSite* src)
   parent_nuclide() = src->parent_nuclide;
   delayed_group() = src->delayed_group;
 
-  // Convert signed surface ID to signed index
+  // If the source has a surface ID assigned, attempt to determine the signed
+  // surface index by comparing the particle's direction with the surface
+  // normal. This only works if the particle is on the surface.
   if (src->surf_id != SURFACE_NONE) {
-    int index_plus_one = model::surface_map[std::abs(src->surf_id)] + 1;
-    surface() = (src->surf_id > 0) ? index_plus_one : -index_plus_one;
+    auto it = model::surface_map.find(std::abs(src->surf_id));
+    if (it != model::surface_map.end()) {
+      int i_surface = it->second;
+      const auto& surf = *model::surfaces[i_surface];
+      if (surf.geom_type() == GeometryType::CSG &&
+          std::abs(surf.evaluate(r())) < FP_COINCIDENT) {
+        int index_plus_one = i_surface + 1;
+        surface() =
+          (u().dot(surf.normal(r())) > 0.0) ? index_plus_one : -index_plus_one;
+      }
+    }
   }
 
   wgt_born() = src->wgt_born;
diff --git a/tests/unit_tests/test_surface_source_write.py b/tests/unit_tests/test_surface_source_write.py
index 6f18d32b71..f8457bd0bd 100644
--- a/tests/unit_tests/test_surface_source_write.py
+++ b/tests/unit_tests/test_surface_source_write.py
@@ -193,6 +193,46 @@ def test_particle_direction(parameter, run_in_tmpdir, model):
                 assert False
 
 
+@pytest.mark.parametrize(
+    "direction",
+    [(1.0, 0.0, 0.0), (-1.0, 0.0, 0.0)],
+)
+@pytest.mark.parametrize(
+    "source_surface_id, reuse_surface_id",
+    [(20, False), (10, False), (10, True)],
+    ids=["matching", "missing", "reused"],
+)
+def test_source_surface_id(
+    direction, source_surface_id, reuse_surface_id, run_in_tmpdir
+):
+    """A source surface is used only when it matches the model surface."""
+    source = openmc.SourceParticle(
+        r=(0.0, 0.0, 0.0), u=direction, surf_id=source_surface_id
+    )
+    openmc.write_source_file([source], "surface_source.h5")
+
+    # Surface 20 is the source surface in this model. Surface 10 is either
+    # absent or reused for an unrelated boundary that does not contain the
+    # source position.
+    left = openmc.XPlane(-1.0, surface_id=1, boundary_type="vacuum")
+    middle = openmc.XPlane(0.0, surface_id=20)
+    right = openmc.XPlane(1.0, surface_id=30, boundary_type="vacuum")
+    bottom = openmc.YPlane(-1.0, surface_id=40, boundary_type="vacuum")
+    top_id = 10 if reuse_surface_id else 50
+    top = openmc.YPlane(1.0, surface_id=top_id, boundary_type="vacuum")
+    cells = [
+        openmc.Cell(region=+left & -middle & +bottom & -top),
+        openmc.Cell(region=+middle & -right & +bottom & -top),
+    ]
+
+    model = openmc.Model(geometry=openmc.Geometry(cells))
+    model.settings.run_mode = "fixed source"
+    model.settings.particles = 10
+    model.settings.batches = 1
+    model.settings.source = openmc.FileSource("surface_source.h5")
+    model.run()
+
+
 @pytest.fixture
 def model_dagmc(request):
     """Model based on the mesh file 'dagmc.h5m' available from