Merge pull request #2256 from RemDelaporteMathurin/centre_for_cylinder_spherical_meshes

Centre for cylinder & spherical meshes

docs/source/io_formats/statepoint.rst

@@ -72,12 +72,23 @@ The current version of the statepoint file format is 17.0.
 
 :Datasets: - **type** (*char[]*) -- Type of mesh.
            - **dimension** (*int*) -- Number of mesh cells in each dimension.
-           - **lower_left** (*double[]*) -- Coordinates of lower-left corner of
-             mesh.
-           - **upper_right** (*double[]*) -- Coordinates of upper-right corner
-             of mesh.
-           - **width** (*double[]*) -- Width of each mesh cell in each
-             dimension.
+           - **Regular Mesh Only:**
+              - **lower_left** (*double[]*) -- Coordinates of lower-left corner of
+                mesh.
+              - **upper_right** (*double[]*) -- Coordinates of upper-right corner
+                of mesh.
+              - **width** (*double[]*) -- Width of each mesh cell in each
+                dimension.
+           - **Rectilinear Mesh Only:**
+              - **x_grid** (*double[]*) -- Mesh divisions along the x-axis.
+              - **y_grid** (*double[]*) -- Mesh divisions along the y-axis.
+              - **z_grid** (*double[]*) -- Mesh divisions along the z-axis.
+           - **Cylindrical & Spherical Mesh Only:**
+              - **r_grid** (*double[]*) -- The mesh divisions along the r-axis.
+              - **phi_grid** (*double[]*) -- The mesh divisions along the phi-axis.
+              - **origin** (*double[]*) -- The origin in cartesian coordinates.
+           - **Spherical Mesh Only:**
+              - **theta_grid** (*double[]*) -- The mesh divisions along the theta-axis.
            - **Unstructured Mesh Only:**
               - **filename** (*char[]*) -- Name of the mesh file.
               - **library** (*char[]*) -- Mesh library used to represent the

docs/source/io_formats/tallies.rst

@@ -351,6 +351,9 @@ attributes/sub-elements:
   :theta_grid:
     The mesh divisions along the theta-axis. (For spherical mesh only.)
 
+  :origin:
+    The origin in cartesian coordinates. (For cylindrical and spherical meshes only.)
+
   :library:
     The mesh library used to represent an unstructured mesh. This can be either
     "moab" or "libmesh". (For unstructured mesh only.)

include/openmc/constants.h

@@ -25,7 +25,7 @@ using double_4dvec = vector<vector<vector<vector<double>>>>;
 constexpr int HDF5_VERSION[] {3, 0};
 
 // Version numbers for binary files
-constexpr array<int, 2> VERSION_STATEPOINT {18, 0};
+constexpr array<int, 2> VERSION_STATEPOINT {18, 1};
 constexpr array<int, 2> VERSION_PARTICLE_RESTART {2, 0};
 constexpr array<int, 2> VERSION_TRACK {3, 0};
 constexpr array<int, 2> VERSION_SUMMARY {6, 0};

include/openmc/mesh.h

@@ -75,6 +75,12 @@ public:
 
   // Methods
 
+  //! 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; };
+
   //! Sample a mesh volume using a certain seed
   //
   //! \param[in] seed Seed to use for random sampling
@@ -262,6 +268,23 @@ public:
 protected:
 };
 
+class PeriodicStructuredMesh : public StructuredMesh {
+
+public:
+  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_;
+  };
+
+  // Data members
+  Position origin_ {0.0, 0.0, 0.0}; //!< Origin of the mesh
+};
+
 //==============================================================================
 //! Tessellation of n-dimensional Euclidean space by congruent squares or cubes
 //==============================================================================
@@ -352,7 +375,7 @@ public:
   int set_grid();
 };
 
-class CylindricalMesh : public StructuredMesh {
+class CylindricalMesh : public PeriodicStructuredMesh {
 public:
   // Constructors
   CylindricalMesh() = default;
@@ -406,7 +429,7 @@ private:
   }
 };
 
-class SphericalMesh : public StructuredMesh {
+class SphericalMesh : public PeriodicStructuredMesh {
 public:
   // Constructors
   SphericalMesh() = default;

include/openmc/xml_interface.h

@@ -9,6 +9,7 @@
 #include "xtensor/xadapt.hpp"
 #include "xtensor/xarray.hpp"
 
+#include "openmc/position.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -49,5 +50,8 @@ xt::xarray<T> get_node_xarray(
   return xt::adapt(v, shape);
 }
 
+Position get_node_position(
+  pugi::xml_node node, const char* name, bool lowercase = false);
+
 } // namespace openmc
 #endif // OPENMC_XML_INTERFACE_H

openmc/mesh.py

