3.1. Balanced system

We obtain a balanced system by using a parametrization $\boldsymbol \varSigma _s$ of $\boldsymbol \varSigma$ with unknowns $\boldsymbol{s}\in \mathbb {R}^{n_s}$, requiring (3.5) to be satisfied at the same number of collocation points, i.e.

(3.7)\begin{equation} \left.\begin{gathered} \boldsymbol{r}_{i,j}(\boldsymbol{s}) := G\boldsymbol{B}_{i,j} -\|\boldsymbol{B}_{i,j}\|^2 \left(\frac{\partial \boldsymbol \varSigma_{s,i,j}}{\partial \varphi}+\iota \frac{\partial \boldsymbol \varSigma_{s,i,j}}{\partial \theta} \right) = 0,\quad \forall\ (i,j) \in \mathcal{C} \\ V(\boldsymbol\varSigma_s)- V_{\text{target}} = 0 \end{gathered}\right\} , \end{equation}

where $\mathcal {C}$ is the set of indices of collocation points on a regular grid, $\boldsymbol \varSigma _{s,i,j} = \boldsymbol \varSigma _{s}(\varphi _{i,j}, \theta _{i,j})$, and $V(\boldsymbol \varSigma _s)$ is the volume enclosed by the surface, computed as

(3.8)\begin{equation} V(\boldsymbol\varSigma_s) = \frac{1}{3 n'_{\varphi} n'_{\theta}} \sum_i \sum_j \boldsymbol \varSigma_{s,i,j} \boldsymbol{\cdot} \left(\frac{\partial \boldsymbol \varSigma_{s,i,j}}{\partial \varphi} \times \frac{\partial \boldsymbol \varSigma_{s,i,j}}{\partial \theta} \right)\quad \text{for } (i,j) \in \mathcal{C}', \end{equation}

where $\mathcal {C}'$ is a set of quadrature points on a regular $n'_{\varphi } \times n'_{\theta }$ grid that is possibly different from the ones used by the collocation method $\mathcal {C}$ to ensure an accurate computation of the volume enclosed by the surface. The equations (3.7) require that the residual at a set of collocation points $(\varphi _{i,j}, \theta _{i,j})$ be zero. We place a grid of $(2n_{\textrm {tor}} +1) \times (2m_{\textrm {pol}}+1)$ equispaced collocation points on $[0,1/n_{\textrm {nfp}})\times [0,1)$ and refer to surfaces that satisfy (3.7) at each point on that grid as BoozerExact surfaces.

Assuming stellarator symmetry, a number of these points are redundant and we choose a subset of them. First, the $x$-component of the residual at the quadrature point associated with the origin, $i = j = 0$, or $\varphi = \theta = 0$, is always zero for all surfaces, so this equation can be ignored. To see why this is the case, note that the $x$-component of stellarator symmetric magnetic fields is zero when $\varphi =\theta = 0$. Similarly, differentiating the first equation in (2.7) with respect to $\varphi$ or $\theta$ shows that the $x$-component of the surface tangents are zero when $\varphi =\theta = 0$. Second, if we include $(\varphi _{i,j}, \theta _{i,j})$ then the collocation point associated with $(-\varphi _{i,j}, -\theta _{i,j})$ is unnecessary. This is because the residual evaluated at $(\varphi _{i,j}, \theta _{i,j})$ is the residual at $(-\varphi _{i,j}, -\theta _{i,j})$ subject to a reflection. This is illustrated in figure 1(a), where the collocation points are grouped by colour. We are free to choose one green, blue, purple and red grouping. Choosing points from two groups of the same colour would result in redundant equations. In all examples, we place collocation points on $[0,1/2n_{\textrm {fp}})\times [0,1)$ from zones 1a, 2a, 3a and 4a, resulting in, respectively, $n_{\varphi } = n_{\textrm {tor}} + 1$ and $m_{\theta } = 2 m_{\textrm {pol}} + 1$ quadrature points in the toroidal and poloidal directions (before removing redundant quadrature points in group 3b and at the origin).

After removing the redundant collocation points, the number of equations that must be satisfied is

(3.9)\begin{equation} 3[2n_{\text{tor}}m_{\text{pol}} + m_{\text{pol}} + n_{\text{tor}} + 1] - 1 + 1 = 3[2n_{\text{tor}}m_{\text{pol}} + m_{\text{pol}} + n_{\text{tor}} + 1]. \end{equation}

This is the same as $n_s+2$, the number of geometric degrees of freedom in the parametrization $n_s$ with $\iota$ and $G$. Since this is a balanced system of nonlinear equations, we can use Newton's method to solve it.

3.2. Overdetermined system

As we will show in numerical examples, the above balanced approach can fail to converge when the magnetic field lines have chaotic regions and islands. For such problems, we propose an alternative approach in which we relax the requirement that the residual be zero everywhere. Instead, we compute a surface that solves (3.5) in a least-squares sense; that is, it minimizes

(3.10)\begin{equation} f(\boldsymbol{s}) := \frac{1}{2} \sum_{i,j} \|\boldsymbol r_{i,j}(\boldsymbol s) \|^2 + \frac{1}{2} w \left( V(\boldsymbol\varSigma_s)- V_{\text{target}} \right)^2. \end{equation}

Here, the nonlinear least squares objective contains the sum of squares of (3.5) at a number of collocation points with a penalty term with weight $w>0$ to approximately enforce the surface volume condition (3.5b). The first-order optimality condition for a minimizer of $f$ is

(3.11)\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{s}}\kern0.09em f = 0 \end{equation}

and the optimization problem can be solved using a Levenberg–Marquardt method or Newton's method. In numerical experiments, we typically first compute a solution with the Levenberg–Marquardt method, which is known to work robustly for nonlinear least squares problems (Nocedal & Wright Reference Nocedal and Wright2006). If the stationarity condition (3.11) is not satisfied to sufficient precision, then we further improve the solution using a few steps of Newton's method. We refer to surfaces that satisfy (3.11) as BoozerLS (‘Boozer least squares’) surfaces. Note that the BoozerLS surface has similarities with the quadratic flux minimizing (QFM) surfaces of Dewar, Gibson & Hudson (Reference Dewar, Gibson and Hudson2010). The QFM surfaces are obtained by finding surfaces so that $(\boldsymbol B\boldsymbol {\cdot }\boldsymbol n)^2$ is minimized on them. This results in surfaces that are approximately tangential to the magnetic field. Surfaces that satisfy (3.5a) have the stricter requirement that the magnetic field must be aligned in the direction ${\partial \boldsymbol \varSigma }/{\partial \varphi }+\iota ({\partial \boldsymbol \varSigma }/{\partial \theta })$ and not just any tangential direction. In all examples, we place equispaced collocation points on $[0,1/2n_{\textrm {fp}})\times [0,1)$, resulting in, respectively, $n_{\varphi } = n_{\textrm {tor}} + 6$ and $m_{\theta } = 2 m_{\textrm {pol}} + 6$ quadrature points in the toroidal and poloidal directions.

3.3. Example: surfaces of an NCSX-like magnetic field

In this section, we use a challenging magnetic field to study the behaviour of the BoozerExact and BoozerLS approaches for computing surfaces. Specifically, we use the 18 non-planar coils and currents, corresponding to the C09R00 boundary configuration, of the national compact stellarator experiment (NCSX) to generate a vacuum magnetic field. The original NCSX design included planar toroidal and poloidal field coils that we do not use here. We have also previously considered this set-up in § 6 of Giuliani et al. (Reference Giuliani, Wechsung, Cerfon, Stadler and Landreman2022). The field is challenging as it does not have nested flux surfaces everywhere, and presents an island chain in the neighbourhood of the low-order rational $\iota =3/7$. The magnetic field $\boldsymbol{B}$ on the surface $\boldsymbol \varSigma$ in (3.5a) is evaluated using the Biot–Savart law since we are working with vacuum fields. Newton's method requires the Jacobian of (3.7) or (3.11), which are functions of the magnetic field and its first spatial derivatives for (3.7), and of its first and second spatial derivatives for (3.11). Since the field is computed with the Biot–Savart law, we have explicit formulae to compute these quantities.

Due to the local convergence properties of the Newton and Levenberg–Marquardt algorithms to compute surfaces, they require initializations that should not be too far from the target surface. We use a continuation approach, i.e. we first compute a surface close to the axis, i.e. with a small value of $V_\textrm {target}$. For that purpose, we generate a toroidal surface centred on the magnetic axis with a fixed minor radius. This does not correspond to a magnetic flux surface. However, in the neighbourhood of the magnetic axis it is typically close enough for Newton's method to converge. We then use previously computed surfaces as initialization for finding the surface for the next larger $V_\textrm {target}$.

In figure 2(a), we show cross-sections of surfaces obtained with the BoozerExact approach, i.e. using the same number of unknowns as equations. In particular, we solve (3.7) with $m_{\textrm {pol}} = n_{\textrm {tor}} = 11$. As can be seen, when the Poincaré plot suggests the existence of magnetic surfaces, the BoozerExact approach converges to surfaces that align well with the Poincaré plot. In chaotic regions or regions containing magnetic islands, the approach may not converge, or it may find self-intersecting surfaces (as shown in red in the figure).

Figure 2.

