2.1. Managing code failures

The equilibrium code may return with an error for a given plasma boundary. There are many reasons for these failures, such as the flux coordinate Jacobian becoming ill-conditioned or the number of iterations exceeding the maximum prescribed value. When such a failure is experienced during the optimization, the objective function is set to an arbitrarily large value (e.g. $10^{12}$) to enforce an effective constraint. The approach of assigning a very large fictitious objective value at unevaluable points is common in the optimization literature (Rasheed, Hirsh & Gelsey Reference Rasheed, Hirsh and Gelsey1997; Emmerich et al. Reference Emmerich, Giotis, Özdemir, Bäck and Giannakoglou2002) and is employed in the STELLOPT code (Lazerson et al. Reference Lazerson, Schmitt and Zhu2020).

As encountering such unevaluable points makes the parameter space non-smooth, it is prudent to try to avoid code failures by placing additional constraints on the parameter space. Such a technique has been employed in the optimization of aircraft by checking that simulation outputs match physical model assumptions, such as the drag coefficient being non-negative (Gelsey Reference Gelsey1995; Gelsey, Schwabacher & Smith Reference Gelsey, Schwabacher and Smith1998). In § 3, we present constraints which prevent the surface from self-intersecting or from obtaining large curvature, leading to an ill-conditioned Jacobian and code failure. Even in the presence of such constraints, unevaluable points may still be present. Here the objective function becomes discontinuous, which is problematic for optimizers which assume function continuity. In this work we employ the Broyden–Fletcher–Goldfarb–Shanno (BFGS) quasi-Newton method with an Armijo–Wolfe line search, which assumes $C^2$ continuity of the objective function (Nocedal & Wright Reference Nocedal and Wright2006). Although convergence of the BFGS algorithm is not guaranteed for non-smooth problems, it has been observed that a reasonable approximation to the optimum is most often achieved (Lemaréchal Reference Lemaréchal1982). With the standard Armijo–Wolfe line search method, the gradient need not be evaluated unless the objective function satisfies the sufficient decrease criterion. Thus the gradient does not need to be evaluated at unevaluable points. In this way, as long as the line search never returns a point where the objective is not differentiable, the BFGS method is well defined. Although there are specialized quasi-Newton methods for non-smooth objective functions (Lewis & Overton Reference Lewis and Overton2013), we have obtained acceptable results with a standard BFGS method.

3.1. Preventing surface self-intersection

We now describe a constraint which prevents self-intersection of the plasma boundary. Given a surface described by the cylindrical coordinates $R(\theta ,\phi )$ and $Z(\theta ,\phi )$, self-intersection of the boundary may occur if either:

1. (a) $R(\theta ,\phi )< 0$ at any point; or

2. (b) the planar curve $\boldsymbol {x}_{\phi _0}(\theta ) = \{ R(\theta ,\phi _0), Z(\theta ,\phi _0)\}$ is self-intersecting for any $\phi _0$.

Condition $(a)$ can be avoided using a penalty objective of the form,

(3.1)\begin{equation} f_R(S_P) = \frac{1}{A_P} \int_{S_P} \,\textrm{d}^2 x \exp \left(- \frac{\left(R - R_\textrm{min}\right)}{w_R} \right), \end{equation}

where $w_R$ is a weight which sets the gradient length scale of the objective, $R_\textrm {min}$ is the minimum allowable major radius and $A_P$ is the area of $S_P$.

Condition $(b)$ can be avoided by introducing a constraint on the global radius of curvature (Gonzalez & Maddocks Reference Gonzalez and Maddocks1999; Walker Reference Walker2016) of each of the curves, $\boldsymbol {x}_{\phi _0}(\theta )$, which will now be defined. The radius of the unique circumcircle containing any three points $\boldsymbol {x}_1$, $\boldsymbol {x}_2$, and $\boldsymbol {x}_3$ (figure 1), can be computed from

(3.2)\begin{equation} \rho(\boldsymbol{x}_1,\boldsymbol{x}_2,\boldsymbol{x}_3) = \frac{|\boldsymbol{x}_1 - \boldsymbol{x}_2\|\boldsymbol{x}_1 - \boldsymbol{x}_3\|\boldsymbol{x}_2-\boldsymbol{x}_3|}{4 \mathcal{A}(\boldsymbol{x}_1,\boldsymbol{x}_2,\boldsymbol{x}_3)}, \end{equation}

where $\mathcal {A}(\boldsymbol {x}_1,\boldsymbol {x}_2,\boldsymbol {x}_3)$ is the area of a triangle with vertices $\boldsymbol {x}_1$, $\boldsymbol {x}_2$, and $\boldsymbol {x}_3$. Under the assumption that these points lie on a non-self-intersecting smooth curve, $\mathcal {C}$ such that $\boldsymbol {x}_1 = \boldsymbol {x}(l_1)$, $\boldsymbol {x}_2 = \boldsymbol {x}(l_2)$, and $\boldsymbol {x}_3 = \boldsymbol {x}(l_3)$, then the radius of the circumcircle satisfies $\rho (\boldsymbol {x}(l_1),\boldsymbol {x}(l_2),\boldsymbol {x}(l_2)) \le \rho (\boldsymbol {x}(l_1),\boldsymbol {x}(l_2),\boldsymbol {x}(l_3))$. This limiting case can be computed from

(3.3)\begin{equation} \rho(\boldsymbol{x}(l_1),\boldsymbol{x}(l_2),\boldsymbol{x}(l_2)) = \lim_{l_2 \rightarrow l_3} \rho(\boldsymbol{x}(l_1),\boldsymbol{x}(l_2),\boldsymbol{x}(l_3)) = \frac{|\boldsymbol{x}(l_1)-\boldsymbol{x}(l_2)|}{2\sqrt{1 - \left(\hat{\boldsymbol{t}}(\boldsymbol{x}(l_2)) \boldsymbol{\cdot} \dfrac{\boldsymbol{x}(l_1)-\boldsymbol{x}(l_2)}{|\boldsymbol{x}(l_1)-\boldsymbol{x}(l_2)|} \right)^2}}, \end{equation}

