Calculation of permanent magnet arrangements for stellarators: A linear least-squares method
CCalculation of permanent magnet arrangements forstellarators: A linear least-squares method
Matt Landreman and Caoxiang Zhu
1. Institute for Research in Electronics and Applied Physics, University of Maryland,College Park, MD 20742, USA2. Princeton Plasma Physics Laboratory, Princeton, NJ 08543, USAE-mail: [email protected]
Abstract.
A problem arising in several engineering areas is to design magnets outsidea volume that produce a desired magnetic field inside it. One instance of this problemis stellarator design, where it has recently been shown that permanent magnets canprovide the required shaping of the magnetic field. Here we demonstrate a robustand efficient algorithm REGCOIL PM to calculate the spatial distribution of thesepermanent magnets. The procedure involves a small number of fixed-point iterations,with a linear least-squares problem solved at each step. The method exploits theBiot-Savart Law’s exact linearity in magnetization density and approximate linearityin magnet size, for magnets far from the target region. No constraint is placed onthe direction of magnetization, so Halbach solutions are found naturally, and themagnitude of the magnetization can be made uniformly equal to a target value.
1. Introduction
Given a desired magnetic field in some region, what arrangement of magnets outside theregion can produce this field? This problem has many applications, including magneticresonance imaging [1, 2, 3], particle accelerators [4, 5], and stellarators for magneticconfinement of plasma [6, 7, 8, 9, 10]. In stellarators, a three-dimensional magneticfield must be carefully shaped in order to provide good confinement of charged particletrajectories and meet other physics objectives. For all these applications, to designmagnets one essentially needs to invert the Biot-Savart law. The Biot-Savart lawprovides a straightforward way to compute the magnetic field B from known currents.However the inverse problem (given B , determine currents) is ill-posed [9] in the sensethat very different currents can give nearly the same B . Therefore there is room forcreativity and innovation in formulating algorithms for these inverse problems such thatthe magnet designs obtained are practical, and the computational cost is low.While a variety of algorithms have been devised to compute the shapes ofelectromagnetic coils [6, 7, 8, 9, 10], in the area of stellarators there is much lessexperience with algorithms for designing permanent magnets. It has recently been a r X i v : . [ phy s i c s . p l a s m - ph ] S e p alculation of permanent magnet arrangements for stellarators M . Allowing the direction of M to be arbitraryis expected to reduce the necessary volume of permanent magnets by approximately afactor of 2 [11]. We do however aim to constraint the magnitude M = | M | , making iteverywhere uniform in the magnet region, since a technical upper limit exists, roughly µ M ≤ . B PM ( r ) = (cid:90) V d r (cid:48) µ π | r − r (cid:48) | (cid:20) r − r (cid:48) )( r − r (cid:48) ) · M ( r (cid:48) ) | r − r (cid:48) | − M ( r (cid:48) ) (cid:21) . (1)Here, the integral is performed over a region V containing the permanent magnets, and B PM is the field at position r due to a magnetization (dipole moment density) M atposition r (cid:48) . Then if the magnet region V is considered fixed, (1) indicates that B ( r )is a linear function of M ( r (cid:48) ). Regularization is required to make the problem well-posed, but Tikhonov-like regularization can be introduced that preserves the linearity.Solving this linear problem for M is naturally faster and more robust than solving anonlinear problem, and it naturally allows the direction of M to be arbitrary. (While M differs slightly from the zero-field magnetization M , the difference is very small forrare-Earth magnets and so will be neglected here. If desired, the material’s relationship M ( B ) could be inverted at each point to obtain M .)However the solution of the linear problem at fixed V will have nonuniform M = | M | , whereas it is preferable to have a solution with M uniformly equal to the limit alculation of permanent magnet arrangements for stellarators d A M B (a) d = xd A M = M / x (b) Permanent magnetEvaluation point B ≈ B Figure 1.
Physical picture of the fixed-point iteration in the REGCOIL PMalgorithm used to achieve a uniform magnetization M = | M | . (a) Suppose a region ofmagnetization M produces a given field B at a distant evaluation region. (b) If onedimension of the magnet is scaled by some factor x and the magnetization is scaled by1 /x , the change to B in the evaluation region is small. of the magnet material. Uniformity can be achieved by noting that the volume integralin (1) makes B PM approximately linear in the magnet thickness. The complicateddependence on r − r (cid:48) in (1) spoils this linearity, but if the magnets are not extremelyclose to the evaluation region, the nonlinearity will be weak. Thus, we can approximatelycorrect the nonuniformity in M by adjusting the magnet thickness. If M is too largeby a factor of two in a given region of V , doubling the thickness of V in this region(at fixed B ) will result in a lowering of M by approximately the same factor of two.This relationship is illustrated in figure 1. While M will not be exactly uniform afterthis change to the shape of V , the procedure can be iterated to improve the uniformity.As we will demonstrate, the number of iterations required can be quite small. While aNewton-type iteration could also be used to account for the nonlinearity, such a methodwould require either analytic derivatives or finite-difference derivatives, increasing thecomputational cost. We will show that the simpler Picard fixed-point iteration is stableand sufficient in practice.A serious experimental design for a permanent magnet stellarator requires a detailedgeometry model with many individual magnet pieces, as in [13]. Here instead wewill make a crude approximation that the magnetization fills a single region withsmooth curved boundaries. This smooth model roughly approximates a large numberof individual magnets, but we acknowledge this approximation is likely insufficient fora serious experimental design. We make this smooth approximation here for severalreasons. The primary reason is that it allows reuse of a significant amount of theREGCOIL code [9]. This made it possible to try out the idea rapidly. Second, thesmooth model is numerically convenient, since it allows us to evaluate integrals withspectral accuracy using uniform grids in periodic coordinates. Having demonstrated theREGCOIL PM algorithm in this paper, it should be straightforward to apply it in thefuture to more realistic geometry models with discrete magnets.In the remainder of this paper, this REGCOIL PM algorithm is defined in greaterdetail and demonstrated for several problems. Section 2 gives further detail about themathematical formulation. Aspects of the discretization and numerical solution arediscussed in section 3, using the NCSX stellarator as an example. In section 4 it is alculation of permanent magnet arrangements for stellarators
2. Mathematical formulation
We consider the common 2-stage approach to stellarator design. In the first stage,the parameter space for optimization is the space of toroidal plasma boundary shapes,and the objective function is a combination of physics figures of merit for the plasmainside this boundary. In the second stage, the shapes of magnets are optimized toproduce the plasma boundary shape resulting from the first stage. Our goal in thispaper is to solve the stage-2 problem. If the stage-2 problem can be approximately solved quickly, the solution can be incorporated into the stage-1 objective function topenalize magnet complexity [15]. In this way, the stage-1 optimization can be made tofind plasma configurations that can be supported by magnets of low complexity, and amore detailed and computationally demanding stage-2 calculation can be done for thefinal magnet design.We thus focus on the problem of finding a permanent magnet arrangement toproduce a desired plasma boundary surface S . To state this problem precisely, firstconsider that S must be a magnetic surface, so we wish to make B · n ≈ S . (Matching the normal component is sufficient to ensure that the full vector B coincides with the target field everywhere inside S .) We then use the linearity ofmagnetostatics to write B = B PM + B f where B PM is the magnetic field (1) producedby the permanent magnets and B f is the field produced by currents that are “fixed”during the permanent magnet design. The quantity B f represents contributions fromthe electromagnets and from current in the plasma, if there is any. Our goal can thenbe stated as achieving f B ≈ f B = (cid:90) S [( B PM + B f ) · n ] d x, (2)an integral over the surface S of the squared normal component of the field.The problem of finding M such that f B ≈ B on and inside S . Here, we will make the problem well-posed in two steps, first constraining themagnet location and then adding a regularization term. In the first step, we restrict thepermanent magnets to lie within a volume V ( d ) with some thickness parameter d . Inthis paper we will take V to be bounded by two toroidal surfaces, a fixed inner surface I and a variable outer surface O , both linking the plasma surface S , with d a function on I that measures the distance to O (figure 2). Specifically, we choose V to be the rangeof position vectors r (cid:48) ( s, θ, ζ ) = r I ( θ, ζ ) + σ s d ( θ, ζ ) n ( θ, ζ ) (3) alculation of permanent magnet arrangements for stellarators Figure 2.
Definitions of geometric quantities. with θ ∈ [0 , π ), ζ ∈ [0 , π ), and s ∈ [0 , r I ( θ, ζ ) is the position vector on thefixed inner surface I of the magnet region, θ and ζ are any poloidal and toroidal angles, d ( θ, ζ ) is a thickness function that will be varied, n = N / | N | is a unit normal of theinner surface, N = ∂ r I ∂ζ × ∂ r I ∂θ (4)is a non-unit-length normal vector, and σ = ± O of V at s = 1 is outside the inner surface I . As long as d does not exceeda ( θ, ζ )-dependent threshold related to the curvature of I , the map r (cid:48) ( s, θ, ζ ) in (3)is invertible. This continuous model for V is an approximation to a set of many smalldiscrete magnets. Other choices for V are possible, such as a set of n discrete hexahedralvolumes, with d a vector of n numbers giving the thickness of each region. However V ( d ) is defined, we will take d to be unknown, to be determined during the permanentmagnet design. Restricting the permanent magnet location to a parameterized volume V ( d ) is a reasonable reflection of practical engineering considerations: we expect thepermanent magnets should be as close to the plasma as possible, limited by the vacuumvessel and any other components that may need to lie in between, but the thickness ofthe permanent magnet region depends on the specific plasma configuration.However, restricting the permanent magnets to be located in a specific V ( d ) like(3) does not fully eliminate the ill-posedness in the problem f B ≈
0. Consider thatnear any point a finite distance from S we could add two oppositely directed magnetswith substantial magnetization, and if these magnets were sufficiently small and closeto each other, the change to B on S would be negligible. This type of ill-posednessis further discussed in [9]. To arrive at a well-posed problem we therefore introduce aregularization term. The most convenient term to introduce is f M = (cid:90) V | M | w ( θ, ζ ) d ( θ, ζ ) d x, (5)a weighted volume integral over the permanent magnet region of the squaredmagnetization density. Here w ( θ, ζ ) is an optional user-supplied weight function that alculation of permanent magnet arrangements for stellarators d ( θ, ζ ) in (5) is motivated by the fixed-point iteration that will be explained shortly. The form of (5) is essentially Tikhonovregularization, but with a physically meaningful weighting. (A similar regularizationterm without the weighting factors was proposed in appendix B of [14].) We can nowdefine a combined objective function f = f B + λf M , (6)where λ is a positive scalar parameter controlling the amount of regularization. Smallvalues of λ correspond to precisely making the target magnetic field (very small f B ) atthe expense of more complicated permanent magnet structures, while large values of λ yield simplified magnet structures at the expense of magnetic field inaccuracies (larger f B ). The problem of finding M that minimizes f at fixed V has the form of a linearleast-squares problem.We can now define the fixed-point iteration to make M uniform. We constrain M to be independent of the radial coordinate s , so M = M ( θ, ζ ). The physicalpicture in figure 1 can then be expressed as d ( θ, ζ ) M ( θ, ζ ) ≈ d ( θ, ζ ) M ( θ, ζ ), wherethe subscripts 1 and 2 refer to a pair of magnet configurations such as panels (a)-(b) offigure 1. If we desire for M to equal a uniform target value M t , then the appropriateupdate rule for d is d j +1 = d j M j M t . (7)The factor of d in (5) can now be explained. The iteration (7) preserves the product M d . In (5), the volume integral contains an explicit M d factor as well as an implicit d factor through the thickness of the integration region V . Therefore (5) is approximatelyconstant during the iterations (7). The explicit d factor in (5) is not critical, but it isconvenient since a “good” value of λ for the first iteration is likely to be a good valuealso for the final iteration.The REGCOIL PM method can now be summarized. First, d ( θ, ζ ) is initializedto a uniform thickness d , and choices of λ and M t are fixed. Then the least-squaresproblem of minimizing f (at fixed d ) is solved for M . An updated thickness d ( θ, ζ ) iscomputed from (7). Using the new V derived from d and (3), the least-squares problemis solved again to yield M . An updated thickness d ( θ, ζ ) is computed from (7), andthe process is repeated until successive iterates are sufficiently close to each other.To evaluate the volume integral in (5), we note (cid:82) V = (cid:82) π dθ (cid:82) π dζ (cid:82) ds |√ g | wherethe Jacobian derived from (3) is √ g = − σdN (cid:2) − σsdH + ( LQ − P ) s d /N (cid:3) , (8)where N = | N | , H = ( LG + QE − P F ) / (2 N ) is the mean curvature of the inner alculation of permanent magnet arrangements for stellarators I , and E = ∂ r I ∂θ · ∂ r I ∂θ , F = ∂ r I ∂θ · ∂ r I ∂ζ , G = ∂ r I ∂ζ · ∂ r I ∂ζ , (9) L = n · ∂ r I ∂θ , P = n · ∂ r I ∂θ∂ζ , Q = n · ∂ r I ∂ζ . This result is derived in Appendix A.
3. Numerical solution
We now discuss the discretization and numerical solution of the equations of the previoussection. The source code for the numerical implementation used here is available onlineat [16], and data for the figures and benchmarks is available at [17].
The magnetization vector is first written as a finite sum of basis functions M ( θ, ζ ) = j max (cid:88) j =1 3 (cid:88) k =1 M j,k p j ( θ, ζ ) e k ( ζ ) , (10)where p j ( θ, ζ ) = (cid:32) sincos (cid:33) j ( m j θ − n j n fp ζ ) (11)are angular basis functions, and e k for k = 1 , , e R , e φ , e Z . The notation in (11) means that either sin or cos is chosen forbasis function j . The number of identical field periods is denoted n fp . The integers m j range from 0 to m max , and n j ranges from 0 to n max (for m j = 0) or − n max to n max (for m j > ζ to be equal to the standard toroidal angle φ on the inner surface, so(3) implies that ζ is not generally the standard toroidal angle off of I . The vectors e k ( ζ )are evaluated at the point on I with the given ζ , meaning that they generally differ fromthe cylindrical basis vectors at any point off I where M is evaluated. The reason forthis choice is so the Cartesian components of M remain constant as you move in thenormal direction from I , reflecting a reasonable engineering constraint.In the common case of stellarator symmetry, only the sin( m j θ − n j ζ ) basis functionsneed to be included in (10) for the e R terms, and only the cos( m j θ − n j ζ ) basis functionsneed to be included for the e φ and e Z terms.If one wished to allow M to vary with s , a sum over basis functions in s (suchas polynomials) could be included in (10). Our numerical implementation allows thispossibility. However if this s dependence is allowed, it is hard to see how to achieve auniform M by the fixed-point iteration proposed here. Therefore for all results in thispaper we do not include s dependence in (10). alculation of permanent magnet arrangements for stellarators f B involves an integral over the plasma boundary surface S .This integral is written as a discrete sum using a uniformly spaced grid of N θ points in θ and a uniform grid of N ζ points in ζ , with θ and ζ angles on S in this case. Moreover,to evaluate f B and f M , integrals over the magnetization volume are required. Theseintegrals are written as discrete sums using a uniform grid of N θ (cid:48) points in θ , a uniformgrid of N ζ (cid:48) points in ζ , and a Gauss-Legendre grid of N s points in s . The magnetthickness d ( θ, ζ ) is stored on the same N θ (cid:48) × N ζ (cid:48) discrete grid points.It could be reasonable to take the independent variables as the components of M on the discrete N θ (cid:48) × N ζ (cid:48) grid points, instead of using the Fourier amplitudes M j,k in(10). We choose the Fourier approach here due to two advantages. First, it is convenientfor imposing stellarator symmetry, which reduces the number of degrees of freedom bya factor of two. Second, it allows the volume integrals in f B and f M to be evaluatedat higher spatial resolution without increasing the number of degrees of freedom forthe least-squares solution, which is numerically efficient in practice. In the Fourierrepresentation used here, one should choose N θ (cid:48) ≥ m max + 1 and N ζ (cid:48) /n fp ≥ n max + 1so there are at least as many degrees of freedom in the grid as in Fourier space. Otherwise f M does not fully regularize every Fourier mode.With M represented by the finite sum (10), and the integrals over S and themagnetization volume in f B and f M approximated by finite sums as described above,minimization of f now has the form of a finite linear-least-squares problem. Suchproblems can be solved by standard methods such as the normal equations, QR decomposition, or singular value decomposition. Before demonstrating the entire REGCOIL PM algorithm, it is valuable to first examinethe behavior of the least-squares solution at fixed d , without the fixed-point iteration.For this discussion we will use the geometry shown in figure 3. This figure showsslices through the geometry at constant φ , where φ is the standard toroidal angle,coinciding with ζ only on the inner magnet surface I . The plasma geometry is that ofthe c09r00 version of NCSX, a free-boundary equilibrium computed using infinitesmallythin approximations of the 18 discrete modular coils. For this paper we will use thec09r00 boundary shape but neglect plasma currents. The contribution to B f fromplasma current could be computed using the same virtual-casing method [18] used forother stellarator coil calculations. The fixed field B f is taken to be a purely toroidalfield, approximating the field from a large number of toroidal field coils. The mean fieldin the plasma region is 0.5 T, one third of the original NCSX design, as this is whatcan be supplied with the array of planar toroidal field coils built for NCSX. The innermagnet surface is taken to be the NCSX vacuum vessel. This vessel is not a uniformdistance from the plasma boundary. For this subsection we consider a uniform magnetthickness d = 0.1 m. The weight w in (5) is set to 1 until section 5.1.Figure 4 shows the trade-off curve (“Pareto frontier”) between f B and f M as the alculation of permanent magnet arrangements for stellarators Z [ m e t e r s ] = 0 Z [ m e t e r s ] = (1/4)(2 /3) Z [ m e t e r s ] = (1/2)(2 /3) Outer magnetization surfaceInner magnetization surfacePlasma boundaryMagnetization regionPlasma region (a) (b) (c)
Figure 3.
Geometry for the discussion of regularization and resolution parametersin section 3.2. Here, φ is the standard toroidal angle, which corresponds to ζ on theinner surface I but not off of I . level of regularization λ is varied. Ideally both f B and f M would both be small, but atrade-off must be made: a small value of one of these quantities requires a large valuefor the other. The trade-off curve plotted is actually 5 curves overlaid, showing thatfactor-of-2 changes in each numerical resolution parameter has negligible effect on thesolutions (Table 1). Three red points indicate solutions that are shown in detail in figure5 At large λ , the trade-off curve extends infinitely far to the left. In this limit, f ≈ λf M , so the solution is M →
0. With no permanent magnets, f B has a nonzerovalue associated with the fixed field B f . At the other limit of small λ , arbitrarily smallvalues of f B and arbitrarily large values of f M are obtained. (The curve eventuallybends to the right but very large numerical resolution is required in this region, so onlythe converged section is displayed.) In this limit, the component of B normal to thetarget plasma surface is made arbitrarily small due to extremely large values of M . Theregularization vanishes in this limit, so very short-scale patterns in M arise. A usermust choose an intermediate value of λ that balances magnet complexity against physicsproperties of the plasma configuration. We now add the fixed-point iteration (7), considering the same 0.5 T NCSX geometryfrom the previous subsection. The iteration converges fastest when the number ofdegrees of freedom in d is close to the number of degrees of freedom in each componentof M , i.e. when N θ (cid:48) = 2 m max + 1 and N ζ (cid:48) /n fp = 2 n max + 1. Otherwise the spatialdependence of d and M in (7) does not match. Therefore for this section we use theparameters of run 5 from table 1. We choose λ = 10 − T / A . We also choose a targetmagnetization M t = 1 . /µ ≈ .
114 MA / m, achievable with rare-Earth magnets.Figure 6 shows the convergence of the fixed-point iterations. It can be seen that the alculation of permanent magnet arrangements for stellarators f M [Amperes meters ]10 f B [ T e s l a m e t e r s ] = 10 = 10 = 10 High regularization Ideal solutionswould be here Low regularization
Figure 4.
Trade-off curve as the regularization parameter λ is varied, at fixed magnetthickness d =0.1 m. Units of λ are Tesla / Ampere . Red dots show the solutions infigure 5. Table 1.
Resolution parameters for the five overlaid blue curves in figure 4. N θ and N ζ : Number of grid points in the poloidal and toroidal angles on the plasmasurface. N θ (cid:48) and N ζ (cid:48) : Number of grid points in the poloidal and toroidal angles inthe magnetization region. m max and n max : Maximum Fourier mode numbers for thecylindrical components of M . N s : Number of Gauss-Legendre points for integrationover the radial coordinate s in the magnetization region. n fp : Number of identicalfield periods.Run N θ = N ζ /n fp N θ (cid:48) = N ζ (cid:48) /n fp m max = n max N s minimum and maximum of M over V both quickly converge to the target M t . Figure6.b shows the difference in d between successive iterates, measured by the maximumover θ and ζ of | d j − − d j | . The difference converges to zero, demonstrating that a fixedpoint has been found.Also shown in figure 6 are results when Anderson acceleration [19, 20] is appliedto the iteration. In Anderson acceleration, a linear combination of the previous fewiterates is used instead of only the previous iterate. The extra computational cost ofthe Anderson step compared to (7) is so small as to be negligible, and here it provides alculation of permanent magnet arrangements for stellarators M [A/m] | B n | [T] M [A/m] | B n | [T] M [A/m] | B n | [T] = 10 T /A = 10 T /A = 10 T /A Figure 5.
Trends as the regularization parameter λ is varied in the least-squaresproblem min f for fixed magnet thickness d . The three solutions shown correspondto the red points in figure 4. As regularization is reduced, the residual normal field B n = ( B PM + B f ) · n on the target plasma surface is reduced, but at the expense ofgreater magnetization magnitude and finer structure in the magnetization. a modest acceleration in convergence.The evolution of the spatial dependence of d and M is shown in figure 7. It canbe seen that both d and M converge rapidly. By eye, M is uniform and equal to M t byiteration 3, and changes to d are hardly visible after iteration 1. The final result for theshape of the magnet region is displayed in figure 8.It is not obvious that for any choice of initial d , the iteration is stable and the fixedpoint obtained is the same. However it appears that a unique solution exists in practice,at least for the examples in this paper. The calculation of this section was repeated withvarious uniform initial d ∈ { . , . , . } meters. (Larger values are not permittedbecause the outer surface begins to self-intersect as √ g crosses through zero.) Thecalculation was also repeated taking the initial d to have a random variation in θ and ζ within [0 , .
1] m. As shown in figure 9, differences between these differently-initializedcalculations converged steadily towards zero as the iterations proceeded. This behavioris in contrast to the formulations in [14, 13] in which dependence on the initial conditionwas observed. The independence of REGCOIL PM results from the initial condition isadvantageous, since a user need not worry about how best to select the initial condition. alculation of permanent magnet arrangements for stellarators Relative departure from target magnetization magnitude M t Picard, max( M )/ M t M )/ M t Anderson2, max( M )/ M t M )/ M t Anderson3, max( M )/ M t M )/ M t m a x , | d j + d j | [ m e t e r s ] Difference in magnet thickness between successive iterates
PicardAnderson2Anderson3 (a) (b)
Figure 6.
Convergence of the fixed-point iterations, showing both the Picard update(7) and its Anderson-accelerated variant. The 2 or 3 after Anderson refers to the‘depth’ of Anderson acceleration. (a) Both the maximum and minimum of M overthe magnetization volume converge to the desired value M t , i.e. M becomes uniform.(b) The difference in the shape of the magnetization region between successive iteratesconverges to zero. d [m] M [MA/m] d [m] M [MA/m] d [m] M [MA/m] d [m] M [MA/m] d [m] M [MA/m] Iteration 0 Iteration 1 Iteration 2 Iteration 3 Iteration 4
Figure 7.
Convergence of the magnet thickness d and magnetization magnitude M during the fixed-point iterations (7). By iteration 3, deviations of M from the targetvalue 1.1 MA/m are invisible on the scale of the figure. alculation of permanent magnet arrangements for stellarators Z [ m e t e r s ] = 0 Z [ m e t e r s ] = (1/4)(2 /3) Z [ m e t e r s ] = (1/2)(2 /3) Outer magnetization surfaceInner magnetization surfacePlasma boundaryMagnetization regionPlasma region (a) (b) (c)
Figure 8.
Cross-sections of the REGCOIL PM solution for the NCSX example insection 3.3. m a x , | d j d . j | [ m e t e r s ] Difference in d compared to d initial = 0.1 m Initial uniform d = 0.15 mInitial uniform d = 0.01 mInitial uniform d = 0.001 mRandom nonuniform initial d Random nonuniform initial d Random nonuniform initial d Figure 9.
For different choices of the initial magnet thickness d , the iterations (7)converge to the same solution. Here, d . j denotes d at iteration j for a calculationinitialized with a uniform d = 0 . | d j − d . j | iscomputed, and the maximum of this difference over θ and ζ is plotted. The differenceconverges towards zero.
4. Verification for Halbach cylinders
A satisfying property of the mathematical formulation of section 2 is that it is consistentwith the analytic solution for cylindrical multipole magnets (“Halbach cylinders”)described by Halbach [21]. A comparison with this analytic result also serves as auseful test of the numerical implementation of section 3. alculation of permanent magnet arrangements for stellarators plasmasurface magnetsurface plasmasurface magnetsurface (a) (b) Figure 10.
The two configurations in cylindrical geometry analyzed in section4.1. (a) Dipoles with uniform magnitude but varying direction, yielding the Halbachcylinder solution. (b) Dipoles oriented normal to the magnet surface with varyingmagnitude.
We first derive the analytic solution by a different method than in [21] to highlightthe parallels with stellarator magnet optimization. We consider two concentric infinitecylinders, an inner one with radius a analogous to the plasma surface, and an outerone with radius b > a analogous to a thin magnet volume. This configuration can beimagined as a high-aspect-ratio limit of an axisymmetric system, so the angle aroundthe cylinder θ is a poloidal angle. Let us try to arrange magnetic dipoles on the outersurface in order to create a normal magnetic field B n = B · n = ¯ B cos( (cid:96)θ ) (12)on the plasma surface, where (cid:96) is a given integer. In other words, suppose there is a fixednormal field B n,f = − ¯ B cos( (cid:96)θ ), and we wish to introduce dipoles to obtain f B = 0. Wewill consider two possible arrangements of dipoles, shown in figure 10: first, dipoles ofuniform magnitude but arbitrary direction, as in a REGCOIL PM solution; and second,dipoles oriented normal to the outer surface but with arbitrary magnitude. This secondcase is considered because dipoles oriented normal to a surface have been considered inrecent papers [12, 13]. We will show the magnetization magnitude in the first approachis half of the maximum magnetization in the second.Outside of the region of dipoles, the magnetic field can be written B = ∇ Φ for ascalar potential Φ. The potential for a single point dipole is Φ = Φ d whereΦ d ( r ) = − µ ( r − r (cid:48) ) · m π | r − r (cid:48) | . (13)Here, m is the magnetic moment, r (cid:48) is the position vector of the dipole, and r is theobservation location. The gradient of (13) gives the expected field B ( r ) = µ π | r − r (cid:48) | (cid:20) r − r (cid:48) )( r − r (cid:48) ) · m | r − r (cid:48) | − m (cid:21) . (14) alculation of permanent magnet arrangements for stellarators r, θ, z ) and Cartesian coordinates ( x, y, z ) with the z axis along the axis of the cylinder, and associated unit vectors ( e x , e y , e z ). Again weuse primes to indicate coordinates on the magnet surface, so the evaluation and sourcepositions are r = r cos θ e x + r sin θ e y + z e z , (15) r (cid:48) = b cos θ (cid:48) e x + b sin θ (cid:48) e y + z (cid:48) e z . Supposing the dipoles cover the surface with a uniform number density η (units of1 / area), then the total potential isΦ( r ) = η (cid:90) d r (cid:48) Φ d = − µ ηb π (cid:90) π dθ (cid:48) (cid:90) ∞−∞ dz (cid:48) ( r − r (cid:48) ) · m ( θ (cid:48) ) | r − r (cid:48) | . (16)For the first of the two configurations, shown in figure 10.a, we consider dipoles m = ˆ m cos( nθ (cid:48) ) e x + ˆ m sin( nθ (cid:48) ) e y , (17)for some integer n and constant ˆ m so | m | is uniform. The integral (16) for this case isevaluated in Appendix B, with the resultΦ uni = µ η ˆ m (cid:16) rb (cid:17) n − cos(( n − θ ) (18)for r < b . We will not need the field for r > b . The magnetic field normal to the innersurface is then B n = (cid:18) ∂ Φ ∂r (cid:19) r = a = µ η ˆ m ( n − b (cid:16) ab (cid:17) n − cos(( n − θ ) . (19)Comparing this result to (12), we see the desired field is produced on the plasma surfaceif we choose n = (cid:96) + 1 and ˆ m = ¯ Bbµ η(cid:96) (cid:18) ba (cid:19) (cid:96) − . (20)Expression (19) for a thin layer of dipoles can be extended to a formula for a finite-thickness magnet with inner radius b and outer radius b by writing η ˆ m = M db andintegrating in b over [ b , b ]. The result is B n = µ M ( n − n − (cid:18) ab (cid:19) n − (cid:34) − (cid:18) b b (cid:19) n − (cid:35) cos(( n − θ ) . (21)This result is equivalent to the radial component of (21a) in [21], Halbach’s multipole,noting the following substitutions: N → n − B r → µ M , r , → b , , and ϕ → θ .Equivalently, the magnetization required to produce the field (12) is M = ¯ B ( (cid:96) − µ (cid:96) (cid:18) b a (cid:19) (cid:96) − (cid:34) − (cid:18) b b (cid:19) (cid:96) − (cid:35) − . (22) alculation of permanent magnet arrangements for stellarators B n (12) exactly, then when it isadded to the aforementioned equal and opposite fixed field B n,f = − ¯ B cos (cid:96)θ one obtains f B = 0. Therefore this dipole configuration is a solution of the RECGOIL PM least-squares step in the limit of small λ . Furthermore, since M is uniform, this configurationis a fixed point of the Picard iteration. Therefore this configuration is a fixed point ofthe overall REGCOIL PM algorithm.We can compare this first configuration of dipoles with the second configuration, inwhich the dipole directions are constrained to lie in the direction normal to the surfaces,now allowing M to vary with θ (cid:48) . This second configuration is illustrated in figure 10.b.We assume the magnitude of the dipoles is | m | = ¯ m cos (cid:96)θ (cid:48) for some constant ¯ m , so m = ¯ m cos( (cid:96)θ (cid:48) ) ( e x cos θ (cid:48) + e y sin θ (cid:48) ) . (23)This expression is substituted into (16), and after evaluating the integrals as shown inAppendix B, one finds Φ = Φ nor forΦ nor = µ η ¯ m (cid:96)θ ) (cid:16) rb (cid:17) (cid:96) (24)for r < b . We will not need the field for r > b . The field normal to the plasma surfaceis then B n = (cid:18) ∂ Φ nor ∂r (cid:19) r = a = µ η ¯ m(cid:96) b cos (cid:96)θ (cid:16) ab (cid:17) (cid:96) − . (25)Comparing this expression to (12) it can be seen that the desired field on the plasmasurface is produced if the maximum dipole magnitude is¯ m = 2 ¯ Bbµ η(cid:96) (cid:18) ba (cid:19) (cid:96) − . (26)The fact that the dipole arrangements (17) and (23) can both produce the same field(12) on the plasma reflects the significant freedom available in choosing the magnets fora given stellarator. Comparing (26) to (20), we see that the required maximum dipolemagnitude is twice as large when the dipoles are constrained to lie normal to the magnetsurface compared to the arbitrary-orientation case. We now compare the analytic result (22) to numerical calculations with REGCOIL PM.Since the code is written for toroidal geometry rather than cylindrical geometry, wechoose a very large but finite aspect ratio, with major radius 30 m, b = 1 m, a = 1 / B = 1 T and, initially, (cid:96) = 3. For this section we neglect the Picarditeration to focus on the behavior of the regularized least-squares problem, fixing d = 1mm. We use the following resolution parameters: 96 grid points poloidally, 512 gridpoints toroidally, 24 Fourier modes poloidally, and 2 grid points radially. alculation of permanent magnet arrangements for stellarators f M [A m ]10 f B [ T m ] [T / A ]10 f B [ T m ] (a) (b) Figure 11.
Scan of the regularization parameter λ for the Halbach cylinderbenchmark problem of section 4.2. (a) Pareto trade-off curve. (b) Normal magneticfield error vs λ . First, the behavior of the regularized least-squares solution is examined as theregularization parameter λ is varied. As shown in figure (11), as λ is decreasedbelow 10 − T /A , the normal field error f B can be made arbitrarily small, indicatingthe permanent magnets exactly produce the desired field. This regime correspondsto the vertical part of the Pareto curve in figure 11.a. For λ above this thresholdvalue, the problem becomes over-regularized, with the regularization term forcing themagnetization to be very small such that the dipoles do not significantly cancel the fixedfield. This regime corresponds to the horizontal part of the Pareto curve in figure 11.a.For the rest of this section we focus on values of λ below the threshold, for which themagnet distribution and f M are insensitive to λ , and f B is very small.Next, figures 12.a-b show a comparison of the magnetization computed byREGCOIL PM to the analytic result (22), as (cid:96) or a are varied. In both figures, errorbars are given for the numerical results, displaying ± M as θ (cid:48) isvaried over [0 , π ] and λ is varied over 10 − − − T /A . The error bars are barelyvisible, indicating that M is found to be uniform and independent of the regularization,as it should be. Extremely close agreement is found between the analytic and numericalresults.Finally, figure 13 displays a 3D rendering of the numerical solution for (cid:96) = 3, a = 1 /
5. NCSX example
We now further develop and analyze the NCSX example. In the following subsections,we demonstrate the ability to remove magnet in regions to make room for ports, free-boundary equilibria using the permanent magnets, and the difficult in raising the field alculation of permanent magnet arrangements for stellarators |M| [Amperes / meter] AnalyticREGCOIL_PM -1 a / b = ratio of plasma / coil surface radius |M| [Amperes / meter] AnalyticREGCOIL_PM
Figure 12.
Comparison of the analytic result (22) for a Halbach cylinder toREGCOIL PM numerical calculations, as described in section 4.2.
Figure 13.
The REGCOIL PM procedure can reproduce the Halbach cylindersolution, as described in section 4.2. Black arrows indicate the magnetization. Here, (cid:96) = 3, the plasma surface is shown in red, and the magnet region is shown in green. magnitude. Finally we present a comparison to the different algorithm of ref [14].
It is infeasible to surround the plasma completely with permanent magnets, since accessis required for heating, diagnostics, and maintenance. We therefore now show howregions of the permanent magnets can be removed for ports. The feasibility of includingports in the 0.5 T NCSX configuration was examined previously in [14, 12, 13].Port regions are selected in REGCOIL PM by increasing the local value of the alculation of permanent magnet arrangements for stellarators P o l o i d a l a n g l e Weight w for excluding port regions Figure 14.
The weight function w ( θ, ζ ) used to exclude permanent magnets at thelocations of ports for the NCSX example. weight w in (5). For the example here we choose the following function for the weight: w ( θ, ζ ) =1 + (cid:88) j A j (cid:20) (cid:18) s j (cid:20) − θ j [1 − cos( θ − θ ,j )] (27) − n fp ∆ ζ j [1 − cos( n fp ζ − n fp ζ ,j )] (cid:35)(cid:33)(cid:35) This function is appropriately periodic in the two angles, and the sum over j allowsmultiple ports to be included. Port j is centered at θ = θ ,j and ζ = ζ ,j , while theextent of the ports in θ and ζ is controlled by ∆ θ j and ∆ ζ j . The parameter s j controlsthe sharpness of the transition from w ≈ w (cid:29)
1. In the example here, we chooseport 1 to have θ , = 0 . ζ , = 1 .
7, ∆ θ = 0 .
4, ∆ ζ = 0 . A = 1000, and s = 5.Additional ports are included with the same parameters but at stellator-symmetric and n fp -symmetric locations. These values are chosen to align the ports with the regions oflowest magnet thickness, which are at the outboard side. The resulting w function isshown in figure 14.The REGCOIL PM solution with ports is displayed in figures 15 and 16. Thesame regularization parameter is used as in section 3.3, λ = 10 − T / A . It can beseen that the change to the magnet geometry is minor. A slight thickening of themagnet volume around the edge of the port is apparent. When ports are included, thevolume of permanent magnets increases only slightly, from 2.012 m to 2.025 m . Themaximum B n also increases only slightly, from 0.00299 T without ports to 0.00303 Twith ports. These results indicate it is likely that ports can be included in permanentmagnet stellarators, at least in some locations.A three-dimensional rendering of the REGCOIL PM solution with ports is shown infigure 17. In the magnetization region, arrows with uniform length are drawn everywhereexcept the ports to show the direction of M . alculation of permanent magnet arrangements for stellarators Magnet thickness d [meters] WITHOUT ports Magnet thickness d [meters] WITH ports Figure 15.
Comparison of the magnet thickness computed by REGCOIL PM forthe NCSX example without and with ports. Z [ m e t e r s ] = 0 Z [ m e t e r s ] = (1/4)(2 /3) Z [ m e t e r s ] = (1/2)(2 /3) Outer magnetization surfaceInner magnetization surfacePlasma boundaryPortMagnetization regionPlasma region (a) (b) (c)
Figure 16.
REGCOIL PM solution for the NCSX example with ports.
To evaluate whether a magnet design is adequate, it is necessary to compute the resultingfree-boundary plasma configuration. To this end, figure 18 shows a comparison of theoriginal c09r00 target configuration with the configurations achieved with permanentmagnets. To compute the latter, our REGCOIL PM implementation saves an MGRIDfile that is used as input to free-boundary VMEC [22, 23]. REGCOIL PM results areshown both with and without ports, corresponding to figures 8 and 16. Panels (a)-(b) offigure 18 show that the magnetic axis and flux surface shapes achieved are very close tothose of the target configuration. Panel (c) shows that the rotational transform profileis reproduced accurately as well. Differences between the REGCOIL PM results withand without ports are barely perceptible, indicating again that it should be possible toinclude ports in the design. More detailed analysis must be done to assess whether the alculation of permanent magnet arrangements for stellarators Figure 17.
The REGCOIL PM solution for the NCSX example with ports, viewedfrom two angles. The red surface is the plasma boundary. The inner and outer magnetboundaries I and O are shown, with their local color indicating d . Everywhere in themagnet region except for the ports, arrows of uniform length display the direction of M . small differences in flux surface shape have a meaningful effect on physics properties.Nonetheless, these preliminary results support the idea that producing the 0.5 T NCSXconfiguration with permanent magnets is feasible. alculation of permanent magnet arrangements for stellarators Z [ m e t e r s ] = 0 Z [ m e t e r s ] = /3 s Rotational transform
Fixed-boundary targetFree boundary, no portsFree boundary, with ports (a) (b) (c)
Figure 18.
The REGCOIL PM solutions reproduce the target flux surface shapesand rotational transform.
Since the NCSX example developed in previous sections has a relatively weak magneticfield ∼ . B would require a doubling of the magnet thickness. In fact the thicknessmust be more than doubled, since the new magnet that is introduced compared to the0.5 Tesla case is farther from the plasma and so has less effect. Figure 19 shows theREGCOIL PM solution for the 1 Tesla case with no ports, and a comparison to figure8 makes clear that a significant increase in magnet thickness is indeed required. Themagnetization volume for the 0.5 Tesla case is 2.0 m , compared to 4.9 m for the 1Tesla case. The magnet thickness for the 1 Tesla case is sufficiently large that thecoordinate system in (3) becomes singular, with √ g crossing zero. This issue is specificto the coordinate system we have chosen in regions where the inner surface is concave,and does not necessarily mean a 1 Tesla solution is impossible. However the significantvolume occupied by the magnets in figure 19 suggests that a ≥ It is interesting to compare the results of REGCOIL PM to the topology optimizationmethod described in [14]. The latter approach is implemented in the code FAMUS.In topology optimization, the presence or absence of a magnet at a given location isrepresented by a continuous variable ρ ∈ [0 , ,
1) so ρ ≈ alculation of permanent magnet arrangements for stellarators Z [ m e t e r s ] = 0 Z [ m e t e r s ] = (1/4)(2 /3) Z [ m e t e r s ] = (1/2)(2 /3) Outer magnetization surfaceInner magnetization surfacePlasma boundaryMagnetization regionPlasma region (a) (b) (c)
Figure 19.
If one attempts to raise the mean field magnitude of the NCSX exampleto 1 Tesla, the magnet region becomes significantly thicker (compare to figure 8).
REGCOIL PM approaches are expected to each have advantages and disadvantages.The potential advantages of REGCOIL PM have already been described. Topologyoptimization is more flexible with respect to the magnet geometry, with no restrictionthat all magnets have one fixed surface specified by the user.We carry out a comparison between the two codes for the 0.5 T NCSX case with noports or plasma current. We first obtain a FAMUS solution, considering dipoles allowedto lie within 14 cm of the NCSX vacuum vessel in the direction away from the plasma.The grid of allowed dipole locations has a resolution of 14 points radially, 64 points in θ ,and 384 points in φ (considering all field periods). The level of regularization in FAMUSis set by hand to achieve a plausible solution, with f B = 2 . × − T m . Then λ inREGCOIL PM is adjusted to match this value of f B , with the result λ = 8 . × − T / A . Both codes achieve the same target magnetization M t = 1 . × A / m .An effective volume of the permanent magnet region can be defined in FAMUS by (cid:80) j | m j | /M t where m j are the discrete dipole moments; the result for this case is 2.32m . The permanent magnet volume of the REGCOIL PM solution is slightly lower,1.96 m .The results of the two codes are shown in figure 20. It can be seen that ρ = 0 or1 nearly everywhere in the FAMUS solution. For both codes, black arrows display themagnetization vector’s projection into the ( R, Z )-plane. While | M | is exactly uniformin the REGCOIL PM solution and very nearly uniform in the FAMUS solution, thearrow lengths vary since a φ component may be present. In panels (a)-(c), the M vectors are shown for 4 of the 14 radial grid locations in FAMUS. In REGCOIL PM,where there is no radial variation, only a single arrow is shown. The dipole locationsin FAMUS are shifted from the symmetry planes by half of the grid spacing (i.e. by2 π/
768 radians), so the figures show the nearest planes of dipoles to the given φ . Thereare many similarities between the solutions from the two codes. Both codes yield athicker magnet layer on the small- R side of the plasma. In these thick regions, the alculation of permanent magnet arrangements for stellarators Z [ m e t e r s ] FAMUS, = 0 Z [ m e t e r s ] FAMUS, = (1/4)(2 /3) Z [ m e t e r s ] FAMUS, = (1/2)(2 /3)
Plasma boundaryPlasma regionMagnetization region Z [ m e t e r s ] REGCOIL_PM, = 0 Z [ m e t e r s ] REGCOIL_PM, = (1/4)(2 /3) Z [ m e t e r s ] REGCOIL_PM, = (1/2)(2 /3)
Plasma boundaryPlasma regionMagnetization region M / M t (a)(d) (b)(e) (c)(f) Figure 20.
Comparison between REGCOIL PM and the topology optimization codeFAMUS at matched f B . Black arrows indicate the direction of the magnetization. direction of the magnetization is very similar between the two codes. At the large- R side, the REGCOIL PM solution has a thin magnet layer, whereas FAMUS eliminatesthe magnets in many of these regions. These two different magnet configurations bothproduce a small field error f B , demonstrating again that there is significant flexibilityin the magnet design.
6. Conclusions
In summary, we have demonstrated an algorithm for computing an arrangement ofpermanent magnets outside of a target volume that produces a desired spatially-dependent magnetic field inside the volume. While the algorithm is applied here tostellarators, the method could be used for other applications as well. The methodhere results in a binary magnetization magnitude: at every point M is either zero orequal to a target value M t . This feature is advantageous since any volume occupiedby magnetization of less than the maximum commercially available magnitude is aninefficient use of space. The method also does not place constraints on the directionof M , meaning that Halbach solutions with rotating M are obtained automatically.While we have not rigorously proved stability or existence of a unique fixed point, the alculation of permanent magnet arrangements for stellarators M . This approximation is likelyinaccurate for a serious experimental design. However it appears straightforward toextend the REGCOIL PM algorithm to a more realistic case of discrete magnet blockswith a uniform direction of M in each block. Each block k would be parameterizedwith a thickness parameter d k . The linear-least-squares solve would have three degreesof freedom per block, one for each coordinate of the block’s M vector. The d k parameterof each block could be updated by applying the same fixed-point iteration used here toeach block. This idea will be explored in future work.Even without this extension to discrete magnet blocks, REGCOIL PM could bevaluable as part of optimization of the plasma shape, i.e. the first stage in thestandard two-stage stellarator design. At each iteration of the plasma optimization,REGCOIL PM could be called, and the resulting magnet thickness could be penalizedin the objective function along with other physics quantities. One could thereby findplasma configurations that can be produced with a relatively low volume of permanentmagnets. Inside this optimization, robustness and speed of a code are more importantthan detailed modeling of all engineering factors, and so the ‘smooth’ REGCOIL PM ofthe present paper would be sufficient and well suited. For this application the numberof fixed-point iterations (eq (7)) could be very small, perhaps one, since the magnetthickness need not be precise. Or, the fixed-point iteration could be avoided entirely,and rather the peak magnitude M from the linear-least-squares solution with uniform d could be penalized. Acknowledgments
Input from Steven Cowley, David Gates, Ken Hammond, Per Helander, TonatiuhS´anchez-Vizuet, and Michael Zarnstorff is gratefully acknowledged. This work wassupported by the U.S. Department of Energy under Contract No. DE-AC02-09CH11466.
Appendix A. Jacobian
Here we derive (8). Applying ∂/∂s , ∂/∂θ , and ∂/∂ζ to (3), one finds √ g = ∂ r (cid:48) ∂s · ∂ r (cid:48) ∂θ × ∂ r (cid:48) ∂ζ (A.1)= − σd (cid:20) N − σsd (cid:18) n · ∂ r I ∂θ × ∂ n ∂ζ + n · ∂ n ∂θ × ∂ r I ∂ζ (cid:19) − s d n · ∂ n ∂θ × ∂ n ∂ζ (cid:21) . alculation of permanent magnet arrangements for stellarators HN , as shown in appendix A of [24]. The last term in(A.1) is evaluated by differentiating n = N /N with (4) to obtain ∂ n ∂θ = 1 N (cid:20) ∂ r I ∂θ∂ζ × ∂ r I ∂θ + ∂ r I ∂ζ × ∂ r I ∂θ − n ∂N∂θ (cid:21) , (A.2) ∂ n ∂ζ = 1 N (cid:20) ∂ r I ∂ζ × ∂ r I ∂θ + ∂ r I ∂ζ × ∂ r I ∂θ∂ζ − n ∂N∂ζ (cid:21) . (A.3)Straightforward manipulation then gives (8). Appendix B. Integrals for section 4
Here we derive expressions (18) and (24). We start by inserting (17) (for uniform-magnitude dipoles) or (23) (for dipoles normal to the magnet surface) into (16). Theresults areΦ uni = − µ ηb ˆ m π (cid:90) π dθ (cid:48) (cid:90) ∞−∞ dz (cid:48) ( r cos θ − b cos θ (cid:48) ) cos nθ (cid:48) + ( r sin θ − b sin θ (cid:48) ) sin nθ (cid:48) [( r cos θ − b cos θ (cid:48) ) + ( r sin θ − b sin θ (cid:48) ) + ( z − z (cid:48) ) ] / (B.1)andΦ nor = − µ ηb ¯ m π (cid:90) π dθ (cid:48) (cid:90) ∞−∞ dz (cid:48) [( r cos θ − b cos θ (cid:48) ) cos θ (cid:48) + ( r sin θ − b sin θ (cid:48) ) sin θ (cid:48) ] cos (cid:96)θ (cid:48) [( r cos θ − b cos θ (cid:48) ) + ( r sin θ − b sin θ (cid:48) ) + ( z − z (cid:48) ) ] / (B.2)respectively. The z (cid:48) integrals are evaluated using (cid:82) ∞−∞ dz (cid:48) [ q + ( z − z (cid:48) ) ] − / = 2 /q for q >
0. Changing the remaining integration variable to γ = θ (cid:48) − θ and using an angle-sumtrigonometric identity, one findsΦ uni = − µ η ˆ m π (cid:90) π dγ A r b − rb cos γ (B.3)with A = rb [cos nγ cos(( n − θ ) − sin nγ sin(( n − θ )] (B.4) − cos(( n − γ ) cos(( n − θ ) + sin(( n − γ ) sin(( n − θ )and Φ nor = − µ η ¯ m π (cid:90) π dγ (cid:2) rb cos γ − (cid:3) [cos (cid:96)γ cos (cid:96)θ − sin (cid:96)γ sin (cid:96)θ ]1 + r b − rb cos γ . (B.5)The contributions from terms ∝ sin (cid:96)γ , sin nγ , and sin(( n − γ ) all vanish. Theremaining integrals can be evaluated using (cid:90) π dγ cos nγ − ρ cos γ + ρ = πρ n − ρ for ρ < , n ≥ , (B.6) (cid:90) π dγ cos nγ cos γ − ρ cos γ + ρ = (cid:40) π ρ − ρ ρ n − for ρ < , n ≥ , πρ − ρ for ρ < , n = 0 . (B.7)The results are (18) and (24). alculation of permanent magnet arrangements for stellarators References [1] R Turner. Gradient coil design: a review of methods.
Magnetic Resonance Imaging , 11:903, 1993.[2] M Poole and R Botwell. Novel gradient coils designed using a boundary element method.
Conceptsin Magnetic Resonance B , 31B:162, 2007.[3] S S Hidalgo-Tobon. Theory of gradient coil design methods for magnetic resonance imaging.
Concepts in Magnetic Resonance A , 36A:223, 2010.[4] L Rossi and E Todesco. Electromagnetic design of superconducting quadrupoles.
Phys. Rev. STAccel. Beams , 9:102401, 2006.[5] S. Russenschuck.
Field Computation for Accelerator Magnets: Analytical and Numerical Methodsfor Electromagnetic Design and Optimization . Wiley, 2011.[6] P Merkel. Solution of stellarator boundary value problems with external currents.
Nucl. Fusion ,27:867, 1987.[7] M Drevlak. Automated optimization of stellarator coils.
Fusion Tech. , 33:106, 1998.[8] D J Strickler, L A Berry, and S P Hirshman. Designing coils for compact stellarators.
Fusion Sci.Tech. , 41:107, 2001.[9] M Landreman. An improved current potential method for fast computation of stellarator coilshapes.
Nucl. Fusion , 57:046003, 2017.[10] C Zhu, S H Hudson, Y Song, and Y Wan. New method to design stellarator coils without thewinding surface.
Nucl. Fusion , 58:016008, 2018.[11] P Helander, M Drevlak, M Zarnstorff, and S C Cowley. Stellarators with permanent magnets.
Phys. Rev. Lett. , 124:095001, 2020.[12] C Zhu, M Zarnstorff, D Gates, and A Brooks. Designing stellarators using perpendicularpermanent magnets.
Nucl. Fusion , 60:076106, 2020.[13] K Hammond, C Zhu, T Brown, K Corrigan, D A Gates, and M Sibilia. Geometric concepts forstellarator permanent magnet arrays.
Nucl. Fusion , 60:106010, 2020.[14] C Zhu, K Hammond, T Brown, D, M Zarnstorff, K Corrigan, M Sibilia, and E Feibush. Topologyoptimization of permanent magnets for stellarators.
Nucl. Fusion , 60:106002, 2020.[15] N Pomphrey, L Berry, A Boozer, A Brooks, R E Hatcher, S P Hirshman, L-P Ku, W H Miner,H E Mynick, W Reiersen, D J Strickler, and P M Valanju. Innovations in compact stellaratorcoil design.
Nucl. Fusion , 41:339, 2001.[16] M Landreman. Dataset on Zenodo, http://doi.org/10.5281/zenodo.4029009. 2020.[17] M Landreman. Dataset on Zenodo, http://doi.org/10.5281/zenodo.4028933. 2020.[18] M Drevlak, C D Beidler, J Geiger, P Helander, and Y Turkin. Optimisation of stellarator equilibriawith ROSE.
Nucl. Fusion , 59:016010, 2019.[19] D G Anderson. Iterative procedures for nonlinear integral equations.
J. Assoc. Comput. Mach. ,12:547, 1965.[20] T S´anchez-Vizuet and M E Solano. A hybridizable discontinuous Galerkin solver for theGradShafranov equation.
Comp. Phys. Comm. , 235:120, 2019.[21] K Halbach. Design of permanent multipole magnets with oriented rare Earth cobalt material.
Nucl. Instrum. Meth. , 169:1, 1980.[22] S P Hirshman and J C Whitson. Steepest-descent moment method for three-dimensionalmagnetohydrodynamic equilibria.
Phys. Fluids , 26:3553, 1983.[23] S P Hirshman, W I van Rij, and P Merkel. Three-dimensional free boundary calculations using aspectral Green’s function method.
Comp. Phys. Comm. , 43:143, 1986.[24] M Landreman and E J Paul. Computing local sensitivity and tolerances for stellarator physicsproperties using shape gradients.