BdIntegral: sandbox DM to prevent shared DM_Plex cache corruption by lmoresi · Pull Request #125 · underworldcode/underworld3

lmoresi · 2026-04-21T06:53:52Z

Summary

DMPlexComputeBdIntegral lazily initialises height-trace FE caches on the coordinate DMField, which lives 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).
BdIntegral.evaluate() now creates a disposable sandbox DM (DMClone + fresh createCoordinateSpace) with its own coordinate field caches. The sandbox is destroyed after each call.
This is a workaround for a PETSc design issue (DMClone_Plex sharing mutable caches). Not optimised -- sandbox is created/destroyed per call, which is correct for mesh deformation and adaptation scenarios.

Reproducer

mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.05, ...)
# ... set up Stokes solver with BCs ...
stokes.solve()
bd = uw.maths.BdIntegral(mesh=mesh, fn=1.0, boundary="Top")
bd.evaluate()
stokes.solve()  # <-- MPI error without this fix

Fails with mpirun -np 2 on PETSc 3.25.0 (and intermittently on 3.24.x depending on mesh partitioning). Mesh size dependent -- only triggers when the partition creates cross-process boundary faces.

Test plan

mpirun -np 2 with cellSize 0.03--0.1 (all pass, previously 0.05/0.09 failed)
Multiple interleaved solve/BdIntegral cycles (5 iterations, 2 boundaries)
Aux-dependent viscosity (exp(T)) with mesh variable in integrand
@jcgraciosa to verify with spherical benchmark script

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)

Copilot

Pull request overview

This PR addresses a PETSc parallel correctness issue where calling DMPlexComputeBdIntegral can mutate shared DM_Plex-associated caches and corrupt other solver DMs cloned from the same mesh, causing subsequent Jacobian assembly/MPI failures. It introduces a disposable “sandbox” DM for BdIntegral.evaluate() to isolate those mutations.

Changes:

Update BdIntegral.evaluate() to run the boundary integral on a disposable sandbox DM and destroy it afterward.
Add a new PETSc C helper UW_DMCreateBdIntegralSandbox() to clone a DM and rebuild its coordinate field/caches for safe BdIntegral evaluation.
Expose the new helper to Cython via petsc_extras.pxi.

Reviewed changes

Copilot reviewed 3 out of 3 changed files in this pull request and generated 6 comments.

File	Description
`src/underworld3/cython/petsc_maths.pyx`	Uses a sandbox DM for BdIntegral evaluation and destroys it after the call.
`src/underworld3/cython/petsc_extras.pxi`	Adds Cython extern declaration for `UW_DMCreateBdIntegralSandbox`.
`src/underworld3/cython/petsc_compat.h`	Implements `UW_DMCreateBdIntegralSandbox()` to rebuild coordinate space/field and copy field setup/aux vec.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Copilot · 2026-04-21T07:00:11Z

+        cdef PetscDMLabel c_dmlabel = NULL
+        cdef DMLabel dmlabel


These Cython locals (c_dmlabel, dmlabel) are now unused after switching to sandbox_label. Please remove them to avoid confusion and keep the evaluate() logic minimal.

Suggested change

cdef PetscDMLabel c_dmlabel = NULL

cdef DMLabel dmlabel

Copilot · 2026-04-21T07:00:11Z

+    // 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));


UW_DMCreateBdIntegralSandbox() hard-codes DMPlexCreateCoordinateSpace(sdm, 1, ...). Underworld meshes can be created with degree>1 (quadratic/triquadratic), and using degree=1 for the sandbox coordinate space can change the geometry mapping and yield incorrect boundary integrals on curved/high-order meshes. Please derive the coordinate degree from src (or pass it in from Python, e.g. mesh.degree) and use that instead of 1.

Copilot · 2026-04-21T07:00:12Z

