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`.