Merge branch 'development' into bugfix/topography-penalty-633

Author: lmoresi

.claude/skills/nonlinear-solver/SKILL.md

@@ -77,7 +77,7 @@ case felt like whack-a-mole. Check them first.
 
 | Trap | Symptom | Fix |
 |---|---|---|
-| Consistent Newton makes the velocity block **non-symmetric**; a Chebyshev/Richardson MG smoother assumes SPD | smoother diverges / stalls → `DIVERGED_LINEAR_SOLVE` or an endless grind | **Now the default** — the FMG bundle ships `mg_levels_ksp_type=gmres` + `pc_type=sor` + `norm_type=none`. Only an issue if you override it, or on GAMG (which uses PETSc's chebyshev default) |
+| **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 |

docs/advanced/fault-networks.md

@@ -15,24 +15,113 @@ net = uw.meshing.FaultNetwork(
     [("Main", main_pts), ("Splay", splay_pts), ("Cross", cross_pts)],
     hierarchy=["Main", "Splay", "Cross"])   # seniority order
 
-mesh = net.prepare(h=0.006).build()          # junctions -> mesh -> split
+mesh = net.prepare(h=0.006).build(width=0.01)  # junctions -> mesh -> split
 
 v = uw.discretisation.MeshVariable("V", mesh, 2, degree=2)
 p = uw.discretisation.MeshVariable("P", mesh, 1, degree=0,
                                    continuous=False)
 stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
-stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel
-stokes.constitutive_model.yield_mode = "min"
-stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0
-stokes.constitutive_model.Parameters.yield_stress = \
-    net.damage_yield(v, dial=0.05)           # the junction glue
-stokes.consistent_jacobian = True
-net.apply_contact(stokes)                    # no-opening pairs, all pieces
+stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
+stokes.constitutive_model.Parameters.shear_viscosity_0 = \
+    net.junction_patch(eta_0=1.0)            # the junction glue (linear)
+net.apply(stokes)                            # no-opening pairs, all pieces
 # ... wall boundary conditions ...
 info = net.solve(stokes)
 print(net.slips(stokes))                     # peak slip per piece
 ```
 
+## One specification, two realisations
+
+A fault is specified once — a trace, its rank in the hierarchy, and the
+properties it carries — and then *realised*. Which realisation you get
+is a keyword on `build`, not a different set of calls:
+
+```python
+net.prepare(h=0.006)
+mesh = net.build(width=0.01)                      # cut, node-pair contact
+mesh = net.build(width=0.002, realisation="ti")   # volumetric weak plane
+net.apply(stokes, eta_1=0.01)                     # eta_1: TI only
+```
+
+Both realisations place the same ribbon band along the same prepared
+pieces, so the cells are identical and results from the two may be
+compared directly. The band is meshed around the trace's own points and
+segments, which become mesh vertices and edges, so **the mesh can be cut
+whatever the width is** — measured complete, with exact vertex
+coincidence, down to a band a tenth of the background element size. The
+realisation is a free choice, not something the mesh grants or refuses.
+
+What differs is what `width` *means*. For the split it is a resolution
+parameter: the band exists to give the cut its own vertices, and its
+thickness is not physics. For the weak plane it is constitutive — the
+layer thickness that sets the slip rate through `V = 2 e_nt w` — so it
+wants two or three elements across it. That is the whole asymmetry, and
+it is about the rheology rather than the mesh.
+
+`slips()` reports each realisation in its own quantity: the tangential
+jump between the two nodes of a cut pair, or the jump in tangential
+velocity across the layer, sampled one half-width plus a cell either
+side of the spine. Both are the fault's own throughput; a probe placed
+further out reads the surrounding flow as well and over-reads short
+strands. The gauge is rank-local: each rank reports only the probe
+pairs it owns and omits a piece it holds no pair of, so a
+max-reduction across ranks recovers the network's answer. This
+matters because `evaluate` answers for any point it is handed,
+extrapolating from the nearest local cell when the point is not in
+the local mesh — a band-less rank would otherwise report a far-field
+extrapolation as the band's slip.
+
+`build(width=None)` keeps the older no-band path — graded refinement
+cut directly. It is split-only, and its mesh is not the one a weak
+plane would use, so do not compare across that choice.
+
+**The band is material, not scaffolding.** It is easy to read the band
+as something the weak plane needs and the split merely tolerates. It is
+not: the band is a meshed region of material *around* the fault, and a
+segmented fault does its interesting work exactly there — at the strand
+tips, and in the ligaments where one cut stops short of the next. Damage
+in those places needs cells to live in, and the band is where they are.
+`net.band` is the mask, `net.footprints` the per-strand ones, in either
+realisation.
+
+`net.band_yield(tau_y)` gives a rheology for the whole band: von
+Mises yield confined to it, everything outside set far too strong to
+yield. Read the next paragraph before using it as glue.
+
+```python
+stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel
+stokes.constitutive_model.Parameters.yield_stress = net.band_yield(4.0)
+stokes.consistent_jacobian = True
+```
+
+A released fault flank sits far below `tau_y` and is untouched, so
+the breakdown appears where the mechanics puts it — but that is tips
+and bends as much as joints. Measured on the S-fault rig: at the
+strength that repairs a stepover, more than half the yielded band cells
+were on strand flanks and free tips, and one step weaker the whole
+main strand had become a weak fault. A uniform threshold cannot pick
+out the welds alone, because the stress concentration at a weld is not
+far enough above the tip and bend concentrations. `band_yield` is a
+damage model for the band; the junction glue is `junction_patch`.
+
+The two realisations' interface parameters correspond, which is worth
+keeping in view when comparing them: the zero-thickness limit of a band
+of viscosity `eta_band` and width `w` is an interface viscosity
+`eta_f = eta_band / w`, which is the `conds` argument of
+`add_fault_bc`. The weak plane's `V = 2 e_nt w` is precisely what the
+contact replaces with a genuine slip rate.
+
+**Properties belong to the fault.** `net.surface(name)` returns the
+retained {class}`~underworld3.meshing.surfaces.Surface` for a piece.
+Friction, accumulated slip, a damage state live there, on the fault,
+and outlive any one realisation of it:
+
+```python
+main = net.surface("Main")
+friction = main.add_variable("mu", size=1)
+friction.data[:] = 0.6
+```
+
 ## The recipe, and why each piece is the way it is
 
 **Hierarchy.** At an X crossing the senior fault runs through and the
@@ -48,17 +137,70 @@ same answer when the junction patch is refined 2x. Make the join as
 small as the mesh allows; buy fidelity with elements, not physical
 size.
 
-**The glue.** `damage_yield` places a compact viscoplastic plug at
-each junction: yield `dial * (1 + 2 * edot_II)` inside, effectively
-rigid outside, sharp `Piecewise` boundaries. The strength and the
-rate-regularisation move together on ONE dial (separating them makes
-the solve harsh without making the zone weaker). Zone stress is
-proportional to the dial down to a ~100x viscosity contrast with
-Newton-from-cold still converging — the compact plug conditions like a
-hole, not like a thin weak layer, so the classic thin-inclusion
-Schur breakdown never appears. `dial=0.05` is near-invisible in the
-stress field at unchanged cost; `dial=0.01` reaches the transmission
-ceiling of an inviscid plug at roughly double cost.
+**The glue, and where it goes.** The split only goes wrong at the
+joints: away from them the cut *is* the target every volumetric
+representation converges to, and adding weakness along a whole strand
+makes the fault over-weak in a way that depends on the band width. So
+the glue is placed, not found. `junction_cells()` reads the places off
+the mesh itself: the ribbon (the band with its extrapolated margins)
+is everything the weak-plane realisation would treat as fault, the cut
+chains are what the split sliced, and a band cell whose nearest spine
+point lies in a piece's margin *and* which sits inside a second
+piece's ribbon is where two pieces meet without being joined — a
+kissing branch, an abutting pair, the intact bridge of a stepover.
+Free tips are excluded on purpose: a margin that runs into intact
+material is a tip, and damage there lengthens the fault instead of
+joining it (measured: with the free tips included, nearly every
+yielded cell was at a tip and the main strand grew 1-6% longer in
+slip). The cells are dilated by one vertex ring, and that ring is
+not optional: the weld's stiffness lives in the intact material
+around the two tips, and the bare junction cells recover only a
+fifth to a quarter of the deficit even when fully plastic.
+
+Two pieces that continue one another along a line — an abutting pair,
+a stepover's continuation — are placed on **one spine**: two ribbons
+laid along the same line interleave their vertices into sliver cells
+(measured: 7800 cells below 1e-6 in area, and the velocity solve
+five times slower). `build()` groups such pieces (end tangents within
+25 degrees, the far start within half a width of the line, within the
+margins' reach), bridges the gap with spine vertices at the local rung,
+and cuts each piece at its own ends; the gap is spine the split does
+not cut, which is exactly what the junction rule reads.
+
+For the rule to see a joint, the ribbons have to meet across it.
+`build()` sees to that: at an end that sits on a prepared junction the
+tip margin is extended until the ribbon reaches the other piece's cut,
+so the whole ligament lies in both ribbons; free tips keep the default
+margin. An abutting pair that `prepare()` did not record as a junction
+(a gap wider than the ligament) is covered as far as the default
+margins overlap — a gap wider than that is two faults, and stays
+welded, which is what a gap of intact rock means.
+
+`junction_patch(eta_0, ratio=0.01)` then makes those cells weak
+isotropic material, `eta = ratio * eta_0`. A viscosity ratio rather
+than a yield stress, because the joint only has to be broken and a
+ratio needs no stress scale — nothing about the block or the loading
+has to be known to set it. Measured against the two end members on
+the S-fault rig (the fault *longer*, one continuous cut, and the fault
+*cut*, abutting cuts, at two resolutions): the patch recovers
+0.8-0.97 of the continuous fault's transmission across the joint; the
+slip crosses on the cut itself (the segment's pair jump reaches the
+continuous fault's); the rest of the network keeps the split's answer
+(main strand within 2.5%); the weak patch reproduces a fully plastic
+patch on the same cells to 1-2% and is insensitive to the ratio from
+0.01 to 0.001; the solve is linear and costs the split's velocity
+iterations, with only the pressure block noticing the contrast
+(hence 0.01, not smaller). Gluing a joint does change the partition
+between the strands that meet there — a reconnected main line takes
+back slip a through-going branch was carrying past the weld — which
+is the junction working, not the patch leaking.
+
+`damage_yield` is the older glue: a viscoplastic plug of radius
+`max(2.5 h, 1.2 pull)` at each *prepared* junction point, yield
+`dial * (1 + 2 * edot_II)` inside, strength and rate-regularisation on
+one dial, sharp `Piecewise` boundaries. It stays available for
+studies of the glue itself, and it does not see stepover bridges,
+which are not prepared junctions.
 
 **No prescribed reconnection.** Nothing tells the network how to link
 up: the stress lobes of the abutting tips decide. A collinear gap
@@ -120,15 +262,45 @@ directly. Place is currently OPT-IN for networks: on graded
 solve is pathological — an open operator-health work item; on uniform
 bases it is healthy.
 
+`build(width=...)` gives the 3-D network the same finite-width
+contract as 2-D: the margin-expanded patches are thickened by
+`±width/2` into ONE fused band (junctions free), and each un-expanded
+patch is embedded in the band as a conforming mid-surface — so the
+same mesh is cut and split (`realisation="split"`) or left whole for
+the volumetric weak plane (`realisation="ti"`), exactly as in 2-D.
+The honoured footprints are exact for planar patches (band cells
+within the patch's own in-plane outline and half a width of its
+plane), the weak-plane director is the patch normal, and `slips()`
+reports the plane form of each gauge: the tangential pair jump for
+the split, the in-plane velocity jump across the layer for the weak
+plane. Junction glue in 3-D remains `damage_yield`'s tubes about the
+junction segments — `junction_cells` is the 2-D ribbon rule and
+refuses in 3-D. Interior networks only: a band that reaches the
+domain boundary is refused loudly (its embedded mid-surface cannot
+yet be clipped against the boundary).
+
 v1 scope, refused loudly outside it: planar patches (the
 `rim_polygon` contract), convex rims, genuine X crossings (a
 near-miss — close but not crossing — is refused rather than guessed
-at); parallel MULTI-fault splitting (the pairing does not yet migrate
-through redistribution — single faults are parallel-validated).
+at). Multi-fault networks split and solve in parallel: `split_faults`
+redistributes ONCE, keyed on the union of the network's facets, and
+every split then runs with serial topology. Either realisation runs
+its velocity block on the geometric multigrid tail the band's base
+mesh owns (`build` adopts it on the final mesh; `net.solve` says so
+when a solve falls back to algebraic multigrid). Placement is
+gather-first, so the cells gmsh fills into the carved cavity — the
+band and its graded surround — live on one rank; only the base's far
+field is balanced. On the crossing-patches fixture that is 8012 of
+8405 cells on one rank at np=4, and the solve is not faster than
+serial: parallel is a correctness mode for this path, not a speed-up,
+until the placed region is rebalanced.
 
 ## Limitations
 
-- 3-D: planar convex patches, X crossings only, serial (above).
+- 3-D: planar convex patches, X crossings only (above). Finite-width
+  bands are interior-only (an outcropping band is refused), and
+  `junction_cells` is the 2-D ribbon rule — 3-D junction glue is
+  `damage_yield`'s tubes.
 - One damage dial per network in `damage_yield` (per-junction values:
   build the expression with `uw.meshing.damage_zone_yield` directly).
 - Time-dependent damage (wear-in/healing) is study-level for now: see