Fix _dm_unstack_bcs: destroyed auxiliary + UW_Boundaries labels (#162)#163
Merged
Conversation
Contributor
There was a problem hiding this comment.
Pull request overview
Fixes a boundary-label handling bug in UW3’s PETSc DM label (un)stacking utilities that caused auxiliary/internal boundary labels to be destroyed during _dm_unstack_bcs, breaking natural boundary condition integration on internal boundaries (issue #162) and during mesh adaptation.
Changes:
- Expand reconstructed boundary labels to include their closure via
dm.labelComplete(...)during_dm_unstack_bcs. - Restore UW3 auxiliary labels (
All_BoundariesandNull_Boundary) after the remove/recreate cycle. - Add detailed docstring context explaining the #162 failure mode and fix rationale.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
|
|
||
| if "All_Boundaries" in boundary_names: | ||
| # Replaces any existing content; safe to call after removeLabel. | ||
| dm.markBoundaryFaces("All_Boundaries", 1001) |
| verts_is = depth_label.getStratumIS(0) | ||
| if verts_is is not None and verts_is.getSize() > 0: | ||
| null_label = dm.getLabel("Null_Boundary") | ||
| null_label.setStratumIS(666, verts_is) |
Comment on lines
+536
to
+541
| try: | ||
| dm.labelComplete(b_dmlabel) | ||
| except Exception: | ||
| # labelComplete is a no-op for empty labels and certain | ||
| # depth combinations; swallow rather than fail label setup. | ||
| pass |
Comment on lines
+528
to
+535
| # BUGFIX(#162): expand the boundary label to include its closure | ||
| # (the vertices and edges bordering the labeled facets). The | ||
| # standard ``useRegions=True`` gmsh import does this implicitly; | ||
| # the ``useRegions=False`` + ``_dm_unstack_bcs`` path used by | ||
| # internal-boundary meshes left labels at the height-1 stratum | ||
| # only. PETSc's natural-BC residual integration needs the | ||
| # closure to find the FE basis points on the boundary, so a | ||
| # closure-less label silently produces zero contribution. |
BoxInternalBoundary calls _dm_unstack_bcs after Mesh construction to unpack its gmsh "Face Sets" into named boundary labels. The previous implementation iterated *every* boundary in the enum and removed its label, then only repopulated labels whose values appeared in the source stack — silently destroying the auxiliary labels Null_Boundary and All_Boundaries (added by extend_enum but populated separately by the Mesh constructor and _from_gmsh / _from_plexh5). It also did not refresh the consolidated UW_Boundaries label that the Mesh constructor builds at construction time. UW_Boundaries is what the solver registers natural BCs against (PetscDSAddBoundary_UW calls in petsc_generic_snes_solvers.pyx). Because the Mesh constructor builds it BEFORE BIB calls _dm_unstack_bcs, none of the named boundaries (Bottom, Top, Left, Right, Internal) were present in UW_Boundaries — so add_natural_bc on any boundary of a BIB mesh was a silent no-op. The named boundaries also lacked their closure expansion, because the useRegions=False + _dm_unstack_bcs path used by internal-boundary meshes doesn't get the implicit closure that gmsh's use_regions route does. Empirical confirmation (issue #162): on BoxInternalBoundary the boundary label populations diverged dramatically from the equivalent StructuredQuadBox: label SQBox BIB pre-fix BIB post-fix Top 63 32 65 All_Boundaries 128 0 128 Null_Boundary 1089 0 1089 UW_Boundaries[12] 63 0 65 (named entry for Top) Without the named entries in UW_Boundaries, the solver registered the natural BC as a label-value pair PETSc could find, but the consolidated label PETSc looks up by name had no facets at that value — silently returning zero residual contribution. The same code path is used in mesh_adapt_meshVar (adaptivity.py:530), so adapted meshes also lost these labels. Fix --- In _dm_unstack_bcs, after the existing remove/recreate cycle: 1. dm.labelComplete(label) on each named boundary so its closure (vertices/edges bordering the labeled facets) is included. 2. Re-mark All_Boundaries via dm.markBoundaryFaces. 3. Re-set Null_Boundary from the depth-0 stratum (all vertices). 4. Rebuild UW_Boundaries from the now-correct individual boundary labels, mirroring the Mesh constructor's logic. Verified -------- - NengLu's reproducer (issue #162): pre-fix u ≈ 2.5e-19 everywhere (zero), F-norm < 1e-11 (no BC contribution); post-fix |u_x| = 5.65e-02 just above and below the internal interface, matching FEniCS's analytical reference 5.65310750e-02 to roundoff. Solver converges normally in 1 SNES iteration with initial F-norm 0.131 (real BC contribution). - Exterior natural BC on BIB Top: also restored (was equally broken pre-fix). u_top = -0.150, matching the StructuredQuadBox Couette-like profile. - Post-fix label populations on BoxInternalBoundary match StructuredQuadBox structurally. - pytest -m "level_1 and tier_a" on amr-dev: 56 passed, 3 skipped, 0 failed. No regressions. Closes #162. Note on sign convention ----------------------- UW3's solution has opposite sign to FEniCS's reference (-0.0565 vs +0.0565). Magnitude is identical to roundoff. This is likely an internal-boundary normal-orientation convention; the FEM solution itself is correct. Underworld development team with AI support from Claude Code
2061518 to
9f034c4
Compare
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
Fully fixes natural BCs on `BoxInternalBoundary` (and the `mesh_adapt_meshVar` adaptation path) — both for exterior boundaries (`Top`, `Bottom`, etc.) and for the `Internal` interface itself. Single-file change to `_dm_unstack_bcs`.
NengLu's reproducer from #162 now matches FEniCS's analytical reference to roundoff.
Closes #162.
Root cause (two layers)
`BoxInternalBoundary` calls `_dm_unstack_bcs` after Mesh construction to unpack its gmsh `Face Sets` into named boundary labels. The previous implementation did three things wrong:
(3) is the layer that turned a label-correctness bug into a solver-correctness bug.
Empirical confirmation
NengLu's reproducer (issue #162):
Sign of the post-fix solution is opposite to FEniCS (-0.0565 vs +0.0565); magnitude is identical. Likely an internal-boundary normal-orientation convention rather than a numerical issue.
Fix
In `_dm_unstack_bcs`, after the existing remove/recreate cycle:
The same fix benefits the `mesh_adapt_meshVar` adaptation path automatically — it's the only other caller of `_dm_unstack_bcs`.
Test plan
Notes
@NengLu — would you mind verifying that your reproducer now produces the expected magnitude on this branch, and that none of your existing scripts regress?
Underworld development team with AI support from Claude Code