From 3aebbaceb6c6bf1a1c65bbedbeae65f831be0e4b Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 9 Sep 2026 19:00:23 +0000 Subject: [PATCH 1/5] Fix coordinate frame of surface normals used by surface tallies Surface::normal() returns the outward normal in the local coordinate frame of the universe that owns the surface, and RectLattice/HexLattice get_normal() likewise return the tile boundary normal in the lattice's own frame. Both were being dotted with Particle::u(), which is the direction at coordinate level zero (the root frame). Whenever the surface or lattice lived under a cell carrying a rotation, the resulting cosine was wrong. The normal was also being evaluated at the root-frame position r() rather than the local position, so for any position-dependent normal (sphere, cylinder, cone, torus, quadric) it was wrong even without a rotation, and it was computed after the crossing had already been carried out, which for a periodic boundary means after the particle had been translated to the partner surface. Concretely, for a plane inside a universe filled into a cell rotated 45 degrees about z, with particles crossing along the lab +x axis: - the surface-crossing flux estimator scored w/|mu| = 1.0 instead of 1/cos(45 deg) = 1.4142, a 41% error; - MuSurfaceFilter binned the crossing at mu = 1.0 instead of 0.7071. Evaluate the normal at the local position, before the crossing, and rotate it up into the root frame with a new rotate_to_root() helper that walks back up the coordinate levels undoing each cell rotation. Net current is unaffected since only the sign of the cosine matters there, but the flux score and MuSurfaceFilter both are. MuSurfaceFilter recomputed the normal itself and could not have applied the same correction, because by the time filters run the coordinate levels no longer identify the surface's universe. It now reads the root-frame normal recorded on the particle by score_surface_tally(). As a side effect this removes an out-of-bounds access to model::surfaces: on a lattice crossing the surface token is SURFACE_NONE, so surface_index() returned -1. Add unit tests covering a rotated universe fill and a rotated lattice, for the flux score and for MuSurfaceFilter. All three fail on the current develop code and pass with this change. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_018gUxmUj7G6Aie4azS9uypa --- include/openmc/geometry.h | 18 ++++ include/openmc/particle_data.h | 11 +++ src/geometry.cpp | 17 ++++ src/particle.cpp | 42 ++++++--- src/tallies/filter_musurface.cpp | 10 +-- src/tallies/tally_scoring.cpp | 5 ++ tests/unit_tests/test_filter_musurface.py | 48 +++++++++++ tests/unit_tests/test_surface_flux.py | 100 ++++++++++++++++++++++ 8 files changed, 235 insertions(+), 16 deletions(-) diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index e8504d48261..5f96751e86b 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -8,6 +8,7 @@ #include "openmc/array.h" #include "openmc/constants.h" +#include "openmc/position.h" #include "openmc/random_ray/source_region.h" // For hash_combine #include "openmc/vector.h" @@ -86,6 +87,23 @@ int check_cell_overlap(GeometryState& p, bool error = true); int cell_instance_at_level(const GeometryState& p, int level); +//============================================================================== +//! Rotate a direction from a local coordinate frame into the root frame +//! +//! Surface and lattice normals are expressed in the local coordinate frame of +//! the universe that contains them. Comparing such a normal against the +//! particle's direction of travel (which is stored in the root frame at +//! coordinate level zero) requires undoing the rotations that were applied +//! while descending to \c level. +//! +//! \param p A particle whose coordinate levels give the chain of rotations +//! \param level The level (zero indexed) that \c u is expressed in +//! \param u A direction in the local frame of \c level +//! \return The same direction expressed in the root coordinate frame +//============================================================================== + +Direction rotate_to_root(const GeometryState& p, int level, Direction u); + //============================================================================== //! Locate a particle in the geometry tree and set its geometry data fields. //! diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 24b7eb53c0a..65266fdb8a2 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -384,6 +384,14 @@ class GeometryState { int& surface() { return surface_; } const int& surface() const { return surface_; } + // Outward unit normal of the surface (or lattice boundary) currently being + // crossed, expressed in the root coordinate frame. Set by + // score_surface_tally() immediately before the tally filters are evaluated, + // so that filters can compare it against the root-frame direction u() + // without having to know which coordinate level the surface lives in. + Direction& surface_normal() { return surface_normal_; } + const Direction& surface_normal() const { return surface_normal_; } + // Surface index based on the current value of the surface_ attribute int surface_index() const { @@ -436,6 +444,9 @@ class GeometryState { int surface_ { SURFACE_NONE}; //!< surface token for surface the particle is currently on + //! Outward normal of the surface being crossed, in the root coordinate frame + Direction surface_normal_ {0.0, 0.0, 1.0}; + BoundaryInfo boundary_; //!< Info about the next intersection int material_ {-1}; //!< index for current material diff --git a/src/geometry.cpp b/src/geometry.cpp index ecacf0bffbf..34a4bc10b53 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -124,6 +124,23 @@ int cell_instance_at_level(const GeometryState& p, int level) //============================================================================== +Direction rotate_to_root(const GeometryState& p, int level, Direction u) +{ + // Each coordinate level below the root was reached by applying the rotation + // matrix of the cell one level above it, so walk back up applying the + // inverse of each rotation in turn. Translations are irrelevant here since + // they do not affect directions. + for (int i = level; i > 0; --i) { + if (p.coord(i).rotated()) { + const auto& c {*model::cells[p.coord(i - 1).cell()]}; + u = u.inverse_rotate(c.rotation_); + } + } + return u; +} + +//============================================================================== + bool find_cell_inner( GeometryState& p, const NeighborList* neighbor_list, bool verbose) { diff --git a/src/particle.cpp b/src/particle.cpp index 98270af56fb..14316b08fac 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -344,32 +344,54 @@ void Particle::event_cross_surface() surface() = boundary().surface(); n_coord() = boundary().coord_level(); + // The surface or lattice being crossed belongs to the universe at the lowest + // coordinate level, so its normal is reported in that level's local frame + // while the particle direction used to score surface tallies lives in the + // root frame. The normal therefore has to be evaluated at the local position + // and rotated up into the root frame, and that has to happen before the + // crossing is carried out, since crossing invalidates the coordinate levels. + int i_surf_level = n_coord() - 1; + if (boundary().lattice_translation()[0] != 0 || boundary().lattice_translation()[1] != 0 || boundary().lattice_translation()[2] != 0) { // Particle crosses lattice boundary - int i_lattice = coord(boundary().coord_level() - 1).lattice(); - bool verbose = settings::verbosity >= 10 || trace(); - cross_lattice(*this, boundary(), verbose); - event() = TallyEvent::LATTICE; + int i_lattice = coord(i_surf_level).lattice(); - // Score cell to cell partial currents + // Determine the lattice boundary normal in the root frame before crossing + bool normal_is_valid = false; + Direction normal; if (!model::active_surface_tallies.empty()) { auto& lat {*model::lattices[i_lattice]}; bool is_valid; - Direction normal = - lat.get_normal(boundary().lattice_translation(), is_valid); + normal = lat.get_normal(boundary().lattice_translation(), is_valid); if (is_valid) { - normal /= normal.norm(); - score_surface_tally(*this, model::active_surface_tallies, normal); + normal = rotate_to_root(*this, i_surf_level, normal / normal.norm()); + normal_is_valid = true; } } + bool verbose = settings::verbosity >= 10 || trace(); + cross_lattice(*this, boundary(), verbose); + event() = TallyEvent::LATTICE; + + // Score cell to cell partial currents + if (normal_is_valid) { + score_surface_tally(*this, model::active_surface_tallies, normal); + } + } else { const auto& surf {*model::surfaces[surface_index()].get()}; + // Determine the surface normal in the root frame before crossing + Direction normal; + if (!model::active_surface_tallies.empty()) { + normal = surf.normal(r_local()); + normal = rotate_to_root(*this, i_surf_level, normal / normal.norm()); + } + // Particle crosses surface // If BC, add particle to surface source before crossing surface if (surf.surf_source_ && surf.bc_) { @@ -387,8 +409,6 @@ void Particle::event_cross_surface() // Score cell to cell partial currents if (!model::active_surface_tallies.empty()) { - Direction normal = surf.normal(r()); - normal /= normal.norm(); score_surface_tally(*this, model::active_surface_tallies, normal); } } diff --git a/src/tallies/filter_musurface.cpp b/src/tallies/filter_musurface.cpp index 340149d4cff..2ab90b0aac3 100644 --- a/src/tallies/filter_musurface.cpp +++ b/src/tallies/filter_musurface.cpp @@ -3,7 +3,6 @@ #include // for abs, copysign #include "openmc/search.h" -#include "openmc/surface.h" #include "openmc/tallies/tally_scoring.h" namespace openmc { @@ -11,10 +10,11 @@ namespace openmc { void MuSurfaceFilter::get_all_bins( const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - // Get surface normal (and make sure it is a unit vector) - const auto surf {model::surfaces[p.surface_index()].get()}; - auto n = surf->normal(p.r()); - n /= n.norm(); + // Use the normal recorded for the crossing being scored. It is already a + // unit vector expressed in the root coordinate frame, which is the frame + // p.u() is in -- recomputing it here from the surface would give the normal + // in the local frame of whichever universe holds the surface. + Direction n = p.surface_normal(); // Determine whether normal should be pointing in or out if (p.surface() < 0) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 27e6671833a..f213ec51a65 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2662,6 +2662,11 @@ void score_surface_tally( { double wgt = p.wgt_last(); + // Make the normal available to filters (e.g. MuSurfaceFilter) that need it. + // The caller is responsible for supplying it in the root coordinate frame so + // that it can be compared directly against p.u(). + p.surface_normal() = normal; + double mu = std::clamp(p.u().dot(normal), -1.0, 1.0); // Sign for net current: +1 if crossing outward (in direction of normal), diff --git a/tests/unit_tests/test_filter_musurface.py b/tests/unit_tests/test_filter_musurface.py index ca0db71f0c6..745b4002615 100644 --- a/tests/unit_tests/test_filter_musurface.py +++ b/tests/unit_tests/test_filter_musurface.py @@ -1,3 +1,5 @@ +import math + import openmc @@ -36,3 +38,49 @@ def test_musurface(run_in_tmpdir): assert element == 0.0 +def test_musurface_rotated_universe(run_in_tmpdir): + """MuSurfaceFilter uses the surface normal in the root coordinate frame. + + The plane lives in a universe filled into a cell rotated 45 degrees about + z, so its normal in the plane's own frame is (1, 0, 0) while the particle + direction is stored in the root frame. Binning must use the root-frame + normal, giving mu = cos(45 deg) rather than 1. + """ + openmc.reset_auto_ids() + + xplane = openmc.XPlane(0.0) + inner1 = openmc.Cell(region=-xplane) + inner2 = openmc.Cell(region=+xplane) + inner_univ = openmc.Universe(cells=[inner1, inner2]) + + sph = openmc.Sphere(r=10.0, boundary_type='vacuum') + root_cell = openmc.Cell(region=-sph, fill=inner_univ) + root_cell.rotation = (0.0, 0.0, 45.0) + + model = openmc.Model() + model.geometry = openmc.Geometry([root_cell]) + + src = openmc.IndependentSource() + src.space = openmc.stats.Point((-5.0, 0.0, 0.0)) + src.angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0)) + + model.settings.run_mode = 'fixed source' + model.settings.batches = 1 + model.settings.particles = 100 + model.settings.source = src + + # 20 equal-width bins from -1 to 1; cos(45 deg) = 0.7071 falls in [0.7, 0.8) + tally = openmc.Tally() + tally.filters = [ + openmc.MuSurfaceFilter(20), + openmc.SurfaceFilter([xplane]), + ] + tally.scores = ['current'] + model.tallies = [tally] + + model.run(apply_tally_results=True) + current_mu = tally.mean.ravel() + + expected_bin = int((math.cos(math.radians(45.0)) + 1.0) / 0.1) + assert current_mu[expected_bin] == 1.0 + assert current_mu.sum() == 1.0 diff --git a/tests/unit_tests/test_surface_flux.py b/tests/unit_tests/test_surface_flux.py index b51c55fc3f4..4a9cff6d77f 100644 --- a/tests/unit_tests/test_surface_flux.py +++ b/tests/unit_tests/test_surface_flux.py @@ -204,3 +204,103 @@ def test_surface_filter_do_not_tally_virtual_surface_crossing(run_in_tmpdir): # Every particle crosses exit the cube with weight 1, so current = 1.0 assert current_mean.sum() == pytest.approx(1.0, rel=1e-8) + + +def test_surface_flux_rotated_universe(run_in_tmpdir): + """Surface flux uses the surface normal in the root coordinate frame. + + The plane is defined inside a universe that is filled into a cell rotated + by 45 degrees about z, so the normal the surface reports in its own frame + is (1, 0, 0) while its normal in the root frame is 45 degrees away from + the particle direction. The surface-crossing estimator must use the + latter, giving w/|mu| = 1/cos(45 deg). + """ + openmc.reset_auto_ids() + + # Universe split by a plane whose local normal is +x + xplane = openmc.XPlane(0.0) + inner1 = openmc.Cell(region=-xplane) + inner2 = openmc.Cell(region=+xplane) + inner_univ = openmc.Universe(cells=[inner1, inner2]) + + sph = openmc.Sphere(r=10.0, boundary_type='vacuum') + root_cell = openmc.Cell(region=-sph, fill=inner_univ) + root_cell.rotation = (0.0, 0.0, 45.0) + + model = openmc.Model() + model.geometry = openmc.Geometry([root_cell]) + + src = openmc.IndependentSource() + src.space = openmc.stats.Point((-5.0, 0.0, 0.0)) + src.angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0)) + + model.settings.run_mode = 'fixed source' + model.settings.batches = 1 + model.settings.particles = 100 + model.settings.source = src + + flux_tally = openmc.Tally() + flux_tally.filters = [openmc.SurfaceFilter([xplane])] + flux_tally.scores = ['flux'] + + current_tally = openmc.Tally() + current_tally.filters = [openmc.SurfaceFilter([xplane])] + current_tally.scores = ['current'] + + model.tallies = [flux_tally, current_tally] + model.run(apply_tally_results=True) + + mu = math.cos(math.radians(45.0)) + assert flux_tally.mean.flat[0] == pytest.approx(1.0 / mu) + # The net current is direction-signed, so it is unaffected by the rotation + assert current_tally.mean.flat[0] == pytest.approx(1.0) + + +def test_lattice_surface_flux_rotated(run_in_tmpdir): + """Lattice boundary normals are also taken in the root coordinate frame. + + Same setup as test_surface_flux_rotated_universe, except the crossing is a + lattice tile boundary rather than a CSG surface. RectLattice reports its + boundary normal in the lattice's own frame, which here is rotated by 45 + degrees relative to the particle direction. + """ + openmc.reset_auto_ids() + + cell1 = openmc.Cell() + cell2 = openmc.Cell() + cell_outer = openmc.Cell() + univ1 = openmc.Universe(cells=[cell1]) + univ2 = openmc.Universe(cells=[cell2]) + univ_outer = openmc.Universe(cells=[cell_outer]) + + lattice = openmc.RectLattice() + lattice.lower_left = (-2.0, -10.0, -10.0) + lattice.pitch = (2.0, 20.0, 20.0) + lattice.universes = [[[univ1, univ2]]] + lattice.outer = univ_outer + + sph = openmc.Sphere(r=10.0, boundary_type='vacuum') + root_cell = openmc.Cell(region=-sph, fill=lattice) + root_cell.rotation = (0.0, 0.0, 45.0) + + model = openmc.Model() + model.geometry = openmc.Geometry([root_cell]) + + src = openmc.IndependentSource() + src.space = openmc.stats.Point((-5.0, 0.0, 0.0)) + src.angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0)) + + model.settings.run_mode = 'fixed source' + model.settings.batches = 1 + model.settings.particles = 100 + model.settings.source = src + + tally = openmc.Tally() + tally.filters = [openmc.CellFromFilter([cell1]), openmc.CellFilter([cell2])] + tally.scores = ['flux'] + model.tallies = [tally] + + model.run(apply_tally_results=True) + + mu = math.cos(math.radians(45.0)) + assert tally.mean.flat[0] == pytest.approx(1.0 / mu) From ff012c9b584b64c4ce6a5d29f7bf026adf5ad5a3 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 9 Sep 2026 19:15:31 +0000 Subject: [PATCH 2/5] Do not recover surface source half-space across coordinate frames Surface source files store unsigned surface IDs, so FileSource::sample() recovers the signed half-space from the site itself. It did that with site.u.dot(surf.normal(site.r)) > 0.0 but site.r and site.u are in the root coordinate frame while Surface::normal() and Surface::evaluate() work in the local frame of the universe holding the surface. For a surface below the root universe the two frames differ, and the coincidence test guarding the computation is not a reliable filter: a surface passing through a point that happens to satisfy the root-frame surface equation is accepted and then signed from mismatched frames. With a plane inside a universe filled into a cell rotated 135 degrees about z, a particle crossing at the origin along the root-frame +x axis was started in the cell on the wrong side of the plane: its direction in the plane's own frame points into the negative half-space, but the root-frame dot product against the unrotated normal is positive. The transform cannot be recovered from the site alone, since a universe may be filled in several places with different rotations, so restrict the recovery to surfaces whose local frame is known to be the root frame. finalize_geometry() now flags surfaces used by cells outside the root universe, and those fall through to the existing SURFACE_NONE path. That is not a loss of accuracy: the cell search resolves an on-surface point with Surface::sense(), which applies the same direction-versus- normal rule but in the correct local frame. Add unit tests for a surface in a rotated universe, which fails on the current develop code, and for a surface in the root universe, which passes both before and after and pins the recovery that is kept. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_018gUxmUj7G6Aie4azS9uypa --- include/openmc/surface.h | 9 +++ src/geometry_aux.cpp | 11 +++ src/source.cpp | 10 ++- tests/unit_tests/test_source_file.py | 105 +++++++++++++++++++++++++++ 4 files changed, 134 insertions(+), 1 deletion(-) diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 2d8580345a4..0a8a759c89e 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -41,6 +41,15 @@ class Surface { unique_ptr bc_; //!< Boundary condition bool surf_source_ {false}; //!< Activate source banking for the surface? + //! Is this surface used only by cells in the root universe? + //! + //! evaluate() and normal() work in the local coordinate frame of the + //! universe holding the surface. That frame coincides with the root (lab) + //! frame only for surfaces in the root universe, so this flag marks the + //! surfaces for which a lab-frame position may be passed to them directly. + //! Set by finalize_geometry(). + bool root_frame_ {true}; + explicit Surface(pugi::xml_node surf_node); Surface(); diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index a740740c1e6..77f82fbf548 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -282,6 +282,17 @@ void finalize_geometry() // Determine number of nested coordinate levels in the geometry model::n_coord_levels = maximum_levels(model::root_universe); + + // Flag surfaces that are used outside of the root universe. Their local + // coordinate frame does not coincide with the root frame, so a lab-frame + // position cannot be handed to Surface::evaluate() or Surface::normal(). + for (const auto& c : model::cells) { + if (c->universe_ == model::root_universe) + continue; + for (auto token : c->surfaces()) { + model::surfaces[std::abs(token) - 1]->root_frame_ = false; + } + } } //============================================================================== diff --git a/src/source.cpp b/src/source.cpp index f951722cc18..d8528d56105 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -554,11 +554,19 @@ SourceSite FileSource::sample(uint64_t* seed) const // surface containing the source site, determine the signed half-space from // the particle direction. Otherwise, ignore the surface ID and allow the // normal cell search to locate the particle. + // + // The site position and direction are in the root coordinate frame, whereas + // evaluate() and normal() work in the local frame of the universe holding + // the surface. The half-space can therefore only be recovered here for + // surfaces in the root universe; for any other surface the frames differ by + // an unknown transform (a universe may be filled in several places with + // different rotations, so the transform cannot be recovered from the site + // alone) and the surface ID is dropped rather than signed incorrectly. if (site.surf_id != SURFACE_NONE) { auto it = model::surface_map.find(std::abs(site.surf_id)); if (it != model::surface_map.end()) { const auto& surf = *model::surfaces[it->second]; - if (surf.geom_type() == GeometryType::CSG && + if (surf.geom_type() == GeometryType::CSG && surf.root_frame_ && std::abs(surf.evaluate(site.r)) < FP_COINCIDENT) { int surf_id = std::abs(site.surf_id); site.surf_id = diff --git a/tests/unit_tests/test_source_file.py b/tests/unit_tests/test_source_file.py index f19ea74a67a..01debdd10a9 100644 --- a/tests/unit_tests/test_source_file.py +++ b/tests/unit_tests/test_source_file.py @@ -169,3 +169,108 @@ def test_source_file_photon_transport(run_in_tmpdir): # Running OpenMC should succeed model.run() + + +def _rotated_universe_model(angle): + """Plane inside a universe filled into a cell rotated by `angle` about z. + + Particles are born at (-5, 0, 0) travelling along +x and cross the plane + at the origin. The plane's normal in its own frame is +x, but in the root + frame it is rotated by `angle`, so for angle > 90 degrees the particle + enters the negative half-space even though its root-frame direction has a + positive dot product with the unrotated normal. + """ + openmc.reset_auto_ids() + + xplane = openmc.XPlane(0.0, surface_id=7) + inner1 = openmc.Cell(cell_id=1, region=-xplane) + inner2 = openmc.Cell(cell_id=2, region=+xplane) + inner_univ = openmc.Universe(cells=[inner1, inner2]) + + sph = openmc.Sphere(r=10.0, boundary_type='vacuum') + root_cell = openmc.Cell(cell_id=3, region=-sph, fill=inner_univ) + root_cell.rotation = (0.0, 0.0, angle) + + model = openmc.Model() + model.geometry = openmc.Geometry([root_cell]) + + src = openmc.IndependentSource() + src.space = openmc.stats.Point((-5.0, 0.0, 0.0)) + src.angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0)) + + model.settings.run_mode = 'fixed source' + model.settings.batches = 1 + model.settings.particles = 20 + model.settings.source = src + return model, inner1, inner2 + + +def test_surface_source_half_space_rotated_universe(run_in_tmpdir): + """A surface source on a surface in a rotated universe starts in the + correct half-space. + + Surface source files store unsigned surface IDs, and the half-space is + recovered when the file is read back. That recovery compares the root-frame + particle direction against Surface::normal(), which is expressed in the + local frame of the universe holding the surface. With a 135 degree + rotation the two disagree in sign, so the particle used to be started in + the cell on the wrong side of the plane. + """ + model, _, _ = _rotated_universe_model(135.0) + model.settings.surf_source_write = { + 'surface_ids': [7], 'max_particles': 50} + model.run() + + # Read the surface source back in and see which cell the particles start in + model, inner1, inner2 = _rotated_universe_model(135.0) + model.settings.source = openmc.FileSource('surface_source.h5') + tally = openmc.Tally() + tally.filters = [openmc.CellFilter([inner1, inner2])] + tally.scores = ['flux'] + model.tallies = [tally] + model.run(apply_tally_results=True) + + # The particle direction in the plane's own frame points into the negative + # half-space, so all of the track length belongs to the -x cell. + assert tally.mean.flat[0] == pytest.approx(10.0) + assert tally.mean.flat[1] == pytest.approx(0.0) + + +def test_surface_source_half_space_root_universe(run_in_tmpdir): + """The half-space is still recovered for a surface in the root universe.""" + openmc.reset_auto_ids() + + xplane = openmc.XPlane(0.0, surface_id=7) + sph = openmc.Sphere(r=10.0, boundary_type='vacuum') + + def build(): + cell1 = openmc.Cell(cell_id=1, region=-sph & -xplane) + cell2 = openmc.Cell(cell_id=2, region=-sph & +xplane) + model = openmc.Model() + model.geometry = openmc.Geometry([cell1, cell2]) + src = openmc.IndependentSource() + src.space = openmc.stats.Point((-5.0, 0.0, 0.0)) + src.angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0)) + model.settings.run_mode = 'fixed source' + model.settings.batches = 1 + model.settings.particles = 20 + model.settings.source = src + return model, cell1, cell2 + + model, _, _ = build() + model.settings.surf_source_write = { + 'surface_ids': [7], 'max_particles': 50} + model.run() + + openmc.reset_auto_ids() + model, cell1, cell2 = build() + model.settings.source = openmc.FileSource('surface_source.h5') + tally = openmc.Tally() + tally.filters = [openmc.CellFilter([cell1, cell2])] + tally.scores = ['flux'] + model.tallies = [tally] + model.run(apply_tally_results=True) + + # Travelling along +x through the plane, the particle enters the +x cell + assert tally.mean.flat[0] == pytest.approx(0.0) + assert tally.mean.flat[1] == pytest.approx(10.0) From 454eb00814347e484b158aace3ed1ee1f5824f5e Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 9 Sep 2026 19:22:17 +0000 Subject: [PATCH 3/5] Use rotate_to_root() for the ray trace plot surface normal PhongRay::on_intersection() carried its own copy of the loop that walks a surface normal from its local coordinate frame back up to the root frame. That is exactly what rotate_to_root() does, so call it instead of repeating the loop. Pure refactor: the two are the same computation, with the helper's index i standing in for the local lev + 1. The plotter has always had this right; it was the surface tally and surface source paths that did not, so this leaves the helper as the single implementation rather than a third copy of the same rotations. Verified by rendering two solid_raytrace plots of a model whose curved surfaces sit in a rotated and translated universe, before and after the change: the PNGs are byte-identical. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_018gUxmUj7G6Aie4azS9uypa --- src/plot.cpp | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/plot.cpp b/src/plot.cpp index 68ce3c7d167..96711191b6e 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -1825,12 +1825,7 @@ void PhongRay::on_intersection() // Need to apply rotations to find the normal vector in // the base level universe's coordinate system. - for (int lev = surf_level - 1; lev >= 0; --lev) { - if (coord(lev + 1).rotated()) { - const Cell& c {*model::cells[coord(lev).cell()]}; - normal = normal.inverse_rotate(c.rotation_); - } - } + normal = rotate_to_root(*this, surf_level, normal); // use the normal opposed to the ray direction if (normal.dot(u()) > 0.0) { From fafc57c91d15d7191056ffb673d73ceec93c5444 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 9 Sep 2026 20:35:27 +0000 Subject: [PATCH 4/5] Document that Surface::root_frame_ does not cover DAGMC surfaces finalize_geometry() computes the flag by walking the surfaces named in each cell's region, and DAGCell inherits Cell::surfaces() rather than overriding it, so it reports none. A DAGMC surface therefore keeps the default root_frame_ = true even when its universe is nested below the root and carries a transform, which the flag's name does not suggest. Nothing acts on that today: the only reader, FileSource::sample(), tests geom_type() == GeometryType::CSG first, so DAGMC surfaces never reach the flag. Record the limitation at both the declaration and the place it is computed so the next reader does not take it at face value. Comments only; no change in behavior. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_018gUxmUj7G6Aie4azS9uypa --- include/openmc/surface.h | 6 ++++++ src/geometry_aux.cpp | 2 ++ 2 files changed, 8 insertions(+) diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 0a8a759c89e..d6b9ce956fe 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -48,6 +48,12 @@ class Surface { //! frame only for surfaces in the root universe, so this flag marks the //! surfaces for which a lab-frame position may be passed to them directly. //! Set by finalize_geometry(). + //! + //! Only meaningful for CSG surfaces. finalize_geometry() determines it by + //! walking the surfaces named in each cell's region, and DAGCell does not + //! report any (it does not override Cell::surfaces()), so a DAGMC surface + //! keeps the default here even when its universe is nested below the root + //! and transformed. Check geom_type() before relying on this flag. bool root_frame_ {true}; explicit Surface(pugi::xml_node surf_node); diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 77f82fbf548..f8139326b03 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -286,6 +286,8 @@ void finalize_geometry() // Flag surfaces that are used outside of the root universe. Their local // coordinate frame does not coincide with the root frame, so a lab-frame // position cannot be handed to Surface::evaluate() or Surface::normal(). + // Note that DAGCell reports no surfaces, so this leaves DAGMC surfaces at + // their default; see the comment on Surface::root_frame_. for (const auto& c : model::cells) { if (c->universe_ == model::root_universe) continue; From 3db070d84cafa2a91004ff7aa580250bea86d41e Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 10 Sep 2026 18:06:43 +0000 Subject: [PATCH 5/5] Take the crossed surface's coordinate level from the boundary event_cross_surface() read the level as n_coord() - 1. That is the same value as boundary().coord_level() - 1, which is where the surrounding code took it from before, because n_coord() is assigned from boundary().coord_level() two lines above. Reviewers should not have to go and check that, and it is fragile: anything later inserted between the assignment and this line would decouple the two silently. Read the level from the boundary instead, and index that level directly when evaluating the normal rather than going through r_local(), which expands to coord_[n_coord_ - 1] and carried the same hidden dependency. No change in behavior; i_lattice and the evaluated position are the same values as before. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_018gUxmUj7G6Aie4azS9uypa --- src/particle.cpp | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index 14316b08fac..6a8e5cd4504 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -344,13 +344,19 @@ void Particle::event_cross_surface() surface() = boundary().surface(); n_coord() = boundary().coord_level(); - // The surface or lattice being crossed belongs to the universe at the lowest - // coordinate level, so its normal is reported in that level's local frame - // while the particle direction used to score surface tallies lives in the - // root frame. The normal therefore has to be evaluated at the local position - // and rotated up into the root frame, and that has to happen before the - // crossing is carried out, since crossing invalidates the coordinate levels. - int i_surf_level = n_coord() - 1; + // The surface or lattice being crossed belongs to the universe at the + // coordinate level the boundary search found it on, so its normal is + // reported in that level's local frame while the particle direction used to + // score surface tallies lives in the root frame. The normal therefore has to + // be evaluated at the local position and rotated up into the root frame, and + // that has to happen before the crossing is carried out, since crossing + // invalidates the coordinate levels. + // + // Take the level from the boundary rather than from n_coord(). The two are + // equal here because of the assignment just above, but reading it from the + // boundary keeps this independent of that, and matches where the level came + // from originally. + int i_surf_level = boundary().coord_level() - 1; if (boundary().lattice_translation()[0] != 0 || boundary().lattice_translation()[1] != 0 || @@ -388,7 +394,7 @@ void Particle::event_cross_surface() // Determine the surface normal in the root frame before crossing Direction normal; if (!model::active_surface_tallies.empty()) { - normal = surf.normal(r_local()); + normal = surf.normal(coord(i_surf_level).r()); normal = rotate_to_root(*this, i_surf_level, normal / normal.norm()); }