where $\hat {\boldsymbol {t}}$ is the unit tangent. This is the radius of the circle that passes through $\boldsymbol {x}(l_1)$ and is tangent to $\mathcal {C}$ at $\boldsymbol {x}(l_2)$, which we define as the self-contact function between these points,

(3.4)\begin{equation} \mathcal{S}_{\mathcal{C}}(\boldsymbol{x}_1,\boldsymbol{x}_2) = \frac{|\boldsymbol{x}_1 - \boldsymbol{x}_2|}{2 \sqrt{1-\left(\hat{\boldsymbol{t}}(\boldsymbol{x}_2) \boldsymbol{\cdot} \dfrac{\left(\boldsymbol{x}_1-\boldsymbol{x}_2 \right)}{|\boldsymbol{x}_1 - \boldsymbol{x}_2|}\right)^2}}, \end{equation}

the radius of the so-called point-tangent circle (figure 1).

Figure 1.

(a) The circumcircle containing $\boldsymbol {x}_1$, $\boldsymbol {x}_2$ and $\boldsymbol {x}_3$ with radius $\rho (\boldsymbol {x}_1,\boldsymbol {x}_2,\boldsymbol {x}_3)$. (b) The points $\boldsymbol {x}_1$, $\boldsymbol {x}_2$, and $\boldsymbol {x}_3$ lie on a closed elliptical curve $\mathcal {C}$ (blue). The point tangent radius for points $\boldsymbol {x}_1$ and $\boldsymbol {x}_2$ and points $\boldsymbol {x}_1$ and $\boldsymbol {x}_3$ are displayed. The global radius of curvature at $\boldsymbol {x}_1$ is $\mathcal {S}_{\mathcal {C}}(\boldsymbol {x}_1,\boldsymbol {x}_2)$.

In the limit that $l_2, l_3 \rightarrow l_1$, we obtain

(3.5)\begin{equation} \lim_{l_2,l_3 \rightarrow l_1} \rho(\boldsymbol{x}(l_1),\boldsymbol{x}(l_2),\boldsymbol{x}(l_3)) = \rho_{\mathcal{C}}(\boldsymbol{x}(l_1)), \end{equation}

where $\rho _{\mathcal {C}}(\boldsymbol {x}(l_1))$ is the local radius of curvature of $\mathcal {C}$ at $\boldsymbol {x}(l_1)$. We define the global radius of curvature as

(3.6)\begin{equation} \rho_{G}(\boldsymbol{x}_1) = \min_{\boldsymbol{x}_2,\boldsymbol{x}_3 \in \mathcal{C}} \rho\left(\boldsymbol{x}_1,\boldsymbol{x}_2,\boldsymbol{x}_3\right) = \min_{\boldsymbol{x}_2 \in \mathcal{C}} \mathcal{S}_{\mathcal{C}}. \end{equation}

This can be thought of as a generalization of the local radius of curvature, as we have the inequality $0 \le \rho _G(\boldsymbol {x}) \le \rho _{\mathcal {C}}(\boldsymbol {x})$. At a given point, the minimizer in (3.6) will either be a point where $(\boldsymbol {x}_1 - \boldsymbol {x}_2) \boldsymbol {\cdot } \hat {\boldsymbol {t}}(\boldsymbol {x}_2) = 0$ or $\rho _{\mathcal {C}}(\boldsymbol {x}_2) = \mathcal {S}_{\mathcal {C}}(\boldsymbol {x}_1,\boldsymbol {x}_2)$ (Smutny Reference Smutny2004). In other words, at a given point $\boldsymbol {x}_1$, the point on the curve that has the smallest self-contact function containing $\boldsymbol {x}_1$ will be where the local radius of curvature is equal to the point-tangent radius or the displacement of the points $(\boldsymbol {x}_1 - \boldsymbol {x}_2)$ is orthogonal to the curve at $\boldsymbol {x}_2$.

If the inequality

(3.7)\begin{equation} \min_{\boldsymbol{x}_1 \in \mathcal{C}} \rho_G(\boldsymbol{x}_1) = \min_{\boldsymbol{x}_1,\boldsymbol{x}_2 \in \mathcal{C}} \mathcal{S}_{\mathcal{C}}(\boldsymbol{x}_1,\boldsymbol{x}_2) \ge \mathcal{R} \end{equation}

is satisfied, this implies that the curve can be ‘thickened’ to a tube without self-contact surrounding $\mathcal {C}$ such that at any given point, the cross-section of the tube is a circle in the plane perpendicular to the local tangent vector of radius $\mathcal {R}$ (Gonzalez & Maddocks Reference Gonzalez and Maddocks1999). For this reason, the concept of the global radius of curvature has been employed for the shape optimization of finite-thickness knots (Gonzalez & Maddocks Reference Gonzalez and Maddocks1999; Carlen et al. Reference Carlen, Laurie, Maddocks and Smutny2005; Walker Reference Walker2016).

We enforce the inequality constraint (3.7) to prevent the self-intersection of every planar curve $\boldsymbol {x}_{\phi _0}(\theta )$ with a penalty function of the form

(3.8)\begin{equation} f_{\mathcal{S}} = \frac{1}{A_P}\int_0^{2{\rm \pi}} \,\textrm{d} \phi \int_0^{2{\rm \pi}} \,\textrm{d}\theta \left|\displaystyle\frac{\partial \boldsymbol{r}}{\partial \theta} \times \displaystyle\frac{\partial \boldsymbol{r}}{\partial \phi}\right|\int_0^{2{\rm \pi}} \,\textrm{d} \theta' \exp\left(-\frac{\left(\mathcal{S}_{\mathcal{C}}\left(\boldsymbol{x}_{\phi}(\theta),\boldsymbol{x}_{\phi}(\theta')\right) - \mathcal{R}\right)}{w_{\mathcal{S}}}\right), \end{equation}

