From ebc8dcecab08f8215150bbf109e06d7e58300eb6 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 21 Apr 2026 16:53:28 +1000 Subject: [PATCH 1/3] BdIntegral: use sandbox DM to avoid corrupting shared DM_Plex caches MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit DMPlexComputeBdIntegral lazily initialises height-trace FE objects on the coordinate DMField, which is cached on the shared DM_Plex struct. Since DMClone shares DM_Plex by refcount, this corrupts all solver DMs cloned from the same mesh — the next Jacobian assembly fails with MPI errors (ranks disagree on off-process matrix entries). Fix: BdIntegral.evaluate() now creates a disposable sandbox DM via DMClone + fresh createCoordinateSpace. The sandbox has its own coordinate DMField with empty height-trace caches, so the lazy init during the integral writes to the sandbox's caches, not the original mesh's. The sandbox is destroyed after each call. Reproducer: mpirun -np 2 with cellSize=0.05 on UnstructuredSimplexBox, solve → BdIntegral → solve. Fails without this fix on PETSc 3.25.0 (and intermittently on 3.24.x depending on mesh partitioning). Related: #96 Underworld development team with AI support from Claude Code (https://claude.com/claude-code) --- src/underworld3/cython/petsc_compat.h | 53 +++++++++++++++++++++++++ src/underworld3/cython/petsc_extras.pxi | 1 + src/underworld3/cython/petsc_maths.pyx | 38 ++++++++++-------- 3 files changed, 75 insertions(+), 17 deletions(-) diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h index c7204da77..b2f13a50a 100644 --- a/src/underworld3/cython/petsc_compat.h +++ b/src/underworld3/cython/petsc_compat.h @@ -258,3 +258,56 @@ 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, const char labelName[], + DM *sandbox) +{ + DM sdm; + PetscInt Nf, i; + PetscSection srcSec; + Vec auxVec; + + PetscFunctionBeginUser; + + // DMClone: shares DM_Plex topology, point SF, labels by reference. + // This is fine — the topology is read-only. The corruption is in + // the coordinate DMField's cached height-trace FEs, which we replace. + PetscCall(DMClone(src, &sdm)); + + // Rebuild coordinate space from scratch — creates a new coordinate + // DMField_DS with empty height-trace caches, independent of src's. + PetscCall(DMPlexCreateCoordinateSpace(sdm, 1, PETSC_FALSE, PETSC_TRUE)); + PetscCall(UW_DMForceCoordinateField(sdm)); + + // Copy fields + DS from src so the section layout matches + 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); +} diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi index 28c6af334..c544e6daf 100644 --- a/src/underworld3/cython/petsc_extras.pxi +++ b/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, const char[], PetscDM*) cdef extern from "petsc.h" nogil: PetscErrorCode PetscDSSetConstants(PetscDS, PetscInt, const PetscScalar[]) diff --git a/src/underworld3/cython/petsc_maths.pyx b/src/underworld3/cython/petsc_maths.pyx index a7d098614..2019a89d2 100644 --- a/src/underworld3/cython/petsc_maths.pyx +++ b/src/underworld3/cython/petsc_maths.pyx @@ -410,43 +410,47 @@ 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 that is topologically identical to mesh.dm + # but has an independent DM_Plex struct. DMPlexComputeBdIntegral + # mutates shared DM_Plex caches (height-trace FEs, closure indices) + # which corrupts solver DMs that were cloned from the same mesh. + # Running the integral on a disposable sandbox avoids this. + cdef PetscDM sandbox_dm = NULL + cdef bytes boundary_bytes = self.boundary.encode('utf-8') + CHKERRQ(UW_DMCreateBdIntegralSandbox( + (mesh.dm).dm, boundary_bytes, &sandbox_dm)) + + # Get the label from the sandbox DM + cdef PetscDMLabel c_dmlabel = NULL + cdef DMLabel dmlabel boundary_enum = mesh.boundaries[self.boundary] cdef PetscInt label_val = boundary_enum.value cdef PetscInt num_vals = 1 - 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 = result From e43a0f62faad9a5597b374c4bd88825c33e65a20 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 22 Apr 2026 10:10:32 +1000 Subject: [PATCH 2/3] Fix CI: version-guard DMPlexCreateCoordinateSpace signature PETSc 3.24.x: DMPlexCreateCoordinateSpace(dm, degree, project, snapFunc) PETSc 3.25.0: DMPlexCreateCoordinateSpace(dm, degree, localized, project) Underworld development team with AI support from Claude Code (https://claude.com/claude-code) --- src/underworld3/cython/petsc_compat.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h index b2f13a50a..fe140ff16 100644 --- a/src/underworld3/cython/petsc_compat.h +++ b/src/underworld3/cython/petsc_compat.h @@ -287,7 +287,13 @@ PetscErrorCode UW_DMCreateBdIntegralSandbox(DM src, const char labelName[], // Rebuild coordinate space from scratch — creates a new coordinate // DMField_DS with empty height-trace caches, independent of src's. + // 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)); // Copy fields + DS from src so the section layout matches From 7d1689ab19e26f252b0554022f107500fd85ce00 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 22 Apr 2026 15:45:18 +1000 Subject: [PATCH 3/3] Address Copilot review: remove unused params, fix comments - Remove unused labelName parameter from UW_DMCreateBdIntegralSandbox - Remove unused c_dmlabel/dmlabel Cython locals - Fix comment: sandbox shares DM_Plex topology via DMClone, the independent part is the coordinate DMField (not the whole DM_Plex) - Clarify DS comment (DMCreateDS recreates from fields, not a copy) Underworld development team with AI support from Claude Code (https://claude.com/claude-code) --- src/underworld3/cython/petsc_compat.h | 14 ++++++++------ src/underworld3/cython/petsc_extras.pxi | 2 +- src/underworld3/cython/petsc_maths.pyx | 19 +++++++++---------- 3 files changed, 18 insertions(+), 17 deletions(-) diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h index fe140ff16..ac4efa517 100644 --- a/src/underworld3/cython/petsc_compat.h +++ b/src/underworld3/cython/petsc_compat.h @@ -270,8 +270,7 @@ PetscErrorCode UW_DMPlexComputeBdIntegral(DM dm, Vec X, // 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, const char labelName[], - DM *sandbox) +PetscErrorCode UW_DMCreateBdIntegralSandbox(DM src, DM *sandbox) { DM sdm; PetscInt Nf, i; @@ -280,13 +279,15 @@ PetscErrorCode UW_DMCreateBdIntegralSandbox(DM src, const char labelName[], PetscFunctionBeginUser; - // DMClone: shares DM_Plex topology, point SF, labels by reference. - // This is fine — the topology is read-only. The corruption is in - // the coordinate DMField's cached height-trace FEs, which we replace. + // 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) @@ -296,7 +297,8 @@ PetscErrorCode UW_DMCreateBdIntegralSandbox(DM src, const char labelName[], #endif PetscCall(UW_DMForceCoordinateField(sdm)); - // Copy fields + DS from src so the section layout matches + // 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) { diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi index c544e6daf..91bcb1763 100644 --- a/src/underworld3/cython/petsc_extras.pxi +++ b/src/underworld3/cython/petsc_extras.pxi @@ -46,7 +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, const char[], PetscDM*) + PetscErrorCode UW_DMCreateBdIntegralSandbox(PetscDM, PetscDM*) cdef extern from "petsc.h" nogil: PetscErrorCode PetscDSSetConstants(PetscDS, PetscInt, const PetscScalar[]) diff --git a/src/underworld3/cython/petsc_maths.pyx b/src/underworld3/cython/petsc_maths.pyx index 2019a89d2..29e26a2ad 100644 --- a/src/underworld3/cython/petsc_maths.pyx +++ b/src/underworld3/cython/petsc_maths.pyx @@ -417,22 +417,21 @@ class BdIntegral: cdef Vec cgvec = a_global - # Create a sandbox DM that is topologically identical to mesh.dm - # but has an independent DM_Plex struct. DMPlexComputeBdIntegral - # mutates shared DM_Plex caches (height-trace FEs, closure indices) - # which corrupts solver DMs that were cloned from the same mesh. - # Running the integral on a disposable sandbox avoids this. + # 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 - cdef bytes boundary_bytes = self.boundary.encode('utf-8') CHKERRQ(UW_DMCreateBdIntegralSandbox( - (mesh.dm).dm, boundary_bytes, &sandbox_dm)) + (mesh.dm).dm, &sandbox_dm)) - # Get the label from the sandbox DM - cdef PetscDMLabel c_dmlabel = NULL - cdef DMLabel dmlabel + # 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') cdef PetscDMLabel sandbox_label = NULL CHKERRQ(DMGetLabel(sandbox_dm, boundary_bytes, &sandbox_label))