ViscoPlasticFlowModel + split-node fault path + consistent tangent: velocity PC setup dies with 'Computed maximum singular value as zero' (error 77)

[lmoresi]

On feature/fault-split-node (PR #502): a Stokes solve with
ViscoPlasticFlowModel, a split-node fault BC (add_fault_bc, so the
solve goes through the rotated path / solve_with_fault), and a
consistent tangent fails during the rotated fieldsplit-Schur KSP setup:

PETSc has generated inconsistent data
Computed maximum singular value as zero
error code 77, in KSP.solve via rotated_bc._solve_rotated_iterative

Facts from the isolation (box + one frictionless split fault, h = 0.03):

• consistent_jacobian = True with picard>0 correctly raises the
  guided NotImplementedError (use "continuation").
• consistent_jacobian = "continuation", picard=3: error 77 above —
  and it fails even with uniform yield_stress = 1e4, i.e. a
  problem that is numerically pure viscous. A spatially varying (masked)
  yield stress makes no difference: both fail identically.
• Default (Picard) tangent with picard=8: both the uniform and the
  masked yield-stress cases solve fine.
• Plain ViscousFlowModel through the same fault path: fine (all the
  fault teaching examples).

So the failing intersection is specifically {ViscoPlasticFlowModel,
rotated fault path, frozen-Picard-then-Newton compilation}. Since the
1e4-yield case is numerically linear, the zero/degenerate operator the
velocity PC sees presumably comes from the continuation machinery's
frozen-tangent compilation of the VP flux rather than from the physics.

Reproducer (worktree env): build a split-fault box, set
ViscoPlasticFlowModel with yield_stress = 1e4,
consistent_jacobian = "continuation", add_fault_bc(0, ...),
solve_with_fault(picard=3).

Workaround in use: the Picard tangent (default) with a few extra
iterations, adequate for mild localised yielding (junction plugs).
The consistent-tangent path matters once yielding is stiff/global
(the notch-class problems on faulted meshes).

Underworld development team with AI support from Claude Code

[lmoresi]

Decomposition run — the split-node formulation is exonerated:

• split mesh + ViscoPlasticFlowModel + plain Dirichlet walls, standard
  SNES path: solves fine.
• ORDINARY box + ViscoPlasticFlowModel + rotated free-slip WALLS +
  consistent_jacobian="continuation": fails with the identical
  "Computed maximum singular value as zero" — no fault in the problem.
• Split-fault case: fails identically (as implied by the above).

So this is a latent rotated-path bug (frozen-tangent / continuation
compilation for the VP model), present for plain rotated free-slip and
merely surfaced by the fault work; VP + rotated walls had not been
exercised together before (the viscoplastic campaigns ran Dirichlet
walls). Likely the same family as the #457/#493 frozen-Picard issues:
the freeze machinery producing a degenerate operator for this model's
flux rather than an unfrozen one — the zero maximum singular value
suggests the velocity block of the frozen Jacobian assembles to zero.
Fix belongs in rotated_bc's tangent staging; one fix covers walls and
faults. Reproducer: box, VP with uniform yield_stress = 1e4 (numerically
pure viscous), add_rotated_freeslip_bc on two walls,
consistent_jacobian = "continuation", solve.

Underworld development team with AI support from Claude Code

[lmoresi]

Root cause found, fixed on feature/fault-split-node (181d618), and the story is bigger than the title: this was never about the rotated path or the faults.

Root cause. Differentiating an unwrapped strain-rate invariant produces half-integer powers of (grad v : grad v) whose value or derivative is 0/0 at a state of rest. Every consistent-tangent assembly with edot = 0 anywhere therefore filled the operator with NaN — measured directly: J(0) norm = nan for a ViscoPlastic model at ANY yield stress (including 1e4, where the problem is numerically pure viscous), healthy at any smooth state. GAMG's eigen-estimator reports the NaN comparisons as "Computed maximum singular value as zero". Standard SNES path, rotated free-slip and split-node faults all affected identically; reproduces on development.

Why the existing protections did not cover it (all measured):

1. The Layer-1 automatic cold-start Picard injection fires correctly and still fails on BOUNDARY-DRIVEN problems: one nrichardson sweep moves boundary data a single element layer, so the deep interior stays at exactly zero strain rate. Body-force-driven problems fill F(0) everywhere — which is why the viscoplastic campaigns validated cleanly and every fault problem (wall-driven) died.
2. The "continuation" tangent's alpha = 0 phase is not Picard-safe: the blended kernel evaluates the Newton branch pointwise, and IEEE 0*NaN = NaN.
3. Decisively: a rigidly-translating stuck region has edot = 0 at the CONVERGED solution. The singular state is physics, not a start-up artifact — no warm-start policy can make it unreachable.

Fix. _jacobian_unwrap now adds 1e-36 under every half-integer power whose argument carries unknowns — the derivative's removable-singularity limit, implemented. Jacobian sources only: the residual never routes through that function and the default Picard tangent never calls it, so both remain bit-identical by construction. Perturbation at any resolvable strain rate: under one part in 1e24. The Layer-1 warm-up is retained as the convergence aid it was designed to be, with its comment corrected.

Verification. All previously-failing reproducers (standard / rotated wall / split fault, True / "continuation") now solve; the cold Jacobian of the never-yielding probe equals its smooth-state Jacobian to all printed digits; 53 tests green across yield-homotopy, rotated free-slip and the fault suites; new regression test test_1067_newton_cold_start.

This should be cherry-picked to development ahead of the PR #502 merge if that merge is delayed — the bug bites any consistent-tangent viscoplastic solve that is boundary-driven or contains rigid regions.

Underworld development team with AI support from Claude Code