Poincaré plots (in grey) and surface cross-sections for the field generated from a set of NCSX modular coils at $\phi = {\rm \pi}/3$. In (a), we show cross-sections of surfaces computed with the BoozerExact approach (in red and black). The red cross-section is self-intersecting as it is close to an island chain. In (b) we show cross-section obtained with the BoozerLS algorithm (in black, green and cyan), which results in nested surfaces that do not self-intersect in the presence of islands. The red and cyan cross-sections in (a,b) have the same target surface label, but are computed using BoozerExact and BoozerLS formulations, respectively. Green cross-sections correspond to additional surfaces not shown in (a). In (c), we show two BoozerExact surfaces that have the same enclosed surface volume.

We next show that the BoozerLS approach is more robust in the presence of islands and chaotic regions. We again use $m_{\textrm {pol}} = n_{\textrm {tor}} = 11$ and the volume constraint penalty $w=1000$, which results in a small relative error in the volume contained by the surface, ranging from $10^{-9}$ to $10^{-5}$. The resulting cross-sections are shown in figure 2(b). It can be seen that when the Poincaré plot indicates the existence of good surfaces, the BoozerLS surface cross-sectionsvisually coincide with the Poincaré plot and with the cross-sections obtained with the balanced BoozerExact approach. In other regions, the BoozerLS approach results in layered smooth surfaces that do not have self-intersections (green curves). The BoozerLS surfaces can either be computed using a standard continuation procedure, or alternatively, using possibly self-intersecting BoozerExact surfaces as initial guesses.

Finally, we numerically show that BoozerExact surfaces, i.e. solutions to (3.7), are not necessarily unique in the presence of chaos. This is more of theoretical interest than a fundamental problem as this is still a practical approach for computing magnetic surfaces when they exist. In addition, we will show in § 4 that it also can be used to optimize coils for quasisymmetry. In figure 2(c), we show two BoozerExact surface cross-sections with the same target volume enclosed, $V_{\textrm {target}}=1$. The red surfaces in figures 2(a) and 2(c) are the same. The cyan cross-section in figure 2(c) is a BoozerExact surface initialized from the cyan BoozerLS surface in figure 2(b). It is still possible that the BoozerLS solver does not converge or for a self-intersecting BoozerLS surface to be found. However, we find that this happens less frequently than it does for BoozerExact surfaces.

To further study the BoozerLS method in the presence of islands or chaotic fields, we examine the behaviour of the residual as the number of Fourier modes ($n_{\textrm {tor}}$, $m_{\textrm {pol}}$) and collocation points used in the parametrization of the surface increases (figure 3). As the number of degrees of freedom and collocation points in $\boldsymbol \varSigma _s$ increase, the residual decreases. However, this is only true for surfaces that do not pass through island chains. In the neighbourhood of the island chain, we observe that the residual does not decrease with increasing surface complexity and appears to stagnate. This hints at an approach to optimize for magnetic fields without islands: adding the square of the residual at the collocation points of the BoozerLS surfaces to the objective function. This approach is not used in the coil optimization discussed next, which relies on BoozerExact surfaces. However, it will be the subject of future work.

Figure 3.

BoozerLS residual in the NCSX-lite configuration with increasing number of Fourier modes used to represent the surface. Panel (a) shows the BoozerLS residuals as a function of the surface label. We observe that it stagnates in the neighbourhoodof $\iota =3/7$, which corresponds to an island chain. Panel (b) shows the rotational transform with respect to surface label for different resolutions. The vertical dashed line approximately indicates the volume for which $\iota =3/7$ is achieved.

4.1. Optimization problem formulation

Now that we are able to compute surfaces and a measure for quasisymmetry, the goal is to formulate and solve optimization problems for coils generating magnetic fields with good quasisymmetry properties on a large number of toroidal surfaces. We seek a design with $N_c$ independent modular coils, to which stellarator and rotational symmetries are applied. After application of these symmetries, the stellarator is made up of $2n_{\textrm {fp}}N_c$ modular coils. Each coil is parametrized with a current and $n_c$ geometric degrees of freedom using a Fourier basis as in Giuliani et al. (Reference Giuliani, Wechsung, Cerfon, Stadler and Landreman2022). All coil degrees of freedom and their currents are summarized into a vector $\boldsymbol{c}\in \mathbb {R}^{N_cn_c + N_c}$. Our guiding principle for designing the coils is to target quasisymmetry of the magnetic field induced by these coils on $N_s$ surfaces. Each surface uses the parametrization from § 2 with $n_s$ geometric degrees of freedom. The resulting surface Fourier coefficients, $\iota$ and $G$ are summarized into a vector $\boldsymbol{s}\in \mathbb {R}^{N_sn_s+2N_s}$. In the optimization formulation below, we use BoozerExact surfaces (§ 3.1). Finally, the optimization also includes various (inequality) constraints for the coils that are motivated by engineering considerations, which are critical in terms of the feasibility and cost-efficiency of a design (Strykowsky et al. Reference Strykowsky, Brown, Chrzanowski, Cole, Heitzenroeder, Neilson, Rej and Viol2009; Neilson et al. Reference Neilson, Gruber, Harris, Rej, Simmons and Strykowsky2010; Klinger et al. Reference Klinger, Baylard, Beidler, Boscary, Bosch, Dinklage, Hartmann, Helander, Maßberg and Peacock2013). These coil constraints are typically not difficult to compute as they explicitly depend on the coil degrees of freedom. This is in contrast to the terms that depend on the surfaces, which have an implicit dependence on the coil degrees of freedom. The optimization problem is

