Determine surface-source starting half-space from particle direction (#4024)

src/particle.cpp

@@ -184,8 +184,11 @@ void Particle::from_source(const SourceSite* src)
 
   // Convert signed surface ID to signed index
   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 index_plus_one = it->second + 1;
+      surface() = (src->surf_id > 0) ? index_plus_one : -index_plus_one;
+    }
   }
 
   wgt_born() = src->wgt_born;

src/source.cpp

@@ -38,6 +38,7 @@
 #include "openmc/simulation.h"
 #include "openmc/state_point.h"
 #include "openmc/string_utils.h"
+#include "openmc/surface.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -547,7 +548,28 @@ SourceSite FileSource::sample(uint64_t* seed) const
 {
   // Sample a particle randomly from list
   size_t i_site = sites_.size() * prn(seed);
-  return sites_[i_site];
+  SourceSite site = sites_[i_site];
+
+  // Surface source files store unsigned surface IDs. If the ID refers to a CSG
+  // surface containing the source site, determine the signed half-space from
+  // the particle direction. Otherwise, ignore the surface ID and allow the
+  // normal cell search to locate the particle.
+  if (site.surf_id != SURFACE_NONE) {
+    auto it = model::surface_map.find(std::abs(site.surf_id));
+    if (it != model::surface_map.end()) {
+      const auto& surf = *model::surfaces[it->second];
+      if (surf.geom_type() == GeometryType::CSG &&
+          std::abs(surf.evaluate(site.r)) < FP_COINCIDENT) {
+        int surf_id = std::abs(site.surf_id);
+        site.surf_id =
+          (site.u.dot(surf.normal(site.r)) > 0.0) ? surf_id : -surf_id;
+        return site;
+      }
+    }
+    site.surf_id = SURFACE_NONE;
+  }
+
+  return site;
 }
 
 //==============================================================================

tests/unit_tests/test_surface_source_write.py

@@ -193,6 +193,40 @@ def test_particle_direction(parameter, run_in_tmpdir, model):
                 assert False
 
 
+def test_surface_source_id(run_in_tmpdir):
+    """Test handling of surface IDs for source particles."""
+
+    # Surface 20 is the surface matching the source in this model. Surface 11 is
+    # absent and surface 10 is reused for an unrelated boundary that does not
+    # contain the source position.
+    sources = []
+    for direction in ((1.0, 0.0, 0.0), (-1.0, 0.0, 0.0)):
+        for source_surface_id in (20, 11, 10):
+            filename = f"surface_source_{len(sources)}.h5"
+            site = openmc.SourceParticle(
+                r=(0.0, 0.0, 0.0), u=direction, surf_id=source_surface_id
+            )
+            openmc.write_source_file([site], filename)
+            sources.append(openmc.FileSource(filename))
+
+    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 = openmc.YPlane(1.0, surface_id=10, 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 = 1000
+    model.settings.batches = 1
+    model.settings.source = sources
+    model.run()
+
+
 @pytest.fixture
 def model_dagmc(request):
     """Model based on the mesh file 'dagmc.h5m' available from