The timestep as a transaction: model.step, a run transcript you can watch and publish, and eight silent defects

[lmoresi]

What this is

A run had no step boundary. Nothing in the library could say what a timestep
did, whether it finished, or what state it started from — so nothing could
check it, replay it, or walk it backwards.

This adds that boundary, a log of it you can watch a run through, a figure of
it you can publish, and fixes eight defects found by building it and then by
using it. Each was silent: a wrong answer or a change that quietly did
nothing, never an error.

The pattern

model = uw.get_default_model()
model.set_reference_quantities(...)      # before the mesh; units are a choice
                                         # you can only make here
mesh = uw.meshing.UnstructuredSimplexBox(...)

model.tracker.time = uw.quantity(0.0, "Myr")
model.record_every = 1

while model.tracker.time < end_time:
    dt = adv_diff.estimate_dt()
    with model.step(dt):
        adv_diff.solve(timestep=dt)
        stokes.solve(zero_init_guess=False)

<step 0 'convect' dt=0.01 solve:SNES_AdvectionDiffusion(T) -> history_shift:EulerianSUPG(T) -> solve:SNES_Stokes(V)>

model.step(dt) owns a time interval. The clock reads the END of that interval
for the whole block, because an implicit residual is centred at t+dt and a
driven boundary must be evaluated there — committing only on exit would put
every implicit coefficient one step late, which is a first-order error that
looks right and converges. The advance commits on clean exit and only then, so
a step that raises or is rejected on a Courant check leaves the clock alone.

Everything the block did is recorded in model.transcript, named by what it
solves. That record answers what a run actually did without the script being
instrumented, which is the question you want to ask of someone else's model.

model.record_every asks each step to keep the state it started from, and
model.rewind() undoes a step: fields, history and clock together. Replaying a
rewound step from its own snapshot reproduces it bit for bit, where
re-running the script does not — warm starts and preconditioner reuse are
solver history outside model state, so two independent runs diverge at 1e-13
from the first step. That asymmetry is why a step that misbehaved can only be
looked at twice this way.

Where a run's transcript lands

On by default. A run that takes a step leaves

transcripts/2026-09-11T14-32-05-my_model/
    my_model.py          the script that launched it, verbatim
    launch.json          argv, interpreter, cwd, version, commit
    transcript.log       one aligned line per step, flushed

The account is only worth having on the run you did not prepare for — requiring
opt-in defeats the case it exists for. The stamp is why it is a stamp: the run
you want is the one from this morning. A directory rather than loose files
because a working directory full of logs and script copies invites mass
deletion.

launch.json is the honest answer to reproducibility. A programmatic launcher
cannot be made reproducible by fiat, but what was actually run can be written
down: command line, interpreter, cwd, package version, and the commit id with a
dirty flag when cwd is a repo. Only the entry script is copied — anything it
imports is not, and the file says so.

Three things make the default tolerable, each tested: nothing is created
until the first step opens, so an import or a mesh-only script leaves no
trace; it is off under pytest; and it is switchable —
model.transcript_file = None, an explicit path (which bypasses the stamped
directory and gets no launch record), or UW_TRANSCRIPT=off /
UW_TRANSCRIPT=/scratch/runs.

One aligned line per step, appended and flushed as the step closes, so
tail -f follows a running job.

# underworld3 run transcript · model 'default' · started 2026-09-10T21:22:40+00:00
# scales: length 2.2e+06 m | time 4.84e+18 s | mass 1.065e+47 kg | temperature 2500 K
# step           t/Myr          dt/Myr    wall/s  outcome    operators, in order
      3         1.51459        0.523655      0.09  ok         [convect] solve:SNES_AdvectionDiffusion_Composed(T) > history_shift:EulerianSUPG(T) > solve:SNES_Stokes(v)
      4         30.5663         29.0517      0.09  ABANDONED  [too big] solve:SNES_AdvectionDiffusion_Composed(T) > history_shift:EulerianSUPG(T) > solve:SNES_Stokes(v)
  -- restore from a snapshot; the clock now reads 1.51459 Myr
  -- rewind to the start of step 3 (t = 0.990939 Myr); 1 step(s) undone
      3         1.51459        0.523655      0.42  ok         [replay] solve:SNES_AdvectionDiffusion_Composed(T) > history_shift:EulerianSUPG(T) > solve:SNES_Stokes(v)

