SNES_Scalar still assembles its tangent with derive_by_array: the #457 transpose survives in the one solver #493 did not convert

[lmoresi]

#493 fixed #457 by replacing derive_by_array + permutedims with explicit-index
loops. It converted SNES_Vector and SNES_Stokes_SaddlePt. SNES_Scalar was
not converted and still carries the same construction, and the same transpose.

It is latent, not active — nothing reachable through the public API hits it
today. But it is latent for the same reason #457 was invisible for months, which
is why it is worth closing rather than leaving.

The transpose, demonstrated

petsc_generic_snes_solvers.pyx ~3732 assigns the result straight to _G3
with no reordering:

G3 = sympy.derive_by_array(F1_jac, L)
...
self._G3 = sympy.ImmutableMatrix(G3)

derive_by_array is dx-first, so for a non-symmetric flux tensor:

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]]

Same root cause as #457: derive_by_array leads with the derivative indices, so
the assembled tangent is the transpose of the one PETSc's [fc, gc, df, dg]
layout expects.

Why nothing fails today

The only anisotropic scalar model, AnisotropicDiffusionModel, builds a
diagonal tensor — sympy.diag(*validated) from a (dim,1)/(1,dim)
vector. A diagonal matrix is its own transpose, so the index order cannot show.
Verified on the supported path: kappa = [1, 7] gives
_G3 = [[1, 0], [0, 7]], correct entry by entry.

That is precisely how #457 hid in Stokes — invisible under frozen C, isotropic
eta(edot) Newton and linear transverse isotropy, wrong exactly when the
tangent lost its major symmetry.

A second failure mode in the same idiom

Feeding a full (non-diagonal) tensor to the scalar DiffusionModel — misuse of
that model, but it is what a user reaches for first — does not raise a useful
error. The diffusivity is carried as one opaque symbol kappa, and
derive_by_array then returns an array that lies about its shape:

F1.shape (2,)   backing len 2
G3.shape (2, 2) backing len 2      <- should be 4

.shape and contents disagree, nothing checks, and it surfaces three frames
later as a bare IndexError: list index out of range inside
sympy.ImmutableMatrix. The traceback names sympy internals and not the
diffusivity, so the cause is not recoverable from the error.

So the one idiom carries three hazards: dx-first index order, shape/content
divergence with no check, and silent success whenever the tangent happens to be
symmetric.

Suggested fix

Convert SNES_Scalar to the explicit-index loops already used by the other two
solvers, finishing #493. The layout contract is
docs/developer/subsystems/petsc-jacobian-layout.md.

The regression test needs a full tensor to exist first. A diagonal one
cannot catch a transpose. That means either AnisotropicDiffusionModel accepts
a full (dim, dim) conductivity, or the test builds F1 directly. The former
is worth having on its own account — a foliation- or fault-aligned conductivity
is a real modelling case and is the direction this project is heading — and it
would make the defect reachable, so it should land with the fix rather than
before it.

Found while checking that #457 was genuinely closed (it is). Related: #457, PR
#493, and Charter §8 on baselines tight enough to fail when a default changes.

Underworld development team with AI support from Claude Code

[lmoresi]

Correction: on the Darcy path this is live, not latent. SNES_Darcy derives
from SNES_Scalar, and dispersivity in Darcy models is exactly the use that
reaches it.

Measured on development (347d7e6), SteadyStateDarcy + DarcyFlowModel:

Parameters.permeability	result
scalar	setup OK, g3 index order correct
diagonal matrix	IndexError: list index out of range
symmetric full matrix	IndexError
non-symmetric full matrix	IndexError

The descriptor accepts the matrix without validation — permeability is
declared as a scalar api_tools.Parameter — and the failure lands later in
_setup_pointwise_functions, as a bare IndexError from sympy.ImmutableMatrix
whose traceback names only sympy internals. Nothing points at the permeability.

What does work

AnisotropicDiffusionModel on a Darcy solver, given a (1, dim) vector:

darcy.constitutive_model = uw.constitutive_models.AnisotropicDiffusionModel
darcy.constitutive_model.Parameters.diffusivity = sympy.Matrix([[2.0, 5.0]])
# -> _G3 = [[2, 0], [0, 5]], index order correct

That builds sympy.diag(...) internally, so it is diagonal only.

Why that matters for dispersivity specifically

The hydrodynamic dispersion tensor

D = alpha_T |v| I + (alpha_L - alpha_T) v (x) v / |v|

is symmetric but not diagonal — the off-diagonal terms vanish only when the
flow is axis-aligned. So the diagonal route cannot represent it for general
flow, and the full-tensor route crashes.

Two consequences, and they are different problems:

1. The crash is the live defect. Any attempt to supply a real dispersion
   tensor fails with an error that does not name its cause.
2. The transpose is NOT the risk here. A dispersion tensor is symmetric, so
   even once the crash is fixed the derive_by_array index order would be
   invisible for this use — which is precisely the trap: it would work, and go
   on working, until someone supplies a non-symmetric conductivity.

So the fix needs both halves, and the ordering in the issue above still holds:
convert SNES_Scalar to explicit-index loops, and let the model accept a full
(dim, dim) tensor so the regression test can be written with a non-symmetric
one. Validating the descriptor (reject a matrix, or accept and handle it) would
turn today's IndexError into a usable message in the meantime.

Underworld development team with AI support from Claude Code

[lmoresi]

Correction on the test requirement, and on what is actually broken.

Earlier comments said the regression test needs a full non-symmetric tensor, and
that closing this waits on an API decision. Both are wrong.

F1 is a read-only property on every concrete scalar solver and F0/F1 on
the base class raise RuntimeError, so there is no route by which a user
supplies an arbitrary flux. Every path to dF1/dL goes through a constitutive
model, and all of them are symmetric — scalar kappa*I, AnisotropicDiffusionModel's
diag(...), and a hydrodynamic dispersion tensor, which is symmetric by
construction. The transpose is therefore unobservable through the public
API, and a test that fabricated a non-symmetric constitutive state would be
guarding a state that cannot occur (Charter §5).

So this issue splits cleanly in two, and only one half is a live defect.

Live, and what a user meets. A matrix supplied to
DarcyFlowModel.Parameters.permeability is accepted without validation — it is
declared a scalar api_tools.Parameter — and then fails inside
_setup_pointwise_functions with a bare IndexError from sympy internals whose
traceback names nothing useful. Diagonal, symmetric-full and non-symmetric-full
all crash; only a scalar works. That is testable today and needs no new API:
the contract is that the parameter either builds the right flux or refuses with
a message naming the permeability.

Latent, and correctness-by-construction. The derive_by_array index order.
PR #748 converts SNES_Scalar to explicit-index loops, finishing #493. It fixes
no answer anybody can currently get wrong, and it should not be described as
though it does. It earns its place because 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 is there for the first such
flux instead of waiting to be rediscovered by it.

Still open here after #748: validating the permeability parameter, and — if
wanted on its own merits rather than as a test fixture — a route to a full
(dim, dim) tensor for dispersivity, which the diagonal AnisotropicDiffusionModel
cannot express for non-axis-aligned flow.

Underworld development team with AI support from Claude Code