2. Overview of the current potential method

When optimizing coils, a balance must be reached between a simple coil design and accuracy in reproducing the desired plasma shape. The REGCOIL code (Landreman Reference Landreman2017) uses a current potential approximation to formulate this optimization as a linear least-squares system with a regularization factor balancing the coil complexity and plasma surface accuracy. The surface current on the coil winding surface is taken to be of the form

(2.1)\begin{equation} \boldsymbol{K}=\boldsymbol{\hat{n}}' \times \boldsymbol{\nabla} \varPhi, \end{equation}

where $\varPhi$ is the current potential and $\boldsymbol {\hat {n}}'$ is the outward unit normal vector to the coil surface. The surface current is then divergence-free and tangential to the coil surface. Since the gradient and level curves of the current potential are perpendicular and on the coil surface, the surface current flows along the level curves of the current potential. The discrete coil shapes are then determined by taking a selection of level curves of the current potential, or streamlines of $\boldsymbol {K}$. These level curves will be closer together in regions of high current density, implying a reduction of the discrete coil–coil separation.

This current potential has secular terms and a single-valued term. The single-valued component of the current potential can be expanded as a sine series over the coil surface under the assumption of stellarator symmetry (Dewar & Hudson Reference Dewar and Hudson1998). The current potential then has the form

(2.2)\begin{equation} \varPhi = G\frac{\zeta'}{2{\rm \pi}} + I\frac{\theta'}{2{\rm \pi}} + \sum_j \varPhi_j \sin(m_j\theta' - n_j\zeta'), \end{equation}

where $G$ and $I$ are currents which link the coil surface poloidally and toroidally, $\theta '$ is a poloidal angle and $\zeta '$ is a toroidal angle.

We consider two figures of merit:

(2.3)\begin{equation} \left. \begin{aligned} \chi^{2}_B &=\int_{S_{\text{plasma}}} (\boldsymbol{B}(\theta,\zeta)\boldsymbol{\cdot} \boldsymbol{\hat{n}}(\theta,\zeta) )^{2} \,\mathrm{d}a ,\\ \chi^{2}_K &=\int_{S_{\text{coil}}} |\boldsymbol{K}(\theta',\zeta') |^{2} \,\mathrm{d}a', \end{aligned} \right\} \end{equation}

where $\boldsymbol {B}$ is the magnetic field and $\boldsymbol {\hat {n}}$ is the outward unit normal vector to the plasma surface. In keeping with Landreman (Reference Landreman2017), primed coordinates refer to the coil surface while unprimed coordinates refer to the plasma surface. Similarly, primed and unprimed quantities are used to refer to those of the coil and plasma surfaces, respectively, when a distinction is necessary. The integral for $\chi ^{2}_B$ is over the desired magnetic surface, $S_{\text {plasma}}$, and the integral for $\chi ^{2}_K$ is over the coil winding surface, $S_{\text {coil}}$. The goal is for $S_{\text {plasma}}$ to be a magnetic surface, so the normal component of the magnetic field will vanish if the plasma surface is properly reconstructed and, thus, so will $\chi ^{2}_B$. Large values of $\chi ^{2}_K$ are correlated with increased coil complexity (Paul et al. Reference Paul, Landreman, Bader and Dorland2018). Regions of higher surface current density require reduced coil–coil spacing; thus, the coil–coil separation is increased by reducing $\chi ^{2}_K$.

We make the approximation that the normal field from the plasma current is negligible such that only the vacuum field from the coils is included in (2.3). In § 6, examples are shown for the vacuum QHS46 and W7-X equilibria. If the derivatives of the normal field from the plasma current with respect to the plasma boundary are computed (Henneberg et al. Reference Henneberg, Hudson, Pfefferlé and Helander2020), the adjoint method we present in § 5 could be modified to include the dependence of $\chi ^{2}_B$ on the plasma boundary through the plasma normal field.

When comparing configurations, it is also useful to consider normalized quantities:

(2.4)\begin{equation} \left. \begin{aligned} \|B_n\|_2 &= \sqrt{\frac{\chi^{2}_B}{B_0^{2} a_{\text{plasma}}}}, \\ \|K\|_2 &= \sqrt{\frac{\chi^{2}_K}{a_{\text{coil}}}}, \end{aligned} \right\} \end{equation}

where $a_{\text {plasma}}$ and $a_{\text {coil}}$ are the plasma and coil surface areas, respectively, and $B_0$ is a constant with the dimensions of magnetic field. In this work, we take it to be the root-mean-squared average over the plasma surface of the magnetic field strength from the equilibrium.

The REGCOIL code minimizes the quantity $\chi ^{2}_B + \lambda \chi ^{2}_K$ by solving a linear least-squares system for the Fourier modes of the single-valued current potential, $\varPhi _j$, defined in (2.2). Here $\lambda$ is the regularization parameter which sets the relative importance of $\chi ^{2}_B$ and $\chi ^{2}_K$ in the optimization. A fixed value of $\lambda$ can be chosen, or the freedom in $\lambda$ can be used to fix a certain target function, such as $\chi ^{2}_K$ or $\|K\|_2$. We employ the notation $\boldsymbol {\varPhi }$ to represent the vector of unknowns, $\{\varPhi _j\}$. The linear least-squares system takes the form

(2.5)\begin{equation} {{{\mathsf{A}}}}\boldsymbol{\varPhi} = \boldsymbol{b}, \end{equation}

with

(2.6)\begin{equation} \left. \begin{gathered} {{{\mathsf{A}}}}={{{\mathsf{A}}}}^{B}+\lambda {{{\mathsf{A}}}}^{K},\\ \boldsymbol{b}=\boldsymbol{b}^{B}+\lambda \boldsymbol{b}^{K}, \end{gathered} \right\} \end{equation}

where the $B$ matrix and vector are associated with the system that optimizes only $\chi ^{2}_B$, and the $K$ matrix and vector are associated with the system that optimizes only $\chi ^{2}_K$. Expressions for these quantities are given in appendix A of Landreman (Reference Landreman2017).

3. Shape gradients

The shape gradient describes the local dependence of a functional on a surface and is independent of the choice of basis or parameterization. Consider a functional $F$ of the plasma surface, such as $\chi ^{2}_B$, $\chi ^{2}_K$, $\|B_n\|_2$ or $\|K\|_2$. The Hadamard–Zolésio structure theorem (Delfour & Zolésio Reference Delfour and Zolésio2011) states that the change in $F$ due to a vector field of infinitesimal displacements $\delta \boldsymbol {r}$ to the plasma surface $S_{\text {plasma}}$ is given by

