Checking growth, transport and shear work

The archive has SHA256 680ce195877487db5520cde1b4233fea89e91e46589daca562045feff24c4ef9. The preceding chapters define the local axes, velocity-unit magnetic disturbances, four channel equations, and centred spatial derivative. Here the density and box length are one in the saved units; the rotating cases use $\Omega=1/T_0$ . We inspect existing records and analytic limits. No new evolution is needed to follow the examples.

Begin with the saved fields and the original result

The archive's manifest.json identifies the other 299 members by byte count and SHA256. A matching fingerprint tells you that a file is the saved file; it does not by itself tell you that the physics is right. The saved records/receipt.json.gz contains eleven cases, forty-five predeclared checks, the original limits, actual values and status. Two of those checks bound resources. Read the actual numerical errors; the final status records which comparisons met their stated limits.

Inside each named case, numbered compressed JSON files contain the time, all four field arrays, energy and stress diagnostics, accumulated shear work, and four fundamental projections. case-returned.json.gz contains the returned case result before its comparisons. The recorded author interpreter was Python 3.14.3. The archive also provides a standard-library reader, source/inspect_saved.py, that verifies the existing members and recomputes energy and stress from all 277 field frames. That is an inspection, not a fresh simulation or a promise of bitwise equality on another computer.

The independent reader sums the actual arrays. This checks that the displayed diagnostics belong to their fields. It does not rerun the derivative or recover every intermediate time stage from the relatively sparse saved frames. Keep that distinction when following the work check below.

Separate spatial error from temporal error

First inspect spatial-64, spatial-128 and spatial-256. They use the same continuum growing seed, $K=0.6/T_0$ , final time $3T_0$ , and step $0.0025T_0$ . Their projected vector errors shrink by factors $3.999838$ and $3.999960$ when the spacing halves. That is the spatial comparison developed in the preceding chapter.

Now change the question: how accurately does the time integrator advance the already differenced spatial equations? Hold the mesh at32 samples. Its fundamental wave has

\begin{aligned} K_hT_0&=0.6\frac{\sin(2\pi/32)}{2\pi/32}\\ &\simeq0.596152111,\\ s_hT_0&\simeq0.661439531. \end{aligned}

The three temporal cases use the eigenvector and exact exponential of this grid equation as their reference. Comparing instead with the original continuum rate would leave a spatial error even as the time step approached zero. That remaining error could hide the time integrator's order. The returned referenceKind is semidiscrete-eigenmode for the three temporal cases and continuum for the other eight. expectedReferenceProjection holds the selected reference vector. Read those fields together; the setup in source/qualify.py specifies the initial vector and exponential.


Time step divided by the time unit	Steps to the common final time	Actual projected vector error
$0.02$	150	$5.00868381\times10^{-10}$
$0.01$	300	$3.14769716\times10^{-11}$
$0.005$	600	$1.97255936\times10^{-12}$

Dividing each coarse error by the next fine error gives $15.912216$ and $15.957427$ . In a fourth-order range, halving the step predicts a factor $2^4=16$ . The observed ratios support that interpretation. The smallest error is already close enough to floating-point effects that further refinement needs a fresh comparison; endless division by sixteen is not a measured result.

Actual time-refinement vector errors and separately accumulated stage-work residuals on one fixed mesh — Both panels use the same three 32-sample temporal cases. The upper points compare saved projections with the exact grid-mode reference; the dashed line is a fourth-order reference anchored to the coarsest result. The lower points compare independently calculated endpoint energy with accumulated stage-stress work, normalized by initial energy. The dotted lower line is the original work-check limit. Both errors are dimensionless; their different definitions set their numerical scales.

Worked example: what fourth order predicts

For a pure growing grid eigenmode, one classical fourth-order Runge–Kutta step multiplies its amplitude by

\begin{aligned} R(z)={}&1+z+\frac{z^2}{2}\\ &+\frac{z^3}{6}+\frac{z^4}{24},\\ z={}&s_h\Delta t. \end{aligned}

The exact multiplier is $e^z$ . Their first missing term is $z^5/120$ . Over a fixed duration, the number of steps grows as the inverse step size; the accumulated relative amplitude error is therefore of fourth order in the step size. The factor sixteen therefore follows from the method's expansion.

The coarsest measured error predicts about $3.1304\times10^{-11}$ after one halving. The actual value is $3.1477\times10^{-11}$ . A leading-order prediction need not match all higher-order terms. It does predict the large reduction and its scale. The comparison uses an analytic multiplier and the existing spatial records, with zero new evolution steps.

Follow the work independently of endpoint energy

Define the total outward stress of a saved field as $T_h=T_{{\rm Re},h}+T_{{\rm Max},h}$ . The continuous-time grid equations obey

\frac{dE_h}{dt}=q\Omega T_h.

The stored work was accumulated from the stress of the four Runge–Kutta stage fields, not assigned from the final energy difference. If $Y_1,\ldots,Y_4$ are those full arrays in one step, its contribution is

