Check an MRI calculation from its saved fields

The useful question is concrete: do the field values, their reported stresses, and the convergence tests describe the same calculation?

This public edition adjusts file-path metadata and reference labels. The original 852,417-byte replay has SHA-256 8e5918ecefd472a171ce7a7a17c2ae44659a71db9e3ac36bf2c4d1aac9e2972c. All 277 compressed field frames here are byte-for-byte identical to it. The metadata edition was also qualified once: 4.24 seconds timed wall and 25,296,896 bytes peak resident memory; its source and a derived timing summary are inside the download.

Open the exact replay bundle

Download linear-channel-replay-public.zip, 849,925 bytes, and the separate saved-field measuring script. Put both files in one new working directory. The ZIP's SHA-256 is:

680ce195877487db5520cde1b4233fea89e91e46589daca562045feff24c4ef9

With Python 3.14.3, the interpreter used for the original qualification, run:

shasum -a 256 linear-channel-replay-public.zip

unzip linear-channel-replay-public.zip -d mri-replay

cd mri-replay

python3 -I -B source/inspect_saved.py --bundle .

python3 -I -B ../measure_saved_fields.py --bundle . --out ../field-measurements.json

Both scripts use the Python standard library. Compare the checksum before extracting. The bundled inspector checks 299 manifest-listed files and 277 frames, then reports PASS_SAVED_FIELD_INSPECTION_ONLY and newTimeSteps: 0. The separate script first verifies those same file hashes, then measures energy, stresses, and Fourier coefficients from the four arrays. Its result is PASS_READ_ONLY_FIELD_MEASUREMENTS. The output filename must be new.

The ZIP has 300 members because manifest.json does not hash itself. Four files under source/ contain the calculation and inspection code. records/receipt.json.gz records the original eleven cases and 45 passed comparisons. A fresh inspection confirms the saved data and that receipt; it does not repeat the original evolution. The original supervised run took 4.67 seconds and reached 26,001,408 bytes of peak resident memory on its author machine. These measured costs belong to that run; timing on another computer can differ.

Read one frame before reading a growth curve

Open records/keplerian-fastest/00600.json.gz with any gzip-capable JSON reader. It contains the frame at code time 1.5. Its field entry contains, in order, radial velocity perturbation, azimuthal velocity perturbation, radial magnetic perturbation, and azimuthal magnetic perturbation. Each array has 256 samples. Position index zero is the start of the periodic vertical interval; index 64 is one-quarter of a wavelength.

The coordinates are radial $x$ , orbital $y$ , and vertical $z$ . The magnetic entries use velocity units,

b_x=\frac{\delta B_x}{\sqrt{\mu_0\rho}},\qquad b_y=\frac{\delta B_y}{\sqrt{\mu_0\rho}}.

Here $\rho$ is the constant density and $\mu_0$ is magnetic permeability. Velocities vary as a sine in height, magnetic perturbations as a cosine. Vertical fluid velocity is zero, and the background vertical field is constant. The spatial solver advances four arrays on a periodic mesh; its vertical derivative uses neighboring samples.

For physical units, choose a box height $L$ , a reference frequency $\Omega_{\rm ref}$ , and a density $\rho$ . The stored interval has height one; time is in $\Omega_{\rm ref}^{-1}$ , and velocity in $L\Omega_{\rm ref}$ . Multiplying stored energy or stress by $\rho L^2\Omega_{\rm ref}^2$ gives an energy density or stress in physical units. The code density is one.

The record's tension is the frequency $K=kv_A$ in these fixed units: $k$ is the vertical wavenumber and $v_A$ the background Alfvén speed. Rotating reference cases have $\Omega=\Omega_{\rm ref}=1$ , so its numerical value also equals $K/\Omega$ . The nonrotating Alfvén case has $\Omega=0$ and tension: 0.8; its clock still uses $\Omega_{\rm ref}$ . Dividing by that case's zero rotation rate would erase the meaning of the saved test.

Four saved fields at time1.5 and the same-height products forming Reynolds and Maxwell stress — Four saved fields and their local stress products, with separate layouts for narrow screens.

At index 64, the two velocities are approximately $3.08033\times10^{-4}$ and $3.07998\times10^{-4}$ . The two magnetic entries are close to zero because the cosine vanishes there. Their large values occur a quarter-wavelength away. This phase difference is why a picture of velocity arrows alone cannot show where magnetic stress is largest. The figure uses the actual saved frame; it adds no evolution.

Form the stresses from products at the same height

For an equally spaced mesh with $N$ samples, calculate

\begin{aligned} E&=\frac{\rho}{2N}\sum_{j=0}^{N-1} \Bigl[u_{x,j}^2+u_{y,j}^2\\ &\qquad\qquad{}+b_{x,j}^2+b_{y,j}^2\Bigr],\\ T_{\rm R}&=\frac{\rho}{N}\sum_j u_{x,j}u_{y,j},\\ T_{\rm M}&=-\frac{\rho}{N}\sum_j b_{x,j}b_{y,j}. \end{aligned}

The first line is the volume-averaged perturbation energy density. The other two lines are Reynolds and Maxwell stresses: correlated motions and magnetic tension transport orbital momentum. With the stated coordinates, positive total stress $T_{\rm R}+T_{\rm M}$ means outward radial transport. Multiplying the fields' separate means would give the wrong result: each sinusoidal component averages to zero, even though its product with another component need not.

The fastest case ends at code time 3 with

\begin{aligned} E(0)&=1.3333333333\times10^{-8},\\ E(3)&=1.2002283886\times10^{-6},\\ T_{\rm R}(3)&=4.5007224026\times10^{-7},\\ T_{\rm M}(3)&=7.5015613910\times10^{-7}. \end{aligned}

