crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Fix negative distances from bins_crossed for CylindricalMesh (#3370)

Browse source
This commit is contained in:
Patrick Shriwise 2025-04-22 23:00:44 -05:00 • committed by GitHub
parent a113440986
commit 5dd6ff6527
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
6 changed files with 29 additions and 33 deletions
Download patch file Download diff file Expand all files Collapse all files

5
include/openmc/mesh.h
View file

@ -138,9 +138,6 @@ public:
//! 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 @@ 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_;

37
src/mesh.cpp
View file

@ -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<MeshDistance, 3> 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();

4
tests/regression_tests/filter_mesh/inputs_true.dat
View file

@ -54,8 +54,8 @@
<mesh id="5" type="cylindrical">
<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>
<z_grid>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</z_grid>
<origin>0.0 0.0 -7.5</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>

3
tests/regression_tests/filter_mesh/test.py
View file

@ -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)

11
tests/unit_tests/test_cylindrical_mesh.py
View file

@ -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()

2
tests/unit_tests/test_spherical_mesh.py
View file

@ -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)))