where $\mathcal {R}$ is the minimum allowable global radius of curvature and $w_{\mathcal {S}}$ is a weight function. In addition to preventing self-intersection of the boundary, this objective improves the regularity of the surface by reducing the curvature of the poloidal cross-sections. A demonstration of this penalty function is presented in § 3.3.

3.2. Surface curvature constraint

While constraining the self-contact function reduces the curvature of the cross-sections of the plasma boundary, we include additional terms which improve the regularity of the plasma boundary. First, we place an effective constraint on the magnitude of the principal curvatures,

(3.9a,b)\begin{equation} |\kappa_1| \le \kappa_{\textrm{max},1} \quad |\kappa_2| \le \kappa_{\textrm{max},2}, \end{equation}

by including a penalty function of the form,

(3.10)\begin{equation} f_{\kappa} = \frac{\int_{S_P} \,\textrm{d}^2 x \left(\exp\left({(|\kappa_1|^2-\kappa_{\textrm{max},1}^2)/w_k^2}\right) + \exp\left({(|\kappa_2|^2-\kappa_{\textrm{max},2}^2)/w_k^2}\right)\right)}{A_P}. \end{equation}

A similar constraint on the principal curvatures is employed in the ROSE code (Drevlak et al. Reference Drevlak, Beidler, Geiger, Helander and Turkin2018). We assume the convention that $\kappa _1 \le \kappa _2$ and $\kappa _{1,2}<0$ indicates concavity.

We additionally add a penalty on the smallest principal curvature,

(3.11)\begin{equation} f_{\bar{\kappa}} = \frac{\int_{S_P} \,\textrm{d}^2 x\, \kappa_1^2}{A_P}. \end{equation}

As large concavity is associated with coil complexity (Paul et al. Reference Paul, Landreman, Bader and Dorland2018), this regularization becomes critical in the two-staged optimization approach. Also, regularity of the boundary improves the convergence of the VMEC code. If the boundary is not close to being a star domain, indicating that there exists a coordinate axis such that the line segment connecting the axis and any point on the boundary is contained within the boundary, then the solver may fail to initialize a guess for the magnetic axis. For highly shaped boundaries, the coordinate surfaces may begin to overlap. Thus, the inclusion of these regularization terms prevents code failures during optimization and improves the convergence toward the optimum. This advantage of curvature penalization is highlighted in the magnetic well optimization in § 4.2.

3.3. Constrained volume optimization

To demonstrate the procedure described in § 3.1 for avoiding self-intersection, we define an objective function,

(3.12)\begin{equation} f(S_P) = V_P + \lambda_{\mathcal{S}} f_{\mathcal{S}} + \lambda_{R} f_R + \lambda_{\kappa} f_{\kappa}, \end{equation}

where $V_P$ is the volume enclosed by the plasma boundary. For this example, we take $\lambda _{\mathcal {S}} = \lambda _{R} = 10^3$, $R_\textrm {min} = \mathcal {R} = 0.2$ m, $w_R = w_{\mathcal {S}} = 10^{-2}$, $\lambda _{\kappa } = 1$, $\kappa _{\textrm {max},1} = \kappa _{\textrm {max},2} = 10\ \textrm {m}^{-1}$ and $w_k = 1$. Here the volume is computed using a surface integral upon application of the divergence theorem,

(3.13)\begin{equation} V_P = \int \,\textrm{d}^3 x = \int_{S_P} \,\textrm{d}^2 x\, Z \hat{\boldsymbol{n}} \boldsymbol{\cdot} \hat{\boldsymbol{Z}}, \end{equation}

to avoid discretizing the volume.

We begin with a stellarator whose boundary is defined by a rotating ellipse,

(3.14a)\begin{gather} R(\theta,\phi) = R_0 + 0.5a\left(\cos \left(\theta - N_{P} \phi\right) + \cos(\theta) \right) \nonumber\\ \quad\quad -\, 0.5 b \left(\cos\left(\theta - N_{P}\phi\right) -\cos(\theta)\right) \end{gather}

(3.14b)\begin{gather} Z(\theta,\phi) = 0.5b\left(\sin(\theta) + \sin\left(\theta - N_{P}\phi\right) \right) \nonumber\\ \quad\quad +\, 0.5 a \left(\sin(\theta) -\sin\left(\theta - N_{P} \phi\right)\right), \end{gather}

where $a$ is the semi-major axis and $b$ is the semi-minor axis. For this example, we take $R_0 = 5$, $a =2$, $b= 1$ and $N_P=3$. We optimize $f(S_P)$ with respect to modes $m \le 3$ and $|n| \le 3$ using the BFGS algorithm provided by SciPy. (This optimization algorithm is used for all demonstrations in this work.)

We compare the surfaces obtained with and without ($\lambda _{\mathcal {S}} = \lambda _{R} = 0$) the constraint terms, shown in figure 2. Without the constraint, the boundary begins to self-intersect. Although the volume is not well-defined for this surface, the discretized surface integral given by (3.13) is reduced from its initial value of $197.39\text{ to }4.18\ \textrm {m}^{3}$ by minimizing the projection of the normal vector in the $\hat {\boldsymbol {Z}}$ direction. With the constraint, the surface is reduced to an axisymmetric torus with a nearly circular cross-section centred at a major radius $R_0 = 0.43\ \textrm {m}$ with an averaged minor radius of $0.24\ \textrm {m}$ and a volume of $0.50\ \textrm {m}^3$.

Figure 2.