(3.1)\begin{equation} \delta F(S_{\text{plasma}};\delta \boldsymbol{r}) = \int_{S_{\text{plasma}}} S_F \delta \boldsymbol{r} \boldsymbol{\cdot} \boldsymbol{\hat{n}} \, \mathrm{d}a, \end{equation}

under the assumption that $\delta F$ exists for all infinitely differentiable and compactly supported $\delta \boldsymbol {r}$. Here $S_F$ is the shape gradient and $\delta F$ is called the shape derivative. In this equation and the following ones, the integral is over the unperturbed plasma surface and $\boldsymbol {\hat {n}}$ is the unit normal outward to this surface. This gives information about how the plasma surface can be deformed to alter a functional of the surface. For the case that the functional $F$ is $\chi ^{2}_K$ or $\|K\|_2$, its shape gradient enables the plasma surface to be deformed to obtain simpler coils.

Given a parameterization of the plasma surface with a set of parameters $\boldsymbol {\varOmega }$, (3.1) can be expressed as

(3.2)\begin{equation} \frac{\partial F}{\partial \varOmega_j} = \int_{S_{\text{plasma}}} S_F \frac{\partial \boldsymbol{r}}{\partial \varOmega_j} \boldsymbol{\cdot} \boldsymbol{\hat{n}} \, \mathrm{d}a. \end{equation}

The shape gradient can be expanded as a Fourier series, assuming stellarator symmetry, as

(3.3)\begin{equation} S_F=\sum_i S_F^{i} \cos(m_i \theta - n_i \zeta), \end{equation}

such that the linear system for the shape gradient Fourier coefficients, $\{S_F^{i}\}$, is

(3.4)\begin{equation} \frac{\partial F}{\partial \varOmega_j} = \sum_i S_F^{i} \int_{S_{\text{plasma}}} \cos(m_i \theta - n_i \zeta) \frac{\partial \boldsymbol{r}}{\partial \varOmega_j} \boldsymbol{\cdot} \boldsymbol{\hat{n}} \, \mathrm{d}a. \end{equation}

If the same number of modes are used to discretize $S_F$ as the number of modes retained in $\{\partial F/\partial \varOmega _j\}$, then the system is square, and the shape gradient can then be constructed from a local parameterization of the surface with a linear solve. If the system is not square, it can still be solved (Landreman & Paul Reference Landreman and Paul2018) using a pseudoinverse or QR factorization. Details of the chosen parameterization of the plasma boundary are provided in the following section.

4. Efficient Fourier representation of the plasma boundary

A toroidal surface such as $S_{\text {plasma}}$ can be specified by two Fourier series in the cylindrical coordinates $R$ and $Z$:

(4.1)\begin{equation} \left. \begin{gathered} R(\theta,\zeta)=\sum_{m,n}R_{m,n}^{c} \cos(m\theta-n\zeta),\\ Z(\theta,\zeta)=\sum_{m,n}Z_{m,n}^{s} \sin(m\theta-n\zeta). \end{gathered} \right\} \end{equation}

Here $\zeta$ is taken to be the cylindrical toroidal angle and $\theta$ is a poloidal angle. As a Fourier series in two separate variables is required, we refer to this as the double Fourier representation. This representation is used in the VMEC (Hirshman & Whitson Reference Hirshman and Whitson1983) and SPEC (Hudson et al. Reference Hudson, Dewar, Dennis, Hole, McGann, Von Nessi and Lazerson2012) equilibrium codes. Here we have assumed stellarator symmetry such that

(4.2)\begin{equation} \left. \begin{gathered} R(-\theta,-\zeta)=R(\theta,\zeta),\\ Z(-\theta,-\zeta)={-}Z(\theta,\zeta). \end{gathered} \right\} \end{equation}

Because this representation does not constrain the poloidal angle, the Fourier amplitudes could be altered corresponding to a redefinition of the poloidal angle, while the surface is left unchanged. This freedom in the choice of poloidal angle has been used to condense the Fourier spectrum using a variational method (Hirshman & Meier Reference Hirshman and Meier1985; Hirshman & Breslau Reference Hirshman and Breslau1998). In the context of fixed-boundary optimization, a uniquely defined poloidal angle may eliminate a null space in the parameter space.

To eliminate such a null space, we consider a representation which requires only one series representation of a radial distance, $l$, from a fixed coordinate axis (Hirshman & Breslau Reference Hirshman and Breslau1998). In this work, we take the coordinate axis to coincide with the magnetic axis from the equilibrium, although this assumption is not necessary. We define a unique poloidal angle $\vartheta$ using the arctangent in the poloidal plane:

(4.3)\begin{equation} \left. \begin{gathered} l=\sqrt{(R-R_0)^{2}+(Z-Z_0)^{2}}, \\ \vartheta=\arctan\left(\frac{Z-Z_0}{R-R_0}\right), \end{gathered} \right\} \end{equation}

where $R_0(\zeta )$ and $Z_0(\zeta )$ are the cylindrical coordinates of the axis. Note that a representation of this form assumes that each cross-section in the poloidal plane is a star domain, indicating that there exists a coordinate axis $(R_0,Z_0)$ such that the line segment connecting the axis and any point on the boundary is contained within the boundary.

The variable $l$ then admits a Fourier representation as

(4.4)\begin{equation} l(\vartheta,\zeta)=\sum_{m,n}l_{m,n}^{c} \cos(m\vartheta-n\zeta). \end{equation}

This representation uses a Fourier expansion in only a single quantity $l$ rather than in both $R$ and $Z$; for this reason, we refer to this as the single Fourier representation. The cylindrical coordinates can now be expressed in the new representation as

(4.5)\begin{equation} \left. \begin{gathered} R(\vartheta,\zeta)=R_0(\zeta)+l(\vartheta,\zeta)\cos(\vartheta),\\ Z(\vartheta,\zeta)=Z_0(\zeta)+l(\vartheta,\zeta)\sin(\vartheta). \end{gathered} \right\} \end{equation}