@@ -1050,6 +1050,9 @@ class CylindricalMesh(StructuredMesh):
         The default value is [0, 2π], i.e. the full phi range.
     z_grid : numpy.ndarray
         1-D array of mesh boundary points along the z-axis.
+    origin : numpy.ndarray
+        1-D array of length 3 the (x,y,z) origin of the mesh in
+        cartesian coordinates
     indices : Iterable of tuple
         An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
         (2, 1, 1), ...]
@@ -1062,6 +1065,7 @@ class CylindricalMesh(StructuredMesh):
         self._r_grid = None
         self._phi_grid = [0.0, 2*pi]
         self._z_grid = None
+        self.origin = (0., 0., 0.)
 
     @property
     def dimension(self):
@@ -1073,6 +1077,10 @@ class CylindricalMesh(StructuredMesh):
     def n_dimension(self):
         return 3
 
+    @property
+    def origin(self):
+        return self._origin
+
     @property
     def r_grid(self):
         return self._r_grid
@@ -1099,6 +1107,12 @@ class CylindricalMesh(StructuredMesh):
                 for p in range(1, np + 1)
                 for r in range(1, nr + 1))
 
+    @origin.setter
+    def origin(self, coords):
+        cv.check_type('mesh origin', coords, Iterable, Real)
+        cv.check_length("mesh origin", coords, 3)
+        self._origin = np.asarray(coords)
+
     @r_grid.setter
     def r_grid(self, grid):
         cv.check_type('mesh r_grid', grid, Iterable, Real)
@@ -1118,6 +1132,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:
@@ -1144,6 +1159,8 @@ class CylindricalMesh(StructuredMesh):
         mesh.r_grid = group['r_grid'][()]
         mesh.phi_grid = group['phi_grid'][()]
         mesh.z_grid = group['z_grid'][()]
+        if 'origin' in group:
+            mesh.origin = group['origin'][()]
 
         return mesh
 
@@ -1240,6 +1257,9 @@ class CylindricalMesh(StructuredMesh):
         subelement = ET.SubElement(element, "z_grid")
         subelement.text = ' '.join(map(str, self.z_grid))
 
+        subelement = ET.SubElement(element, "origin")
+        subelement.text = ' '.join(map(str, self.origin))
+
         return element
 
     @classmethod
@@ -1263,6 +1283,8 @@ class CylindricalMesh(StructuredMesh):
         mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()]
         mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()]
         mesh.z_grid = [float(x) for x in get_text(elem, "z_grid").split()]
+        mesh.origin = [float(x) for x in get_text(elem, "origin", default=[0., 0., 0.]).split()]
+
         return mesh
 
     @property
@@ -1304,12 +1326,13 @@ class CylindricalMesh(StructuredMesh):
         """
         # create points
         pts_cylindrical = self.vertices.T.reshape(-1, 3)
-
         pts_cartesian = np.copy(pts_cylindrical)
-        r, phi = pts_cylindrical[:, 0], pts_cylindrical[:, 1]
-        pts_cartesian[:, 0] = r * np.cos(phi)
-        pts_cartesian[:, 1] = r * np.sin(phi)
 
+        r, phi = pts_cylindrical[:, 0], pts_cylindrical[:, 1]
+        pts_cartesian[:, 0] = r * np.cos(phi) + self.origin[0]
+        pts_cartesian[:, 1] = r * np.sin(phi) + self.origin[1]
+        pts_cartesian[:, 2] += self.origin[2]
+        
         return super().write_data_to_vtk(
             points=pts_cartesian,
             filename=filename,
@@ -1346,6 +1369,9 @@ class SphericalMesh(StructuredMesh):
     phi_grid : numpy.ndarray
         1-D array of mesh boundary points along the phi-axis in radians.
         The default value is [0, 2π], i.e. the full phi range.
+    origin : numpy.ndarray
+        1-D array of length 3 the (x,y,z) origin of the mesh in
+        cartesian coordinates
     indices : Iterable of tuple
         An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
         (2, 1, 1), ...]
@@ -1358,6 +1384,7 @@ class SphericalMesh(StructuredMesh):
         self._r_grid = None
         self._theta_grid = [0, pi]
         self._phi_grid = [0, 2*pi]
+        self.origin = (0., 0., 0.)
 
     @property
     def dimension(self):
@@ -1369,6 +1396,10 @@ class SphericalMesh(StructuredMesh):
     def n_dimension(self):
         return 3
 
+    @property
+    def origin(self):
+        return self._origin
+
     @property
     def r_grid(self):
         return self._r_grid
@@ -1395,6 +1426,12 @@ class SphericalMesh(StructuredMesh):
                 for t in range(1, nt + 1)
                 for r in range(1, nr + 1))
 
+    @origin.setter
+    def origin(self, coords):
+        cv.check_type('mesh origin', coords, Iterable, Real)
+        cv.check_length("mesh origin", coords, 3)
+        self._origin = np.asarray(coords)
+
     @r_grid.setter
     def r_grid(self, grid):
         cv.check_type('mesh r_grid', grid, Iterable, Real)
@@ -1414,6 +1451,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:
@@ -1440,6 +1478,8 @@ class SphericalMesh(StructuredMesh):
         mesh.r_grid = group['r_grid'][()]
         mesh.theta_grid = group['theta_grid'][()]
         mesh.phi_grid = group['phi_grid'][()]
+        if 'origin' in group:
+            mesh.origin = group['origin'][()]
 
         return mesh
 
@@ -1466,6 +1506,9 @@ class SphericalMesh(StructuredMesh):
         subelement = ET.SubElement(element, "phi_grid")
         subelement.text = ' '.join(map(str, self.phi_grid))
 
+        subelement = ET.SubElement(element, "origin")
+        subelement.text = ' '.join(map(str, self.origin))
+
         return element
 
     @classmethod
@@ -1489,6 +1532,8 @@ class SphericalMesh(StructuredMesh):
         mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()]
         mesh.theta_grid = [float(x) for x in get_text(elem, "theta_grid").split()]
         mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()]
+        mesh.origin = [float(x) for x in get_text(elem, "origin", default=[0., 0., 0.]).split()]
+
         return mesh
 
     @property
@@ -1530,12 +1575,13 @@ class SphericalMesh(StructuredMesh):
         """
         # create points
         pts_spherical = self.vertices.T.reshape(-1, 3)
