Near the edge of tokamak plasmas, strong particle and energy sources and sinks (e.g. radiation, contact with a material wall,...) drive the plasma *away* from thermal equilibrium, thus invalidating the assumptions underpinning fluid theories.
The XGC, gyrokinetic, particle-in-cell (PIC) suite-of-codes, XGC1, XGCa and XGC0, were developed to provide a *comprehensive* description of kinetic transport phenomena in this complicated region; including: heating and cooling, radiation losses, neutral particle recycling and impurity transport.

A particle's “state” is described by the position, ${\bf x}$, and velocity, ${\bf v}$; constituting a a six-dimensional phase space. Gyrokinetic codes average over the very-fast, “gyromotion” of charged particles in strong magnetic fields, ${\bf B}$, and phase space becomes five-dimensional, ${\bf X} \equiv ({\bf x},v_\parallel,\mu)$, where ${\bf x}$ is the position of the “guiding center”, $v_\parallel$ is the velocity parallel to ${\bf B}$, and $\mu$ is the “magnetic moment”. The density of particles is given by a distribution function, $f({\bf X},t)$; the evolution of which, including collisons, is formally described by the “Vlasov-Maxwell” equation: $f$ evolves along “characteristics”, which are the dynamical trajectories of the guiding centers [1], \begin{eqnarray} \dot {\bf x} & = & \left[ v_\parallel {\bf b} + v_\parallel^2 \nabla B \times {\bf b} + {\bf B} \times ( \mu \nabla B - {\bf E}) / B^2 \right] / D,\\ \dot v_\parallel & = & - ( {\bf B} + v_\parallel \nabla B \times {\bf b} ) \cdot ( \mu \nabla B - {\bf E}), \end{eqnarray} where ${\bf E}$ is the electric field, ${\bf b}={\bf B}/|B|$, and $D \equiv 1 + v_\parallel {\bf b} \cdot \nabla \times {\bf b}/B$ ensures the conservation of phase-space volume (Liouville theorem). The non-thermal-equilibrium demands that gyrokinetic codes must evolve the full distribution function, by applying classical “full-f” [2,3] and noise-reducing, “total-f” techniques [4]. (This is in contrast to so-called “$\delta f$” methods, which evolve only a small perturbation to an assumed-static, usually-Maxwellian distribution.) As full-f codes, XGC can include heat and torque input, radiation cooling; and neutral particle recycling [5].

Multiple particle-species (e.g. ions and electrons, ions and impurities) are included; and XGC uses a field-alligned, unstructured mesh in cylindrical coordinates, and so can easily accommodate the irregular magnetic fields in the plasma edge (e.g. the “X-point”, “separatrix”).
XGC calculates transport in the *entire* plasma volume; from the “closed-flux-surface”, good-confinement region (near the magnetic axis) to the “scrape-off layer” (where magnetic fieldlines intersect the wall and confinement is lost).
Collisions between ions, electrons and impurities are evaluated using either (i) a linear, Monte-Carlo operator [6] for test-particles, and the Hirshman-Sigmar operator [7] for field-particle collisions; or (ii) a fully-nonlinear, Fokker-Planck-Landau collision operator [8,9].
XGC codes efficiently exploit massively-parallel computing architecture.

[1] Robert G. Littlejohn, Phys. Fluids

[2] C.S. Chang, Seunghoe Ku & H. Weitzner, Phys. Plasmas

[3] S. Ku, C.S. Chang & P.H. Diamond, Nucl. Fusion

[4] S. Ku, R. Hager

[5] D.P. Stotler, C.S. Chang

[6] Allen H. Boozer & Gioietta Kuo‐Petravic, Phys. Fluids

[7] S.P. Hirshman & D.J. Sigmar, Nucl. Fusion

[8] E. S. Yoon & C.S. Chang, Phys. Plasmas

[9] Robert Hager, E.S. Yoon

[10] Robert Hager & C.S. Chang, Phys. Plasmas

[11] D.J. Battaglia, K.H. Burrell

