From 884789f31ec8461ef4e5f79318fa8819fcee44f1 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Fri, 10 Jul 2026 07:22:29 +1000 Subject: [PATCH] fix(mesh): add mesh.canonical_normal accessor for partition-safe internal-boundary normals (#327) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `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 --- .../discretisation/discretisation_mesh.py | 55 ++++++++++++++++ tests/test_0502_boundary_integrals.py | 63 ++++++++++++------- 2 files changed, 95 insertions(+), 23 deletions(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 8859464a..1c499d2b 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -2709,6 +2709,61 @@ def Gamma_P1(self): self._update_projected_normals() return self._projected_normals.sym + def canonical_normal(self, boundary_name): + 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. + + The primary caller is code that needs a **partition-safe** normal + on an *internal* boundary — see :issue:`327`. On an internal + boundary at a partition seam, PETSc's per-quadrature ``petsc_n[]`` + (surfacing as :attr:`Gamma` / :attr:`Gamma_N`) is derived from + ``support[0]`` of the DMPlex facet closure, which is + partition-dependent; different ranks disagree on which cell is + "support[0]" for the one seam facet, and the outward normal of + that facet flips sign. A signed integral of ``Gamma[k]`` is then + wrong by O(seam-facets / total-facets). The analytic normal + returned here is partition-independent and does not touch + ``petsc_n``, so it sidesteps the defect entirely for the mesh + classes that know their internal-boundary geometry + (:class:`BoxInternalBoundary`, :class:`AnnulusInternalBoundary`, + :class:`SphericalShellInternalBoundary`). + + Parameters + ---------- + boundary_name : str + Name of the boundary label to look up (case-sensitive; must + match one of :attr:`boundaries`). + + Returns + ------- + 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``. + + See Also + -------- + Gamma : the raw per-quadrature normal — use for external + boundaries; may be partition-dependent on internal seams. + Gamma_N : normalised :attr:`Gamma`. + Gamma_P1 : projected P1 normals, useful for curved external + boundaries. + """ + 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 + # =================================================================== # Bounding surfaces — per-surface tangent-slip + restore. # See docs/developer/design/boundary-slip-strategy.md. SEPARATE from diff --git a/tests/test_0502_boundary_integrals.py b/tests/test_0502_boundary_integrals.py index 5f7d323f..a67e040f 100644 --- a/tests/test_0502_boundary_integrals.py +++ b/tests/test_0502_boundary_integrals.py @@ -159,24 +159,22 @@ def test_bd_integral_internal_coordinate_fn(): assert abs(value - 0.5) < 0.01, f"Expected 0.5, got {value}" -# TODO(BUG): internal-boundary facet-normal orientation is rank-dependent at -# partition seams: at np2 one seam facet contributes with flipped sign, so -# |integral of n_y| = 1 - 2/32 = 0.9375. Scalar integrands (length, -# coordinate functions) are exact in parallel; only signed-normal integrands -# are affected. Found while unskipping after the BF-13 constructor fix -# (2026-07 audit) — separate defect, not covered by BF-13. -@pytest.mark.skipif( - uw.mpi.size > 1, - reason="Internal-boundary normal orientation is rank-dependent at partition seams (see TODO(BUG) above)", -) +# 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. def test_bd_integral_internal_normal_ny(): - """Integrate n_y along internal boundary at y=0.5. - The internal boundary has normals pointing in +y or -y direction, - so integrating n_y should give +1 or -1 (length 1 boundary).""" + """Integrate n_y along internal boundary at y=0.5 using the analytic + canonical normal. The internal boundary has normals pointing in +y or + -y direction, so integrating n_y should give +1 or -1 (length 1 + boundary).""" mesh_internal, _, _ = _get_internal_mesh() - Gamma = mesh_internal.Gamma - n_y = Gamma[1] + normal = mesh_internal.canonical_normal("Internal") + n_y = normal[1] bd_int = uw.maths.BdIntegral(mesh_internal, fn=n_y, boundary="Internal") value = bd_int.evaluate() @@ -200,17 +198,14 @@ def test_bd_integral_internal_normal_nx(): assert abs(value) < 0.01, f"Expected ~0, got {value}" -@pytest.mark.skipif( - uw.mpi.size > 1, - reason="Internal-boundary normal orientation is rank-dependent at partition seams (see TODO(BUG) above)", -) def test_bd_integral_internal_normal_weighted(): - """Integrate x * n_y along internal boundary at y=0.5. - int_0^1 x * n_y dx = n_y * 0.5. Since |n_y| = 1, result should be ~0.5.""" + """Integrate x * n_y along internal boundary at y=0.5 using the analytic + canonical normal. int_0^1 x * n_y dx = n_y * 0.5. Since |n_y| = 1, + result should be ~0.5.""" mesh_internal, x_i, _ = _get_internal_mesh() - Gamma = mesh_internal.Gamma - n_y = Gamma[1] + normal = mesh_internal.canonical_normal("Internal") + n_y = normal[1] bd_int = uw.maths.BdIntegral(mesh_internal, fn=x_i * n_y, boundary="Internal") value = bd_int.evaluate() @@ -218,6 +213,28 @@ def test_bd_integral_internal_normal_weighted(): assert abs(abs(value) - 0.5) < 0.01, f"Expected |value| = 0.5, got {value}" +def test_bd_integral_internal_canonical_normal_off_grid_zint(): + """Regression for the failing partition-through-boundary case from #327. + + With ``zintCoord=0.55`` (off-grid), the mpirun -n 2 partition seam runs + through the internal boundary and one seam facet's ``petsc_n[]`` flips + sign: ``∫ Gamma[1] dS`` returns 0.9375 = 1 − 2/32 instead of 1.0. The + analytic ``canonical_normal("Internal")`` returned by the mesh factory + does not touch ``petsc_n[]`` and gives the exact value regardless of + partition.""" + mesh_off = BoxInternalBoundary( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + cellSize=1.0/32.0, zintCoord=0.55, simplex=True, + ) + # Need at least one variable so BdIntegral has a section to integrate against + uw.discretisation.MeshVariable("T_off", mesh_off, 1, degree=2) + + n_y = mesh_off.canonical_normal("Internal")[1] + val = uw.maths.BdIntegral(mesh_off, fn=n_y, boundary="Internal").evaluate() + assert abs(abs(val) - 1.0) < 1e-6, ( + f"canonical_normal internal integral should be exactly ±1, got {val}") + + def test_bd_integral_internal_does_not_affect_external(): """External boundaries should still work on the internal-boundary mesh."""