A true log records the backtracks, so rewind() and a bare load_state() each
write their own line — a log that shows step 3, then step 3 again with nothing
in between, is not a log of what happened. Four other things reach the file that
the bounded in-memory transcript does not keep: an abandoned step, a step aged out
by transcript_limit, and everything up to a kill.
wall/s is how long the block took; a step that suddenly takes ten times as
long is the first sign of a solver in trouble.

A .jsonl suffix (or transcript_format = "jsonl") writes the same record as one
JSON object per line, read back with uw.read_transcript. JSON lines rather than
YAML because one self-contained record per line is the point: a killed run
leaves a truncated final line that fails to parse, so the reader drops it and
keeps everything before, where a half-written YAML mapping frequently still
parses as a real record with its last key missing. YAML stays the right format
for a document written once and edited by hand, which is what Model.to_yaml
uses it for. The two formats differ deliberately: text is a report (one time
column, one unit, named in the header), JSON is a record (every value in the
units the run actually held it in).

The figure

A green-on-black terminal is not where a run belongs in a paper.

uw.transcript_diagram(model, out="figures/run.pdf")     # or a .jsonl log
uw.transcript_diagram(model, out="figures/run.svg")     # same figure, SVG
uw.transcript_flowchart(model)                          # Mermaid, for docs

Time runs down the page: one row per step, A4 portrait, paginated, so the
figure drops into a document column. The PDF is written directly — base-14
fonts, Flate-compressed content streams, a real cross-reference table — so a
run becomes a figure with nothing installed. .svg writes the same layout as
one continuous page.

The letter is the layout. Each distinct operator sequence gets one, defined
once at the foot:

step   t/Myr   dt/Myr   seq   dt
  11   9.856    1.920    A    ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇
  12     145    135.2    A    ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄  abandoned
  11   9.856    1.920    A    ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇
  12   12.56    2.704    B    ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇

A   AdvectionDiffusion(T)  >  shift EulerianSUPG(T)  >  Stokes(v)       14 steps
B   AdvectionDiffusion(T)  >  shift EulerianSUPG(T)  >  Stokes(v)  >
    AdvectionDiffusion(T)  >  shift EulerianSUPG(T)  >  Stokes(v)        1 step

A column of A with a single B in it says at a glance that one step did
something different. A hundred spelled-out sequences say nothing and hide the
one that matters. Backtracks are drawn in the left gutter as the path the run
took — a dashed arrow up from the step it bailed out of to the step it returned
to, then a solid arrow down to the row that takes that step again — and the
dt axis goes
logarithmic when the range exceeds 20x and says so, since a rejected step is
often tens of times the accepted ones — which is why it was rejected — and on
a linear axis it flattens everything else to nothing.

model.clear_transcript() starts a new run's account. A driver that runs the same
model many times — an inversion, a parameter sweep, a restart — otherwise gets
the concatenation of every run the process has done, and rewind() walks back
into the previous one.

Nothing is compulsory. A script that never opens a step behaves exactly as
before and the recording calls are no-ops.

Naming

A transcript is what was actually played, false starts and re-takes
included, which is precisely what this holds — and it carries its own contrast
with the score, the structure a run is supposed to repeat. "Journal" said
only "a log of some kind", and collides with the publishing sense in a
scientific codebase. record is kept as a verb: a step records what it
did; the thing it produces is the transcript.

