Skip to content

Commit 28a06bb

Browse files
lmoresilmoresi
committed
Add PETSc patch for ghost facet double-counting in boundary assembly
DMPlexComputeBdResidual_Internal and DMPlexComputeBdIntegral do not exclude ghost facets (SF leaves) from the facet IS, unlike the volume residual code which checks the ghost label. For internal boundaries at partition junctions, this causes shared facets to be integrated on multiple ranks, producing O(1) errors after MPI summation. The patch filters ghost facets using ISDifference, following the same pattern used by DMPlexComputeBdFaceGeomFVM for FVM face flux. It applies cleanly to PETSc v3.18.0 through v3.24.5. Includes: - petsc-custom/patches/plexfem-ghost-facet-fix.patch - build-petsc.sh integration (auto-applies after clone) - MR description for upstream PETSc submission Upstream MR branch pushed to gitlab.com/lmoresi/petsc (fix/plexfem-ghost-facet-boundary-residual) Underworld development team with AI support from Claude Code
1 parent c91fedd commit 28a06bb

File tree

3 files changed

+164
-0
lines changed

3 files changed

+164
-0
lines changed

‎petsc-custom/build-petsc.sh‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -57,6 +57,27 @@ clone_petsc() {
5757
echo "Clone complete."
5858
}
5959

60+
apply_patches() {
61+
echo "Applying UW3 patches to PETSc..."
62+
cd "$PETSC_DIR"
63+
64+
# Fix ghost facet double-counting in boundary residual/integral assembly.
65+
# Without this, internal boundary natural BCs and BdIntegral produce
66+
# incorrect results in parallel (shared facets integrated on multiple ranks).
67+
# Submitted upstream: [TODO: add PETSc MR link when available]
68+
local patch="${SCRIPT_DIR}/patches/plexfem-ghost-facet-fix.patch"
69+
if [ -f "$patch" ]; then
70+
if git apply --check "$patch" 2>/dev/null; then
71+
git apply "$patch"
72+
echo " Applied: plexfem-ghost-facet-fix.patch"
73+
else
74+
echo " Skipped: plexfem-ghost-facet-fix.patch (already applied or conflict)"
75+
fi
76+
fi
77+
78+
echo "Patches complete."
79+
}
80+
6081
configure_petsc() {
6182
echo "Configuring PETSc with AMR tools..."
6283
cd "$PETSC_DIR"
@@ -147,6 +168,7 @@ show_help() {
147168
echo " build Build PETSc"
148169
echo " test Run PETSc tests"
149170
echo " petsc4py Build and install petsc4py"
171+
echo " patch Apply UW3 patches to PETSc source"
150172
echo " clean Remove PETSc directory"
151173
echo " help Show this help"
152174
}
@@ -155,6 +177,7 @@ show_help() {
155177
case "${1:-all}" in
156178
all)
157179
clone_petsc
180+
apply_patches
158181
configure_petsc
159182
build_petsc
160183
build_petsc4py
@@ -175,6 +198,9 @@ case "${1:-all}" in
175198
build)
176199
build_petsc
177200
;;
201+
patch)
202+
apply_patches
203+
;;
178204
test)
179205
test_petsc
180206
;;

