The question
Scientific software is only useful if its numbers can be trusted. This study asks a precise, testable question:
Does the browser equilibrium solver produce the same magnetic flux field as a separately written solver — and how closely?
Three things are measured: whether the solver converges at its designed rate towards an exact solution; whether it agrees point by point with a separate Python implementation; and whether the results reproduce identically on different software versions.
Method
The problem
The Grad–Shafranov equation is solved on a rectangular domain R ∈ [1.5, 2.5] m, Z ∈ [−0.5, 0.5] m with zero flux on the boundary, using second-order finite differences on an n × n grid. The cylindrical term −ψR/R is included explicitly.
ψRR − ψR/R + ψZZ = S(R, Z)
Two test cases
- Exact (manufactured) solution. An exact solution ψ = 0.1·sin(π(R−1.5))·sin(π(Z+0.5)) is chosen first and the matching source term derived from it, so the numerical error can be measured exactly.
- Pressure and current profiles. Polynomial pressure and current profiles are prescribed and the resulting nonlinear problem is solved — the same type of problem the laboratory solves interactively.
The reference code
The reference is a separate Python program using SciPy sparse solvers and a Newton iteration. It shares no code with the browser solver and was written from the mathematical problem statement, not from the browser code. The browser solver uses a different algorithm (nonlinear relaxation) to reach the same discrete solution, so agreement between them is a strong check on both.
Agreement is measured as the largest point-by-point difference in ψ over the whole grid, together with the position of the magnetic axis. Accuracy targets were fixed before any comparison was run.
Result 1 · Convergence to the exact solution
Halving the grid spacing should reduce the error by a factor of four for a second-order scheme. It does.
| Grid | Max residual (T) | Max |ψ − ψexact| (Wb/rad) | Ratio |
|---|---|---|---|
| n = 17 | 3.130829 × 10⁻¹⁴ | 3.241862 × 10⁻⁴ | — |
| n = 33 | 1.394440 × 10⁻¹³ | 8.093579 × 10⁻⁵ | 4.0055 |
observed order of accuracy = log₂(3.241862 × 10⁻⁴ / 8.093579 × 10⁻⁵) = 2.0020
The browser solver is held to the same standard. Its accuracy target — an observed order between 1.95 and 2.05 — was fixed before the run, and it measures 2.0020. Its relative residual on the fine grid (n = 33) is 9.03 × 10⁻¹², inside the convergence target of 1 × 10⁻¹¹, and the boundary condition is met exactly.
Result 2 · The two codes agree to round-off
The browser solution is compared with the Python solution at every interior grid point, for four cases. The accuracy targets were fixed in advance: a relative residual of at most 1 × 10⁻¹¹ and a maximum point-by-point difference below 2 × 10⁻¹¹ Wb/rad, with an exact match of the magnetic-axis position.
| Case | Grid | Reference residual (T) | Browser relative residual | Max difference Δψ (Wb/rad) | Result | Axis |
|---|---|---|---|---|---|---|
| Exact solution | n = 17 | 3.130829 × 10⁻¹⁴ | 4.616059 × 10⁻¹² | 8.100326 × 10⁻¹³ | pass | exact |
| Exact solution | n = 33 | 1.394440 × 10⁻¹³ | 9.026586 × 10⁻¹² | 1.821057 × 10⁻¹² | pass | exact |
| Pressure and current profiles | n = 17 | 1.075529 × 10⁻¹⁶ | 4.014248 × 10⁻¹² | 7.137086 × 10⁻¹⁵ | pass | exact |
| Pressure and current profiles | n = 33 | 5.342948 × 10⁻¹⁶ | 9.503926 × 10⁻¹² | 2.081885 × 10⁻¹⁴ | pass | exact |
The largest disagreement anywhere in the comparison is 1.82 × 10⁻¹² Wb/rad, eleven times smaller than the accuracy target and comparable to the rounding error 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 · Reproducible across software versions
The complete benchmark suite was re-run on two different Node.js versions and compared field by field. Both equilibrium benchmarks — the subject of this study — reproduced bit for bit, as did the large majority of the suite. The few benchmarks that showed tiny differences use adaptive step control, and their differences were orders of magnitude below their accuracy targets; every benchmark passed on both versions.
The Python reference was likewise re-run with newer versions of NumPy and SciPy. Every value in all four cases came back bitwise identical.
Scope of the result
- Fixed boundary. The study covers fixed-boundary equilibria. Free-boundary equilibria, where external coils determine the plasma shape, are a planned extension.
- Numerical verification. The study confirms that the equation is solved correctly. Comparison with measured plasma data is part of the validation work planned with research partners.
- Next step: external reproduction. Both codes were written within HBF. An independent reproduction by another research group is the natural next step and is one of the areas open in the collaboration programme.
- Grid-based axis. The magnetic axis is located on the grid without sub-grid refinement, so both codes report it at the same grid point.
Reproduce it yourself
Python reference
Download the reference solver and its recorded output. It needs Python 3 with NumPy and SciPy and runs in under a second.
reference.py— the Python/SciPy reference solverreference.json— its recorded output for all four cases
python3 reference.py
Compare the psi_WbPerRad arrays in the new output with the recorded file.
In the browser
Open Torus Lab and use the Observatory equilibrium controls to solve the same problem interactively, then export the run to inspect the full numerical record.
The complete benchmark suite and its records are available to research partners on request.
Compare with your code
Run it against your solver
If your group maintains an equilibrium code, we would welcome an independent comparison. The problem is fully specified above, and the reference files are small.
Hardip Bhesaniya · info@hbfnetwork.de