PPPL

Research Codes

Overview

A peeling-ballooning eigenfunction calculated using M3D-C$^1$. The red and blue show the perturbed pressure, and the magenta curve the location of the last closed flux surface of the plasma equilibrium. The green triangles show the finite element mesh, which has resolution packed near the edge of the plasma to efficiently resolve the eigenfunction. The peeling-ballooning instability leads to Edge Localized Modes (ELMs) in tokamaks.

The following research codes are actively maintained:

The Department of Energy (DOE) Princeton Plasma Physics Laboratory (PPPL) Theory Department Code Licence Release form is here.

XGC: the X-point Gyrokinetic Code for Transport in Tokamaks

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 28, 2015 (1985)
[2] C.S. Chang, Seunghoe Ku & H. Weitzner, Phys. Plasmas 11, 2649 (2004)
[3] S. Ku, C.S. Chang & P.H. Diamond, Nucl. Fusion 49, 115021 (2009)
[4] S. Ku, R. Hager et al., J. Comp. Phys. 315, 467 (2016)
[5] D.P. Stotler, C.S. Chang et al., J. Nucl. Mater. 438 S1275 (2013)
[6] Allen H. Boozer & Gioietta Kuo‐Petravic, Phys. Fluids 24, 851 (1981)
[7] S.P. Hirshman & D.J. Sigmar, Nucl. Fusion 21, 1079 (1981)
[8] E. S. Yoon & C.S. Chang, Phys. Plasmas 21, 032503 (2014)
[9] Robert Hager, E.S. Yoon et al., J. Comp. Phys. 315 644 (2016)
[10] Robert Hager & C.S. Chang, Phys. Plasmas 23, 042503 (2016)
[11] D.J. Battaglia, K.H. Burrell et al., Phys. Plasmas 21, 072508 (2014)
[#h34: C-S. Chang, R. Hager, S-H. Ku, 2016-08-05]

Gkeyll: Energy-conserving, discontinuous, high-order discretizations for gyro-kinetic simulations

Discontinuities at cell boundaries are allowed and used to compute a “numerical flux” needed to update the solution. Shown is a DG fit (in the least-square sense) of $x^4+\sin(5x)$ onto constant (left), linear (middle) and quadratic (right) basis functions.

Fusion energy gain in tokamaks depends sensitively on the plasma edge [1,2]; but, because of open- and closed-fieldlines, interaction with divertor plates, neutral particles, and large electromagnetic fluctuations, the edge region is, understandably, difficult to treat analytically. Large-scale, kinetic, numerical simulations are required. The “GKEYLL” code [3] is a flexible, robust, powerful, algorithm that provides numerical calculations of gyrokinetic turbulence, which importantly preserves the conservation laws of gyrokinetics.

Gyrokinetics [4] describes how a distribution of particles, described by a density-distribution function, $f(t,{\bf z})$, evolves in time, $t$; where ${\bf z}=({\bf x},{\bf v})$ describes the position and velocity, i.e. a point in “phase space”, of the guiding-center. Elegant theoretical and numerical descriptions of this dynamical system exploit the “Hamiltonian properties”, i.e. $\partial f / \partial t + \{f,H\} = 0$, where $H({\bf z})$ is the Hamiltonian (i.e. “energy”, e.g. for a Vlasov system, $H(x, v, t)$ $=$ $mv^2/2$ $+$ $q \, \phi(x,t)$, where $\phi$ is the electro-static potential) and $\{g,f\}$ is the Poisson bracket operator [5]. For reliable simulations, numerical discretizations must preserve so-called “quadratic invariants”, $\int H \{ f,H\} \, d{\bf z} $ $=$ $ \int f \{ f,H\} \, d{\bf z} = 0$.

Discontinuous Galerkin (DG) algorithms represent the “state-of-art” discretizations of hyperbolic, partial-differential equations [6]. DG combines the key advantages of finite-elements (e.g. low phase error, high accuracy, flexible geometries) with those of finite-volumes (e.g. up-winding, guaranteed positivity/monotonicity, . . . ), and makes efficient use of parallel computing architectures. DG is inherently super-convergent: e.g., whereas finite-volume methods interpolate $p$ points to get $p$-th order accuracy, DG methods in contrast interpolate $p$ points to get $(2p−1)$-th order accuracy, a significant advantage for $p>1$! Use of DG schemes may lead to significant advances towards a production-quality, edge gyrokinetic simulation software with reasonable computational costs.


[1] A. Loarte, et al. Nucl. Fusion 47, S203 (2007)
[2] P.C. Stangeby, The Plasma Boundary of Magnetic Fusion Devices, Institute of Physics Publishing (2000).
[3] A. Hakim, 26th IAEA Fusion Energy Conference (2016)
[4] E.A. Frieman & Liu Chen, Phys. Fluids 25, 502 (1982)
[5] John R. Cary & Alain J. Brizard, Rev. Mod. Phys. 81, 693 (2009)
[6] Bernardo Cockburn & Chi-Wang Shu, J. Comput. Phys. 141, 199 (1998)
[#h44: A. Hakim, 2016-08-05]

NOVA and NOVA-K

Hot beam ion orbits in the Tokamak Fusion Test Reactor (TFTF) shot #103101 [2,4]. Shown (click to enlarge) are passing (a) and trapped (b) orbits represented in both $\psi,R$ and $Z,R$ planes at a time when TAEs were observed.

NOVA is a suite of codes including the linear ideal eigenmode solver to find the solutions of ideal magnetohydrodynamic (MHD) system of equations [1], including such effects as plasma compressibility and realistic tokamak geometry. The kinetic post-processor of the suite, NOVA-K [2], analyses those solutions to find their stability properties. NOVA-K evaluates the wave particle interaction of the eigenmodes such as Toroidal Alfvén Eigenmodes (TAEs) or Reversed Shear Alfvén Eigenmodes (RSAEs) by employing the quadratic form with the perturbed distribution function coming from the drift kinetic equations [2,4]. The hot particle growth rate of ideal MHD eigenmode is expressed via the equation \begin{eqnarray} \frac{\gamma_{h}}{\omega_{AE}}=\frac{\Im\delta W_{k}}{2\delta K}, \end{eqnarray} where $\omega_{AE}$ is Alfvén eigenmode frequency, $\delta W_{k}$ is the potential energy of the mode, and $\delta K$ is the inertial energy. In computations of the hot ion contribution to the growth rate NOVA-K makes use of the constant of motion to represent the hot particle drift orbits. Example of such representation is shown in the figure.

NOVA is routinely used for AE structure computations and comparisons with the experimentally observed instabilities [5,6]. The main limitations of NOVA code are caused by neglecting thermal ion finite Larmor radius (FLR), toroidal rotation, and drift effects for the eigenmode computations. Therefore it can not describe accurately radiative damping for example. Finite element methods are used in radial direction and Fourier harmonics are used in poloidal and toroidal directions.

NOVA-K is able to predict various kinetic growth and damping rates perturbatively, such as the phase space gradient drive from energetic particles, continuum damping, radiative damping, ion/electron Landau damping and trapped electron collisional damping.

More information can be found on Dr. N. N. Gorelenkov's NOVA page.


[1] C. Z. Cheng, M. C. Chance, Phys. Fluids 29, 3695 (1986)
[2] C. Z. Cheng, Phys. Reports 211, 1 (1992)
[3] G. Y. Fu, C. Z. Cheng, K. L. Wong, Phys. Fluids B 5, 4040 (1993)
[4] N. N. Gorelenkov, C. Z. Cheng, G. Y. Fu, Phys. Plasmas 6, 2802 (1999)
[5] M. A. Van Zeeland, G. J. Kramer et al., Phys. Rev. Lett. 97, 135001 (2006)
[6] G. J. Kramer, R. Nazikian et al., Phys. Plasmas 13, 056104 (2006)
[#h45: N. N. Gorelenkov, 2018-07-05]

Orbit

The guiding center code Orbit, using guiding center equations first derived in White and Chance [1] and more completely in The Theory of Toroidally Confined Plasmas [2], can use equal arc, PEST, or any other straight field line coordinates $\psi_p, \theta, \zeta$, which along with the parallel velocity and energy completely specify the particle position and velocity. The guiding center equations depend only on the magnitude of the $B$ field and functions $I$ and $g$ with $\vec{B} = g\nabla \zeta + I\nabla \theta + \delta \nabla \psi_p$. In simplest form, in an axisymmetric configuration, in straight field line coordinates, and without field perturbation or magnetic field ripple, the equations are \begin{eqnarray} \dot{\theta} = \frac{ \rho_{\parallel}B^2}{D} (1 - \rho_{\parallel}g^{\prime}) + \frac{g}{D}\left[ (\mu + \rho_{\parallel}^2B)\frac{\partial B}{\partial\psi_p} + \frac{\partial \Phi}{\partial\psi_p}\right], \label{tdot2} \end{eqnarray} \begin{eqnarray} \dot{\psi_p} = -\frac{g}{D}\left[(\mu +\rho_{\parallel}^2B)\frac{\partial B}{\partial\theta} + \frac{\partial \Phi}{\partial\theta}\right], \label{psdot2} \end{eqnarray} \begin{eqnarray} \dot{\rho_{\parallel}} = -\frac{(1 -\rho_{\parallel}g^{\prime})}{D} \left[(\mu +\rho_{\parallel}^2B)\frac{\partial B}{\partial\theta} + \frac{\partial\Phi}{\partial\theta}\right], \label{rdot2} \end{eqnarray} \begin{eqnarray} \dot{\zeta} = \frac{ \rho_{\parallel}B^2}{D} (q +\rho_{\parallel}I^{\prime}_{\psi_p}) - \frac{I}{D} \left[(\mu +\rho_{\parallel}^2B)\frac{\partial B}{\partial\psi_p} + \frac{\partial \Phi}{\partial\psi_p}\right], \label{zdot2} \end{eqnarray} with $ D = gq + I + \rho_\parallel(g I'_{\psi_p} - I g'_{\psi_p})$. The terms in $\partial \Phi/\partial\psi_p$, $\partial \Phi/\partial\theta$, $\partial \Phi/\partial\zeta$ are easily recognized as describing $\vec{E}\times \vec{B}$ drift.

Orbit can read numerical equilibria developed by TRANSP or other routines, particle distributions produced by TRANSP, and mode eigenfuncitons produced by NOVA.

The code uses a fourth order Runge Cutta integration routine. It is divided into a main program Orbit.F, which is essentially a heavily commented name list and a set of switches for choosing the type of run, the diagnostics, the data storage, and output.

The code has been used extensively since its inception at PPPL, General Atomics, RFX (Padova) and Kiev for the analysis of mode-induced particle loss for fishbones and TAE modes, induced ripple loss, the modification of particle distributions, local particle transport, etc.

To obtain a copy take all files from /u/ftp/pub/white/Orbit. To submit a job modify the script batch, modify Orbit.F to choose the run desired, type make, and type qsub -q sque batch


[1] R. B. White & M. S. Chance, Phys. Fluids 27, 2455 (1984)
[1] R. B. White, The theory of Toroidally Confined Plasmas, Imperial College Press, 3rd. ed. (2014)
[#h46: R. B. White, 2018-07-05]

GTS: Gyrokinetic Tokamak Simulation Code

Snapshot of a GTS ion-temperature-gradient instability simulation showing the field-line-following mesh and the quasi-2D structure of the electrostatic potential in the presence of microturbulence. Notice the fine structure on the two poloidal planes perpendicular to the toroidal direction, while the potential along the field lines changes very little as you go around the torus. (Image generated from a GTS simulation by Kwan-Liu Ma and his group at the University of California, Davis.)

The Gyrokinetic Tokamak Simulation (GTS) code is a full-geometry, massively parallel $\delta f$ particle-in-cell code [1,2] developed at PPPL to simulate plasma turbulence and transport in practical fusion experiments. The GTS code solves the modern gyrokinetic equation in conservative form [3]: \begin{equation} \frac{\partial f_a}{\partial t}+\frac{1}{B^{*}}\nabla_{\bf Z}\cdot (\dot{{\bf Z}}B^{*}f_a)=\sum\limits_{b} C[f_a, f_b]. \end{equation} for a gyro-center distribution function $f ({\bf Z}, t)$ in 5-dimension phase space ${\bf Z}$, along with the gyrokinetic Poisson equation and Ampere's law for potentials using a $\delta f$ method.

GTS features high robustness at treating globally consistent, shaped cross-section tokamaks; in particular, the highly challenging spherical tokamak geometry such as NSTX and its upgrade NSTX-U. GTS simulations directly import plasma profiles of temperature, density and toroidal rotation from the TRANSP experimental database, along with the related numerical MHD equilibria, including perturbed equilibria. General magnetic coordinates and a field-line-following mesh are employed [1]. The particle gyro-center motion is calculated by Lagrangian equations in the flux coordinates, which allows for accurate particle orbit integration thanks to the separation between fast parallel motion and slow perpendicular drifts. The field-line-following mesh accounts for the nature of the quasi-2D mode structure of drift-wave turbulence in toroidal systems, and hence offers a highly efficient spatial resolution for strongly anisotropic fluctuations in fusion plasmas. Fully-kinetic electron physics is included. In particular, both trapped and untrapped electrons are included in the non-adiabatic response. GTS solves the field equations in configuration space for the turbulence potentials using finite element method on unstructured mesh, which is carried out by PETSc. The real space, global field solver, in principle, retains all toroidal modes from $(m,n) = (0,0)$ all the way up to a limit which is set by grid resolution, and therefore retains full-channel nonlinear energy couplings.

One remarkable feature in GTS, which distinguishes it from the other $\delta f$ particle simulations, is that the weight equations satisfy the incompressibility condition in extended phase space $({\bf Z}, w)$ [4]. Satisfying the incompressibility is actually required in order to correctly solve the $\delta f$ kinetic equation using simulation markers whose distribution function $F({\bf Z}, w, t)$ is advanced along the marker trajectory in the extended phase space according to $F({\bf Z},w,t)=const.$.

In GTS, Coulomb collisions between like particles are implemented via a linearized Fokker-Plank operator with particle, momentum and energy conservation. Electron-ion collisions are simulated by the Lorentz operator. By modeling the same gyrokinetic-Poisson system, GTS is extended to performing global neoclassical simulation in additional to traditional turbulence simulation. More importantly, GTS now is able to do a global gyrokinetic simulation with self-consistent turbulence and neoclassical dynamics coupled together. This remarkable capability is shown to lead to significant new features regarding nonlinear turbulence dynamics, impacting a number of important transport issues in tokamak plasmas. In particular, this capability is critical for the proposed study of nonlinear NTM physics. For example, it allows to calculate bootstrap current in the presence of turbulence [5, 6], which plays a key role for NTM evolution. Currently, GTS has been extended to simulating finite beta physics. This includes low-n shear-Alfvén modes, current driven tearing modes, kinetic ballooning modes and micro-tearing modes.

A state-of-the-art electromagnetic algorithm has been developed and implemented into GTS [7,8] with the goal to achieve a robust, global electromagnetic simulation capability to attack the highly challenging electron transport problem in high-$\beta$ NSTX-U plasmas and be used as a first-principles-based module for integrated whole device modeling of turbulence/neoclassical/MHD physics.

On the front of physics study, GTS simulations have been applied to wide experiments for various problems. Recent applications include discovering new turbulence sources responsible for plasma transport and understanding underlying physics behind the confinement scaling in spherical tokamak experiments [4,9], validating the physics of turbulence-driven plasma flow and first-principles-based model prediction of intrinsic rotation profile against profile against experiments [10], and plasma self-generated non-inductive current in turbulent fusion plasmas [5].

GTS has three levels of parallelism: a one-dimensional domain decomposition in the toroidal direction, dividing both the grid and the particles, a particle distribution within each domain, which further divides the particles between processors, and a loop-level multi-threading method. The domain decomposition and the particle distribution are implemented with MPI, while the loop-level multi-threading is implemented with OpenMP directives.


[1] W. X. Wang, Z. Lin et al., Phys. Plasmas 13, 092505 (2006)
[2] W. X. Wang, P. H. Diamond et al., Phys. Plasmas 17, 072511 (2010)
[3] A. J. Brizard & T. S. Hahm, Rev. Mod. Phys. 79, 421 (2007)
[4] W. X. Wang, S. Ethier et al., Phys. Plasmas 22, 102509 (2015)
[5] W. X. Wang et al., Proc. 24th IAEA Fusion Energy Conference, (2012), TH/P7-14
[6] T. S. Hahm, Nucl. Fusion 53 104026 (2013)
[7] E. A. Startsev & W. W. Lee, Phys. Plasmas 21, 022505 (2014)
[8] E. A. Startsev et al., Paper BM9.00002, APS-DPP Conference, San Jose, CA (2016)
[9] W. X. Wang, S. Ethier et al., Nucl. Fusion 55, 122001 (2015)
[10] B. A. Grierson, W. X. Wang et al., Phys. Rev. Lett. 118, 015002 (2017)
[#h47: W. Wang, 2018-07-11]

HMHD: A 3D Extended MHD Code

A snapshot of plasmoid-mediated turbulent reconnection simulation showing the parallel current density and samples of magnetic field lines.

HMHD is a massively parallel, general purpose 3D extended MHD code that solves the fluid equations of particle density $n$ and momentum density $n\boldsymbol{u}$: \[ \partial_{t}n+\nabla\cdot\left(n\boldsymbol{u}\right)=0, \] \[ \partial_{t}\left(n\boldsymbol{u}\right)+\nabla\cdot\left(n\boldsymbol{uu}\right)=-\nabla\left(p_{i}+p_{e}\right)+\boldsymbol{J}\times\boldsymbol{B}-\nabla\cdot\boldsymbol{\Pi}+\boldsymbol{F}, \] where $\boldsymbol{J}=\nabla\times\boldsymbol{B}$ is the electric current, $p_{e}$ and $p_{i}$ are the electron and ion pressures, $\boldsymbol{\Pi}$ is the viscous stress tensor, and $\boldsymbol{F}$ is an external force. The magnetic field ${\bf B}$ is stepped by the Faraday's law \[ \partial_{t}\boldsymbol{B}=-\nabla\times\boldsymbol{E}, \] where the electric field $\boldsymbol{E}$ is determined by a generalized Ohm's law that incorporates the Hall term and the electron pressure term in the following form: \[ \boldsymbol{E}=-\boldsymbol{u}\times\boldsymbol{B}+d_{i}\frac{\boldsymbol{J}\times\boldsymbol{B}-\nabla p_{e}}{n}+\eta\boldsymbol{J}, \] with $\eta$ the resistivity and $d_{i}$ the ion skin depth. The set of equations is completed by additional equations for electron and ion pressures $p_{e}$ and $p_{i}$, where several options are available. In the simplest level an isothermal equation of state $p_{i}=p_{e}=nT$ is assumed; in a sophisticated level, the ion and electron pressure can be evolved individually with viscous and resistive heating, anisotropic thermal conduction, and thermal exchange between the two species; various additional options are available between these two levels. HMHD is flexible to switch on and off various effects in the governing equations.

HMHD uses a single Cartesian grid, with the capability of variable grid spacing. The numerical algorithm [1] approximates spatial derivatives by finite differences with a five-point stencil in each direction. The time-stepping scheme has several options including a second-order accurate trapezoidal leapfrog method as well as three-stage or four-stage strong stability preserving Runge-Kutta methods [2,3]. HMHD is written in Fortran 90 and parallelized with MPI for domain decomposition, augmented with OpenMP for multi-threading in each domain. HMHD has been employed to carry out large-scale 2D simulations plasmoid-mediated reconnection in resistive MHD [4,5,6] and Hall MHD [7,8]. It has been used to carry out 3D self-generated turbulent reconnection simulations [9,10].


[1] P. N. Guzdar, J. F. Drake et al., Phys. Fluids B 5, 3712 (1993)
[2] S. Gottlieb, C.-W. Shu & E. Tadmor, SIAM Review 43, 89 (2001)
[3] R. J. Spiteri & S. J. Ruuth, SIAM J. Numer. Anal. 40, 469 (2002)
[4] Y.-M. Huang & A. Bhattacharjee, Phys. Plasmas 17, 062104 (2010)
[5] Y.-M. Huang & A. Bhattacharjee, Phys. Rev. Lett. 109, 265002 (2012)
[6] Y.-M. Huang, L. Comisso & A. Bhattacharjee, Astrophys. J. 849, 75 (2017)
[7] Y.-M. Huang, A. Bhattacharjee & B. P. Sullivan, Phys. Plasmas 18, 072109 (2011)
[8] J. Ng, Y.-M. Huang et al., Phys. Plasmas 22, 112104 (2015)
[9] Y.-M. Huang & A. Bhattacharjee, Astrophys. J. 818, 20 (2016)
[10] D. Hu, A. Bhattacharjee & Y.-M. Huang, Phys. Plasmas 25, 062305 (2018)
[#h48: Y-M. Huang, 2018-07-11]

FOCUS: Flexible Optimized Coils Using Space curves

Modular coils for a conventional rotating elliptical stellarator. The color on the internal plasma boundary indicates the strength of mean curvature.

Finding an easy-to-build coils set has been a critical issue for stellarator design for decades. The construction of the coils is only one component of modern fusion experiments; but, realizing that it is the currents in the coils that produce the “magnetic bottle” that confines the plasma, it is easy to understand that designing and accurately constructing suitable coils is paramount.

Conventional approaches simplify this problem by assuming that coils are lying on a defined toroidal “winding” surface [1]. The Flexible Optimized Coils using Space Curves (FOCUS) code [2] represents coils as one-dimensional curves embedded in three-dimensional space. A curve is described directly, and completely generally, in Cartesian coordinates as ${\bf x}(t) = x(t) \, {\bf i} + y(t) \, {\bf j} + z(t) \, {\bf k}$. The coil parameters, ${\bf X}$, are then varied to minimize a target function consisting of multiple objective functions, \begin{equation} \chi^2({\bf X}) \equiv w_{B}\int_S \frac{1}{2} \left ( { \bf B} \cdot {\bf n} \right )^2 ds \ + \ w_{\Psi}\int_0^{2\pi} \frac{1}{2} \left( \Psi_\zeta \; - \; \Psi_o \right)^2 d\zeta \ + \ w_{L} \; \sum_{i=1}^{N_C} L_i + \cdots \end{equation} These objective functions cover both “physics” requirements and “engineering” constraints, such as the normal magnetic field, the toroidal flux, the resonant magnetic harmonics, coil length, coil-to-coil separation and coil-to-plasma separation.

With analytically calculated derivatives, FOCUS computes the gradient and Hessian with respect to coil parameters fast and accurately. It allows FOCUS to employ numerous powerful optimization algorithms, like the Newton method [3]. The figure (click to enlarge) shows using FOCUS to design modular coils for a rotating elliptical stellerator. FOCUS is also applied to analyze the error field sensitivity to coil deviations [4], vary the shape of plasma surface in order to simplify the coil geometry [5] and design non-axisymmetric RMP coils for tokamaks.


[1] P. Merkel, Nucl. Fusion, 27, 867 (1987)
[2] Caoxiang Zhu, Stuart R. Hudson et al., Nucl. Fusion, 58, 016008 (2017)
[3] Caoxiang Zhu, Stuart R. Hudson et al., Plasma Phys. Control. Fusion 60, 065008 (2018)
[4] Caoxiang Zhu, Stuart R. Hudson et al., Plasma Phys. Control. Fusion 60, 054016 (2018)
[5] S. R. Hudson, C. Zhu et al., Phys. Lett. A, in press
[#h49: C. Zhu, 祝曹祥, 2016-08-05]

RBQ1D

(a) Mode evolution for different levels of collisionality featuring intermitency and steady saturation within RBQ1D depending on the effective collisionality and (b) scaling of saturation amplitude as a function of collisionality within RBQ1D. Calculations are done for the global TAE shown in panel (c).
The interaction between fast ions and Alfvénic eigenmodes has proved to be numerically expensive to be modeled in realistic large tokamak configurations such as ITER. Therefore, it is motivating to explore reduced, numerically efficient models such as the Resonance Broadened Quasilinear model, used to build the code in its one-dimensional version (RBQ1D). The code is capable of modeling the fast ion distribution function in the direction of the canonical toroidal momentum while self-consistently evolving the amplitude of interacting Alfvénic modes. The theoretical approach is based on the model proposed by Berk et al [1] to address the resonant energetic particle interaction in both regimes of isolated and overlapping modes. RBQ1D is written by using the same structure of the conventional quasilinear equations for the fast ion distribution function but with the resonant delta function broadened primarily in the radial direction. The model was designed to reproduce the expected saturation level for non-overlapping modes from analytic theory. In the model, the nonlinear trapping (bounce) frequency is the fundamental variable for the dynamics in the vicinity at a resonance. The diffusion equation of RBQ1D, derived using action and angle variables, is [2,3] \begin{eqnarray} \frac{\partial f}{\partial t}=\pi\sum_{l,M}\frac{\partial}{\partial P_{\varphi}}C_{M}^{2}\mathcal{E}^{2}\frac{G_{m^{\prime}l}^{*}G_{ml}}{\left|\partial\Omega_{l}/\partial P_{\varphi}\right|_{res}}\mathcal{F}_{lM}\frac{\partial}{\partial P_{\varphi}}f_{lM}+\nu_{eff}^{3}\sum_{l,M}\frac{\partial^{2}}{\partial P_{\varphi}^{2}}\left(f_{lM}-f_{0}\right), \end{eqnarray} and the amplitude evolution satisfies \begin{eqnarray} C_{M}\left(t\right)\sim e^{\left(\gamma_{L}+\gamma_{d}\right)t}\Rightarrow\frac{dC_{M}^{2}}{dt}=2\left(\gamma_{L}+\gamma_{d}\right)C_{M}^{2}. \end{eqnarray} RBQ1D employs a finite-difference scheme used for numerical integration of the distribution function. It recovers several scenarios for amplitude evolution, such as pulsations, intermittency and quasi-steady saturation, see figure (a). The code is interfaced with linear ideal MHD code NOVA, which provides eigenstructures, and the stability code NOVA-K, which provides damping rates and wave-particle interaction matrices for resonances in 3D constant of motion space. RBQ1D employs an iterative procedure to account for mode structure non-uniformities within the resonant island [4]. RBQ1D has been subject to rigorous verification exercises [4]. Both the wave-particle resonant diffusion and Coulomb collisional scattering diffusion operator are thoroughly verified against analytical expressions in limiting cases. In addition, the collisional scattering frequency dependence of the modes saturation level is similar to the value expected by the analytical theory, as shown in figure (b).
[1] H. L. Berk, B. N. Breizman et al., Nucl. Fusion 35, 1661 (1995)
[2] V. N. Duarte, Ph.D. Thesis, U. São Paulo (2017)
[3] N. Gorelenkov, V. Duarte et al., Nucl. Fusion 58, 082016 (2018)
[4] N. N. Gorelenkov, Invited Talk, APS-DPP (2018)
[#h50: V. Duarte, 2018-09-25]

M3D-C$^1$ Extended MHD code

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 17, 102508 (2010)
[2] S.C. Jardin, N. Ferraro et al., Com. Sci. Disc. 5, 014002 (2012)
[3] N.M. Ferraro, T.E. Evans et al., Nucl. Fusion 53, 073042 (2013)
[#h4: N. Ferraro, 2016-05-23]

Stepped Pressure Equilibrium Code (SPEC)

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. 49, 717 (1996)
[2] H. Grad, Phys. Fluids 10, 137 (1967)
[3] S.R. Hudson, R.L. Dewar et al., Phys. Plasmas 19, 112502 (2012)
[4] M.J. Hole, S.R. Hudson & R.L. Dewar, J. Plas. Physics 72, 1167 (2006)
[5] S.R. Hudson, M.J. Hole & R.L. Dewar, Phys. Plasmas 14, 052505 (2007)
[6] J.B. Taylor, Rev. Modern Phys. 58, 741 (1986)
[7] G.R. Dennis, S.R. Hudson et al.,Phys. Plasmas 20, 032509 (2013)
[#h5: S.R. Hudson, 2016-05-23]

Designing Stellarator Coils with the COILOPT++, Coil-Optimization Code

Concept design for a quasi-axisymmetric stellarator fusion reactor with the modular coils straightened and spaced for ease of access. (Figure courtesy of Tom Brown.)

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 1, 253 (1958)
[2] Dennis J. Strickler, Lee A. Berry & Steven P. Hirshman, Fusion Sci. & Technol. 41, 107 (2002)
[3] J.A. Breslau, N. Pomphrey et al., in preparation (2016)
[4] George H. Neilson, David A. Gates et al., IEEE Trans. Plasma Sci. 42, 489 (2014)
[5] Wikipedia, Biot-Savart Law
[6] Wikipedia, Levenberg–Marquardt Algorithm
[7] Wikipedia, Differential Evolution
[#h12: J. Breslau, 2016-05-23]

STELLOPT: Stellarator Optimization and Equilibrium Reconstruction Code

Design of the National Compact Stellarator eXperiment (NCSX) optimized by STELLOPT. The shape of three-dimensional nested flux surfaces is optimized to have good MHD stability, plasma confinement and turbulence transport. External non-planar coils with relatively simple geometries, large coil-plasma space and coil-coil separation are also obtained.

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 1, 253 (1958)
[2] S.P. Hirshman, D.A. Spong et al., Phys. Rev. Lett. 80, 528 (1998)
[3] D.A. Spong, S.P. Hirshman et al., Nucl. Fusion 40, 563 (2000)
[4] S.P. Hirshman, J.C. Whitson, Phys. Fluids 26, 3553 (1983)
[5] Donald W. Marquardt, SIAM J. Appl. Math. 11, 431 (1963)
[6] D.E. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning (Reading, Addison-Wesley, 1989)
[7] Riccardo Poli, James Kennedy & Tim Blackwell, Swarm Intell. 1, 33 (2007)
[8] S. Lazerson & I.T. Chapman, Plasma Phys. Control. Fusion 55, 084004 (2013)
[9] J.C. Schmitt, J. Bialek et al., Rev. Sci. Instrum. 85, 11E817 (2014)
[10] S. Lazerson S and the DIII-D Team, Nucl. Fusion 55, 1 (2015)
[#h35: S. Lazerson, 2016-08-18]