diff --git a/include/openmc/constants.h b/include/openmc/constants.h index a1d94e5819e..1c13bc071d7 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -58,6 +58,11 @@ constexpr double FP_COINCIDENT {1e-12}; constexpr double TORUS_TOL {1e-10}; constexpr double RADIAL_MESH_TOL {1e-10}; +// Tolerance on the normalized normal of a general plane for treating that +// plane as axis-aligned when computing a bounding box. Matches the value of +// Surface._atol used by PlaneMixin.bounding_box in openmc/surface.py. +constexpr double PLANE_ALIGNMENT_TOL {1e-12}; + // Maximum number of random samples per history constexpr int MAX_SAMPLE {100000}; diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 2d8580345a4..b459d847dc1 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -171,6 +171,7 @@ class SurfacePlane : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + BoundingBox bounding_box(bool pos_side) const override; double A_, B_, C_, D_; }; @@ -337,6 +338,7 @@ class SurfaceXTorus : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + BoundingBox bounding_box(bool pos_side) const override; double x0_, y0_, z0_, A_, B_, C_; }; @@ -354,6 +356,7 @@ class SurfaceYTorus : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + BoundingBox bounding_box(bool pos_side) const override; double x0_, y0_, z0_, A_, B_, C_; }; @@ -371,6 +374,7 @@ class SurfaceZTorus : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + BoundingBox bounding_box(bool pos_side) const override; double x0_, y0_, z0_, A_, B_, C_; }; diff --git a/openmc/surface.py b/openmc/surface.py index c2afeb613ac..441d9486ab2 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -561,22 +561,18 @@ def bounding_box(self, side): nhat = self._get_normal() ll = np.array([-np.inf, -np.inf, -np.inf]) ur = np.array([np.inf, np.inf, np.inf]) - # If the plane is axis aligned, find the proper bounding box - if np.any(np.isclose(np.abs(nhat), 1., rtol=0., atol=self._atol)): - sign = nhat.sum() - a, b, c, d = self._get_base_coeffs() - vals = [d/val if not np.isclose(val, 0., rtol=0., atol=self._atol) - else np.nan for val in (a, b, c)] - if side == '-': - if sign > 0: - ur = np.array([v if not np.isnan(v) else np.inf for v in vals]) - else: - ll = np.array([v if not np.isnan(v) else -np.inf for v in vals]) - elif side == '+': - if sign > 0: - ll = np.array([v if not np.isnan(v) else -np.inf for v in vals]) - else: - ur = np.array([v if not np.isnan(v) else np.inf for v in vals]) + # A plane only bounds a half-space when its normal is parallel to a + # coordinate axis, in which case it bounds it along that axis alone. + axis = int(np.argmax(np.abs(nhat))) + on_axis = np.isclose(abs(nhat[axis]), 1., rtol=0., atol=self._atol) + off_axis = np.delete(nhat, axis) + if on_axis and np.all(np.isclose(off_axis, 0., rtol=0., atol=self._atol)): + coeffs = self._get_base_coeffs() + intercept = coeffs[3]/coeffs[axis] + if (side == '+') == (coeffs[axis] > 0): + ll[axis] = intercept + else: + ur[axis] = intercept return BoundingBox(ll, ur) diff --git a/src/surface.cpp b/src/surface.cpp index 81b756deae7..06f0f88b3c4 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -352,6 +352,50 @@ double SurfacePlane::evaluate(Position r) const return A_ * r.x + B_ * r.y + C_ * r.z - D_; } +BoundingBox SurfacePlane::bounding_box(bool pos_side) const +{ + // A general plane bounds a half-space in one direction only when its normal + // is parallel to a coordinate axis; otherwise both half-spaces are unbounded + // along every axis. This mirrors PlaneMixin.bounding_box on the Python side, + // so that a plane whose off-axis coefficients are rotation-matrix roundoff + // (e.g. B = 1 with A = C = 6.1e-17) yields the same box through both APIs. + const Direction n = normal({}); + const double norm = n.norm(); + if (norm == 0.0) + return {}; + + int axis = -1; + for (int i = 0; i < 3; ++i) { + if (std::abs(std::abs(n[i] / norm) - 1.0) > PLANE_ALIGNMENT_TOL) + continue; + + bool aligned = true; + for (int j = 0; j < 3; ++j) { + if (j != i && std::abs(n[j] / norm) > PLANE_ALIGNMENT_TOL) { + aligned = false; + break; + } + } + if (aligned) { + axis = i; + break; + } + } + if (axis == -1) + return {}; + + // The half-space is bounded below when the outward normal points along the + // positive axis direction and we are on the positive side, or vice versa. + BoundingBox bbox; + const double intercept = D_ / n[axis]; + if (pos_side == (n[axis] > 0.0)) { + bbox.min[axis] = intercept; + } else { + bbox.max[axis] = intercept; + } + return bbox; +} + double SurfacePlane::distance(Position r, Direction u, bool coincident) const { const double f = A_ * r.x + B_ * r.y + C_ * r.z - D_; @@ -1034,6 +1078,17 @@ double SurfaceXTorus::evaluate(Position r) const std::pow(std::sqrt(y * y + z * z) - A_, 2) / (C_ * C_) - 1.; } +BoundingBox SurfaceXTorus::bounding_box(bool pos_side) const +{ + // The torus interior is compact: it extends +/-B_ along the axis of + // revolution and +/-(A_ + C_) in the two perpendicular directions. Mirrors + // XTorus.bounding_box on the Python side. + if (pos_side) + return {}; + return {{x0_ - B_, y0_ - A_ - C_, z0_ - A_ - C_}, + {x0_ + B_, y0_ + A_ + C_, z0_ + A_ + C_}}; +} + double SurfaceXTorus::distance(Position r, Direction u, bool coincident) const { double x = r.x - x0_; @@ -1087,6 +1142,14 @@ double SurfaceYTorus::evaluate(Position r) const std::pow(std::sqrt(x * x + z * z) - A_, 2) / (C_ * C_) - 1.; } +BoundingBox SurfaceYTorus::bounding_box(bool pos_side) const +{ + if (pos_side) + return {}; + return {{x0_ - A_ - C_, y0_ - B_, z0_ - A_ - C_}, + {x0_ + A_ + C_, y0_ + B_, z0_ + A_ + C_}}; +} + double SurfaceYTorus::distance(Position r, Direction u, bool coincident) const { double x = r.x - x0_; @@ -1140,6 +1203,14 @@ double SurfaceZTorus::evaluate(Position r) const std::pow(std::sqrt(x * x + y * y) - A_, 2) / (C_ * C_) - 1.; } +BoundingBox SurfaceZTorus::bounding_box(bool pos_side) const +{ + if (pos_side) + return {}; + return {{x0_ - A_ - C_, y0_ - A_ - C_, z0_ - B_}, + {x0_ + A_ + C_, y0_ + A_ + C_, z0_ + B_}}; +} + double SurfaceZTorus::distance(Position r, Direction u, bool coincident) const { double x = r.x - x0_; diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt index 991f219f528..b6892549bde 100644 --- a/tests/cpp_unit_tests/CMakeLists.txt +++ b/tests/cpp_unit_tests/CMakeLists.txt @@ -9,6 +9,7 @@ set(TEST_NAMES test_photon test_ray test_region + test_surface test_tensor test_geometry # Add additional unit test files here diff --git a/tests/cpp_unit_tests/test_surface.cpp b/tests/cpp_unit_tests/test_surface.cpp new file mode 100644 index 00000000000..1529cf21031 --- /dev/null +++ b/tests/cpp_unit_tests/test_surface.cpp @@ -0,0 +1,186 @@ +#include + +#include +#include + +#include + +#include "openmc/cell.h" +#include "openmc/constants.h" +#include "openmc/surface.h" + +using namespace openmc; + +namespace { + +template +std::unique_ptr make_surface( + pugi::xml_document& doc, int id, const char* type, const char* coeffs) +{ + pugi::xml_node n = doc.append_child("surface"); + n.append_attribute("id") = id; + n.append_attribute("type") = type; + n.append_attribute("coeffs") = coeffs; + return std::make_unique(n); +} + +} // anonymous namespace + +TEST_CASE("General plane bounding box") +{ + pugi::xml_document doc; + + SECTION("Exactly axis-aligned planes bound one axis only") + { + // +x normal: the positive half-space starts at x = D/A + auto px = make_surface(doc, 1, "plane", "1.0 0.0 0.0 5.0"); + BoundingBox pos = px->bounding_box(true); + CHECK(pos.min.x == Catch::Approx(5.0)); + CHECK(pos.min.y == -INFTY); + CHECK(pos.min.z == -INFTY); + CHECK(pos.max.x == INFTY); + + BoundingBox neg = px->bounding_box(false); + CHECK(neg.max.x == Catch::Approx(5.0)); + CHECK(neg.min.x == -INFTY); + CHECK(neg.max.y == INFTY); + + // -y normal: the sense of the bound flips with the sign of the normal + auto ny = make_surface(doc, 2, "plane", "0.0 -1.0 0.0 3.0"); + BoundingBox ny_pos = ny->bounding_box(true); + CHECK(ny_pos.max.y == Catch::Approx(-3.0)); + CHECK(ny_pos.min.y == -INFTY); + CHECK(ny_pos.min.x == -INFTY); + + BoundingBox ny_neg = ny->bounding_box(false); + CHECK(ny_neg.min.y == Catch::Approx(-3.0)); + CHECK(ny_neg.max.y == INFTY); + } + + SECTION("Coefficients need not be normalized") + { + // 4z - 10 = 0 is the same plane as z = 2.5 + auto pz = make_surface(doc, 3, "plane", "0.0 0.0 4.0 10.0"); + CHECK(pz->bounding_box(true).min.z == Catch::Approx(2.5)); + CHECK(pz->bounding_box(false).max.z == Catch::Approx(2.5)); + } + + SECTION("Rotation roundoff is still axis aligned (issue #2632)") + { + // The surface s64 from the reported model: a y-plane at -1.3 written with + // cos(pi/2) in the x and z slots. + auto p = make_surface( + doc, 4, "plane", "6.123233995736766e-17 1.0 6.123233995736766e-17 -1.3"); + BoundingBox pos = p->bounding_box(true); + CHECK(pos.min.y == Catch::Approx(-1.3)); + // The two off-axis directions must stay unbounded + CHECK(pos.min.x == -INFTY); + CHECK(pos.min.z == -INFTY); + CHECK(pos.max.x == INFTY); + CHECK(pos.max.z == INFTY); + } + + SECTION("Oblique planes bound nothing") + { + auto p = make_surface( + doc, 5, "plane", "0.7071067811865476 0.0 0.7071067811865476 11.45"); + for (bool side : {false, true}) { + BoundingBox bb = p->bounding_box(side); + CHECK(bb.min.x == -INFTY); + CHECK(bb.min.y == -INFTY); + CHECK(bb.min.z == -INFTY); + CHECK(bb.max.x == INFTY); + CHECK(bb.max.y == INFTY); + CHECK(bb.max.z == INFTY); + } + } + + SECTION("A plane tilted well beyond tolerance bounds nothing") + { + // 1e-5 is far outside PLANE_ALIGNMENT_TOL, so this must not be treated as + // an x-plane, and in particular must not acquire a bound on y. + auto p = make_surface(doc, 6, "plane", "1.0 1e-5 0.0 5.0"); + BoundingBox bb = p->bounding_box(true); + CHECK(bb.min.x == -INFTY); + CHECK(bb.min.y == -INFTY); + } + + SECTION("An off-axis coefficient above tolerance bounds nothing") + { + // Although the x component of the normalized normal is within + // PLANE_ALIGNMENT_TOL of one, the y component is above the tolerance. + auto p = make_surface(doc, 7, "plane", "1.0 1e-11 0.0 5.0"); + for (bool side : {false, true}) { + BoundingBox bb = p->bounding_box(side); + CHECK(bb.min.x == -INFTY); + CHECK(bb.min.y == -INFTY); + CHECK(bb.min.z == -INFTY); + CHECK(bb.max.x == INFTY); + CHECK(bb.max.y == INFTY); + CHECK(bb.max.z == INFTY); + } + } + + SECTION("A degenerate plane bounds nothing") + { + auto p = make_surface(doc, 8, "plane", "0.0 0.0 0.0 1.0"); + BoundingBox bb = p->bounding_box(true); + CHECK(bb.min.x == -INFTY); + CHECK(bb.max.x == INFTY); + } +} + +TEST_CASE("Torus bounding box") +{ + pugi::xml_document doc; + + // x0 y0 z0 A B C, so the interior spans +/-B along the axis of revolution + // and +/-(A + C) in the perpendicular directions. + SECTION("x-torus") + { + auto t = make_surface( + doc, 1, "x-torus", "1.0 2.0 3.0 5.0 0.5 0.25"); + BoundingBox in = t->bounding_box(false); + CHECK(in.min.x == Catch::Approx(0.5)); + CHECK(in.max.x == Catch::Approx(1.5)); + CHECK(in.min.y == Catch::Approx(-3.25)); + CHECK(in.max.y == Catch::Approx(7.25)); + CHECK(in.min.z == Catch::Approx(-2.25)); + CHECK(in.max.z == Catch::Approx(8.25)); + + // The exterior of a torus is unbounded + BoundingBox out = t->bounding_box(true); + CHECK(out.min.x == -INFTY); + CHECK(out.max.z == INFTY); + } + + SECTION("y-torus") + { + auto t = make_surface( + doc, 2, "y-torus", "1.0 2.0 3.0 5.0 0.5 0.25"); + BoundingBox in = t->bounding_box(false); + CHECK(in.min.y == Catch::Approx(1.5)); + CHECK(in.max.y == Catch::Approx(2.5)); + CHECK(in.min.x == Catch::Approx(-4.25)); + CHECK(in.max.x == Catch::Approx(6.25)); + CHECK(in.min.z == Catch::Approx(-2.25)); + CHECK(in.max.z == Catch::Approx(8.25)); + + CHECK(t->bounding_box(true).max.y == INFTY); + } + + SECTION("z-torus") + { + auto t = make_surface( + doc, 3, "z-torus", "1.0 2.0 3.0 5.0 0.5 0.25"); + BoundingBox in = t->bounding_box(false); + CHECK(in.min.z == Catch::Approx(2.5)); + CHECK(in.max.z == Catch::Approx(3.5)); + CHECK(in.min.x == Catch::Approx(-4.25)); + CHECK(in.max.x == Catch::Approx(6.25)); + CHECK(in.min.y == Catch::Approx(-3.25)); + CHECK(in.max.y == Catch::Approx(7.25)); + + CHECK(t->bounding_box(true).min.z == -INFTY); + } +} diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index f5c0b8f8b62..887ac319bfc 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -4,6 +4,7 @@ import numpy as np import math import openmc +import openmc.lib import pytest @@ -72,6 +73,60 @@ def test_plane_from_points(): assert s.d == 1.0 +def test_plane_bounding_box(): + # A plane written by a rotation has roundoff in the off-axis slots but is + # still a y-plane and should bound the half-space along y alone (#2632) + eps = math.cos(math.pi/2) + s = openmc.Plane(eps, 1., eps, -1.3) + ll, ur = (+s).bounding_box + assert ll == pytest.approx((-np.inf, -1.3, -np.inf)) + assert np.all(np.isinf(ur)) + ll, ur = (-s).bounding_box + assert ur == pytest.approx((np.inf, -1.3, np.inf)) + assert np.all(np.isinf(ll)) + + # The intercept comes from the raw coefficients, which need not be + # normalized: 4z - 10 = 0 is the plane z = 2.5 + s = openmc.Plane(0., 0., 4., 10.) + assert (+s).bounding_box[0] == pytest.approx((-np.inf, -np.inf, 2.5)) + assert (-s).bounding_box[1] == pytest.approx((np.inf, np.inf, 2.5)) + + # The sense of the bound flips with the sign of the normal + s = openmc.Plane(0., -1., 0., 3.) + assert (+s).bounding_box[1] == pytest.approx((np.inf, -3., np.inf)) + assert (-s).bounding_box[0] == pytest.approx((-np.inf, -3., -np.inf)) + + # An off-axis normalized coefficient above the alignment tolerance makes + # the half-space unbounded along every axis + s = openmc.Plane(1., 1.e-11, 0., 5.) + assert_infinite_bb(s) + p = (4., 2.e11, 0.) + assert p in +s + + # A plane tilted well beyond the alignment tolerance bounds nothing + assert_infinite_bb(openmc.Plane(1., 1.e-5, 0., 5.)) + + +def test_plane_bounding_box_lib(mpi_intracomm): + """Check Python and C++ bounding consistency.""" + openmc.reset_auto_ids() + + cyl = openmc.YCylinder(r=17.7) + y_min = openmc.Plane(6.123233995736766e-17, 1., 6.123233995736766e-17, -1.3) + y_max = openmc.YPlane(5.6, boundary_type='vacuum') + cell = openmc.Cell(region=+y_min & -y_max & -cyl) + model = openmc.Model(geometry=openmc.Geometry([cell])) + model.settings.batches = 1 + model.settings.particles = 100 + + python_box = cell.bounding_box + with openmc.lib.TemporarySession(model, intracomm=mpi_intracomm): + lib_box = openmc.lib.cells[cell.id].bounding_box + + np.testing.assert_array_equal(python_box.lower_left, lib_box.lower_left) + np.testing.assert_array_equal(python_box.upper_right, lib_box.upper_right) + + def test_xplane(): s = openmc.XPlane(3., boundary_type='reflective') assert s.x0 == 3. @@ -195,12 +250,12 @@ def test_cylinder(): assert s.dy == -1 assert s.dz == 1 assert s.r == 2 - + # Check radius must be positive with pytest.raises(ValueError): openmc.Cylinder(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r=0.0) with pytest.raises(ValueError): - openmc.Cylinder(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r=-1.0) + openmc.Cylinder(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r=-1.0) # Check bounding box assert_infinite_bb(s) @@ -250,7 +305,7 @@ def test_xcylinder(): assert s.y0 == y assert s.z0 == z assert s.r == r - + # Check radius must be positive with pytest.raises(ValueError): openmc.XCylinder(y0=y, z0=z, r=0.0) @@ -302,7 +357,7 @@ def test_ycylinder(): assert s.x0 == x assert s.z0 == z assert s.r == r - + # Check radius must be positive with pytest.raises(ValueError): openmc.YCylinder(x0=x, z0=z, r=0.0) @@ -345,7 +400,7 @@ def test_zcylinder(): assert s.x0 == x assert s.y0 == y assert s.r == r - + # Check radius must be positive with pytest.raises(ValueError): openmc.ZCylinder(x0=x, y0=y, r=0.0) @@ -389,7 +444,7 @@ def test_sphere(): assert s.y0 == y assert s.z0 == z assert s.r == r - + # Check radius must be positive with pytest.raises(ValueError): openmc.Sphere(x0=x, y0=y, z0=z, r=0.0) @@ -434,7 +489,7 @@ def cone_common(apex, r2, cls): assert s.y0 == y assert s.z0 == z assert s.r2 == r2 - + # Check radius must be positive with pytest.raises(ValueError): cls(x0=x, y0=y, z0=z, r2=0.0) @@ -478,12 +533,12 @@ def test_cone(): assert s.dy == -1 assert s.dz == 1 assert s.r2 == 4 - + # Check radius must be positive with pytest.raises(ValueError): openmc.Cone(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r2=0.0) with pytest.raises(ValueError): - openmc.Cone(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r2=-1.0) + openmc.Cone(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r2=-1.0) # Check bounding box assert_infinite_bb(s) @@ -664,9 +719,9 @@ def torus_common(center, R, r1, r2, cls): assert s.a == R assert s.b == r1 assert s.c == r2 - + # Check radius must be positive - params = [(0.0, r1, r2), (R, 0.0, r2), (R, r1, 0.0), + params = [(0.0, r1, r2), (R, 0.0, r2), (R, r1, 0.0), (-1.0, r1, r2), (R, -1.0, r2), (R, r1, -1.0)] for a,b,c in params: with pytest.raises(ValueError):