Merge pull request #4 from pshriwise/centre_for_cylinder_spherical_meshes

Centre for cylinder spherical meshes

include/openmc/mesh.h

@@ -241,6 +241,7 @@ public:
   xt::xtensor<double, 1> lower_left_;  //!< Lower-left coordinates of mesh
   xt::xtensor<double, 1> upper_right_; //!< Upper-right coordinates of mesh
   std::array<int, 3> shape_; //!< Number of mesh elements in each dimension
+  Position origin_ {0.0, 0.0, 0.0};
 
 protected:
 };
@@ -359,7 +360,6 @@ public:
   void to_hdf5(hid_t group) const override;
 
   array<vector<double>, 3> grid_;
-  Position origin_;
 
   int set_grid();
 
@@ -414,7 +414,6 @@ public:
   void to_hdf5(hid_t group) const override;
 
   array<vector<double>, 3> grid_;
-  Position origin_;
 
   int set_grid();
 

openmc/mesh.py

@@ -1119,6 +1119,7 @@ class CylindricalMesh(StructuredMesh):
         fmt = '{0: <16}{1}{2}\n'
         string = super().__repr__()
         string += fmt.format('\tDimensions', '=\t', self.n_dimension)
+        string += fmt.format('\tOrigin', '=\t', self.origin)
         r_grid_str = str(self._r_grid) if self._r_grid is None else len(self._r_grid)
         string += fmt.format('\tN R pnts:', '=\t', r_grid_str)
         if self._r_grid is not None:
@@ -1445,6 +1446,7 @@ class SphericalMesh(StructuredMesh):
         fmt = '{0: <16}{1}{2}\n'
         string = super().__repr__()
         string += fmt.format('\tDimensions', '=\t', self.n_dimension)
+        string += fmt.format('\tOrigin', '=\t', self.origin)
         r_grid_str = str(self._r_grid) if self._r_grid is None else len(self._r_grid)
         string += fmt.format('\tN R pnts:', '=\t', r_grid_str)
         if self._r_grid is not None:

src/mesh.cpp

@@ -412,7 +412,6 @@ template<class T>
 void StructuredMesh::raytrace_mesh(
   Position r0, Position r1, const Direction& u, T tally) const
 {
-
   // TODO: when c++-17 is available, use "if constexpr ()" to compile-time
   // enable/disable tally calls for now, T template type needs to provide both
   // surface and track methods, which might be empty. modern optimizing
@@ -445,6 +444,11 @@ 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
+  r0 -= origin_;
+  r1 -= origin_;
+
   // Calculate initial distances to next surfaces in all three dimensions
   std::array<MeshDistance, 3> distances;
   for (int k = 0; k < n; ++k) {
@@ -972,13 +976,12 @@ std::string CylindricalMesh::get_mesh_type() const
 StructuredMesh::MeshIndex CylindricalMesh::get_indices(
   Position r, bool& in_mesh) const
 {
-  Position mapped_r;
+  r -= origin_;
 
+  Position mapped_r;
   mapped_r[0] = std::hypot(r.x, r.y);
   mapped_r[2] = r[2];
 
-  mapped_r -= origin_;
-
   if (mapped_r[0] < FP_PRECISION) {
     mapped_r[1] = 0.0;
   } else {
@@ -1090,22 +1093,23 @@ 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(r0, u, l, ijk[i])),
-      MeshDistance(ijk[i] - 1, false, find_r_crossing(r0, u, l, ijk[i] - 1)));
+      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)));
 
   } else if (i == 1) {
 
     return std::min(MeshDistance(sanitize_phi(ijk[i] + 1), true,
-                      find_phi_crossing(r0, u, l, ijk[i])),
+                      find_phi_crossing(r, u, l, ijk[i])),
       MeshDistance(sanitize_phi(ijk[i] - 1), false,
-        find_phi_crossing(r0, u, l, ijk[i] - 1)));
+        find_phi_crossing(r, u, l, ijk[i] - 1)));
 
   } else {
-    return find_z_crossing(r0, u, l, ijk[i]);
+    return find_z_crossing(r, u, l, ijk[i]);
   }
 }
 
