diff --git a/.claude/skills/adapt-on-top-faults/SKILL.md b/.claude/skills/adapt-on-top-faults/SKILL.md index 207bf3bd2..48cd277fa 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 @@ -146,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 @@ -219,6 +294,53 @@ 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 — 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) | symptom / edge | cause & handling | @@ -230,8 +352,13 @@ for step in range(nsteps): | 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. | -| NVB at np>1 raises NotImplementedError | native `_nvb_transform` extension not built (needs the custom-PETSc/amr env). | +| 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. | +| `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..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 @@ -70,9 +63,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.) --- @@ -127,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. --- @@ -167,6 +164,63 @@ 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** — 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") +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: + +- **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 + 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`). + 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. +- **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** 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 +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 +237,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`. 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