‎petsc-custom/patches/PETSC_MR_DESCRIPTION.md‎

Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
# PETSc Merge Request: Ghost Facet Fix
2+
3+
## Submission details
4+
5+
- **Fork**: https://gitlab.com/lmoresi/petsc
6+
- **Branch**: `fix/plexfem-ghost-facet-boundary-residual`
7+
- **Target**: `petsc/petsc` → `release`
8+
- **Create MR**: https://gitlab.com/lmoresi/petsc/-/merge_requests/new?merge_request%5Bsource_branch%5D=fix%2Fplexfem-ghost-facet-boundary-residual
9+
10+
## Title
11+
12+
DMPlex: filter ghost facets from boundary residual and integral assembly
13+
14+
## Description
15+
16+
### Problem
17+
18+
`DMPlexComputeBdResidual_Internal` and `DMPlexComputeBdIntegral` construct
19+
`facetIS` as all facets at depth `dim-1` but do not exclude ghost facets
20+
(point SF leaves). The volume residual code in `DMPlexComputeResidualByKey`
21+
checks the `"ghost"` label before calling `DMPlexVecSetClosure` (line 5355),
22+
but the boundary residual path has no equivalent check.
23+
24+
For **external** boundaries this is benign: ghost facets on the domain
25+
boundary still have `support[0]` pointing to a valid local cell, and
26+
`DMLocalToGlobal(ADD_VALUES)` correctly resolves ownership. But for
27+
**internal** boundaries — facets in the mesh interior that carry a
28+
`DMLabel` value and a `DM_BC_NATURAL` condition — facets at partition
29+
junctions are present on multiple ranks. Each rank independently
30+
integrates the boundary flux via `PetscFEIntegrateBdResidual`, and
31+
`DMLocalToGlobal(ADD_VALUES)` sums all contributions, causing
32+
double-counting.
33+
34+
The same issue affects `DMPlexComputeBdIntegral`: ghost boundary facets
35+
are included in the per-face integral sum, which is then `MPI_Allreduce`d,
36+
again double-counting shared facets.
37+
38+
### Use case: internal boundaries in geodynamics
39+
40+
[Underworld3](https://github.com/underworldcode/underworld3) is a
41+
Python/PETSc geodynamics framework that uses DMPlex FEM throughout.
42+
Geophysical models commonly require traction (natural) boundary conditions
43+
on **internal surfaces** — for example, a fault plane or a material
44+
interface embedded within the mesh. These are set up by labeling interior
45+
facets with a `DMLabel` and registering `DM_BC_NATURAL` via
46+
`PetscDSAddBoundary` / `PetscWeakFormAddBdResidual`.
47+
48+
This works correctly in serial but produces O(1) errors whenever the
49+
labeled internal boundary is split across partition boundaries, because
50+
shared facets are integrated on multiple ranks.
51+
52+
### Fix
53+
54+
After obtaining `facetIS` from `DMLabelGetStratumIS(depthLabel, dim-1, ...)`,
55+
filter out SF leaves using `ISDifference`. This ensures each boundary facet
56+
is processed by exactly one rank (the owner).
57+
58+
This follows the same pattern already used in `DMPlexComputeBdFaceGeomFVM`
59+
(plexfem.c, around line 1087) which calls `PetscFindInt(face, nleaves, leaves, &loc)`
60+
to skip ghost faces during FVM face flux computation.
61+
62+
### Testing
63+
64+
Tested with Underworld3 Stokes solver using `DM_BC_NATURAL` on an internal
65+
boundary (annular mesh with labeled interior circle). Results match serial
66+
to 10+ significant figures across all processor counts tested (1, 2, 4, 8).
67+
Before the fix, any partition that splits the internal boundary produces
68+
O(1) errors from double-counted boundary contributions.
69+
70+
The patch applies cleanly to all PETSc releases from v3.18.0 through
71+
v3.24.5 — the affected code has been structurally unchanged across these
72+
versions.
73+
74+
## Reference examples
75+
76+
- Stokes Green's function with internal boundary natural BCs:
77+
https://github.com/underworldcode/underworld3/blob/development/docs/examples/fluid_mechanics/advanced/Ex_Stokes_Annulus_Kernels.py
78+
- Cartesian variant:
79+
https://github.com/underworldcode/underworld3/blob/development/docs/examples/fluid_mechanics/advanced/Ex_Stokes_Cartesian_Kernels.py

‎petsc-custom/patches/plexfem-ghost-facet-fix.patch‎

Lines changed: 59 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,59 @@
1+
diff --git a/src/dm/impls/plex/plexfem.c b/src/dm/impls/plex/plexfem.c
2+
index 0299f569070..26e6432585d 100644
3+
--- a/src/dm/impls/plex/plexfem.c
4+
+++ b/src/dm/impls/plex/plexfem.c
5+
@@ -2874,6 +2874,26 @@ PetscErrorCode DMPlexComputeBdIntegral(DM dm, Vec X, DMLabel label, PetscInt num
6+
PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
7+
PetscCall(DMGetDimension(dm, &dim));
8+
PetscCall(DMLabelGetStratumIS(depthLabel, dim - 1, &facetIS));
9+
+ /* Filter out ghost facets (SF leaves) so that each boundary facet is only
10+
+ counted on one rank. Without this, shared facets at partition boundaries
11+
+ are integrated on multiple ranks, causing double-counting after MPI sum. */
12+
+ if (facetIS) {
13+
+ PetscSF sf;
14+
+ PetscInt nleaves;
15+
+ const PetscInt *leaves;
16+
+
17+
+ PetscCall(DMGetPointSF(dm, &sf));
18+
+ PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL));
19+
+ if (nleaves > 0 && leaves) {
20+
+ IS leafIS, ownedFacetIS;
21+
+
22+
+ PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nleaves, leaves, PETSC_USE_POINTER, &leafIS));
23+
+ PetscCall(ISDifference(facetIS, leafIS, &ownedFacetIS));
24+
+ PetscCall(ISDestroy(&leafIS));
25+
+ PetscCall(ISDestroy(&facetIS));
26+
+ facetIS = ownedFacetIS;
27+
+ }
28+
+ }
29+
PetscCall(DMGetLocalSection(dm, &section));
30+
PetscCall(PetscSectionGetNumFields(section, &Nf));
31+
/* Get local solution with boundary values */
32+
@@ -5113,6 +5133,27 @@ static PetscErrorCode DMPlexComputeBdResidual_Internal(DM dm, Vec locX, Vec locX
33+
PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
34+
PetscCall(DMGetDimension(dm, &dim));
35+
PetscCall(DMLabelGetStratumIS(depthLabel, dim - 1, &facetIS));
36+
+ /* Filter out ghost facets (SF leaves) so that boundary residual contributions
37+
+ from shared facets are only assembled on the owning rank. Without this,
38+
+ internal boundary natural BCs at partition junctions get double-counted
39+
+ because LocalToGlobal with ADD_VALUES sums contributions from all ranks. */
40+
+ if (facetIS) {
41+
+ PetscSF sf;
42+
+ PetscInt nleaves;
43+
+ const PetscInt *leaves;
44+
+
45+
+ PetscCall(DMGetPointSF(dm, &sf));
46+
+ PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL));
47+
+ if (nleaves > 0 && leaves) {
48+
+ IS leafIS, ownedFacetIS;
49+
+
50+
+ PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nleaves, leaves, PETSC_USE_POINTER, &leafIS));
51+
+ PetscCall(ISDifference(facetIS, leafIS, &ownedFacetIS));
52+
+ PetscCall(ISDestroy(&leafIS));
53+
+ PetscCall(ISDestroy(&facetIS));
54+
+ facetIS = ownedFacetIS;
55+
+ }
56+
+ }
57+
PetscCall(PetscDSGetNumBoundary(prob, &numBd));
58+
for (bd = 0; bd < numBd; ++bd) {
59+
PetscWeakForm wf;

0 commit comments

Comments
 (0)