diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h index c7204da77..ac4efa517 100644 --- a/src/underworld3/cython/petsc_compat.h +++ b/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); +} diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi index 28c6af334..91bcb1763 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, 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..29e26a2ad 100644 --- a/src/underworld3/cython/petsc_maths.pyx +++ b/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( + (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 = result