In practice, the Fourier series in (4.4) may decay rather slowly for relevant stellarator boundaries with widely varying curvature, as equally spaced grid points in $\vartheta$ accumulate in regions of the surface where $l$ varies slowly with respect to $\vartheta$ and become sparse where $l$ varies rapidly with respect to $\vartheta$. In figure 1(b), this accumulation occurs on the inboard and outboard sides of the plasma surface.

Figure 1. The bean-shaped cross-section of the W7-X configuration with a 0.25 m uniformly offset surface sampled at 32 values of the poloidal using (a) the initial VMEC representation, (b) the single Fourier representation with an arctangent angle, (c) the single Fourier representation with an arclength angle on the plasma surface and (d) the single Fourier representation with an arclength angle on a surface uniformly offset from the plasma surface by 0.25 m.

To improve the efficiency of this representation, a second unique poloidal angle can be defined by the arclength along the plasma surface at constant toroidal angle:

(4.6)\begin{equation} \bar{\theta} (\vartheta,\zeta)=2{\rm \pi} \dfrac{\displaystyle\int_0^{\vartheta} \sqrt{\left(\dfrac{\partial R}{\partial \vartheta'}\right)^{2} + \left(\dfrac{\partial Z}{\partial \vartheta'}\right)^{2}} \,\mathrm{d}\vartheta' }{\displaystyle\int_0^{2{\rm \pi}} \sqrt{\left(\dfrac{\partial R}{\partial \vartheta'}\right)^{2} + \left(\dfrac{\partial Z}{\partial \vartheta'}\right)^{2}} \,\mathrm{d}\vartheta' }. \end{equation}

This choice of poloidal angle has been shown to correspond to the minimum of a spectral width functional (Hirshman & Breslau Reference Hirshman and Breslau1998). The Hirshman–Breslau technique can be used to obtain improved spectral convergence over the arclength angle by minimizing an energy functional with a nonlinear conjugate gradient method. In contrast, we present a method for improving the spectral convergence by computing an arclength angle on a uniform offset surface. This simply requires a Fourier transform to compute the coefficients in the new representation rather than the solution of a nonlinear optimization problem. We define a uniform offset surface by

(4.7)\begin{equation} \boldsymbol{r}_{\text{arclength}}=\boldsymbol{r}+a_{\text{arclength}}\boldsymbol{\hat{n}}, \end{equation}

where $a_{\text {arclength}}$ is a constant offset distance, taken normally to the plasma surface $\boldsymbol {r}$. This parameterization is particularly important for constraining the coil winding surface as described in § 5.3. It is critical that the offset distance $a_{\text {arclength}}$ is not too large; otherwise, the offset surface will self-intersect. The maximum offset distance is constrained by the principal curvatures, $\kappa _1$ and $\kappa _2$, of $S_{\text {plasma}}$. Here we use the convention that $\kappa _{1,2}<0$ indicates concavity. The maximum curvature in the outwardly concave regions of the plasma surface ($\kappa _{1,2}<0$) determines the minimum outward radius of curvature, and thus the maximum outward ($a_{\text {arclength}}>0$) offset distance. Conversely, the maximum curvature in the outwardly convex regions of the plasma surface ($\kappa _{1,2}>0$) determines the maximum inward ($a_{\text {arclength}}<0$) offset distance (Farouki Reference Farouki1986). The constraint on the offset distance is, therefore, given by

(4.8)\begin{equation} -\frac{1}{\max\left\{\kappa_1, \kappa_2\right\}} < a_{\text{arclength}} < \frac{1}{|\min\left\{\kappa_1, \kappa_2\right\}|}. \end{equation}

Note that for any closed toroidal surface, the surface integral of the Gaussian curvature ($\int \kappa _1 \kappa _2 \, d^{2} a$) must vanish, indicating that the minimum of the principal curvatures must always be negative and the maximum of the principal curvatures must be positive.

In this case, rather than using the arclength on the plasma surface to define $\bar {\theta }$, the arclength is computed on a surface uniformly offset from the plasma boundary as defined in (4.7) with the restriction on outward offsets as in (4.8). In this way, the poloidal angle can be tuned to accommodate a particular offset of the coil surface from the plasma surface or to increase resolution in regions of large (positive or negative) curvature.

We define $\omega$ to be the difference between the arclength and arctangent angles:

(4.9)\begin{equation} \omega(\vartheta,\zeta)=\bar{\theta}(\vartheta,\zeta)-\vartheta. \end{equation}

The quantities $l$ and $\omega$ can now be expressed as Fourier series in $\bar {\theta }$:

(4.10)\begin{equation} \left. \begin{gathered} l(\bar{\theta},\zeta)=\sum_{m,n}l_{m,n}^{c} \cos(m\bar{\theta}-n\zeta),\\ \omega(\bar{\theta},\zeta)=\sum_{m,n}\omega_{m,n}^{s} \sin(m\bar{\theta}-n\zeta). \end{gathered} \right\} \end{equation}

Note that while we express $\omega$ in a Fourier series to describe the relationship between the arctangent angle $\vartheta$ and the arclength angle $\bar {\theta }$, it is only necessary to differentiate with respect to $l_{m,n}^{c}$ when constructing shape gradients. In this way, the poloidal angle is fixed as the plasma boundary is perturbed, and the corresponding null space is eliminated. The cylindrical coordinates can now be expressed in the new representation as

(4.11)\begin{equation} \left. \begin{gathered} R(\bar{\theta},\zeta)=R_0(\zeta)+l(\bar{\theta},\zeta)\cos(\bar{\theta}-\omega(\bar{\theta},\zeta)),\\ Z(\bar{\theta},\zeta)=Z_0(\zeta)+l(\bar{\theta},\zeta)\sin(\bar{\theta}-\omega(\bar{\theta},\zeta)). \end{gathered} \right\} \end{equation}

The arclength angle is advantageous over the arctangent angle because it allows for a more efficient representation of the plasma surface with a uniform sampling of the poloidal angle. This is particularly useful for representing plasma boundaries with regions of large curvature. In figure 1, we compare the distribution of poloidal grid points on the bean-shaped cross-section of the W7-X boundary, which features both concave and convex regions. The grid points in the plasma arclength angle are equally spaced along the plasma surface cross-section while the points in the arctangent angle accumulate on the inboard and outboard sides. Observe that in highly convex regions of the plasma surface, the normal vector is changing rapidly, causing points on the offset surface to diverge from each other. By fixing the points on the offset surface to be spaced by uniform arclength increments, the points on the plasma surface tend to accumulate in the highly convex regions. Extra resolution in these regions is beneficial for improved convergence of REGCOIL calculations as well as a more efficient Fourier series representation. Alternatively, for surfaces which feature highly concave regions, it may be more advantageous to use an inward offset when computing the arclength angle. In this case, the maximum inward offset is set by the maximum curvature in the outwardly convex regions as in (4.8).