(4.3a)\begin{gather} \min_{\boldsymbol{c} \in \mathbb{R}^{N_{\!c} n_{c} + N_c},\quad \boldsymbol{s}\in \mathbb{R}^{N_{\!s} n_s + 2 N_s}} \quad \hat J (\boldsymbol{c}, \boldsymbol{s}) \end{gather}

(4.3b)\begin{gather}\text{subject to }\, (3.7)\, \text{on each of the}\, N_s\, \text{surfaces,} \end{gather}

(4.3c)\begin{gather}\frac{1}{N_s}\sum_{k=1}^{N_s}\iota_k = \bar{\iota}, \end{gather}

(4.3d)\begin{gather}\varPsi = \varPsi_0, \end{gather}

(4.3e)\begin{gather}R_{\text{major}} = R_0, \end{gather}

(4.3f)\begin{gather}\sum_{i = 1}^{N_c} L^{(i)}_{c} \leq L_{\max}, \end{gather}

(4.3g)\begin{gather}\kappa_i \leq \kappa_{\max},\quad i = 1, \ldots, N_c, \end{gather}

(4.3h)\begin{gather}\frac{1}{L^{(i)}_{c}}\int_{\boldsymbol\varGamma^{(i)}} \kappa_i^2\,{\rm d} l \leq \kappa_{\text{msc}},\quad i = 1, \ldots, N_c, \end{gather}

(4.3i)\begin{gather}\| \boldsymbol\varGamma^{(i)}- \boldsymbol\varGamma^{(j)} \|\geq d_{\min} \quad \text{for}\ i \neq j, \end{gather}

where the objective function is the average (normalized) non-quasiaxisymmetry on the surfaces $\boldsymbol \varSigma _{s,k}\ k= 1, \ldots, N_s$:

(4.4)\begin{equation} \hat{J}(\boldsymbol{c},\boldsymbol{s}) = \frac{1}{N_s}\sum^{N_s}_{k=1} \frac{\displaystyle\int_{\boldsymbol\varSigma_{s,k} }B_{\text{non-QA}, k}(\boldsymbol{c}, \boldsymbol{s})^2\,{\rm d} S}{\displaystyle\int_{\boldsymbol\varSigma_{s,k} }B_{\text{QA}, k}(\boldsymbol{c}, \boldsymbol{s})^2\,{\rm d} S}.\end{equation}

The above objective is the average ratio of the squared 2-norm of the field magnitude's non-quasisymmetric and quasisymmetric components. The objective is subject to the surface constraints (4.3b)–(4.3e) and the coil constraints (4.3f)–(4.3i). In particular, (4.3c) ensures that the average rotational transform across the surfaces is $\bar {\iota }$. Equation (4.3d) fixes the toroidal flux on the innermost surface to a given value $\varPsi _0$ and prevents the currents from going to zero. Equation (4.3e) fixes the major radius on the innermost surface to a given $R_0$ and prevents the length scale of the stellarator from changing. Equation (4.3f) prevents the sum of the independent modular coil lengths $\sum _{i = 1}^{N_c} L^{(i)}_{c}(\boldsymbol {c})$ from exceeding a given value $L_{\max }>0$. Equations (4.3g) and (4.3h), respectively, prevent the curvature and mean squared curvature on each coil from exceeding values $\kappa _{\max }$ and $\kappa _{\textrm {msc}}$. Finally, (4.3i) ensures that the coils stay at least $d_{\min }>0$ away from one another.

4.2. Computational aspects

Our computational approach to finding minimizers for (4.3) uses a gradient-based descent algorithm to minimize the objective. This requires appropriate gradients of the objective, which take into account the equality PDE constraints (4.3b) describing the surfaces, as well as the remaining equality and inequality constraints (4.3c)–(4.3i). To enforce the PDE constraints (4.3b), we use the method of Lagrange multipliers (Tröltzsch Reference Tröltzsch2010; Giuliani et al. Reference Giuliani, Wechsung, Cerfon, Stadler and Landreman2022). For that purpose, we consider the objective as a function just of the coil degrees of freedom, i.e. $J(\boldsymbol {c}) := \hat J(\boldsymbol {c}, \boldsymbol {s}(\boldsymbol {c}))$ by considering the surfaces as functions of the coils, i.e. $\boldsymbol{s} = \boldsymbol{s}(\boldsymbol{c})$ implicitly defined through the solution of (3.7). Computing the gradient of $J$ with respect to $\boldsymbol{c}$ can be done efficiently by introducing Lagrange multipliers and using the adjoint method. Computing gradients then requires solving two systems of equations per surface, both with $n_{s}+2$ equations: the nonlinear system of (3.7) and a linear system for the adjoint variable. This method compares favourably with finite difference approaches for computing gradients, where the number of times (3.7) must be solved scales with $n_c$, the number of degrees of freedom used to represent the coils.