Cross-sections of the initial (blue) and optimized surfaces obtained by minimizing the objective function given by (3.12). The surface obtained without the constraints ($\lambda _R = \lambda _{\mathcal {S}} = 0$) becomes self-intersecting (red), and the surfaced obtained with the constraints remains regular (green).

In figure 3 we display the value of the global radius of curvature (3.6) for the initial and final surfaces. On the initial surface, $\rho _G$ is minimized at the endpoints of the semi-major axis where the global radius of curvature matches the local radius of curvature. The maximum value of $\rho _G$ is obtained in the region near the endpoint of the ellipse's semi-minor axis, where the point-tangent radius is minimized for $\boldsymbol {x}_2$ on the opposite side of the semi-minor axis. The optimized boundary features decreased values of $\rho _G$ on the outboard side, where $\rho _G$ matches the local radius of curvature. We note that the minimum value of $\mathcal {S}_{\mathcal {C}}$ on the optimized slightly violates the constraint at 0.15 m owing to the penalty formulation.

Figure 3.

The global radius of curvature (3.6) of the (a) initial surface (3.14) and (b) surface optimized for minimal volume with regularity constraints (3.12).

4.1. Target rotational transform

We now discuss a demonstration of optimization to obtain a target rotational transform profile. We consider a target function which quantifies the difference between the rotational transform and a desired profile, $\iota _{\text {target}}(\psi )$,

(4.1)\begin{equation} f_{\iota} = \tfrac{1}{2} \int_{V_P}\, \textrm{d}\psi \left(\iota(\psi)-\iota_{\text{target}}(\psi) \right)^2, \end{equation}

in order to obtain a low shear stellarator. The gradient of this objective function is obtained by computing an equilibrium with a perturbed toroidal current profile,

(4.2)\begin{equation} \delta I_T(\psi) = w(\psi) (\iota(\psi) - \iota_{\text{target}}(\psi)), \end{equation}

as described in Antonsen et al. (Reference Antonsen, Paul and Landreman2019). We present a benchmark problem in appendix A to demonstrate that the gradient-based optimizer can converge to a known minimum of an objective involving $f_{\iota }$ in a 2D space.

Although the VMEC code assumes nested surfaces such that the magnetic field is integrable, the solution to Laplace's equation with the prescribed boundary $S_P$ will not generally have continuously nested surfaces. We can, however, aim to obtain a vacuum field with a value of the rotational transform that improves the integrability of the ‘real’ vacuum field. (Note that although we can run VMEC with prescribed $\iota (\psi )$ rather than $I_T(\psi )$, it will not generally be a vacuum field.) We choose $\iota _{\text {target}}$ to be the ‘most irrational’ noble between $p/q = 0/1$ and $p'/q'=1/2$ (Greene, MacKay & Stark Reference Greene, MacKay and Stark1986),

(4.3)\begin{equation} \iota_{\text{target}} = \frac{p + \gamma p'}{q + \gamma q'} \approx 0.381966, \end{equation}

where $\gamma = (\sqrt {5} + 1)/2$ is the golden mean. This is a noble irrational, indicating that its path in the Farey tree is eventually alternating. At a given level in the Farey tree, the path that changes direction is said to be ‘most irrational’ if a cantorus with this rotation number has the smallest flux of field line trajectories (Meiss Reference Meiss1992). Thus, the value of $\iota _{\text {target}}$ is chosen to have the smallest flux for rotation numbers between $p/q = 0/1$ and $p'/q' = 1/2$. According to the Kolmogorov–Arnold–Moser (KAM) theorem, invariant circles with sufficiently irrational frequencies persist under small perturbations. Thus, choosing such a value for the rotational transform is likely to result in a large volume of magnetic surfaces. We remark that moderate shear may be desirable for certain stellarator design studies. Our adjoint formalism is quite flexible, and the choice of $\iota _{\text {target}}$ could be modified accordingly.

We define our objective function to be

(4.4)\begin{equation} f(S_P) = f_{\iota} + \lambda_R f_R + \lambda_{\mathcal{S}} f_{\mathcal{S}}, \end{equation}

where $\lambda _{\mathcal {S}} = \lambda _{R} = 10^3$, $R_\textrm {min} = \mathcal {R} = 0.2\ \textrm {m}$ and $w_R = w_{\mathcal {S}} = 10^{-2}$. We begin with a surface given by a rotating ellipse with a non-planar axis,

(4.5a)\begin{gather} R(\theta,\phi) = R_0 + R_{0,1} \cos(\phi) + 0.5a\left(\cos \left(\theta - N_{P} \phi\right) + \cos(\theta) \right) \nonumber\\ -\, 0.5 b \left(\cos\left(\theta - N_{P}\phi\right) -\cos(\theta)\right) \end{gather}

(4.5b)\begin{gather} \hskip8pt Z(\theta,\phi) ={-}\,R_{0,1} \sin(\phi) + 0.5b\left(\sin(\theta) + \sin\left(\theta - N_{P}\phi\right) \right) \nonumber\\ +\, 0.5 a \left(\sin(\theta) -\sin\left(\theta - N_{P} \phi\right)\right), \end{gather}

with $a = 2$, $b = 1$, $N_P = 3$, $R_0 = 5$ and $R_{0,1} = -0.5$. We consider a vacuum field with $p(\psi ) = I_T(\psi ) = 0$.

To investigate the benefit of increasing the dimensionality of the optimization space, we optimize with respect to the boundary harmonics $m \le 3$, $|n| \le 3$. We then use this result to optimize with respect to $m \le 4$, $|n| \le 4$ and then with respect to $m \le 5$, $|n| \le 5$. In comparison with the result of the initial optimization in the low-dimensional space, we are able to reduce the objective function by 70 % (figure 4). We obtain a rotational transform profile which very closely matches the target value with an objective value of $f_{\iota } = 4.93 \times 10^{-10}$. If the optimization space is further increased, the optimum is reduced by less than 1 %. For the further analysis of the optimization in this section, we present results from the optimization with respect to the modes $m \le 5$ and $|n| \le 5$.