To further evaluate these representations, we compare the convergence of the Fourier series for the W7-X surface. The double Fourier series representation (4.1) is compact in both $m$ and $n$ as shown in figure 2, as the surface was optimized with a truncated spectrum. (The surface is from a fixed-boundary equilibrium which predated the coil optimization.) For the other representations, the modes form a tilted band which is thin in the $n$ direction and decays slowly in the $m$ direction as shown in figure 3. With the arctangent representation (figure 3a) many modes must be retained in the Fourier series to accurately reconstruct the plasma surface. With the arclength angle defined on the plasma boundary (figure 3b), the magnitude of the Fourier modes decreases more rapidly in $m$. This decay is enhanced when the arclength angle is computed on a surface displaced from the plasma by 0.25 m (figure 3d). A similar uniformly offset surface is used in § 6.3 for computation of the shape gradients.

Figure 2. Magnitudes of (a) $R_{m,n}^{c}$ and (b) $Z_{m,n}^{s}$ for the W7-X boundary on a grid of mode numbers $m$ and $n/N_p$, where $N_P$ is the number of field periods.

Figure 3. (a) Magnitudes of $l_{m,n}^{c}$ for the arctangent poloidal angle on a grid of mode numbers $m$ and $n/N_P$. The dominant $l_{m,n}^{c}$ lie on a long, tilted band. Magnitudes of (b) $l_{m,n}^{c}$ and (c) $\omega _{m,n}^{s}$ for the arclength poloidal angle on the plasma surface. The amplitudes are more localized to the smaller mode numbers than with the arctangent angle, but still lie on a long, tilted band. Magnitudes of (d) $l_{m,n}^{c}$ and (e) $\omega _{m,n}^{s}$ for the arclength poloidal angle on a surface uniformly offset from the plasma surface by 0.25 m. The amplitudes are even more localized to the smaller mode numbers.

A spectrum which is condensed in one basis is, in general, not condensed in another given basis. We observe that the Fourier spectrum of the chosen surface is, by design, condensed in the VMEC representation (figure 2), but it is more diffuse for the various single Fourier representations (figure 3). To demonstrate that this is a consequence of changing basis in general rather than a weakness of the single Fourier representations specifically, we present a transformation from a condensed spectrum in a single Fourier representation to the VMEC representation in figure 4. Starting from the VMEC representation of the W7-X surface, we compute the single Fourier representation with the arclength poloidal angle on a surface uniformly offset from the plasma surface by 0.25 m. The original VMEC spectrum has a total of 43 non-zero modes while this single Fourier spectrum has 382 modes with magnitudes greater than $10^{-5}$, so we truncate this single Fourier series to only include the 43 modes with the largest magnitudes. Starting from the truncated single Fourier spectrum, we then compute the double Fourier representation in (4.1) with the initial VMEC poloidal angle. We find that the VMEC spectrum calculated in this way now has a total of 734 modes with magnitudes greater than $10^{-5}$. We also present figure 4(d) showing the magnitudes of the single and double Fourier spectra sorted in descending order. We observe that the single Fourier series decays faster than the double Fourier representation. Thus, we can conclude that the single Fourier representation is more ‘efficient’, as fewer modes are required to represent the surface.

Figure 4. Magnitudes of (a) $l_{m,n}^{c}$ for the arclength poloidal angle on a surface uniformly offset from the plasma surface by 0.25 m as in figure 3(d), but retaining only the 43 largest-amplitude modes (the number of non-zero modes in the VMEC representation). Magnitudes of (b) $R_{m,n}^{c}$ and (c) $Z_{m,n}^{s}$ obtained by starting with the truncated single Fourier spectrum in (a). A dispersion similar to that of figure 3(d) is observed for these spectra. (d) Magnitudes of the full set of $l_{m,n}^{c}$ (figure 3d) in blue and $R_{m,n}^{c}$ and $Z_{m,n}^{s}$ computed from the truncated single Fourier series (b,c) in red, sorted in descending order.

5. Differentiation with respect to the Fourier coefficients

To take the derivative of $\chi ^{2}_I$, where $I\in \{B,K\}$, with respect to some parameter of the plasma surface, one must consider the explicit and implicit dependence of the functional on the parameters. We call the set of parameters of the plasma surface $\boldsymbol {\varOmega }$ and a particular one $\varOmega _j$. The derivative of $\chi ^{2}_I$ with respect to $\varOmega _j$ is computed using the chain rule:

(5.1)\begin{equation} \frac{\partial \chi^{2}_I (\boldsymbol{\varOmega},\boldsymbol{\varPhi}(\boldsymbol{\varOmega}))}{\partial \varOmega_j}=\left.\frac{\partial \chi^{2}_I }{\partial \varOmega_j}\right|_{\mathrm{d}\boldsymbol{\varPhi}=0}+ \frac{\partial \chi^{2}_I}{\partial \boldsymbol{\varPhi}}\boldsymbol{\cdot} \frac{\partial \boldsymbol{\varPhi}}{\partial \varOmega_j}. \end{equation}

The bar notation here means that the first derivative on the right-hand side of (5.1) is taken with $\mathrm {d}\boldsymbol {\varPhi }=0$. The derivatives involving $\boldsymbol {\varPhi }$ indicate element-wise operation; thus, $\partial \chi ^{2}_I/\partial \boldsymbol {\varPhi }$ and $\partial \boldsymbol {\varPhi }/ \partial \varOmega _j$ are vectors of the same dimension as $\boldsymbol {\varPhi }$. The dot product between these vectors in the second term on the right-hand side of (5.1) indicates contraction over $\boldsymbol {\varPhi }$. The second term on the right-hand side of (5.1) accounts for the implicit dependence of $\chi ^{2}_I$ on $\varOmega _j$ through $\boldsymbol {\varPhi }$. In §§ 5.1 and 5.2, $\boldsymbol {\varPhi }$ will gain dependence on $\boldsymbol {\varOmega }$ due to constraints imposed on the total derivative.

5.1. Optimal solution constraint

