Fix bug in cylindrical and spherical meshes when grid boundary is coincident with cell boundary

src/mesh.cpp

@@ -500,7 +500,8 @@ void StructuredMesh::raytrace_mesh(
   }
 
   // translate start and end positions,
-  // this needs to come after the get_indices call because it does its own translation
+  // this needs to come after the get_indices call because it does its own
+  // translation
   local_coords(r0);
   local_coords(r1);
 
@@ -1075,15 +1076,16 @@ double CylindricalMesh::find_r_crossing(
   const double inv_denominator = 1.0 / denominator;
 
   const double p = (u.x * r.x + u.y * r.y) * inv_denominator;
-  double D = p * p + (r0 * r0 - r.x * r.x - r.y * r.y) * inv_denominator;
+  double c = r.x * r.x + r.y * r.y - r0 * r0;
+  double D = p * p - c * inv_denominator;
 
   if (D < 0.0)
     return INFTY;
 
   D = std::sqrt(D);
 
   // the solution -p - D is always smaller as -p + D : Check this one first
-  if (-p - D > l)
+  if (-p - D > l && std::abs(c) > 1e-10)
     return -p - D;
   if (-p + D > l)
     return -p + D;
@@ -1303,12 +1305,13 @@ double SphericalMesh::find_r_crossing(
   // |r+s*u| = |r| + 2*s*r*u + s^2 (|u|==1 !)
   const double r0 = grid_[0][shell];
   const double p = r.dot(u);
-  double D = p * p - r.dot(r) + r0 * r0;
+  double c = r.dot(r) - r0 * r0;
+  double D = p * p - c;
 
   if (D >= 0.0) {
     D = std::sqrt(D);
     // the solution -p - D is always smaller as -p + D : Check this one first
-    if (-p - D > l)
+    if (-p - D > l && std::abs(c) > 1e-10)
       return -p - D;
     if (-p + D > l)
       return -p + D;

tests/unit_tests/test_filter_mesh.py

@@ -1,6 +1,7 @@
 import math
 
 import numpy as np
+import pytest
 from uncertainties import unumpy
 
 import openmc
@@ -112,6 +113,105 @@ def test_cylindrical_mesh_estimators(run_in_tmpdir):
     std_dev = unumpy.std_devs(delta)
     assert np.all(diff < 3*std_dev)
 
+
+@pytest.mark.parametrize("scale", [1.0, 1e2, 1e4, 1e5])
+def test_cylindrical_mesh_coincident(scale, run_in_tmpdir):
+    """Test for cylindrical mesh boundary being coincident with a cell boundary"""
+
+    fuel = openmc.Material()
+    fuel.add_nuclide('U235', 1.)
+    fuel.set_density('g/cm3', 4.5)
+
+    zcyl = openmc.ZCylinder(r=1.25*scale)
+    box = openmc.rectangular_prism(4*scale, 4*scale, boundary_type='reflective')
+    cell1 = openmc.Cell(fill=fuel, region=-zcyl)
+    cell2 = openmc.Cell(fill=None, region=+zcyl & box)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell1, cell2])
+
+    model.settings.particles = 100
+    model.settings.batches = 10
+    model.settings.inactive = 0
+
+    cyl_mesh = openmc.CylindricalMesh()
+    cyl_mesh.r_grid = [0., 1.25*scale]
+    cyl_mesh.phi_grid = [0., 2*math.pi]
+    cyl_mesh.z_grid = [-1e10, 1e10]
+    cyl_mesh_filter = openmc.MeshFilter(cyl_mesh)
+    cell_filter = openmc.CellFilter([cell1])
+
+    tally1 = openmc.Tally()
+    tally1.filters = [cyl_mesh_filter]
+    tally1.scores = ['flux']
+    tally2 = openmc.Tally()
+    tally2.filters = [cell_filter]
+    tally2.scores = ['flux']
+    model.tallies = openmc.Tallies([tally1, tally2])
+
+    # Run OpenMC
+    sp_filename = model.run()
+
+    # Get flux for each of the two tallies
+    with openmc.StatePoint(sp_filename) as sp:
+        t1 = sp.tallies[tally1.id]
+        t2 = sp.tallies[tally2.id]
+        mean1 = t1.mean.ravel()[0]
+        mean2 = t2.mean.ravel()[0]
+
+    # The two tallies should be exactly the same
+    assert mean1 == pytest.approx(mean2)
+
+
+@pytest.mark.parametrize("scale", [1.0, 1e2, 1e4, 1e5])
+def test_spherical_mesh_coincident(scale, run_in_tmpdir):
+    """Test for spherical mesh boundary being coincident with a cell boundary"""
+
+    fuel = openmc.Material()
+    fuel.add_nuclide('U235', 1.)
+    fuel.set_density('g/cm3', 4.5)
+
+    sph = openmc.Sphere(r=1.25*scale)
+    rcc = openmc.model.RectangularParallelepiped(
+        -2*scale, 2*scale, -2*scale, 2*scale, -2*scale, 2*scale,
+        boundary_type='reflective')
+    cell1 = openmc.Cell(fill=fuel, region=-sph)
+    cell2 = openmc.Cell(fill=None, region=+sph & -rcc)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell1, cell2])
+
+    model.settings.particles = 100
+    model.settings.batches = 10
+    model.settings.inactive = 0
+
+    sph_mesh = openmc.SphericalMesh()
+    sph_mesh.r_grid = [0., 1.25*scale]
+    sph_mesh.phi_grid = [0., 2*math.pi]
+    sph_mesh.theta_grid = [0., math.pi]
+    sph_mesh_filter = openmc.MeshFilter(sph_mesh)
+    cell_filter = openmc.CellFilter([cell1])
+
+    tally1 = openmc.Tally()
+    tally1.filters = [sph_mesh_filter]
+    tally1.scores = ['flux']
+    tally2 = openmc.Tally()
+    tally2.filters = [cell_filter]
+    tally2.scores = ['flux']
+    model.tallies = openmc.Tallies([tally1, tally2])
+
+    # Run OpenMC
+    sp_filename = model.run()
+
+    # Get flux for each of the two tallies
+    with openmc.StatePoint(sp_filename) as sp:
+        t1 = sp.tallies[tally1.id]
+        t2 = sp.tallies[tally2.id]
+        mean1 = t1.mean.ravel()[0]
+        mean2 = t2.mean.ravel()[0]
+
+    # The two tallies should be exactly the same
+    assert mean1 == pytest.approx(mean2)
+
+
 def test_get_reshaped_data(run_in_tmpdir):
     """Test that expanding MeshFilter dimensions works as expected"""