Unify boundary-normal interface behind mesh.Gamma (#327) + np=4 label deadlock fix#361
Open
lmoresi wants to merge 4 commits into
Open
Unify boundary-normal interface behind mesh.Gamma (#327) + np=4 label deadlock fix#361lmoresi wants to merge 4 commits into
lmoresi wants to merge 4 commits 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
…match) getNonEmptyStratumValuesIS() is rank-local, but the per-value loop makes collective calls (labelComplete). When a rank's partition touches no facet of some boundary — routine at np>=4 on BoxInternalBoundary — ranks make different numbers of collective calls and the run deadlocks in mesh construction (stacks: one rank in labelComplete, others in markBoundaryFaces). np=2 escaped because both halves touch every boundary. Iterate the allgathered union of stratum values so every rank makes the same sequence of collective calls; guard the local setStratumIS to values live on this rank. Verified: BoxInternalBoundary construction and the full 0502 suite now pass at np=4 (24/24, previously hung indefinitely). Underworld development team with AI support from Claude Code
mesh.Gamma is now the single user-facing boundary-normal symbol. The consumers that know their boundary (uw.maths.BdIntegral, natural BCs via add_condition) resolve it per boundary through the new Mesh._resolve_boundary_normals: - external boundaries: unchanged — Gamma compiles to PETSc's exact per-quadrature outward unit normal petsc_n[] (zero solver risk); - internal boundaries: petsc_n[] orientation follows the partition-dependent DMPlex support[0] (issue #327), so Gamma components are substituted with the mesh factory's declared analytic normal (canonical_normal), partition-independent by construction; - internal boundary with no declared normal: loud RuntimeError — there is no orientation convention for an arbitrary internal surface until DMPlexOrientLabel supports open surfaces (pending upstream). Internal-boundary detection (_boundary_is_internal): any labeled facet with support size 2, MAX-allreduced and cached per boundary. Supporting changes: - Convention: declared internal normals point from region Inner to region Outer. 2D BoxInternalBoundary Internal flipped [0,-1] -> [0,1] to match 3D (+z) and the radial (unit_e_0) annulus/spherical convention. - AnnulusInternalBoundary: attach the existing boundary_normals Enum (was defined but never attached); SphericalShellInternalBoundary: add and attach one. TODO(BUG) noted on Box external entries, which are inward-pointing. - Gamma_N docstring: deprecated alias of Gamma (numerically identical). Gamma_P1 docstring: deprecated, with diagnosis (off-kernel petsc_n fallback + sign cancellation on internal facets -> sub-unit vectors). - SNES_Stokes_Constrained.add_constraint_bc default normal now boundary_normal(boundary) (oriented-then-averaged per-boundary field) instead of the deprecated global Gamma_P1. - tests(0502): internal-normal tests use plain mesh.Gamma with tight 1e-6 assertions and the Inner->Outer (+y) sign; back-compat test for the canonical_normal accessor; #327 off-grid zint regression kept. Verified: 0502 suite 24/24 at np=1, np=2, np=4; natural-BC/Nitsche/ constrained-freeslip regression trio 23 passed. Underworld development team with AI support from Claude Code
Contributor
There was a problem hiding this comment.
Pull request overview
This PR makes mesh.Gamma the single user-facing boundary-normal symbol across the codebase by resolving it per-boundary in boundary-integrand consumers (exact PETSc quadrature normals on external boundaries; declared analytic normals on internal boundaries to avoid rank-dependent orientation at partition seams, #327). It also fixes a pre-existing MPI deadlock in _dm_unstack_bcs at np>=4 by ensuring every rank performs the same sequence of collective calls.
Changes:
- Add
Mesh.canonical_normal()plus_boundary_is_internal()and_resolve_boundary_normals()to transparently substitute analytic normals for internal-boundarymesh.Gammausage in BdIntegral and natural BCs. - Attach/adjust mesh-factory
boundary_normalsdeclarations (including updatedBoxInternalBoundaryinternal normal convention) and update solver defaults to use per-boundary assembled normals for constraint BCs. - Fix
_dm_unstack_bcsdeadlock by iterating the globally unioned set of label values; update/extend regression tests around internal-boundary normals in parallel.
Reviewed changes
Copilot reviewed 9 out of 9 changed files in this pull request and generated 3 comments.
Show a summary per file
| File | Description |
|---|---|
| tests/test_0502_boundary_integrals.py | Unskips and strengthens internal-boundary normal integral tests; adds regression for off-grid zintCoord seam case. |
| src/underworld3/systems/solvers.py | Changes default constraint-BC normal to mesh.boundary_normal(boundary) (per-boundary assembled normal). |
| src/underworld3/meshing/spherical.py | Adds boundary_normals Enum attachment to spherical internal-boundary meshes for canonical/resolved normals. |
| src/underworld3/meshing/cartesian.py | Flips 2D internal boundary declared normal to +y (Inner→Outer convention) and adds an explicit TODO note about external normal direction. |
| src/underworld3/meshing/annulus.py | Attaches boundary_normals Enum to annulus internal-boundary meshes for canonical/resolved normals. |
| src/underworld3/discretisation/discretisation_mesh.py | Introduces canonical_normal, internal-boundary detection/cache, and Gamma resolution logic; refactors boundary-facet enumeration helper. |
| src/underworld3/cython/petsc_maths.pyx | Resolves mesh.Gamma per boundary before JIT-compiling boundary integrands. |
| src/underworld3/cython/petsc_generic_snes_solvers.pyx | Resolves mesh.Gamma per boundary for Neumann/natural-BC expressions. |
| src/underworld3/adaptivity.py | Fixes np>=4 deadlock in _dm_unstack_bcs by iterating the allgathered union of label values. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+2727
to
+2728
| r"""Analytic outward-pointing normal for a boundary, or ``None`` | ||
| if no analytic normal was declared for that boundary. |
Comment on lines
+2758
to
+2761
| 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``. |
|
|
||
| The integrand may reference the outward unit normal via ``mesh.Gamma_N`` | ||
| (components map to ``petsc_n[0], petsc_n[1], ...`` in the generated C code). | ||
| The integrand may reference the outward unit normal via ``mesh.Gamma`` |
The factory-declared analytic normal (canonical_normal, consumed when mesh.Gamma resolves on an internal boundary) describes the ORIGINAL factory geometry. After mesh.deform() / adapt() the internal surface may have warped, and silently integrating the stale declaration would be wrong. - nuke_coords_and_rebuild bumps a _geometry_version counter (0 at construction; every later call is a genuine coordinate change — deform, adapt, re-extract and direct coordinate writes all funnel through it). - boundary_normals is now a property whose setter stamps the declaration with the current geometry version. - canonical_normal raises RuntimeError when the declaration predates the current geometry, telling the user to re-declare (mesh.boundary_normals = mesh.boundary_normals when the expression is still valid, e.g. radial normals after a radius-preserving remesh) or to use an explicit normal expression. External boundaries are unaffected: Gamma -> petsc_n[] is exact on the deformed geometry and never consults the declaration. Test: deform a BoxInternalBoundary with a displacement vanishing on the internal plane; Gamma resolution raises; re-declaration restores the exact +1 integral. Serial-only until #360 (deform crashes at np>1). Verified: 0502 25 passed serial, 24 passed + 1 skip at np=2 and np=4; natural-BC/Nitsche/constrained regression trio 23 passed. Underworld development team with AI support from Claude Code
Member
Author
|
Added 2f98dec: declared analytic normals are invalidated by mesh deformation. The declaration describes the factory geometry, so after Mechanism: New regression test (serial-only until #360: Underworld development team with AI support from Claude Code |
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.
Summary
One user-facing normal symbol.
mesh.Gammanow works on any boundary; the four alternative interfaces (Gamma_N,Gamma_P1,boundary_normal(b),canonical_normal(b)) are no longer needed in user code (all retained for back-compat; the first two are documented as deprecated).Consumers that know their boundary —
uw.maths.BdIntegraland natural BCs (add_condition/add_natural_bc) — resolveGammaper boundary via the newMesh._resolve_boundary_normals:Gammacompiles to PETSc's exact per-quadrature outward unit normalpetsc_n[]. This path is byte-identical to before, so zero solver risk.petsc_n[]orientation follows the partition-dependent DMPlexsupport[0](Internal-facet normal orientation is rank-dependent at partition seams (vector surface integrals wrong at np>1) #327:mpirun -n 2gave 0.9375 instead of 1.0 for a unit signed-normal integral).Gammacomponents are substituted with the mesh factory's declared analytic normal (canonical_normal), which is partition-independent by construction.RuntimeError. There is no orientation convention for an arbitrary internal surface untilDMPlexOrientLabelhandles open surfaces (pending upstream with PETSc); failing loudly beats integrating a partition-dependent sign. (Noteboundary_normal()silently returns zeros there — it skips support-size-2 facets — which is why it is not a fallback.)Internal-boundary detection: any labeled facet with support size 2, MAX-allreduced across ranks and cached per boundary.
Orientation convention (behaviour change)
Declared internal-boundary normals point from region
Innerto regionOuter:BoxInternalBoundary:Internalflipped[0,-1]→[0,1](now consistent with 3D[0,0,1])AnnulusInternalBoundary/SphericalShellInternalBoundary: radially outwardunit_e_0(the annulus Enum existed but was never attached to the mesh; the spherical factory had none — both now attached)Code that consumed the old 2D
[0,-1]declaration viacanonical_normalwill see a sign flip.Second commit: np≥4 deadlock in
_dm_unstack_bcs(pre-existing)While validating at np=4, every
BoxInternalBoundaryconstruction deadlocked — ondevelopmenttoo, i.e. pre-existing and independent of this change._dm_unstack_bcsiterated the rank-localgetNonEmptyStratumValuesIS()while making collectivelabelCompletecalls per value; ranks whose partition misses a boundary make fewer collective calls → deadlock (np=2 escaped because both halves touch all boundaries). Fixed by iterating the allgathered union of values, with local stratum access guarded. Kept as its own commit for cherry-pick ease.Validation
test_0502_boundary_integrals.py: 24/24 at np=1, np=2, np=4 (np=4 previously hung indefinitely; now 5s)test_1000_poissonNaturalBC(external natural BCs through the new resolver),test_1060_nitsche_freeslip(mesh.Gammainadd_natural_bc),test_1061_constrained_freeslip(newboundary_normal(boundary)default) — 23 passedmesh.Gammainternal signed-normal integrals with 1e-6 assertions and the Inner→Outer sign;canonical_normalaccessor back-compat; the Internal-facet normal orientation is rank-dependent at partition seams (vector surface integrals wrong at np>1) #327 off-gridzintCoord=0.55partition-seam regressionRelationship to PR #359
This branch contains #359's commit (
canonical_normalaccessor) and supersedes it — the accessor alone still required users to pick the right interface. Suggest closing #359 in favour of this PR (leaving that to the maintainers).Deferred (flagged, not in this PR)
Gamma_N/Gamma_P1(not mechanical — curved-boundary guides should point atboundary_normal)TODO(BUG)incartesian.py: Box externalboundary_normalsentries are inward-pointing whilecanonical_normaldocuments outward — needs a consumer audit before flipping_assemble_boundary_normalparallel seam ADD-reduce TODODMPlexOrientLabelopen-surface support still pending upstreamUnderworld development team with AI support from Claude Code