To take this derivative such that the current potential satisfies the linear least-squares system given in (2.5), we differentiate the constraint given by (2.5) with respect to the Fourier coefficients to obtain

(5.2)\begin{equation} \frac{\partial \boldsymbol{\varPhi}}{\partial \varOmega_j} = {{{\mathsf{A}}}}^{{-}1}\left(\frac{\partial \boldsymbol{b}}{\partial \varOmega_j}-\frac{\partial {{{\mathsf{A}}}}}{\partial \varOmega_j}\boldsymbol{\varPhi}\right), \end{equation}

which can then be substituted into (5.1). Ordinarily, determining $\partial \boldsymbol {\varPhi } / \partial \varOmega _j$ would require solving a linear system involving ${\mathsf{A}}$ for each $\varOmega _j$; however, this can be avoided here by using an adjoint method. The second term on the right-hand side of (5.1) is

(5.3)\begin{equation} \frac{\partial \chi^{2}_I}{\partial \boldsymbol{\varPhi}}\boldsymbol{\cdot} {{{\mathsf{A}}}}^{{-}1}\left[\frac{\partial \boldsymbol{b}}{\partial \varOmega_j}-\frac{\partial {{{\mathsf{A}}}}}{\partial \varOmega_j}\boldsymbol{\varPhi}\right]=\left(\frac{\partial \chi^{2}_I}{\partial \boldsymbol{\varPhi}}\right)^{\textrm{T}} \left({{{\mathsf{A}}}}^{{-}1}\left[\frac{\partial \boldsymbol{b}}{\partial \varOmega_j}-\frac{\partial {{{\mathsf{A}}}}}{\partial \varOmega_j}\boldsymbol{\varPhi}\right]\right), \end{equation}

where the second expression uses the transpose to perform the inner product between the vectors of the first expression as a matrix multiplication. The quantity on the left-hand side of this product is a row vector while the quantity in parentheses on the right is a column vector. Note that

(5.4)\begin{equation} \left(\frac{\partial \chi^{2}_I}{\partial \boldsymbol{\varPhi}}\right)^{\textrm{T}} {{{\mathsf{A}}}}^{{-}1} = \left[({{{\mathsf{A}}}}^{{-}1})^{\textrm{T}} \frac{\partial \chi^{2}_I}{\partial \boldsymbol{\varPhi}}\right]^{\textrm{T}}, \end{equation}

meaning that the operator ${\mathsf{A}}^{-1}$ can be moved to the left-hand side of the inner product by taking its transpose. By defining an adjoint variable $\boldsymbol {q}_I$ through

(5.5)\begin{equation} {{{\mathsf{A}}}}^{\textrm{T}} \boldsymbol{q}_I=\frac{\partial \chi^{2}_I}{\partial \boldsymbol{\varPhi}}, \end{equation}

the inner product term (5.3) becomes

(5.6)\begin{equation} \boldsymbol{q}_I^{\textrm{T}} \left(\frac{\partial \boldsymbol{b}}{\partial \varOmega_j}-\frac{\partial {{{\mathsf{A}}}}}{\partial \varOmega_j}\boldsymbol{\varPhi}\right) = \boldsymbol{q}_I \boldsymbol{\cdot} \left(\frac{\partial \boldsymbol{b}}{\partial \varOmega_j}-\frac{\partial {{{\mathsf{A}}}}}{\partial \varOmega_j}\boldsymbol{\varPhi}\right); \end{equation}

thus, a linear system involving the matrix ${\mathsf{A}}$ only needs to be solved once for all $\varOmega _j$ to get the adjoint variable since $\partial \chi ^{2}_I / \partial \boldsymbol {\varPhi }$ is not specific to a particular $\varOmega _j$. The dependence of $\chi ^{2}_I$ on $\boldsymbol {\varPhi }$ is given by the corresponding matrices and vectors:

(5.7)\begin{equation} \frac{\partial \chi^{2}_I}{\partial \boldsymbol{\varPhi}} = 2 N_p ( {{{\mathsf{A}}}}^{I} \boldsymbol{\varPhi} - \boldsymbol{b}^{I}), \end{equation}

where the definitions from appendix A of Landreman (Reference Landreman2017) have been used. Then (5.1) becomes

(5.8)\begin{equation} \left.\frac{\partial \chi^{2}_I(\boldsymbol{\varOmega},\boldsymbol{\varPhi}(\boldsymbol{\varOmega}))}{\partial \varOmega_j}\right|_{{{{\mathsf{A}}}}\boldsymbol{\varPhi}=\boldsymbol{b}}=\left.\frac{\partial \chi^{2}_I }{\partial \varOmega_j}\right|_{\mathrm{d}\boldsymbol{\varPhi}=0}+\boldsymbol{q}_I\boldsymbol{\cdot}\left(\frac{\partial \boldsymbol{b}}{\partial \varOmega_j}-\frac{\partial {{{\mathsf{A}}}}}{\partial \varOmega_j}\boldsymbol{\varPhi}\right). \end{equation}

This equation applies to both $\chi ^{2}_B$ and $\chi ^{2}_K$, so both adjoint variables $\boldsymbol {q}_B$ and $\boldsymbol {q}_K$ are computed from (5.5). Since the quantity $\chi ^{2}_B + \lambda \chi ^{2}_K$ is minimized with respect to $\boldsymbol {\varPhi }$ by REGCOIL,

(5.9)\begin{equation} \frac{\partial \chi_B^{2}}{\partial \boldsymbol{\varPhi}} + \lambda \frac{\partial \chi_K^{2}}{\partial \boldsymbol{\varPhi}} = 0. \end{equation}

This means that these adjoint variables are related by the equation

(5.10)\begin{equation} \boldsymbol{q}_B + \lambda \boldsymbol{q}_K = 0, \end{equation}

so it is only necessary to determine one of them by solving a linear system.

5.2. Fixed norm constraint

In the previous section, differentiation was performed at fixed $\lambda$. In order to evaluate different magnetic configurations on the same footing, we use the freedom in $\lambda$ to fix a quantity such as $\|K\|_2$ (2.4) which is related to the engineering complexity of the coils. To account for this constraint, we define the following quantity:

(5.11)\begin{equation} G=\|K\|_2-\|K\|^{\text{target}}_2, \end{equation}

