Revert "Fix ghost facet double-counting in solver natural BCs under MPI"

lmoresi · lmoresi · commit c91fedd7c972 · 2026-03-04T11:41:19.000+11:00
This reverts commit d998732.
diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h
@@ -139,106 +139,6 @@ PetscErrorCode UW_DMPlexSetSNESLocalFEM(DM dm, PetscBool flag, void *ctx)
 #endif
 }
 
-// Create a ghost-filtered copy of a DMLabel on a DM.
-//
-// PETSc's boundary integration and residual routines iterate over ALL stratum
-// points in a label, including ghost facets.  For internal boundaries whose
-// facets sit between cells on different ranks, this leads to double-counting.
-//
-// This helper creates a new label (dstName) that contains only the owned
-// (non-ghost) points from srcLabel.  Ghost points are identified via the
-// point SF: every leaf in the SF is a ghost point.
-//
-// Usage:
-//   UW_CreateOwnedLabel(dm, srcLabel, "UW_Boundaries_no_ghost");
-//   // The new label is now attached to dm under the given name.
-PetscErrorCode UW_CreateOwnedLabel(DM dm, DMLabel srcLabel, const char dstName[])
-{
-    PetscInt  pStart, pEnd;
-
-    PetscFunctionBeginUser;
-
-    // --- Build ghost-point bitset from the point SF ---
-    PetscSF       sf;
-    PetscInt      nleaves;
-    const PetscInt *ilocal;
-    PetscCall(DMGetPointSF(dm, &sf));
-    PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &ilocal, NULL));
-    PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
-
-    PetscBT ghostBT;
-    PetscCall(PetscBTCreate(pEnd - pStart, &ghostBT));
-    if (ilocal) {
-        for (PetscInt i = 0; i < nleaves; i++) {
-            PetscCall(PetscBTSet(ghostBT, ilocal[i] - pStart));
-        }
-    }
-
-    // --- Remove any pre-existing label with this name, then create fresh ---
-    DMLabel existing;
-    PetscCall(DMGetLabel(dm, dstName, &existing));
-    if (existing) {
-        PetscCall(DMRemoveLabel(dm, dstName, NULL));
-    }
-    PetscCall(DMCreateLabel(dm, dstName));
-
-    DMLabel dstLabel;
-    PetscCall(DMGetLabel(dm, dstName, &dstLabel));
-
-    // --- Copy each stratum, keeping only owned points ---
-    PetscInt numValues;
-    IS       valueIS;
-    PetscCall(DMLabelGetNumValues(srcLabel, &numValues));
-    PetscCall(DMLabelGetValueIS(srcLabel, &valueIS));
-
-    const PetscInt *values;
-    PetscCall(ISGetIndices(valueIS, &values));
-
-    for (PetscInt v = 0; v < numValues; v++) {
-        PetscInt val = values[v];
-
-        IS origIS;
-        PetscCall(DMLabelGetStratumIS(srcLabel, val, &origIS));
-        if (!origIS) continue;
-
-        PetscInt        n;
-        const PetscInt *indices;
-        PetscCall(ISGetLocalSize(origIS, &n));
-        PetscCall(ISGetIndices(origIS, &indices));
-
-        // Count owned points
-        PetscInt nOwned = 0;
-        for (PetscInt i = 0; i < n; i++) {
-            if (!PetscBTLookup(ghostBT, indices[i] - pStart)) nOwned++;
-        }
-
-        // Build owned-only index array
-        PetscInt *owned;
-        PetscCall(PetscMalloc1(nOwned, &owned));
-        PetscInt j = 0;
-        for (PetscInt i = 0; i < n; i++) {
-            if (!PetscBTLookup(ghostBT, indices[i] - pStart)) {
-                owned[j++] = indices[i];
-            }
-        }
-
-        PetscCall(ISRestoreIndices(origIS, &indices));
-        PetscCall(ISDestroy(&origIS));
-
-        IS ownedIS;
-        PetscCall(ISCreateGeneral(PETSC_COMM_SELF,
-                                  nOwned, owned, PETSC_OWN_POINTER, &ownedIS));
-        PetscCall(DMLabelSetStratumIS(dstLabel, val, ownedIS));
-        PetscCall(ISDestroy(&ownedIS));
-    }
-
-    PetscCall(ISRestoreIndices(valueIS, &values));
-    PetscCall(ISDestroy(&valueIS));
-    PetscCall(PetscBTDestroy(&ghostBT));
-
-    PetscFunctionReturn(PETSC_SUCCESS);
-}
-
 // Simplified wrapper for DMPlexComputeBdIntegral.
 // Takes a single boundary pointwise function (for field 0) instead of an Nf-element array.
 //
diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi
@@ -43,7 +43,6 @@ cdef extern from "petsc_compat.h":
     PetscErrorCode UW_PetscDSViewWF(PetscDS)     
     PetscErrorCode UW_PetscDSViewBdWF(PetscDS, PetscInt)
     PetscErrorCode UW_DMPlexSetSNESLocalFEM( PetscDM, PetscBool, void *)
-    PetscErrorCode UW_CreateOwnedLabel( PetscDM, PetscDMLabel, const char[])
     PetscErrorCode UW_DMPlexComputeBdIntegral( PetscDM, PetscVec, PetscDMLabel, PetscInt, const PetscInt*, void*, PetscScalar*, void*)
 
 cdef extern from "petsc.h" nogil:
diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -1317,7 +1317,7 @@ class SNES_Scalar(SolverBaseClass):
             bc = PetscDSAddBoundary_UW(cdm.dm,
                                 bc_type,
                                 str(boundary+f"{bc.components}").encode('utf8'),
-                                "UW_Boundaries_no_ghost".encode('utf8'),   # ghost-filtered for natural BCs
+                                "UW_Boundaries".encode('utf8'),   # consolidated boundary label
                                 bc.f_id,  # field ID in the DM
                                 num_constrained_components,
                                 <const PetscInt *> &comps_view[0],
@@ -1535,7 +1535,6 @@ class SNES_Scalar(SolverBaseClass):
                 )
 
         ## Now add the boundary residual / jacobian terms
-        ## Use ghost-filtered label to prevent double-counting on internal boundaries
 
         cdef DMLabel c_label
 
@@ -1545,7 +1544,7 @@ class SNES_Scalar(SolverBaseClass):
             boundary_id = bc.PETScID
 
             value = self.mesh.boundaries[bc.boundary].value
-            bc_label = self.dm.getLabel("UW_Boundaries_no_ghost")
+            bc_label = self.dm.getLabel("UW_Boundaries")
 
             label_val = value
 
@@ -2038,7 +2037,7 @@ class SNES_Vector(SolverBaseClass):
             bc = PetscDSAddBoundary_UW(cdm.dm,
                                 bc_type,
                                 str(boundary+f"{bc.components}").encode('utf8'),
-                                "UW_Boundaries_no_ghost".encode('utf8'),   # ghost-filtered for natural BCs
+                                "UW_Boundaries".encode('utf8'),   # was: str(boundary)
                                 bc.f_id,  # field ID in the DM
                                 num_constrained_components,
                                 <const PetscInt *> &comps_view[0],
@@ -2073,7 +2072,7 @@ class SNES_Vector(SolverBaseClass):
             bc = PetscDSAddBoundary_UW(cdm.dm,
                                 bc_type,
                                 str(boundary+f"{bc.components}").encode('utf8'),
-                                "UW_Boundaries".encode('utf8'),
+                                "UW_Boundaries".encode('utf8'),   # was: str(boundary)
                                 bc.f_id,  # field ID in the DM
                                 num_constrained_components,
                                 <const PetscInt *> &comps_view[0],
@@ -2278,7 +2277,6 @@ class SNES_Vector(SolverBaseClass):
                 )
 
         ## SNES VECTOR ADD Boundary terms
-        ## Use ghost-filtered label to prevent double-counting on internal boundaries
 
         cdef DMLabel c_label
 
@@ -2288,7 +2286,8 @@ class SNES_Vector(SolverBaseClass):
             boundary_id = bc.PETScID
 
             value = self.mesh.boundaries[bc.boundary].value
-            bc_label = self.dm.getLabel("UW_Boundaries_no_ghost")
+            bc_label = self.dm.getLabel("UW_Boundaries")
+            #bc_label = self.dm.getLabel(boundary)
 
             label_val = value
 
@@ -3455,7 +3454,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             bc = PetscDSAddBoundary_UW(cdm.dm,
                                 bc_type,
                                 str(boundary+f"{bc.components}").encode('utf8'),
-                                str("UW_Boundaries_no_ghost").encode('utf8'),  # ghost-filtered for natural BCs
+                                str("UW_Boundaries").encode('utf8'),
                                 bc.f_id,  # field ID in the DM
                                 num_constrained_components,
                                 <const PetscInt *> &comps_view[0],
@@ -3553,7 +3552,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             boundary_id = bc.PETScID
 
             value = self.mesh.boundaries[bc.boundary].value
-            bc_label = self.dm.getLabel("UW_Boundaries_no_ghost")
+            bc_label = self.dm.getLabel("UW_Boundaries")
+            # bc_label = self.dm.getLabel(boundary)
 
             label_val = value
 
diff --git a/src/underworld3/cython/petsc_maths.pyx b/src/underworld3/cython/petsc_maths.pyx
@@ -15,32 +15,6 @@ cdef extern from "petsc.h" nogil:
     PetscErrorCode DMPlexComputeCellwiseIntegralFEM( PetscDM, PetscVec, PetscVec, void* )
 
 
-def create_owned_label(dm, src_label_name, dst_label_name):
-    """Create a ghost-filtered copy of a DMLabel on a DM.
-
-    Ghost facets (identified via the point SF) are excluded from the
-    destination label.  This prevents double-counting when PETSc
-    integrates over boundary label strata in parallel.
-
-    Parameters
-    ----------
-    dm : PETSc.DM
-        The distributed mesh DM.
-    src_label_name : str
-        Name of the source label (e.g. ``"UW_Boundaries"``).
-    dst_label_name : str
-        Name for the new ghost-filtered label (e.g. ``"UW_Boundaries_no_ghost"``).
-    """
-
-    cdef DM c_dm = dm
-    src_label = dm.getLabel(src_label_name)
-    if src_label is None:
-        return
-
-    cdef DMLabel c_src = src_label
-    UW_CreateOwnedLabel(c_dm.dm, c_src.dmlabel, dst_label_name.encode('utf8'))
-
-
 class Integral:
     """
     The `Integral` class constructs the volume integral
diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py
@@ -438,14 +438,6 @@ class replacement_boundaries(Enum):
 
             uw.mpi.barrier()
 
-            # Ghost-filtered copy for natural BC residual integration.
-            # PETSc boundary integration iterates ALL label stratum points
-            # including ghost facets, causing double-counting on internal
-            # boundaries.  The _no_ghost variant is used by solvers for
-            # natural BC registration.
-            from underworld3.cython.petsc_maths import create_owned_label
-            create_owned_label(self.dm, "UW_Boundaries", "UW_Boundaries_no_ghost")
-
         ## ---
         ## Note - coarsening callback is tricky because the coarse meshes do not have the labels
         ##
diff --git a/tests/test_0503_solver_natural_bc.py b/tests/test_0503_solver_natural_bc.py
