Parallel mesh-adaptation fixes: MUMPS heap-corruption, patch volumes, collective remesh, field remap

[lmoresi]

Summary

Makes parallel mesh adaptation (np>=3) correct for the OT/MA mesh-motion movers. Four parallel-correctness fixes, uncovered in one investigation (the last three surfaced by running the adaptive convection harness at np=4). Serial behaviour is unchanged/bit-identical throughout; tier-A is green (177 passed / 0 failed) at every stage.

1. MUMPS-in-parallel — the heap corruption (root cause)

MUMPS (parallel LU) re-solved in a loop at np>=3 corrupts the heap (probabilistic SEGV/SIGBUS/MPI-deadlock). NOT the singular pure-Neumann solve, n_components>dim, or DMClone-shared caches (all disproved). The default GMRES+GAMG never touches MUMPS and is clean. The movers crashed because every sub-solve was wired through _use_direct_solver (lagged MUMPS LU), and _use_iterative_solver reintroduced MUMPS as the GAMG coarse solver.

• _use_direct_solver: iterative fallback under MPI; serial keeps lagged MUMPS.
• _use_iterative_solver: GAMG coarse uses redundant+svd (not MUMPS) under MPI.
• _wire: forwards elliptic.

2. Parallel-correct _patch_volumes

The per-vertex dual area (= lumped P1 mass diagonal) is assembled via the FE mass matrix (DM.createMassMatrix + M.1) in parallel so PETSc does the cross-rank localToGlobal(ADD); the hand-rolled local sum under-counted partition-boundary vertices. TODO to switch to DMCreateMassMatrixLumped once petsc4py binds it.

3. Collective remesh decision

mesh_metric_mismatch computed misalignment with np.corrcoef on rank-local cells → ranks disagreed across skip_threshold and the collective mover deadlocked. Now: alignment from globally-reduced moment sums (serial bit-identical) + the skip/adapt decision OR-reduced — if any rank remeshes, all do.

4. global_evaluate for the post-adapt field remap

The field remap (evaluate old field at new DOF coords) used the local uw.function.evaluate; DOFs that move across a rank-partition seam land off-rank, where local evaluate returns stale/garbage → a growing field anomaly at the seams. Fixed to uw.function.global_evaluate (serial-identical) in core _ot_adapt.py (reference-mesh and adapted-position remaps — affects all mesh.OT_adapt users) and the harness.

Verification

• np=5: OT/MA movers + mesh.OT_adapt(field) crash/hang-free; solvers converge and match serial to ~1e-10.
• np=4 adaptive convection (8 steps): clean run, unanimous adapt/skip across ranks, valid non-tangled adapted mesh; the field anomaly at the partition seams is removed (smoother T, correct boundary-layer values).
• tier-A: 177 passed, 0 failed (throughout).

Follow-ups (not in this PR)

• A general guard/warning for pc_factor_mat_solver_type=mumps / pc_type=lu re-solved at np>=3 (UW3-wide; this PR covers the movers).
• Boundary slivers at the partition seams (anisotropic mover slip handling) — lower priority, value-correct.

Design doc: docs/developer/design/parallel-repeated-solve-corruption.md.

Underworld development team with AI support from Claude Code

[Copilot]

Pull request overview

This PR addresses a reproducible parallel crash/hang in the adaptive mesh movers by avoiding repeated parallel MUMPS usage (direct LU) and by fixing a parallel correctness issue in per-vertex patch-volume computation used by equidistribution.

Changes:

• Disable the MUMPS-based “direct” solver path under MPI by falling back to the iterative configuration; also remove MUMPS from the GAMG coarse solve in parallel (use redundant+svd).
• Make _patch_volumes parallel-correct by computing the lumped P1 mass diagonal via FE assembly (M·1) when running under MPI.
• Add a design note documenting the investigation and fix.

Reviewed changes

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

File	Description
src/underworld3/meshing/smoothing.py	Gates MUMPS usage under MPI, adjusts GAMG coarse options, and fixes parallel patch-volume assembly for equidistribution movers.
docs/developer/design/parallel-repeated-solve-corruption.md	Documents the root-cause investigation and the mover-local mitigation strategy.

Comments suppressed due to low confidence (2)

src/underworld3/meshing/smoothing.py:1414

• In parallel (uw.mpi.size > 1), normalising patch with the local np.mean(patch) makes the scaling partition-dependent (each rank will use a different normalisation constant). Since patch is then used to define vol_field and to form b = rho_t * patch, this can introduce MPI-layout-dependent results. Prefer a partition-invariant normalisation (global mean over the MeshVariable’s global Vec) when running under MPI, while keeping the current numpy mean in serial for bit-identical behaviour.

            patch = _patch_volumes(tris, old_coords, n_verts, vol_field)
            patch /= float(np.mean(patch))
        else:
            patch = np.ones(n_verts, dtype=np.double)
        _va = vol_field.array
        _va[...] = patch.reshape(_va.shape)

src/underworld3/meshing/smoothing.py:1761

• patch_mean is currently computed as an average of per-rank means (allreduce(mean)/size), which is not the mean “over the domain” unless every rank owns the same number of vertices (and it can also be sensitive to ghost duplication depending on the local array layout). Since patch_mean is used to normalise both vol_field and the source term, this can make results depend on the partitioning. Prefer using the MeshVariable’s global mean under MPI (via vol_field.mean()), while keeping the numpy mean in serial for bit-identical behaviour.

            patch = _patch_volumes(tris, old_coords, n_verts, vol_field)
        # Normalise so the mean over the domain is the cell mean.
        patch_mean = float(np.mean(patch))
        if uw.mpi.size > 1:
            patch_mean = uw.mpi.comm.allreduce(patch_mean) / uw.mpi.size

---

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