diff --git a/docs/developer/design/run-score-and-transcript.md b/docs/developer/design/run-score-and-transcript.md new file mode 100644 index 000000000..926c57128 --- /dev/null +++ b/docs/developer/design/run-score-and-transcript.md @@ -0,0 +1,171 @@ +# The run score and the run transcript + +*Status: design note. Vocabulary and the model it implies. No implementation.* + +## Why this note exists + +Underworld3 records what a run did — `model.transcript`, the on-disk log, the +figure. Building those made it clear that the hard part is not the recording +but the **view**: the record has to make a *relation* visible, and until we +knew which relation, each addition was a patch on the last. + +The relation is this. A value stored by one action in one cycle is picked up by +a read in the next. Every silent defect found while building the timestepping +machinery was a read taking the wrong write: + +- an Eulerian history that initialises only on its first solve, so a solver + reused for a second run silently reads the **previous run's** store +- the history-before-advection ordering trap +- `old_frame_traceback` on a deforming mesh (#423), where a semi-Lagrangian + history recorded in the old frame and read after the mesh moved amplifies + ~10% per cycle, with no error and no symptom other than the growth rate +- reading a solver's `F0`/`F1` templates before the solve configured them + +None of those is visible in a single cycle. All of them are structure across +cycles. + +## Two documents, not one + +**The score** is what the run is *supposed* to do: which parts play, in what +order, with what meter. It is the same for every cycle of a well-behaved run. + +**The transcript** is what the run *actually did*, including the things that +are not in the piece — a cycle abandoned, a jump back to an earlier cycle, a +re-take. + +Both already exist in the current implementation without those names. The +figure's operator-sequence legend (`A`, `B`, ...) is an **inferred score**; the +rows are the **transcript**. A run in which every bar is `A` played the score. +A `B` is a bar played differently. + +That naming makes the original motivating question mechanical: *is the model +doing the scientific task you say it is* becomes **a diff of the transcript +against the score**. + +## Vocabulary + +The terms below are where the ideas came from, and they are deliberately NOT +the words the code and the figures use. "Bar", "note", "rest" and "simile" were +useful for getting the model right and are forced as public vocabulary: what +ships says **step**, **part**, **did nothing**, and **unchanged**. The mapping +is kept here because the reasoning depends on it — each musical term carries a +convention that is the reason the corresponding decision was made. + + +| term | meaning here | +|---|---| +| **bar** → *step* | one turn of the orchestrating loop. Numbered monotonically, never reused, and the thing you refer to. NOT necessarily a physical time interval — it may be a task. | +| **beat** | a position inside a step at which alignment is required. **Barriers** sit on beats: a mesh deform, an adapt, a migration, a remesh. | +| **part** (kept) | a participant with its own column. Two kinds: *actors* (solvers, swarm pushes, mesh movers) and *state-holders* (DDt histories, fields, particle coordinates). | +| **note** → *a mark* | what a part did in a step, carrying its own duration — its `dt`, which need not be the bar's. | +| **rest** → *did nothing* | notated absence. Distinguishes *did nothing this bar* from *was not being watched*. | +| **tuplet** | n notes in the space of the bar, bracketed with the ratio. Sub-cycling, notated as ordinary rather than flagged as anomalous. | +| **tie** | a value written in one step and read in the next, drawn as an arc across the barline. | +| **tempo** | how step numbers map to real time. Deliberately separate from the meter: `dt` varies, wall clock varies more, and the vertical axis is ordinal. | +| **performance event** | not part of the piece: a bar abandoned, a jump back, a re-take. Transcript only. | + +Two conventions borrowed with the vocabulary and worth keeping: + +- **Score order.** Parts appear in a stable, conventional order (actors, then + histories, then swarms, then mesh), not order of first appearance, so a + reader finds the same part in the same place in every run's transcript. +- **Rests are compulsory.** In a score every part accounts for every beat. A + part that does nothing is written as a rest, never left blank. + +## Why two dimensions + +Different parts advance on different `dt`. A swarm may sub-cycle twice for one +Stokes solve; two histories on the same field may be at different orders. Those +cannot be laid on one axis without pretending they share a clock — which is the +failure being looked for. + +So: **parts across the page, bars down the page.** Vertical is ordinal, not +time. A part whose accumulated time drifts away from its neighbours' is then a +visible misalignment rather than something you would have to instrument for. + +## What the notation makes checkable + +Each check is a reading of the notation rather than a separate assertion: + +| reading | defect | +|---|---| +| a tie with no note at its head | a read of uninitialised history | +| two notes tied into one read | the physical step taken twice | +| a tuplet whose ratio does not fill its bar | sub-cycling that fails to tile the interval | +| a tie crossing a beat that carries a barrier, with nothing re-expressing it | the `old_frame_traceback` class (#423) | +| transcript ≠ score | the model is not doing what the script says | + +**None of these belong in the loop.** An earlier version asserted one of them +there — "a history must advance exactly once per bar" — and it was wrong twice +over: it would fire on legitimate sub-cycling, and it could not see the check +that matters most, since #423's signature is a growth rate across bars. It has +been removed. Derived from the notation the rule is the honest one — **the +notes in a bar tile its interval exactly once** — and sub-cycling satisfies it. + +The separation that follows: **execution records, analysis judges.** A step +raises only on structural failures, where the transcript could not be +well-formed — a step opened inside another step, a rewind to a bar that kept no +snapshot. Everything above is a *finding*, produced by a pass over a finished +transcript, which can look across bars, 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. + +## Where reads and writes come from + +Feasibility, because this is the part that decides whether the model is +buildable: + +- **Writes** come from the runtime hooks already in place: the unknown after a + solve, `psi_star[i]` in `update_post_solve`, particle coordinates after an + advection. +- **Reads** are derivable *symbolically*. A solve's residual is a SymPy form, + so the variables it depends on are in `F0` / `F1`'s atoms. This is how the + discrete adjoint obtained `∂F/∂ψ*` at all. No kernel instrumentation. + +Within one solve, reads and writes are not ordered — the kernel reads +`psi_star` throughout the Newton iteration — so the edge is *write → solve*, +not write → instant. The cell is atomic; the meaning is in the edges between +cells. That is exactly where the cross-cycle relation lives, so the limitation +does not bite. + +## What is missing today + +Everything above except two items is a renaming of something already captured, +which is a reasonable sign the vocabulary fits rather than being imposed. + +1. **Barriers are not events.** `_deform_mesh` does not declare itself, so + there is no beat to draw the rule at — and the check that most wants the + barrier (#423) has no anchor without it. +2. **Rests are not recorded.** A part that does nothing in a bar simply does + not appear, so silence and absence are indistinguishable. + +A third, smaller: **a bar is not necessarily an interval, and `model.step(dt)` +insists that it is.** The signature requires a `dt`, so there is no container +for "the next task". If the event clock is the general thing, the timestep is +the common case rather than the definition. + +## Inferred score, then declared score + +The score is **inferred** today — the figure takes the most common bar as the +norm. That costs nothing and can only ever say *this bar differs from its +neighbours*. + +A **declared** score — the script stating what a bar is supposed to contain — +turns the diff into *this run disagrees with its own description*, which is the +stronger claim and the one that motivated the work. The path is to work +towards the declarative model from the inferred one rather than to require it +up front; nothing in the vocabulary above depends on which we have. + +## A note on names + +"Transcript" rather than tape or record: a transcript is what was actually +played, including the false starts and the re-takes, which is precisely the +thing being kept. The word also carries its own contrast with the score, so +the pair names itself. + +The API followed: what was `model.journal` is `model.transcript`, and the +renderers follow from the pair — `transcript_diagram` draws the transcript, +`transcript_flowchart` draws the score it implies. + +"Record" is kept as a **verb**. A step records what it did; the thing it +produces is the transcript. diff --git a/docs/developer/guides/HOW-TO-WRITE-UW3-SCRIPTS.md b/docs/developer/guides/HOW-TO-WRITE-UW3-SCRIPTS.md index 1646bd499..875f255ea 100644 --- a/docs/developer/guides/HOW-TO-WRITE-UW3-SCRIPTS.md +++ b/docs/developer/guides/HOW-TO-WRITE-UW3-SCRIPTS.md @@ -14,9 +14,10 @@ This guide captures critical lessons learned from writing and debugging Underwor 3. [Units System Integration](#units-system-integration) 4. [Mesh and Variable Creation](#mesh-and-variable-creation) 5. [Solver Setup and Execution](#solver-setup-and-execution) -6. [Common Pitfalls and Anti-Patterns](#common-pitfalls-and-anti-patterns) -7. [Testing Best Practices](#testing-best-practices) -8. [Debugging Techniques](#debugging-techniques) +6. [The Timestepping Pattern](#the-timestepping-pattern) +7. [Common Pitfalls and Anti-Patterns](#common-pitfalls-and-anti-patterns) +8. [Testing Best Practices](#testing-best-practices) +9. [Debugging Techniques](#debugging-techniques) --- @@ -414,8 +415,478 @@ model = stokes.constitutive_model # Confusing with uw.Model --- +## The Timestepping Pattern + +### ⚠️ RULE: the clock lives on `model.tracker`, never in a loose variable + +Every time-dependent script needs a clock. Declare the model and its reference +quantities first (see [RULE #2](#rule-2-reference-quantities-before-mesh-creation) +— they must precede mesh creation), then keep the clock on the tracker: + +```python +uw.reset_default_model() +model = uw.get_default_model() +model.set_reference_quantities( + domain_depth=uw.quantity(500, "km"), + material_density=uw.quantity(3300, "kg/m**3"), + material_viscosity=uw.quantity(1e21, "Pa*s"), +) + +mesh = uw.meshing.UnstructuredSimplexBox(...) # inherits the reference quantities +# ... variables, solvers ... + +model.tracker.time = uw.quantity(0.0, "Myr") +model.tracker.step = 0 +model.tracker.dt = None + +while model.tracker.time < end_time: + dt = adv_diff.estimate_dt() # a dimensional quantity when units are active + + adv_diff.solve(timestep=dt) + stokes.solve(zero_init_guess=False) + + model.tracker.time = model.tracker.time + dt + model.tracker.step += 1 + model.tracker.dt = dt +``` + +**Start from the model, not from the mesh.** A default model is created for you +if you never ask for one, so it is easy to write a whole script without noticing +it exists — and then reference quantities, which must be set before the mesh, are +already too late. Declaring the model on the first line makes the units decision +explicit at the only point where it can still be made. `estimate_dt()` then +returns a dimensional quantity and the clock carries units with no extra work. + +Neither the model nor the units are enforced today; both arrived after much of +the surrounding code. Treat this ordering as the pattern regardless, because +retrofitting units to a script written without them means rebuilding the mesh. + +**Why this and not `t = 0.0; t += dt`.** The tracker is captured by +`model.save_state()` and reverted by `model.load_state()`. A loose Python +variable is not — the tracker's own docstring says so: + +> Everything on the tracker is captured by `snapshot` and reverted by +> `restore`; loose Python variables are not. + +So a script that keeps its clock in a local and then backsteps gets its +**fields** restored and its **time** left in the future. Nothing raises; the +run simply carries a clock that disagrees with the state it is describing, and +every output written from that point is mislabelled. Measured on a real +five-step run: after restoring a snapshot taken at t = 1140, a loose `t` still +read 2280 while the fields were correctly back at 1140. + +This holds for both snapshot flavours — the in-memory token and the on-disk +snapshot used for restart. + +`time`, `step` and `dt` are pre-seeded on the tracker as a convention. Anything +else you assign to it (`model.tracker.rms_velocity = ...`) is captured and +restored the same way, so a diagnostic you want to survive a backstep belongs +there too. + +### Wrap the step + +`model.step(dt)` makes the loop a transaction: + +```python +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) +``` + +Three things follow, and none of them requires anything else in the script to +change. The clock reads the END of the interval for the whole block, which is +where an implicit scheme centres its residual, so a time-dependent coefficient +is evaluated at the right time. The advance commits only on clean exit, so a +step that raises — or one abandoned because the Courant number came out too +large — leaves the clock exactly as it was. And everything the block did is +recorded: + +```python +>>> for entry in model.transcript[-3:]: +... print(entry) + solve:SNES_Stokes(V)> + solve:SNES_Stokes(V)> + solve:SNES_Stokes(V)> +``` + +That record is worth having on its own. It answers what a run actually did, +in order, without the script being instrumented for it — which is the question +you want to ask of someone else's model, or your own six months later. + +Opening a step is optional. A script that never does behaves exactly as before. + +### Recording a run + +Ask the step to keep the state it started from and the transcript becomes a +restorable record: + +```python +model.record_every = 1 # keep every step; None (default) keeps none +model.record_limit = 8 # how many snapshots to retain + +while model.tracker.time < end_time: + with model.step(dt): + adv_diff.solve(timestep=dt) + stokes.solve(zero_init_guess=False) + +model.rewind() # undo the last step: fields, history and clock +``` + +The snapshot is taken before the operators run, which is the only correct +point — a `DDt` shifts its history in its post-solve hook, so a snapshot taken +afterwards holds the shifted history rather than the step's input. + +Two things this buys beyond backstepping. Replaying a step from its own +snapshot reproduces it exactly, where re-running the script does not, so a step +that misbehaved can be looked at twice. And an adjoint needs precisely this: the +state at each step and the order the operators were applied in. + +Snapshots cost roughly 13 bytes per primary degree of freedom per step. Older +steps lose their snapshot and keep their transcript record, so the account of what +happened outlives the state it happened to. + +A driver that runs the same model more than once — an inversion, a parameter +sweep, a restart — should start each run with a clean account: + +```python +model.clear_transcript() +model.tracker.time = uw.quantity(0.0, "Myr") +model.tracker.step = 0 +``` + +Without it the transcript is the concatenation of every run the process has done, +and `rewind()` will walk back into the previous one. + +On a mesh that deforms or adapts the snapshot cannot be taken yet; the run +warns once, keeps recording, and `rewind()` will not reach those steps. + +### Writing the transcript down + +`model.transcript` is what the run can still undo. It lives in memory, it is +bounded, and it dies with the process. The transcript on disk is what the run +*did* — and it is **on by default**, because the account is only worth having +on the run you did not prepare for: + +``` +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 stamp is the point: the run you want is the one from this morning, and a +fixed filename would have overwritten it. A directory rather than loose files +because a working directory full of logs and script copies invites mass +deletion, which loses the one you needed. + +`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: the command line, the interpreter, the working directory, the package +version, and the commit id with a dirty flag if the work is under version +control. Only the entry script is copied — anything it imports is not, which is +what the commit id is there to cover. + +Three things keep the default tolerable. **Nothing is created until the first +step opens**, so an import, or a script that only builds a mesh, leaves no +trace. **It is off under pytest**, because 1800 tests should not each leave a +directory. And it can be turned off or sent elsewhere: + +```python +model.transcript_file = "output/run.log" # somewhere else, no launch record +model.transcript_file = "output/run.jsonl" # JSON lines instead +model.transcript_file = None # off +``` + +```bash +UW_TRANSCRIPT=off # off for the session +UW_TRANSCRIPT=/scratch/runs # put the stamped directories there +``` + +The file itself is one aligned line per step, appended and flushed as it +closes, so `tail -f` follows a running job: + +``` +# underworld3 step log · 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 + 0 0.175907 0.175907 0.44 ok [convect] solve:SNES_AdvectionDiffusion_Composed(T) > history_shift:EulerianSUPG(T) > solve:SNES_Stokes(v) + 1 0.501546 0.325639 0.09 ok [convect] solve:SNES_AdvectionDiffusion_Composed(T) > history_shift:EulerianSUPG(T) > solve:SNES_Stokes(v) + 2 0.990939 0.489393 0.09 ok [convect] solve:SNES_AdvectionDiffusion_Composed(T) > history_shift:EulerianSUPG(T) > solve:SNES_Stokes(v) + 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) +``` + +`wall/s` is how long the block took. It is not physics, but it is the number +you want when watching: a step that suddenly takes ten times as long is the +first sign of a solver in trouble. + +**A true log records the backtracks.** `rewind()` and a bare `load_state()` +each write their own line, because a log that shows step 3, then step 3 again +with nothing in between, is not a log of what happened. `rewind` writes the +more specific note and suppresses the generic one. + +Four other things go in the file that are not in `model.transcript`, all +deliberate: + +- **An abandoned step.** A rejected step is the part of a run's history that is + otherwise invisible, and it is usually what you want when asking why a run + went the way it did. +- **A step aged out by `transcript_limit`.** The account of what happened outlives + both the state and the bounded in-memory list. +- **Everything up to a kill.** The file is flushed per step. + +### For parsing: JSON lines + +A path ending `.jsonl`, `.ndjson` or `.json` — or `model.transcript_format = +"jsonl"` — writes the same record as one JSON object per line: + +```json +{"kind": "run", "model": "default", "started": "2026-09-10T21:22:40+00:00", "scales": {"length": {"magnitude": 2200000.0, "units": "meter"}, ...}} +{"kind": "step", "index": 0, "label": "convect", + "t0": {"magnitude": 0.0, "units": "megayear"}, + "t1": {"magnitude": 0.1759, "units": "megayear"}, + "dt": {"magnitude": 5551210433127.7, "units": "second"}, + "completed": true, "restorable": true, "wall": 0.44, + "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)"}, + {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 1.1469e-06}, + {"kind": "solve", "name": "SNES_Stokes(v)"}]} +{"kind": "rewind", "message": "...", "to_step": 3, "steps_undone": 1, "t": {"magnitude": 0.9909, "units": "megayear"}} +``` + +Read it back with `uw.read_transcript(path)`, which returns one entry per run — an +inversion driver that ran the forward model thirteen times leaves thirteen runs +in one file, delimited by the header `clear_transcript()` writes. + +**Why JSON lines and not YAML.** One self-contained record per line is the +whole point. A killed run leaves a truncated final line that *fails* to parse, +so `read_transcript` drops it and keeps everything before; a half-written YAML +mapping frequently still parses, as a real record with its last key missing. +Line-oriented also means `grep`, `wc -l` and `jq -c` work without a parser, and +`json` is stdlib with predictable float round-tripping. YAML is the right +format for a whole document written once and edited by hand — which is what +`Model.to_yaml` uses it for — but a log is a stream. + +The two formats differ in one more way. The text log is a **report**: the time +column is converted into one unit, named in the header. The JSON log is a +**record**: every value keeps the units the run actually held it in, which is +why `t0` may read in Myr beside a `dt` in seconds — the clock came from the +tracker and the interval from `estimate_dt()`. + +Rank 0 writes; the other ranks record in memory as usual. + +### The same account, as a figure + +A terminal is not where a run belongs in a paper. + +```python +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 +``` + +`transcript_diagram` puts **time down the page**: one row per step, A4 portrait, +paginated, so it drops into a document column and opens anywhere. Each row +carries the step index, the clock, `dt` as a number and as a bar, a wall-clock +tick, and one letter. + +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 whose state it returned +to, then a solid arrow **down** from there to the row that takes that step +again. The pair is what makes a repeated step index read as a repeat rather +than a typo. Two calls that make the same jump — a `load_state` and then a +`rewind` to the same place — are one backtrack in the run's story and one arrow +on the page. + +The PDF and the SVG are both written directly — no plotting library, no +rasterisation, nothing fetched at render time, and a print-safe palette that +separates in greyscale. The `dt` axis goes logarithmic when the range exceeds +20x and says so: 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. + +**The letter is the layout.** Each distinct operator sequence gets one, defined +once at the foot of the figure: + +``` +step t/Myr dt/Myr seq dt + 9 6.472 1.176 A ▇▇▇▇▇▇▇▇▇▇▇▇ + 10 7.936 1.464 A ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ + 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. + +The two columns are independent, which is worth reading carefully: **`seq` is +what the step ran; `ok` / `abandoned` is whether it was kept.** In the figure +above the abandoned step ran the ordinary sequence `A` and was then rejected by +a check in the script — nothing failed. `B` is the same three operators run +twice inside one step block. The figure reports that it differs and makes no +claim about whether it is wrong. + +`transcript_flowchart` renders one step's operator flow as Mermaid. When a run has +more than one distinct sequence, each becomes its own subgraph labelled with +the steps that took it, so an anomalous step is visible rather than averaged +away. + +Both accept a live model, a `.jsonl` log, or the list `read_transcript` returns. +Not a text log: that one is a report, and reading it back is refused with the +one line that fixes it. + +### Recording, not judging + +The step's job is to record faithfully. It does not decide whether what it +recorded was a mistake. + +Two kinds of check are easy to confuse, and only one belongs in the loop. + +**Structural checks** ask whether the transcript is well-formed — a step cannot +be opened inside another step; a rewind cannot reach a step that kept no +snapshot. These cannot legitimately fail, so they raise, immediately. + +**Findings** ask whether what was recorded looks wrong. They belong to a pass +over the transcript, after the run. That is not a deferral for convenience; it +is where they can actually be computed: + +- A finding may need to look **across bars**. The free-surface instability in + `#423` is a history recorded in the old frame and read after the mesh moved, + growing about 10% per cycle. Its signature *is* the growth rate, so no + per-step check can see it at all. +- A finding may be **wrong about what is legitimate**. "A history advanced + twice in this bar" is a mistake when a step was taken twice and perfectly + correct when a swarm sub-cycles. Inside the loop that has to be guessed; + over the transcript it is a question about whether the operations tile the + bar's interval. +- A finding can be **re-run on an old transcript** when a new pathology is + learned. A warning fired at run time cannot. + +So the transcript records that a history shifted twice, with both shifts in +order, and says nothing about it. Reading that is +`docs/developer/design/run-score-and-transcript.md`'s subject, and the analysis +pass it describes is not built yet. + +One thing worth knowing while it is not: a history advances on **every** solve, +whether or not that call passed a timestep — omitting it reuses the last value. +There is no "solve without advancing" switch, so a corrector or a Picard +iteration on a coupled system has to put the history back between passes: + +```python +saved = copy.deepcopy(adv_diff.Unknowns.DuDt.state) +adv_diff.solve(timestep=dt) # the extra pass +adv_diff.Unknowns.DuDt.state = saved +``` + +That a coupled iteration inside one step needs this is a gap in the library, +not a rule the user broke. + +### Backstepping + +The pattern above is what makes speculative stepping safe: + +```python +snap = model.save_state() # before the step, not after + +try: + with model.step(big_dt): + adv_diff.solve(timestep=big_dt) + stokes.solve(zero_init_guess=False) + if courant_number() > courant_limit: + raise StepRejected # abandons the step; the clock stays put +except StepRejected: + model.load_state(snap) # fields go back; the clock never moved + for sub_dt in substeps(big_dt): + with model.step(sub_dt): + ... +``` + +Take the snapshot **before** the operator, not after. A `DDt` history plugin +shifts its history in its post-solve hook, so a snapshot taken after a solve +holds the shifted history, which is not the state that step ran from. + +Restoring a snapshot is bit-exact and repeatable. Re-running the same script is +not: warm starts and preconditioner reuse are solver history that is not part +of model state, so two independent runs of the same problem on the same solver +objects diverge at the 1e-13 level from the first step. If you need to look at a +step twice, restore it rather than re-run it. + +### Two worked cases + +**`docs/examples/convection/intermediate/Ex_Convection_Annulus_Recorded.py`** — +Boussinesq convection in an annulus. Four reference quantities, a body force +written as a force (Ra falls out of the nondimensionalisation rather than being +typed in), rotated free-slip on the curved boundaries, and a varying +`estimate_dt()`. It then demonstrates the four things the transcript buys, in +order: the transcript, a rejected step, a bit-exact replay, and a step taken +twice showing up 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.** The backward pass of a discrete adjoint +needs exactly what the transcript holds: the state at each step and the order the +operators were applied in. Walking `model.transcript` backwards — +`load_state(entry.snapshot)`, replay, transpose-solve — replaces the +hand-written checkpoint dictionary that an adjoint normally carries, and +removes its dependence on knowing in advance which arrays the backward pass +will want. + +### Time-dependent expressions + +`mesh.t` is the model clock as a symbol. It is repacked from +`model.tracker.time` before every solve, so a time-dependent source or +boundary condition follows the loop above with no recompilation per step: + +```python +omega = 2 * sympy.pi / period +stokes.add_dirichlet_bc((V0 * sympy.sin(omega * mesh.t), 0.0), "Top") +``` + +Two things to know. A script that never advances `model.tracker.time` leaves +`mesh.t` at zero, so the clock and the pattern above are the same subject. And +`mesh.t` should appear inside an expression rather than be handed bare to a +scalar setter — `poisson.f = mesh.t` stores a value, `poisson.f = 1.0 * mesh.t` +keeps the symbol. + +--- + ## Common Pitfalls and Anti-Patterns +### ❌ Rebinding the name instead of setting `.sym` + +To change the value of an expression, set `.sym`. It is the only settable +property — `.value` and `.data` are derived, read-only views. + +```python +# ✅ CORRECT - a value change; the container keeps its identity +viscosity.sym = sympy.Integer(0) +solver._update_constants() # only if you are not about to solve + +# ❌ WRONG - rebinds a Python name and changes nothing +viscosity = 0 +``` + +The second line leaves every expression that already references the atom +pointing at the old object with its old value, and nothing complains. The +identity is the point: because the container is unchanged, a ramped value +reaches every residual that mentions it with no rebuild. + +`expr.copy(other)` does the same job from another expression, and assigning to +a constitutive parameter slot (`Parameters.diffusivity = 0.0`) is also a value +change rather than a replacement. + ### ❌ Swarm Variable Creation After Population ```python @@ -694,6 +1165,15 @@ TypeError: unsupported operand type(s) for *: 'UnitAwareDerivativeMatrix' and 'N - [ ] If using units: Set reference quantities BEFORE mesh creation - [ ] If using units with [M]: Provide material_density or equivalent +### Writing a Timestepping Loop + +- [ ] Declare the model and its reference quantities BEFORE creating the mesh +- [ ] Keep `time`, `step` and `dt` on `model.tracker`, not in local variables +- [ ] Wrap each step in `with model.step(dt):` +- [ ] To change an expression's value set `.sym`, never rebind the name +- [ ] Take snapshots BEFORE the operator you might want to undo +- [ ] Use `mesh.t` inside an expression for time dependence, never bare + ### Creating a Swarm - [ ] Create mesh first @@ -727,6 +1207,17 @@ TypeError: unsupported operand type(s) for *: 'UnitAwareDerivativeMatrix' and 'N ## Version History +- **2026-09-09**: The timestepping pattern + - Start from the model and its reference quantities, not from the mesh + - Clock on `model.tracker`, not loose variables (snapshot consistency) + - Disk snapshots now carry dimensional values (magnitude + units) + - `mesh.t` now resolves to the model clock (#410) + - `model.step(dt)` — the step as a transaction, and the step transcript + - `model.record_every` / `model.rewind()` — the transcript as a restorable record + - A step warns when a history advances more than once + - Set `.sym` to change a value; rebinding the name changes nothing + - Backstepping recipe; snapshot before the operator + - `mesh.t` is not the model clock and is silently zero in a solve - **2025-11-15**: Initial version - Swarm ordering rules from test_0850/0851 debugging - Units everywhere-or-nowhere principle @@ -744,3 +1235,5 @@ TypeError: unsupported operand type(s) for *: 'UnitAwareDerivativeMatrix' and 'N - `docs/developer/COORDINATE-UNITS-TECHNICAL-NOTE.md`: Coordinate units implementation - `docs/beginner/tutorials/12-Units_System.ipynb`: Units system tutorial - `docs/beginner/tutorials/13-Non_Dimensional_Scaling.ipynb`: Dimensional analysis +- `docs/advanced/snapshot-restore.md`: Snapshot and restore semantics +- `tests/test_0009_model_tracker.py`: The pattern, enforced diff --git a/docs/developer/index.md b/docs/developer/index.md index a8c1f5123..938131827 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -145,6 +145,7 @@ guides/mpi-hang-supervision :hidden: :caption: Design Documents +design/run-score-and-transcript design/UNITS_SIMPLIFIED_DESIGN_2025-11 design/ND_UNITS_BOUNDARY_CONTRACT design/WHY_UNITS_NOT_DIMENSIONALITY diff --git a/docs/examples/convection/README.md b/docs/examples/convection/README.md index 72717925f..95d5621ac 100644 --- a/docs/examples/convection/README.md +++ b/docs/examples/convection/README.md @@ -42,6 +42,15 @@ Thermal convection combines heat transfer and fluid mechanics to model buoyancy- - Multiple convection cells and interactions - Introduces: domain geometry effects, cell interactions +7. **Annulus Convection, Recorded** - `Ex_Convection_Annulus_Recorded.py` + - Boussinesq convection in a 2D annulus, written in the timestepping pattern + - Reference quantities first, so the buoyancy is written as a force and the + Rayleigh number falls out of the nondimensionalisation + - Rotated free-slip on the curved boundaries; a varying `estimate_dt()` + - Demonstrates the run's own record: the transcript, a rejected step, a + bit-exact replay, and a step taken twice showing in the record + - See `docs/developer/guides/HOW-TO-WRITE-UW3-SCRIPTS.md` + ### 🎓 Advanced Examples (`advanced/`) **Complex convection systems.** diff --git a/docs/examples/convection/intermediate/Ex_Convection_Annulus_Recorded.py b/docs/examples/convection/intermediate/Ex_Convection_Annulus_Recorded.py new file mode 100644 index 000000000..3411eddc4 --- /dev/null +++ b/docs/examples/convection/intermediate/Ex_Convection_Annulus_Recorded.py @@ -0,0 +1,518 @@ +# --- +# jupyter: +# jupytext: +# formats: py:percent +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# kernelspec: +# display_name: Python 3 +# language: python +# name: python3 +# --- + +# %% [markdown] +""" +# Convection in an Annulus — a recorded run + +**PHYSICS:** convection +**DIFFICULTY:** intermediate + +## Description + +Boussinesq thermal convection in a 2D annulus, written in the timestepping +pattern: the model and its reference quantities come first, the clock lives on +`model.tracker`, and each timestep is a `model.step(dt)` block. + +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. The physics here is +unchanged. What the pattern adds is that the run keeps an account of itself, and +that account is worth four things this script demonstrates in turn: + +1. **what ran** — an ordered transcript, named by what each solver solves +2. **a rejected step** — the clock does not move when a step is abandoned +3. **playback** — a recorded step replays bit-for-bit, where a re-run does not +4. **recording without judging** — a step taken twice is visible in the + transcript, and the transcript makes no claim about whether that is wrong +5. **a log on disk** — the same account in aligned columns, flushed as each + step closes, so a run that dies keeps its history and a run in progress can + be watched with `tail -f` +6. **a figure** — the same account as a portrait PDF (or SVG), and the step's + operator flow as Mermaid + +## Key concepts + +- **Units first.** Four reference quantities fix the scaling, and the body + force is then written as the physics — `-rho0 alpha T g rhat` — rather than + as a Rayleigh number. Ra falls out of the nondimensionalisation; the script + prints it so you can check. +- **Rotated free-slip on curved boundaries.** `add_rotated_freeslip_bc` + enforces `v.n = 0` to machine precision on a circle, where a penalty or + Nitsche condition leaks at ~1e-3. +- **A varying timestep.** `estimate_dt()` returns a dimensional quantity that + goes straight into `model.step(dt)` and `adv.solve(timestep=dt)`. With an + implicit SUPG transport this is an accuracy choice, not a stability limit. + +## Parameters + +Override from the command line, e.g. `-uw_n_steps 20 -uw_cell_size 0.075`. +""" + +# %% +import os + +import numpy as np +import sympy + +import underworld3 as uw + + +def say(*args): + """Rank-safe print that keeps its own formatting. + + `uw.pprint`'s default `clean_display=True` rewrites the string it is given + — it strips braces and collapses runs of whitespace — so an aligned table + printed through it loses its columns. Pass `clean_display=False` whenever + the layout is yours rather than SymPy's. + """ + uw.pprint(*args, clean_display=False) + + +params = uw.Params( + uw_cell_size=0.1, # mesh resolution, as a fraction of the outer radius + uw_n_steps=8, # timesteps in the recorded run + uw_dt_fraction=0.5, # accuracy factor on estimate_dt() + uw_demos=1, # run the transcript demonstrations after the loop +) + +# %% [markdown] +""" +## The model comes first + +Reference quantities must be set BEFORE the mesh is created, so the model is +the first thing the script declares rather than something the mesh conjures for +you. Four quantities fix all four dimensions this problem uses: + +| quantity | fixes | +|---|---| +| `shell_thickness` | length | +| `thermal_diffusivity` | time, as `d^2 / kappa` | +| `mantle_viscosity` | mass | +| `temperature_contrast` | temperature | +""" + +# %% +uw.reset_default_model() +model = uw.get_default_model() + +SHELL_THICKNESS = uw.quantity(2200, "km") +KAPPA = uw.quantity(1e-6, "m**2/s") +ETA = uw.quantity(1e22, "Pa*s") +DELTA_T = uw.quantity(2500, "K") + +RHO0 = uw.quantity(3300, "kg/m**3") +ALPHA = uw.quantity(3e-5, "1/K") +GRAVITY = uw.quantity(9.81, "m/s**2") + +model.set_reference_quantities( + shell_thickness=SHELL_THICKNESS, + thermal_diffusivity=KAPPA, + mantle_viscosity=ETA, + temperature_contrast=DELTA_T, +) + +RAYLEIGH = (RHO0 * ALPHA * DELTA_T * GRAVITY * SHELL_THICKNESS**3 / (KAPPA * ETA)) +say(f"Ra = {float(RAYLEIGH.to('dimensionless').magnitude):.3e}") +say(f"diffusion time d^2/kappa = " + f"{model.get_fundamental_scales()['time'].to('Gyr')}") + +# %% [markdown] +""" +## Mesh and variables + +An annulus with `radiusInner / radiusOuter = 0.55`, roughly Earth's +core-mantle ratio. The mesh is built in model units; `mesh.X.coords` reads back +in metres because the model declares a length scale. +""" + +# %% +R_OUTER = 1.0 +R_INNER = 0.55 + +mesh = uw.meshing.Annulus( + radiusInner=R_INNER, + radiusOuter=R_OUTER, + cellSize=params.uw_cell_size, + degree=1, + qdegree=3, +) + +v = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=2) +p = uw.discretisation.MeshVariable("p", mesh, 1, degree=1) +T = uw.discretisation.MeshVariable("T", mesh, 1, degree=3) + +# %% [markdown] +""" +## Stokes: rotated free-slip, and the body force as physics + +`add_rotated_freeslip_bc(0.0, boundary)` rotates each boundary node into its +own normal / tangent frame and constrains the normal component strongly. On a +circle that is exact to machine precision; a penalty or Nitsche condition +leaks at around 1e-3, which on a convection run shows up as spurious radial +flow at the boundary. + +The buoyancy is written as the force it is, in the units it has. The +nondimensionalisation turns it into `Ra T rhat` — that is where the Rayleigh +number printed above comes from, and writing it this way means the script +never has to be told what Ra is. +""" + +# %% +stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) +stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel +stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA +stokes.tolerance = 1.0e-8 +stokes.petsc_options.delValue("ksp_monitor") + +stokes.add_rotated_freeslip_bc(0.0, "Upper") +stokes.add_rotated_freeslip_bc(0.0, "Lower") + +radius = sympy.sqrt(mesh.X.dot(mesh.X)) +rhat = mesh.X / radius +# Name the coefficient rather than letting the product collapse into an +# anonymous number. Python multiplies the three quantities at assignment, so +# without this the run records a bare -0.97119 kg/(K m^2 s^2) and nothing +# saying where it came from. +BUOYANCY = uw.expression( + r"\rho_0 \alpha g", + RHO0 * ALPHA * GRAVITY, + "buoyancy coefficient: reference density x thermal expansivity x gravity", +) +stokes.bodyforce = -BUOYANCY * T.sym[0] * rhat + +# %% [markdown] +""" +## Transport + +`uw.systems.AdvDiffusion` composes an implicit Eulerian SUPG transport step +(Crank-Nicolson by default). Hot inner boundary, cold outer. +""" + +# %% +adv = uw.systems.AdvDiffusion(mesh, u_Field=T, V_fn=v.sym) +adv.constitutive_model = uw.constitutive_models.DiffusionModel +adv.constitutive_model.Parameters.diffusivity = KAPPA +adv.add_dirichlet_bc(1.0, "Lower") +adv.add_dirichlet_bc(0.0, "Upper") +adv.tolerance = 1.0e-8 +adv.petsc_options.delValue("ksp_monitor") + +# %% [markdown] +""" +## Initial condition + +A conductive profile with a mode-5 perturbation. Built from the nodal +coordinates in model units, which is what `beta`-style level sets and initial +conditions generally want: `T.coords` reads in metres, so divide by the length +scale once and work in the box's own units. +""" + +# %% +LENGTH_SCALE = float(model.get_fundamental_scales()["length"].to("m").magnitude) + + +def myr(q): + """A time quantity as a Myr string. `UWQuantity.__format__` delegates to + the bare float, so pint's `~` format specs do not apply to it.""" + return f"{float(q.to('Myr').magnitude):.4f} Myr" + +Xn = np.asarray(T.coords)[:, :2] / LENGTH_SCALE +rn = np.sqrt((Xn**2).sum(axis=1)) +thn = np.arctan2(Xn[:, 1], Xn[:, 0]) +shell = (rn - R_INNER) / (R_OUTER - R_INNER) + +T.array[:, 0, 0] = (1.0 - shell) + 0.1 * np.sin(5.0 * thn) * np.sin(np.pi * shell) + +adv.Unknowns.DuDt.initialise_history() +stokes.solve(zero_init_guess=True) + +# %% [markdown] +""" +## Diagnostics on the tracker + +`model.tracker.time`, `.step` and `.dt` are pre-seeded by convention. Anything +else you assign to it is captured by a snapshot and restored by a rewind, in +the same breath as the fields — which is exactly what a diagnostic wants, and +what a loose Python variable cannot give you. +""" + +# %% +v_rms_fn = sympy.sqrt(v.sym.dot(v.sym)) +area = float(uw.maths.Integral(mesh, sympy.sympify(1.0)).evaluate()) + + +def v_rms(): + return float(uw.maths.Integral(mesh, v_rms_fn).evaluate()) / area + + +# %% [markdown] +""" +## The loop + +Everything above is ordinary. This is the pattern: + +```python +while ...: + dt = adv.estimate_dt() + with model.step(dt, label="convect"): + adv.solve(timestep=dt, zero_init_guess=False) + stokes.solve(zero_init_guess=False) +``` + +`record_every = 1` asks each step to keep the state it started from — every +field, the transport history, the clock and the tracker diagnostics together — +captured before the operators run, which is the only correct point. +""" + +# %% +model.tracker.time = uw.quantity(0.0, "Myr") +model.tracker.step = 0 +model.tracker.dt = None +model.tracker.v_rms = v_rms() + +model.record_every = 1 +model.record_limit = params.uw_n_steps + +# The on-disk transcript needs no setting up: it is on by default, and this run +# will leave one under `transcripts/` beside a copy of this script. Section 5 +# shows what landed. To send it elsewhere, or to turn it off: +# +# model.transcript_file = "output/annulus.log" # somewhere else +# model.transcript_file = "output/annulus.jsonl" # JSON lines, for parsing +# model.transcript_file = None # off +say(f"transcript: {model.transcript_file or '(off)'}") + +for _ in range(int(params.uw_n_steps)): + dt = params.uw_dt_fraction * adv.estimate_dt() + + with model.step(dt, label="convect"): + adv.solve(timestep=dt, zero_init_guess=False) + stokes.solve(zero_init_guess=False) + model.tracker.v_rms = v_rms() + + say(f"step {model.tracker.step:>3d} " + f"t = {myr(model.tracker.time)} " + f"dt = {myr(dt)} " + f"v_rms = {model.tracker.v_rms:.4e}") + +# %% [markdown] +""" +## 1. What ran + +The transcript is an ordered account of each step: the interval it covered and +the operators it applied, named by what they solve. It answers "is this model +doing the thing the write-up says it does" without the script being +instrumented for it — which is the question you want to ask of someone else's +model, or of your own six months later. + +Note the `history_shift` between the two solves. That is the transport history +advancing — the thing that makes a step taken twice visible at all. +""" + +# %% +if params.uw_demos: + say("") + say("--- 1. the transcript " + "-" * 55) + for entry in model.transcript: + say(f" step {entry.index:>2d} dt = {myr(entry.dt):>12s} " + + " -> ".join(f"{e['kind']}:{e['name']}" for e in entry.events)) + say(f" {len(model.restore_points)} of {len(model.transcript)} steps " + f"are restorable") + +# %% [markdown] +""" +## 2. A rejected step + +A `model.step` block is a transaction. If it does not exit cleanly — an +exception, or a step abandoned because a diagnostic came out wrong — the clock +and the step counter are left exactly as they were, and nothing is added to the +transcript. Backstepping no longer has to remember to unwind a counter. + +The fields are yours to restore: take a snapshot before the block, and load it +in the handler. The clock never moved, so the two stay consistent. +""" + + +# %% +class StepRejected(Exception): + """Raised inside a step block to abandon it.""" + + +if params.uw_demos: + say("") + say("--- 2. a rejected step " + "-" * 51) + + before = (myr(model.tracker.time), model.tracker.step, len(model.transcript)) + snap = model.save_state() # BEFORE the step, not after + + reckless_dt = 50.0 * params.uw_dt_fraction * adv.estimate_dt() + try: + with model.step(reckless_dt, label="too big"): + adv.solve(timestep=reckless_dt, zero_init_guess=False) + stokes.solve(zero_init_guess=False) + if v_rms() > 4.0 * model.tracker.v_rms: + raise StepRejected("v_rms jumped; the step is not resolved") + except StepRejected as why: + model.load_state(snap) + say(f" rejected: {why}") + + after = (myr(model.tracker.time), model.tracker.step, len(model.transcript)) + say(f" clock/step/transcript before : {before}") + say(f" clock/step/transcript after : {after}") + say(f" unchanged: {before == after}") + +# %% [markdown] +""" +## 3. Playback + +`model.rewind()` puts the run back to the start of a completed step — fields, +transport history, clock and tracker diagnostics together — and truncates the +transcript to match, so it continues to describe the run that actually happened. + +Replaying the step from there reproduces it exactly. Re-*running* the script +does not: warm starts and preconditioner reuse are solver history rather than +model state, so two independent runs of the same problem diverge at the 1e-13 +level from the first step. If you need to look at a step twice, restore it +rather than re-run it. +""" + +# %% +if params.uw_demos: + say("") + say("--- 3. playback " + "-" * 58) + + T_end = np.asarray(T.array)[:, 0, 0].copy() + v_rms_end = model.tracker.v_rms + + target = model.rewind() + say(f" rewound to the start of step {target.index}: " + f"t = {myr(model.tracker.time)}, " + f"v_rms = {model.tracker.v_rms:.4e} " + f"(was {v_rms_end:.4e})") + + with model.step(target.dt, label="replay"): + adv.solve(timestep=target.dt, zero_init_guess=False) + stokes.solve(zero_init_guess=False) + model.tracker.v_rms = v_rms() + + T_replay = np.asarray(T.array)[:, 0, 0] + say(f" replayed: identical to the original step: " + f"{np.array_equal(T_replay, T_end)} " + f"max |dT| = {np.abs(T_replay - T_end).max():.3e}") + +# %% [markdown] +""" +## 4. Recording, not judging + +Call a solver twice inside one step — a predictor/corrector, a Picard iteration +on the coupled system, a retry — and its history advances twice, so the +physical step is taken twice. The timestep history and the solve counter look +identical to a single step, so the transcript is the only place it shows. + +The step records that and says nothing about it. Whether two shifts in one bar +are a mistake or legitimate sub-cycling is a reading of the transcript, made by +a later pass that can look across bars; inside the loop it would have to be +guessed. See `docs/developer/design/run-score-and-transcript.md`. + +What you see below is the bar's operator sequence with everything in it twice — +which is also why the figure gives that bar its own letter. +""" + +# %% +if params.uw_demos: + say("") + say("--- 4. a step taken twice " + "-" * 47) + + dt = params.uw_dt_fraction * adv.estimate_dt() + with model.step(dt, label="taken twice"): + adv.solve(timestep=dt, zero_init_guess=False) # a "predictor" + stokes.solve(zero_init_guess=False) + adv.solve(timestep=dt, zero_init_guess=False) # and a "corrector" + stokes.solve(zero_init_guess=False) + + entry = model.transcript[-1] + shifts = [e for e in entry.events if e["kind"] == "history_shift"] + say(f" history shifts in this one step: {len(shifts)}") + say(f" the step as recorded: {entry}") + +# %% [markdown] +""" +## 5. The log on disk + +The transcript lands on disk without being asked, in a stamped directory under +`transcripts/` beside a copy of the script that launched it and a `launch.json` +recording what invoked it. One line per step, flushed as it closes — so +`tail -f` follows a running job, and a run that is killed keeps everything up +to the moment it died. Nothing is created for a script that never takes a step. + +Three differences from `model.transcript`, all deliberate. An **abandoned** step +appears in the file and not in memory. A step aged out by `transcript_limit` +leaves memory but stays in the file. And a **backtrack** — `rewind()` or a +bare `load_state()` — writes its own line, because a log that shows step 7 and +then step 7 again, with nothing in between, is not a log of what happened. +""" + +# %% +if params.uw_demos: + say("") + say("--- 5. the log on disk " + "-" * 51) + say(f" {model.transcript_file}") + + run_dir = os.path.dirname(model.transcript_file) if model.transcript_file else "" + if run_dir: + say(f" the run directory holds: " + f"{', '.join(sorted(os.listdir(run_dir)))}") + say("") + with open(model.transcript_file, encoding="utf-8") as handle: + for line in handle.read().splitlines(): + say(" " + line) + +# %% [markdown] +""" +## 6. The same account, as a figure + +A terminal is not where a run belongs in a paper. `uw.transcript_diagram` renders +the record with **time running down the page** — one row per step, A4 portrait, +paginated — and writes it as a PDF, which opens anywhere, or an SVG if the +suffix says so. Both are written directly: no plotting library, no +rasterisation, no theme to fight with. `uw.transcript_flowchart` renders one +step's operator flow as Mermaid, for dropping into documentation. + +The layout decision worth knowing about: each distinct operator sequence gets a +**letter**, defined once at the foot of the figure. A column of `A` with a +single `B` in it says at a glance that one step did something different, where +a hundred spelled-out sequences say nothing and hide the one that matters. The +doubled step below is found that way rather than by reading. + +Both take a live model, which matters here because the text log is a report and +cannot be read back — pass a `.jsonl` log or the model itself. +""" + +# %% +if params.uw_demos: + say("") + say("--- 6. the figure " + "-" * 56) + + stem = model.transcript_file.rsplit(".", 1)[0] # beside the transcript + say(f" {uw.transcript_diagram(model, out=stem + '.pdf', title='Annulus convection - run log')}") + say(f" {uw.transcript_diagram(model, out=stem + '.svg', title='Annulus convection - run log')}") + say("") + for line in uw.transcript_flowchart(model).splitlines(): + say(" " + line) + +# %% +say("") +say(f"final: t = {myr(model.tracker.time)}, " + f"{model.tracker.step} steps, " + f"v_rms = {model.tracker.v_rms:.4e}") diff --git a/publications/blog-posts/figures/run-transcripts/make_figures.py b/publications/blog-posts/figures/run-transcripts/make_figures.py new file mode 100644 index 000000000..810580ac7 --- /dev/null +++ b/publications/blog-posts/figures/run-transcripts/make_figures.py @@ -0,0 +1,189 @@ +"""Figures for "A transcript for every run". + +Runs a short annulus convection model, deliberately rejects one step and +replays another, and renders the transcript it leaves as a figure. Run it from +this directory: + + python make_figures.py + +It writes, beside itself: + + run-transcript.svg the run as a figure (web) + run-transcript.pdf the same, for print + run-score.svg the same run as a score (web) + run-score.pdf the same, for print + transcript.log the text transcript the run wrote + transcript.jsonl the machine-readable transcript + score.txt the score rendered from the transcript + +Everything here is Underworld3's own machinery: the transcript is written +without being asked, and the figure is rendered from it afterwards. +""" + +import os + +import numpy as np +import sympy + +import underworld3 as uw + +params = uw.Params( + uw_cell_size=0.1, # mesh resolution, as a fraction of the outer radius + uw_n_steps=14, # timesteps before the demonstrations + uw_dt_fraction=0.5, # accuracy factor on estimate_dt() +) + +HERE = os.path.dirname(os.path.abspath(__file__)) + +# --- the model, in the units it is quoted in ------------------------------- +SHELL_THICKNESS = uw.quantity(2200, "km") +KAPPA = uw.quantity(1e-6, "m**2/s") +ETA = uw.quantity(1e22, "Pa*s") +DELTA_T = uw.quantity(2500, "K") +RHO0 = uw.quantity(3300, "kg/m**3") +ALPHA = uw.quantity(3e-5, "1/K") +GRAVITY = uw.quantity(9.81, "m/s**2") + +R_OUTER, R_INNER = 1.0, 0.55 + +uw.reset_default_model() +model = uw.get_default_model() +model.set_reference_quantities( + shell_thickness=SHELL_THICKNESS, + thermal_diffusivity=KAPPA, + mantle_viscosity=ETA, + temperature_contrast=DELTA_T, +) +# The transcript is on by default and lands in a stamped run directory; this +# script copies the two renderings out beside itself afterwards. +os.environ.setdefault("UW_TRANSCRIPT", os.path.join(HERE, "transcripts")) + +mesh = uw.meshing.Annulus( + radiusInner=R_INNER, radiusOuter=R_OUTER, + cellSize=params.uw_cell_size, degree=1, qdegree=3, +) +v = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=2) +p = uw.discretisation.MeshVariable("p", mesh, 1, degree=1) +T = uw.discretisation.MeshVariable("T", mesh, 1, degree=3) + +stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) +stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel +stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA +stokes.tolerance = 1.0e-8 +stokes.petsc_options.delValue("ksp_monitor") +# Rotated free-slip: exact on a circle, where a penalty condition leaks. +stokes.add_rotated_freeslip_bc(0.0, "Upper") +stokes.add_rotated_freeslip_bc(0.0, "Lower") + +radius = sympy.sqrt(mesh.X.dot(mesh.X)) +# The buoyancy as the force it is; the Rayleigh number falls out of the +# non-dimensionalisation rather than being typed in. +# Name the coefficient rather than letting the product collapse into an +# anonymous number. Python multiplies the three quantities at assignment, so +# without this the run records a bare -0.97119 kg/(K m^2 s^2) and nothing +# saying where it came from. +BUOYANCY = uw.expression( + r"\rho_0 \alpha g", + RHO0 * ALPHA * GRAVITY, + "buoyancy coefficient: reference density x thermal expansivity x gravity", +) +stokes.bodyforce = -BUOYANCY * T.sym[0] * mesh.X / radius + +adv = uw.systems.AdvDiffusion(mesh, u_Field=T, V_fn=v.sym) +adv.constitutive_model = uw.constitutive_models.DiffusionModel +adv.constitutive_model.Parameters.diffusivity = KAPPA +adv.add_dirichlet_bc(1.0, "Lower") +adv.add_dirichlet_bc(0.0, "Upper") +adv.tolerance = 1.0e-8 +adv.petsc_options.delValue("ksp_monitor") + +# Conductive profile with a mode-5 perturbation. +scale = float(model.get_fundamental_scales()["length"].to("m").magnitude) +X = np.asarray(T.coords)[:, :2] / scale +r = np.sqrt((X**2).sum(axis=1)) +th = np.arctan2(X[:, 1], X[:, 0]) +shell = (r - R_INNER) / (R_OUTER - R_INNER) +T.array[:, 0, 0] = (1.0 - shell) + 0.1 * np.sin(5.0 * th) * np.sin(np.pi * shell) +adv.Unknowns.DuDt.initialise_history() +stokes.solve(zero_init_guess=True) + +v_rms_fn = sympy.sqrt(v.sym.dot(v.sym)) +area = float(uw.maths.Integral(mesh, sympy.sympify(1.0)).evaluate()) + + +def v_rms(): + return float(uw.maths.Integral(mesh, v_rms_fn).evaluate()) / area + + +# --- the run --------------------------------------------------------------- +model.tracker.time = uw.quantity(0.0, "Myr") +model.tracker.step = 0 +model.tracker.v_rms = v_rms() +model.record_every = 1 +model.record_limit = int(params.uw_n_steps) + 4 + +for _ in range(int(params.uw_n_steps)): + dt = params.uw_dt_fraction * adv.estimate_dt() + with model.step(dt, label="convect"): + adv.solve(timestep=dt, zero_init_guess=False) + stokes.solve(zero_init_guess=False) + model.tracker.v_rms = v_rms() + + +class StepRejected(Exception): + """Raised inside a step block to abandon it.""" + + +# A step fifty times too large, rejected on a diagnostic after it ran. +snapshot = model.save_state() +reckless = 50.0 * params.uw_dt_fraction * adv.estimate_dt() +try: + with model.step(reckless, label="too big"): + adv.solve(timestep=reckless, zero_init_guess=False) + stokes.solve(zero_init_guess=False) + if v_rms() > 4.0 * model.tracker.v_rms: + raise StepRejected("v_rms jumped") +except StepRejected: + model.load_state(snapshot) + +# Back one step, and take it again: a replay reproduces it exactly. +target = model.rewind() +with model.step(target.dt, label="replay"): + adv.solve(timestep=target.dt, zero_init_guess=False) + stokes.solve(zero_init_guess=False) + model.tracker.v_rms = v_rms() + +# One step with the transport solved twice — a predictor/corrector written +# without noticing that the history advances on every solve. +dt = params.uw_dt_fraction * adv.estimate_dt() +with model.step(dt, label="taken twice"): + adv.solve(timestep=dt, zero_init_guess=False) + stokes.solve(zero_init_guess=False) + adv.solve(timestep=dt, zero_init_guess=False) + stokes.solve(zero_init_guess=False) + +# --- the renderings -------------------------------------------------------- +# From the RECORD on disk rather than from the live model: the in-memory +# transcript holds what the run can still undo, so the abandoned step and the +# backtracks — the rows worth looking at — are only in the file. +import shutil + +run_dir = os.path.dirname(model.transcript_file) +record = os.path.join(run_dir, "transcript.jsonl") +for name in ("transcript.log", "transcript.jsonl"): + shutil.copyfile(os.path.join(run_dir, name), os.path.join(HERE, name)) + +title = "Annulus convection — run transcript" +uw.transcript_diagram(record, out=os.path.join(HERE, "run-transcript.svg"), title=title) +uw.transcript_diagram(record, out=os.path.join(HERE, "run-transcript.pdf"), title=title) + +score_title = "Annulus convection — score" +uw.transcript_score_figure(record, out=os.path.join(HERE, "run-score.svg"), title=score_title) +uw.transcript_score_figure(record, out=os.path.join(HERE, "run-score.pdf"), title=score_title) + +score = uw.transcript_score(record, width=15) +with open(os.path.join(HERE, "score.txt"), "w", encoding="utf-8") as handle: + handle.write(score + "\n") + +uw.pprint(score, clean_display=False) +uw.pprint(f"figures and transcript written to {HERE}", clean_display=False) diff --git a/publications/blog-posts/figures/run-transcripts/run-score.pdf b/publications/blog-posts/figures/run-transcripts/run-score.pdf new file mode 100644 index 000000000..3102b1522 Binary files /dev/null and b/publications/blog-posts/figures/run-transcripts/run-score.pdf differ diff --git a/publications/blog-posts/figures/run-transcripts/run-score.svg b/publications/blog-posts/figures/run-transcripts/run-score.svg new file mode 100644 index 000000000..f7f5ad03d --- /dev/null +++ b/publications/blog-posts/figures/run-transcripts/run-score.svg @@ -0,0 +1,103 @@ + + +Annulus convection — score +started 2026-09-12T17:25:37+00:00 · no terminator: still running, or interrupted · 3 parts +step +t/Myr +dt/Myr +AdvectionDiffusion(T) +Stokes(v) +EulerianSUPG(T) + +0 +0.175907 +0.175907 + +1 + +3 + +2 + +1 +12 +×12 +0.5015 +12.56 + + +0.3256 +2.704 + + + + + + + + +13 +16.8708 +4.31057 + +1 + +3 + +2 +14 +446.031 +429.161 + +1 + +3 + +2 +abandoned +13 +16.8708 +4.31057 + +1 + +3 + +2 +14 +25.454 +8.58321 + +1 + +4 + +3 + +6 + +2 + +5 + + + + + + + +rewind 1 (+1) + + + +again + +ran; the digit is the order it ran within the step + +ran in a step that was then abandoned + +did nothing in this step + + +the steps between did exactly this, unchanged; first and last values are shown where they differ + diff --git a/publications/blog-posts/figures/run-transcripts/run-transcript.pdf b/publications/blog-posts/figures/run-transcripts/run-transcript.pdf new file mode 100644 index 000000000..e759b370a Binary files /dev/null and b/publications/blog-posts/figures/run-transcripts/run-transcript.pdf differ diff --git a/publications/blog-posts/figures/run-transcripts/run-transcript.svg b/publications/blog-posts/figures/run-transcripts/run-transcript.svg new file mode 100644 index 000000000..d1800a320 --- /dev/null +++ b/publications/blog-posts/figures/run-transcripts/run-transcript.svg @@ -0,0 +1,133 @@ + + +Annulus convection — run transcript +started 2026-09-12T17:25:37+00:00 · t = 0 to 25.45 Myr · 17 steps · 1 abandoned · 2 backtrack(s) +scales: length 2.2e+06 m time 4.84e+18 s mass 1.065e+47 kg temperature 2500 K +step +t/Myr +dt/Myr +seq +dt (log scale) +wall + +0 +0.1759 +0.1759 +A + + +1 +0.5015 +0.3256 +A + + +2 +0.9909 +0.4894 +A + + +3 +1.515 +0.5237 +A + + +4 +2.096 +0.581 +A + + +5 +2.746 +0.6503 +A + + +6 +3.478 +0.7323 +A + + +7 +4.317 +0.8383 +A + + +8 +5.296 +0.9797 +A + + +9 +6.472 +1.176 +A + + +10 +7.936 +1.464 +A + + +11 +9.856 +1.92 +A + + +12 +12.56 +2.704 +A + + +13 +16.87 +4.311 +A + + +14 +446 +429.2 +A + +abandoned + +13 +16.87 +4.311 +A + + +14 +25.45 +8.583 +B + + + + + + +rewind 1 (+1) + + + +again +Operator sequences — what the seq letter on each row stands for +A +16 steps +AdvectionDiffusion(T) > shift EulerianSUPG(T) > Stokes(v) +B +1 step +AdvectionDiffusion(T) > shift EulerianSUPG(T) > Stokes(v) > AdvectionDiffusion(T) > +shift EulerianSUPG(T) > Stokes(v) + diff --git a/publications/blog-posts/figures/run-transcripts/score.txt b/publications/blog-posts/figures/run-transcripts/score.txt new file mode 100644 index 000000000..d3a1cca0f --- /dev/null +++ b/publications/blog-posts/figures/run-transcripts/score.txt @@ -0,0 +1,18 @@ +score · model 'default' +started 2026-09-12T17:25:37+00:00 +no terminator: this run is still going, or it was interrupted. What follows is the score of a prefix. + + step t/Myr dt/Myr │ AdvectionDiffus │ Stokes(v) │ EulerianSUPG(T) │ +─────────────────────────────────────────────────────────────────────────────────── + 0 0.175907 0.175907 │ 1 │ 3 │ 2 │ + …13 16.8708 4.31057 │ ↓ │ ↓ │ ↓ │ ×13 unchanged, dt 0.1759 → 4.311 + 14 446.031 429.161 │ 1 │ 3 │ 2 │ ABANDONED + │ restore from a snapshot; the clock now reads 16.8708 Myr + │ rewind to the start of step 13 (t = 12.5602 Myr); 1 step(s) undone + 13 16.8708 4.31057 │ 1 │ 3 │ 2 │ + 14 25.454 8.58321 │ 1,4 │ 3,6 │ 2,5 │ + +17 step(s), 3 part(s): AdvectionDiffusion(T), Stokes(v), EulerianSUPG(T) +↓ the steps between did exactly this, unchanged +· this part did nothing in that step +digits are the order the parts ran within the step diff --git a/publications/blog-posts/figures/run-transcripts/transcript.jsonl b/publications/blog-posts/figures/run-transcripts/transcript.jsonl new file mode 100644 index 000000000..7b92458ab --- /dev/null +++ b/publications/blog-posts/figures/run-transcripts/transcript.jsonl @@ -0,0 +1,21 @@ +{"kind": "run", "model": "default", "started": "2026-09-12T17:25:37+00:00", "scales": {"length": {"magnitude": 2200000.00000002, "units": "meter"}, "time": {"magnitude": 4.839999999999929e+18, "units": "second"}, "mass": {"magnitude": 1.0648000000000114e+47, "units": "kilogram"}, "temperature": {"magnitude": 2500.000000000002, "units": "kelvin"}}} +{"kind": "part", "part": "SNES_Stokes#10", "label": "SNES_Stokes(v)", "at_step": 0, "fingerprint": "Matrix([[N.x*\\rho_0 \\alpha g*{T}(N.x, N.y)/sqrt(N.x**2 + N.y**2)], [N.y*\\rho_0 \\alpha g*{T}(N.x, N.y)/sqrt(N.x**2 + N.y**2)]])Matrix([[\\eta*\\uplambda*({v}_{ 0,0}(N.x, N.y) + {v}_{ 1,1}(N.x, N.y)) + 2*\\eta*{v}_{ 0,0}(N.x, N.y) - {p}(N.x, N.y), 2*\\eta*({v}_{ 0,1}(N.x, N.y)/2 + {v}_{ 1,0}(N.x, N.y)/2)], [2*\\eta*({v}_{ 0,1}(N.x, N.y)/2 + {v}_{ 1,0}(N.x, N.y)/2), \\eta*\\uplambda*({v}_{ 0,0}(N.x, N.y) + {v}_{ 1,1}(N.x, N.y)) + 2*\\eta*{v}_{ 1,1}(N.x, N.y) - {p}(N.x, N.y)]])Matrix([[{v}_{ 0,0}(N.x, N.y) + {v}_{ 1,1}(N.x, N.y)]])", "solver": "SNES_Stokes", "unknown": "v", "dim": 2, "cdim": 2, "forms": {"F0": {"symbol": "\\mathbf{f}_0\\left( \\mathbf{u} \\right)", "description": "Velocity equation body force term (pointwise).", "latex": "\\left[\\begin{matrix}\\frac{\\mathrm{x} \\rho_0 \\alpha g {T}(\\mathbf{x})}{\\sqrt{\\mathrm{x}^{2} + \\mathrm{y}^{2}}}\\\\\\frac{\\mathrm{y} \\rho_0 \\alpha g {T}(\\mathbf{x})}{\\sqrt{\\mathrm{x}^{2} + \\mathrm{y}^{2}}}\\end{matrix}\\right]", "text": "Matrix([[N.x*\\rho_0 \\alpha g*{T}(N.x, N.y)/sqrt(N.x**2 + N.y**2)], [N.y*\\rho_0 \\alpha g*{T}(N.x, N.y)/sqrt(N.x**2 + N.y**2)]])", "where": [{"symbol": "\\rho_0 \\alpha g", "latex": "\\mathtt{\\text{0.9711900000000001 [kilogram / kelvin / meter ** 2 / second ** 2]}}", "value": "0.9711900000000001 [kilogram / kelvin / meter ** 2 / second ** 2]", "units": "kilogram / kelvin / meter ** 2 / second ** 2", "description": "buoyancy coefficient: reference density x thermal expansivity x gravity", "where": []}]}, "F1": {"symbol": "\\mathbf{F}_1\\left( \\mathbf{u} \\right)", "description": "Velocity equation flux/stress term (pointwise).", "latex": "\\left[\\begin{matrix}\\eta \\uplambda \\left({v}_{ 0,0}(\\mathbf{x}) + {v}_{ 1,1}(\\mathbf{x})\\right) + 2 \\eta {v}_{ 0,0}(\\mathbf{x}) - {p}(\\mathbf{x}) & 2 \\eta \\left(\\frac{{v}_{ 0,1}(\\mathbf{x})}{2} + \\frac{{v}_{ 1,0}(\\mathbf{x})}{2}\\right)\\\\2 \\eta \\left(\\frac{{v}_{ 0,1}(\\mathbf{x})}{2} + \\frac{{v}_{ 1,0}(\\mathbf{x})}{2}\\right) & \\eta \\uplambda \\left({v}_{ 0,0}(\\mathbf{x}) + {v}_{ 1,1}(\\mathbf{x})\\right) + 2 \\eta {v}_{ 1,1}(\\mathbf{x}) - {p}(\\mathbf{x})\\end{matrix}\\right]", "text": "Matrix([[\\eta*\\uplambda*({v}_{ 0,0}(N.x, N.y) + {v}_{ 1,1}(N.x, N.y)) + 2*\\eta*{v}_{ 0,0}(N.x, N.y) - {p}(N.x, N.y), 2*\\eta*({v}_{ 0,1}(N.x, N.y)/2 + {v}_{ 1,0}(N.x, N.y)/2)], [2*\\eta*({v}_{ 0,1}(N.x, N.y)/2 + {v}_{ 1,0}(N.x, N.y)/2), \\eta*\\uplambda*({v}_{ 0,0}(N.x, N.y) + {v}_{ 1,1}(N.x, N.y)) + 2*\\eta*{v}_{ 1,1}(N.x, N.y) - {p}(N.x, N.y)]])", "where": [{"symbol": "\\eta", "latex": "\\mathtt{\\text{1e+22 [pascal * second]}}", "value": "1e+22 [pascal * second]", "units": "pascal * second", "description": "Shear viscosity", "where": []}, {"symbol": "\\uplambda", "latex": "0", "value": "0", "units": null, "description": "Numerical Penalty", "where": []}]}, "PF0": {"symbol": "\\mathbf{h}_0\\left( \\mathbf{p} \\right)", "description": "Pressure equation constraint term (continuity).", "latex": "\\left[\\begin{matrix}{v}_{ 0,0}(\\mathbf{x}) + {v}_{ 1,1}(\\mathbf{x})\\end{matrix}\\right]", "text": "Matrix([[{v}_{ 0,0}(N.x, N.y) + {v}_{ 1,1}(N.x, N.y)]])", "where": []}}, "boundary_conditions": [{"mechanism": "rotated_freeslip", "type": "rotated free-slip", "boundary": "Upper", "latex": "\\mathbf{u}\\cdot\\hat{\\mathbf{n}} = 0", "text": "u . n = 0", "normal": "mesh"}, {"mechanism": "rotated_freeslip", "type": "rotated free-slip", "boundary": "Lower", "latex": "\\mathbf{u}\\cdot\\hat{\\mathbf{n}} = 0", "text": "u . n = 0", "normal": "mesh"}], "terms": [{"name": "bodyforce", "description": "body force per unit volume; F0 is its negative", "latex": "\\left[\\begin{matrix}- \\frac{\\mathrm{x} \\rho_0 \\alpha g {T}(\\mathbf{x})}{\\sqrt{\\mathrm{x}^{2} + \\mathrm{y}^{2}}}\\\\- \\frac{\\mathrm{y} \\rho_0 \\alpha g {T}(\\mathbf{x})}{\\sqrt{\\mathrm{x}^{2} + \\mathrm{y}^{2}}}\\end{matrix}\\right]", "text": "Matrix([[-N.x*\\rho_0 \\alpha g*{T}(N.x, N.y)/sqrt(N.x**2 + N.y**2)], [-N.y*\\rho_0 \\alpha g*{T}(N.x, N.y)/sqrt(N.x**2 + N.y**2)]])", "where": [{"symbol": "\\rho_0 \\alpha g", "latex": "\\mathtt{\\text{0.9711900000000001 [kilogram / kelvin / meter ** 2 / second ** 2]}}", "value": "0.9711900000000001 [kilogram / kelvin / meter ** 2 / second ** 2]", "units": "kilogram / kelvin / meter ** 2 / second ** 2", "description": "buoyancy coefficient: reference density x thermal expansivity x gravity", "where": []}]}, {"name": "penalty", "description": "augmented-Lagrangian grad-div penalty (0 = off)", "latex": "\\uplambda", "text": "0", "where": []}, {"name": "ViscousFlowModel.shear_viscosity_0", "description": "constitutive parameter", "latex": "\\mathtt{\\text{1e+22 [pascal * second]}}", "text": "1e+22 [pascal * second]", "where": []}], "terms_declared": true} +{"kind": "step", "index": 0, "label": "convect", "t0": {"magnitude": 0.0, "units": "megayear"}, "t1": {"magnitude": 0.1759072436791037, "units": "megayear"}, "dt": {"magnitude": 5551210433127.684, "units": "second"}, "completed": true, "restorable": true, "wall": 0.2239170828834176, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 1.1469443043652406e-06, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "step", "index": 1, "label": "convect", "t0": {"magnitude": 0.1759072436791037, "units": "megayear"}, "t1": {"magnitude": 0.5015458223417135, "units": "megayear"}, "dt": {"magnitude": 10276372010003.176, "units": "second"}, "completed": true, "restorable": true, "wall": 0.06961754243820906, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 2.123217357438704e-06, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "step", "index": 2, "label": "convect", "t0": {"magnitude": 0.5015458223417135, "units": "megayear"}, "t1": {"magnitude": 0.9909387375900567, "units": "megayear"}, "dt": {"magnitude": 15444065862241.113, "units": "second"}, "completed": true, "restorable": true, "wall": 0.07091937493532896, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 3.190922698810194e-06, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "step", "index": 3, "label": "convect", "t0": {"magnitude": 0.9909387375900567, "units": "megayear"}, "t1": {"magnitude": 1.5145938870536368, "units": "megayear"}, "dt": {"magnitude": 16525299744711.88, "units": "second"}, "completed": true, "restorable": true, "wall": 0.06876895809546113, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 3.4143181290727526e-06, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "step", "index": 4, "label": "convect", "t0": {"magnitude": 1.5145938870536368, "units": "megayear"}, "t1": {"magnitude": 2.0956275212136184, "units": "megayear"}, "dt": {"magnitude": 18336027013367.04, "units": "second"}, "completed": true, "restorable": true, "wall": 0.06908804224804044, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 3.788435333340353e-06, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "step", "index": 5, "label": "convect", "t0": {"magnitude": 2.0956275212136184, "units": "megayear"}, "t1": {"magnitude": 2.7459250179997823, "units": "megayear"}, "dt": {"magnitude": 20521828284579.05, "units": "second"}, "completed": true, "restorable": true, "wall": 0.06877437513321638, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 4.240047166235403e-06, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "step", "index": 6, "label": "convect", "t0": {"magnitude": 2.7459250179997823, "units": "megayear"}, "t1": {"magnitude": 3.4782572379243297, "units": "megayear"}, "dt": {"magnitude": 23110647263490.9, "units": "second"}, "completed": true, "restorable": true, "wall": 0.06792695820331573, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 4.774927120556041e-06, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "step", "index": 7, "label": "convect", "t0": {"magnitude": 3.4782572379243297, "units": "megayear"}, "t1": {"magnitude": 4.316530568513284, "units": "megayear"}, "dt": {"magnitude": 26453894457393.99, "units": "second"}, "completed": true, "restorable": true, "wall": 0.06781033333390951, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 5.4656806730153666e-06, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "step", "index": 8, "label": "convect", "t0": {"magnitude": 4.316530568513284, "units": "megayear"}, "t1": {"magnitude": 5.296186199263291, "units": "megayear"}, "dt": {"magnitude": 30915580532956.42, "units": "second"}, "completed": true, "restorable": true, "wall": 0.06854437477886677, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 6.387516639040676e-06, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "step", "index": 9, "label": "convect", "t0": {"magnitude": 5.296186199263291, "units": "megayear"}, "t1": {"magnitude": 6.472194257077156, "units": "megayear"}, "dt": {"magnitude": 37111991885266.82, "units": "second"}, "completed": true, "restorable": true, "wall": 0.06765404203906655, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 7.66776691844367e-06, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "step", "index": 10, "label": "convect", "t0": {"magnitude": 6.472194257077156, "units": "megayear"}, "t1": {"magnitude": 7.936157563583196, "units": "megayear"}, "dt": {"magnitude": 46199168441395.016, "units": "second"}, "completed": true, "restorable": true, "wall": 0.06849033338949084, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 9.545282735825555e-06, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "step", "index": 11, "label": "convect", "t0": {"magnitude": 7.936157563583196, "units": "megayear"}, "t1": {"magnitude": 9.856349768374267, "units": "megayear"}, "dt": {"magnitude": 60596657521914.71, "units": "second"}, "completed": true, "restorable": true, "wall": 0.06813162518665195, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 1.2519970562379255e-05, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "step", "index": 12, "label": "convect", "t0": {"magnitude": 9.856349768374267, "units": "megayear"}, "t1": {"magnitude": 12.560193831583604, "units": "megayear"}, "dt": {"magnitude": 85326829409134.98, "units": "second"}, "completed": true, "restorable": true, "wall": 0.06861650012433529, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 1.7629510208499222e-05, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "step", "index": 13, "label": "convect", "t0": {"magnitude": 12.560193831583604, "units": "megayear"}, "t1": {"magnitude": 16.870760189578128, "units": "megayear"}, "dt": {"magnitude": 136031128899048.02, "units": "second"}, "completed": true, "restorable": true, "wall": 0.06833104230463505, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 2.8105605144431817e-05, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "step", "index": 14, "label": "too big", "t0": {"magnitude": 16.870760189578128, "units": "megayear"}, "t1": {"magnitude": 446.03133916692144, "units": "megayear"}, "dt": {"magnitude": 1.354327788713541e+16, "units": "second"}, "completed": false, "restorable": false, "wall": 0.0689316252246499, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 0.0027981979105652084, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "restore", "message": "restore from a snapshot; the clock now reads 16.8708 Myr", "source": "memory", "t": {"magnitude": 16.870760189578128, "units": "megayear"}} +{"kind": "rewind", "message": "rewind to the start of step 13 (t = 12.5602 Myr); 1 step(s) undone", "to_step": 13, "steps_undone": 1, "t": {"magnitude": 12.560193831583604, "units": "megayear"}} +{"kind": "step", "index": 13, "label": "replay", "t0": {"magnitude": 12.560193831583604, "units": "megayear"}, "t1": {"magnitude": 16.870760189578128, "units": "megayear"}, "dt": {"magnitude": 136031128899048.02, "units": "second"}, "completed": true, "restorable": true, "wall": 0.2741635418497026, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 2.8105605144431817e-05, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} +{"kind": "step", "index": 14, "label": "taken twice", "t0": {"magnitude": 16.870760189578128, "units": "megayear"}, "t1": {"magnitude": 25.453971769124998, "units": "megayear"}, "dt": {"magnitude": 270865557742708.22, "units": "second"}, "completed": true, "restorable": true, "wall": 0.1339148748666048, "events": [{"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 5.5963958211304166e-05, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}, {"kind": "solve", "name": "SNES_AdvectionDiffusion_Composed(T)", "part": "SNES_AdvectionDiffusion_Composed#17"}, {"kind": "history_shift", "name": "EulerianSUPG(T)", "dt": 5.5963958211304166e-05, "part": "EulerianSUPG#18"}, {"kind": "solve", "name": "SNES_Stokes(v)", "part": "SNES_Stokes#10"}]} diff --git a/publications/blog-posts/run-transcripts.md b/publications/blog-posts/run-transcripts.md new file mode 100644 index 000000000..9c8903f4f --- /dev/null +++ b/publications/blog-posts/run-transcripts.md @@ -0,0 +1,219 @@ +--- +title: "A Transcript for Every Run" +status: draft +feeds_into: [paper-2, release-post] +target: underworldcode.org (Ghost) +figures: figures/run-transcripts +--- + +Six months after a result, the questions are always the same: which version of +the code produced it, what the parameters were, and whether the run did what +the caption says it did. Underworld 1 could answer all three from a single +file. Underworld3 could answer none of them, until recently, and the reason is +worth setting out, because it is a consequence of what makes the current code +usable at all. + +## The Underworld 1 answer + +A model in Underworld 1 was an XML document. The launcher read it, StGermain +assembled the components it named, and the run followed. Provenance came for +free: the XML *was* the model, so keeping it kept everything. Two people with +the same XML and the same binary were running the same experiment, and could +say so without qualification. + +The cost was that the XML was a programming language without being one. It had +no expressions worth the name, no control flow a reader could follow, and no +debugger. Anything the schema had not anticipated could not be said, and a +model that needed something new needed C. A great deal of scientific intent +ended up encoded as combinations of components that happened to compose, which +is a poor medium for explaining what a model does. + +## Why that answer is not available now + +A model in Underworld3 is a Python program. The library takes SymPy expressions +for constitutive behaviour, boundary conditions and source terms; it puts no +constraint on what happens between solves; and the timestep loop is written by +the user. That openness is deliberate, and it is the same property that makes +the library drivable by a colleague who has never used it, or by a language +model: the surface is compositional, so a reader can predict what a call will +do from what the pieces mean. + +It also means no serialisable model document exists, even in principle. The +model is an arbitrary program, and what it does is a sequence of calls — +including calls whose arguments depend on runtime state. Nothing can be +recovered from the objects afterwards, because the objects do not record the +order in which anyone touched them. + +That is a real trade. Underworld 1 had a complete statement of intent and very +little expressive power. Underworld3 has the expressive power and, until +recently, kept no statement of anything at all. + +## What a run leaves now + +Every Underworld3 run that takes a timestep writes a **transcript**: a record +of what it did, appended and flushed as each step closes. No configuration is +involved. A run that takes no step writes nothing, and a run under `pytest` +writes nothing, so the default costs nothing where it would only be noise. + +``` +transcripts/2026-09-12T03-31-34-make_figures/ + make_figures.py the script that launched it, verbatim + launch.json argv, interpreter, cwd, versions, commit + transcript.log one aligned line per step, flushed + transcript.jsonl the same record, machine-readable +``` + +The directory is stamped with the time the run started, so a second run does +not overwrite the first. The path is printed when the run begins and again +when it ends, and a `latest` symlink points at the most recent one. + +`launch.json` is the replacement for the XML, and it is a weaker thing +honestly labelled. It cannot state what the model *is*; it states what was +*run*: + +```json +{ + "started": "2026-09-12T03:31:34+00:00", + "argv": ["make_figures.py"], + "python": "3.12.12", + "underworld3": "0.0.0", + "mpi_size": 1, + "git_commit": "638a6c3c5a3859fc40e2a87f4f918ba9f3aebcb1", + "git_dirty": true, + "script": "make_figures.py" +} +``` + +The entry script is copied beside it. Modules it imports are not, which is +what the commit id is there to cover, and the file says so rather than leaving +a reader to find out. `git_dirty` is the field that earns its place: a commit +id with uncommitted changes on top of it identifies nothing, and recording the +flag is the difference between provenance and the appearance of it. + +## Where this lands against FAIR + +**Findable** is the stamped directory, the `latest` pointer and the announced +path. A run that produced a figure can be located without anyone having +remembered to write down where it went. + +**Accessible** is the two file formats. The text transcript is read with +`tail -f` while a job is running; the JSON Lines record is read with `jq`, or +with `uw.read_transcript`. Neither needs Underworld3 installed to open. + +**Interoperable** is the way dimensional values are stored. Each carries its +own magnitude and unit string — `{"magnitude": 4.31, "units": "megayear"}` — +rather than a bespoke convention that a reader has to be told about. + +**Reusable** is the launch record together with the transcript. The two +together let someone decide what to change, which is what reuse requires. + +## Provenance is not reproducibility + +The transcript identifies a run. It does not promise that running the script +again produces the same numbers, and on this code it will not. Warm starts and +preconditioner reuse are solver history that sits outside model state, so two +independent runs of the same script on the same machine diverge at the 1e-13 +level from the first step. + +Restoring is exact where re-running is not. A step restored from the snapshot +it began with and taken again reproduces its temperature field bit for bit — +`max |ΔT| = 0` across every step of the annulus run below. That asymmetry is +the practical reason a run keeps its own restore points: a step that +misbehaved can be looked at twice, which re-running cannot give you. + +## Reading a run + +The transcript is complete, which creates the problem a profiler has: every +step is in it, and most steps are identical. Removing the repetitive ones by +hand would be a judgement about what mattered. A **score** removes them by +rule instead — it groups consecutive steps that did exactly the same thing, +and carries the first and last value of anything that changed across the +group, so a timestep that grew by a factor of eight survives the grouping. + +``` +score · model 'default' +started 2026-09-12T15:22:22+00:00 +no terminator: this run is still going, or it was interrupted. What follows is the score of a prefix. + + step t/Myr dt/Myr │ AdvectionDiffus │ Stokes(v) │ EulerianSUPG(T) │ +─────────────────────────────────────────────────────────────────────────────────── + 0 0.175907 0.175907 │ 1 │ 3 │ 2 │ + …13 16.8708 4.31057 │ ↓ │ ↓ │ ↓ │ ×13 unchanged, dt 0.1759 → 4.311 + 14 446.031 429.161 │ 1 │ 3 │ 2 │ ABANDONED + │ restore from a snapshot; the clock now reads 16.8708 Myr + │ rewind to the start of step 13 (t = 12.5602 Myr); 1 step(s) undone + 13 16.8708 4.31057 │ 1 │ 3 │ 2 │ + 14 25.454 8.58321 │ 1,4 │ 3,6 │ 2,5 │ + +17 step(s), 3 part(s): AdvectionDiffusion(T), Stokes(v), EulerianSUPG(T) +↓ the steps between did exactly this, unchanged +· this part did nothing in that step +digits are the order the parts ran within the step +``` + +Each column is a participant: the advection-diffusion solver, the Stokes +solver, and the transport history the first of them holds. The digits give the +order they ran within the step. The grouped line asserts that steps 1 to 13 +were *identical* to step 0 in what ran and in what order, so nothing can hide +behind it. Seventeen steps become six rows, and the rows that remain are the +ones a reader would have picked out. Grouping is optional, because a run whose +timestep is itself the thing under examination is easier to read one row at a +time. + +Three things are visible in those six rows, and in @fig-score, that no print +statement was written to report. Step 14 attempted 429 Myr and was rejected on a velocity +diagnostic, so the run went back to step 13 and took it again. The replayed +step 13 took 0.30 s of wall clock against the original's 0.06, because the +restore discarded the warm start. And the last step ran everything twice — +`1,4 │ 3,6 │ 2,5` — which is a predictor-corrector written without noticing +that the transport history advances on every solve, so the temperature +advanced two intervals while the clock advanced one. + +```{figure} figures/run-transcripts/run-score.svg +:label: fig-score +:alt: A score of a 17-step annulus convection run. Three columns — AdvectionDiffusion(T), Stokes(v) and EulerianSUPG(T) — each carry a filled mark per step with the order it ran: 1, 3 and 2. Steps 1 to 12 group into one shaded band with a downward arrow in each column, labelled x12, showing t running 0.5015 to 12.56 Myr and dt 0.3256 to 2.704 Myr. Step 14 at 446.031 Myr is drawn with hollow dashed marks and labelled abandoned. A dashed red arrow in the left gutter labelled "rewind 1 (+1)" runs back from it to step 13, and a blue arrow labelled "again" runs down to the replayed step 13. The final step carries two marks in every column, numbered 1 and 4, 3 and 6, 2 and 5. + +The same run as a score. Each column is a participant and each row a step; a +filled mark carries the order that participant ran within the step. Steps 1 to +12 are grouped into one band, which asserts that they did exactly what step 0 +did, and reports the range of `t` and `dt` across them. The final step carries +two marks in every column, which is the predictor-corrector advancing the +transport history twice. +``` + +```{figure} figures/run-transcripts/run-transcript.svg +:label: fig-run-transcript +:alt: A portrait figure of a 17-step annulus convection run. Steps run down the page with dt as a horizontal bar; dt grows from 0.18 to 4.31 Myr over the first fourteen. Step 14 is drawn hatched in red at 429 Myr and marked abandoned. A dashed arrow in the left gutter runs back from it to step 13, and a solid arrow labelled "again" runs down to the replayed step 13. The final step carries the letter B where every other step carries A. + +The same run as a figure. Time runs down the page, `dt` across it, and the two +backtracks are drawn in the left gutter: back out of the abandoned step, then +down to the step that was taken again. Each distinct operator sequence gets a +letter, defined once at the foot, so the one step that did something different +is the one letter that differs. +``` + +The figure and the score are rendered from the JSON record after the run, or +during it — the score above is of a run still in progress, which is why it +says so. Both are produced by +[`make_figures.py`](figures/run-transcripts/make_figures.py), which is also +the run they describe. + +## What it costs, and what it does not do + +A step appends two lines, one to each format, and flushes. The snapshots that +make `rewind` possible cost about 13 bytes per primary degree of freedom per +step and are off unless asked for; the transcript itself is bytes. + +Three things are not in it yet. A part that does nothing in a step does not +appear, so silence and absence are indistinguishable. Mesh deformation and +adaptation are not recorded as events, so a history stored before the mesh +moved and read after it has nothing in the record to flag it. And the score is +inferred from what repeated rather than declared by the script, so it can say +that one step differs from its neighbours and cannot yet say that a run +disagrees with its own description. + +The last of those is where the Underworld 1 comparison ends up. A declared +score would be the XML's descendant — a statement of what a step is supposed +to contain — checked against the transcript rather than executed from it. The +document would describe the model without having to be the only way to express +it. diff --git a/scripts/test.sh b/scripts/test.sh index f58ca36f9..eff3eebc6 100755 --- a/scripts/test.sh +++ b/scripts/test.sh @@ -134,6 +134,7 @@ if [ $PARALLEL_ONLY -eq 0 ]; then # unbatched file would have closed the issue without closing the gap. # level_2/tier_b, ~55s serial; passes at np=1 and np=2. $PYTEST tests/test_1072_free_surface_spherical.py || status=1 + $PYTEST tests/test_1074_free_surface_config_drift.py || status=1 # Diffusion / Advection tests $PYTEST tests/test_1100*py || status=1 diff --git a/src/underworld3/__init__.py b/src/underworld3/__init__.py index 4ed69e815..50bfd7c5a 100644 --- a/src/underworld3/__init__.py +++ b/src/underworld3/__init__.py @@ -218,9 +218,16 @@ def view(): create_model, get_default_model, reset_default_model, + read_transcript, ThermalConvectionConfig, create_thermal_convection_model, ) +from .utilities.transcript_report import ( + transcript_diagram, + transcript_flowchart, + transcript_score, + transcript_score_figure, +) from .parameters import ParameterRegistry, ParameterType from .materials import MaterialRegistry, MaterialProperty from .constitutive_models import MultiMaterialConstitutiveModel diff --git a/src/underworld3/checkpoint/disk_snapshot.py b/src/underworld3/checkpoint/disk_snapshot.py index 9cbf7c22a..6a61a7a38 100644 --- a/src/underworld3/checkpoint/disk_snapshot.py +++ b/src/underworld3/checkpoint/disk_snapshot.py @@ -770,7 +770,8 @@ def _read_swarm_from_sidecar(swarm, sidecar_path: str) -> None: # # Serialisation is *generic over dataclass fields* — no per-class # special code. Handled value types: None, bool, int, float, str, -# numpy.ndarray, list (JSON-encoded), dict (recursive subgroup). Other +# numpy.ndarray, list (JSON-encoded), dict (recursive subgroup), +# dimensional quantities (magnitude + unit string). Other # types (notably sympy expressions in DDtSymbolicState.psi_star) are # marked with `_skipped` and not round-tripped — documented as # a v1.x limitation; consumers either use a non-Symbolic DDt flavor @@ -788,6 +789,40 @@ def _is_h5_attr_scalar(value: Any) -> bool: ) or isinstance(value, (bool, str)) +_MAGNITUDE_SUFFIX = "__magnitude" +_UNITS_SUFFIX = "__units" + + +def _is_quantity(value: Any) -> bool: + """A dimensional value, duck-typed. + + ``uw.quantity`` returns a ``UWQuantity``, which is NOT a + ``pint.Quantity`` subclass, so an isinstance test against either would + miss one of them. Both carry ``magnitude`` and ``units``. + """ + return hasattr(value, "magnitude") and hasattr(value, "units") + + +def _write_quantity(h5group, name: str, value: Any) -> None: + """Store a dimensional value as magnitude + unit string. + + The magnitude goes through the ordinary dispatch (attr for a scalar, + dataset for an array), so an array-valued quantity round-trips too. + """ + h5group.attrs[name + _UNITS_SUFFIX] = str(value.units) + _serialise_field(h5group, name + _MAGNITUDE_SUFFIX, np.asarray(value.magnitude).item() + if np.ndim(value.magnitude) == 0 else np.asarray(value.magnitude)) + + +def _read_quantity(h5group, name: str) -> Any: + """Inverse of :func:`_write_quantity`.""" + units = h5group.attrs[name + _UNITS_SUFFIX] + if isinstance(units, bytes): + units = units.decode() + magnitude = _deserialise_field(h5group, name + _MAGNITUDE_SUFFIX, None) + return uw.quantity(magnitude, str(units)) + + def _serialise_field(h5group, name: str, value: Any) -> None: """Write a Python value into an HDF5 group as attr/dataset/subgroup. @@ -797,6 +832,7 @@ def _serialise_field(h5group, name: str, value: Any) -> None: - attr `__json` for JSON-encodable lists / nested simple structures - dataset `` for numpy arrays - subgroup `` for dict values, recursing + - attrs `__magnitude` + `__units` for dimensional values - attr `__skipped` = '' for anything else """ if value is None: @@ -808,6 +844,9 @@ def _serialise_field(h5group, name: str, value: Any) -> None: if isinstance(value, str): h5group.attrs[name] = value return + if _is_quantity(value): + _write_quantity(h5group, name, value) + return if isinstance(value, np.ndarray): if name in h5group: del h5group[name] @@ -860,6 +899,21 @@ def _group_to_dict(h5group) -> dict: out[k] = _group_to_dict(item) else: out[k] = np.asarray(item[...]) + + # Fold `__magnitude` + `__units` back into one dimensional + # value. Done as a post-pass because the two halves may arrive from + # different loops above (a scalar magnitude is an attr, an array + # magnitude a dataset). + for units_key in [k for k in out if k.endswith(_UNITS_SUFFIX)]: + base = units_key[: -len(_UNITS_SUFFIX)] + magnitude_key = base + _MAGNITUDE_SUFFIX + if magnitude_key not in out: + continue + units = out.pop(units_key) + magnitude = out.pop(magnitude_key) + if isinstance(units, bytes): + units = units.decode() + out[base] = uw.quantity(magnitude, str(units)) return out @@ -886,6 +940,9 @@ def _deserialise_field(h5group, name: str, fallback: Any) -> Any: if (name + "__json") in h5group.attrs: return json.loads(h5group.attrs[name + "__json"]) + if (name + _UNITS_SUFFIX) in h5group.attrs: + return _read_quantity(h5group, name) + if (name + "__skipped") in h5group.attrs: # Skipped at write time — keep the current value rather than # clobber it with a placeholder. diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 1c807d4a7..a74bf0bf1 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -270,6 +270,40 @@ class Constitutive_Model(uw_object): _class_instance_counts = {} @timing.routine_timer_decorator + def _declared_terms(self): + """The parameters this model was given, by the name they were set under. + + Enumerated from ``Parameters`` so every constitutive model satisfies + the contract without writing it out; a model whose terms need a better + account overrides this. + """ + import types + + parameters = getattr(self, "Parameters", None) + if parameters is None: + return [] + + terms, seen = [], set() + for name in sorted(a for a in dir(parameters) if not a.startswith("_")): + try: + value = getattr(parameters, name) + except Exception: + continue + # `dir()` sees anything bound into the namespace, including the + # module imports that leak into it. + if isinstance(value, types.ModuleType) or callable(value): + continue + key = str(value) + if key in seen: + continue # `viscosity` and `shear_viscosity_0` alias + seen.add(key) + terms.append({ + "name": f"{type(self).__name__}.{name}", + "value": getattr(value, "sym", value), + "description": "constitutive parameter", + }) + return terms + def __init__(self, unknowns, material_name: str = None): """ Initialize a constitutive model. diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index aa22e699e..314ff81b6 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -1351,6 +1351,179 @@ class SolverBaseClass(uw_object): """Coordinate system of the underlying mesh.""" return inner_self._owning_solver.mesh.CoordinateSystem + def _constraint_mechanisms(self): + """Every way a constraint can have been put on this solver. + + ONE enumeration, used by the mixed-mechanism guard and by + :meth:`describe`. A mechanism added later and not registered here + breaks the guard first, which is a loud failure — where a second list + kept for reporting would simply omit it from the description and say + nothing. + """ + return { + "essential": list(getattr(self, "essential_bcs", None) or []), + "natural": list(getattr(self, "natural_bcs", None) or []), + "rotated_freeslip": list(getattr(self, "_rotated_freeslip_bcs", None) or []), + "fault_contact": list(getattr(self, "_fault_contact_faults", None) or []), + "multipliers": list(getattr(self, "_multipliers", None) or []), + } + + def _declared_terms(self): + """The terms this solver was GIVEN, by the name they were given under. + + The contract a solver implements so its description can say where a + residual came from. ``F0`` for Stokes is ``-bodyforce``; without this, + a description can show the assembled product and not that the user + wrote ``-rho0 * alpha * g * T * rhat``, nor which name to change. + + Return a list of ``{"name", "value", "description"}``. ``None`` means + the solver has not adopted the contract, and :meth:`describe` reports + that rather than passing it off as "no terms". + """ + return None + + def describe(self, depth=4): + """What this solver solves, as data. + + The residual templates with their symbols and descriptions, the named + expressions they contain — expanded RECURSIVELY, so a constitutive + model written in terms of further named quantities is followed rather + than printed as one opaque value — and the boundary conditions. + + One description, two consumers. :meth:`view` renders it for a reader + and the run transcript serialises it, so the equation a note quotes and + the equation the run recorded cannot drift apart. + + Parameters + ---------- + depth : int, default 4 + How far to follow named expressions into each other. A cycle stops + at the symbol that repeats, whatever the depth. + + Returns + ------- + dict + """ + import sympy + + def unpack(expression, level, seen): + """Named expressions inside ``expression``, and inside those.""" + out = [] + if level > depth: + return out + try: + found = uw.function.fn_extract_expressions(expression) + except Exception: + return out + for named in sorted(found, key=lambda e: str(getattr(e, "symbol", e))): + symbol = str(getattr(named, "symbol", named)) + if symbol in seen: + continue + seen.add(symbol) + value = getattr(named, "sym", None) + description = str(getattr(named, "description", "") or "") + out.append({ + "symbol": symbol, + "latex": sympy.latex(value) if value is not None else None, + "value": str(value) if value is not None else None, + "units": (str(named.units) + if getattr(named, "units", None) else None), + "description": ("" if description == "No description provided" + else description), + "where": unpack(value, level + 1, seen) if value is not None else [], + }) + return out + + forms, seen = {}, set() + for name in ("F0", "F1", "PF0"): + template = getattr(self, name, None) + if template is None: + continue + expression = getattr(template, "sym", None) + if expression is None: + continue + # The template's own symbol and docstring are the equation's + # published names; they belong beside its value. + declared = getattr(type(self), name, None) + forms[name] = { + "symbol": getattr(declared, "name", None), + "description": (getattr(declared, "__doc__", "") or "").strip().split("\n")[0], + "latex": sympy.latex(expression), + "text": str(expression), + } + forms[name]["where"] = unpack(expression, 1, seen) + + mechanisms = self._constraint_mechanisms() + conditions = [] + for kind in ("essential", "natural"): + for bc in mechanisms[kind]: + function = getattr(bc, "fn", None) + if function is None: + function = getattr(bc, "fn_f", None) + conditions.append({ + "mechanism": kind, + "type": str(getattr(bc, "type", kind)), + "boundary": str(getattr(bc, "boundary", "?")), + "latex": sympy.latex(function) if function is not None else None, + "text": str(function) if function is not None else None, + }) + # Rotated free-slip is applied by machinery outside the solver, but the + # solver holds what was asked for — so a description that skipped it + # would report "no boundary conditions" for a model whose entire + # boundary treatment is rotated. + datum = getattr(self, "_rotated_freeslip_datum", None) or {} + for boundary, normal in mechanisms["rotated_freeslip"]: + value = datum.get(boundary) + conditions.append({ + "mechanism": "rotated_freeslip", + "type": "rotated free-slip" if value is None + else "rotated normal datum", + "boundary": str(boundary), + "latex": sympy.latex(value) if value is not None else r"\mathbf{u}\cdot\hat{\mathbf{n}} = 0", + "text": str(value) if value is not None else "u . n = 0", + "normal": "mesh" if normal is None else str(normal), + }) + for fault in mechanisms["fault_contact"]: + conditions.append({ + "mechanism": "fault_contact", + "type": "fault contact", + "boundary": str(getattr(fault, "name", fault)), + "latex": None, "text": None, + }) + + terms = self._declared_terms() + described_terms = None + if terms is not None: + described_terms = [] + for term in terms: + value = term.get("value") + described_terms.append({ + "name": term.get("name"), + "description": term.get("description", ""), + "latex": sympy.latex(value) if value is not None else None, + "text": str(value) if value is not None else None, + "where": unpack(value, 1, set()) if value is not None else [], + }) + + return { + "solver": type(self).__name__, + "unknown": getattr(getattr(self, "u", None), "name", None), + "dim": getattr(self.mesh, "dim", None), + "cdim": getattr(self.mesh, "cdim", None), + "forms": forms, + "boundary_conditions": conditions, + "terms": described_terms, + "terms_declared": terms is not None, + } + + def _describe_where(self, entries, display, Latex, level=0): + """Render the "Where:" tree from :meth:`describe`.""" + for entry in entries: + indent = "\\quad " * (level + 1) + tail = f" \\quad ({entry['description']})" if entry["description"] else "" + display(Latex(f"${indent}{entry['symbol']} = {entry['latex']}${tail}")) + self._describe_where(entry.get("where", []), display, Latex, level + 1) + def _object_viewer(self): '''This will add specific information about this object to the generic class viewer ''' @@ -2260,12 +2433,55 @@ class SolverBaseClass(uw_object): cdef double[::1] vals_view = np.ascontiguousarray(values, dtype=np.float64) CHKERRQ(PetscDSSetConstants(cds.ds, n_constants, &vals_view[0])) - def _update_constants(self): + def _update_constants(self, record=True): """Re-pack current UWexpression values and call PetscDSSetConstants. Called before each solve() to ensure constants are current without requiring JIT recompilation. + + ``record=False`` suppresses the step-transcript entry. Pass it from any + site that pushes constants for its OWN assembly rather than to + dispatch a solve — otherwise the transcript reports one operator as two. + The rotated free-slip loop is such a site: it re-attaches the + auxiliary vector and re-packs before running its own manual Krylov + loop, after the public ``solve()`` has already announced itself. """ + # Refresh mesh.t from the model clock first, so a time-dependent + # expression is repacked with the rest of the constants rather than + # needing its own hook (or a recompile) per timestep. + try: + self.mesh._sync_time_from_model() + except AttributeError: + pass + + # Note the solve in the model's step transcript, if a step is open. This + # is the one place every solver passes through before solving, so one + # hook records them all, in order. A no-op outside a model.step block. + if record: + try: + # Name it by what it SOLVES, not by its auto-generated instance + # id: a transcript reading "Stokes(V) -> AdvDiffusion(T)" is + # auditable, one reading "Solver_8_ -> Solver_14_" is not. + try: + unknown = self.u.name + except Exception: + unknown = "?" + # `name` is what gets printed; `part` is what a column keys on. + # A rendered label is not an identity — two solvers that happen + # to render the same would collapse into one part, and changing + # how the label is built would silently re-partition every + # transcript ever written. + part = f"{type(self).__name__}#{self.instance_number}" + label = f"{type(self).__name__}({unknown})" + model = uw.get_default_model() + # What it solves, not only that it solved: the residual is + # SymPy, so the weak form can be written into the transcript + # exactly as implemented. + model._describe_part(self, part, label) + model._record_step_event("solve", label, part=part) + except Exception: + pass + if not self.constants_manifest or self.dm is None: return @@ -4101,14 +4317,14 @@ class SNES_Scalar(SolverBaseClass): ) - exprs = uw.function.fn_extract_expressions(self.F0) - exprs = exprs.union(uw.function.fn_extract_expressions(self.F1)) - - if len(exprs) != 0: + # Rendered from describe(), so the "Where:" a reader sees and the + # equation the run transcript records come from one description. + where = [] + for form in self.describe()["forms"].values(): + where.extend(form.get("where", [])) + if where: display(Markdown("*Where:*")) - - for expr in exprs: - expr._object_viewer() + self._describe_where(where, display, Latex) display( @@ -5150,14 +5366,14 @@ class SNES_Vector(SolverBaseClass): Latex(eqF1), Latex(eqf0), ) - exprs = uw.function.fn_extract_expressions(self.F0) - exprs = exprs.union(uw.function.fn_extract_expressions(self.F1)) - - if len(exprs) != 0: + # Rendered from describe(), so the "Where:" a reader sees and the + # equation the run transcript records come from one description. + where = [] + for form in self.describe()["forms"].values(): + where.extend(form.get("where", [])) + if where: display(Markdown("*Where:*")) - - for expr in exprs: - expr._object_viewer() + self._describe_where(where, display, Latex) display( Markdown(fr"# Boundary Conditions"),) @@ -6216,10 +6432,9 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): mechanism is already in place and which one was refused. """ - rotated = list(getattr(self, "_rotated_freeslip_bcs", None) or []) + list( - getattr(self, "_fault_contact_faults", None) or [] - ) - multipliers = list(getattr(self, "_multipliers", None) or []) + mechanisms = self._constraint_mechanisms() + rotated = mechanisms["rotated_freeslip"] + mechanisms["fault_contact"] + multipliers = mechanisms["multipliers"] if adding == "solve": # The dispatch reads both lists, so it can only report the pair. diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 7f045997e..4ee75bf84 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -1385,13 +1385,16 @@ def _setup_symbolic_coordinates(self, coordinate_system_type): self._Gamma.y._ccodestr = "petsc_n[1]" self._Gamma.z._ccodestr = "petsc_n[2]" - # Time coordinate — PETSc passes this as petsc_t to all pointwise - # functions. Solvers set dm.time before each solve via solve(time=t). - # Users reference it as mesh.t in expressions (e.g. V0 * sympy.sin(omega * mesh.t)) - from ..utilities.unit_aware_coordinates import TimeSymbol - - self._t = TimeSymbol("t") - self._t._units = None # patched below by _patch_time_units + # Time coordinate. This is a live-rampable ``constants[]`` atom, NOT + # PETSc's ``petsc_t``: the high-level solve() wrappers never set + # petsc_t, so an expression built on it evaluated to zero inside every + # solve (silently — a time-dependent BC was identically zero). Time is + # owned by the orchestration model; ``mesh.t`` reads that clock. + # ``_sync_time_from_model`` repacks it before each solve, from the + # solver's ``_update_constants``, so no kernel is recompiled per step. + self._t = uw.expression( + r"t", 0.0, "model time — the clock on uw.get_default_model().tracker" + ) # Add unit awareness to coordinate symbols if mesh has units or model has scales from ..utilities.unit_aware_coordinates import patch_coordinate_units @@ -4516,30 +4519,60 @@ def CoordinateSystem(self) -> CoordinateSystem: @property def t(self): - r"""Symbolic time coordinate. + r"""Symbolic model time. - PETSc passes a time value (``petsc_t``) to all pointwise residual - and Jacobian functions. Use ``mesh.t`` in expressions to reference - this time without forcing JIT recompilation each timestep. + A live-rampable ``constants[]`` atom carrying the clock owned by the + orchestration model, ``uw.get_default_model().tracker.time``. Every + solver repacks it from that clock immediately before solving, so an + expression built on ``mesh.t`` follows time with no JIT recompilation + per step. - The low-level PETSc solver accepts ``time=t`` to set the value - of ``petsc_t`` for pointwise functions. If not provided, ``petsc_t`` - defaults to 0. Note: the high-level Python ``solve()`` wrappers - do not yet pass ``time=`` through — set it directly via - ``UW_DMSetTime`` at the Cython level if needed. + Maintain the clock as part of the timestepping loop (see + ``docs/developer/guides/HOW-TO-WRITE-UW3-SCRIPTS.md``). A script that + never advances it leaves ``mesh.t`` at zero. - When the scaling system is active, ``mesh.t`` carries time units - (derived from the model's time scale) so that dimensional analysis - works correctly in expressions. + A dimensional clock is non-dimensionalised on the way in, so the value + the kernels see is always in solver units. + + .. note:: + Assign it as part of an expression rather than bare. A boundary + condition takes a Matrix / array form, and a bare atom handed to a + scalar setter is stored by value. Examples -------- >>> omega = 2 * np.pi / period >>> stokes.add_dirichlet_bc((V0 * sympy.sin(omega * mesh.t), 0.0), "Top") - >>> stokes.solve(time=current_time) # sets petsc_t before SNES + >>> model.tracker.time = 1.5 * uw.quantity(1, "Myr") + >>> stokes.solve() # mesh.t picks the clock up """ return self._t + def _sync_time_from_model(self): + """Repack ``mesh.t`` from the model clock. Called by every solver's + ``_update_constants`` immediately before a solve, so an expression + containing ``mesh.t`` sees the current time without a rebuild. + + Silent no-op when the model has no clock: ``mesh.t`` then stays at + whatever it was last set to (0.0 for a fresh mesh), which is the + behaviour a script that never advances a clock already expects. + """ + try: + model = uw.get_default_model() + time = model.tracker.time + except Exception: + return + if time is None: + return + try: + if hasattr(time, "magnitude") or hasattr(time, "_pint_qty"): + time = float(uw.non_dimensionalise(time)) + self._t.sym = sympy.sympify(float(time)) + except Exception: + # A clock we cannot reduce to a number is not worth failing a + # solve over; leave mesh.t as it stands. + return + @property def nullspace_rotations(self): """Symbolic velocity fields for rigid-body rotation null modes. @@ -5208,7 +5241,12 @@ def snapshot_payload(self) -> dict: - ``name``: stable string identifier for the mesh. - ``mesh_version``: current ``_mesh_version`` integer. - - ``coords``: deformed mesh coordinates (numpy array). + - ``coords``: deformed mesh coordinates, in MODEL UNITS — the + representation :meth:`_deform_mesh` writes back. ``mesh.X.coords`` + is the unit-aware view and returns metres when a model declares a + length scale; capturing that and restoring it through + ``_deform_mesh`` would multiply the mesh by the length scale on + every restore, silently and without changing any array's shape. - ``vars``: ``{var.clean_name: gvec_array.copy()}`` for every mesh variable on this mesh. @@ -5216,7 +5254,7 @@ def snapshot_payload(self) -> dict: section / DM-topology data sufficient to rebuild the DM on restore. """ - coords = numpy.asarray(self.X.coords).copy() + coords = numpy.asarray(self._coords).copy() var_arrays: dict[str, numpy.ndarray] = {} for var in self.vars.values(): var._sync_lvec_to_gvec() @@ -5261,7 +5299,7 @@ def apply_snapshot_payload(self, payload: dict) -> None: ) coords = numpy.asarray(payload["coords"]) - expected_shape = numpy.asarray(self.X.coords).shape + expected_shape = numpy.asarray(self._coords).shape if coords.shape != expected_shape: raise SnapshotInvalidatedError( f"mesh {self.name!r}: coordinate shape changed " diff --git a/src/underworld3/model.py b/src/underworld3/model.py index 1a5619ae8..3b0a34981 100644 --- a/src/underworld3/model.py +++ b/src/underworld3/model.py @@ -43,6 +43,106 @@ class PintNativeModelMixin: pass +class _AutoSentinel: + """Sentinel for "the user has not said where the transcript goes". + + Distinct from None, which means "off": the two have to be told apart, + because a default that cannot be switched off is worse than no default. + + Copy-stable on purpose. ``PrivateAttr`` deep-copies its default, and a bare + ``object()`` would come back as a DIFFERENT object per model, so every + identity check against it would silently fail. + """ + + __slots__ = () + + def __repr__(self): + return "" + + def __copy__(self): + return self + + def __deepcopy__(self, memo): + return self + + +_AUTO = _AutoSentinel() + +TRANSCRIPTS_DIR = "transcripts" + + +def _transcript_disabled(): + """Whether the automatic transcript should stay off. + + Off under pytest — 1800 tests should not each leave a directory — and off + when ``UW_TRANSCRIPT`` says so, which is the switch for CI and for anyone + who does not want the files. + """ + setting = os.environ.get("UW_TRANSCRIPT", "").strip().lower() + if setting in ("off", "0", "no", "none", "false"): + return True + if setting: + return False + return "PYTEST_CURRENT_TEST" in os.environ or "pytest" in sys.modules + + +def _launch_stem(): + """A short name for the run, taken from the script that started it.""" + entry = sys.argv[0] if sys.argv else "" + if not entry or entry == "-c": + return "interactive" + stem = os.path.splitext(os.path.basename(entry))[0] + return "".join(c if (c.isalnum() or c in "-_") else "-" for c in stem) or "run" + + +def _launch_manifest(): + """What was invoked, as far as it can be known. + + A programmatic launcher cannot be made reproducible by fiat, but what was + actually run CAN be written down: the command line, the interpreter, the + working directory, the package version, and the commit if there is one. + That is the difference between "I cannot reproduce this" and "I know + exactly what produced it and can decide what to change". + """ + from datetime import datetime, timezone + + import underworld3 as uw + + manifest = { + "started": datetime.now(timezone.utc).isoformat(timespec="seconds"), + "argv": list(sys.argv), + "executable": sys.executable, + "cwd": os.getcwd(), + "underworld3": getattr(uw, "__version__", "unknown"), + "underworld3_path": os.path.dirname(getattr(uw, "__file__", "") or ""), + "python": sys.version.split()[0], + "mpi_size": int(uw.mpi.size), + } + try: + manifest["host"] = os.uname().nodename + except Exception: + pass + # The entry script is copied beside this; imported modules are NOT, so a + # commit id is what covers the rest when the work is under version control. + try: + import subprocess + + sha = subprocess.run( + ["git", "rev-parse", "HEAD"], capture_output=True, text=True, + timeout=5, cwd=os.getcwd(), + ) + if sha.returncode == 0: + manifest["git_commit"] = sha.stdout.strip() + dirty = subprocess.run( + ["git", "status", "--porcelain"], capture_output=True, text=True, + timeout=5, cwd=os.getcwd(), + ) + manifest["git_dirty"] = bool(dirty.stdout.strip()) + except Exception: + pass + return manifest + + class ModelState(Enum): """Model lifecycle states""" @@ -54,6 +154,182 @@ class ModelState(Enum): ERROR = "error" +class ModelStep: + """What one timestep did — the transcript entry for a ``model.step`` block. + + Ordered, so ``[e.name for e in step.events]`` is the sequence of operators + the step actually applied. That sequence is what makes a step auditable + (did this run do what the write-up says?) and what a replay or an adjoint + needs in order to walk the run backwards. + """ + + __slots__ = ("index", "t0", "dt", "label", "events", "completed", "snapshot", + "wall") + + def __init__(self, index, t0, dt, label=None): + self.index = index + self.t0 = t0 + self.dt = dt + self.label = label + self.events = [] + self.completed = False + # Seconds of wall clock the block took. Not physics, but the number you + # want when watching a run: a step that suddenly takes ten times as + # long is the first sign of a solver in trouble. + self.wall = None + # The state this step STARTED from, when the recording policy kept one. + # Taken before the operators ran, which is the only correct point: a + # DDt shifts its history in its post-solve hook, so a snapshot taken + # afterwards holds the shifted history rather than the step's input. + self.snapshot = None + + @property + def restorable(self): + """Whether this step kept the state it started from.""" + return self.snapshot is not None + + @property + def t1(self): + """The end of the interval this step covers.""" + return self.t0 + self.dt + + def _record(self, kind, name, **detail): + self.events.append({"kind": kind, "name": name, **detail}) + + def as_dict(self): + """This step as plain JSON-able data — the on-disk log's line format. + + Dimensional values become ``{"magnitude": ..., "units": ...}``, the + same split the on-disk snapshot uses, so a log written by a run with + units is readable without a live model to interpret it. + + ``snapshot`` is deliberately absent: it is megabytes of field data and + does not survive the process. ``restorable`` records whether one was + held, which is what a reader of the log can act on. + + Each value keeps the units the run actually held it in, which is why + ``t0`` may read in Myr beside a ``dt`` in seconds: the clock came from + the tracker and the interval from ``estimate_dt()``. The log is a + transcript of the run, not a tidied report of it — convert on the way out. + """ + return { + "kind": "step", + "index": self.index, + "label": self.label, + "t0": _jsonable_quantity(self.t0), + "t1": _jsonable_quantity(self.t1), + "dt": _jsonable_quantity(self.dt), + "completed": bool(self.completed), + "restorable": bool(self.restorable), + "wall": None if self.wall is None else float(self.wall), + "events": [dict(e) for e in self.events], + } + + def __repr__(self): + state = "" if self.completed else " ABANDONED" + seq = " -> ".join(f"{e['kind']}:{e['name']}" for e in self.events) or "(nothing)" + tag = f" {self.label!r}" if self.label else "" + return f"" + + +def _quantity_parts(value): + """``(magnitude, unit string or None)`` for a value that may be dimensional.""" + if hasattr(value, "magnitude") and hasattr(value, "units"): + try: + return float(value.magnitude), str(value.units) + except (TypeError, ValueError): + return None, str(value.units) + try: + return float(value), None + except (TypeError, ValueError): + return None, None + + +def _abbreviate_unit(unit): + """A short unit name for a column header. Falls back to the full name.""" + if unit is None: + return "" + return { + "second": "s", "minute": "min", "hour": "hr", "day": "d", + "year": "yr", "kiloyear": "kyr", "megayear": "Myr", "gigayear": "Gyr", + "meter": "m", "kilometer": "km", "kelvin": "K", "kilogram": "kg", + }.get(str(unit), str(unit)) + + +def _in_units_of(value, unit): + """``value`` as a bare number in ``unit``, or its own magnitude if it cannot + be converted. A text log is a report: one time column, one unit.""" + if unit is None: + magnitude, _ = _quantity_parts(value) + return magnitude + try: + return float(value.to(unit).magnitude) + except Exception: + magnitude, _ = _quantity_parts(value) + return magnitude + + +def _bare(value, unit): + """A serialised value as a bare number, converted to ``unit`` if it can be. + + ``value`` is what :meth:`ModelStep.as_dict` produced: a float, or a + ``{"magnitude", "units"}`` pair. The text log shows one time column in one + unit, so a ``dt`` in seconds beside a clock in Myr is converted rather than + printed as it stands. + """ + if not isinstance(value, dict): + try: + return float(value) + except (TypeError, ValueError): + return float("nan") + magnitude = value.get("magnitude") + units = value.get("units") + if unit is None or units is None or str(units) == str(unit): + try: + return float(magnitude) + except (TypeError, ValueError): + return float("nan") + try: + import underworld3 as uw + + return float(uw.quantity(float(magnitude), str(units)).to(unit).magnitude) + except Exception: + try: + return float(magnitude) + except (TypeError, ValueError): + return float("nan") + + +def _pretty_time(value): + """A compact, readable rendering of a clock value for a log note.""" + magnitude, unit = _quantity_parts(value) + if magnitude is None: + return str(value) + if unit is None: + return f"{magnitude:.6g}" + return f"{magnitude:.6g} {_abbreviate_unit(unit)}" + + +def _jsonable_quantity(value): + """A number, or a dimensional value split into magnitude and units. + + Duck-typed, because ``uw.quantity`` returns a ``UWQuantity``, which is not + a ``pint.Quantity`` subclass — an isinstance test against either would + miss one of them. Both carry ``magnitude`` and ``units``. + """ + if hasattr(value, "magnitude") and hasattr(value, "units"): + magnitude = value.magnitude + try: + magnitude = float(magnitude) + except (TypeError, ValueError): + magnitude = str(magnitude) + return {"magnitude": magnitude, "units": str(value.units)} + try: + return float(value) + except (TypeError, ValueError): + return str(value) + + class Model(PintNativeModelMixin, BaseModel): """ Central orchestrator for Underworld3 simulations. @@ -141,6 +417,45 @@ class Model(PintNativeModelMixin, BaseModel): # src/underworld3/checkpoint/tracker.py. _tracker: Any = PrivateAttr(default=None) + # The step transcript: an ordered record of what each timestep actually did. + # ``_open_step`` is the ModelStep currently in progress (None outside a + # ``with model.step(dt):`` block); ``_transcript`` is the bounded history of + # completed steps. See :meth:`step`. + _open_step: Any = PrivateAttr(default=None) + _transcript: Any = PrivateAttr(default_factory=list) + _transcript_limit: Any = PrivateAttr(default=512) + + # Optional on-disk log of the transcript: one JSON object per line, appended + # and flushed as each step closes. See :attr:`transcript_file`. + # ``_AUTO`` until the user says otherwise: a transcript lands in + # ``transcripts/-