From 7551e12ef6b0d142ea279d6557b3e81d7ff46225 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 1 Aug 2026 08:53:21 +1000 Subject: [PATCH] skills: parallel adapt engines, band pinning, and FMG on adapt children adapt-on-top-faults - engines section: nvb vs edge_split, both parallel in 2-D and 3-D and bit-confluent; edge_split has no conforming closure so refinement cannot escape the marked region, and marks on the DIAMETER (the volume proxy reported the target met while the mesh was 3.2x coarser across the feature). - repair=True: gates on reducing the largest angle, NOT on Delaunay. Delaunay maximises the minimum angle while P1 depends on the maximum, and flipping a gmsh mesh toward Delaunay raised the 99th-percentile max angle. Worth it on a poor base (156 -> 115 degrees, slivers 3.84% -> 0.00%), marginal on a clean one, and it gives up bit-confluence, so it is opt-in. - relax on a mesh refined onto an interface makes things WORSE (+77% leak); pin_bands is the fix. - new section on sizing the band and representing the fault margin: the -2 Cov(eta, edot) leak metric, why a within-cell marking rule loses to the plain distance size field, why the optimal band width depends on which objective you pick, and what a step-edged margin buys and costs. - gotchas: Mesh(dm) takes the DM over (bare SIGSEGV if you keep using the old handle); Mesh(dm) without boundaries= loses the boundary enum; evaluate() "Total components 8 != 6" on a variable-heavy mesh. adaptive-meshing - PIN THE INTERFACE section for relax(pin_bands=...), including the signed-vs-unsigned distance rule and the pinned_labels merge trap. - cross-reference to nonlinear-solver for the FMG setup. nonlinear-solver - new section: FMG on an adapt-on-top child. The child carries its own graded custom-P tail and solvers pick it up automatically; the base must have refinement>=1; a base-only tail triples the V-cycle count; V-cycle counts are insensitive to element quality (a pass, not a failed measurement) so use GAMG as the quality probe; relax can trip #424 into the dense RBF fallback; repair invalidates the any-degree transfer but not the vertex prolongation. - cross-references to adapt-on-top-faults and adaptive-meshing. Underworld development team with AI support from Claude Code --- .claude/skills/adapt-on-top-faults/SKILL.md | 123 +++++++++++++++++++- .claude/skills/adaptive-meshing/SKILL.md | 36 +++++- .claude/skills/nonlinear-solver/SKILL.md | 47 +++++++- 3 files changed, 202 insertions(+), 4 deletions(-) diff --git a/.claude/skills/adapt-on-top-faults/SKILL.md b/.claude/skills/adapt-on-top-faults/SKILL.md index 207bf3bd2..a9c90d324 100644 --- a/.claude/skills/adapt-on-top-faults/SKILL.md +++ b/.claude/skills/adapt-on-top-faults/SKILL.md @@ -75,6 +75,77 @@ director), `fault.refinement_metric_function(...)`. --- +## Engines: `nvb` vs `edge_split` — and what runs in parallel + +`engine="nvb"` (default) is graded newest-vertex bisection: bounded conforming +closure, a similarity-class bound that keeps child quality tied to the base, and +**partition-independent** output. `engine="edge_split"` splits the **longest edge** +of every cell coarser than the metric asks for and needs **no conforming closure +at all**, because splitting an edge divides every incident cell at the same new +vertex. Consequences: + +- refinement **cannot escape the marked region** — the band hugs the feature + instead of a halo around it; +- it marks on the cell **DIAMETER**, not `(dim!·vol)^(1/dim)`. The volume proxy + reported the target met while the mesh was **3.2× coarser** across the feature; +- it gives up the similarity-class bound, so quality at depth is not guaranteed + the way bisection's is — that is what `repair=` and `relax()` are for. + +**Both run in parallel, 2-D and 3-D, and both are bit-confluent** (identical mesh +at any communicator size). `edge_split` drives the same compiled `uwnvb_bisect` +transform as NVB, so it inherits star-forest propagation, co-partitioning, labels +and coordinates. Verified at np=1/2/3/4 up to 56k cells. + +```python +child = base.adapt(metric, max_levels=3, engine="edge_split") +child = base.adapt(metric, max_levels=3, engine="edge_split", repair=True) +``` + +`repair=True` runs a **reconnection (Lawson flip) pass** after each generation — +2-D and `edge_split` only; it raises rather than silently doing nothing otherwise. +It gates on **reducing the largest angle**, NOT on Delaunay: Delaunay maximises the +*minimum* angle while P1 interpolation depends on the *maximum* (Babuška–Aziz), and +flipping a gmsh mesh toward Delaunay was measured to RAISE the 99th-percentile max +angle 126.8° → 129.3°. gmsh optimises shape, not the empty-circle property. + +- **worth it on a POOR base** — anisotropic, graded, relaxed, or read from a file: + 99th-pct max angle 156° → 115°, slivers below q=0.1 3.84 % → 0.00 %. On a clean + gmsh base it moves 124.7° → 120.5° and the error not at all. +- ⚠️ **it gives up bit-confluence.** Which cavities may be flipped depends on where + the partitioner cut (no cavity may contain a cell incident on a shared point). + Conformity, orientation, volume, labels and the SF stay exact at every rank + count; only the choice of flips near a seam differs. Hence opt-in. +- Seam cost is small and **shrinks with resolution**: frozen repair sites 0.9–3.5 % + at 56k cells, np=2..8, halving with every halving of the target size. In a fault + band specifically, 5.5 % at np=2 and 13 % at np=4 on a 4k-cell mesh. +- ⚠️ the 99th-pct angle recovers under a frozen seam but the **absolute max does + not** — a few worst cells sit on the seam (148° vs 123° serial). +- it invalidates the cell-parent map for the any-degree MG transfer (a flipped + cell can straddle two coarse cells), so degree ≥ 2 falls back to the geometric + prolongation builder. The exact vertex prolongation survives — flips move no + vertex. + +## Relaxing an adapted fault mesh — PIN THE BAND + +`child.relax()` on a mesh refined onto an interface **makes things worse**. The +MMPDE mover optimises element shape against an equilateral reference and knows +nothing about where the material changes, so it slides the small cells that +refinement placed on the interface *off* it. Measured on a step-edged fault: +manufactured stress across the interface **+77 %**, and it stopped being confined +to the fault. Counter-intuitively it *reduces* the number of straddling cells +(1343 → 965) and is still worse, because the survivors are bigger. + +```python +child.relax(pin_bands=[fault]) # interface = the surface itself +child.relax(pin_bands=[(fault, 0.02)], pin_halo=2) # weak zone of half-width 0.02 +``` + +Leak unchanged to five decimal places (0.03075 → 0.03076), confinement preserved, +straddling count identical — while the mover still reshapes the rest of the +domain. `pin_halo` (default 1) pins extra rings; pinning only the cut cells lets +the mover pull on them from outside. `pin_bands` **merges** with the auto-pinned +boundaries, so it cannot silently release the domain edge. + ## Fault as a constitutive weak zone (iso and TI) The metric only needs the fault GEOMETRY (pure distance). The constitutive weak zone @@ -219,6 +290,51 @@ for step in range(nsteps): --- +## Sizing the band, and how the fault margin is represented + +The artefact that matters for a fault is **stress manufactured by elements that +straddle the weak-zone margin** — high strain rate at one end, high viscosity at +the other. It is exactly + +```python +leak = 2 * (eta.mean(axis=1) * edot.mean(axis=1) - (eta * edot).mean(axis=1)) +``` + +per cell (vertex values), i.e. `−2 Cov(η, ε̇)`: **zero** for any cell wholly inside +or wholly outside the weak zone, positive only across the transition. It lives +strictly *inside* elements — plotting nodal `2ηε̇` cannot show it, because at a node +the two fields are sampled at the same point and are consistent by construction. + +Measured guidance, all at matched cell count: + +- **Band width.** The answer depends on what you are minimising, and the two + objectives disagree. *Total* leak: narrower is better (concentrate cells where + ∇η is steepest). Leak **into the matrix** (what usually matters): an optimum at + core half-width ≈ the **influence width**, 2.6× better than a narrow band. + Straddling-cell count: wider is monotonically better. +- **Don't invent a marking rule.** Marking on within-cell η variation is + intuitive and measurably *worse* per DOF than the plain distance size field: + N^-0.37 (absolute jump) or a complete stall (log ratio) against **N^-1.04**. The + leak is spread across the whole transition, not concentrated in a few cells, so + there is nothing for a targeting rule to target. ⚠️ The log ratio is largest + where η is *smallest* — it refines the fault core, the opposite end from the + problem. +- **A step-edged margin confines it.** `influence_function(profile="step")` (the + DEFAULT profile) plus marking on the distance level set puts essentially **0 %** + of the leak beyond d=0.03, against 11.4 % for a smooth blend, and converges + slightly faster (N^-1.32). The price: total leak 2.5× higher and the **worst + single cell 20× worse** (21.4 vs 1.04) — concentrated into a one-cell collar + welded to the interface rather than spread. For a viscous solve that is a clear + win; for a yielding model the worst cell is what reaches yield first, so weigh + it. Mark geometrically on the level set: once the edge is sharp, sampled η + depends on which side a vertex happens to fall. +- **Exact fixes.** An element-wise constant (P0) viscosity makes `Cov(η, ε̇) ≡ 0` + on any mesh — not reduced, zero. So does aligning the interface with element + boundaries. Both move the error from *inside* elements to *where the element + boundaries fall*, which makes `relax(pin_bands=...)` the lever rather than shape + repair. ⚠️ P0 also breaks any within-cell marking rule (contrast is identically + zero) — it would have to be reposed on the facet jump. + ## Gotchas / rough edges (candidates to fix as we go) | symptom / edge | cause & handling | @@ -231,7 +347,12 @@ for step in range(nsteps): | rotated free-slip Stokes "not converged" | `s.snes` isn't the solving object (manual loop). Read `stokes._rotated_freeslip_info['ksp_reason']` / `['nonlinear_iterations']` (PR #298); sanity-check v·n leakage. | | nonlinear TI/VEP + rotated free-slip + timestepping gives a wrong (frozen) answer | pre-#298 the rotated path is ONE linear solve (one Newton step from u=0). PR #298 runs it inside a Newton/Picard loop — ensure it's merged for nonlinear/warm-start runs. | | FMG not used under rotated free-slip | rotated_bc uses a self-contained KSP and reads `solver._custom_mg`, not the native `dm_hierarchy` — `set_custom_fmg(..., field_id=0)`. FUNDAMENTAL (DM-less rotated operator can't use the DM-coupled fieldsplit; #298 keeps it), not a quick fix. | -| NVB at np>1 raises NotImplementedError | native `_nvb_transform` extension not built (needs the custom-PETSc/amr env). | +| NVB at np>1 raises NotImplementedError | native `_nvb_transform` extension not built (needs the custom-PETSc/amr env). Both `nvb` and `edge_split` are otherwise fully parallel, 2-D and 3-D. | +| high stress appears in the matrix beside the fault | elements STRADDLING the weak-zone margin: one end sees high strain rate, the other high viscosity. The FE forms `mean(η)·mean(ε̇)`; the honest cell average is `mean(η ε̇)`, and the difference is `−2 Cov(η, ε̇)` across the cell. Zero for any cell wholly in or wholly out. See the band-width section below. | +| refinement band narrower than the fault's INFLUENCE | measured: η still 0.07 at d=0.06 while the mesh has already coarsened 4×, so the artefact peaks on the transition flank, not on the fault. **85 % of it sits at d>0.01.** Size the flat core from the *influence* width, not the fault. | +| `uw.function.evaluate` fails "Total components 8 != 6" | cached-interpolation mismatch on a mesh already carrying several solver variables. Sample the field numerically from `surface.unsigned_distance` instead. | +| bare SIGSEGV, no traceback, after building a Mesh from a raw DM | `uw.discretisation.Mesh(dm, ...)` TAKES THE DM OVER. Read geometry from `child.dm`, never the handle you passed in. | +| `KeyError: 'Left'` from a Mesh built on a refined DM | `Mesh(dm)` without `boundaries=` loses the boundary ENUM even though the labels are on the DM. Pass `boundaries=base.boundaries`. | **Build/run**: this lives on the `feature/adapt-on-top` worktree; env `.pixi/envs/amr-dev/bin/{python,mpirun}`; `./uw build` after source changes. diff --git a/.claude/skills/adaptive-meshing/SKILL.md b/.claude/skills/adaptive-meshing/SKILL.md index d54697f19..80883d345 100644 --- a/.claude/skills/adaptive-meshing/SKILL.md +++ b/.claude/skills/adaptive-meshing/SKILL.md @@ -21,10 +21,42 @@ Companion: the `uw-visualisation` skill for rendering results. **Choosing the paradigm:** THIS skill is the **mover** (node movement / equidistribution, `smooth_mesh_interior`) — the mesh deforms to follow a field. For -**local refinement** instead (`mesh.adapt(engine="nvb")` returns a refined CHILD; a +**local refinement** instead (`mesh.adapt(...)` returns a refined CHILD; a fault resolved by a fine band + custom-P FMG + rotated free-slip + dynamic topography + advection-diffusion), use the **`adapt-on-top-faults`** skill. Different tools — -don't mix them. +don't mix them. For the FMG setup that consumes an adapt child's hierarchy, see the +**`nonlinear-solver`** skill. + +--- + +## PIN THE INTERFACE when you relax a mesh that was refined onto one + +The two operations fight. The mover optimises element **shape** against an +equilateral reference and knows nothing about where the material changes, so it +slides the small cells that refinement placed on an interface *off* it. Measured +on a step-edged fault: manufactured stress across the interface **+77 %**, and it +stopped being confined to the fault. It even *reduces* the number of straddling +cells (1343 → 965) while making things worse, because the survivors are bigger — +leak per straddling cell up 2.5×. + +```python +child.relax(pin_bands=[fault]) # interface = the surface +child.relax(pin_bands=[(fault, 0.02)], pin_halo=2) # weak zone, half-width 0.02 +``` + +Leak unchanged to five decimals, confinement preserved, straddling count identical +— and the mover still reshapes everywhere else. Notes: + +- `pin_halo` (default 1) pins extra rings. Pinning only the cut cells lets the + mover pull on them from outside and drag the pinned ring out of shape anyway. +- `pin_bands` **merges** with `pinned_labels`. Passing `pinned_labels` yourself + REPLACES the default of "pin every named boundary", so a hand-rolled version + that substitutes the band label silently lets the mover deform the domain. +- `mesh.label_interface_band(surface, offset, halo)` is the underlying helper if + you want the label for something else. It uses the SIGNED distance at offset 0 + and the UNSIGNED distance at a non-zero offset — the unsigned distance is never + negative, so a straddle test against it at offset 0 can never fire, and a weak + zone has two margins that the unsigned form catches at once. --- diff --git a/.claude/skills/nonlinear-solver/SKILL.md b/.claude/skills/nonlinear-solver/SKILL.md index 7bd3965a0..0051eb92c 100644 --- a/.claude/skills/nonlinear-solver/SKILL.md +++ b/.claude/skills/nonlinear-solver/SKILL.md @@ -167,6 +167,48 @@ conditioning — the coarsest grid cannot represent the viscosity contrast — a levels *every* smoother fails there (richardson outright, gmres with ρ>1). Use the δ/ξ continuation to stay in the solvable region. +## FMG on an ADAPT-ON-TOP child (locally refined meshes) + +An `adapt()` child carries its **own custom-P geometric MG tail** — one level per +refinement generation — on `child._custom_mg_coarse_meshes`, and solvers built on +it pick it up automatically. So the usual advice above ("never use a non-nested +hierarchy") is satisfied without you assembling anything: + +```python +child = base.adapt(metric, max_levels=3, engine="edge_split") +stokes = uw.systems.Stokes(child, velocityField=v, pressureField=p) +stokes.solve() # pc=mg auto-attached off the child's tail +``` + +Requirements and traps, all measured: + +- **The base must be built with `refinement>=1`.** That uniform tail is what the + adapt levels extend; without it the hierarchy starts at the adapted mesh and + there is no coarse grid. +- **Keep the GRADED tail** (one MG level per generation), which is what + `_adapt_nested` stores. Handing the solver a base-only tail instead — coarse + base straight to the fully adapted mesh — **triples the V-cycle count**. +- **V-cycle counts are insensitive to element quality here, and that is a PASS not + a failed measurement.** On a fault child the velocity block takes 2 iterations + (iso) or 2–3 (TI) across meshes ranging from 156° to 105° max angle. The + geometric hierarchy's coarse spaces come from the mesh hierarchy, not from the + fine operator, so shape does not move it — which is exactly what makes + adapt-on-top viable. **If you want a solver-side probe of mesh quality, use + GAMG**, which does respond (iso 79 → 64 velocity iterations with `repair=True`). +- **`relax()` can trip #424.** On a relaxed, unrepaired child the barycentric + transfer hit 22 zero columns and fell back to the DENSE global RBF builder — a + performance cliff, not just a warning. Watch for the `custom_mg: barycentric + transfer build failed ... retrying with 'rbf'` message. +- **`repair=True` invalidates the any-degree nested transfer** (a flipped cell can + straddle two coarse cells), so degree ≥ 2 falls back to the geometric builder. + The exact ½,½ vertex prolongation survives, because flips move no vertex. +- Under **rotated free-slip** none of this is automatic — that path builds its own + KSP and reads `solver._custom_mg`. See the `adapt-on-top-faults` skill. + +Companion skills: **`adapt-on-top-faults`** (building the child, engines, repair, +band sizing), **`adaptive-meshing`** (the mover, and `relax(pin_bands=...)` for +relaxing a mesh that was refined onto an interface). + ## Gotchas - **`./uw build` → `amr-dev` env**; verify `uw.__file__` is the worktree site-packages. @@ -183,4 +225,7 @@ continuation to stay in the solvable region. - Continuation driver: `underworld3.systems.yield_continuation`. - Diagnostics: `SNES_*.get_snes_diagnostics()` / `solve_with_diagnostics()`. - Related skills: `plasticity-solvers` (yield law + tangent per model), - `free-surface-convection`, `adaptive-meshing`. + `free-surface-convection`, `adaptive-meshing` (mover + `relax(pin_bands=...)`), + `adapt-on-top-faults` (locally refined children and their MG tail). +- Reconnection / refinement engines: + `docs/developer/design/mesh-reconnection-and-delaunay-adapt.md`.