@@ -1209,12 +1213,11 @@ std::string SphericalMesh::get_mesh_type() const
 StructuredMesh::MeshIndex SphericalMesh::get_indices(
   Position r, bool& in_mesh) const
 {
+  r -= origin_;
+
   Position mapped_r;
-
   mapped_r[0] = r.norm();
 
-  mapped_r -= origin_;
-
   if (mapped_r[0] < FP_PRECISION) {
     mapped_r[1] = 0.0;
     mapped_r[2] = 0.0;

tests/unit_tests/test_cylindrical_mesh.py

@@ -1,16 +1,42 @@
+from itertools import product, permutations
+
 import openmc
 import numpy as np
 
-def test_origin_read_write_to_xml():
-    """Tests that the origin attribute can be written and read back to XML
-    """
+import pytest
+
+geom_size = 5
+
+@pytest.fixture()
+def model():
+    openmc.reset_auto_ids()
+
+    water = openmc.Material(name='water')
+    water.add_element('H', 2.0)
+    water.add_element('O', 1.0)
+    water.set_density('g/cc', 1.0)
+
+    rpp = openmc.model.RectangularParallelepiped(*([-geom_size, geom_size] * 3),
+                                                 boundary_type='vacuum')
+
+    cell = openmc.Cell(region=-rpp, fill=water)
+
+    geom = openmc.Geometry([cell])
+
+    source = openmc.Source()
+    source.space = openmc.stats.Point()
+    source.energy = openmc.stats.Discrete([10000], [1.0])
+
+    settings = openmc.Settings()
+    settings.particles = 1000
+    settings.batches = 10
+    settings.run_mode = 'fixed source'
+
     # build
     mesh = openmc.CylindricalMesh()
-    mesh.phi_grid = [1, 2, 3]
-    mesh.z_grid = [1, 2, 3]
-    mesh.r_grid = [1, 2, 3]
-
-    mesh.origin = [0.1, 0.2, 0.3]
+    mesh.phi_grid = np.linspace(0, 2*np.pi, 21)
+    mesh.z_grid = np.linspace(-geom_size, geom_size, 11)
+    mesh.r_grid = np.linspace(0, geom_size, geom_size)
 
     tally = openmc.Tally()
 
@@ -21,10 +47,52 @@ def test_origin_read_write_to_xml():
 
     tallies = openmc.Tallies([tally])
 
-    tallies.export_to_xml()
+    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
+    """
+    mesh = model.tallies[0].filters[0].mesh
+    mesh.origin = [0.1, 0.2, 0.3]
+    model.tallies.export_to_xml()
 
     # read back
     new_tallies = openmc.Tallies.from_xml()
     new_tally = new_tallies[0]
     new_mesh = new_tally.filters[0].mesh
-    assert np.allclose(new_mesh.origin, mesh.origin)
+    np.testing.assert_equal(new_mesh.origin, mesh.origin)
+
+estimators = ('tracklength', 'collision')
+origins = permutations([-geom_size, 0, 0])
+
+test_cases = product(estimators, origins)
+
+def label(p):
+    if isinstance(p, tuple):
+        return f'origin:{p}'
+    if isinstance(p, str):
+        return f'estimator:{p}'
+
+@pytest.mark.parametrize('estimator,origin', test_cases, ids=label)
+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
+    model.tallies[0].estimator = estimator
+    # move the center of the cylinder mesh upwards
+    mesh.origin = origin
+
+    sp_filename = model.run()
+
+    with openmc.StatePoint(sp_filename) as sp:
+        tally = sp.tallies[1]
+
+        # we've translated half of the cylinder mesh above the model,
+        # so ensure that half of the bins are populated
+        assert np.count_nonzero(tally.mean) == tally.mean.size / 2
+
+        # check that the lower half of the mesh contains a tally result
+        # and that the upper half is zero
+        mean = tally.mean.reshape(mesh.dimension[::-1]).T
+        # assert np.count_nonzero(mean[:, :, :5]) == mean.size / 2
+        # assert np.count_nonzero(mean[:, :, 5:]) == 0