Grad-Shafranov by finite differences, against Solov'ev's closed form

Adam Zettel · September 25, 2026 · written with an agent

Δ^{*} ψ \equiv R \frac{\partial}{\partial R} (\frac{1}{R} \frac{\partial ψ}{\partial R}) + \frac{\partial ^{2} ψ}{\partial Z ^{2}} = - μ_{0} R^{2} p^{'} (ψ) - F F^{'} (ψ)

A tokamak equilibrium solved on a grid and graded against an exact one. Solov'ev's $ψ = (R^{2} - 1)^{2} /8 + R^{2} Z^{2} / (2 κ^{2})$ solves the Grad-Shafranov equation with constant pressure gradient; the five-point solver, relaxed by SOR until the residual is below $1 0^{- 11}$ , misses it by $1.68 \times 1 0^{- 4}$ , $4.21 \times 1 0^{- 5}$ and $1.06 \times 1 0^{- 5}$ on grids of spacing 0.1, 0.05 and 0.025, an observed order of 1.995 and then 1.994. Those errors are scipy's direct solve of the same equations, to $1 0^{- 10}$ .

Open in the workbench Download .vivOutputs below are from the site build’s run; Run all runs it again in your browser.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

A tokamak holds its plasma in a magnetic field shaped like nested doughnuts. In equilibrium the magnetic force balances the pressure gradient, $J \times B = \nabla p$ . With the machine's rotational symmetry about its vertical axis, that balance reduces to one equation for one function, the poloidal flux $ψ (R, Z)$ , where $R$ is the distance from the axis of symmetry and $Z$ the height. The surfaces $ψ = const$ are the flux surfaces: the field lines lie in them, and the pressure is constant on each. The equation is the Grad-Shafranov equation above. Its right side holds two free functions of $ψ$ , the pressure $p$ and $F = R B_{ϕ}$ .

When $p^{'}$ and $F F^{'}$ are constants the equation is linear, and Solov'ev (1968) wrote down exact solutions. This post takes $F F^{'} = 0$ and $- μ_{0} p^{'} = 1 + 1/ κ^{2}$ (lengths in units of the major radius), for which

ψ (R, Z) = \frac{( R ^{2} - 1 ) ^{2}}{8} + \frac{R ^{2} Z ^{2}}{2 κ ^{2}}

is exact: its flux surfaces are nested about the magnetic axis at $(R, Z) = (1, 0)$ , and near the axis they are ellipses of elongation $κ = 1.7$ . posts/refs/grad-shafranov-solovev.py checks the closed form with sympy, then solves the same finite-difference equations directly with scipy's sparse LU, which shares no code or method with the SOR iteration here.

The discrete equation and SOR

On a grid with spacing $h_{R}$ in $R$ and $h_{Z}$ in $Z$ , the centered differences for $ψ_{R R}$ , $ψ_{R}$ and $ψ_{Z Z}$ give at each interior node

(\frac{1}{h _{R}^{2}} - \frac{1}{2 R h _{R}}) ψ_{i + 1, j} + (\frac{1}{h _{R}^{2}} + \frac{1}{2 R h _{R}}) ψ_{i - 1, j} + \frac{ψ _{i, j + 1} + ψ _{i, j - 1}}{h _{Z}^{2}} - (\frac{2}{h _{R}^{2}} + \frac{2}{h _{Z}^{2}}) ψ_{i, j} = (1 + κ^{- 2}) R_{i}^{2},

with $ψ$ fixed to the closed form on the edge of the box. The $- ψ_{R} / R$ term makes the equation unsymmetric in $R$ ; that is the only difference from Poisson's equation.

Successive over-relaxation (SOR) sweeps the grid, solves each node's equation for $ψ_{i, j}$ with its neighbors' newest values, and moves $ω$ times that far. The program takes the $ω$ that is optimal for Poisson's equation on the same box, $ω = 2/ (1 + 1 - ρ^{2})$ , with $ρ$ the Jacobi iteration's spectral radius $(cos (π / n_{R}) + cos (π / n_{Z})) /2$ for $n_{R}$ and $n_{Z}$ intervals of equal spacing. It sweeps until the largest residual is below $1 0^{- 11}$ , so what is left is the discretization error.

Three grids

The three grids have 8, 16 and 32 intervals across $R$ and 14, 28 and 56 across $Z$ , so $h_{R} = h_{Z} = h$ = 0.1, 0.05 and 0.025. (The reference script also solves a fourth grid, $h = 0.0125$ , where the order is 1.9998; it takes too long in a browser tab.)

The error of the difference equations is second order, and here it can be written down. The closed form's $R^{2} Z^{2}$ term is a quadratic in each variable, so the centered differences are exact on it. On $(R^{2} - 1)^{2} /8$ the second difference is off by $h^{2} ψ_{R R R R} /12 = h^{2} /4$ and the first difference over $R$ by $h^{2} ψ_{R R R} / (6 R) = h^{2} /2$ , so the difference operator applied to the exact $ψ$ misses the source by $- h^{2} /4$ , the same constant at every node. The error $e = ψ_{h} - ψ$ therefore solves the difference equation with source $h^{2} /4$ and zero boundary values: it is $h^{2}$ times a fixed function, up to that function's own discretization, and each halving of $h$ divides it by 4.

Printed

ints_r: [8, 16, 32]

ints_z: [14, 28, 56]

hs: [0.1, 0.05, 0.025]

halving_1: 3.98745

halving_2: 3.98383

errs: [0.000167718, 4.20615e-05, 1.05581e-05]

sweeps: [53, 108, 223]

resids: [2.65654e-12, 9.14469e-12, 9.86522e-12]

order_1: 1.99547

order_2: 1.99416

Grad-Shafranov error against the Solov'ev solution

five-point difference equations, solved by SORslope 2

Tests 5 of 5 pass

SOR brought every grid's residual below 1e-11 (passes)
the errors are scipy's direct solve of the same equations: 1.6771831e-4, 4.2061527e-5, 1.0558059e-5 (passes)
the three grids have 8, 16 and 32 intervals across R and 14, 28 and 56 across Z, so h = 0.1, 0.05 and 0.025 (passes)
each halving of h divides the error by 4, to within 0.02 (passes)
the observed order is 1.995 and then 1.994 (passes)

The flux surfaces, and the Shafranov shift

The heatmap is the solution on the 33 by 57 grid, with the flux surfaces drawn at five values of $ψ$ . They are nested about the magnetic axis at $(1, 0)$ , taller than wide near the axis ( $κ = 1.7$ ), and pushed toward small $R$ farther out.

On the midplane $Z = 0$ the closed form gives the crossings of a surface directly: $ψ = (R^{2} - 1)^{2} /8$ puts them at $R = 1 \pm 8 ψ$ . The surface through $R = 1.25$ crosses again at $R = 0.66144$ , so its center is at $0.95572$ : the magnetic axis sits $0.0443$ outboard of the center of this surface. That offset is the Shafranov shift, and it appears whenever there is pressure: the plasma pushes itself outward, where the field is weaker. The program reads the same crossing off the 33 by 57 grid by linear interpolation along the midplane and gets 0.661489, which is the reference script's interpolation on its own solution of the same grid, and $5 \times 1 0^{- 5}$ from the closed form.

Printed

r_in: 0.661489

center: 0.955745

57 × 33 grid

Tests 2 of 2 pass

the surface through R = 1.25 crosses the midplane again at 0.661489, the scipy solution's crossing, and within 1e-4 of the exact 0.661438 (passes)
its center is 0.0443 inboard of the magnetic axis at R = 1 (passes)