docs/developer/design/run-score-and-transcript.md is the design note behind
that — the vocabulary (bar, beat, part, note, rest, tuplet, tie, tempo), what
each makes checkable, and the three things still missing: barriers are not
events, rests are not recorded, and a bar is not necessarily an interval
although model.step(dt) insists that it is. No implementation in it.

Two worked cases

The pattern is written up in HOW-TO-WRITE-UW3-SCRIPTS.md, and two runnable
cases exercise it in different geometries.

docs/examples/convection/intermediate/Ex_Convection_Annulus_Recorded.py —
Boussinesq convection in a 2D annulus. Four reference quantities fix all four
dimensions, and the buoyancy is then written as the force it is
(-RHO0 * ALPHA * GRAVITY * T * rhat) so the Rayleigh number falls out of the
nondimensionalisation rather than being typed in; the script prints
Ra = 2.585e+06 as a check. Rotated free-slip on the curved boundaries, and a
varying estimate_dt() that goes straight into model.step(dt). After the
loop it demonstrates the four things the record buys, in order: the transcript, a
rejected step that leaves the clock where it was, a bit-exact replay
(max |dT| = 0.000e+00), and a step taken twice showing as its own
operator sequence.

Compare ../advanced/Ex_Convection_Cylinder.py, which solves the same physics
with a bare for step in range(n) loop and no clock at all.

An adjoint driven from the transcript (a sinking blob under Stokes flow, not
in this repo). The backward pass of a discrete adjoint needs exactly what the
record holds: the state at each step and the order the operators ran in.
Walking model.transcript backwards — load_state(entry.snapshot), replay,
transpose-solve — replaced the hand-written ck = {"B", "V", "P", "dt"}
checkpoint dictionary that adjoint carried, and with it the dictionary's
dependence on knowing in advance which arrays the backward pass would want.
Taylor test 1.00000 / 1.00003, the same quality as the hand-taped version. It
is also what surfaced three of the defects below.

The defects