where $\|K\|^{\text {target}}_2$ is the fixed value desired for $\|K\|_2$. In effect, this constraint allows the coil complexity to be fixed to an acceptable level. Let the optimal solution constraint be described by

(5.12)\begin{equation} \boldsymbol{F}={{{\mathsf{A}}}}\boldsymbol{\varPhi} - \boldsymbol{b}. \end{equation}

Now $\lambda$ must be allowed to vary unlike in the previous section. By taking the total differentials of $\boldsymbol {F}$ and $G$, we obtain

(5.13)\begin{equation} \left. \begin{gathered} \mathrm{d}\boldsymbol{F} = \sum_j \left( \frac{\partial {{{\mathsf{A}}}}}{\partial \varOmega_j}\boldsymbol{\varPhi} - \frac{\partial \boldsymbol{b}}{\partial \varOmega_j} \right)\mathrm{d}\varOmega_j + {{{\mathsf{A}}}} \mathrm{d}\boldsymbol{\varPhi} +( {{{\mathsf{A}}}}^{K} \boldsymbol{\varPhi} - \boldsymbol{b}^{K} )\,\mathrm{d}\lambda = 0,\\ \mathrm{d}G = \left.\sum_j \frac{\partial G}{\partial \varOmega_j}\right|_{\mathrm{d}\boldsymbol{\varPhi}=0} \mathrm{d}\varOmega_j + \frac{\partial G}{\partial \boldsymbol{\varPhi}} \boldsymbol{\cdot} \mathrm{d}\boldsymbol{\varPhi} = 0. \end{gathered} \right\} \end{equation}

The differential of $\lambda$ can be eliminated from these equations, as described in § A.1, and $\partial \boldsymbol {\varPhi } / \partial \varOmega _j$ subject to the constraints $\boldsymbol {F}=0$ and $G=0$ is obtained:

(5.14)\begin{align} \left.\frac{\partial \boldsymbol{\varPhi} }{\partial \varOmega_j}\right|_{\boldsymbol{F}=0,\ G=0} &={-}{{{\mathsf{A}}}}^{{-}1}\left[ \left( \frac{\partial {{{\mathsf{A}}}}}{\partial \varOmega_j}\boldsymbol{\varPhi} - \frac{\partial \boldsymbol{b}}{\partial \varOmega_j} \right)\right.\nonumber\\ &\quad \left.+\, \frac{({{{\mathsf{A}}}}^{K}\boldsymbol{\varPhi} - \boldsymbol{b}^{K})}{\tilde{\boldsymbol{q}}\boldsymbol{\cdot}({{{\mathsf{A}}}}^{K}\boldsymbol{\varPhi} - \boldsymbol{b}^{K})} \left[\left. \frac{\partial G}{\partial \varOmega_j}\right|_{\mathrm{d}\boldsymbol{\varPhi}=0} - \tilde{\boldsymbol{q}}\boldsymbol{\cdot} \left( \frac{\partial {{{\mathsf{A}}}}}{\partial \varOmega_j}\boldsymbol{\varPhi} - \frac{\partial \boldsymbol{b}}{\partial \varOmega_j} \right) \right] \right], \end{align}

where an additional adjoint variable $\tilde {\boldsymbol {q}}$ is chosen to satisfy

(5.15)\begin{equation} {{{\mathsf{A}}}}^{\textrm{T}} \tilde{\boldsymbol{q}}=\frac{\partial G}{\partial \boldsymbol{\varPhi}}. \end{equation}

For the case of $G$ given by (5.11), the explicit derivative $\partial G / \partial \varOmega _j = 0$, and $\partial G /\partial \boldsymbol {\varPhi }$ is proportional to $\partial \chi _K^{2} /\partial \boldsymbol {\varPhi }$. These conditions cause the implicit dependence in (5.1) to cancel, leaving only the explicit dependence:

(5.16)\begin{equation} \left.\frac{\partial \chi^{2}_I(\boldsymbol{\varOmega},\boldsymbol{\varPhi}(\boldsymbol{\varOmega}))}{\partial \varOmega_j}\right|_{ {{{\mathsf{A}}}}\boldsymbol{\varPhi}=\boldsymbol{b}, \ \|K\|_2=\|K\|^{\text{target}}_2 }=\left.\frac{\partial \chi^{2}_I }{\partial \varOmega_j}\right|_{\mathrm{d}\boldsymbol{\varPhi}=0}, \end{equation}

as discussed in § A.3.

5.3. Uniform offset constraint

Up to this point, the derivatives have been computed with a fixed winding surface. Without placing an additional constraint on the winding surface, the plasma boundary may move uniformly towards the winding surface in order to minimize the coil–plasma distance, as the shaping components of the magnetic field decay with distance from the coils (Boozer Reference Boozer2000; Landreman & Boozer Reference Landreman and Boozer2016). For this reason, we instead compute derivatives with respect to the plasma surface parameters while maintaining the constraint that the coil winding surface is uniformly offset a distance $a$ from the plasma surface. The coil surface is then defined by

(5.17)\begin{equation} \boldsymbol{r}'=\boldsymbol{r}+a\boldsymbol{\hat{n}}. \end{equation}

The coil offset distance must similarly be constrained by the principal curvatures as

(5.18)\begin{equation} 0 < a < \frac{1}{|\min\left\{\kappa_1, \kappa_2\right\}|}. \end{equation}

Unlike the restriction on the offset surface that defines the poloidal angle (4.8), the coil offset distance is constrained by (5.18) to be greater than zero since the coils must lie in the vacuum region. We remark that as the plasma surface evolves during the optimization, a constraint must be placed on the principal curvatures such that the above inequality is satisfied. As concavity of the plasma boundary has been shown to be correlated with coil complexity (Paul et al. Reference Paul, Landreman, Bader and Dorland2018), placing such a constraint on the curvature of the boundary is reasonable and can be performed with a penalty formulation (Paul, Landreman & Antonsen Reference Paul, Landreman and Antonsen2020b). If desired, the winding surface could be offset by a non-uniform function on the surface, and instead, constraints would be imposed on the minimum and maximum offset distance.

This constraint couples the geometry of the winding surface to that of the plasma surface and, therefore, introduces additional derivatives with respect to the plasma surface Fourier amplitudes. For example, ${\mathsf{A}}$ gains additional dependence through the ${\mathsf{A}}^{K}$ term (2.6). Some equations are given in § A.2 showing how the derivatives of the coil winding surface position and normal vectors become coupled to the plasma surface through (5.17).

