SNES_Scalar: explicit-index tangent loops, finishing the #457 conversion (#747)

[lmoresi]

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
g3[i][j] should be dF1[i]/dL[j]         = [[1, 7], [3, 2]]

derive_by_array is dx-first, so it yields [dg][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 κ and
AnisotropicDiffusionModel builds a diagonal one — and wrong the moment it is
not.

Not academic. SNES_Darcy derives from SNES_Scalar, and a hydrodynamic
dispersion tensor

D = α_T |v| I + (α_L − α_T) v⊗v / |v|

is exactly the case that reaches it. It is symmetric but not diagonal unless
the flow is axis-aligned.

Scope

Only G3 changes. G1 and G2 flatten to unambiguous sequences (one free
index each), and 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 it, and it surfaces frames later as a bare IndexError from sympy
internals naming nothing useful.

Verified as a no-op where the code was already right

case	shapes	g3[i][j] == dF1[i]/dL[j]
scalar κ	(1,1) (2,1) (1,2) (2,2) — unchanged	CORRECT
diagonal κ	(1,1) (2,1) (1,2) (2,2) — unchanged	CORRECT

35 passed, 1 xpassed across Poisson, Poisson-natural-BC, Poisson-spherical,
Poisson-constants, Darcy, Darcy sign regression, Transient Darcy, Richards,
AdvDiffusion Cartesian, AdvDiffusion annulus, and Poisson-with-units. The xpass
is a pre-existing fragile AdvDiff test whose own marker states either outcome is
acceptable.

This 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 — and a
symmetric tensor would pass whether or not the bug were present, so the test
must use a non-symmetric one.

Closing #747 needs an API decision on where a full tensor belongs:
DarcyFlowModel.permeability (smallest change, but permeability and
dispersivity are different physics sharing a slot), AnisotropicDiffusionModel
widened from diagonal to full (honest to its name, already works on the Darcy
solver), or a new model. Left for the maintainers.

Underworld development team with AI support from Claude Code

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[lmoresi]

Correction to this PR's framing, after maintainer review. Two claims in the
description overreach and should not stand as written.

The transpose is not observable through the public API

F1 is a read-only property on every concrete scalar solver — derived from the
constitutive model or the DDt flux — and F0/F1 on the base class raise
RuntimeError. There is no setter, so a user cannot hand the solver an
arbitrary flux. Every route to dF1/dL therefore runs through a constitutive
model, and every shipped one is symmetric:

model	dF1/dL
DiffusionModel (scalar kappa)	kappa * I — symmetric
AnisotropicDiffusionModel	diag(kappa_0 ...) — symmetric
hydrodynamic dispersion a_T|v|I + (a_L-a_T) v(x)v/|v|	symmetric by construction

So this PR is not fixing a wrong answer that anyone can currently produce.
It makes the assembly correct by construction rather than by accident, and it
removes the derive_by_array shape-lying hazard. Both are worth having; neither
is "the Darcy tangent is wrong today", and the description implied otherwise.

The regression test I proposed would have tested an unreachable state

Asserting the index order needs a non-symmetric dF1/dL, and nothing in the API
can produce one. Writing a test that fabricates that state through a
constitutive model would be guarding a state that cannot occur — Charter §5.
Dropping that requirement, and with it the claim that this PR waits on a
tensor-API decision. It does not.

What IS worth testing is reachable today

darcy.constitutive_model.Parameters.permeability = sympy.Matrix([[2, 3], [3, 5]])
# accepted unvalidated -> bare IndexError three frames down in sympy internals

permeability is declared a scalar api_tools.Parameter; a matrix is accepted
and fails later with a traceback naming nothing useful. That is the failure a
user actually meets. The test is that the API either builds the right flux or
refuses with a message naming the permeability.

Why the explicit loops still earn their place

The flux law is general. A non-symmetric dF1/dL is legitimate outside
diffusion even though no diffusive tensor is asymmetric — so the correct
construction has forward value for the first such flux, rather than waiting to
be rediscovered by it.

Underworld development team with AI support from Claude Code