fix(mesh): add mesh.canonical_normal accessor for partition-safe internal-boundary normals (#327)#359
fix(mesh): add mesh.canonical_normal accessor for partition-safe internal-boundary normals (#327)#359lmoresi wants to merge 1 commit into
Conversation
…rnal-boundary normals (#327) `mesh.Gamma[i]` JIT-substitutes `petsc_n[i]`, PETSc's per-quadrature outward-normal-from-`support[0]` in `DMPlexComputeBdIntegral`. On an internal boundary at a partition seam, different ranks disagree on which cell is `support[0]` for the one seam facet, so the outward normal flips sign there and a signed integral like `∫ n_y dS` is silently off by O(seam facets / total facets) in parallel — the exact `1 − 2/32 = 0.9375` signature #327 reports. The mesh factory Enums (`BoxInternalBoundary.boundary_normals_{2,3}D`, `AnnulusInternalBoundary.boundary_normals`, `SphericalShellInternalBoundary.boundary_normals`) already declare the analytic outward-pointing normal for each boundary, but nothing consumed them. Wire that up: `Mesh.canonical_normal(boundary_name)` returns the declared sympy Matrix (or None if not declared), giving callers a partition-independent alternative to `Gamma` for meshes that know their internal-boundary geometry analytically. The `test_0502` "internal_normal_ny" / "internal_normal_weighted" tests are un-skipped (they were serial-only with a TODO(BUG) marker citing #327) and now use the canonical accessor. A new `test_bd_integral_internal_canonical_normal_off_grid_zint` regression constructs the failing `zint=0.55` configuration that reliably picks a partition seam through the internal boundary at np=2 (from the sweep in the /remote-control session's #327 investigation). 23/23 pass in `tests/test_0502_boundary_integrals.py` at both np=1 and np=2. Not fixed: `Gamma[i]` itself remains partition-dependent on internal boundaries — users writing code that reads `Gamma[i]` on an internal boundary still get the wrong answer at np>1 unless they switch to `canonical_normal`. The full fix (canonicalise `support[0]` at construction, or expose `DMPlexOrientLabel` via a Cython wrapper) is larger and out of scope for this PR — see the #327 investigation notes. Underworld development team with AI support from Claude Code
There was a problem hiding this comment.
Pull request overview
Adds a partition-safe way to obtain consistent internal-boundary normals by exposing mesh-factory analytic normals via a new Mesh.canonical_normal(boundary_name) accessor, and updates boundary-integral tests to use it for signed-normal integrals affected by #327.
Changes:
- Add
Mesh.canonical_normal(boundary_name)to return mesh-factory declared analytic normals (avoidspetsc_n[]/ partition seam sign flips). - Unskip and update internal-boundary signed-normal tests to use
canonical_normal("Internal"). - Add a regression test for the off-grid
zintCoord=0.55seam-crossing configuration.
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated 4 comments.
| File | Description |
|---|---|
src/underworld3/discretisation/discretisation_mesh.py |
Introduces the canonical_normal() accessor and documents intended use for partition-safe internal-boundary normals. |
tests/test_0502_boundary_integrals.py |
Updates tests to use canonical_normal() for signed internal-normal integrals and adds an off-grid regression reproducer. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| # The `mesh.Gamma` normal on an *internal* boundary is derived from petsc_n[] | ||
| # which uses DMPlex support[0] — partition-dependent at seam facets, so a | ||
| # signed-normal integral like ∫ n_y dS is silently off by O(seam facets / | ||
| # total facets) in parallel (issue #327). These tests use the analytic | ||
| # `canonical_normal` accessor, which returns the mesh factory's declared | ||
| # outward normal for the boundary (a partition-independent sympy expression); | ||
| # see the accessor's docstring in discretisation_mesh.py. |
| r"""Analytic outward-pointing normal for a boundary, or ``None`` | ||
| if no analytic normal was declared for that boundary. | ||
|
|
||
| Sourced from the mesh factory's ``boundary_normals`` Enum: for | ||
| axis-aligned box boundaries this is a constant sympy Matrix, for | ||
| annulus / spherical-shell radial boundaries it is the analytic | ||
| radial unit vector, and so on. |
| sympy.Matrix or None | ||
| Row matrix of length :attr:`cdim` giving the outward-pointing | ||
| normal, or ``None`` if this mesh factory did not declare | ||
| an analytic normal for ``boundary_name``. |
| bn = getattr(self, "boundary_normals", None) | ||
| if bn is None: | ||
| return None | ||
| try: | ||
| member = bn[boundary_name] | ||
| except (KeyError, AttributeError): | ||
| return None | ||
| value = getattr(member, "value", member) | ||
| return value |
|
Superseded by #361, which contains this branch's commit (the canonical_normal accessor) and builds the full unified interface on top: mesh.Gamma resolves per boundary (external → exact petsc_n[]; internal → the declared analytic normal), with deformation invalidation and an np≥4 label-completion deadlock fix. See #361 for details. Underworld development team with AI support from Claude Code |
Partial fix for #327. Wires up the mesh factories' declared analytic normals as a partition-safe alternative to `Gamma[i]` on internal boundaries.
Background
`mesh.Gamma[i]` JIT-substitutes `petsc_n[i]`, PETSc's per-quadrature
outward-normal-from-`support[0]` in `DMPlexComputeBdIntegral`. On an internal
boundary at a partition seam, different ranks disagree on which cell is
`support[0]` for the one seam facet, so the outward normal flips sign there and
a signed integral like `∫ n_y dS` is silently off by O(seam facets / total facets)
in parallel.
Reproducer (from the /remote-control session's #327 investigation):
```python
mesh = BoxInternalBoundary(minCoords=(0,0), maxCoords=(1,1),
cellSize=1/32, zintCoord=0.55, simplex=True)
np=1: ∫ Gamma[1] dS = 1.000000
np=2: ∫ Gamma[1] dS = 0.937500 ← wrong, one seam facet flipped
```
At `zintCoord=0.5` (grid-aligned) gmsh happens to route the np=2 seam around
the internal boundary and the bug is invisible; move `zint` off-grid and it
reappears reliably.
Change
The mesh factory Enums (`BoxInternalBoundary.boundary_normals_{2,3}D`,
`AnnulusInternalBoundary.boundary_normals`, `SphericalShellInternalBoundary`,
etc.) already declare the analytic outward-pointing normal for each boundary,
but nothing consumed them. This PR wires them up:
`Mesh.canonical_normal(boundary_name)` — returns the declared sympy Matrix
outward normal for a boundary, or `None` if this mesh factory did not declare
one. Because it is a symbolic expression built from mesh coordinates it does
not touch `petsc_n[]`, so any signed-normal integral evaluated with it is
partition-independent for the meshes whose internal-boundary geometry the
factory knows (Box, Annulus, SphericalShell, ...).
Tests
in `tests/test_0502_boundary_integrals.py` were serial-only with a
`TODO(BUG)` marker citing Internal-facet normal orientation is rank-dependent at partition seams (vector surface integrals wrong at np>1) #327. Un-skipped and switched to use
`mesh.canonical_normal("Internal")`.
`zint=0.55` seam-crossing configuration (which is the one that reliably
reproduces the bug on Linux gmsh) and asserts the canonical normal gives
exactly `|·| = 1.0` at np≥1.
What this does not fix
`Gamma[i]` itself remains partition-dependent on internal boundaries. Any
user code that reads `Gamma[i]` on an internal boundary at np>1 still gets
the wrong answer unless they switch to `canonical_normal`. The full fix is
larger — either canonicalising `support[0]` at construction (Cython for
`dm.setCone`) or exposing `DMPlexOrientLabel` with a proper depth-encoded
label constructor. Both routes explored in the #327 investigation notes;
neither is a session slice.
Underworld development team with AI support from Claude Code