Repository navigation
Viscoelastic stress history in a units model: the particle store takes dimensional magnitudes, the nodal flavour crashes, and no history store carries stress units #788
Description
Activity
Labelled
fixed-in-PR. The fix is in #789, which is not merged: it sits two deep in the stack rooted onfeature/forward-parallel, which has no PR.#789 is green and mergeable against its base after today's refresh, and its
developmentmerge was validated locally atlevel_1 tier_a— 1397 passed, 0 failed. The two CI failures it carried for ten days were testsdevelopmenthad already renamed, not its own.The label means the fix is written, not released.
- added a commit that references this issue
on Oct 8, 2026 Fixed by #789, merged to
development(3a9069d9).The stress history stores are now built with their units at construction,
.datastays non-dimensional, and there is one write rule — so the particle store no longer takes dimensional magnitudes and the nodal flavour no longer crashes.Two things from this work worth keeping beside the fix:
- swarm advection was unit-blind, which is the same class of defect one layer down;
v.coordsis metres in a units model, which is the trap that makes this easy to reintroduce.
The template for testing it is
test_1064: a units run against a plain run,dtin kyr, a moving start, order 2.Closing manually:
Closeslines here fire whendevelopmentreachesmain. Removing thefixed-in-PRlabel.
Measured on the Maxwell shear box with reference quantities (length 1 km, viscosity 1e21 Pa s, time 1 Myr), `ViscoElasticPlasticFlowModel` order 1, dt 0.1 Myr.
Particle flavours write dimensional magnitudes into non-dimensional storage. `uw.function.evaluate(flux, particle_coords)` returns a `UnitAwareArray` whose magnitude is in Pa. `Lagrangian.update_post_solve`, `Lagrangian.initialise_history` and their `Lagrangian_Swarm` counterparts write that magnitude straight into `psi_star[k].data`, the non-dimensional work array the proxy and the constitutive model read. Each step multiplies by the stress scale (3.2e7 here):
| step | stored max |psi*| | velocity storage max |
|---|---|---|
| 0 | 2.9e6 | 0.5 |
| 1 | 8.3e13 | 0.62 |
| 2 | 2.4e21 | 1.6e6 |
The nodal and integration-point flavours reduce through `_to_nondim_ndarray` at the same point, so the storage is right there. Non-dimensional storage is the design; the missing piece is the reduction on write.
The nodal flavour crashes in the same model: `commit_flux_to_history` does `np.copy(self.psi_star[0].array[...])` on a `UnitAwareArray` (`TypeError: no implementation found for 'numpy.copy'`).
No history store carries stress units. Every flavour resolves `psi_units` from `psi_fn` at construction, and the solver constructs the history with `sympy.Matrix.zeros` as the placeholder `psi_fn`, so `_psi_units` is "dimensionless" although the flux is in Pa. The particle stores are created with no `units` at all, so `psi_star[k].array` reads back non-dimensional numbers with no units. Reading `.array` in Pa needs the units attached when the solver assigns the real flux, for every flavour.
There is no units-model test of any stress history; the shear-box tests all run without reference quantities.
Found while reviewing #785 (marked `TODO(BUG)` there at the particle write).