No new adjoint variables are needed to impose this constraint. With the combination of this and the fixed $\|K\|_2$ constraint, $G$ given by (5.11) now has explicit dependence on the plasma surface parameters; however, $\partial G /\partial \boldsymbol {\varPhi }$ is still proportional to $\partial \chi _K^{2} /\partial \boldsymbol {\varPhi }$. The total derivatives of $\chi ^{2}_B$ and $\chi ^{2}_K$ are then given by

(5.19)\begin{align} \left. \begin{aligned} \left.\frac{\partial \chi^{2}_B(\boldsymbol{\varOmega},\boldsymbol{\varPhi}(\boldsymbol{\varOmega}))}{\partial \varOmega_j}\right|_{ {{{\mathsf{A}}}}\boldsymbol{\varPhi}=\boldsymbol{b}, \ \|K\|_2=\|K\|^{\text{target}}_2, \ \boldsymbol{r}'=\boldsymbol{r}+a\boldsymbol{\hat{n}} }&=\left.\frac{\partial \chi^{2}_B }{\partial \varOmega_j}\right|_{\mathrm{d}\boldsymbol{\varPhi}=0} -\lambda \left( \frac{\chi^{2}_K}{a_{\text{coil}}}\frac{\partial a_{\text{coil}}}{\partial \varOmega_j} - \left.\frac{\partial \chi^{2}_K }{\partial \varOmega_j}\right|_{\mathrm{d}\boldsymbol{\varPhi}=0} \right),\\ \left. \frac{\partial \chi^{2}_K(\boldsymbol{\varOmega},\boldsymbol{\varPhi}(\boldsymbol{\varOmega}))}{\partial \varOmega_j}\right|_{ {{{\mathsf{A}}}}\boldsymbol{\varPhi}=\boldsymbol{b}, \ \|K\|_2=\|K\|^{\text{target}}_2, \ \boldsymbol{r}'=\boldsymbol{r}+a\boldsymbol{\hat{n}} }&=\frac{\chi^{2}_K}{a_{\text{coil}}}\frac{\partial a_{\text{coil}}}{\partial \varOmega_j}. \end{aligned} \right\} \end{align}

A detailed derivation of these equations is given in § A.3. Note that the implicit dependence of $\chi ^{2}_B$ is $-\lambda$ times the implicit dependence of $\chi ^{2}_K$, and the total dependence of $\chi ^{2}_K$ is such that $\|K\|_2$ vanishes.

Since higher-order derivatives of the position vector are needed to compute the normal vector and its derivatives, any high-frequency noise in the plasma surface representation is amplified when taking derivatives with the constraint given by (5.17). For this reason, the poloidal angle defined with respect to the arclength on the winding surface as described in § 4 greatly improves the numerical performance with this constraint.

6. Benchmarks and demonstrations

In this section, we present several benchmark calculations to verify the derivatives and constraints from the previous section.

6.1. Finite-difference derivatives

The derivatives of $\chi ^{2}_I$ are implemented analytically as described in the previous section rather than with a finite-difference approximation. The former is advantageous over the latter in that it takes far less computation time for a high-dimensional derivative; however, the finite-difference result should converge to the analytic result as the step size is reduced. This is used to check the accuracy of the analytic derivatives. Consider a numerical forward-difference approximation of $\partial \chi ^{2}_I / \partial \varOmega _j$:

(6.1)\begin{equation} \left(\frac{\partial \chi^{2}_{I}}{\partial \varOmega_j}\right)_{\text{num.}}=\frac{\chi^{2}_I(\varOmega_j+\varDelta\varOmega)-\chi^{2}_I(\varOmega_j)}{\varDelta\varOmega}. \end{equation}

When doing this calculation, the perturbation of the current potential corresponding to the perturbation of $\varOmega _j$ is included in order to account for the implicit dependence of $\chi ^{2}_I$ on the plasma surface. This expression is a first-order approximation, so it should agree with the analytic expression to ${O}(\varDelta \varOmega )$.

We now define the fractional difference between the analytic and numerical derivatives as

(6.2)\begin{equation} \text{fractional difference}=\left(\frac{\partial \chi^{2}_I}{\partial \varOmega_j} - \left(\frac{\partial \chi^{2}_{I}}{\partial \varOmega_j}\right)_{\text{num.}}\right) \left(\frac{\partial \chi^{2}_I}{\partial \varOmega_j}\right)^{{-}1}. \end{equation}

The following results were generated using the QHS46 equilibrium (Shohet et al. Reference Shohet, Anderson, Anderson and Talmadge1991) with a $0.25\ \text {m}$ uniformly offset coil winding surface which was held fixed during the calculation of derivatives with respect to the plasma parameters. A plot of the fractional difference for the derivatives of $\|B\|_2$ with only the optimal solution constraint and uniform offset constraint imposed is given in figure 5(a). The surface is described by the single Fourier representation given in (4.10) with the arclength angle defined on a surface chosen to have the same $0.25\ \text {m}$ uniform offset as the coil winding surface. A similar fractional difference plot is given for $\|K\|_2$ in figure 5(b). We find that the fractional difference scales linearly with $\varDelta \varOmega$ for $10^{-7} \lesssim \varDelta \varOmega \lesssim 10^{-2}$. For $\varDelta \varOmega \lesssim 10^{-7}$, the round-off error begins to dominate, and the linear relationship is no longer observed (Sauer Reference Sauer2012). Similar trends are also observed in the fractional difference for $\|B\|_2$ with the addition of the fixed $\|K\|_2$ constraint, presented in figure 5(c). For this calculation with the fixed $\|K\|_2$ constraint, $\|K\|_2^{\text {target}}=1.6\ \text {MA}\,\text {m}^{-1}$ was used.

Figure 5. Modes with $m$ up to 2 and $n$ up to $\pm 2 N_p$ are displayed with lines on the red end of the spectrum corresponding to larger magnitudes of the derivative and lines on the blue end of the spectrum corresponding to smaller magnitudes. We present the fractional difference between analytic and numerical calculations of (a) $\partial \|B_n\|_2 / \partial l^{c}_{m,n}$, (b) $\partial \|K\|_2 / \partial l^{c}_{m,n}$ and (c) $\partial \|B_n\| / \partial l^{c}_{m,n}$ while holding $\|K\|_2$ fixed. For each of these, the offset distance is additionally fixed.

