Skip to content

Fix coordinate frame of surface normals in surface tallies and surface sources - #4119

Open
GuySten wants to merge 4 commits into
openmc-dev:developfrom
GuySten:claude/elegant-mayer-nn4gpe
Open

Fix coordinate frame of surface normals in surface tallies and surface sources#4119
GuySten wants to merge 4 commits into
openmc-dev:developfrom
GuySten:claude/elegant-mayer-nn4gpe

Conversation

@GuySten

@GuySten GuySten commented Sep 9, 2026

Copy link
Copy Markdown
Contributor

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, and RectLattice/HexLattice::get_normal() likewise return tile boundary normals in the lattice's frame. Both were dotted with Particle::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:

quantity before correct
surface-crossing flux estimator 1.0000 1.4142 (41% low)
MuSurfaceFilter bin [0.9, 1.0) [0.7, 0.8)
lattice-crossing flux 1.0000 1.4142

Net 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 new rotate_to_root() helper that walks up the coordinate levels undoing each cell rotation. score_surface_tally() records that root-frame normal on the particle so MuSurfaceFilter reads it instead of recomputing: the filter could not have applied the correction itself, because by the time filters run cross_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-bounds model::surfaces[-1] access, since on a lattice crossing the surface token is SURFACE_NONE and surface_index() returned -1.

This also corrects nested DAGMC geometry, where a DAGMCUniverse sits below the root (as in tests/regression_tests/dagmc/universes). DAGSurface::normal() hands the position to MOAB's get_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, since r_local() and r() are then the same and rotate_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 from site.u.dot(surf.normal(site.r)). Again site.r/site.u are root-frame while normal() and evaluate() are local-frame, and the |evaluate(site.r)| < FP_COINCIDENT guard 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):

-x cell +x cell
before 0.0 10.0
after 10.0 0.0

This one is not fixed by rotating the normal. There is no coordinate stack inside sample(), and building one is circular: exhaustive_find_cell() consumes p.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 existing SURFACE_NONE path. That costs no accuracy: the cell search resolves an on-surface point with Surface::sense(), which applies the identical direction-versus-normal rule but in the correct local frame, because find_cell_inner() passes it r_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 calls rotate_to_root(), so the helper is the single implementation rather than a third copy of the same rotations. Pure refactor, 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.

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 DAGCell inherits Cell::surfaces() rather than overriding it, so it reports none. A DAGMC surface therefore keeps root_frame_ = true even when nested and transformed, which the name does not suggest. Nothing acts on that — the only reader tests geom_type() == GeometryType::CSG first — 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 develop as well as this branch:

test on develop here
test_surface_flux_rotated_universe fails passes
test_lattice_surface_flux_rotated fails passes
test_musurface_rotated_universe fails passes
test_surface_source_half_space_rotated_universe fails passes
test_surface_source_half_space_root_universe passes passes

The 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

  • I have performed a self-review of my own code
  • I have run clang-format (version 18) on any C++ source files (if applicable)
  • I have followed the style guidelines for Python source files (if applicable)
  • I have made corresponding changes to the documentation (if applicable) — not applicable, no user-facing behaviour is documented differently
  • I have added tests that prove my fix is effective or that my feature works (if applicable)

🤖 Generated with Claude Code

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 GuySten added the Bugs label Sep 9, 2026
@GuySten GuySten changed the title Fix coordinate frame of surface normals in surface tallies and surface sources- #47 Fix coordinate frame of surface normals in surface tallies and surface sources Sep 9, 2026
@GuySten
GuySten marked this pull request as ready for review September 9, 2026 21:59
@GuySten
GuySten requested a review from pshriwise as a code owner September 9, 2026 21:59

@pshriwise pshriwise left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for this @GuySten! I'll have more of a look later, but one thing jumped out at me for now.

Comment thread src/particle.cpp
// 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;

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Something to think about: The previous evaluation of i_lattice is set using the boundary property, but here we're now using n_coord() -1. I'd have to look to see if these are consistent in the context of when this code is called.

@pshriwise pshriwise left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Correcting review status.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants