$p=2$. One could also remove the redundancy in the definition of $\theta$ by
using the unique and spectrally condensed representation of Hirshman and
Breslau [109] or by solving the nonlinear constraint equation of Hirshman and
Meier [110] once the optimal surface has been obtained.

The quantity
$\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}=\sqrt{\chi^{2}_{J}/A_{\text{coil}}}$
is the 2-norm of the current density, where $A_{\text{coil}}$ is the winding
surface area,

$\displaystyle A_{\text{coil}}=\int_{\text{coil}}d^{2}x\,.$ (3.13)

Although we are using a current potential approach rather than directly
optimizing coil shapes, including
$\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$ in the objective
function allows us to obtain coils with good engineering properties.
Derivatives of coil-specific metrics (such as curvature) could be computed
from the current potential if desired. For example, consider $N$ contours
beginning at equally-spaced toroidal angles $\phi_{0}^{i}$ and $\theta_{0}=0$.
The $i^{\text{th}}$ contour is defined by functions $\theta_{i}(s)$ and
$\phi_{i}(s)$ for parameter $s$, where $\partial\Phi/\partial s=0$. The
derivatives of coil metrics which depend on
$\textbf{x}(\theta_{i}(s),\phi_{i}(s))$, could be computed with the adjoint
method which will be described in Section 3.4. As the direct targeting of coil
metrics introduces additional arbitrary weights in the objective function and
the solution to another adjoint equation must be obtained to compute its
gradient, we instead include
$\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$ in our objective
function.

To demonstrate this correlation between
$\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$ and coil shape
complexity, we compute the coil set on the actual W7-X winding surface using
REGCOIL. The regularization parameter $\lambda$ is varied to achieve several
values of $\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$. Coil shapes
are obtained from the contours of $\Phi$. In Figure 3.1, two of the W7-X non-
planar coils computed in this way are shown, and the corresponding coil
metrics are given in Table 3.1. (These correspond to the two leftmost coils in
Figure 3.5.) We consider the average and maximum length $l$, toroidal extent
$\Delta\phi$, curvature $\kappa$, and the minimum coil-coil distance
$d_{\text{coil-coil}}^{\text{min}}$. The average, maximum, and minimum are
taken over the set of 5 unique coils. The coil shapes become more complex as
$\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$ increases, quantified
by increasing $\kappa$ and $\Delta\phi$ and decreasing $d_{\text{coil-
coil}}^{\text{min}}$. Here the curvature, $\kappa$, of a three-dimensional
parameterized curve, $\textbf{x}(t)$, is,

$\displaystyle\kappa=\frac{\left\rvert\textbf{x}^{\prime}(t)\times\textbf{x}^{\prime\prime}(t)\right\rvert}{\left\rvert\textbf{x}^{\prime}(t)\right\rvert^{3}}.$
(3.14)

We have compared coil shapes on a single winding surface, finding them to
become simpler as $\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$
decreases. As
$\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}=\left(\chi^{2}_{J}/A_{\text{coil}}\right)^{1/2}$,
we would find similar trends with $\chi^{2}_{J}$. We have chosen to include
$\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$ in the objective
function as it is normalized by $A_{\text{coil}}$, so it is a more useful
quantity for comparison of coil shapes on different winding surfaces.

(a)

(b)

(c)

Figure 3.1: Two non-planar W7-X coils (corresponding to the two leftmost coils in Figure 3.5) computed with REGCOIL using the actual W7-X winding surface. The regularization parameter $\lambda$ is chosen to achieve the shown values of $\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$. As $\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$ increases, the average length, toroidal extent, and curvature increase. Figure adapted from [185] with permission. $\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$ [MA/m] | 2.20 | 2.70 | 3.20
---|---|---|---
$J_{\text{max}}$ [MA/m] | 4.55 | 9.50 | 29.1
$\chi^{2}_{B}$ [T2 m2] | 1.89 | $5.25\times 10^{-3}$ | $2.10\times 10^{-3}$
Average $l$ [m] | 8.03 | 9.18 | 9.81
Max $l$ [m] | 8.26 | 10.5 | 11.8
Average $\Delta\phi$ [rad.] | 0.146 | 0.222 | 0.253
Max $\Delta\phi$ [rad.] | 0.161 | 0.282 | 0.372
Average $\kappa$ [m-1] | 1.04 | 1.29 | 1.32
Max $\kappa$ [m-1] | 2.54 | 20.3 | 56.1
$d_{\text{coil-coil}}^{\text{min}}$ [m] | 0.353 | 0.182 | 0.0758
Table 3.1: Comparison of metrics for coils computed with REGCOIL using the
actual W7-X winding surface. Average and max are evaluated for the set of 5
unique coils. The regularization parameter $\lambda$ is varied to achieve
these values of $\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$. Table
adapted from [185] with permission.

To minimize $f$, the relative weights in (3.10) ($\alpha_{V}$,
$\alpha_{\mathcal{S}}$, and $\alpha_{J}$) are chosen such that each of the
terms in the objective function have similar magnitudes, though much tuning of
these parameters is required to obtain results which simultaneously improve
the physics properties (decrease $\chi^{2}_{B}$) and engineering properties
(increase $V_{\text{coil}}$ and $d_{\text{coil-coil}}^{\text{min}}$, decrease
$\kappa$ and $\Delta\phi$).

####  3.3.2 Optimization constraints

Minimization of $f$ is performed subject to the inequality constraint
$d_{\text{min}}\geq d_{\text{min}}^{\text{target}}$. Here $d_{\text{min}}$ is
the minimum distance between the coil-winding surface and the plasma surface,

$\displaystyle d_{\text{min}}=\min_{\theta,\phi}\left(d_{\text{coil-
plasma}}\right)=\min_{\theta,\phi}\left(\min_{\theta_{p},\phi_{p}}\,\mathinner{\\!\left\lvert\textbf{x}_{C}-\textbf{x}_{P}\right\rvert}\right),$
(3.15)

and $d_{\text{min}}^{\text{target}}$ is the minimum tolerable coil-plasma
separation. The quantities $\theta_{p}$ and $\phi_{p}$ are poloidal and
toroidal angles on the plasma surface, $\textbf{x}_{P}$ and $\textbf{x}_{C}$
are the position vectors on the plasma and winding surface, and
$d_{\text{coil-plasma}}$ is the coil-plasma distance as a function of $\theta$
and $\phi$.

The maximum current density $J_{\text{max}}$ is also constrained,

$\displaystyle J_{\text{max}}=\max_{\theta,\phi}\,J.$ (3.16)

This roughly corresponds to a fixed minimum coil-coil spacing. This constraint
is enforced by fixing $J_{\text{max}}$ to obtain the regularization parameter
$\lambda$ in the REGCOIL solve, so we avoid the need for an equality
constraint or the inclusion of $J_{\text{max}}$ in the objective function.
Rather, $\overrightarrow{\bm{\Phi}}(\Omega)$ is determined such that
$J_{\text{max}}$ is fixed. The inequality-constrained nonlinear optimization
is performed using the NLOPT [125] software package using a conservative
convex separable quadratic approximation (CCSAQ) [224]. While there are
several gradient-based inequality-constrained algorithms available, we choose
to use CCSAQ as it is relatively insensitive to the bound constraints imposed
on the optimization parameters. We recognize that there are many possible
combinations of constraints, objective functions, and regularization
conditions that could be used. For example,
$\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$ could be fixed to
determine $\lambda$ while $J_{\text{max}}$ could be included in the objective
function. We found that the formulation we have presented produces the best
coil shapes.

###  3.4 Derivatives of $f$ and the adjoint method

We must compute derivatives of $f$ with respect to the geometric parameters
$\Omega$ in order to use gradient-based optimization methods. The spectral
width $\mathcal{S}$ and the volume $V_{\text{coil}}$ are explicit functions of
$\Omega$, so their analytic derivatives can be obtained. On the other hand,
$\chi^{2}_{B}$ and $\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$
depend both explicitly on coil geometry and on $\bm{\Phi}(\Omega)$. One
approach to obtain the derivatives of these quantities could be to solve the
REGCOIL linear system $N_{\Omega}+1$ times, taking a finite-difference step in
each Fourier coefficient. However, if $N_{\Omega}$ is large, the computational
cost of this method could be prohibitively expensive. Instead, we will apply
the adjoint method to compute derivatives. This technique will be demonstrated
below.

The derivative of $\chi^{2}_{B}$ can be computed using the chain rule,

$\displaystyle\frac{\partial\chi^{2}_{B}(\Omega,\overrightarrow{\bm{\Phi}}(\Omega))}{\partial\Omega_{m,n}}=\frac{\partial\chi^{2}_{B}(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\Omega_{m,n}}+\frac{\partial\chi^{2}_{B}(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\overrightarrow{\bm{\Phi}}}\cdot\frac{\partial\overrightarrow{\bm{\Phi}}(\Omega)}{\partial\Omega_{m,n}},$
(3.17)

where $\overrightarrow{\bm{\Phi}}(\Omega)$ is understood to vary with $\Omega$
such that (3.8) is satisfied. The dot product is a contraction over the
current potential basis functions, $\\{\Phi_{m,n}\\}$. We can compute
$\partial\overrightarrow{\bm{\Phi}}(\Omega)/\partial\Omega_{m,n}$ by
differentiating the linear system (3.8) with respect to $\Omega_{m,n}$,

$\displaystyle\frac{\partial\overleftrightarrow{\textbf{A}}(\Omega)}{\partial\Omega_{m,n}}\overrightarrow{\bm{\Phi}}+\overleftrightarrow{\textbf{A}}\frac{\partial\overrightarrow{\bm{\Phi}}(\Omega)}{\partial\Omega_{m,n}}=\frac{\partial\overrightarrow{\textbf{b}}(\Omega)}{\partial\Omega_{m,n}},$
(3.18)

and formally solving this equation to obtain,

$\displaystyle\frac{\partial\bm{\Phi}(\Omega)}{\partial\Omega_{m,n}}=\overleftrightarrow{\textbf{A}}^{-1}\left(\frac{\partial\overrightarrow{\textbf{b}}(\Omega)}{\partial\Omega_{m,n}}-\frac{\partial\overleftrightarrow{\textbf{A}}(\Omega)}{\partial\Omega_{m,n}}\overrightarrow{\bm{\Phi}}\right).$
(3.19)

Equation (3.19) is inserted into (3.17),

$\displaystyle\frac{\partial\chi^{2}_{B}(\Omega,\overrightarrow{\bm{\Phi}}(\Omega))}{\partial\Omega_{m,n}}=\frac{\partial\chi^{2}_{B}(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\Omega_{m,n}}+\frac{\partial\chi^{2}_{B}(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\overrightarrow{\bm{\Phi}}}\cdot\left[\overleftrightarrow{\textbf{A}}^{-1}\left(\frac{\partial\overrightarrow{\textbf{b}}(\Omega)}{\partial\Omega_{m,n}}-\frac{\partial\overleftrightarrow{\textbf{A}}(\Omega)}{\partial\Omega_{m,n}}\overrightarrow{\bm{\Phi}}\right)\right].$
(3.20)

This expression could be evaluated by solving the linear system (3.18) for
$\partial\overrightarrow{\bm{\Phi}}/\partial\Omega_{m,n}$ and performing the
inner product with $\partial\chi^{2}_{B}/\partial\overrightarrow{\bm{\Phi}}$.
However, the computational cost of this method scales similarly to that of
finite differencing, as described in Section 2.2.1. Instead, we can exploit
the adjoint property of the operator to obtain,

$\displaystyle\frac{\partial\chi^{2}_{B}(\Omega,\overrightarrow{\bm{\Phi}}(\Omega))}{\partial\Omega_{m,n}}=\frac{\partial\chi^{2}_{B}(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\Omega_{m,n}}+\left[\left(\overleftrightarrow{\textbf{A}}^{-1}\right)^{T}\frac{\partial\chi^{2}_{B}(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\overrightarrow{\bm{\Phi}}}\right]\cdot\left(\frac{\partial\overrightarrow{\textbf{b}}(\Omega)}{\partial\Omega_{m,n}}-\frac{\partial\overleftrightarrow{\textbf{A}}(\Omega)}{\partial\Omega_{m,n}}\overrightarrow{\bm{\Phi}}\right).$
(3.21)

For any invertible matrix,
$\left(\overleftrightarrow{\textbf{A}}^{-1}\right)^{T}=\left(\overleftrightarrow{\textbf{A}}^{T}\right)^{-1}$.
Hence we can instead solve a linear system involving the matrix
$\overleftrightarrow{\textbf{A}}^{T}$ to compute an adjoint variable
$\overrightarrow{\textbf{q}}$, defined as the solution of

$\displaystyle\overleftrightarrow{\textbf{A}}^{T}\overrightarrow{\textbf{q}}=\frac{\partial\chi^{2}_{B}(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\overrightarrow{\bm{\Phi}}}.$
(3.22)

Rather than compute a finite-difference derivative for each $\Omega_{m,n}$ or
solve a linear system to compute each
$\partial\overrightarrow{\bm{\Phi}}/\partial\Omega_{m,n}$ as in (3.19), we
solve two linear systems: the forward (3.8) and adjoint (3.22). The adjoint
equation is similar to the forward equation
($\overleftrightarrow{\textbf{A}}^{T}$ has the same dimensions and
eigenspectrum as $\overleftrightarrow{\textbf{A}}$), so the same computational
tools can be used to solve the adjoint problem. We then perform an inner
product with $\overrightarrow{\textbf{q}}$ to obtain the derivatives with
respect to each $\Omega_{m,n}$,

$\displaystyle\frac{\partial\chi^{2}_{B}(\Omega,\overrightarrow{\bm{\Phi}}(\Omega))}{\partial\Omega_{m,n}}=\frac{\partial\chi^{2}_{B}(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\Omega_{m,n}}+\overrightarrow{\textbf{q}}\cdot\left(\frac{\partial\overrightarrow{\textbf{b}}(\Omega)}{\partial\Omega_{m,n}}-\frac{\partial\overleftrightarrow{\textbf{A}}(\Omega)}{\partial\Omega_{m,n}}\bm{\Phi}\right).$
(3.23)

The derivatives $\partial\overrightarrow{\textbf{b}}/\partial\Omega_{m,n}$,
$\partial\overleftrightarrow{\textbf{A}}/\partial\Omega_{m,n}$,
$\partial\chi^{2}_{B}/\partial\Omega_{m,n}$, and
$\partial\chi^{2}_{B}/\partial\overrightarrow{\bm{\Phi}}$ can be computed
analytically. In the above discussion, the regularization parameter $\lambda$
has been assumed to be fixed. A similar method can be used if a $\lambda$
search is performed to obtain a target $J_{\text{max}}$ (see Appendix C). The
same method is used to compute derivatives of
$\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$.

We note that adjoint methods provide the most significant reduction in
computational cost when the linear solve is expensive. For the REGCOIL system,
this is not the case, as the cost of constructing
$\overleftrightarrow{\textbf{A}}$ and $\overrightarrow{\textbf{b}}$ exceeds
that of the solve. We have implemented OpenMP multithreading for the
construction of $\partial\overleftrightarrow{\textbf{A}}/\partial\Omega$ and
$\partial\overrightarrow{\textbf{b}}/\partial\Omega$ such that the cost of
computing the gradients via the adjoint method is cheaper than computing
finite-difference derivatives serially.

The constraint functions, $d_{\text{min}}$ and $J_{\text{max}}$, must also be
differentiated with respect to $\Omega_{m,n}$. As $d_{\text{min}}$ is defined
in terms of the minimum function, we approximate it using the smooth log-sum-
exponent function [29],

$\displaystyle d_{\text{min,
lse}}=-\frac{1}{q}\log\left(\frac{\int_{S_{C}}d^{2}x_{C}\,\int_{S_{P}}d^{2}x_{P}\,\exp\left(-q\mathinner{\\!\left\lvert\textbf{x}_{C}-\textbf{x}_{P}\right\rvert}\right)}{\int_{S_{C}}d^{2}x_{C}\,\int_{S_{P}}d^{2}x_{P}\,}\right).$
(3.24)

This function can be analytically differentiated with respect to
$\Omega_{m,n}$. As $q$ approaches infinity, $d_{\text{min, lse}}$ approaches
$d_{\text{min}}$. For $q$ very large, the function obtains very sharp
gradients. A typical value of $q=10^{4}$ m-1 was used. The log-sum-exponent
function is also used to approximate $J_{\text{max}}$, as described in
Appendix C.

###  3.5 Winding surface optimization results

####  3.5.1 Trends with optimization parameters

Beginning with the actual W7-X winding surface, we perform scans over the
coefficients $\alpha_{V}$ and $\alpha_{\mathcal{S}}$ in the objective function
(3.10). The plasma surface was obtained from a fixed-boundary VMEC solution
that predated the coil design and is free from modular coil ripple. The
constraint target is set to be the minimum coil-plasma distance on the initial
winding surface, $d_{\text{min}}^{\text{target}}=0.37$ m. The cross-sections
of the optimized surfaces in the poloidal plane are shown in Figures 3.2 and
3.3 along with the last-closed flux surface (red), a constant offset surface
at $d_{\text{min}}^{\text{target}}$ (black solid), and the initial winding
surface (black dashed).

We perform a scan over $\alpha_{\mathcal{S}}$ with $\alpha_{V}=\alpha_{J}=0$.
For optimal values of $\alpha_{\mathcal{S}}$, the addition of the spectral
width term should simply reparameterize the surface, eliminating the zero-
gradient direction in parameter space. Thus we expect that when $\chi^{2}_{B}$
is the only other term in the objective function, the winding surface should
collapse to a constant offset surface. When $\alpha_{\mathcal{S}}$ is too
large, the surface shape changes to favor a condensed Fourier series. When
$\alpha_{\mathcal{S}}$ is too small, the optimization may terminate
prematurely in a local minimum due to the non-uniqueness of the
representation. Indeed we find that with increasing $\alpha_{\mathcal{S}}$,
the winding surface approaches a torus with a circular cross-section, which
has a minimal Fourier spectrum. At moderately small values of
$\alpha_{\mathcal{S}}$ ($\sim 0.3$) the surface approaches a constant offset
surface at $d_{\text{min}}^{\text{target}}$, as $\chi^{2}_{B}$ is dominant in
objective function. For very small values of $\alpha_{\mathcal{S}}$ ($\sim
0.003$), we find that the optimization terminates at a point relatively close
to the initial surface, and the resulting winding surface deviates from a
constant offset surface. An intermediate value of $\alpha_{\mathcal{S}}=0.3$
was chosen for the following optimizations of the W7-X winding surface.

A scan over $\alpha_{V}$ is performed at fixed $\alpha_{\mathcal{S}}=0.3$ and
$\alpha_{J}=0$ such that the spectral width does not greatly increase. As
$\alpha_{V}$ increases, $d_{\text{coil-plasma}}$ increases significantly on
the outboard side while it remains fixed in the inboard concave regions. This
trend is not surprising, as concave plasma shapes have been shown to be
inefficient to produce with coils [137]. Interestingly, the winding surface
obtains a somewhat pointed shape at the triangle cross-section ($\phi=0.5$
$2\pi/N_{p}$), becoming elongated at the tip of the triangle and “pinching”
toward the plasma surface at the edges.

Figure 3.2: Optimized winding surfaces obtained with $\alpha_{V}=\alpha_{J}=0$
and the values of $\alpha_{\mathcal{S}}$ shown. The actual W7-X winding
surface is used as the initial surface in the optimization (black dashed). As
$\alpha_{\mathcal{S}}$ increases, the magnitude of the spectral-width term in
the objective function increases, and the winding surface approaches a
cylindrical torus with a minimal Fourier spectrum. For moderately small values
of $\alpha_{\mathcal{S}}$, the winding surface approaches a uniform offset
surface from the plasma surface (black solid). Figure adapted from [185] with
permission. Figure 3.3: Optimized winding surfaces obtained with
$\alpha_{\mathcal{S}}=0.3$, $\alpha_{J}=0$, and the values of $\alpha_{V}$
shown. The actual W7-X winding surface is used as the initial surface in the
optimization (black dashed). As $\alpha_{V}$ increases, $d_{\text{coil-
plasma}}$ increases on the outboard side while it remains fixed in the concave
region. Figure adapted from [185] with permission.

####  3.5.2 Optimal W7-X winding surface

We now include nonzero $\alpha_{J}$ and attempt a comprehensive optimization.
The $J_{\text{max}}$ constraint is selected such that the metrics ($l$,
$\kappa$, and $\Delta\phi$) of the coils computed on the initial surface
roughly match those of the actual non-planar coil set. The coil-plasma
distance constraint $d_{\text{min}}^{\text{target}}$ is set to be the minimum
$d_{\text{coil-plasma}}$ on the initial winding surface. Parameters
$\alpha_{V}=0.5$, $\alpha_{\mathcal{S}}=0.24$, and $\alpha_{J}=1.6\times
10^{-6}$ were used in the objective function. Optimization was performed over
118 Fourier coefficients $\big{(}\mathinner{\\!\left\lvert n\right\rvert}\leq
4$ and $m\leq 6$ in (3.9)$\big{)}$ and the objective function was evaluated a
total of 5165 times to reach the optimum ($1.5\times 10^{4}$ linear solves
rather than $6.1\times 10^{5}$ required for finite-difference derivatives).
The optimal surface and coil set are shown in Figures 3.4 and 3.5, and the
corresponding metrics are shown in Table 3.2. We find a solution which
increases $V_{\text{coil}}$ by 22% and decreases $\chi^{2}_{B}$ by 52% over
the initial winding surface. (Note that it is numerically impossible to obtain
a current distribution that exactly reproduces the plasma surface, so
$\chi^{2}_{B}$ is nonzero when computed from the REGCOIL solution on the
initial winding surface.) In addition, the optimized coil set features a
smaller average and maximum $\Delta\phi$ and $\kappa$ and larger
$d_{\text{coil-coil}}^{\text{min}}$. The length of the coils increases to
accommodate for the increase in $V_{\text{coil}}$. Again we find that the
increase in $V_{\text{coil}}$ is most pronounced in the outboard convex
regions while $d_{\text{coil-plasma}}$ is maintained in the concave regions of
the bean-shaped cross-sections. The “pinching” feature of the winding surface
is again present in the triangle cross-section ($\phi=0.5\,2\pi/N_{p}$).

It should be noted that the decrease in $d_{\text{coil-plasma}}$ at the bottom
and top of the bean cross-section ($\phi=0$) might interfere with the current
W7-X divertor baffles. However, the increase in volume on the outboard side
would allow for increased flexibility for the neutral beam injection duct
[200]. We have performed this optimization to show that a winding surface
could be constructed that increases $V_{\text{coil}}$ (and thus the average
$d_{\text{coil-plasma}}$), improves coil shapes, and decreases $\chi^{2}_{B}$.
If further engineering considerations were necessary, these could be
implemented. The surface we have obtained is optimal with respect to the
engineering considerations and constraints we have imposed, which differ from
those of the W7-X team [15]. Thus the direct comparison between our method and
those of [15] cannot be made based on these results.

Figure 3.4: The actual W7-X coil-winding surface and plasma surface are shown with our optimized winding surface. In comparison with the actual surface, the optimized surface reduces $\chi^{2}_{B}$ by 52% and increases $V_{\text{coil}}$ by 22%. Figure adapted from [185] with permission. Figure 3.5: Comparisons of coil set computed with REGCOIL using the actual W7-X winding surface (dark blue) and the optimized surface (light blue). Figure reproduced from [185] with permission. | Initial | Optimized | Actual coil set
---|---|---|---
$\chi^{2}_{B}$ [T2m2] | 0.115 | 0.0711 |
$V_{\text{coil}}$[m3] | 156 | 190 |
$\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$ [MA/m] | 2.21 | 2.16 |
$J_{\text{max}}$ [MA/m] | 7.70 | 7.70 |
Average $l$ [m] | 8.51 | 8.95 | 8.69
Max $l$ [m] | 8.84 | 9.14 | 8.74
Average $\Delta\phi$ [rad.] | 0.190 | 0.179 | 0.198
Max $\Delta\phi$ [rad.] | 0.222 | 0.197 | 0.208
Average $\kappa$ [m-1] | 1.21 | 1.10 | 1.20
Max $\kappa$ [m-1] | 9.01 | 4.84 | 2.59
$d_{\text{coil-coil}}^{\text{min}}$ [m] | 0.223 | 0.271 | 0.261
Table 3.2: Comparison of metrics of the actual W7-X winding surface and our
optimized surface. We also show metrics of the coil set computed on the
winding surfaces using REGCOIL and the metrics for the actual W7-X nonplanar
coils. Regularization in REGCOIL is chosen such that the coil metrics computed
on the initial surface roughly match those of the actual coil set. Coil
complexity improves from the initial to the final surface (decreased average
and max $\Delta\phi$ and $\kappa$, increased $d_{\text{coil-
coil}}^{\text{min}}$). The average and max $l$ increases to allow for the
increase in $V_{\text{coil}}$. Table adapted from [185] with permission.

####  3.5.3 Optimal HSX winding surface

We perform the same procedure for the optimization of the HSX winding surface.
Parameters $\alpha_{V}=3.13\times 10^{-4}$, $\alpha_{\mathcal{S}}=0$, and
$\alpha_{J}=3\times 10^{-10}$ were used in the objective function. We found
that the spectral width term was not necessary to obtain a satisfying optimum
in this case. The initial winding surface was taken to be a toroidal surface
on which the actual modular coils lie. The plasma equilibrium used is a fixed-
boundary VMEC solution without coil ripple. Optimization was performed over
100 Fourier coefficients $\big{(}\mathinner{\\!\left\lvert n\right\rvert}\leq
5$ and $m\leq 4$ in (3.9)$\big{)}$ and the objective function was evaluated a
total of 560 times to reach the optimum ($1.7\times 10^{3}$ linear solves
rather than $5.7\times 10^{4}$ required for forward-difference derivatives).
The coil-plasma distance constraint was set to be
$d_{\text{min}}^{\text{target}}=0.14$ m, the minimum coil-plasma distance on
the actual winding surface. The optimal surface and coil set are shown in
Figures 3.6 and 3.7, and the corresponding coil metrics are shown in Table
3.3. We find a solution that increases $V_{\text{coil}}$ by 18% and decreases
$\chi^{2}_{B}$ by 4% over the initial winding surface. The coil set computed
with REGCOIL using the optimized surface appears qualitatively similar to that
computed with the initial surface but with increased $d_{\text{coil-plasma}}$
on the outboard side. The average and maximum $\Delta\phi$ and $\kappa$
decreased while $d_{\text{coil-coil}}^{\text{min}}$ was increased for the coil
set computed on the optimal surface in comparison to that of the initial
surface. As was observed in the W7-X optimization (Figure 3.4), the optimized
HSX winding surface obtains a somewhat pinched shape near the triangle cross-
section ($\phi=0.5\,2\pi/N_{p}$).

Figure 3.6: The actual HSX coil-winding surface and plasma surface are shown with our optimized winding surface. In comparison with the actual surface, the optimized surface has decreased $\chi^{2}_{B}$ by 4% and increased $V_{\text{coil}}$ by 18%. Figure adapted from [185] with permission. Figure 3.7: The coils obtained from REGCOIL using the actual HSX winding surface (dark blue) and optimized surface (light blue). Figure reproduced from [185] with permission. | Initial | Optimized | Actual coil set
---|---|---|---
$\chi^{2}_{B}$ [T2m2] | $1.53\times 10^{-5}$ | $1.47\times 10^{-5}$ |
$V_{\text{coil}}$[m3] | 2.60 | 3.07 |
$\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$ [MA/m] | 0.956 | 0.891 |
$J_{\text{max}}$ [MA/m] | 1.84 | 1.84 |
Average $l$ [m] | 2.26 | 2.39 | 2.24
Max $l$ [m] | 2.49 | 2.46 | 2.33
Average $\Delta\phi$ [rad.] | 0.372 | 0.365 | 0.362
Max $\Delta\phi$ [rad.] | 0.530 | 0.505 | 0.478
Average $\kappa$ [m-1] | 5.15 | 4.80 | 5.05
Max $\kappa$ [m-1] | 33.4 | 25.8 | 11.7
$d_{\text{coil-coil}}^{\text{min}}$ [m] | 0.0850 | 0.0853 | 0.0930
Table 3.3: Comparison of metrics of the actual HSX winding surface and our
optimized surface. We also show metrics of the coil set computed on the
winding surfaces using REGCOIL and the metrics for the actual HSX modular
coils. Regularization in REGCOIL is chosen such that the coil metrics computed
on the initial surface roughly match those of the actual coil set. Coil
complexity improves from the initial to the final surface (decreased average
and max $\Delta\phi$ and $\kappa$, increased $d_{\text{coil-
coil}}^{\text{min}}$). The average and max $l$ increases to allow for the
increase in $V_{\text{coil}}$. Table adapted from [185] with permission.

###  3.6 Local winding surface sensitivity

With the adjoint method we have computed derivatives of the objective function
with respect to Fourier components of the winding surface, $\partial
f/\partial\Omega$. While this representation of derivatives is convenient for
gradient-based optimization, the sensitivity to local displacements of the
surface is obscured. Alternatively, it is possible to represent the
sensitivity of $f$ with respect to normal displacements of surface area
elements of a given winding surface $S_{C}$,

$\displaystyle\delta
f(S_{C};\delta\textbf{x})=\int_{S_{C}}d^{2}x\,\,\mathcal{G}\delta\textbf{x}\cdot\hat{\textbf{n}}.$
(3.25)

The shape gradient and shape derivatives are described in detail in Section
2.1. As both $\chi^{2}_{B}$ and
$\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$ are defined in terms of
surface integrals over the winding surface, it can be shown that the shape
derivative of these functions can be written in the Hadamard form [171]. The
shape gradients $\mathcal{G}_{\chi^{2}_{B}}$ and
$\mathcal{G}_{\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}}$ can be
computed from the Fourier derivatives ($\partial\chi^{2}_{B}/\partial\Omega$
and
$\partial\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}/\partial\Omega$)
using a singular value decomposition method [138]. Here the perturbations
$\delta f$ and $\delta\textbf{x}$ are written in terms of the Fourier
derivatives, and $\mathcal{G}$ is also represented in a finite Fourier series,

$\displaystyle\frac{\partial
f(\Omega)}{\partial\Omega_{m,n}}=\int_{S_{C}}d^{2}x\,\left(\sum_{m,n}\mathcal{G}_{m,n}\cos(m\theta+nN_{p}\phi)\right)\frac{\partial\textbf{x}(\Omega)}{\partial\Omega_{m,n}}\cdot\hat{\textbf{n}}.$
(3.26)

After discretizing in $\theta$ and $\phi$, (3.26) takes the form of a
(generally not square) matrix equation which can be solved using the Moore-
Penrose pseudoinverse to obtain $\mathcal{G}_{m,n}$.

We compute $\mathcal{G}_{\chi^{2}_{B}}$ and
$\mathcal{G}_{\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}}$ (Figure
3.9) at fixed $\lambda$. These quantities are computed on the actual W7-X
winding surface and a surface uniformly offset from the plasma surface with
$d_{\text{coil-plasma}}=0.61$ m (the area-averaged $d_{\text{coil-plasma}}$
over the actual surface). We consider surfaces that are equidistant from the
plasma surface on average as $\mathcal{G}$ scales inversely with
$A_{\text{coil}}$. The poloidal cross-sections of these surfaces are shown in
Figure 3.8. For each surface $\lambda$ is chosen to achieve
$J_{\text{max}}=7.7$ MA/m as was used in Section 3.5.2. On both surfaces we
observe a narrow region featuring a large positive
$\mathcal{G}_{\chi^{2}_{B}}$, indicating that $d_{\text{coil-plasma}}$ should
decrease at that location in order that $\chi^{2}_{B}$ decreases. This
corresponds to locations on the plasma surface with significant concavity
(Figure 3.11(b)). The maximum $\mathcal{G}_{\chi^{2}_{B}}$ occurs at
$\phi=0.15\,2\pi/N_{p}$ on both surfaces (Figure 3.4). In comparison with this
region, the magnitude of $\mathcal{G}_{\chi^{2}_{B}}$ is relatively small over
the majority of the area of the surfaces shown, demonstrating that engineering
tolerances might be more relaxed in these locations. There is also a region of
negative $\mathcal{G}_{\chi^{2}_{B}}$ near $\phi=\pi/N_{p}$ and $\theta=0$.
This is the “tip” of the triangle-shaped cross-section, where $d_{\text{coil-
plasma}}$ was increased over the course of the optimization (Figures 3.2, 3.3,
and 3.4). We find that $\mathcal{G}_{\chi^{2}_{B}}$ computed on the actual
winding surface has similar trends to that computed on the surface uniformly
offset from the plasma. This indicates that the shape gradient depends on the
specific geometry of the winding surface. We have computed
$\mathcal{G}_{\chi^{2}_{B}}$ for several other winding surfaces with varying
$d_{\text{coil-plasma}}$. Regardless of the winding surface chosen, we observe
increased sensitivity in the concave regions.

The quantity
$\mathcal{G}_{\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}}$ roughly
quantifies how coil complexity changes with normal displacements of the coil
surface. In view of Figure 3.10, the locations of large
$\mathcal{G}_{\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}}$ overlap
with areas of increased $J$. On the actual winding surface, the maximum of
$\mathcal{G}_{\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}}$ occurs
near the location of the closest approach between coils (two rightmost coils
in Figure 3.5(a)). The shape gradients
$\mathcal{G}_{\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}}$ and
$\mathcal{G}_{\chi^{2}_{B}}$ have very similar trends. The concave regions of
the plasma surface are difficult to produce with external coils, resulting in
increased coil complexity and $J$. Therefore,
$\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$ is most sensitive to
displacements of the coil-winding surface in these regions.

Figure 3.8: The cross-sections of the two winding surfaces used to compute
$\mathcal{G}_{\chi^{2}_{B}}$ and
$\mathcal{G}_{\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}}$ are shown
in the poloidal plane. Figure adapted from [185] with permission.

(a) Offset from plasma

(b) Actual

(c) Offset from plasma

(d) Actual

Figure 3.9: Shape gradient for $\chi^{2}_{B}$ ((a) and (b))
$||\textbf{J}||_{2}$ ((c) and (d)). These functions are computed using the
W7-X plasma surface and a uniform offset winding surface from the plasma
surface with $d_{\text{coil-plasma}}=0.61$ m ((a) and (c)) and the actual
winding surface ((b) and (d)). The region of increased
$\mathcal{G}_{\chi^{2}_{B}}$ corresponds with concave regions of the plasma
surface (Figure 3.11(b)). Regions of large positive
$\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$ correspond to regions
with increased $J$ (Figure 3.10). Figure adapted from [185] with permission.

(a) Offset from plasma

(b) Actual

Figure 3.10: Current density magnitude, $J$, computed from REGCOIL using the
W7-X plasma surface and (a) a uniform offset winding surface from the plasma
surface with $d_{\text{coil-plasma}}=0.61$ m and (b) the actual winding
surface. Figure adapted from [185] with permission.

We recognize several ways that the shape gradient technique could be improved
to provide more relevant diagnostics for experimental design. With a winding
surface representation, the shape gradient does not allow for calculation of
the sensitivity to lateral coil displacements. Also, our analysis does not
account for field ripple due to the finite number of coils. Although Figure
3.9 indicates that the coils should move toward the plasma to reduce the field
error, the ripple fields might be significant with a filamentary model. A
similar calculation could be performed using the filamentary coil sensitivity
techniques presented in Section 2.1 and discussed further in Chapter 5.
Finally, $\chi^{2}_{B}$ does not account for the sensitivity to resonant
fields that could cause the formation of islands, though there is ongoing work
toward computing the shape gradient for such a metric [76].

Sensitivity studies on NCSX similarly found that coil errors on the inboard
side in regions of small $d_{\text{coil-plasma}}$ had a significant effect on
flux surface quality [236]. The necessity of small $d_{\text{coil-plasma}}$
for bean-shaped plasmas has been noted in many coil optimization efforts [216,
60] and has been demonstrated by evaluating the singular value decomposition
of the discretized Biot-Savart integral operator [137]. We can identify these
regions where the fidelity of the plasma surface requires tighter tolerance on
coil positions using the shape gradient.

###  3.7 Metrics for configuration optimization

The results presented here and in [137] indicate that the concave regions of
the surface are both the regions where a small coil-plasma distance is
required and the sensitivity to the winding surface position is highest. The
regions of concavity can be determined by considering the principal curvatures
of the plasma surface. Let $\hat{\textbf{n}}(\textbf{x}_{0})$ represent the
normal vector at the plasma surface at some point $\textbf{x}_{0}$, and let
$A_{n}$ represent a plane that includes this normal vector. The intersection
of the plane and the surface makes a curve $\textbf{x}(l)$, which has
curvature $\kappa_{0}$ at the point $\textbf{x}_{0}$, as calculated from
(3.14). The two principal curvatures $\kappa_{1}$ and $\kappa_{2}$ represent
the maximum and minimum curvatures, $\kappa_{0}$, from all possible planes
$A_{n}$. We choose the convention for the principal curvatures such that
convex curves have positive curvature and concave curves have negative
curvatures. Therefore, small values of the second principal curvature,
$\kappa_{2}$, represent regions on the surface where the concavity is
increased.

The second principal curvature for the W7-X plasma surface is shown in Figure
3.11(b). Although $\kappa_{2}$ and the shape gradients are evaluated on
different surfaces, we note that regions of high concavity (negative
$\kappa_{2}$) coincide with regions of large, positive $\mathcal{G}$ (Figure
3.9). The regions of high concavity also correspond to the regions where the
optimization procedure tends to place the winding surface closest to the
plasma (Figure 3.11). We recognize that our winding surface optimization
accounts for several engineering considerations in addition to reproducing the
desired plasma surface. However, for a wide range of parameters the winding
surfaces we obtain feature small $d_{\text{coil-plasma}}$ in the bean-shaped
cross-sections (Figures 3.2 and 3.3). Thus $\kappa_{2}$, which is exceedingly
fast to compute, may serve as a target for optimization of the plasma
configuration. By minimizing the regions of high concavity, it may be possible
to find stellarator equilibria that are more amenable to coils that are
positioned farther from the plasma. Any increase in the minimal distance
between the plasma and the coils has implications for the size of a reactor,
where $d_{\text{coil-plasma}}$ is set by the required blanket width. Similar
metrics are considered in the ROSE code, such as the integrated absolute value
of the Gaussian curvature and integrated absolute value of the maximum
curvature [59].

(a)

(b)

Figure 3.11: (a) The minimum distance between the W7-X plasma surface and the
optimized winding surface obtained in Section 3.5.2 and (b) the second
principle curvature $\kappa_{2}$ are shown as a function of location on the
plasma surface. Locations of large negative $\kappa_{2}$ coincide with regions
where the optimization resulted in small $d_{\text{coil-plasma}}$. Figure
adapted from [185] with permission.

###  3.8 Conclusions

We have outlined a new method for the optimization of the stellarator coil-
winding surface using a continuous current potential approach. Rather than
evolving filamentary coil shapes, we use REGCOIL to obtain the current density
on a winding surface and optimize the winding surface using analytic gradients
of the objective function. We have shown that we can indirectly improve the
coil curvature and toroidal extent by targeting the root-mean-squared current
density in our objective function (Figure 3.1). This approach offers several
potential advantages over other nonlinear coil optimization tools.

  1. 1.

The difficulty of the optimization is reduced by the application of the
REGCOIL method, which takes the form of a linear least-squares system. The
optimal coil shapes on a given winding surface can thus be efficiently and
robustly computed.

  2. 2.

By fixing the maximum current density to obtain the regularization in REGCOIL,
we eliminate the need to implement an additional equality constraint or
arbitrary weight in the objective function.

  3. 3.

By using REGCOIL to compute coil shapes on a given surface, we can apply the
adjoint method for computing derivatives (Section 3.4). This allows us to
reduce the number of function evaluations required during the nonlinear
optimization by a factor of $\approx 50$.

  4. 4.

Given the critical role coil design plays in the stellarator optimization
process, it is important to have many tools that approach the problem from
different angles. Our approach differs from the other available nonlinear coil
optimization applications [154, 216, 218, 32, 243] as we optimize a continuous
current potential.

We have demonstrated this method by optimizing coils for W7-X and HSX
(Sections 3.5.2 and 3.5.3). We find that we can simultaneously decrease the
integral-squared error in reproducing the plasma surface, increase the volume
contained within the winding surface, maintain the minimum coil-plasma
distance, and improve the coil metrics over REGCOIL solutions computed on the
initial winding surfaces (Tables 3.2 and 3.3). Several features of these
optimized winding surfaces are noteworthy. While the coil-plasma distance must
be small in concave regions, it can increase greatly on the outboard, convex
side of the bean cross-section. At triangle-shaped cross-sections, the winding
surface obtains a somewhat “pinched” appearance (Figures 3.3, 3.4, and 3.6). A
similar W7-X winding surface shape has been obtained with the ONSET code (see
ref. [154], Figure 5). Further work is required to understand this behavior.

There are several limitations to this approach that should be noted. First, we
have applied a local nonlinear optimization algorithm. This is a reasonable
choice if the initial condition is close to a global optimum. Second, we
currently have not added coil-specific metrics to our objective function (for
example, curvature or length). This could be implemented if necessary for
engineering purposes.

We should also note that this application does not allow for the full benefits
of adjoint methods. While adjoint methods significantly reduce CPU time if the
solve is the computational bottleneck, this is not the case for the REGCOIL
system. Other applications that are dominated by the linear solve CPU time
would see increased benefits from the implementation of an adjoint method, as
will be seen in the following Chapters.

We demonstrate a technique for visualization of shape derivatives in real
space rather than Fourier space. This shape gradient describes how an
objective function changes with respect to normal displacements of the winding
surface. We apply this technique to visualize the derivatives of the integral-
squared normal field on the plasma surface and the root-mean-squared current
density for the W7-X plasma surface and two winding surfaces (Figure 3.9).
This diagnostic identifies the concave regions as being very sensitive to the
positions of coils, as has been observed from previous coil optimization
efforts. We will continue to gain insight from the shape gradient concept in
Chapters 4, 5, and 6.

## Chapter 4: Adjoint-based optimization of neoclassical properties

Several critical quantities for stellarator design arise from neoclassical
physics, the kinetic theory of collisional transport in the presence of
magnetic field gradients and curvature. This so-called neoclassical transport
results from the random-walk of charged particles as they exhibit guiding
center motion. Due to the complicated guiding center orbits present in a 3D
field, neoclassical transport is generally enhanced in a stellarator. One of
the primary goals of stellarator optimization is to reduce this transport.
Furthermore, the bootstrap current, driven by collisional processes, should be
minimized in low-shear designs or if an island divertor system is to be used.
These neoclassical properties are described by solutions of the drift-kinetic
equation (DKE),

$\displaystyle\left(v_{||}\hat{\textbf{b}}+\textbf{v}_{\text{d}}\right)\cdot\nabla
f=C(f),$ (4.1)

where $f$ is the distribution function,
$v_{||}=\textbf{v}\cdot\hat{\textbf{b}}$ is the parallel component of the
velocity, $\textbf{v}_{\text{d}}$ is the guiding center drift velocity, and
$C$ is the collision operator. The DKE is obtained from the Fokker-Planck
equation under the assumption that the plasma is strongly magnetized such that
(4.1) describes length scales much longer than the gyroradius and frequencies
much smaller than the gyrofrequency. We have taken the equilibrium limit,
assuming time scales longer than the gyroperiod but shorter than the transport
time scale on which the profiles relax. In this Chapter we make an additional
assumption of local thermodynamic equilibrium, such that $f\approx f_{M}$, a
Maxwellian distribution (defined in Section 4.2), to lowest order. This
assumption is valid in stellarator configurations, provided that the
collisionless orbits are sufficiently confined and the collision frequency is
not too low [227, 33]. The departure from a Maxwellian, $f_{1}$, is driven by
gradients in $f_{M}$ due to variations in the density, temperature, and
electrostatic potential. The drift-kinetic equation is described in many
references, including Chapter 7 in [99] and [94, 97].

In this Chapter, we will apply both the discrete and continuous adjoint
methods described in Chapter 2 to efficiently compute derivatives of functions
that depend on such solutions of the drift kinetic equation. This analysis
will allow us to efficiently optimize the local magnetic field for several
neoclassical quantities in addition to analyzing their sensitivity to changes
in the magnetic field.

The material in this Chapter has been adapted from [186].

###  4.1 Introduction

Neoclassical transport is governed by solutions of the drift kinetic equation
(DKE) (5.131) from which moments (e.g., radial fluxes and bootstrap current)
are computed. The DKE local to a flux surface can be solved numerically [140,
18]. However, this four-dimensional problem is expensive to solve within an
optimization loop, especially in low-collisionality regimes for which
increased pitch-angle resolution is required to resolve the collisional
boundary layer.

Therefore, it is sometimes desirable to consider an analytic reduction of the
DKE. Under the assumption of low collisionality, a bounce-averaged DKE can be
considered [17, 34]. While bounce-averaging can significantly reduce the
computational cost by decreasing the spatial dimensionality, this approach
typically requires restrictions on the geometry, such as closeness to
omnigeneity or a model magnetic field. Additional reduction of the DKE can be
made in low-collisionality regimes, resulting in semi-analytic expressions.
For example the effective ripple, $\epsilon_{\text{eff}}$ [168], quantifies
the geometric dependence of the $1/\nu$ radial transport ($\nu$ is the
collision frequency) and has been widely used during optimization studies
[242, 134, 106]. (The effective ripple will be discussed further in Chapter 5
and Appendix M.) The $1/\nu$ regime, though, is only relevant when $E_{r}$ is
small enough that the typical poloidal rotation frequency is much smaller than
the typical collision frequency [116], which is not always an experimentally-
relevant regime. A low-collisionality semi-analytic bootstrap current model
[205] is also commonly adopted for stellarator design [15, 114]. However, this
analytic expression is known to be ill-behaved near rational surfaces.
Furthermore, benchmarks with numerical solutions of the DKE in the low-
collisionality limit have been shown to differ significantly from the semi-
analytic model [16, 127]. Any analytic reduction of the DKE implies additional
assumptions, such as on the collisionality, size of $E_{r}$, or on the
magnetic geometry.

Due to the limitations of bounce-averaged and semi-analytic models, there are
benefits to computing neoclassical quantities using numerical solutions to the
DKE without approximation. With the numerical methods currently used for
stellarator optimization, this approach becomes computationally challenging
within an optimization loop. Due to their fully three-dimensional nature,
optimization of stellarator geometry requires navigation through high-
dimensional spaces, such as the space of the shape of the outer boundary of
the plasma or the shapes of electromagnetic coils. The number of parameters
required to describe these spaces, $N$, is often quite large
($\mathcal{O}(10^{2})$). Knowledge of the gradient of the objective function
with respect to these parameters can significantly improve the convergence to
a local minimum. Once a descent direction is identified, each iteration
reduces to a one-dimensional line search. Gradient-based optimization with the
Levenberg-Marquardt algorithm in the STELLOPT code [218] has been widely used
in the stellarator community and led to the design of NCSX [197].

Although derivative information is valuable, numerically computing the
derivative of a figure of merit $f$ (for example, with finite-difference
derivatives) can be prohibitively expensive, as $f$ must be evaluated
$\mathcal{O}(N)$ times. For neoclassical optimization, this implies solving
the DKE $\mathcal{O}(N)$ times; thus including finite-collisionality
neoclassical quantities in the objective function is often impractical. In
this Chapter, we describe an adjoint method for neoclassical optimization.
With this method, the computation of the derivatives of $f$ with respect to
$N$ parameters has cost comparable to solving the DKE twice, thus making the
inclusion of these quantities possible within an optimization loop. In this
Chapter, we obtain derivatives of neoclassical figures of merit with respect
to local geometric parameters on a surface rather than the outer boundary or
coil shapes. However, the geometric derivatives we compute provide an
important step toward adjoint-based optimization of MHD equilibria, as
discussed in Section 4.5.2 and Chapter 5.

In Section 4.2, we provide an overview of the numerical solution of the DKE
local to a flux surface. In Section 4.3 the adjoint neoclassical method is
described. The continuous and discrete approaches for this problem are
presented, and their implementation and benchmarks are discussed in Section
4.4. The adjoint method is used to compute derivatives of moments of the
neoclassical distribution function with respect to local geometric quantities.
The derivative information can be used to identify regions of increased
sensitivity to magnetic perturbations, as discussed in Section 4.5.1. We
demonstrate adjoint-based optimization in Section 4.5.2 by locally modifying
the field strength on a flux surface. A discussion of the application of this
method for optimization of MHD equilibria is presented in 4.5.2. Finally, the
adjoint method is applied to accelerate the calculation of the ambipolar
electric field in Section 4.5.3.

###  4.2 Drift kinetic equation

The local drift kinetic equation is,

$\left(v_{||}\hat{\textbf{b}}+\textbf{v}_{E}\right)\cdot\nabla
f_{1s}-C_{s}(f_{1s})=-\textbf{v}_{\text{m}s}\cdot\nabla\psi\frac{\partial
f_{Ms}}{\partial\psi},$ (4.2)

Here $\hat{\textbf{b}}=\textbf{B}/B$ is a unit vector in the direction of the
magnetic field, $v_{||}=\textbf{v}\cdot\hat{\textbf{b}}$ is the parallel
component of the velocity, and $2\pi\psi$ is the toroidal flux. The Fokker-
Planck collision operator is $C_{s}(f_{1s})$, linearized about a Maxwellian
$f_{Ms}=n_{s}v_{ts}^{-3}\pi^{-3/2}e^{-v^{2}/v_{ts}^{2}}$ where
$v_{ts}=\sqrt{2T_{s}/m_{s}}$ is the thermal speed, $n_{s}$ is the density,
$T_{s}$ is the temperature, $m_{s}$ is the mass, and the subscript indicates
species. In (4.2), derivatives are performed holding
$W_{s}=m_{s}v^{2}/2+q_{s}\Phi$ and $\mu=v_{\perp}^{2}/2B$ fixed, where
$v=\sqrt{\textbf{v}\cdot\textbf{v}}$ is the magnitude of velocity, $\Phi$ is
the electrostatic potential, $v_{\perp}=\sqrt{v^{2}-v_{||}^{2}}$ is the
perpendicular velocity, and $q_{s}$ is the charge. The radial magnetic drift
is,

$\textbf{v}_{\text{m}s}\cdot\nabla\psi=\frac{m_{s}}{q_{s}B^{2}}\left(v_{||}^{2}+\frac{v_{\perp}^{2}}{2}\right)\hat{\textbf{b}}\times\nabla
B\cdot\nabla\psi,$ (4.3)

assuming a magnetic field in MHD force balance, and $\textbf{v}_{E}$ is the
$\textbf{E}\times\textbf{B}$ velocity,

$\textbf{v}_{E}=\frac{\textbf{B}\times\nabla\Phi}{B^{2}}.$ (4.4)

Throughout we assume $\Phi=\Phi(\psi)$ such that (4.2) is linear. In (4.2) we
will not consider the effect of inductive electric fields, as these can be
assumed to be small for stellarators without inductive current drive. We also
do not consider the effects of magnetic drifts tangential to the flux surface
in (4.2), as these only become important when $E_{r}$ is small [184]. We can
assume radial locality, manifested by the absence of any radial derivatives of
$f_{1s}$ in (4.2), when $\nu_{*}\gg\rho_{*}$ [33], where
$\nu_{*}=\nu/(v_{t}/L)\ll 1$ is the normalized collision frequency for
macroscopic scale length $L$ and $\rho_{*}=v_{t}m/(LqB)$ is the normalized
gyrofrequency. Numerical solutions to (4.2) are computed with the Stellarator
Fokker-Planck Iterative Neoclassical Solver (SFINCS) [140] code which allows
for general stellarator geometry with flux surfaces.

SFINCS solves (4.2) locally on a flux surface $\psi$, a four-dimensional
system. The SFINCS coordinates include two angles (poloidal angle $\theta$ and
toroidal angle $\phi$), speed $X_{s}=v/v_{ts}$, and pitch angle
$\xi_{s}=v_{||}/v$. Specifics about the implementation of (4.2) in the SFINCS
code are described in Appendix D. We will refer to two choices of
implementation: the full trajectory model and the DKES trajectory model. The
full trajectory model maintains $\mu$ conservation as radial coupling (terms
involving $\partial f_{1s}/\partial\psi$) is dropped. While the DKES model
does not conserve $\mu$ when $E_{r}\neq 0$, the adjoint operator under the
DKES model takes a particularly simple form, as discussed in Section 4.3.1.
This model also does not introduce any unphysical constraints on the
distribution function when $E_{r}=0$, as occurs for the full trajectory model
[140]. These constraints motivate the introduction of particle and heat
sources, which are discussed in the following Section. We will discuss details
of the implementation of the DKE in the SFINCS code, as these need to be
considered in arriving at the adjoint equation. However, the adjoint
neoclassical approach is quite general and could be implemented in other
drift-kinetic codes with slight modification.

From solutions of (4.2), several neoclassical quantities are computed,
including the flux-surface averaged parallel flow,

$\displaystyle V_{||,s}=\frac{\left\langle B\int
d^{3}v\,f_{1s}v_{||}\right\rangle_{\psi}}{n_{s}\langle
B^{2}\rangle_{\psi}^{1/2}},$ (4.5)

the radial particle flux,

$\displaystyle\Gamma_{s}=\left\langle\int
d^{3}v\,\left(\textbf{v}_{\text{m}s}\cdot\nabla\rho\right)f_{1s}\right\rangle_{\psi},$
(4.6)

and the radial heat flux (sometimes referred to as an energy flux),

$\displaystyle Q_{s}=\left\langle\int
d^{3}v\,\frac{m_{s}v^{2}}{2}\left(\textbf{v}_{\text{m}s}\cdot\nabla\rho\right)f_{1s}\right\rangle_{\psi}.$
(4.7)

Here the flux-surface average of a quantity $A$ is,

$\displaystyle\langle A\rangle_{\psi}$
$\displaystyle=\frac{\int_{0}^{2\pi}d\theta\int_{0}^{2\pi}d\phi\,\sqrt{g}A}{V^{\prime}(\psi)}$
(4.8a) $\displaystyle V^{\prime}(\psi)$
$\displaystyle=\int_{0}^{2\pi}d\theta\int_{0}^{2\pi}d\phi\,\sqrt{g},$ (4.8b)

and $\sqrt{g}=\left(\nabla\psi\times\nabla\theta\cdot\nabla\phi\right)^{-1}$
is the Jacobian. We will also consider species-summed quantities including the
bootstrap current, $J_{b}=\sum_{s}q_{s}n_{s}V_{||,s}$, the radial current,
$J_{r}=\sum_{s}q_{s}\Gamma_{s}$, and the total heat flux,
$Q_{\text{tot}}=\sum_{s}Q_{s}$. Here the effective normalized radius is
$\rho=\sqrt{\psi/\psi_{0}}$, where $2\pi\psi_{0}$ is the toroidal flux at the
boundary.

####  4.2.1 Sources and constraints

To avoid unphysical constraints on $f_{1s}$ implied by the moment equations of
(4.2) in the presence of a non-zero $E_{r}$ [140], particle and heat sources
are added to the DKE (D.1),

$\displaystyle\mathbb{L}_{0s}f_{1s}-C_{s}(f_{1s})-f_{Ms}\left(X_{s}^{2}-\frac{5}{2}\right)S_{1s}^{f}(\psi)-f_{Ms}\left(X_{s}^{2}-\frac{3}{2}\right)S_{2s}^{f}(\psi)=\mathbb{S}_{0s},$
(4.9)

where $S_{1s}^{f}(\psi)$ and $S_{2s}^{f}(\psi)$ are unknowns such that
$S_{1s}^{f}$ provides a particle source and $S_{2s}^{f}$ provides a heat
source. The collisionless trajectory operator in SFINCS coordinates is,

$\displaystyle\mathbb{L}_{0s}=\dot{\textbf{x}}\cdot\nabla+\dot{X}_{s}\frac{\partial}{\partial
X_{s}}+\dot{\xi}_{s}\frac{\partial}{\partial\xi_{s}},$ (4.10)

and the inhomogeneous drive term is
$\mathbb{S}_{0s}=-(\textbf{v}_{\text{m}s}\cdot\nabla\psi)\partial
f_{Ms}/\partial\psi$. The source functions are determined via the requirement
that $\langle\int d^{3}v\,f_{1s}\rangle_{\psi}=0$ and $\langle\int
d^{3}v\,X_{s}^{2}f_{1s}\rangle_{\psi}=0$ (i.e. $f_{1s}$ does not provide net
density or pressure). So, the following system of equations is solved,

$\displaystyle\underbrace{\left[\begin{array}[]{ccc}\mathbb{L}_{0s}-C_{s}&-f_{Ms}(X_{s}^{2}-\frac{5}{2})&-f_{Ms}(X_{s}^{2}-\frac{3}{2})\\\
\mathbb{L}_{1s}&0&0\\\
\mathbb{L}_{2s}&0&0\end{array}\right]}_{\mathbb{L}_{s}}\underbrace{\left[\begin{array}[]{c}f_{1s}\\\
S_{1s}^{f}\\\
S_{2s}^{f}\end{array}\right]}_{F_{s}}=\underbrace{\left[\begin{array}[]{c}\mathbb{S}_{0s}\\\
0\\\ 0\end{array}\right]}_{\mathbb{S}_{s}}.$ (4.20)

The velocity-space averaging operations are denoted
$\mathbb{L}_{1s}f_{1s}=\langle\int d^{3}v\,f_{1s}\rangle_{\psi}$ and
$\mathbb{L}_{2s}f_{1s}=\langle\int d^{3}v\,f_{1s}X_{s}^{2}\rangle_{\psi}$. The
full multi-species system can be written as,

$\displaystyle\left[\begin{array}[]{c}\mathbb{L}_{1}\\\ \vdots\\\
\mathbb{L}_{N_{\text{species}}}\end{array}\right]\left[\begin{array}[]{c}F_{1}\\\
\vdots\\\
F_{N_{\text{species}}}\end{array}\right]=\left[\begin{array}[]{c}\mathbb{S}_{1}\\\
\vdots\\\ \mathbb{S}_{N_{\text{species}}}\end{array}\right].$ (4.30)

Here the linear systems corresponding to each species as in (4.20) are coupled
through the collision operator. We use the following notation to refer to the
above system,

$\displaystyle\mathbb{L}F=\mathbb{S}.$ (4.31)

###  4.3 Adjoint approach

The goal of the adjoint neoclassical approach is to compute derivatives of a
moment of the distribution function efficiently, $\mathcal{R}$ (e.g.,
$V_{||,s},\Gamma_{s},Q_{s},J_{b},J_{r},Q_{\text{tot}})$, with respect to many
parameters. Consider a set of parameters,
$\Omega=\\{\Omega_{i}\\}_{i=1}^{N_{\Omega}}$, on which $\mathcal{R}$ depends.
Computing a forward-difference derivative with respect to $\Omega$ requires
$N_{\Omega}+1$ solutions of (4.31). With the adjoint approach,
$\partial\mathcal{R}/\partial\Omega$ can be computed with one solution of
(4.31) and one solution of a linear adjoint equation of the same size as
(4.31). Thus if $N_{\Omega}$ is very large and the solution to (4.31) is
computationally expensive to obtain, the adjoint approach can reduce the cost
by $N_{\Omega}$. For stellarator optimization, it is desirable to compute
derivatives with respect to parameters that describe the magnetic geometry. In
fully three-dimensional geometry, $N_{\Omega}$ is $\mathcal{O}(10^{2})$ and
solving (4.31) is the most expensive part of computing $\mathcal{R}$ (rather
than constructing the linear system or taking a moment of the distribution
function). The discretized linear system is typically very large ($N\sim
10^{5}-10^{6}$ for the calculations shown in the Chapter) and sparse. Thus
matrix-matrix products are significantly less expensive than the linear solve,
which is performed with a preconditioned Krylov iterative method.
Consequently, the adjoint method provides a factor of $N_{\Omega}\sim 10^{2}$
savings over both the forward sensitivity and finite-difference methods, as
described in Section 2.2.1. The adjoint method also allows us to avoid
additional round-off or truncation error arising from finite-difference
derivatives. In what follows, we consider $\Omega$ to be a set of parameters
describing the magnetic geometry, which will be specified in Section 4.4.

We compute the derivatives of $\mathcal{R}$ using two approaches. In the first
approach, we define an inner product that involves integrals over the
distribution function, and an adjoint operator is obtained with respect to
this inner product. This is the continuous approach introduced in Section
2.2.2. In the second approach, we consider the DKE after discretization,
defining an adjoint operator with respect to the Euclidean dot product. This
is the discrete approach introduced in Section 2.2.1. While these approaches
should provide identical results within discretization error, the advantages
and drawbacks of each method will be discussed at the end of Section 4.3.2.

####  4.3.1 Continuous approach

Let $F=\\{F_{s}\\}_{s=1}^{N_{\text{species}}}$ be the set of unknowns computed
with SFINCS before discretization, denoted by the column vector in (4.30) with
$F_{s}$ given by (4.20). That is, $F$ consists of a set of
$N_{\text{species}}$ distribution functions over $(\theta,\phi,X_{s},\xi_{s})$
and their associated source functions. We define an inner product between two
such quantities in the following way,

$\displaystyle\langle F,G\rangle=\sum_{s}\left\langle\int
d^{3}v\,\frac{f_{1s}g_{1s}}{f_{Ms}}\right\rangle_{\psi}+S_{1s}^{f}S_{1s}^{g}+S_{2s}^{f}S_{2s}^{g}.$
(4.32)

Here the superscript on $S_{1s}$ and $S_{2s}$ denotes the distribution
function with which the source functions are associated and the sum is over
species. The space of continuous functions, $F$, of this form such that
$\langle F,F\rangle$ is bounded will be denoted by $\mathcal{H}$. It can be
seen that (4.32) is indeed an inner product, as it satisfies conjugate
symmetry ($\langle G,F\rangle=\langle F,G\rangle$ $\forall
F,G\in\mathcal{H}$), linearity ($\langle F+G,H\rangle=\langle
F,H\rangle+\langle G,H\rangle$ $\forall F,G,H\in\mathcal{H}$ and $\langle
F,aG\rangle=a\langle F,G\rangle$ $\forall F,G\in\mathcal{H}$,
$a\in\mathbb{R}$), and positive definiteness ($\langle F,F\rangle\geq 0$ and
$\langle F,F\rangle=0$ only if $F=0$ $\forall F\in\mathcal{H}$) [199]. This
implies that if $\mathcal{H}$ is finite-dimensional, then for any linear
operator $L$ there exists a unique adjoint operator $L^{\dagger}$ such that
$\langle LF,G\rangle=\langle F,L^{\dagger}G\rangle$ for all
$F,G\in\mathcal{H}$. While here $\mathcal{H}$ is not finite-dimensional, we
will show that such an adjoint operator exists for this inner product.

Note that the norm associated with this inner product $||F||=\sqrt{\langle
F,F\rangle}$ is similar to the free energy norm,

$\displaystyle W=\sum_{s}\left\langle\int
d^{3}v\,\frac{T_{s}f_{1s}^{2}}{2f_{Ms}}\right\rangle_{\psi},$ (4.33)

which obeys a conservation equation in gyrokinetic theory [132, 2, 141]. The
choice of inner product (4.32) is advantageous, as the linearized Fokker-
Planck collision operator becomes self-adjoint for species linearized about
Maxwellians with the same temperature. In what follows, we assume that all
included species are of the same temperature. This assumption could be lifted,
with a modification to the collision operator that appears in the adjoint
equation (Appendix E). This assumption is not necessary when using the
discrete approach (Section 4.3.2).

Consider a moment of the distribution function
$\mathcal{R}\in\\{V_{||,s},\Gamma_{s},Q_{s},J_{b},J_{r},Q_{\text{tot}}\\}$,
which can be written as an inner product with a vector
$\widetilde{\mathcal{R}}\in\mathcal{H}$,

$\displaystyle\mathcal{R}=\langle F,\widetilde{\mathcal{R}}\rangle,$ (4.34)

according to (4.32). For example,

$\displaystyle\widetilde{J_{r}}=\left[\begin{array}[]{c}q_{s}\textbf{v}_{\text{m}s}\cdot\nabla\psi
f_{Ms}\\\ 0\\\ 0\end{array}\right]_{s=1}^{N_{\text{species}}},$ (4.38)

where the column structure corresponds with that in (4.20) and (4.30).

We are interested in computing the derivative of $\mathcal{R}$ with respect to
a set of parameters, $\Omega=\\{\Omega_{i}\\}_{i=1}^{N_{\Omega}}$ such that
the DKE is satisfied. Computing such a derivative with the forward sensitivity
method requires that we compute $\partial F(\Omega)/\partial\Omega_{i}$ from
the linearized DKE,

$\displaystyle\frac{\partial\mathbb{L}(\Omega)}{\partial\Omega_{i}}F+\mathbb{L}\frac{\partial
F(\Omega)}{\partial\Omega_{i}}=\frac{\partial\mathbb{S}(\Omega)}{\partial\Omega_{i}},$
(4.39)

for each $\Omega_{i}$ and evaluate the derivative using the chain rule,

$\displaystyle\frac{\partial\mathcal{R}(\Omega,F(\Omega))}{\partial\Omega_{i}}=\frac{\partial\mathcal{R}(\Omega,F)}{\partial\Omega_{i}}+\left\langle\widetilde{\mathcal{R}},\frac{\partial
F(\Omega)}{\partial\Omega_{i}}\right\rangle.$ (4.40)

We see that the forward sensitivity method requires solutions of $N_{\Omega}$
linear systems of the same dimension as the DKE (4.31).

To avoid this additional computational cost, we instead apply the adjoint
method by constructing the Lagrangian functional, enforcing (4.31) as a
constraint,

$\displaystyle\mathcal{L}(\Omega,F,\lambda^{\mathcal{R}})=\mathcal{R}(\Omega,F)+\left\langle\lambda^{\mathcal{R}},\mathbb{L}F-\mathbb{S}\right\rangle.$
(4.41)

Here $\lambda^{\mathcal{R}}$ is the Lagrange multiplier. We obtain the adjoint
equation by finding a stationary point of $\mathcal{L}$ with respect to $F$,

$\displaystyle\delta\mathcal{L}(\Omega,F,\lambda^{\mathcal{R}};\delta
F)=\langle\delta
F,\widetilde{\mathcal{R}}\rangle+\left\langle\lambda^{\mathcal{R}},\mathbb{L}\delta
F\right\rangle=0.$ (4.42)

We can now use the adjoint property to express the above as,

$\displaystyle\delta\mathcal{L}(\Omega,F,\lambda^{\mathcal{R}};\delta
F)=\langle\delta
F,\widetilde{\mathcal{R}}+\mathbb{L}^{\dagger}\lambda^{\mathcal{R}}\rangle.$
(4.43)

A stationary point of $\mathcal{L}$ with respect to $F$ corresponds to
$\lambda^{\mathcal{R}}$ which satisfies the weak form of the adjoint equation,

$\displaystyle\mathbb{L}^{\dagger}\lambda^{\mathcal{R}}+\widetilde{\mathcal{R}}=0.$
(4.44)

With this adjoint variable, we can now compute derivatives of $\mathcal{R}$
with respect to any parameter by computing the corresponding perturbations of
$\mathcal{L}$,

$\displaystyle\frac{\partial\mathcal{R}(\Omega,F(\Omega))}{\partial\Omega_{i}}=\frac{\partial\mathcal{L}(\Omega,F,\lambda^{\mathcal{R}})}{\partial\Omega_{i}}=\frac{\partial\mathcal{R}(\Omega,F)}{\partial\Omega_{i}}+\left\langle\lambda^{\mathcal{R}},\frac{\partial\mathbb{L}(\Omega)}{\partial\Omega_{i}}F-\frac{\partial\mathbb{S}(\Omega)}{\partial\Omega_{i}}\right\rangle.$
(4.45)

The first term on the right hand side accounts for the explicit dependence on
$\Omega_{i}$ while the second accounts for the implicit dependence on
$\Omega_{i}$ through $F$. Thus, using (4.45), the derivative with respect to
$\Omega$ can be computed with the solution to two linear systems, (4.31) and
(4.44). The partial derivatives on the right hand side of (4.45) can be
computed analytically by considering the explicit geometric dependence of
$\mathcal{R}$, $\mathbb{L}$, and $\mathbb{S}$.

When $N_{\Omega}$ is large, the cost of computing
$\partial\mathcal{R}/\partial\Omega$ using (4.45) is dominated not by the
linear solve but by constructing $\partial\mathbb{S}/\partial\Omega$ and
$\partial\mathbb{L}/\partial\Omega$ and computing the inner product. Thus the
cost still scales with $N_{\Omega}$. However, we obtain a significant savings
in comparison with forward-difference derivatives, as shown in Section 4.4.

The adjoint operator for each species takes the following form,

$\displaystyle\mathbb{L}_{s}^{\dagger}=\left[\begin{array}[]{c c
c}\mathbb{L}_{0s}^{\dagger}-C_{s}&f_{Ms}&f_{Ms}X_{s}^{2}\\\
\mathbb{L}_{1s}^{\dagger}&0&0\\\
\mathbb{L}_{2s}^{\dagger}&0&0\end{array}\right],$ (4.49)

where $\mathbb{L}_{1s}^{\dagger}=5/2\mathbb{L}_{1s}-\mathbb{L}_{2s}$ and
$\mathbb{L}_{2s}^{\dagger}=3/2\mathbb{L}_{1s}-\mathbb{L}_{2s}$. The same
column structure is used as for the forward operator (4.30),
$\mathbb{L}^{\dagger}=\\{\mathbb{L}_{s}^{\dagger}\\}_{i=1}^{N_{\text{species}}}$.
The quantity $\mathbb{L}_{0s}^{\dagger}$ satisfies $\langle\int
d^{3}v\,g_{1s}\mathbb{L}_{0s}f_{1s}/f_{Ms}\rangle_{\psi}=\langle\int
d^{3}v\,f_{1s}\mathbb{L}_{0s}^{\dagger}g_{1s}/f_{Ms}\rangle_{\psi}$ and
depends on which trajectory model is applied. The expression (4.49) can be
verified by noting that

$\displaystyle\langle\mathbb{L}F,G\rangle$
$\displaystyle=\sum_{s}\left\langle\frac{f_{1s}\left((\mathbb{L}^{\dagger}_{0s}-C_{s})g_{1s}+f_{Ms}\left(S_{1s}^{g}+S_{2s}^{g}X_{s}^{2}\right)\right)}{f_{Ms}}\right\rangle_{\psi}+S_{1s}^{f}\mathbb{L}_{1s}^{\dagger}g_{1s}+S_{2s}^{f}\mathbb{L}_{2s}^{\dagger}g_{1s}$
$\displaystyle=\langle F,\mathbb{L}^{\dagger}G\rangle.$ (4.50)

For the DKES trajectories the adjoint operator is,

$\displaystyle\mathbb{L}_{0s}^{\dagger}=-\mathbb{L}_{0s}.$ (4.51)

This anti-self-adjoint property is used in obtaining the variational principle
which provides bounds on neoclassical transport coefficients in the DKES code
[230]. For full trajectories it is,

$\displaystyle\mathbb{L}_{0s}^{\dagger}=-\mathbb{L}_{0s}+\frac{q_{s}}{T_{s}}\Phi^{\prime}(\psi)\textbf{v}_{\text{m}s}\cdot\nabla\psi.$
(4.52)

The anti-self-adjoint property does not hold for this trajectory model as the
$\textbf{E}\times\textbf{B}$ drift (F.9) is no longer divergenceless. Appendix
F contains details on obtaining these adjoint operators.

####  4.3.2 Discrete approach

Next, we consider the discrete adjoint approach. Let
$\overrightarrow{\textbf{F}}$ be the set of unknowns computed with SFINCS
after discretization of $F$. The linear DKE (4.31) upon discretization can
then be written schematically as,

$\displaystyle\overleftrightarrow{\textbf{L}}\overrightarrow{\textbf{F}}=\overrightarrow{\textbf{S}}.$
(4.53)

In this case, we can define an inner product as the vector dot product,

$\displaystyle\langle\overrightarrow{\textbf{F}},\overrightarrow{\textbf{G}}\rangle=\overrightarrow{\textbf{F}}\cdot\overrightarrow{\textbf{G}}.$
(4.54)

In real Euclidean space, the adjoint operator,
$\left(\overleftrightarrow{\textbf{L}}\right)^{\dagger}$, which satisfies,

$\displaystyle\left\langle\overleftrightarrow{\textbf{L}}\overrightarrow{\textbf{F}},\overrightarrow{\textbf{G}}\right\rangle=\left\langle\overrightarrow{\textbf{F}},\left(\overleftrightarrow{\textbf{L}}\right)^{\dagger}\overrightarrow{\textbf{G}}\right\rangle$
(4.55)

is simply the transpose of the matrix,
$\left(\overleftrightarrow{\textbf{L}}\right)^{T}$. Again, the moments of the
distribution function, $\mathcal{R}$ can be expressed as an inner product with
a vector $\overrightarrow{\textbf{R}}$,

$\displaystyle\mathcal{R}=\langle\overrightarrow{\textbf{F}},\overrightarrow{\textbf{R}}\rangle.$
(4.56)

Using the discrete approach, the following adjoint equation must be solved

$\displaystyle\left(\overleftrightarrow{\textbf{L}}\right)^{T}\overrightarrow{\bm{\lambda}}^{\mathcal{R}}=\overrightarrow{\textbf{R}}.$
(4.57)

The adjoint variable, $\overrightarrow{\bm{\lambda}}^{\mathcal{R}}$, can again
be used to compute the derivative of $\mathcal{R}$ with respect to $\Omega$,

$\displaystyle\frac{\partial\mathcal{R}\left(\Omega,\overrightarrow{\textbf{F}}(\Omega)\right)}{\partial\Omega_{i}}=\frac{\partial\mathcal{R}\left(\Omega,\overrightarrow{\textbf{F}}\right)}{\partial\Omega_{i}}+\left\langle\overrightarrow{\bm{\lambda}}^{\mathcal{R}},\left(\frac{\partial\overrightarrow{\textbf{S}}(\Omega)}{\partial\Omega_{i}}-\frac{\partial\overleftrightarrow{\textbf{L}}(\Omega)}{\partial\Omega_{i}}\overrightarrow{\textbf{F}}\right)\right\rangle.$
(4.58)

As with the continuous approach, the partial derivatives on the right hand
side can be computed analytically. In this way, the derivative of
$\mathcal{R}$ with respect to $\Omega$ can be computed with only two linear
solves, (4.53) and (4.57).

In the SFINCS implementation, the DKE is typically solved with the
preconditioned GMRES algorithm. In the continuous approach, a preconditioner
matrix for both the forward and adjoint operator must be $LU$-factorized. Here
the preconditioner matrix is the same as the full matrix but without cross-
species or speed coupling. As the adjoint matrix is sufficiently different
from the forward matrix, we do not obtain convergence when the same
preconditioner is used for both problems. However, in the discrete approach,
the $LU$-factorization for the preconditioner of the forward matrix can be
reused for the preconditioner of the adjoint matrix. (If a matrix $A$ has been
factorized as $A=LU$ then $A^{T}=U^{T}L^{T}$ where $U^{T}$ is lower triangular
and $L^{T}$ is upper triangular). This provides a significant reduction in
memory and computational cost for the discrete approach.

Furthermore, the discrete adjoint approach provides the exact derivatives for
the discretized problem. With this method, the adjoint equation is obtained
using the vector dot product and matrix transpose, which can be computed
without any numerical approximation. The error in the derivatives obtained by
the adjoint method is therefore only limited by the tolerance to which the
linear solve is performed with GMRES. On the other hand, the continuous
adjoint approach relies on a continuous inner product that must ultimately be
approximated numerically. Thus the continuous approach provides the exact
derivatives only in the limit that the discrete approximation of the inner
product exactly reproduces the continuous inner product. Therefore we expect
the results of the discrete and adjoint approaches to agree within
discretization error, as will be demonstrated in Section 4.4.

The continuous approach can be advantageous in that an adjoint equation may be
prescribed independently of the discretization scheme. Note that in the
discrete approach, the adjoint operator is obtained from the matrix transpose
of the discretized forward operator, which implies that the same spatial and
velocity resolution parameters must be used for both the forward and adjoint
solutions. In this Chapter, we will employ the same discretization parameters
for both the adjoint and forward problems, but this restriction is not
required for the continuous approach.

###  4.4 Implementation and benchmarks

The adjoint method has been implemented in the SFINCS code111The adjoint
method is implemented in the main branch of the SFINCS code
https://github.com/landreman/sfincs. using both the discrete and continuous
approaches. The magnetic geometry is specified in Boozer coordinates (Appendix
A.4) such that the covariant form of the magnetic field is,

$\displaystyle\textbf{B}=I(\psi)\nabla\vartheta_{B}+G(\psi)\nabla\varphi_{B}+K(\psi,\vartheta_{B},\varphi_{B})\nabla\psi,$
(4.59)

where $I(\psi)=\mu_{0}I_{T}(\psi)/2\pi$ and $G(\psi)=\mu_{0}I_{P}(\psi)/2\pi$,
$I_{T}(\psi)$ is the toroidal current enclosed by $\psi$, and $I_{P}(\psi)$ is
the poloidal current outside of $\psi$. The contravariant form is,

$\displaystyle\textbf{B}=\nabla\psi\times\nabla\vartheta_{B}-\iota(\psi)\nabla\psi\times\nabla\varphi_{B},$
(4.60)

where $\iota(\psi)$ is the rotational transform. The Jacobian is obtained from
dotting (4.59) with (4.60),

$\displaystyle\sqrt{g}=\frac{G(\psi)+\iota(\psi)I(\psi)}{B^{2}}.$ (4.61)

As $K(\psi,\vartheta_{B},\varphi_{B})$ does not appear in any of the
trajectory coefficients ((D.2) and (D.4)), in the drive term in (D.1), or in
the geometric factors used to define the moments of the distribution function
((4.5), (4.6), and (4.7)), all the geometric dependence enters through
$B(\psi,\vartheta_{B},\varphi_{B})$, $G(\psi)$, $I(\psi)$, and $\iota(\psi)$.
We choose to use Boozer coordinates for these computations as it reduces the
number of geometric parameters that must be considered, but the neoclassical
adjoint method is not limited to this choice of coordinate system.

We approximate $B$ by a truncated Fourier series,

$\displaystyle B=\sum_{m,n}B_{m,n}^{c}\cos(m\vartheta_{B}-nN_{P}\varphi_{B}),$
(4.62)

where the sum is taken over Fourier modes $m\leq m_{\max}$ and $|n|\leq
n_{\max}$ and $N_{P}$ is the number of periods. In (4.62), we have assumed
stellarator symmetry such that
$B(-\vartheta_{B},-\varphi_{B})=B(\vartheta_{B},\varphi_{B})$, and $N_{p}$
symmetry such that
$B(\vartheta_{B},\varphi_{B}+2\pi/N_{P})=B(\vartheta_{B},\varphi_{B})$. Thus
we compute derivatives with respect to the parameters
$\Omega=\\{B_{m,n}^{c},I(\psi),G(\psi),\iota(\psi)$}. Additionally,
derivatives with respect to $E_{r}$ are computed, which are used for efficient
ambipolar solutions and computing derivatives of geometric quantities at
ambipolarity (Section 4.5.3) rather than at fixed $E_{r}$.

To demonstrate, we compute $\partial\mathcal{R}/\partial B_{0,0}^{c}$ for
moments of the ion distribution function using the discrete and continuous
adjoint methods. A 3-mode model of the standard configuration W7-X geometry at
$\rho=\sqrt{\psi/\psi_{0}}=0.5$ is used (Table 1 in [16]),

$\displaystyle
B=B_{0,0}^{c}+B_{0,1}^{c}\cos(N_{P}\varphi_{B})+B_{1,1}^{c}\cos(\vartheta_{B}-N_{P}\varphi_{B})+B_{1,0}^{c}\cos(\vartheta_{B}),$
(4.63)

where $B_{0,1}^{c}=0.04645B_{0,0}^{c}$, $B_{1,1}^{c}=-0.04351B_{0,0}^{c}$, and
$B_{1,0}^{c}=-0.01902B_{0,0}^{c}$. Electron and ion ($q_{i}=e$) species are
included, and the derivatives are computed at the ambipolar $E_{r}$ with the
full trajectory model. The derivatives are also computed with a forward-
difference approach with varying step size $\Delta B_{0,0}^{c}$. In Figure 4.1
we show the fractional-difference between $\partial\mathcal{R}/\partial
B_{0,0}^{c}$ computed using the adjoint method and with forward-difference
derivatives. We see that at large values of $\Delta B_{0,0}^{c}$, the adjoint
and numerical derivatives begin to differ significantly due to discretization
error from the forward-difference approximation. The fractional error
decreases proportional to $\Delta B_{0,0}^{c}$ as expected until the rounding
error begins to dominate [203] when $\Delta B_{0,0}^{c}/B_{0,0}^{c}$ is
approximately $10^{-4}$, where $B_{0,0}^{c}$ is the value of the unperturbed
mode. The discrete and continuous approaches show qualitatively similar
trends. However, the minimum fractional difference is lower in the discrete
approach due to the additional discretization error that arises with the
continuous approach. With sufficient resolution parameters (41 $\theta$ grid
points, 61 $\phi$ grid points, 85 $\xi$ basis functions, and 7 $X$ basis
functions), the fractional error of the continuous approach is $\leq 0.1\%$
and should not be significant for most applications. We find similar agreement
for other derivatives and with the DKES trajectory model.

To demonstrate that the discrete and continuous methods indeed produce the
same derivative information, we compute the fractional difference between the
derivatives computed with the two methods as a function of the resolution
parameters. As an example, in Figure 4.2(a) we show the fractional difference
in $\partial Q_{i}/\partial\iota$, where $Q_{i}$ is the radial ion heat flux,
as a function of the number of Legendre polynomials used for the pitch angle
discretization, $N_{\xi}$, keeping the other resolution parameters fixed. As
$N_{\xi}$ is increased, the fractional differences converge to a finite value,
approximately $10^{-4}$, due to the discretization error in the other
resolution parameters. Similar resolution parameters are required for the
convergence of the moment itself, $Q_{i}$, and its derivative computed with
the continuous method, $\partial Q_{i}/\partial\iota$. Convergence of $Q_{i}$
within 5% is obtained with $N_{\xi}=38$, similar to that required for the
convergence of $\partial Q/\partial\iota$, as can be seen in Figure 4.2(a).

In Figure 4.2(b), we compare the cost of calculating derivatives of one moment
with respect to $N_{\Omega}$ parameters using the continuous and discrete
adjoint methods and forward-difference derivatives. All computations are
performed on the Edison computer at NERSC using 48 processors, and the elapsed
wall time is reported. Here we include the cost of solving the linear system
and computing diagnostics $N_{\Omega}+1$ times for the forward-difference
approach, and the cost of solving the forward and adjoint linear systems and
computing diagnostics for the adjoint approaches. The cost of the continuous
approach is slightly more than that of the discrete approach due to the cost
of factorizing the adjoint preconditioner. However, at large $N_{\Omega}$ the
cost of computing diagnostics for the adjoint approach (e.g., computing
$\partial\mathbb{S}/\partial\Omega$ and $\partial\mathbb{L}/\partial\Omega$
and performing the inner product in (4.45)) dominates that of solving the
adjoint linear system; thus the discrete and continuous approaches become
comparable in cost. In this regime, the adjoint approach provides speed-up by
a factor of approximately $50$.

(a) Discrete approach

(b) Continuous approach

Figure 4.1: Fractional difference between derivatives with respect to
$B_{0,0}^{c}$ computed with the adjoint method and with a forward-difference
derivative with step size $\Delta B_{0,0}^{c}$. The full trajectory model was
used with (a) the discrete and (b) the continuous adjoint approaches. Figure
adapted from [186] with permission.

(a)

(b)

Figure 4.2: (a) The fractional difference between $\partial
Q_{i}/\partial\iota$ computed with the continuous and discrete approaches
converges with the number of pitch angle Legendre modes, $N_{\xi}$. (b)
Comparison of the computational cost of computing
$\partial\mathcal{R}/\partial\Omega$ with forward-difference derivatives and
the adjoint approach as a function of $N_{\Omega}$, the number of parameters
in the gradient. Figure reproduced from [186] with permission.

###  4.5 Applications of the adjoint method

####  4.5.1 Local magnetic sensitivity analysis

With the adjoint method, it is possible to compute derivatives of a moment of
the distribution function with respect to the Fourier amplitudes of the field
strength, $\\{\partial\mathcal{R}/\partial B_{m,n}^{c}\\}$. Rather than
consider sensitivity in Fourier space, we would like to compute the
sensitivity to local perturbations of the field strength. We now quantify the
relationship between these two representations of sensitivity information.

Consider the Gateaux functional derivative [52] of $\mathcal{R}$ with respect
to $B$,

$\displaystyle\delta\mathcal{R}(B(\textbf{x});\delta
B)=\lim_{\epsilon\rightarrow 0}\frac{\mathcal{R}(B(\textbf{x})+\epsilon\delta
B(\textbf{x}))-\mathcal{R}(B(\textbf{x}))}{\epsilon}.$ (4.64)

Here the field strength is perturbed at fixed $I(\psi)$, $G(\psi)$, and
$\iota(\psi)$. As $\delta\mathcal{R}(B(\textbf{x});\delta B)$ is a linear
functional of $\delta B$, by the Riesz representation theorem [199],
$\delta\mathcal{R}$ can be expressed as an inner product with $\delta B$ and
some element of the appropriate space. The function $\delta B$ is defined on a
flux surface, $\psi$; thus it is sensible to express $\delta\mathcal{R}$ in
the following way,

$\displaystyle\delta\mathcal{R}(B(\textbf{x});\delta B)=\left\langle
S_{\mathcal{R}}\delta B(\textbf{x})\right\rangle_{\psi}.$ (4.65)

Here $\delta\mathcal{R}$ quantifies the change in the moment $\mathcal{R}$
associated with a local perturbation to the field strength, $\delta
B(\textbf{x})$. The function $S_{\mathcal{R}}$ is analogous to the shape
gradient introduced in Section 2.1, which will be discussed further in Section
4.5.2.

Suppose that $B$ is stellarator symmetric and $N_{P}$ symmetric. If $E_{r}=0$,
then $S_{\mathcal{R}}$ must also possess stellarator and $N_{P}$ symmetry
(Appendix G). However, when $E_{r}\neq 0$, $S_{\mathcal{R}}$ is no longer
guaranteed to have stellarator symmetry. Nonetheless, it may be desirable to
ignore the stellarator-asymmetric part of $S_{\mathcal{R}}$ if an optimized
stellarator-symmetric configuration is desired. For the remainder of this
Chapter, we will make this assumption, though the analysis could be extended
to consider the effect of breaking of stellarator symmetry. A truncated
Fourier series can approximate the quantity $S_{\mathcal{R}}$ under these
assumptions,

$\displaystyle
S_{\mathcal{R}}=\sum_{m,n}S_{m,n}\cos(m\vartheta_{B}-nN_{P}\varphi_{B}),$
(4.66)

where the sum is taken over $m\leq m_{\max}$ and $|n|\leq n_{\max}$. The
quantity $\delta B(\textbf{x})$ can be written in terms of perturbations to
the Fourier coefficients,

$\displaystyle\delta B(\textbf{x})=\sum_{m,n}\delta
B_{m,n}^{c}\cos(m\vartheta_{B}-nN_{P}\varphi_{B}),$ (4.67)

and now $\delta\mathcal{R}$ can be written in terms of these perturbations to
the Fourier coefficients,

$\displaystyle\delta\mathcal{R}=\sum_{m,n}\frac{\partial\mathcal{R}}{\partial
B_{m,n}^{c}}\delta B_{m,n}^{c}.$ (4.68)

In this way, (4.65) can be expressed as a linear system,

$\displaystyle\frac{\partial\mathcal{R}}{\partial
B_{m,n}^{c}}=\sum_{m^{\prime},n^{\prime}}D_{m,n;m^{\prime},n^{\prime}}S_{m^{\prime},n^{\prime}},$
(4.69)

where,

$D_{m,n;m^{\prime},n^{\prime}}=\\\
V^{\prime}(\psi)^{-1}\int_{0}^{2\pi}d\vartheta_{B}\int_{0}^{2\pi}d\varphi_{B}\,\sqrt{g}\cos(m\vartheta_{B}-nN_{P}\varphi_{B})\cos(m^{\prime}\vartheta_{B}-n^{\prime}N_{P}\varphi_{B}).$
(4.70)

If the same number of modes is used to discretize $\delta\mathcal{R}$ and
$S_{\mathcal{R}}$, then the linear system is square.

(a)

(b)

(c)

Figure 4.3: (a) The local magnetic sensitivity function for the bootstrap
current, $S_{J_{b}}$, is shown for the W7-X standard configuration. Positive
values indicate that increasing the field strength at a given location will
increase $J_{b}$ through (4.65). (b) The local sensitivity function for the
ion particle flux, $S_{\Gamma_{i}}$. (c) The magnetic field strength on the
$\rho=0.7$ surface. Figure adapted from [186] with permission.

In contrast to derivatives with respect to the Fourier modes of $B$, the
sensitivity function, $S_{\mathcal{R}}$, is a spatially local quantity,
quantifying the change in a figure of merit resulting from a local
perturbation of the field strength. In this way, $S_{\mathcal{R}}$ can inform
where perturbations to the magnetic field strength can be tolerated. The
sensitivity function could be related directly to a local magnetic tolerance,
as described in Section 2.1.3. In contrast with the work in [138], here we are
considering perturbations to the field strength on any flux surface rather
than at the plasma boundary. However, $S_{\mathcal{R}}$ still provides insight
into where trim coils should be placed or coil displacements can be tolerated
without sacrificing desired neoclassical properties. The sensitivity function
can also be used for gradient-based optimization in the space of the field
strength on a flux surface, as demonstrated in Section 4.5.2.

We compute $S_{J_{b}}$ for the W7-X standard configuration at $\rho=0.70$,
shown in Figure 4.3(a). We use a fixed-boundary equilibrium that preceded the
coil design and does not include coil ripple, and the full equilibrium is used
rather than the truncated Fourier series considered in Section 4.4. The same
resolution parameters are used as in Section 4.4, and derivatives with respect
to $B_{m,n}^{c}$ are computed for $m_{\max}=n_{\max}=20$. The largest modes
for this configuration are the helical curvature $B_{1,1}^{c}$, the toroidal
curvature $B_{1,0}^{c}$, and the toroidal mirror $B_{0,1}^{c}$. We find that
$S_{J_{b}}$ is large and negative on the inboard side, indicating that
increasing the magnitude of the toroidal curvature component of $B$ would lead
to an increase in $J_{b}$. This result is in agreement with previous analysis
[155], which found that at low collisionality, the bootstrap current
coefficients depend strongly on the toroidal curvature. Additionally, we note
a localized region of strong sensitivity on the inboard side near the bean-
shaped cross-section. Experimental [55] and numerical [75] evidence indicates
that the magnitude of the bootstrap current is increased in the lower mirror-
ratio configuration of W7-X, where the mirror-ratio is defined as
$(B_{\text{max}}-B_{\text{min}})/(B_{\text{max}}+B_{\text{min}})$. Our result
appears to be consistent with these observations: we note that the localized
region of strongly positive $S_{J_{b}}$ is near the maximum of the magnetic
field strength (Figure 4.3(c)), indicating that increasing the mirror-ratio
would lead to a decrease in the magnitude of bootstrap current, as $J_{b}<0$
for this configuration.

In Figure 4.3(b) is the sensitivity function for the ion particle flux,
$S_{\Gamma_{i}}$, computed for the same configuration using $m_{\max}=20$ and
$n_{\max}=20$. We find that the particle flux is more sensitive to
perturbations on the outboard side in localized regions, while on the inboard
side the sensitivity is relatively small in magnitude.

####  4.5.2 Gradient-based optimization

##### Optimization of the magnetic field strength

As a second demonstration of the adjoint neoclassical method, we consider
optimizing in the space of the field strength on a surface, taking
$\Omega=\\{B_{m,n}^{c}\\}$. As Boozer coordinates are used, the covariant form
(4.59) satisfies $(\nabla\times\textbf{B})\cdot\nabla\psi=0$ and the
contravariant form (4.60) satisfies $\nabla\cdot\textbf{B}=0$. As we will
artificially modify the field strength while keeping other geometry parameters
fixed, the resulting field will not necessarily satisfy both of these
conditions with both the covariant and contravariant forms. While there is no
guarantee that the resulting field strength will be consistent with a global
equilibrium solution, it provides insight into how local changes to the field
strength can impact neoclassical properties. As a second step, the outer
boundary could be optimized to match the desired field strength on a single
surface. In Section 4.5.2, we discuss how the derivatives computed in this
Chapter could be coupled to the optimization of an MHD equilibrium.

We perform optimization with a BFGS quasi-Newton method (Chapter 6 in [170])
using an objective function $\chi^{2}=J_{b}^{2}$, implemented in the
sfincs_adjoint branch of the STELLOPT code. A backtracking line search is used
at each iteration to find a step size that satisfies a condition of sufficient
decrease of $\chi^{2}$. We use the same equilibrium as in Section 4.5.1,
retaining modes $m\leq 12$ and $|n|\leq 12$, and compute derivatives with
respect to these modes. Convergence to $\chi^{2}\leq 10^{-10}$ was obtained
within 8 BFGS iterations (28 function evaluations), as shown in Figure 4.4(a).
The difference in field strength between the initial and optimized
configuration, $B_{\text{opt}}-B_{\text{init}}$, is shown in Figure 4.4(b). As
expected from the analysis in Section 4.5.1, the field strength increased on
the outboard side and decreased on the inboard side in comparison with
$B_{\text{init}}$. (Note that $J_{b}<0$.)

(a)

(b)

Figure 4.4: (a) Convergence of $\chi^{2}=J_{b}^{2}$ for optimization over
$\Omega=\\{B_{m,n}^{c}\\}$ with an adjoint-based BFGS method. (b) The change
in field strength from the initial to optimized configuration. Figure adapted
from [186] with permission.

##### Optimization of MHD equilibria

The local sensitivity function, $S_{\mathcal{R}}$, along with
$\partial\mathcal{R}/\partial I$, $\partial\mathcal{R}/\partial G$, and
$\partial\mathcal{R}/\partial\iota$, can be used to determine how
perturbations to the outer boundary of the plasma, $S_{P}$, result in
perturbations to $\mathcal{R}$. This is quantified through the idea of the
shape gradient, introduced in Section 2.1. The partial derivatives of
$\mathcal{R}$ can be computed with the adjoint method outlined in Section 4.3,
and the shape gradient can be obtained with only one additional MHD
equilibrium solution through the application of another adjoint method.

Consider a figure of merit which is integrated over the toroidal confinement
volume, $V_{P}$,

$\displaystyle
f_{\mathcal{R}}(S_{P})=\int_{V_{P}}d^{3}x\,w(\psi)\mathcal{R}(\psi),$ (4.71)

where $w(\psi)$ is a weighting function. That is, SFINCS is run on a set of
$\psi$ surfaces within $V_{P}$ and the volume integral is computed
numerically. Here we consider $S_{P}$ to be the plasma boundary used for a
fixed-boundary MHD equilibrium calculation. From the Hadamard-Zolesio
structure theorem (Section 2.1), the perturbation to $f_{\mathcal{R}}$
resulting from normal perturbation to $S_{P}$ can be written in the following
form,

$\displaystyle\delta
f_{\mathcal{R}}(S_{P};\delta\textbf{x})=\int_{S_{P}}d^{2}x\,\left(\delta\textbf{x}\cdot\hat{\textbf{n}}\right)\mathcal{G},$
(4.72)

under certain assumptions of smoothness [52]. This can be thought of as
another instance of the Riesz representation theorem, as $\delta
f_{\mathcal{R}}$ is a linear functional of $\delta\textbf{x}$. Here
$\hat{\textbf{n}}$ is the outward unit normal on $S_{P}$ and
$\delta\textbf{x}$ is a vector field describing the perturbation to the
surface. Intuitively, only normal perturbations to $S_{P}$ result in a change
to $f_{\mathcal{R}}$. The shape gradient is $\mathcal{G}$, which quantifies
the contribution of a local normal perturbation of the boundary to the change
in $f_{\mathcal{R}}$. The shape gradient can be used for fixed-boundary
optimization of equilibria or analysis of sensitivity to perturbations of
magnetic surfaces. It can be computed using a second adjoint method, where a
perturbed MHD force balance equation is solved with the addition of a bulk
force that depends on derivatives computed from the neoclassical adjoint
method. This will be described in detail in Chapter 5. While the continuous
neoclassical adjoint method described in this Chapter arises from the self-
adjointness of the linearized Fokker-Planck operator, the adjoint method for
MHD equilibria arises from the self-adjointness of the MHD force operator. In
practice, these two adjoint methods could be coupled by first computing an MHD
equilibrium solution, computing neoclassical transport and its geometric
derivatives from this equilibrium with the neoclassical adjoint method, and
passing these derivatives back to the equilibrium code to compute the shape
gradient with the perturbed MHD adjoint method. In this way, derivatives of
neoclassical quantities with respect to the shape of the outer boundary are
computed with only two equilibrium solutions and two DKE solutions.

Rather than solve an additional adjoint equation, the outer boundary could be
optimized by numerically computing derivatives of
$\\{B_{m,n}^{c}(\psi),G(\psi),I(\psi)\\}$ with respect to the double Fourier
series describing the outer boundary shape in cylindrical coordinates,
$\\{R_{m,n}^{c},Z_{m,n}^{s}\\}$, using a finite-difference method. This could
be done using the STELLOPT code [213, 197] with BOOZ_XFORM [202] to perform
the coordinate transformation. For example, if the rotational transform is
held fixed in the VMEC equilibrium calculation [111], the derivative of a
moment, $\mathcal{R}$, with respect to a boundary coefficient, $R_{m,n}^{c}$,
can be computed as,

$\displaystyle\frac{\partial\mathcal{R}(\psi)}{\partial
R_{m,n}^{c}(\psi)}=\sum_{m^{\prime},n^{\prime}}\frac{\partial\mathcal{R}(\psi)}{\partial
B_{m^{\prime},n^{\prime}}^{c}(\psi)}\frac{\partial
B_{m^{\prime},n^{\prime}}^{c}(\psi)}{\partial
R_{m,n}^{c}(\psi)}+\frac{\partial\mathcal{R}(\psi)}{\partial
G(\psi)}\frac{\partial G(\psi)}{\partial
R_{m,n}^{c}(\psi)}+\frac{\partial\mathcal{R}(\psi)}{\partial
I(\psi)}\frac{\partial I(\psi)}{\partial R_{m,n}^{c}(\psi)},$ (4.73)

where $\partial\mathcal{R}(\psi)/\partial B_{m,n}^{c}(\psi)$,
$\partial\mathcal{R}(\psi)/\partial G(\psi)$, and
$\partial\mathcal{R}(\psi)/\partial I(\psi)$ are computed with the
neoclassical adjoint method and $\partial B_{m,n}^{c}(\psi)/\partial
R_{m,n}^{c}(\psi)$, $\partial G(\psi)/\partial R_{m,n}^{c}(\psi)$, and
$\partial I(\psi)/\partial R_{m,n}^{c}(\psi)$ are computed with finite-
difference derivatives using STELLOPT. Similarly, derivatives of
$\\{B_{m,n}^{c}(\psi),G(\psi),I(\psi)\\}$ could be computed with respect to
coil parameters using a free-boundary equilibrium solution, allowing for
direct optimization of neoclassical quantities with respect to coil shapes.
The neoclassical calculation with SFINCS is typically significantly more
expensive than the equilibrium calculation (for the geometry discussed in
Section 4.5.1 fixed-boundary VMEC took 54 seconds while SFINCS took 157
seconds on 4 processors of the NERSC Edison computer). As such, combining
adjoint-based with finite-difference derivatives can still result in a
significant computational savings.

####  4.5.3 Ambipolarity

As stellarators are not intrinsically ambipolar, the radial electric field is
not truly an independent parameter. The ambipolar $E_{r}$ must be obtained
which satisfies the condition $J_{r}(E_{r})=0$. The application of adjoint-
based derivatives for computing the ambipolar solution is discussed in Section
4.5.3. An adjoint method to compute derivatives with respect to geometric
parameters at fixed ambipolarity is discussed in Section 4.5.3.

##### Accelerating ambipolar solve

A nonlinear root-finding algorithm must be used to compute the ambipolar
$E_{r}$. This root-finding can be accelerated with derivative information,
such as with a Newton-Raphson method [195]. The derivative required, $\partial
J_{r}/\partial E_{r}$, can be computed with the discrete or continuous adjoint
method as described in Section 4.3 with the replacement $\Omega_{i}\rightarrow
E_{r}$, considering $\mathcal{R}=J_{r}$.

We implement three nonlinear root finding methods: Brent’s method [30], the
Newton-Raphson method, and a hybrid between the bisection and Newton-Raphson
methods [195]. Brent’s method guarantees at least linear convergence by
combining quadratic interpolation with bisection and does not require
derivatives. The Newton-Raphson method can provide quadratic convergence under
certain assumptions but in general is not guaranteed to converge. If an
iterate lies near a stationary point or a poor initial guess is given, the
method can fail. For this reason, we implement the hybrid method, which
combines the possible quadratic convergence properties of Newton-Raphson with
the guaranteed linear convergence of the bisection method. Both Brent’s method
and the hybrid method require the root to be bracketed and therefore may
require additional function evaluations to obtain the bracket.

We compare these methods in Figure 4.5, using the W7-X standard configuration
considered in Section 4.5.1 with the full trajectory model and the discrete
adjoint approach, beginning with an initial guess of $E_{r}=-10$ kV/m with
bounds at $E_{r}^{\min}=-100$ kV/m and $E_{r}^{\max}=$ 100 kV/m. The root is
located at $E_{r}=-3.84$ kV/m. For this example, the hybrid and Newton methods
had nearly identical convergence properties. However, the Newton method is
less expensive as it does not require $J_{r}$ to be evaluated at the bounds of
the interval. The Newton method provides a 22% savings in wall clock time over
Brent’s method to obtain the root within the same tolerance.

In the above discussion, we have assumed that there is only one stable root of
interest. Of course, a given configuration may possess several roots,
especially if the ions and electrons are in different collisionality regimes
[92]. Multiple roots can be obtained by performing several root solves with
different initial values and brackets, which could be trivially parallelized.
Thus the adjoint method could still provide an acceleration in this more
general case.

Figure 4.5: The ambipolar root is obtained with Brent, Newton-Raphson, and
Newton hybrid nonlinear root solvers. The derivatives obtained with the
adjoint method provide better convergence properties for the Newton methods.
Figure adapted from [186] with permission.

##### Derivatives at ambipolarity

The adjoint method described in Section 4.3 assumes that $E_{r}$ is held
constant when computing derivatives with respect to $\Omega$. However, $E_{r}$
cannot truly be determined independently from geometric quantities, as the
ambipolar solution should be recomputed as the geometry is altered. It is
therefore desirable to compute derivatives at fixed ambipolarity (fixed
$J_{r}=0$) rather than at fixed $E_{r}$. This is performed by solving an
additional adjoint equation,

$\displaystyle\mathbb{L}^{\dagger}\lambda^{J_{r}}+\widetilde{J_{r}}=0,$ (4.74)

in the continuous approach or,

$\displaystyle\left(\overleftrightarrow{\textbf{L}}\right)^{T}\overrightarrow{\bm{\lambda}}^{J_{r}}=\overrightarrow{\textbf{J}_{r}},$
(4.75)

in the discrete approach. Details are described in Appendix H.

It should be noted that by computing derivatives at ambipolarity, we assume
that a given moment $\mathcal{R}$ is a differentiable function of the geometry
at fixed $J_{r}=0$. That is, this method cannot be applied to cases in which a
stable root disappears as the geometry varies. As this will occur at a
stationary point of $J_{r}(E_{r})$, this situation could be avoided within an
optimization loop by computing derivatives at constant $E_{r}$ rather than
constant $J_{r}$ if $|\partial J_{r}/\partial E_{r}|$ falls below a given
threshold at ambipolarity.

Although an additional adjoint solve is required, this method of computing
derivatives at ambipolarity is advantageous as several linear solves are
typically needed to obtain the ambipolar root. A comparison of the
computational cost between the adjoint method and the forward-difference
method for derivatives at ambipolarity is shown in Figure 4.6(a). Here the
full trajectory model is used, and the results for both the discrete and
continuous adjoint methods are shown. For the finite-difference derivative,
the ambipolar solve is performed with Brent’s method at each step in $\Omega$.
As in Figure 4.2(b), we find that for large $N_{\Omega}$, the cost of the
continuous and discrete approaches are essentially the same, as the cost is no
longer dominated by the linear solve. When computing the derivatives at
ambipolarity, both adjoint methods decrease the cost by a factor of
approximately $200$ for large $N_{\Omega}$.

In Figure 4.6(b) we show a benchmark between derivatives at ambipolarity,
$(\partial\mathcal{R}/\partial B_{0,0}^{c})_{J_{r}}$, computed with the
discrete adjoint method and with forward-difference derivatives. For the
forward-difference method, the Newton solver is used to obtain the ambipolar
$E_{r}$ as $B_{0,0}^{c}$ is varied. As the forward difference step size
$\Delta B_{0,0}^{c}$ decreases, the fractional difference again decreases
proportional to $\Delta B_{0,0}^{c}$ until it reaches a minimum when $\Delta
B_{0,0}^{c}/B_{0,0}^{c}$ is approximately $10^{-4}$. In comparison with Figure
4.1, we see that the minimum fractional difference is slightly larger at fixed
ambipolarity than at fixed $E_{r}$, as the tolerance parameters associated
with the Newton solver introduce an additional source of error to the forward-
difference approach.

(a)

(b)

Figure 4.6: (a) The cost of computing the gradient
$\partial\mathcal{R}/\partial\Omega$ at ambipolarity scales with $N_{\Omega}$,
the number of parameters in $\Omega$. (b) The fractional difference between
$\partial\mathcal{R}/\partial B_{0,0}^{c}$ at constant ambipolarity obtained
with the adjoint method and with finite-difference derivatives. Figure adapted
from [186] with permission.

In Figures 4.7(a) and 4.7(b) we compare the sensitivity function for the
particle flux, $S_{\Gamma_{i}}$, computed using derivatives at constant
$E_{r}$ with that computed at constant $J_{r}$. Here derivatives are computed
using the discrete adjoint method with full trajectories, and the sensitivity
function is constructed as described in Section 4.5.1. The configuration and
numerical parameters are the same as described in Section 4.5.1. At constant
$J_{r}$, the large region of increased sensitivity on the outboard side that
appears at constant $E_{r}$ remains, though the overall magnitude of the
sensitivity decreases. Thus it may be important to account for the effect of
the ambipolar $E_{r}$ when optimizing for radial transport. In Figures 4.7(c)
and 4.7(d) we perform the same comparison for $S_{J_{b}}$, finding the
derivatives at fixed $E_{r}$ and at fixed $J_{r}$ to be virtually identical.
This is to be expected, as numerical calculations of neoclassical transport
coefficients for W7-X have found that the bootstrap coefficients are much less
sensitive to $E_{r}$ than those for the radial transport (Figures 18 and 26 in
[16]). Furthermore, the bootstrap current in the $1/\nu$ regime is independent
of $E_{r}$, and the finite-collisionality correction is small for optimized
stellarators, such as W7-X [102]. Therefore, the ambipolarity corrections to
the derivatives are less important for $J_{b}$ than for the radial transport.

(a)

(b)

(c)

(d)

Figure 4.7: The sensitivity function for the ion particle flux,
$S_{\Gamma_{i}}$, is computed at (a) constant $E_{r}$ and (b) constant
$J_{r}$. Similarly, $S_{J_{b}}$ is computed at (c) constant $E_{r}$ and (d)
constant $J_{r}$. Figure adapted from [186] with permission.

###  4.6 Conclusions

We have described a method by which moments $\mathcal{R}$ of the neoclassical
distribution function can be differentiated efficiently with respect to many
parameters. The adjoint approach requires defining an inner product from which
the adjoint operator is obtained. We consider two choices for this inner
product. One choice corresponds with computing the adjoint of the linear
operator after discretization, and the other corresponds with computing it
before discretization. In the case of the former, the Euclidean dot product
can be used, and in the case of the latter, an inner product whose
corresponding norm is similar to the free energy norm (4.32) is defined. In
Section 4.4, we show that these approaches provide the same derivative
information within discretization error, as expected. Both methods provide a
reduction in computational cost by a factor of approximately $50$ in
comparison with forward-difference derivatives when differentiating with
respect to many ($\mathcal{O}(10^{2})$) parameters. In Section 4.5.3 the
adjoint method is extended to compute derivatives at ambipolarity. This method
provides a reduction in cost by a factor of approximately $200$ over a
forward-difference approach. We have implemented this method in the SFINCS
code, and similar methods could be applied to other drift kinetic solvers.

In this Chapter, we consider derivatives with respect to geometric quantities
that enter the DKE through Boozer coordinates. However, the adjoint
neoclassical method we have described is much more general, allowing for many
possible applications. For example, derivatives of the radial fluxes with
respect to the temperature and density profiles could be used to accelerate
the solution of the transport equations using a Newton method [13]. The
transport solution could furthermore be incorporated into the optimization
loop to self-consistently evolve the macroscopic profiles in the presence of
neoclassical fluxes. Rather than simply optimizing for minimal fluxes, an
objective function such as the total fusion power could be considered [107],
with optimization accelerated by adjoint-based derivatives.

Another application of the continuous adjoint formulation is the correction of
discretization error. The same solution obtained in Section 4.3.1 can be used
to quantify and correct for the error in a moment, $\mathcal{R}$, providing
similar accuracy to that computed with a higher-order stencil or finer mesh
without the associated cost. This method has been applied in the field of
computational fluid dynamics by solving adjoint Euler equations [231, 189] and
could prove useful for efficiently obtaining solutions of the DKE in low-
collisionality regimes.

In Section 4.5.2, we have shown an example of adjoint-based neoclassical
optimization, where the optimization space is taken to be the Fourier modes of
the field strength on a surface, $\\{B_{m,n}^{c}\\}$. While optimization
within this space is not necessarily consistent with a global equilibrium
solution, it demonstrates the adjoint neoclassical method for efficient
optimization. In Section 4.5.2, two approaches to self-consistently optimize
MHD equilibria are discussed. Further discussion and demonstration will be
provided in Chapter 5.

In Appendix G we show that when $E_{r}=0$ and the unperturbed geometry is
stellarator symmetric, the sensitivity functions for moments of the
distribution function are also stellarator symmetric. However, when $E_{r}\neq
0$ this is no longer true. This implies that obtaining minimal neoclassical
transport in the $\sqrt{\nu}$ regime may require breaking of stellarator
symmetry. In this Chapter, we have ignored the effects of stellarator
symmetry-breaking, though we hope to extend this work to study these effects
in the future.

## Chapter 5: Adjoint shape gradient for MHD equilibria

Most stellarator optimization to date has assumed that the magnetic field
satisfies the MHD equilibrium equations with either a fixed or free-boundary
approach, as detailed in Section 1.4.2. If a gradient-based optimization
approach is applied, derivatives of quantities that depend on the equilibrium
solutions must be computed with respect to the shapes of the filamentary coils
or plasma boundary. In this Chapter, we demonstrate an adjoint approach for
obtaining the coil or surface shape gradient of such functions. With the shape
gradient efficiently computed, shape derivatives with respect to any shape
perturbation can be calculated.

The material in this Chapter has been adapted with permission from [10] and
[187].

###  5.1 Introduction

Several figures of merit quantifying confinement must be considered in the
numerical optimization of stellarator MHD equilibrium. These figures of merit
describing a configuration depend on the shape of the outer plasma boundary or
the shape of the electro-magnetic coils. It is thus desirable to obtain
derivatives with respect to these shapes for optimization of equilibria or
identification of sensitivity information. These so-called shape derivatives
can be computed by directly perturbing the shape, recomputing the equilibrium,
and computing the resulting change to a figure of merit that depends on the
equilibrium solution. However, this direct finite-difference approach requires
recomputing the equilibrium for each possible perturbation of the shape. For
stellarators whose geometry is described by a set of $N_{\Omega}\sim 10^{2}$
parameters, this requires $N_{\Omega}$ solutions to the MHD equilibrium
equations. Despite this computational complexity, gradient-based optimization
of stellarators has proceeded with the direct approach (e.g. [197, 134, 196]).

As the target optimized configuration can never be realized exactly, an
analysis of the sensitivity to perturbations, such as errors in coil
fabrication or assembly, is central to the success of a stellarator. Tight
tolerances have proven to be a significant driver of the cost of stellarator
experiments [220, 130]; thus an improvement to the algorithms used to conduct
sensitivity studies can have a substantial impact on the field. In studies of
the coil tolerances for flux surface quality of LHD [240] and NCSX [31, 236],
perturbations of several distributions were manually applied to each coil.
Sensitivity analysis can also be performed with analytic derivatives.
Numerical derivatives with respect to tilt angle and coil translation of the
CNT coils have been used to compute the sensitivity of the rotational
transform on axis [88]. Analytic derivatives have recently been applied to
study coil sensitivities of the CNT stellarator by considering the
eigenvectors of the Hessian matrix [243]. Thus, in addition to gradient-based
optimization, derivatives with respect to shape can be applied to sensitivity
analysis.