The constraints (4.3c)–(4.3i) are enforced approximately using a penalty method, with penalties that only become active when the constraint is violated. We next detail the penalty terms used, and discuss a practical strategy to choose and adjust the weights for these penalties at the end of this section.

For the penalty on the sum of the independent modular coil lengths, we use the standard quadratic penalty

(4.5)\begin{equation} \frac{1}{2}\left[\max \left( 0, \sum_{i = 1}^{N_c} L^{(i)}_{c}(\boldsymbol{c}) - L_{\max} \right) \right]^2. \end{equation}

The penalty on the maximum curvature for coil $\boldsymbol \varGamma ^{(i)}$ is

(4.6)\begin{equation} \sum_{i=1}^{N_c}\frac{1}{2} \int_{\boldsymbol \varGamma^{(i)}} \left[\max \left( 0, \kappa_i(\boldsymbol{c}) - \kappa_{\text{max}} \right) \right]^2{\rm d} l_i, \end{equation}

and, as in Giuliani et al. (Reference Giuliani, Wechsung, Cerfon, Stadler and Landreman2022), the penalty on the minimum pairwise coil distance is

(4.7)\begin{equation} \frac{1}2\sum_{i\neq j} \int_{\boldsymbol \varGamma^{(i)}} \int_{\boldsymbol \varGamma^{(j)}} [\max(0, d_\text{min} - \| \boldsymbol \varGamma^{(i)}-\boldsymbol \varGamma^{(j)} \|_2)]^2 {\rm d} l_j\,{\rm d} l_i. \end{equation}

For the equality constraints, we add quadratic penalties to the objective. For instance, to prevent the currents in the coils from becoming zero, for the toroidal flux $\varPsi$ of the innermost surface we use

(4.8)\begin{equation} \tfrac{1}{2}\left( \varPsi(\boldsymbol{c}, \boldsymbol{s})-\varPsi_0 \right)^2. \end{equation}

Here, the toroidal flux through a toroidal cross-section of constant Boozer angle $\varphi =0$, $S_{\varphi =0}$, of a surface is given by

(4.9)\begin{equation} \varPsi(\boldsymbol{c}, \boldsymbol{s}) = \int_{S_{\varphi=0}} \boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{n}\,{\rm d} S = \int_0^1 \boldsymbol{A} \boldsymbol{\cdot} \frac{\partial}{\partial \theta} \boldsymbol\varSigma_s(0,\theta)\,{\rm d} \theta, \end{equation}

where $\boldsymbol{A}$ is the vector potential, such that $\boldsymbol{B} = \boldsymbol {\nabla } \times \boldsymbol {A}$, and where we have used Stokes’ theorem to turn the surface integral for the flux into a line integral. Additionally, we include a penalty term on the major radius $R_\textrm {major}$ of the innermost surface:

(4.10)\begin{equation} \tfrac{1}{2} (R_{\text{major}}(\boldsymbol{s})-R_0)^2. \end{equation}

This provides a characteristic length scale to the stellarator configurations we optimize, and maintains the aspect ratio of the surfaces in the volume. Details of how we compute this quantity are given in Appendix B. To constrain the mean rotational transform of the surfaces to be close to a given value $\overline \iota$, we use the penalty term

(4.11)\begin{equation} \frac{1}{2}\left(\frac{1}{N_s}\sum_{k=1}^{N_s}\iota_k(\boldsymbol{s}) - \bar{\iota}\right)^2, \end{equation}

where $\iota _k$ is the rotational transform on the $k$th surface. While here we target the average rotational transform, other choices are possible in our framework, such as constraining any number of surface rotational transforms to a given value or interval. Finally, we include a penalty that favours coils with an approximately uniform incremental arclength

