The Grad-Shafranov Solver

How the axisymmetric equilibrium equation is discretized and solved to find nested flux surfaces.

The Equation

The Grad-Shafranov equation is the force-balance condition for an axisymmetric toroidal plasma, a nonlinear elliptic PDE for the poloidal flux psi(R, Z). It reads Delta-star psi = -mu0 R^2 dp/dpsi - F dF/dpsi, where Delta-star is a modified Laplacian, p(psi) is the pressure profile and F(psi) = R B_t is the poloidal current function. The two free functions p(psi) and F(psi) are chosen inputs; psi is the unknown.

Discretization

A fixed-boundary solver imposes the plasma shape as the boundary where psi is constant and solves on a grid inside it. A free-boundary solver instead includes external coil currents and solves for the boundary self-consistently, requiring the plasma-vacuum interface to be found as part of the solution. Finite-difference and finite-element schemes both discretize Delta-star; the source term is nonlinear in psi through the profile functions.

Iterative Solution

Because the right side depends on psi, the equation is solved by iteration: guess psi, evaluate the source, invert the elliptic operator to get a new psi, and repeat under-relaxed until convergence. A Picard or Newton iteration handles the nonlinearity. Free-boundary solvers add an outer loop that adjusts to the coil currents and enforces the flux at the boundary and X-points.

import numpy as np
# One Picard sweep of a fixed-boundary Grad-Shafranov solve (schematic)
def gs_sweep(psi, R, dpsi_grid, invert_delta_star, p_prime, ff_prime):
    # source S = -mu0 R^2 p'(psi) - F F'(psi)
    mu0 = 4e-7 * np.pi
    S = -mu0 * R**2 * p_prime(psi) - ff_prime(psi)
    psi_new = invert_delta_star(S)   # solve Delta* psi_new = S with fixed BC
    return 0.5 * psi + 0.5 * psi_new  # under-relaxed update

Relevance

Every tokamak equilibrium, stability input, and transport grid begins with a Grad-Shafranov solution; it produces the flux surfaces, the Shafranov shift, and the q profile. For the Hyperion breeder concept, the equilibrium reconstruction underlying its design point, including negative triangularity of -0.30, is a Grad-Shafranov solution for a machine in simulation, not a fit to measured hardware.