\begin{aligned} \Delta W={}&\frac{q\Omega\Delta t}{6}\\ &\times\bigl[T_h(Y_1)+2T_h(Y_2)\\ &+2T_h(Y_3)+T_h(Y_4)\bigr]. \end{aligned}

Energy per volume is read separately from the squared actual field samples. The comparison is independent of an endpoint subtraction used to define work, though both quantities use the same numerical field history. The original residual is

r_W=\frac{|E_h(t_f)-E_h(0)-W(t_f)|}{E_h(0)}.

For the coarse temporal case, the saved dimensionless energy and work are approximately

\begin{aligned} \widehat E_h(0)&\simeq1.10819394\times10^{-8},\\ \widehat E_h(t_f)&\simeq5.86371701\times10^{-7},\\ \widehat W(t_f)&\simeq5.75289762\times10^{-7}. \end{aligned}

The hat means division by the reference energy density $\rho_0L_0^2/T_0^2$ , with $\rho_0$ the chosen reference density. The unrounded saved numbers give $r_W=3.28304\times10^{-8}$ , below the predeclared $10^{-7}$ limit. Do not subtract the displayed rounded values and expect to recover the last digits of such a small residual. The fine temporal case gives $1.29096\times10^{-10}$ . These results test finite-step work consistency; they are not exact energy conservation to round-off.

Rigid rotation provides a different check. Its late saved stress is nonzero, but $q=0$ , so shear power and accumulated shear work vanish. The ordinary perturbation energy stays constant. The zero-field Keplerian test instead checks the epicyclic invariant $u_x^2+4u_y^2$ , not constant ordinary kinetic energy. Keeping those two limits separate helps expose a wrong rotational source term. Stone and Gardiner discuss artificial epicyclic growth in Section III of their pinned paper; their conserving time scheme is different from this channel's Runge–Kutta method.

Read zero divergence in the geometry actually tested

At uniform density, conversion from velocity units to physical magnetic units multiplies every component by a constant, so a zero divergence is unchanged. Magnetic divergence sums derivatives of matching components,

\begin{aligned} \nabla\cdot\mathbf B_{\rm phys} &=\partial_x B_{{\rm phys},x}\\ &+\partial_y B_{{\rm phys},y}\\ &+\partial_z B_{{\rm phys},z}. \end{aligned}

Here the horizontal components vary only with height, the horizontal directions are homogeneous, and the vertical background field is constant. Each derivative in that expression vanishes. A vertical variation of the radial field does not enter as a radial derivative. The audit copies this same channel into a periodic $3\times3\times N$ array and differences the represented field. The maximum is exactly zero in every tested case.

That result checks the representation of these assumptions and periodic indices. It does not exercise a general multidimensional induction equation or qualify constrained transport. Athena's treatment of shearing radial boundaries and magnetic flux is a separate problem described in Stone and Gardiner's Sections IV–V and the official manual. No such boundary remapping is tested by a horizontally homogeneous vertical channel.

Questions to carry to the next calculation

1. A finer time step stops reducing the error against the continuum wave. What should you test? Hold the mesh fixed and compare with its own exact grid-mode solution. A spatial error can remain after the temporal error becomes small. Then refine the mesh separately while holding the physical initial seed fixed.

2. Someone stores work by setting it equal to final energy minus initial energy. What would a tiny work residual demonstrate? Only the assignment itself. An informative comparison needs a separately accumulated source, here the stage stresses, and an independently evaluated energy from field samples.

3. A divergence check is zero for every saved channel. May it certify arbitrary three-dimensional MHD? No. Ask which derivatives were nonzero in the test. This construction has no horizontal gradients and no varying vertical magnetic component; it does not challenge the cancellations a general induction solver must preserve.

4. What does the complete original result let you reuse? Eleven specified linear periodic cases, their actual arrays, fixed comparisons and declared tolerances. To evaluate another calculation, record its units and geometry, identify its reference, vary space and time separately, compare transport with energy gain, and choose limits that remove the mechanism. Nonlinear saturation and a global accretion rate require additional models and evidence.

Sources

The primary local equations and linear growth relation are in Balbus and Hawley 1991, ApJ 376, 214–222, printed 216–218, equations 2.2, 2.5–2.9 and 2.13–2.14b. Primary scan. Stone and Gardiner, arXiv:1006.0139v1, Sections II–III, especially equations 11–18 and Section III.2, discuss local source terms, epicyclic energy and work; Sections IV–V treat the additional boundary and magnetic-flux problems. Pinned primary source. The official Athena++ shearing-box manual distinguishes radial shearing-periodic, orbital periodic and independently chosen vertical boundaries.

All numerical values here are observations of the pinned linear-channel-replay-public.zip, its original returned cases and 277 field frames. The checks described here reread saved fields and evaluate analytic expressions. They run no Athena calculation, add no evolution step and make no nonlinear-MHD acceptance claim.

Try the prediction

Separate time error from spatial error

Predict the error ratio when the time step halves on one fixed mesh.
Compare the three saved temporal cases against their own grid-mode reference.
Read the vector errors, endpoint energies and separately accumulated stage-stress work.
Explain the fourth-order ratios and identify what the homogeneous-channel divergence check tests.