Figure 4.

(a) The $f_{\iota }$ objective value (4.1), (b) the L2 norm of the $f(S_P)$ objective (4.4) gradient and (c) the initial and final profiles of the rotational transform. In (a,b) the optimization was performed in a staged approach: first with respect to modes $m \le 3$ and $|n| \le 3$, then with respect to $m \le 4$ and $|n| \le 4$ and, finally, with respect to $m \le 5$ and $|n| \le 5$.

We note that we were not able to reduce the gradient norm to the requested tolerance of $10^{-8}$. This result can be attributed to approximations made in computing the gradient with the adjoint method. As discussed in Antonsen et al. (Reference Antonsen, Paul and Landreman2019) and Paul et al. (Reference Paul, Antonsen, Landreman and Cooper2020), the adjoint equations are derived under the assumption that the VMEC code satisfies MHD force balance (1.1). However, there is always some residual error in force balance due to discretization error and the assumption of continuously nested flux surfaces. For these calculations, we converged to a force balance tolerance of $10^{-11}$ with 99 flux surfaces and mode numbers $|m| \le 11$, $|n| \le 11$. These resolution parameters were chosen to strike a balance between an accurate adjoint solution and efficiency of the optimization. Preliminary calculations indicate that increasing the resolution parameters over what was used in this work does not cause a significant change in the optimum. Furthermore, the linear adjoint solution is approximated by adding a small perturbation to the nonlinear MHD force balance. This technique is effectively a forward-difference approximation of the adjoint equation, which introduces an error that scales with the magnitude of the perturbation. When there are small errors in the gradient, the computed gradient may no longer be a descent direction, making convergence difficult near the optimum (Dekeyser Reference Dekeyser2014).

Nonetheless, we effectively eliminate the magnetic shear throughout the volume. In figure 5 we display the initial and optimized boundaries (solid) along with an interior surface at $\psi /\psi _0 = 0.06$ (dashed) and the magnetic axis (star). We can consider some of the features of the optimized surface in consideration of the expression for the rotational transform near the magnetic axis (Mercier Reference Mercier1964; Helander Reference Helander2014), which indicates that rotating ellipticity and torsion of the axis contribute to the on-axis transform in a vacuum field. We see that the interior surfaces become slightly more elliptical in order to increase $\iota$ near the axis, whereas the ellipticity of the boundary is slightly decreased, becoming more square. However, the torsion of the magnetic axis is maintained.

Figure 5.

Poloidal cross-sections of the initial boundary (4.5) and boundary which optimizes (4.4) (solid) along with the surface at $\psi /\psi _0 = 0.06$ (dashed) and magnetic axis (star).

To analyse the effect of low shear on the integrability of the field, we compute the vacuum field using the SPEC code (Hudson et al. Reference Hudson, Dewar, Dennis, Hole, McGann, Von Nessi and Lazerson2012). The SPEC calculations are performed with a single volume with Beltrami parameter $\mu = 0$ such that the magnetic field satisfies $\boldsymbol {\nabla } \times \boldsymbol {B} = 0$ with a Neumann boundary condition (1.1b). In figure 6 we show a Poincaré section for the initial field, which has a small island chain at the $\iota = 3/8$ resonance. With the optimized boundary, we eliminate this island chain and obtain a large volume of nested surfaces.

Figure 6.

Poincaré section computed from the SPEC vacuum field using the initial boundary (4.5) and boundary which minimizes (4.4): (a) initial; (b) optimized.

4.2. Magnetic well

We next consider an objective function which aims to obtain a magnetic well,

(4.6)\begin{equation} f_w = \frac{\int_{V_P} \,\textrm{d} \psi \left(w_1(\psi) V'(\psi) - w_2(\psi) V'(\psi)\right)}{\int_{V_P} \,\textrm{d} \psi \left(w_1(\psi) V'(\psi) + w_2(\psi) V'(\psi)\right)}, \end{equation}

with $w_1(\psi ) = \exp \left ({-(\psi /\psi _0-1.0)^2/0.1^2}\right )$ and $w_2(\psi ) = \exp \left ({-(\psi /\psi _0)^2/0.1^2}\right )$. When $V''(\psi ) < 0$, a magnetic well is said to be present, which provides a stabilizing term in the Mercier criterion for interchange modes (Mercier & Luc Reference Mercier and Luc1974). We can consider $f_w$ to be a normalized finite-difference approximation of $V''(\psi )$; thus minimization of $f_w$ is performed in order to achieve a magnetic well. A similar objective is employed in the ROSE code (Drevlak et al. Reference Drevlak, Beidler, Geiger, Helander and Turkin2018), computed from integration along a field line. Some experimental observations point toward a correlation between the existence of a magnetic hill and the onset of fluctuations (Castellano et al. Reference Castellano, Jiménez, Hidalgo, Pedrosa, Fraguas, Pastor, Herranz and Alejaldre2002). Nonetheless, we remark that the existence of a magnetic well does not always imply ideal MHD stability (Watanabe et al. Reference Watanabe, Sakakibara, Narushima, Funaba, Narihara, Tanaka, Yamaguchi, Toi, Ohdachi and Kaneko2005).

The gradient of this objective function is obtained by computing an equilibrium with a perturbed pressure profile,

(4.7)\begin{equation} \delta p(\psi) = \frac{w_1(\psi) - w_2(\psi) - f_w\left(w_1(\psi) + w_2(\psi)\right)}{\int_{V_P} \,\textrm{d}^3 x \left(w_1(\psi) + w_2(\psi)\right)}, \end{equation}

as described in Paul et al. (Reference Paul, Antonsen, Landreman and Cooper2020).

