From 1cae2f712a947e863fa87328d92e9fc668dd1b94 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 1 Sep 2023 08:52:51 -0500
Subject: [PATCH] Prevent segfault in distance to boundary calculation (#2659)

---
 src/geometry.cpp                             |  3 +-
 tests/unit_tests/test_no_visible_boundary.py | 30 ++++++++++++++++++++
 2 files changed, 32 insertions(+), 1 deletion(-)
 create mode 100644 tests/unit_tests/test_no_visible_boundary.py

diff --git a/src/geometry.cpp b/src/geometry.cpp
index 5e015c1bc0f..47152e17295 100644
--- a/src/geometry.cpp
+++ b/src/geometry.cpp
@@ -419,13 +419,14 @@ BoundaryInfo distance_to_boundary(Particle& p)
     double& d = info.distance;
     if (d_surf < d_lat - FP_COINCIDENT) {
       if (d == INFINITY || (d - d_surf) / d >= FP_REL_PRECISION) {
+        // Update closest distance
         d = d_surf;
 
         // If the cell is not simple, it is possible that both the negative and
         // positive half-space were given in the region specification. Thus, we
         // have to explicitly check which half-space the particle would be
         // traveling into if the surface is crossed
-        if (c.is_simple()) {
+        if (c.is_simple() || d == INFTY) {
           info.surface_index = level_surf_cross;
         } else {
           Position r_hit = r + d_surf * u;
diff --git a/tests/unit_tests/test_no_visible_boundary.py b/tests/unit_tests/test_no_visible_boundary.py
new file mode 100644
index 00000000000..8f07b8da9e7
--- /dev/null
+++ b/tests/unit_tests/test_no_visible_boundary.py
@@ -0,0 +1,30 @@
+import openmc
+
+
+def test_no_visible_boundary(run_in_tmpdir):
+    copper = openmc.Material()
+    copper.add_nuclide('Cu63', 1.0)
+    copper.set_density('g/cm3', 0.3)
+    air = openmc.Material()
+    air.add_nuclide('N14', 1.0)
+    air.set_density('g/cm3', 0.0012)
+
+    # Create a simple model of a neutron source directly impinging on a thin
+    # disc of copper. Neutrons leaving the back of the disc see no surfaces in
+    # front of them.
+    disc = openmc.model.RightCircularCylinder((0., 0., 1.), 0.1, 1.2)
+    box = openmc.rectangular_prism(width=10, height=10, boundary_type='vacuum')
+    c1 = openmc.Cell(fill=copper, region=-disc)
+    c2 = openmc.Cell(fill=air, region=+disc & box)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([c1, c2])
+    model.settings.run_mode = 'fixed source'
+    model.settings.particles = 1000
+    model.settings.batches = 5
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Point(),
+        angle=openmc.stats.Monodirectional((0., 0., 1.))
+    )
+
+    # Run model to ensure it doesn't segfault
+    model.run()