From de71b7f7d7abd862d7fd69cc178363aee5ceeff9 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 16 Sep 2026 10:33:04 +1000 Subject: [PATCH] SNES_Scalar: explicit-index tangent loops, finishing the #457 conversion (#747) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit PR #493 fixed #457 by replacing `derive_by_array` + `permutedims` with explicit-index loops, and converted `SNES_Vector` and `SNES_Stokes_SaddlePt`. `SNES_Scalar` was missed, and still carried the same transpose: F1 = [L0 + 7 L1, 3 L0 + 2 L1] k = [[1,7],[3,2]] ImmutableMatrix(derive_by_array(F1, L)) = [[1,3],[7,2]] <- k TRANSPOSED `derive_by_array` is dx-FIRST, so it yields [dg][df] = dF1[dg]/dL[df] where PETSc's g3 wants [df][dg]. Invisible whenever the flux tensor is symmetric — which covers every shipped scalar model, since `DiffusionModel` takes a scalar kappa and `AnisotropicDiffusionModel` builds a diagonal one — and wrong the moment it is not. **This is not academic.** `SNES_Darcy` derives from `SNES_Scalar`, and a hydrodynamic dispersion tensor `D = a_T|v| I + (a_L - a_T) v(x)v / |v|` is exactly the case that reaches it. Only G3 changes. G1 and G2 flatten to unambiguous sequences (one free index each), and the shapes are held identical to the previous form — G0 (1,1), G1 (cdim,1), G2 (1,cdim), G3 (cdim,cdim) — so the JIT sees the same flat layout and only the index order moves. The loops also avoid the second hazard in that idiom. Handed a flux built from a matrix-valued symbol, `derive_by_array` returns an array whose `.shape` disagrees with its backing store (claims (2,2), holds 2 entries); nothing checks, and it surfaces frames later as a bare `IndexError` from sympy internals naming nothing useful. Verified a no-op where the code was already right: scalar and diagonal kappa give identical shapes and `g3[i][j] == dF1[i]/dL[j]` in both. 35 passed, 1 xpassed (a pre-existing fragile AdvDiff test whose marker states either outcome is acceptable) across Poisson, Darcy, Transient Darcy, Richards, AdvDiffusion and the units Poisson. Does NOT close #747: a full (dim, dim) tensor is still unreachable through any shipped model, so the regression test that would pin this cannot be written yet. That needs an API decision — see the issue. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 44 ++++++++++++++++--- 1 file changed, 38 insertions(+), 6 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 8cbe26c5..022ee14b 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -3786,12 +3786,44 @@ class SNES_Scalar(SolverBaseClass): f0_jac = self._jacobian_source(f0) F1_jac = self._jacobian_source(F1, self._newton_flux(F1)) - G0 = sympy.derive_by_array(f0_jac, U) - G1 = sympy.derive_by_array(f0_jac, L) - G2 = sympy.derive_by_array(F1_jac, U) - G3 = sympy.derive_by_array(F1_jac, L) - - # Re-organise if needed / make hashable + # Explicit-index Jacobian construction. `sympy.derive_by_array` is + # dx-FIRST (the derivative indices lead), so `derive_by_array(F1, L)` + # yields [dg][df] = dF1[dg]/dL[df] — the TRANSPOSE of the g3 PETSc + # wants. For a scalar unknown that is invisible whenever the flux + # tensor is symmetric (every shipped model: scalar kappa, and + # AnisotropicDiffusionModel's diagonal kappa), and wrong as soon as it + # is not. Same defect as issue #457 in the Stokes/vector tangents, + # fixed there by PR #493; this is the solver that conversion missed + # (#747). Layout contract: + # docs/developer/subsystems/petsc-jacobian-layout.md + # + # derive_by_array carries a second hazard these loops avoid: handed a + # flux built from a matrix-valued symbol it returns an array whose + # .shape disagrees with its backing store, which surfaces later as a + # bare IndexError from sympy internals naming nothing useful. + # + # Shapes are unchanged from the previous form, so the JIT sees the same + # flat sequences: G0 (1,1), G1 (cdim,1), G2 (1,cdim), G3 (cdim,cdim). + # Only G3's index order moves. + Uc = U[0] + G0 = sympy.zeros(1, 1) + G0[0, 0] = sympy.diff(f0_jac[0], Uc) + + # G1[df, 0] = d f0 / d L[df] + G1 = sympy.zeros(cdim, 1) + for df in range(cdim): + G1[df, 0] = sympy.diff(f0_jac[0], L[df]) + + # G2[0, df] = d F1[df] / d U + G2 = sympy.zeros(1, cdim) + for df in range(cdim): + G2[0, df] = sympy.diff(F1_jac[df], Uc) + + # G3[df, dg] = d F1[df] / d L[dg] + G3 = sympy.zeros(cdim, cdim) + for df in range(cdim): + for dg in range(cdim): + G3[df, dg] = sympy.diff(F1_jac[df], L[dg]) self._G0 = sympy.ImmutableMatrix(G0) self._G1 = sympy.ImmutableMatrix(G1)