We define our objective function to be

(4.8)\begin{equation} f(S_P) = f_w + \frac{\lambda_V}{2}\left(V_P - V_P^{\text{init}}\right)^2 + \lambda_{\mathcal{S}} f_{\mathcal{S}} + \lambda_R f_R + \lambda_{\kappa} f_{\kappa} + \lambda_{\bar{\kappa}} f_{\bar{\kappa}}, \end{equation}

with $\lambda _{\mathcal {S}} = 10^{-3}$, $\lambda _{R} = 10^{2}$, $w_{\mathcal {S}} = w_{R} = 10^{-2}$, $\kappa _{\textrm {max}} \equiv \kappa _{\textrm {max},1} = \kappa _{\textrm {max},2} = 7\ \textrm {m}^{-1}$, $w_{\kappa } = 1$, $\lambda _{\bar {\kappa }}= \lambda _V = 1$ and $\lambda _{\kappa } = 10^{2}$. Here $V_{P}^{\text {init}}$ is the volume of the initial surface. To determine the importance of the maximum curvature regularization term, we optimize with four sets of parameters as described in table 1. The values $\lambda _{\kappa } = 10^2$ and $\kappa _{\textrm {max}} = 7\ \textrm {m}^{-1}$ were chosen to balance the curvature and well metrics.

Table 1.

Parameters used for the optimization of the magnetic well and the resulting values of the optimum objective function (4.8).

In addition to the regularization terms, we include a term in the objective function which penalizes a change in the volume, as an increase in the inverse aspect ratio can yield a Shafranov shift (§ 3.7 in Wesson & Campbell Reference Wesson and Campbell2011). This shift in the flux surfaces toward a smaller major radius implies that the volume of a flux surface increases less rapidly than its cross-sectional area. Furthermore, the flux through a surface increases more rapidly than its cross-sectional area because the geometric centre is moving into a region of increased field strength (assuming the field is mostly toroidal). Thus, the volume increases less rapidly than the flux, and a negative value of $V''(\psi )$ can be achieved (Taylor Reference Taylor1965).

We begin with a boundary given by a rotating ellipse (A 2) with $a = 2$, $b = 1$, $R_0 = 5$ and $R_{0,1} = 0$. We optimize the boundary with respect to the modes $m \le 10$ and $|n| \le 10$. The equilibrium is computed with $p(\psi ) = 0$ and $I_T(\psi ) = 0$ such that a vacuum field is considered. We arrive at the boundary given in figure 7. With the addition of the volume constraint, the aspect ratio remains roughly constant (3.54 for the initial boundary and 3.46 for the optimized boundary), so that the well is provided by the shaping of the boundary rather than the Shafranov shift that arises due to the inverse aspect ratio. We note that the optimized surface features triangularity that rotates from outward-pointing with horizontal elongation to inward-pointing with vertical elongation.

Figure 7.

Cross-sections of the boundary (orange) optimized for the magnetic well (4.8) with regularization on the curvature $(\lambda _{\kappa } = 10^2, \lambda _{\bar {\kappa }} = 1)$. Several interior surfaces are shown, along with the location of the magnetic axis (blue star) and axis normal vector (black arrow). The optimized boundary features rotating triangularity and elongation.

To demonstrate the effect of the curvature regularization terms, we examine the optimization with $\lambda _{\kappa } = \lambda _{\bar {\kappa }} = 0$. We converge to the boundary shown in figure 8. As can be seen, regions of large curvature are exhibited, including a dimple-like feature on the inboard side and a concave ‘pinching’ feature. Again the aspect ratio remains roughly constant (3.46).

Figure 8.

Cross-sections of the boundary optimized for the magnetic well (4.8) without regularization on the principal curvature $(\lambda _{\kappa } = \lambda _{\bar {\kappa }} = 0)$. Several interior surfaces are shown, along with the location of the magnetic axis (blue star) and axis normal vector (black arrow). The optimized boundary features a ‘dimple’ feature at the $\phi = 0$ plane and regions of increased concavity.

We can understand the geometric dependence of the magnetic well by considering the expression for $V''(\psi )$ on the axis that arises from the near-axis expansion in the inverse coordinate representation (Landreman & Jorge Reference Landreman and Jorge2020),

(4.9)\begin{equation} V''(\psi) \propto \int_0^{2{\rm \pi}}\,\textrm{d} \varphi_B \frac{1}{B_0^4}\left[3 (B_{1s}^2 + B_{1c}^2) - 4 B_0 B_{20}\right], \end{equation}

where $\varphi _B$ is the Boozer toroidal angle and we have made the assumption of a vacuum field. The field strength near the axis is expanded as

(4.10)\begin{align} B &= B_0(\varphi_B) + r \left[ B_{1c}(\varphi_B) \cos(\vartheta_B) + B_{1s}(\varphi_B) \sin(\vartheta_B)\right] \nonumber\\ &\quad + r^2 \left[B_{20}(\varphi_B) + B_{2c} \cos(2 \vartheta_B) + B_{2s} \sin(2 \vartheta_B) \right] + {O}\left(r^3\right), \end{align}

where $r \propto \sqrt {\psi }$ is the effective minor radius and $\vartheta _B$ is the Boozer poloidal angle.