6.2. Area shape gradient benchmark

As a benchmark for calculating shape gradients with the modified REGCOIL code, the shape gradient of the plasma surface area is computed. There is an analytic expression for this shape gradient (Landreman & Paul Reference Landreman and Paul2018), which is given by

(6.3)\begin{equation} S_{a_{\text{plasma}}}=2H, \end{equation}

where $H = (\kappa _1 + \kappa _2)/2$ is the mean curvature of the plasma surface. Again, we assume the convention that positive $H$ indicates convexity. The shape gradient of the plasma area calculated using (3.4) is given in figure 6(a). The calculation is performed for the W7-X boundary. (This is the boundary of a fixed-boundary equilibria that preceded the coil design and does not include coil ripple.) The shape gradient is computed on a grid of $N_{\bar {\theta }}=200$ grid points in the poloidal angle and $N_{\zeta } = 200$ grid points in the toroidal angle. The derivatives of the area are computed for $m \le 35$ and $|n/N_P| \le 35$, and the shape gradient is discretized with a Fourier series with the same set of modes, so (3.4) is a square linear system. We compare this to the expected result, given in figure 6(b). The average error in the shape gradient,

(6.4)\begin{equation} \text{error} = \frac{\int_{S_{\text{plasma}}} \, |S_{a_{\text{plasma}}} - 2 H| \, \mathrm{d}a }{\int_{S_{\text{plasma}}} \, |2 H| \, \mathrm{d}a}, \end{equation}

is computed to be $3.69 \times 10^{-4}$. In figure 6(d), we display the convergence of the error with respect to the modes retained in the derivatives and in computing the shape gradient. We see that the expected and computed shape gradients converge to each other with increasing numerical resolution.

Figure 6. Plot in the $\vartheta$ and $\zeta$ plane of (a) the shape gradient of the plasma area calculated in the REGCOIL code using a solution of (3.4) for a square system with $m \le 35$ and $|n/N_P| \le 35$ and (b) the shape gradient of the plasma area calculated using the explicit form involving the mean curvature (6.3). In (c), the same quantity plotted in (a) is displayed on the W7-X plasma boundary, and in (d) the error defined in (6.4) is shown as a function of the maximum poloidal and toroidal mode number used to compute the shape gradient.

6.3. Shape gradients of $\|B_n\|_2$ and $\|K\|_2$

In this section, (3.4) is used to obtain the shape gradients of $\|B_n\|_2$ and $\|K\|_2$ for the QHS46 and W7-X configurations using the derivatives obtained from the adjoint method. For both configurations, we choose the winding surface to be uniformly offset from the plasma boundary (5.17). In the case of W7-X we choose $a = 0.37$ m to match the minimum coil–plasma distance of the actual W7-X winding surface. For QSH46 we choose $a = 0.21$ m such that the offset distance is scaled according to the minor radius of the plasma boundary. For the QHS46 calculations, a fixed value of $\|K\|_2^{\text {target}}=1.7\ \text {MA}\,\text {m}^{-1}$ was used, while for W7-X it was taken to be ${2.4}\ \text {MA}\,\text {m}^{-1}$. The computation of the derivatives necessary for the evaluation of these shape gradients required approximately 7 h on a single AMD Opteron 6136 CPU. The cost of construction of the shape gradient from the parameter derivatives is negligible in comparison (${<}1\ \textrm {min}$ on a laptop).

In figure 7, results are shown for the three-dimensional shape of the QHS46 plasma surface. The parameter derivatives with respect to modes with $m \le 30$ and $|n/N_P| \le 30$ are computed on a grid with $N_{\bar {\theta }} = N_{\zeta } = 100$, and the shape gradient is discretized with a Fourier series with the same set of modes. The shape gradient of $\|B_n\|_2$ for only the optimal solution constraint is shown in figure 7(a) and similarly for $\|K\|_2$ in figure 7(b). We see that for both of these quantities, the shape gradient is large and positive in the regions of convex curvature and strongly negative in the nearby regions where the surface curvature is locally reduced. This indicates that in order to reduce both of these figures of merit, the sharp ridge features on the boundary must be locally rounded. Similar features are present when the shape gradient of $\|B_n\|_2$ is computed with the fixed $\|K\|_2$ constraint as shown in figure 7(c), though the magnitude is slightly increased. In figure 7(d), the shape gradient with the additional uniform offset constraint is qualitatively similar, but an additional feature is present near the triangle-shaped cross-section.

Figure 7. The shape gradient of (a) $\|B_n\|_2$, (b) $\|K\|_2$, (c) $\|B_n\|_2$ with fixed $\|K\|_2$ and (d) $\|B_n\|_2$ with fixed $\|K\|_2$ and fixed offset distance is computed for the QHS46 boundary with a uniform offset winding surface of 0.21 m. The shape gradient is computed using a linear solve of (3.4) for a square system with modes $m \le 30$ and $|n/N_P| \le 30$.

In figure 8, we present the results for the W7-X boundary. The parameter derivatives with respect to modes with $m \le 30$ and $|n/N_P| \le 30$ are computed on a grid with $N_{\bar {\theta }} = N_{\zeta } = 100$, and the shape gradient is discretized with a Fourier series with the same set of modes. We see that the shape gradient for $\|B_n\|_2$ is again peaked in the regions of large convexity and slightly negative in the regions of concavity. This indicates that the normal field could be reduced by pushing the surface inward in the sharply convex regions and outward in the concave regions. We see that the shape gradient for $\|K\|_2$ similarly has a large magnitude in the convex regions. Considering the shape gradient of $\|B_n\|_2$ with the additional constraints, we see that the trends are qualitatively similar.

Figure 8. The shape gradient of (a) $\|B_n\|_2$, (b) $\|K\|_2$, (c) $\|B_n\|_2$ with fixed $\|K\|_2$ and (d) $\|B_n\|_2$ with fixed $\|K\|_2$ and fixed offset distance is computed for the W7-X boundary with a uniform offset winding surface of 0.37 m. The shape gradient is computed using a linear solve of (3.4) for a square system with modes $m \le 30$ and $|n/N_P| \le 30$.

Finally, we mention that the calculations presented in this section rely on the choice of the offset distance, $a$. We find that the trends presented here do not strongly vary with this choice. This is consistent with the observations in Paul et al. (Reference Paul, Landreman, Bader and Dorland2018).