diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 8faf45f1b99..5a727c2b6f1 100644 --- a/include/openmc/mesh.h +++ b/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_; diff --git a/src/mesh.cpp b/src/mesh.cpp index 101da746849..8280177ac13 100644 --- a/src/mesh.cpp +++ b/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 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(); diff --git a/tests/regression_tests/filter_mesh/inputs_true.dat b/tests/regression_tests/filter_mesh/inputs_true.dat index 0667c034148..4995749bbd3 100644 --- a/tests/regression_tests/filter_mesh/inputs_true.dat +++ b/tests/regression_tests/filter_mesh/inputs_true.dat @@ -54,8 +54,8 @@ 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 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 - -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 - 0.0 0.0 0.0 + 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 + 0.0 0.0 -7.5 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 diff --git a/tests/regression_tests/filter_mesh/test.py b/tests/regression_tests/filter_mesh/test.py index 6214e0486c0..165ba2a0c09 100644 --- a/tests/regression_tests/filter_mesh/test.py +++ b/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) diff --git a/tests/unit_tests/test_cylindrical_mesh.py b/tests/unit_tests/test_cylindrical_mesh.py index b408bc0e917..645269825a8 100644 --- a/tests/unit_tests/test_cylindrical_mesh.py +++ b/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() diff --git a/tests/unit_tests/test_spherical_mesh.py b/tests/unit_tests/test_spherical_mesh.py index 0b579be35f7..a95a4151a2e 100644 --- a/tests/unit_tests/test_spherical_mesh.py +++ b/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)))