(4.12)\begin{equation} \sum^{N_c}_{i=1}\int^1_0 (\|\boldsymbol \varGamma^{(i)'}\|_2 - t^{(i)})^2\,{\rm d} \xi, \end{equation}

with $t^{(i)} = \int ^1_0 \|\boldsymbol \varGamma ^{(i)'}\|_2\,\textrm {d}\xi$, where $\xi \in [0,1)$ is the angle parameterization of the coil. This prevents the incremental arclengths from getting close to zero, where the length penalty term is non-differentiable. Note that it might not be possible to attain a uniform incremental arclength at all points on the coil. However, this is not necessary as the purpose of this penalty is to avoid the above issue. Numerical experiments confirm that this term does not substantially affect the value of the objective over a wide range of penalty weights.

We compute analytical gradients of the (reduced) objective $J(\boldsymbol{c})$, and use the Broyden–Fletcher–Goldfarb–Shanno (BFGS) method (Nocedal & Wright Reference Nocedal and Wright2006) with linesearch to find a minimizer. Compared with steepest descent, BFGS has the advantage that it approximates second-order information of the objective, which typically results in faster convergence or, in some cases, avoids that the iteration progress stalls. To choose appropriate penalty weights, we follow a simple strategy that appropriately increases these weights throughout the iteration. Namely, if any of the constraints are not satisfied to within a $0.1\,\%$ relative tolerance after 3000 iterations of BFGS, we increase the weights on the associated penalty terms by a factor of 10. Additionally, we allow for a maximum number of 45 000 BFGS iterations. However, since we choose a judicious set of initial weights informed from previous runs, weights are increased rarely.

4.3. Coil optimization for precise quasiaxisymmetry on surfaces

In this example, we revisit a coil set obtained in Wechsung et al. (Reference Wechsung, Landreman, Giuliani, Cerfon and Stadler2022b) with the second stage of the standard two-stage optimization method, where the goal was to produce the precise quasiaxisymmetry (QA) magnetic field presented in Landreman & Paul (Reference Landreman and Paul2022) using a realistic set of coils. In the first stage, that field was optimized for precise QA on a toroidal volume with aspect ratio 6, and the average rotational transform in the volume was 0.42. The stage II optimization found $N_c=4$ coils with a minimum pairwise coil distance of $d_{\min } = 0.1\,\textrm {m}$, and maximum curvature and mean squared curvature of $\kappa _{\max } = 5\,\textrm {m}^{-1}$ and $\kappa _{\textrm {msc}} = 5\,\textrm {m}^{-2}$. Various coil designs are found where the sum of the lengths of the four base coils is constrained to be below $L_{\max } = 18, 20, 22, 24\,\textrm {m}$. Here, the major radius of the device is approximately 1 m. In the following optimization, we use the same values for $d_{\min }$, $\kappa _{\max }$, $\kappa _{\textrm {msc}}$ as in the stage II optimization from Wechsung et al. (Reference Wechsung, Landreman, Giuliani, Cerfon and Stadler2022b).

Building on the results from this two-stage coil design, we aim at improving the quasisymmetry on multiple surfaces uniformly distributed throughout the volume using the optimization approach presented in the previous section. Since the fields resulting from the coils in Wechsung et al. (Reference Wechsung, Landreman, Giuliani, Cerfon and Stadler2022b) appear to avoid islands, we use BoozerExact (rather than BoozerLS) surfaces in this optimization. Should islands appear in the vicinity of the surfaces over the course of the optimization, the BoozerExact surfaces might begin to self-intersect and no longer correspond to flux surfaces of the magnetic field that we are trying to optimize. However, as will be shown, this does not happen in the final coil designs. Generalizing the formulation (4.3) to BoozerLS surfaces, in order to handle magnetic fields with islands in a robust manner, and to optimize the coils to eliminate such islands, will be part of future work.

We first optimize coils using five BoozerExact surfaces uniformly spaced in terms of the surface label $\sqrt V$, which is proportional to the minor radius. We choose $m_{\textrm {pol}} = n_{\textrm {tor}} = 10$, and the outermost surface in the initial configuration has aspect ratio approximately equal to 6, similar to the optimization in Landreman & Paul (Reference Landreman and Paul2022). We set the target toroidal flux (4.3d) and major radius (4.3e) of the innermost surface to be the same as the one in the stage II coil set. In addition, we use the same coil engineering requirements (4.3f)–(4.3i) as the ones used in Wechsung et al. (Reference Wechsung, Landreman, Giuliani, Cerfon and Stadler2022b). Each of the four base modular coils is described by a Fourier representation with 16 Fourier modes, which results in an optimization problem with 400 degrees of freedom corresponding to the geometry and coil currents.

Unlike in Landreman & Paul (Reference Landreman and Paul2022), we do not directly fix the aspect ratio of the outermost surface. Instead, we fix the major radius of the innermost surface and the volume of the outermost surface. We empirically observe that this effectively limits the variation of the aspect ratio of the outermost surface over the course of the optimization. Once the coils are optimized on five surfaces, we increase the number of BoozerExact surfaces to nine by keeping the original five surfaces and introducing four new surfaces in-between, and again run at most 45 000 BFGS iterations.

In the following description of the results, we will refer to the stage I precise QA magnetic field of Landreman & Paul (Reference Landreman and Paul2022) as QA-I, the configurations obtained from the classical stage II coil optimization in Wechsung et al. (Reference Wechsung, Landreman, Giuliani, Cerfon and Stadler2022b) for different maximum coil lengths as QA-II[$L_{\max }$], and the stage III coils obtained from solving the optimization problem described in this manuscript as QA-III[$L_{\max }$].

When optimizing on nine surfaces, all coil sets considered in this manuscript reached optimality, ($\| \boldsymbol {\nabla } J \|_{\infty } < 10^{-10}$), within the maximum number of iterations and the gradient was reduced by at least six orders of magnitude. More iterations are required for longer coil lengths since the optimization problem becomes more ill-conditioned (see figure 4). In figure 5, we provide two similar measures of quasisymmetry: the norm of the non-quasisymmetric component of the field magnitude and the maximum symmetry-breaking Fourier modes. Significant improvements in quasisymmetry are observed after optimization on five surfaces. Additional minor improvement is observed when increasing from five to nine surfaces (figure 5a). This suggests that for this particular set-up, five surfaces were sufficient to optimize for QA throughout the volume and not just in the neighbourhood of the surfaces.

Figure 4.

Norm of gradient in each iteration where five surfaces are used. Gradients at the initial configuration were $O(10^{-4})$, but we truncate the $y$-axis in order to better observe convergence to optimality. Every 3000 iterations (indicated by dotted vertical lines), the penalty weights are increased until the (in)equality constraints in (4.3) are satisfied to sufficient precision. Since a suitable set of initial weights were chosen, this does not occur frequently.

Figure 5.

Two measures of the field's departure from quasisymmetry: the ratio $\|B_{\textrm {non-QA}}\|/\sqrt {\textrm {area}}$ (a), and the maximum symmetry-breaking Fourier mode (b). In both plots, the magnetic field is normalized to be 1 T on the magnetic axis. The curves corresponding to the QA-III coils in panel (a) have two line styles, dotted ($\cdot \cdot \cdot \cdot \cdot$) and solid (—–), which correspond to configurations resulting from optimization on five and nine surfaces, respectively. The maximum symmetry-breaking Fourier mode is calculated by computing a two-dimensional fast Fourier transform of the field strength, and finding the maximum $|B_{m,n}|$ for $n\neq 0$, corresponding to the Fourier modes introducing variation in $\varphi$. The thick and thin vertical black lines correspond, respectively, to the five surfaces used in the initial optimization and the four surfaces added in the nine surface optimization run. The shaded grey region of the plot corresponds to field strengths below $50\,\mathrm {\mu }$T, corresponding to the Earth's background magnetic field.

In our implementation, we use MPI-based (message-passing-interface-based) parallelism to launch $N_s$ tasks that concurrently compute surfaces and adjoint variables, meaning that the time taken to evaluate the objective and its gradient does not substantially increase with the number of surfaces on which quasisymmetry is optimized. Since five surfaces might not be sufficient in all configurations, the possibility of parallelizing the surface and adjoint solves, as available in our implementation, means that we do not have to compromise on the number of surfaces used. Our implementation relies heavily on the open-source stellarator optimization package SIMSOPT (Landreman et al. Reference Landreman, Medasani, Wechsung, Giuliani, Jorge and Zhu2021) where many aspects of this work are contributed. We also use SIMSOPT's ‘optimizable’ framework for efficient representation of our objective (Appendix C). In the above example, a typical iteration of BFGS takes approximately two seconds on Intel Xeon Platinum 8268 processors when optimizing on five or nine surfaces.

In figure 5(b), we provide the magnitude of the maximum symmetry-breaking mode of the field magnitude and compare with eight previous configurations taken from figure 6(a) in Landreman & Paul (Reference Landreman and Paul2022). The shaded area of the figure indicates field magnitudes that are smaller than the Earth's background magnetic field ($50\,\mathrm {\mu }$T). We highlight that after optimization, the symmetry-breaking errors of the QA-III[24] coils improve not only on the QA-II[24] field, but even on the original QA-I field. This might have a number of explanations. For example, the dimension of the optimization problem here is 400 and thus just over three times larger than the dimension of the problem used to design QA-I, which is 120. This also might be because we use discretely exact gradients in this work, and all coils in the nine surface optimization reached optimality, to very close precision. This is in contrast to the QA-I field computations, where only finite difference gradients were used, thus potentially limiting the accuracy of gradients and as a result how much the objective could be decreased.

In the physics results discussed below, we now only analyse the coils optimized for QA using nine surfaces. Geometric properties of the final coils after optimization for quasisymmetry on nine surfaces are provided in table 1 and the coils are shown in figure 6. The full coil set for the QA-III[24] configuration is shown in figure 8(a). As $L_{\max }$ increases, the coil-to-surface separation increases as coils can move farther away from the outermost surface, thereby reducing the effects of the discrete nature of the coils and producing better quasisymmetry in the volume. As the coil length increases, the maximum curvatures, mean-squared curvatures, and the coil–coil separation approach and attain $\kappa _{\max }$, $\kappa _{\textrm {msc}}$, $d_{\min }$, respectively. We observe that although the surfaces do not change substantially in the optimization (figure 9), the coils do (figure 6).

Figure 6.

From (a–d), inboard view of the 18, 20, 22, 24 m coils before and after optimization. Grey and gold correspond, respectively,to the stage II coils and stage III coils (nine surface optimization). The colour on the surface corresponds to the magnetic field strength.

Table 1.

Comparison of geometric coil properties obtained from stage II and stage III optimization for different modular coil lengths. Stage III results are optimized using nine surfaces. Coil-to-surface distance is computed with respect to the surface used for coil design in the stage II rows, and the outermost surface on which quasisymmetry is optimized for the stage III rows.

Cross-sections of surfaces in the magnetic field after stage III optimization on nine surfaces and Poincaré plots are provided in figure 9(a), where there are no visible islands. The rotational transform is approximately constant and equal to $\bar {\iota }$ throughout the volume on which we are optimizing, varying between $0.415$ and $0.428$ for the different coil lengths. We also compare cross-sections of surfaces in the stage II and stage III fields in figure 9(b). We observe that although the coils change substantially (figure 6), the surfaces do not.

We compute two measures to quantify the confinement properties of our optimized magnetic configurations: alpha particle losses (figure 7) and the thermal collisional transport magnitude $\epsilon _{\textrm {eff}}^{3/2}$ (figure 8b). The alpha particle losses are calculated as follows. Each configuration is first scaled to match the minor radius and average field strength of the ARIES-CS reactor, 1.70 m and 5.9 T. Alpha particles are initialized isotropically in velocity space with 3.5 MeV energy, and uniformly on the flux surface with toroidal flux $1/4$ that of the outermost computational surface, corresponding to half of the effective minor radius. The collisionless guiding-centre trajectories are followed using the SIMPLE code (Albert, Kasilov & Kernbichler Reference Albert, Kasilov and Kernbichler2020a,Reference Albert, Kasilov and Kernbichlerb), and particles are considered lost if they exit the outermost surface. When considering the confinement of alpha particles, we see the most significant improvements for the two shorter coil sets. For QA-III[18], losses are reduced from 17.7 % to 6.6 %, for QA-III[20], losses are reduced from 1.7 % to 0.7 %. For the longer QA-III[22] and QA-III[24] coils, the (small) particle losses are comparable to those for QA-II[22] and QA-II[24], respectively, despite significantly better quasisymmetry throughout the volume. We note that this imperfect correlation between quasisymmetry and energetic particle confinement is consistent with the findings in Bader et al. (Reference Bader, Anderson, Drevlak, Faber, Hegna, Henneberg, Landreman, Schmitt, Suzuki and Ware2021), and in Landreman & Paul (Reference Landreman and Paul2022) where the QA+Well configuration had lower particle losses than the QA configuration despite worse quasisymmetry. For $\epsilon _{\textrm {eff}}^{3/2}$ (Nemov et al. Reference Nemov, Kasilov, Kernbichler and Heyn1999) we observe significant improvement compared with the coils of Wechsung et al. (Reference Wechsung, Landreman, Giuliani, Cerfon and Stadler2022b), and our procedure is even able to improve the results from the target field of Landreman & Paul (Reference Landreman and Paul2022).

Figure 7.

Losses of alpha particles spawned on the surface with normalized toroidal flux $0.25$. We observe significant improvements in confinement for the shorter coils. For $L_{\max }=18\,\textrm {m}$, losses are reduced from 17.7 % to 6.6 %, for $L_{\max }=20\,\textrm {m}$, losses are reduced from 1.7 % to 0.7 %. The stage II and III coils with lengths $L_{\max }=22, 24\,\textrm {m}$ present comparable particle losses, despite significantly better quasisymmetry for the stage III coil sets. The previous configurations, shown in grey, are detailed in Landreman & Paul (Reference Landreman and Paul2022).

Figure 8.

(a) The QA-III[24] stellarator and the outermost surface on which quasisymmetry is optimized. (b) The thermal collisional transport magnitude $\epsilon _{\textrm {eff}}^{3/2}$ is significantly reduced; the QA-III[20], QA-III[22] and QA-III[24] configurations even improve on the QA-I configuration. The previous configurations, shown in grey, are detailed in Landreman & Paul (Reference Landreman and Paul2022).

Figure 9.

Shown in (a) and (b) with solid lines (—–) are cross-sections of the nine surfaces used in the QA-III[18] coil set optimization at $\phi = 0, {\rm \pi}/4, {\rm \pi}/2$ in black, red and blue. Grey lines on the left correspond to Poincaré plots. Shown in (b)are cross-sections of surfaces in the QA-II[18] field using dashed lines (- - -). The longer coil sets result in comparable cross-sections, but the difference between stage II and III cross-sections become smaller.