An adjoint method for gradient-based optimization of stellarator coil shapes
AAn adjoint method for gradient-based optimization of stellarator coil shapes
E. J. Paul
Department of Physics, University of Maryland, College Park, MD 20742, USA ∗ M. Landreman
Institute for Research in Electronics and Applied Physics,University of Maryland, College Park, MD 20742, USA
A. Bader
Department of Engineering Physics, University of Wisconsin, Madison, WI 53706, USA
W. Dorland
Department of Physics, University of Maryland, College Park, MD 20742, USA
We present a method for stellarator coil design via gradient-based optimization of the coil-windingsurface. The REGCOIL (Landreman 2017
Nucl. Fusion ∗ [email protected] a r X i v : . [ phy s i c s . p l a s m - ph ] J a n I. INTRODUCTION
Stellarators confine particles by generating rotational transform with external coils. The 3-dimensionalnature of a stellarator presents great opportunity, allowing a large space within which to find optimalplasma configurations. However, designing coils to produce the necessary non-axisymmetric magnetic fieldis a significant challenge for the stellarator program. The design of simple coils which can be reasonablyengineered and produce a plasma with optimal physics properties is required in order for the steady-state,disruption-free confinement of optimized stellarators to be realized.Stellarator coils are usually designed to produce a target outer plasma boundary. The plasma boundaryis separately optimized for various physics quantities, including magnetohydrodynamic (MHD) stability,neoclassical confinement, and profiles of rotational transform and pressure [1]. The coil shapes are thenoptimized such that one of the magnetic surfaces approximately matches the desired plasma surface. Ingeneral the desired plasma configuration can not be produced exactly due to engineering constraints on thecoil complexity.In addition to minimization of the magnetic field error, there are several factors that should be consideredin the design of coils shapes. The winding surface upon which the currents lie should be sufficiently separatedfrom the plasma surface to allow for neutron shielding to protect the coils, the vacuum vessel, and a divertorsystem. In a reactor, the coil-plasma distance is closely tied to the tritium breeding ratio and overall costof electricity as it determines the allowable blanket thickness. The coil-plasma distance was targeted in theARIES-CS study to reduce machine size [2]. In practice the minimum feasible coil-plasma separation is afunction of the desired plasma shape. Concave regions (such as the bean W7-X cross section) are especiallydifficult to produce [3] and require the winding surface to be near to the plasma surface. While decreasingthe inter-coil spacing minimizes ripple fields, increasing coil-coil spacing allows adequate space for removalof blanket modules, heat transport plumbing, diagnostics, and support structures. The curvature of a coilshould be below a certain threshold to allow for the finite thickness of the conducting material and to avoidprohibitively high manufacturing costs. The length of each coil should also be considered, as expense willgrow with the amount of conducting material that needs to be produced. For these reasons, identifying coilswith suitable engineering properties can impact the size and cost of a stellarator device.Most coil design codes have assumed the coils to lie on a closed toroidal winding surface enclosing thedesired plasma surface. In
NESCOIL [4], the currents on this surface are determined by minimizing the integral-squared normal magnetic field on the target plasma surface. Using a stream function approach, the currentpotential on the winding surface is decomposed in Fourier harmonics. This takes the form of a least-squaresproblem which can be solved with a single linear system. The coil filament shapes can be obtained from thecontours of the current potential. Because it is guaranteed to find a global minimum,
NESCOIL is often usedin the preliminary stages of the design process [5–7]. It was used for the initial coil configuration studies forNCSX [8]. The W7-X coils were designed using an extension of
NESCOIL which modified the winding surfacegeometry for quality of magnetic surfaces and engineering properties of the coils [9]. However, the inversionof the Biot-Savart integral by
NESCOIL is fundamentally ill-posed, resulting in solutions with amplifiednoise. The
REGCOIL [10] approach addresses this problem with Tikhonov regularization. Here the surface-average-squared current density, corresponding to the squared-inverse distance between coils, is added tothe objective function. With the addition of this regularization term,
REGCOIL is able to simultaneouslyincrease the minimum coil-coil distances and improve reconstruction of the desired plasma surface over
NESCOIL solutions. In this work we build on the
REGCOIL method to optimize the current distributionin 3 dimensions. The current distribution on a single winding surface is computed with
REGCOIL , and thewinding surface geometry is optimized to reproduce the plasma surface with fidelity and improve engineeringproperties of the coil shapes.Other nonlinear coil optimization tools exist which evolve discrete coil shapes rather than continuoussurface current distributions. Drevlak’s
ONSET code [11] optimizes coils within limiting inner and outercoil surfaces. The
COILOPT [12, 13] code, developed for the design of the NCSX coil set [14], optimizescoil filaments on a winding surface which is allowed to vary.
COILOPT++ [15] improved upon
COILOPT , bydefining coils using splines, which allows one to straighten modular coils in order to improve access to theplasma. The need for a winding surface was eliminated with the
FOCUS [16] code, which represents coilsas 3-dimensional space curves. The
FOCUS approach employs analytic differentiation for gradient-basedoptimization, as we do in this work. As the design of optimal coils is central to the development of aneconomical stellarator, it is important to have several approaches. The current potential method could haveseveral possible advantages, including the possible implementation of adjoint methods. Furthermore, thecomplexity of the nonlinear optimization is reduced over other approaches, as the current distribution on thewinding surface is efficiently and robustly computed by solving a linear system. By optimizing the windingsurface it is possible to gain insight into what features of plasma surfaces require coils to be close to theplasma, and what features allow coils to be placed farther away [3].Many engineering design problems can be formulated in terms of the minimization of an objective functionwith respect to some free parameters. A powerful tool for such problems is gradient-based optimization,which requires knowledge of the sensitivity of the objective function with respect to design parameters. Thesegradients can be computed by finite differencing the objective function, but the finite step size introduceserrors and the step size must be chosen carefully. Also, if the optimization space is very large, finite differ-encing can be computationally expensive. Although derivative-free optimization techniques exist, they areless efficient than gradient based algorithms, are limited in the types of constraints that can be implemented,and are typically effective only for small problems [17]. Adjoint methods allow for efficient computationof gradients of the objective function with respect to a large number of design parameters. The cost ofcomputing the derivatives in this way scales independently of the number of design parameters and linearlywith the number of objective functions. In addition to gradient-based optimization, these derivatives canalso be used for uncertainty quantification in scientific computation [18] or to construct sensitivity maps forvisualization of how an objective function changes with respect to normal displacements of a surface [19, 20].Adjoint methods were developed in the 1970s for sensitivity analysis of drag and flow dynamics [21] andhave been widely used for shape optimization in the field of aerodynamics and computational fluid dynamics(CFD) [19, 20, 22–24]. Only recently have these methods been used for tokamak physics in the contextof fitting model parameters with experimental edge data on ASDEX-Upgrade [25] and advanced divertordesign with plasma edge simulations [26]. As stellarator design requires many more geometric parametersthan tokamak design, adjoint-based optimization could provide a significant reduction to computational costto this field.The design of magnetic resonance imaging (MRI) coils has also benefited from adjoint methods [27]. MRIgradient coils which lie on a cylindrical winding surface must provide a specified spatial variation in themagnetic field within a region of interest. This inverse problem is often solved with a linear least-squaressystem by minimizing the squared departure from the desired field at specified points with respect to thecurrent in differential surface elements [28]. This method is comparable to the
NESCOIL [4] approach forstellarator coil design. Gradient coil design was improved by the addition of a regularization term related tothe integral-squared current density [29] or the integral-squared curvature [30], comparable to the
REGCOIL approach. The adjoint method is applied to compute the sensitivity of an objective function with respectto the current potential on the winding surface. Here the Biot-Savart law is written in terms of a matrixequation using the least-squares finite element method, and the adjoint of this matrix is inverted to computethe derivatives [27]. As the adjoint formalism has proven fruitful in this field, we anticipate that it couldhave similar applications in the closely-related field of stellarator coil design.In the sections that follow, we present a new method for design of the coil-winding surface using adjoint-based optimization. An adjoint solve is performed to obtain gradients of several figures of merit, the integral-squared normal magnetic field on the plasma surface and root-mean-squared current density on the windingsurface, with respect to the Fourier components describing the coil surface. A brief overview of the
REGCOIL approach is given in II. The optimization method and objective function are described in section III. Theadjoint method for computing gradients of the objective function is outlined in section IV. Optimizationresults for the W7-X and HSX winding surfaces are presented in section V. In section VI we demonstratea method for visualization of shape derivatives on the winding surface. We discuss properties of optimizedwinding surface configurations in section VII. In section VIII we summarize our results and conclude.
II. OVERVIEW OF THE REGCOIL SYSTEM
First, we review the problem of determining coil shapes once the plasma boundary and coil winding surfacehave been specified. Given the winding surface geometry, our task is to obtain the surface current density, K . The divergence-free surface current density can be related to a scalar current potential Φ, the streamfunction for K , K = n × ∇ Φ . (1)Here n is the unit normal on the winding surface. The current potential Φ can be decomposed into single-valued and secular terms, Φ( θ, ζ ) = Φ sv ( θ, ζ ) + Gζ π + Iθ π . (2)Here ζ is the usual toroidal angle, and θ is a poloidal angle. The quantities G and I are the currents linkingthe surface poloidally and toroidally, respectively. The single-valued term (Φ sv ) is determined by solving the REGCOIL system. It is chosen to minimize the primary objective function, χ = χ B + λχ K . (3)Here χ B is the surface-integrated-squared normal magnetic field on the desired plasma surface, χ B = (cid:90) plasma d A B n . (4)The normal component of the magnetic field on the plasma surface B n includes contributions from currentsin the plasma, current density K on the winding surface, and currents in other external coils. The quantity χ K is the surface-integrated-squared current density on the winding surface, χ K = (cid:90) coil d A K . (5)Here K = | K | . Minimization of χ B by itself ( λ = 0) is fundamentally ill-posed, as very different coil shapescan provide almost identical B n on the plasma surface (for example, oppositely directed currents cancel inthe Biot-Savart integral). The addition of χ K to the objective function is a form of Tikhonov regularization.As we will show, minimization of χ K also simplifies coil shapes. The formulation in REGCOIL allows for finercontrol of regularization while improving engineering properties of the coil set over the
NESCOIL formulation,which relies on Fourier series truncation for regularization.The regularization parameter λ can be chosen to obtain a target maximum current density K max , corre-sponding to a minimum tolerable inter-coil spacing. A 1D nonlinear root finding algorithm is typically usedfor this process.The single-valued part of the current potential Φ sv is represented using a finite Fourier series,Φ sv ( θ, ζ ) = (cid:88) j Φ j sin( m j θ − n j ζ ) . (6)Only a sine series is needed if stellarator symmetry is imposed on the current density ( K ( − θ, − ζ ) = K ( θ, ζ )).As the minimization of χ with respect to Φ j is a linear least-squares problem, it can be solved via the normalequations to obtain a unique solution. The Fourier amplitudes Φ j are determined by the minimization of χ , ∂χ ∂ Φ j = ∂χ B ∂ Φ j + λ ∂χ K ∂ Φ j = 0 , (7)which takes the form of a linear system, (cid:88) j A k,j Φ j = b k . (8)We will use the notation A Φ = b . Throughout bold-faced type will denote the vector space of basis functionsfor Φ sv unless otherwise noted. For additional details see [10]. III. WINDING SURFACE OPTIMIZATION
We use
REGCOIL to compute the distribution of current on a fixed, two-dimensional winding surface.To design coil shapes in 3-dimensional space, we modify the winding surface geometry by minimizing anobjective function (12). This objective function quantifies key physics and engineering properties and iseasy to calculate from the
REGCOIL solution. Optimal coil geometries are obtained by nonlinear, constrainedoptimization.
A. Objective function
The Cartesian components of the winding surface can be decomposed in Fourier harmonics. x = (cid:88) m,n r cmn cos( mθ + nN p ζ ) cos( ζ ) , (9) y = (cid:88) m,n r cmn cos( mθ + nN p ζ ) sin( ζ ) , (10) z = (cid:88) m,n z smn sin( mθ + nN p ζ ) . (11)Here N p is the number of toroidal periods. Stellarator symmetry of the winding surface is assumed( R ( − θ, − ζ ) = R ( θ, ζ ) and z ( − θ, − ζ ) = − z ( θ, ζ ), where R = x + y ). We take the Fourier components ofthe winding surface, Ω = ( r cmn , z smn ), as our optimization parameters and assume a desired plasma surface tobe held fixed. Throughout Ω displayed with a subscript index will refer to a single Fourier component, whilein the absence of a subscript it refers to the set of Fourier components. For a given winding surface geometry,Ω, and desired plasma surface, the current potential Φ(Ω) can be determined by solving the REGCOIL systemto obtain a solution which both reproduces the desired plasma surface with fidelity and maximizes coil-coildistance, as described in section II.We define an objective function, f , which will be minimized with respect to Ω, f (Ω , Φ (Ω)) = χ B (Ω , Φ (Ω)) − α V V / (Ω) + α S S p (Ω) + α K (cid:107) K (cid:107) (Ω , Φ (Ω)) . (12)The coefficients α V , α S , and α K weigh the relative importance of the terms in f . We take χ B (4) as ourproxy for the desired physics properties of the plasma surface. The normal magnetic field depends on Φ , thesingle-valued current potential on the surface, and Ω, the geometric properties of the coil-winding surface.The quantity V coil is the total volume enclosed by the coil-winding surface, V coil = (cid:90) coil d V. (13)We use V / as a proxy for the coil-plasma separation. The quantity S p is a measure of the spectral widthof the Fourier series describing the coil-winding surface [31], S p = (cid:88) m,n m p (cid:16) ( r cmn ) + ( z smn ) (cid:17) . (14)Smaller values of S p correspond to Fourier spectra which decay rapidly with increasing m . We take advantageof the non-uniqueness of the representation in (11) to obtain surface parameterization which are more efficient.There is no unique definition of θ , and minimization of S p removes this redundancy. We use a typical valueof p = 2. The quantity (cid:107) K (cid:107) is the 2-norm of the current density, defined in terms of an area integral overthe surface, (cid:107) K (cid:107) = (cid:32) (cid:82) coil d A | K | A coil (cid:33) / , (15)where A coil is the winding surface area, A coil = (cid:90) coil d A . (16)Although we are using a current potential approach rather than directly optimizing coil shapes, including (cid:107) K (cid:107) in the objective function allows us to obtain coils with good engineering properties. The direct targetingof coil metrics (such as the curvature) introduces additional arbitrary weights in the objective function, andthe solution to another adjoint equation must be obtained to compute its gradient. This will be left forfuture work.FIG. 1: Two non-planar W7-X coils (corresponding to the two leftmost coils in figure 5) computed with REGCOIL using the actual W7-X winding surface. The regularization parameter λ is chosen to achieve theshown values of (cid:107) K (cid:107) . As (cid:107) K (cid:107) increases, the average length, toroidal extent, and curvature increase.To demonstrate this correlation between (cid:107) K (cid:107) and coil shape complexity, we compute the coil set on theactual W7-X winding surface using REGCOIL . The regularization parameter λ is varied to achieve severalvalues of (cid:107) K (cid:107) . Coil shapes are obtained from the contours of Φ. In figure 1, two of the W7-X non-planarcomputed in this way are shown, and the corresponding coil metrics are given in table I. These correspondto the two leftmost coils in figure 5. We consider the average and maximum length l , toroidal extent ∆ ζ ,and curvature κ and the minimum coil-coil distance d mincoil-coil . The average, maximum, and minimum aretaken over the set of 5 unique coils. The coil shapes become more complex as (cid:107) K (cid:107) increases, quantified byincreasing κ and ∆ ζ and decreasing d mincoil-coil . Here the curvature, κ , of a 3-dimensional parameterized curve, r ( t ), is κ = (cid:12)(cid:12)(cid:12)(cid:12) d r dt × d r dt (cid:12)(cid:12)(cid:12)(cid:12)(cid:44)(cid:12)(cid:12)(cid:12)(cid:12) d r dt (cid:12)(cid:12)(cid:12)(cid:12) . (17)We have compared coil shapes on a single winding surface, finding them to become simpler as (cid:107) K (cid:107) decreases.As (cid:107) K (cid:107) = (cid:0) χ K /A coil (cid:1) / , we would find similar trends with χ K . We have chosen to include (cid:107) K (cid:107) in theobjective function as it is normalized by A coil , so it is a more useful quantity for comparison of coil shapeson different winding surfaces.To minimize f , the relative weights in (12) ( α V , α S , and α K ) are chosen such that each of the termsin the objective function have similar magnitudes, though much tuning of these parameters is required toobtain results which simultaneously improve the physics properties (decrease χ B ) and engineering properties(increase V coil and d mincoil-coil , decrease κ and ∆ ζ ). B. Optimization constraints
Minimization of f is performed subject to the inequality constraint d min ≥ d targetmin . Here d min is theminimum distance between the coil-winding surface and the plasma surface, d min = min θ,ζ (cid:0) d coil-plasma (cid:1) = min θ,ζ (cid:18) min θ p ,ζ p (cid:12)(cid:12) r coil − r plasma (cid:12)(cid:12)(cid:19) , (18)and d targetmin is the minimum tolerable coil-plasma separation. The quantities θ p and ζ p are poloidal andtoroidal angles on the plasma surface, r plasma and r coil are the position vectors on the plasma and winding (cid:107) K (cid:107) [MA/m] 2.20 2.70 3.20 K max [MA/m] 4.55 9.50 29.1Average l [m] 8.03 9.18 9.81Max l [m] 8.26 10.5 11.8Average ∆ ζ [rad.] 0.146 0.222 0.253Max ∆ ζ [rad.] 0.161 0.282 0.372Average κ [m − ] 1.04 1.29 1.32Max κ [m − ] 2.54 20.3 56.1 d mincoil-coil [m] 0.353 0.182 0.0758 TABLE I: Comparison of metrics for coils computed with
REGCOIL using the actual W7-X winding surface.Average and max are evaluated for the set of 5 unique coils. The regularization parameter λ is varied toachieve these values of (cid:107) K (cid:107) .surface, and d coil-plasma is the coil-plasma distance as a function of θ and ζ .The maximum current density K max is also constrained, K max = max θ,ζ K. (19)This roughly corresponds to a fixed minimum coil-coil spacing. This constraint is enforced by fixing K max toobtain the regularization parameter λ in the REGCOIL solve, so we avoid the need for an equality constraintor the inclusion of K max in the objective function. Rather, Φ(Ω) is determined such that K max is fixed. Theinequality-constrained nonlinear optimization is performed using the NLOPT [32] software package usinga conservative convex separable quadratic approximation (CCSAQ) [33]. While there are several gradient-based inequality-constrained algorithms available, we chose to use CCSAQ as it is relatively insensitive tothe bound constraints imposed on the optimization parameters. We recognize that there are many possiblecombinations of constraints, objective functions, and regularization conditions that could be used. Forexample, (cid:107) K (cid:107) could be fixed to determine λ while K max could be included in the objective function. Wefound that the formulation we have presented produces the best coil shapes. IV. DERIVATIVES OF f AND THE ADJOINT METHOD
We must compute derivatives of f with respect to the geometric parameters Ω in order to use gradient-based optimization methods. The spectral width S p and volume V coil are explicit functions of Ω, so theiranalytic derivatives can be obtained. On the other hand, χ B and (cid:107) K (cid:107) depend both explicitly on coilgeometry and on Φ (Ω). One approach to obtain the derivatives of these quantities could be to solve the REGCOIL linear system N Ω +1 times, taking a finite difference step in each Fourier coefficient. However, if N Ω (number of Fourier modes) is large, the computational cost of this method could be prohibitively expensive.Instead we will apply the adjoint method to compute derivatives. This technique will be demonstrated below.The derivative of χ B can be computed using the chain rule, ∂χ B (Ω , Φ (Ω)) ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) A Φ = b = ∂χ B ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) Φ + ∂χ B ∂ Φ · ∂ Φ ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) A Φ = b . (20)The subscript A Φ = b indicates that Φ varies with Ω according to (8), with A and b denoting the matrixand right hand side of the linear system in (8). The dot product is a contraction over the current potentialbasis functions, Φ j . We can compute ∂ Φ /∂ Ω j by differentiating the linear system (8) with respect to Ω j , ∂ A ∂ Ω j Φ + A ∂ Φ ∂ Ω j = ∂ b ∂ Ω j , (21)and formally solving this equation to obtain ∂ Φ ∂ Ω j = A − (cid:32) ∂ b ∂ Ω j − ∂ A ∂ Ω j Φ (cid:33) . (22)Equation (22) is inserted into (20), ∂χ B (Ω , Φ (Ω)) ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) A Φ = b = ∂χ B ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) Φ + ∂χ B ∂ Φ · A − (cid:32) ∂ b ∂ Ω j − ∂ A ∂ Ω j Φ (cid:33) . (23)This expression could be evaluated by inverting A for each of the geometric components Ω j and performingthe inner product with ∂χ B /∂ Φ for each Ω j . However, the computational cost of this method scales similarlyto that of finite differencing. Instead, we can exploit the adjoint property of the operator. For a given innerproduct ( , ), the adjoint of an operator, A , is defined as the operator A † satisfying ( b, Ac ) = ( A † b, c ). Aswe are working in R n , the adjoint operator corresponds to the matrix transpose, so ∂χ B (Ω , Φ (Ω)) ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) A Φ = b = ∂χ B ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) Φ + (cid:34)(cid:16) A − (cid:17) T ∂χ B ∂ Φ (cid:35) · (cid:32) ∂ b ∂ Ω j − ∂ A ∂ Ω j Φ (cid:33) . (24)For any invertible matrix, (cid:0) A − (cid:1) T = (cid:0) A T (cid:1) − . Hence we can instead invert the operator A T to compute anadjoint variable q , defined as the solution of A T q = ∂χ B ∂ Φ . (25)Rather than finite difference in each Ω j or invert A for each ∂ Φ /∂ Ω j as in (22), we solve two linear systems:the forward (8) and adjoint (25). The adjoint equation is similar to the forward equation ( A T has the samedimensions and eigenspectrum as A ), so the same computational tools can be used to solve the adjointproblem. We then perform an inner product with q to obtain the derivatives with respect to each Ω j , ∂χ B (Ω , Φ (Ω)) ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) A Φ = b = ∂χ B ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) Φ + q · (cid:32) ∂ b ∂ Ω j − ∂ A ∂ Ω j Φ (cid:33) . (26)The derivatives ∂ b /∂ Ω j , ∂ A /∂ Ω j , (cid:0) ∂χ B /∂ Ω j (cid:1) Φ , and ∂χ B /∂ Φ can be computed analytically. In the abovediscussion, the regularization parameter λ has been assumed to be fixed. A similar method can be used ifa λ search is performed to obtain a target K max (see appendix A). The same method is used to computederivatives of (cid:107) K (cid:107) .We note that adjoint methods provide the most significant reduction in computational cost when thelinear solve is expensive. For the REGCOIL system this is not the case, as the cost of constructing A and b exceeds that of the solve. We have implemented OpenMP multithreading for the construction of ∂ A /∂ Ωand ∂ b /∂ Ω such that the cost of computing the gradients via the adjoint method is cheaper than computingfinite differences serially.The constraint functions, d min and K max , must also be differentiated with respect to Ω j . As d min is definedin terms of the minimum function, we approximate it using the smooth log-sum-exponent function [34]. d min, lse = − q log (cid:82) coil d A (cid:82) plasma d A exp (cid:16) − q (cid:12)(cid:12) r coil − r plasma (cid:12)(cid:12)(cid:17)(cid:82) coil d A (cid:82) plasma d A (27)This function can be analytically differentiated with respect to Ω j . As q approaches infinity, d min, lse ap-proaches d min . For q very large, the function obtains very sharp gradients. A typical value of q = 10 m − was used. The log-sum-exponent function is also used to approximate K max . V. WINDING SURFACE OPTIMIZATION RESULTSA. Trends with optimization parameters
Beginning with the actual W7-X winding surface, we perform scans over the coefficients α V and α S inthe objective function (12). The plasma surface was obtained from a fixed-boundary VMEC solution thatpredated the coil design and is free from modular coil ripple. The constraint target is set to be the minimumcoil-plasma distance on the initial winding surface, d targetmin = 0 .
37 m. The cross sections of the optimizedsurfaces in the poloidal plane are shown in figures 2 and 3 along with the last-closed flux surface (red), aconstant offset surface at d targetmin (black solid), and the initial winding surface (black dashed).With increasing α S at fixed α V = α K = 0, the winding surface approaches a cylindrical torus whichhas a minimal Fourier spectra. At moderately small values of α S (0.3) the surface approaches a constantoffset surface at d targetmin , as χ B is dominant in objective function. For very small values of α S (0.003), wefind that the optimization terminates at a point relatively close to the initial surface, and the resultingwinding surface deviates from a constant offset surface. An intermediate value of α S = 0 . α V is performed at fixed α S = 0 . α K = 0 such that the spectral width does not greatlyincrease. As α V increases, d coil-plasma increases significantly on the outboard side while it remains fixed inthe inboard concave regions. This trend is not surprising, as concave plasma shapes have been shown to beinefficient to produce with coils [3]. Interestingly, the winding surface obtains a somewhat pointed shapeat the triangle cross-section ( ζ = 0 . π/N p ), becoming elongated at the tip of the triangle and ‘pinching’toward the plasma surface at the edges. B. Optimal W7-X winding surface
We now include nonzero α K and attempt a comprehensive optimization. The K max constraint is selectedsuch that the metrics ( l , κ , and ∆ ζ ) of the coils computed on the initial surface roughly match those of theactual non-planar coil set. The coil-plasma distance constraint d targetmin is set to be the minimum d coil-plasma onthe initial winding surface. Parameters α V = 0 . α S = 0 .
24, and α K = 1 . × − were used in the objectivefunction. Optimization was performed over 118 Fourier coefficients (cid:0) | n | ≤ m ≤ (cid:1) and theobjective function was evaluated a total of 5165 times to reach the optimum (1 . × linear solves ratherthan 6 . × required for finite difference derivatives). The optimal surface and coil set are shown in figures4 and 5, and the corresponding metrics are shown in table II. We find a solution which increases V coil by 22%and decreases χ B by 52% over the initial winding surface (note that it is numerically impossible to obtaina current distribution that exactly reproduces the plasma surface, so χ B is nonzero when computed fromthe REGCOIL solution on the initial winding surface). In addition, the optimized coil set features a smalleraverage and maximum ∆ ζ and κ and larger d mincoil-coil . The length of the coils increases to accommodate for theincrease in V coil . Again we find that the increase in V coil is most pronounced in the outboard convex regionswhile d coil-plasma is maintained in the concave regions of the bean-shaped cross-sections. The ‘pinching’feature of the winding surface is again present in the triangle cross-section ( ζ = 0 . π/N p ).It should be noted that the decrease in d coil-plasma at the bottom and top of the bean cross section ( ζ = 0)might interfere with the current W7-X divertor baffles. However, the increase in volume on the outboardside would allow for increased flexibility for the neutral beam injection duct [35]. We have performed thisoptimization to show that a winding surface could be constructed which increases V coil (and thus the average d coil-plasma ), improves coil shapes, and decreases χ B . If further engineering considerations were necessarythese could be implemented. C. Optimal HSX winding surface
We perform the same procedure for optimization of the HSX winding surface. Parameters α V = 3 . × − , α S = 0, and α K = 3 × − were used in the objective function. We found that the spectral widthterm was not necessary to obtain a satisfying optimum in this case. The initial winding surface was taken tobe a toroidal surface on which the actual modular coils lie. The plasma equilibrium used is a fixed-boundary0FIG. 2: Optimized winding surfaces obtained with α V = α K = 0 and the values of α S shown. The actualW7-X winding surface is used as the initial surface in the optimization (black dashed). As α S increases,the magnitude of the spectral-width term in the objective function increases, and the winding surfaceapproaches a cylindrical torus with a minimal Fourier spectra. For moderately small values of α S , thewinding surface approaches a uniform offset surface from the plasma surface (black solid). Initial Optimized Actual coil set χ B [T m ] 0.115 0.0711 V coil [m ] 156 190 (cid:107) K (cid:107) [MA/m] 2.21 2.16 K max [MA/m] 7.70 7.70Average l [m] 8.51 8.95 8.69Max l [m] 8.84 9.14 8.74Average ∆ ζ [rad.] 0.190 0.179 0.198Max ∆ ζ [rad.] 0.222 0.197 0.208Average κ [m − ] 1.21 1.10 1.20Max κ [m − ] 9.01 4.84 2.59 d mincoil-coil [m] 0.223 0.271 0.261 TABLE II: Comparison of metrics of the actual W7-X winding surface and our optimized surface. We alsoshow metrics of the coil set computed on the winding surfaces using
REGCOIL and the metrics for the actualW7-X nonplanar coils. Regularization in
REGCOIL is chosen such that the coil metrics computed on theinitial surface roughly match those of the actual coil set. Coil complexity improves from the initial to thefinal surface (decreased average and max ∆ ζ and κ , increased d mincoil-coil ). The average and max l increases toallow for the increase in V coil .1FIG. 3: Optimized winding surfaces obtained with α S = 0 . α K = 0, and the values of α V shown. Theactual W7-X winding surface is used as the initial surface in the optimization (black dashed). As α V increases, d coil-plasma increases on the outboard side while it remains fixed in the concave region.VMEC solution without coil ripple. Optimization was performed over 100 Fourier coefficients (cid:0) | n | ≤ m ≤ (cid:1) and the objective function was evaluated a total of 560 times to reach the optimum(1 . × linear solves rather than 5 . × required for finite difference derivatives). The coil-plasmadistance constraint was set to be d targetmin = 0 .
14 m, the minimum coil-plasma distance on the actual windingsurface. The optimal surface and coil set are shown in figures 6 and 7, and the corresponding coil metrics areshown in table III. We find a solution which increases V coil by 18% and decreases χ B by 4% over the initialwinding surface. The coil set computed with REGCOIL using the optimized surface appears qualitativelysimilar to that computed with the initial surface but with increased d coil-plasma on the outboard side. Theaverage and maximum ∆ ζ and κ decreased while d mincoil-coil was increased for the coil set computed on theoptimal surface in comparison to that of the initial surface. As was observed in the W7-X optimization (figure4), the optimized HSX winding surface obtains a somewhat pinched shape near the triangle cross-section( ζ = 0 . π/N p ).2 ζ =0.0 π N p ζ =0.25 π N p Actual SurfaceOptimizedPlasma ζ =0.5 π N p ζ =0.75 π N p FIG. 4: The actual W7-X coil-winding surface and plasma surface are shown with our optimized windingsurface. In comparison with the actual surface, the optimized surface reduced χ B by 52% and increased V coil by 22%.3 (a)(b) FIG. 5: Comparisons of coil set computed with
REGCOIL using the actual W7-X winding surface (darkblue) and the optimized surface (light blue).4 ζ =0.0 π N p ζ =0.25 π N p Actual SurfaceOptimizedPlasma ζ =0.5 π N p ζ =0.75 π N p FIG. 6: The actual HSX coil-winding surface and plasma surface are shown with our optimized windingsurface. In comparison with the actual surface, the optimized surface has decreased χ B by 4% andincreased V coil by 18%.5 (a)(b) FIG. 7: The coils obtained from
REGCOIL using the actual HSX winding surface (dark blue) and optimizedsurface (light blue).6
Initial Optimized Actual coil set χ B [T m ] 1 . × − . × − V coil [m ] 2.60 3.07 (cid:107) K (cid:107) [MA/m] 0.956 0.891 K max [MA/m] 1.84 1.84Average l [m] 2.26 2.39 2.24Max l [m] 2.49 2.46 2.33Average ∆ ζ [rad.] 0.372 0.365 0.362Max ∆ ζ [rad.] 0.530 0.505 0.478Average κ [m − ] 5.15 4.80 5.05Max κ [m − ] 33.4 25.8 11.7 d mincoil-coil [m] 0.0850 0.0853 0.0930 TABLE III: Comparison of metrics of the actual HSX winding surface and our optimized surface. We alsoshow metrics of the coil set computed on the winding surfaces using
REGCOIL and the metrics for the actualHSX modular coils. Regularization in
REGCOIL is chosen such that the coil metrics computed on the initialsurface roughly match those of the actual coil set. Coil complexity improves from the initial to the finalsurface (decreased average and max ∆ ζ and κ , increased d mincoil-coil ). The average and max l increases toallow for the increase in V coil .7 VI. WINDING SURFACE SENSITIVITY MAPS
With the adjoint method we have computed derivatives of the objective function with respect to Fouriercomponents of the winding surface, ∂f /∂
Ω. While this representation of derivatives is convenient for gradient-based optimization, visualization of the surface sensitivity in real space is obscured. Alternatively, it ispossible to represent the sensitivity of f with respect to normal displacements of surface area elements of agiven winding surface Ω, δf (Ω , δ r ) = (cid:90) coil d A S δ r · n . (28)Here, S ( θ, ζ ) is a scalar function that will be called the sensitivity. The form (28) implies f is unchangedby tangential displacements of the surface. The shape derivative δf can be formally defined as follows [36].Consider a vector field, δ r , which describes displacements of the surface, Ω. The surface varies smoothlyfrom Ω to Ω (cid:15) , where each point on Ω undergoes transformation T (cid:15) .Ω (cid:15) = (cid:8) T (cid:15) ( r ) : r ∈ Ω (cid:9) , (29)and T (cid:15) is the displacement of each point on the surface by the vector field (cid:15)δ r , T (cid:15) ( r ) = r + (cid:15)δ r ( r ) . (30)The shape derivative, δf (Ω , δ r ), of a functional of the surface geometry, f (Ω), is then defined as δf (Ω , δ r ) = lim (cid:15) → f (Ω (cid:15) ) − f (Ω) (cid:15) . (31)Note that the definition of δf only depends on the direction of δ r , not its magnitude. The shape derivativeis a Gˆateaux derivative, a directional derivative defined for a functional of a vector space. At each point onthe winding surface δf (Ω) is defined for each direction δ r , corresponding to perturbations of the surface atthat location in the specified direction. Under some assumptions, the shape derivative can be representedin the form of (28) (called the Hadamard-Zol`esio structure theorem by some authors) [36]. This so-calledHadamard form for shape derivatives is convenient for computation and has been applied to constructsensitivity maps of Navier-Stokes flows for car aerodynamic design [19, 20]. This representation could havepotential applications for stellarator design, allowing for visualization of regions on the winding surface whichrequire tight engineering tolerances for a given figure of merit.As both χ B and (cid:107) K (cid:107) are defined in terms of surface integrals over the winding surface, it can be shownthat the shape derivative of these functions can be written the Hadamard form [37]. The surface sensitivityfunctions S χ B and S (cid:107) K (cid:107) can be computed from the Fourier derivatives ( ∂χ B /∂ Ω and ∂ (cid:107) K (cid:107) /∂ Ω) usinga singular value decomposition method [38]. Here the perturbations δf and δ r are written in terms of theFourier derivatives, and S is also represented in a finite Fourier series, ∂f∂ Ω j = (cid:90) coil d A (cid:32)(cid:88) mn S mn cos( mθ + nN p ζ ) (cid:33) ∂ r ∂ Ω j · n . (32)After discretizing in θ and ζ , (32) takes the form of a (generally not square) matrix equation which can besolved using the Moore-Penrose pseudoinverse to obtain S mn .We compute S χ B and S (cid:107) K (cid:107) (figure 9) at fixed λ . These quantities are computed on the actual W7-Xwinding surface and a surface uniformly offset from the plasma surface with d coil-plasma = 0 .
61 m (the area-averaged d coil-plasma over the actual surface). We consider surfaces that are equidistant from the plasmasurface on average as S scales inversely with A coil . The poloidal cross-sections of these surfaces are shownin figure 8. For each surface λ is chosen to achieve K max = 7 . S χ B , indicating that d coil-plasma shoulddecrease at that location in order that χ B decreases. This corresponds to locations on the plasma surfacewith significant concavity (see figure 11(b)). The maximum S χ B occurs at ζ = 0 .
15 2 π/N p on both surfaces(see figure 4). In comparison with this region, the magnitude of S χ B is relatively small over the majority ofthe area of the surfaces shown, demonstrating that engineering tolerances might be more relaxed in these8 π N p ζ =0.25 π N p Offset from plasmaActualPlasma π N p ζ =0.75 π N p FIG. 8: The cross sections of the two winding surfaces used to compute S χ B and S (cid:107) K (cid:107) are shown in thepoloidal plane.locations. There is also a region of negative S χ B near ζ = 0 . π/N p and θ = 0. This is the ‘tip’ of thetriangle-shaped cross-section, where d coil-plasma was increased over the course of the optimization (figures 2,3, and 4). We find that S χ B computed on the actual winding surface has similar trends to that computed onthe surface uniformly offset from the plasma. Although on average these surfaces are equidistant from theplasma surface, the magnitude of S χ B is higher on the actual winding surface over much of the area. Thisindicates that the surface sensitivity function depends on the specific geometry of the winding surface. Wehave computed S χ B for several other winding surfaces with varying d coil-plasma . Regardless of the windingsurface chosen, we observe increased sensitivity in the concave regions.The quantity S (cid:107) K (cid:107) roughly quantifies how coil complexity changes with normal displacements of thecoil surface. In view of figure 10, the locations of large S (cid:107) K (cid:107) overlap with areas of increased K . On theactual winding surface, the maximum of S (cid:107) K (cid:107) occurs near the location of closest approach between coils(two rightmost coils in figure 5(a)). The sensitivity functions S (cid:107) K (cid:107) and S χ B have very similar trends. Theconcave regions of the plasma surface are difficult to produce with external coils, resulting in increased coilcomplexity and K . Therefore, (cid:107) K (cid:107) is most sensitive to displacements of the coil-winding surface in theseregions.Studies of the plasma magnetic field sensitivity to perturbations of the coil placement on NCSX similarlyfound that coil errors on the inboard side in regions of small d coil-plasma had a significant effect on fluxsurface quality [39]. The necessity of small d coil-plasma for bean-shaped plasmas has been noted in many coiloptimization efforts [2, 12] and has been demonstrated by evaluating the singular value decomposition ofthe discretized Biot-Savart integral operator [3]. We are able to identify these regions where fidelity of theplasma surface requires tighter tolerance on coil positions using the surface sensitivity function. VII. METRICS FOR CONFIGURATION OPTIMIZATION
Historically, stellarator design has proceeded by first optimizing an equilibrium based on various desiredproperties, such as neoclassical transport and MHD stability. Calculating the coils is a second step, done9 (c) Offset from plasma (d) Actual
FIG. 9: Surface sensitivity functions for χ K (upper subplots) and χ B (lower subplots). These functions arecomputed using the W7-X plasma surface and a uniform offset winding surface from the plasma surfacewith d coil-plasma = 0 .
61 m ((a) and (c)) and the actual winding surface ((b) and (d)). The region ofincreased S χ B corresponds with concave regions of the plasma surface (see figure 11(b)). Regions of largepositive (cid:107) K (cid:107) correspond to regions with increased K (see figure 10)only after the equilibrium has been determined. The results presented here and in [3] indicate that theconcave regions of the surface are both the areas where the optimizing routine chooses winding surfaces thatlie close to the plasma and where the sensitivity to winding surface position is highest.The regions of concavity can be determined by considering the principal curvatures of the plasma surface.Let n represent the normal vector at the plasma surface at some point r , then let A n represent a planethat includes this normal vector. The intersection of the plane and the surface makes a curve r , which hascurvature κ at the point r , as calculated from (17). Then the two principal curvatures P and P representthe maximum and minimum curvatures, κ , from all possible planes A n . The signs of P and P depend onthe convention chosen for the normal vector n . We choose the convention such that convex curves, r , havepositive curvature and concave curves have negative curvatures. Therefore, minima of the second principalcurvature, P , represent regions on the surface where the concavity is maximum.The second principal curvature for the W7-X plasma is shown in figure 11(b). The regions of high concavityare represented by negative values of the second principal curvature. Although P and the sensitivityfunctions are evaluated on different surfaces, we note that regions of high concavity (negative P ) coincidewith regions of high sensitivity (figure 9). The regions of high concavity also correspond to the regions wherethe optimization procedure tends to place the winding surface closest to the plasma (see figure 11). We0 (a) Offset from plasma FIG. 10: Current density magnitude, K , computed from REGCOIL using the W7-X plasma surface and (a) auniform offset winding surface from the plasma surface with d coil-plasma = 0 .
61 m and (b) the actualwinding surface. (a) (b)
FIG. 11: (a) The minimum distance between the W7-X plasma surface and the optimized winding surfaceobtained in section V B and (b) the second principle curvature P are shown as a function of location onthe plasma surface. Locations of large negative P coincide with regions where the optimization resulted insmall d coil-plasma .recognize that our winding surface optimization accounts for several engineering consideration in additionto reproducing the desired plasma surface. However, for a wide range of parameters the winding surfaceswe obtain feature small d coil-plasma in the bean-shaped cross-sections (figures 2 and 3). Thus P , whichis exceedingly fast to compute, may serve as a target for optimization of the plasma configuration. Byminimizing the regions of high concavity, it may be possible to find stellarator equilibria which are moreamenable to coils that are positioned farther from the plasma. Any increase in the minimal distance betweenthe plasma and the coils has implications for the size of a reactor, where the d coil-plasma is set by the blanketwidth.1 VIII. CONCLUSIONS
We have outlined a new method for optimization of the stellarator coil-winding surface using a continuouscurrent potential approach. Rather than evolving filamentary coil shapes, we use
REGCOIL to obtain thecurrent density on a winding surface, and optimize the winding surface using analytic gradients of theobjective function. We have shown that we can indirectly improve the coil curvature and toroidal extent bytargeting the root-mean-squared current density in our objective function (figure 1). This approach offersseveral potential advantages over other nonlinear coil optimization tools.1. The difficulty of the optimization is reduced by the application of the
REGCOIL method, which takesthe form of a linear least-squares system. The optimal coil shapes on a given winding surface can thusbe efficiently and robustly computed.2. By fixing the maximum current density in order to obtain the regularization in
REGCOIL , we eliminatethe need to implement an additional equality constraint or arbitrary weight in the objective function.3. By using
REGCOIL to compute coil shapes on a given surface, we are able to apply the adjoint methodfor computing derivatives (section IV). This allows us to reduce the number of function evaluationsrequired during the nonlinear optimization by a factor of ≈ REGCOIL solutions computed on the initial winding surfaces (tables II and III).Several features of these optimized winding surfaces are noteworthy. While the coil-plasma distance must besmall in concave regions, it can increase greatly on the outboard, convex side of the bean cross-section. Attriangle-shaped cross-sections, the winding surface obtains a somewhat ‘pinched’ appearance (figures 3, 4,and 6). A similar W7-X winding surface shape has been obtained with the
ONSET code (see ref. [11], figure5). Further work is required to understand this behavior.There are several limitations of this approach that should be noted. First, we have applied a local nonlinearoptimization algorithm. This is a reasonable choice if the initial condition is close to a global optimum. Wenote that several global gradient-based optimization algorithms exist, which could be used if a global searchis desired. Second, we currently have not added coil-specific metrics to our objective function (for example,curvature or length). This could be implemented if necessary for engineering purposes.We should also note that this application does not allow for the full benefits of adjoint methods. Whileadjoint methods significantly reduce CPU time if the solve is the computational bottleneck, this is not thecase for the
REGCOIL system, Other applications that are dominated by the linear solve CPU time wouldsee increased benefits from the implementation of an adjoint method. In particular, the field of stellaratordesign could benefit from further incorporation of these methods in other aspects of the design process, suchas computation of neoclassical transport and magnetic equilibria, as stellarators feature complex geometrywith many free parameters describing a given configuration,We demonstrate a technique for visualization of shape derivatives in real space rather than Fourier space.This surface sensitivity function describes how an objective function changes with respect to normal displace-ments of the winding surface. We apply this technique to visualize the derivatives of the integral-squared B n on the plasma surface and the root-mean-squared current density for the W7-X plasma surface and threewinding surfaces (figure 9). This diagnostic identifies the concave regions as being very sensitive to thepositions of coils, as has been observed from previous coil optimization efforts. This visualization techniquecould have potential applications for quantification of engineering tolerances in stellarator design and couldbe extended to the analysis of the sensitivity of other physics properties.2 Appendix A: Adjoint derivative at fixed K max We enforce K max = constant in the REGCOIL solve in order to obtain the regularization parameter λ byrequiring that the following constraint be satisfied within a given tolerance: G (Ω , Φ (Ω , λ )) = K max (Ω , Φ (Ω , λ )) − K targetmax = 0 . (A1)Here K targetmax is the target maximum current density and Φ is chosen to satisfy the forward equation (8), F (Ω , Φ , λ ) = A (Ω , λ ) Φ − b (Ω , λ ) = 0 . (A2)A log-sum-exponent function is used to approximate the maximum function, similar to that used to approx-imate d coil-plasma (27). K max ≈ K max , lse = 1 p log (cid:32) (cid:82) coil d A exp ( pK ) A coil (cid:33) (A3)We compute the total differential of F , d F = (cid:88) j (cid:32) ∂ A ∂ Ω j Φ − ∂ b ∂ Ω j (cid:33) d Ω j + A d Φ + (cid:16) A K Φ − b K (cid:17) dλ = 0 . (A4)Here A K = ∂ A /∂λ and b K = ∂ b /∂λ . We left multiply by A − and solve for d Φ . d Φ = − (cid:88) j A − (cid:32) ∂ A ∂ Ω j Φ − ∂ b ∂ Ω j (cid:33) d Ω j − A − (cid:16) A K Φ − b K (cid:17) dλ (A5)We also compute the total differential of G , dG = (cid:88) j ∂G∂ Ω j d Ω j + ∂G∂ Φ · d Φ = 0 . (A6)Using the form for d Φ (A5), we compute dλ in terms of d Ω j , dλ = (cid:32) ∂G∂ Φ · (cid:20) A − (cid:16) A K Φ − b K (cid:17)(cid:21)(cid:33) − (cid:88) j ∂G∂ Ω j − ∂G∂ Φ · A − (cid:32) ∂ A ∂ Ω j Φ − ∂ b ∂ Ω j (cid:33) d Ω j . (A7)Using (A5) and (A7), the derivative of Φ with respect to Ω j subject to equations (A1) and (A2) is given bythe following expression: ∂ Φ ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) F =0 , G =0 = − A − (cid:32) ∂ A ∂ Ω j Φ − ∂ b ∂ Ω j (cid:33) − A − (cid:0) A K Φ − b K (cid:1) ∂G∂ Φ · (cid:104) A − (cid:0) A K Φ − b K (cid:1)(cid:105) ∂G∂ Ω j − ∂G∂ Φ · A − (cid:32) ∂ A ∂ Ω j Φ − ∂ b ∂ Ω j (cid:33) . (A8)We use the adjoint method to avoid inverting the operator A for each Ω j , ∂ Φ ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) F =0 , G =0 = − A − (cid:32) ∂ A ∂ Ω j Φ − ∂ b ∂ Ω j (cid:33) − A − (cid:0) A K Φ − b K (cid:1) ∂G∂ Φ · (cid:104) A − (cid:0) A K Φ − b K (cid:1)(cid:105) ∂G∂ Ω j − (cid:20)(cid:16) A T (cid:17) − ∂G∂ Φ (cid:21) · (cid:32) ∂ A ∂ Ω j Φ − ∂ b ∂ Ω j (cid:33) . (A9)3We introduce a new adjoint vector (cid:101) q , defined to be the solution of A T (cid:101) q = ∂G∂ Φ . (A10)Equation (A9) is then used to compute the derivatives of χ B with respect to Ω j : ∂χ B ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) F =0 , G =0 = ∂χ B ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) Φ ,λ + ∂χ B ∂ Φ · ∂ Φ ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) F =0 , G =0 . (A11)This result can be written in terms of both adjoint variables, q and (cid:101) q : ∂χ B ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) F =0 , G =0 = ∂χ B ∂ Ω j (cid:12)(cid:12)(cid:12)(cid:12) Φ ,λ − q · (cid:32) ∂ A ∂ Ω j Φ − ∂ b ∂ Ω j (cid:33) − q · (cid:0) A K Φ − b K (cid:1)(cid:101) q · (cid:0) A K Φ − b K (cid:1) ∂G∂ Ω j − (cid:101) q · (cid:32) ∂ A ∂ Ω j Φ − ∂ b ∂ Ω j (cid:33) . (A12)The same method is used to compute derivatives of (cid:107) K (cid:107) . So, to obtain the derivatives at fixed K max , wecompute a solution to the two adjoint equations, (25) and (A10), in addition to the forward equation, (8). ACKNOWLEDGEMENTS
The authors would like to thank I. Abel and T. Antonsen for helpful input and discussions. This workwas supported by the US Department of Energy through grants DE-FG02-93ER-54197 and DE-FC02-08ER-54964. The computations presented in this paper have used resources at the National Energy ResearchScientific Computing Center (NERSC). [1] J. N¨uhrenberg and R. Zille, Physics Letters A , 113 (1988).[2] L. El-Guebaly, P. Wilson, D. Henderson, M. Sawan, G. Sviatoslavsky, R. Slaybaugh, B. Kiedrowski, A. Ibrahim,C. Martin, R. Raffray, S. Malang, J. Lyon, L. P. Ku, X. Wang, L. Bromberg, B. Merrill, L. Waganer, F. Na-jmabadi, and the Aries-CS Team, Fusion Science and Technology , 747 (2008).[3] M. Landreman and A. H. Boozer, Physics of Plasmas , 032506 (2016).[4] P. Merkel, Nuclear Fusion , 867 (1987).[5] D. A. Spong and J. H. Harris, Plasma and Fusion Research , S2039 (2010).[6] L. P. Ku and A. H. Boozer, Nuclear Fusion , 013004 (2011).[7] M. Drevlak, F. Brochard, P. Helander, J. Kisslinger, M. Mikhailov, C. N¨uhrenberg, J. N¨uhrenberg, andY. Turkin, Contributions to Plasma Physics , 459 (2013).[8] N. Pomphrey, L. Berry, A. Boozer, A. Brooks, R. Hatcher, S. Hirshman, L.-P. Ku, W. Miner, H. Mynick,W. Reiersen, D. Strickler, and P. Valanju, Nuclear Fusion , 339 (2001).[9] C. Beidler, G. Grieger, F. Herrnegger, E. Harmeyer, J. Kisslinger, , W. Lotz, H. Maassberg, P. Merkel,J. N¨uhrenberg, F. Rau, J. Sapper, F. Sardei, R. Scardovelli, A. Schl¨uter, and H. Wobig, Fusion Technology ,148 (1990).[10] M. Landreman, Nuclear Fusion , 046003 (2017).[11] M. Drevlak, Fusion Technol. , 106 (1998).[12] D. J. Strickler, L. A. Berry, and S. P. Hirshman, Fusion Science and Technology , 107 (2002).[13] D. J. Strickler, S. P. Hirshman, D. A. Spong, M. J. Cole, J. F. Lyon, B. E. Nelson, D. E. Williamson, and A. S.Ware, Fusion Science and Technology , 15 (2004).[14] M. Zarnstorff, L. Berry, A. Brooks, E. Fredrickson, G.-Y. Fu, S. Hirshman, S. Hudson, L.-P. Ku, E. Lazarus,D. Mikkelsen, D. Monticello, G. H. Neilson, N. Pomphrey, A. Reiman, D. Spong, D. Strickler, A. Boozer, W. A.Cooper, R. J. Goldston, R. Hatcher, M. Y. Isaev, C. Kessel, J. Lewandowski, J. Lyon, P. Merkel, H. Mynick,B. Nelson, C. N¨uhrenberg, M. Redi, W. Reiersen, P. Rutherford, R. Sanchez, J. Schmidt, and R. White, PlasmaPhysics and Controlled Fusion , A237 (2001).[15] T. Brown, J. Breslau, D. Gates, N. Pomphrey, and A. Zolfaghari, in IEEE 26th Symposium on Fusion Engineering(SOFE) (Austin, Texas, 2015).[16] C. Zhu, S. R. Hudson, Y. Song, and Y. Wan, Nuclear Fusion , 016008 (2018). [17] J. Nocedal and S. J. Wright, Numerical Optimization (Springer, 2006) pp. 193–194.[18] C. J. Roy and W. L. Oberkampf, Computer Methods in Applied Mechanics and Engineering , 2131 (2011).[19] C. Othmer, International Journal for Numerical Methods in Fluids , 861 (2008).[20] C. Othmer, Journal of Mathematics in Industry , 6 (2014).[21] O. Pironneau, Journal of Fluid Mechanics , 97 (1974).[22] G. Kuruvila, S. Ta’asan, and M. Salas, in (Reno, Nevada, 1995).[23] A. Jameson, L. Martinelli, and N. Pierce, Theoretical and Computational Fluid Dynamics , 213 (1998).[24] W. K. Anderson and V. Venkatakrishnan, Computers and Fluids , 443 (1999).[25] J. Kim, D. P. Coster, J. Neuhauser, R. Schneider, and the ASDEX-Upgrade Team, Journal of Nuclear Materials , 644 (2001).[26] M. Baelmans, M. Blommaert, W. Dekeyser, and T. Van Oevelen, Nuclear Fusion , 036022 (2017).[27] F. Jia, Z. Liu, M. Zaitsev, J. Hennig, and J. G. Korvink, Structural and Multidisciplinary Optimization , 523(2014).[28] R. Turner, Magnetic Resonance Imaging , 903 (1993).[29] L. K. Forbes, M. A. Brideson, and S. Crozier, IEEE Transactions on Magnetics , 2134 (2005).[30] L. K. Forbes and S. Crozier, Journal of Physics D: Applied Physics , 3447 (2001).[31] S. P. Hirshman and H. K. Meier, Physics of Fluids , 1387 (1985).[32] S. G. Johnson, “The NLopt nonlinear-optimization package,” (2014).[33] K. Svanberg, SIAM Journal of Optimization , 555 (2002).[34] S. Boyd and L. Vandenberghe, Convex Optimization (Cambridge University Press, 2004) p. 72.[35] N. Rust, B. Heinemann, B. Mendelevitch, A. Peacock, and M. Smirnow, Fusion Engineering and Design , 728(2011).[36] M. C. Delfour and J.-P. Zol´esio, Shapes and Geometries (Society for Industrial and Applied Mathematics, 2011)pp. 480–481.[37] A. A. Novotny and J. Sokolowski,
Topological Derivatives in Shape Optimization (Springer, 2013) pp. 35–38.[38] M. Landreman, (in preparation).[39] D. Williamson, A. Brooks, T. Brown, J. Chrzanowski, M. Cole, H.-M. Fan, K. Freudenberg, P. Fogarty, T. Har-grove, P. Heitzenroeder, G. Lovett, P. Miller, R. Myatt, B. Nelson, W. Reiersen, and D. Strickler, FusionEngineering and Design75-79