These numbers use code energy-density units. Their magnetic-to-velocity stress ratio is approximately 1.666746, close to the continuum fastest-mode value $5/3$ . The small difference is a finite-mesh error. This ratio follows from the particular linear eigenmode. A turbulent disk needs its own stress-ratio calculation.

Compare work with energy, using the right sampling

Background shear supplies energy at the rate

\frac{dE}{dt}=q\Omega\left(T_{\rm R}+T_{\rm M}\right),

where $q=-d\ln\Omega/d\ln R$ measures differential rotation. The saved integratedShearWork was accumulated from stresses at the four Runge–Kutta stages of every time step. It was computed independently of the final energy difference. For the fastest case it is $1.1868950553\times10^{-6}$ . Subtract it from $E(3)-E(0)$ , then divide the absolute residual by $E(0)$ . The result is about $2.29\times10^{-11}$ .

The measuring script also applies the trapezoid rule to stresses in the 25 saved frames. That gives $1.1903702511\times10^{-6}$ , about 0.293% above the energy increase. The wider spacing between saved frames explains this difference. A trapezoid drawn over a rising exponential lies above the curve; the original evolution uses many smaller steps and four stage values per step. The tight original work comparison uses the stage integral. A later integral over sparse snapshots has its own sampling error.

Separate spatial and temporal error

The three spatial cases start with the same continuum eigenvector and end at the same time. Their relative four-component errors are:


Vertical samples	Relative field error
64	0.001287717890
128	0.000321942474
256	0.000080486418

The script projects each final field onto the declared sine or cosine and compares that vector with the saved analytic target. Doubling the mesh reduces the error by factors of 3.99984 and 3.99996, as expected for the centered derivative's second-order spatial accuracy.

Saved-field relative errors on three spatial meshes and three time steps, with reduction ratios — Spatial and temporal errors measured again from saved fields.

The three temporal cases instead keep 32 samples and use time steps 0.02, 0.01, and 0.005. Their initial vector and target use the centered derivative's effective wavenumber. This removes the continuum spatial mismatch from that comparison. The errors are approximately $5.009\times10^{-10}$ , $3.148\times10^{-11}$ , and $1.973\times10^{-12}$ ; halving the time step improves them by factors 15.9122 and 15.9574. A fourth-order method approaches a factor of 16. Mixing a continuum target into this test would leave a spatial error floor and hide the time-order result.

The saved target is named expectedReferenceProjection; referenceKind is semidiscrete-eigenmode for temporal cases and continuum for the others. The analytic reference supplies the original seed and comparison targets to the controller. The spatial stepper receives the four fields, rotation, shear, Alfvén speed, and time step; it receives no expected final array. Remeasuring products checks the saved diagnostics. The continuum limits and two convergence sequences test different scientific properties. A file checksum establishes file identity. Physical accuracy comes from the separate scientific comparisons.

Optional fresh evolution and the first failure

To advance the same eleven cases again, use a new output directory from the extracted root:

python3 -S -B source/qualify.py --out rerun-01

This command evolves fields. It uses only standard-library Python and leaves the original records/ untouched. The author qualification used an outside 120-second supervisor; the script has an internal 110-second elapsed-time check, a two-gigabyte starting free-space floor, and recorded checks for 128 MB peak resident memory and 15 MB new output. The resident-memory comparison uses macOS reporting units. Preserve the printed output and newly generated directory if any comparison fails; another platform's rounding or resource accounting deserves investigation before any threshold is changed.

Each completed case is saved before its comparisons. A failed comparison records its name, measured value, and limit; failure in the compressed receipt gives the exception. Read the first failed comparison first. A partial case can have saved frames even if it never returned a complete case row. Existing output directories are refused. A hash failure in the inspection commands stops inspection before numerical records are decoded.

The frozen calculations are ideal, incompressible, horizontally homogeneous linear channels with periodic height. They test growth, oscillatory limits, work, and resolution dependence over a short, small-amplitude interval. Their stress has no thermal-pressure normalization, so it is not an accretion $\alpha$ .

The browser Lab evaluates the continuum linear channel and exports its own sampled fields. Its initial amplitude, representative starting time, and amplitude stopping rule differ from this offline bundle, whose seed scale is fixed at $10^{-4}$ and whose cases end at time 3. Compare growth rate, stress ratio, or $E(t)/E(0)$ at matching parameters and elapsed times; the absolute energies need not match. The Lab does not import these offline frames. Its export and this download are two separately identified datasets.

Try two checks yourself

Predict the result of averaging each field before taking its product. At a quarter-wavelength, read the two velocities, multiply them, then compare with the volume average. Explain why these are different operations.

Next compare the two work estimates in field-measurements.json. Identify which one used every integration stage and which one used saved frames. Explain why changing the snapshot cadence could improve the latter without changing the evolved trajectory.

For another replay dataset, keep the same sequence: identify the exact files; read coordinates and units; calculate a diagnostic from its raw fields; choose a physically independent comparison; separate mesh error from time error; preserve the first failed result. That sequence gives each number a question it can answer.

Sources

The archived model notes identify Balbus and Hawley (1991), equations2.2 and2.5–2.9, and Stone and Gardiner (2010), SectionsII–III, as the physical and numerical references. Balbus–Hawley primary scan. Stone–Gardiner pinned version.

Try the prediction

Compare a Lab export with the replay

Predict which comparisons remain the same when the seed amplitude changes.
At Keplerian shear, compare the reference and fastest-growing tension settings; export each displayed linear state.
Compare growth rate, magnetic-to-velocity stress ratio, and energy normalized to that state's initial energy.
Explain why absolute Lab export energies need not equal the separately seeded offline records, and use the downloaded replay for spatial convergence.