From 7dbfd2f24c1f18a32b7d02ed0a07c18815feffd6 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 20 Jun 2026 12:56:42 +1000 Subject: [PATCH 1/2] fix(meshing): SPD-floor the mmpde metric to stop silent NaN-bail on deformed meshes _winslow_mmpde's energy/gradient use fractional powers (sqrt(detM), detM**((1-p)/2), and S**q with S = tr(J M^-1 J^T)) that are defined only for an SPD metric. The metric is a guide field FE-evaluated at a FIXED reference cloud; once the interior deforms, a reference point can fall outside the current mesh and the P1 density is then evaluated by FE extrapolation, going negative. A scalar density rho<0 gives M = rho*I with det = rho^2 > 0 (so a detM>0 check passes) but M negative-definite, so S<0 and S**q = NaN. The energy is then non-finite and the mover bails with zero displacement -> adaptation silently stops on roughly half the adapts of an adaptive convection loop (mesh under-resolves; the symptom was previously misread as adapt/solver coupling or "passive fault damping"). Project every evaluated metric tensor onto SPD with a small relative eigenvalue floor at the _eval_M chokepoint. A genuinely SPD metric is returned unchanged (short-circuit no-op); extrapolation garbage becomes a benign "coarsen here" instead of a NaN. No behaviour change on valid metrics (bit-identical), no API change. Tests: test_0762/0850/0750/0830 green (33 passed, 5 skipped MMG, 4 pre-existing xfail). Validated on the stagnant-lid adaptive convection driver: all forced adapts move, T stays in [0,1], vrms bounded over 52 steps (no-fault) and 32 steps (passive fault). Underworld development team with AI support from Claude Code --- src/underworld3/meshing/smoothing.py | 34 ++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/src/underworld3/meshing/smoothing.py b/src/underworld3/meshing/smoothing.py index b364dfe8..1061a653 100644 --- a/src/underworld3/meshing/smoothing.py +++ b/src/underworld3/meshing/smoothing.py @@ -3164,6 +3164,40 @@ def _eval_M(pts): else: _eval_M = _eval_M_analytic + # --- SPD sanitiser on the evaluated metric ------------------------- + # The MMPDE functional G uses fractional powers that are defined ONLY + # for an SPD metric: sqrt(detM), detM**((1-p)/2), and S**q with + # S = tr(J M⁻¹ Jᵀ). The metric is a guide field FE-evaluated at the + # FIXED reference cloud; once the interior has deformed, a reference + # point can fall OUTSIDE the current mesh and the P1 metric field is + # then evaluated by FE EXTRAPOLATION (out-of-cell basis functions go + # negative), yielding a non-SPD tensor — e.g. a scalar density ρ·I + # with ρ<0. Its determinant ρ² stays positive (so a detM>0 test + # passes) but M is negative-definite, so S<0 and S**q = NaN → the + # energy is non-finite and the mover bails with zero displacement + # (no adaptation). Project every evaluated tensor onto SPD with a + # small RELATIVE eigenvalue floor: a genuine SPD metric is returned + # unchanged (no-op), while extrapolation garbage becomes a benign + # "coarsen here" (tiny positive eigenvalues) instead of a NaN. + _eval_M_raw = _eval_M + + def _spd_sanitise(M): + Ms = 0.5 * (M + np.swapaxes(M, -1, -2)) + w, Vc = np.linalg.eigh(Ms) + if w.size: + wmax = float(np.nanmax(w)) + else: + wmax = 1.0 + floor = max(wmax, 1.0) * 1.0e-8 + if np.all(np.isfinite(w)) and float(w.min()) >= floor: + return Ms # already SPD → no-op + w = np.clip(np.nan_to_num(w, nan=floor, posinf=wmax, neginf=floor), + floor, None) + return np.einsum('nij,nj,nkj->nik', Vc, w, Vc) + + def _eval_M(pts): + return _spd_sanitise(_eval_M_raw(pts)) + # Mesh-owned boundary slip is applied per outer iter via mesh.boundary_slip # (below). Pre-touch Gamma_P1 here so the projected-normal MeshVariable # exists before any DM snapshot (footgun-safe; redundant with the central From 732557ea1ff0a649c11e141736dcaf860dd2ab2e Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 20 Jun 2026 15:35:12 +1000 Subject: [PATCH 2/2] fix(meshing): per-tensor SPD masking + empty-rank guard (review) Address Copilot review on _spd_sanitise: - Guard the empty batch (a parallel rank that owns no cells): return early before the eigenvalue reductions, which previously raised on a zero-size array. - Make the SPD test per-tensor: only the cells whose OWN eigenvalues are non-finite or below the floor are projected; every already-SPD tensor is returned bit-identical to the symmetrised input. One bad point can no longer perturb otherwise-valid cells (matches the "no-op for SPD" contract). Verified: empty input is safe, a single garbage cell leaves the rest bit-identical, all output tensors SPD. Underworld development team with AI support from Claude Code --- src/underworld3/meshing/smoothing.py | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/src/underworld3/meshing/smoothing.py b/src/underworld3/meshing/smoothing.py index 1061a653..5272cfbb 100644 --- a/src/underworld3/meshing/smoothing.py +++ b/src/underworld3/meshing/smoothing.py @@ -3182,18 +3182,26 @@ def _eval_M(pts): _eval_M_raw = _eval_M def _spd_sanitise(M): + # Symmetrise: the metric is symmetric by construction, so for a valid + # tensor this is an exact no-op (M_ij == M_ji bit-for-bit). Ms = 0.5 * (M + np.swapaxes(M, -1, -2)) + if Ms.shape[0] == 0: + return Ms # rank owns no cells w, Vc = np.linalg.eigh(Ms) - if w.size: - wmax = float(np.nanmax(w)) - else: - wmax = 1.0 + wmax = float(np.nanmax(w)) floor = max(wmax, 1.0) * 1.0e-8 - if np.all(np.isfinite(w)) and float(w.min()) >= floor: - return Ms # already SPD → no-op - w = np.clip(np.nan_to_num(w, nan=floor, posinf=wmax, neginf=floor), - floor, None) - return np.einsum('nij,nj,nkj->nik', Vc, w, Vc) + # Per-tensor SPD test: a cell is "bad" only if one of its OWN + # eigenvalues is non-finite or below the floor. Project just those + # cells; every already-SPD tensor is returned untouched (bit-identical + # to the symmetrised input), so one bad point cannot perturb the rest. + bad = ~np.isfinite(w).all(axis=1) | (w.min(axis=1) < floor) + if not bad.any(): + return Ms + out = Ms.copy() + wf = np.clip(np.nan_to_num(w[bad], nan=floor, posinf=wmax, neginf=floor), + floor, None) + out[bad] = np.einsum('nij,nj,nkj->nik', Vc[bad], wf, Vc[bad]) + return out def _eval_M(pts): return _spd_sanitise(_eval_M_raw(pts))