Free boundary equilibrium calculations with ISOLVER

This option is selected with LEVGEO=12

The ISolver TRANSP equilibrium module is derived from Jon Menard’s original standalone IDL version of ISolver. It performs free boundary equilibrium calculations incorporating information on the evolving currents in PF coils and passive structures surrounding the plasma and optionally any feedback circuits controlling them, which allows it to evolve the $q$ profile self-consistently for predictive calculations. As of early 2023, TRANSP has ISolver device data available for 17 different tokamaks and spherical tori, including DIII-D, NSTX-U, SPARC, MAST, Asdex-U, EAST, KSTAR, and ITER. New devices can be added by providing information on coil geometry and composition, circuits, and feedback controllers in a specified format and pre-processing it through the make_eqtok utility; this task is normally performed by a TRANSP developer at users’ request.

Internally, ISolver uses Picard iteration to converge on a new equilibrium solution at each geometry time step. It begins with the existing toroidal current density distribution on an $(R, z)$ grid and the $q$ profile either input or, if NLISODIF=.TRUE., evolved according to the flux diffusion equation on a flux surface,

$$\left.\frac{\partial\psi}{\partial t}\right|_\phi = -\frac{V_{loop}}{2\pi} = -\eta_{\parallel} J^{OH} = -\eta_{\parallel} \frac{\left<\vec{J}\cdot\vec{B}\right> - \left<\vec{J}^{drive}\cdot\vec{B}\right>}{\left<\vec{B}\cdot\nabla\phi\right>} $$
where $\eta_{\parallel}$ is the plasma resistivity and $\left<\vec{J}^{drive}\cdot\vec{B}\right>$ represents driven currents.

Boundary conditions for this equation can be chosen either to match the plasma current (NLPCUR=.TRUE.) or the surface voltage (NLVSUR=.TRUE.) in the input data. Because of the relatively coarse resolution of the typical $(R, z)$ mesh, ISolver is poor at handling equilibria containing current sheets.

With $q$ updated, the toroidal field flux function can be evaluated according to
$$FV'\left<\frac{1}{R^2}\right> = 4\pi^2 q \Delta\psi $$
where $V' = dV/d\bar\psi $ is the differential volume and $\Delta\psi$ is the enclosed poloidal flux. Next the plasma current is updated according to
$$J_\phi^{p,F}(R,z)=\frac{R}{\Delta\psi}\frac{dp}{d\bar\psi} + \frac{F}{\mu_0 R \Delta\psi}\frac{dF}{d\bar\psi},$$
$$J_\phi^{new}(R,z)=\omega J_\phi^{p,F}(R,z) + (1-\omega)J_\phi^{old}(R,z),$$
where $\omega$ is the relaxation parameter for the iteration. This step is followed by solving for the plasma poloidal flux, given by the elliptic equation
$$\Delta^\star\psi^{plasma} = -\mu_0 R J_\phi^{new}(R,z)$$
The total poloidal flux is then the sum of this plasma flux and that due to the coils. The coil currents may be set from the input data, but this usually produces poor results. The default behavior is to determine the currents based on a least-squares fit of the resulting plasma boundary shape to a reference boundary provided. The free-boundary alternative is to solve the circuit equation to advance these currents. ISolver will search for the magnetic axis, X points, and flux surfaces based on the new plasma+coil flux solution, and will repeat the plasma current update step with successively smaller relaxation parameters until good surfaces are found.

Once an updated current distribution has been found to produce good flux surfaces, the solution is checked for convergence. If this has been been achieved, ISolver returns the updated equilibrium to TRANSP; if not, it loops back to the $q$ profile update step with the new current distribution and the iteration continues.