Merge pull request #2053 from paulromano/torus-bugfix

Improve robustness of torus distance calculation when particle is coincident

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;

tests/regression_tests/torus/large_major/inputs_true.dat

@@ -0,0 +1,37 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+  <cell id="1" material="void" region="-1" universe="1" />
+  <cell id="2" material="1" region="1 -2" universe="1" />
+  <cell id="3" material="2" region="2 -3" universe="1" />
+  <surface coeffs="0.0 0.0 0.0 1000.0 30.0 30.0" id="1" type="z-torus" />
+  <surface coeffs="0.0 0.0 0.0 1000.0 35.0 35.0" id="2" type="z-torus" />
+  <surface boundary="vacuum" coeffs="0.0 0.0 0.0 1000.0 40.0 40.0" id="3" type="z-torus" />
+</geometry>
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+  <material id="1">
+    <density units="g/cm3" value="1.0" />
+    <nuclide ao="1.0" name="W184" />
+  </material>
+  <material id="2">
+    <density units="g/cm3" value="5.0" />
+    <nuclide ao="1.0" name="Fe56" />
+  </material>
+</materials>
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+  <run_mode>fixed source</run_mode>
+  <particles>1000</particles>
+  <batches>10</batches>
+  <source strength="1.0">
+    <space type="point">
+      <parameters>-1000.0 0 0</parameters>
+    </space>
+  </source>
+</settings>
+<?xml version='1.0' encoding='utf-8'?>
+<tallies>
+  <tally id="1">
+    <scores>flux</scores>
+  </tally>
+</tallies>

tests/regression_tests/torus/large_major/results_true.dat

@@ -0,0 +1,3 @@
+tally 1:
+9.675396E+02
+9.363406E+04

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