Stokes uu_G3 transposed tangent + TI un-frozen Picard: fix issue #457 at source

[lmoresi]

Fixes #457 at source. Two independent defects, both hidden by major-symmetric tangents until the transverse-isotropic Newton tangent — the first in-repo tangent without major symmetry — exposed them.

1. The Stokes velocity-block Jacobian was the major transpose of the true tangent

sympy.derive_by_array is dx-first, but SNES_Stokes_SaddlePt translated it with permutedims((0,2,1,3)) — correct only under the dx-last convention that petsc-jacobian-layout.md wrongly recorded. The assembled uu_G3 therefore held dF1[gc,dg]/dL[fc,df] in PETSc's [fc,gc,df,dg] slot. Frozen constitutive tensors and isotropic Newton terms are major-symmetric, so nothing in the suite could see it; the 2026-04-21 layout audit signed the Stokes path off against the mis-stated convention at the same time SNES_Vector was being fixed for the identical bug.

All uu, up, and natural-BC (including Nitsche G2/G3) blocks now use the explicit-index loops SNES_Vector and SNES_MultiComponent already use. The block-constraint boundary G0 keeps its form — shape-protected, and r(n⊗n) is symmetric.

2. The TI c-tensor builds silently un-froze the Picard tangent

TransverseIsotropicFlowModel._build_c_tensor baked the unwrapped .sym contents of η₀/η₁ into the tensor, exposing their strain-rate dependence to sympy.diff — so "Picard" for TI was the (transposed) Newton derivative, never a frozen-coefficient tangent. It now bakes the wrapped Parameter atoms exactly like the isotropic path. TransverseIsotropicVEPFlowModel gets the same discipline plus a persistent _eta1_yield_eff UWexpression container (the ViscoPlasticFlowModel pattern) so the yield-limited weak-plane viscosity freezes under Picard and is recovered by _jacobian_unwrap under Newton.

