New method to design stellarator coils without the winding surface

For convenience, we work within the filamentary approximation, by which we mean that each coil has zero cross-sectional area. (This approximation is not fundamental, and can easily be relaxed by adding multiple parallel curves.) We represent the geometry of the coils as closed, one-dimensional curves embedded in three-dimensional space.

The geometry of such curves is elegantly described by the fundamental theorem of curves [18], which states that any regular curves in three-dimensional space with non-zero curvature can be determined (up to rigid displacements and rotations) by the curvature, $\kappa(s)$, and torsion, $\tau(s)$, where s parameterizes arc-length. Given κ and τ, and boundary conditions (such as the starting point), the geometry, ${\bf x}(s)$, unit tangent vector, ${\bf t}(s)$, unit normal and bi-normal, ${\bf n}(s)$ and ${\bf b}(s)$, of the curve can be obtained by integration of the Frenet–Serret formulas [19],

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle \label{eq:curve} \left[\begin{array}{@{}c@{}} \tvec{x} \nonumber \\ \tvec{t} \nonumber \\ \tvec{n} \nonumber \\ \tvec{b} \nonumber \\ \end{array}\right]^\prime = \left[\begin{array}{@{}ccc@{}} 1 & 0 & 0 \nonumber \\ 0 & \kappa & 0 \nonumber \\ - \kappa & 0 & \tau \nonumber \\ 0 & \tau & 0 \nonumber \\ \end{array}\right] \cdot \left[\begin{array}{@{}c@{}} \tvec{t} \nonumber \\ \tvec{n} \nonumber \\ \tvec{b} \nonumber \\ \end{array}\right], \nonumber \end{align} \tag{ 1 }$$

where the 'prime' denotes derivative with respect to the parameter s, i.e. $\newcommand{\e}{{\rm e}} {\bf x}^\prime \equiv {\rm d}{\bf x} / {\rm d}s$.

Each curve is described by two independent functions of arc length, namely the curvature and torsion. Using the curvature function as an independent degree of freedom to implicitly define the coil geometry has its appeal, as the curvature of the coils is frequently an engineering constraint and any curvature constraints can potentially be enforced before the nonlinear optimization proceeds. However, solving 12 ODEs for each coil at each iteration of the nonlinear optimization turns out to be somewhat cumbersome. Furthermore, to exploit efficient nonlinear optimization algorithms, we shall later construct the derivatives of the coil geometry with respect to the independent degrees of freedom, and this requires integration of the derivatives of equation (1) with respect to the parameters describing the curvature and torsion, and this becomes even more cumbersome. Not all curvatures and torsion result in closed curves, and additional constraints must be included to ensure that the curves are both closed and smooth.