As can be seen from (A34) in Landreman & Sengupta (Reference Landreman and Sengupta2019), the poloidally independent shift in the flux surface in the normal direction, $X_{20}$, increases proportional to $B_{20}$. Here, the normal vector is $\hat {\boldsymbol {n}} = \boldsymbol {r}''(l)/|\boldsymbol {r}''(l)|$, where $\boldsymbol {r}(l)$ is the position vector along the magnetic axis parameterized by the length. In the axisymmetric limit, a positive value of $X_{20}$ indicates a net shift of the surfaces’ geometric centre toward a smaller major radius. This correlation between $X_{20}$ and $B_{20}$ is consistent with the statement that a Shafranov shift provides a magnetic well. In figures 7 and 8 we display the magnetic axis (blue star) and normal vector (black arrow) along with the shapes of several magnetic surfaces for the optimized configurations. While both the unconstrained and constrained optima feature a net shift in the geometric centre of the flux surfaces in the normal direction, this is achieved with vastly different shaping of the boundary. Interestingly, both configurations feature negative triangularity at the $\phi = 0.5$ $(2{\rm \pi} /N_P)$ plane.

We note that $B_{20}$ arises owing to other shaping components of the surface (Landreman & Sengupta Reference Landreman and Sengupta2019). In axisymmetry, negative values of $B_{2c}$ contribute to positive (outward-pointing) triangularity ($X_{2c}<0$) ((B11) in Landreman Reference Landreman2021). Assuming stellarator symmetry, negative values of $B_{2c}$ contribute to positive values of $B_{20}$ if the surface is vertically elongated, and positive values of $B_{2c}$ contribute to positive values of $B_{20}$ if the surface is horizontally elongated. This implies that positive triangularity coupled with vertical elongation or negative triangularity coupled with horizontal elongation contributes to the magnetic well in axisymmetry. These trends are consistent with the stability analysis of oblate plasmas with negative triangularity (Pogutse & Yurchenko Reference Pogutse and Yurchenko1982; Kesner, Ramos & Gang Reference Kesner, Ramos and Gang1995; Medvedev et al. Reference Medvedev, Kikuchi, Villard, Takizuka, Diamond, Zushi, Nagasaki, Duan, Wu and Ivanov2015). In three dimensions, the connection between $B_{20}$ and $B_{2c}$ is more complicated (Landreman & Sengupta Reference Landreman and Sengupta2019, (A41) and (A42)). Nonetheless, we note that the magnetic well near the axis arises at second order in the expansion parameter, which includes the effects of ellipticity and triangularity. Thus, it is not surprising that the optimized boundary with curvature constraints features rotating ellipticity and triangularity.

In comparing the convergence of the optimization with and without the curvature constraints (figure 9), we see that the presence of the curvature objective prevents some failures of the VMEC code. Each function evaluation that resulted in a VMEC failure is visualized as a ‘spike’ that extends above $|f_w - f^{\text {opt}}_w| = 1$, as we artificially set the value of the objective function to $10^{12}$ at these points. Although the constrained optimization still features some VMEC failures (12 versus 19), the inclusion of the curvature constraints improves the convergence toward the optimum, and a deeper magnetic well is achieved throughout the volume.

Figure 9.

Convergence of the optimization of the objective function (4.8) (a) with the inclusion of the curvature constraints ($\lambda _{\kappa } = 10^2, \lambda _{\bar {\kappa }} = 1$) and (b) without the curvature constraints ($\lambda _{\kappa } = \lambda _{\bar {\kappa }} = 0$). (c) With the inclusion of the curvature constraint, we achieve a deeper magnetic well throughout the volume.

4.3. Quasi-symmetry

We next consider optimization for quasi-symmetry near the magnetic axis. This is quantified through the objective,

(4.11)\begin{equation} f_{QS} = \frac{\int_{V_P} \,\textrm{d}^3 x\, w(\psi) (B - \bar{B})^2}{\int_{V_P} \,\textrm{d}^3 x\, B^2}, \end{equation}

where $\bar {B} = \int _{V_P} \,\textrm {d}^3 x \, w(\psi ) B$ and we take $w(\psi ) = \exp \left ({-(\psi /\psi _0)^2/0.1^2}\right )$. This objective aims to make the field strength constant on the magnetic axis, which is a feature of quasi-axisymmetric and quasi-helically symmetric equilibria. A similar objective function has also been included in optimization for energetic particle confinement (Drevlak et al. Reference Drevlak, Geiger, Helander and Turkin2014). The gradient of this objective function is obtained by computing an equilibrium with the addition of an anisotropic pressure tensor, $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {P} = p_{\|} \hat {\boldsymbol {b}}\hat {\boldsymbol {b}} + p_{\perp } (\boldsymbol {I} -\hat {\boldsymbol {b}}\hat {\boldsymbol {b}})$ with

(4.12a)\begin{gather} p_{\|}(\psi,B) = \frac{w(\psi) (B - \bar{B})^2 - f_{QS}w(\psi) B^2}{\int_{V_P} \,\textrm{d}^3 x\, w(\psi) B^2}, \end{gather}

(4.12b)\begin{gather}p_{{\perp}} = p_{\|} - B \displaystyle\frac{\partial p_{\|}(\psi,B)}{\partial B}, \end{gather}

as described in Paul et al. (Reference Paul, Antonsen, Landreman and Cooper2020).

We take our objective function to be

(4.13)\begin{equation} f(S_P) = f_{QS} + \lambda_{\iota} f_{\iota} + \lambda_{\bar{\kappa}} f_{\bar{\kappa}} + \lambda_{\kappa} f_{\kappa} + \lambda_{\mathcal{S}} f_{\mathcal{S}} + \lambda_R f_R, \end{equation}

where $\iota _{\text {target}}(\psi )$ is taken to be the initial rotational transform profile. We also take $\lambda _{\iota } = 1$, $\lambda _{\bar {\kappa }} = 10^{-3}$, $\lambda _{\kappa } = 5 \times 10^{-2}$, $\kappa _{\textrm {max},1} = \kappa _{\textrm {max},2} = 10\ \textrm {m}^{-1}$, $w_{\kappa } = 1$, $w_{\mathcal {S}} = w_R = 10^{-1}$, $\lambda _R = 10^3$, $\lambda _{\mathcal {S}} = 10^{-1}$, $\mathcal {R} = 0.5\ \textrm {m}$ and $R_\textrm {min} = 0.2\ \textrm {m}$. The additional constraint on the rotational transform is required to prevent the surface from becoming axisymmetric to reduce the toroidal ripple on the axis.