Measured (unpatched, this branch, on the #457 reproduction)

quantity	before	after
Picard velocity-block asymmetry	1.089e-01	4.0e-17
Newton -snes_test_jacobian rel	1.09e-01	1.3e-4 (FD-limited; isotropic control 1.06e-4)
Newton convergence	Picard pace / stall	9 iterations

Tests

tests/test_1066_stokes_jacobian_layout.py, validated both ways (all three fail on the pre-fix build at their primary assertions — 0.313 rel FD mismatch, velocity-gradient dependence throughout — and pass on this branch):

1. Assembled uu_G3 vs a central-difference oracle built from the residual flux, with a built-in guard that the oracle tangent is genuinely major-asymmetric (the test cannot pass by symmetry).
2. TI Picard uu_G3 has zero velocity-gradient dependence (frozen) and exact major symmetry; Newton positive control.
3. TI-VEP c-tensor coefficient freezing with a Newton-unwrap positive control.

Full level_1 and tier_a gate: 554 passed, 0 failed.

Docs: petsc-jacobian-layout.md corrected — dx-first convention, a warning recording the earlier mis-audit, the Stokes migration, and a ban on permutedims for new Jacobian code.

Out of scope, noted for the record: #463 (TI-VEP viscosity property discards the yield-limited η₁ — orthogonal property-level defect) and the isotropic VEP _build_c_tensor's own .sym use (its flux path runs through stress(), not this tensor, and its tangents are symmetric either way; it should be tidied with #463).

Underworld development team with AI support from Claude Code

[Copilot]

🟡 Not ready to approve

The natural-BC gradient-flux up_G2/up_G3 Jacobian currently differentiates the raw boundary flux instead of the Newton-consistent _jacobian_source() version, which can reintroduce the very kink/unwrap inconsistency the new path is meant to avoid.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

This review doesn't count toward merge requirements. Sign up for the private preview to control whether Copilot approvals count.

Pull request overview

Fixes the root causes of issue #457 by correcting Stokes Jacobian block assembly to match PETSc’s [fc,gc,df,dg] layout under SymPy’s dx-first derive_by_array convention, and by restoring the intended Picard “frozen-coefficient” contract for transverse-isotropic (TI / TI-VEP) tangents via wrapped-atom coefficient handling.

Changes:

• Replaces derive_by_array + permutedims in SNES_Stokes_SaddlePt with explicit-index Jacobian construction for uu, up, and natural-BC blocks to avoid major-transpose assembly.
• Updates TI and TI-VEP c-tensor coefficient construction to keep viscosity/yield coefficients inside wrapped atoms under Picard, with explicit Newton unwrapping.
• Adds targeted regression tests and corrects/strengthens the Jacobian layout documentation to match the dx-first convention.

File summaries

File	Description
tests/test_1066_stokes_jacobian_layout.py	Adds regression tests for Stokes uu_G3 layout vs FD oracle and for TI/TI-VEP Picard coefficient-freezing behavior.
src/underworld3/cython/petsc_generic_snes_solvers.pyx	Migrates Stokes + natural-BC Jacobian assembly to explicit-index loops matching PETSc layout (removes permutedims path for these blocks).
src/underworld3/constitutive_models.py	Fixes TI/TI-VEP c-tensor construction to use wrapped Parameter atoms (and adds persistent yield-limited η₁ container) so Picard tangents remain frozen.
docs/developer/subsystems/petsc-jacobian-layout.md	Corrects SymPy convention documentation (dx-first) and records the Stokes migration/ban on new permutedims Jacobian code.

Review details

Suppressed comments (1)

src/underworld3/cython/petsc_generic_snes_solvers.pyx:7583

• Similarly, up_G3 for the natural-BC gradient term differentiates bd_F1 instead of the Newton-consistent bd_F1_jac that was just constructed. If _jacobian_source() is smoothing/unwrapping bd_F1 for Newton, up_G3 should use the same source to keep the block Jacobian consistent with the residual linearisation policy.

                                G3[fc * dim + df, dg] = sympy.diff(
                                    bd_F1[fc, df], Gp[0, dg]
                                )

• Files reviewed: 4/4 changed files
• Comments generated: 1
• Review effort level: Lite

We're testing this review assessment. Please use 👍 or 👎 to tell us if it's correct.

[lmoresi]

Adversarial review — live probes at the PR head (isolated worktree, own build)

We attempted to break both fixes independently. Verdict: clean, two nits (response commit below).

Findings

1. MINOR — petsc-jacobian-layout.md claimed "no remaining in-repo derive_by_array + permutedims Jacobian path" while the block-constrained boundary uu_G0 still used it (cbc_permutation, content r·(n⊗n) — symmetric and shape-protected, so numerically harmless). Fixed in response commit b0cc240 by converting that site to the explicit loops; the claim is now literally true (0 live permutedims in the file). Retest: layout + constrained free-slip + constrained SolCx + Stokes, 25 passed.
2. MINOR — two G1 loops named the trial-side derivative index df where the layout doc uses dg. Cosmetic; renamed in response commit b0cc240.
3. NOT-PR — tests/parallel/test_1063_constrained_freeslip_parallel.py fails at np=2 on BOTH pre-fix development and the PR build with byte-identical iso values (0.59828... vs golden 0.61945...; TI differs between builds only at 6e-13). Pre-existing golden/environment mismatch — filed as ptest 1063 (constrained free-slip np2): iso golden value mismatch is environment-dependent and PRE-existing #495.

Attacks that failed (independently measured)

• 3D TI Newton, tilted director (1,2,3)/√14: symbolic FD-oracle on uu_G3 rel 6.7e-11 (the transposed form scores 1.63e-1); PETSc -snes_test_jacobian 3D TI 7.9e-08 vs isotropic 3D control 6.9e-08 — same order.
• 2D reproduction: pre-fix -snes_test_jacobian 0.16–0.20 → PR build 5.5e-07 (sharper than the PR body's 1.3e-4).
• sympy convention independently re-verified dx-FIRST — the doc correction and transpose diagnosis stand.
• Flatten orders of the reshaped up_*/Nitsche blocks verified against the fe.c flat-index formula; no in-repo solver puts p/grad-p in F0, so the old wrong-ordered-but-zero blocks never fired; pu_* blocks are correct because Nc_f=1 makes dx-first coincide.
• Navier-Stokes: F0's u-dependence is diagonal (ρc₀/Δt·I) — the old transpose was invisible there; no in-repo asymmetric G0; test_0610 passes. A user-supplied asymmetric u-dependent body force previously got a transposed G0 — this fix improves that case; converged answers unchanged.
• TI-VEP container: hybrid mode keeps the container out of the ETD tensor; per-instance _uw_id (no cross-model collision); yield_stress updates flow through without a stale tangent; is_constant_expr = False (the JIT will not freeze it as a constant); .units behavior identical to the established _plastic_eff_viscosity pattern.
• Test vacuity: the FD oracle is lambdify + central differences (no shared sympy.diff machinery); the asymmetry guard has 16x margin (0.163 vs 0.01); 3/3 pass on the PR build, 3/3 fail at their primary assertions on pre-fix development (uu_G3 rel 3.13e-01).
• Charter sweep clean; np=2: test_1010 6 passed; Nitsche/rotated/constrained 38 passed; VEP regression 8 passed; full level_1+tier_a gate on the PR build: 531 passed, 0 failed.