The accurate calculation of the equilibrium, stability and dynamical evolution of magnetically-confined plasma is fundamental for fusion research. The most suitable, macroscopic model to address some of the most critical challenges confronting tokamak plasmas is given by the extended-magnetohydrodynamic (MHD) equations, which describe plasmas as electrically conducting fluids of ions and electrons. The M3D-C$^1$ code [?] solves the fluid equations: for example, the “single-fluid” model, in which the ions and electrons are considered to have the same fluid velocity and temperature, the dynamical equations for the particle number density $n$, the fluid velocity $\vec{u}$, the total pressure $p$ are \begin{eqnarray} \frac{\partial n}{\partial t} + \nabla \cdot (n \vec{u}) & = & 0 \\ n m_i \left( \frac{\partial \vec{u}}{\partial t} + \vec{u} \cdot \nabla \vec{u} \right) & = & \vec{J} \times \vec{B} - \nabla p \color{blue}{ - \nabla \cdot \Pi + \vec{F}} \\ \frac{\partial p}{\partial t} + \vec{u} \cdot \nabla p + \Gamma p \nabla \cdot \vec{u} & = & \color{blue}{(\Gamma - 1) \left[Q - \nabla \cdot \vec{q} + \eta J^2 - \vec{u} \cdot \vec{F} - \Pi : \nabla u \right]} \end{eqnarray} together with a “generalized Ohm's law”, $\vec{E} = - \vec{u} \times \vec{B} \color{blue}{ + \eta \vec{J}}$; and with a reduced set of Maxwell's equations for the electrical-current density, $\vec{J} = \nabla \times \vec{B} / \mu_0$, and for the time evolution of the magnetic field, $\partial_t \vec{B} = - \nabla \times \vec{E}$. The manifestly divergence-free magnetic field, $\vec{B}$, and the fluid velocity, $\vec{u}$, are represented using stream functions and potentials, $\vec{B} = \nabla \psi \times \nabla \varphi + F \nabla \varphi$ and $\vec{u} = R^2 \nabla U \times \nabla \varphi + R^2 \omega \nabla \varphi + R^{-2} \nabla_\perp \chi$.

The “ideal-MHD” model is described by the above equations with the terms in blue set to zero. M3D-C$^1$ is also capable of computing the “two-fluid model”, which accommodates differences in the ion and electron fluid velocities: the generalized Ohm's law is augmented with $\color{red}{( \vec{J} \times \vec{B} - \nabla p_e - \nabla \cdot \Pi_e + \vec{F}_e ) / n_e}$, additional terms must be added to $\partial_t p$, and an equation for the “electron pressure”, $p_e$, must be included.

Physically-meaningful, “reduced” models provide accurate solutions in certain physical limits and are obtained, at a fraction of the computational cost, by restricting the scalar fields that are evolved: e.g., the two-field, reduced model is obtained by only evolving $\psi$ and $U$, and the “four-field”, reduced model by only evolving $\psi$, $U$, $F$, and $\omega$.

To obtain accurate solutions efficiently, over a broad range of temporal and spatial scales, M3D-C1 employs advanced numerical methods, such as: high-order, $C^1$, finite elements, an unstructured geometrical mesh, fully-implicit and semi-implicit time integration, physics-based preconditioning, domain-decomposition parallelization, and the use of scalable, parallel, sparse, linear algebra solvers.

M3D-C$^1$ is used to model numerous tokamak phenomena, including: edge localized modes (ELMs) [1]; sawtooth cycles [2]; and vertical displacement events (VDEs), resistive wall modes (RWMs), and perturbed (i.e. three-dimensional, 3D) equilibria [3].

[1] N.M. Ferraro, S.C. Jardin & P.B. Snyder, Phys. Plasmas

[2] S.C. Jardin, N. Ferraro

[3] N.M. Ferraro, T.E. Evans

Building on the theoretical foundations of Bruno & Laurence [1], that three-dimensional, (3D) magnetohydrodynamic (MHD) equilibria with “stepped”-pressure profiles are well-defined and guaranteed to exist, whereas 3D equilibria with integrable magnetic-fields and smooth pressure (or with non-integrable fields and continuous-but-fractal pressure) are pathological [2], Dr. S.R. Hudson wrote the “Stepped Pressure Equilibrium Code”, (SPEC) [3].
SPEC finds minimal-plasma-energy states, subject to the constraints of conserved helicity and fluxes in a collection, $i=1, N_V$, of nested sub-volumes, ${\cal R}_i$, by extremizing the multi-region, relaxed-MHD (MRxMHD), energy functional, ${\cal F}$, introduced by Dr. M.J. Hole, Dr. S.R. Hudson & Prof. R.L. Dewar [4,5],
\begin{eqnarray}
{\cal F} \equiv \sum_{i=1}^{N_V}
\left\{
\int_{{\cal R}_i} \! \left( \frac{p}{\gamma-1} + \frac{B^2}{2} \right) dv -
\frac{\mu_i}{2}
\left( \int_{{\cal R}_i} \!\! {\bf A}\cdot{\bf B} \, dv - H_i \right)
\right\}.
\end{eqnarray}
Relaxation is allowed in each ${\cal R}_i$, so unphysical, parallel currents are avoided and magnetic reconnection is allowed; and the ideal-MHD constraints are enforced at selected “ideal interfaces”, ${\cal I}_i$, on which ${\bf B}\cdot{\bf n}=0$.
The Euler-Lagrange equations derived from $\delta {\cal F}=0$ are: $\nabla \times {\bf B} = \mu_i {\bf B}$ in each ${\cal R}_i$;
and continuity of total pressure, $[[p+B^2/2]]=0$, across each ${\cal I}_i$, so that non-trivial, *stepped * pressure profiles may be sustained.
If $N_V=1$, MRxMHD equilibria reduce to so-called “Taylor states” [6]; and as $N_V \rightarrow \infty$, MRxMHD equilibria approach ideal-MHD equilibria [7].
Discontinuous solutions are admitted.
The figure (click to enlarge) shows an $N_V=4$ equilibrium with magnetic islands, chaotic fieldlines and a non-trivial pressure.

[1] Oscar P. Bruno & Peter Laurence, Commun. Pure Appl. Math.

[2] H. Grad, Phys. Fluids

[3] S.R. Hudson, R.L. Dewar

[4] M.J. Hole, S.R. Hudson & R.L. Dewar, J. Plas. Physics

[5] S.R. Hudson, M.J. Hole & R.L. Dewar, Phys. Plasmas

[6] J.B. Taylor, Rev. Modern Phys.

[7] G.R. Dennis, S.R. Hudson

Electrical currents flowing *inside* magnetically-confined plasmas cannot produce the magnetic field required for the confinement of the plasma itself.
Quite aside from understanding the physics of plasmas, the task of designing external, current-carrying coils that produce the confining magnetic field, ${\bf B}_C$, remains a fundamental problem; particularly for the geometrically-complicated, non-axisymmetric, “stellarator” class of confinement device [1].
The coils are subject to severe, engineering constraints: the coils must be “buildable”, and at a reasonable cost; and the coils must be supported against the forces they exert upon each other.
For reactor maintenance, the coils must allow access to internal structures, such as the vacuum vessel, and must allow room for diagnostics; many, closely-packed coils that give precise control over the external magnetic field might *not* be satisfactory.

The COILOPT code [2] and its successor, COILOPT++ [3], vary the geometrical degrees-of-freedom, ${\bf x}$, of a discrete set of coils to minimize a physics+engineering, “cost-function”, $\chi({\bf x})$, defined as the surface integral over a given, “target”, plasma boundary, ${\cal S}$, of the squared normal component of the “error” field,
\begin{eqnarray}
\chi^2 \equiv \oint_{\cal S} \left( \delta {\bf B} \cdot {\bf n} \right)^2 da + \mbox{engineering constraints},
\end{eqnarray}
where $\delta {\bf B} \equiv {\bf B}_C - {\bf B}_P$ is the difference between the externally-produced magnetic field (as computed using the Biot-Savart law [5]) and the magnetic field, ${\bf B}_P$, produced by the plasma currents (determined by an equilibrium calculation).
Using mathematical optimization algorithms, by finding an ${\bf x}$ that minimizes $\chi^2$ we find a coil configuration that minimizes the *total*, normal magnetic field at the boundary; thereby producing a “good flux surface”.