-
         pts_cartesian = np.copy(pts_spherical)
+
         r, theta, phi = pts_spherical[:, 0], pts_spherical[:, 1], pts_spherical[:, 2]
-        pts_cartesian[:, 0] = r * np.sin(theta) * np.cos(phi)
-        pts_cartesian[:, 1] = r * np.sin(theta) * np.sin(phi)
-        pts_cartesian[:, 2] = r * np.cos(theta)
+
+        pts_cartesian[:, 0] = r * np.sin(phi) * np.cos(theta) + self.origin[0]
+        pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta) + self.origin[1]
+        pts_cartesian[:, 2] = r * np.cos(phi) + self.origin[2]
 
         return super().write_data_to_vtk(
             points=pts_cartesian,

src/mesh.cpp

@@ -467,7 +467,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
@@ -500,6 +499,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
+  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) {
@@ -1003,13 +1007,15 @@ void RectilinearMesh::to_hdf5(hid_t group) const
 // CylindricalMesh implementation
 //==============================================================================
 
-CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node}
+CylindricalMesh::CylindricalMesh(pugi::xml_node node)
+  : PeriodicStructuredMesh {node}
 {
   n_dimension_ = 3;
 
   grid_[0] = get_node_array<double>(node, "r_grid");
   grid_[1] = get_node_array<double>(node, "phi_grid");
   grid_[2] = get_node_array<double>(node, "z_grid");
+  origin_ = get_node_position(node, "origin");
 
   if (int err = set_grid()) {
     fatal_error(openmc_err_msg);
@@ -1026,10 +1032,12 @@ std::string CylindricalMesh::get_mesh_type() const
 StructuredMesh::MeshIndex CylindricalMesh::get_indices(
   Position r, bool& in_mesh) const
 {
-  Position mapped_r;
+  local_coords(r);
 
+  Position mapped_r;
   mapped_r[0] = std::hypot(r.x, r.y);
   mapped_r[2] = r[2];
+
   if (mapped_r[0] < FP_PRECISION) {
     mapped_r[1] = 0.0;
   } else {
@@ -1141,22 +1149,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]);
   }
 }
 
@@ -1227,6 +1236,7 @@ void CylindricalMesh::to_hdf5(hid_t group) const
   write_dataset(mesh_group, "r_grid", grid_[0]);
   write_dataset(mesh_group, "phi_grid", grid_[1]);
   write_dataset(mesh_group, "z_grid", grid_[2]);
+  write_dataset(mesh_group, "origin", origin_);
 
   close_group(mesh_group);
 }
@@ -1235,13 +1245,15 @@ void CylindricalMesh::to_hdf5(hid_t group) const
 // SphericalMesh implementation
 //==============================================================================
 
-SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node}
+SphericalMesh::SphericalMesh(pugi::xml_node node)
+  : PeriodicStructuredMesh {node}
 {
   n_dimension_ = 3;
 
   grid_[0] = get_node_array<double>(node, "r_grid");
   grid_[1] = get_node_array<double>(node, "theta_grid");
   grid_[2] = get_node_array<double>(node, "phi_grid");
+  origin_ = get_node_position(node, "origin");
 
   if (int err = set_grid()) {
     fatal_error(openmc_err_msg);
@@ -1258,9 +1270,11 @@ std::string SphericalMesh::get_mesh_type() const
 StructuredMesh::MeshIndex SphericalMesh::get_indices(
   Position r, bool& in_mesh) const
 {
-  Position mapped_r;
+  local_coords(r);
 
+  Position mapped_r;
   mapped_r[0] = r.norm();
+
   if (mapped_r[0] < FP_PRECISION) {
     mapped_r[1] = 0.0;
     mapped_r[2] = 0.0;
@@ -1490,6 +1504,7 @@ void SphericalMesh::to_hdf5(hid_t group) const
   write_dataset(mesh_group, "r_grid", grid_[0]);
   write_dataset(mesh_group, "theta_grid", grid_[1]);
   write_dataset(mesh_group, "phi_grid", grid_[2]);
+  write_dataset(mesh_group, "origin", origin_);
 
   close_group(mesh_group);
 }

src/xml_interface.cpp

@@ -48,4 +48,11 @@ bool get_node_value_bool(pugi::xml_node node, const char* name)
   return false;
 }
 
+Position get_node_position(
+  pugi::xml_node node, const char* name, bool lowercase)
+{
+  vector<double> arr = get_node_array<double>(node, name, lowercase);
+  return Position(arr);
+}
+
 } // namespace openmc

tests/regression_tests/filter_mesh/inputs_true.dat

@@ -57,11 +57,13 @@
     <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>
   </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>
     <theta_grid>0.0 0.39269908169872414 0.7853981633974483 1.1780972450961724 1.5707963267948966 1.9634954084936207 2.356194490192345 2.748893571891069 3.141592653589793</theta_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>
+    <origin>0.0 0.0 0.0</origin>
   </mesh>
   <filter id="1" type="mesh">
     <bins>1</bins>

tests/unit_tests/test_cylindrical_mesh.py

@@ -0,0 +1,105 @@
+from itertools import product, permutations
+
+import openmc
+import numpy as np
+
+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 = 2000
+    settings.batches = 10
+    settings.run_mode = 'fixed source'
+
+    # build
+    mesh = openmc.CylindricalMesh()
+    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()
+
+    mesh_filter = openmc.MeshFilter(mesh)
+    tally.filters.append(mesh_filter)
+
+    tally.scores.append("flux")
+
+    tallies = openmc.Tallies([tally])
+
+    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
+    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)))
+
+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 half of the mesh that is outside of the geometry
+        # contains the zero values
+        mean = tally.get_reshaped_data('mean', expand_dims=True)
+        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

tests/unit_tests/test_deplete_activation.py

@@ -33,7 +33,7 @@ def model():
     )
     model.settings.run_mode = 'fixed source'
 
-    rx_tally = openmc.Tally()
+    rx_tally = openmc.Tally(name='activation tally')
     rx_tally.scores = ['(n,gamma)']
     model.tallies.append(rx_tally)
 
@@ -53,7 +53,7 @@ def test_activation(run_in_tmpdir, model, reaction_rate_mode, reaction_rate_opts
     # Determine (n.gamma) reaction rate using initial run
     sp = model.run()
     with openmc.StatePoint(sp) as sp:
-        tally = sp.tallies[1]
+        tally = sp.get_tally(name='activation tally')
         capture_rate = tally.mean.flat[0]
 
     # Create one-nuclide depletion chain

tests/unit_tests/test_spherical_mesh.py

@@ -0,0 +1,109 @@
+from itertools import product, permutations
+
+import openmc
+import numpy as np
+
+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 = 2000
+    settings.batches = 10
+    settings.run_mode = 'fixed source'
+
+    # build
+    mesh = openmc.SphericalMesh()
+    mesh.phi_grid = np.linspace(0, 2*np.pi, 13)
+    mesh.theta_grid = np.linspace(0, np.pi, 7)
+    mesh.r_grid = np.linspace(0, geom_size, geom_size)
+
+    tally = openmc.Tally()
+
+    mesh_filter = openmc.MeshFilter(mesh)
+    tally.filters.append(mesh_filter)
+
+    tally.scores.append("flux")
+
+    tallies = openmc.Tallies([tally])
+
+    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
+    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)))
+origins |= set(permutations((offset, 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 spherical mesh
+    mesh.origin = origin
+
+    sp_filename = model.run()
+
+    with openmc.StatePoint(sp_filename) as sp:
+        tally = sp.tallies[1]
+
+        # we've translated half of the spherical 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 half of the mesh that is outside of the geometry
+        # contains the zero values
+        mean = tally.get_reshaped_data('mean', expand_dims=True)
+        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