Skip to content

Let a mesh surface crossing change more than one index - #4096

Open
GuySten wants to merge 1 commit into
openmc-dev:developfrom
GuySten:mesh-multiaxis-crossing
Open

Let a mesh surface crossing change more than one index#4096
GuySten wants to merge 1 commit into
openmc-dev:developfrom
GuySten:mesh-multiaxis-crossing

Conversation

@GuySten

@GuySten GuySten commented Aug 31, 2026

Copy link
Copy Markdown
Contributor

Let a mesh surface crossing change more than one index

Pure refactor of StructuredMesh. No behaviour change, no results change, no
regolds. Second of the pieces split out of #3857, and independent of it.

Problem

StructuredMesh::raytrace_mesh assumes that crossing a surface changes exactly
one mesh index, and that the distance to the next surface on the other axes is
unaffected:

// Update cell and calculate distance to next surface in k-direction.
// The two other directions are still valid!
ijk[k] = distances[k].next_index;
distances[k] = distance_to_grid_boundary(ijk, k, local_r, u, traveled_distance);

in_mesh = ((ijk[k] >= 1) && (ijk[k] <= shape_[k]));

That is true of every mesh currently in the code, but it is a property of those
meshes rather than of structured meshes in general, and the assumption is spread
across three separate places: MeshDistance::next_index being a single int,
the hard-coded refresh of only distances[k], and the inline bounds check.

Two more things are entangled with it. Periodic index wrapping is done inside
each mesh's distance_to_grid_boundary, so next_index is pre-sanitized:

MeshDistance(sanitize_phi(ijk[i] + 1), true, find_phi_crossing(...))

which means the value stored is not "where I am going" but "where I am going,
already folded". And the logic for re-entering the mesh from outside sits inline
in the middle of raytrace_mesh.

Changes

Three changes, each exercised by code in this PR.

  • MeshDistance::next_index (an int) becomes offset (a MeshIndex). A
    crossing now says how each index changes rather than naming a single
    destination. The existing meshes each set one component, so ijk[i] + 1
    becomes {1, 0, 0} and so on.
  • virtual void sanitize_index(MeshIndex&), called once after the index is
    advanced. CylindricalMesh and SphericalMesh override it with the
    sanitize_phi / sanitize_theta calls they previously folded into
    distance_to_grid_boundary. Default is a no-op.
  • distance_to_mesh(...) extracts the re-entry block, moved verbatim, out of
    the middle of raytrace_mesh.

Plus MESH_MAX_AXES, replacing the hard-coded 3 in MeshIndex, shape_ and
the per-axis distance array -- three independent literals that all had to agree.

Why this is behaviour-preserving

  • offset defaults to {0, 0, 0} and the old code set next_index = ijk[i]
    before any early return, so a crossing that does not happen leaves the index
    unchanged either way.
  • sanitize_index applies the same sanitize_phi / sanitize_theta to the
    same axis, just after the index is advanced instead of before it is stored.
  • Exactly distances[k] is still refreshed after a crossing.
  • distance_to_mesh is the previous expression, unmoved.

Value independent of #3857

The periodic wrapping that cylindrical and spherical meshes do gets a name and
one place to live rather than being folded into each mesh's distance
calculation, where it made next_index mean two things at once. And
distance_to_mesh is now a separately readable and testable function rather
than a block in the middle of the longest function in the file.

Performance

Benchmarked against a synthetic stress case -- a 50³ mesh surface tally, 2000 particles, single-threaded -- chosen so that mesh tallying dominates transport time. Cylindrical and spherical meshes show no resolvable change. A regular mesh shows a small increase, bounded at roughly 12% of transport in that configuration by the least favourable reading of the data and under 3% by the robust one; four repeats per arm is not enough to pin it down further. Real models, where transport is not dominated by a fine mesh surface tally, would see a fraction of that.

Testing

No results change, so the existing mesh regression and unit tests cover this
unchanged. Nothing new is executed and no allocation is added on the transport
path -- MeshIndex is std::array<int, 3>, so MeshDistance grows by two ints
and stays a plain aggregate.

No new data members and no allocation. MeshDistance grows by two ints and
stays a plain aggregate.

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)
  • I have added tests that prove my fix is effective or that my feature works (if applicable)

@GuySten
GuySten marked this pull request as ready for review August 31, 2026 21:38
@GuySten
GuySten requested a review from paulromano August 31, 2026 21:38
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant