# The DESC Stellarator Code Suite Part III: Quasi-symmetry optimization

D. W. Dudt 1  R. Conlin 1  D. Panici 1  E. Kolemen 1<EMAIL_ADDRESS>1Princeton University, Princeton, New Jersey 08544

###### Abstract

The DESC stellarator optimization code takes advantage of advanced numerical
methods to search the full parameter space much faster than conventional
tools. Only a single equilibrium solution is needed at each optimization step
thanks to automatic differentiation, which efficiently provides exact
derivative information. A Gauss-Newton trust-region optimization method uses
second-order derivative information to take large steps in parameter space and
converges rapidly. With just-in-time compilation and GPU portability, high-
dimensional stellarator optimization runs take orders of magnitude less
computation time with DESC compared to other approaches. This paper presents
the theory of the DESC fixed-boundary local optimization algorithm along with
demonstrations of how to easily implement it in the code. Example quasi-
symmetry optimizations are shown and compared to results from conventional
tools. Three different forms of quasi-symmetry objectives are available in
DESC, and their relative advantages are discussed in detail. In the examples
presented, the triple product formulation yields the best optimization results
in terms of minimized computation time and particle transport. This paper
concludes with an explanation of how the modular code suite can be extended to
accommodate other types of optimization problems.

##  1 Introduction

###  1.1 Motivation

Stellarators are an attractive candidate for fusion energy generation due to
their potential for steady-state operation and reduced susceptibility to major
disruptions. Unlike axisymmetric configurations such as tokamaks, the three-
dimensional (3D) magnetic geometry of stellarators unfortunately does not
guarantee that all charged particles will be confined by these devices.
Trapped particles, which do not sample full magnetic field lines, are subject
to radial drifts that can cause them to cross magnetic flux surfaces and leave
the plasma volume. This necessitates stellarator optimization: the search for
magnetic geometries that maximize the confinement of charged particles. There
are many other physics and engineering objectives that are also important
qualities of an “optimal” stellarator design, including the stability of the
equilibrium magnetic field and the complexity of the coils needed to generate
that field. The design space of stellarators is very large – their 3D external
magnetic fields posses about an order of magnitude more degrees of freedom
than tokamaks (Boozer, 2015) – and this high-dimensional optimization problem
is computationally challenging. Developing efficient methods to search this
large parameter space and find optimal solutions is essential for the
continued improvement of stellarator performance towards practical energy
production.

An ideal magnetohydrodynamic (MHD) static equilibrium magnetic field
$\mathbf{B}$ and current density $\mathbf{J}$ satisfies

$\displaystyle\mathbf{J}\times\mathbf{B}$ $\displaystyle=\nabla p$ (1a)
$\displaystyle\nabla\times\mathbf{B}$ $\displaystyle=\mu_{0}\mathbf{J}$ (1b)
$\displaystyle\nabla\cdot\mathbf{B}$ $\displaystyle=0$ (1c)

throughout the plasma volume, where $\mu_{0}$ is the magnetic constant and $p$
is the pressure. Magnetic field lines are generally chaotic in three-
dimensional geometries, but nested flux surfaces are desirable for achieving
higher plasma pressures. This additional requirement of magnetic surfaces is
given by

$\mathbf{B}\cdot\nabla\psi=0,$ (2)

where $2\pi\psi$ is the toroidal magnetic flux through a surface of constant
toroidal angle $\phi$. Equilibrium solutions that meet the criteria of
Equations 1 and 2 are uniquely defined by the shape of their boundary surface
along with the pressure and rotational transform (or plasma current density)
profiles (Kruskal & Kulsrud, 1958). Stellarator optimization seeks to find
configurations with the most desirable properties in the landscape of these
free variables.

###  1.2 Literature Review

Attempts have been made to avoid searching this high dimensional space by
directly constructing optimized stellarators using expansions in the distance
from the magnetic axis (Landreman & Sengupta, 2018; Landreman et al., 2019;
Plunk et al., 2019). While these methods may provide useful starting points
for further optimization, they are inadequate for designing realistic
stellarators with acceptable confinement near the plasma edge. The traditional
approach to this problem is through fixed-boundary optimization, which was
pioneered by Nührenberg and Zille (Nührenberg & Zille, 1988). The boundary
shape is typically assumed to posses stellarator symmetry and parameterized in
a finite double Fourier series of the form

$\displaystyle R^{b}(\theta,\phi)$
$\displaystyle=\sum_{n=-N}^{N}\sum_{m=0}^{M}R^{b}_{mn}\cos(m\theta-
nN_{FP}\phi)$ (3a) $\displaystyle Z^{b}(\theta,\phi)$
$\displaystyle=\sum_{n=-N}^{N}\sum_{m=0}^{M}Z^{b}_{mn}\sin(m\theta-
nN_{FP}\phi)$ (3b)

where $\theta$ is an arbitrary poloidal angle, $(R,\phi,Z)$ are the toroidal
coordinates, and $N_{FP}$ is the number of (toroidal) field periods of the
device. For a given set of desired profiles, the Fourier coefficients
$\\{R^{b}_{mn},Z^{b}_{mn}\\}$ are then used as the optimization parameters to
search for an equilibrium with target properties. A set of coils that can
produce this desired plasma boundary are then found through a second
optimization stage, which must consider manufacturing constraints. Single-
stage methods that combine the equilibrium and coil optimization problems have
also been explored recently and promise faster computation times (Giuliani et
al., 2020).

Stellarator optimization is therefore cast as a conventional multidimensional
optimization problem, and a plethora of algorithms are available for problems
of this type. Derivative-free routines have been used in the past, such as
Brent’s method as implemented in the ROSE code (Drevlak et al., 2019). These
approaches are inexpensive per iteration, but can suffer from slow convergence
and are best suited for smaller problems. Large scale non-linear optimization
is significantly aided by gradient information, but analytic derivatives are
difficult to obtain for complex objectives. In the Levenberg-Mardquart
algorithm used by the STELLOPT code (Spong et al., 1998), the required
derivatives are computed through finite differences. The cost of this approach
is that it requires numerous equilibrium evaluations, and since this scales
with dimensionality it can put a practical limit on the boundary resolution
$(M,N)$ available to study. Furthermore, finite differences approximations are
very sensitive to the step size, which can make this approach inaccurate and
unreliable. Recent progress has been made through the use of adjoint methods,
as in the ALPOpt code (Paul et al., 2021). The adjoint approach reduces the
computation burden to only two equilibrium solutions and avoids the noise of
numerical derivatives (Antonsen et al., 2019; Paul et al., 2020). Adjoint
methods are labor-intensive to implement, however, and are not applicable for
all desired objectives and derivate orders.

Another method for obtaining gradient information is automatic
differentiation: a programmatic technique that exploits the primitive
operations underlying all mathematical instructions to efficiently compute
exact derivatives of arbitrary (differentiable) functions to any order.
Automatic differentiation (AD) is trivial to implement in modern programming
languages and is already being used in the new stellarator coil design code
FOCUSADD (McGreivy et al., 2021), but has not yet been incorporated into other
optimization objectives. The bottleneck process for stellarator optimization
is the equilibrium calculations that most physics and engineering objectives
depend on, and all of the existing codes – STELLOPT, ROSE, ALPOpt, SIMSOPT
(Landreman et al., 2021), etc. – primarily rely on the equilibrium solver VMEC
(Hirshman & Whitson, 1983). Since VMEC was written in FORTRAN without support
for automatic differentiation, this approach is not available for these
optimization codes and they are limited to using inaccurate or expensive
gradient information.

###  1.3 Optimization with DESC

DESC (Dudt et al., 2022) is a pseudo-spectral equilibrium solver written in
Python with the latest numerical techniques, and has now developed into a full
stellarator optimization code with this equilibrium solver at its core. Unlike
the decentralized approach of other stellarator optimization suites that wrap
their algorithms around calls to external codes that evaluate each objective,
DESC’s philosophy is to compute everything “in-house”. This requires re-
writing existing programs as part of the DESC package, but it allows
derivative information to freely propagate throughout the code and results in
significant performance improvements. The centralized design of DESC offers
many significant advantages over existing stellarator optimization tools.
First and foremost, the use of AD through the JAX package (Bradbury et al.,
2018) provides exact derivative information from a single equilibrium
solution. The computation time is therefore independent of the size of the
optimization space, in contrast to finite difference gradients which require a
number of equilibrium calculations that scales linearly with the number of
optimization parameters. Since AD is exact it can also improve robustness and
convergence over numerical derivatives. In comparison to adjoint methods, AD
can theoretically be used with any deterministic optimization objective, is
much easier to implement, and can compute higher orders of derivatives with no
additional coding. Furthermore, the single equilibrium required at each
optimization step is found using a perturbation of the previous solution as
the initial guess, which converges much faster than a “cold start”. Computing
equilibrium solutions has historically been the bottleneck of the stellarator
optimization process, and this approach minimizes that burden. Altogether,
DESC provides a tool to efficiently search the stellarator optimization
landscape and locate configurations of interest.

This is the final paper in a three-part series about the DESC stellarator code
suite. Part I (Panici et al., 2022) highlights the advantages of the DESC
equilibrium solver through a comparison to VMEC. Part II (Conlin et al., 2022)
presents the novel perturbation and continuation methods available in DESC.
This paper, Part III, explains the theory and practice of how these tools can
be combined to solve stellarator optimization problems with DESC. Section 2
details how the optimization problem is formulated and solved in conjunction
with the equilibrium problem. Different metrics of quasi-symmetry are then
defined as examples of optimization objectives, and code snippets are given to
demonstrate how this optimization can be performed with DESC. In Section 3,
the results of a quasi-symmetry optimization problem are presented for a low-
dimensional example that is easy to visualize. The effects of using different
quasi-symmetry objectives and numerical options are compared and contrasted.
These results are intended to demonstrate the tools that DESC provides, not to
show a novel quasi-symmetric result. Speed comparisons are also made between
DESC and STELLOPT to reveal how the computation times of both codes scale for
larger problems. This paper only covers fixed-boundary optimization and not
coil optimization, which is a natural extension and still under development.
Global and stochastic optimization methods, such as particle swarm and genetic
algorithms both used by STELLOPT and ROSE, are not discussed but could be
included in a future version of the code.

##  2 Methods

###  2.1 Equilibrium Solver

As explained in a previous publication (Dudt & Kolemen, 2020), DESC solves the
inverse equilibrium problem by minimizing the MHD force balance errors at a
series of collocation points. The inverse coordinate mapping is given by
$R=R(\rho,\theta,\zeta)$, $Z=Z(\rho,\theta,\zeta)$, and
$\lambda=\lambda(\rho,\theta,\zeta)$, where $(R,\phi,Z)$ are the toroidal
coordinates, $(\rho,\vartheta,\zeta)$ are straight field-line coordinates
corresponding to the choice of toroidal angle $\zeta=\phi$, and the poloidal
angle is related to the coordinate used in the boundary parameterization (3)
through: $\vartheta=\theta+\lambda(\rho,\theta,\zeta)$. Nested flux surfaces
are assumed, and the flux surface label is chosen to be
$\rho=\sqrt{\psi/\psi_{\text{edge}}}$. This map between the flux and toroidal
coordinate systems is discretized with global Fourier-Zernike basis sets, and
the spectral coefficients are denoted as $R_{lmn}$, $Z_{lmn}$, and
$\lambda_{lmn}$ ($l$, $m$, and $n$ are the radial, poloidal, and toroidal mode
numbers, respectively, of the basis functions).

Let the MHD equilibrium force balance error be defined as
$\mathbf{F}=\mathbf{J}\times\mathbf{B}-\nabla p$, which has a radial and
helical component:

$\displaystyle f_{\rho}$
$\displaystyle=\left(\sqrt{g}\left(B^{\zeta}J^{\theta}-B^{\theta}J^{\zeta}\right)-p^{\prime}\right)|\nabla\rho|$
(4a) $\displaystyle f_{\beta}$
$\displaystyle=\sqrt{g}J^{\rho}|\mathbf{\beta}|$ (4b)

where $\mathbf{\beta}=B^{\zeta}\nabla\theta-B^{\theta}\nabla\zeta$ and
$\sqrt{g}$ is the Jacobian of the $(\rho,\theta,\zeta)$ computational
coordinate system. Evaluating each function at a series of coordinates results
in a system of nonlinear equations of the form

$\mathbf{f}(\mathbf{x},\mathbf{c})=\begin{bmatrix}f_{\rho,i}\\\
f_{\beta,j}\end{bmatrix}$ (5)

where $i=0,\ldots,I$ and $j=0,\ldots,J$ denote the $i$th and $j$th collocation
points for each function. Here $\mathbf{x}$ represents the independent
variables, which are a subset of the spectral coefficients $R_{lmn}$,
$Z_{lmn}$, and $\lambda_{lmn}$ such that the boundary conditions
$R(\rho=1)=R^{b}(\theta,\phi)$ and $Z(\rho=1)=Z^{b}(\theta,\phi)$ are
satisfied. The vector $\mathbf{c}$ is all of the parameters that define a
unique equilibrium problem, which includes the boundary coefficients
$R^{b}_{mn}$ and $Z^{b}_{mn}$, coefficients for the profiles $p(\rho)$ and
$\iota(\rho)$, and the total flux through the boundary surface $\psi_{edge}$.
An equilibrium is defined as the least-squares solution

$\mathbf{x}^{*}=\text{arg
min}_{\mathbf{x}}|\mathbf{f}(\mathbf{x},\mathbf{c})|^{2}$ (6)

for a fixed set of parameters $\mathbf{c}$, and is typically solved with a
quasi-Newton method.

###  2.2 Optimization Approach

Stellarator optimization desires to find equilibrium solutions that also
satisfy secondary objectives. Let the system of equations
$\mathbf{g}(\mathbf{x},\mathbf{c})$ represent a set of costs that are desired
to be minimized. Examples of quasi-symmetry objective functions are given in
Section 2.3, but these could be any general engineering or physics target. The
optimal configuration is defined as the least-squares solution

$\mathbf{c}^{*}=\text{arg
min}_{\mathbf{c}}|\mathbf{g}(\mathbf{x}^{*},\mathbf{c})|^{2}$ (7)

where $\mathbf{x}^{*}$ is the equilibrium solution given by (6) for the set of
parameters $\mathbf{c}^{*}$.

DESC employs a constrained optimization approach similar to the prediction-
correction continuation method for solving systems of equations (Yu & Dutton,
1998). Starting from an initial equilibrium solution
$(\mathbf{x}_{k}^{*},\mathbf{c}_{k})$ that it not optimal, the parameters are
perturbed to a new state
$(\mathbf{x}_{k+1},\mathbf{c}_{k+1})=(\mathbf{x}_{k}^{*}+\Delta\mathbf{x},\mathbf{c}_{k}+\Delta\mathbf{c})$
that better satisfies the objective:
$|\mathbf{g}(\mathbf{x}_{k+1},\mathbf{c}_{k+1})|^{2}<|\mathbf{g}(\mathbf{x}_{k},\mathbf{c}_{k})|^{2}$.
The perturbation is restricted to a subspace that maintains the equilibrium
constraint, and is performed using a trust region Newton method as detailed in
the remainder of this section. Equilibrium force balance is never satisfied
exactly, however, and this “prediction” step may violate the equilibrium
constraints slightly. Therefore the equilibrium problem is then re-solved with
the new parameters $\mathbf{c}_{k+1}$ to get the corresponding equilibrium
state $\mathbf{x}_{k+1}^{*}$, but this “correction” step should be a minor
adjustment from the perturbed state $\mathbf{x}_{k+1}$. This process is then
iterated until a set of parameters $\mathbf{c}^{*}$ are found that satisfy the
objectives to the desired tolerance, and the Newton methods are expected to
yield quadratic convergence. These iterations can also be included in a
broader optimization loop that changes the resolution and dimensionality of
the problem at each step.

The optimization algorithm implemented in DESC can be summarized by the
pseudocode in Listing 1. Stopping criteria are determined by relative
tolerances on the cost function values and step sizes, in addition to maximum
numbers of iterations.

⬇

1while optimization stopping criteria are not met:

2

3 perturb equilibrium solution to improve objective

4

5 while equilibrium stopping criteria are not met:

6

7 solve equilibrium force balance

Listing 1: Pseudocode outlining the DESC optimization algorithm.

The derivation of the optimal perturbations begins by Taylor expanding both
the constraint (equilibrium) and objective functions about the current state:

$\displaystyle\mathbf{f}(\mathbf{x}+\Delta\mathbf{x},\mathbf{c}$
$\displaystyle+\Delta\mathbf{c})=\mathbf{f}(\mathbf{x},\mathbf{c})+\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\Delta\mathbf{x}+\frac{\partial\mathbf{f}}{\partial\mathbf{c}}\Delta\mathbf{c}$
(8a)
$\displaystyle+\frac{1}{2}\frac{\partial^{2}\mathbf{f}}{\partial\mathbf{x}^{2}}\Delta\mathbf{x}\Delta\mathbf{x}^{T}+\frac{1}{2}\frac{\partial^{2}\mathbf{f}}{\partial\mathbf{c}^{2}}\Delta\mathbf{c}\Delta\mathbf{c}^{T}+\frac{\partial^{2}\mathbf{f}}{\partial\mathbf{x}\partial\mathbf{c}}\Delta\mathbf{x}\Delta\mathbf{c}^{T}$
$\displaystyle\mathbf{g}(\mathbf{x}+\Delta\mathbf{x},\mathbf{c}$
$\displaystyle+\Delta\mathbf{c})=\mathbf{g}(\mathbf{x},\mathbf{c})+\frac{\partial\mathbf{g}}{\partial\mathbf{x}}\Delta\mathbf{x}+\frac{\partial\mathbf{g}}{\partial\mathbf{c}}\Delta\mathbf{c}$
(8b)
$\displaystyle+\frac{1}{2}\frac{\partial^{2}\mathbf{g}}{\partial\mathbf{x}^{2}}\Delta\mathbf{x}\Delta\mathbf{x}^{T}+\frac{1}{2}\frac{\partial^{2}\mathbf{g}}{\partial\mathbf{c}^{2}}\Delta\mathbf{c}\Delta\mathbf{c}^{T}+\frac{\partial^{2}\mathbf{g}}{\partial\mathbf{x}\partial\mathbf{c}}\Delta\mathbf{x}\Delta\mathbf{c}^{T}$

The perturbations themselves are assumed to be a combination of first and
second-order terms

$\displaystyle\Delta\mathbf{x}$
$\displaystyle=\epsilon\mathbf{x}_{1}+\epsilon^{2}\mathbf{x}_{2}$ (9a)
$\displaystyle\Delta\mathbf{c}$
$\displaystyle=\epsilon\mathbf{c}_{1}+\epsilon^{2}\mathbf{c}_{2}$ (9b)

where $\epsilon\ll 1$ is some small parameter.

This paper refers to the “order” of a method by the number of terms included
in the perturbation expansion (9). However, the terminology typically used in
the optimization literature would be a degree higher than this. The
discrepancy is because the Jacobian matrices of first derivatives such as
$\frac{\partial\mathbf{g}}{\partial\mathbf{c}}$ give the Hessian matrices of
second derivatives for the scalar least-squares problem in (7).

####  2.2.1 First Order

Substituting (9) into (8) and collecting only the terms up to order $\epsilon$
yields:

$\displaystyle\mathbf{0}$
$\displaystyle=\mathbf{f}(\mathbf{x},\mathbf{c})+\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\epsilon\mathbf{x}_{1}+\frac{\partial\mathbf{f}}{\partial\mathbf{c}}\epsilon\mathbf{c}_{1}$
(10a) $\displaystyle\mathbf{0}$
$\displaystyle=\mathbf{g}(\mathbf{x},\mathbf{c})+\frac{\partial\mathbf{g}}{\partial\mathbf{x}}\epsilon\mathbf{x}_{1}+\frac{\partial\mathbf{g}}{\partial\mathbf{c}}\epsilon\mathbf{c}_{1}.$
(10b)

Equation (10a) can be solved for $\epsilon\mathbf{x}_{1}$ in terms of
$\epsilon\mathbf{c}_{1}$

$\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\epsilon\mathbf{x}_{1}=-\mathbf{f}(\mathbf{x},\mathbf{c})-\frac{\partial\mathbf{f}}{\partial\mathbf{c}}\epsilon\mathbf{c}_{1},$
(11)

and substituting this expression into (10b) gives an equation for the first-
order parameter perturbation:

$\left[\frac{\partial\mathbf{g}}{\partial\mathbf{x}}\left(\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\right)^{-1}\frac{\partial\mathbf{f}}{\partial\mathbf{c}}-\frac{\partial\mathbf{g}}{\partial\mathbf{c}}\right]\epsilon\mathbf{c}_{1}=\mathbf{g}(\mathbf{x},\mathbf{c})-\frac{\partial\mathbf{g}}{\partial\mathbf{x}}\left(\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\right)^{-1}\mathbf{f}(\mathbf{x},\mathbf{c}).$
(12)

Equation (12) is a linear system that can be solved for the optimal
perturbation $\epsilon\mathbf{c}_{1}$ using a trust-region step, as explained
in Section 2.2.3. The corresponding change to the independent variables,
$\epsilon\mathbf{x}_{1}$, is then determined from (11). All of the Jacobian
matrices are computed through automatic differentiation, and
$\frac{\partial\mathbf{f}}{\partial\mathbf{x}}$ is typically already know from
solving the equilibrium problem (6) at the previous optimization step. A
pseudo-inverse is taken for the term
$\left(\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\right)^{-1}$.

####  2.2.2 Second Order

Collecting the terms proportional to $\epsilon^{2}$ that were omitted from (8)
yields:

$\displaystyle\mathbf{0}=\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\epsilon^{2}\mathbf{x}_{2}$
$\displaystyle+\frac{\partial\mathbf{f}}{\partial\mathbf{c}}\epsilon^{2}\mathbf{c}_{2}+\frac{1}{2}\frac{\partial^{2}\mathbf{f}}{\partial\mathbf{x}^{2}}\epsilon^{2}\mathbf{x}_{1}\mathbf{x}_{1}^{T}$
(13a)
$\displaystyle+\frac{1}{2}\frac{\partial^{2}\mathbf{f}}{\partial\mathbf{c}^{2}}\epsilon^{2}\mathbf{c}_{1}\mathbf{c}_{1}^{T}+\frac{\partial^{2}\mathbf{f}}{\partial\mathbf{x}\partial\mathbf{c}}\epsilon^{2}\mathbf{x}_{1}\mathbf{c}_{1}^{T}$
$\displaystyle\mathbf{0}=\frac{\partial\mathbf{g}}{\partial\mathbf{x}}\epsilon^{2}\mathbf{x}_{2}$
$\displaystyle+\frac{\partial\mathbf{g}}{\partial\mathbf{c}}\epsilon^{2}\mathbf{c}_{2}+\frac{1}{2}\frac{\partial^{2}\mathbf{g}}{\partial\mathbf{x}^{2}}\epsilon^{2}\mathbf{x}_{1}\mathbf{x}_{1}^{T}$
(13b)
$\displaystyle+\frac{1}{2}\frac{\partial^{2}\mathbf{g}}{\partial\mathbf{c}^{2}}\epsilon^{2}\mathbf{c}_{1}\mathbf{c}_{1}^{T}+\frac{\partial^{2}\mathbf{g}}{\partial\mathbf{x}\partial\mathbf{c}}\epsilon^{2}\mathbf{x}_{1}\mathbf{c}_{1}^{T}$

In a similar process, equation (13a) can be solved for
$\epsilon\mathbf{x}_{2}$ in terms of $\epsilon\mathbf{c}_{1}$ and
$\epsilon\mathbf{c}_{2}$

$\displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\epsilon^{2}\mathbf{x}_{2}$
$\displaystyle=-\frac{\partial\mathbf{f}}{\partial\mathbf{c}}\epsilon^{2}\mathbf{c}_{2}-\frac{1}{2}\frac{\partial^{2}\mathbf{f}}{\partial\mathbf{x}^{2}}\epsilon^{2}\mathbf{x}_{1}\mathbf{x}_{1}^{T}$
(14)
$\displaystyle-\frac{1}{2}\frac{\partial^{2}\mathbf{f}}{\partial\mathbf{c}^{2}}\epsilon^{2}\mathbf{c}_{1}\mathbf{c}_{1}^{T}-\frac{\partial^{2}\mathbf{f}}{\partial\mathbf{x}\partial\mathbf{c}}\epsilon^{2}\mathbf{x}_{1}\mathbf{c}_{1}^{T},$

and substituting this expression into (13b) gives an equation for the second-
order parameter perturbation:

$\displaystyle\left[\left(\frac{\partial\mathbf{g}}{\partial\mathbf{x}}\right)\right.$
$\displaystyle\left.\left(\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\right)^{-1}\frac{\partial\mathbf{f}}{\partial\mathbf{c}}-\frac{\partial\mathbf{g}}{\partial\mathbf{c}}\right]\epsilon^{2}\mathbf{c}_{2}=$
(15) $\displaystyle\frac{1}{2}$
$\displaystyle\left[\frac{\partial^{2}\mathbf{g}}{\partial\mathbf{x}^{2}}-\left(\frac{\partial\mathbf{g}}{\partial\mathbf{x}}\right)\left(\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\right)^{-1}\frac{\partial^{2}\mathbf{f}}{\partial\mathbf{x}^{2}}\right]\epsilon^{2}\mathbf{x}_{1}\mathbf{x}_{1}^{T}$
$\displaystyle+\frac{1}{2}$
$\displaystyle\left[\frac{\partial^{2}\mathbf{g}}{\partial\mathbf{c}^{2}}-\left(\frac{\partial\mathbf{g}}{\partial\mathbf{x}}\right)\left(\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\right)^{-1}\frac{\partial^{2}\mathbf{f}}{\partial\mathbf{c}^{2}}\right]\epsilon^{2}\mathbf{c}_{1}\mathbf{c}_{1}^{T}$
$\displaystyle+$
$\displaystyle\left[\frac{\partial^{2}\mathbf{g}}{\partial\mathbf{x}\partial\mathbf{c}}-\left(\frac{\partial\mathbf{g}}{\partial\mathbf{x}}\right)\left(\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\right)^{-1}\frac{\partial^{2}\mathbf{f}}{\partial\mathbf{x}\partial\mathbf{c}}\right]\epsilon^{2}\mathbf{x}_{1}\mathbf{c}_{1}^{T}.$

Trust-region steps are again used to solve for the second order correction to
the optimal perturbation, $\epsilon\mathbf{c}_{2}$, and then the independent
variables $\epsilon\mathbf{x}_{2}$, using equations (15) and (14),
respectively. One difference is that the trust-region radius used for the
second order terms is restricted to be a fraction of the magnitude of the
first order step. Note that the same matrix appears on the left-hand side of
Equations (12) and (15), so once the single value decomposition is known from
the first order solution it does not have to be recomputed for the second
order solution. Also note that the Hessian matrices on the right-hand side of
(15) are not needed in full, but can be efficiently computed with Jacobian-
vector products. Therefore this second order approximation requires relatively
little additional computation time after the first order solution is found,
and the extra accuracy in representing the parameter space is usually worth
the trade-off. This adaptive procedure helps to ensure the optimization is
robust to different initial conditions.

####  2.2.3 Trust-Region Step

Equations (11), (12), (14), and (15) are all linear systems of the form
$\mathbf{A}\mathbf{y}=\mathbf{b}$. A trust-region method is used instead of
taking the full Newton step, which improves convergence while far from the
optimum (Nocedal & Wright, 2006). This involves casting these equations into a
problem of the form

$\text{min}_{\mathbf{y}}|\mathbf{A}\mathbf{y}-\mathbf{b}|^{2},\hskip
2.84526pt|\mathbf{y}|\leq\delta$ (16)

where $\delta$ is the trust region radius. The solution requires a single-
value decomposition of the matrix $\mathbf{A}$, but this computation time is
typically insignificant compared to the automatic differentiation used to
evaluate the matrix. The trust region radius is chosen through an adaptive
procedure based on the ratio of the cost reduction predicted by the model to
the actual reduction achieved by the new equilibrium solution. If this ratio
is above $0.75$ and the trust region radius was restricting the step size, the
radius is doubled for the next iteration. The trust region radius is reduced
to a quarter of its former size if the ratio is below $0.25$ – this includes
prospective steps that fail to actually reduce the cost function, which are
rejected.

###  2.3 Quasi-symmetry Objective Functions

Charged particle confinement is not guaranteed in three-dimensional toroidal
geometries, and this is typically the primary objective of stellarator
optimization. Since neoclassical confinement calculations are historically too
expensive to run within an optimization loop, quasi-symmetry is a popular
proxy function for particle confinement that is much cheaper to compute.
Quasi-symmetry ensures particles remain close to flux surfaces by
approximately conserving some canonical momentum (Rodriguez et al., 2020). The
formal definition of a quasi-symmetric magnetic field is the existence of a
vector field $\mathbf{u}$ such that:

$\displaystyle\mathbf{B}\times\mathbf{u}$ $\displaystyle=\nabla\psi$ (17a)
$\displaystyle\mathbf{u}\cdot\nabla B$ $\displaystyle=0$ (17b)
$\displaystyle\nabla\cdot\mathbf{u}$ $\displaystyle=0.$ (17c)

However, this definition is not amenable to computational analysis. Three
other quasi-symmetry objective functions are implemented in DESC, and although
they are theoretically equivalent, each one can have distinct numerical
consequences in practice.

####  2.3.1 Boozer Coordinates

A magnetic field is quasi-symmetric if its magnitude can be written in the
form

$|\mathbf{B}|=B(\psi,M\vartheta_{B}-N\zeta_{B})$ (18)

where $\vartheta_{B}$ and $\zeta_{B}$ are the poloidal and toroidal angles in
Boozer coordinates, respectively, and $M$, $N$ $\in\mathbb{Z}$ determine the
type of quasi-symmetry. The optimization objective function derived from this
definition is

$\mathbf{f}_{B}=\\{B_{mn}~{}|~{}m/n\neq M/N\\}$ (19)

where $B_{mn}$ are the Fourier coefficients of $|\mathbf{B}|$ in Boozer
coordinates on a particular surface, and the $m=n=0$ mode is always excluded.
This corresponds to the set of non-symmetric Boozer modes, which must vanish
in quasi-symmetry. An associated scalar, dimensionless metric is

$\hat{f}_{B}=\frac{|\mathbf{f}_{B}|}{\sqrt{\sum_{m,n}B_{mn}^{2}}}$ (20)

which is commonly used to quantify the departure from true quasi-symmetry.
This is the traditional definition of quasi-symmetry, and the normalized
scalar has the convenient property of always being in the range
$\hat{f}_{B}\in[0,1)$. However, the transformation to Boozer coordinates can
be computationally expensive. DESC uses the same Boozer coordinate
transformation algorithm as BOOZ_XFORM, but was rewritten within the python
code to work with automatic differentiation. This implementation does not take
advantage of “generalized” Boozer coordinates and therefore assumes that an
equilibrium exists (Rodriguez & Bhattacharjee, 2021), but the force balance is
generally never satisfied to machine precision for numerical solutions.
Furthermore, this definition only quantifies the “global” deviation from
quasi-symmetry on each flux surface rather than providing a “local” deviation
throughout the surface.

####  2.3.2 Two-Term

A magnetic field is also quasi-symmetric if the quantity

$C=\frac{\left(\mathbf{B}\times\nabla\psi\right)\cdot\nabla
B}{\mathbf{B}\cdot\nabla B}$ (21)

is a flux function: $C=C(\psi)$. This definition assumes
$\mathbf{J}\cdot\nabla\rho=J^{\rho}=0$, which is less restrictive than the
requirement of force balance used in the true Boozer form, but must hold in
equilibrium by (4b) with the same caveat about numerical resolution. In Boozer
coordinates, this flux function is also equivalent to

$C=\frac{MG+NI}{M\mbox{$\,\iota\\!\\!$-}-N}.$ (22)

The covariant components of the magnetic field in Boozer coordinates can be
computed in any flux coordinate system by

$\displaystyle G$ $\displaystyle=\langle B_{\zeta}\rangle$ (23a)
$\displaystyle I$ $\displaystyle=\langle B_{\theta}\rangle$ (23b)

where the flux surface average of a quantity $Q$ is defined as

$\langle Q\rangle\equiv\frac{\int_{0}^{2\pi}\int_{0}^{2\pi}Q\sqrt{g}d\theta
d\zeta}{\int_{0}^{2\pi}\int_{0}^{2\pi}\sqrt{g}d\theta d\zeta}.$ (24)

A more useful form that avoids computational singularities is

$f_{C}=\left(M\mbox{$\,\iota\\!\\!$-}-N\right)\left(\mathbf{B}\times\nabla\psi\right)\cdot\nabla
B-\left(MG+NI\right)\mathbf{B}\cdot\nabla B.$ (25)

Evaluating this function at a series of collocation points on a given surface
yields the vector form

$\mathbf{f}_{C}=\\{f_{C}(\theta_{i},\zeta_{j})~{}|~{}i\in[0,2\pi),j\in[0,2\pi/N_{FP})\\}.$
(26)

These residuals have units of $\mathrm{T^{3}}$; a dimensionless scalar metric
is chosen to be

$\hat{f}_{C}=\frac{\langle|f_{C}|\rangle}{\langle B\rangle^{3}}.$ (27)

This quasi-symmetry objective does not rely on a transformation to Boozer
coordinates, but it still requires specifying the helicity of the quasi-
symmetry. Unlike the Boozer form, this two-term definition can reveal the
local quasi-symmetry errors within each surface.

####  2.3.3 Triple Product

Finally, a magnetic field is also quasi-symmetric if the function

$f_{T}=\nabla\psi\times\nabla B\cdot\nabla\left(\mathbf{B}\cdot\nabla
B\right)$ (28)

vanish throughout the region of interest. Similarly, this function is
evaluated at a series of collocation points to yield the vector form

$\mathbf{f}_{T}=\\{f_{T}(\theta_{i},\zeta_{j})~{}|~{}i\in[0,2\pi),j\in[0,2\pi/N_{FP})\\}.$
(29)

This function has units of $\frac{\text{T}^{4}}{\text{m}^{2}}$; a
dimensionless scalar metric is chosen to be

$\hat{f}_{T}(\psi)=\frac{\langle R\rangle^{2}\langle|f_{T}|\rangle}{\langle
B\rangle^{4}}.$ (30)

The triple product form benefits from the same advantages as the two-term
definition: it does not rely on a transformation to Boozer coordinates, and
provides local errors rather than a flux-surface quantity. Furthermore, it
does not assume any equilibrium conditions or require the helicity to be
specified a-priori, which could be useful for certain optimization scenarios.
This form is also more amenable to optimizing for quasi-symmetry throughout a
volume rather than a single flux surface, which may not be physically
achievable but could still have practical applications (Garren & Boozer,
1991a, b).

##  3 Results

###  3.1 Code Implementation

DESC is designed to have an approachable user-interface through python
scripts, similar to the SIMSOPT framework. The code snippet in Listing 2
demonstrates how to run DESC as it was used to generate the results in Section
3.2. This example shows the QuasisymmetryTwoTerm function corresponding to
$\mathbf{f}_{C}$, targeting quasi-helical symmetry. It is also normalized
using the same denominator as (27). The syntax for the objective functions
QuasisymmetryBoozer and QuasisymmetryTripleProduct corresponding to
$\mathbf{f}_{B}$ and $\mathbf{f}_{T}$, respectively, are similar.

The "opt_subspace" argument is a transformation matrix that can be used to
restrict the optimization parameters to a custom subspace. In this example, it
is used to optimize over two boundary coefficients given in the double-angle
Fourier series basis of (3). DESC uses a different but equivalent basis from
this VMEC convention, and the vmec_boundary_subspace function relates them
through the Ptolemy identities. The results presented in this paper were run
with the default options in DESC, but the numerical resolution, target
surfaces, stopping criteria, and all other options can be fully customized.

Although it is not shown in this example, multiple optimization objectives can
be targeted together with relative weighting between them. Examples of other
objectives that are available in DESC include geometric quantities such as the
plasma volume and aspect ratio. The code is designed in a modular structure so
that new objective functions can easily be added to the existing framework.

⬇

1from desc import set_device

2set_device("gpu") # if running on a GPU

3

4from desc.equilibrium import Equilibrium

5from desc.objectives import (

6 ObjectiveFunction,

7 QuasisymmetryTwoTerm,

8 get_fixed_boundary_constraints,

9)

10from desc.vmec_utils import vmec_boundary_subspace

11

12# load an initial equilibrium solution

13eq = Equilibrium.load("path/to/equilibrium.h5")

14

15QH = (1, eq.NFP) # type of quasi-symmetry

16

17# objective function to optimize

18objective = ObjectiveFunction(

19 QuasisymmetryTwoTerm(helicity=QH, norm=True),

20 get_fixed_boundary_constraints(),

21)

22

23# optimization variables

24perturb_options = {

25 "dRb": True, # optimize R boundary mode RBC(2,1)

26 "dZb": True, # optimize Z boundary mode ZBS(2,1)

27 "opt_subspace": vmec_boundary_subspace(

28 eq, RBC=[2, 1], ZBS=[2, 1]

29 ),

30}

31

32# run optimization

33eq.optimize(objective, perturb_options=perturb_options)

Listing 2: Example code to run quasi-helical symmetry optimization using the
two-term objective function.

###  3.2 Quasi-Symmetry Optimization

This section demonstrates the capability of DESC to optimize stellarator
equilibria, using quasi-symmetry as an example objective function. As the
initial state, the $m=1$, $n=2$ boundary modes for both $R^{b}$ and $Z^{b}$ in
(3) of a quasi-helically symmetric STELLOPT benchmark solution (Landreman,
2022) were modified to degrade the quality of quasi-symmetry. The equilibrium
was then optimized for quasi-symmetry on the last closed flux surface
($\rho=1$) in this two-dimensional parameter space with all other boundary
modes and profiles held fixed at the benchmark solution values. The
optimization was performed using each of the three objective functions
described previously, and with both first and second-order optimization
methods. This simple design space was chosen for comparison to a previous
study of this problem (Rodriguez et al., 2022), and it lends itself well to
visualization. Note that in the example used here (in contrast to the previous
study), quasi-symmetry is only being targeted on a single flux surface
(instead of multiple surfaces) and the rotational transform profile (instead
of the current profile) is held fixed during optimization.

Figure 1: Quasi-symmetry optimization paths using the Boozer coordinate
objective in DESC. The STELLOPT optimization path is shown for comparison.
Figure 2: Quasi-symmetry optimization paths using the two-term objective in
DESC. The optimal STELLOPT solution is shown for comparison. Figure 3: Quasi-
symmetry optimization paths using the triple product objective in DESC. The
optimal STELLOPT solution is shown for comparison.

Figures 1, 2, and 3 show the optimization paths of the DESC first and second-
order methods for each of the three quasi-symmetry objectives. The STELLOPT
optimization path for its Boozer form objective is included in Figure 1, and
the final solution is also indicated in Figures 2 and 3 for reference. The
optimization landscape shown in gray scale in these plots are the respective
dimensionless scalar error metrics given in Section 2.3, which are
proportional but not equivalent to the least-squares cost functions that are
actually being minimized. Comparing across all three figures reveals that each
of the objective functions has a similar narrow valley of quasi-symmetric
solutions, although the global minima are not identical and are still far from
perfect quasi-symmetry. The DESC optimization runs targeting the Boozer
coordinate and two-term objectives converge to results very close to the
STELLOPT solution, while the triple product minimum is farther along the
valley.

Although these two optimization codes converge to similar results, their paths
to arrive there from the initial configuration are very different. The
Levenberg-Marquadt algorithm in STELLOPT rarely permits the full Gauss-Newton
step in favor of a more conservative approach that is often closer to the
steepest descent direction. DESC’s algorithm is conceptually similar with an
adaptive trust-region radius to improve robustness, but allows the full Gauss-
Newton step whenever possible. As a consequence, DESC converges to the optimum
in a few large steps while STELLOPT requires nearly a hundred smaller steps to
traverse the optimization landscape in this example (see Figure 1). In a
process analogous to geodesic acceleration in the Levenberg-Marquadt
algorithm, the optimization step size in DESC is further increased by
including higher-order approximations. Contrasting the first and second-order
results reveals that the additional derivative information enables the
optimizer to take slightly larger steps, which can cause faster convergence in
fewer iterations.

Figure 4: Contours of magnetic field magnitude in Boozer coordinates for the
initial equilibrium (top) and optimized solution (bottom) targeting the Boozer
measure of quasi-symmetry. Figure 5: Comparison of quasi-symmetry and
neoclassical confinement errors for the initial and optimized configurations.
Only the second-order DESC results are shown. The three quasi-symmetry metrics
were computed in DESC; the effective ripple was computed with NEO from VMEC
equilibria; the particle and heat fluxes were computed in SFINCS also using
VMEC equilibria.

Figure 4 gives a qualitative picture of the quasi-helical symmetry
optimization though a comparison of $|B|$ in Boozer coordinates between the
initial and final configurations. The contours of the optimized solution
clearly have the desired slope of $4/1$ that the original stellarator was
lacking, although the maximum field strength was also changed in the process.
However, all of these quasi-symmetry metrics are only used as proxy functions
for the true goal of particle confinement. In order to validate that the
quasi-symmetry optimization does improve the confinement properties, the NEO
(Nemov et al., 1999) and SFINCS codes (Landreman et al., 2014) were run on the
solutions to calculate their neoclassical confinement, and the results for
each equilibrium are included in Figure 5. The quasi-symmetry errors were all
computed in DESC, including for the STELLOPT solution which was re-computed
with the DESC equilibrium solver. The optimal boundary coefficients from each
solution were also used to compute high-resolution VMEC equilibria for inputs
to NEO and SFINCS, including for the DESC solutions. Since the first and
second-order results all converged to within $1.3~{}\text{mm}$ of each other,
the errors were only evaluated for the second-order DESC results.

A surprising result is that all four of the optimized solutions have nearly
equivalent values for the three different quasi-symmetry errors, despite each
solution targeting a different form during optimization. This suggests that
the optimization landscape is flat near the global minimum: all of the
equilibria in the phase space “valley” have similar levels of quasi-symmetry.
The confinement results tell a different story, however. $\epsilon_{eff}$ is
the effective ripple in the $1/\nu$ regime and is commonly used as a measure
of particle transport (Nemov et al., 1999). $\hat{f}_{P}$ and $\hat{f}_{Q}$
are the normalized particle and heat fluxes of the ions:

$\displaystyle\hat{f}_{P}=\frac{\bar{R}}{\bar{n}\bar{v}}\left\langle\int
d^{3}vf_{i}\mathbf{v}_{mi}\cdot\nabla\rho\right\rangle$ (31a)
$\displaystyle\hat{f}_{Q}=\frac{\bar{R}}{\bar{n}\bar{v}^{3}\bar{m}}\left\langle\int
d^{3}vf_{i}\frac{m_{i}v^{2}}{2}\mathbf{v}_{mi}\cdot\nabla\rho\right\rangle$
(31b)

where $\bar{R}=1~{}\text{m}$, $\bar{n}=1\times 10^{20}~{}\text{m}^{-3}$,
$\bar{v}=\sqrt{2\bar{T}/\bar{m}}$, $\bar{T}=1~{}\text{keV}$, and $\bar{m}$ is
the proton mass. All of the optimized solutions substantially reduced the
effective ripple and both fluxes, confirming that quasi-symmetry is correlated
with neoclassical confinement; but the relationship is complex. The DESC
optimization targeting the triple product metric, which converged to a result
farther away in parameter space from the other solutions, achieved nearly an
order of magnitude better confinement than the others according to these
measures. This lack of correlation supports previous claims that quasi-
symmetry is not an accurate indicator of transport levels (Martin & Landreman,
2020; Rodriguez et al., 2022). The triple product quasi-symmetry measure is
the only one that does not depend on equilibrium conditions, and this
independence from the force balance constraints could aid the optimization
process. Unlike the other two metrics which depend on flux surface averaged
quantities, the truly local nature of the triple product definition could also
explain why it was the most effective proxy for confinement in this
optimization example.

###  3.3 Computation Speed

Computation times for running realistic quasi-symmetry optimization problems
in DESC were benchmarked against STELLOPT. STELLOPT is a parallelized FORTRAN
code that runs on multiple CPUs, while DESC is a python code aided by just-in-
time compilation that is currently compatible with a single CPU or GPU. All
tests were run on the Traverse computing cluster at Princeton University using
AMD EPYC 7281 model CPUs, which have 16 cores, and NVIDIA V100 Tensor Core
GPUs, wich have 32G of memory. The same numerical resolutions were used
between the two codes for a fair comparison.

Table 1: STELLOPT computation times (in seconds). The VMEC equilibrium inputs used were: MPOL = 9, NTOR = 8, NS_ARRAY = 17 33 65, and FTOL_ARRAY = 1E-8 1E-10 1E-12. The BOOZ_XFORM optimization inputs used were: MBOZ = 17, NBOZ = 16, and TARGET_HELICITY(65) = 0. | Number of Optimization Parameters
---|---
CPUs | 8 | 24 | 48 | 80
2 | 362 | 939 | 1757 | 2831
4 | 188 | 399 | 665 | 1074
8 | 150 | 246 | 346 | 550
16 | 118 | 180 | 234 | 322
Table 2: DESC computation times (in seconds) for 560 optimization parameters. The numerical resolution used was $L=M=N=8$ for the equilibria and $M=N=16$ for the Boozer spectrum with the $f_{B}$ objective. Quasi-symmetry was targeted on the $\rho=1$ surface. |  | Initial Iteration | Subsequent Iterations
---|---|---|---
Objective |  | Perturb | Solve | Total | Perturb | Solve | Total
$f_{C}$ | 1 GPU | 157 | 33 | 190 | 2 | 3 | 5
$f_{T}$ | 1 GPU | 89 | 33 | 122 | 2 | 3 | 5
$f_{B}$ | 1 CPU | 763 | 174 | 937 | 409 | 24 | 434
$f_{C}$ | 1 CPU | 165 | 168 | 333 | 12 | 25 | 37
$f_{T}$ | 1 CPU | 250 | 169 | 419 | 16 | 23 | 39

Table 1 gives the computation times of a single step of the STELLOPT
Levenberg-Marquadt optimization algorithm for a range of phase space
dimensions and at different levels of parallelization. The number of
optimization parameters used corresponds to optimizing all of the stellarator-
symmetric boundary modes within the resolution $M,N\leq 1,2,3,4$, excluding
the major radius. Since the bottleneck process of the STELLOPT algorithm is
computing gradients through finite differences at each optimization step, its
computation time is expected to increase linearly with the number of
optimization parameters, and decrease linearly with the number of parallel
CPUs that are sharing the workload. The timing data in this table confirm that
scaling – the longer run times of larger optimization problems can be reduced
through greater parallelization. However, for a moderately sized problem on a
demanding number of state-of-the-art CPU cores, STELLOPT still takes several
minutes per iteration, which results in hours of total computation time.

In contrast, table 2 gives the computation times for a single Gauss-Newton
trust-region optimization step in DESC, broken down between the perturbation
and equilibrium solve involved at each iteration. Only the first-order
optimization method was tested to give the closest comparison to STELLOPT, but
including the second-order terms does not add significantly to the
perturbation times as explained previously. The equilibrium solve time is
heavily dependent on the number of Gauss-Newton steps required to solve (6) to
the desired tolerances. The average times shown in this table are for ten
iterations of the equilibrium subproblem, but in practice the number will vary
between optimizer iterations (larger optimization steps typically require more
force balance correction). Thanks to the use of automatic differentiation, the
computation times are independent of the number of optimization parameters –
these results included all $560$ boundary modes in the $M,N\leq 8$ equilibrium
resolution. Because the objective functions get compiled to work with
automatic differentiation the first time they are called, the initial
iteration takes substantially longer than subsequent ones. The benefit is that
calls to the compiled functions can then be evaluated much faster, yielding
large savings in the long run when many optimization steps are required. These
functions only need to be recompiled if the resolution of the equilibrium
changes, which typically does not occur in a traditional optimization
procedure. Even with a high level of parallelization in STELLOPT, the reliance
on finite differencing prevents it from approaching the speed of automatic
differentiation for large numbers of optimization parameters.

This table also reveals a large speedup of running on a GPU: the compile times
are roughly halved and iterations are about $8$ times faster than on the CPU
for this problem. It is important to note that the same code can run on either
a CPU or GPU without any modifications required by the user. Special attention
is needed for discussion of the quasi-symmetry objective $f_{B}$: because the
conversion to Boozer coordinates requires re-evaluating Fourier basis
functions in the new coordinates, differentiating through this function with
automatic differentiation is very memory intensive. The numerical resolution
used in these tests was too large to run the Boozer form on the GPU with
limited memory capacity. The perturbation time of the Boozer form on the CPU
is also much slower than the other quasi-symmetry objectives, but it is still
a significant improvement over STELLOPT especially when considering the number
of CPUs and optimization parameters involved. Whether running on a single CPU
or GPU, the total run time is orders of magnitude faster using DESC compared
to STELLOPT.

##  4 Conclusions

DESC is designed to be the next-generation stellarator optimization code. Its
unique architecture allows it to take advantage of automatic differentiation
and seamlessly pass that derivative information from its equilibrium solver to
the optimizer, resulting in dramatic speed improvements over existing tools.
This paper has explained the DESC optimization approach and demonstrated the
simplicity of using the python code. Through examples of quasi-symmetry
optimization, it has shown that DESC can achieve similar if not better quality
results as conventional tools such as STELLOPT in orders of magnitude less
computation time, and even other modern tools like SIMSOPT are not designed
for this level of speed performance. The triple product metric is believed to
be the best quasi-symmetry objective due to its lack of physics assumptions
and its numerical advantages, but more test cases are needed to confirm if
this result generalizes to other problems. The purpose has been to showcase
the potential of DESC, and discovering reactor-relevant stellarator designs is
left for future work.

As an open source software with a modular structure, the DESC code suite is
continually improving upon existing performance and expanding to provide new
functionality. The example quasi-symmetry optimizations used default options,
but there is great flexibility in specifying alternatives to better
accommodate certain problems. Quasi-symmetry was provided as an example
optimization objective due to its popularity in the literature, but other
physics and engineering objectives should be added and used in tandem. Thanks
to automatic differentiation, the process of incorporating new optimization
targets with exact derivative information is relatively simple for the
developer. There is also room to add alternatives to the Gauss-Newton trust-
region optimization method, such as global optimization algorithms. More work
is needed to better understand the trade-offs between robustness and
convergence speed among the various possibilities. Plans to extend DESC to
handle free-boundary equilibria and single-stage coil optimization are already
under development. Anyone who is interested in taking advantage of these
capabilities is strongly encouraged to install and use the publicly available
code, and recommendations for additional features are always welcome.

