Cross-code validation of the equilibrium solver

The question

A numerical kernel that runs inside a web page is easy to distrust, and reasonably so. The question this study answers is narrow and falsifiable:

Does the browser fixed-boundary Grad–Shafranov solver produce the same flux field as an independently written solver, to what tolerance, and does that agreement survive a change of JavaScript engine?

Three things are measured. First, whether the discretisation converges at its designed order against an exact manufactured solution. Second, whether the converged field agrees node-for-node with a separate Python implementation of the same problem. Third, whether the whole benchmark suite reproduces when the engine underneath it changes — because a result that only holds on one runtime is not a result.

Method

Problem

The Grad–Shafranov operator is solved on a rectangular domain R ∈ [1.5, 2.5] m, Z ∈ [−0.5, 0.5] m with constant Dirichlet flux ψ = 0 on the complete boundary, discretised with second-order centred differences on an n × n grid. The cylindrical −ψ_R/R term is included explicitly rather than absorbed.

ψ_RR − ψ_R/R + ψ_ZZ = S(R, Z)

Two source terms

• Manufactured. An exact solution ψ = 0.1·sin(π(R−1.5))·sin(π(Z+0.5)) is chosen first and the source term derived from it analytically, so the discretisation error can be measured against a known answer rather than estimated.
• Profile polynomial. A pressure polynomial p′ and F² polynomial are prescribed and the resulting nonlinear problem is solved, exercising the same operator under the nonlinearity the laboratory actually uses.

The independent implementation

The reference is a separate Python program using SciPy sparse direct solves and a Newton iteration with an analytic Jacobian. It shares no code with the browser kernel, is written against the problem statement rather than against the JavaScript, and reads no JavaScript output — its first line records that constraint. The browser kernel uses nonlinear point relaxation, a different algorithm reaching the same discrete system.

Agreement is measured as the maximum absolute nodal difference in ψ over the whole grid, plus an exact-match test on the located magnetic axis. Both are compared against declared tolerances fixed before the comparison ran.

Result 1 · Convergence against the manufactured solution

Halving the grid spacing should reduce the error by a factor of four for a second-order scheme. It does.

observed order of accuracy = log₂(3.241862 × 10⁻⁴ / 8.093579 × 10⁻⁵) = 2.0020

The browser kernel is held to the same standard by benchmark NUM-11, whose acceptance interval for the observed order was fixed at [1.95, 2.05] before the run. Its measured relative residual on the fine grid is 4.616059 × 10⁻¹², against a tolerance of 1 × 10⁻¹¹, and its boundary error is exactly 0 Wb/rad.

Result 2 · Node-for-node agreement between the two codes

Benchmark NUM-12 compares the browser field against the Python field at every interior node, for all four cases. The declared tolerance is a maximum nodal difference below 2 × 10⁻¹¹ Wb/rad with an exact magnetic-axis match.

The largest disagreement anywhere in the comparison is 1.82 × 10⁻¹² Wb/rad, an order of magnitude inside the declared tolerance and comparable to the accumulated round-off of the two different algorithms. The located magnetic axis matches exactly in every case, at R = 2.0 m, Z = 0.0 m for the n = 33 manufactured grid.

Result 3 · Reproducibility across engines

A benchmark that only reproduces on the machine that produced it is not evidence. The full 14-case suite was regenerated on a different JavaScript engine — Node v22.22.2 against the v24.19.0 that produced the shipped record — with every source file verified byte-identical first.

Nine of fourteen benchmarks reproduce bit-identically. All fourteen still pass. Both equilibrium benchmarks — the subject of this study — are in the bit-identical set.

What the five drifting cases actually are

• NUM-05 and NUM-06 use adaptive step control along field lines. Different engine builds take marginally different branch decisions, so the step sequence differs. The quantity that moves by 3.7% is the reported error metric itself, which changes from 1.34 × 10⁻¹² to 1.29 × 10⁻¹². The acceptance tolerance is 1 × 10⁻⁷ — the drift is five orders of magnitude below the level that matters.
• NUM-13 localises Poincaré crossings. Its arclength error moves from 1.26 × 10⁻¹⁰ to 1.71 × 10⁻¹⁰ m against a 2 × 10⁻⁸ m tolerance, and the crossing count is unchanged at 3.
• NUM-08 moves in the ninth significant figure of a filament quadrature.
• NUM-25 differs by 57 bytes in an exported archive, because the embedded run timestamp has a different string length. The tool's own source comments anticipate exactly this. No numerical record changed: the canonical mismatch count is 0.