The shape gradient quantifies the change in a figure of merit associated with
a local perturbation to a shape. Thus, if the shape gradient can be obtained,
the shape derivative with respect to any perturbation is known (more precise
definitions of the shape derivative and gradient are given in Sections 2.1 and
5.2). The shape gradient representation can be computed from parameter
derivatives by solving a small linear system (Sections 2.1.2). However,
computing parameter derivatives can often be computationally expensive, as
numerical derivatives require evaluating the objective function at least
$N_{\Omega}+1$ times if one-sided finite-difference derivatives are used, or
$2N_{\Omega}$ times for centered differences. As computing the objective
function often involves solving a linear or nonlinear system, such as the MHD
equilibrium equations, this implies solving the system of equations $\geq
N_{\Omega}+1$ times. Numerical derivatives also introduce additional noise,
and the finite-difference step size must be chosen carefully.

Rather than use parameter derivatives, in this Chapter we will use an adjoint
method to compute the shape gradient. This is sometimes termed adjoint shape
sensitivity or adjoint shape optimization, which has its origins in
aerodynamic engineering and computational fluid dynamics [190, 82]. As with
adjoint methods for parameter derivatives, this technique only requires the
solution of two linear or nonlinear systems of equations. This technique has
been applied to magnetic confinement fusion for the design of tokamak divertor
shapes by solving forward and adjoint fluid equations [48, 49, 50]. As
stellarators require many parameters to describe their shape, adjoint shape
sensitivity could significantly decrease the cost of computing the shape
gradient. If one is optimizing in the space of parameters describing the
boundary of the plasma or the shape of coils, the shape gradient
representation obtained from the adjoint method can be converted to parameter
derivatives upon multiplication with a small matrix (Section 2.1).

We begin in Section 5.2 with a brief review of shape calculus concepts in the
context of MHD equilibria. In Section 5.3, the fundamental adjoint relations
for perturbations to MHD equilibria are derived and discussed. These relations
take a form that is similar to that of transport coefficients that are related
by Onsager symmetry [177, 178]. Specifically, perturbations to the equilibrium
are characterized as a set of generalized responses to a complementary set of
generalized forces. The responses and forces can be thought of as being
related by a matrix operator, which is symmetric. The resulting relations
among forces and responses can be used to compute the shape gradient of
functions of the equilibria with respect to displacements of the plasma
boundary or the coil shapes. In Section 5.4, the continuous adjoint method
that takes advantage of the generalized self-adjointness relations is
discussed. Several applications to stellarator figures of merit will be
demonstrated in Section 5.5.

Although the adjoint relations are based on the equations of linearized MHD,
we perform numerical calculations in this Chapter with nonlinear MHD solutions
with the addition of a small perturbation. Demonstration is performed using
nonlinear stellarator MHD equilibrium codes based on a variational principle,
VMEC [111] and ANIMEC [43]. We obtain expressions for the shape gradients of
the volume-averaged $\beta$ (Section 5.5.1), rotational transform (Section
5.5.2), vacuum magnetic well (Section 5.5.3), magnetic ripple (Section 5.5.4),
effective ripple in the $1/\nu$ neoclassical regime [168] where $\nu$ is the
collision frequency (Section 5.5.5), and departure from quasi-symmetry
(Section 5.5.6). Finally, we demonstrate that the adjoint method for
neoclassical optimization outlined in Chapter 4 can be coupled with a
linearized adjoint MHD solution to compute derivatives of several neoclassical
quantities with respect to the shape of the plasma boundary (Section 5.5.7).
We present calculations of the shape gradient with the adjoint approach for
the volume-averaged $\beta$, rotational transform, and vacuum magnetic well
figures of merit, which do not require modification to VMEC. The calculation
for the magnetic ripple is computed with a minor modification of the ANIMEC
code. The adjoint force balance equations needed to compute the shape gradient
for the other figures of merit require the addition of a bulk force that will
necessitate further modification of an equilibrium or linearized MHD code.
Numerical calculations for these figures of merit will, therefore, not be
presented in this Chapter.

###  5.2 Shape calculus review

We now review shape calculus fundamentals introduced in Chapter 2 in the
context of functions that depend on MHD equilibrium quantities. Consider a
functional, $F(S_{P})$, that depends implicitly on the plasma boundary,
$S_{P}$, through the solution to the fixed-boundary MHD equilibrium equations
(Section 1.4.1) with boundary condition
$\textbf{B}\cdot\hat{\textbf{n}}|_{S_{P}}=0$ where $\hat{\textbf{n}}$ is the
outward unit normal on $S_{P}$. We define a functional integrated over the
plasma volume, $V_{P}$,

$\displaystyle f(S_{P})=\int_{V_{P}}d^{3}x\,F(S_{P}),$ (5.1)

where $S_{P}$ is the boundary of $V_{P}$. Consider a vector field describing
displacements of the surface, $\delta\textbf{x}$, and a displaced surface
$S_{P,\epsilon}=\\{\textbf{x}_{0}+\epsilon\delta\textbf{x}\mathrel{\mathop{\ordinarycolon}}\textbf{x}_{0}\in
S_{P}\\}$. The shape derivative of $F$ is defined as,

$\displaystyle\delta F(S_{P};\delta\textbf{x})=\lim_{\epsilon\rightarrow
0}\frac{F(S_{P,\epsilon})-F(S_{P})}{\epsilon}.$ (5.2)

The shape derivative of $f$ is defined by the same expression with $F\to f$.
Under certain assumptions of smoothness of $\delta F$ with respect to
$\delta\textbf{x}$, the shape derivative of the volume-integrated quantity,
$f$, can be written in the following way (Section 2.1),

$\displaystyle\delta f(S_{P};\delta\textbf{x})=\int_{V_{P}}d^{3}x\,\delta
F(S_{P};\delta\textbf{x})+\int_{S_{P}}d^{2}x\,\delta\textbf{x}\cdot\hat{\textbf{n}}F.$
(5.3)

The first term accounts for the Eulerian perturbation to $F$ while the second
accounts for the motion of the boundary. This is referred to as the transport
theorem for domain functionals and will be used throughout this Chapter to
compute the shape derivatives of figures of merit of interest.

According to the Hadamard-Zolesio structure theorem [52], the shape derivative
of a functional of $S_{P}$ (not restricted to the form of (5.1)) can be
written in the following form,

$\displaystyle\delta
f(S_{P};\delta\textbf{x})=\int_{S_{P}}d^{2}x\,\delta\textbf{x}\cdot\hat{\textbf{n}}\mathcal{G},$
(5.4)

assuming $\delta f$ exists for all $\delta\textbf{x}$ and is sufficiently
smooth. In the above expression, $\mathcal{G}$ is the shape gradient. This is
an instance of the Riesz representation theorem, which states that any linear
functional can be expressed as an inner product with an element of the
appropriate space [199]. As the shape derivative of $f$ is linear in
$\delta\textbf{x}$, it can be written in the form of (5.4). Intuitively, the
shape derivative does not depend on tangential perturbations to the surface.
The shape gradient can be computed from derivatives with respect to the set of
parameters, $\Omega$, used to discretize $S_{P}$,

$\displaystyle\frac{\partial
f(\Omega)}{\partial\Omega_{i}}=\int_{S_{P}}d^{2}x\,\frac{\partial\textbf{x}(\Omega)}{\partial\Omega_{i}}\cdot\hat{\textbf{n}}\mathcal{G}.$
(5.5)

For example, $\Omega=\\{R_{m,n}^{c},Z_{m,n}^{s}\\}$ could be assumed, where
these are the Fourier coefficients (5.70) in a cosine and sine representation
of the cylindrical coordinates $(R,Z)$ of $S_{P}$. Upon discretization of the
right-hand side on a surface, the above takes the form of a linear system that
can be solved for $\mathcal{G}$ [138]. However, this approach requires
performing at least one additional equilibrium calculation for each parameter
with a finite-difference approach.

The shape gradient can also be computed with respect to perturbations of
currents in the vacuum region. We now consider $f$ to depend on the shape of a
set of filamentary coils, $C=\\{C_{k}\\}$, through a free-boundary solution to
the MHD equilibrium equations (Section 1.4.1). We consider a vector field of
displacements to the coils, $\delta\textbf{x}_{C}$. The shape derivative of
$f$ can also be written in shape gradient form,

$\displaystyle\delta
f(C;\delta\textbf{x}_{C})=\sum_{k}\oint_{C_{k}}dl\,\delta\textbf{x}_{C_{k}}\cdot\widetilde{\bm{\mathcal{G}}}_{k},$
(5.6)

where $\widetilde{\bm{\mathcal{G}}}_{k}$ is the shape gradient for coil $k$,
$C_{k}$ is the line integral along coil $k$, and the sum is taken over coils.
Again, $\widetilde{\bm{\mathcal{G}}}_{k}$ can be computed from derivatives
with respect to a set of a parameters describing coil shapes (5.84), analogous
to (5.5). Note that we have defined the shape gradient in a slightly different
way here than that introduced in Chapter 2 (2.15) (without the cross with
$\hat{\textbf{t}}$), although we will find in this Chapter that
$\widetilde{\bm{\mathcal{G}}}_{k}$ is perpendicular to $\hat{\textbf{t}}$ for
the functionals under consideration. We distinguish the shape gradient as
defined in (5.6) from that defined in (2.15) with a tilde.

To avoid the cost of direct computation of the shape gradient, we apply an
adjoint approach. The shape gradient is thus obtained without perturbing the
plasma surface or coil shapes directly, but instead by solving an additional
adjoint equation that depends on the figure of merit of interest. We perform
the calculation with the direct approach to demonstrate that the same
derivative information is computed with either method.

###  5.3 Adjoint relations for MHD equilibria

The goal of this Section is to generalize the well-known self-adjointness [20]
of the MHD force operator,

$\displaystyle\int_{V_{P}}d^{3}x\,\left(\bm{\xi}_{1}\cdot\textbf{F}[\bm{\xi}_{2}]-\bm{\xi}_{2}\cdot\textbf{F}[\bm{\xi}_{1}]\right)-\frac{1}{\mu_{0}}\int_{S_{P}}d^{2}x\,\hat{\textbf{n}}\cdot\left(\bm{\xi}_{2}\delta\textbf{B}[\bm{\xi}_{1}]\cdot\textbf{B}-\bm{\xi}_{1}\delta\textbf{B}[\bm{\xi}_{2}]\cdot\textbf{B}\right)=0,$
(5.7)

to allow for perturbations of interest for stellarator optimization. In this
expression, the perturbed magnetic field is expressed in terms of the
displacement vector,

$\displaystyle\delta\textbf{B}[\bm{\xi}_{1,2}]=\nabla\times\left(\bm{\xi}_{1,2}\times\textbf{B}\right),$
(5.8)

which follows from the assumption that the rotational transform is fixed by
the perturbation (flux-freezing). The MHD force operator,

$\displaystyle\textbf{F}[\bm{\xi}_{1,2}]=\frac{\left(\nabla\times\delta\textbf{B}[\bm{\xi}_{1,2}]\right)\times\textbf{B}}{\mu_{0}}+\frac{\left(\nabla\times\textbf{B}\right)\times\delta\textbf{B}[\bm{\xi}_{1,2}]}{\mu_{0}}-\nabla\left(\delta
p[\bm{\xi}_{1,2}]\right),$ (5.9)

is a linearization of the MHD equilibrium equation,

$\displaystyle\frac{\left(\nabla\times\textbf{B}\right)\times\textbf{B}}{\mu_{0}}=\nabla
p,$ (5.10)

with boundary condition,

$\displaystyle\textbf{B}\cdot\hat{\textbf{n}}\rvert_{S_{P}}=0,$ (5.11)

under the assumption that the magnetic field is perturbed according to (5.8)
and the pressure is perturbed according to,

$\displaystyle\delta p[\bm{\xi}_{1,2}]=-\bm{\xi}_{1,2}\cdot\nabla p-\gamma
p\nabla\cdot\bm{\xi}_{1,2},$ (5.12)

where $\gamma$ is the adiabatic index. As $\bm{\xi}$ describes the motion of
field lines, modes which perturb the plasma boundary exhibit non-zero
$\bm{\xi}\cdot\hat{\textbf{n}}\rvert_{S_{P}}$. The self-adjointness provides a
relationship between two perturbations about an MHD equilibrium state
described by (5.10)-(5.11). This relation is incredibly valuable for ideal MHD
stability analysis, forming the basis for the energy principle.

As described in Section 2.2.2, when formulating a continuous adjoint approach,
the adjoint of the linearized operator appearing in the forward PDE must be
obtained. However, we cannot directly apply the self-adjointness relation from
MHD stability theory (5.7) for the stellarator optimization problem. While MHD
perturbations assume fixed rotational transform, stellarator optimization is
often performed instead at fixed toroidal current. While the MHD self-
adjointness relation allows for perturbations of the plasma boundary, we would
also like to consider linearized equilibrium states corresponding to
perturbations of coils in the vacuum region. We now form the appropriate
generalized self-adjointness relations corresponding to fixed-boundary
perturbations (applied perturbations to the plasma boundary) and free-boundary
perturbations (applied perturbations to electro-magnetic coils). Even though
the boundary shape changes in the former case, we refer to it as “fixed
boundary” since the equilibrium code is run in fixed-boundary mode, and since
the associated adjoint problem will turn out to have no boundary perturbation.

The resulting expressions will allow us to relate the “direct perturbations,”
those corresponding to a linearized equilibrium state associated with the
direct perturbation of the plasma boundary or coil shapes, and “adjoint
perturbations,” with which we can compute the shape gradient efficiently. The
adjoint perturbation will correspond to the change in the equilibrium when an
additional bulk force acts on the plasma or the toroidal current profile is
changed. For the adjoint perturbation, there is no change to the outer flux
surface in the fixed-boundary case or to the coil currents in the free-
boundary case. In this Section, we will show that aspects of the direct and
adjoint changes are related to each other in a manner similar to Onsager
symmetry. Thus, it will be shown that by calculating the adjoint perturbation,
with a judiciously chosen added force or change in the toroidal current
profile, the solution to the direct problem can be determined.

We consider equilibria in which the magnetic field in the plasma can be
expressed in terms of scalar functions
$\psi(\textbf{x}),\chi(\psi),\vartheta(\textbf{x})$, and
$\varphi(\textbf{x})$,

$\displaystyle\textbf{B}=\nabla\psi\times\nabla\vartheta-\nabla\chi\times\nabla\varphi=\nabla\psi\times\nabla\alpha,$
(5.13)

where ($\psi$,$\vartheta$,$\varphi$) form any magnetic coordinate system
(Appendix A.3). We will regard $\psi$ as labeling the flux surfaces and
consider toroidal geometries for which,

$\displaystyle\alpha=\vartheta-\iota(\psi)\varphi,$ (5.14)

label field lines in a flux surface, where $\vartheta$ is a poloidal angle,
$\varphi$ is a toroidal angle, and $\iota(\psi)=\chi^{\prime}(\psi)$ is the
rotational transform, with $\chi(\psi)$ being the poloidal flux function. With
these definitions, the magnetic flux passing toroidally through a poloidally
closed curve of constant $\psi$ is $2\pi\psi$, and the flux passing poloidally
between the magnetic axis and the surface of constant $\psi$ is
$2\pi\chi(\psi)$. Thus, we assume that good flux surfaces exist and leave
aside the issues of islands and chaotic field lines. In addition to the
representation of the magnetic field, we assume that MHD force balance (5.10)
is satisfied with a scalar pressure, $p(\psi)$.

As mentioned, we will consider two cases, a fixed-boundary case in which the
shape of the outer flux surface is prescribed, and a free-boundary case for
which outside the plasma, whose surface is defined by a particular value of
toroidal flux, the force balance equation (5.10) does not apply, but rather,
the magnetic field is determined by Ampere’s law,

$\displaystyle\nabla\times\textbf{B}=\mu_{0}\textbf{J},$ (5.15)

with a given current density $\textbf{J}_{C}$, representing current flowing
outside the confinement region. The fixed-boundary and free-boundary equations
are discussed in detail in Section 1.4.1.

From (5.10) it follows that current density stream-lines also lie in the
$\psi=$ constant surfaces. The toroidal current passing through a surface,
$\mathcal{S}_{T}(\psi)$ (Figure A.2), whose perimeter is a closed poloidal
loop at constant $\psi$ is given by,

$\displaystyle
I_{T}(\psi)=\int_{\mathcal{S}_{T}(\psi)}d^{2}x\,\hat{\textbf{n}}\cdot\textbf{J}=\int_{\mathcal{S}_{T}(\psi)}d\psi\,d\vartheta\,\sqrt{g}\nabla\varphi\cdot\textbf{J},$
(5.16)

where $\sqrt{g}^{-1}=\nabla\psi\times\nabla\vartheta\cdot\nabla\phi$.

Equations 5.13, 5.14, 5.14, 5.10, 5.15 and 5.16 describe our base equilibrium