+PetscErrorCode UW_DMCreateBdIntegralSandbox(DM src, const char labelName[],
+                                            DM *sandbox)
+{


The labelName parameter is currently unused in UW_DMCreateBdIntegralSandbox(), but is threaded through from Cython. If it’s not needed, remove it from the signature to avoid confusion (or explicitly mark it unused to silence compiler warnings).

Copilot · 2026-04-21T07:00:12Z

+    // 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));


This comment says “Copy fields + DS from src…”, but the implementation copies fields and then calls DMCreateDS(sdm) (which recreates a DS from options rather than copying src’s DS). Consider correcting the comment to reflect what the code actually does (copy fields/section, then recreate DS) to avoid future maintenance mistakes.

Copilot · 2026-04-21T07:00:12Z

+        # 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.


The comment says the sandbox DM has an “independent DM_Plex struct”, but UW_DMCreateBdIntegralSandbox() currently starts with DMClone(), which (as noted in the C comment) shares the underlying DM_Plex by refcount. This is misleading; consider rewording to state that the sandbox rebuilds the coordinate field/coordinate DMField caches (not the DM_Plex) to avoid cache contamination.

Suggested change

# 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 disposable sandbox DM that is topologically identical to

# mesh.dm. The sandbox is used to rebuild coordinate field /

# coordinate DMField-related caches needed by boundary integration

# without contaminating cached state on the solver DM. This does

# not imply an independent underlying DM_Plex struct.

Copilot · 2026-04-21T07:00:12Z

+        cdef PetscDM sandbox_dm = NULL
+        cdef bytes boundary_bytes = self.boundary.encode('utf-8')
+        CHKERRQ(UW_DMCreateBdIntegralSandbox(
+            (<DM>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)



sandbox_dm and the borrowed global vector (a_global) are only cleaned up on the success path. If UW_DMCreateBdIntegralSandbox(), DMGetLabel(), or UW_DMPlexComputeBdIntegral() raises via CHKERRQ, DMDestroy(&sandbox_dm) and mesh.dm.restoreGlobalVec(a_global) will be skipped, leaking PETSc objects and potentially leaving the DM’s global Vec checked-out. Please wrap the sandbox/global-vec lifetime in a try/finally (or equivalent) to guarantee cleanup on errors.

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)

- 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)

Copilot AI review requested due to automatic review settings April 21, 2026 06:53

lmoresi mentioned this pull request Apr 21, 2026

Fix BdIntegral parallel hang: coordinate DM label contamination #115

Merged

5 tasks

Copilot started reviewing on behalf of lmoresi April 21, 2026 06:54 View session

Copilot AI reviewed Apr 21, 2026

View reviewed changes

lmoresi added 2 commits April 22, 2026 10:10

lmoresi merged commit b838e0f into development Apr 22, 2026
1 check passed

lmoresi deleted the bugfix/bdintegral-sandbox-dm branch June 13, 2026 00:47

Navigation Menu

Uh oh!

BdIntegral: sandbox DM to prevent shared DM_Plex cache corruption - #125

BdIntegral: sandbox DM to prevent shared DM_Plex cache corruption#125
lmoresi merged 3 commits into
developmentfrom
bugfix/bdintegral-sandbox-dm

lmoresi commented Apr 21, 2026 •

edited by jcgraciosa

Loading

Uh oh!

Copilot AI left a comment

Uh oh!

Copilot AI Apr 21, 2026

Uh oh!

Copilot AI Apr 21, 2026

Uh oh!

Copilot AI Apr 21, 2026

Uh oh!

Copilot AI Apr 21, 2026

Uh oh!

Copilot AI Apr 21, 2026

Uh oh!

Copilot AI Apr 21, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

Conversation

lmoresi commented Apr 21, 2026 • edited by jcgraciosa Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

Reproducer

Test plan

Related

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Copilot AI Apr 21, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 21, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 21, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 21, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 21, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 21, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

lmoresi commented Apr 21, 2026 •

edited by jcgraciosa

Loading