Stress histories in a units-aware model: reduce on the way in, read back in Pa (#788)

[lmoresi]

Stacked on #785. Makes the viscoelastic stress histories, and swarm advection, work in a model with reference quantities.

Design. A history's stores hold non-dimensional values in .data, which is what the solver reads and what copies between stores go through. A solver builds its stress history with the stress's units, so .array and evaluate of the store read in Pa, as for any MeshVariable. What enters is reduced on the way in: the timestep where each flavour receives it, and anything written from evaluate (one rule, _write_evaluated: a result that carries units is reduced, one that does not is already non-dimensional).

What failed before, Maxwell shear box with reference quantities (1 km, 1e21 Pa s, 1 Myr):

	before
semi_lagrangian	TypeError on np.copy of a unit-aware array
integration_point	TypeError: unhashable UWQuantity (timestep in the trace cache key)
lagrangian	1.5e148 after 20 steps; and the particles did not move at all (swarm advection used the dimensional velocity)
forward, eulerian	correct non-dimensional stress, nothing read back in Pa
order 2, any flavour	silently BDF-1 whenever the timestep unit differs from the reference time unit
a history started from a moving flow	initial level written dimensional (Eulerian, Lagrangian_Swarm)

Test test_1064_stress_history_units: every flavour, in a units model with the timestep in kyr (not the reference unit) and as the same problem in plain numbers, from a moving start, at orders 1 and 2. The stores agree to 1e-6, and the stress reads back in Pa as the non-dimensional value times the stress scale.

The stress-history, swarm and advection test files pass (324, plus one known fragile test that now passes); the style gate is clean.

History of this PR. The first version stamped the flux's units onto every store; the second made the stores unitless with a separate dimensional read (carried). Adversarial review found the second reimplemented MeshVariable units and that its test was blind to the timestep, the order-2 guard, the initial writes and the swarm motion; this version replaces both.

Fixes #788.

Underworld development team with AI support from Claude Code

🤖 Generated with Claude Code

https://claude.ai/code/session_017kSkAq7oJ5J3XisuvLBovo

[lmoresi]

Adversarial review before the push, one pass over the diff with probes run in the worktree. Found and fixed in the pushed commit:

• The stamp read the flux's units with get_units, and the flux contains the store's own symbol, so the store's units fed back: pascal * second ** 2 / megayear ** 2 on the first stamp, one more s/Myr factor on every re-assignment (a parameter change re-assigns psi_fn), and .array read 2.7e34. The stamp now uses reduced units, a fixed point under re-assignment; the test re-assigns the flux and checks the units and the value again.
• The integration-point flavour's projection scratch was in its store list and read back without reduction, so one re-assignment put dimensional values into storage (8.7e14 Pa after one step). It is scratch and is no longer stamped.
• The per-flavour store lists had been inserted between a timing decorator and __init__, timing the getter instead of the constructor.
• _units_of_components raised on a scalar flux; it coerces now, and refuses components of different dimensionality rather than taking the first.
• The particle flavour computed its own units at construction, before its slots existed, and nothing stamped after; every flavour that assigns psi_fn before building its slots now stamps once the slots exist. Symbolic lists no stores.
• The test now has a .array baseline in Pa and drops a tautological units check.

Recorded, not fixed: the strain-rate forcing stores keep units=None; get_units returns None for some fluxes that mix a unitless symbol in; _dt keeps the caller's form in the trace-back flavours. Checked and fine: the component layout against the old .array[:, i, j] semantics for symmetric and full tensors, the timestep reduction on floats, quantities and expressions, snapshot restore by name, the JIT reading non-dimensional storage regardless of the stamp.

[lmoresi]

Why this is red, and what it needs

The two CI failures are not this PR's. They are byte-identical across #785, #789, #795 and #800:

tests/test_1060_nitsche_freeslip.py::test_nitsche_normal_velocity_zero
    assert np.float64(0.0001008010278702676) < 0.0001
tests/test_1070_free_surface_plume.py::test_freesurface_strong_constraint_beats_penalty
    AssertionError: strong constraint not better: 1.73e-02 vs penalty 3.25e-02

Both tests were already fixed on development, for exactly these numbers, and neither name exists there any more. development's CI is green at 1391e4aa. The replacements name the failures in their own docstrings:

• test_1060 → test_nitsche_constrains_the_wall_normal_velocity, a relative bound of 0.1: "the previous version of this test asserted an ABSOLUTE 1e-4, which scaled with the buoyancy forcing rather than with anything about Nitsche, and sat a fraction of a percent from failing for unrelated reasons."
• test_1070 → test_freesurface_strong_constraint_tracks_the_prescribed_rate, an absolute bound of 5e-2 against the datum: "It replaces a comparison against the penalty path (errors["strong"] < 0.5 * errors["penalty"]) … measured 2026-09-14: 1.06e-2 on macOS, 1.73e-2 on the Linux CI runner, where the annulus triangulates differently."

So this branch is stale, not broken — 58 commits behind development.

The merge is not clean, and the conflicts are the same on all four

git merge origin/development conflicts in exactly two files on every one of the four branches (#800 adds a third, tests/test_0066_integration_point_slcn.py):

scripts/test.sh — one hunk, take THIS branch's side. development narrowed the glob to tests/test_1100*py; this side has tests/test_110*py. test_1101_advdiff_swarm_rotating_gaussian.py matches none of development's globs (test_1100*, test_1110*, test_1120*), so the narrow side would leave a test file dark — #721 recurring. scripts/check_test_coverage.py verifies the resolution.

src/underworld3/systems/ddt.py — five hunks, and four of them are mechanical. Four have a zero-line development side: they are this branch's own additions (applies_inflow_value, commits_flux_in_post_solve, the inflow_value setter — the #745/#783 inflow work), flagged only because surrounding context moved. None of those three names exists on development. Take this branch's side.

The fifth, at the top of _DDtBase, is a real design conflict and needs the author: this branch adds update_exp_coefficients / _exp_alpha / _exp_phi to the base class; development has since put update_exp_coefficients in three per-flavour classes (ddt.py:1299, 1734, 3612) with _exp_alpha/_exp_phi at :3623, 3628, and adds _note_history_shift in that same place. Keeping both gives a base-class definition under three overrides. development's per-flavour layout is the later design, but which one this branch's ETD work is written against is the author's call, not a merge-tool decision — and this is the DDt core, where CLAUDE.md asks for benchmarking rather than a judgement.

Order matters

All four touch ddt.py in the same region, so they conflict with each other as well as with development. They want merging one at a time, smallest first (#789, #785, #795, then #800), with a rebuild and the transport tests between each — the base moves for the rest after every merge.

Found in the backlog sweep (#819). No commits pushed to this branch.