From 0561eaee21deaf114bd34646e79061e98ebcf657 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Thu, 12 Mar 2026 12:16:53 +1100 Subject: [PATCH 1/2] petsc-custom: add internal-boundary ownership patch, mpi test, and docs --- CHANGES.md | 11 ++ petsc-custom/README.md | 9 + petsc-custom/build-petsc.sh | 23 ++- ...DESCRIPTION_INTERNAL_BOUNDARY_OWNERSHIP.md | 134 +++++++++++++++ ...xfem-internal-boundary-ownership-fix.patch | 162 ++++++++++++++++++ ...est_0765_internal_boundary_integral_mpi.py | 61 +++++++ 6 files changed, 392 insertions(+), 8 deletions(-) create mode 100644 petsc-custom/patches/PETSC_MR_DESCRIPTION_INTERNAL_BOUNDARY_OWNERSHIP.md create mode 100644 petsc-custom/patches/plexfem-internal-boundary-ownership-fix.patch create mode 100644 tests/parallel/test_0765_internal_boundary_integral_mpi.py diff --git a/CHANGES.md b/CHANGES.md index 6c7c2132e..936e5efe2 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,5 +1,16 @@ # CHANGES: Underworld3 +## 2026-03-12 + + - Fixed PETSc DMPlex internal-boundary MPI rank-dependence in custom PETSc patch workflow. + - Added `plexfem-internal-boundary-ownership-fix.patch`: + - ghost-facet ownership filtering in boundary integral / residual / Jacobian paths + - part-consistent residual assembly (`support[key.part]`) with support-size guards + - Added MPI regression test: + - `tests/parallel/test_0765_internal_boundary_integral_mpi.py` + - Validation benchmark: + - `Ex_Stokes_Kramer_latest.py` (`case1`, natural BC) now shows stable velocity L2 across `np=1,2,4,8` (no `np=8` branch split). + ## 2025-12-21 - PETSc 3.24 compatibility verified (conda-forge petsc 3.24.2 works correctly) diff --git a/petsc-custom/README.md b/petsc-custom/README.md index 9598e1022..9ac447749 100644 --- a/petsc-custom/README.md +++ b/petsc-custom/README.md @@ -39,6 +39,15 @@ cd petsc-custom ./build-petsc.sh clean # Remove everything ``` +## PETSc Patch Selection + +- `build-petsc.sh` applies PETSc boundary-assembly patches from `petsc-custom/patches`. +- Preferred patch: + - `plexfem-internal-boundary-ownership-fix.patch` +- Legacy fallback (only if preferred patch is absent): + - `plexfem-ghost-facet-fix.patch` +- Do not apply both patches together. + ## What Gets Installed The build downloads and compiles: diff --git a/petsc-custom/build-petsc.sh b/petsc-custom/build-petsc.sh index ff8c94f7e..f2ea51802 100755 --- a/petsc-custom/build-petsc.sh +++ b/petsc-custom/build-petsc.sh @@ -61,14 +61,21 @@ apply_patches() { echo "Applying UW3 patches to PETSc..." cd "$PETSC_DIR" - # Fix ghost facet double-counting in boundary residual/integral assembly. - # Without this, internal boundary natural BCs and BdIntegral produce - # incorrect results in parallel (shared facets integrated on multiple ranks). - # Submitted upstream: [TODO: add PETSc MR link when available] - local patch="${SCRIPT_DIR}/patches/plexfem-ghost-facet-fix.patch" - if [ -f "$patch" ]; then - if git apply --check "$patch" 2>/dev/null; then - git apply "$patch" + # Internal-boundary ownership + part-consistent assembly fix in plexfem.c. + # Supersedes the older ghost-facet-only patch. + local patch_new="${SCRIPT_DIR}/patches/plexfem-internal-boundary-ownership-fix.patch" + local patch_old="${SCRIPT_DIR}/patches/plexfem-ghost-facet-fix.patch" + + if [ -f "$patch_new" ]; then + if git apply --check "$patch_new" 2>/dev/null; then + git apply "$patch_new" + echo " Applied: plexfem-internal-boundary-ownership-fix.patch" + else + echo " Skipped: plexfem-internal-boundary-ownership-fix.patch (already applied or conflict)" + fi + elif [ -f "$patch_old" ]; then + if git apply --check "$patch_old" 2>/dev/null; then + git apply "$patch_old" echo " Applied: plexfem-ghost-facet-fix.patch" else echo " Skipped: plexfem-ghost-facet-fix.patch (already applied or conflict)" diff --git a/petsc-custom/patches/PETSC_MR_DESCRIPTION_INTERNAL_BOUNDARY_OWNERSHIP.md b/petsc-custom/patches/PETSC_MR_DESCRIPTION_INTERNAL_BOUNDARY_OWNERSHIP.md new file mode 100644 index 000000000..fb529fc47 --- /dev/null +++ b/petsc-custom/patches/PETSC_MR_DESCRIPTION_INTERNAL_BOUNDARY_OWNERSHIP.md @@ -0,0 +1,134 @@ +# PETSc Patch: Internal Boundary Ownership Consistency + +## Patch file + +- `plexfem-internal-boundary-ownership-fix.patch` + +## Scope + +File touched in PETSc: + +- `src/dm/impls/plex/plexfem.c` + +## Problem + +Internal-boundary contributions in parallel were rank-dependent for natural BC / boundary integral workflows. +Symptoms included: + +- Different `v_l2`/`p_l2` across MPI sizes for the same problem setup. +- Different boundary callback totals between ranks/partitions. + +Root causes: + +1. Boundary facet sets were not consistently filtered to owned facets in all boundary assembly paths. +2. `DMPlexComputeBdResidualSingleByKey` always used `support[0]` during closure gather/scatter, which is not part-consistent when `key.part != 0`. +3. Part filtering for boundary points needed to guard against insufficient support for `part > 0`. + +## Fix Summary + +The patch makes boundary assembly ownership-consistent and part-consistent: + +1. Filter SF leaf (ghost) facets from `facetIS` using `ISDifference(...)` in: + - `DMPlexComputeBdIntegral` + - `DMPlexComputeBdResidual_Internal` + - `DMPlexComputeBdJacobian_Internal` + +2. In `DMPlexComputeBdResidualSingleByKey`: + - For `key.part > 0`, filter points with `supportSize <= key.part`. + - Use `support[key.part]` (not hardcoded `support[0]`) for: + - local closure reads (`locX`, `locX_t`) + - auxiliary closure lookup + - local residual insertion (`DMPlexVecSetClosure`) + - Add `PetscCheck(supportSize > key.part, ...)` guards before use. + +## Why this resolves the bug + +Each physical boundary facet is assembled once (owned rank only), and each part uses the correct adjacent cell side. This removes over/under-assembly and residual/Jacobian side mismatches that previously depended on partition topology. + +## Validation + +Validated with Underworld3 annulus internal-boundary cases: + +- `Ex_Stokes_Kramer_latest.py` with natural BC, `np=1,2,4,8`: + - stable velocity L2 around `2.32345982e-03` + - no branch split at `np=8` + +- Boundary callback totals (`bdres`/`bdjac`) match between `np=7` and `np=8` after patch. + +- A UW MPI regression test was added: + - `tests/parallel/test_0765_internal_boundary_integral_mpi.py` + - checks annulus internal-boundary circumference in parallel against analytic value. + +## Benchmark Script And Case + +- Benchmark script (GitHub): + - https://github.com/gthyagi/UW3_Annulus_Spherical_Benchmarks/blob/main/benchmarks/annulus/Ex_Stokes_Kramer_latest.py + +- Recommended case for reproducing/validating this fix: + - `-uw_case case1 -uw_bc_type natural` + +- Example runs: + +```bash +# from UW3_Annulus_Spherical_Benchmarks/benchmarks/annulus +PY=/path/to/uw/.pixi/envs/amr-dev/bin/python +SCRIPT=Ex_Stokes_Kramer_latest.py + +for n in 1 2 4 8; do + mpirun -np "$n" "$PY" "$SCRIPT" -uw_case case1 -uw_bc_type natural +done +``` + +## Patch Application Matrix + +Use `git apply --check` before applying: + +```bash +git apply --check /path/to/patch.patch +``` + +### Case A: old patch is not applied + +- Recommended and simplest path. +- Apply only: + - `plexfem-internal-boundary-ownership-fix.patch` + +### Case B: old patch is already applied + +- Do not stack both full patches. +- `plexfem-internal-boundary-ownership-fix.patch` usually overlaps with the old patch and will fail `--check`. +- Recommended options: + 1. Revert/restore `plexfem.c` to clean PETSc source, then apply only the new patch. + 2. Or create/use an incremental patch from old -> new (if you need sequential layering). + +## Verification Commands + +### 1) UW MPI regression test + +```bash +mpirun -np 4 python -m pytest --with-mpi tests/parallel/test_0765_internal_boundary_integral_mpi.py -q +``` + +Expected: +- `2 passed` + +### 2) Benchmark validation + +```bash +# from UW3_Annulus_Spherical_Benchmarks/benchmarks/annulus +for n in 1 2 4 8; do + mpirun -np "$n" "$PY" "$SCRIPT" -uw_case case1 -uw_bc_type natural +done +``` + +Expected: +- `Relative velocity L2 error` should be stable across `np=1,2,4,8` +- Typical value around `2.32345982e-03` for this case + +## Notes + +- This patch supersedes the narrower `plexfem-ghost-facet-fix.patch` by also addressing part-consistent residual assembly and Jacobian ownership filtering. +- Apply **only one** of these two patches. + - Preferred: `plexfem-internal-boundary-ownership-fix.patch` + - Legacy fallback: `plexfem-ghost-facet-fix.patch` +- Do **not** apply both patches together (overlapping hunks / duplicate logic). diff --git a/petsc-custom/patches/plexfem-internal-boundary-ownership-fix.patch b/petsc-custom/patches/plexfem-internal-boundary-ownership-fix.patch new file mode 100644 index 000000000..10d7c429a --- /dev/null +++ b/petsc-custom/patches/plexfem-internal-boundary-ownership-fix.patch @@ -0,0 +1,162 @@ +diff --git a/src/dm/impls/plex/plexfem.c b/src/dm/impls/plex/plexfem.c +index 0299f569070..1f07f49f60f 100644 +--- a/src/dm/impls/plex/plexfem.c ++++ b/src/dm/impls/plex/plexfem.c +@@ -2874,6 +2874,24 @@ PetscErrorCode DMPlexComputeBdIntegral(DM dm, Vec X, DMLabel label, PetscInt num + PetscCall(DMPlexGetDepthLabel(dm, &depthLabel)); + PetscCall(DMGetDimension(dm, &dim)); + PetscCall(DMLabelGetStratumIS(depthLabel, dim - 1, &facetIS)); ++ /* Filter out ghost facets (SF leaves) so each facet contributes once */ ++ if (facetIS) { ++ PetscSF sf; ++ PetscInt nleaves; ++ const PetscInt *leaves; ++ ++ PetscCall(DMGetPointSF(dm, &sf)); ++ PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL)); ++ if (nleaves > 0 && leaves) { ++ IS leafIS, ownedFacetIS; ++ ++ PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nleaves, leaves, PETSC_USE_POINTER, &leafIS)); ++ PetscCall(ISDifference(facetIS, leafIS, &ownedFacetIS)); ++ PetscCall(ISDestroy(&leafIS)); ++ PetscCall(ISDestroy(&facetIS)); ++ facetIS = ownedFacetIS; ++ } ++ } + PetscCall(DMGetLocalSection(dm, §ion)); + PetscCall(PetscSectionGetNumFields(section, &Nf)); + /* Get local solution with boundary values */ +@@ -4961,6 +4979,28 @@ PetscErrorCode DMPlexComputeBdResidualSingleByKey(DM dm, PetscWeakForm wf, Petsc + PetscCall(ISDestroy(&pointIS)); + pointIS = isectIS; + } ++ if (key.part > 0) { ++ IS filteredIS = NULL; ++ const PetscInt *allPoints; ++ PetscInt *keptPoints = NULL; ++ PetscInt nAllFaces, nKept = 0, fidx; ++ ++ PetscCall(ISGetLocalSize(pointIS, &nAllFaces)); ++ PetscCall(ISGetIndices(pointIS, &allPoints)); ++ PetscCall(PetscMalloc1(nAllFaces, &keptPoints)); ++ for (fidx = 0; fidx < nAllFaces; ++fidx) { ++ const PetscInt point = allPoints[fidx]; ++ PetscInt supportSize; ++ ++ PetscCall(DMPlexGetSupportSize(dm, point, &supportSize)); ++ if (supportSize <= key.part) continue; ++ keptPoints[nKept++] = point; ++ } ++ PetscCall(ISRestoreIndices(pointIS, &allPoints)); ++ PetscCall(ISDestroy(&pointIS)); ++ PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nKept, keptPoints, PETSC_OWN_POINTER, &filteredIS)); ++ pointIS = filteredIS; ++ } + PetscCall(ISGetLocalSize(pointIS, &numFaces)); + PetscCall(ISGetIndices(pointIS, &points)); + PetscCall(PetscMalloc4(numFaces * totDim, &u, (locX_t ? (size_t)numFaces * totDim : 0), &u_t, numFaces * totDim, &elemVec, (locA ? (size_t)numFaces * totDimAux : 0), &a)); +@@ -4978,21 +5018,24 @@ PetscErrorCode DMPlexComputeBdResidualSingleByKey(DM dm, PetscWeakForm wf, Petsc + for (face = 0; face < numFaces; ++face) { + const PetscInt point = points[face], *support; + PetscScalar *x = NULL; +- PetscInt i; ++ PetscInt i, supportSize, cell; + + PetscCall(DMPlexGetSupport(dm, point, &support)); +- PetscCall(DMPlexVecGetClosure(plex, section, locX, support[0], NULL, &x)); ++ PetscCall(DMPlexGetSupportSize(dm, point, &supportSize)); ++ PetscCheck(supportSize > key.part, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Boundary point %" PetscInt_FMT " has support size %" PetscInt_FMT " but requested part %" PetscInt_FMT, point, supportSize, key.part); ++ cell = support[key.part]; ++ PetscCall(DMPlexVecGetClosure(plex, section, locX, cell, NULL, &x)); + for (i = 0; i < totDim; ++i) u[face * totDim + i] = x[i]; +- PetscCall(DMPlexVecRestoreClosure(plex, section, locX, support[0], NULL, &x)); ++ PetscCall(DMPlexVecRestoreClosure(plex, section, locX, cell, NULL, &x)); + if (locX_t) { +- PetscCall(DMPlexVecGetClosure(plex, section, locX_t, support[0], NULL, &x)); ++ PetscCall(DMPlexVecGetClosure(plex, section, locX_t, cell, NULL, &x)); + for (i = 0; i < totDim; ++i) u_t[face * totDim + i] = x[i]; +- PetscCall(DMPlexVecRestoreClosure(plex, section, locX_t, support[0], NULL, &x)); ++ PetscCall(DMPlexVecRestoreClosure(plex, section, locX_t, cell, NULL, &x)); + } + if (locA) { + PetscInt subp; + +- PetscCall(DMGetEnclosurePoint(plexA, dm, encAux, support[0], &subp)); ++ PetscCall(DMGetEnclosurePoint(plexA, dm, encAux, cell, &subp)); + PetscCall(DMPlexVecGetClosure(plexA, sectionAux, locA, subp, NULL, &x)); + for (i = 0; i < totDimAux; ++i) a[face * totDimAux + i] = x[i]; + PetscCall(DMPlexVecRestoreClosure(plexA, sectionAux, locA, subp, NULL, &x)); +@@ -5028,10 +5071,13 @@ PetscErrorCode DMPlexComputeBdResidualSingleByKey(DM dm, PetscWeakForm wf, Petsc + } + for (face = 0; face < numFaces; ++face) { + const PetscInt point = points[face], *support; ++ PetscInt supportSize; + + if (mesh->printFEM > 1) PetscCall(DMPrintCellVector(point, name, totDim, &elemVec[face * totDim])); + PetscCall(DMPlexGetSupport(plex, point, &support)); +- PetscCall(DMPlexVecSetClosure(plex, NULL, locF, support[0], &elemVec[face * totDim], ADD_ALL_VALUES)); ++ PetscCall(DMPlexGetSupportSize(dm, point, &supportSize)); ++ PetscCheck(supportSize > key.part, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Boundary point %" PetscInt_FMT " has support size %" PetscInt_FMT " but requested part %" PetscInt_FMT, point, supportSize, key.part); ++ PetscCall(DMPlexVecSetClosure(plex, NULL, locF, support[key.part], &elemVec[face * totDim], ADD_ALL_VALUES)); + } + PetscCall(DMSNESRestoreFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom)); + PetscCall(PetscQuadratureDestroy(&qGeom)); +@@ -5113,6 +5159,24 @@ static PetscErrorCode DMPlexComputeBdResidual_Internal(DM dm, Vec locX, Vec locX + PetscCall(DMPlexGetDepthLabel(dm, &depthLabel)); + PetscCall(DMGetDimension(dm, &dim)); + PetscCall(DMLabelGetStratumIS(depthLabel, dim - 1, &facetIS)); ++ /* Filter out ghost facets (SF leaves) so each facet contributes once */ ++ if (facetIS) { ++ PetscSF sf; ++ PetscInt nleaves; ++ const PetscInt *leaves; ++ ++ PetscCall(DMGetPointSF(dm, &sf)); ++ PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL)); ++ if (nleaves > 0 && leaves) { ++ IS leafIS, ownedFacetIS; ++ ++ PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nleaves, leaves, PETSC_USE_POINTER, &leafIS)); ++ PetscCall(ISDifference(facetIS, leafIS, &ownedFacetIS)); ++ PetscCall(ISDestroy(&leafIS)); ++ PetscCall(ISDestroy(&facetIS)); ++ facetIS = ownedFacetIS; ++ } ++ } + PetscCall(PetscDSGetNumBoundary(prob, &numBd)); + for (bd = 0; bd < numBd; ++bd) { + PetscWeakForm wf; +@@ -6117,13 +6181,31 @@ static PetscErrorCode DMPlexComputeBdJacobian_Internal(DM dm, Vec locX, Vec locX + PetscInt dim, numBd, bd; + DMLabel depthLabel; + DMField coordField = NULL; +- IS facetIS; ++ IS facetIS = NULL; + + PetscFunctionBegin; + PetscCall(DMGetDS(dm, &prob)); + PetscCall(DMPlexGetDepthLabel(dm, &depthLabel)); + PetscCall(DMGetDimension(dm, &dim)); + PetscCall(DMLabelGetStratumIS(depthLabel, dim - 1, &facetIS)); ++ /* Filter out ghost facets (SF leaves) so each facet contributes once */ ++ if (facetIS) { ++ PetscSF sf; ++ PetscInt nleaves; ++ const PetscInt *leaves; ++ ++ PetscCall(DMGetPointSF(dm, &sf)); ++ PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL)); ++ if (nleaves > 0 && leaves) { ++ IS leafIS, ownedFacetIS; ++ ++ PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nleaves, leaves, PETSC_USE_POINTER, &leafIS)); ++ PetscCall(ISDifference(facetIS, leafIS, &ownedFacetIS)); ++ PetscCall(ISDestroy(&leafIS)); ++ PetscCall(ISDestroy(&facetIS)); ++ facetIS = ownedFacetIS; ++ } ++ } + PetscCall(PetscDSGetNumBoundary(prob, &numBd)); + PetscCall(DMGetCoordinateField(dm, &coordField)); + for (bd = 0; bd < numBd; ++bd) { diff --git a/tests/parallel/test_0765_internal_boundary_integral_mpi.py b/tests/parallel/test_0765_internal_boundary_integral_mpi.py new file mode 100644 index 000000000..17ed49e24 --- /dev/null +++ b/tests/parallel/test_0765_internal_boundary_integral_mpi.py @@ -0,0 +1,61 @@ +""" +MPI regression test for internal-boundary boundary-integral ownership. + +This guards against rank-dependent over/under-assembly caused by +ghost/internal facet handling in PETSc boundary assembly paths. +""" + +import numpy as np +import pytest +import underworld3 as uw + + +pytestmark = [ + pytest.mark.level_2, + pytest.mark.tier_a, + pytest.mark.mpi(min_size=2), + pytest.mark.timeout(120), +] + + +def _annulus_mesh(): + return uw.meshing.AnnulusInternalBoundary( + radiusOuter=2.22, + radiusInternal=2.0, + radiusInner=1.22, + cellSize_Inner=1.0 / 32.0, + cellSize_Internal=(1.0 / 32.0) / 2.0, + cellSize_Outer=1.0 / 32.0, + ) + + +mesh_annulus = _annulus_mesh() +# PETSc integration path requires at least one variable on the mesh. +_dummy_var = uw.discretisation.MeshVariable("T_annulus_mpi_bd", mesh_annulus, 1, degree=1) + + +@pytest.mark.mpi(min_size=2) +def test_internal_boundary_circumference_parallel(): + """ + Internal boundary circumference should match 2*pi*R in parallel. + """ + value = float(uw.maths.BdIntegral(mesh=mesh_annulus, fn=1.0, boundary="Internal").evaluate()) + expected = 2.0 * np.pi * 2.0 + + rel_err = abs(value - expected) / expected + assert rel_err < 2.0e-2, f"Internal circumference rel_err={rel_err:.3e}, value={value}, expected={expected}" + + gathered = uw.mpi.comm.allgather(value) + assert max(gathered) - min(gathered) < 1.0e-12, f"Rank mismatch in integral values: {gathered}" + + +@pytest.mark.mpi(min_size=2) +def test_outer_boundary_circumference_parallel(): + """ + External boundary circumference remains correct in parallel. + """ + value = float(uw.maths.BdIntegral(mesh=mesh_annulus, fn=1.0, boundary="Upper").evaluate()) + expected = 2.0 * np.pi * 2.22 + + rel_err = abs(value - expected) / expected + assert rel_err < 2.0e-2, f"Outer circumference rel_err={rel_err:.3e}, value={value}, expected={expected}" From 861390b41043e3baed63711ec8e2e501fdc91d64 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Thu, 12 Mar 2026 12:37:04 +1100 Subject: [PATCH 2/2] docs(petsc): expand MR notes with reproducer, metrics, risk, and backport guidance --- ...DESCRIPTION_INTERNAL_BOUNDARY_OWNERSHIP.md | 130 ++++++++++++++++++ 1 file changed, 130 insertions(+) diff --git a/petsc-custom/patches/PETSC_MR_DESCRIPTION_INTERNAL_BOUNDARY_OWNERSHIP.md b/petsc-custom/patches/PETSC_MR_DESCRIPTION_INTERNAL_BOUNDARY_OWNERSHIP.md index fb529fc47..c82033429 100644 --- a/petsc-custom/patches/PETSC_MR_DESCRIPTION_INTERNAL_BOUNDARY_OWNERSHIP.md +++ b/petsc-custom/patches/PETSC_MR_DESCRIPTION_INTERNAL_BOUNDARY_OWNERSHIP.md @@ -10,6 +10,17 @@ File touched in PETSc: - `src/dm/impls/plex/plexfem.c` +## PETSc Version Details + +- Upstream baseline used for patch generation: + - PETSc repository branch: `release` + - Local baseline commit: `5a03b95372f` (before applying this patch) +- Runtime/tested configuration: + - PETSc: `3.24.5` (release build) + - petsc4py: `3.24.5` + - `PETSC_ARCH`: `petsc-4-uw-openmpi` + - `PETSC_DIR`: `.../underworld3/petsc-custom/petsc` + ## Problem Internal-boundary contributions in parallel were rank-dependent for natural BC / boundary integral workflows. @@ -59,6 +70,125 @@ Validated with Underworld3 annulus internal-boundary cases: - `tests/parallel/test_0765_internal_boundary_integral_mpi.py` - checks annulus internal-boundary circumference in parallel against analytic value. +## PETSc-Only Reproducer (For Upstream MR) + +Goal: provide a PETSc-native (no UW dependency) reproducer under +`src/dm/impls/plex/tests` that demonstrates rank-dependent internal-boundary +assembly before this patch and rank-invariant behavior after. + +### Proposed test shape + +- New PETSc test source: + - `src/dm/impls/plex/tests/ex_internalbd_ownership.c` (name illustrative) +- Setup: + - create a distributed 2D DMPlex box mesh + - define a facet `DMLabel` for an interior line (internal boundary) + - add a simple natural-boundary contribution over that internal label + - compute a deterministic scalar diagnostic (boundary integral or derived norm) +- Runs: + - `nsize: 1` and `nsize: 4` (or `8`) + - compare scalar outputs with strict tolerance + +### Why C test (not petsc4py) + +`petsc4py` does not currently expose the needed internal DMPlex boundary +assembly entry points directly enough for a clean reproducer of this path. +Upstream PETSc C tests are the correct home for long-term regression coverage. + +## Exact Before / After Numbers + +The following were measured on the same mesh/problem configuration +(`internal_bc_rank_reproducer.py`, cellsize `0.03125`): + +### Before fix (rank-dependent branch split) + +- `np=7`: `v_l2_int=2.3234597202763355e-03` +- `np=8`: `v_l2_int=6.2729732160772709e-03` + +### After fix (rank-invariant) + +- `np=7`: `v_l2_int=2.3234597202763251e-03` +- `np=8`: `v_l2_int=2.3234597202740370e-03` + +### Tolerance criteria + +- Rank-invariance target for this case: + - relative difference between MPI sizes `< 1e-10` for `v_l2_int` in this benchmark family +- Practical acceptance in PETSc test harness: + - choose tolerance robust to platform/compiler variation (e.g. `1e-9` to `1e-8`) + +## Changed-Function Map + +Patch updates are limited to: + +- `DMPlexComputeBdIntegral` +- `DMPlexComputeBdResidual_Internal` +- `DMPlexComputeBdJacobian_Internal` +- `DMPlexComputeBdResidualSingleByKey` + +No public API signatures were changed. + +## Rationale: `support[key.part]` vs `support[0]` + +Boundary weak-form callbacks can be registered with a nonzero `part` to target +a specific side of a facet. Using `support[0]` unconditionally can pull +closures from the wrong adjacent cell and produce side-inconsistent assembly. +Therefore `support[key.part]` must be used, with explicit support-size guards, +and points invalid for `part > 0` must be filtered. + +## Ownership Model Note (SF leaves) + +In parallel DMPlex, ghost facets appear as SF leaves on non-owning ranks. +If these are assembled in boundary paths, shared internal facets can +contribute multiple times and then be summed by MPI reductions / ADD_VALUES. +Filtering SF leaves ensures each physical facet contributes exactly once +(owner rank only), removing duplicate assembly. + +## PETSc Regression Test Plan (Upstream) + +Add a regression test under `src/dm/impls/plex/tests` with MPI coverage: + +- `nsize: 1` +- `nsize: 4` (or `8`) +- check scalar diagnostic equality within tolerance +- include in PETSc test harness (`TEST` block) so CI exercises the path + +Suggested command pattern once merged in PETSc tree: + +```bash +make -C $PETSC_DIR check TESTDIR=src/dm/impls/plex/tests +``` + +Or direct: + +```bash +mpirun -np 4 ./ex_internalbd_ownership [args] +``` + +## Version Scope + +- Baseline for patch generation: + - PETSc `release` branch, local baseline commit `5a03b95372f` +- Runtime validated in this workflow: + - PETSc `3.24.5` + - petsc4py `3.24.5` + - OpenMPI `5.0.10` +- Affected code sections are structurally stable across recent PETSc releases, + so backport applicability is expected to be broad with minimal adaptation. + +## Risk / Performance + +- No API changes. +- Correctness-only fix in boundary assembly ownership/part handling. +- Additional operations are lightweight (`ISDifference`, support checks) and + expected overhead is negligible relative to FE assembly costs. + +## Backport Recommendation + +- Recommended target: PETSc `release` branch. +- Also suitable for older maintained branches where the same boundary assembly + code is present (confirm with `git apply --check` and test run). + ## Benchmark Script And Case - Benchmark script (GitHub):