diff --git a/src/surface.cpp b/src/surface.cpp
index dbfa0bf3dc..c0b3f0a593 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -1026,9 +1026,11 @@ double torus_distance(double x1, double x2, double x3, double u1, double u2,
double c1p = 2 * four_A2 * (u1 * x1 + u2 * x2);
double c0p = four_A2 * (x1 * x1 + x2 * x2);
- // Coefficient for equation: a t^4 + b t^3 + c t^2 + d t + e = 0
+ // Coefficient for equation: a t^4 + b t^3 + c t^2 + d t + e = 0. If the point
+ // is coincident, the 'e' coefficient should be zero. Explicitly setting it to
+ // zero helps avoid numerical issues below with root finding.
double coeff[5];
- coeff[0] = c0 * c0 - c0p;
+ coeff[0] = coincident ? 0.0 : c0 * c0 - c0p;
coeff[1] = 2 * c0 * c1 - c1p;
coeff[2] = c1 * c1 + 2 * c0 * c2 - c2p;
coeff[3] = 2 * c1 * c2;
diff --git a/tests/regression_tests/torus/large_major/__init__.py b/tests/regression_tests/torus/large_major/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/tests/regression_tests/torus/large_major/inputs_true.dat b/tests/regression_tests/torus/large_major/inputs_true.dat
new file mode 100644
index 0000000000..ae3f72d0fc
--- /dev/null
+++ b/tests/regression_tests/torus/large_major/inputs_true.dat
@@ -0,0 +1,37 @@
+
+
+ |
+ |
+ |
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ fixed source
+ 1000
+ 10
+
+
+ -1000.0 0 0
+
+
+
+
+
+
+ flux
+
+
diff --git a/tests/regression_tests/torus/large_major/results_true.dat b/tests/regression_tests/torus/large_major/results_true.dat
new file mode 100644
index 0000000000..cb9f54970c
--- /dev/null
+++ b/tests/regression_tests/torus/large_major/results_true.dat
@@ -0,0 +1,3 @@
+tally 1:
+9.675396E+02
+9.363406E+04
diff --git a/tests/regression_tests/torus/large_major/test.py b/tests/regression_tests/torus/large_major/test.py
new file mode 100644
index 0000000000..256716f4e1
--- /dev/null
+++ b/tests/regression_tests/torus/large_major/test.py
@@ -0,0 +1,41 @@
+import openmc
+import pytest
+
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.fixture
+def model():
+ model = openmc.Model()
+ tungsten = openmc.Material()
+ tungsten.set_density('g/cm3', 1.0)
+ tungsten.add_nuclide('W184', 1.0)
+ ss = openmc.Material()
+ ss.set_density('g/cm3', 5.0)
+ ss.add_nuclide('Fe56', 1.0)
+ model.materials.extend([tungsten, ss])
+
+ # Create nested torii with very large major radii
+ R = 1000.0
+ vacuum = openmc.ZTorus(a=R, b=30.0, c=30.0)
+ first_wall = openmc.ZTorus(a=R, b=35.0, c=35.0)
+ vessel = openmc.ZTorus(a=R, b=40.0, c=40.0, boundary_type='vacuum')
+ cell1 = openmc.Cell(region=-vacuum)
+ cell2 = openmc.Cell(fill=tungsten, region=+vacuum & -first_wall)
+ cell3 = openmc.Cell(fill=ss, region=+first_wall & -vessel)
+ model.geometry = openmc.Geometry([cell1, cell2, cell3])
+
+ model.settings.run_mode ='fixed source'
+ model.settings.particles = 1000
+ model.settings.batches = 10
+ model.settings.source = openmc.Source(space=openmc.stats.Point((-R, 0, 0,)))
+
+ tally = openmc.Tally()
+ tally.scores = ['flux']
+ model.tallies.append(tally)
+ return model
+
+
+def test_torus_large_major(model):
+ harness = PyAPITestHarness('statepoint.10.h5', model)
+ harness.main()