Release v3.1.0#349
Open
jcgraciosa wants to merge 648 commits into
Open
Conversation
smooth_mesh_interior(metric=, method='spring'|'ma', boundary_slip=) — two topology-preserving node-redistribution solvers toward a target size/density field: - volumetric elastic spring (default): equal edge springs (shape, equant cells, no slivers) + per-cell area constraint A0 proportional to 1/rho_tgt (size); PCG truss-energy minimisation with an inversion-rejecting line search. ~0.3s, robust, scales with amplitude. - Benamou-Froese-Oberman Monge-Ampere: convex-branch det(I+D2phi)=g damped Picard on a variationally-recovered Hessian, pure-Neumann + the SNES_Scalar.constant_nullspace hook. Robust isotropic refinement, ~60x costlier. Optional tangential boundary slip (per-ring radius projection; nodes provably stay on the surface). Honest grading ~1.3-1.8x at fixed node count (extreme OT grading needs more nodes / a topology change). Docs: docs/developer/subsystems/mesh-metric-redistribution.md. Underworld development team with AI support from Claude Code
Validation and benchmark scripts for the metric smoother and the upcoming Monge-Ampere efficiency / convergence-rate work: honest per-node grading metric with mesh plots, exact radial equidistribution ground truth, setup sanity probe, cost / convergence / preconditioner benchmarks, boundary-slip A/B, interior-feature refinement test. Underworld development team with AI support from Claude Code
The Monge-Ampere path's cost was dominated by re-running a full GMRES+GAMG setup on every Picard iteration even though the mesh (hence the phi-Poisson Laplacian and the SPD Hessian-recovery mass matrix) is constant within the loop -- only the RHS changes. The constant_nullspace near-nullspace re-attach forced the GAMG re-setup, and newtonls added a redundant 2nd iterate + line search. _use_direct_solver puts the cached phi / Hessian / grad-phi sub-solvers on ksponly + MUMPS LU with snes_lag_jacobian=-2 and snes_lag_preconditioner=-2 (factor once, reuse: inner iters become back-substitutions), plus mat_mumps_icntl_24=1 for the singular pure-Neumann phi. n_picard default 40->25 (deep/near grading is flat from iter ~20). A direct solve is exact (tighter than the GMRES rtol) so the Picard fixed point -- hence grading/quality -- is bit-for-bit unchanged: d/n 1.02/1.43/1.71/1.54 (AMP 0/2/8/20), AMP=0 exact no-op, no tangle. Cost: cold 12-18s->1-2s, warm <=34s->1-2.5s; the warm>>cold GAMG-resetup degradation is eliminated. Serial MUMPS is the expedient that proved the diagnosis; the portable lever is the factor/setup-once-reuse pattern (port to GAMG for parallel) -- see the doc Implementation note + Open items. Validator: scripts/ma_cost_grading.py; profiler: ma_profile_phases.py. Underworld development team with AI support from Claude Code
Quantified the user-proposed quasi-Newton on the cofactor form cof(I+D2phi):D2(dphi)=g-det(I+D2phi) offline (fixed geometry, res-16 Annulus, transport-map metric), vs the shipped BFO-Picard. Decisive negative: (1) same recovered Hessian => same large irreducible det(I+H_rec)-g floor, so it cannot beat the settled fixed-node grading cap (predicted, not a grading lever); (2) it is LESS robust than BFO at the same Hessian quality -- all three standard convexity safeguards fail to reach BFO's fixed point (residual line-search freezes lambda->8e-3; det>0 backtrack stalls under-deformed d/n 1.58; PD-projection of eig(I+H) creeps, never contracts, overshoots into a cell-inverting map d/n 1.49) while BFO reaches d/n 1.713 in 16 iters. Root cause: BFO's sqrt(...)-2 is a closed-form convex-branch solve using only the deviatoric recovered Hessian + g, side-stepping the noisy under-estimated full det; cofactor-Newton feeds the full noisy recovered H into both C_k and the residual. UW3 forbids 2nd derivatives of mesh-var functions so a sharp D2phi is unavailable. Do not re-explore. Parallel path is therefore NOT Newton: keep BFO (robust convex-branch) and port the proven factor/setup-once-reuse pattern to GAMG (snes_lag_jacobian=-2 / KSPSetReusePreconditioner + constant near-nullspace + Krylov warm-start). Design doc + ma_newton_phase0.py. Underworld development team with AI support from Claude Code
…-reuse prototype Two findings from challenging the GAMG-reuse parallel prototype: 1. phi order, not the solver, sets the grading. P2 == P3 to ~3 dp across AMP 0/2/8/20 on the direct path (matches the recorded baseline, AMP=0 no-op exact, no tangle) at ~2x lower cost -- smaller matrices also help the direct factorisation scale. P1 is NOT grading-equivalent (~1.40 vs 1.71 at AMP=8, ~18% weaker), so P2 is the floor. phi_degree is now a parameter, default 3 -> 2. Canonical cost_compare at P2: MA cold ~0.7-0.9 s (was ~12-18 s originally) -- ~15-20x combined, grading bit-for-bit. 2. The BFO+GAMG-reuse parallel path (linear_solver="gamg", selectable; "direct" stays default). FGMRES+GAMG(SOR) for the elliptic phi-Poisson (CG unjustified -- no exact-symmetry guarantee + non-symmetric SOR smoother), CG+Jacobi for the SPD mass systems; snes_lag_*=-2 + Krylov warm-start. The constant nullspace is verified correctly wired (not the failure mode); P3 was a major GAMG confound (P2+gamg converges where P3+gamg diverges to a no-op) but even at P2 the warm post-_deform_mesh re-solve stays erratic, and this build has no alternative AMG (hypre/ML absent). Accepted position: MUMPS direct (itself MPI-parallel) for now; the P2 size cut only helps it. A robust iterative path would need a single-Dirichlet-pin (not the constant nullspace) and/or hypre, and is gated behind parallel-exact assembly + 3D (the solver is not the parallel bottleneck yet). Spring-as-MA-initial-guess: settled-rejected (geometric mechanism, unchanged by phi-order/solver speed) -- not re-run. Design doc + subsystem doc updated; validators ma_phi_order.py, ma_phi2_validate.py, ma_gamg_vs_direct.py, ma_solver_scaling.py. Underworld development team with AI support from Claude Code
…t warm-cost scaling P1 vs P2 × GAMG/direct, RES 16→64 (ma_p1_gamg_scaling.py, AMP=8): - P1 does NOT make GAMG robust. When P1+GAMG converges it is textbook-good (18-22 iters, N-independent — genuinely more AMG-friendly than P2's 77→103) but it still fails erratically (P1+gamg diverges res-32 & res-64; P2+gamg fails res-16 & res-32). The pure-Neumann + warm-resolve breakdown is order-independent and resolution-erratic; MUMPS direct is robust at every (res, order). - P1 grading ~1.40 vs P2 ~1.71 at every resolution (~18% weaker everywhere) — not a grading option regardless of solver. - Key actionable side-finding: the DIRECT warm cost scales badly with N — P2-direct warm 1.3 s (res-16) → 17.8 s (res-48) → 46.4 s (res-64), far above cold (9.5 s). The res-16 warm≈cold result does not extrapolate; per-call post-_deform_mesh rebuild + MUMPS refactorisation + cache-invalidated evaluate() re-interp are O(N)-growing. This is the next per-timestep-scaling work item. Design doc + memory updated. Underworld development team with AI support from Claude Code
Visual confirmation of the phi_degree decision: undeformed vs P1/P2/P3
MA redistribution (res-16, AMP=8) + outer-rim zoom. Shows P2 and P3
are visually identical strong clean grading (d/n 1.707 vs 1.713),
P1 a visibly weaker nudge (1.397), all meshes valid (no slivers /
inversion). Figures: /tmp/metric_mesh/ma_grids{,_zoom}.png.
Underworld development team with AI support from Claude Code
User caught that the P2 rim cells are far tighter than the nominal 1/3. ma_radial_anisotropy.py (res-16, AMP=8, vs undeformed): the deep/near ≈1.71 is a per-node mean of ALL incident edges and averages the collapsed radial edges with the frozen/expanded tangential ones. Band-mean radial at the rim ≈0.38× (~1/3, the isotropic-edge-criterion prediction) but the MIN radial edge ≈0.14× (~1/7) and minA/meanA ≈0.019 (thinnest tri ~1/52 mean area). All collapse is radial; tangential ~frozen. Mechanism: the outer ring is pinned (boundary) and the metric peaks exactly at r=R_O, so equidistribution jams the next ring(s) against the fixed wall into one near-degenerate sliver layer, independent of AMP. The isotropic AMP=1/s^2-1 rule is wrong here (annulus transport is purely radial; boundary-peaked-vs-pinned over-collapses the wall layer). d/n is OK as a regression guard but does NOT certify mesh quality near a boundary-peaked feature -- use minA/meanA or a radial/tangential split. Levers documented (offset metric peak inward / cap AMP to a quality floor / pinned-boundary 1-D radial OT design). Design doc + memory updated; figure /tmp/metric_mesh/ma_radial_profile.png. Underworld development team with AI support from Claude Code
…pendent
ma_show_grids_hires.py renders undeformed/P1/P2 at res 32 & 48
(AMP=8), full + outer-rim zoom, annotated with d/n + minA/meanA +
min-radial-vs-undeformed. Confirms the P2 rim sliver does NOT
improve with resolution: min radial 0.14x (res-16) -> 0.10x
(res-32) -> 0.11x (res-48); minA/meanA stays ~0.02-0.03. It is a
structural pinned-boundary + boundary-peaked-metric artifact, not a
discretisation error. P1 stays clean at all resolutions
(minA/meanA ~0.2). Figures /tmp/metric_mesh/ma_grids_hires{,_zoom}.png.
Underworld development team with AI support from Claude Code
…the snuggle fix Interior blob (away from the pinned boundary), AMP=8: - GAMG is ROBUST for localised cases (revises the blanket "GAMG fragile"). Every metric shape x width x resolution converged in ~27-54 iters, cost competitive with direct, zero failures. The earlier fragility was SPECIFICALLY the boundary-peaked-metric- against-pinned-boundary pathology (metric spiking where the operator is pinned). For the realistic localised-feature use case the parallel GAMG path is viable. - "Snuggle" fix: a wider Gaussian is wrong (one width couples the resolution scale and the reach: narrow = sharp but isolated pucker, bulk idle; broad = global motion but feature washed out, far/near 2.4->1.5). A heavy-tailed Lorentzian monitor 1+AMP/(1+d^2/Wc^2) gives the best feature resolution (far/near 2.74 > narrow 2.42) AND ~3x more inward migration of distant nodes (+0.025 vs +0.008): the whole mesh rakes coherently toward the feature. Mild quality cost (minA 0.089 vs 0.105), no tangle. Standard r-adaptation lesson (monitor needs global reach). Diagnostic-only (no src/default change). Design doc + memory updated; figures /tmp/metric_mesh/ma_heavytail.png, ma_localised_reach.png. Underworld development team with AI support from Claude Code
Tested defining the metric in (r,theta) so it "pulls in theta" plus boundary sliding (reference: compact Cartesian Lorentzian at an interior point gives far/near 2.74). 1. A separable (r,theta) Lorentzian is the wrong shape -- an anisotropic spoke, not a blob: the chord 2(1-cos dtheta) saturates at the antipode (no angular reach, far/near 1.38); true-wrapped-angle + balanced cores is a low-gradient radial ridge the smoother washes out (far/near 1.12, ~no-op). Use a compact |X-P|^2 Lorentzian about the feature point -- it has the correct combined radial+angular extent and pulls in theta inherently. 2. boundary_slip works mechanically (rim radial drift ~1e-16, nodes provably stay on the ring; GAMG robust ~31 it) but is NOT a concentrator: near a boundary feature it relaxes the mesh (far/near 1.21->1.12, minA 0.32->0.48, rim-count near theta0 16->18) -- it removes the hard pin so the rim equalises, it does not drag rim nodes tangentially toward the feature. Slip buys boundary quality, not feature concentration. 3. Boundary-proximal features are choked: same compact Lorentzian gives far/near 2.74 at r0=0.78 (interior) vs 1.21 at r0=0.88 (near rim) -- the fixed-node + pinned-boundary cap, feature side. Recommendation: compact Cartesian |X-P|^2 Lorentzian about the feature point, keep features with interior room, slip only for boundary quality; drop the polar-separable idea. Diagnostic-only. Design doc + memory updated; figs /tmp/metric_mesh/ma_polar_slip*, ma_lorentzian_slip.png. Underworld development team with AI support from Claude Code
…al dead end User's reframing: exploit the abundant tangential node budget (slide spare angular nodes toward the feature) not the scarce pinned radial one. Built (1) the exact 1-D angular OT map as the target and (2) a new opt-in _winslow_elliptic(move_anisotropy= (w_r,w_theta)) that rescales the realised displacement in the local radial/tangential frame (default None = byte-identical). Angle-only feature rho(theta), AMP=8, res-24: - (1) exact angular OT: far/near 2.21, frac@theta0 0.159->0.415, radial drift 1e-16, no tangle. Exact, trivial, robust -- the right tool for separable/structured features, used directly. - (2) scalar BFO: ~ZERO angular concentration (far/near 0.98, frac 0.158 ~ uniform) for ANY weighting. move_anisotropy works as designed (suppresses spurious radial drift 6.8e-2 -> 8e-3) but there is no tangential transport to reweight -- the scalar potential never generates it. Root cause = the foundational ~1.7 cap, both directions: a scalar equidistribution potential with fixed topology cannot do large coherent bulk transport, radial OR tangential (hoop/fixed-topology stiffness cuts both ways). "(1) as a target for (2)" proved (2) cannot reach it. Generalisable path: explicit dimensional-split 1-D OT for separable features, or a true anisotropic metric-tensor adaptation (heavy) -- NOT anisotropic diffusivity / move-weighting on the scalar potential. move_anisotropy is kept as an opt-in quality knob (off-direction drift suppressor), not a concentrator. Design doc + memory updated; fig /tmp/metric_mesh/ma_angular_ot.png. Underworld development team with AI support from Claude Code
ma_metric_tensor_viz.py: scalar density -> M = (1/h0^2)[I + beta ghat ghat^T (|grad rho|/ref)^2], eigen-clamped, desired-cell ellipses on a polar sample grid for a radial rho(r) and an angular rho(theta) feature. Correct and confirms the design: - radial feature -> tangentially elongated ellipses (short across r), angular feature -> radially elongated (short across theta), each concentrated at the feature and isotropic where grad rho->0. - The eigenframe auto-aligns to r-hat / theta-hat with NO (r,theta) frame specified anywhere (M sees only the Cartesian gradient) -- resolves the user's (r,theta) puzzle; the scalar API is preserved. - max anisotropy = the eigen-clamp band (8.3:1), as designed. Honest nuance: a gradient-based metric refines where rho CHANGES (flanks), isotropic at a smooth peak / far field; for small cells at the feature core use smoothed |grad rho| or Hessian-based M=|H(rho)| (needs the recovered-Hessian path). Gradient form is the first-derivative UW3-clean first cut. The metric construction (the ~1-day half) is verified and cheap; remaining for (3) is the anisotropic mover (metric-Winslow / M-weighted displacement solve), improving alignment/quality not the fixed-node cap. Design doc updated; fig /tmp/metric_mesh/ma_metric_tensor.png. Underworld development team with AI support from Claude Code
Self-contained handoff so a fresh session can start building the metric-Winslow / M-weighted displacement solve without re-deriving: what to read (memory dead-ends, design doc, the _CofDiff tensor- constitutive pattern + cache/lag/MUMPS infra), the concrete plan, the validation targets (anisotropy-aware diagnostics + 1-D OT reference), the accepted quality-not-cap caveat, and scope (~1-2 weeks, new branch). Metric construction is done/verified; only the mover remains. Underworld development team with AI support from Claude Code
_winslow_anisotropic: M-weighted Laplace (Winslow) smooth of the
coordinate map with an eigen-clamped, gradient-derived anisotropic
metric tensor D. Displacement form ∇·(D∇u_c)=-Σ_j∂_j D_{jc} solved
per coordinate component, sharing the _CofDiff-style DiffusionModel
tensor operator (_c=D) and the factor-once-reuse direct solver;
homogeneous Dirichlet on the pinned boundary makes it non-singular
(no constant_nullspace — side-steps the GAMG-pure-Neumann
fragility). Linear (one solve/component/outer step, no Picard).
Reuses _winslow_elliptic's signed-area backtrack, boundary_slip and
move_anisotropy. Wired into smooth_mesh_interior as the opt-in
method="anisotropic"; spring default + ma path untouched.
Critical no-op guard: the L2 projection of a uniform-ρ zero
gradient leaves ~1e-18 round-off; normalising by that noisy max
fabricated O(1) anisotropy. A scale-aware g_eps floor makes AMP=0
an exact isotropic no-op while AMP>0 is bit-identical to the
verified ma_metric_tensor_viz construction.
Smoke (res-16 Annulus): AMP=0 exact no-op (max|Δx|=0), AMP=8 moves
nodes, mesh stays valid via the backtrack, ~3s/5 outer steps.
Underworld development team with AI support from Claude Code
Two formulation fixes after the radial-feature validation arc: 1. Build the eigen-clamped metric tensor D ONCE on the undeformed mesh and hold it fixed & Lagrangian (rides material points), exactly as _winslow_spring computes its rest-lengths/A0 once. Re-projecting ∇ρ on the progressively distorted mesh each outer step is a positive feedback — D blew up on squashed cells → catastrophic over-collapse (verified). With D fixed the outer loop is a stable damped fixed-point iteration of one linear operator toward the M-harmonic map. 2. Under-relax the per-step displacement (relax). The decoupled direct Winslow form has no Rado–Kneser–Choquet non-folding guarantee, so a single un-damped elliptic jump folds; its stable regime is bounded by the metric anisotropy. Characterised the Pareto frontier (scripts/aniso_param_sweep.py): aniso_cap=2 + relax≈0.1-0.2 is the robust default; higher cap needs gentler relax + more n_outer; cap≳6 folds regardless (would need the coupled/inverse Winslow — out of prototype scope). Defaults aniso_cap 8→2, relax 0.4→0.2. Validated (scripts/aniso_validate_radial.py, anisotropy-aware radial/tangential split + minA/meanA, NOT d/n, vs MA + spring + exact 1-D OT; grids rendered). Interior radial feature (r=0.70, res-16 AMP=8): (3) minA/meanA 0.47 vs isotropic MA 0.18 / spring 0.25 — 2.6x cleaner, no slivers, still redistributes toward the feature. Even on the pathological boundary-peaked case (r=R_O) (3) keeps minA/meanA 0.24 vs MA 0.019 (12x less degenerate). (3) trades grading MAGNITUDE for clean anisotropic cell ALIGNMENT — exactly its intended role; does not beat the fixed node-count cap. Underworld development team with AI support from Claude Code
Completes the (3) validation arc with the anisotropy-aware diagnostics on the remaining model problems. Angular-only feature (res-24 AMP=8, scripts/aniso_validate_ angular.py, vs the exact 1-D angular-OT target): (3) far/near 1.02, frac@θ0 0.160 ≈ uniform — like the scalar MA it CANNOT do separable angular concentration (expected: separable ⇒ the explicit 1-D OT is exact + cheap; the settled fixed-topology cap). But (3) keeps minA/meanA 0.243 vs scalar MA 0.144 — cleaner. Non-separable blob P=(0.78,0) (the case (3) is for, scripts/aniso_validate_nonsep.py): (3) minA/meanA 0.295 vs MA 0.109 / spring 0.119 — 2.7x cleaner, no slivers; the zoomed grid shows MA/spring pull a degenerate slivered knot into the blob while (3) gives a clean, well-shaped, blob-aligned densification. Concentration milder (far/near 1.10 vs MA 1.37) — (3) trades magnitude for clean alignment, exactly its intended role. Verdict across the arc: (3) is the cleanest method everywhere (2.6-12x better minA/meanA than isotropic MA), never slivers, is linear/cheap; it does NOT beat the fixed node-count cap (for separable features the explicit 1-D OT remains exact + cheaper). A validated prototype matching the kickoff brief. Underworld development team with AI support from Claude Code
Design doc: append the "(3) anisotropic mover — IMPLEMENTED & VALIDATED" section (formulation, the two settled formulation findings — build-D-once feedback + the decoupled-form folding / aniso_cap Pareto, the validation table across all model problems, and the verdict + open follow-ups). Subsystem doc: extend "the two solvers" → three (method="anisotropic" row + recommendation + Open items entry + example). Project memory (project-anisotropic-mover) records the non-obvious settled findings so future sessions don't re-derive them: build D once (re-projection is a feedback collapse), aniso_cap is the binding lever (cap=2 robust, ≥6 folds — coupled/ inverse Winslow needed, out of scope), the g_eps no-op guard, and the honest "cleanest-not-strongest, doesn't beat the node-count cap" verdict. Underworld development team with AI support from Claude Code
Answers a review question: a GRADIENT metric resolves a Gaussian blob's EDGE, not its CORE — by design. Renders the desired-cell ellipse field (exactly the shipped mover's construction: β=200, aniso_cap=2, h0 from mesh) beside the realised method="anisotropic" mesh and an edge-length/cell-aspect profile vs d=|X-P|. Confirms: target = fine anisotropic shell at the |∇ρ|-max ring, isotropic-coarse core + far. Realised mesh matches — core (d≈0) mean edge 0.079 vs undeformed 0.041 (~1.9× COARSENED: nodes evacuated from the ∇ρ=0 core into the edge shell), edge-length minimum + aspect maximum on the flank ring, far field untouched (0.0416 vs 0.0412). If the blob CORE must be resolved, the Hessian metric M=|H(ρ)| is the tool (documented follow-up, out of prototype scope). Underworld development team with AI support from Claude Code
"Run 20 steps non-adaptive, then refine." Reuses the fixed_mesh_convection setup (Annulus res-16, free-slip Stokes + AdvDiffusionSLCN, Ra=1e5, 20 steps, FIXED mesh; state cached to npz). Then builds the mover-side analogue of UW3's adaptation metric (adaptivity.metric_from_gradient): a Lagrangian density ρ = 1 + AMP·t, t = clip((|∇T|-g_lo)/(g_hi-g_lo),0,1) with the same percentile window idea — fed as method="anisotropic" metric=. |∇T| is pre-projected to a scalar field so the mover's internal metric.diff(X) stays a first derivative (UW3-legal) and the metric rides material points. Result (fig /tmp/metric_mesh/aniso_convection.png): the metric correctly targets the inner thermal BL + plume edges; the refined mesh gathers cells there (anisotropic, aligned along the BL/plume, short across ∇T) and stays clean — minA/meanA 0.57→0.25, no slivers. Fixed node budget ⇒ relative redistribution (the documented node-count-cap caveat): same metric *form* as adaptation, budget-limited *magnitude*. Underworld development team with AI support from Claude Code
Fold the pipeline walkthrough (cache-build components 1-4, per-call build-D-once + damped MMPDE loop steps 5-6), the solver limitations (2D-tri-only, decoupled-direct folding bound, node-count cap, gradient-resolves-edges-not-cores, Lagrangian-fixed metric, serial-exact assembly, linear/decoupled), and the unexplored corners (no solution-accuracy proof yet, no dynamic-adaptive loop, coupled/inverse Winslow, Hessian metric, the metric_from_gradient-style ρ helper, GAMG path, auto-tuning, free-surface slip) into the design doc's "(3) … IMPLEMENTED & VALIDATED" section as the canonical reference. Underworld development team with AI support from Claude Code
scripts/aniso_cost_and_gamg.py: interior radial feature, res 16/24/32/48, direct vs gamg, decomposing cold (one-off setup) / warm (per-timestep) / per-outer-step / D-build. Underworld development team with AI support from Claude Code
GAMG: bit-parity with direct at res 16–48 (the mover is non-singular — homogeneous Dirichlet — so no pure-Neumann fragility; first metric method with a working parallel-capable solver path). Cost: cold≈warm (no MA warm-rebuild pathology), ~O(N) in #cells, per-step = a fixed number of SPD elliptic solves + local work — structurally parallel-cheap (no nonlinear solve / global transport). Lever is n_outer. Folded a "GAMG parity + cost per step (measured)" subsection into the design doc; updated the Solver-limitations + Corners (GAMG validated → remaining is parallel-exact assembly; concrete 3D scope: solver core already dim-general, the 2D-specific piece is the tet signed-volume backtrack). Subsystem-doc Open-items entry updated. Underworld development team with AI support from Claude Code
aniso_convection_res32.py: rebuild the res-16 P3 T from the cached 20-step npz, interpolate onto res-32 T nodes (uw.function.evaluate across meshes), re-solve Stokes for a consistent V, run 5 warm steps (res-32 state itself cached), then refine on ρ∝|∇T| with method="anisotropic". Warm start preserves+sharpens the res-16 plume structure; the mover snuggles points into the inner BL + plume conduits, mesh stays valid (minA/meanA 0.55→0.19, no slivers — finer |∇T| ⇒ a more aggressive but clean redistribution than res-16's 0.25). Fig /tmp/metric_mesh/aniso_convection_res32.png. Underworld development team with AI support from Claude Code
- smooth_mesh_interior gains method_kwargs: a clean pass-through for the per-method tuning (aniso_cap/relax/n_outer/linear_solver /beta/move_anisotropy) without bloating the shared signature; documented with the validated anisotropic defaults. - New public uw.meshing.metric_density_from_gradient(mesh, field, amp=, lo/hi_percentile=, ...): encapsulates the demonstrated pattern (project ∇field — first-derivative, UW3-clean — percentile -normalise, freeze as a Lagrangian scalar) and returns the ρ density expression. The relative, fixed-node-budget analogue of adaptivity.metric_from_gradient. Cached per (mesh,degree,name, topology) so it is cheap + leak-free called every step in an adaptive loop. Scale-aware g_eps floor ⇒ a (near-)uniform field yields ρ≡1 exactly (the smoke test caught the percentile-of- round-off fabrication — same class as the mover's own guard). - Exported from underworld3.meshing. aniso_api_smoke.py: public import, helper cache, method_kwargs pass-through (bad kwarg → TypeError), valid clean move, exact uniform no-op — all pass. Underworld development team with AI support from Claude Code
mesh-adaptation.md: add a "Two families of mesh adaptation" table (MMG re-mesh / absolute h / variables reset vs node-redistribution / relative density / topology preserved / anisotropic / cheap- parallel) in Core Concepts, and a "Node redistribution — the snuggling mover" section with the metric_density_from_gradient + smooth_mesh_interior(method="anisotropic") recipe, the gradient-resolves-edges-not-cores + node-count-cap caveats, the method_kwargs knobs, and cross-links to the subsystem + design docs. Frames it as a peer adaptive strategy with the SAME intent API as adaptivity.metric_from_gradient. Underworld development team with AI support from Claude Code
- scripts/adaptive_convection_harness.py: Ra=1e5, uniform res-24 reference vs res-16 + periodic anisotropic adaptation (locked-in API), matched-time Nu(t)/vrms(t) rms comparison + figure, both histories cached. The node-update/ALE correction is an explicit documented hook (apply_adaptation_correction): --correction none = uncorrected baseline (quantifies the drift the next phase must remove); --correction ale raises with the precise spec (V_fn = v - v_mesh post-adapt, per the settled free-surface ALE precedent). - Design note: "NEXT-PHASE KICKOFF BRIEF — dynamic adaptive convection": phase-closed summary (API/docs/harness entry points), the ALE correction as the headline open piece + acceptance test, prioritised follow-ups (3D scope, parallel, Hessian metric, auto-tune), and the resume recipe. Records that the metric is 1/h^2 per principal direction in 2D AND 3D (NOT 1/h^dim) — verified adaptivity.py:154. Underworld development team with AI support from Claude Code
adapt_with_correction(mode="interp"): exploit that the anisotropic mover is topology-preserving (vector size / DOF order / partition invariant) and the pinned boundary leaves the domain unchanged — so the correct remap is the OLD P3 field evaluated at the NEW node positions via the local FE evaluate (uw.function.evaluate, true basis, max fidelity, no cross-rank migration), done by a deform-back -> evaluate -> deform-forward around the mover, then re-solve Stokes. NOT read_timestep's kd-tree path (that exists for the decomposition-changing mesh.adapt case and was diffusive — degraded the solution). Result (Ra=1e5, res-16 adapt-every-5 vs uniform res-24 ref, n15): rms ΔNu 0.028 vs uncorrected 0.097 (~3.5x better); the per-adapt Nu kink is eliminated and the trajectory tracks the reference. vrms shows a small expected interpolation-diffusion under- prediction (the remap tax — motivates the future ALE option). Also: nusselt() sign fixed (was -Nu: divided by the signed conductive gradient instead of the conductive flux); harness overlays the cached uncorrected baseline; --correction gains "interp"; "ale" raises as the next-phase adapter-bound option. Underworld development team with AI support from Claude Code
Second submesh flavour alongside the existing rock/air subdomain (extract_region): pull any level out of a mesh's nested refinement hierarchy as a standalone solver-ready uw.Mesh with full parent lineage, and transfer fields back and forth. Contract: refine-DM mode only. Available only when a genuine nested refinement hierarchy exists; transfer uses PETSc's nested interpolator/injector (exact, parallel-local, no point location). No geometric/General fallback, no KDTree -- a non-refined mesh raises rather than silently degrading. Key construction: build the transfer pair by refining a single-field clone of the coarse level (dm_f = dm_c.refine()), so refine() itself establishes the setCoarseDM linkage and regular-refinement flag and createInterpolation/createInjection take the nested exact branch. Independently cloning two hierarchy levels breaks this and petsc4py exposes no setter to restore it. Includes: - refined_pair_prototype.py: coarsened_companion + prolongate/restrict/ inject/sample, refine-DM contract guard - example_refined_companion.py: sibling of test_region_ds_submesh.py, same annulus+radial-buoyancy Stokes so the two submesh flavours are directly comparable - test_refined_pair_solver.py: gating test (extract -> labels -> solve -> map back), round-trip recovers coarse to ~6e-16 - test_refined_pair_contract.py: enforces no-refinement -> raises - probe_hierarchy_labels.py: boundary labels survive every hierarchy level (box, annulus-with-internal-boundary, spherical shell) - submesh-solver-architecture.md: "Two Submesh Flavours" section No source changes; investigation lives under docs/examples/. Underworld development team with AI support from Claude Code
- Drop unused numpy import - Fix stale cache-value-shape comments to match the stored tuples - Key interp/inject caches by variable objects via nested WeakKeyDictionary instead of id() (id reuse after GC could alias a stale entry with a mismatched FE layout) - Invalidate _canonical_data AND _cached_data_array on direct PETSc writes (matches the core _re_extract_from_parent invalidation; .array /.data are copies that otherwise go stale) - Destroy scratch global Vecs after each transfer (try/finally) so repeated per-timestep transfers don't accumulate PETSc objects; cached Mat/DM are built once and kept Verified: contract, gating (round-trip 5.9e-16) and annulus example all still pass unchanged. Underworld development team with AI support from Claude Code
…ngs (#214) (#255) The genuine bug: SNES_TransientDarcy.solve projected +darcy_flux for the velocity field, where darcy_flux is the ASSEMBLY flux F = kappa(grad(h)-s). The physical Darcy velocity is v = -F = -kappa(grad(h)-s) (steady SNES_Darcy already used -darcy_flux), so transient runs reported a sign-flipped velocity. Fixed to -darcy_flux to match steady. Docstrings (issue #214): - SNES_Poisson: equation written +div(kappa grad u)=f but the solver assembles -div(kappa grad u)=f -> corrected both docstring instances. - SNES_Darcy / DarcyFlowModel: stated the velocity as +kappa(grad h - s); the code (correctly) returns v = -kappa(grad h - s). Clarified that the model flux property is the assembly flux F and the physical velocity is -F. The assembly flux F1 = kappa(grad(h)-s) is UNCHANGED (it already solves -div(kappa(grad h - s)) = f correctly for all f, including f != 0). Tests (test_1004b): v.grad(h) < 0 (down-gradient) for steady AND transient (the transient assertion fails on the pre-fix +darcy_flux, verified); plus an f != 0 manufactured-solution test locking the source/assembly sign (the existing Darcy tests all use f=0). Closes #214. Supersedes #243 (which flipped the flux instead, introducing an untested f != 0 source-sign regression). Underworld development team with AI support from Claude Code
…ng + drop dead cell_dim read (#252) (#256) Follow-ups from the #218 review: - disk_snapshot._write_disk_snapshot is a legitimate internal user of the user-deprecated write_checkpoint (it depends on the .mesh.00000.h5 / .<var>.00000.h5 filename convention its reload path reads). Wrap the call in warnings.catch_warnings()+ignore FutureWarning so persistent snapshots don't spam the deprecation warning. (Full migration to write_timestep tracked in #252.) - checkpoint_xdmf read topo["cells"].attrs["cell_dim"] into cellDim but never used it (topology_precision now derives from cells.dtype.itemsize). Removed — it was a needless dependency on a PETSc-internal attr that could KeyError. Verified: test_0007/0008 (snapshot) pass with -W error::FutureWarning (26/26); test_0003/0005 (save_load/xdmf) green (39/39 total). Item 3 (unify swarmVar HyperSlab XDMF block with the flattened meshVar form) deferred: cosmetic, ParaView-sensitive, no easy local validation — left in #252. Underworld development team with AI support from Claude Code
…eformed meshes (#259) _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)) defined only for an SPD metric. The metric is FE-evaluated at a fixed reference cloud; after the interior deforms a reference point can fall outside the mesh and the P1 density is evaluated by FE extrapolation, going negative — M=rho*I with rho<0 has det=rho^2>0 (passes a detM>0 check) but is negative-definite, so S<0 and S**q=NaN, the energy is non-finite, and the mover bails with zero displacement (adaptation silently stops on ~half the adapts). Project each evaluated tensor onto SPD with a small relative eigenvalue floor at the _eval_M chokepoint: per-tensor test, only the failing cells are projected, already-SPD tensors returned bit-identical, empty-rank safe. No API change, no behaviour change on valid metrics. Tests: test_0762/0850/0750/0830 green. Validated on the stagnant-lid adaptive convection driver (all forced adapts move, T in [0,1], vrms bounded). Underworld development team with AI support from Claude Code
…ed mesh (#264) * fix(mesh): boundary normals and domain membership must track a deformed mesh Two related bugs where analytic/cached UNDEFORMED geometry survived a mesh.deform(), corrupting Nitsche/penalty free-slip BCs and semi-Lagrangian advection on free surfaces. (1) Projected boundary normals stayed radial on a deformed surface. mesh.Gamma's base scalars carry the C-code petsc_n[] (PETSc facet normal), defined only inside surface-integral kernels; _update_projected_normals point-evaluated it, so it fell back to the coordinate -> (x,y)/r = radial. Gamma_P1 therefore never tracked deformation and add_nitsche_bc constrained v.r_hat instead of v.n_hat_true, leaking throughflow proportional to surface tilt. Fix: new Mesh.boundary_normal(label) assembles the exact PETSc facet normals (computeCellGeometryFVM, outward, area-weighted) from ONLY that boundary's faces -> per-boundary (a corner shared with another boundary keeps its own one-sided normal, never an average across the discontinuity), deformation-tracking, smooth on a smooth boundary. add_nitsche_bc (Stokes and Scalar) now uses it; global Gamma_P1 is unchanged for back-compat. The field is refreshed on every deform so BCs that captured its .sym at setup read the new geometry. Parallel-safe (guards getStratumIS against dead label values on ranks holding no faces of a boundary). (2) points_in_domain()/function.evaluate() used the original boundary. The boundary-skeleton kd-tree was built from mesh._nav_coords, captured as a reference to the ORIGINAL coords in __init__ and never refreshed in nuke_coords_and_rebuild (only the per-cell control-point arrays were, on adapt). On a volume mesh _nav_coords stayed at the undeformed boundary, so a region that bulged past the original surface read as EXTERIOR -- stranding SL trace-back feet there and mis-locating evaluations, injecting the cold boundary value at topographic highs. Fix: invalidate the boundary kd-tree and refresh _nav_coords from the current DM coordinates on every deform. Tests: test_0056 (boundary_normal tracks a deformed annulus + a Cartesian free-surface corner) and test_0057 (deformed-domain membership + evaluation), both tier_a. test_1060 Nitsche free-slip passes serial and np2. Follow-up: ADD-reduce the unnormalised facet contributions at partition-seam boundary vertices before normalising (those few vertices currently use a this-rank-only stencil; parallel-safe and correct elsewhere). Underworld development team with AI support from Claude Code * Address Copilot review on #264 - boundary_normal(): drop unused cKDTree import (tree is built in _assemble_boundary_normal). - _assemble_boundary_normal(): drop unused `import underworld3 as uw`; count P1 vertices PER FACET from its closure (vStart..vEnd) instead of assuming nverts == cdim, so non-simplex facets (3D quad/hex boundaries) are handled correctly. - add_nitsche_bc docstring: `normal` default is now the per-boundary, deformation-tracking mesh.boundary_normal(boundary), not mesh.Gamma_N. - tests: drop unused surf_now; tidy an editorial comment fragment. test_0056 + test_0057 pass (3 passed). Underworld development team with AI support from Claude Code
Restores the intended "RBF metric eval" robustness. The mmpde mover's RBF/Shepard metric path baked the metric via an FE evaluation at the FIXED reference cloud `ref`; on a deformed mesh that reference can mis-locate / drift outside the deformed interior and the FE evaluation returns garbage (the P1 density — strictly positive by construction — comes back NEGATIVE, even at a field's own DOFs). A negative density is a non-SPD metric, which the mover then either NaN-bails on (a hidden stall) or, with the #259 SPD-floor, acts on as "coarsen hard here" → giant cells / holes / divergence in the adapted mesh. Fix: bake the metric at the CURRENT mesh NODES (its own DOF locations) and Shepard-interpolate from there. A Shepard (positive-weight, convex) average of the positive nodal values is GUARANTEED >= 0 (monotone) — the metric can never go non-SPD from the eval — and nodes are always inside the mesh (no out-of-domain / drift) and need no per-step cell location (fast). `ref` is kept for the _edge_mats reference frame. Best paired with #264 (deformed-mesh point-location fix), which makes the nodal FE bake exact; together the adapted mesh stays clean under forced every-step adaptation at R=5 (folded=0, cell-area-ratio flat) where it previously tore holes. Validated in a sibling worktree (this is a one-line source change; CI runs the full mover suite). Underworld development team with AI support from Claude Code
SemiLagrangian.update_pre_solve crashed for any unit-bearing model in the time
loop:
ValueError: Cannot subtract arrays with incompatible units:
'meter' and 'meter / second'
Root cause: the trace-back is performed in the mesh's NON-DIMENSIONAL (DM)
coordinate space (evaluate()/global_evaluate treat plain arrays as DM coords,
and DM point-location uses DM values 0..L_model, NOT dimensional metres). But
the has_units branch kept dimensional coords (mesh.X.coords ~ metres) and
velocity (m/s) and left dt a unitless float, so:
- coords[m] - v[m/s]*dt[dimensionless] -> unit subtraction crash, and
- even patched, dimensional coords (~1e9) fed to a DM that locates in 0..1000
space would mislocate.
Fix: route the whole trace-back through ND/DM space regardless of units —
departure point from psi_star[i].coords_nd (== .coords for a non-units model;
the ND reduction of the dimensional coords when units are active), velocity
non-dimensionalised, dt non-dimensional. Both the node-velocity and
midpoint-velocity legs are corrected. The non-units path is unchanged (it
already used exactly this ND logic).
Verified: units-active AdvDiffusionSLCN now runs and tracks the equivalent
non-dimensional run to ~1e-3 (the small residual is the constitutive diffusivity
scaling under units, a separate concern, not the trace-back). Non-units SL
regression unchanged (test_0855/0610/1054/1100 pass; test_1110 fails identically
with/without this change — a pre-existing MG DMCreateInjection issue). New
regression: tests/test_1056_units_slcn_traceback.py.
Unblocks the time-loop half of #263 (Tutorial_Thermal_Convection_Units).
Underworld development team with AI support from Claude Code
Companion to the units design doc. States the rule that the DM/PETSc/solver layer (incl. evaluate/global_evaluate, point-location, JIT residuals) is non-dimensional-only, units live above it in the MeshVariable API (.array/.coords dimensional vs .data/.coords_nd ND), the crossing functions are non_dimensionalise/dimensionalise, and the trap that evaluate strips the UnitAwareArray label but keeps the value (silent mislocation) — the #267 class of bug. Documents the canonical correct pattern (the #267 trace-back fix) and notes the enforcement guard as in-progress/tracked. Registered in the developer design toctree next to the units design doc. Underworld development team with AI support from Claude Code
…280) uw.Params reads overrides from the PETSc options database, but petsc4py only auto-populates that DB from sys.argv on some platforms (e.g. macOS) and not others (e.g. Gadi). On Gadi, Params(...) therefore silently fell back to its defaults even when the CLI argument was present — a hard-to-debug, environment-dependent failure with no warning. Params.__init__ now calls uw.parse_cmd_line_options() (idempotent) before reading the options DB, so a '-uw_<name> <value>' override is applied consistently. Regression: tests/test_0821_params_cli_override.py. Closes #111. Underworld development team with AI support from Claude Code
…, #158) (#281) #156: Projection (and its vector/tensor/multi-component variants) does a linear SPD L2 projection, but inherits the heavy SNES_Scalar newtonls/gmres/gamg stack — GAMG setup/repartition is the memory/communication bottleneck for repeated post-processing projections (OOM-killed mid-sequence on Gadi). Add an opt-in SNES_Projection.linear_solver(pc='jacobi', rtol=1e-10) that switches to ksponly+CG+cheap-PC and drops the unused GAMG options. The global default is unchanged (internal paths that rely on it are unaffected). Verified the lightweight solve matches the default projection; tests/test_0506. #158: document the 'project components, compose analytically' rule for boundary stress recovery — projecting the composed nᵀσn mixes derived-τ and discontinuous p at lower accuracy (~2x error on the Thieulot benchmark); project the τ components and form σ_rr = nᵀτ_proj n − p instead. New how-to under docs/developer/subsystems/ covering both rules. Underworld development team with AI support from Claude Code
* Fix Darcy Cartesian example script run * Ex_Darcy_Cartesian: use public mesh.deform() not internal _deform_mesh The example called mesh._deform_mesh(new_coords=...), which the coord-mutation capability gate (#246) now rejects once the mesh carries variables (the internal primitive moves nodes without transferring fields). Switched to the public mesh.deform(new_coords) API, which performs the sanctioned field-carrying deform. Verified the example runs end-to-end. Underworld development team with AI support from Claude Code --------- Co-authored-by: lmoresi <louis.moresi@anu.edu.au>
A MeshVariable's .sym lives in non-dimensional (model-unit) space, so the MathematicalMixin arithmetic operators could not compose with a unit-bearing scalar operand: T - uw.quantity(500, 'K') raised TypeError: Cannot subtract UWQuantity from EnhancedMeshVariable even for compatible units (the Pint quantity reached sympy and failed to sympify). This blocks the fundamental buoyancy form T - T_ref. Add a coercion branch to the operand handling in __add__/__sub__/__rsub__/ __mul__/__rmul__/__truediv__/__rtruediv__: a UWQuantity (anything with .magnitude and .units) is reduced to its non-dimensional value via uw.non_dimensionalise before composing with the variable's .sym. Compatible units (incl. prefixes like kK) reduce consistently via the model scaling; the change is behaviour-preserving for all non-quantity operands (sympy / variable / plain scalar hit the same branches as before — verified the sympy-scalar matrix quirk is pre-existing, not introduced here). Scope: variable on the LEFT (meshvar ± / * quantity). 'quantity - meshvar' (quantity on the left) is handled by the UWQuantity class and is out of scope. Verified: T ±/* quantity compose with exact ND values (T - 500K == T.sym - 0.5), kK prefix consistent, no regression (test_0700/0750/0754 unchanged). Regression: tests/test_0756_meshvar_quantity_arithmetic.py. Helps unblock the buoyancy expression in PR #263. Closes #282. Underworld development team with AI support from Claude Code
The Stokes bodyforce setter tried to handle UWQuantity components via item._sympify_() — a method that does not exist on UWQuantity (_sympy_ also raises on a dimensional quantity). So any units-active buoyancy assignment (stokes.bodyforce = [0, ρ₀ α g (T − T_ref)]) crashed with AttributeError: 'UWQuantity' object has no attribute '_sympify_'. The body force feeds the non-dimensional solver, so non_dimensionalise the UWQuantity component (consistent with the ND<->units boundary contract and the #282 fix). The component may be symbolic (a buoyancy carries the T field), so non_dimensionalise yields a 1x1 array/Matrix whose element is extracted (the existing shape-handling is preserved). Plain non-quantity components are unchanged. Verified: a symbolic UWQuantity buoyancy assigns to bodyforce as the correct ND expression; test_1010 (Stokes) unchanged. Regression: tests/test_1011_bodyforce_units.py. The second blocker (after #267/#282) for the thermal-convection units tutorial (#263). Underworld development team with AI support from Claude Code
…, rotation-gauge (#265) * Stokes_Constrained parallel correctness: gauge, convergence, knockout audit Three fixes from a focused parallel-correctness investigation of the in-saddle Lagrange-multiplier free-slip solver (Stokes_Constrained), validated against the Zhong (2008) spherical-shell response (3-D, SphericalShellInternalBoundary). Item 1 — consistent pressure gauge (raw-field reproducibility). On an enclosed constrained problem the constant pressure and constant multiplier are independent gauge freedoms; the solver lands on a partition-dependent level for each, so raw pressure/multiplier are not reproducible across ranks. Added a guarded automatic pressure gauge (auto_pressure_gauge, default on) that pins the surface-mean pressure on the first constraint boundary when a pressure null space is active and the user has registered no gauge of their own. It is physics- neutral (velocity bit-identical) and fixes the raw mean pressure (10.5% -> 0.4% serial vs np=8). The pressure pin does NOT fix the raw multiplier (separate gauge freedom): dynamic topography is read gauge-invariantly via topography(reference="mean"). New parallel regression in test_1063. Item 2 — the grouped-Schur solve was false-converging. The grouped u | [p,h] Schur preconditioner uses inner iterative sub-solves (a variable/nonlinear preconditioner). A non-flexible Krylov method (the base's gmres) is invalid for that and false-converged: it reported CONVERGED_RTOL in ~2 iterations on the preconditioned norm while the TRUE residual ||r||/||b|| blew up to ~1e3, landing on a wrong, partition-dependent answer (the ~0.4% velocity spread and the failure to reproduce the Zhong response). Fixed with four defaults in Stokes_Constrained: ksp_type=fgmres (flexible Krylov), ksp_norm_type= unpreconditioned (honest true-residual stopping test), Eisenstat-Walker rtol0=rtolmax=tolerance*0.1 (EW's default 0.3 capped the linear solve at one iteration), and a bounded ksp_max_it. The solve now genuinely converges (true residual ~1e-13) and reproduces Zhong (Ut within 1%, Ub within 0.1%); new test_1064 default-constrained test asserts this. Velocity partition spread 0.36% -> 0.12%. The residual 0.12% is a separate, pre-existing bug: add_natural_bc on an INTERNAL surface assembles a partition-dependent load (the force *function* integrates identically to 1e-16, but the assembled vector differs 0.027% because interior facets at partition seams lack their neighbour support cell on the non-overlapped assembly DM). Marked with TODO(BUG) at the natural-BC assembly; filed in the planning file. Universal (affects Nitsche too), amplified here by the augmented conditioning. Item 3 — interior-multiplier knockout audit. Confirmed the knockout is genuinely lossless (a converged solve with it disabled moves the answer ~1e-8 in velocity, no interior-h blow-up); the previously observed "2-5% when off" is an iterative-conditioning artifact, not lost physics. Rewrote the docstrings accordingly. Added a guard that warns a monolithic direct solve (pc_type lu/cholesky) of the constrained saddle point is a serial diagnostic only (wrong response + segfaults in parallel). Underworld development team with AI support from Claude Code * fix(stokes): project rigid-rotation gauge out of constrained free-slip solution For a free-slip (non-Dirichlet) velocity the rigid-body rotations are a true nullspace of the velocity block (A_uu·rotation = 0). The monolithic nullspace attached for the solve is not removed inside the fieldsplit/Schur iteration: the inner velocity KSP has no rotation nullspace, so it leaves an unconstrained rigid rotation in the velocity. The operator is blind to it (true residual still converges to machine precision), but the rotation amplitude is partition- dependent — the tangential velocity differed serial-vs-parallel by ~0.1–1% even at a converged residual, while the physical (radial / gauge-invariant) flow was already partition-clean. SNES_Stokes_SaddlePt.solve() now projects the velocity rotation modes out of the converged solution via MatNullSpaceRemove (_remove_velocity_rotation_gauge). Because the rotation is a genuine nullspace, removing it does not change the residual and yields the same rotation-free solution on any decomposition. This is the velocity analogue of the constant-pressure gauge (set_pressure_gauge). Surgical: only the velocity rotation modes (zero in the pressure/multiplier blocks), so the existing pressure/multiplier gauge handling is untouched; a no-op when there are no rotation modes (e.g. Dirichlet velocity BCs). Validated on the SphericalShellInternalBoundary constrained free-slip problem: the tangential surface velocity serial-vs-np8 relative difference drops from 1.2% to 2.5e-14, with the iteration count unchanged. Regression clean: test_1064 + test_1010 (9 passed, 2 pre-existing xfails) and test_1063 parallel (np=2) all pass. Also refines the TODO(BUG) at the natural-BC assembly loop to document the underlying PETSc interior-facet support[0] partition-dependence and the workarounds that were tested and ruled out. Underworld development team with AI support from Claude Code * test(stokes): guard the free-slip rotation-gauge fix Adds a fast serial tier-a/level-1 test that solves a fully free-slip annulus (both boundaries constraint BCs, no Dirichlet velocity BC -> rigid-rotation nullspace active) driven by a purely radial body force, and asserts the converged velocity carries no rigid-rotation component. The existing test_1063 cases all pin the rotation nullspace with a Dirichlet BC, so they did not cover this path; this test fails without _remove_velocity_rotation_gauge (the rotation coefficient is ~1e-3) and passes with it (~4e-15). Underworld development team with AI support from Claude Code * chore(stokes): invalidate rotation-nullspace cache + document in-solve finding Two clean-ups to the free-slip rotation-gauge handling: 1. Fix a latent staleness bug: the cached velocity rotation nullspace (_velocity_rotation_nullspace) is built from node coordinates but was never invalidated, so it would go stale after a solver rebuild / mesh deformation. It is now initialised in __init__ and reset in _reset_stokes_nullspace and the _build teardown path, alongside the monolithic Stokes nullspace. 2. Record, in _remove_velocity_rotation_gauge, the empirical result of trying the in-solve alternative: a rotations-only DMSetNullSpaceConstructor on the velocity field DOES attach the rotation nullspace to the fieldsplit velocity sub-block (verified), but the converged GLOBAL velocity still retained the rotation gauge (coefficient ~0.18, unchanged) — because removing a nullspace gauge is a projection on the converged global solution, and attaching it to a sub-block operator only constrains the inner Krylov solves (fgmres + a variable fieldsplit/Schur preconditioner reintroduces the rotation). So the explicit post-solve projection is the correct and necessary mechanism, not a stopgap. No behaviour change. Underworld development team with AI support from Claude Code
The Nitsche free-slip/constraint penalty term gamma*mu/h previously used a
single GLOBAL scalar h = mesh.get_min_radius() (the smallest cell anywhere).
The viscosity mu was already local; only h was wrong. On a non-uniform,
adaptive, or deforming mesh this mis-scales the boundary penalty: the global
minimum cell size is applied on every facet (over-stiffening coarser boundary
cells), and it drifts as refinement or distortion changes the global minimum.
Replace it with a LOCAL, per-cell size via a new deformation-tracking field
Mesh.cell_size():
- a cell-constant (degree-0, discontinuous) scalar MeshVariable holding each
cell's characteristic length (self._radii), so the boundary kernel sees the
adjacent cell's size;
- refreshed on mesh.deform() and mesh.adapt() (and via a reinit callback in
the remesh transaction), mirroring the boundary_normal machinery, so it is
never stale on a moved/re-refined mesh;
- the fill is purely rank-local (no collective), so it is parallel-safe.
Wired into add_nitsche_bc on both SNES_Stokes_SaddlePt and SNES_Vector behind
local_h=True (default). local_h=False restores the exact legacy global-h path.
On a uniform mesh local == global, so existing tuned behaviour is preserved.
Validation:
- Existing Nitsche suite green serial + np2 (test_1060/0056/1011/1014/1061;
the only np2 failure, test_1061, segfaults on unmodified development too —
its LU/mumps block-constrained path uses no Nitsche).
- New test_1065_nitsche_local_h (serial + np2): the penalty size is local
(~6x the global minimum at a coarse free-slip boundary), tracks deformation,
and still solves free-slip correctly.
- Held-lid adaptive free-surface convection (fs_convection_goal4, res 24,
rho_g 7.5e5, gamma 10): global-h spikes vhmax to 2873 (2 glitches, spurious
h_inf mountain 2.7e-2); local-h holds vhmax at 342 (0 glitches, h_inf 7.5e-3).
Medians identical => same cost, spikes removed.
Underworld development team with AI support from Claude Code
Contributor
Author
|
@lmoresi — notes for review:
|
lmoresi
approved these changes
Jul 7, 2026
lmoresi
left a comment
lmoresi
left a comment
Member
There was a problem hiding this comment.
Thanks @jcgraciosa - looks good to me.
Let's cherry pick that one IS fix back onto development so it is not re-invented later.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
5 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Underworld3 v3.1.0
Highlights
solver stack: write models in physical units (Pa·s, km, K) and let the
framework non-dimensionalise and re-dimensionalise automatically.
mesh mover are now formally supported, providing the foundation for
anisotropic adaptive remeshing workflows.
deformation capability gate, semi-Lagrangian improvements, projection
enhancements, and Darcy sign-convention corrections.
Supported (validated)
Preview (present, unguaranteed)
preconditionerproperty).Bug Fixes
-uw_*parameter overrides now apply regardless of platform (fix(params): apply -uw_* CLI overrides regardless of platform (#111) #280)New Features
mesh.deform()API with safecache invalidation, replacing the internal
mesh._deform_mesh()(Invalidate evaluate/DMInterp/topology caches on Mesh._deform_mesh #188)Improvements
Breaking Changes
mesh._deform_mesh()now raises if the mesh carries live variables — usethe public
mesh.deform()instead (Invalidate evaluate/DMInterp/topology caches on Mesh._deform_mesh #188)Contributors
Thanks to everyone who contributed to this release:
Louis Moresi, Thyagarajulu Gollapalli, Ben Knight, Neng Lu, Julian Giordani, Juan Carlos Graciosa, Saurabh Shukla, @tiannh7.