From 28a802780bc45e92bdc82f783bffbf6635090f19 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 27 Sep 2026 10:52:17 -0700 Subject: [PATCH] Capability guides live in docs; skills are symlinks; a family names its guides MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Curated guidance the code cannot state about itself — which transport scheme, which boundary treatment, how to make a hard solve converge — was split between the AI skills in .claude/skills and the rulings in CLAUDE.md, two paths with no index and no reader outside an assistant. Each guide is now one MyST page under docs/developer/guides with front matter naming the families it applies to; the seven skills are symlinks to those pages, so a skill cannot drift from its guide. Two guides are new: transport schemes, drafted from what the tests established this month with the rulings still open marked as such, and the boundary-condition rulings, gathered from CLAUDE.md and the issues. uw.capabilities() reads the front matter and lists a guide beside its family; the class-level view does the same, so uw.systems.Stokes.view() ends with the guides that apply to Stokes; the server serves them with uw_guides and uw_guide. The developer index carries the table. The adversarial review contract gains the rule, and CLAUDE.md the pointer: a change to a family is reviewed against every guide that names it, in the same change. test_0030 fails on a copied skill, a guide without front matter, or a family no class carries. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01Na7qBenCp67rDTZhFGTh5V --- .claude/skills/adapt-on-top-faults/SKILL.md | 370 +--------------- .claude/skills/adaptive-meshing/SKILL.md | 406 +---------------- .claude/skills/cetz-figures/SKILL.md | 144 +------ .../skills/free-surface-convection/SKILL.md | 225 +--------- .claude/skills/nonlinear-solver/SKILL.md | 322 +------------- .claude/skills/plasticity-solvers/SKILL.md | 227 +--------- .claude/skills/uw-visualisation/SKILL.md | 126 +----- CLAUDE.md | 8 + docs/developer/guides/adapt-on-top-faults.md | 371 ++++++++++++++++ docs/developer/guides/adaptive-meshing.md | 407 ++++++++++++++++++ docs/developer/guides/adversarial-review.md | 10 + .../guides/boundary-condition-rulings.md | 89 ++++ docs/developer/guides/cetz-figures.md | 145 +++++++ .../guides/free-surface-convection.md | 226 ++++++++++ docs/developer/guides/nonlinear-solver.md | 323 ++++++++++++++ docs/developer/guides/plasticity-solvers.md | 228 ++++++++++ docs/developer/guides/transport-schemes.md | 77 ++++ docs/developer/guides/uw-visualisation.md | 127 ++++++ docs/developer/index.md | 29 ++ src/underworld3/constitutive_models.py | 10 +- .../cython/petsc_generic_snes_solvers.pyx | 7 + src/underworld3/mcp/__init__.py | 25 ++ src/underworld3/systems/ddt.py | 7 + src/underworld3/utilities/capabilities.py | 79 +++- tests/test_0030_capability_guides.py | 62 +++ 25 files changed, 2235 insertions(+), 1815 deletions(-) mode change 100644 => 120000 .claude/skills/adapt-on-top-faults/SKILL.md mode change 100644 => 120000 .claude/skills/adaptive-meshing/SKILL.md mode change 100644 => 120000 .claude/skills/cetz-figures/SKILL.md mode change 100644 => 120000 .claude/skills/free-surface-convection/SKILL.md mode change 100644 => 120000 .claude/skills/nonlinear-solver/SKILL.md mode change 100644 => 120000 .claude/skills/plasticity-solvers/SKILL.md mode change 100644 => 120000 .claude/skills/uw-visualisation/SKILL.md create mode 100644 docs/developer/guides/adapt-on-top-faults.md create mode 100644 docs/developer/guides/adaptive-meshing.md create mode 100644 docs/developer/guides/boundary-condition-rulings.md create mode 100644 docs/developer/guides/cetz-figures.md create mode 100644 docs/developer/guides/free-surface-convection.md create mode 100644 docs/developer/guides/nonlinear-solver.md create mode 100644 docs/developer/guides/plasticity-solvers.md create mode 100644 docs/developer/guides/transport-schemes.md create mode 100644 docs/developer/guides/uw-visualisation.md create mode 100644 tests/test_0030_capability_guides.py diff --git a/.claude/skills/adapt-on-top-faults/SKILL.md b/.claude/skills/adapt-on-top-faults/SKILL.md deleted file mode 100644 index 9fec8795a..000000000 --- a/.claude/skills/adapt-on-top-faults/SKILL.md +++ /dev/null @@ -1,369 +0,0 @@ ---- -name: adapt-on-top-faults -description: Recipe for Underworld3 FAULT models on an NVB adapt-on-top mesh — resolve a fault Surface by LOCAL refinement (mesh.adapt(metric, max_levels=...) returns a child; NVB is the 2D default engine), drive it from the fault's EXACT signed distance, use rotated strong free-slip (composes with transverse-isotropy where Nitsche does not), recover dynamic topography from the constraint reaction, and run advection-diffusion on the adapted mesh with field transfer across re-adaptation. Reach for THIS for instantaneous/coupled fault-flow problems. For MMPDE node-movement convection use the `adaptive-meshing` skill instead; for rendering use `uw-visualisation`. ---- - -# adapt-on-top-faults - -The validated recipe for **fault problems on a locally-refined (adapt-on-top) mesh**. -Distilled from the annulus fault study (2026-07, `feature/adapt-on-top`). - -**This is the REFINEMENT paradigm**, not the mover one: -- `mesh.adapt(metric, max_levels=...)` bisects the base finest **locally** and returns - a **new child mesh** (`child.parent is mesh`). It is *adapt / re-adapt*, NOT node - movement — non-cumulative (each call re-marks from the static base). The child owns - a custom-P geometric-MG (FMG) tail so solvers on it get multigrid for free. -- For **MMPDE / equidistribution node movement** (deforming the same mesh to a field) - use the **`adaptive-meshing`** skill instead. Different tool; don't mix them up. - -Reference implementations (copy from these — all run np1/2/4): -`~/+Simulations/nvb_parallel_fault_study/` (weak-fault Stokes + FMG), -`~/+Simulations/shear_box_fault_study/` (iso vs TI, orientation sweeps), -`~/+Simulations/annulus_fault_study/` (rotated free-slip + topography + moving fault + -advection-diffusion). Companion skills: `uw-visualisation`, `adaptive-meshing`. - ---- - -## The core loop (fault → metric → adapt → child) - -```python -import underworld3 as uw, numpy as np, sympy - -# Base mesh MUST be built with refinement>=1 (supplies the coarse MG tail that NVB -# extends). Base cellSize ~ 2x the target h_near is the SWEET SPOT (see gotchas). -base = uw.meshing.Annulus(radiusInner=0.55, radiusOuter=1.0, cellSize=0.08, - refinement=1, qdegree=3) # or UnstructuredSimplexBox(..., refinement=2) - -# fault as a Surface (polyline control points, N x 3 with z=0 in 2D) -fault = uw.meshing.Surface("fault", base, fault_pts, symbol="F") -fault.discretize() - -# METRIC = a CALLABLE built from the fault's EXACT signed distance. This is the key -# to clean, non-patchy grading: it is evaluated at each refined level's centroids, -# so it resolves itself at the new resolution (no P1-field aliasing). -metric = fault.refinement_metric_function(h_near=0.02, h_far=0.08, width=0.05, - profile="linear") -child = base.adapt(metric, max_levels=3) # -> graded child (NVB is the 2D default) -``` - -- NVB (the 2D default engine) = graded newest-vertex bisection (bounded closure, - parallel via the native `uwnvb` transform; bit-confluent serial↔parallel); - `engine=` is the advanced selector. `engine="sbr"` = uniform - patch (the default; not graded). NVB is 2D only for now. -- `max_levels` is the isotropic-equivalent depth (NVB runs `2*max_levels` bisection - passes). The metric shape decides the grading; `max_levels` just caps it. - -### Metric options (all accepted by `adapt`) -1. **callable** `metric(centroids)->M` — **preferred for faults**. Evaluated per level. -2. MeshVariable / sympy expression — sampled via `uw.function.evaluate` from the BASE - mesh → a peaked `M=1/h²` aliases → *patchy* levels. Avoid for thin features. - -**Custom refinement shape** — pass any callable. For a *fat, uniformly-fine* band -(not just a thin line at the fault), a flat-core metric: -```python -def metric(pts, _f=fault): - d = _f.unsigned_distance(pts) # EXACT distance at arbitrary points - core, ramp, hn, hf = 0.02, 0.05, 0.01, 0.08 - h = np.where(d < core, hn, np.minimum(hn + (hf-hn)*(d-core)/ramp, hf)) - return 1.0 / h**2 -``` - -`Surface` distance API (all exact, arbitrary query points): -`fault.signed_distance(coords)`, `fault.unsigned_distance(coords)`, -`fault.director` (unit normal = normalised ∇(signed distance); the TI weak-plane -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 -needs the fault's distance/normal ON THE CHILD. Two ways: - -```python -# (A) re-home the SAME fault onto the child (cleanest; distance recomputes on child) -fault.remap_to(child) -eta = fault.influence_function(width=0.04, value_near=1e-3, value_far=1.0, - profile="smoothstep") # isotropic weak zone - -# (B) or build a child-side Surface (needed if the base fault is still in use for a -# base-mesh metric — symbol disambiguation refuses a base-mesh symbol in a child -# solver). fault_c = uw.meshing.Surface("fault_c", child, fault_pts); fault_c.discretize() -``` - -Isotropic weak zone: -```python -stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel -stokes.constitutive_model.Parameters.shear_viscosity_0 = eta # drops near fault -``` - -Transverse-isotropic weak PLANE (the physical fault; low fault-parallel shear): -```python -stokes.constitutive_model = uw.constitutive_models.TransverseIsotropicFlowModel -stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_bulk # normal viscosity (constant) -stokes.constitutive_model.Parameters.shear_viscosity_1 = eta # weak near fault -> bulk far -stokes.constitutive_model.Parameters.director = fault.director # unit fault normal -``` -The **TI Jacobian is the full consistent tangent** (not isotropic + defect -correction — that framing is WRONG). The TI velocity FMG V-cycle needs a few more -iters than iso (~8 vs ~2) because of a directional near-null mode the isotropic -point-smoother doesn't damp — bounded and contrast-independent, not a bug. TI needs -~2 elements across the weak zone to resolve (iso ~1). - ---- - -## Rotated strong free-slip (the reason to use this, not Nitsche) - -Nitsche free-slip is INCOMPATIBLE with the TI model (its penalty scales by an -isotropic viscosity). Rotated strong free-slip imposes `u·n̂=0` as an ESSENTIAL -constraint in a per-node (n,t) frame → machine-zero leakage AND composes with TI. - -```python -# normal=None (the default) is measure-weighted and consistent with the assembly — -# prefer it. An analytic nhat is exact for the TRUE circle but keeps a consistency -# error against the faceted integral (#560); use it only when the constraint must -# follow the geometry rather than the mesh. -nhat = mesh.CoordinateSystem.unit_e_0 # exact radial normal (annulus/sphere) -stokes.add_rotated_freeslip_bc(0, "Upper", normal=nhat) -stokes.add_rotated_freeslip_bc(0, "Lower", normal=nhat) -stokes.petsc_use_pressure_nullspace = True # enclosed -> pressure gauge -stokes.solve() # rigid-rotation gauge auto-removed -# Convergence status: read stokes._rotated_freeslip_info = {ksp_reason, -# nonlinear_iterations, rotation_gauge_removed, reaction} — NOT s.snes.getConvergedReason() -# (the rotated solve is a manual loop, not snes.solve). Also sanity-check v·n leakage (~1e-16). -``` - -**Nonlinear rheology / warm-start / timestepping (PR #298, `feature/rotated-snes`):** -rotated free-slip now works *inside* the nonlinear iteration — a nonlinearity probe -auto-dispatches power-law / VEP / TI-with-yield models to a Newton/Picard loop -(`solve_rotated_freeslip_nonlinear`), so warm-started time loops are correct. On the -worktree state *before* #298 lands, the rotated path is a SINGLE linear solve — it -silently returns one Newton linearisation from `u=0` for a nonlinear model. If you -run nonlinear TI + timestepping with rotated free-slip, make sure #298 is in. - -**Dynamic topography** from the constraint reaction (the reason to bother): -```python -h = uw.discretisation.MeshVariable("h", mesh, 1, degree=1) -stokes.dynamic_topography("Upper", h, buoyancy_scale=rho_g) # h = -(σ_nn - mean)/ρg -xs, sig = stokes.boundary_normal_traction("Upper") # or the raw σ_nn (lumped-mass) -``` - ---- - -## FMG under rotated free-slip - -`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 - 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 - ``` - (~4 velocity iters vs ~26 GAMG). Otherwise it falls back to GAMG (fine, just slower). - -The default velocity-block preconditioner and the pressure Schur (`1/η`) are already -near-optimal for TI — do **not** hand-roll a "TI-aware Schur"; measured, the default -`1/η₀` beats every alternative (the weak TI mode is fault-parallel shear, which is -volume-preserving, so pressure sees η₀). - ---- - -## Moving fault + re-adaptation + field transfer - -`adapt()` is non-cumulative and deterministic. Move the fault, re-adapt from the -static base, carry any field by interpolation: - -```python -for step in range(N): - fault_pts = move(fault_pts, step) # kinematics - fault = uw.meshing.Surface(f"fault{step}", base, fault_pts, symbol="F"); fault.discretize() - child = base.adapt(fault.refinement_metric_function(...), max_levels=3) - # verify: folded=0 (all cell |vol|>0), base unchanged (non-cumulative) - # carry a field child_{k-1} -> child_k by interpolation: - T = uw.discretisation.MeshVariable(f"T{step}", child, 1, degree=1) - if T_prev is None: - T.data[:,0] = uw.function.evaluate(T0_expr, T.coords) - else: - T.data[:,0] = uw.function.evaluate(T_prev.sym, T.coords) # mesh->mesh interp - T_prev = T -``` -Transfer error is **non-accumulating** for a smooth field (bounded by the per-mesh P1 -representation floor, ~0.2% on a fine base, ~3% on a coarse base) — repeated -re-meshing does not diffuse a smooth field away. Sharp fronts lose more per transfer. - ---- - -## Advection-diffusion on the adapted mesh - -```python -T = uw.discretisation.MeshVariable("T", child, 1, degree=2) -adv = uw.systems.AdvDiffusionSLCN(child, u_Field=T, V_fn=stokes.u.sym, order=1, - monotone_mode="clamp") # clamp = bounded, no overshoot -adv.constitutive_model = uw.constitutive_models.DiffusionModel -adv.constitutive_model.Parameters.diffusivity = 2e-4 -adv.f = 0.0 -dt = 0.015 # FIXED dt — SLCN is semi-Lagrangian (unconditionally - # stable). estimate_dt() reports the tiny fault-band - # Courant limit; do NOT use it to size the step. -for step in range(nsteps): - adv.solve(timestep=dt) -``` -- **SLCN dt is NOT Courant-limited** by the fine fault-band cells — pick dt from the - coarse-region advection, verify accuracy. -- The scalar AD auto-FMG-injection bug (velocity-block custom-P mismatching a scalar - operator on adapt children → PtAP error 60) is **FIXED**: `auto_inject_custom_mg` - now calls `snes.setUp()` before reading the finest reduced map (so the DM section - is the finalized space the operator lives on) and validates that map against the - assembled operator. Custom-P installs successfully on scalar SLCN AdvDiffusion on - NVB adapt children (no skip-guard, no workaround) — `bugfix/custom-mg-parallel`. - ---- - -## 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 | -|---|---| -| patchy along-fault refinement (level 4 here, 2 there) | P1-interpolated `M=1/h²` aliasing. Use the **callable exact-distance** metric. | -| child solver rejects `fault.distance.sym` (foreign-mesh error) | symbol disambiguation. `fault.remap_to(child)` or build a child-side `Surface`. **Rough edge**: needing two Surface objects. | -| adapt added-cell count jitters ±30% step-to-step | small-number geometry on a thin band. Deterministic; quality (near-fault h) is constant. Base ≈ **2× target** minimises it (CoV ~10% at 0.08 base vs 19% at 0.14, 24% at 0.05). | -| tiny AD timestep / slow advection | `estimate_dt` returns the fault-band Courant limit. Use a fixed dt (SLCN is unconditionally stable). **Rough edge**: `estimate_dt` not adapt-aware. | -| 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 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. | -| scattered 1-level refinement across the WHOLE domain; far-field quality drops | the metric's far clip ``h_far`` sits below the base mesh's cell DIAMETERS (gmsh ``cellSize`` is a target edge length; diameters run 1.2–2.5x it — measured 0.108–0.223 for cellSize 0.18+ref 1). ``edge_split`` marks on diameter, so 100% of the domain refines once and the unrequested bisection DE-CONDITIONS the grid (far-field median q 0.372 -> 0.295, measured). Set ``h_far >= 1.05 * cell_diameters(base.dm).max()``. | -| 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/adapt-on-top-faults/SKILL.md b/.claude/skills/adapt-on-top-faults/SKILL.md new file mode 120000 index 000000000..e81a9ec22 --- /dev/null +++ b/.claude/skills/adapt-on-top-faults/SKILL.md @@ -0,0 +1 @@ +../../../docs/developer/guides/adapt-on-top-faults.md \ No newline at end of file diff --git a/.claude/skills/adaptive-meshing/SKILL.md b/.claude/skills/adaptive-meshing/SKILL.md deleted file mode 100644 index c2c98da26..000000000 --- a/.claude/skills/adaptive-meshing/SKILL.md +++ /dev/null @@ -1,405 +0,0 @@ ---- -name: adaptive-meshing -description: The canonical, workable recipe for Underworld3 moving-mesh / adaptive-mesh convection (annulus stagnant-lid, faults, free surface). Reach for THIS first when setting up any model with a deforming or adapted mesh — it encodes the combination that does not blow up, tangle, or inject spurious energy, and explains the failure modes so you don't re-derive them. Use before choosing movers, free-slip BCs, restart, or field-transfer options. ---- - -# adaptive-meshing - -The one workable combination for UW3 moving/adaptive-mesh convection, distilled -from many sessions that each re-picked options and stomped on each other's -defaults. **Start from this recipe; change one thing at a time and verify.** - -Reference implementation (current, validated): the **`underworld3.workflows` -adaptive-convection example** — -`docs/examples/workflows/adaptive_convection/` on the `feature/adaptive-convection` -worktree (`config.py`+`simulate.py` no-fault; `fault_config.py`+`fault_simulate.py` -fault; `diagnostics.py`, `render.py`, `compare.py`). Express adaptive runs as a -WORKFLOW (a `WorkflowConfig` + `@workflow_step` DAG + `Run`), NOT a monolithic -driver. The older `scripts/fault_convection_adapt_loop.py` (feature/fault-convection) -is superseded — its ideas are folded into the workflow + this skill. -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(...)` 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. 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. - ---- - -## Mover quick-start (copy-paste — this is the hard-to-discover bit) - -The user entry is `uw.meshing.node_redistribution(mesh, metric, ...)` (the -purposeful spelling; it dispatches to `mesh.redistribute_nodes`, which drives -the MMPDE mover on 2D simplex meshes — `smooth_mesh_interior` is the -machinery underneath and takes the same kwargs). Minimal correct setup to -adapt a mesh to a field `T` each step: - -```python -import underworld3 as uw - -# metric from |grad T|: refinement=R is a factor on the BACKGROUND spacing h0, -# not a finest:coarsest ratio. The envelope is h in [h0/R, h0*coarsening], and -# coarsening="auto" is R**(1/d) — so R=5 in 2-D spans h0/5 to 2.2*h0, a ratio -# of R**(1+1/d) ~ 11. Use refinement=R, NOT strategy= (caps at ~2, under-grades). -rho = uw.meshing.metric_density_from_gradient( - mesh, T, refinement=5, coarsening="auto", metric_choice="front-following") - -# move the mesh — the mover (Huang-Kamenski MMPDE) is variational, -# non-folding, clusters AND aligns cells. It OWNS field transfer -# (remaps T + SLCN history, fires on_remesh hooks). -uw.meshing.node_redistribution( - mesh, rho, - method_kwargs=dict(step_frac=0.2, accel="cg", momentum=0.0), # mmpde's OWN kwargs - slip_surfaces=True, # boundary nodes slide tangentially (parallel-safe) - skip_threshold=0.9) # skip the move when the mesh is already aligned -``` - -For a **sharp feature (fault)** pass an anisotropic SPD TENSOR metric instead of -the scalar `rho` (thin ACROSS the feature normal n) and bake a gmsh base: - -```python -import sympy -n = sympy.Matrix([nx, ny]) # constant fault-normal unit vector -d = dfac.sym[0] # DIRECT unsigned distance field (P1) -M = rho * sympy.eye(2) + (Rf**2 - 1.0) * sympy.exp(-(d/w)**2) * (n * n.T) -# mesh built with: uw.meshing.Annulus(..., refine_lines=[xy], refine_size_min=smin) -uw.meshing.node_redistribution(mesh, M, - method_kwargs=dict(step_frac=0.2, accel="cg", momentum=0.0), - slip_surfaces=True, skip_threshold=None) # tensor metric: do the skip check yourself -``` - -Pitfalls that make it "not work": `method="anisotropic"`/`"ot"`/`"spring"`/`"ma"` -(RETIRED 2026-07 — they now raise ValueError; mmpde is the default and only -metric mover); injecting `relax`/`n_outer` (starves mmpde's CG); `strategy=` instead of -`refinement=R` (under-grades); a scalar bump for a fault (refines a fat corridor, -leaves the centre coarse); signed `Surface.distance.sym` for `d` (bleeds along the -line extension — use a direct unsigned distance). Full rationale + the rest of the -recipe (BCs, restart, field transfer, cadence) below. - ---- - -## The canonical recipe (defaults that work) - -### 1. Mover — mmpde -`uw.meshing.smooth_mesh_interior(mesh, metric=..., method="mmpde", -method_kwargs=dict(step_frac=0.2, accel="cg", momentum=0.0), slip_surfaces=True)`. -- mmpde = Huang–Kamenski variational, **non-folding** (energy → ∞ as detJ → 0), - clusters AND aligns cells. It is the only clean mover. -- **NOT** the `anisotropic` mover (shreds/freezes on static features) or `OT` - (slivers — OT is optimal transport, sliver-prone, not a route around the cap). -- Do **not** inject `relax`/`n_outer` into mmpde (those are the anisotropic - mover's knobs and starve mmpde's internal CG). - -### 2. Metric -- Thermal: `metric_density_from_gradient(mesh, T, refinement=R, - metric_choice="front-following")`. `refinement=R` (≈5) is the maximum local - refinement **on the background cell size h0**, not the finest:coarsest ratio: - the metric targets `h ∈ [h0/R, h0·coarsening]`, and `coarsening="auto"` takes - the budget-conserving `R**(1/d)`. So R=5 in 2-D asks for h0/5 up to 2.2·h0 — - a finest:coarsest ratio of `R**(1+1/d)` ≈ 11, and ≈ 8.5 in 3-D. Named - `strategy=` caps at ~2 and under-grades. R≈5 extracts ~all the grading the - node budget/layout allows; don't over-tune R (benign no-op above budget). - Passing `refinement` takes the **envelope branch**, which ignores `amp`, - `lo/hi_percentile`, `mode` and `power`. -- Fault / sharp feature: a **hand-built anisotropic SPD tensor** - `M = ρ·I + (Rf²−1)·exp(−(d/w)²)·n nᵀ` (thin ACROSS the feature normal n). A - scalar bump refines a fat isotropic corridor and leaves the centre-line coarse. -- Use a DIRECT unsigned distance field (geometry_tools) for `d`, NOT the Surface's - signed `.distance.sym` (its zero-contour bleeds the metric along the line - extension). - -### 3. Creation vs maintenance (the cap) — and the gmsh base COMPOUNDS -mmpde **cannot create** strong refinement from a uniform mesh — it saturates at a -fixed-topology cap (~1.8× on the fault, re-measured), because the mover has a FIXED -node budget (it redistributes, never adds nodes). To go finer you need MORE NODES, -which only gmsh can add at construction. **Bake refinement into the gmsh base** -(`Annulus(refine_lines=[xy], refine_size_min=...)` — now real, see the Faults -section) and the mover doesn't just MAINTAIN it: the extra gmsh nodes lift it off -its budget cap so it **compounds** (gmsh f2 base 0.44 → mover 0.29; f3 base 0.30 → -mover 0.19 ≈ 5× finer). Measured (Ra1e6 Rf8 res24, all folded=0): -uniform 0.55 (~1.8×) → gmsh-f2 0.29 (~3.4×) → gmsh-f3 0.19 (~5×). Judge by -fault/bulk nearest-neighbour spacing RATIO, never by global misalignment. - -### 4. Cadence — adapt as often as you like (forced every step is FINE) -Adapt every step or every few — on a CORRECT build (see §5) the mover converges -dead-flat and forced every-step adaptation under vigorous convection is stable -(validated: 39/39 forced adapts, mesh folded=0, area-ratio ~14 flat). Skipping -when aligned just saves cost. (An earlier claim that forcing adaptation -"tangles / over-injects energy" was WRONG — that was the deformed-eval bug in §5.) - -### 5. THE bug that wrecked adaptation (holes) — and the fix -SYMPTOM: giant empty cells / holes in the adapted mesh, intermittent area-ratio -spikes, convection wrecked. ROOT CAUSE (proven): `uw.function.evaluate` -**mis-locates points on a deformed mesh** (the nav kd-tree `mesh._nav_coords` was -captured from the ORIGINAL coords and never refreshed), so the metric — built with -strictly-positive nodal values — evaluates to **NEGATIVE garbage** (even at its own -DOFs) → non-SPD → the mover wrecks the mesh. FIX (both): -- **`8a9d2ff2`** (refresh `_nav_coords` + projected normals on every deform) — - the real fix; makes `function.evaluate`/`points_in_domain` track deformation. -- **Monotone RBF metric bake** (in `_mmpde_mover`, formerly `_winslow_mmpde`): Shepard-interpolate the - metric from its **positive nodal values** — a convex average is guaranteed ≥0 - (monotone) + fast (no cell-location). Use RBF for the metric; it doesn't need - high-precision eval. -With these, the metric stays positive/SPD; #259's SPD-floor never fires (harmless). -The mover, `accel`, and `refinement=R` (e.g. R=5) are all FINE — they were -red-herring symptoms of the eval bug. Always confirm a CLEAN BUILD first -(`./uw build`; `md5 site-packages/.../smoothing.py == src`); a stale build is its -own cause of holes. - -### 6. Stokes free-slip -Penalty (`add_natural_bc(KFS·v·n·n)`, KFS≈1e6) or Nitsche γ=10 both work on a -UNIFORM mesh. Nitsche preferred for sharp fault corners. Do NOT diagnose free-slip -with nodal v·n: Nitsche enforces v·n=0 weakly, so nodal v·n is large even when -correct — use vrms (kinetic energy). (An earlier "warm-restart × Nitsche → blow-up, -use cold restart" diagnosis was WRONG / confounded by the §5 eval bug + a stale build.) - -**On an ADAPTIVE mesh, Nitsche free-slip γ=10 is TOO SOFT → intermittent vrms -spikes at adapt steps** (under-enforcement, NOT a blow-up: single-step velocity -garbage that T survives because the solve is cold-restarted — e.g. vrms -144→8887→144). The mover pins the boundary and refines just beneath it, creating -high-aspect-ratio near-boundary cells whose Nitsche inverse-estimate constant needs -a bigger penalty. Two fixes, BOTH good (corroborated across sessions): -- **Nitsche γ=100** (not 10) for a free-slip top on an adaptive mesh — clean, smooth - signal. This is INDEPENDENT of the penalty's mesh-size scaling: local vs global `h` - are equivalent here (both garbage at γ=10, both clean at γ=100). -- **Don't fully pin the boundary in the mover** — let it tangential-slip - (`slip_surfaces=True`, the §1 default), NOT `pinned_labels=[Upper,...]`. Avoiding - the distorted near-boundary cells lets γ=10 work. Pin a boundary only when you must - hold a prescribed shape (e.g. a free-surface height the integrator just set). - -Aside — free-SURFACE held-lid is the OPPOSITE problem: the GLOBAL `h = -get_min_radius` drifts BELOW the surface cell as the interior refines → the Nitsche -penalty OVER-stiffens → a spurious one-step surface "mountain" (vhmax spike). Fix = -a LOCAL per-cell `h` via `mesh.cell_size()` (deformation/adaptation-tracking); -`add_nitsche_bc(..., local_h=True)` is the default (PR #275). So held-lid wants LESS -penalty (local-h), free-slip-adaptive wants MORE (γ=100) — different knobs. - -### 7. Advection + timestep -`AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V.sym, theta=1.0, monotone_mode='clamp')`. -theta=1.0 (backward Euler) for stability; monotone clamp bounds SL overshoot; -`V_fn=V.sym` is the PHYSICAL velocity. **dt can be LARGE — SLCN is unconditionally -stable; the smallest/median adapted cell does NOT have to govern dt.** Use -`estimate_dt(percentile=50)` × a multiplier (e.g. dt_mult 3–5) to advance physical -time faster (needed to develop convection in stiff stagnant-lid regimes). - -### 8. Rheology -- **isotropic** (linear) ⇒ `snes_type=ksponly` (one exact KSP solve; default - newtonls rejects steps on vigorous flow). Cheap. Use for resolved features. -- **TI / anisotropic** weak fault (`TransverseIsotropicFlowModel`: - `shear_viscosity_0=η_FK`, `shear_viscosity_1=η_weak`, `director`=fault normal): - keep the **default newtonls** (NOT ksponly — ksponly converges to the WRONG - answer); use **penalty** free-slip (Nitsche trips the anisotropic-Jacobian bug); - **GAMG** (FMG ~7× slower). Measured on the gmsh-resolved + one-sided base - (Ra1e6 Δη1e3): ran clean to t=0.06, folded=0, vrms→20, **~10×** the isotropic - cost per step (GAMG eats the 1000× anisotropic contrast — well under the feared - 20×). The TI fault visibly steers the flow (persistent recirculation at the trace). -- Viscosity-bearing fields are **P0/P1 only** (positivity; higher order overshoots). - FK viscosity `η = exp(θ(1−T))`, θ = ln(Δη). Floor a weak zone, don't multiply → 0. - -### 9. Build discipline -- `./uw build` after ANY source change; verify `uw.__file__` is in site-packages - and (when debugging) that `site-packages/.../smoothing.py` matches `src/...`. A - **stale mover build is a top cause of "giant empty elements / holes"** in the - adapted mesh. -- **NEVER** `pip install -e .` (contaminates all envs). Run from inside the worktree - (`pixi run -e amr-dev`). - ---- - -## Faults — the gmsh-resolved, on-fault, one-sided recipe (validated 2026-06-23) - -The full fault recipe, hard-won. Reference implementation: the -`underworld3.workflows` adaptive-convection example -(`docs/examples/workflows/adaptive_convection/{config,fault_config}.py`, on the -`feature/adaptive-convection` worktree — built on the FIXED mover base). Use the -WORKFLOW system, not a monolithic driver. - -### gmsh line refinement (`refine_lines`) — NOW IMPLEMENTED in `Annulus` -```python -uw.meshing.Annulus(radiusOuter=1, radiusInner=0.5, cellSize=1/24, - refine_lines=[xy], # list of (N,2) polylines (model coords) - refine_size_min=cellSize/3, # cell size ON the line (factor 2-3 is plenty) - refine_dist_min=0.02, refine_dist_max=0.12) # size ramps back to cellSize -``` -A gmsh **Distance + Threshold** field along the polyline; INTERIOR points are -embedded so nodes land ON the line. Backward-compatible (default `None`). It is a -core meshing change → lands as its own small meshing PR, separate from the workflow -example. (Before this session it was only *called* in old scripts and never existed -— don't trust `refine_lines` on any branch but the one carrying this commit.) - -### Keeping refinement ON the fault (the metric-composition trap) -The fault metric is `M = ρ·I + (Rf²−1)·exp(−(d/wₙ)²)·n nᵀ` with the isotropic SIZE -density `ρ`. **Do NOT fuse the fault density into ρ_T by PRODUCT** — the cold -surface thermal BL (ρ_T ~ R^d ~ 20–25) out-competes the fault near its top and the -refinement drifts ABOVE the fault, starving the deep fault ("seems to repel"): -- Use `ρ = max(ρ_T, fault_ρ)` (NOT `ρ_T · fault_ρ`). -- Make `fault_ρ = 1 + amp·gauss` with **amp > ρ_T** (~25 at R=5) so the fault wins - the max along its whole length. -- DIAGNOSE this by comparing the gmsh BASE (step 0, refinement centered on the - fault — correct) vs the DEVELOPED mesh (drifted above) → it's the MOVER's metric, - not gmsh. Render the mesh with the **fault trace overlaid** (`render.py --fault`). - -### One-sided fault influence (the clean control) -Even with max+amp, the symmetric metric DEMANDS both flanks while realized nodes -drift to the hanging wall. Make it one-sided: -- Store the **SIGNED** distance in `dfac` (the gaussians square it, so magnitude is - unchanged — the sign only feeds a `0.5(1+tanh(m·d/w))` gate). Probe the - radially-outward side once to define "upper" regardless of the distance tool's - orientation convention. -- `fault_metric_side` (both/upper/lower) gates the refinement; `fault_rheology_side` - gates the weak zone. **`both=upper`** is the physical recipe: a one-sided - hanging-wall damage zone with the mesh refined on the same side (refinement and - rheology coincide; gmsh-f3 → fault/bulk ~0.19, folded=0). `metric_side=lower` - instead pulls refinement onto the footwall to counter the upward drift. - -### The wedge fill (anti-collision) -The fault pull and the surface-BL pull compete for the coarse cells in the radial -sliver BETWEEN them. `fault_wedge=True` gmsh-fills that wedge (sample radial -segments from each fault point up to the surface, add as a second `refine_lines` -point set) so both pulls have their own budget and merge into one coherent fine -wedge instead of colliding. - -### Weak zone (rheology) — geometric blend + gaussian PEAKED on the fault (2026-06-24) -The weak-zone viscosity blends `η_FK` (background) and `floor` (fault) by the -influence `f`∈[0,1] (P1, positive). Get it right with TWO rules (verify by -reconstructing + rendering the REALIZED η field — `render_fields.py` — and the -combined `η_FK(T)^(1−f)` field; never assume the floor is reached): -- **GEOMETRIC blend `η_weak = η_FK^(1−f)·floor^f`** (NOT arithmetic - `η_FK·(1−f)+floor·f`, whose `(1−f)` term leaks the stiff-lid background through - → η≈8 at f=0.97 in a 1000× lid). Geometric reaches the floor genuinely (η≈1.2 at - f=0.97) AND ties the contrast to the LOCAL η_FK — so the fault automatically - bites hardest where it cuts the cold stiff lid (physically correct), nothing in - the hot interior. -- **`f` must be PEAKED on the fault (gaussian), NOT a TOPHAT block.** A top-hat - makes a uniform weak BLOCK; its sharp edges mean strain follows ∇f (TWO parallel - lines at the block edges, not on the fault), and a one-sided block is offset to - the hanging wall (η_1=1 core sits ABOVE the drawn line). Use a **gaussian - `f=exp(−(d/w)²)`** (peak f=1 ON the fault → η_1=1 on the drawn line, strain - localizes INTO the slot as a single feature). side=both = symmetric; side=upper = - hanging-wall halo (gaussian taper up, sharp footwall recovery — NO halving gate). - Centre the METRIC too (`metric_side=both`) so nodes refine on the line. -THERMAL CONTRAST (Louis's insight, confirmed): even a symmetric gaussian gives a -much bigger η-contrast on the COLD (upper/surface) side than the warm (footwall) -side, because η_FK rises ~1000× toward the surface — so the fault's dynamical -prominence naturally concentrates in the cold lid. This is a feature, not a bug. -Verified (gmsh-f3, gaussian width=0.025, side+metric=both): weak zone centred on -the drawn fault on the adapted mesh, folded=0. TUNE AT STEP 0 (build mesh+fields, -no solve — fast; plot the 1D η_1(d) profile + the 2D field). COST: genuine 1000× -TI contrast ~25–145 s/step (cold-start steps slow, ~25 s once developed; GAMG). -For TI see §8. (`fault_config.py`: `fault_profile=gaussian` (default still tophat — -pass gaussian), geometric blend in `create_solvers`; `render_fields.py` light maps.) - ---- - -## Failure modes — symptom → cause → fix - -| Symptom | Cause | Fix | -|---|---|---| -| Giant empty cells / holes in adapted mesh; intermittent area-ratio spikes | `function.evaluate` mis-locates on deformed mesh → metric → negative/non-SPD (the real bug). OR stale build. | `8a9d2ff2` + monotone RBF metric bake; `./uw build` + verify md5 site-packages==src | -| Decays when it should convect | over-diffusion OR under-resolution OR dt too small to develop | check a resolved arbiter (finer space+time); raise node budget; **larger dt** (dt_mult 3–5) | -| Fault won't refine under convection | mmpde creation cap; field gives no signal in cold lid | gmsh `refine_lines` base — the gmsh nodes let mmpde COMPOUND past the cap (~5×) | -| Refinement drifts ABOVE the fault / "repels", deep fault starved | fault density fused by PRODUCT with ρ_T → thermal BL out-competes the fault near the surface | `ρ = max(ρ_T, fault_ρ)`; `amp > ρ_T` (~25); or one-sided `fault_metric_side`; +`fault_wedge` | -| Refinement on both flanks but you want one side | symmetric (unsigned-distance) metric | signed `dfac` + `fault_metric_side`/`fault_rheology_side` tanh gate (`both=upper` = physical) | -| TI fault solve: ksponly gives wrong answer | ksponly skips the Picard the inexact GAMG inner solve needs | keep default newtonls; penalty free-slip; GAMG (~10× isotropic cost, fine) | -| free-slip ADAPTIVE: vrms spikes at adapt steps (e.g. 144→8887→144), T stays bounded | Nitsche γ=10 too soft on the distorted near-boundary cells the pinned-top+interior-refine creates (under-enforcement) | **Nitsche γ=100** (not 10); OR don't pin the boundary (tangential-slip `slip_surfaces=True`). h-scaling (local vs global) is moot here | -| free-SURFACE held-lid: spurious one-step surface "mountain" / vhmax spike | global `h=get_min_radius` drifts below the surface cell as interior refines → Nitsche OVER-stiff | LOCAL per-cell h: `add_nitsche_bc(local_h=True)` (default, PR #275) = `mesh.cell_size()` | - -**Diagnose by:** vrms (KE), Nu (`BdIntegral` surface flux), fault/bulk NN-spacing -RATIO (cKDTree), folded-element count + min cell area. Compare runs **at matched -physical time t** (dt differs between meshes), not by step number. When unsure -whether behaviour is physical, build a **resolved arbiter** (uniform mesh finer in -BOTH space and time) — if it agrees with one candidate, that's the truth. - -## Diagnostics (in the workflow example, reusable) -`diagnostics.py`: `mesh_quality` (folded / area-ratio / aspect), `nn_spacing_ratios` -(BL + fault/bulk), `NusseltSurface`, `vrms`, `History`. `render.py`: T+mesh+ -streamlines on the `Run` layout, **`--fault`** overlays the fault trace (read from -the run manifest) + **`--focus-fault`** auto-crops on it + **`--mesh-only`** for the -clean mesh, **`--all`** for every frame. `compare.py`/`fault_refine_plot.py`: -matched-physical-time comparison + fault/bulk-ratio time series. -**Rendering long runs as they go:** a completion-only Monitor is NOT enough — arm a -Monitor that POLLS for new `run.mesh.NNNNN.xdmf` checkpoints and emits the index, so -you render each step as it lands. -**Checkpoint an ADAPTIVE mesh with `meshUpdates=True`** (per-step geometry) or the -saved frames pair deformed fields with stale step-0 geometry. - -## The verified canonical command - -Requires the fix (§5): `8a9d2ff2` (deformed-mesh point-location) + the monotone RBF -metric bake in `smoothing.py`. Validated 2026-06-22 at vigorous Ra1e6/Δη1e3 -stagnant-lid: forced adapt EVERY step (39/39), larger dt, → vrms→18, Nu→1.86, -|v|max 61, mesh CLEAN every frame (folded=0, area-ratio ~14 flat), no abort. - -```bash -# no-fault baseline (the workflow CLI; one flag per config field) -pixi run -e amr-dev python docs/examples/workflows/adaptive_convection/simulate.py \ - --output-dir ~/+Simulations//baseline \ - --rayleigh 1e6 --delta-eta 1e3 --cellsize 0.0417 \ - --resolution-ratio 5 --adapt-every 1 --dt-mult 4 --max-steps 80 --max-t 0.06 - -# resolved fault: gmsh base (factor 3) + on-fault + one-sided hanging wall + TI -pixi run -e amr-dev python docs/examples/workflows/adaptive_convection/fault_simulate.py \ - --output-dir ~/+Simulations//fault \ - --rayleigh 1e6 --delta-eta 1e3 --cellsize 0.0417 --resolution-ratio 5 \ - --fault-base-smin 0.0139 --fault-anisotropy 8 \ - --metric-combine max --fault-refine-amp 25 \ - --fault-rheology-side upper --fault-metric-side upper \ - --rheology ti --dt-mult 4 --max-steps 40 --max-t 0.06 -``` - -Key choices, verified: -- **`--resolution-ratio 5`** (R=5) is fine — the mover handles it on a correct build. -- **`--adapt-every 1`** force adapt every step (the strongest mesh test). -- **`--dt-mult 4`** — larger dt for STABILITY is fine (SLCN unconditional); but it - costs transient ACCURACY (over-diffusive backward-Euler DELAYS the convective - onset vs a resolved arbiter — dt×1.5 recovers it). Use small dt_mult for faithful - transients, large to reach quasi-steady fast. -- **`--freeslip penalty`** (default) — REQUIRED for TI (Nitsche trips the - anisotropic-Jacobian bug). It's the raw velocity penalty `kfs·(v·n)n`. -- Fault: `--fault-base-smin` (gmsh resolve), `--metric-combine max` + - `--fault-refine-amp 25` (keep refinement ON the fault), `--fault-*-side` (one-sided), - `--rheology ti` (real fault). Drop the fault flags for the no-fault control. - -Render with `render.py` (`--fault --focus-fault` to see the trace + refinement -coincide; `--mesh-only`; `--all`). Judge the mesh by folded/area-ratio, the physics -by vrms/Nu, ALWAYS at matched physical time. - -## Related memory -`project_adaptive_convection_as_workflow` (THIS session: workflow port, gmsh -refine_lines, on-fault/one-sided/wedge, TI, dt-accuracy), `project_mmpde_holes_real_root_cause`, -`project_fault_refine_fixed_topology_cap`, `project_uw_workflow_landing`, -`feedback_debug_adaptive_solver_method`, `project_fault_convection_working_settings`. diff --git a/.claude/skills/adaptive-meshing/SKILL.md b/.claude/skills/adaptive-meshing/SKILL.md new file mode 120000 index 000000000..13942c4af --- /dev/null +++ b/.claude/skills/adaptive-meshing/SKILL.md @@ -0,0 +1 @@ +../../../docs/developer/guides/adaptive-meshing.md \ No newline at end of file diff --git a/.claude/skills/cetz-figures/SKILL.md b/.claude/skills/cetz-figures/SKILL.md deleted file mode 100644 index f65e81a3a..000000000 --- a/.claude/skills/cetz-figures/SKILL.md +++ /dev/null @@ -1,143 +0,0 @@ ---- -name: cetz-figures -description: Build schematic / labelled-geometry figures for underworld3 papers using Typst + cetz. Use when the figure is primarily about topology, annotation, and math-typeset labels (meshes, solver diagrams, flow charts). Prefer Python → SVG → `#image()` for data-heavy figures (fields, colormaps, arrow plots) instead. ---- - -# cetz-figures - -Scaffold for Typst/cetz figures in `publications/**/figures/` alongside the -existing `arrays-sync-flow.typ`. This skill exists because upstream Claude -sessions draft cetz blind — this one actually compiles. - -## When to use cetz - -- Mesh schematics with a few labelled triangles / vertices / control points. -- Solver / data-flow diagrams (see `arrays-sync-flow.typ` in this repo). -- Anything where labels should render in the paper's math/text fonts. -- Anything that benefits from recompiling with the paper. - -## When to use something else - -- **Data-heavy plots** (scalar fields, colormaps, quiver plots, anything with - dense per-pixel or per-cell data) — generate SVG from Python/matplotlib, - include via `#image("foo.svg")`. cetz will fight you here. -- **Geometry computation** (Delaunay, intersections, interpolation) — do it in - Python offline, emit JSON with the shape `{"vertices": [...], - "triangles": [...], ...}`, let Typst just draw. - -## When NOT to use TikZ - -Evaluated and rejected: -- Slower compile than Typst. -- Drags in a LaTeX toolchain that isn't otherwise required by the project. -- No Typst math-font advantage over cetz for labels in our paper context. - -Keep TikZ in back pocket only if a co-author insists on TikZ source. - -## Before you draw (thesis-first discipline) - -When a user hands you a figure request — especially one replacing an -existing ASCII sketch, whiteboard photo, or reference figure — **do not -start transcribing**. The source artefact is a hypothesis about what -to communicate, not a specification. Before opening cetz: - -1. **State the thesis in one sentence.** What is the figure arguing? - If you can't write it plainly, you don't yet understand the figure. - Ask the user to articulate it. - -2. **Honestly audit the source.** If the original ASCII / sketch is - being replaced, it's often replaced *because it doesn't land well*. - Say what's broken about it before proposing the replacement — the - user often agrees and the new figure can do more than the old. - -3. **Enumerate design decisions as explicit questions, not assumptions.** - For a curved-boundary normals figure that's typically: - - Geometry (circle / ellipse / arc span / zoom level) - - Sampling density (number of facets, quadrature points per facet) - - Overlay vs. side-by-side - - Whether to show error quantitatively (arcs, annotations) or leave - it as a visible angle - - Where the figure lives in the repo (which doc / which branch) - Present a proposed interpretation with the decisions flagged; let - the user resolve them before you compile. - -4. **Only then open cetz.** Iterate visually — the discipline above is - about not committing to a design prematurely, not about planning - exhaustively. Once you start, compile often. - -The user's phrase "I'm not quite sure what this is intended to -illustrate" is the canonical trigger for this discipline. If you hear -it (or catch yourself about to transcribe without checking), stop and -do the four steps above. - -## Key gotchas (hit during iteration — not hypothetical) - -1. **Don't give a helper a parameter named after a `cetz.draw` export** - — `anchor`, `fill`, `stroke`. The cause is the `import cetz.draw: *` the - helper needs (gotcha 6): it runs *inside* the function body and shadows the - parameter, so the parameter name resolves to cetz's function rather than to - the value you passed. `anchor` panics with `"Unknown anchor 'anchor' for - element 'none'"`; `fill` gives `"expected color, gradient, tiling, or none, - found function"` pointing into `canvas.typ`, nowhere near your code. Rename - to `align-to`, `bg`, `edge`. See `cetz-cheatsheet.md`. - -2. **Clipping is a Typst concern, not a cetz one.** cetz has no `\clip`. - Wrap the canvas in `#box(clip: true, width: ..., height: ..., ...)` and - draw slightly oversized inside the canvas — the box clips the overflow. - -3. **Painter's algorithm — order matters.** Draw the background first, the - highlight second. No z-index exists. Verified in `mesh-demo.typ`. - -4. **Semi-transparency via `rgb(r, g, b, a)`** (alpha 0–255) or - `color.transparentize(col, 50%)`. Both work for fill and stroke. - -5. **Math in labels just works.** `content(pos, $v_1$)` renders in the - document math font. No escape hatch needed. This is a real cetz win over - SVG. - -6. **`import cetz.draw: *` inside the canvas closure.** Without it, `line`, - `circle`, `content` aren't in scope. Helper functions that draw need - their own `import cetz.draw: *` line inside. - -## Project layout pattern - -Each blog post or paper section gets its own subdirectory under -`figures/`, so a post's figures travel together: - -``` -publications/blog-posts/figures/ -└── / - ├── .typ # cetz drawing - ├── .png # committed output - ├── -data.json # (optional) precomputed geometry - └── generate-.py # (optional) Python that writes the JSON -``` - -Concrete example: `publications/blog-posts/figures/finding-particles/` -holds `mesh-demo.*` and `domain-demo.*` for the post -`finding-particles.md`. - -The JSON intermediate is the forward bridge to underworld3 — see -`underworld-bridge.md`. - -## Reference files - -- `cetz-cheatsheet.md` — what worked from memory vs. needed lookup. -- `underworld-bridge.md` — JSON schema for future `uw.meshing` export. -- `examples/` — self-contained copies (each with `.typ`, `.png`, `.json`, - and generator `.py`) of the figures this skill is scaffolded from. - These are snapshots; the live versions may have drifted if a post or - doc was iterated on further. - - `mesh-demo.*` — element-level point-in-cell test. - Live: `publications/blog-posts/figures/finding-particles/`. - - `domain-demo.*` — parallel domain centroid ambiguity. - Live: `publications/blog-posts/figures/finding-particles/`. - - `facet-vs-true-normals.*` — facet normal vs. smooth-surface normal - on a curved boundary. Live: - `docs/advanced/figures/curved-bc/`. - -## Canonical reference in the repo - -`publications/blog-posts/figures/arrays-sync-flow.typ` — prior cetz figure -in the repo, established version (0.3.4) and house style (hex colours, -helper-function pattern). Follow its conventions. diff --git a/.claude/skills/cetz-figures/SKILL.md b/.claude/skills/cetz-figures/SKILL.md new file mode 120000 index 000000000..3ea8ffa2a --- /dev/null +++ b/.claude/skills/cetz-figures/SKILL.md @@ -0,0 +1 @@ +../../../docs/developer/guides/cetz-figures.md \ No newline at end of file diff --git a/.claude/skills/free-surface-convection/SKILL.md b/.claude/skills/free-surface-convection/SKILL.md deleted file mode 100644 index b4e4c4bfb..000000000 --- a/.claude/skills/free-surface-convection/SKILL.md +++ /dev/null @@ -1,224 +0,0 @@ ---- -name: free-surface-convection -description: The Underworld3 free-surface convection method we are hardening — the THREE-NUMBER pointwise topography integrator (held-lid stress equilibrium h_∞ + L-stable exponential relaxation), NOT FSSA. Reach for THIS before touching any free-surface / dynamic-topography convection run, choosing a surface-update scheme, or "stabilising" a free surface. It records the method, why FSSA is explicitly rejected, and the failure modes. ---- - -# free-surface-convection - -The free-surface scheme used in `~/+Simulations/FreeSurface/convection/fs4_compare.py` -and the design doc `docs/developer/design/FREESLIP_DYNAMIC_TOPOGRAPHY_FREESURFACE.md`. -**This is the method we are HARDENING — do not replace it; do not add FSSA.** - -> The 3-number integrator is necessary but NOT sufficient — see **Hardening strategies -> (2026-06)** below for the material-surface advection, tangential topography term, -> free-slip-inner nullspace, and graded/higher-order-mesh fixes that make it actually -> work. Reference impl + diagnostic tools live in `~/+Simulations/FreeSurface/convection/`. - -## ⚠️ NOT FSSA - -The `docs/examples/free_surface/advanced/Annulus*FS.py` examples use **FSSA** -(`add_natural_bc(δt·(Γ·v)Γ/2, "Upper")`). **That is NOT our method.** FSSA buys -stability by adding an implicit surface traction that **UNDER-deforms the surface** -— it trades accuracy for stability. Our scheme is designed to be stable **and** -accurate. If you find yourself adding `FSSA`, an `add_natural_bc` traction on the -free surface, or a `Gamma.dot(v)` stabiliser — STOP, you have the wrong method. -(Those example files are a template for a *different* approach, not this one.) - -## The three-number pointwise integrator (THE method) - -Two Stokes solves per step on the SAME mesh, then a pointwise surface update: - -1. **Free solve** — stress-free top (NO velocity BC on `Upper`; pressure datum is - pinned by the stress-free condition → no pressure nullspace). The surface - normal velocity `u_n` of this solve IS the kinematic rate `ḣ`. -2. **Held-lid solve** — a second Stokes solve with a RIGID free-slip held lid - (`u_n = 0`, via `add_nitsche_bc(0.0, "Upper", local_h=True)` — see [[project_nitsche_local_h_pr275]]) - and a DRIVING-ONLY body force. Its surface normal stress `σ_nn` gives the - equilibrium topography `h_∞ = -(σ_nn - mean)/ρg`. (The free solve forces - `σ_nn = 0`, so the equilibrium MUST come from the held-lid stress.) -3. **Pointwise exp step**, per surface node, from THREE numbers (`h`, `ḣ=u_n`, - `h_∞`): - ``` - γ = ḣ / (h_∞ − h) # local relaxation rate, clamp γ ≥ 0 - h ← h_∞ + (h − h_∞)·exp(−γ·dt) - ``` - L-stable: the step is bounded between `h` and `h_∞`, so it **cannot overshoot** - regardless of a noisy local `γ` (no "drunken sailor"). 1 extra solve/step; - beats RK4 at large dt. NO per-node freeze-clamp (that was the old `relax` bug). -4. The nodal surface increment is **carried inward by a Laplacian diffuser** - (smooth, minimal mesh deformation — NOT full mmpde adaptation), then - `mesh.deform()`. Uniform meshes are fine; adaptivity is NOT required. - -Reference impl: `fs4_compare.py` → `_surface_step`, `_h_inf` (held-lid σ_nn via a -`Projection`), `_surf_un`, `_carry_diffuser`. The free-slip RIGID-top run (no -surface motion) is the reference; the free-surface run is the same driving solve -PLUS this surface update. See [[project_fs4_adaptive_2x2]], -[[project_stress_equilibrium_freesurface]], [[project_freeslip_topo_freesurface]]. - -## Performance - -- Free-slip (rigid top) stagnant-lid runs are FAST. The free-surface cost is the - extra held-lid solve **plus** that moving the surface forces a COLD-START Stokes - each step (can't warm-start across a deformed mesh). -- Use **uniform** meshes for this problem — the diffuser gives minimal deformation, - no mmpde needed. (FMG works on a uniform `refinement=N` hierarchy; scalar solvers - must avoid FMG — PETSc err62, issue underworldcode/underworld3#276.) - -## Hardening strategies (2026-06) — the integrator alone is not enough - -The 3-number integrator moves the surface correctly, but several *other* things must -be right or it runs away / tangles. All implemented in `fs4_compare.py` (flags noted). - -### 1. Material-surface advection — THE key fix (`--advect-velocity`) -The runaway (`u_n` 42→125→285→445, cold lid leaking in, plumes punching through) was -NOT an `h_∞`/BC bug (`h_∞` is verified correct, even in the stagnant FK lid — held-lid -free-slip is the EASY case there). The bug: the surface moves by the L-stable relaxed -rate `ũ_n = Δh/Δt ≤ u_n`, but T was advected with the **stress-free solve velocity** -(surface-normal = full `u_n`). Net material then crosses the surface. A free surface is -a MATERIAL boundary: advect T with a velocity whose surface-normal = `ũ_n`. Modes: -- `consistent` (the right way): a THIRD Stokes solve, same buoyancy, `v·n̂ = ũ_n` - PRESCRIBED at the surface (penalty), tangential stress-free. `ũ_n = (shape_new−shape0)/dt` - = the full ∂h/∂t at fixed θ (correct ALE target). -- `blend`: `α·v_free + (1−α)·v_held`, `α = φ1(γΔt) = (1−e^{−γΔt})/(γΔt)` (the exp-decay - time-average). By Stokes LINEARITY this *is* the prescribed-`ũ_n` solve for UNIFORM α - (and free for FK, which is linear in v). BUT the single mean-α collapse is NOT close - enough once γ varies per surface node — the planform diverges (mode-3 vs mode-1/2), - throughflow ~23 vs ~0.08. Per-node α breaks div-free (∇α·(v_free−v_held)). So the - per-node `consistent` 3rd solve is REQUIRED for structured planforms. -- `free`: advect with stress-free v (the inconsistent baseline — the runaway). - -### 2. Tangential topography advection (`--no-tangent-advect` to disable; default ON) -The pointwise relaxation omits the `v_t·∂_s h` term — a surface rotation/convergence -should carry the topography pattern along the surface; without it you get edge artefacts -where ∂_s h is large (plume-bulge edges). Fix = operator split per step: (1) departure- -point semi-Lagrangian transport of the surface shape in θ by `ω = v_t/r`, then (2) the -L-stable normal relaxation. Lowers throughflow + improves mesh quality. - -### 3. Free-slip inner boundary — rotation nullspace (`--inner freeslip`) -The rigid rotation `[-y,x]` is a velocity nullspace ONLY while the boundary is CIRCULAR. -Once the free surface DEFORMS, do NOT attach it to the held/consistent solves (→ held -22 s/`DIVERGED_LINEAR_SOLVE`, throughflow blow-up). Keep `petsc_use_pressure_nullspace`; -strip the gauge with the exact post-solve projection `_project_out_rotation` on `v`, -`v_cons` (drives advection) AND `v_h` (one consistent non-rotating frame). The undeformed -free-slip *reference* (`--surface freeslip`) is fine WITH the nullspace attached. - -### 4. Graded / higher-order meshes (drive node movement consistently) -- **Surface-ring detection**: tie the tolerance to the FINEST cell - (`0.5·mesh.get_min_radius()`), NOT the nominal `cellsize`. On a gmsh-graded mesh - (`cellSizeOuter`) the old tolerance scoops the first interior ring → a 2%-thick - "surface band" → tangling (looks like the surface "destroying itself"; it isn't — - the diffuser was fed a corrupt surface). BETTER (TODO): build the ring from the DMPlex - `Upper` label (`dm.getLabel("Upper").getStratumIS`), removing the tolerance entirely. -- **Node movement**: the solve velocity is P2; the mesh geometry is P1. Drive `u_n` and - the tangential transport from a P1 length-smoothed `Vector_Projection` of V (`v_p1`), - NOT a point-evaluation of the P2 field. -- **Stress smoothing**: `topo_proj.smoothing_length` = a fixed PHYSICAL length - (`--smooth-length`), not cell-count, so `h_∞` is mesh/order-independent. - -### 5. Cost — there is no acceleration win (don't chase it) -3 Stokes solves/step (free→u_n, held→h_∞, consistent→advect). Warm-start does NOT help -(outer KSP already 1 iter; FMG supplies its own nested guess — measured SLOWER). Blend- -skip rarely fires (α-spread always large). Operator/PC reuse: already reused across -`solve()`s (first 5.5 s setup, steady ~945 ms = irreducible FMG solve; RHS-only resolve -same cost). The per-step cost is the geometric FMG hierarchy REBUILD on the deforming -mesh — intrinsic to moving meshes, "live with it." The UNIFIED-PENALTY single solver -(`penalty·(v·n̂ − V₁·n̂)·n̂`; penalty=0→free, V₁=0→held, V₁=ũ_n→consistent — held & -consistent share the matrix) is the cleanest formulation (no recompile on the constant) -but doesn't cut the irreducible solve. - -### Elastic-plate flexure `h_∞` — IMPLEMENTED (`--flexure-D`) -Generalizes the LOCAL Airy `h_∞ = −σ_nn/ρg` (the D=0 limit) to a flexed plate -`(D ∂_s^4 + ρg) h_∞ = −σ_nn`, solved SPECTRALLY on the ring (serial Fourier — the -feasible substitute for UW3's blocked 1D-manifold FE solve): per mode -`h_∞(m) = −σ_nn(m)/(ρg + D(m/r_o)⁴)`. `D` sets the flexural wavelength `(D/ρg)^{1/4}` -and damps short-wavelength loads — the physically-grounded, mesh-independent length- -smoothing. In `_h_inf` (h_∞-ONLY — the stable form). **PROTOTYPE — amplitude response correct -(stiffer plate → less deflection) but it does NOT low-pass the SURFACE**: filtering `h_∞` only -sets a smooth set-point; the surface still picks up short-wavelength content from the SL -tangential transport + partial relaxation. "Filter every surface number (h, ḣ, h_∞)" was -TRIED and REJECTED — filtering the GEOMETRY `h` injects a spurious smooth-the-mesh motion into -`ũ_n` → flow runs away (Vrms 50→345); filtering `ḣ` alone is stable but elevates Vrms with no -benefit. So making flexure a TRUE surface low-pass without destabilising is OPEN/hard -(`_flex_filter` helper is in place). Examples: `stagnant_lid_mode1_study/{figures/flexure_*.png, -runs/flexure_D*}`. - -### Open / next (not yet done) -- **Label-based surface ring** (replace the radial heuristic with the `Upper` stratum). -- **Flexure D calibration** to a realistic lithospheric flexural wavelength. - -### Diagnostic tools (`~/+Simulations/FreeSurface/convection/stagnant_lid_mode1_study/scripts/`) -- `heldlid_hinf_check.py` — verify `h_∞` via 4 independent free-slip enforcements × Δη sweep -- `stitch_compare.py` — side-by-side montage of per-run dirs (`--dirs a,b,c`) -- `unified_penalty_solver.py` — the one-solver penalty formulation probe -- `resolve_timing_probe.py` — repeated-solve / lag-Jacobian / reuse-PC timing - -## ★★ Body force must be FULL Boussinesq on a deforming mesh (2026-07-26) - -**On any run where the mesh surface actually moves, the body force must retain the -ρ₀ background: `bodyforce = (thermal_buoyancy − rho_0_g)·r̂` (i.e. ρ = ρ₀(1−αΔT)).** -The reduced (driving-only) form has **NO restoring force for surface deformation -anywhere in the momentum system** — `buoyancy_scale` enters only the kinematic target -h∞ = −σ_nn/ρg, which exerts zero force on the flow. Consequences and evidence (all -measured, `~/+Simulations/FreeSurface/annulus_fs_convection/teaching/`): - -- **Relaxation A/B (`relaxation_test.py`)** — imposed 5% mode-4 topography, no thermal - driving: reduced form = surface FROZEN (velocities are round-off); full density = - **262→5 km in 10 steps at the Cathles rate** (measured 2.1e4 vs ρ₀g/(2ηk) = 2.5e4, - within 20% at res 0.06). -- **Convection A/B (`restoring_force_demo.png`, gap-Ra 3e4, ρ₀g 1e6)** — reduced h - grows monotonically without limit; full density *rings* about its supported amplitude - early (damped), then tracks the developing flow ~35% lower with HIGHER Nu (1.93 vs - 1.65). Without the term, soft surfaces run away entirely (measured to 50% of radius). -- **The FreeSurface manager needs NO change** — held/consistent solves inherit - `stokes.bodyforce`; h∞ then self-consistently includes the self-load (a moving target - the integrator follows cleanly — verified in the relaxation test). - -Rules that come with it: -1. **ρ₀g and Ra are COUPLED: αΔT = (thermal buoyancy coeff)/ρ₀g must be ≲ 0.3.** - Softness is not a free knob — the old "soft" runs at ρg=2e5 were αΔT = 1.2–4 (no - such fluid) and their runaway was partly parameter nonsense. Soft surfaces require - weak driving. -2. **Tighten the solver tolerance (~1e-8)**: the hydrostatic RHS dominates the dynamic - signal by 1e3–1e4, and a relative tolerance judges the total. -3. Prefer `snes_type=ksponly` for the isoviscous solves; the huge RHS makes `newtonls` - thrash worse. -4. Driver switch: `fs_convection.py -uw_full 1`. -5. **`FreeSurface(background_buoyancy="analytic")` is REQUIRED with full density** - (98961147): the recovered reaction contains the self-load +h_current and the - reduced-form negation otherwise flips it -> h_inf = -h + drive/rho_g (parks at HALF - equilibrium; -1 eigenvalue = period-2 ringing; steady flow THROUGH the stationary - surface). "analytic" subtracts the geometric height - no extra solve. The CBF - recovery itself is exact (probe lesson: select boundary DOFs by LABEL, never a - radius mask, on a deformed mesh). `background_buoyancy=` = exact two-reaction - reference mode. - -Related transport fact (same campaign): the serial T-blow-up on deforming FS meshes was -the **old-frame SL reach-back** amplifying per loop cycle (issue #423; smaller dt makes -it WORSE) — retired in `b507aca1`; the manager now uses the standard ALE path + clamp + -deform-aware foot restore. The parallel datum defect is #421. - -## Failure modes — symptom → cause - -| Symptom | Cause | -|---|---| -| Surface deforms but `u_n` RUNS AWAY (e.g. u_n 42→125→285→445), cold lid leaks in, plumes punch through | **RESOLVED**: material-surface advection inconsistency — T advected with stress-free `u_n` while surface moves by relaxed `ũ_n`. Fix = `--advect-velocity consistent` (Hardening §1). NOT an h_∞ bug, NOT fixed by FSSA. | -| `held` solve 22 s / `DIVERGED_LINEAR_SOLVE`, throughflow blows up, with `--inner freeslip` | rigid-rotation `[-y,x]` attached as a nullspace on the DEFORMED (non-circular) surface — invalid. Don't attach it on the moving surface; use the post-solve projection (Hardening §3). | -| Graded-mesh surface "destroys itself" (q→0.2, h_max 2% at step 1) | surface-detection tolerance scooped the first interior ring → 2%-thick band, NOT real deformation. Tie tolerance to finest cell (Hardening §4). The diffuser is innocent. | -| Topography GROWS without saturating; flow-through persists even at huge deformation | **missing ρ₀ background** — reduced body force has no surface restoring force (see the Full-Boussinesq section above). Fix = ρ = ρ₀(1−αΔT); check αΔT ≲ 0.3 | -| T leaves [0,1] on a deforming mesh, mesh-locked hot/cold spikes in the squeezed band, worse at SMALLER dt | old-frame SL reach-back amplifying per loop cycle (issue #423) — use the standard ALE path (`old_frame_traceback=False`, the manager default since b507aca1) | -| Stress-free top but surface not updated each step | nothing stops throughflow (the stress-free top is an open boundary unless the integrator moves the surface to track `u_n`) | -| Nu decays when it should be steady (kinematic free surface) | LAG: the SL foot reaches beyond an under-moved surface → cold pump. Fix = the h_∞ relaxation, not more smoothing | -| Surface "mountain" / one-step spike on adaptive mesh | held-lid Nitsche penalty over-stiffened by GLOBAL h; use `local_h=True` (default, PR #275) = `mesh.cell_size()` | - -## Dead ends (already tried — do NOT repeat) - -- **FSSA** signed-traction free-surface: diverges / under-deforms — rejected. -- **High-k post-smoothing** of the surface: the instability is low-m, smoothing - the wrong band. -- Per-node freeze-clamp in the relaxation (the old `relax` fatal bug). - -## Diagnose by - -`h_max` (deflection, as % of r_o), `u_n` / `vhmax` (surface throughflow — should NOT -grow unbounded), `hinf_max` (the equilibrium target), `vrms`, `Nu`. Compare the -free-surface run against the free-slip RIGID-top reference at matched physical time. diff --git a/.claude/skills/free-surface-convection/SKILL.md b/.claude/skills/free-surface-convection/SKILL.md new file mode 120000 index 000000000..7166340ad --- /dev/null +++ b/.claude/skills/free-surface-convection/SKILL.md @@ -0,0 +1 @@ +../../../docs/developer/guides/free-surface-convection.md \ No newline at end of file diff --git a/.claude/skills/nonlinear-solver/SKILL.md b/.claude/skills/nonlinear-solver/SKILL.md deleted file mode 100644 index 297b3df85..000000000 --- a/.claude/skills/nonlinear-solver/SKILL.md +++ /dev/null @@ -1,321 +0,0 @@ ---- -name: nonlinear-solver -description: How to make a hard nonlinear Stokes solve (Drucker-Prager / yield-stress viscoplastic) CONVERGE reliably in Underworld3 the way the working recipe actually does it — automatic warm-start (one Picard step on a cold start) plus a MULTI-SOLVE δ-continuation (constant δ per solve, warm-start the next, sharper δ), the consistent-Newton tangent, and a non-symmetry-safe multigrid smoother. Reach for THIS when a viscoplastic solve stalls / diverges and you are about to hand-tune PETSc options, ramp δ, or "just add a monitor". It carries the CONFIG TRAP LIST — the setup mistakes that each produce a different failure a few steps in — and the one thing you must NOT do (ramp δ inside a single SNES solve). For the yield-law maths and which tangent per model, see `plasticity-solvers`. ---- - -# nonlinear-solver - -The recipe that gets a **hard viscoplastic (Drucker–Prager) Stokes** problem to -converge, and — more importantly — the list of setup mistakes that stop it. The -central lesson from the Spiegelman hard-case study (`η_bg=1e26`, `V=10`): every -failure was a **solver-configuration** error, not a bad Jacobian. If the correct -setup is a minefield for an expert, that is an API regression — so the goal is to -make the correct path the default path. - -Design of record: `docs/developer/design/nonlinear-solver-homotopy-warmstart.md`. -Yield-law maths, tangent-per-model, quadratic-convergence check: `plasticity-solvers`. - ---- - -## The recipe (what actually converges) - -1. **Warm start.** Start the continuation at **large δ**, where the yield surface is - smooth and the problem is easy, and take **one Picard step** into the Newton - basin. One Picard step is defect-correction iteration 1 — contractive, cheap. From - a *warm* iterate, take **no** Picard step (it wastes the good quadratic start). - A cold `v=0` start is safe on its own terms: `ε̇=0` makes `η_pl` infinite, which - the soft-min carries to the viscous branch (see the trap list for the one form - that must be written carefully). - -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. - ---- - -## DO NOT ramp δ inside one SNES solve - -Ramping δ **inside a single SNES solve** (a `SNESSetUpdate` callback that sharpens -the yield surface between Newton iterations) is **proven dead** — it diverges -`DIVERGED_LINEAR_SOLVE` after ~2 iterations and grinds for ~2 hours, **even on the -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*.** (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.) - ---- - -## CONFIG TRAP LIST - -Each of these produces a *different* failure a few steps in — that is why the hard -case felt like whack-a-mole. Check them first. - -| Trap | Symptom | Fix | -|---|---|---| -| **Perfect plasticity's consistent tangent is SINGULAR along the flow**: on the hard-`Min` plastic branch η = τ_y/2ε̇_II, so 2η + 2η′ε̇_II = 0 — the velocity block is symmetric but semi-definite in every yielded cell. (An earlier version of this row blamed *asymmetry*; that is wrong for any η(ε̇_II) law — the rank-one term η′ ε̇⊗ε̇/ε̇_II is symmetric. Pressure-dependent yield adds a non-symmetric v–p coupling, not a non-symmetric velocity block. Corrected 2026-08-26, maintainer review.) | benign while yielded cells are few (the viscous neighbours regularise); with a large yielded fraction the velocity sub-solve caps out and Newton stalls at ~1e-3, no failure reason | give the plastic branch a positive tangent: a small δ soft-min (`yield_mode="softmin"`, powermean, `yield_anchor="yield"`), a rounded viscosity floor, or rate-strengthening ξ; Picard converges regardless (full 2η stiffness) but is linear-rate. The FMG bundle's `gmres`+`sor` smoother is Newton-safe either way | -| `preconditioner="fmg"` (vs explicit `pc_type=mg` + manual mg opts) | outer KSP "converges" in **1 iteration** → no real Newton correction → stall → `DIVERGED_LINE_SEARCH` | use explicit `pc_type=mg` with the smoother opts above; bound the outer KSP (`ksp_max_it`~80) so a hostile step fails fast | -| Cold plastic start `v=0`, or any rigid/unyielded point | `DIVERGED_FNORM_NAN` at iteration 0 | **Not** a div/0: `ε̇=0` gives `η_pl=+inf`, which `Min` and the sqrt soft-min carry correctly to the viscous branch. Only a soft-min form that computes `η_ve·η_pl/(η_ve+η_pl)` breaks (`inf/inf`). Fixed in the power-mean; if you hand-roll a blend, write the harmonic mean as `η_ve/(1+η_ve/η_pl)`. **Do not reach for a strain-rate floor** — it hides this rather than fixing it | -| LU velocity block with all-Dirichlet-ish BC | pressure nullspace singular | attach the Stokes nullspace / avoid a bare LU there | -| Hand-rolled `snes_monitor` to "see what's happening" | you read residuals but miss the tell | use `solve_with_diagnostics` / `get_snes_diagnostics` instead (below) | - -**The diagnostic tell:** `solver.get_snes_diagnostics()["linear_iterations"] ≈ 1` -per Newton step means the linear solve is doing **no real work** (the FMG-1-iteration -trap). A healthy consistent-Newton solve does real Krylov work each step and -converges quadratically. Use `solve_with_diagnostics()`, not a hand-rolled monitor. - ---- - -## Automatic warm-start (Layer 1 — landed) - -`solver.has_solution` is a **public, read-only** status flag: `True` only after a -solve whose SNES converged; reset on a structural rebuild (remesh / adapt / -mesh-mover — the `is_setup=False` hook); kept through coefficient changes (viscosity, -δ, BC values, time step). A **diverged** solve leaves it `False`, so the next solve -auto-cold-starts rather than warming off a corrupted iterate. - -On a **cold** (`zero_init_guess=True`) Stokes solve under the **consistent-Newton -tangent**, a single Picard step is now taken automatically (reusing the existing -`picard=1` machinery). The default (frozen) tangent path is bit-identical. - -```python -stokes.consistent_jacobian = True -stokes.solve() # cold → one automatic Picard step, then Newton -if stokes.has_solution: - ... -``` - ---- - -## Implementation status (this line of work) - -- **Layer 1a — DONE:** `has_solution` + cold consistent-Newton Picard warm-up - (`petsc_generic_snes_solvers.pyx`; test `test_0201`). -- **Layer 1b — DONE:** `zero_init_guess` is tri-state — `None` (default) auto-detects - from `has_solution`, `True` forces fresh, `False` insists on warm. Note warm and cold - agree only to the *convergence tolerance*, not bitwise. -- **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 — 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. 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. - ---- - -## Multigrid depth — how to measure a smoother honestly - -**A two-level hierarchy is a coarse-grid correction, not a V-cycle.** Smoother -comparisons made on one are misleading: the gmres-over-richardson margin measured on -the Spiegelman notch is only 5 % at 3 levels but **25 % at 4** (ρ per V-cycle 0.746 → -0.560), because a deeper cycle applies the smoother on more coarse operators. Judge a -smoother at depth or not at all. - -To get depth without a monster problem, refine a **deliberately ultra-coarse NESTED -base**: `make_notch_mesh.py 1` (492 cells) + uniform `refinement=N` gives 3 levels / -7,872 cells at `N=2` and 4 levels / 31,488 at `N=3` — deeper *and* smaller than the old -2-level 38,580-cell setup. In MG you want the coarsest grid as coarse as it can be -before the problem breaks down. - -- **Never use a non-nested hierarchy** here — it does not give strong MG convergence - (maintainer ruling). Uniform refinement nests by construction. -- Accepted tradeoff: uniform refinement does **not** snap new boundary nodes back to - the analytic notch arcs (no CAD/EGADS model attached), so the corner geometry is - frozen at the coarse mesh's chords on every level. - -Measure with `fmg_contraction_probe.py` (ρ_MG per V-cycle; `<0.5` healthy, `0.8–0.95` -struggling, `≥0.98` hangs) or `smoother_depth_sweep.py` (pays the mesh build + viscous -seed once, sweeps smoothers in-process) in the Spiegelman study. - -**`solve_report` cannot see the smoother.** It records the *Newton* contraction; the -outer KSP is Eisenstat–Walker-collapsed to ~1 iteration/step, so the smoother's work -hides inside the velocity sub-block. Probe the `fieldsplit_velocity_` sub-KSP directly. - -**A smoother will not rescue small ξ.** At the hard corner the failure is operator -conditioning — the coarsest grid cannot represent the viscosity contrast — and at 4 -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). - -## The Schur complement: pair the penalty with FMG, never with GAMG - -**Symptom this is for**: the velocity block's iteration count is rock solid but -the pressure sub-solve wanders into the hundreds and eventually stops -converging. - -**First: it is probably not the pressure block.** `S = -B A^-1 B^T` is applied -*through* the velocity solve, so a velocity solve that exits at its iteration -cap makes the Schur operator inconsistent between applications — and no Krylov -method converges against an operator that moves under it. The pressure block -then caps too, and the outer flounders. Measured on SolCx (eta 1e6, P2-P0disc, -h=1/30), changing **only** `fieldsplit_velocity_ksp_max_it`: - -| velocity cap | sec | outer | pressure/app | velocity/app | -|---|---|---|---|---| -| 200 (default) | 976.0 | 44 | **200.0** | **200.0** | -| 5000 | **25.6** | **2** | **30.0** | 618.0 | - -**38x from a number that is not in the pressure block**, and the velocity error -is identical in both rows. Before tuning the Schur solve, check whether either -block sat at exactly its cap — `solve_report.sub` gives iterations and -applications per block, and a per-application count equal to the cap to the -digit is the tell. - -**Then: the penalty is the lever on the Schur count, and it needs FMG.** -`stokes.penalty = lambda` adds `lambda*mu*(div u)(div v)`, which makes the -eta-scaled mass matrix a better approximation to S. Matched on one mesh -(2592 cells), same discrete solve, only the velocity preconditioner differs: - -| lambda | velocity PC | sec | outer | Schur/app | velocity/app | velocity total | -|---|---|---|---|---|---|---| -| 0 | GAMG | 15.49 | 2 | 125.5 | 94.7 | 24802 | -| 0 | **FMG** | **3.88** | 1 | **59.0** | **8.8** | **546** | -| 10 | GAMG | 20.68 | 7 | 22.3 | **199.9 capped** | 33976 | -| 10 | **FMG** | **3.06** | 1 | **18.0** | **13.5** | **270** | - -- **With FMG, `penalty = 10` improves every axis at once**: 21% faster, Schur - count 3.3x smaller, total velocity work halved. FMG absorbs grad-div - augmentation (8.8 -> 13.5 iterations per application); GAMG does not - (94.7 -> capped). -- **With GAMG, do not use it at all.** The same `penalty = 10` makes the solve - *slower* (15.49 -> 20.68 s), because augmentation is exactly what drives GAMG - into its cap. Uncapping rescues it to 11.03 s but it still needs **833** - iterations per application, and FMG is 3.6x faster on the same mesh. - Feasible is not competitive. - -**The accuracy cost is consistent, so it is safe to pair by default.** The -penalty is grad-div, not a true augmented Lagrangian — `div(P2)` is not inside -`P0`, so the term does not vanish at the discrete solution and it does perturb -the answer. But the perturbation converges away: same rate, and the gap shrinks -under refinement. - -| cells | lambda=0 v err | rate | lambda=10 v err | rate | gap | -|---|---|---|---|---|---| -| 648 | 2.112e-1 | — | 2.327e-1 | — | 1.102 | -| 2592 | 1.266e-1 | 1.67 | 1.376e-1 | 1.69 | 1.087 | -| 10368 | 8.727e-2 | 1.45 | 9.305e-2 | 1.48 | **1.066** | - -For a pressure-dependent constitutive law use the mechanical pressure, -`p_mech = p - lambda*mu*(div u)`; the raw `p` is the multiplier. - -**Traps.** - -- **FMG needs a refined base or you silently get GAMG.** Measured: - `refinement=0` -> one hierarchy level -> default velocity PC is `gamg`; - `refinement=2` -> `mg`. So `penalty` set "with FMG" on an unrefined mesh is - actually the harmful GAMG pairing. Check - `snes.getKSP().getPC().getFieldSplitSubKSP()[0].getPC().getType()`, or read - `solver.pc_fallbacks`. -- **Scaling `saddle_preconditioner` by a constant does nothing** — it does not - change the Krylov subspace. `1/eta` and `101/eta` both give 28 iterations, - identical to every digit, so an "AL-matched" `1/(eta*(1+lambda))` cannot help. - The 1/eta *weighting* itself is worth 1.9x (28 vs 52 with a flat `1`). -- **Eisenstat-Walker is inert under `snes_type=ksponly`** — identical iterations - and error on or off. And `outer 1` is not an EW artefact: it is what a full - Schur factorisation gives when the Schur complement is solved well. -- Measurements: `~/+Simulations/pressure_schur_625/` (#625). - -## Gotchas - -- **`./uw build` → `amr-dev` env**; verify `uw.__file__` is the worktree site-packages. -- **Run VEP/consistent-Newton tests UNFORKED** — `pytest --forked` SIGABRTs (fork of - multithreaded PETSc). -- Benchmark **every** default change — "Solver Stability is Paramount". -- ξ (rate-strengthening) is a **non-homotopic** regularisation: put a user loop - *around* `solve()`, never inside the δ-march. - -## Reference - -- Design: `docs/developer/design/nonlinear-solver-homotopy-warmstart.md`, - `jacobian-consistent-tangent.md`, `solver-strategies-catalogue.md`. -- 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` (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/nonlinear-solver/SKILL.md b/.claude/skills/nonlinear-solver/SKILL.md new file mode 120000 index 000000000..2c57852aa --- /dev/null +++ b/.claude/skills/nonlinear-solver/SKILL.md @@ -0,0 +1 @@ +../../../docs/developer/guides/nonlinear-solver.md \ No newline at end of file diff --git a/.claude/skills/plasticity-solvers/SKILL.md b/.claude/skills/plasticity-solvers/SKILL.md deleted file mode 100644 index d5f83bfd0..000000000 --- a/.claude/skills/plasticity-solvers/SKILL.md +++ /dev/null @@ -1,226 +0,0 @@ ---- -name: plasticity-solvers -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 what the yield campaigns measured actually works — -and records what was retired. - -**The default call is now just:** - -```python -stokes.constitutive_model = cm # ViscoPlastic / ViscoElasticPlastic / TI-VEP -cm.Parameters.yield_stress = tau_y # finite -> plasticity active -stokes.consistent_jacobian = True # Newton tangent (non-elastic DP; see table) -stokes.solve() # cold start takes ONE Picard step automatically -``` - ---- - -## 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. **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 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 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) -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) -# a nonzero difference == the Newton form is present -``` - ---- - -## δ 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. - ---- - -## 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)`. - ---- - -## Failure modes → fixes - -| Symptom | Cause | Fix | -|---------|-------|-----| -| `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) | - ---- - -## Gotchas - -- **`./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). -- `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. - ---- - -## Reference - -- 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 -tangent — the origin of the "~20 iterations is intrinsic" folklore. - -## SNESFAS — do not reach for it - -Nonlinear multigrid (SNESFAS) looks tempting for hard viscoplastic solves but is -**not a viable option** at present (maintainer ruling 2026-07-17): there are no -good preconditioners for the nonlinear hierarchy, and it abandons the robust -linear-solver path (consistent tangent / continuation + fieldsplit + MG) that -this skill is built around. It stays options-only for experiments; treat it as a -future investigation. See `docs/developer/design/solver-strategies-catalogue.md` -and `MULTIGRID_MINIMAL_CONTROL_2026-07.md` (ruling 6). diff --git a/.claude/skills/plasticity-solvers/SKILL.md b/.claude/skills/plasticity-solvers/SKILL.md new file mode 120000 index 000000000..73d85daee --- /dev/null +++ b/.claude/skills/plasticity-solvers/SKILL.md @@ -0,0 +1 @@ +../../../docs/developer/guides/plasticity-solvers.md \ No newline at end of file diff --git a/.claude/skills/uw-visualisation/SKILL.md b/.claude/skills/uw-visualisation/SKILL.md deleted file mode 100644 index 1a43f6ea5..000000000 --- a/.claude/skills/uw-visualisation/SKILL.md +++ /dev/null @@ -1,125 +0,0 @@ ---- -name: uw-visualisation -description: Render Underworld3 mesh fields (T, V, viscosity, the adapted mesh) correctly with PyVista. Use whenever you need to SEE a UW3 result — a field colormap, the moving/adapted mesh, streamlines, or compare runs. Reach for THIS before hand-rolling a renderer; getting the four cosmetic settings wrong makes renders look grey/patchy/blocky and wastes a round-trip with Louis. ---- - -# uw-visualisation - -Canonical PyVista recipe for Underworld3 fields. This exists because every fresh -Claude session re-derives the renderer and gets the colormap / background / -lighting / DOF-sampling wrong, producing "grey/patchy/weird" images Louis -rejects. The settings below match his reference renders exactly. - -**Use PyVista (`underworld3.visualisation`), NOT matplotlib.** Louis reaffirmed -this even after seeing the legacy matplotlib renderer -(`scripts/fault_convection_frames.py`) — that one is NOT preferred. - -## Hard rules (artifacts + output location) - -- **Outputs go under `~/+Simulations/...`, NEVER `/tmp`** (Louis can't view /tmp - or harness task paths). Mirror the run's `--sim-dir`; write `T_.png` into - the run directory, comparison figures into the sim-dir root. -- `pv.OFF_SCREEN = True` at import; finish with `pl.screenshot(path); pl.close()`. - -## The field+mesh pattern (copy this exactly) - -```python -import numpy as np, underworld3 as uw, underworld3.visualisation as vis, pyvista as pv -pv.OFF_SCREEN = True - -mesh = uw.discretisation.Mesh(f"{label}.mesh.00000.h5") # or the live mesh -T = uw.discretisation.MeshVariable("T_v2p1", mesh, 1, degree=3, continuous=True) -T.read_timestep(label, "T_v2p1", 0, outputPath=D) # or use the live var - -pv_T = vis.meshVariable_to_pv_mesh_object(T) # Delaunay through T's OWN DOFs -pv_T.point_data["T"] = np.asarray(T.data[:, 0]) # attach DOF values DIRECTLY (P3-faithful) -edges = vis.mesh_to_pv_mesh(mesh).extract_all_edges() - -pl = pv.Plotter(off_screen=True, window_size=(1000, 1000)) -pl.set_background("white") # rule 2 -pl.add_mesh(pv_T, scalars="T", cmap="RdBu_r", clim=(0, 1), # rules 1 + clim required - show_edges=False, lighting=False) # rule 3 -pl.add_mesh(edges, color="black", line_width=0.5, lighting=False) # mesh overlay -pl.view_xy(); pl.camera.zoom(1.3) -pl.screenshot(out); pl.close() -``` - -## The four things that make renders look bad (all COSMETIC) - -1. `cmap="coolwarm"` → muddy grey-lavender midtone — this IS the "blue/grey/red" - Louis rejects. **Use `cmap="RdBu_r"`** (clean blue→white→red). -2. PyVista's default grey background bleeds through RdBu_r's white (T≈0.5) → dirty - grey. **Always `pl.set_background("white")`.** -3. Default lighting darkens the colormap. **Always `lighting=False`** on every - `add_mesh`. -4. Re-evaluating via `scalar_fn_to_pv_points` / vertex-only sampling drops the - high-order DOFs → blocky. **Attach `T.data[:,0]` directly** to the DOF-cloud - mesh from `meshVariable_to_pv_mesh_object` (it is correct for annulus/box/disc - — do NOT avoid it). `clim` MUST be passed (default `clim=""` trips `np.any`). -5. Resampling ANY field (even P1) onto a regular pixel grid via - `uw.function.evaluate` **dapples at element boundaries** — grid points that - straddle a facet get located into a neighbouring cell with slightly-off - reference coords (Louis: "artefacts across the elements", S-fault rig). - Render derived fields NODALLY on the mesh's own triangulation instead: - evaluate at `mesh_to_pv_mesh(mesh).points` (exact at vertices for P1, - whichever cell the locator picks), attach as point_data, let VTK - interpolate WITHIN elements. On a SPLIT mesh never Delaunay the DOF cloud - (it re-triangulates across the slit) — use the mesh's own cells. - -## Seeing the MESH (adaptation / moving mesh) - -The full-annulus T colormap **washes out mesh detail** — at whole-domain zoom the -grading is invisible. To judge adaptation you MUST crop: - -- Zoom the feature region with a parallel camera: - `pl.camera.parallel_projection = True; pl.camera.parallel_scale = half_width; - pl.camera.focal_point = (cx, cy, 0)`. -- For mesh-only views, drop the field and draw `edges` on white, `line_width≈0.7`. -- Real corruption vs render artifact: apparent "holes / lumps" are often a - mesh-overlay/low-res artifact. Before calling adaptation broken, CHECK the - field's value range is bounded and count folded elements (negative cell area) - programmatically — do NOT diagnose from a render alone. -- **Overlay the feature you're refining to** (a fault trace, an interface): draw it - as a red `pv.PolyData` line over the mesh. Without it you cannot tell whether the - refinement sits ON the feature or has drifted off it (a real failure mode — see - the `adaptive-meshing` skill). Read the geometry from the run manifest so any run - renders the same way. - -## Adaptive / long runs - -- **Render each checkpoint as it lands**, not just the last frame: arm a Monitor that - polls for new `run.mesh.NNNNN.{xdmf,h5}` and emits the index → render on each - event. A completion-only watch leaves you blind for a multi-hour (e.g. TI) run. -- The per-step mesh GEOMETRY must have been written (`write_timestep(..., - meshUpdates=True)`) or you'll render deformed fields on the stale step-0 mesh. - Load the per-step `run.mesh.NNNNN.h5` as the mesh, then `read_timestep` the vars. - -## Velocity - -Same pattern; use **streamlines, not glyphs**. Build a pv mesh for V, add -`pv_mesh.streamlines(...)` or evaluate V on a line seed. Magnitude with the same -white-bg / lighting=False rules. - -## Quantities to judge a convection run (not just pretty pictures) - -- `vrms` from `uw.function.evaluate(V.sym.dot(V.sym), mesh.X.coords)` → the clean - kinetic-energy indicator (more reliable than nodal boundary metrics). -- Surface heat flux Nu via `uw.maths.BdIntegral` on the Upper boundary. -- Mesh quality: fault/bulk nearest-neighbour spacing RATIO (cKDTree) for refinement; - folded-element count + min cell area for tangling. - -## Templates in this skill - -- `render_field.py` — single/`--all`-steps T+mesh render of a run directory. -- `render_field_streamlines.py` — T colormap + mesh + **V streamlines** (sparse - seeds, thin lines, short integration so weak/closed cells read clearly, not - black spiral-blobs). Use for convection. `--tag --all`. -- `zoom_compare.py` — side-by-side cropped mesh+field for N runs at one step. - -Copy these into the run's `scripts/` (or run in place), point `--sim-dir` at the -run, and adjust the field/variable names. They already encode every rule above. - -## Related memory - -`feedback_use_uw_pyvista_visualisation.md`, `feedback_pyvista_viz_pattern.md`, -`feedback_render_all_steps.md`, `project_adaptation_corruption_was_render_artifact.md`. diff --git a/.claude/skills/uw-visualisation/SKILL.md b/.claude/skills/uw-visualisation/SKILL.md new file mode 120000 index 000000000..b0ba15a52 --- /dev/null +++ b/.claude/skills/uw-visualisation/SKILL.md @@ -0,0 +1 @@ +../../../docs/developer/guides/uw-visualisation.md \ No newline at end of file diff --git a/CLAUDE.md b/CLAUDE.md index 9f53ad445..505231528 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -103,6 +103,14 @@ they are about (`mesh-adaptation-architecture.md`), never whimsical ones. ## Rulings a session needs in hand +**Capability guides are documentation, and a change to a family updates +them.** Curated guidance lives in `docs/developer/guides/` with front matter +naming the families it applies to (see the developer index, "Capability +guides"); `uw.capabilities()` and the class-level `view()` list a guide beside +its family, and the skills in `.claude/skills` are symlinks to those pages. +Touching a solver, constitutive model, history scheme or boundary mechanism +means checking every guide that names it, in the same change. + **Adversarial review before the PR opens.** Every branch gets one before it becomes a PR, and again after any substantial post-review commit. The checklist — parallel rank asymmetry, frame and unit boundaries, determinism, tests that cannot diff --git a/docs/developer/guides/adapt-on-top-faults.md b/docs/developer/guides/adapt-on-top-faults.md new file mode 100644 index 000000000..43b130479 --- /dev/null +++ b/docs/developer/guides/adapt-on-top-faults.md @@ -0,0 +1,371 @@ +--- +name: adapt-on-top-faults +description: Recipe for Underworld3 FAULT models on an NVB adapt-on-top mesh — resolve a fault Surface by LOCAL refinement (mesh.adapt(metric, max_levels=...) returns a child; NVB is the 2D default engine), drive it from the fault's EXACT signed distance, use rotated strong free-slip (composes with transverse-isotropy where Nitsche does not), recover dynamic topography from the constraint reaction, and run advection-diffusion on the adapted mesh with field transfer across re-adaptation. Reach for THIS for instantaneous/coupled fault-flow problems. For MMPDE node-movement convection use the `adaptive-meshing` skill instead; for rendering use `uw-visualisation`. +families: [Stokes] +kind: recipe +--- + +# adapt-on-top-faults + +The validated recipe for **fault problems on a locally-refined (adapt-on-top) mesh**. +Distilled from the annulus fault study (2026-07, `feature/adapt-on-top`). + +**This is the REFINEMENT paradigm**, not the mover one: +- `mesh.adapt(metric, max_levels=...)` bisects the base finest **locally** and returns + a **new child mesh** (`child.parent is mesh`). It is *adapt / re-adapt*, NOT node + movement — non-cumulative (each call re-marks from the static base). The child owns + a custom-P geometric-MG (FMG) tail so solvers on it get multigrid for free. +- For **MMPDE / equidistribution node movement** (deforming the same mesh to a field) + use the **`adaptive-meshing`** skill instead. Different tool; don't mix them up. + +Reference implementations (copy from these — all run np1/2/4): +`~/+Simulations/nvb_parallel_fault_study/` (weak-fault Stokes + FMG), +`~/+Simulations/shear_box_fault_study/` (iso vs TI, orientation sweeps), +`~/+Simulations/annulus_fault_study/` (rotated free-slip + topography + moving fault + +advection-diffusion). Companion skills: `uw-visualisation`, `adaptive-meshing`. + +--- + +## The core loop (fault → metric → adapt → child) + +```python +import underworld3 as uw, numpy as np, sympy + +# Base mesh MUST be built with refinement>=1 (supplies the coarse MG tail that NVB +# extends). Base cellSize ~ 2x the target h_near is the SWEET SPOT (see gotchas). +base = uw.meshing.Annulus(radiusInner=0.55, radiusOuter=1.0, cellSize=0.08, + refinement=1, qdegree=3) # or UnstructuredSimplexBox(..., refinement=2) + +# fault as a Surface (polyline control points, N x 3 with z=0 in 2D) +fault = uw.meshing.Surface("fault", base, fault_pts, symbol="F") +fault.discretize() + +# METRIC = a CALLABLE built from the fault's EXACT signed distance. This is the key +# to clean, non-patchy grading: it is evaluated at each refined level's centroids, +# so it resolves itself at the new resolution (no P1-field aliasing). +metric = fault.refinement_metric_function(h_near=0.02, h_far=0.08, width=0.05, + profile="linear") +child = base.adapt(metric, max_levels=3) # -> graded child (NVB is the 2D default) +``` + +- NVB (the 2D default engine) = graded newest-vertex bisection (bounded closure, + parallel via the native `uwnvb` transform; bit-confluent serial↔parallel); + `engine=` is the advanced selector. `engine="sbr"` = uniform + patch (the default; not graded). NVB is 2D only for now. +- `max_levels` is the isotropic-equivalent depth (NVB runs `2*max_levels` bisection + passes). The metric shape decides the grading; `max_levels` just caps it. + +### Metric options (all accepted by `adapt`) +1. **callable** `metric(centroids)->M` — **preferred for faults**. Evaluated per level. +2. MeshVariable / sympy expression — sampled via `uw.function.evaluate` from the BASE + mesh → a peaked `M=1/h²` aliases → *patchy* levels. Avoid for thin features. + +**Custom refinement shape** — pass any callable. For a *fat, uniformly-fine* band +(not just a thin line at the fault), a flat-core metric: +```python +def metric(pts, _f=fault): + d = _f.unsigned_distance(pts) # EXACT distance at arbitrary points + core, ramp, hn, hf = 0.02, 0.05, 0.01, 0.08 + h = np.where(d < core, hn, np.minimum(hn + (hf-hn)*(d-core)/ramp, hf)) + return 1.0 / h**2 +``` + +`Surface` distance API (all exact, arbitrary query points): +`fault.signed_distance(coords)`, `fault.unsigned_distance(coords)`, +`fault.director` (unit normal = normalised ∇(signed distance); the TI weak-plane +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 +needs the fault's distance/normal ON THE CHILD. Two ways: + +```python +# (A) re-home the SAME fault onto the child (cleanest; distance recomputes on child) +fault.remap_to(child) +eta = fault.influence_function(width=0.04, value_near=1e-3, value_far=1.0, + profile="smoothstep") # isotropic weak zone + +# (B) or build a child-side Surface (needed if the base fault is still in use for a +# base-mesh metric — symbol disambiguation refuses a base-mesh symbol in a child +# solver). fault_c = uw.meshing.Surface("fault_c", child, fault_pts); fault_c.discretize() +``` + +Isotropic weak zone: +```python +stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel +stokes.constitutive_model.Parameters.shear_viscosity_0 = eta # drops near fault +``` + +Transverse-isotropic weak PLANE (the physical fault; low fault-parallel shear): +```python +stokes.constitutive_model = uw.constitutive_models.TransverseIsotropicFlowModel +stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_bulk # normal viscosity (constant) +stokes.constitutive_model.Parameters.shear_viscosity_1 = eta # weak near fault -> bulk far +stokes.constitutive_model.Parameters.director = fault.director # unit fault normal +``` +The **TI Jacobian is the full consistent tangent** (not isotropic + defect +correction — that framing is WRONG). The TI velocity FMG V-cycle needs a few more +iters than iso (~8 vs ~2) because of a directional near-null mode the isotropic +point-smoother doesn't damp — bounded and contrast-independent, not a bug. TI needs +~2 elements across the weak zone to resolve (iso ~1). + +--- + +## Rotated strong free-slip (the reason to use this, not Nitsche) + +Nitsche free-slip is INCOMPATIBLE with the TI model (its penalty scales by an +isotropic viscosity). Rotated strong free-slip imposes `u·n̂=0` as an ESSENTIAL +constraint in a per-node (n,t) frame → machine-zero leakage AND composes with TI. + +```python +# normal=None (the default) is measure-weighted and consistent with the assembly — +# prefer it. An analytic nhat is exact for the TRUE circle but keeps a consistency +# error against the faceted integral (#560); use it only when the constraint must +# follow the geometry rather than the mesh. +nhat = mesh.CoordinateSystem.unit_e_0 # exact radial normal (annulus/sphere) +stokes.add_rotated_freeslip_bc(0, "Upper", normal=nhat) +stokes.add_rotated_freeslip_bc(0, "Lower", normal=nhat) +stokes.petsc_use_pressure_nullspace = True # enclosed -> pressure gauge +stokes.solve() # rigid-rotation gauge auto-removed +# Convergence status: read stokes._rotated_freeslip_info = {ksp_reason, +# nonlinear_iterations, rotation_gauge_removed, reaction} — NOT s.snes.getConvergedReason() +# (the rotated solve is a manual loop, not snes.solve). Also sanity-check v·n leakage (~1e-16). +``` + +**Nonlinear rheology / warm-start / timestepping (PR #298, `feature/rotated-snes`):** +rotated free-slip now works *inside* the nonlinear iteration — a nonlinearity probe +auto-dispatches power-law / VEP / TI-with-yield models to a Newton/Picard loop +(`solve_rotated_freeslip_nonlinear`), so warm-started time loops are correct. On the +worktree state *before* #298 lands, the rotated path is a SINGLE linear solve — it +silently returns one Newton linearisation from `u=0` for a nonlinear model. If you +run nonlinear TI + timestepping with rotated free-slip, make sure #298 is in. + +**Dynamic topography** from the constraint reaction (the reason to bother): +```python +h = uw.discretisation.MeshVariable("h", mesh, 1, degree=1) +stokes.dynamic_topography("Upper", h, buoyancy_scale=rho_g) # h = -(σ_nn - mean)/ρg +xs, sig = stokes.boundary_normal_traction("Upper") # or the raw σ_nn (lumped-mass) +``` + +--- + +## FMG under rotated free-slip + +`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 + 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 + ``` + (~4 velocity iters vs ~26 GAMG). Otherwise it falls back to GAMG (fine, just slower). + +The default velocity-block preconditioner and the pressure Schur (`1/η`) are already +near-optimal for TI — do **not** hand-roll a "TI-aware Schur"; measured, the default +`1/η₀` beats every alternative (the weak TI mode is fault-parallel shear, which is +volume-preserving, so pressure sees η₀). + +--- + +## Moving fault + re-adaptation + field transfer + +`adapt()` is non-cumulative and deterministic. Move the fault, re-adapt from the +static base, carry any field by interpolation: + +```python +for step in range(N): + fault_pts = move(fault_pts, step) # kinematics + fault = uw.meshing.Surface(f"fault{step}", base, fault_pts, symbol="F"); fault.discretize() + child = base.adapt(fault.refinement_metric_function(...), max_levels=3) + # verify: folded=0 (all cell |vol|>0), base unchanged (non-cumulative) + # carry a field child_{k-1} -> child_k by interpolation: + T = uw.discretisation.MeshVariable(f"T{step}", child, 1, degree=1) + if T_prev is None: + T.data[:,0] = uw.function.evaluate(T0_expr, T.coords) + else: + T.data[:,0] = uw.function.evaluate(T_prev.sym, T.coords) # mesh->mesh interp + T_prev = T +``` +Transfer error is **non-accumulating** for a smooth field (bounded by the per-mesh P1 +representation floor, ~0.2% on a fine base, ~3% on a coarse base) — repeated +re-meshing does not diffuse a smooth field away. Sharp fronts lose more per transfer. + +--- + +## Advection-diffusion on the adapted mesh + +```python +T = uw.discretisation.MeshVariable("T", child, 1, degree=2) +adv = uw.systems.AdvDiffusionSLCN(child, u_Field=T, V_fn=stokes.u.sym, order=1, + monotone_mode="clamp") # clamp = bounded, no overshoot +adv.constitutive_model = uw.constitutive_models.DiffusionModel +adv.constitutive_model.Parameters.diffusivity = 2e-4 +adv.f = 0.0 +dt = 0.015 # FIXED dt — SLCN is semi-Lagrangian (unconditionally + # stable). estimate_dt() reports the tiny fault-band + # Courant limit; do NOT use it to size the step. +for step in range(nsteps): + adv.solve(timestep=dt) +``` +- **SLCN dt is NOT Courant-limited** by the fine fault-band cells — pick dt from the + coarse-region advection, verify accuracy. +- The scalar AD auto-FMG-injection bug (velocity-block custom-P mismatching a scalar + operator on adapt children → PtAP error 60) is **FIXED**: `auto_inject_custom_mg` + now calls `snes.setUp()` before reading the finest reduced map (so the DM section + is the finalized space the operator lives on) and validates that map against the + assembled operator. Custom-P installs successfully on scalar SLCN AdvDiffusion on + NVB adapt children (no skip-guard, no workaround) — `bugfix/custom-mg-parallel`. + +--- + +## 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 | +|---|---| +| patchy along-fault refinement (level 4 here, 2 there) | P1-interpolated `M=1/h²` aliasing. Use the **callable exact-distance** metric. | +| child solver rejects `fault.distance.sym` (foreign-mesh error) | symbol disambiguation. `fault.remap_to(child)` or build a child-side `Surface`. **Rough edge**: needing two Surface objects. | +| adapt added-cell count jitters ±30% step-to-step | small-number geometry on a thin band. Deterministic; quality (near-fault h) is constant. Base ≈ **2× target** minimises it (CoV ~10% at 0.08 base vs 19% at 0.14, 24% at 0.05). | +| tiny AD timestep / slow advection | `estimate_dt` returns the fault-band Courant limit. Use a fixed dt (SLCN is unconditionally stable). **Rough edge**: `estimate_dt` not adapt-aware. | +| 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 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. | +| scattered 1-level refinement across the WHOLE domain; far-field quality drops | the metric's far clip ``h_far`` sits below the base mesh's cell DIAMETERS (gmsh ``cellSize`` is a target edge length; diameters run 1.2–2.5x it — measured 0.108–0.223 for cellSize 0.18+ref 1). ``edge_split`` marks on diameter, so 100% of the domain refines once and the unrequested bisection DE-CONDITIONS the grid (far-field median q 0.372 -> 0.295, measured). Set ``h_far >= 1.05 * cell_diameters(base.dm).max()``. | +| 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/docs/developer/guides/adaptive-meshing.md b/docs/developer/guides/adaptive-meshing.md new file mode 100644 index 000000000..4738d3b37 --- /dev/null +++ b/docs/developer/guides/adaptive-meshing.md @@ -0,0 +1,407 @@ +--- +name: adaptive-meshing +description: The canonical, workable recipe for Underworld3 moving-mesh / adaptive-mesh convection (annulus stagnant-lid, faults, free surface). Reach for THIS first when setting up any model with a deforming or adapted mesh — it encodes the combination that does not blow up, tangle, or inject spurious energy, and explains the failure modes so you don't re-derive them. Use before choosing movers, free-slip BCs, restart, or field-transfer options. +families: [AdvDiffusion, Stokes] +kind: recipe +--- + +# adaptive-meshing + +The one workable combination for UW3 moving/adaptive-mesh convection, distilled +from many sessions that each re-picked options and stomped on each other's +defaults. **Start from this recipe; change one thing at a time and verify.** + +Reference implementation (current, validated): the **`underworld3.workflows` +adaptive-convection example** — +`docs/examples/workflows/adaptive_convection/` on the `feature/adaptive-convection` +worktree (`config.py`+`simulate.py` no-fault; `fault_config.py`+`fault_simulate.py` +fault; `diagnostics.py`, `render.py`, `compare.py`). Express adaptive runs as a +WORKFLOW (a `WorkflowConfig` + `@workflow_step` DAG + `Run`), NOT a monolithic +driver. The older `scripts/fault_convection_adapt_loop.py` (feature/fault-convection) +is superseded — its ideas are folded into the workflow + this skill. +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(...)` 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. 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. + +--- + +## Mover quick-start (copy-paste — this is the hard-to-discover bit) + +The user entry is `uw.meshing.node_redistribution(mesh, metric, ...)` (the +purposeful spelling; it dispatches to `mesh.redistribute_nodes`, which drives +the MMPDE mover on 2D simplex meshes — `smooth_mesh_interior` is the +machinery underneath and takes the same kwargs). Minimal correct setup to +adapt a mesh to a field `T` each step: + +```python +import underworld3 as uw + +# metric from |grad T|: refinement=R is a factor on the BACKGROUND spacing h0, +# not a finest:coarsest ratio. The envelope is h in [h0/R, h0*coarsening], and +# coarsening="auto" is R**(1/d) — so R=5 in 2-D spans h0/5 to 2.2*h0, a ratio +# of R**(1+1/d) ~ 11. Use refinement=R, NOT strategy= (caps at ~2, under-grades). +rho = uw.meshing.metric_density_from_gradient( + mesh, T, refinement=5, coarsening="auto", metric_choice="front-following") + +# move the mesh — the mover (Huang-Kamenski MMPDE) is variational, +# non-folding, clusters AND aligns cells. It OWNS field transfer +# (remaps T + SLCN history, fires on_remesh hooks). +uw.meshing.node_redistribution( + mesh, rho, + method_kwargs=dict(step_frac=0.2, accel="cg", momentum=0.0), # mmpde's OWN kwargs + slip_surfaces=True, # boundary nodes slide tangentially (parallel-safe) + skip_threshold=0.9) # skip the move when the mesh is already aligned +``` + +For a **sharp feature (fault)** pass an anisotropic SPD TENSOR metric instead of +the scalar `rho` (thin ACROSS the feature normal n) and bake a gmsh base: + +```python +import sympy +n = sympy.Matrix([nx, ny]) # constant fault-normal unit vector +d = dfac.sym[0] # DIRECT unsigned distance field (P1) +M = rho * sympy.eye(2) + (Rf**2 - 1.0) * sympy.exp(-(d/w)**2) * (n * n.T) +# mesh built with: uw.meshing.Annulus(..., refine_lines=[xy], refine_size_min=smin) +uw.meshing.node_redistribution(mesh, M, + method_kwargs=dict(step_frac=0.2, accel="cg", momentum=0.0), + slip_surfaces=True, skip_threshold=None) # tensor metric: do the skip check yourself +``` + +Pitfalls that make it "not work": `method="anisotropic"`/`"ot"`/`"spring"`/`"ma"` +(RETIRED 2026-07 — they now raise ValueError; mmpde is the default and only +metric mover); injecting `relax`/`n_outer` (starves mmpde's CG); `strategy=` instead of +`refinement=R` (under-grades); a scalar bump for a fault (refines a fat corridor, +leaves the centre coarse); signed `Surface.distance.sym` for `d` (bleeds along the +line extension — use a direct unsigned distance). Full rationale + the rest of the +recipe (BCs, restart, field transfer, cadence) below. + +--- + +## The canonical recipe (defaults that work) + +### 1. Mover — mmpde +`uw.meshing.smooth_mesh_interior(mesh, metric=..., method="mmpde", +method_kwargs=dict(step_frac=0.2, accel="cg", momentum=0.0), slip_surfaces=True)`. +- mmpde = Huang–Kamenski variational, **non-folding** (energy → ∞ as detJ → 0), + clusters AND aligns cells. It is the only clean mover. +- **NOT** the `anisotropic` mover (shreds/freezes on static features) or `OT` + (slivers — OT is optimal transport, sliver-prone, not a route around the cap). +- Do **not** inject `relax`/`n_outer` into mmpde (those are the anisotropic + mover's knobs and starve mmpde's internal CG). + +### 2. Metric +- Thermal: `metric_density_from_gradient(mesh, T, refinement=R, + metric_choice="front-following")`. `refinement=R` (≈5) is the maximum local + refinement **on the background cell size h0**, not the finest:coarsest ratio: + the metric targets `h ∈ [h0/R, h0·coarsening]`, and `coarsening="auto"` takes + the budget-conserving `R**(1/d)`. So R=5 in 2-D asks for h0/5 up to 2.2·h0 — + a finest:coarsest ratio of `R**(1+1/d)` ≈ 11, and ≈ 8.5 in 3-D. Named + `strategy=` caps at ~2 and under-grades. R≈5 extracts ~all the grading the + node budget/layout allows; don't over-tune R (benign no-op above budget). + Passing `refinement` takes the **envelope branch**, which ignores `amp`, + `lo/hi_percentile`, `mode` and `power`. +- Fault / sharp feature: a **hand-built anisotropic SPD tensor** + `M = ρ·I + (Rf²−1)·exp(−(d/w)²)·n nᵀ` (thin ACROSS the feature normal n). A + scalar bump refines a fat isotropic corridor and leaves the centre-line coarse. +- Use a DIRECT unsigned distance field (geometry_tools) for `d`, NOT the Surface's + signed `.distance.sym` (its zero-contour bleeds the metric along the line + extension). + +### 3. Creation vs maintenance (the cap) — and the gmsh base COMPOUNDS +mmpde **cannot create** strong refinement from a uniform mesh — it saturates at a +fixed-topology cap (~1.8× on the fault, re-measured), because the mover has a FIXED +node budget (it redistributes, never adds nodes). To go finer you need MORE NODES, +which only gmsh can add at construction. **Bake refinement into the gmsh base** +(`Annulus(refine_lines=[xy], refine_size_min=...)` — now real, see the Faults +section) and the mover doesn't just MAINTAIN it: the extra gmsh nodes lift it off +its budget cap so it **compounds** (gmsh f2 base 0.44 → mover 0.29; f3 base 0.30 → +mover 0.19 ≈ 5× finer). Measured (Ra1e6 Rf8 res24, all folded=0): +uniform 0.55 (~1.8×) → gmsh-f2 0.29 (~3.4×) → gmsh-f3 0.19 (~5×). Judge by +fault/bulk nearest-neighbour spacing RATIO, never by global misalignment. + +### 4. Cadence — adapt as often as you like (forced every step is FINE) +Adapt every step or every few — on a CORRECT build (see §5) the mover converges +dead-flat and forced every-step adaptation under vigorous convection is stable +(validated: 39/39 forced adapts, mesh folded=0, area-ratio ~14 flat). Skipping +when aligned just saves cost. (An earlier claim that forcing adaptation +"tangles / over-injects energy" was WRONG — that was the deformed-eval bug in §5.) + +### 5. THE bug that wrecked adaptation (holes) — and the fix +SYMPTOM: giant empty cells / holes in the adapted mesh, intermittent area-ratio +spikes, convection wrecked. ROOT CAUSE (proven): `uw.function.evaluate` +**mis-locates points on a deformed mesh** (the nav kd-tree `mesh._nav_coords` was +captured from the ORIGINAL coords and never refreshed), so the metric — built with +strictly-positive nodal values — evaluates to **NEGATIVE garbage** (even at its own +DOFs) → non-SPD → the mover wrecks the mesh. FIX (both): +- **`8a9d2ff2`** (refresh `_nav_coords` + projected normals on every deform) — + the real fix; makes `function.evaluate`/`points_in_domain` track deformation. +- **Monotone RBF metric bake** (in `_mmpde_mover`, formerly `_winslow_mmpde`): Shepard-interpolate the + metric from its **positive nodal values** — a convex average is guaranteed ≥0 + (monotone) + fast (no cell-location). Use RBF for the metric; it doesn't need + high-precision eval. +With these, the metric stays positive/SPD; #259's SPD-floor never fires (harmless). +The mover, `accel`, and `refinement=R` (e.g. R=5) are all FINE — they were +red-herring symptoms of the eval bug. Always confirm a CLEAN BUILD first +(`./uw build`; `md5 site-packages/.../smoothing.py == src`); a stale build is its +own cause of holes. + +### 6. Stokes free-slip +Penalty (`add_natural_bc(KFS·v·n·n)`, KFS≈1e6) or Nitsche γ=10 both work on a +UNIFORM mesh. Nitsche preferred for sharp fault corners. Do NOT diagnose free-slip +with nodal v·n: Nitsche enforces v·n=0 weakly, so nodal v·n is large even when +correct — use vrms (kinetic energy). (An earlier "warm-restart × Nitsche → blow-up, +use cold restart" diagnosis was WRONG / confounded by the §5 eval bug + a stale build.) + +**On an ADAPTIVE mesh, Nitsche free-slip γ=10 is TOO SOFT → intermittent vrms +spikes at adapt steps** (under-enforcement, NOT a blow-up: single-step velocity +garbage that T survives because the solve is cold-restarted — e.g. vrms +144→8887→144). The mover pins the boundary and refines just beneath it, creating +high-aspect-ratio near-boundary cells whose Nitsche inverse-estimate constant needs +a bigger penalty. Two fixes, BOTH good (corroborated across sessions): +- **Nitsche γ=100** (not 10) for a free-slip top on an adaptive mesh — clean, smooth + signal. This is INDEPENDENT of the penalty's mesh-size scaling: local vs global `h` + are equivalent here (both garbage at γ=10, both clean at γ=100). +- **Don't fully pin the boundary in the mover** — let it tangential-slip + (`slip_surfaces=True`, the §1 default), NOT `pinned_labels=[Upper,...]`. Avoiding + the distorted near-boundary cells lets γ=10 work. Pin a boundary only when you must + hold a prescribed shape (e.g. a free-surface height the integrator just set). + +Aside — free-SURFACE held-lid is the OPPOSITE problem: the GLOBAL `h = +get_min_radius` drifts BELOW the surface cell as the interior refines → the Nitsche +penalty OVER-stiffens → a spurious one-step surface "mountain" (vhmax spike). Fix = +a LOCAL per-cell `h` via `mesh.cell_size()` (deformation/adaptation-tracking); +`add_nitsche_bc(..., local_h=True)` is the default (PR #275). So held-lid wants LESS +penalty (local-h), free-slip-adaptive wants MORE (γ=100) — different knobs. + +### 7. Advection + timestep +`AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V.sym, theta=1.0, monotone_mode='clamp')`. +theta=1.0 (backward Euler) for stability; monotone clamp bounds SL overshoot; +`V_fn=V.sym` is the PHYSICAL velocity. **dt can be LARGE — SLCN is unconditionally +stable; the smallest/median adapted cell does NOT have to govern dt.** Use +`estimate_dt(percentile=50)` × a multiplier (e.g. dt_mult 3–5) to advance physical +time faster (needed to develop convection in stiff stagnant-lid regimes). + +### 8. Rheology +- **isotropic** (linear) ⇒ `snes_type=ksponly` (one exact KSP solve; default + newtonls rejects steps on vigorous flow). Cheap. Use for resolved features. +- **TI / anisotropic** weak fault (`TransverseIsotropicFlowModel`: + `shear_viscosity_0=η_FK`, `shear_viscosity_1=η_weak`, `director`=fault normal): + keep the **default newtonls** (NOT ksponly — ksponly converges to the WRONG + answer); use **penalty** free-slip (Nitsche trips the anisotropic-Jacobian bug); + **GAMG** (FMG ~7× slower). Measured on the gmsh-resolved + one-sided base + (Ra1e6 Δη1e3): ran clean to t=0.06, folded=0, vrms→20, **~10×** the isotropic + cost per step (GAMG eats the 1000× anisotropic contrast — well under the feared + 20×). The TI fault visibly steers the flow (persistent recirculation at the trace). +- Viscosity-bearing fields are **P0/P1 only** (positivity; higher order overshoots). + FK viscosity `η = exp(θ(1−T))`, θ = ln(Δη). Floor a weak zone, don't multiply → 0. + +### 9. Build discipline +- `./uw build` after ANY source change; verify `uw.__file__` is in site-packages + and (when debugging) that `site-packages/.../smoothing.py` matches `src/...`. A + **stale mover build is a top cause of "giant empty elements / holes"** in the + adapted mesh. +- **NEVER** `pip install -e .` (contaminates all envs). Run from inside the worktree + (`pixi run -e amr-dev`). + +--- + +## Faults — the gmsh-resolved, on-fault, one-sided recipe (validated 2026-06-23) + +The full fault recipe, hard-won. Reference implementation: the +`underworld3.workflows` adaptive-convection example +(`docs/examples/workflows/adaptive_convection/{config,fault_config}.py`, on the +`feature/adaptive-convection` worktree — built on the FIXED mover base). Use the +WORKFLOW system, not a monolithic driver. + +### gmsh line refinement (`refine_lines`) — NOW IMPLEMENTED in `Annulus` +```python +uw.meshing.Annulus(radiusOuter=1, radiusInner=0.5, cellSize=1/24, + refine_lines=[xy], # list of (N,2) polylines (model coords) + refine_size_min=cellSize/3, # cell size ON the line (factor 2-3 is plenty) + refine_dist_min=0.02, refine_dist_max=0.12) # size ramps back to cellSize +``` +A gmsh **Distance + Threshold** field along the polyline; INTERIOR points are +embedded so nodes land ON the line. Backward-compatible (default `None`). It is a +core meshing change → lands as its own small meshing PR, separate from the workflow +example. (Before this session it was only *called* in old scripts and never existed +— don't trust `refine_lines` on any branch but the one carrying this commit.) + +### Keeping refinement ON the fault (the metric-composition trap) +The fault metric is `M = ρ·I + (Rf²−1)·exp(−(d/wₙ)²)·n nᵀ` with the isotropic SIZE +density `ρ`. **Do NOT fuse the fault density into ρ_T by PRODUCT** — the cold +surface thermal BL (ρ_T ~ R^d ~ 20–25) out-competes the fault near its top and the +refinement drifts ABOVE the fault, starving the deep fault ("seems to repel"): +- Use `ρ = max(ρ_T, fault_ρ)` (NOT `ρ_T · fault_ρ`). +- Make `fault_ρ = 1 + amp·gauss` with **amp > ρ_T** (~25 at R=5) so the fault wins + the max along its whole length. +- DIAGNOSE this by comparing the gmsh BASE (step 0, refinement centered on the + fault — correct) vs the DEVELOPED mesh (drifted above) → it's the MOVER's metric, + not gmsh. Render the mesh with the **fault trace overlaid** (`render.py --fault`). + +### One-sided fault influence (the clean control) +Even with max+amp, the symmetric metric DEMANDS both flanks while realized nodes +drift to the hanging wall. Make it one-sided: +- Store the **SIGNED** distance in `dfac` (the gaussians square it, so magnitude is + unchanged — the sign only feeds a `0.5(1+tanh(m·d/w))` gate). Probe the + radially-outward side once to define "upper" regardless of the distance tool's + orientation convention. +- `fault_metric_side` (both/upper/lower) gates the refinement; `fault_rheology_side` + gates the weak zone. **`both=upper`** is the physical recipe: a one-sided + hanging-wall damage zone with the mesh refined on the same side (refinement and + rheology coincide; gmsh-f3 → fault/bulk ~0.19, folded=0). `metric_side=lower` + instead pulls refinement onto the footwall to counter the upward drift. + +### The wedge fill (anti-collision) +The fault pull and the surface-BL pull compete for the coarse cells in the radial +sliver BETWEEN them. `fault_wedge=True` gmsh-fills that wedge (sample radial +segments from each fault point up to the surface, add as a second `refine_lines` +point set) so both pulls have their own budget and merge into one coherent fine +wedge instead of colliding. + +### Weak zone (rheology) — geometric blend + gaussian PEAKED on the fault (2026-06-24) +The weak-zone viscosity blends `η_FK` (background) and `floor` (fault) by the +influence `f`∈[0,1] (P1, positive). Get it right with TWO rules (verify by +reconstructing + rendering the REALIZED η field — `render_fields.py` — and the +combined `η_FK(T)^(1−f)` field; never assume the floor is reached): +- **GEOMETRIC blend `η_weak = η_FK^(1−f)·floor^f`** (NOT arithmetic + `η_FK·(1−f)+floor·f`, whose `(1−f)` term leaks the stiff-lid background through + → η≈8 at f=0.97 in a 1000× lid). Geometric reaches the floor genuinely (η≈1.2 at + f=0.97) AND ties the contrast to the LOCAL η_FK — so the fault automatically + bites hardest where it cuts the cold stiff lid (physically correct), nothing in + the hot interior. +- **`f` must be PEAKED on the fault (gaussian), NOT a TOPHAT block.** A top-hat + makes a uniform weak BLOCK; its sharp edges mean strain follows ∇f (TWO parallel + lines at the block edges, not on the fault), and a one-sided block is offset to + the hanging wall (η_1=1 core sits ABOVE the drawn line). Use a **gaussian + `f=exp(−(d/w)²)`** (peak f=1 ON the fault → η_1=1 on the drawn line, strain + localizes INTO the slot as a single feature). side=both = symmetric; side=upper = + hanging-wall halo (gaussian taper up, sharp footwall recovery — NO halving gate). + Centre the METRIC too (`metric_side=both`) so nodes refine on the line. +THERMAL CONTRAST (Louis's insight, confirmed): even a symmetric gaussian gives a +much bigger η-contrast on the COLD (upper/surface) side than the warm (footwall) +side, because η_FK rises ~1000× toward the surface — so the fault's dynamical +prominence naturally concentrates in the cold lid. This is a feature, not a bug. +Verified (gmsh-f3, gaussian width=0.025, side+metric=both): weak zone centred on +the drawn fault on the adapted mesh, folded=0. TUNE AT STEP 0 (build mesh+fields, +no solve — fast; plot the 1D η_1(d) profile + the 2D field). COST: genuine 1000× +TI contrast ~25–145 s/step (cold-start steps slow, ~25 s once developed; GAMG). +For TI see §8. (`fault_config.py`: `fault_profile=gaussian` (default still tophat — +pass gaussian), geometric blend in `create_solvers`; `render_fields.py` light maps.) + +--- + +## Failure modes — symptom → cause → fix + +| Symptom | Cause | Fix | +|---|---|---| +| Giant empty cells / holes in adapted mesh; intermittent area-ratio spikes | `function.evaluate` mis-locates on deformed mesh → metric → negative/non-SPD (the real bug). OR stale build. | `8a9d2ff2` + monotone RBF metric bake; `./uw build` + verify md5 site-packages==src | +| Decays when it should convect | over-diffusion OR under-resolution OR dt too small to develop | check a resolved arbiter (finer space+time); raise node budget; **larger dt** (dt_mult 3–5) | +| Fault won't refine under convection | mmpde creation cap; field gives no signal in cold lid | gmsh `refine_lines` base — the gmsh nodes let mmpde COMPOUND past the cap (~5×) | +| Refinement drifts ABOVE the fault / "repels", deep fault starved | fault density fused by PRODUCT with ρ_T → thermal BL out-competes the fault near the surface | `ρ = max(ρ_T, fault_ρ)`; `amp > ρ_T` (~25); or one-sided `fault_metric_side`; +`fault_wedge` | +| Refinement on both flanks but you want one side | symmetric (unsigned-distance) metric | signed `dfac` + `fault_metric_side`/`fault_rheology_side` tanh gate (`both=upper` = physical) | +| TI fault solve: ksponly gives wrong answer | ksponly skips the Picard the inexact GAMG inner solve needs | keep default newtonls; penalty free-slip; GAMG (~10× isotropic cost, fine) | +| free-slip ADAPTIVE: vrms spikes at adapt steps (e.g. 144→8887→144), T stays bounded | Nitsche γ=10 too soft on the distorted near-boundary cells the pinned-top+interior-refine creates (under-enforcement) | **Nitsche γ=100** (not 10); OR don't pin the boundary (tangential-slip `slip_surfaces=True`). h-scaling (local vs global) is moot here | +| free-SURFACE held-lid: spurious one-step surface "mountain" / vhmax spike | global `h=get_min_radius` drifts below the surface cell as interior refines → Nitsche OVER-stiff | LOCAL per-cell h: `add_nitsche_bc(local_h=True)` (default, PR #275) = `mesh.cell_size()` | + +**Diagnose by:** vrms (KE), Nu (`BdIntegral` surface flux), fault/bulk NN-spacing +RATIO (cKDTree), folded-element count + min cell area. Compare runs **at matched +physical time t** (dt differs between meshes), not by step number. When unsure +whether behaviour is physical, build a **resolved arbiter** (uniform mesh finer in +BOTH space and time) — if it agrees with one candidate, that's the truth. + +## Diagnostics (in the workflow example, reusable) +`diagnostics.py`: `mesh_quality` (folded / area-ratio / aspect), `nn_spacing_ratios` +(BL + fault/bulk), `NusseltSurface`, `vrms`, `History`. `render.py`: T+mesh+ +streamlines on the `Run` layout, **`--fault`** overlays the fault trace (read from +the run manifest) + **`--focus-fault`** auto-crops on it + **`--mesh-only`** for the +clean mesh, **`--all`** for every frame. `compare.py`/`fault_refine_plot.py`: +matched-physical-time comparison + fault/bulk-ratio time series. +**Rendering long runs as they go:** a completion-only Monitor is NOT enough — arm a +Monitor that POLLS for new `run.mesh.NNNNN.xdmf` checkpoints and emits the index, so +you render each step as it lands. +**Checkpoint an ADAPTIVE mesh with `meshUpdates=True`** (per-step geometry) or the +saved frames pair deformed fields with stale step-0 geometry. + +## The verified canonical command + +Requires the fix (§5): `8a9d2ff2` (deformed-mesh point-location) + the monotone RBF +metric bake in `smoothing.py`. Validated 2026-06-22 at vigorous Ra1e6/Δη1e3 +stagnant-lid: forced adapt EVERY step (39/39), larger dt, → vrms→18, Nu→1.86, +|v|max 61, mesh CLEAN every frame (folded=0, area-ratio ~14 flat), no abort. + +```bash +# no-fault baseline (the workflow CLI; one flag per config field) +pixi run -e amr-dev python docs/examples/workflows/adaptive_convection/simulate.py \ + --output-dir ~/+Simulations//baseline \ + --rayleigh 1e6 --delta-eta 1e3 --cellsize 0.0417 \ + --resolution-ratio 5 --adapt-every 1 --dt-mult 4 --max-steps 80 --max-t 0.06 + +# resolved fault: gmsh base (factor 3) + on-fault + one-sided hanging wall + TI +pixi run -e amr-dev python docs/examples/workflows/adaptive_convection/fault_simulate.py \ + --output-dir ~/+Simulations//fault \ + --rayleigh 1e6 --delta-eta 1e3 --cellsize 0.0417 --resolution-ratio 5 \ + --fault-base-smin 0.0139 --fault-anisotropy 8 \ + --metric-combine max --fault-refine-amp 25 \ + --fault-rheology-side upper --fault-metric-side upper \ + --rheology ti --dt-mult 4 --max-steps 40 --max-t 0.06 +``` + +Key choices, verified: +- **`--resolution-ratio 5`** (R=5) is fine — the mover handles it on a correct build. +- **`--adapt-every 1`** force adapt every step (the strongest mesh test). +- **`--dt-mult 4`** — larger dt for STABILITY is fine (SLCN unconditional); but it + costs transient ACCURACY (over-diffusive backward-Euler DELAYS the convective + onset vs a resolved arbiter — dt×1.5 recovers it). Use small dt_mult for faithful + transients, large to reach quasi-steady fast. +- **`--freeslip penalty`** (default) — REQUIRED for TI (Nitsche trips the + anisotropic-Jacobian bug). It's the raw velocity penalty `kfs·(v·n)n`. +- Fault: `--fault-base-smin` (gmsh resolve), `--metric-combine max` + + `--fault-refine-amp 25` (keep refinement ON the fault), `--fault-*-side` (one-sided), + `--rheology ti` (real fault). Drop the fault flags for the no-fault control. + +Render with `render.py` (`--fault --focus-fault` to see the trace + refinement +coincide; `--mesh-only`; `--all`). Judge the mesh by folded/area-ratio, the physics +by vrms/Nu, ALWAYS at matched physical time. + +## Related memory +`project_adaptive_convection_as_workflow` (THIS session: workflow port, gmsh +refine_lines, on-fault/one-sided/wedge, TI, dt-accuracy), `project_mmpde_holes_real_root_cause`, +`project_fault_refine_fixed_topology_cap`, `project_uw_workflow_landing`, +`feedback_debug_adaptive_solver_method`, `project_fault_convection_working_settings`. diff --git a/docs/developer/guides/adversarial-review.md b/docs/developer/guides/adversarial-review.md index c9e45368b..14a13813f 100644 --- a/docs/developer/guides/adversarial-review.md +++ b/docs/developer/guides/adversarial-review.md @@ -88,6 +88,16 @@ name. An anonymous float collapses into the assembled product and the transcript can only show the number. Examples in `docs/` name their coefficients. +**A family's guide says what the family can do now.** The curated guides in +`docs/developer/guides/` carry front matter naming the families they apply +to, and `uw.capabilities()` lists them beside each family. A change to a +solver family, a constitutive model, a history scheme or a boundary +mechanism is reviewed against every guide that names it: the guide is +updated in the same change, or the review says why it still holds. The AI +skills in `.claude/skills` are symlinks to those pages; +`tests/test_0030_capability_guides.py` fails on a copy, on a guide without +front matter, and on a family name no class carries. + ## Where the reviews live `docs/reviews/[YYYY-MM]/`, indexed by `docs/reviews/README.md`, and posted on diff --git a/docs/developer/guides/boundary-condition-rulings.md b/docs/developer/guides/boundary-condition-rulings.md new file mode 100644 index 000000000..b5b7f9605 --- /dev/null +++ b/docs/developer/guides/boundary-condition-rulings.md @@ -0,0 +1,89 @@ +--- +name: boundary-condition-rulings +description: Which boundary treatment to use in Underworld3 and why — rotated strong free-slip over Nitsche or penalty, Nitsche only for a condition that must evolve, natural tractions, essential data as controls in an adjoint, what happens at corners, and the mesh-side traps that leave a condition silently empty. Read with the class-level view of a solver, which lists the mechanisms it accepts. +families: [Stokes, Stokes_Constrained, VE_Stokes, NavierStokes, Poisson, Diffusion, AdvDiffusion, SteadyStateDarcy, TransientDarcy, Richards] +kind: guide +--- + +# Boundary conditions: the rulings + +A solver lists the mechanisms it accepts (`uw.systems.Stokes.view()`: +essential, natural, Nitsche, rotated free-slip, fault contact, and for a +saddle point the multiplier constraint). This page says which to choose, +and what each one commits you to. + +## Free-slip: rotated strong free-slip, not Nitsche or penalty + +`solver.add_rotated_freeslip_bc(conds, boundary, normal=None)`, value first: +`conds=0` is free-slip, a scalar or expression prescribes the wall-normal +datum strongly. + +- It enforces $v\cdot\hat n = 0$ to machine precision; Nitsche and penalty + leak at about $10^{-3}$. +- It is correct on curved, tilted and deformed boundaries: the normal is + taken per node, measure-weighted to match the facet integral the + assembler evaluates (#560). Leave `normal=None` unless the constraint + must follow the true surface rather than the mesh; an analytic normal is + exact for the geometry but carries a consistency error against the + faceted assembly. +- It works inside the nonlinear SNES and with geometric multigrid, and is + transparent to the tangent (`consistent_jacobian` behaves as it would + without it). +- Its reaction is the boundary normal traction, read with + `solver.boundary_normal_traction(boundary)`; no augmented-Lagrangian + splitting is involved. + +Governing document: [rotated free-slip](../subsystems/rotated-freeslip.md). + +## Nitsche: only for a condition that must evolve + +A hard rotated constraint cannot morph. A condition that changes character +in time, a Dirichlet-to-Neumann ramp or a traction that switches on, is a +Nitsche condition (`add_nitsche_bc`). Two rulings come with it: + +- The penalty scale must be recalibrated whenever the cell size is + redefined (#734). A value of 10 was a cliff, a hundred times worse than + 12.5, after the cell size changed under it; the scale is tied to the + cell size by #697. +- Nitsche is the one mechanism whose boundary Jacobian reads the unknown, + so a solver carrying it takes the matrix route for its adjoint and + refuses a parameter in a Dirichlet datum until the boundary tangent is + in the reaction term. + +## Natural conditions: tractions and fluxes + +`add_natural_bc(value, boundary)` is a facet load. `mesh.Gamma` in the +expression resolves to the facet normal: exact per quadrature point on an +external boundary, the declared analytic normal on an internal one, where +PETSc's normal is orientation-ambiguous (#327). A parameter in a natural +condition is differentiated by the adjoint as a facet part with nothing +further from the user. + +## Essential data + +`add_essential_bc` and `add_dirichlet_bc` are the same call. Components +given as `None` (or `sympy.oo`) are left free. Three things to know: + +- PETSc constrains the closure of each labelled boundary, corner vertices + included, and inserts the data boundary by boundary in the order they + were registered, so where two boundaries meet the later one's datum is + the one applied. Register the lid last if the lid's value is the one the + corner should carry. +- A named parameter in a datum is a control: `solver.gradient(...)` + carries its reaction term (#762). A bare number is not a control, since + it has no name. +- A datum given as a plain number is scaled by the model's units at + registration; a symbolic datum is scaled when it is compiled. Either way + the transcript's key shows the value as applied. + +## The mesh side: conditions that come out empty + +- Gmsh physical groups must be numbered in the order of the boundary + enumeration, or the conditions attach to the wrong stratum and are + silently empty. Check `mesh.view()`, which lists the boundaries with + their sizes at level 1. +- A patch on a fault must be at least three cells across to be split; + smaller, the split refuses. +- A boundary condition on a field that is not the solver's unknown (a + pressure datum on a Stokes solver) is accepted by the solver but not by + the adjoint, which refuses rather than drop it. diff --git a/docs/developer/guides/cetz-figures.md b/docs/developer/guides/cetz-figures.md new file mode 100644 index 000000000..18a0ec934 --- /dev/null +++ b/docs/developer/guides/cetz-figures.md @@ -0,0 +1,145 @@ +--- +name: cetz-figures +description: Build schematic / labelled-geometry figures for underworld3 papers using Typst + cetz. Use when the figure is primarily about topology, annotation, and math-typeset labels (meshes, solver diagrams, flow charts). Prefer Python → SVG → `#image()` for data-heavy figures (fields, colormaps, arrow plots) instead. +families: [] +kind: recipe +--- + +# cetz-figures + +Scaffold for Typst/cetz figures in `publications/**/figures/` alongside the +existing `arrays-sync-flow.typ`. This skill exists because upstream Claude +sessions draft cetz blind — this one actually compiles. + +## When to use cetz + +- Mesh schematics with a few labelled triangles / vertices / control points. +- Solver / data-flow diagrams (see `arrays-sync-flow.typ` in this repo). +- Anything where labels should render in the paper's math/text fonts. +- Anything that benefits from recompiling with the paper. + +## When to use something else + +- **Data-heavy plots** (scalar fields, colormaps, quiver plots, anything with + dense per-pixel or per-cell data) — generate SVG from Python/matplotlib, + include via `#image("foo.svg")`. cetz will fight you here. +- **Geometry computation** (Delaunay, intersections, interpolation) — do it in + Python offline, emit JSON with the shape `{"vertices": [...], + "triangles": [...], ...}`, let Typst just draw. + +## When NOT to use TikZ + +Evaluated and rejected: +- Slower compile than Typst. +- Drags in a LaTeX toolchain that isn't otherwise required by the project. +- No Typst math-font advantage over cetz for labels in our paper context. + +Keep TikZ in back pocket only if a co-author insists on TikZ source. + +## Before you draw (thesis-first discipline) + +When a user hands you a figure request — especially one replacing an +existing ASCII sketch, whiteboard photo, or reference figure — **do not +start transcribing**. The source artefact is a hypothesis about what +to communicate, not a specification. Before opening cetz: + +1. **State the thesis in one sentence.** What is the figure arguing? + If you can't write it plainly, you don't yet understand the figure. + Ask the user to articulate it. + +2. **Honestly audit the source.** If the original ASCII / sketch is + being replaced, it's often replaced *because it doesn't land well*. + Say what's broken about it before proposing the replacement — the + user often agrees and the new figure can do more than the old. + +3. **Enumerate design decisions as explicit questions, not assumptions.** + For a curved-boundary normals figure that's typically: + - Geometry (circle / ellipse / arc span / zoom level) + - Sampling density (number of facets, quadrature points per facet) + - Overlay vs. side-by-side + - Whether to show error quantitatively (arcs, annotations) or leave + it as a visible angle + - Where the figure lives in the repo (which doc / which branch) + Present a proposed interpretation with the decisions flagged; let + the user resolve them before you compile. + +4. **Only then open cetz.** Iterate visually — the discipline above is + about not committing to a design prematurely, not about planning + exhaustively. Once you start, compile often. + +The user's phrase "I'm not quite sure what this is intended to +illustrate" is the canonical trigger for this discipline. If you hear +it (or catch yourself about to transcribe without checking), stop and +do the four steps above. + +## Key gotchas (hit during iteration — not hypothetical) + +1. **Don't give a helper a parameter named after a `cetz.draw` export** + — `anchor`, `fill`, `stroke`. The cause is the `import cetz.draw: *` the + helper needs (gotcha 6): it runs *inside* the function body and shadows the + parameter, so the parameter name resolves to cetz's function rather than to + the value you passed. `anchor` panics with `"Unknown anchor 'anchor' for + element 'none'"`; `fill` gives `"expected color, gradient, tiling, or none, + found function"` pointing into `canvas.typ`, nowhere near your code. Rename + to `align-to`, `bg`, `edge`. See `cetz-cheatsheet.md`. + +2. **Clipping is a Typst concern, not a cetz one.** cetz has no `\clip`. + Wrap the canvas in `#box(clip: true, width: ..., height: ..., ...)` and + draw slightly oversized inside the canvas — the box clips the overflow. + +3. **Painter's algorithm — order matters.** Draw the background first, the + highlight second. No z-index exists. Verified in `mesh-demo.typ`. + +4. **Semi-transparency via `rgb(r, g, b, a)`** (alpha 0–255) or + `color.transparentize(col, 50%)`. Both work for fill and stroke. + +5. **Math in labels just works.** `content(pos, $v_1$)` renders in the + document math font. No escape hatch needed. This is a real cetz win over + SVG. + +6. **`import cetz.draw: *` inside the canvas closure.** Without it, `line`, + `circle`, `content` aren't in scope. Helper functions that draw need + their own `import cetz.draw: *` line inside. + +## Project layout pattern + +Each blog post or paper section gets its own subdirectory under +`figures/`, so a post's figures travel together: + +``` +publications/blog-posts/figures/ +└── / + ├── .typ # cetz drawing + ├── .png # committed output + ├── -data.json # (optional) precomputed geometry + └── generate-.py # (optional) Python that writes the JSON +``` + +Concrete example: `publications/blog-posts/figures/finding-particles/` +holds `mesh-demo.*` and `domain-demo.*` for the post +`finding-particles.md`. + +The JSON intermediate is the forward bridge to underworld3 — see +`underworld-bridge.md`. + +## Reference files + +- `cetz-cheatsheet.md` — what worked from memory vs. needed lookup. +- `underworld-bridge.md` — JSON schema for future `uw.meshing` export. +- `examples/` — self-contained copies (each with `.typ`, `.png`, `.json`, + and generator `.py`) of the figures this skill is scaffolded from. + These are snapshots; the live versions may have drifted if a post or + doc was iterated on further. + - `mesh-demo.*` — element-level point-in-cell test. + Live: `publications/blog-posts/figures/finding-particles/`. + - `domain-demo.*` — parallel domain centroid ambiguity. + Live: `publications/blog-posts/figures/finding-particles/`. + - `facet-vs-true-normals.*` — facet normal vs. smooth-surface normal + on a curved boundary. Live: + `docs/advanced/figures/curved-bc/`. + +## Canonical reference in the repo + +`publications/blog-posts/figures/arrays-sync-flow.typ` — prior cetz figure +in the repo, established version (0.3.4) and house style (hex colours, +helper-function pattern). Follow its conventions. diff --git a/docs/developer/guides/free-surface-convection.md b/docs/developer/guides/free-surface-convection.md new file mode 100644 index 000000000..4647fffa6 --- /dev/null +++ b/docs/developer/guides/free-surface-convection.md @@ -0,0 +1,226 @@ +--- +name: free-surface-convection +description: The Underworld3 free-surface convection method we are hardening — the THREE-NUMBER pointwise topography integrator (held-lid stress equilibrium h_∞ + L-stable exponential relaxation), NOT FSSA. Reach for THIS before touching any free-surface / dynamic-topography convection run, choosing a surface-update scheme, or "stabilising" a free surface. It records the method, why FSSA is explicitly rejected, and the failure modes. +families: [Stokes, AdvDiffusion] +kind: recipe +--- + +# free-surface-convection + +The free-surface scheme used in `~/+Simulations/FreeSurface/convection/fs4_compare.py` +and the design doc `docs/developer/design/FREESLIP_DYNAMIC_TOPOGRAPHY_FREESURFACE.md`. +**This is the method we are HARDENING — do not replace it; do not add FSSA.** + +> The 3-number integrator is necessary but NOT sufficient — see **Hardening strategies +> (2026-06)** below for the material-surface advection, tangential topography term, +> free-slip-inner nullspace, and graded/higher-order-mesh fixes that make it actually +> work. Reference impl + diagnostic tools live in `~/+Simulations/FreeSurface/convection/`. + +## ⚠️ NOT FSSA + +The `docs/examples/free_surface/advanced/Annulus*FS.py` examples use **FSSA** +(`add_natural_bc(δt·(Γ·v)Γ/2, "Upper")`). **That is NOT our method.** FSSA buys +stability by adding an implicit surface traction that **UNDER-deforms the surface** +— it trades accuracy for stability. Our scheme is designed to be stable **and** +accurate. If you find yourself adding `FSSA`, an `add_natural_bc` traction on the +free surface, or a `Gamma.dot(v)` stabiliser — STOP, you have the wrong method. +(Those example files are a template for a *different* approach, not this one.) + +## The three-number pointwise integrator (THE method) + +Two Stokes solves per step on the SAME mesh, then a pointwise surface update: + +1. **Free solve** — stress-free top (NO velocity BC on `Upper`; pressure datum is + pinned by the stress-free condition → no pressure nullspace). The surface + normal velocity `u_n` of this solve IS the kinematic rate `ḣ`. +2. **Held-lid solve** — a second Stokes solve with a RIGID free-slip held lid + (`u_n = 0`, via `add_nitsche_bc(0.0, "Upper", local_h=True)` — see [[project_nitsche_local_h_pr275]]) + and a DRIVING-ONLY body force. Its surface normal stress `σ_nn` gives the + equilibrium topography `h_∞ = -(σ_nn - mean)/ρg`. (The free solve forces + `σ_nn = 0`, so the equilibrium MUST come from the held-lid stress.) +3. **Pointwise exp step**, per surface node, from THREE numbers (`h`, `ḣ=u_n`, + `h_∞`): + ``` + γ = ḣ / (h_∞ − h) # local relaxation rate, clamp γ ≥ 0 + h ← h_∞ + (h − h_∞)·exp(−γ·dt) + ``` + L-stable: the step is bounded between `h` and `h_∞`, so it **cannot overshoot** + regardless of a noisy local `γ` (no "drunken sailor"). 1 extra solve/step; + beats RK4 at large dt. NO per-node freeze-clamp (that was the old `relax` bug). +4. The nodal surface increment is **carried inward by a Laplacian diffuser** + (smooth, minimal mesh deformation — NOT full mmpde adaptation), then + `mesh.deform()`. Uniform meshes are fine; adaptivity is NOT required. + +Reference impl: `fs4_compare.py` → `_surface_step`, `_h_inf` (held-lid σ_nn via a +`Projection`), `_surf_un`, `_carry_diffuser`. The free-slip RIGID-top run (no +surface motion) is the reference; the free-surface run is the same driving solve +PLUS this surface update. See [[project_fs4_adaptive_2x2]], +[[project_stress_equilibrium_freesurface]], [[project_freeslip_topo_freesurface]]. + +## Performance + +- Free-slip (rigid top) stagnant-lid runs are FAST. The free-surface cost is the + extra held-lid solve **plus** that moving the surface forces a COLD-START Stokes + each step (can't warm-start across a deformed mesh). +- Use **uniform** meshes for this problem — the diffuser gives minimal deformation, + no mmpde needed. (FMG works on a uniform `refinement=N` hierarchy; scalar solvers + must avoid FMG — PETSc err62, issue underworldcode/underworld3#276.) + +## Hardening strategies (2026-06) — the integrator alone is not enough + +The 3-number integrator moves the surface correctly, but several *other* things must +be right or it runs away / tangles. All implemented in `fs4_compare.py` (flags noted). + +### 1. Material-surface advection — THE key fix (`--advect-velocity`) +The runaway (`u_n` 42→125→285→445, cold lid leaking in, plumes punching through) was +NOT an `h_∞`/BC bug (`h_∞` is verified correct, even in the stagnant FK lid — held-lid +free-slip is the EASY case there). The bug: the surface moves by the L-stable relaxed +rate `ũ_n = Δh/Δt ≤ u_n`, but T was advected with the **stress-free solve velocity** +(surface-normal = full `u_n`). Net material then crosses the surface. A free surface is +a MATERIAL boundary: advect T with a velocity whose surface-normal = `ũ_n`. Modes: +- `consistent` (the right way): a THIRD Stokes solve, same buoyancy, `v·n̂ = ũ_n` + PRESCRIBED at the surface (penalty), tangential stress-free. `ũ_n = (shape_new−shape0)/dt` + = the full ∂h/∂t at fixed θ (correct ALE target). +- `blend`: `α·v_free + (1−α)·v_held`, `α = φ1(γΔt) = (1−e^{−γΔt})/(γΔt)` (the exp-decay + time-average). By Stokes LINEARITY this *is* the prescribed-`ũ_n` solve for UNIFORM α + (and free for FK, which is linear in v). BUT the single mean-α collapse is NOT close + enough once γ varies per surface node — the planform diverges (mode-3 vs mode-1/2), + throughflow ~23 vs ~0.08. Per-node α breaks div-free (∇α·(v_free−v_held)). So the + per-node `consistent` 3rd solve is REQUIRED for structured planforms. +- `free`: advect with stress-free v (the inconsistent baseline — the runaway). + +### 2. Tangential topography advection (`--no-tangent-advect` to disable; default ON) +The pointwise relaxation omits the `v_t·∂_s h` term — a surface rotation/convergence +should carry the topography pattern along the surface; without it you get edge artefacts +where ∂_s h is large (plume-bulge edges). Fix = operator split per step: (1) departure- +point semi-Lagrangian transport of the surface shape in θ by `ω = v_t/r`, then (2) the +L-stable normal relaxation. Lowers throughflow + improves mesh quality. + +### 3. Free-slip inner boundary — rotation nullspace (`--inner freeslip`) +The rigid rotation `[-y,x]` is a velocity nullspace ONLY while the boundary is CIRCULAR. +Once the free surface DEFORMS, do NOT attach it to the held/consistent solves (→ held +22 s/`DIVERGED_LINEAR_SOLVE`, throughflow blow-up). Keep `petsc_use_pressure_nullspace`; +strip the gauge with the exact post-solve projection `_project_out_rotation` on `v`, +`v_cons` (drives advection) AND `v_h` (one consistent non-rotating frame). The undeformed +free-slip *reference* (`--surface freeslip`) is fine WITH the nullspace attached. + +### 4. Graded / higher-order meshes (drive node movement consistently) +- **Surface-ring detection**: tie the tolerance to the FINEST cell + (`0.5·mesh.get_min_radius()`), NOT the nominal `cellsize`. On a gmsh-graded mesh + (`cellSizeOuter`) the old tolerance scoops the first interior ring → a 2%-thick + "surface band" → tangling (looks like the surface "destroying itself"; it isn't — + the diffuser was fed a corrupt surface). BETTER (TODO): build the ring from the DMPlex + `Upper` label (`dm.getLabel("Upper").getStratumIS`), removing the tolerance entirely. +- **Node movement**: the solve velocity is P2; the mesh geometry is P1. Drive `u_n` and + the tangential transport from a P1 length-smoothed `Vector_Projection` of V (`v_p1`), + NOT a point-evaluation of the P2 field. +- **Stress smoothing**: `topo_proj.smoothing_length` = a fixed PHYSICAL length + (`--smooth-length`), not cell-count, so `h_∞` is mesh/order-independent. + +### 5. Cost — there is no acceleration win (don't chase it) +3 Stokes solves/step (free→u_n, held→h_∞, consistent→advect). Warm-start does NOT help +(outer KSP already 1 iter; FMG supplies its own nested guess — measured SLOWER). Blend- +skip rarely fires (α-spread always large). Operator/PC reuse: already reused across +`solve()`s (first 5.5 s setup, steady ~945 ms = irreducible FMG solve; RHS-only resolve +same cost). The per-step cost is the geometric FMG hierarchy REBUILD on the deforming +mesh — intrinsic to moving meshes, "live with it." The UNIFIED-PENALTY single solver +(`penalty·(v·n̂ − V₁·n̂)·n̂`; penalty=0→free, V₁=0→held, V₁=ũ_n→consistent — held & +consistent share the matrix) is the cleanest formulation (no recompile on the constant) +but doesn't cut the irreducible solve. + +### Elastic-plate flexure `h_∞` — IMPLEMENTED (`--flexure-D`) +Generalizes the LOCAL Airy `h_∞ = −σ_nn/ρg` (the D=0 limit) to a flexed plate +`(D ∂_s^4 + ρg) h_∞ = −σ_nn`, solved SPECTRALLY on the ring (serial Fourier — the +feasible substitute for UW3's blocked 1D-manifold FE solve): per mode +`h_∞(m) = −σ_nn(m)/(ρg + D(m/r_o)⁴)`. `D` sets the flexural wavelength `(D/ρg)^{1/4}` +and damps short-wavelength loads — the physically-grounded, mesh-independent length- +smoothing. In `_h_inf` (h_∞-ONLY — the stable form). **PROTOTYPE — amplitude response correct +(stiffer plate → less deflection) but it does NOT low-pass the SURFACE**: filtering `h_∞` only +sets a smooth set-point; the surface still picks up short-wavelength content from the SL +tangential transport + partial relaxation. "Filter every surface number (h, ḣ, h_∞)" was +TRIED and REJECTED — filtering the GEOMETRY `h` injects a spurious smooth-the-mesh motion into +`ũ_n` → flow runs away (Vrms 50→345); filtering `ḣ` alone is stable but elevates Vrms with no +benefit. So making flexure a TRUE surface low-pass without destabilising is OPEN/hard +(`_flex_filter` helper is in place). Examples: `stagnant_lid_mode1_study/{figures/flexure_*.png, +runs/flexure_D*}`. + +### Open / next (not yet done) +- **Label-based surface ring** (replace the radial heuristic with the `Upper` stratum). +- **Flexure D calibration** to a realistic lithospheric flexural wavelength. + +### Diagnostic tools (`~/+Simulations/FreeSurface/convection/stagnant_lid_mode1_study/scripts/`) +- `heldlid_hinf_check.py` — verify `h_∞` via 4 independent free-slip enforcements × Δη sweep +- `stitch_compare.py` — side-by-side montage of per-run dirs (`--dirs a,b,c`) +- `unified_penalty_solver.py` — the one-solver penalty formulation probe +- `resolve_timing_probe.py` — repeated-solve / lag-Jacobian / reuse-PC timing + +## ★★ Body force must be FULL Boussinesq on a deforming mesh (2026-07-26) + +**On any run where the mesh surface actually moves, the body force must retain the +ρ₀ background: `bodyforce = (thermal_buoyancy − rho_0_g)·r̂` (i.e. ρ = ρ₀(1−αΔT)).** +The reduced (driving-only) form has **NO restoring force for surface deformation +anywhere in the momentum system** — `buoyancy_scale` enters only the kinematic target +h∞ = −σ_nn/ρg, which exerts zero force on the flow. Consequences and evidence (all +measured, `~/+Simulations/FreeSurface/annulus_fs_convection/teaching/`): + +- **Relaxation A/B (`relaxation_test.py`)** — imposed 5% mode-4 topography, no thermal + driving: reduced form = surface FROZEN (velocities are round-off); full density = + **262→5 km in 10 steps at the Cathles rate** (measured 2.1e4 vs ρ₀g/(2ηk) = 2.5e4, + within 20% at res 0.06). +- **Convection A/B (`restoring_force_demo.png`, gap-Ra 3e4, ρ₀g 1e6)** — reduced h + grows monotonically without limit; full density *rings* about its supported amplitude + early (damped), then tracks the developing flow ~35% lower with HIGHER Nu (1.93 vs + 1.65). Without the term, soft surfaces run away entirely (measured to 50% of radius). +- **The FreeSurface manager needs NO change** — held/consistent solves inherit + `stokes.bodyforce`; h∞ then self-consistently includes the self-load (a moving target + the integrator follows cleanly — verified in the relaxation test). + +Rules that come with it: +1. **ρ₀g and Ra are COUPLED: αΔT = (thermal buoyancy coeff)/ρ₀g must be ≲ 0.3.** + Softness is not a free knob — the old "soft" runs at ρg=2e5 were αΔT = 1.2–4 (no + such fluid) and their runaway was partly parameter nonsense. Soft surfaces require + weak driving. +2. **Tighten the solver tolerance (~1e-8)**: the hydrostatic RHS dominates the dynamic + signal by 1e3–1e4, and a relative tolerance judges the total. +3. Prefer `snes_type=ksponly` for the isoviscous solves; the huge RHS makes `newtonls` + thrash worse. +4. Driver switch: `fs_convection.py -uw_full 1`. +5. **`FreeSurface(background_buoyancy="analytic")` is REQUIRED with full density** + (98961147): the recovered reaction contains the self-load +h_current and the + reduced-form negation otherwise flips it -> h_inf = -h + drive/rho_g (parks at HALF + equilibrium; -1 eigenvalue = period-2 ringing; steady flow THROUGH the stationary + surface). "analytic" subtracts the geometric height - no extra solve. The CBF + recovery itself is exact (probe lesson: select boundary DOFs by LABEL, never a + radius mask, on a deformed mesh). `background_buoyancy=` = exact two-reaction + reference mode. + +Related transport fact (same campaign): the serial T-blow-up on deforming FS meshes was +the **old-frame SL reach-back** amplifying per loop cycle (issue #423; smaller dt makes +it WORSE) — retired in `b507aca1`; the manager now uses the standard ALE path + clamp + +deform-aware foot restore. The parallel datum defect is #421. + +## Failure modes — symptom → cause + +| Symptom | Cause | +|---|---| +| Surface deforms but `u_n` RUNS AWAY (e.g. u_n 42→125→285→445), cold lid leaks in, plumes punch through | **RESOLVED**: material-surface advection inconsistency — T advected with stress-free `u_n` while surface moves by relaxed `ũ_n`. Fix = `--advect-velocity consistent` (Hardening §1). NOT an h_∞ bug, NOT fixed by FSSA. | +| `held` solve 22 s / `DIVERGED_LINEAR_SOLVE`, throughflow blows up, with `--inner freeslip` | rigid-rotation `[-y,x]` attached as a nullspace on the DEFORMED (non-circular) surface — invalid. Don't attach it on the moving surface; use the post-solve projection (Hardening §3). | +| Graded-mesh surface "destroys itself" (q→0.2, h_max 2% at step 1) | surface-detection tolerance scooped the first interior ring → 2%-thick band, NOT real deformation. Tie tolerance to finest cell (Hardening §4). The diffuser is innocent. | +| Topography GROWS without saturating; flow-through persists even at huge deformation | **missing ρ₀ background** — reduced body force has no surface restoring force (see the Full-Boussinesq section above). Fix = ρ = ρ₀(1−αΔT); check αΔT ≲ 0.3 | +| T leaves [0,1] on a deforming mesh, mesh-locked hot/cold spikes in the squeezed band, worse at SMALLER dt | old-frame SL reach-back amplifying per loop cycle (issue #423) — use the standard ALE path (`old_frame_traceback=False`, the manager default since b507aca1) | +| Stress-free top but surface not updated each step | nothing stops throughflow (the stress-free top is an open boundary unless the integrator moves the surface to track `u_n`) | +| Nu decays when it should be steady (kinematic free surface) | LAG: the SL foot reaches beyond an under-moved surface → cold pump. Fix = the h_∞ relaxation, not more smoothing | +| Surface "mountain" / one-step spike on adaptive mesh | held-lid Nitsche penalty over-stiffened by GLOBAL h; use `local_h=True` (default, PR #275) = `mesh.cell_size()` | + +## Dead ends (already tried — do NOT repeat) + +- **FSSA** signed-traction free-surface: diverges / under-deforms — rejected. +- **High-k post-smoothing** of the surface: the instability is low-m, smoothing + the wrong band. +- Per-node freeze-clamp in the relaxation (the old `relax` fatal bug). + +## Diagnose by + +`h_max` (deflection, as % of r_o), `u_n` / `vhmax` (surface throughflow — should NOT +grow unbounded), `hinf_max` (the equilibrium target), `vrms`, `Nu`. Compare the +free-surface run against the free-slip RIGID-top reference at matched physical time. diff --git a/docs/developer/guides/nonlinear-solver.md b/docs/developer/guides/nonlinear-solver.md new file mode 100644 index 000000000..1674d0bf7 --- /dev/null +++ b/docs/developer/guides/nonlinear-solver.md @@ -0,0 +1,323 @@ +--- +name: nonlinear-solver +description: How to make a hard nonlinear Stokes solve (Drucker-Prager / yield-stress viscoplastic) CONVERGE reliably in Underworld3 the way the working recipe actually does it — automatic warm-start (one Picard step on a cold start) plus a MULTI-SOLVE δ-continuation (constant δ per solve, warm-start the next, sharper δ), the consistent-Newton tangent, and a non-symmetry-safe multigrid smoother. Reach for THIS when a viscoplastic solve stalls / diverges and you are about to hand-tune PETSc options, ramp δ, or "just add a monitor". It carries the CONFIG TRAP LIST — the setup mistakes that each produce a different failure a few steps in — and the one thing you must NOT do (ramp δ inside a single SNES solve). For the yield-law maths and which tangent per model, see `plasticity-solvers`. +families: [Stokes, ViscoPlasticFlowModel] +kind: recipe +--- + +# nonlinear-solver + +The recipe that gets a **hard viscoplastic (Drucker–Prager) Stokes** problem to +converge, and — more importantly — the list of setup mistakes that stop it. The +central lesson from the Spiegelman hard-case study (`η_bg=1e26`, `V=10`): every +failure was a **solver-configuration** error, not a bad Jacobian. If the correct +setup is a minefield for an expert, that is an API regression — so the goal is to +make the correct path the default path. + +Design of record: `docs/developer/design/nonlinear-solver-homotopy-warmstart.md`. +Yield-law maths, tangent-per-model, quadratic-convergence check: `plasticity-solvers`. + +--- + +## The recipe (what actually converges) + +1. **Warm start.** Start the continuation at **large δ**, where the yield surface is + smooth and the problem is easy, and take **one Picard step** into the Newton + basin. One Picard step is defect-correction iteration 1 — contractive, cheap. From + a *warm* iterate, take **no** Picard step (it wastes the good quadratic start). + A cold `v=0` start is safe on its own terms: `ε̇=0` makes `η_pl` infinite, which + the soft-min carries to the viscous branch (see the trap list for the one form + that must be written carefully). + +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. + +--- + +## DO NOT ramp δ inside one SNES solve + +Ramping δ **inside a single SNES solve** (a `SNESSetUpdate` callback that sharpens +the yield surface between Newton iterations) is **proven dead** — it diverges +`DIVERGED_LINEAR_SOLVE` after ~2 iterations and grinds for ~2 hours, **even on the +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*.** (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.) + +--- + +## CONFIG TRAP LIST + +Each of these produces a *different* failure a few steps in — that is why the hard +case felt like whack-a-mole. Check them first. + +| Trap | Symptom | Fix | +|---|---|---| +| **Perfect plasticity's consistent tangent is SINGULAR along the flow**: on the hard-`Min` plastic branch η = τ_y/2ε̇_II, so 2η + 2η′ε̇_II = 0 — the velocity block is symmetric but semi-definite in every yielded cell. (An earlier version of this row blamed *asymmetry*; that is wrong for any η(ε̇_II) law — the rank-one term η′ ε̇⊗ε̇/ε̇_II is symmetric. Pressure-dependent yield adds a non-symmetric v–p coupling, not a non-symmetric velocity block. Corrected 2026-08-26, maintainer review.) | benign while yielded cells are few (the viscous neighbours regularise); with a large yielded fraction the velocity sub-solve caps out and Newton stalls at ~1e-3, no failure reason | give the plastic branch a positive tangent: a small δ soft-min (`yield_mode="softmin"`, powermean, `yield_anchor="yield"`), a rounded viscosity floor, or rate-strengthening ξ; Picard converges regardless (full 2η stiffness) but is linear-rate. The FMG bundle's `gmres`+`sor` smoother is Newton-safe either way | +| `preconditioner="fmg"` (vs explicit `pc_type=mg` + manual mg opts) | outer KSP "converges" in **1 iteration** → no real Newton correction → stall → `DIVERGED_LINE_SEARCH` | use explicit `pc_type=mg` with the smoother opts above; bound the outer KSP (`ksp_max_it`~80) so a hostile step fails fast | +| Cold plastic start `v=0`, or any rigid/unyielded point | `DIVERGED_FNORM_NAN` at iteration 0 | **Not** a div/0: `ε̇=0` gives `η_pl=+inf`, which `Min` and the sqrt soft-min carry correctly to the viscous branch. Only a soft-min form that computes `η_ve·η_pl/(η_ve+η_pl)` breaks (`inf/inf`). Fixed in the power-mean; if you hand-roll a blend, write the harmonic mean as `η_ve/(1+η_ve/η_pl)`. **Do not reach for a strain-rate floor** — it hides this rather than fixing it | +| LU velocity block with all-Dirichlet-ish BC | pressure nullspace singular | attach the Stokes nullspace / avoid a bare LU there | +| Hand-rolled `snes_monitor` to "see what's happening" | you read residuals but miss the tell | use `solve_with_diagnostics` / `get_snes_diagnostics` instead (below) | + +**The diagnostic tell:** `solver.get_snes_diagnostics()["linear_iterations"] ≈ 1` +per Newton step means the linear solve is doing **no real work** (the FMG-1-iteration +trap). A healthy consistent-Newton solve does real Krylov work each step and +converges quadratically. Use `solve_with_diagnostics()`, not a hand-rolled monitor. + +--- + +## Automatic warm-start (Layer 1 — landed) + +`solver.has_solution` is a **public, read-only** status flag: `True` only after a +solve whose SNES converged; reset on a structural rebuild (remesh / adapt / +mesh-mover — the `is_setup=False` hook); kept through coefficient changes (viscosity, +δ, BC values, time step). A **diverged** solve leaves it `False`, so the next solve +auto-cold-starts rather than warming off a corrupted iterate. + +On a **cold** (`zero_init_guess=True`) Stokes solve under the **consistent-Newton +tangent**, a single Picard step is now taken automatically (reusing the existing +`picard=1` machinery). The default (frozen) tangent path is bit-identical. + +```python +stokes.consistent_jacobian = True +stokes.solve() # cold → one automatic Picard step, then Newton +if stokes.has_solution: + ... +``` + +--- + +## Implementation status (this line of work) + +- **Layer 1a — DONE:** `has_solution` + cold consistent-Newton Picard warm-up + (`petsc_generic_snes_solvers.pyx`; test `test_0201`). +- **Layer 1b — DONE:** `zero_init_guess` is tri-state — `None` (default) auto-detects + from `has_solution`, `True` forces fresh, `False` insists on warm. Note warm and cold + agree only to the *convergence tolerance*, not bitwise. +- **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 — 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. 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. + +--- + +## Multigrid depth — how to measure a smoother honestly + +**A two-level hierarchy is a coarse-grid correction, not a V-cycle.** Smoother +comparisons made on one are misleading: the gmres-over-richardson margin measured on +the Spiegelman notch is only 5 % at 3 levels but **25 % at 4** (ρ per V-cycle 0.746 → +0.560), because a deeper cycle applies the smoother on more coarse operators. Judge a +smoother at depth or not at all. + +To get depth without a monster problem, refine a **deliberately ultra-coarse NESTED +base**: `make_notch_mesh.py 1` (492 cells) + uniform `refinement=N` gives 3 levels / +7,872 cells at `N=2` and 4 levels / 31,488 at `N=3` — deeper *and* smaller than the old +2-level 38,580-cell setup. In MG you want the coarsest grid as coarse as it can be +before the problem breaks down. + +- **Never use a non-nested hierarchy** here — it does not give strong MG convergence + (maintainer ruling). Uniform refinement nests by construction. +- Accepted tradeoff: uniform refinement does **not** snap new boundary nodes back to + the analytic notch arcs (no CAD/EGADS model attached), so the corner geometry is + frozen at the coarse mesh's chords on every level. + +Measure with `fmg_contraction_probe.py` (ρ_MG per V-cycle; `<0.5` healthy, `0.8–0.95` +struggling, `≥0.98` hangs) or `smoother_depth_sweep.py` (pays the mesh build + viscous +seed once, sweeps smoothers in-process) in the Spiegelman study. + +**`solve_report` cannot see the smoother.** It records the *Newton* contraction; the +outer KSP is Eisenstat–Walker-collapsed to ~1 iteration/step, so the smoother's work +hides inside the velocity sub-block. Probe the `fieldsplit_velocity_` sub-KSP directly. + +**A smoother will not rescue small ξ.** At the hard corner the failure is operator +conditioning — the coarsest grid cannot represent the viscosity contrast — and at 4 +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). + +## The Schur complement: pair the penalty with FMG, never with GAMG + +**Symptom this is for**: the velocity block's iteration count is rock solid but +the pressure sub-solve wanders into the hundreds and eventually stops +converging. + +**First: it is probably not the pressure block.** `S = -B A^-1 B^T` is applied +*through* the velocity solve, so a velocity solve that exits at its iteration +cap makes the Schur operator inconsistent between applications — and no Krylov +method converges against an operator that moves under it. The pressure block +then caps too, and the outer flounders. Measured on SolCx (eta 1e6, P2-P0disc, +h=1/30), changing **only** `fieldsplit_velocity_ksp_max_it`: + +| velocity cap | sec | outer | pressure/app | velocity/app | +|---|---|---|---|---| +| 200 (default) | 976.0 | 44 | **200.0** | **200.0** | +| 5000 | **25.6** | **2** | **30.0** | 618.0 | + +**38x from a number that is not in the pressure block**, and the velocity error +is identical in both rows. Before tuning the Schur solve, check whether either +block sat at exactly its cap — `solve_report.sub` gives iterations and +applications per block, and a per-application count equal to the cap to the +digit is the tell. + +**Then: the penalty is the lever on the Schur count, and it needs FMG.** +`stokes.penalty = lambda` adds `lambda*mu*(div u)(div v)`, which makes the +eta-scaled mass matrix a better approximation to S. Matched on one mesh +(2592 cells), same discrete solve, only the velocity preconditioner differs: + +| lambda | velocity PC | sec | outer | Schur/app | velocity/app | velocity total | +|---|---|---|---|---|---|---| +| 0 | GAMG | 15.49 | 2 | 125.5 | 94.7 | 24802 | +| 0 | **FMG** | **3.88** | 1 | **59.0** | **8.8** | **546** | +| 10 | GAMG | 20.68 | 7 | 22.3 | **199.9 capped** | 33976 | +| 10 | **FMG** | **3.06** | 1 | **18.0** | **13.5** | **270** | + +- **With FMG, `penalty = 10` improves every axis at once**: 21% faster, Schur + count 3.3x smaller, total velocity work halved. FMG absorbs grad-div + augmentation (8.8 -> 13.5 iterations per application); GAMG does not + (94.7 -> capped). +- **With GAMG, do not use it at all.** The same `penalty = 10` makes the solve + *slower* (15.49 -> 20.68 s), because augmentation is exactly what drives GAMG + into its cap. Uncapping rescues it to 11.03 s but it still needs **833** + iterations per application, and FMG is 3.6x faster on the same mesh. + Feasible is not competitive. + +**The accuracy cost is consistent, so it is safe to pair by default.** The +penalty is grad-div, not a true augmented Lagrangian — `div(P2)` is not inside +`P0`, so the term does not vanish at the discrete solution and it does perturb +the answer. But the perturbation converges away: same rate, and the gap shrinks +under refinement. + +| cells | lambda=0 v err | rate | lambda=10 v err | rate | gap | +|---|---|---|---|---|---| +| 648 | 2.112e-1 | — | 2.327e-1 | — | 1.102 | +| 2592 | 1.266e-1 | 1.67 | 1.376e-1 | 1.69 | 1.087 | +| 10368 | 8.727e-2 | 1.45 | 9.305e-2 | 1.48 | **1.066** | + +For a pressure-dependent constitutive law use the mechanical pressure, +`p_mech = p - lambda*mu*(div u)`; the raw `p` is the multiplier. + +**Traps.** + +- **FMG needs a refined base or you silently get GAMG.** Measured: + `refinement=0` -> one hierarchy level -> default velocity PC is `gamg`; + `refinement=2` -> `mg`. So `penalty` set "with FMG" on an unrefined mesh is + actually the harmful GAMG pairing. Check + `snes.getKSP().getPC().getFieldSplitSubKSP()[0].getPC().getType()`, or read + `solver.pc_fallbacks`. +- **Scaling `saddle_preconditioner` by a constant does nothing** — it does not + change the Krylov subspace. `1/eta` and `101/eta` both give 28 iterations, + identical to every digit, so an "AL-matched" `1/(eta*(1+lambda))` cannot help. + The 1/eta *weighting* itself is worth 1.9x (28 vs 52 with a flat `1`). +- **Eisenstat-Walker is inert under `snes_type=ksponly`** — identical iterations + and error on or off. And `outer 1` is not an EW artefact: it is what a full + Schur factorisation gives when the Schur complement is solved well. +- Measurements: `~/+Simulations/pressure_schur_625/` (#625). + +## Gotchas + +- **`./uw build` → `amr-dev` env**; verify `uw.__file__` is the worktree site-packages. +- **Run VEP/consistent-Newton tests UNFORKED** — `pytest --forked` SIGABRTs (fork of + multithreaded PETSc). +- Benchmark **every** default change — "Solver Stability is Paramount". +- ξ (rate-strengthening) is a **non-homotopic** regularisation: put a user loop + *around* `solve()`, never inside the δ-march. + +## Reference + +- Design: `docs/developer/design/nonlinear-solver-homotopy-warmstart.md`, + `jacobian-consistent-tangent.md`, `solver-strategies-catalogue.md`. +- 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` (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/docs/developer/guides/plasticity-solvers.md b/docs/developer/guides/plasticity-solvers.md new file mode 100644 index 000000000..6d0e3c56c --- /dev/null +++ b/docs/developer/guides/plasticity-solvers.md @@ -0,0 +1,228 @@ +--- +name: plasticity-solvers +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`. +families: [Stokes, ViscoPlasticFlowModel, ViscoElasticPlasticFlowModel] +kind: recipe +--- + +# 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 what the yield campaigns measured actually works — +and records what was retired. + +**The default call is now just:** + +```python +stokes.constitutive_model = cm # ViscoPlastic / ViscoElasticPlastic / TI-VEP +cm.Parameters.yield_stress = tau_y # finite -> plasticity active +stokes.consistent_jacobian = True # Newton tangent (non-elastic DP; see table) +stokes.solve() # cold start takes ONE Picard step automatically +``` + +--- + +## 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. **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 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 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) +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) +# a nonzero difference == the Newton form is present +``` + +--- + +## δ 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. + +--- + +## 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)`. + +--- + +## Failure modes → fixes + +| Symptom | Cause | Fix | +|---------|-------|-----| +| `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) | + +--- + +## Gotchas + +- **`./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). +- `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. + +--- + +## Reference + +- 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 +tangent — the origin of the "~20 iterations is intrinsic" folklore. + +## SNESFAS — do not reach for it + +Nonlinear multigrid (SNESFAS) looks tempting for hard viscoplastic solves but is +**not a viable option** at present (maintainer ruling 2026-07-17): there are no +good preconditioners for the nonlinear hierarchy, and it abandons the robust +linear-solver path (consistent tangent / continuation + fieldsplit + MG) that +this skill is built around. It stays options-only for experiments; treat it as a +future investigation. See `docs/developer/design/solver-strategies-catalogue.md` +and `MULTIGRID_MINIMAL_CONTROL_2026-07.md` (ruling 6). diff --git a/docs/developer/guides/transport-schemes.md b/docs/developer/guides/transport-schemes.md new file mode 100644 index 000000000..d6fa772ba --- /dev/null +++ b/docs/developer/guides/transport-schemes.md @@ -0,0 +1,77 @@ +--- +name: transport-schemes +description: Which transport scheme and which time history to use in Underworld3, and why — nodal, integration-point or grid histories; semi-Lagrangian, Eulerian SUPG or Lagrangian swarm transport; the Courant number to run at; what each choice does to a settled state and to a peak. The evidence is the tests and notes named beside each ruling. +families: [AdvDiffusion, NavierStokes, Eulerian, EulerianSUPG, SemiLagrangian, Lagrangian, Lagrangian_Swarm, IntegrationPointSemiLagrangian] +kind: guide +status: draft, rulings to be confirmed +--- + +# Transport schemes: which one, when + +A time-dependent solve in Underworld3 is a residual plus a history: the +history (`DuDt`, `DFDt`) says where the quantity was at the previous levels +and how it got there. The scheme is the choice of where that history lives +and how it is carried. The classes describe themselves +(`uw.capabilities("histories")`); this page says which to choose. + +## The choices + +| scheme | history lives on | carried by | use for | +|---|---|---|---| +| `Eulerian` | mesh nodes | nothing moves; the transport term is in the residual | diffusion-dominated fields, or with SUPG below | +| `EulerianSUPG` | mesh nodes | implicit advection with streamline-upwind stabilisation, assembled in the residual | the momentum equation of `NavierStokes`; a field advected by a resolved velocity | +| `SemiLagrangian` | mesh nodes | traced back along the flow to a departure point and interpolated there | advection-diffusion at moderate Courant number; the transport the adjoint can differentiate | +| `IntegrationPointSemiLagrangian` | integration points | the same trace, from the quadrature points | stress and other flux histories that must not be smoothed through the nodes | +| `Lagrangian_Swarm` | particles | the particles move; the mesh reads a proxy | material identity, and any history that must follow the material exactly | +| `Symbolic` | nowhere | the user's own expression | a history the script supplies itself | + +## What the tests established + +**Advection of a step (Waters and King, 2026-09-16, #749).** With the +same time integrator, a nodal history converges as the timestep falls; an +integration-point history degrades as the timestep falls, ringing and then +diverging; a grid history diverges at every timestep. The peak of the +profile is set by the integrator, the settled state by the transport. Rule: +for a scalar field carried by the flow, put the history on the nodes. + +**Stress histories (2026-09-11, #735).** For a viscoelastic flux history +the grid and integration-point histories converge together with +resolution; the nodal history converges to a different answer, with about +half again too much stress in the first two cells off a no-slip wall, +because the nodal projection smooths the history through the wall. Rule: +a flux history lives on the integration points (or on the grid, with DEVSS +opted in); a scalar field's history lives on the nodes. The two rules are +not in conflict: they are different quantities. + +**Courant number for the integration-point trace (#703, #737).** The +integration-point semi-Lagrangian scheme runs at Courant number about +one, not below it. Damping it to run at smaller steps was tried and +disliked; the memory term amplifies the low-Courant mode. Refine the mesh +and the timestep together. + +**The momentum equation (#687).** `NavierStokes` carries its momentum +transport as `EulerianSUPG`: implicit advection in the residual, with a +partition-independent cell size in the stabilisation. + +**What the adjoint can differentiate.** A semi-Lagrangian trace is +differentiable in the velocity, and the interpolation at the departure +points is materialised, so a run built on it is adjointable end to end. A +particle step is adjointable exactly when the particle set is fixed across +it, which `swarm.advection` checks by counting. + +## The time integrator is a separate choice + +`order` sets the depth of the history and the scheme's order in time; +`theta` sets the weighting (`0.5` is Crank-Nicolson, `1` backward Euler). +`model.step(dt)` carries the clock, and the coefficients of the scheme are +exact rationals in the residual, so the transcript's key shows the +integrator a part used. A fixed timestep keeps an objective from depending +on the control through the schedule; an adaptive one records its decisions +through `model.rewind(reason=...)`. + +## Rulings still open + +- Whether a nodal history should ever be offered for a flux quantity, or + refused. +- A default Courant target for the nodal semi-Lagrangian scheme, and + whether `estimate_dt()` should report it. diff --git a/docs/developer/guides/uw-visualisation.md b/docs/developer/guides/uw-visualisation.md new file mode 100644 index 000000000..31843102f --- /dev/null +++ b/docs/developer/guides/uw-visualisation.md @@ -0,0 +1,127 @@ +--- +name: uw-visualisation +description: Render Underworld3 mesh fields (T, V, viscosity, the adapted mesh) correctly with PyVista. Use whenever you need to SEE a UW3 result — a field colormap, the moving/adapted mesh, streamlines, or compare runs. Reach for THIS before hand-rolling a renderer; getting the four cosmetic settings wrong makes renders look grey/patchy/blocky and wastes a round-trip with Louis. +families: [] +kind: recipe +--- + +# uw-visualisation + +Canonical PyVista recipe for Underworld3 fields. This exists because every fresh +Claude session re-derives the renderer and gets the colormap / background / +lighting / DOF-sampling wrong, producing "grey/patchy/weird" images Louis +rejects. The settings below match his reference renders exactly. + +**Use PyVista (`underworld3.visualisation`), NOT matplotlib.** Louis reaffirmed +this even after seeing the legacy matplotlib renderer +(`scripts/fault_convection_frames.py`) — that one is NOT preferred. + +## Hard rules (artifacts + output location) + +- **Outputs go under `~/+Simulations/...`, NEVER `/tmp`** (Louis can't view /tmp + or harness task paths). Mirror the run's `--sim-dir`; write `T_.png` into + the run directory, comparison figures into the sim-dir root. +- `pv.OFF_SCREEN = True` at import; finish with `pl.screenshot(path); pl.close()`. + +## The field+mesh pattern (copy this exactly) + +```python +import numpy as np, underworld3 as uw, underworld3.visualisation as vis, pyvista as pv +pv.OFF_SCREEN = True + +mesh = uw.discretisation.Mesh(f"{label}.mesh.00000.h5") # or the live mesh +T = uw.discretisation.MeshVariable("T_v2p1", mesh, 1, degree=3, continuous=True) +T.read_timestep(label, "T_v2p1", 0, outputPath=D) # or use the live var + +pv_T = vis.meshVariable_to_pv_mesh_object(T) # Delaunay through T's OWN DOFs +pv_T.point_data["T"] = np.asarray(T.data[:, 0]) # attach DOF values DIRECTLY (P3-faithful) +edges = vis.mesh_to_pv_mesh(mesh).extract_all_edges() + +pl = pv.Plotter(off_screen=True, window_size=(1000, 1000)) +pl.set_background("white") # rule 2 +pl.add_mesh(pv_T, scalars="T", cmap="RdBu_r", clim=(0, 1), # rules 1 + clim required + show_edges=False, lighting=False) # rule 3 +pl.add_mesh(edges, color="black", line_width=0.5, lighting=False) # mesh overlay +pl.view_xy(); pl.camera.zoom(1.3) +pl.screenshot(out); pl.close() +``` + +## The four things that make renders look bad (all COSMETIC) + +1. `cmap="coolwarm"` → muddy grey-lavender midtone — this IS the "blue/grey/red" + Louis rejects. **Use `cmap="RdBu_r"`** (clean blue→white→red). +2. PyVista's default grey background bleeds through RdBu_r's white (T≈0.5) → dirty + grey. **Always `pl.set_background("white")`.** +3. Default lighting darkens the colormap. **Always `lighting=False`** on every + `add_mesh`. +4. Re-evaluating via `scalar_fn_to_pv_points` / vertex-only sampling drops the + high-order DOFs → blocky. **Attach `T.data[:,0]` directly** to the DOF-cloud + mesh from `meshVariable_to_pv_mesh_object` (it is correct for annulus/box/disc + — do NOT avoid it). `clim` MUST be passed (default `clim=""` trips `np.any`). +5. Resampling ANY field (even P1) onto a regular pixel grid via + `uw.function.evaluate` **dapples at element boundaries** — grid points that + straddle a facet get located into a neighbouring cell with slightly-off + reference coords (Louis: "artefacts across the elements", S-fault rig). + Render derived fields NODALLY on the mesh's own triangulation instead: + evaluate at `mesh_to_pv_mesh(mesh).points` (exact at vertices for P1, + whichever cell the locator picks), attach as point_data, let VTK + interpolate WITHIN elements. On a SPLIT mesh never Delaunay the DOF cloud + (it re-triangulates across the slit) — use the mesh's own cells. + +## Seeing the MESH (adaptation / moving mesh) + +The full-annulus T colormap **washes out mesh detail** — at whole-domain zoom the +grading is invisible. To judge adaptation you MUST crop: + +- Zoom the feature region with a parallel camera: + `pl.camera.parallel_projection = True; pl.camera.parallel_scale = half_width; + pl.camera.focal_point = (cx, cy, 0)`. +- For mesh-only views, drop the field and draw `edges` on white, `line_width≈0.7`. +- Real corruption vs render artifact: apparent "holes / lumps" are often a + mesh-overlay/low-res artifact. Before calling adaptation broken, CHECK the + field's value range is bounded and count folded elements (negative cell area) + programmatically — do NOT diagnose from a render alone. +- **Overlay the feature you're refining to** (a fault trace, an interface): draw it + as a red `pv.PolyData` line over the mesh. Without it you cannot tell whether the + refinement sits ON the feature or has drifted off it (a real failure mode — see + the `adaptive-meshing` skill). Read the geometry from the run manifest so any run + renders the same way. + +## Adaptive / long runs + +- **Render each checkpoint as it lands**, not just the last frame: arm a Monitor that + polls for new `run.mesh.NNNNN.{xdmf,h5}` and emits the index → render on each + event. A completion-only watch leaves you blind for a multi-hour (e.g. TI) run. +- The per-step mesh GEOMETRY must have been written (`write_timestep(..., + meshUpdates=True)`) or you'll render deformed fields on the stale step-0 mesh. + Load the per-step `run.mesh.NNNNN.h5` as the mesh, then `read_timestep` the vars. + +## Velocity + +Same pattern; use **streamlines, not glyphs**. Build a pv mesh for V, add +`pv_mesh.streamlines(...)` or evaluate V on a line seed. Magnitude with the same +white-bg / lighting=False rules. + +## Quantities to judge a convection run (not just pretty pictures) + +- `vrms` from `uw.function.evaluate(V.sym.dot(V.sym), mesh.X.coords)` → the clean + kinetic-energy indicator (more reliable than nodal boundary metrics). +- Surface heat flux Nu via `uw.maths.BdIntegral` on the Upper boundary. +- Mesh quality: fault/bulk nearest-neighbour spacing RATIO (cKDTree) for refinement; + folded-element count + min cell area for tangling. + +## Templates in this skill + +- `render_field.py` — single/`--all`-steps T+mesh render of a run directory. +- `render_field_streamlines.py` — T colormap + mesh + **V streamlines** (sparse + seeds, thin lines, short integration so weak/closed cells read clearly, not + black spiral-blobs). Use for convection. `--tag --all`. +- `zoom_compare.py` — side-by-side cropped mesh+field for N runs at one step. + +Copy these into the run's `scripts/` (or run in place), point `--sim-dir` at the +run, and adjust the field/variable names. They already encode every rule above. + +## Related memory + +`feedback_use_uw_pyvista_visualisation.md`, `feedback_pyvista_viz_pattern.md`, +`feedback_render_all_steps.md`, `project_adaptation_corruption_was_render_artifact.md`. diff --git a/docs/developer/index.md b/docs/developer/index.md index d74b816d4..08b180d12 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -44,6 +44,26 @@ same topic are reference or historical material subordinate to the governing doc | Docstring format | NumPy/Sphinx with RST `:math:` — Charter §6 and the [Style Guide docstring section](UW3_Style_and_Patterns_Guide.md) | | Documentation file format | MyST Markdown (`.md`) for Sphinx — CLAUDE.md "Documentation Requests" section | +## Capability guides + +Curated guidance that the code cannot state about itself: which scheme to +choose, which boundary treatment to prefer, how to make a hard solve +converge. Each is one page with front matter naming the families it applies +to, so `uw.capabilities()`, `uw.systems.Stokes.view()` and the MCP server list +it beside the family; the AI skills in `.claude/skills` are symlinks to these +pages, never copies. A change to a family is reviewed against its guides. + +| Guide | Applies to | +|---|---| +| [Transport schemes](guides/transport-schemes.md) | advection-diffusion, Navier-Stokes, the history schemes | +| [Boundary-condition rulings](guides/boundary-condition-rulings.md) | every solver | +| [Nonlinear solver recipe](guides/nonlinear-solver.md) | viscoplastic Stokes | +| [Plasticity solvers](guides/plasticity-solvers.md) | viscoplastic and VEP Stokes | +| [Adaptive meshing](guides/adaptive-meshing.md) | moving-mesh convection | +| [Adapt-on-top faults](guides/adapt-on-top-faults.md) | fault models on adapted meshes | +| [Free-surface convection](guides/free-surface-convection.md) | free-surface Stokes | +| [Visualisation](guides/uw-visualisation.md), [CeTZ figures](guides/cetz-figures.md) | figures | + ## Documentation Structure This documentation is organized into focused sections: @@ -140,6 +160,15 @@ guides/state-as-dataclass guides/BINDER_CONTAINER_SETUP guides/hpc-cluster-setup guides/mpi-hang-supervision +guides/transport-schemes +guides/boundary-condition-rulings +guides/nonlinear-solver +guides/plasticity-solvers +guides/adaptive-meshing +guides/adapt-on-top-faults +guides/free-surface-convection +guides/uw-visualisation +guides/cetz-figures ``` ```{toctree} diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 6b2499650..44424922a 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -759,8 +759,16 @@ def describe_class(cls, depth=4): "text": None, "units": getattr(attr, "units", None), "description": (getattr(attr, "description", "") or "").strip(), "where": []}) + facts = {} + try: + from underworld3.utilities.capabilities import guides_for + linked = guides_for(cls.__name__) + if linked: + facts["guides"] = linked + except Exception: + pass return record("constitutive_model_family", cls.__name__, doc.split("\n")[0], - documentation=doc or None, terms=terms or None) + documentation=doc or None, facts=facts or None, terms=terms or None) def describe(self, depth=4): """What this constitutive model is, as data: its parameters as terms, diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 6e1255330..7994f2e9f 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -1533,6 +1533,13 @@ class SolverBaseClass(uw_object): conditions.append({"mechanism": method, "type": method[4:-3].replace("_", " "), "boundary": "any", "latex": None, "text": (getattr(fn, "__doc__", "") or "").strip().split("\n")[0] or None}) + try: + from underworld3.utilities.capabilities import guides_for + linked = guides_for(cls.__name__, *public) + if linked: + facts["guides"] = linked + except Exception: + pass return record("solver_family", cls.__name__, doc.split("\n")[0], documentation=doc or None, facts=facts, forms=forms or None, terms=terms or None, conditions=conditions or None, terms_declared=bool(terms)) diff --git a/src/underworld3/mcp/__init__.py b/src/underworld3/mcp/__init__.py index 6832b7789..2a4bd460c 100644 --- a/src/underworld3/mcp/__init__.py +++ b/src/underworld3/mcp/__init__.py @@ -334,5 +334,30 @@ def uw_capability(name: str, format: str = "markdown") -> str: return f"error: {exc}" +@server.tool(name="uw_guides", annotations=_READ_ONLY) +def uw_guides() -> str: + """The capability guides in this checkout: curated guidance the code + cannot state about itself (which transport scheme, which boundary + treatment, how to make a hard solve converge), each with the families + it applies to. uw_guide reads one.""" + from ..utilities.capabilities import guides + found = guides() + if not found: + return "no guides found: the server is not running inside an Underworld3 checkout (set UW_DOCS to its docs directory)" + return _yaml([{k: v for k, v in g.items() if k != "path"} for g in found.values()]) + + +@server.tool(name="uw_guide", annotations=_READ_ONLY) +def uw_guide(name: str) -> str: + """One capability guide in full, as Markdown. name is from uw_guides, + such as transport-schemes, boundary-condition-rulings or + nonlinear-solver.""" + from ..utilities.capabilities import guide_text, guides + text = guide_text(name) + if text is None: + return f"error: no guide named {name!r}; guides are {sorted(guides())}" + return text + + def main(): server.run(transport="stdio") diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index fa6645ca6..baf24de11 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -572,6 +572,13 @@ def describe_class(cls, depth=4): facts[f"default {key}"] = params[key].default except (TypeError, ValueError): pass + try: + from underworld3.utilities.capabilities import guides_for + linked = guides_for(cls.__name__) + if linked: + facts["guides"] = linked + except Exception: + pass return record("history_family", cls.__name__, doc.split("\n")[0], documentation=doc or None, facts=facts) def describe(self, depth=4): diff --git a/src/underworld3/utilities/capabilities.py b/src/underworld3/utilities/capabilities.py index 2ef7a5d49..073875fa8 100644 --- a/src/underworld3/utilities/capabilities.py +++ b/src/underworld3/utilities/capabilities.py @@ -40,6 +40,76 @@ def families(): return out +def guides_directory(): + """The checkout's ``docs/developer/guides``, found upward from the + working directory or from ``UW_DOCS``; ``None`` outside a checkout.""" + import os + candidates = [] + if os.environ.get("UW_DOCS"): + candidates.append(os.path.join(os.environ["UW_DOCS"], "developer", "guides")) + here = os.path.abspath(os.getcwd()) + while True: + candidates.append(os.path.join(here, "docs", "developer", "guides")) + parent = os.path.dirname(here) + if parent == here: + break + here = parent + for c in candidates: + if os.path.isdir(c): + return c + return None + + +def guides(): + """The capability guides in the checkout: ``{name: {name, description, + families, kind, path}}`` read from each page's front matter. Empty + outside a checkout.""" + import glob + import os + import yaml + directory = guides_directory() + out = {} + if directory is None: + return out + for path in sorted(glob.glob(os.path.join(directory, "*.md"))): + with open(path, encoding="utf-8") as handle: + head = handle.read(4000) + if not head.startswith("---"): + continue + parts = head.split("---", 2) + if len(parts) < 3: + continue + try: + meta = yaml.safe_load(parts[1]) or {} + except yaml.YAMLError: + continue + if not isinstance(meta, dict) or "families" not in meta: + continue + name = str(meta.get("name") or os.path.splitext(os.path.basename(path))[0]) + out[name] = {"name": name, "description": str(meta.get("description") or ""), + "families": [str(f) for f in (meta.get("families") or [])], + "kind": str(meta.get("kind") or "guide"), "path": path} + return out + + +def guides_for(*names): + """The names of the guides whose ``families`` include any of ``names`` + (a public name or a class name).""" + wanted = {str(n) for n in names if n} + return [g["name"] for g in guides().values() if wanted & set(g["families"])] + + +def guide_text(name): + """The body of one guide, front matter removed, or ``None``.""" + g = guides().get(name) + if g is None: + return None + with open(g["path"], encoding="utf-8") as handle: + text = handle.read() + parts = text.split("---", 2) + return parts[2].lstrip("\n") if text.startswith("---") and len(parts) == 3 else text + + def _summary_row(name, description): """A family reduced to what a catalogue line needs.""" facts = dict(description.get("facts") or {}) @@ -51,6 +121,9 @@ def _summary_row(name, description): facts["given"] = [t["name"] for t in description["terms"]] if description.get("conditions"): facts["conditions"] = [c.get("mechanism") for c in description["conditions"]] + linked = guides_for(name, description.get("name")) + if linked: + facts["guides"] = linked return record(description.get("kind", "family"), name, description.get("summary", ""), facts={"class": description.get("name"), **facts}) @@ -91,5 +164,9 @@ def family(name): if cls is None: cls = next((c for c in members.values() if c.__name__ == name), None) if cls is not None: - return cls.describe_class() + d = cls.describe_class() + linked = guides_for(name, cls.__name__) + if linked: + d.setdefault("facts", {})["guides"] = linked + return d return None diff --git a/tests/test_0030_capability_guides.py b/tests/test_0030_capability_guides.py new file mode 100644 index 000000000..b727c9ff3 --- /dev/null +++ b/tests/test_0030_capability_guides.py @@ -0,0 +1,62 @@ +"""Capability guides are documentation, and nothing else copies them. + +Each curated guide is one page under docs/developer/guides with front +matter naming the families it applies to. The AI skills in .claude/skills +are symlinks to those pages, uw.capabilities() lists a guide beside its +family, and a family's class-level description names its guides. This +file enforces the single source: a copied skill, a guide without front +matter, or a family no class carries, fails here. +""" +import os +import pathlib + +import pytest +import yaml + +import underworld3 as uw +from underworld3.utilities.capabilities import families, guides, guide_text, guides_for + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + +ROOT = pathlib.Path(__file__).resolve().parent.parent +GUIDES = ROOT / "docs" / "developer" / "guides" +SKILLS = ROOT / ".claude" / "skills" + + +def _front_matter(path): + text = path.read_text(encoding="utf-8") + assert text.startswith("---"), f"{path.name}: no front matter" + return yaml.safe_load(text.split("---", 2)[1]) or {} + + +def test_every_skill_is_a_symlink_into_the_guides(): + skills = sorted(p for p in SKILLS.glob("*/SKILL.md")) + assert skills, "no skills found" + for skill in skills: + assert skill.is_symlink(), f"{skill} is a copy; make it a symlink into docs/developer/guides" + target = (skill.parent / os.readlink(skill)).resolve() + assert target.parent == GUIDES.resolve() and target.exists(), (skill, target) + + +def test_every_guide_names_real_families(): + known = {name for group in families().values() for name in group} + known |= {cls.__name__ for group in families().values() for cls in group.values()} + found = guides() + assert {"transport-schemes", "boundary-condition-rulings", "nonlinear-solver"} <= set(found) + for name, g in found.items(): + meta = _front_matter(pathlib.Path(g["path"])) + assert meta.get("name") == name and meta.get("description"), name + unknown = set(g["families"]) - known + assert not unknown, f"guide {name} names families no class carries: {sorted(unknown)}" + + +def test_families_list_their_guides(): + assert "boundary-condition-rulings" in guides_for("Stokes") + assert "transport-schemes" in guides_for("SemiLagrangian") + stokes = uw.systems.Stokes.describe_class() + assert "nonlinear-solver" in stokes["facts"]["guides"] + cat = uw.capabilities("solvers") + row = next(c for c in cat["children"][0]["children"] if c["name"] == "Stokes") + assert "boundary-condition-rulings" in row["facts"]["guides"] + text = guide_text("transport-schemes") + assert text.startswith("# Transport schemes") and guide_text("nothing") is None