Fix GenericFluxModel flux storage — preserve symbolic field dependencies

[bknight1]

Summary

GenericFluxModel.Parameters.flux setter (and __init__) wrapped each flux component in a UWexpression via validate_parameters(), creating intermediate symbols (q_0, q_1) that the JIT treated as independent constants. This lost the dependency on the unknown field variable and produced a singular matrix.

Fix

Store raw sympy expressions directly in the flux setter, bypassing the UWexpression wrapper. The __init__ now also sets a raw sympy.zeros vector instead of validated components.

Tests

New test file tests/test_1060_generic_flux_model_update.py with three tier_a, level_1 tests:

• test_flux_update_via_force_setup: Changing flux via _force_setup=True works (solution scales correctly with k)
• test_flux_update_via_solver_link: Changing flux via _solver back-link triggers automatic rewire (no _force_setup needed)
• test_flux_stores_raw_sympy: Verifies flux stores raw sympy, not UWexpression-wrapped symbols (no q_0, q_1 artifacts)

Verification

All three tests pass with the fix; before the fix they would produce singular matrices.

[Copilot]

Pull request overview

This PR fixes GenericFluxModel.Parameters.flux storage so flux components are preserved as raw SymPy expressions (instead of being wrapped into UWexpression placeholders), preventing loss of symbolic dependencies on the unknown field during JIT assembly (which previously could yield singular systems).

Changes:

• Updated GenericFluxModel._Parameters.__init__ and flux setter to store raw SymPy vectors directly (no validate_parameters() wrapping into q_i symbols).
• Added a new tier_a / level_1 regression test module covering flux updates across solves and ensuring flux is stored as raw SymPy.

Reviewed changes

Copilot reviewed 2 out of 2 changed files in this pull request and generated 1 comment.

File	Description
src/underworld3/constitutive_models.py	Stores GenericFluxModel flux as raw SymPy to preserve field-variable dependencies for JIT assembly.
tests/test_1060_generic_flux_model_update.py	New regression tests validating flux update propagation and raw-SymPy storage.

Comments suppressed due to low confidence (2)

tests/test_1060_generic_flux_model_update.py:86

• This test asserts internal/private solver state (_solver_is_setup, _needs_function_rewire). Those flags are implementation details and may change without breaking the external behaviour this test is trying to validate (that changing Parameters.flux triggers a rewire and the next solve() uses the new expression). It’s safer to validate behaviour via the second solve + ratio checks only.

    assert not solver.constitutive_model._solver_is_setup
    assert solver._needs_function_rewire

tests/test_1060_generic_flux_model_update.py:116

• The raw-sympy check relies on "q_" not in str(flux), which is brittle (string representations can change, and legitimate symbols could contain q_). Prefer a structural check that elements are not UWexpression instances.

    # Check it's a plain sympy Matrix, not a UWexpression
    assert isinstance(flux, sympy.Matrix), f"Expected sympy.Matrix, got {type(flux)}"
    # The elements should be raw sympy, not UWexpressions
    # UWexpressions show as "q_0 = ..." in str representation
    flux_str = str(flux)
    assert "q_" not in flux_str, (

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[lmoresi]

Adversarial review — live probes on the PR head and its exact parent commit

Verdict: the fix is correct and merge-worthy; the stated mechanism is wrong in a way that matters for the record. No change to the fix itself is needed — a wording correction and one test cleanup follow in a response commit, then this merges.

The true mechanism (probed, decisive)

Reproduced pre-fix on the PR's exact parent (967fb30b): Poisson + GenericFluxModel, flux = grad(u) → DIVERGED_LINEAR_SOLVE after 0 iterations, solution identically zero. Diagnostics:

Evidence	Value (pre-fix)
constants_manifest	empty — q_0/q_1 never routed to constants[]
_is_truly_constant(q_0)	False (it contains the unknown's gradient)
_G3 (dF1/d∇u)	zero matrix
Same pre-fix code, consistent_jacobian=True	G3 = Identity, solves correctly (peak 0.07365, the textbook value)

So this is not "the JIT treated the components as independent constants" — the manifest was empty and the residual was never corrupted (PR #416's reveal machinery handled the wrappers correctly, which is why the Newton-mode run converges with the identical residual). The defect is in the default Picard tangent: _jacobian_source deliberately differentiates the RAW wrapped flux (frozen-coefficient semantics; unwrap-before-differentiate runs only on the Newton path). For every other model that's safe — the wrapper holds only a coefficient, and d(κ·∇u)/d∇u = κ ≠ 0. GenericFluxModel wrapped the entire flux, so the opaque q_i swallowed the ∇u structure and "freeze the coefficient" froze everything: G0–G3 = 0 → structurally singular operator.

This distinction matters: the current docstring/comment wording implicitly indicts the #416 constants machinery, which is exonerated, and would misdirect the next person who hits a singular matrix.

Findings

1. Fix works — PROBED. All 3 new tests pass; independent probe: GenericFluxModel q = 3∇u vs DiffusionModel κ = 3 agree to 0.0 relative difference; flux updates between solves work via both _force_setup and the _solver-link reset.
2. Mechanism wording — PROBED, correction required (above). Applied in the response commit.
3. Units survive raw storage — PROBED. Dimensional κ as uw.expression under active reference quantities: bit-identical to DiffusionModel, and κ still lands in constants_manifest (coefficient-level rampability retained). Nothing in the units chain depended on the whole-flux wrapper.
4. Behavioral note — PROBED. With raw storage, GenericFluxModel's default tangent is now effectively full Newton: probe q = (1+u²)∇u converges in 2 SNES iterations under the Picard default (there is no wrapper left to freeze, so the Picard/Newton switch is a no-op for this model). Usually a win; documented in the response commit.
5. Collateral clean — PROBED. Diff touches only GenericFluxModel._Parameters; validate_parameters and all other models untouched; test_0103 (rampable constants), test_0200, test_0812 (units) + the new file: 14/14 on the PR build.
6. Test cleanups (minor). Parameters._solver = solver in the test setup is redundant (the constitutive_model setter wires it) and masks a regression in the automatic wiring — deleted in the response commit. Tests hard-code MUMPS; house preference is plain lu (all review probes used it) — left as-is to keep the response commit minimal, noted here.

The trap that remains (scoped precisely)

The general hazard is narrower than "wrapped field-dependent parameters break Jacobians": coefficient-level wrapping under the Picard default is intentional frozen-coefficient semantics and stays nonsingular. The dangerous pattern is exactly one — a wrapped Parameter that swallows the entire ∇u dependence of a flux term — and GenericFluxModel was the only in-tree instance. Future models (or user-built constitutive models) can re-create it; a warn-if-wrapped-sym-contains-the-unknown's-gradient guard is conceivable but out of scope here.

Adversarial review per project policy — Underworld development team with AI support from Claude Code