Due to its efficient computations, DESC is well positioned to perform a large-
scale exploration of the stellarator design space. The perturbation and
continuation methods could provide physical insight on the bifurcation of
tokamaks into quasi-axisymmetric stellarators (Plunk, 2020). DESC could also
be used to generate a large database of equilibria that could then be analyzed
with the aid of machine learning techniques to identify regions in the
parameter landscape that deserve closer attention. Another important
application is equilibrium reconstruction, which is an optimization problem to
find the equilibrium parameters that best match experimental diagnostics. The
rate of generating full equilibria reconstructions with conventional tools is
disappointingly slow (Hanson et al., 2009), but DESC could be used to provide
timely analysis of experiments. Whether improving the usefulness of existing
devices or discovering better designs for future reactors, the optimization
capabilities of DESC open new possibilities for stellarator research that were
previously unavailable.

Thank you to Eduardo Rodriguez for many insightful conversations about the
nature of quasi-symmetry and its potential for stellarator optimization, and
to Albert Mollen and Elizabeth Paul for their instruction on the SFINCS code.
Thank you to the SIMSOPT development team, whose python version of BOOZ_XFORM
was the basis for the implementation in DESC.

This work was supported by the U.S. Department of Energy under contract
numbers DE-AC02-09CH11466, DE-SC0022005 and Field Work Proposal No. 1019. The
United States Government retains a non-exclusive, paid-up, irrevocable, world-
wide license to publish or reproduce the published form of this manuscript, or
allow others to do so, for United States Government purposes.

The data that support the findings of this study are available from the
corresponding author upon reasonable request. The latest version of the DESC
code can be accessed at the public repository:
https://github.com/PlasmaControl/DESC. The results in this paper were
generated from the git commit hash 0fcc708, and the input files used are
available with the publication.

## References

  * Antonsen et al. (2019) Antonsen, T., Paul, E. J. & Landreman, M. 2019 J. Plasma Phys. 85, 905850207.
  * Boozer (2015) Boozer, A. H. 2015 J. Plasma Phys. 81, 515810606\.
  * Bradbury et al. (2018) Bradbury, James, Frostig, Roy, Hawkins, Peter, Johnson, Matthew James, Leary, Chris, Maclaurin, Dougal, Necula, George, Paszke, Adam, VanderPlas, Jake, Wanderman-Milne, Skye & Zhang, Qiao 2018 JAX: composable transformations of Python+NumPy programs.
  * Conlin et al. (2022) Conlin, Rory, Dudt, Daniel W., Panici, Dario & Kolemen, Egemen 2022 The desc stellarator code suite part ii: Perturbation and continuation methods .
  * Drevlak et al. (2019) Drevlak, M., Beidler, C.D., Geiger, J., Helander, P. & Turkin, Y. 2019 Optimisation of stellarator equilibria with ROSE. Nuclear Fusion 59 (1), 016010.
  * Dudt & Kolemen (2020) Dudt, D. W. & Kolemen, E. 2020 Phys. Plasmas 27, 102513.
  * Dudt et al. (2022) Dudt, D. W., W., Conlin, D., Panici, Unalmis, K., Kim, P. & Kolemen, E. 2022 Desc. https://github.com/PlasmaControl/DESC.
  * Garren & Boozer (1991a) Garren, D. A. & Boozer, A. H. 1991a Phys. Fluids B 3, 2805–2821.
  * Garren & Boozer (1991b) Garren, D. A. & Boozer, A. H. 1991b Phys. Fluids B 3, 2822–2834.
  * Giuliani et al. (2020) Giuliani, Andrew, Wechsung, Florian, Cerfon, Antoine, Stadler, Georg & Landreman, Matt 2020 Single-stage gradient-based stellarator coil design: Optimization for near-axis quasi-symmetry. arXiv:2010.02033 [physics] , arXiv: 2010.02033.
  * Hanson et al. (2009) Hanson, J. D., Hirshman, S. P., Knowlton, S. F., Lao, L. L., Lazarus, E. A. & Shields, J. M. 2009 Nucl. Fusion 49, 075031.
  * Hirshman & Whitson (1983) Hirshman, S. P. & Whitson, J. C. 1983 Phys. Fluids 26, 3553–3568.
  * Kruskal & Kulsrud (1958) Kruskal, M. D. & Kulsrud, R. M. 1958 Phys. Fluids 1, 265–274.
  * Landreman (2022) Landreman, M. 2022 stellopt_scenarios. https://github.com/landreman/stellopt_scenarios/tree/master/34DOF_varyBoundary_targetQuasisymmetry.
  * Landreman et al. (2021) Landreman, M., Medasani, B., Wechsung, F., Giuliani, A., Jorge, R. & Zhu, C. 2021 Simsopt: A flexible framework for stellarator optimization. Journal of Open Source Software 6 (65), 3525.
  * Landreman & Sengupta (2018) Landreman, M. & Sengupta, W. 2018 J. Plasma Phys. 84, 905840616.
  * Landreman et al. (2019) Landreman, M., Sengupta, W. & Plunk, G. G. 2019 J. Plasma Phys. 85, 905850103.
  * Landreman et al. (2014) Landreman, M., Smith, H. M., Mollen, A. & Helander, P. 2014 Phys. Plasmas 21, 042503.
  * Martin & Landreman (2020) Martin, M. F. & Landreman, M. 2020 J. Plasma Phys. 86, 905860317.
  * McGreivy et al. (2021) McGreivy, N., Hudson, S. R. & Zhu, C. 2021 Nucl. Fusion 61, 026020.
  * Nemov et al. (1999) Nemov, V. V., Kasilov, S. V., Kernbichler, W. & Heyn, M. F. 1999 Phys. Plasmas 6, 4622–4632.
  * Nocedal & Wright (2006) Nocedal, J. & Wright, S. J. 2006 Numerical Optimization, 2nd edn. Springer.
  * Nührenberg & Zille (1988) Nührenberg, J. & Zille, R. 1988 Phys. Lett. A 129, 113–117.
  * Panici et al. (2022) Panici, Dario, Conlin, Rory, Dudt, Daniel W. & Kolemen, Egemen 2022 The desc stellarator code suite part i: Quick and accurate equilibria computations .
  * Paul et al. (2020) Paul, E. J., Antonsen, T., Landreman, M. & Cooper, W. A. 2020 J. Plasma Phys. 86, 905850207.
  * Paul et al. (2021) Paul, E. J., Landreman, M. & Antonsen, T. 2021 J. Plasma Phys. 87, 905870214.
  * Plunk (2020) Plunk, G. G 2020 J. Plasma Phys. 86, 905860409\.
  * Plunk et al. (2019) Plunk, G. G., Landreman, M. & Helander, P. 2019 J. Plasma Phys. 85, 905850602.
  * Rodriguez & Bhattacharjee (2021) Rodriguez, E. & Bhattacharjee, A. 2021 Phys. Plasmas 28, 012508.
  * Rodriguez et al. (2020) Rodriguez, E., Helander, P. & Bhattacharjee, A. 2020 Phys. Plasmas 27, 062501.
  * Rodriguez et al. (2022) Rodriguez, E., Paul, E. J. & Bhattacharjee, A. 2022 J. Plasma Phys. 88, 905880109.
  * Spong et al. (1998) Spong, D. A., Hirshman, S. P., Whitson, J. C., Batchelor, D. B., Carreras, B. A., Lynch, V. E. & Rome, J. A. 1998 J* optimization of small aspect ratio stellarator/tokamak hybrid devices. Physics of Plasmas 5 (5), 1752–1758.
  * Yu & Dutton (1998) Yu, Z. & Dutton, R. W. 1998 VLSI Design 6, 141–145.