Geometric aggregation for placed-hierarchy FMG: coarse levels as decimated point clouds, not band-only meshes

[lmoresi]

The #622 closing measurement located the placed-hierarchy FMG gap in the levels themselves: band-only levels barely coarsen in DOF count (the background rides every level; a cycle smooths ~175% of the fine level's nonzeros), and scattered transfers between non-nested clouds fatten the Galerkin product (265-481 nnz/row against 90 assembled) even at 4 entries per transfer row. GAMG wins by SELECTION: aggressive, well-separated coarse points keep both the level sizes and the product sparsity down.

Proposal — do the same job geometrically:

1. Build each coarse level's point set by DECIMATION of the fine level's node cloud: a maximal independent set at radius alpha*h_local (alpha ~ 2-3), kd-tree greedy — the geometric analogue of aggregation, and it coarsens the background and the band TOGETHER (true ~8x DOF reduction per level, tunable by alpha).
2. Delaunay-triangulate each decimated cloud into a scaffold mesh so it carries a DM space — the existing CustomMGHierarchy machinery (BC-per-level reduction, coordinate-based transfers) then consumes it unchanged. The mesh is a scaffold for the section; the operator is Galerkin from the fine level regardless.
3. Transfers by the standard sparse local RBF (custom_mg.rbf_prolongation: no conditioning path on an indefinite saddle solve, and its only test asserts convergence rather than accuracy #429; nnn small — the MIS guarantees every fine point has a nearby coarse point, and the zero-column repair covers stragglers) or barycentric.

This keeps everything the placed-ribbon workflow wants — no constraint on how the finest mesh was made, split-at-finest-only per the Galerkin-degradation pattern — while giving the hierarchy GAMG-like level economics. Expected first test: the split-fault contact benchmark of ~/+Simulations/place_route_health/ (GAMG 18.6 s warm; band-only custom-P 77-97 s at BETTER iteration counts), then the damage-contrast regime where GAMG degrades and the chase actually pays (#576/#579's 13x).

Underworld development team with AI support from Claude Code

[lmoresi]

Immediate supporting measurement (the fewer-levels half of the program, same benchmark): dropping the band-only intermediate levels closes the whole gap — full tail 97.0 s warm (vel 6 its), no-mid 30.7 s (18), TWO-LEVEL (coarsest 1,689 DOFs only) 19.6 s against GAMG's 18.6 s, at 22 vs 46 velocity iterations. The intermediates were negative-value at this scale: each bought a few iterations and charged a near-full extra smoothing pass per cycle. The decimated-cloud middle level this issue proposes is the piece that should keep two-level cycle cost while restoring the low iteration count — the headroom that matters in the contrast regime where GAMG's iterations degrade.

Underworld development team with AI support from Claude Code

[lmoresi]

Prototype measured (place_route_health/ribbon_cloud.py). The decimated-cloud level WORKS: MIS decimation of the split fine mesh's vertex cloud (alpha=2.2, hull-seeded per wall), scipy Delaunay scaffold (slivers filtered; element quality is irrelevant — the level's operator is Galerkin PtAP, only point positions and the DOF section enter), box walls labelled geometrically for the BC-per-level reduction. As a TWO-LEVEL tail (cloud + split finest): converged, velocity 8 iterations against 22 for the base-only two-level and 6 for the full band hierarchy — the selection restores the iteration count, as proposed. Warm 27.6 s at this toy scale (base-only 19.6, GAMG 18.6): the iteration-vs-cycle-cost dial is now fully in hand, and the 8-iteration configuration is the one that matters for the contrast regime where GAMG's counts blow up.

Open defect, precisely located: the THREE-level chain [base L0, cloud, fine] DIVERGES (KSP dies instantly, zero block iterations) while both two-level chains converge — the L0 -> cloud transfer is the culprit, the one non-nested pair that needed 39 zero-column repairs. Candidates: the nearest-fine-DOF injection producing an inconsistent Galerkin coarse operator at that density of repairs (39 of 6,423 — the repair was designed for a handful), or the pair simply coarsening too unevenly (base L0 is uniform 0.35-spacing; the cloud is graded 0.03-0.26). Next probes: decimate a SECOND cloud from the first (nested-ish gradation, few repairs) instead of reusing base L0; or widen the transfer stencil (local RBF, larger nnn) for that pair only; or verify the repair rows on that transfer directly.

Underworld development team with AI support from Claude Code

[lmoresi]

Design reframe (Louis, 2026-08-21): the localised-refinement hierarchy is a COMPOSITE / OVERLAY grid — levels differ only in the refined patch, the background is bit-identical between levels and its nodes fall through (the coordinate transfers already produce identity rows there, observed). The measured pathology (a cycle smooths ~175% of the fine level's nonzeros; band-only levels condemned as negative value) is not the levels — it is WHOLE-LEVEL SMOOTHING on overlay grids. The classical answer is FAC/MLAT: smooth each overlay level ONLY on its refined patch plus an interface halo; the background is smoothed once, at the level that owns it. Patch sizes halve per level, so the tail costs a geometric series (~2x the finest patch) instead of levels-times-whole-mesh — the band-only hierarchy becomes near-optimal under patch smoothing.

Implementation path in the existing machinery: a PCMG level smoother is any KSP/PC, so mask it to the patch with ASM/GASM over one user-supplied subdomain IS. The active set is recoverable two ways, both in hand: the adapt engine records marked cells per generation (child._adapt_markers), and geometrically a node whose transfer row is identity (coincident coordinates) is background — one kd-tree pass.

Revised program shape: (1) exact nested transfers for the base ladder (the P2 node clouds do not nest even where meshes do — L1's Galerkin operator measured 481 nnz/row where a native coarse operator runs ~90-150; use the recorded/native any-degree nested prolongations for dm_hierarchy pairs); (2) FAC patch smoothing for overlay levels; (3) decimated clouds only where genuinely non-nested placed structure needs bridging; (4) split at the finest only, unsplit below (validated). Items (1) and (2) subsume most of what the cloud chain was compensating for.

Underworld development team with AI support from Claude Code

[lmoresi]

Items 1 and 2 landed and measured (ec21348 on feature/fault-outcrop-3d-cap; benchmark = the split-contact ribbon of ~/+Simulations/place_route_health, warm repeats, same machine-state controls).

Item 1 — the fat Galerkin chain was mostly STRUCTURAL junk, not scattered geometry. The point-located builders emit ~1e-16 weights that are structural nonzeros, and PtAP fills by structure, compounding level over level: that is the 90 -> 265 -> 481 nnz/row cascade. Dropping numerically-zero weights collapses the full-tail chain to 138/119/173/90 (the item-1 acceptance band) and the warm contact solve from 97 s to 52 s at bit-identical answers (slip 0.1509, leak 0). On top of that, native refine() pairs now get the EXACT any-degree nested prolongation (parent cells recovered topologically, weights by the dual-basis identity under the parent's affine pullback; verified exact to 1e-15 at P1/P2, 2-D and 3-D). One measured warning for the record: PETSc's DMCreateInterpolation general path is NOT the FE embedding (row sums up to 1.375, quadratic reproduction error 1e-2) — we build our own. Full-tail iterations reproduce the recorded 6-7 under the robust smoother; the same-day dial is GAMG 27.4 s / 56 its, two-level 24.7 / 22, no-mid 27.8 / 18, full tail 51.9 / 12 — the full-tail premium fell from 5x to 1.9x of the two-level cost.

Item 2 — FAC patch smoothing is implemented and correct, with three load-bearing configuration findings: (1) BASIC ASM, never restricted — discarding the halo correction stalls the outer KSP at 80-375 iterations where basic runs 6 against a whole-level baseline of 4; (2) SOR subdomain solves, never the ASM default ILU, which takes a numeric zero pivot on the rotated Galerkin patch block (min diag 8e-5 near the constraint; PC_FAILED before the first iteration); (3) the rotated path must force the FAC levels' smoother setup before its options-DB cleanup, because PCASM creates its sub-KSP lazily at first apply. The patch is read off each level's own transfer (identity row + bit-coincident node = background, decided per node so the per-node rotation cannot disturb it; slit-duplicated nodes count as patch), so placed levels work with no adapt bookkeeping.

The honest FAC verdict at this toy's proportions: a net cost. The ribbon band is 66-82% of every refined level's DOFs (the toy is deliberately fat — width 0.06 in a unit box), so there is no background to save and the ASM gather/scatter charges ~15% (60.3 vs 51.9 s; forcing the finest level in as well: 68.1 s, iterations unchanged). The geometric-series premise needs patches that SHRINK relative to their levels — production proportions. That is exactly where the item-3 contrast test goes next, with the A/B dial (UW_FAC_* / UW_CUSTOM_MG_* knobs, TODO-marked) already in place.

Underworld development team with AI support from Claude Code

[lmoresi]

Item 3 first measurements (contrast sweep on the split-contact ribbon; weak zone = the cells supporting the slit, eta = 1/contrast as a per-cell P0 field; genuine solves forced by perturbing eta between solves).

The predicted iteration headroom is REAL and holds through two decades: velocity block, GAMG vs full-tail custom-P FMG (robust gmres smoother) — contrast 1: 56 vs 7 iterations; 1e2: 63 vs 10; 1e4: 178 vs 36. GAMG's count roughly triples per two decades of contrast; FMG's stays within ~5x of its flat-viscosity count. The smoother variant matters exactly as the recorded mechanism says it should: the fast (richardson) smoother degrades to 23 iterations already at 1e2 where gmres holds 10 — the contrast regime needs the Krylov smoother.

But the sweep is BLOCKED at 1e4 by the pressure block, not the velocity block: both arms' pressure sub-solve sits at its 200-iteration gasm cap every outer application — silently unconverged (the repeating 200 is a cap, not convergence; #625's regime). The 1/mu pressure-mass Schur preconditioner spans the full 1e4 coefficient range once the viscosity is zonal, and beyond ~1e2 the wall clock measures cap-priced pressure applications rather than the velocity preconditioner: GAMG 449 s vs FMG 3487 s at 1e4, where FMG is paying its per-V-cycle cost times those applications. 1e6 was not run — it would measure the cap.

So the order of battle for converting the iteration headroom into wall time is: (a) Schur/pressure preconditioning under zonal contrast — a viscosity-weighted pressure mass or a better Schur preconditioner, the #625 follow-on; then (b) per-V-cycle cost, which is the FAC lever and needs thin-patch proportions to bite. Answers agree between arms at every rung (peak slip 0.1859 at 1e2, 0.680 at 1e4, leak 0 throughout).

Underworld development team with AI support from Claude Code

[lmoresi]

Contrast diagnosis closed (836bc47): the FMG's V-cycle count was NOT constant under contrast — 7 -> 10 -> 36 velocity iterations over 1e0 -> 1e4 where constancy is the expectation — and the mechanism is now pinned by two discriminating measurements. Level ablation at 1e4: full / no_mid / two-level tails converge in 20 / 19 / 24 iterations — tail depth is irrelevant, so the coarse corrections contribute nothing to the weak-band modes; even the band-carrying mid level buys nothing, which implicates the coefficient-blind transfers and coarse spaces rather than the level geometry. Doubling the smoother work (gmres/8 for /4) exactly halves the count (20 -> 10): convergence is proportional to total fine-level smoothing, the signature of smoother-limited band modes.

The fix falls out of the FAC machinery already in place: make the patch smoother STRONG. With a direct subdomain solve on the FAC patch (which contains the weak zone) the count collapses to 2 iterations at contrast 1e4 and 3 at 1e6 — constant across six decades and below the flat-viscosity count of 7. One prerequisite worth recording: the factorization must be shifted (sub_pc_factor_shift_type nonzero) — the rotated Galerkin patch block carries near-zero pivots inherited from the constraint-zeroed transfer rows, and an unshifted factorization takes NUMERIC_ZEROPIVOT even as exact LU.

So the composite picture: FAC patch smoothing was a wash as a COST device at toy proportions, but as a ROBUSTNESS device it is the answer to the contrast regime — the patch is where the hard physics lives, and solving it directly each cycle is what restores the multigrid contract. Remaining engineering to make the cycle cheap: shrink the subdomain from the whole refined patch to the weak zone plus halo (a much smaller factorization), reuse factors across cycles, and the #625 pressure-Schur fix, which is still the wall-clock gate for every arm at contrast.

Underworld development team with AI support from Claude Code

[lmoresi]

Switched the analysis to 2-D (the 3-D slit-subdomain experiment was not converging in useful time; 3-D resumes when the skeleton is settled). The rig is line_mg.py in the place_route_health campaign: 2-D box, interior line cut, split fault, frictionless contact, 2.5k cells — under a second per warm solve, so the design space is explorable interactively.

The clean result: with ordinary multigrid for the smooth problem and one strong subdomain solve on the fault patch — basic ASM, subdomain = the cut-affected rows (508 of 9,848, 5%), shifted LU — the velocity block converges in 6 / 6 / 7 / 7 iterations at contrast 1 / 1e2 / 1e4 / 1e6. Constant across six decades. The controls on the same rig: GAMG degrades 45 / 29 / 68 / 200 (its cap); plain FMG without the strong patch drifts 8 / 6 / 16 / 50. This confirms the design separation on proportions that actually resemble the target: the smooth refined bulk needs nothing special, and the strong solve needs only a subdomain proportional to the fault trace — keyed on the physics (split-duplicated nodes plus an operator-layer halo, UW_FAC_PATCH=slit, 910c636), never on the refinement structure, which in a global subduction mesh is most of the node budget.

The known blocker is unchanged and now cleanly isolated: the pressure gasm sub-solve caps at 200 iterations from contrast 1e4 (all the wall-time inflation from 1.2 to 7.3 s is pressure), so the #625 Schur fix under zonal contrast is the next structural piece. After that, the rig's next rungs are multi-segment faults with junction merging, and an artificially-segmented straight fault to test whether the along-fault coarse space is needed at realistic segment counts.

Underworld development team with AI support from Claude Code

[lmoresi]

Two more design questions answered on the 2-D rig (49d01f0), both in the direction that simplifies the California-scale architecture.

Segmentation is free. The fault trace's strong-solve subdomain split into k along-strike LU blocks (one ASM subdomain each, overlap 1): 7 velocity iterations at k = 1, 2, 4, 8, 16 at contrast 1e4, and still 7 at k = 16, contrast 1e6. Block-wise segment solving costs nothing at any segment count we can test on this rig. The mechanism reading: the stiff coupling is ACROSS the fault and block-local, while the along-fault long-range mode is smooth along the trace, so the ordinary coarse ladder carries it despite coefficient-blind transfers. The fault-network coarse space stays in reserve rather than becoming a prerequisite. (One instrumentation lesson for the record: the first k-sweep was flat at plain-FMG counts because the slit detection keyed on transfer identity rows onto coarse nodes — empty for a cut mesh, whose slit vertices are new points at edge crossings — and the mode silently declined. Detection is now coincident fine pairs, the split's plus/minus nodes at identical coordinates. Check the probe, not just the iteration count.)

The strong lid is benign. The subduction-shaped fixture — a cold strong layer (eta 1e4, as a sharp jump or a tanh profile) crossed mid-trace by the 1e-4 gouge, i.e. 1e8 across the fault walls inside the lid: slit-FAC holds 7 iterations for BOTH lid profiles at 7.2 s, with the lid boundary lying entirely outside the patch; GAMG runs 88 (jump) / 80 (smooth); patchless FMG 16/16. So the extended, co-dimension-1 lid jump is the benign kind of discontinuity for the Galerkin ladder plus a Krylov smoother, and the malignant kind — the thin cross-cutting weak inclusion — is exactly what the patch already owns. GAMG's jump-vs-smooth separation is mild here because the gouge, present in every arm, dominates its degradation; the discontinuity-not-contrast effect of the old SolCx control is being carried by the fault itself.

Cumulative picture on the rig: ordinary multigrid + a fault-trace strong solve (in as many blocks as convenient) is constant at 6-8 velocity iterations against every structural aggravation tried so far — contrast to 1e6, sixteen-way segmentation, a strong lid of either profile — and is also the fastest arm in wall clock at contrast >= 1e4 despite the still-unfixed pressure cap (#625), because fewer outer applications mean fewer cap-priced pressure solves. Next: the junction case (two segments meeting at a point), then the #625 Schur fix.

Underworld development team with AI support from Claude Code

[lmoresi]

Clarification measured (d4b86ea): the strong-patch mechanism is not tied to split-node faults. The patch solve — ASM blocks with shifted factorizations over a row set — never sees the constitutive model; only the automatic detection (coincident split pairs) was split-specific. For painted fault models detection is unnecessary: the modeler knows the painted cells, and a boolean cell mask on the solver (solver._fac_zone_cells) keys the finest level's patch to them directly.

Measured on the 2-D rig with an UNSPLIT painted weak band (the weak-fault model, standard Stokes path, no contact machinery): the zone-patched FMG holds 9 velocity iterations at contrast 1e4 and 7 at 1e6, against GAMG's 40-and-degrading — the same contrast-independence the split-contact configuration shows, with no split anywhere in the problem. The subdomain solver barely matters on the unsplit band at this size (SOR and LU blocks both give 9): the split-contact patch is the more singular one and is where the shifted LU earns its keep.

Consequence: the choice among split-node, painted-weak, and TI fault representations returns to being a physics decision — branch behaviour, gouge-width resolution, whether topological slip is required — with the solver architecture indifferent. TI is covered by the same mechanism by construction (the block factorization is exactly the treatment for the band's stiff anisotropic coupling, which is what defeats point smoothers); a painted-TI confirmation run on the rig is the remaining extension.

The segmentation ladder also completes the scaling statement: the patch reduces the "direct solve" from N volume DOFs to the trace m ~ N^(2/3), and since iterations are flat in the number of along-strike blocks, fixed-size segment blocks make the total factorization cost linear in the trace and embarrassingly parallel — no superlinear direct-solver scaling survives anywhere in the design.

Underworld development team with AI support from Claude Code