mesh.t was silently zero inside every solve (#410). It was bound to
PETSc's petsc_t, which the high-level solve() wrappers never set, and
solve(time=...) was accepted and ignored. A time-dependent boundary condition
written as sin(omega * mesh.t) — the usage mesh.t's own docstring
advertises — was identically zero, and nothing warned. Now a live-rampable
constants[] atom carrying model.tracker.time, repacked in
_update_constants, which is the single choke point every solver already
passes before solving. No kernel is recompiled per timestep and no solver
needed its own hook. Follows the #410 ruling rather than repairing the PETSc
plumbing that ruling records as tried and failed.

A dimensional value on the tracker was dropped by the on-disk snapshot.
A pint quantity fell through _serialise_field's unserialisable branch, was
written as __skipped, and was absent after load_state. With reference
quantities set, estimate_dt() is dimensional and so is the natural clock — so
the entry the pattern asks for was exactly the one a restart lost, while plain
floats beside it survived. Now stored as magnitude plus unit string.

A constants[] slot that stopped being constant was packed as 0.0. This is
the "rampable constant in exponent position does not ramp" report, and the
mechanism is not what that report assumed. The atom is not compiled out.
(1 + T**2)**(-m) + 1 is the NUMBER 2 while m is zero, so the whole
diffusivity banks as one constant and the collector stops there. Ramp m and
it depends on T again, the slot can no longer be reduced, and the solve
received a zero diffusivity: measured 2.0, then 0.0, then 0.0, with a
DIVERGED_LINEAR_SOLVE and nothing saying why. Now raises, naming the slot,
showing its content, and giving the two lines that force a rebuild.

A free surface whose derived lids had drifted solved anyway. held and
consistent are separate Stokes solvers configured once, at construction.
Change the free solve afterwards and they keep the old values: measured
viscosity 1000 against 1, body force -5 against -1, tolerance 1e-11 against
1e-6. h_inf, the equilibrium the surface relaxes toward, comes from the HELD
solve — so the surface relaxed toward an equilibrium computed with stale
physics while the free solve used the new. solve() now refuses and names what
drifted.

An abandoned step held on to its snapshot. Nothing could reach it — the
record never joins the transcript — so it was field-sized memory pinned until the
exception's traceback was collected.

A snapshot rescaled the mesh whenever units were active. mesh.X.coords is
the unit-aware view and returns metres once a model declares a length scale;
the DM coordinate vector _deform_mesh writes back into holds model units.
Capture read the first and restore wrote the second, so every restore
multiplied the mesh by the length scale — a 500 km box came back
250,000,000 km across, compounding on each round trip. Nothing raised: shapes
matched and every field was restored correctly, so only the geometry was wrong.
The visible symptom is uw.function.evaluate returning the value at one corner
for every sample point, because every sample point is now outside the domain.
model.rewind() goes straight through that path, which is how it surfaced. The
swarm path was unaffected — it captures the raw DMSwarmPIC_coor.

The transcript counted one operator as two. The hook lives in
_update_constants, which every solver passes on its way to a solve — except
that solve_rotated_freeslip pushes constants a second time for its own
assembly, after the public solve() has already announced the solver. So a
Stokes solve on a curved boundary recorded twice, in a record whose whole value
is that it says what ran. _update_constants now takes record=False for a
setup push.

estimate_dt() lost its units under the pattern's own idiom.
dt = fraction * solver.estimate_dt() came back as a bare float whenever the
estimate was a Python float rather than a numpy scalar: np.squeeze promotes
it to a 0-d array, uw.dimensionalise maps an array to a UnitAwareArray, and
that drops its units under arithmetic. The guard for exactly this already
existed in _dimensionalise_dt but sat on the no-units branch, so
Stokes.estimate_dt (numpy scalar) was fine and AdvDiffusion.estimate_dt
(Python float, accuracy basis) was not — the two silently disagreed.

Recording, not judging

The step's job is to record faithfully. It does not decide whether what it
recorded was a mistake.

An earlier revision of this branch asserted one rule in the loop — a history
must advance exactly once per step — and warned when it did not. It has been
removed. It claimed the authority of an invariant for a heuristic about usage;
it was already known to be wrong for a swarm that sub-cycles; it sat at the
wrong layer, since what knows a step was taken twice is the dataflow rather
than the container; and it could only ever see one step, whereas the check that
matters most cannot be seen from one — #423's signature is a growth rate of
~10% per cycle.

So the line now falls between structural checks — is the transcript
well-formed? a step cannot nest; a rewind cannot reach a step that kept no
snapshot — which stay in the loop and raise; and findings — does what was
recorded look wrong? — which belong to a pass over a finished transcript. That
pass can look across steps, can be re-run on an old transcript when a new
pathology is learned, and never has to decide mid-run whether something was
deliberate.

The recording is unchanged: both shifts of a doubled step are in the
transcript, in order. What went is the judgement.

Notes for review

The guards all have negative controls, which matter more than the positive
cases: an ordinary one-solve-per-step loop is asserted silent; an untouched
free-surface manager solves; an ordinary constant still ramps; and a parameter
whose value is an expression over mesh variables still tracks when the FIELD
changes, since symbolic sharing is correct and only re-assignment drifts.

The free-surface guard compares against what copying the parameter today
would produce, applying the same velocity rebinding _copy_constitutive_model
applies. A nonlinear rheology is deliberately rebound onto each derived
solver's own unknowns, so those expressions are meant to differ textually — the
first version of the guard failed test_1070's nonlinear-viscosity case for
exactly that reason.

test_1074 is added to scripts/test.sh beside test_1072; the test_107*
group is not otherwise batched and it would never have run.

Related: #708 shares a trigger with the constants defect — a zero at compile
time — but not the mechanism. There the victim is a raw zero history matrix, so
no slot exists to lose and the term is arithmetically absent. That PR's fix,
populating before compiling, is the right one for it.

Known gaps, not addressed here

A snapshot still cannot restore across a mesh deform or adapt, so
model.rewind() cannot reach a step on a moving or adapting mesh. The run
warns once and keeps recording. That is the mesh-rebuild-on-restore path
already planned for v1.2.

FreeSurface is not a snapshot state-bearer: _h_inf and _conserve_target
are not captured, and _pending_v_mesh_disp is not a field of
DDtSemiLagrangianState. Filed rather than fixed.

uw.pprint's default clean_display=True rewrites the string it is given —
three passes of a brace-stripping regex plus \s+ collapsed to a single space
— so the project's recommended rank-safe print destroys indentation, column
alignment, and any literal braces. The annulus example works around it with
clean_display=False. Narrowing the cleaner is not free, since the documented
use is uw.pprint(f"Expression: {expr}") and so it has to run on f-strings
too. Filed, not fixed here.

🤖 Generated with Claude Code

https://claude.ai/code/session_01Na7qBenCp67rDTZhFGTh5V

[Copilot]

🟡 Changes recommended

There are a few correctness and reliability issues in newly-added guards/docs (string-based drift comparison, a brittle error path, and a doc inconsistency) that should be fixed before merging.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Pull request overview

This PR introduces an explicit timestep boundary (with model.step(dt):) as a transactional unit of work, records step activity in a journal (optionally with restorable snapshots for rewind/replay), and addresses several previously-silent defects discovered while building that mechanism (time-dependent mesh.t, disk snapshot of dimensional tracker values, constants-slot collapse, and FreeSurface derived-solver drift).

Changes:

• Add ModelStep journaling + model.step(dt) context manager, plus optional record_every snapshots and model.rewind() support.
• Make mesh.t a live rampable constant synced from model.tracker.time before each solve; hook solve + history-shift events into the step journal.
• Extend disk snapshot serialization to round-trip dimensional quantities (magnitude + unit string), add regressions, and wire new FreeSurface drift guard + constant-slot guard tests into the test runner.

File summaries

File	Description
tests/test_1074_free_surface_config_drift.py	Regression coverage for FreeSurface derived-solver drift refusal + negative controls
tests/test_0104_constant_slot_still_constant.py	Regression coverage for constants[] slot “stops being constant” error + recovery path
tests/test_0011_model_step_journal.py	Regression coverage for model.step semantics, journaling, recording, and rewind/replay invariants
tests/test_0009_model_tracker.py	Regression coverage for mesh.t tracking model clock and dimensional tracker round-trip on disk snapshots
src/underworld3/utilities/unit_aware_coordinates.py	Make time-units patching tolerant of new mesh._t implementation
src/underworld3/utilities/_jitextension.py	Refuse to silently pack 0.0 into a stale constants[] slot; raise diagnostic error
src/underworld3/systems/solvers.py	Fix source-term setter to avoid baking UWexpression atoms by mistaken quantity duck-typing
src/underworld3/systems/free_surface.py	Add derived-solver drift detection and refusal before solving
src/underworld3/systems/ddt.py	Emit journal events when histories shift to detect “advanced twice per step” cases
src/underworld3/model.py	Add step journaling/recording/rewind API and model.step(dt) transactional context manager
src/underworld3/discretisation/discretisation_mesh.py	Implement mesh.t as model-time-backed constants[] atom; sync it from model before solves
src/underworld3/cython/petsc_generic_snes_solvers.pyx	Central hook: sync mesh.t + record solve events in the open step journal
src/underworld3/checkpoint/disk_snapshot.py	Add dimensional-quantity serialization (magnitude + unit string) for disk snapshots
scripts/test.sh	Ensure new test_1074 is actually executed by the script runner
docs/developer/guides/HOW-TO-WRITE-UW3-SCRIPTS.md	Document the new timestepping pattern, journaling/recording/rewind usage, and mesh.t behavior

Review details

• Files reviewed: 15/15 changed files
• Comments generated: 4
• Review effort level: Lite

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.