We begin with a rotating ellipse boundary with torsion, given by

(4.14a)\begin{align} &R(\theta,\phi) = R_0 + R_{0,1} \cos(\phi) + R_{0,2}\cos(2\phi) + 0.5a\left(\cos \left(\theta - N_{P} \phi\right) + \cos(\theta) \right) \nonumber\\ & \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad -\, 0.5 b \left(\cos\left(\theta - N_{P}\phi\right) -\cos(\theta)\right) \end{align}

(4.14b)\begin{align} &Z(\theta,\phi) = R_{0,1} \sin(\phi) + Z_{0,2}\cos(2\phi) + 0.5b\left(\sin(\theta) + \sin\left(\theta - N_{P}\phi\right) \right) \nonumber\\ & \quad \quad \quad \quad \quad \quad \quad \quad +\, 0.5 a \left(\sin(\theta) -\sin\left(\theta - N_{P} \phi\right)\right), \end{align}

with $R_0 = 5$, $R_{0,1} = R_{0,2} = 0.3$, $a = 2$, $b = 1$ and $N_P = 3$. We perform optimization at $\beta = 1.2\,\%$ with the profiles shown in figure 10. We optimize with respect to the modes $m \le 10$ and $|n| \le 10$ in a staged approach: first optimizing with respect to $m \le 3$ and $|n| \le 3$, then with respect to $m \le 6$ and $|n| \le 6$, then with respect to the full set of modes.

Figure 10.

(a) The pressure and (b) the integrated toroidal current profiles used for the optimization of quasi-symmetry on the axis (4.13).

We display the magnetic field strength on the initial and optimized boundaries in figure 11. We note the initial configuration features a large toroidal variation of the field strength, with increased field strength near the ‘corners’ that arise owing to the large torsion of the axis. Although the axis ripple figure of merit (4.11) reduces the toroidal variation of the field strength on the axis, which must vanish in both quasi-axisymmetry and quasi-helical symmetry, we find that the optimized magnetic field is driven closer to quasi-axisymmetry.

Figure 11.

The magnetic field strength on (a) the initial boundary (4.14) and (b) the boundary optimized for quasi-symmetry on the axis (4.13).

We remark that a design based solely on the objective (4.11) may be limited in that quasi-symmetry is not targeted away from the axis. We therefore quantify the departure from quasi-symmetry in the two configurations in two additional ways. The first is through the figure of merit,

(4.15)\begin{equation} f_{TP}(\psi) =\frac{\left \langle \left \rvert \displaystyle\frac{\partial B}{\partial \theta} \displaystyle\frac{\partial }{\partial \phi} \left(\boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} B \right) - \displaystyle\frac{\partial B}{\partial \phi} \displaystyle\frac{\partial }{\partial \theta}\left(\boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} B\right) \right \rvert \right \rangle_{\psi}}{\langle B^3 \rangle_{\psi}/r_{\text{minor}}(\psi)}, \end{equation}

where $r_{\text {minor}}(\psi ) = \sqrt {A(\psi )/{\rm \pi} }$ is the averaged minor radius of a flux surface where $A(\psi )$ is the toroidally averaged cross-sectional area. (This is the same definition of the minor radius used in the VMEC code.) The flux-surface average of a quantity $Q$ is

(4.16)\begin{equation} \langle Q \rangle_{\psi} = \frac{\int_0^{2{\rm \pi}} \,\textrm{d}\theta \int_0^{2{\rm \pi}} \,\textrm{d} \phi Q \sqrt{g}}{\int_0^{2{\rm \pi}} \,\textrm{d} \theta \int_0^{2{\rm \pi}} \,\textrm{d} \phi \sqrt{g}}, \end{equation}

where $\sqrt {g} = (\boldsymbol {\nabla } \psi \times \boldsymbol {\nabla } \theta \boldsymbol {\cdot } \boldsymbol {\nabla } \phi )^{-1}$ is the flux coordinate Jacobian. Equation (4.15) is a normalized figure of merit which employs the triple product form for quasi-symmetry (Helander Reference Helander2014; Rodriguez, Helander & Bhattacharjee Reference Rodriguez, Helander and Bhattacharjee2020),

(4.17)\begin{equation} \boldsymbol{\nabla} \psi \times \boldsymbol{\nabla} B \boldsymbol{\cdot} \boldsymbol{\nabla} \left(\boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} B \right) = 0, \end{equation}

and allows us to quantify the quasi-symmetry error without specifying the helicity of the symmetry. We also perform a Boozer coordinate ($\psi$,$\vartheta _B$,$\varphi _B$) transformation (Sanchez et al. Reference Sanchez, Hirshman, Ware, Berry and Spong2000) to obtain the Fourier harmonics of the field strength,

(4.18)\begin{equation} B(\psi,\vartheta_B,\varphi_B) = \sum_{m,n} B_{m,n} \cos(m \vartheta_B - n \varphi_B). \end{equation}

The deviation from quasi-axisymmetry can then be quantified as

(4.19)\begin{equation} f_{QA}(\psi) = \sqrt{\frac{\sum_{m,n\ne 0} B_{m,n}(\psi)^2}{\sum_{m,n}B_{m,n}(\psi)^2}}. \end{equation}

Although quasi-symmetry was only targeted on the axis, we reduce the quasi-symmetry error throughout the volume with respect to both metrics without introducing a reduction in the rotational transform (figure 12).

Figure 12.

(a) The rotational transform, (b) the quasi-axisymmetric error (4.19) and (c) the triple product error (4.15) for the initial boundary (4.14) and boundary optimized for small magnetic ripple on the axis (4.13).