A new design methodology [3] has been developed for “modular” coils for so-called “quasi-axisymmetric” stellarators (stellarators for which the field strength *appears* axisymmetric in appropriate coordinates), leading to coils that are both simpler to construct and that allow easier access [4].
By straightening the coils on the outboard side of the device, additional space is created for insertion and removal of toroidal vessel segments, blanket modules, and so forth.
COILOPT++ employs a “spline+straight” representation, with fast, parallel optimization algorithms (e.g. the Levenberg-Marquardt algorithm [6], differential evolution [7], ...), to quickly generate coil designs that produce the target equilibrium.
COILOPT++ has allowed design of modular coils for a moderate-aspect-ratio, quasi-axisymmetric, stellarator reactor configuration, an example of which is shown in the figure (click to enlarge).

[1] Lyman Spitzer Jr., Phys. Fluids

[2] Dennis J. Strickler, Lee A. Berry & Steven P. Hirshman, Fusion Sci. & Technol.

[3] J.A. Breslau, N. Pomphrey

[4] George H. Neilson, David A. Gates

[5] Wikipedia, Biot-Savart Law

[6] Wikipedia, Levenberg–Marquardt Algorithm

[7] Wikipedia, Differential Evolution

One of the defining characteristics of the “stellarator” class of magnetic confinement device [1] is that the confining magnetic field is, for the most part, generated by external, current-carrying coils; stellarators are consequently more stable than their axisymmetric, “tokamak” cousins, for which an essential component of the confining field is produced by internal plasma currents.
Stellarators have historically had degraded confinement, as compared to tokamaks, and require more-complicated, “three-dimensional” geometry; however, by *exploiting* the three-dimensional shaping of the plasma boundary (which in turn determines the global, plasma equilibrium), stellarators may be designed to provide *optimized* plasma confinement.
This is certainly easier said-than-done: the plasma equilibrium is a *nonlinear* function of the boundary, and the particle and heat transport are nonlinear functions of the plasma equilibrium*!*

STELLOPT [2,3] is a versatile, *optimization* code that constructs suitable, magnetohydrodynamic equilibrium states via minimization of a “cost-function”, $\chi^2$, that quantifies how “attractive” an equilibrium is; for example, $\chi^2$ quantifies the stability of the plasma to small perturbations (including “ballooning” stability, “kink” stability) and the particle and heat transport (including “neoclassical” transport, “turbulent” transport, “energetic particle” confinement).
The independent, degrees-of-freedom, ${\bf x}$, describe the plasma boundary, which is input for the VMEC equilibrium code [4].
The boundary shape is varied (using either a Levenberg-Marquardt [5], Differential-Evolution [6] or Particle-Swarm [7] algorithms) to find minima of $\chi^2$; thus constructing an optimal, plasma equilibrium with both satisfactory stability and confinement properties.

A recent extension of STELLOPT enables “equilibrium reconstruction” [8,9]. By minimizing $\chi^2({\bf x})$ $\equiv$ $\sum_{i} [ y_i - f_i({\bf x}) ]^2 / \sigma_i^2$, where the $y_i$ and the $f_i({\bf x})$ are, respectively, the experimental measurements and calculated “synthetic diagnostics” (i.e. numerical calculations that mimic signals measured by magnetic sensors) of Thomson scattering, charge exchange, interferometry, Faraday rotation, motional Stark effect, and electro-cyclotron emission reflectometry, to name just a few, STELLOPT solves the “inverse” problem of inferring the plasma state given the experimental measurements. (The $\sigma_i$ are user-adjustable “weights”, which may reflect experimental uncertainties.) Equilibrium reconstruction is invaluable for understanding the properties of confined plasmas in present-day experimental devices, such as the DIIID device [10].

[1] Lyman Spitzer Jr., Phys. Fluids

[2] S.P. Hirshman, D.A. Spong

[3] D.A. Spong, S.P. Hirshman

[4] S.P. Hirshman, J.C. Whitson, Phys. Fluids

[5] Donald W. Marquardt, SIAM J. Appl. Math.

[6] D.E. Goldberg,

[7] Riccardo Poli, James Kennedy & Tim Blackwell, Swarm Intell.

[8] S. Lazerson & I.T. Chapman, Plasma Phys. Control. Fusion

[9] J.C. Schmitt, J. Bialek

[10] S. Lazerson S and the DIII-D Team, Nucl. Fusion