From 7551e12ef6b0d142ea279d6557b3e81d7ff46225 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 1 Aug 2026 08:53:21 +1000 Subject: [PATCH 1/4] 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`. From 661e847fefa3e1ef4f11e9481d5d86b0262208f3 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 13 Aug 2026 11:22:10 +1000 Subject: [PATCH 2/4] skills: bring the two adapt/solver skills current with development MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit adapt-on-top-faults: the rotated free-slip path now picks up an adapt child's mesh-owned MG tail automatically (custom_mg.build_transfers, the #467 fix) — the "FUNDAMENTAL, not a quick fix" gotcha row was describing a bug that is gone. Point the band-sizing section at the interface-alignment primitives that now exist (place_sheet / place_thin_volume / remove_embedded, #517-#526). nonlinear-solver: three stale claims corrected — rotated free-slip is no longer the exception to the automatic tail pickup (#467); the adapt tail is one MG level per DOUBLING of h, not per refinement generation (mg_coarsening_ratio=2.0, #515); refinement=0 still yields a hierarchy that starts at the base, so "no coarse grid" overstated it. Three capabilities that landed since the branch: preconditioner="gamg" is respected on adapt children (#530), single-field FMG (#478/#534), and solver.pc_fallbacks as the observability hook (#534). The supersession note no longer names the removed in-SNES ramp API. Underworld development team with AI support from Claude Code --- .claude/skills/adapt-on-top-faults/SKILL.md | 30 ++++++++----- .claude/skills/nonlinear-solver/SKILL.md | 50 ++++++++++++++------- 2 files changed, 51 insertions(+), 29 deletions(-) diff --git a/.claude/skills/adapt-on-top-faults/SKILL.md b/.claude/skills/adapt-on-top-faults/SKILL.md index a9c90d324..48cd277fa 100644 --- a/.claude/skills/adapt-on-top-faults/SKILL.md +++ b/.claude/skills/adapt-on-top-faults/SKILL.md @@ -217,15 +217,19 @@ xs, sig = stokes.boundary_normal_traction("Upper") # or the raw σ_nn --- -## FMG under rotated free-slip — READ THIS (rough edge) +## FMG under rotated free-slip -`rotated_bc.solve_rotated_freeslip` builds its OWN fieldsplit KSP and uses custom-P -FMG **only if `solver._custom_mg` is set** (via `set_custom_fmg`). It does **not** -read the native `dm_hierarchy`. So: +`rotated_bc.solve_rotated_freeslip` builds its OWN fieldsplit KSP, but it resolves +the multigrid hierarchy through the same `custom_mg.build_transfers` rule as the +standard path: an explicit `set_custom_fmg` registration wins, otherwise a +mesh-owned adapt tail is picked up opportunistically. So: - On an `adapt()` **child**, the mesh-owned custom-P tail is auto-picked-up for the - Stokes velocity block → FMG for free (standard, non-rotated solves). -- On a plain **refined `Annulus`** with **rotated** free-slip, native FMG is ignored; - to get FMG you must build a coarse-mesh tail and call: + velocity block → FMG for free, **rotated free-slip included** (the old + unreachability — rotated solves silently falling back to GAMG on adapt + children — was #467, fixed). +- On a plain **refined `Annulus`** with **rotated** free-slip, the native + `dm_hierarchy` is still not read; to get FMG you must build a coarse-mesh tail + and call: ```python from underworld3.utilities.custom_mg import set_custom_fmg set_custom_fmg(stokes, [Annulus(cs=0.16), Annulus(cs=0.08)], field_id=0) # velocity block @@ -330,10 +334,12 @@ Measured guidance, all at matched cell count: 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. + boundaries — and there are now primitives that do exactly that: `place_sheet` / + `place_thin_volume` / `remove_embedded` in `utilities/place_surface.py` + (#517–#526). Both fixes 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) @@ -346,7 +352,7 @@ Measured guidance, all at matched cell count: | surface-breaking fault: huge local vmax; topo peak keeps growing | real stress singularity at the outcrop. Topo peak **saturates** (integrable/log, bounded by finite buoyancy) — far field is fine; only the pointwise outcrop value is mesh-dependent. | | 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. | +| FMG under rotated free-slip on a plain (non-adapt) refined mesh | the rotated KSP resolves hierarchies via `custom_mg.build_transfers`: an adapt child's mesh-owned tail is picked up AUTOMATICALLY (#467 fixed the old silent GAMG fallback), but the native `dm_hierarchy` is still not read — on a plain refined mesh, `set_custom_fmg(..., field_id=0)`. | | 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. | diff --git a/.claude/skills/nonlinear-solver/SKILL.md b/.claude/skills/nonlinear-solver/SKILL.md index 0051eb92c..4f915e845 100644 --- a/.claude/skills/nonlinear-solver/SKILL.md +++ b/.claude/skills/nonlinear-solver/SKILL.md @@ -70,9 +70,10 @@ proven solver config**. Mechanism: the continuation only sharpens δ from a **converged**, well-conditioned iterate; the in-SNES ramp sharpens δ **mid-solve** at a far-from-solution iterate where the consistent-Newton Jacobian on a sharpening surface is ill-conditioned and the linear solve fails. **Hide the *continuation*, -not the *ramp*.** (This supersedes the `enable_yield_homotopy()` in-SNES ramp still -described in `plasticity-solvers`; that path is retained only as a dead-experiment -record — use the multi-solve continuation above.) +not the *ramp*.** (An in-SNES ramp API once shipped and has been removed from the +source entirely — use the multi-solve continuation above; `plasticity-solvers` +carries the yield-law substrate and the evidence on when a δ-march is worth it +at all.) --- @@ -169,10 +170,11 @@ 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: +An `adapt()` child carries its **own custom-P geometric MG tail** — subsampled to +one level per **DOUBLING of h** (`mg_coarsening_ratio=2.0`, the `adapt()` default) +— 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") @@ -182,12 +184,15 @@ 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**. +- **Build the base with `refinement>=1` for a deeper tail.** The custom-P tail + always starts at the BASE mesh — with `refinement=0` it is + `[base] + the intermediate doubling levels`, so there IS a coarse grid — but + the uniform base levels extend it downward, and in MG you want the coarsest + grid as coarse as it can be. +- **Keep the GRADED tail.** `_adapt_nested` stores one MG level per doubling of + resolution (`_subsample_mg_levels`; per-bisection-pass levels were measured + 2.3–7.3× slower). 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 @@ -195,15 +200,26 @@ Requirements and traps, all measured: 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`). + That is now actionable: `solver.preconditioner = "gamg"` is **respected** on an + adapt child (#530) — before that guard the opportunistic pickup silently + clobbered it back to `pc=mg`, so any FMG-vs-GAMG comparison was vacuous. +- **Single-field solvers get FMG too** (#478/#534): `preconditioner = "fmg"` on a + Poisson/projection-class solver builds the custom-P tail over the mesh's own + `dm_hierarchy` — the section is not Stokes-or-adapt-child only. - **`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. + performance cliff, not just a warning. +- **Every PC degradation is recorded in `solver.pc_fallbacks`** (#534) — the + requested/installed/reason record for the #424 barycentric→rbf retry, a + collapsed hierarchy, a declined pickup. Read that, don't scrape warnings. - **`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. +- Under **rotated free-slip** the mesh-owned adapt tail is picked up automatically + too — the rotated KSP resolves hierarchies through the same + `custom_mg.build_transfers` rule (#467 fixed the old silent GAMG fallback). See + the `adapt-on-top-faults` skill for the plain-refined-mesh case, which still + needs `set_custom_fmg`. Companion skills: **`adapt-on-top-faults`** (building the child, engines, repair, band sizing), **`adaptive-meshing`** (the mover, and `relax(pin_bands=...)` for From 66c33a617a51f7e1bb9ad92e357438de24ea70a6 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 13 Aug 2026 11:24:13 +1000 Subject: [PATCH 3/4] skills: rewrite plasticity-solvers around what the yield campaigns measured Fixes #454. The skill taught the retired yield-homotopy doctrine in ~8 places, headlined by an enable method that no longer exists in src/ and an in-SNES delta-ramp that is separately proven to diverge. The doctrine rested on a unit-scaling error: re-measured on the corrected problem, the delta-march never succeeded where a direct hard-Min solve failed, and the ruling was to regularise the problem, not the solver. Rewritten around the evidence: Newton with the automatic Picard entry (Picard is an entry requirement, not an accelerator), rescue on failure OR stagnation, grid sequencing as the validated warm start, the per-model tangent table, and the #475 yield_mode / yield_smoother / yield_anchor substrate presented as a modelling choice with the multi-solve-only discipline for any delta march. yield_continuation is described honestly per open issue #473 (the cold-start guarantee does not hold on a Piecewise yield stress; the step control is effectively one-shot). Floors updated to the post-#475 semantics (viscosity_min_rounding). Kept: the Newton-confirmation check, the VEP-indefinite ruling, the SNESFAS ruling, the Picard-folklore footnote. Underworld development team with AI support from Claude Code --- .claude/skills/plasticity-solvers/SKILL.md | 215 ++++++++++++++------- 1 file changed, 150 insertions(+), 65 deletions(-) diff --git a/.claude/skills/plasticity-solvers/SKILL.md b/.claude/skills/plasticity-solvers/SKILL.md index 5b468ef63..d5f83bfd0 100644 --- a/.claude/skills/plasticity-solvers/SKILL.md +++ b/.claude/skills/plasticity-solvers/SKILL.md @@ -1,96 +1,178 @@ --- name: plasticity-solvers -description: How to get hard-Min viscoplastic / visco-elastic-plastic (VEP) Stokes solves to CONVERGE in Underworld3 — the consistent Newton tangent (solver.consistent_jacobian), the δ-soft-min yield law, and the yield homotopy that pairs them. Reach for THIS first when a Drucker-Prager / yield-stress Stokes solve stalls, diverges (DIVERGED_LINEAR_SOLVE / line-search fail), or grinds through ~20+ nonlinear iterations. Tells you which tangent to use per model, how to confirm you are actually running Newton, and the measured failure modes. +description: How to get hard-Min viscoplastic / visco-elastic-plastic (VEP) Stokes solves to CONVERGE in Underworld3 — Newton with the automatic Picard entry (solver.consistent_jacobian), which tangent per model, grid sequencing for the hard cases, and the δ-soft-min substrate (yield_mode / yield_smoother / yield_anchor) as a modelling choice. Reach for THIS first when a Drucker-Prager / yield-stress Stokes solve stalls, diverges (DIVERGED_LINEAR_SOLVE / line-search fail), or grinds through ~20+ nonlinear iterations. Tells you which tangent to use per model, how to confirm you are actually running Newton, and the measured failure modes. For the solver-config trap list and multigrid, see `nonlinear-solver`. --- # plasticity-solvers The workable recipe for **nonlinear convergence of yielding (viscoplastic / VEP) Stokes** in Underworld3. Hard-`Min` yield laws have a non-differentiable kink that -breaks naive solvers; this encodes the combination that converges. +breaks naive solvers; this encodes what the yield campaigns measured actually works — +and records what was retired. -**One call does it:** +**The default call is now just:** ```python stokes.constitutive_model = cm # ViscoPlastic / ViscoElasticPlastic / TI-VEP cm.Parameters.yield_stress = tau_y # finite -> plasticity active -cm.enable_yield_homotopy() # <- the strategy: δ-ramp + the right tangent -stokes.solve(zero_init_guess=False) +stokes.consistent_jacobian = True # Newton tangent (non-elastic DP; see table) +stokes.solve() # cold start takes ONE Picard step automatically ``` -`enable_yield_homotopy()` is the recommended default for any hard-Min solve. It picks -the Jacobian tangent per model (below) and ramps the soft-min δ→0 within one solve. - --- -## The two ingredients +## The doctrine (measured, 2026-07 campaigns) Yielding viscoplasticity is `η_eff = Min(η_visc, η_yield)`, `η_yield = τ_y/(2·ε̇_II)`. The `Min` kink is what makes it hard. -1. **Consistent Newton tangent** — `solver.consistent_jacobian`: - - `False` (default) — frozen-viscosity **Picard** tangent: contractive, globally - stable, **linear** (slow). Bit-identical to long-standing UW3. - - `True` — full **Newton** tangent (the assembled `dF1/dL` carries `∂η/∂(grad v)`): - **quadratic** near the solution; the kink can break it far from it. - - `"continuation"` — Picard→Newton α-blend (α a `constants[]` atom ramped 0→1). - -2. **δ-soft-min yield law + homotopy** — `g = 1 + ½(f-1+√((f-1)²+δ²)) − offset`, - `η_eff = η_ve/g`, `f = η_ve/η_pl`: - - **δ = 0 (default) ≡ exact `Min`** to machine precision. - - δ and the onset offset are `constants[]` atoms → δ is **runtime-rampable with no - JIT recompile**. `enable_yield_homotopy()` ramps δ from `delta_start`→0 within one - solve (residual-paced absolute schedule, via the `SNESSetUpdate` hook). The smooth - (δ>0) problem warm-starts the sharper one; δ ends at 0 so the **converged answer is - on the exact yield surface**. - -This is **problem-space** continuation (ramp the residual smooth→sharp). A smooth -Jacobian on a sharp `Min` residual is the consistent tangent of a *different* (harmonic) -problem and diverges worse than Picard — don't do that. +1. **Picard is an ENTRY requirement, not an accelerator.** On a cold start under + the consistent tangent, `solve()` takes one automatic Picard (frozen-tangent) + step and then runs Newton (fires only when `picard==0`, + `consistent_jacobian is True`, and the start is cold — from a warm iterate, + 0 Picard really is 0). Do NOT front-load Picard where Newton works: measured, + opening with 5 / 25 Picard steps cost 12 / 30 total iterations against pure + Newton's 7. + +2. **Newton-first; spend Picard only to rescue.** When Newton fails + (`DIVERGED_LINE_SEARCH` / `DIVERGED_LINEAR_SOLVE`) **or stalls admissibly** + (steps accepted, residual flat — no FAIL reason ever fires), revert to the best + iterate and buy a Picard block: `solve(picard=N)`, or + `consistent_jacobian="continuation"` (staged Picard→Newton α-blend; α is a + `constants[]` atom, no recompile). Rescue on *stagnation*, not only on a + failure reason — the failure-only trigger measured byte-identical to doing + nothing at the cliff. + +3. **Grid sequencing is the validated warm start for hard problems.** Solve + coarse (it finds the localisation structure cheaply), transfer the state up + (the linear-exact local RBF, #430), warm-start the fine solve. Measured on the + notch: 2–3× deeper residual, more localised, fewer iterations than any cold + fine strategy. No packaged API yet — hand-roll the cascade with + `uw.function.evaluate` per level; PETSc's `-snes_grid_sequence` does NOT work + on UW3 meshes. See `docs/developer/design/multilevel-nonlinear-stokes-strategy.md`. + +4. **`solver.has_solution` / tri-state `zero_init_guess`** make warm-started + campaigns safe: `None` (default) auto-detects; a diverged solve or a remesh + clears the flag so the next solve cold-starts (with its Picard entry) instead + of warming off a corrupted iterate. + +--- + +## The retired doctrine — do not resurrect it + +An earlier line of work paired the δ-soft-min with a **yield homotopy** and shipped +a model-level enable method for an in-SNES δ-ramp. That API **has been removed from +the source**, and the doctrine it taught rested on a unit-scaling error in the +campaign that motivated it. +Re-measured on the correctly-scaled problem (13 points across two parameter axes): + +- the δ-march **never succeeded where a direct hard-Min solve failed**, and where + both work the direct solve is 4–5× faster with better residuals; +- homotopy rescues **Picard**, not Newton — under the consistent tangent it adds + nothing; +- ramping δ **inside** a single SNES solve is separately proven dead (diverges + `DIVERGED_LINEAR_SOLVE` within ~2 iterations even on the proven config). + +The ruling that closed the campaign: **regularise the PROBLEM (give the shear band +a physical length scale), not the solver.** Where a hard-Min solve will not +converge, sharpening δ is not the missing lever — a viscous seed, the Picard +rescue, and grid sequencing are. + +--- + +## Which tangent for which model (measured) + +`solver.consistent_jacobian` takes `False` | `True` | `"continuation"`: + +| Model | Use | Why | +|-------|-----|-----| +| `ViscoPlasticFlowModel` (non-elastic) | **`True`** (Newton) | Quadratic near the solution; the automatic Picard entry handles the cold start. | +| `ViscoElasticPlasticFlowModel` (VEP) | **`False`** (Picard) | The consistent yield tangent over the elastic stress-history block makes the Jacobian **indefinite → `DIVERGED_LINEAR_SOLVE`**. Picard is contractive. | +| `TransverseIsotropicVEPFlowModel` (TI-VEP) | **`False`** (Picard) | Same as VEP (elastic). | +| Any, far from the solution | **`"continuation"`** | Staged Picard→Newton; Picard locates the basin, Newton finishes. Beat pure Newton at every notch point measured — but its stage switch is a residual LEVEL and one-way, so it can overspend Picard on easy problems. | + +> Measured: VEP loading-through-yield — Picard converges (σ locks at τ_y), +> Newton diverges every step (`DIVERGED_LINEAR_SOLVE`). --- ## Confirm you are actually running Newton -A consistent-Newton solve should converge **quadratically** — the residual roughly -squares each iteration and reaches ~1e-12 in 3–6 nonlinear steps. A **linear** tail -(residual dropping by a roughly constant factor over ~15–25 steps) means you are on the -Picard tangent — check `solver.consistent_jacobian is True` and that the viscosity is a -field of the unknowns, not a constant. +A consistent-Newton solve on a smooth-enough problem converges **quadratically** — +the residual roughly squares each iteration and reaches ~1e-12 in 3–6 nonlinear +steps. A **linear** tail (a roughly constant reduction factor over ~15–25 steps) +means you are on the Picard tangent — check `solver.consistent_jacobian is True` +and that the viscosity is a function of the unknowns, not a constant. (On genuinely +hard localising problems the quadratic phase may never be reached — that is the +problem, not the tangent; see the doctrine above.) -Direct symbolic check that the Newton term is in the Jacobian (`dF1/dL` differs between -the frozen and unwrapped flux by exactly the `∂η/∂(grad v)` term): +Direct symbolic check that the Newton term is present (`dF1/dL` differs between the +frozen and unwrapped flux by exactly the `∂η/∂(grad v)` term): ```python import sympy from underworld3.function.expressions import unwrap_expression -F1 = sympy.Array(stokes.F1.sym) # residual flux (η wrapped) -L = sympy.Array(stokes.Unknowns.L) # velocity-gradient symbols +F1 = sympy.Array(stokes.F1.sym) +L = sympy.Array(stokes.Unknowns.L) G_picard = sympy.derive_by_array(F1, L) F1_unwrapped = sympy.Array( [unwrap_expression(e, mode="symbolic_keep_constants") for e in F1], F1.shape) G_newton = sympy.derive_by_array(F1_unwrapped, L) -assert sympy.simplify(sympy.Array(G_newton) - sympy.Array(G_picard)) != \ - sympy.Array([0]*len(list(sympy.flatten(F1)))*len(list(sympy.flatten(L)))) -# nonzero difference == the Newton form is present +# a nonzero difference == the Newton form is present ``` --- -## Which tangent for which model (measured) - -`enable_yield_homotopy(consistent_tangent="auto")` is **model-aware**: +## δ smoothing — a modelling choice, not a convergence strategy (#475 substrate) + +If you want a *rounded* yield law at all (as physics or as a formulation choice), +the substrate is three model properties; δ is a `constants[]` atom, so changing it +never recompiles: + +- **`yield_mode`**: `"min"` (default — exact hard `Min`), `"softmin"` (the + δ-parameterised family below), `"harmonic"` (a **distinct physical model**, a + parallel blend — not an approximation to `Min`). +- **`yield_smoother`**: `"sqrt"` or `"powermean"`. **δ is NOT the same parameter + in the two families**: the power mean's sharpness is `s = 1/(δ + 0.001)`, so + δ ≤ 1 and δ = 1 IS the harmonic mean; the sqrt family's δ is a percentage stress + deviation, generous entry O(10), and δ = 0 is exactly `Min`. The power mean at + δ = 0 lands within 0.07 % of `Min` — an order of magnitude inside a 1e-8 solver + tolerance. +- **`yield_anchor`**: which point is pinned to the exact law — the SIDE of `Min` + belongs to the anchor, not the family. `"onset"` (default, historical) is exact + on the unyielded branch but sits BELOW `Min` at and above yield — a *weaker* + problem than the sharp one. `"yield"` pins τ/τ_y = 1 exactly and sits on-or-above + `Min` everywhere; the cost is stiffer unyielded material (bounded ×2 sqrt, + ×2^δ powermean, both → 1 as δ → 0). + +**If you march δ toward the sharp law, the only sound discipline is multi-solve:** +hold δ constant for a full solve to tolerance, warm-start the next smaller δ, +sharpen only between converged solves. Never ramp δ inside one SNES solve. The +packaged march is `stokes.solve(homotopy=True)` / +`underworld3.systems.yield_continuation` — usable, with two open caveats (#473): +its documented cold-start guarantee does NOT hold on a multi-material +(`Piecewise`) yield stress, so give it a viscous pre-solve anyway; and its +adaptive step control is effectively one-shot (one early decision pins the step +for the whole march). Do not expect it to cross a cliff the direct solve cannot — +measured, it never has. -| Model | `"auto"` picks | Why | -|-------|----------------|-----| -| `ViscoPlasticFlowModel` (non-elastic) | **Newton** + δ-ramp | δ-ramp keeps the residual smooth while Newton finds the basin; both sharpen to exact Min together. From a viscous guess the consistent tangent is robust (≈3–4 quadratic iters). | -| `ViscoElasticPlasticFlowModel` (VEP) | **Picard** + δ-ramp | The consistent yield tangent over the elastic stress-history block makes the Jacobian **indefinite → `DIVERGED_LINEAR_SOLVE`**. Picard is contractive; with the δ-ramp it converges to the exact yield surface. | -| `TransverseIsotropicVEPFlowModel` (TI-VEP) | **Picard** + δ-ramp | same as VEP (elastic). | - -Override with `consistent_tangent=True` / `False` / `"continuation"`. +--- -> Measured: VEP loading-through-yield — Picard+δ-ramp converges (σ locks at τ_y), -> Newton+δ-ramp diverges every step (`DIVERGED_LINEAR_SOLVE`). +## Floors + +- **`shear_viscosity_min`** (default `-oo` = off) is applied through + `uw.maths.smooth_max`, but the default rounding scale is zero under + `yield_mode="min"` and `δ·|floor|` under the smooth modes — so it **vanishes as + δ → 0**, leaving an exact `Max` corner that kills the consistent tangent + (`nl=0, DIVERGED_LINEAR_SOLVE`). Set **`viscosity_min_rounding`** (a few per + cent of the floor) and the cutoff is differentiable at any δ, including 0. +- A viscosity floor bounds the viscosity contrast and therefore how localised the + solution can be — relaxing it toward zero is a solution-SELECTION continuation, + independent of δ. Use it deliberately. +- **Do not add a strain-rate floor for the cold start.** At `ε̇=0`, `η_pl=+inf` + is carried correctly to the viscous branch by `Min` and by both smooth families; + only a hand-rolled product-over-sum harmonic blend breaks (`inf/inf`) — write it + as `η_ve/(1+f)`. --- @@ -98,10 +180,11 @@ Override with `consistent_tangent=True` / `False` / `"continuation"`. | Symptom | Cause | Fix | |---------|-------|-----| -| `DIVERGED_LINEAR_SOLVE`, 0 iters, VEP | consistent Newton tangent on the elastic block → indefinite | use Picard (`consistent_tangent=False`, or `"auto"` on a VEP model) | -| Line-search failure at first step | `bt` line search trips on the kink as δ steps down | `enable_yield_homotopy` sets `snes_linesearch_type="basic"` (full step); keep it | -| Converges but σ sits **below** τ_y | solving a *smoothed* problem (fixed δ>0, or a smooth surrogate) | use the δ-ramp (δ→0) so the endpoint is exact Min | -| Linear (~20-iter) convergence | running Picard when you wanted Newton | `consistent_jacobian=True` on a non-elastic model (see "Confirm" above) | +| `DIVERGED_LINEAR_SOLVE`, 0 iters, VEP | consistent Newton over the elastic block → indefinite | Picard (`consistent_jacobian=False`) | +| `DIVERGED_LINEAR_SOLVE` at nl=0 with a viscosity floor set | δ→0 leaves the floor's `Max` corner exact | set `viscosity_min_rounding` | +| Newton stalls with no divergence reason | admissible uselessness — steps accepted, residual flat | revert to best iterate, Picard block (`picard=N` / `"continuation"`); consider grid sequencing | +| Converges but σ sits **below** τ_y | a fixed δ>0 soft-min under the default `"onset"` anchor is a WEAKER law | that is the modelling choice you made — use `yield_anchor="yield"`, or δ→0 / `yield_mode="min"` for the exact surface | +| Linear (~20-iter) convergence | Picard tangent when you wanted Newton | `consistent_jacobian=True` on a non-elastic model (see "Confirm" above) | --- @@ -109,10 +192,7 @@ Override with `consistent_tangent=True` / `False` / `"continuation"`. - **`./uw build` → `amr-dev` env.** Verify `uw.__file__` is the worktree site-packages. - **Run VEP tests UNFORKED** — `pytest --forked` SIGABRTs here (fork of multithreaded PETSc). -- `enable_yield_homotopy` raises the `snes_max_it` floor to 100 and needs the model - **attached to a solver first** (`solver.constitutive_model = model`). -- `harmonic` yield mode is a **distinct physical model** (parallel blend), not an - approximation to Min — the homotopy does not apply to it. +- `harmonic` yield mode is a **distinct physical model**, not an approximation to Min. - If you project η, use a **low-order** field (P0/P1) — higher order overshoots and η is not guaranteed positive. @@ -120,11 +200,16 @@ Override with `consistent_tangent=True` / `False` / `"continuation"`. ## Reference -- `ViscousFlowModel._combine_yield`, `enable_yield_homotopy` / `_yield_homotopy_step` - in `constitutive_models.py`; `solver.consistent_jacobian` / `_jacobian_source`. -- Design: `docs/developer/design/jacobian-consistent-tangent.md`. -- Tests: `tests/test_1053_yield_homotopy.py`. -- Benchmark: `docs/examples/WIP/Benchmark/Ex_VP_Spiegelman_Benchmark.py`. +- Yield law: `ViscousFlowModel._combine_yield`, `yield_anchor`, `yield_smoother`, + `viscosity_min_rounding` in `constitutive_models.py`. +- Tangent: `solver.consistent_jacobian` / `_jacobian_source` in + `petsc_generic_snes_solvers.pyx`; design + `docs/developer/design/jacobian-consistent-tangent.md`. +- Warm start / continuation: `docs/developer/design/nonlinear-solver-homotopy-warmstart.md`; + grid sequencing: `docs/developer/design/multilevel-nonlinear-stokes-strategy.md`. +- Tests: `test_0201_solver_has_solution_warmstart.py`, `test_1055_yield_smoother.py`, + `test_1057_yield_homotopy_solve.py`, `test_1059_yield_anchor.py`. +- Solver-config traps, smoother, FMG/multigrid: the `nonlinear-solver` skill. Footnote: before this work UW3 differentiated the flux with the viscosity still wrapped, so `∂η/∂(grad v)` was dropped and viscoplastic solves silently ran the Picard From 59d02276b2646f39c655c17ef1894a36e25ce4c9 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 13 Aug 2026 11:27:09 +1000 Subject: [PATCH 4/4] The nonlinear-solver recipe demotes the delta-march to a rescue The recipe's step 2 still presented solve(homotopy=True) as the default entry point ("one call - automatic"), which is the retracted doctrine in its post-API form: the evidence that recommended a delta-march first rested on a unit-scaling error, and the driver's cold-start guarantee is broken (issue #473). The step now escalates honestly - grid sequencing first, the delta-continuation as rescue of last resort with the #473 caveat stated - and the Layer-2 status note carries the same demotion. The delta-discipline itself (constant per solve, never in-SNES) is unchanged; plasticity-solvers holds the ruling and evidence. Underworld development team with AI support from Claude Code --- .claude/skills/nonlinear-solver/SKILL.md | 52 +++++++++++------------- 1 file changed, 24 insertions(+), 28 deletions(-) diff --git a/.claude/skills/nonlinear-solver/SKILL.md b/.claude/skills/nonlinear-solver/SKILL.md index 4f915e845..33b70a5ae 100644 --- a/.claude/skills/nonlinear-solver/SKILL.md +++ b/.claude/skills/nonlinear-solver/SKILL.md @@ -27,38 +27,31 @@ Yield-law maths, tangent-per-model, quadratic-convergence check: `plasticity-sol the soft-min carries to the viscous branch (see the trap list for the one form that must be written carefully). -2. **Multi-solve δ-continuation** (NOT an in-solve ramp). Hold δ **constant** for a - full nonlinear solve to tolerance; warm-start the next, smaller δ from that - converged state; march δ down to the sharp surface. δ is a `constants[]` atom, - so each step is a recompile-free `PetscDSSetConstants` update. - - **This is now one call** — the model advertises the homotopy and `solve()` marches it: - - ```python - stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel - cm.Parameters.shear_viscosity_0 = ... - cm.Parameters.yield_stress = ... # a plain pressure-dependent yield - report = stokes.solve(homotopy=True) # smooth mode, tangent, march: automatic - report["settled_delta"] # smallest δ that converged - ``` - - `solve(homotopy=True)` sets the smooth mode (softmin + power-mean), picks the - tangent the model asks for (Newton for DP, Picard for elastic VEP), and runs a - residual-guided march that accelerates on easy steps and reverts + retries a - failed δ more gently. Tune with - `homotopy_options=dict(delta0=…, down=…, dmin=…, entry_maxit=…, step_maxit=…)`; - the driver is also callable directly as - `underworld3.systems.yield_continuation`. +2. **If a single solve at the sharp surface fails**, escalate in this order: + **grid sequencing first** (solve coarse, transfer, re-solve fine — the + measured 2-3x win at the notch), and only then a **multi-solve + δ-continuation** as the rescue of last resort. The δ-discipline, when you do + reach for it: hold δ **constant** for a full nonlinear solve to tolerance; + warm-start the next, smaller δ from that converged state; march down to the + sharp surface. δ is a `constants[]` atom, so each step is a recompile-free + `PetscDSSetConstants` update. + + The packaged driver is `stokes.solve(homotopy=True)` (also callable directly + as `underworld3.systems.yield_continuation`; tune with + `homotopy_options=dict(delta0=…, down=…, dmin=…, entry_maxit=…, + step_maxit=…)`). **Treat it as a rescue, not the default**: the evidence + that once made a δ-march the recommended entry point was retracted (it + rested on a unit-scaling error — `plasticity-solvers` carries the ruling + and the surviving evidence), and the driver's documented cold-start + guarantee does not currently hold (issue #473: entry can fail on a + pressure-dependent yield, and the step control is effectively one-shot). + Newton + the automatic Picard entry handles the standard cases without it. 3. **Consistent-Newton tangent** for non-elastic DP (`consistent_jacobian=True`); **Picard** for elastic VEP — see `plasticity-solvers` for the per-model table. 4. **`bt` line search** with the consistent tangent on a smooth (δ>0) surface. -The δ-march is cheap: with the power-mean smoother a converged δ warm-starts every -sharper δ in ≈0 Newton iterations, so a residual-guided auto-descent costs almost -nothing. - --- ## DO NOT ramp δ inside one SNES solve @@ -128,10 +121,13 @@ if stokes.has_solution: - **Layer 3 — DONE:** the FMG velocity smoother defaults to `gmres`+`sor` with `mg_levels_ksp_norm_type=none` (fixed-cost V-cycle), unconditionally — see "Multigrid depth" below. -- **Layer 2 — DONE:** the model advertises the homotopy +- **Layer 2 — SHIPPED, DEMOTED TO RESCUE:** the model advertises the homotopy (`supports_yield_homotopy` / `_yield_homotopy_control`) and `stokes.solve(homotopy=True, homotopy_options=...)` runs the residual-guided - continuation, returning the march summary. + continuation, returning the march summary. The doctrine that made this the + recommended entry point was retracted (unit-scaling error — see + `plasticity-solvers`), and its cold-start guarantee is broken (issue #473); + use it after Newton + Picard entry and grid sequencing have failed. ---