From c4e39ae0cb5cce455f46d47a53fd19768d972c98 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 9 Apr 2025 17:01:03 -0500 Subject: [PATCH 1/7] Use local coordinates to determine initial mesh indices --- src/mesh.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 101da746849..c1979689eb2 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -925,6 +925,7 @@ void StructuredMesh::raytrace_mesh( // Calculate index of current cell. Offset the position a tiny bit in // direction of flight + local_coords(r0); MeshIndex ijk = get_indices(r0 + TINY_BIT * u, in_mesh); // if track is very short, assume that it is completely inside one cell. @@ -939,7 +940,7 @@ void StructuredMesh::raytrace_mesh( // 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(r0); local_coords(r1); // Calculate initial distances to next surfaces in all three dimensions From 17f6a82cbc9df2ea48683a88529113bedcdc13c7 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 9 Apr 2025 23:05:58 -0500 Subject: [PATCH 2/7] Adjust how global/local positions are handled in raytrace_mesh --- src/mesh.cpp | 29 +++++++++++++---------------- 1 file changed, 13 insertions(+), 16 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index c1979689eb2..e03119dabdd 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -910,6 +910,12 @@ void StructuredMesh::raytrace_mesh( // compilers will (hopefully) eliminate the complete code (including // calculation of parameters) but for the future: be explicit + // 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; + local_coords(r0); + local_coords(r1); + // Compute the length of the entire track. double total_distance = (r1 - r0).norm(); if (total_distance == 0.0 && settings::solver_type != SolverType::RANDOM_RAY) @@ -925,8 +931,7 @@ void StructuredMesh::raytrace_mesh( // Calculate index of current cell. Offset the position a tiny bit in // direction of flight - local_coords(r0); - 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 @@ -937,12 +942,6 @@ 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) { @@ -1006,7 +1005,7 @@ 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); @@ -1631,23 +1630,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]); } } From 0443b2d924abe0f5b74a2320f64977c03dc02f57 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 9 Apr 2025 23:13:28 -0500 Subject: [PATCH 3/7] Set cylinder origin and adjust Z grid to represnt equivalent cylindrical mesh --- tests/regression_tests/filter_mesh/inputs_true.dat | 4 ++-- tests/regression_tests/filter_mesh/test.py | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) 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) From 993e23d39364bb603cc9f87f723f77b2a8bf4036 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 10 Apr 2025 12:30:07 -0500 Subject: [PATCH 4/7] Rework local variable naming of coordinates in raytrace_mesh. Update test values for cylindrical mesh. --- include/openmc/mesh.h | 4 ++-- src/mesh.cpp | 19 ++++++++++--------- tests/unit_tests/test_cylindrical_mesh.py | 11 ++++++++--- tests/unit_tests/test_spherical_mesh.py | 2 -- 4 files changed, 20 insertions(+), 16 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 8faf45f1b99..d94af3ee6b1 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -139,7 +139,7 @@ class Mesh { virtual void prepare_for_point_location() {}; //! Update a position to the local coordinates of the mesh - virtual void local_coords(Position& r) const {}; + virtual Position& local_coords(Position& r) const { return r; }; //! Return a position in the local coordinates of the mesh virtual Position local_coords(const Position& r) const { return r; }; @@ -427,7 +427,7 @@ 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(Position& r) const override { return r -= origin_; }; Position local_coords(const Position& r) const override { diff --git a/src/mesh.cpp b/src/mesh.cpp index e03119dabdd..96a4967a629 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -910,16 +910,17 @@ void StructuredMesh::raytrace_mesh( // compilers will (hopefully) eliminate the complete code (including // calculation of parameters) but for the future: be explicit - // 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; - local_coords(r0); - local_coords(r1); // Compute the length of the entire track. double total_distance = (r1 - r0).norm(); if (total_distance == 0.0 && settings::solver_type != SolverType::RANDOM_RAY) - return; + 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; + + const Position& local_r = local_coords(r0); const int n = n_dimension_; @@ -945,7 +946,7 @@ void StructuredMesh::raytrace_mesh( // 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 +976,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])); @@ -1008,7 +1009,7 @@ void StructuredMesh::raytrace_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 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))) From a01b659979f1384a2d5e2d8b98d277e23a27385c Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 10 Apr 2025 12:34:12 -0500 Subject: [PATCH 5/7] Clang format --- src/mesh.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 96a4967a629..d1f8f818453 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -910,11 +910,10 @@ void StructuredMesh::raytrace_mesh( // compilers will (hopefully) eliminate the complete code (including // calculation of parameters) but for the future: be explicit - // Compute the length of the entire track. double total_distance = (r1 - r0).norm(); if (total_distance == 0.0 && settings::solver_type != SolverType::RANDOM_RAY) - return; + return; // keep a copy of the original global position to pass to get_indices, // which performs its own transformation to local coordinates From 0b761b7e5966f2e2d1f899abf613f8d1d8f3dc39 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 11 Apr 2025 09:46:51 -0500 Subject: [PATCH 6/7] Rm whitespace --- src/mesh.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index d1f8f818453..ed9aa0e2503 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -918,7 +918,6 @@ void StructuredMesh::raytrace_mesh( // 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; - const Position& local_r = local_coords(r0); const int n = n_dimension_; From 3957796423c5951c34e9ac6456ed028b823d5eda Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 14 Apr 2025 15:39:42 -0500 Subject: [PATCH 7/7] Single local_coords implementation --- include/openmc/mesh.h | 5 ----- src/mesh.cpp | 6 +++--- 2 files changed, 3 insertions(+), 8 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index d94af3ee6b1..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 Position& local_coords(Position& r) const { return r; }; - //! 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} {}; - Position& local_coords(Position& r) const override { return r -= origin_; }; - Position local_coords(const Position& r) const override { return r - origin_; diff --git a/src/mesh.cpp b/src/mesh.cpp index ed9aa0e2503..8280177ac13 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -918,7 +918,7 @@ void StructuredMesh::raytrace_mesh( // 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; - const Position& local_r = local_coords(r0); + Position local_r = local_coords(r0); const int n = n_dimension_; @@ -1481,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); @@ -1759,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();