Instead, we employ a very easy way of describing one-dimensional curves embedded in three-dimensional space. A curve is described directly, and completely generally, in Cartesian coordinates as $\newcommand{\tvec}[1]{{\bf #1}} \tvec{x}(t) = x(t) \, {\bf i} + y(t) \, {\bf j} + z(t) \, {\bf k}$. Three functions are required to specify the geometry, namely $x(t)$, $y(t)$ and $z(t)$, with the constraints that each function be periodic, e.g. ${\bf x}(t+T) = {\bf x}(t)$ for some T. The curve parameter, t, at this point is arbitrary.

A variety of mathematical representations are possible. For purpose of illustration, and because our initial interest is in smooth coils, we use a Fourier representation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:fseries} x(t) = x_{c,0} + \sum_ {n=1}^{N_{\rm F}} \left [x_{c,n} \cos(nt) + x_{s,n} \sin(nt) \right ], \nonumber \end{align} \tag{ 2 }$$

with t varying between $[0, 2\pi]$, and similarly for $y(t)$ and $z(t)$. The shape of a coil is then fully determined by the $3 \times (2N_{\rm F}+1)$ Fourier coefficients. No additional assumptions are made here, so this representation fits for all kinds of closed smooth coils, such as helical, modular, saddle, etc. (For coils with straight segments, for example rectangular and so-called 'window pane' coils, the Fourier representation is not suitable. In future upgrades to FOCUS, additional representations for the curves, such as cubic splines and piecewise-linear, will be implemented.)

Hereafter, $\newcommand{\tvec}[1]{{\bf #1}} \tvec{X}$ will be used to represent the geometrical degrees of freedom of a set of N_C coils, ${\bf x}_i$, and the coil currents, I_i, which are also treated as independent degrees of freedom. The total number of degrees of freedom is $N_{D} = N_C \times [ 3 \times (2N_{\rm F}+1) + 1 ]$.

The coil parameters are to be varied to minimize a target function consisting of both 'physics' and 'engineering' objective functions,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle \label{eq:target} \chi^2(\tvec{X}) = \sum_j w_j \left( \frac{f_j({\tvec X}) - f_{j,o}}{f_{j,o}} \right)^2, \nonumber \end{align} \tag{ 3 }$$

where $\newcommand{\tvec}[1]{{\bf #1}} f_j({\tvec X})$ is the value of the $j{\rm th}$ objective function, to be defined below, for a given set of coil parameters, $f_{j, o}$ denotes the desired value, and w_j is a user-prescribed weight.

The key criterion is to satisfy the condition of $\newcommand{\tvec}[1]{{\bf #1}} B_n = \tvec{B} \cdot \tvec{n} = 0$ on the pre-defined plasma boundary, as much as is possible with a finite number of coils. A smaller B_n error usually implies a better approximation to the target equilibrium properties (but so-called 'resonant' normal field errors may be particularly problematic). This can be achieved by minimizing the squared total normal field integrated over the boundary,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}[1]{{\rm d} #1} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle \label{eq:bnormal} f_B(\tvec{X}) \equiv \int_S \frac{1}{2} \left ( \tvec{B} \cdot \tvec{n} \right)^2 \dd{s}. \nonumber \end{align} \tag{ 4 }$$

The magnetic field at position $\boldsymbol {\bar{\bf x}}$ produced by a coils set is calculated using the Biot–Savart law,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}[1]{{\rm d} #1} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle \label{eq:biotsavart} \tvec{B}_V(\boldsymbol {\bar{\bf x}}) = \frac{\mu_0}{4\pi}\sum_{i=1}^{N_C} I_i \ \int_{C_i} \frac{\dd{\tvec{l}_i} \times \tvec{r} }{r^3}. \nonumber \end{align} \tag{ 5 }$$

Here, I_i is the current in $i{\rm th}$ coil, $\newcommand{\tvec}[1]{{\bf #1}} \newcommand{\dd}[1]{{\rm d} #1} \newcommand{\e}{{\rm e}} \dd{\tvec{l}_i} \equiv {\bf x}_i^\prime \, {\rm d}t$ is the differential line element, and $\newcommand{\tvec}[1]{{\bf #1}} \tvec{r} = \boldsymbol {\bar {\bf x}}- \tvec{x}_i$ is the displacement vector between the evaluate point on the surface and the differential element.

The plasma field, $\newcommand{\tvec}[1]{{\bf #1}} \tvec{B}_{P}$, produced by the reference plasma must be pre-computed using a suitable equilibrium code, and the following will assume that the plasma field does not change during the optimization of the coil geometry. The FOCUS code can accommodate an arbitrary $\newcommand{\tvec}[1]{{\bf #1}} \tvec{B}_{P}$, but for simplicity of illustration in this paper we restrict attention to vacuum fields, so that $\newcommand{\tvec}[1]{{\bf #1}} \tvec{B} = \tvec{B}_{V}$, and we hereafter suppress the subscript V for notational clarity.

If a small geometrical deformation, $\newcommand{\tvec}[1]{{\bf #1}} \delta \tvec{x}_i$, is applied to the $i{\rm th}$ coil, the corresponding changes in $\newcommand{\tvec}[1]{{\bf #1}} \tvec{r}$ and $\newcommand{\e}{{\rm e}} r\equiv \vert {\bf r}\vert$ are $\newcommand{\tvec}[1]{{\bf #1}} \delta \tvec{r} = - \delta \tvec{x}_i$ and $\newcommand{\tvec}[1]{{\bf #1}} \delta r = - \tvec{r} \cdot {\delta \tvec{x}_i}/r$, and the change in the magnetic field is

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}[1]{{\rm d} #1} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle \delta \tvec{B}({\boldsymbol {\bar {\bf x}}}) = \frac{\mu_0}{4\pi} I_i \int_0^{2\pi} \left [ \frac{3\tvec{r} \cdot \tvec{x}_i^\prime}{r^5} \tvec{r}\times \delta \tvec{x}_i \, + \, \frac{2}{r^3}\delta \tvec{x}_i \times \tvec{x}_i^\prime \, + \, \frac{3\tvec{r} \cdot \delta \tvec{x}_i}{r^5} \tvec{x}_i^\prime \times \tvec{r} \right ] \dd{t}. \nonumber \end{align} \tag{ 6 }$$

The change in f_B is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}[1]{{\rm d} #1} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle \label{eq:Obderiv} \delta f_B(\tvec{X}) =\,& \int_S (\tvec{B} \cdot \tvec{n}) \, \frac{\mu_0}{4\pi} I_i \int_0^{2\pi} \left [ \frac{3\tvec{r} \cdot \tvec{x}_i^\prime}{r^5} \tvec{n} \times \tvec{r} \, \right.\nonumber\\&\left.+ \, \frac{2}{r^3} \tvec{x}_i^\prime \times \tvec{n} \, + \, \frac{3 \tvec{x}_i^\prime \times \tvec{r} \cdot \tvec{n}}{r^5} \tvec{r}\right ] \cdot \delta \tvec{x}_i \; \dd{t} \; \dd{s}. \nonumber \end{align} \tag{ 7 }$$

The variation of f_B with respect to variations in the coil currents $\delta I_i$ is similar (and somewhat simpler).

To avoid trivial solutions, for example $f_B \rightarrow 0$ when $I_i \rightarrow 0$, $\forall i$, it is sufficient to constrain the enclosed toroidal flux. If ${\bf B}\cdot{\bf n}=0$ on the boundary, then the toroidal flux through any poloidal cross-sectional surface is constant. We include an objective function defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}[1]{{\rm d} #1} \newcommand{\tvec}[1]{{\bf #1}} \newcommand{\z}{\zeta} \displaystyle f_\Psi(\tvec{X}) \equiv \frac{1}{2\pi} \int_0^{2\pi} \frac{1}{2} \left( \frac{\Psi_\z \; - \; \Psi_o}{\Psi_o} \right)^2 \dd{\z}, \nonumber \end{align} \tag{ 8 }$$

where the flux through a poloidal surface, ${{\mathcal T}}$, produced by cutting the boundary with plane $\newcommand{\z}{\zeta} \z=const.$ is computed using Stokes' theorem,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}[1]{{\rm d} #1} \newcommand{\tvec}[1]{{\bf #1}} \newcommand{\z}{\zeta} \displaystyle \Psi_\z({\bf X}) \equiv \int_{{\mathcal T}} \! \tvec{B} \cdot \dd{\tvec{S}} = \oint_{\partial {{\mathcal T}}} \!\! \tvec{A} \cdot \dd{\tvec{l}}. \nonumber \end{align} \tag{ 9 }$$

Here $\newcommand{\tvec}[1]{{\bf #1}} \newcommand{\dd}[1]{{\rm d} #1} \dd \tvec{l}$ is a differential line segment on the boundary curve of the poloidal surface, and the magnetic vector potential, $\newcommand{\tvec}[1]{{\bf #1}} \tvec{A}$, is

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}[1]{{\rm d} #1} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle \tvec{A}({\boldsymbol {\bar {\bf x}}}) = \frac{\mu_0}{4\pi}\sum_{i=1}^{N_C} I_i \ \int_{C_i} \frac{\dd{\tvec{l}_i}}{r}. \nonumber \end{align} \tag{ 10 }$$

The variation of $f_\Psi$ resulting from $\newcommand{\tvec}[1]{{\bf #1}} \delta \tvec{x}_i$ is

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}[1]{{\rm d} #1} \newcommand{\tvec}[1]{{\bf #1}} \newcommand{\z}{\zeta} \displaystyle \label{eq:Otderiv} \delta f_\Psi(\tvec{X}) = \frac{1}{2\pi} \int_0^{2\pi} \left (\frac{\Psi_\z \; - \; \Psi_o}{\Psi_o} \right) \ \frac{\delta \Psi_\z}{\Psi_o}\dd{\z}, \nonumber \end{align} \tag{ 11 }$$

where $\newcommand{\z}{\zeta} \newcommand{\tvec}[1]{{\bf #1}} \newcommand{\dd}[1]{{\rm d} #1} \delta \Psi_\z = \int_{\partial {{\mathcal T}}} \delta \tvec{A} \cdot \dd{\tvec{l}}$ and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}[1]{{\rm d} #1} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle \delta \tvec{A}({\boldsymbol {\bar {\bf x}}}) = \frac{\mu_0}{4\pi} I_i \ \int_0^{2\pi} \left [ -\frac{\tvec{r} \cdot \tvec{x}_i^\prime}{r^3} \ \delta \tvec{x}_i \, + \, \frac{\tvec{r} \cdot \delta \tvec{x}_i}{r^3} \ \tvec{x}_i^\prime \right ] \dd{t}. \nonumber \end{align} \tag{ 12 }$$

In addition to the physical constraints, namely that the coils produce the required magnetic field, engineering constraints should also be considered, such as the coil–coil separation, coil-plasma separation, minimum radius of curvature, electromagnetic forces, stored energy, etc. The engineering constraints may or may not compete with the physics constraints. For example, the 'ripple' in ${\bf B}\cdot{\bf n}$ on the reference boundary, for a fixed number of coils, is reduced as the coils are moved further away from the plasma. This also increases the coil-plasma separation, thus increasing the space available for divertors etc. Similarly, the ripple in ${\bf B}\cdot{\bf n}$ is reduced by positioning the coils approximately equally spaced in toroidal angle, and this increases the coil–coil separation. These engineering constraints are compatible with the physics constraint of minimizing ${\bf B}\cdot{\bf n}$ on the reference boundary. In fact, in the following we will show that we can obtain acceptable coil sets without including the engineering constraints on coil-plasma and coil–coil separation.

There are also engineering constraints that do compete with the physics, most notably that the total number of discrete coils is less than infinity and the constraint that the length of the coils be reasonable. Without a constraint on the lengths of each coil, we found that the coils can become both arbitrarily long and can develop more 'wiggles' so as to better produce the required magnetic field. We found that it was required to include a penalty on the coil length. This penalizes coils with large curvature and leads to shorter coils. This is attractive, as to a good approximation the cost of building a coil set is proportional to the total length of the coils.

We include an objective function of the form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle f_L(\tvec{X}) = \frac{1}{N_C} \; \sum_{i=1}^{N_C} \frac{{\rm e}^{L_i}}{{\rm e}^{L_{i,o}}}, \label{eq:lengthpenalty} \nonumber \end{align} \tag{ 13 }$$

where $L_i({\bf X})$ is the length of ith coil,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\abs}[1]{\vert#1\vert} \newcommand{\dd}[1]{{\rm d} #1} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle L_i(\tvec{X}) = \int_0^{2\pi} \!\! \abs{\tvec{x}_i^\prime} \; \dd{t}, \nonumber \end{align} \tag{ 14 }$$

and $L_{i, 0}$ is a user-specified normalization. The variations in f_L resulting from a variation $\newcommand{\tvec}[1]{{\bf #1}} \delta \tvec{x}_i$ is

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}[1]{{\rm d} #1} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle \label{eq:Olderiv} \delta f_L(\tvec{X}) = \frac{1}{N_C} \; \frac{{\rm e}^{L_i}}{{\rm e}^{L_{i,0}}} \int_0^{2\pi} \frac{({\tvec{x}^\prime_i} \cdot {\tvec{x}^{\prime\prime}_i}){\tvec{x}^\prime_i} - ({\tvec{x}^\prime_i} \cdot {\tvec{x}^\prime_i}){\tvec{x}^{\prime\prime}_i}}{({\tvec{x}_i^\prime} \cdot {\tvec{x}_i^\prime})^{3/2}} \cdot \delta \tvec{x}_i \ \dd{t}. \nonumber \end{align} \tag{ 15 }$$

The length penalty, equation (13), will encourage the coils to shorter. This can be replaced with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:Lquad} f_L({\bf X}) = \frac{1}{N_C} \; \sum_{i=1}^{N_C} \frac{1}{2} \frac{(L_i - L{i,o})^2}{L_{i,o}^2}, \nonumber \end{align} \tag{ 16 }$$

to make each coil length close to a target value.

The parameterization given in equation (2) is not unique. Consider the curve variation $\delta {\bf x}_i = {\bf x}^\prime_i(t) \, \delta u$. Inserting this into equations (7), (11) and (15), we can easily get that $\delta f_B$, $\delta f_\Psi$ and $\delta f_L$ are all zero. This 'variation' is tangential to the curve and leaves the geometry of curve unchanged and thus leaves the magnetic field unchanged, but it does effectively change the parameterization of the curve and thereby changes the Fourier harmonics.

This non-uniqueness of the curve parameters can be exploited for numerical efficiency. For each coil, a spectral width of the Fourier harmonics is defined as [20]

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle M_i(\tvec{X}) = \sum_ {n=1}^{N_{\rm F}} n^{\,p} ( {X^{i}_{c,n}}^{2} + {X^{i}_{s,n}}^{2} + {Y^{i}_{c,n}}^{2} + {Y^{i}_{s,n}}^{2} + {Z^{i}_{c,n}}^{2} + {Z^{i}_{s,n}}^{2}), \nonumber \end{align} \tag{ 17 }$$

where $p \geqslant 1$. Curve parameterizations that minimize M_i are attractive for numerical purposes, because curves with complicated geometry can be described with fewer Fourier harmonics. Minimizing the spectral width, which is a purely 'numerical' constraint, should not compete with the physics or engineering properties. So when minimizing M, we only allow tangential variations in the geometry of coils. In the Fourier parameter space, the tangential variation is defined as $\newcommand{\e}{{\rm e}} \delta {\bf X} \equiv {\bf X}^\prime \delta u$, where ${\bf X}^\prime$ is the tangential direction in parameter space and $\delta u$ is an arbitrary function. If we choose $\newcommand{\e}{{\rm e}} \delta u \equiv - \partial_{\bf X} M \cdot {\bf X}^\prime$, the spectral width will be forced to decrease subject without changing the coil geometry. This shall be used below.

Our target function is the weighted sum of the normal field error, the toroidal flux error, and the length penalty, and is represented as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle \chi^2(\tvec{X}) = w_B \;f_B(\tvec{X}) + w_{\Psi} \;f_\Psi(\tvec{X}) + w_L \; f_L(\tvec{X}). \label{eq:minimalobjectivefunction} \nonumber \end{align} \tag{ 18 }$$

The inclusion of each of these objective functions has been justified in the above discussion: reducing the normal field error is required so that the correct magnetic bottle is produced, the toroidal flux constraint is required to avoid trivial solutions, and the length penalty is required to avoid infinitely long coils.

What is most remarkable is that geometry of the coils is completely free and the only so-called engineering constraint is a penalty on the length, and perhaps the implicit constraint that the number of coils is finite and fixed. The results to be shown below support the hypothesis that this is the minimal objective function. Other constraints, such as including a penalty on the coil–coil or coil-plasma separation can be included but they are not obligatory. We call equation (18) the minimal objective function as not one of f_B, $f_\Psi$ or f_L can be removed, except perhaps for the qualification that if trivial solutions could be guaranteed to be avoided (by not allowing the currents to vary, for example) then perhaps including the toroidal flux constraint may be removed.

The task of designing optimal coils is now reduced to a single-object minimization problem. A variety of minimization algorithms can be applied, such as the Levenberg-Marquardt (LM) method [21], which was successfully used in COILOPT and COILOPT++. Unlike other nonlinear coil-designing codes, where the derivatives are computed numerically, we construct the analytically calculated derivatives from equations (7), (11) and (15).

Calculating the objective functions and their derivatives involves computing line and surface integrals. Each coil is treated as a single filament with $N_{s}$ segments. The predefined plasma surface is discretized into $\newcommand{\z}{\zeta} N_{\theta} \times N_{\z}$ surface elements. We shall illustrate how we implement the numerical calculations in the code by one example. The derivative of the normal field error f_B with respect to an arbitrary cosine harmonic of the $i{\rm th}$ coil denoted as $x_{c, n}^i$ is computed as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\pdv}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\dd}[1]{{\rm d} #1} \newcommand{\tvec}[1]{{\bf #1}} \newcommand{\z}{\zeta} \displaystyle \pdv{f_B}{x_{c,n}^i} & = \int_S (\tvec{B} \cdot \tvec{n}) \; \pdv{\tvec{B} \cdot \tvec{n}}{x_{c,n}^i} \dd{s} \nonumber \\ & \approx \sum_{j=1}^{N_\theta} \sum_{k=1}^{N_\z} {\bigg [} (B_x n_x + B_y n_y + B_z n_z) \nonumber \\ & {\bigg (}\pdv{B_x}{x_{c,n}^i} n_x + \pdv{B_y}{x_{c,n}^i} n_y + \pdv{B_z}{x_{c,n}^i} n_z {\bigg)} \sqrt g \; {\bigg ]}_{j+\frac{1}{2},k+\frac{1}{2}} \frac{2\pi}{N_{\theta}} \frac{2\pi}{N_{\z}}. \nonumber \end{align} \tag{ 19 }$$

The integral is evaluated by the mid-point rule. The calculation of $\partial{B_x} / \partial{x_{c, n}^i}$, which involves a line integral, can be approximated

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\pdv}[2]{\frac{\partial #1}{\partial #2}} \displaystyle \pdv{B_x}{x_{c,n}^i} \approx \frac{\mu_0}{4\pi} I_i \ \sum_{m=1}^{N_s} \left . - \frac{3 (\Delta y \ {z^\prime_t} - \Delta z \ {y^\prime_t}) \ \Delta x \ \cos(nt)}{r^5} \right \vert _{m} \frac{2\pi}{N_{s}}. \nonumber \end{align} \tag{ 20 }$$

The derivatives of toroidal-flux error and the coil-length penalty are calculated similarly.

Analytically calculated derivatives offer several advantages in both accelerating the computation speed and accessing alternative algorithms. For the proof of concept of the new coil representation and the cost function employed in FOCUS, we use the very simple, steepest descent algorithm. Defining an artificial 'time', τ, the descent direction is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\pdv}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle \label{eq:descent} \pdv{\tvec{X}}{\tau} \equiv - \pdv{\chi^2}{\bf X} - {\bf X}^\prime \; \left( \pdv{M}{\bf X} \cdot {\bf X}^\prime \right). \nonumber \end{align} \tag{ 21 }$$

The first term in the right hand side is the usual gradient descent, and the second is in the tangential direction and is included for the spectral condensation. The relationship between the target function and time is

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\pdv}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle \pdv{\chi^2}{\tau} =\,& \pdv{\chi^2}{\tvec{X}} \cdot \pdv{\tvec{X}}{\tau} = - \pdv{\chi^2}{\tvec{X}} \cdot \pdv{\chi^2}{\bf X} \; \nonumber\\&- \; \left ( \pdv{\chi^2}{\tvec{X}} \cdot {\bf X}^\prime \right) \; \left( \pdv{M}{\bf X} \cdot {\bf X}^\prime \right) = - \left (\pdv{\chi^2}{\tvec{X}} \right)^2, \nonumber \end{align} \tag{ 22 }$$

where we have used $\partial_{\bf X} \chi^2 \cdot {\bf X}^\prime = 0$, namely that tangential variations do not change the target function. The first derivative, $\partial \chi^2 / \partial \tau$, is always non-positive. After reaching the 'geometrical' minimum, $\partial \chi^2 / \partial {\bf X} = 0$, the inclusion of the tangential term in equation (21) advances the Fourier harmonics towards the 'spectral minimum',

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\pdv}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle \pdv{M}{\tau} = \pdv{M}{\tvec{X}} \cdot \pdv{\tvec{X}}{\tau} = - \left ( \pdv{M}{\bf X} \cdot {\bf X}^\prime \right)^2. \nonumber \end{align} \tag{ 23 }$$

The ODEs in equation (21) are integrated using a fixed order Runge–Kutta method with the subroutine D02BJF in the NAG Library [22]. We emphasize that this article does not attempt to illustrate the attractiveness of the steepest descent method but rather to illustrate that the objective function described in equation (18) is sufficient to define suitable coils.

The steepest descent method is continuous. This means that the coils cannot pass through the plasma, as this would result in a singularity in the normal field. We can define the linking number of each coil and the plasma using the Gauss linking integral [23] with respect to the coil and the magnetic axis:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dd}[1]{{\rm d} #1} \newcommand{\tvec}[1]{{\bf #1}} \displaystyle {\rm link}(C_i, {\rm Plasma}) = \frac{1}{4\pi} \ \oint_{C_i} \oint_{\rm axis} \frac{\tvec{r}_1 - \tvec{r}_2}{\vert \tvec{r}_1 - \tvec{r}_2\vert ^3} \cdot \dd{\tvec{r}_1} \times \dd{\tvec{r}_2}. \nonumber \end{align} \tag{ 24 }$$

With the steepest descent algorithm, this linking number is implicitly conserved during the optimization. This amounts to an additional, implied constraint.

A two-field period, rotating elliptical plasma boundary is specified as

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\z}{\zeta} \displaystyle R = & 3.0 + 0.3 \cos(\theta) \quad \quad \quad \quad \quad \quad \quad \quad - 0.06 \cos(\theta - N_P \z), \nonumber \\ Z = & \quad \ - 0.3 \sin(\theta) - 0.06 \sin(- N_P \z) - 0.06 \sin(\theta - N_P \z), \nonumber \end{align*}$$

where $N_P = 2$ and $\newcommand{\z}{\zeta} \z$ is the geometrical toroidal angle. As one of the simplest configurations for stellarators, this rotating elliptical boundary can be produced with conventional helical coils, such as in LHD [24]. Helical coils, with certain pitch angles, can efficiently produce desired rotational transforms, particularly for classical stellarators with helix structures. However, the access to the plasma and ease of maintenance can be badly restricted. Modular coils, on the other hand, can overcome these difficulties, but differences in the local structures of the fields produced by helical coils and modular coils will exist [25]. For illustration, we use FOCUS to construct modular coils for this boundary.

To initialize the coil-design calculation in FOCUS using the steepest-descent method, we must supply a suitable initial guess for the coil geometry. We anticipate that, generally, codes such as NESCOIL and REGCOIL will be used for this, but here a very simple initial guess for the coils is made by placing 16 circular coils equally spaced in the toroidal angle, as shown in figure 1(a). The radius of the coils is $R_C=0.75$ m, and the chosen Fourier resolution is $N_{\rm F}=4$.

Zoom In Zoom Out Reset image size

To test the different objective functions separately, we first include only the physical constraints, i.e. we set the length penalty weight to zero, $w_L=0$. The limit on the time integration is set to an arbitrary value, e.g. $\tau_{\rm end}=10^3$. As illustrated in figure 2, f_B decreases from the initial value of $3.63\times10^{-2}$ to $1.26\times10^{-5}$ at $\tau = 10^3$, a significant reduction. The final shape of the optimized coils are shown in figure 1(b). We can see that all the coils are getting longer and further away from the plasma to better reproduce the given boundary and to reduce the ripple. A Poincaré plot of the flux surfaces produced by the optimized coils, figure 3, indicates that a satisfactory magnetic field has been obtained.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this optimization, only physics targets were included in the total target function, i.e. $\chi^2=$ $w_B \;f_B + w_{\Psi} \;f_\Psi$, with $w_B / w_{\Psi}=100/1$. Engineering constraints such as coil–coil separation, coil-plasma separation, even the coil length penalty were not included. The coils are not constrained to lie on a winding surface and are completely free to move in space. The fact that we get reasonable coils is remarkable.

Figure 2 also indicates that the target function is not converged at $\tau = 10^3$. Even at $\tau = 10^6$, the target function is still not converged. We do not give any physical meaning to the artificial time parameter; the only salient point here is that the coil geometry seems to not converge as $\tau \rightarrow \infty$ for this case where there is no penalty on the length. We now explore the convergence properties of FOCUS by varying different numerical parameters.

Increasing the number of Fourier harmonics, N_F, will presumably allow access to a greater variety of coil geometries and in turn decrease the normal field error. As shown in figure 4(a), when the number of Fourier harmonics is restricted to $N_{\rm F}=1$, which means that the coils can only be planar, the normal field error cannot decrease as much as compared to the cases with larger N_F. For $N_{\rm F}=1$, the optimized coils, as illustrated in figure 6(a), are elliptical and rotate to follow the pattern of the plasma shape. For the $N_{\rm F}\geqslant 4$ cases, there are no distinct differences between the coils, showing that the coil geometry has converged with respect to the Fourier resolution. In other words, for this elliptical configuration, the coils can be adequately represented by $N_{\rm F} = 4$.

Zoom In Zoom Out Reset image size

With more coils, the more accurately a given reference boundary can be recovered. However, more coils imply higher cost and less separation between each coil, which is vital for diagnostic ports. In practice, a compromise between accuracy and space must be found. Figure 4(b) shows the convergence with respect to N_C. As expected, with more coils lower f_B has been achieved. The toroidal Fourier harmonics of the magnetic field magnitude, $\newcommand{\z}{\zeta} B(\z) = \int_{0}^{2\pi} \vert B\vert \, {\rm d}{\theta} / 2\pi$, on the plasma boundary are shown in figure 5. Typically there will be a local maximum in the Fourier amplitude (i.e. the ripple) that corresponds with the number of coils. This implies that the modular ripples are reduced with more coils. Without sufficient coils, the shape of the optimized coils would be too complicated, like the coils shown in figure 6(b). The four coils have large twists to produce the desired magnetic field.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the calculations shown in figures 4(a) and (b), the initial coils are circular coils with an arbitrarily chosen radius of $R_C = 0.75$ m. For different initial radii, differing in a small range as shown in figure 4(c), we found that the final B_n errors tend to be identical. The optimized coils have similar geometries (details not shown here), which indicates that they are at the same minimum. Thus, the optimizing process is not very sensitive to the initial radius, at least not for this equilibrium.

Now we explore the effect of varying the length penalty weight. We assign w_L with different values ranging from $10^{-6}$ to $1.0$. The results in figure 4(d) show that increasing the length penalty has a deleterious impact on the normal field error. And, the optimization converges faster (with respect to τ) as the length weight increases. This can be used to control the convergence of the code, and to make a balance between the physical accuracy and engineering feasibility.

By choosing suitable weights, a suitable coils set for the rotating ellipse configuration is obtained. Choosing the weights such that, for the initial configuration, we have $w_Bf_B = 1.0$, $w_{\Psi}f_\Psi = 0.5$ and $w_L \,f_L = 0.2$, and with 16 circular coils with an initial radius of $R_C=0.75$ m, we get an attractive coils set at $\tau=10^3$, as shown in figure 6(c).

We now consider a more complicated geometry, that of W7-X. The coils system of W7-X mainly consists 50 non-planar modular coils, which produce the dominant component of the confining magnetic field, and 20 planar modular coils, which are primarily used for flexible control of the rotational transform. The coils are arranged in five field periods [26]. The last closed flux surface (LCFS) in the 'standard configuration' [27] is used here as the target plasma boundary. For this configuration, only the modular coils are used.

With the given plasma boundary, we start the optimizations with a set of circular coils that are positioned surrounding the plasma surface at equal toroidal intervals, i.e. the $i{\rm th}$ coil is placed in the plane $\newcommand{\z}{\zeta} \zeta = 2\pi(i-1)/{N_C}$. The first calculation uses 50 coils ($R_{\rm coil}=1.25$ m). The normal field error, toroidal flux, $\Psi_o = 2.133$ Wb, and length penalty, $L_o = 8.5$ m, and using the quadratic form equation (16), are all included, with artificially chosen weights. By $\tau=5.0$, the normal field error f_B has been reduced from $1.44\times10^1$ to $2.95\times10^{-4}$, while $f_{\Psi}$ from $1.07\times 10^{-2}$ to $3.70\times 10^{-8}$ and f_L from $2.89\times 10^{-3}$ to $3.26\times 10^{-7}$. The optimized coils are shown in figure 7. Poincaré plots of the produced magnetic field, as shown in figure 8, indicate that the target boundary is reproduced very well.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An interesting exploration is to reduce the number of coils N_C, while keeping all the other parameters unchanged. The Poincaré plots of FOCUS optimized coils with different N_C are shown in figure 9. It's encouraging that, even with as less as 30 coils, FOCUS can approximate the target boundary very well. The coils shapes, as shown in figure 10 for $N_C=30$, appear 'reasonable', by which we mean that the coils are sufficiently smooth, sufficiently well-separated, etc. We understand that with less coils the magnetic ripple would presumably increase, as we studied in the rotating elliptical case. We emphasize that this paper is primarily concerned with the introduction of a new coil optimization algorithm, and the W7-X calculation is used for an initial illustration. We have not performed a comprehensive experimental design study of W7-X, and that we have not estimated the effect that the small changes in the magnetic field have on plasma performance, and we have not estimated the coil–coil electromagnetic forces, for example. The results show that FOCUS is able to design coils for stellarators as complicated as W7-X with a high precision.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Besides the W7-X device, another existing large stellarator/heliotron machine is the large helical device (LHD) [28] in Japan. Unlike the modular coils used in W7-X, the main coils of LHD are a pair of helical windings. The approximated LHD coils, including two helical windings and three pairs of vertical field coils, are shown in figure 11.

Zoom In Zoom Out Reset image size

Helical curves on a torus can be represented by the Fourier series as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle x & = [ R + r \cos(\theta) ] \cos(N\theta),\nonumber \\ y & = [ R + r \cos(\theta) ] \sin(N\theta), \nonumber \\ z & = r \sin(\theta), \nonumber \end{align} \tag{ 25 }$$

where R, r are the major and minor radius of the torus, and $N=10/2$ for the LHD case.

When fitting the actual LHD coils with an $N_{\rm F}=6$ Fourier series, truncation errors are introduced. Consequently, the fitted coils have slightly different geometries from the actual coils, which ultimately results in the relatively large errors for the produced magnetic field ($f_B= 1.062\times 10^{-1}$). Figure 12(a) shows that the initial fitted coils have deteriorated magnetic flux surfaces. There are no closed flux surfaces at the target boundary and the magnetic axis is significantly shifted outwards. When we use FOCUS to optimize the fitted coils, the ${\bf B}\cdot{\bf n}$ error f_B is decreased to $7.112\times 10^{-3}$ at $\tau=10^2$. The Poincaré plots in figure 12(b) indicate that the optimized coils fix the errors and can produce remarkably good flux surfaces that match the target plasma boundary.

Zoom In Zoom Out Reset image size