Solving the Grad-Shafranov Equation Numerically

Compute a tokamak plasma equilibrium by iterating the Grad-Shafranov elliptic PDE on a grid until the flux converges.

The equation

Axisymmetric MHD equilibrium is governed by the Grad-Shafranov equation for the poloidal flux psi(R,Z): the elliptic operator Delta* psi = -mu0 R^2 dp/dpsi - F dF/dpsi, where p(psi) is pressure and F(psi) = R B_toroidal. The right-hand side depends on psi, so the problem is nonlinear and solved by iteration.

Discretization

On a uniform (R,Z) grid the operator Delta* psi = R d/dR( (1/R) dpsi/dR ) + d2psi/dZ2 becomes a five-point finite-difference stencil, with the 1/R factor breaking the usual symmetry. Boundary values of psi are fixed by the plasma-facing conductors.

Picard iteration

• Guess an initial psi and normalize it between axis and boundary.
• Evaluate the source term using the current psi through the profiles p'(psi) and FF'(psi).
• Solve the linear elliptic system for a new psi.
• Relax and repeat until psi stops changing.

import numpy as np
nR,nZ=65,65; R=np.linspace(1.0,3.0,nR); Z=np.linspace(-1,1,nZ)
dR=R[1]-R[0]; dZ=Z[1]-Z[0]
psi=np.zeros((nR,nZ))
def source(psi): return 1.0*np.ones_like(psi)  # placeholder p',FF'
for it in range(500):                # Gauss-Seidel sweep
    old=psi.copy()
    for i in range(1,nR-1):
        for j in range(1,nZ-1):
            rhs=-source(psi)[i,j]*R[i]**2
            psi[i,j]=0.25*(psi[i+1,j]+psi[i-1,j]+psi[i,j+1]+psi[i,j-1]-dR*dR*rhs)
    if np.max(np.abs(psi-old))<1e-6: break
print('converged in',it,'sweeps')

Reading the solution

Contours of psi are the nested flux surfaces; the innermost closed contour is the magnetic axis and the last closed one is the separatrix. The equilibrium fixes the plasma shape - including the negative triangularity used in the Hyperion breeder design - which then feeds stability and transport analysis. Real solvers add free-boundary coil currents and profile constraints, but the iterate-the-elliptic-solve loop shown here is the core.