From 723695de8e9992562e09fb6e5e9b3da99f17bb47 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 3 Aug 2026 14:56:05 +1000 Subject: [PATCH 1/2] fix(solver): Stokes uu_G3 was the major transpose of the true tangent; TI c-tensors un-froze the Picard tangent (#457) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two independent defects, both invisible behind major-symmetric tangents until the transverse-isotropic Newton tangent (the first in-repo tangent WITHOUT major symmetry) exposed them via -snes_test_jacobian and a velocity-block symmetry probe. 1. SNES_Stokes_SaddlePt Jacobian layout. sympy.derive_by_array is dx-FIRST, but the uu/up/natural-BC blocks translated it with permutedims((0,2,1,3)) — valid only for the dx-LAST convention this repo's layout doc wrongly recorded. The result: dF1[gc,dg]/dL[fc,df] landed in PETSc's [fc,gc,df,dg] slot — the major transpose. Every isotropic rheology has a major-symmetric tangent (frozen C is symmetric; the isotropic Newton term is a rank-one e^ (x) e^), so the whole test suite passed over it; the 2026-04-21 audit "verified" the Stokes path against the mis-stated convention while SNES_Vector was being fixed for the same bug. All blocks now use the same explicit-index loops as SNES_Vector/SNES_MultiComponent, writing each entry directly into PETSc's [fc,gc,df,dg] layout. The block-constraint boundary G0 keeps its form (shape-protected and r*(n(x)n) symmetric). F0/F1 sources are shape-normalised (spherical templates deliver (dim,1) F0). 2. TI constitutive freezing. TransverseIsotropicFlowModel._build_c_tensor baked the UNWRAPPED .sym contents of eta_0/eta_1 into the c-tensor, so any strain-rate dependence was visible to sympy.diff and the default (Picard) tangent was silently NOT frozen — it ran the (transposed) Newton derivative instead. Now bakes the wrapped Parameter atoms exactly like the isotropic path. TransverseIsotropicVEPFlowModel gets the same treatment plus a persistent _eta1_yield_eff UWexpression container (the ViscoPlasticFlowModel._plastic_eff_viscosity pattern) so the yield-limited weak-plane viscosity freezes under Picard and is recovered by _jacobian_unwrap under Newton. Measured on the #457 reproduction (annulus-class TI, eta_1/eta_0 = 0.01, shear-thinning eta_0): - Picard velocity-block asymmetry: 1.089e-01 -> 4.0e-17 (machine zero) - Newton -snes_test_jacobian: 1.09e-01 -> 1.3e-4 (FD-limited; isotropic control 1.06e-4); Newton converges in 9 its where the transposed tangent stalled at Picard pace. Docs: petsc-jacobian-layout.md corrected (dx-first convention, warning about the earlier mis-audit, Stokes migration recorded, permutedims banned for new Jacobian code). Tests: test_1066_stokes_jacobian_layout.py — (1) assembled uu_G3 vs an FD oracle from the residual flux with a built-in guard that the oracle tangent is genuinely major-asymmetric; (2) TI Picard G3 has zero velocity-gradient dependence (frozen) + major symmetry, with a Newton positive control; (3) TI-VEP c-tensor coefficient freezing with a Newton-unwrap positive control. All three fail on the pre-fix build at their primary assertions (0.313 rel FD mismatch; L-dependence everywhere) and pass on this build. Full level_1 tier_a gate: 554 passed, 0 failed. Closes #457. Underworld development team with AI support from Claude Code --- .../subsystems/petsc-jacobian-layout.md | 83 +++-- src/underworld3/constitutive_models.py | 41 ++- .../cython/petsc_generic_snes_solvers.pyx | 205 ++++++++++--- tests/test_1066_stokes_jacobian_layout.py | 283 ++++++++++++++++++ 4 files changed, 528 insertions(+), 84 deletions(-) create mode 100644 tests/test_1066_stokes_jacobian_layout.py diff --git a/docs/developer/subsystems/petsc-jacobian-layout.md b/docs/developer/subsystems/petsc-jacobian-layout.md index 5f22636eb..09cd9a539 100644 --- a/docs/developer/subsystems/petsc-jacobian-layout.md +++ b/docs/developer/subsystems/petsc-jacobian-layout.md @@ -61,33 +61,38 @@ pairs like-kind indices. ## Sympy's native convention -`sympy.derive_by_array(F, x)` returns an array with **`F`'s indices -first, then `x`'s**: +```{warning} +An earlier revision of this page stated the opposite convention +("F's indices first, then x's") and audited the Stokes `permutedims` +path as correct on that basis. That was wrong — it is how the issue +#457 transposed-tangent bug survived the 2026-04-21 audit. +``` + +`sympy.derive_by_array(F, x)` is **dx-FIRST**: the returned array has +**`x`'s indices first, then `F`'s** (see the sympy documentation for +`derive_by_array`; verified numerically against finite differences in +the #457 investigation): ```python G3 = sympy.derive_by_array(F1, L) # F1 shape (Nc, dim), L shape (Nc, dim) # → G3 shape (Nc, dim, Nc, dim) -# → G3[i, j, k, l] = ∂F1[i, j] / ∂L[k, l] -# → index order [fc, df, gc, dg] +# → G3[k, l, i, j] = ∂F1[i, j] / ∂L[k, l] +# → index order [gc, dg, fc, df] ``` -This is the same mathematical object PETSc wants, but with the **two -middle axes swapped** relative to PETSc's flat layout. Translating -sympy → PETSc therefore needs the permutation **`(0, 2, 1, 3)`** — swap -axes 1 and 2 — before flattening. - -For `g1` (3D, sympy order `[fc, gc, df]`) and `g2` (3D, sympy order -`[fc, df, gc]`): - -| Array | Sympy shape after `derive_by_array` | Permutation | Final 2D shape | -|-------|--------------------------------------|-------------|----------------| -| `g0` | `(Nc, Nc)` — already `[fc, gc]` | none | `(Nc, Nc)` | -| `g1` | `(Nc, Nc, dim)` — `[fc, gc, df]` | none | `(Nc·Nc, dim)` | -| `g2` | `(Nc, dim, Nc)` — `[fc, df, gc]` | `(0, 2, 1)` | `(Nc·Nc, dim)` | -| `g3` | `(Nc, dim, Nc, dim)` — `[fc, df, gc, dg]` | `(0, 2, 1, 3)` | `(Nc·Nc, dim·dim)` | +Translating sympy → PETSc's `[fc, gc, df, dg]` therefore needs the +permutation **`(2, 0, 3, 1)`** before flattening. The historically-used +`(0, 2, 1, 3)` yields `[gc, fc, dg, df]` — the **major transpose** +`∂F1[gc, dg]/∂L[fc, df]` in the `[fc, gc, df, dg]` slot. That is +identical to the true tangent whenever the tangent has major symmetry +(any frozen major-symmetric `C`, the isotropic Newton term +`η'·(ε̂⊗ε̂)`, linear transverse isotropy) and silently wrong when it +does not (the transverse-isotropic Newton term is `B⊗a` with `B ∦ a` — +issue #457). -The row-major flatten of the 2D form matches PETSc's flat index exactly. +Do not reach for `permutedims` at all — use the explicit construction +below. ## Alternative: explicit construction @@ -142,13 +147,28 @@ An asymmetric F1 like raw `smoothing * L` exposes it immediately. Results of walking every solver in `petsc_generic_snes_solvers.pyx` as of 2026-04-21: -| Solver | File:line | Construction | Status | -|--------|-----------|--------------|--------| -| `SNES_Scalar` | `petsc_generic_snes_solvers.pyx:1255` | no permutation (Nc=1) | ✅ correct — single component makes the swap a no-op | -| `SNES_Vector` | `petsc_generic_snes_solvers.pyx:2013` | explicit per-entry construction | ✅ correct (2026-04-21 migration) | -| `SNES_MultiComponent` | `petsc_generic_snes_solvers.pyx:2919` | explicit per-entry construction | ✅ correct (2026-04-20 fix) | -| `SNES_Stokes_SaddlePt` | `petsc_generic_snes_solvers.pyx:3552` | `permutedims` with `(0, 2, 1, 3)` | ✅ correct — matches PETSc layout | -| `SNES_NavierStokes` | inherits Stokes | inherited | ✅ correct | +| Solver | Construction | Status | +|--------|--------------|--------| +| `SNES_Scalar` | no permutation (Nc=1) | ✅ correct — single component makes the swap a no-op | +| `SNES_Vector` | explicit per-entry construction | ✅ correct (2026-04-21 migration) | +| `SNES_MultiComponent` | explicit per-entry construction | ✅ correct (2026-04-20 fix) | +| `SNES_Stokes_SaddlePt` | explicit per-entry construction | ✅ correct (2026-08 migration, issue #457 — the earlier `permutedims (0, 2, 1, 3)` form assembled the **major transpose** of `uu_G3` and was mis-audited as correct on 2026-04-21 because this page then stated the sympy convention backwards) | +| `SNES_NavierStokes` | inherits Stokes | ✅ correct (inherited) | + +### `SNES_Stokes_SaddlePt` migration (2026-08, issue #457) + +The saddle-point solver used `derive_by_array + permutedims((0,2,1,3))` +for the `uu`, `up`, and natural-BC Jacobian blocks. Under the true +(dx-first) sympy convention this places `∂F1[gc,dg]/∂L[fc,df]` in +PETSc's `[fc,gc,df,dg]` slot — the major transpose. Every isotropic +rheology (and linear transverse isotropy) has a major-symmetric tangent, +so the transpose was invisible to the whole test suite; the +transverse-isotropic Newton tangent is the first in-repo tangent without +major symmetry, and `-snes_test_jacobian` flagged it at ~1e-1 relative +(vs FD-limited ~1e-4 after the fix). All blocks now use the same +explicit-index loops as `SNES_Vector`/`SNES_MultiComponent`. The +scalar-`p` blocks (`up_*`, `pu_*`) were layout-safe only where a size-1 +axis made the swap a no-op; they are explicit now too. ### `SNES_Vector` migration (2026-04-21) @@ -190,11 +210,12 @@ consumers. When writing a new class that registers a `PetscDSSetJacobian` callback: -1. **Decide on construction style.** For multi-field or novel residual - shapes, prefer the explicit-index pattern from `SNES_MultiComponent` - — it reads like the PETSc documentation and is robust against sympy - convention drift. Reserve `derive_by_array + permutedims` for cases - that match an already-validated solver pattern. +1. **Use the explicit-index pattern from `SNES_MultiComponent` — always.** + It reads like the PETSc documentation and is robust against sympy + convention drift. There is no remaining in-repo + `derive_by_array + permutedims` Jacobian path to copy from, by design: + the last one (Stokes) carried a transposed tangent for months because + every test tangent happened to be major-symmetric (issue #457). 2. **Write a validation test.** For every solver that can be reached at `Nc > 1` with `smoothing > 0`, include a test that: diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index d7db6e05b..3eaf24843 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -3109,9 +3109,16 @@ def _build_c_tensor(self): d = self.dim dv = uw.maths.tensor.idxmap[d][0] - # Use .sym to get sympy expressions from Parameters - eta_0 = self.Parameters.shear_viscosity_0.sym - eta_1 = self.Parameters.shear_viscosity_1.sym + # Bake the WRAPPED Parameter atoms into the tensor, exactly like the + # isotropic ViscousFlowModel path. sympy.diff treats a UWexpression + # atom as a constant, so the default (Picard) tangent stays genuinely + # frozen even when eta_0/eta_1 depend on strain rate; only the Newton + # path (_jacobian_unwrap) exposes that dependence. Baking `.sym` + # (unwrapped contents) here silently un-froze the TI Picard tangent — + # issue #457. The director is `.sym` only for component indexing; + # it carries no velocity dependence. + eta_0 = self.Parameters.shear_viscosity_0 + eta_1 = self.Parameters.shear_viscosity_1 n = self.Parameters.director.sym Delta = eta_0 - eta_1 @@ -3306,6 +3313,17 @@ def __init__(self, unknowns, order=1, integrator: str = "bdf", "Equivalent value of strain rate 2nd invariant (accounting for stress history)", ) + # Persistent container for the yield-limited weak-plane viscosity + # (the ViscoPlasticFlowModel._plastic_eff_viscosity pattern): the + # combined eta_1_eff is stored INSIDE this atom so the Picard tangent + # freezes it, and only _jacobian_unwrap (Newton) sees the strain-rate + # dependence of the yield law. See issue #457. + self._eta1_yield_eff = expression( + R"{\eta_{1,\textrm{eff,p}}}", + 1, + "Yield-limited weak-plane viscosity (effective)", + ) + self._order = order self._yield_mode = "softmin" self._yield_softness = 0.1 @@ -3784,8 +3802,12 @@ def _eta_for_tensor(self, integrator_mode, apply_yield): branch of the hybrid integrator, where the bulk is structurally non-yieldable so clipping is a no-op anyway). """ + # NOTE: bake the WRAPPED Parameter atoms (never `.sym`) — the Picard + # tangent freezes coefficients only while they stay inside UWexpression + # atoms; unwrapped contents expose grad-v to sympy.diff and silently + # un-freeze the default tangent (issue #457). if integrator_mode == "etd": - eta_0 = self.Parameters.shear_viscosity_0.sym + eta_0 = self.Parameters.shear_viscosity_0 eta_1_eff = self.Parameters.shear_viscosity_1 else: # bdf eta_0_raw = self.Parameters.shear_viscosity_0 @@ -3794,7 +3816,7 @@ def _eta_for_tensor(self, integrator_mode, apply_yield): c0 = self._bdf_c0 mu_val = mu.sym if hasattr(mu, 'sym') else mu if mu_val is sympy.oo: - eta_0 = eta_0_raw.sym if hasattr(eta_0_raw, 'sym') else eta_0_raw + eta_0 = eta_0_raw else: eta_0 = eta_0_raw * mu * dt_e / (c0 * eta_0_raw + mu * dt_e) eta_1_eff = self.Parameters.ve_effective_viscosity @@ -3806,7 +3828,11 @@ def _eta_for_tensor(self, integrator_mode, apply_yield): # constants[] atom (so a homotopy can ramp it without a recompile) and # honours yield_smoother. Previously inlined here, which pinned this # model to the sqrt family and to a baked float delta. - eta_1_eff = self._combine_yield(eta_1_eff, vp_eff) + # The combined coefficient is stored INSIDE the persistent container + # so the Picard tangent freezes the yield law (the raw composite + # carries grad-v through the resolved fault-plane shear rate). + self._eta1_yield_eff._sym = self._combine_yield(eta_1_eff, vp_eff) + eta_1_eff = self._eta1_yield_eff return eta_0, eta_1_eff def _assemble_c_tensor(self, eta_0, eta_1_eff): @@ -3896,7 +3922,8 @@ def stress_projection(self): def _build_c_tensor_ve(self): """Build anisotropic tensor with VE η₁ only (no yield).""" d = self.dim - eta_0 = self.Parameters.shear_viscosity_0.sym + # Wrapped atom, not `.sym` — same freezing contract as _eta_for_tensor. + eta_0 = self.Parameters.shear_viscosity_0 eta_1_ve = self.Parameters.ve_effective_viscosity n = self.Parameters.director.sym Delta = eta_0 - eta_1_ve diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index fb3e7f820..4ea0a030c 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -7343,42 +7343,100 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): F1_for_jac = self._jacobian_source(sympy.Array(F1_jac_src)) else: F1_for_jac = self._jacobian_source(F1, self._newton_flux(F1)) + # Normalise to strict (dim, dim) Array indexing for the explicit + # Jacobian loops below. + F1_for_jac = sympy.Array(F1_for_jac).reshape(dim, dim) - G0 = sympy.derive_by_array(F0_jac, self.u.sym) - G1 = sympy.derive_by_array(F0_jac, self.Unknowns.L) - G2 = sympy.derive_by_array(F1_for_jac, self.u.sym) - G3 = sympy.derive_by_array(F1_for_jac, self.Unknowns.L) + # Explicit-index Jacobian construction — writes each entry directly + # into PETSc's flat [fc, gc, df, dg] layout via row-major matrices. + # sympy.derive_by_array is dx-FIRST (derivative indices lead), so the + # previous derive_by_array + permutedims((0,2,1,3)) form assembled the + # MAJOR TRANSPOSE of uu_G3: it placed dF1[gc,dg]/dL[fc,df] in the + # [fc,gc,df,dg] slot. Invisible whenever the tangent has major symmetry + # (frozen C, isotropic eta(edot) Newton, linear transverse isotropy), + # wrong exactly when it does not (transverse-isotropic Newton — issue + # #457). Same explicit-loop construction as SNES_Vector above; the + # layout contract is docs/developer/subsystems/petsc-jacobian-layout.md. + U_list = [self.u.sym[0, c] for c in range(dim)] + L = self.Unknowns.L + Nc = dim - # reorganise indices from sympy to petsc orssdering / reshape to Matrix form - # ijkl -> LJKI (hence 3120) - # ij k -> KJ I (hence 210) - # i jk -> J KI (hence 201) + # Normalise the F0 Jacobian source to a flat list of dim scalar + # entries — F0 arrives as (1, dim) or (dim, 1) depending on how the + # template/bodyforce was written (Array indexing is strict). + _f0_flat = sympy.Array(F0_jac).reshape(dim) + f0_jac_list = [_f0_flat[c] for c in range(dim)] - # The indices need to be interleaved, but for symmetric problems - # there are lots of symmetries. This means we can find it hard to debug - # the required permutation for a non-symmetric problem - permutation = (0,2,1,3) # ? same symmetry as I_ijkl ? # OK - # permutation = (0,2,3,1) # ? same symmetry as I_ijkl ? # OK - # permutation = (3,1,2,0) # ? same symmetry as I_ijkl ? # OK + # uu_G0[fc, gc] = dF0[fc] / dU[gc] + G0 = sympy.zeros(Nc, Nc) + for fc in range(Nc): + for gc in range(Nc): + G0[fc, gc] = sympy.diff(f0_jac_list[fc], U_list[gc]) - self._uu_G0 = sympy.ImmutableMatrix(sympy.permutedims(G0, permutation).reshape(dim,dim)) - self._uu_G1 = sympy.ImmutableMatrix(sympy.permutedims(G1, permutation).reshape(dim,dim*dim)) - self._uu_G2 = sympy.ImmutableMatrix(sympy.permutedims(G2, permutation).reshape(dim*dim,dim)) - self._uu_G3 = sympy.ImmutableMatrix(sympy.permutedims(G3, permutation).reshape(dim*dim,dim*dim)) + # uu_G1[fc*Nc + gc, df] = dF0[fc] / dL[gc, df] + G1 = sympy.zeros(Nc * Nc, dim) + for fc in range(Nc): + for gc in range(Nc): + for df in range(dim): + G1[fc * Nc + gc, df] = sympy.diff(f0_jac_list[fc], L[gc, df]) - fns_jacobian += [self._uu_G0, self._uu_G1, self._uu_G2, self._uu_G3] + # uu_G2[fc*Nc + gc, df] = dF1[fc, df] / dU[gc] + G2 = sympy.zeros(Nc * Nc, dim) + for fc in range(Nc): + for gc in range(Nc): + for df in range(dim): + G2[fc * Nc + gc, df] = sympy.diff(F1_for_jac[fc, df], U_list[gc]) + + # uu_G3[fc*Nc + gc, df*dim + dg] = dF1[fc, df] / dL[gc, dg] + G3 = sympy.zeros(Nc * Nc, dim * dim) + for fc in range(Nc): + for gc in range(Nc): + for df in range(dim): + for dg in range(dim): + G3[fc * Nc + gc, df * dim + dg] = sympy.diff( + F1_for_jac[fc, df], L[gc, dg] + ) - # U/P block (check permutations - hard to validate without a full collection of examples) + self._uu_G0 = sympy.ImmutableMatrix(G0) + self._uu_G1 = sympy.ImmutableMatrix(G1) + self._uu_G2 = sympy.ImmutableMatrix(G2) + self._uu_G3 = sympy.ImmutableMatrix(G3) - G0 = sympy.derive_by_array(F0_jac, self.p.sym) - G1 = sympy.derive_by_array(F0_jac, self._G) - G2 = sympy.derive_by_array(F1_for_jac, self.p.sym) - G3 = sympy.derive_by_array(F1_for_jac, self._G) + fns_jacobian += [self._uu_G0, self._uu_G1, self._uu_G2, self._uu_G3] - self._up_G0 = sympy.ImmutableMatrix(G0.reshape(dim)) # zero in tests - self._up_G1 = sympy.ImmutableMatrix(sympy.permutedims(G1, permutation).reshape(dim,dim)) # zero in stokes tests - self._up_G2 = sympy.ImmutableMatrix(sympy.permutedims(G2, permutation).reshape(dim,dim)) # ? - self._up_G3 = sympy.ImmutableMatrix(sympy.permutedims(G3, permutation).reshape(dim*dim,dim)) # zeros + # U/P block. The constraint field is scalar (Nc_p == 1), so the g-index + # is size 1 and only the derivative indices need explicit placement. + p_scalar = self.p.sym[0] + Gp = self._G # (1, dim) row of dp/dx_dg symbols + + # up_G0[fc, 0] = dF0[fc] / dp + G0 = sympy.zeros(dim, 1) + for fc in range(dim): + G0[fc, 0] = sympy.diff(f0_jac_list[fc], p_scalar) + + # up_G1[fc, dg] = dF0[fc] / d(dp/dx_dg) + G1 = sympy.zeros(dim, dim) + for fc in range(dim): + for dg in range(dim): + G1[fc, dg] = sympy.diff(f0_jac_list[fc], Gp[0, dg]) + + # up_G2[fc, df] = dF1[fc, df] / dp + G2 = sympy.zeros(dim, dim) + for fc in range(dim): + for df in range(dim): + G2[fc, df] = sympy.diff(F1_for_jac[fc, df], p_scalar) + + # up_G3[fc*dim + df, dg] = dF1[fc, df] / d(dp/dx_dg) + G3 = sympy.zeros(dim * dim, dim) + for fc in range(dim): + for df in range(dim): + for dg in range(dim): + G3[fc * dim + df, dg] = sympy.diff(F1_for_jac[fc, df], Gp[0, dg]) + + self._up_G0 = sympy.ImmutableMatrix(G0) # zero in stokes tests + self._up_G1 = sympy.ImmutableMatrix(G1) # zero in stokes tests + self._up_G2 = sympy.ImmutableMatrix(G2) # pressure coupling + self._up_G3 = sympy.ImmutableMatrix(G3) # zero in stokes tests fns_jacobian += [self._up_G0, self._up_G1, self._up_G2, self._up_G3] @@ -7430,24 +7488,48 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): if bc.fn_f is not None: - permutation = (0,2,1,3) # ? same symmetry as I_ijkl ? # OK - bd_F0 = sympy.Array(bc.fn_f) bc.fns["u_f0"] = sympy.ImmutableDenseMatrix(bd_F0) fns_bd_residual += [bc.fns["u_f0"]] - G0 = sympy.derive_by_array(bd_F0, self.Unknowns.u.sym) - G1 = sympy.derive_by_array(bd_F0, self.Unknowns.L) - bc.fns["uu_G0"] = sympy.ImmutableMatrix(sympy.permutedims(G0, permutation).reshape(dim,dim)) # sympy.ImmutableMatrix(sympy.permutedims(G0, permutation).reshape(dim,dim)) - bc.fns["uu_G1"] = sympy.ImmutableMatrix(sympy.permutedims(G1, permutation).reshape(dim,dim*dim)) # sympy.ImmutableMatrix(sympy.permutedims(G1, permutation).reshape(dim,dim*dim)) + # Boundary Jacobians follow the same PETSc [fc, gc, df, dg] + # layout as the bulk blocks — build them with the same + # explicit-index loops (the old permutedims form transposed + # fc/gc here too; harmless only while every natural-BC + # tangent happened to be symmetric). + bd_f0_list = list(sympy.Array(bc.fn_f).reshape(dim)) + + # uu_G0[fc, gc] = d bd_F0[fc] / dU[gc] + G0 = sympy.zeros(dim, dim) + for fc in range(dim): + for gc in range(dim): + G0[fc, gc] = sympy.diff(bd_f0_list[fc], U_list[gc]) + + # uu_G1[fc*dim + gc, df] = d bd_F0[fc] / dL[gc, df] + G1 = sympy.zeros(dim * dim, dim) + for fc in range(dim): + for gc in range(dim): + for df in range(dim): + G1[fc * dim + gc, df] = sympy.diff(bd_f0_list[fc], L[gc, df]) + + bc.fns["uu_G0"] = sympy.ImmutableMatrix(G0) + bc.fns["uu_G1"] = sympy.ImmutableMatrix(G1) fns_bd_jacobian += [bc.fns["uu_G0"], bc.fns["uu_G1"]] - G0 = sympy.derive_by_array(bc.fns["u_f0"], P) - G1 = sympy.derive_by_array(bc.fns["u_f0"], self._G) + # up_G0[fc, 0] = d bd_F0[fc] / dp + G0 = sympy.zeros(dim, 1) + for fc in range(dim): + G0[fc, 0] = sympy.diff(bd_f0_list[fc], p_scalar) + + # up_G1[fc, dg] = d bd_F0[fc] / d(dp/dx_dg) + G1 = sympy.zeros(dim, dim) + for fc in range(dim): + for dg in range(dim): + G1[fc, dg] = sympy.diff(bd_f0_list[fc], Gp[0, dg]) - bc.fns["up_G0"] = sympy.ImmutableMatrix(G0.reshape(dim)) - bc.fns["up_G1"] = sympy.ImmutableMatrix(sympy.permutedims(G1, permutation).reshape(dim,dim)) + bc.fns["up_G0"] = sympy.ImmutableMatrix(G0) + bc.fns["up_G1"] = sympy.ImmutableMatrix(G1) fns_bd_jacobian += [bc.fns["up_G0"], bc.fns["up_G1"]] # Gradient boundary residual (f1_bd) and its Jacobians (g2, g3) @@ -7461,16 +7543,47 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # smooth-kink the Jacobian source so its tangent is Newton- # consistent (same fix as the bulk). Residual u_F1 stays exact. bd_F1_jac = self._jacobian_source(bd_F1) - G2 = sympy.derive_by_array(bd_F1_jac, self.Unknowns.u.sym) - G3 = sympy.derive_by_array(bd_F1_jac, self.Unknowns.L) - bc.fns["uu_G2"] = sympy.ImmutableMatrix(sympy.permutedims(G2, permutation).reshape(dim*dim, dim)) - bc.fns["uu_G3"] = sympy.ImmutableMatrix(sympy.permutedims(G3, permutation).reshape(dim*dim, dim*dim)) + + # uu_G2[fc*dim + gc, df] = d bd_F1[fc, df] / dU[gc] + G2 = sympy.zeros(dim * dim, dim) + for fc in range(dim): + for gc in range(dim): + for df in range(dim): + G2[fc * dim + gc, df] = sympy.diff( + bd_F1_jac[fc, df], U_list[gc] + ) + + # uu_G3[fc*dim + gc, df*dim + dg] = d bd_F1[fc, df] / dL[gc, dg] + G3 = sympy.zeros(dim * dim, dim * dim) + for fc in range(dim): + for gc in range(dim): + for df in range(dim): + for dg in range(dim): + G3[fc * dim + gc, df * dim + dg] = sympy.diff( + bd_F1_jac[fc, df], L[gc, dg] + ) + + bc.fns["uu_G2"] = sympy.ImmutableMatrix(G2) + bc.fns["uu_G3"] = sympy.ImmutableMatrix(G3) fns_bd_jacobian += [bc.fns["uu_G2"], bc.fns["uu_G3"]] - G2 = sympy.derive_by_array(bc.fns["u_F1"], P) - G3 = sympy.derive_by_array(bc.fns["u_F1"], self._G) - bc.fns["up_G2"] = sympy.ImmutableMatrix(G2.reshape(dim, dim)) - bc.fns["up_G3"] = sympy.ImmutableMatrix(G3.reshape(dim, dim*dim)) + # up_G2[fc, df] = d bd_F1[fc, df] / dp + G2 = sympy.zeros(dim, dim) + for fc in range(dim): + for df in range(dim): + G2[fc, df] = sympy.diff(bd_F1[fc, df], p_scalar) + + # up_G3[fc*dim + df, dg] = d bd_F1[fc, df] / d(dp/dx_dg) + G3 = sympy.zeros(dim * dim, dim) + for fc in range(dim): + for df in range(dim): + for dg in range(dim): + G3[fc * dim + df, dg] = sympy.diff( + bd_F1[fc, df], Gp[0, dg] + ) + + bc.fns["up_G2"] = sympy.ImmutableMatrix(G2) + bc.fns["up_G3"] = sympy.ImmutableMatrix(G3) fns_bd_jacobian += [bc.fns["up_G2"], bc.fns["up_G3"]] # Pressure boundary residual and Jacobians (pu, pp blocks) diff --git a/tests/test_1066_stokes_jacobian_layout.py b/tests/test_1066_stokes_jacobian_layout.py new file mode 100644 index 000000000..0f98fc3cb --- /dev/null +++ b/tests/test_1066_stokes_jacobian_layout.py @@ -0,0 +1,283 @@ +"""Stokes velocity-block Jacobian layout and Picard-freezing contracts (issue #457). + +Two independent defects hid behind major-symmetric tangents until the +transverse-isotropic Newton tangent (the first in-repo tangent WITHOUT major +symmetry) exposed them: + +1. ``SNES_Stokes_SaddlePt`` assembled ``uu_G3`` as the MAJOR TRANSPOSE of the + true tangent (``derive_by_array`` is dx-first; the old + ``permutedims((0,2,1,3))`` translation assumed dx-last). Test 1 pins the + assembled ``_uu_G3`` against a finite-difference oracle built from the + RESIDUAL flux — the same comparison ``-snes_test_jacobian`` makes — using a + tangent whose major-asymmetric part is provably non-trivial, so the test + cannot pass by symmetry. + +2. ``TransverseIsotropicFlowModel._build_c_tensor`` baked the UNWRAPPED + ``.sym`` contents of its viscosity Parameters into the c-tensor, exposing + the strain-rate dependence to ``sympy.diff`` and silently un-freezing the + default (Picard) tangent. Tests 2/3 pin the freezing contract: coefficients + hide inside UWexpression atoms under Picard; only the Newton unwrap sees + through them. + +See docs/developer/subsystems/petsc-jacobian-layout.md. +""" + +import numpy as np +import pytest +import sympy + +import underworld3 as uw +from underworld3.function.expressions import UWexpression, unwrap_expression + + +ETA1_FACTOR = 0.01 # weak plane: anisotropy strongly active + + +def _ti_stokes(consistent_jacobian): + """Tiny TI Stokes with a strain-rate-dependent eta_0 (assigned as raw + sympy — the Parameter setter wraps it into the Parameter atom).""" + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) + v = uw.discretisation.MeshVariable("V_jl", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P_jl", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + + x, y = mesh.X + grad = sympy.Matrix( + [[v.sym[0].diff(x), v.sym[0].diff(y)], + [v.sym[1].diff(x), v.sym[1].diff(y)]] + ) + edot = (grad + grad.T) / 2 + eII = sympy.sqrt( + sympy.Rational(1, 2) * (edot[0, 0] ** 2 + edot[1, 1] ** 2) + + edot[0, 1] ** 2 + ) + eta0 = (sympy.Float(0.001) + eII) ** sympy.Rational(-2, 3) + + stokes.constitutive_model = uw.constitutive_models.TransverseIsotropicFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta0 + stokes.constitutive_model.Parameters.shear_viscosity_1 = ETA1_FACTOR * eta0 + rt2 = sympy.sqrt(2) + stokes.constitutive_model.Parameters.director = sympy.Matrix([[1 / rt2, 1 / rt2]]) + + stokes.add_dirichlet_bc((1.0, 0.0), "Top") + stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") + stokes.consistent_jacobian = consistent_jacobian + return stokes, v, p, mesh + + +def _numeric_fns(exprs, stokes, v, p, mesh, Lsyms): + """Fully unwrap each expression, replace gradient/field/coordinate objects + with plain symbols/numbers, lambdify over the 4 velocity-gradient symbols.""" + L = stokes.Unknowns.L + N = mesh.CoordinateSystem.N + flat_L = [Lsyms[k][l] for k in range(2) for l in range(2)] + rep = {L[k, l]: Lsyms[k][l] for k in range(2) for l in range(2)} + rep[v.sym[0]] = sympy.Float(0.1) + rep[v.sym[1]] = sympy.Float(-0.2) + rep[p.sym[0]] = sympy.Float(0.7) + rep[N[0]] = sympy.Float(0.4) + rep[N[1]] = sympy.Float(0.55) + + fns = [] + for e in exprs: + eu = unwrap_expression(e, mode="nondimensional") + eu = sympy.sympify(eu).xreplace(rep) + leftovers = [a for a in eu.free_symbols if a not in flat_L] + assert not leftovers, f"unresolved symbols {leftovers} in {e}" + fns.append(sympy.lambdify(flat_L, eu, "numpy")) + return fns + + +@pytest.mark.level_1 +@pytest.mark.tier_a +def test_stokes_uu_g3_matches_fd_oracle(): + """The assembled Newton uu_G3, read back through PETSc's flat + ((fc*d+gc)*d+df)*d+dg convention, must match a central-difference + derivative of the residual flux. The TI Newton tangent has no major + symmetry, so a transposed assembly fails this at O(1e-1).""" + d = 2 + stokes, v, p, mesh = _ti_stokes(consistent_jacobian=True) + stokes._setup_pointwise_functions() + + Lsyms = [[sympy.Symbol(f"l{k}{l}") for l in range(2)] for k in range(2)] + Lval = np.array([[0.31, -1.17], [0.73, 0.11]]) + flat_val = Lval.flatten() + + # FD ground truth from the residual flux (what -snes_test_jacobian FDs). + res_fns = _numeric_fns(list(stokes._u_F1), stokes, v, p, mesh, Lsyms) + FD = np.zeros((d, d, d, d)) # [fc, gc, df, dg] = dF1[fc, df]/dL[gc, dg] + h = 1.0e-6 + for gc in range(d): + for dg in range(d): + up = flat_val.copy() + dn = flat_val.copy() + up[gc * d + dg] += h + dn[gc * d + dg] -= h + for fc in range(d): + for df in range(d): + idx = fc * d + df + FD[fc, gc, df, dg] = (res_fns[idx](*up) - res_fns[idx](*dn)) / (2 * h) + + scale = np.abs(FD).max() + + # Guard against a true-by-construction pass: the oracle tangent must have + # a non-trivial major-asymmetric part, else this test cannot see a + # transposed assembly at all. + major_asym = max( + abs(FD[a, b, c, e] - FD[b, a, e, c]) + for a in range(d) for b in range(d) for c in range(d) for e in range(d) + ) + assert major_asym / scale > 1.0e-2, ( + "oracle tangent is (near-)major-symmetric — test cannot detect a " + "transposed uu_G3; strengthen the anisotropy in the setup" + ) + + # The solver's assembled uu_G3, read exactly as PETSc reads it. + G3M = stokes._uu_G3 + g3_fns = {} + entries = [] + order = [] + for fc in range(d): + for gc in range(d): + for df in range(d): + for dg in range(d): + entries.append(G3M[fc * d + gc, df * d + dg]) + order.append((fc, gc, df, dg)) + fns = _numeric_fns(entries, stokes, v, p, mesh, Lsyms) + G3 = np.zeros((d, d, d, d)) + for (fc, gc, df, dg), f in zip(order, fns): + G3[fc, gc, df, dg] = f(*flat_val) + + rel = np.abs(G3 - FD).max() / scale + assert rel < 1.0e-5, ( + f"assembled uu_G3 differs from the FD tangent by {rel:.3e} rel — " + "layout transposition or missing tangent terms (issue #457)" + ) + + +@pytest.mark.level_1 +@pytest.mark.tier_a +def test_ti_picard_tangent_is_frozen(): + """Under the default (Picard) tangent the TI viscosity coefficients must + stay hidden inside UWexpression atoms: no raw velocity-gradient + (Derivative) atoms may appear in the assembled uu_G3. The Newton build of + the same problem MUST expose them (positive control that the dependence + exists and the Newton unwrap sees it).""" + d = 2 + + stokes, _, _, _ = _ti_stokes(consistent_jacobian=False) + stokes._setup_pointwise_functions() + # NOTE: velocity gradients are UnderworldFunction symbols, not + # sympy.Derivative atoms — detect dependence the way the tangent + # machinery does: differentiate w.r.t. the L entries. + L = stokes.Unknowns.L + picard_live = [ + (i, k, l) + for i, e in enumerate(stokes._uu_G3) + for k in range(d) + for l in range(d) + if sympy.diff(sympy.sympify(e), L[k, l]) != 0 + ] + assert not picard_live, ( + "Picard uu_G3 depends on the velocity gradient — the TI c-tensor " + "build has un-frozen the default tangent (issue #457): first hits " + f"{picard_live[:4]}" + ) + + # Frozen also means major-symmetric: C has minor+major symmetry, so the + # Picard block must satisfy G3[fc,gc,df,dg] == G3[gc,fc,dg,df] exactly. + G3M = stokes._uu_G3 + for fc in range(d): + for gc in range(d): + for df in range(d): + for dg in range(d): + a = G3M[fc * d + gc, df * d + dg] + b = G3M[gc * d + fc, dg * d + df] + residue = sympy.expand(a - b) + if residue != 0: + # expand() is the cheap canonicaliser; fall back to + # simplify only on the (rare) unproven entry. + residue = sympy.simplify(residue) + assert residue == 0, ( + f"Picard uu_G3 not major-symmetric at " + f"[{fc},{gc},{df},{dg}]" + ) + + stokes_n, _, _, _ = _ti_stokes(consistent_jacobian=True) + stokes_n._setup_pointwise_functions() + Ln = stokes_n.Unknowns.L + newton_live = any( + sympy.diff(sympy.sympify(e), Ln[k, l]) != 0 + for e in stokes_n._uu_G3 + for k in range(d) + for l in range(d) + ) + assert newton_live, ( + "Newton uu_G3 has no velocity-gradient dependence — the unwrap " + "path is not seeing the strain-rate-dependent viscosity, so this " + "test's Picard assertion is vacuous" + ) + + +@pytest.mark.level_1 +@pytest.mark.tier_a +def test_ti_vep_c_tensor_coefficients_are_frozen(): + """TI-VEP with yield active: the c-tensor COEFFICIENTS (including the + yield-limited eta_1_eff, whose raw form carries grad-v through the + resolved fault-plane shear rate) must live inside UWexpression atoms. + E_eff's strain-rate content is residual structure and is exercised + elsewhere; here only the coefficient freezing is pinned.""" + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) + v = uw.discretisation.MeshVariable("V_tv", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P_tv", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + cm = uw.constitutive_models.TransverseIsotropicVEPFlowModel(stokes.Unknowns) + stokes.constitutive_model = cm + cm.Parameters.shear_viscosity_0 = 1.0 + cm.Parameters.shear_viscosity_1 = 0.01 + cm.Parameters.shear_modulus = 100.0 + cm.Parameters.dt_elastic = 0.1 # BDF eta_0 composite is nan at the oo default + cm.Parameters.yield_stress = 0.5 + cm.Parameters.director = sympy.Matrix([0.0, 1.0]) + cm.Parameters.strainrate_inv_II_min = 1.0e-6 + + assert cm.is_viscoplastic + cm._build_c_tensor() + + # Velocity gradients are UnderworldFunction symbols (not sympy.Derivative + # atoms) — detect dependence via d/dL, as the tangent machinery does. + # NOTE: iterating a rank-4 sympy Array yields SUB-ARRAYS — index scalars + # explicitly. + import itertools + + L = stokes.Unknowns.L + d = mesh.dim + c = sympy.Array(cm.c) + live = [ + (idx, k, l) + for idx in itertools.product(range(d), repeat=4) + for k in range(d) + for l in range(d) + if sympy.diff(sympy.sympify(c[idx]), L[k, l]) != 0 + ] + assert not live, ( + "TI-VEP c-tensor depends on the velocity gradient (yield-limited " + "eta_1_eff baked unwrapped) — Picard tangent un-frozen (issue #457): " + f"first hits {live[:4]}" + ) + + # Positive control: the dependence must exist INSIDE the atoms — the + # Newton unwrap of the same tensor must reveal it. + from underworld3.cython.generic_solvers import _jacobian_unwrap + + c_unwrapped = _jacobian_unwrap(c) + revealed = any( + sympy.diff(sympy.sympify(c_unwrapped[idx]), L[k, l]) != 0 + for idx in itertools.product(range(d), repeat=4) + for k in range(d) + for l in range(d) + ) + assert revealed, ( + "Newton unwrap of the TI-VEP c-tensor reveals no strain-rate " + "dependence — the yield law has been lost from the tangent chain" + ) From b0cc24058889a259bc3452f7ceffa307846f5ac6 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 3 Aug 2026 15:32:15 +1000 Subject: [PATCH 2/2] =?UTF-8?q?review:=20response=20commit=20=E2=80=94=20r?= =?UTF-8?q?etire=20the=20last=20permutedims=20(block-constraint=20uu=5FG0)?= =?UTF-8?q?,=20align=20G1=20index=20naming=20to=20dg?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The layout doc's 'no remaining permutedims' claim is now literally true: the block-constrained boundary uu_G0 (content r*(n(x)n), symmetric and shape-protected — numerically unaffected) moves to the same explicit [fc, gc] loops as every other block. Trial-side derivative index in the two G1 loops renamed df -> dg to match the layout doc's convention. Retest: layout + constrained free-slip + constrained SolCx + Stokes suites, 25 passed. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 27 ++++++++++--------- 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 4ea0a030c..32f46d5a6 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -7373,12 +7373,12 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): for gc in range(Nc): G0[fc, gc] = sympy.diff(f0_jac_list[fc], U_list[gc]) - # uu_G1[fc*Nc + gc, df] = dF0[fc] / dL[gc, df] + # uu_G1[fc*Nc + gc, dg] = dF0[fc] / dL[gc, dg] G1 = sympy.zeros(Nc * Nc, dim) for fc in range(Nc): for gc in range(Nc): - for df in range(dim): - G1[fc * Nc + gc, df] = sympy.diff(f0_jac_list[fc], L[gc, df]) + for dg in range(dim): + G1[fc * Nc + gc, dg] = sympy.diff(f0_jac_list[fc], L[gc, dg]) # uu_G2[fc*Nc + gc, df] = dF1[fc, df] / dU[gc] G2 = sympy.zeros(Nc * Nc, dim) @@ -7506,12 +7506,12 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): for gc in range(dim): G0[fc, gc] = sympy.diff(bd_f0_list[fc], U_list[gc]) - # uu_G1[fc*dim + gc, df] = d bd_F0[fc] / dL[gc, df] + # uu_G1[fc*dim + gc, dg] = d bd_F0[fc] / dL[gc, dg] G1 = sympy.zeros(dim * dim, dim) for fc in range(dim): for gc in range(dim): - for df in range(dim): - G1[fc * dim + gc, df] = sympy.diff(bd_f0_list[fc], L[gc, df]) + for dg in range(dim): + G1[fc * dim + gc, dg] = sympy.diff(bd_f0_list[fc], L[gc, dg]) bc.fns["uu_G0"] = sympy.ImmutableMatrix(G0) bc.fns["uu_G1"] = sympy.ImmutableMatrix(G1) @@ -7617,7 +7617,6 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): ## stiffness uu = ∂fn_f/∂u = r·(n⊗n) (0, 0) which conditions the [p,h] ## Schur complement (r=0 ⇒ bare KKT, uu=0). Guarded: no-op for ordinary ## Stokes. - cbc_permutation = (0, 2, 1, 3) for cbc in self._block_constraint_bcs: n_row = cbc.normal # sympy 1×dim Matrix g_sym = cbc.g @@ -7643,11 +7642,15 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): cbc.fns["h_f0"] = sympy.ImmutableDenseMatrix(fn_h) fns_bd_residual += [cbc.fns["h_f0"]] - # uu (0, 0): ∂fn_f/∂u = r·(n⊗n) — AL stiffness (mirror Nitsche shape) - G0 = sympy.derive_by_array(sympy.Array(fn_f), self.Unknowns.u.sym) - cbc.fns["uu_G0"] = sympy.ImmutableMatrix( - sympy.permutedims(G0, cbc_permutation).reshape(dim, dim) - ) + # uu (0, 0): ∂fn_f/∂u = r·(n⊗n) — AL stiffness (mirror Nitsche shape). + # Explicit [fc, gc] placement like every other Jacobian block (the + # content is symmetric, but no permutedims survives on principle — + # see petsc-jacobian-layout.md). + G0 = sympy.zeros(dim, dim) + for fc in range(dim): + for gc in range(dim): + G0[fc, gc] = sympy.diff(fn_f[fc], U_list[gc]) + cbc.fns["uu_G0"] = sympy.ImmutableMatrix(G0) fns_bd_jacobian += [cbc.fns["uu_G0"]] # uh (0, h): ∂fn_f/∂h = n — mirror the up_G0 (velocity,scalar) shape