Merge pull request #2443 from pshriwise/sph-cyl-mesh-fix

Improvements to spherical/cylindrical mesh radial boundary coincidence with geometry

include/openmc/constants.h

@@ -51,6 +51,10 @@ constexpr double FP_PRECISION {1e-14};
 constexpr double FP_REL_PRECISION {1e-5};
 constexpr double FP_COINCIDENT {1e-12};
 
+// Coincidence tolerances
+constexpr double TORUS_TOL {1e-10};
+constexpr double RADIAL_MESH_TOL {1e-10};
+
 // Maximum number of random samples per history
 constexpr int MAX_SAMPLE {100000};
 

src/mesh.cpp

@@ -1066,6 +1066,9 @@ double CylindricalMesh::find_r_crossing(
   // s^2 * (u^2 + v^2) + 2*s*(u*x+v*y) + x^2+y^2-r0^2 = 0
 
   const double r0 = grid_[0][shell];
+  if (r0 == 0.0)
+    return INFTY;
+
   const double denominator = u.x * u.x + u.y * u.y;
 
   // Direction of flight is in z-direction. Will never intersect r.
@@ -1085,7 +1088,10 @@ double CylindricalMesh::find_r_crossing(
   D = std::sqrt(D);
 
   // the solution -p - D is always smaller as -p + D : Check this one first
-  if (-p - D > l && std::abs(c) > 1e-10)
+  if (std::abs(c) <= RADIAL_MESH_TOL)
+    return INFTY;
+
+  if (-p - D > l)
     return -p - D;
   if (-p + D > l)
     return -p + D;
@@ -1304,14 +1310,19 @@ double SphericalMesh::find_r_crossing(
   // solve |r+s*u| = r0
   // |r+s*u| = |r| + 2*s*r*u + s^2 (|u|==1 !)
   const double r0 = grid_[0][shell];
+  if (r0 == 0.0)
+    return INFTY;
   const double p = r.dot(u);
   double c = r.dot(r) - r0 * r0;
   double D = p * p - c;
 
+  if (std::abs(c) <= RADIAL_MESH_TOL)
+    return INFTY;
+
   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 && std::abs(c) > 1e-10)
+    if (-p - D > l)
       return -p - D;
     if (-p + D > l)
       return -p + D;

src/surface.cpp

@@ -1044,7 +1044,7 @@ double torus_distance(double x1, double x2, double x3, double u1, double u2,
   // zero but possibly small and positive. A tolerance is set to discard that
   // zero.
   double distance = INFTY;
-  double cutoff = coincident ? 1e-10 : 0.0;
+  double cutoff = coincident ? TORUS_TOL : 0.0;
   for (int i = 0; i < 4; ++i) {
     if (roots[i].imag() == 0) {
       double root = roots[i].real();

tests/unit_tests/test_cylindrical_mesh.py

@@ -75,7 +75,7 @@ def label(p):
         return f'estimator:{p}'
 
 @pytest.mark.parametrize('estimator,origin', test_cases, ids=label)
-def test_offset_mesh(model, estimator, origin):
+def test_offset_mesh(run_in_tmpdir, model, estimator, origin):
     """Tests that the mesh has been moved based on tally results
     """
     mesh = model.tallies[0].filters[0].mesh
@@ -98,8 +98,99 @@ def test_offset_mesh(model, estimator, origin):
         centroids = mesh.centroids
         for ijk in mesh.indices:
             i, j, k = np.array(ijk) - 1
-            print(centroids[:, i, j, k])
             if model.geometry.find(centroids[:, i, j, k]):
                 mean[i, j, k] == 0.0
             else:
                 mean[i, j, k] != 0.0
+
+@pytest.fixture()
+def void_coincident_geom_model():
+    """A model with many geometric boundaries coincident with mesh boundaries
+       across many scales
+    """
+    openmc.reset_auto_ids()
+    model = openmc.model.Model()
+
+    model.materials = openmc.Materials()
+    radii = [0.1,1, 5, 50, 100, 150, 250]
+    cylinders = [openmc.ZCylinder(r=ri) for ri in radii]
+    cylinders[-1].boundary_type = 'vacuum'
+
+    regions = openmc.model.subdivide(cylinders)[:-1]
+    cells = [openmc.Cell(region=r, fill=None) for r in regions]
+    geom = openmc.Geometry(cells)
+
+    model.geometry = geom
+
+    settings = openmc.Settings(run_mode='fixed source')
+    settings.batches = 2
+    settings.particles = 1000
+    model.settings = settings
+
+    mesh = openmc.CylindricalMesh()
+    mesh.r_grid = np.linspace(0, 250, 501)
+    mesh.z_grid = [-250, 250]
+    mesh.phi_grid = np.linspace(0, 2*np.pi, 2)
+    mesh_filter = openmc.MeshFilter(mesh)
+
+    tally = openmc.Tally()
+    tally.scores = ['flux']
+    tally.filters = [mesh_filter]
+
+    model.tallies = openmc.Tallies([tally])
+
+    return model
+
+
+# convenience function for checking tally results
+# in the following tests
+def _check_void_cylindrical_tally(statepoint_filename):
+    with openmc.StatePoint(statepoint_filename) as sp:
+        flux_tally = sp.tallies[1]
+        mesh = flux_tally.find_filter(openmc.MeshFilter).mesh
+        neutron_flux = flux_tally.get_reshaped_data().squeeze()
+        # we expect the tally results to be the same as the mesh grid width
+        # for these cases
+        d_r = mesh.r_grid[1] - mesh.r_grid[0]
+        assert neutron_flux == pytest.approx(d_r)
+
+def test_void_geom_pnt_src(run_in_tmpdir, void_coincident_geom_model):
+    src = openmc.Source()
+    src.space = openmc.stats.Point()
+    src.angle = openmc.stats.PolarAzimuthal(mu=openmc.stats.Discrete([0.0], [1.0]))
+    src.energy = openmc.stats.Discrete([14.06e6], [1])
+    void_coincident_geom_model.settings.source = src
+
+    sp_filename = void_coincident_geom_model.run()
+    _check_void_cylindrical_tally(sp_filename)
+
+
+def test_void_geom_boundary_src(run_in_tmpdir, void_coincident_geom_model):
+    # update source to a number of points on the outside of the cylinder
+    # with directions pointing toward the origin
+    bbox = void_coincident_geom_model.geometry.bounding_box
+
+    # can't source particle directly on the geometry boundary
+    outer_r = bbox[1][0] - 1e-08
+
+    n_sources = 100
+    radial_vals = np.linspace(0.0, 2.0*np.pi, n_sources)
+
+    sources = []
+
+    energy = openmc.stats.Discrete([14.06e6], [1])
+    for val in radial_vals:
+        src = openmc.Source()
+        src.energy = energy
+
+        pnt = np.array([np.cos(val), np.sin(val), 0.0])
+        u = -pnt
+        src.space = openmc.stats.Point(outer_r*pnt)
+        src.angle = openmc.stats.Monodirectional(u)
+        src.strength = 0.5/n_sources
+        sources.append(src)
+
+    void_coincident_geom_model.settings.source = sources
+    sp_filename = void_coincident_geom_model.run()
+
+    _check_void_cylindrical_tally(sp_filename)

tests/unit_tests/test_spherical_mesh.py

@@ -79,7 +79,7 @@ def label(p):
         return f'estimator:{p}'
 
 @pytest.mark.parametrize('estimator,origin', test_cases, ids=label)
-def test_offset_mesh(model, estimator, origin):
+def test_offset_mesh(run_in_tmpdir, model, estimator, origin):
     """Tests that the mesh has been moved based on tally results
     """
     mesh = model.tallies[0].filters[0].mesh
@@ -102,8 +102,108 @@ def test_offset_mesh(model, estimator, origin):
         centroids = mesh.centroids
         for ijk in mesh.indices:
             i, j, k = np.array(ijk) - 1
-            print(centroids[:, i, j, k])
             if model.geometry.find(centroids[:, i, j, k]):
                 mean[i, j, k] == 0.0
             else:
                 mean[i, j, k] != 0.0
+
+# Some void geometry tests to check our radial intersection methods on
+# spherical and cylindrical meshes
+
+@pytest.fixture()
+def void_coincident_geom_model():
+    """A model with many geometric boundaries coincident with mesh boundaries
+    across many scales
+    """
+    openmc.reset_auto_ids()
+
+    model = openmc.Model()
+
+    model.materials = openmc.Materials()
+    radii = [0.1, 1, 5, 50, 100, 150, 250]
+    spheres = [openmc.Sphere(r=ri) for ri in radii]
+    spheres[-1].boundary_type = 'vacuum'
+
+    regions = openmc.model.subdivide(spheres)[:-1]
+    cells = [openmc.Cell(region=r, fill=None) for r in regions]
+    geom = openmc.Geometry(cells)
+
+    model.geometry = geom
+
+    settings = openmc.Settings(run_mode='fixed source')
+    settings.batches = 2
+    settings.particles = 5000
+    model.settings = settings
+
+    mesh = openmc.SphericalMesh()
+    mesh.r_grid = np.linspace(0, 250, 501)
+    mesh_filter = openmc.MeshFilter(mesh)
+
+    tally = openmc.Tally()
+    tally.scores = ['flux']
+    tally.filters = [mesh_filter]
+
+    model.tallies = openmc.Tallies([tally])
+
+    return model
+
+
+# convenience function for checking tally results
+# in the following tests
+def _check_void_spherical_tally(statepoint_filename):
+    with openmc.StatePoint(statepoint_filename) as sp:
+        flux_tally = sp.tallies[1]
+        mesh = flux_tally.find_filter(openmc.MeshFilter).mesh
+        neutron_flux = flux_tally.get_reshaped_data().squeeze()
+        # the flux values for each bin should equal the width
+        # width of the mesh bins
+        d_r = mesh.r_grid[1] - mesh.r_grid[0]
+        assert neutron_flux == pytest.approx(d_r)
+
+
+def test_void_geom_pnt_src(run_in_tmpdir, void_coincident_geom_model):
+    # add isotropic point source
+    src = openmc.Source()
+    src.space = openmc.stats.Point()
+    src.energy = openmc.stats.Discrete([14.06e6], [1])
+    void_coincident_geom_model.settings.source = src
+
+    # run model and check tally results
+    sp_filename = void_coincident_geom_model.run()
+    _check_void_spherical_tally(sp_filename)
+
+
+def test_void_geom_boundary_src(run_in_tmpdir, void_coincident_geom_model):
+    # update source to a number of points on the outside of the sphere
+    # with directions pointing toward the origin
+    n_sources = 20
+    phi_vals = np.linspace(0, np.pi, n_sources)
+    theta_vals = np.linspace(0, 2.0*np.pi, n_sources)
+
+    bbox = void_coincident_geom_model.geometry.bounding_box
+    # can't source particles directly on the geometry boundary
+    outer_r = bbox[1][0] - 1e-08
+
+    sources = []
+
+    energy = openmc.stats.Discrete([14.06e6], [1])
+
+    for phi, theta in zip(phi_vals, theta_vals):
+
+        src = openmc.Source()
+        src.energy = energy
+
+        pnt = np.array([np.sin(phi)*np.cos(theta), np.sin(phi)*np.sin(theta), np.cos(phi)])
+        u = -pnt
+        src.space = openmc.stats.Point(outer_r*pnt)
+        src.angle = openmc.stats.Monodirectional(u)
+        # set source strengths so that we can still expect
+        # a tally value of 0.5
+        src.strength = 0.5/n_sources
+
+        sources.append(src)
+
+    void_coincident_geom_model.settings.source = sources
+
+    sp_filename = void_coincident_geom_model.run()
+    _check_void_spherical_tally(sp_filename)