diff --git a/.gitignore b/.gitignore index cd3992068..c3235d615 100644 --- a/.gitignore +++ b/.gitignore @@ -155,6 +155,7 @@ uw_*.msh # Visualization outputs (exclude generated, allow curated assets) docs/examples/**/output/** +docs/examples/**/transcripts/** docs/examples/**/*.png docs/examples/**/*.jpg docs/examples/**/*.gif diff --git a/CLAUDE.md b/CLAUDE.md index 0f7868a1d..82745c019 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -141,6 +141,13 @@ This keeps feature branches independent and makes cross-pollination of fixes str - Cherry-pick to `main` if critical → tag patch release - Cherry-pick to active feature branches (underworld-claude handles this) +### Adversarial Review Before the PR Opens +**Every branch gets an adversarial review before it becomes a PR**, and again +after any substantial post-review commit. The checklist — parallel rank +asymmetry, frame/unit boundaries, determinism, tests that cannot fail, CI +reach, and the solver/transcript contracts — is in +`docs/developer/guides/adversarial-review.md`. Post the findings on the PR. + ### Git Worktrees for Session Isolation **Use a worktree for any multi-file change** (docs cleanup, refactoring, features). Multiple Claude sessions sharing one working directory will overwrite each other's work. diff --git a/docs/developer/UW3_STYLE_CHARTER.md b/docs/developer/UW3_STYLE_CHARTER.md index cc470ee86..17bb5823d 100644 --- a/docs/developer/UW3_STYLE_CHARTER.md +++ b/docs/developer/UW3_STYLE_CHARTER.md @@ -88,6 +88,7 @@ These settle the June-2026 drift (see `docs/reviews/2026-07/API-CONSISTENCY-REVI | Solver capabilities | Anything that configures or reads one solver is a METHOD on that solver, lazily importing its `utilities/*` implementation (the `boundary_flux` pattern) — never a free function as the documented entry point. | | Namespaces | Every user-facing module is exported from its subpackage `__init__`/`__all__` in the PR that creates it. No deep-import-only features. | | Docstrings | NumPy/Sphinx style with RST `:math:`. This SUPERSEDES the Markdown-for-pdoc prescription still printed in `UW3_Style_and_Patterns_Guide.md` — that section is wrong; do not follow it. | +| Constants in symbolic forms | A number that enters a residual, a scheme or a constitutive law is EXACT where it is exact: `sympy.Rational(1, 2)` for a theta of one half, `sympy.Integer(4)` for a weight of four, never `0.5` or `4.0`. A float cannot cancel, and prints as `0.5` or `64.0` in the form the run records. A regulariser that keeps a denominator finite is `uw.maths.functions.vanishing`, which prints as :math:`\varepsilon`, never a literal `1e-30`. User-supplied floats are made exact at the boundary with `sympy.nsimplify(x, rational=True)`. | ## 7. Data Access diff --git a/docs/developer/design/run-plan-and-transcript.md b/docs/developer/design/run-plan-and-transcript.md new file mode 100644 index 000000000..58ad60264 --- /dev/null +++ b/docs/developer/design/run-plan-and-transcript.md @@ -0,0 +1,174 @@ +# The run plan 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 plan** is what the run is *supposed* to do: which parts run, in what +order, on what alignment. It is the same for every cycle of a well-behaved +run. + +**The transcript** is what the run *actually did*, including the things that +were never planned — 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 plan**; the +rows are the **transcript**. A run in which every step is `A` followed the +plan. A `B` is a step that did something else. + +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 plan**. + +## Vocabulary + +The terms below are where the ideas came from, and they are deliberately NOT +the words the code and the figures use. "Score", "bar", "note", "rest" and +"simile" were useful for getting the model right and are forced as public +vocabulary: what ships says **transcript**, **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: + +- **Stable part order.** Parts appear in a fixed 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. +- **Absence is written down.** Every part accounts for every step. A part + that did nothing is marked as having done nothing, 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 ≠ plan | 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. 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 plan, then declared plan + +The plan is **inferred** today — the figure takes the most common step as the +norm. That costs nothing and can only ever say *this step differs from its +neighbours*. + +A **declared** plan — the script stating what a step 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 log or record: a transcript is what actually +happened, including the false starts and the re-takes, which is precisely the +thing being kept. + +The API followed: what was `model.journal` is `model.transcript`, and the +renderers are named for it — `transcript_table` sets it in text, +`transcript_figure` draws it, `transcript_flowchart` draws the sequence it +implies. Nothing user-facing is called a score: the word belongs to the +derivation above, not to the shipped vocabulary. + +"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/CODE-REVIEW-PROCESS.md b/docs/developer/guides/CODE-REVIEW-PROCESS.md index 3d269df70..6227488fc 100644 --- a/docs/developer/guides/CODE-REVIEW-PROCESS.md +++ b/docs/developer/guides/CODE-REVIEW-PROCESS.md @@ -76,6 +76,17 @@ The code review process serves to: ## Review Process Workflow +### Phase 0: Adversarial review (before the PR opens) + +Run the pass described in [adversarial-review.md](adversarial-review.md): an +independent attempt to refute the change, aimed at the failure modes this +codebase produces — rank asymmetry, frame and unit boundaries, determinism, +tests that cannot fail, CI reach, and the contracts a change must not quietly +leave. Fix what it finds, then open the PR and post the findings on it. + +Every PR, including one-line bug fixes and docs-only changes; again after any +substantial post-review commit. + ### Phase 1: Preparation 1. **Author**: Prepare change materials diff --git a/docs/developer/guides/HOW-TO-WRITE-UW3-SCRIPTS.md b/docs/developer/guides/HOW-TO-WRITE-UW3-SCRIPTS.md index 1646bd499..c5fcd6629 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,500 @@ 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) + Stokes(V)> + Stokes(V)> + 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. + +What 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 — with the state it started from and the +order the operators were applied in both on record. + +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] AdvDiffusion(T) > shift EulerianSUPG(T) > Stokes(v) + 1 0.501546 0.325639 0.09 ok [convect] AdvDiffusion(T) > shift EulerianSUPG(T) > Stokes(v) + 2 0.990939 0.489393 0.09 ok [convect] AdvDiffusion(T) > shift EulerianSUPG(T) > Stokes(v) + 3 1.51459 0.523655 0.09 ok [convect] AdvDiffusion(T) > shift EulerianSUPG(T) > Stokes(v) + 4 30.5663 29.0517 0.09 ABANDONED [too big] AdvDiffusion(T) > shift EulerianSUPG(T) > 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] AdvDiffusion(T) > shift EulerianSUPG(T) > 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. + +**The outcome column says how the step went, not only that it ran.** A step +whose solves all converged reads `ok`; one where a solve did not reads +`DIVERGED`, and the solve is named underneath with its reason and its work: + +``` + 7 2.63841 0.523655 4.81 DIVERGED [convect] Stokes(v) > ... + !! SNES_Stokes(v): DIVERGED_LINEAR_SOLVE after 6 its (1200 ksp), |F| 3.11e-04 + ~~ RuntimeWarning: Stokes: the velocity block fell back to 'gamg' — no mesh hierarchy was available... +``` + +Warnings raised inside the block are recorded the same way, with where they +came from. The transcript gets a copy, not the only copy: the warning is still +shown, and your own filters still apply. "The velocity block fell back to +gamg" changes what the numbers mean, and a record that kept the residual norms +but not that line would be an account of the run with the explanation removed. + +The figure marks the same three states per solve — converged, converged with a +fieldsplit block that hit its iteration cap, and diverged. The middle one is +worth the separate mark: a capped block did not solve, so the Schur operator +was applied through a velocity solve that was still moving, and the outer SNES +can still report CONVERGED (#625). + +**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-plan-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. + +**Reproducing a run from its record.** `read_transcript()` gives the steps +back in order, each with the interval it covered and the operators it applied, +and a run kept with `record_every` can be re-entered at any recorded step with +`load_state(entry.snapshot)`. That is what makes a result reproducible from +the record rather than from a re-run of the script: the same state, the same +operators, in the same order, without anyone having decided in advance which +arrays would be wanted. + +### 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 +1187,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 +1229,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 +1257,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/guides/adversarial-review.md b/docs/developer/guides/adversarial-review.md new file mode 100644 index 000000000..c9e45368b --- /dev/null +++ b/docs/developer/guides/adversarial-review.md @@ -0,0 +1,94 @@ +# Adversarial review: refute the change before it opens as a PR + +**Audience**: anyone (human or AI) preparing an Underworld3 branch for review. + +The [code review process](CODE-REVIEW-PROCESS.md) describes how a change is +reviewed once it is a PR. This document is the pass that runs *before* that: an +independent attempt to refute the change, aimed at the failure modes this +codebase actually produces. Four parallel reviewers on three already-pushed +branches found twenty-odd defects in 2026-09, including a SEGV reachable from +the API's own error message and a JIT change that aborted half of all `np>=2` +runs. Every one of them was on work already pushed, which is why the ordering +matters. + +## Running it + +1. Read `docs/developer/UW3_STYLE_CHARTER.md`. It is two pages and it is + normative; the style findings that reach a PR are almost always in it. +2. Run the review on the branch diff — an independent pass whose goal is to + break the change, not to confirm it. +3. Fix what it finds, then open the PR. +4. Post the findings on the PR, including the attacks that failed. +5. Run it **again** after any substantial post-review commit. A refactor + landing after the review is when a second pass is worth most. + +The posted review is terse: findings and evidence (numbers, `file:line`, probe +name), one line each; a short list of attacks that failed, with numbers; merge +conditions as bullets. Voice is "we". Cite a rule at a defect +("Charter §4: uncommented swallow"), never to justify the review itself. + +## What we attack + +**Parallel.** Rank asymmetry, collective ordering, empty-rank edges. A +collective call inside a rank-dependent branch hangs the job rather than +failing it, so it survives every serial test. + +**Frames and units.** Does the value cross a boundary in the units the other +side expects? Snapshot capture read `mesh.X.coords` (metres) while restore +wrote model units, and silently rescaled the mesh on every restore. + +**Determinism.** Anything whose result can depend on hash seed, dict order, or +partition. A JIT change made the generated C hash-seed dependent and aborted +roughly half of multi-rank runs. + +**Tests that cannot fail.** Every regression test must be shown to fail on its +bug. A ramp of 1000 → 1e-6 where `"1.0e-6"` sorts as a prefix of `"1.0"` moves +nothing; an MPI property that is seed-flaky caught its bug zero times in six. +Delete a test that cannot fail. A flaky test is the test's fault until proven +otherwise. + +**Diagnostics that do not discriminate.** Before building on a measure, check +it gives a different answer with the feature on and off. Three field-level +measures of population control did not. + +**CI reach.** Does the new test file match the CI glob? Eleven test files — +an entire subsystem's suite — never ran because they fell outside the pattern. + +## Contracts a change must not quietly leave + +These are enforced by tests. The review checks the test is still the one doing +the enforcing, and that no exemption list has grown a new entry without a +reason beside it. + +**A solver says what it solves.** Every solver class declares `_solver_terms`, +a tuple of `(attribute, description)` naming the terms it is given, from which +`_declared_terms()` is built. `describe()` then returns the residual templates, +the named expressions inside them (followed recursively into the constitutive +model), the boundary conditions and those terms. `view()` renders that +description and the run transcript serialises it, so a note and a run cannot +quote different equations. A new solver that skips the roster fails +`tests/test_0016_solver_description_contract.py`; adding its name to that +file's exemption list is not the fix. + +**Constraints have one enumeration.** `_constraint_mechanisms()` lists every +way a constraint can be on a solver, and both the mixed-mechanism guard +(#464) and `describe()` read it. A mechanism added without registering there +breaks the guard — loudly — rather than vanishing from every description. + +**A run records what it did.** `model.step(dt)` is the transaction that carries +the transcript. A new operation that changes model state inside a step — a +solve, a history shift, a mesh deformation — records itself, or the transcript +silently understates the run. See +[HOW-TO-WRITE-UW3-SCRIPTS](HOW-TO-WRITE-UW3-SCRIPTS.md) and +`docs/developer/design/run-plan-and-transcript.md`. + +**Named quantities keep their names.** A coefficient written as +`uw.expression(r"\rho_0 \alpha g", ...)` appears in the description under that +name. An anonymous float collapses into the assembled product and the +transcript can only show the number. Examples in `docs/` name their +coefficients. + +## Where the reviews live + +`docs/reviews/[YYYY-MM]/`, indexed by `docs/reviews/README.md`, and posted on +the PR itself. diff --git a/docs/developer/index.md b/docs/developer/index.md index a8c1f5123..16eff5b9d 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -130,6 +130,7 @@ guides/HOW-TO-WRITE-UW3-SCRIPTS guides/notebook-style-guide guides/GMSH_INTEGRATION_GUIDE guides/CODE-REVIEW-PROCESS +guides/adversarial-review guides/style-gates guides/SPELLING_CONVENTION guides/version-management @@ -145,6 +146,7 @@ guides/mpi-hang-supervision :hidden: :caption: Design Documents +design/run-plan-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..19b1627b4 --- /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-plan-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/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..bf3faf592 100644 --- a/src/underworld3/__init__.py +++ b/src/underworld3/__init__.py @@ -218,9 +218,17 @@ 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_table, + transcript_figure, + transcript_key, +) 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..2647afb16 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -1351,6 +1351,277 @@ class SolverBaseClass(uw_object): """Coordinate system of the underlying mesh.""" return inner_self._owning_solver.mesh.CoordinateSystem + def _transcript_identity(self): + """``(part, label)`` for this solver in a run transcript. + + 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. + + ``label`` 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. + """ + try: + unknown = self.u.name + except Exception: + unknown = "?" + return (f"{type(self).__name__}#{self.instance_number}", + f"{type(self).__name__}({unknown})") + + def _record_solve_outcome(self, report): + """Write the outcome of the solve just finished onto its transcript event. + + The ``solve`` event is recorded BEFORE the solve, because the order + operators ran in is the thing the transcript exists to preserve. The + outcome is only known afterwards, so it is attached to that same event + rather than appended as a second one: one solve, one row, carrying both + where it came in the sequence and how it went. + + Without this a transcript says a run solved 300 times and nothing about + the fifty that diverged — which is the difference between a record of a + run and an account of it. + """ + if report is None: + return + try: + model = uw.get_default_model() + part, _ = self._transcript_identity() + model._record_solve_outcome(part, report) + except Exception: + pass + + 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 []), + } + + #: The terms this solver is given, as ``(attribute, description)`` pairs, + #: in the order a reader should meet them. A subclass that sets this gets + #: :meth:`_declared_terms` for free; ``None`` means the solver has not + #: adopted the contract, which :meth:`describe` reports as such rather + #: than passing it off as "no terms". Inherited, so a derived solver + #: extends its parent's roster rather than restating it: + #: ``_solver_terms = SNES_Stokes._solver_terms + (("rho", "..."),)``. + #: Enforced by ``tests/test_0016_solver_description_contract.py``. + _solver_terms = None + + def _declared_terms(self): + """The terms this solver was GIVEN, by the name they were given under. + + The contract a solver satisfies 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. + + Returns a list of ``{"name", "value", "description"}`` built from + :attr:`_solver_terms`, with the constitutive model's own terms + appended — the model is where most of the named physics lives, and a + solver that reported only its own attributes would stop at + ``constitutive_model`` as an opaque object. + + Override only where the roster cannot express what was given; the + declaration is the intended path, so that adopting the contract is a + line of data rather than a method to keep in step with :meth:`describe`. + """ + roster = type(self)._solver_terms + if roster is None: + return None + + terms = [] + for attribute, description in roster: + # A term can be a property that is not answerable yet — a solver + # described before it is configured, at collection time or in a + # notebook. Report that in place of the value; raising here would + # take out view() and the transcript for a solver that is merely + # incomplete. + try: + value = getattr(self, attribute, None) + value = getattr(value, "sym", value) + except Exception as exc: + value, description = None, f"{description} (unavailable: {exc})" + terms.append({ + "name": attribute, + "value": value, + "description": description, + }) + + model = getattr(self, "constitutive_model", None) + if model is not None: + declared = getattr(model, "_declared_terms", None) + if callable(declared): + terms.extend(declared() or []) + else: + terms.append({ + "name": "constitutive_model", + "value": None, + "description": f"{type(model).__name__} " + f"(does not declare its terms)", + }) + return terms + + 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 = list(uw.function.fn_extract_expressions(expression)) + except Exception: + return out + # A value that IS a named expression — Parameters.diffusivity set + # to the user's own uw.expression — contains no sub-expressions, + # so extraction returns nothing and the user's name, units and + # description would never appear. It is the child. + if (hasattr(expression, "symbol") and hasattr(expression, "sym") + and not any(e is expression for e in found)): + found.append(expression) + 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 ''' @@ -1504,6 +1775,7 @@ class SolverBaseClass(uw_object): ) self._solve_report = report self._solve_history.append(report) + self._record_solve_outcome(report) return report def _capture_rotated_report(self, info): @@ -1547,6 +1819,7 @@ class SolverBaseClass(uw_object): ) self._solve_report = report self._solve_history.append(report) + self._record_solve_outcome(report) return report def guard(self, *, wall_per_step): @@ -2260,12 +2533,47 @@ 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=False): """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 — but it is ALSO called by a + # Parameter change, the continuation alpha toggle, reaction assembly + # and residual-field evaluation, none of which is a solve. Only the + # solve() bodies pass record=True; a call from anywhere else must not + # write a solve event, or one solve reads as several in the record + # (found in review: 2-4 events per solve under continuation). + if record: + try: + part, label = self._transcript_identity() + 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 @@ -4038,7 +4346,7 @@ class SNES_Scalar(SolverBaseClass): ierr = DMSetAuxiliaryVec_UW(dm.dm, NULL, 0, 0, cmesh_lvec.vec); CHKERRQ(ierr) # Update constants (e.g. changed material params) before solve - self._update_constants() + self._update_constants(record=True) # Pure-Neumann scalar problems: attach a constant nullspace # to the (now set-up) Jacobian. No-op unless @@ -4101,14 +4409,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( @@ -5086,7 +5394,7 @@ class SNES_Vector(SolverBaseClass): ierr = DMSetAuxiliaryVec_UW(dm.dm, NULL, 0, 0, cmesh_lvec.vec); CHKERRQ(ierr) # Update constants (e.g. changed material params) before solve - self._update_constants() + self._update_constants(record=True) # Custom geometric-MG prolongation on the (top-level vector) PC, if # registered via set_custom_fmg or owned by an adapt() mesh. Mirrors the @@ -5150,14 +5458,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"),) @@ -5802,7 +6110,7 @@ class SNES_MultiComponent(SolverBaseClass): cmesh_lvec = self.mesh.lvec ierr = DMSetAuxiliaryVec_UW(dm.dm, NULL, 0, 0, cmesh_lvec.vec); CHKERRQ(ierr) - self._update_constants() + self._update_constants(record=True) self._snes_solve_with_retries(gvec, divergence_retries, verbose) @@ -6216,10 +6524,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. @@ -9616,7 +9923,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): UW_DMSetTime(_time_dm_stokes.dm, t_nd) self.mesh.update_lvec() self.dm.setAuxiliaryVec(self.mesh.lvec, None) - self._update_constants() + self._update_constants(record=True) # guard() refuses rotated free-slip, but the BC can be added AFTER arming. # Re-check here: this path never reaches the instrumentation, so an armed @@ -9669,7 +9976,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.dm.setAuxiliaryVec(self.mesh.lvec, None) # Update constants (e.g. changed material params) before solve - self._update_constants() + self._update_constants(record=True) gvec = self.dm.getGlobalVec() gvec.setArray(0.0) 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..8d0845db1 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().astimezone().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,196 @@ 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 needs in + order to reproduce it. + """ + + __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(_operator_text(e) for e in self.events) or "(nothing)" + tag = f" {self.label!r}" if self.label else "" + return f"" + + +def _operator_text(event): + """One operator as a reader sees it, in every view: ``Stokes(v)`` — the + name the user wrote, not ``solve:SNES_Stokes(v)`` — and a history shift + as ``shift EulerianSUPG(T)``. Anything else keeps its kind as a prefix.""" + from underworld3.utilities.transcript_report import _short_operator + + kind, name = event.get("kind"), _short_operator(event.get("name", "?")) + if kind == "solve": + return name + if kind == "history_shift": + return f"shift {name}" + return f"{kind}:{name}" + + +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 +431,48 @@ 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/-