Fix negative distances from bins_crossed for CylindricalMesh (#3370)

Author: pshriwise

include/openmc/mesh.h

@@ -138,9 +138,6 @@ class Mesh {
   //! Perform any preparation needed to support point location within the mesh
   virtual void prepare_for_point_location() {};
 
-  //! Update a position to the local coordinates of the mesh
-  virtual void local_coords(Position& r) const {};
-
   //! Return a position in the local coordinates of the mesh
   virtual Position local_coords(const Position& r) const { return r; };
 
@@ -427,8 +424,6 @@ class PeriodicStructuredMesh : public StructuredMesh {
   PeriodicStructuredMesh() = default;
   PeriodicStructuredMesh(pugi::xml_node node) : StructuredMesh {node} {};
 
-  void local_coords(Position& r) const override { r -= origin_; };
-
   Position local_coords(const Position& r) const override
   {
     return r - origin_;

src/mesh.cpp

@@ -915,6 +915,11 @@ void StructuredMesh::raytrace_mesh(
   if (total_distance == 0.0 && settings::solver_type != SolverType::RANDOM_RAY)
     return;
 
+  // keep a copy of the original global position to pass to get_indices,
+  // which performs its own transformation to local coordinates
+  Position global_r = r0;
+  Position local_r = local_coords(r0);
+
   const int n = n_dimension_;
 
   // Flag if position is inside the mesh
@@ -925,7 +930,7 @@ void StructuredMesh::raytrace_mesh(
 
   // Calculate index of current cell. Offset the position a tiny bit in
   // direction of flight
-  MeshIndex ijk = get_indices(r0 + TINY_BIT * u, in_mesh);
+  MeshIndex ijk = get_indices(global_r + TINY_BIT * u, in_mesh);
 
   // if track is very short, assume that it is completely inside one cell.
   // Only the current cell will score and no surfaces
@@ -936,16 +941,10 @@ void StructuredMesh::raytrace_mesh(
     return;
   }
 
-  // translate start and end positions,
-  // this needs to come after the get_indices call because it does its own
-  // translation
-  local_coords(r0);
-  local_coords(r1);
-
   // Calculate initial distances to next surfaces in all three dimensions
   std::array<MeshDistance, 3> distances;
   for (int k = 0; k < n; ++k) {
-    distances[k] = distance_to_grid_boundary(ijk, k, r0, u, 0.0);
+    distances[k] = distance_to_grid_boundary(ijk, k, local_r, u, 0.0);
   }
 
   // Loop until r = r1 is eventually reached
@@ -975,7 +974,7 @@ void StructuredMesh::raytrace_mesh(
       // The two other directions are still valid!
       ijk[k] = distances[k].next_index;
       distances[k] =
-        distance_to_grid_boundary(ijk, k, r0, u, traveled_distance);
+        distance_to_grid_boundary(ijk, k, local_r, u, traveled_distance);
 
       // Check if we have left the interior of the mesh
       in_mesh = ((ijk[k] >= 1) && (ijk[k] <= shape_[k]));
@@ -1005,10 +1004,10 @@ void StructuredMesh::raytrace_mesh(
 
       // Calculate the new cell index and update all distances to next
       // surfaces.
-      ijk = get_indices(r0 + (traveled_distance + TINY_BIT) * u, in_mesh);
+      ijk = get_indices(global_r + (traveled_distance + TINY_BIT) * u, in_mesh);
       for (int k = 0; k < n; ++k) {
         distances[k] =
-          distance_to_grid_boundary(ijk, k, r0, u, traveled_distance);
+          distance_to_grid_boundary(ijk, k, local_r, u, traveled_distance);
       }
 
       // If inside the mesh, Tally inward current
@@ -1482,7 +1481,7 @@ std::string CylindricalMesh::get_mesh_type() const
 StructuredMesh::MeshIndex CylindricalMesh::get_indices(
   Position r, bool& in_mesh) const
 {
-  local_coords(r);
+  r = local_coords(r);
 
   Position mapped_r;
   mapped_r[0] = std::hypot(r.x, r.y);
@@ -1630,23 +1629,21 @@ StructuredMesh::MeshDistance CylindricalMesh::distance_to_grid_boundary(
   const MeshIndex& ijk, int i, const Position& r0, const Direction& u,
   double l) const
 {
-  Position r = r0 - origin_;
-
   if (i == 0) {
 
     return std::min(
-      MeshDistance(ijk[i] + 1, true, find_r_crossing(r, u, l, ijk[i])),
-      MeshDistance(ijk[i] - 1, false, find_r_crossing(r, u, l, ijk[i] - 1)));
+      MeshDistance(ijk[i] + 1, true, find_r_crossing(r0, u, l, ijk[i])),
+      MeshDistance(ijk[i] - 1, false, find_r_crossing(r0, u, l, ijk[i] - 1)));
 
   } else if (i == 1) {
 
     return std::min(MeshDistance(sanitize_phi(ijk[i] + 1), true,
-                      find_phi_crossing(r, u, l, ijk[i])),
+                      find_phi_crossing(r0, u, l, ijk[i])),
       MeshDistance(sanitize_phi(ijk[i] - 1), false,
-        find_phi_crossing(r, u, l, ijk[i] - 1)));
+        find_phi_crossing(r0, u, l, ijk[i] - 1)));
 
   } else {
-    return find_z_crossing(r, u, l, ijk[i]);
+    return find_z_crossing(r0, u, l, ijk[i]);
   }
 }
 
@@ -1762,7 +1759,7 @@ std::string SphericalMesh::get_mesh_type() const
 StructuredMesh::MeshIndex SphericalMesh::get_indices(
   Position r, bool& in_mesh) const
 {
-  local_coords(r);
+  r = local_coords(r);
 
   Position mapped_r;
   mapped_r[0] = r.norm();

tests/regression_tests/filter_mesh/inputs_true.dat

@@ -54,8 +54,8 @@
     <mesh id="5" type="cylindrical">
       <r_grid>0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5</r_grid>
       <phi_grid>0.0 0.3490658503988659 0.6981317007977318 1.0471975511965976 1.3962634015954636 1.7453292519943295 2.0943951023931953 2.443460952792061 2.792526803190927 3.141592653589793 3.490658503988659 3.839724354387525 4.1887902047863905 4.537856055185257 4.886921905584122 5.235987755982989 5.585053606381854 5.93411945678072 6.283185307179586</phi_grid>
-      <z_grid>-7.5 -6.5625 -5.625 -4.6875 -3.75 -2.8125 -1.875 -0.9375 0.0 0.9375 1.875 2.8125 3.75 4.6875 5.625 6.5625 7.5</z_grid>
-      <origin>0.0 0.0 0.0</origin>
+      <z_grid>0.0 0.9375 1.875 2.8125 3.75 4.6875 5.625 6.5625 7.5 8.4375 9.375 10.3125 11.25 12.1875 13.125 14.0625 15.0</z_grid>
+      <origin>0.0 0.0 -7.5</origin>
     </mesh>
     <mesh id="6" type="spherical">
       <r_grid>0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5</r_grid>

tests/regression_tests/filter_mesh/test.py

@@ -61,9 +61,10 @@ def model():
     np.testing.assert_allclose(recti_mesh.volumes, recti_mesh_exp_vols)
 
     cyl_mesh = openmc.CylindricalMesh(
+        origin=(0, 0, -7.5),
         r_grid=np.linspace(0, 7.5, 18),
         phi_grid=np.linspace(0, 2*pi, 19),
-        z_grid=np.linspace(-7.5, 7.5, 17),
+        z_grid=np.linspace(0, 15, 17),
     )
     dr = 0.5 * np.diff(np.linspace(0, 7.5, 18)**2)
     dp = np.full(cyl_mesh.dimension[1], 2*pi / 18)

tests/unit_tests/test_cylindrical_mesh.py

@@ -49,6 +49,7 @@ def model():
 
     return openmc.Model(geometry=geom, settings=settings, tallies=tallies)
 
+
 def test_origin_read_write_to_xml(run_in_tmpdir, model):
     """Tests that the origin attribute can be written and read back to XML
     """
@@ -63,8 +64,9 @@ def test_origin_read_write_to_xml(run_in_tmpdir, model):
     np.testing.assert_equal(new_mesh.origin, mesh.origin)
 
 estimators = ('tracklength', 'collision')
-origins = set(permutations((-geom_size, 0, 0)))
-origins |= set(permutations((geom_size, 0, 0)))
+offset = geom_size + 0.001
+origins = set(permutations((-offset , 0, 0)))
+origins |= set(permutations((offset, 0, 0)))
 
 test_cases = product(estimators, origins)
 
@@ -74,8 +76,9 @@ def label(p):
     if isinstance(p, str):
         return f'estimator:{p}'
 
+
 @pytest.mark.parametrize('estimator,origin', test_cases, ids=label)
-def test_offset_mesh(run_in_tmpdir, model, estimator, origin):
+def test_offset_mesh(model, estimator, origin):
     """Tests that the mesh has been moved based on tally results
     """
     mesh = model.tallies[0].filters[0].mesh
@@ -103,6 +106,7 @@ def test_offset_mesh(run_in_tmpdir, model, estimator, origin):
             else:
                 mean[i, j, k] != 0.0
 
+
 @pytest.fixture()
 def void_coincident_geom_model():
     """A model with many geometric boundaries coincident with mesh boundaries
@@ -155,6 +159,7 @@ def _check_void_cylindrical_tally(statepoint_filename):
         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.IndependentSource()
     src.space = openmc.stats.Point()

tests/unit_tests/test_spherical_mesh.py

@@ -63,8 +63,6 @@ def test_origin_read_write_to_xml(run_in_tmpdir, model):
     np.testing.assert_equal(new_mesh.origin, mesh.origin)
 
 estimators = ('tracklength', 'collision')
-# TODO: determine why this is needed for spherical mesh
-# but not cylindrical mesh
 offset = geom_size + 0.001
 
 origins = set(permutations((-offset, 0, 0)))