Fix coordinate frame of surface normals in surface tallies and surface sources - #4119
Open
GuySten wants to merge 4 commits into
Open
Fix coordinate frame of surface normals in surface tallies and surface sources#4119GuySten wants to merge 4 commits into
GuySten wants to merge 4 commits into
Conversation
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 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018gUxmUj7G6Aie4azS9uypa
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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018gUxmUj7G6Aie4azS9uypa
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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018gUxmUj7G6Aie4azS9uypa
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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018gUxmUj7G6Aie4azS9uypa
GuySten
marked this pull request as ready for review
September 9, 2026 21:59
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Description
Surface and lattice boundary normals are reported in the local coordinate frame of the universe that owns them, but particle directions are stored in the root (lab) frame. Several places compared the two directly, which is wrong whenever the surface sits under a cell carrying a rotation.
1. Surface tally normals (
3aebbac)Surface::normal()returns the outward normal in the surface's own frame, andRectLattice/HexLattice::get_normal()likewise return tile boundary normals in the lattice's frame. Both were dotted withParticle::u(), the direction at coordinate level 0.Two further problems at the same call sites: the normal was evaluated at
r()(root-frame position) rather than the local position, which is wrong even without a rotation for any position-dependent normal (sphere, cylinder, cone, torus, quadric); and it was computed after the crossing was carried out, so for a periodic boundary it used the position after translation to the partner surface.Measured with a plane inside a universe filled into a cell rotated 45 degrees about z, particles crossing along the lab +x axis:
MuSurfaceFilterbinNet current is unaffected, since only the sign of the cosine matters there.
The fix evaluates the normal at
r_local(), before the crossing, and rotates it into the root frame with a newrotate_to_root()helper that walks up the coordinate levels undoing each cell rotation.score_surface_tally()records that root-frame normal on the particle soMuSurfaceFilterreads it instead of recomputing: the filter could not have applied the correction itself, because by the time filters runcross_surface()has re-descended the coordinate stack and the surface's level is no longer identifiable. As a side effect this removes an out-of-boundsmodel::surfaces[-1]access, since on a lattice crossing the surface token isSURFACE_NONEandsurface_index()returned -1.This also corrects nested DAGMC geometry, where a
DAGMCUniversesits below the root (as intests/regression_tests/dagmc/universes).DAGSurface::normal()hands the position to MOAB'sget_angle()to find the facet normal, so a root-frame position there was not merely mis-rotated — it was a lookup at a point that need not lie on the surface at all. For a root-level DAGMC universe the behavior is unchanged, sincer_local()andr()are then the same androtate_to_root()is a no-op.2. Surface source half-space (
ff012c9)Surface source files store unsigned surface IDs, so
FileSource::sample()recovers the signed half-space fromsite.u.dot(surf.normal(site.r)). Againsite.r/site.uare root-frame whilenormal()andevaluate()are local-frame, and the|evaluate(site.r)| < FP_COINCIDENTguard does not filter this out: 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 in a universe rotated 135 degrees about z, a particle crossing the origin along root-frame +x was started in the cell on the wrong side of the plane (track-length flux, cm per source particle):
This one is not fixed by rotating the normal. There is no coordinate stack inside
sample(), and building one is circular:exhaustive_find_cell()consumesp.surface()to disambiguate on-surface points, which is the value being computed. Nor can the transform be recovered from the site alone, since a universe may be filled in several places with different rotations.So the recovery is restricted to surfaces whose local frame is the root frame.
finalize_geometry()flags surfaces used by cells outside the root universe, and those fall through to the existingSURFACE_NONEpath. That costs no accuracy: the cell search resolves an on-surface point withSurface::sense(), which applies the identical direction-versus-normal rule but in the correct local frame, becausefind_cell_inner()passes itr_local()/u_local().3. Plotter refactor (
454eb00)PhongRay::on_intersection()already carried its own copy of the walk-up loop; the plotter has always had this right. It now callsrotate_to_root(), so the helper is the single implementation rather than a third copy of the same rotations. Pure refactor, verified by rendering twosolid_raytraceplots of a model whose curved surfaces sit in a rotated and translated universe, before and after the change: the PNGs are byte-identical.4. Comment on the scope of
root_frame_(fafc57c)The flag added in (2) is computed by walking the surfaces named in each cell's region, and
DAGCellinheritsCell::surfaces()rather than overriding it, so it reports none. A DAGMC surface therefore keepsroot_frame_ = trueeven when nested and transformed, which the name does not suggest. Nothing acts on that — the only reader testsgeom_type() == GeometryType::CSGfirst — but the limitation is now recorded at the declaration and at the place it is computed. Comments only.Tests
Five unit tests were added. Each was checked against
developas well as this branch:developtest_surface_flux_rotated_universetest_lattice_surface_flux_rotatedtest_musurface_rotated_universetest_surface_source_half_space_rotated_universetest_surface_source_half_space_root_universeThe last one passes both ways by design: it pins the root-universe recovery path that is deliberately kept, so this cannot later be "simplified" by deleting the reconstruction outright.
Also run locally: all 12 C++ Catch2 suites, and 17 geometry/tally/filter/source unit-test files compared before and after, with identical pass/fail counts apart from the new tests.
One coverage gap worth flagging for reviewers: no test in the repository combines DAGMC with surface tallies, so the DAGMC reasoning in (1) rests on reading the code rather than on a passing test. The DAGMC CI job confirms nothing existing regressed, but it does not exercise the changed path.
No tracking issue is associated with this change.
Checklist
🤖 Generated with Claude Code