BdIntegral: sandbox DM to prevent shared DM_Plex cache corruption

src/underworld3/cython/petsc_compat.h

@@ -258,3 +258,64 @@ PetscErrorCode UW_DMPlexComputeBdIntegral(DM dm, Vec X,
 
     PetscFunctionReturn(PETSC_SUCCESS);
 }
+
+// Create a sandbox DM for BdIntegral that won't contaminate the
+// original DM's shared DM_Plex caches.
+//
+// DMClone shares DM_Plex by refcount, and DMPlexComputeBdIntegral
+// lazily initialises height-trace FE caches on the shared struct.
+// This corrupts any solver DM cloned from the same mesh.
+//
+// Strategy: DMClone (cheap, shares topology) then rebuild the
+// coordinate field from scratch.  The fresh coordinate field has its
+// own DMField_DS with empty height-trace caches, so the lazy init
+// during BdIntegral writes to the sandbox's caches, not the original's.
+PetscErrorCode UW_DMCreateBdIntegralSandbox(DM src, DM *sandbox)
+{
+    DM            sdm;
+    PetscInt      Nf, i;
+    PetscSection  srcSec;
+    Vec           auxVec;
+
+    PetscFunctionBeginUser;
+
+    // DMClone shares DM_Plex topology, point SF, and labels by reference.
+    // The topology itself is read-only — the problem is the coordinate
+    // DMField's lazily-cached height-trace FEs.  We replace the coordinate
+    // field below so those caches are independent.
+    PetscCall(DMClone(src, &sdm));
+
+    // Rebuild coordinate space from scratch — creates a new coordinate
+    // DMField_DS with empty height-trace caches, independent of src's.
+    // Degree 1 (P1/Q1) matches UW3's mesh coordinate convention.
+    // Signature changed in PETSc 3.25: (dm, degree, localized, project)
+    // vs 3.24: (dm, degree, project, snapFunc)
+#if PETSC_VERSION_GE(3, 25, 0)
+    PetscCall(DMPlexCreateCoordinateSpace(sdm, 1, PETSC_FALSE, PETSC_TRUE));
+#else
+    PetscCall(DMPlexCreateCoordinateSpace(sdm, 1, PETSC_FALSE, NULL));
+#endif
+    PetscCall(UW_DMForceCoordinateField(sdm));
+
+    // Replicate the field layout and section from src so that the
+    // solution vector is compatible with the sandbox's DS.
+    PetscCall(DMGetLocalSection(src, &srcSec));
+    PetscCall(PetscSectionGetNumFields(srcSec, &Nf));
+    for (i = 0; i < Nf; ++i) {
+        PetscObject obj;
+        DMLabel     label;
+        PetscCall(DMGetField(src, i, &label, &obj));
+        PetscCall(DMSetField(sdm, i, label, obj));
+    }
+    PetscCall(DMCreateDS(sdm));
+    PetscCall(DMSetLocalSection(sdm, srcSec));
+
+    // Copy auxiliary data (viscosity, etc.) — just a pointer, no mutation
+    PetscCall(DMGetAuxiliaryVec(src, NULL, 0, 0, &auxVec));
+    if (auxVec) {
+        PetscCall(DMSetAuxiliaryVec(sdm, NULL, 0, 0, auxVec));
+    }
+
+    *sandbox = sdm;
+    PetscFunctionReturn(PETSC_SUCCESS);
+}

src/underworld3/cython/petsc_extras.pxi

@@ -46,6 +46,7 @@ cdef extern from "petsc_compat.h":
     PetscErrorCode UW_DMPlexSetSNESLocalFEM( PetscDM, PetscBool, void *)
     PetscErrorCode UW_DMForceCoordinateField(PetscDM)
     PetscErrorCode UW_DMPlexComputeBdIntegral( PetscDM, PetscVec, PetscDMLabel, PetscInt, const PetscInt*, void*, PetscScalar*, void*)
+    PetscErrorCode UW_DMCreateBdIntegralSandbox(PetscDM, PetscDM*)
 
 cdef extern from "petsc.h" nogil:
     PetscErrorCode PetscDSSetConstants(PetscDS, PetscInt, const PetscScalar[])

src/underworld3/cython/petsc_maths.pyx

@@ -410,43 +410,46 @@ class BdIntegral:
             self.mesh.vars.values(), verbose=verbose)
         cdef PtrContainer ext = _getext_result.ptrobj
 
-        # Prepare the solution vector
+        # Prepare the solution vector on the original mesh DM
         self.mesh.update_lvec()
         a_global = mesh.dm.getGlobalVec()
         mesh.dm.localToGlobal(self.mesh.lvec, a_global)
 
         cdef Vec cgvec = a_global
-        cdef DM dm_c = mesh.dm
 
-        # Get the DMLabel and label value for the named boundary
+        # Create a sandbox DM via DMClone + fresh coordinate space.
+        # DMPlexComputeBdIntegral lazily initialises height-trace FE caches
+        # on the coordinate DMField, which is shared via DMClone's refcount.
+        # This corrupts solver DMs cloned from the same mesh.  The sandbox
+        # gets its own coordinate field with empty caches, so the lazy init
+        # writes there instead of on the original mesh DM.
+        cdef PetscDM sandbox_dm = NULL
+        CHKERRQ(UW_DMCreateBdIntegralSandbox(
+            (<DM>mesh.dm).dm, &sandbox_dm))
+
+        # Get the boundary label from the sandbox (shared with mesh.dm via DMClone)
         boundary_enum = mesh.boundaries[self.boundary]
         cdef PetscInt label_val = boundary_enum.value
         cdef PetscInt num_vals = 1
+        cdef bytes boundary_bytes = self.boundary.encode('utf-8')
 
-        c_label = mesh.dm.getLabel(self.boundary)
+        cdef PetscDMLabel sandbox_label = NULL
+        CHKERRQ(DMGetLabel(sandbox_dm, boundary_bytes, &sandbox_label))
 
         # Output value
         cdef PetscScalar result = 0.0
 
-        # Label can be None if this boundary has no facets on this process
-        # (normal in parallel) or if the DM doesn't carry the label at all.
-        # The C wrapper handles NULL labels gracefully (contributes 0 to the
-        # MPI reduction so all ranks still participate).
-        cdef PetscDMLabel c_dmlabel = NULL
-        cdef DMLabel dmlabel
-        if c_label:
-            dmlabel = c_label
-            c_dmlabel = dmlabel.dmlabel
-
-        # Call the boundary integral
+        # Call the boundary integral on the SANDBOX DM
         ierr = UW_DMPlexComputeBdIntegral(
-            dm_c.dm, cgvec.vec,
-            c_dmlabel, num_vals, &label_val,
+            sandbox_dm, cgvec.vec,
+            sandbox_label, num_vals, &label_val,
             ext.fns_bd_residual[0],
             &result, NULL
         )
         CHKERRQ(ierr)
 
+        # Clean up sandbox
+        CHKERRQ(DMDestroy(&sandbox_dm))
         mesh.dm.restoreGlobalVec(a_global)
 
         cdef double vald = <double> result