diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index c690dc4fee..9554a51436 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -241,6 +241,7 @@ public: xt::xtensor lower_left_; //!< Lower-left coordinates of mesh xt::xtensor upper_right_; //!< Upper-right coordinates of mesh std::array 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, 3> grid_; - Position origin_; int set_grid(); @@ -414,7 +414,6 @@ public: void to_hdf5(hid_t group) const override; array, 3> grid_; - Position origin_; int set_grid(); diff --git a/openmc/mesh.py b/openmc/mesh.py index 307d9d952c..47d34ebe2e 100644 --- a/openmc/mesh.py +++ b/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: diff --git a/src/mesh.cpp b/src/mesh.cpp index 82644db08b..e0fa305edd 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -412,7 +412,6 @@ template 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 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; diff --git a/tests/unit_tests/test_cylindrical_mesh.py b/tests/unit_tests/test_cylindrical_mesh.py index 952e8cf5a2..e67f21b15c 100644 --- a/tests/unit_tests/test_cylindrical_mesh.py +++ b/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) \ No newline at end of file + 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