The strict --check mode fails on a different engine by design, because it hashes the complete report including the environment block. That is the correct behaviour for a release gate and the reason this study reports the per-record comparison instead.

The Python side reproduces exactly

The independent reference was also regenerated on newer libraries — numpy 2.4.4 and scipy 1.17.1 against the bundled numpy 2.3.5 and scipy 1.17.0. Every ψ value in all four cases came back bitwise identical, maximum |Δψ| = 0. Runtime was 0.7 s.

Limitations

These are the boundaries of the claim, stated so they can be checked rather than assumed.

• Fixed boundary only. The comparison covers the stated fixed rectangular boundary cases. Free-boundary equilibrium, where the plasma shape is determined by external coils, is not solved and not supported.
• Not a physical equilibrium. The manufactured case is an elliptic-kernel test with an invented source term. It verifies the operator and its discretisation, not that any resulting configuration is physically realisable, force-balanced in a real device, or MHD stable.
• Two implementations, one specification. Independent implementation guards against coding error. It does not guard against a shared misreading of the problem. A third implementation by a group that did not write either of these is the next meaningful test. Model validation of this kind is one of the areas open in the collaboration programme.
• Grid-extremum axis. The magnetic axis is located as a grid extremum with no subgrid critical-point refinement, so exact axis agreement is partly a consequence of both codes using the same grid.
• No experimental content. Nothing here involves measured plasma data. This is numerical verification, not physical validation.
• Engine drift is characterised, not eliminated. Five benchmarks are demonstrably engine-sensitive at levels far below tolerance. Anyone relying on bitwise reproducibility across runtimes should read Result 3 rather than assume it.

Code and data

Everything needed to reproduce this study ships inside the package you are reading. No download, account or network access is required.

• docs/observatory-equilibrium-reference/reference.py — the independent Python/SciPy solver, 2 433 bytes sha256 abaf9498447e02b0a8db8bb76ef3fe13afd5216ff9c44fc2ce937cb2df60d870
• docs/observatory-equilibrium-reference/reference.json — its recorded output, 79 697 bytes sha256 7946a6d17e90692128330cc725eba6bb2bbc04c5a5f637ce621d51913d8ccf9d
• docs/v748-verification/numerical-benchmarks.json — all 14 benchmark records, 83 830 bytes sha256 6da09fbbe0c1dd7d027ca37287d07d226a9d4bb7e669c43c71b18d540e0f0d4f
• tools/verify_observatory.mjs — the benchmark runner (Node built-ins only, no packages)
• torus-lab/observatory/equilibrium.mjs — the browser Grad–Shafranov kernel under test
• SHA256SUMS — hashes for every file in the package

The benchmark record carries its own sourceInventory: sixteen files with their SHA-256 hashes, recorded at the moment the benchmarks ran. Verifying those hashes proves the numbers on this page belong to the code in this package.

Reproduce it yourself

From the extracted package directory. Nothing writes outside it.

1 · Confirm the code under test is unmodified

Every file listed in the benchmark's source inventory should match its recorded hash.

sha256sum -c SHA256SUMS --quiet

2 · Re-run the independent Python solver

Requires Python 3 with numpy and scipy. It writes reference.json beside itself, so copy it elsewhere first if you want to keep the original for comparison.

cp docs/observatory-equilibrium-reference/reference.py /tmp/
cd /tmp && python3 reference.py

Compare the psi_WbPerRad arrays in your output against the bundled file. Result 3 above reports zero difference across a numpy and scipy version change.

3 · Re-run the browser kernel benchmarks

Requires Node.js. No npm install — the runner imports only Node built-ins and the package's own modules.

node tools/verify_observatory.mjs

This rewrites docs/v748-verification/numerical-benchmarks.json. Copy the original aside first, then compare the two files field by field to reproduce the table in Result 3. Expect the generatedAt, executionTimestamp and environment fields to differ.

4 · Strict gate

node tools/verify_observatory.mjs --check

Passes only on the recording environment. On any other engine it fails by design, which is what Result 3 measures.

5 · In the browser

Open Torus Lab and use the Observatory equilibrium controls to solve the same fixed-boundary problem interactively, then export the run package to inspect the numerical record directly.