Design of quasi-axisymmetric stellarators with variable-thickness perpendicular permanent magnets based on a two-step magnet design strategy

A new method for designing stellarators with variable-thickness perpendicular permanent magnets has been developed, based on a two-step magnet design strategy that we previously proposed for designing stellarators with standardized magnet blocks. Our design strategy uses a ‘local compensation’ method to obtain an initial magnet design and a ‘global fine-tuning’ method to further optimize the initial design. The new method is compared to the previously proposed Fourier decomposition method, and the results indicate that this new method can yield a similar magnet design for an N fp = 2 quasi-axisymmetric stellarator at a lower computational cost. High accuracy is achieved, as demonstrated by a normal magnetic field square on the plasma surface of χB2=3.43×10−8T2m2 , a maximum residual of Bn=3.44Gs for a ∼1 T total field, and a flux-surface-averaged residual B n of Bn/B=3.16×10−5 relative to the total field. Furthermore, the new method can automatically yield a magnet design with large areas of continuous vacancies for ports and plasma access, which is more flexible than the Fourier decomposition method. These results indicate that the new method is robust and effective when used to design stellarators with variable-thickness perpendicular permanent magnets, and, most importantly, that the two-step magnet design strategy has great potential to develop permanent magnet design methods suitable for various magnet layouts.


Introduction
Stellarators continuously attract attention because of their many desirable inherent characteristics, for example, their steady-state operation, absence of a requirement for an external current drive, low recirculating power, freedom from disruption, etc [1,2]. Traditionally, stellarators have relied on external coils to confine and shape the plasma. However, the engineering complexities, in particular, the fabrication and assembly of the extremely complicated three-dimensional non-planar coils [3,4] have significantly restricted the development of stellarators. Recently, Helander et al [5] pointed out that permanent magnets can help shape the plasma and thus simplify stellarator coils. With the use of permanent magnets to create most of the required rotational transform and identical planar coils to create a strong toroidal magnetic field, the engineering difficulties and the construction costs of stellarators can be significantly reduced.
The objective of the permanent magnet design is to reduce the normal magnetic field on the plasma surface to a satisfactory level. The magnetic field, generally created by the coils and the plasma current, will not be too large if an appropriate equilibrium configuration and coil design are used. For example, for a β = 0 ESTELL (Evolutive STELlarator of Lorraine) configuration [6] using 12 identical circular crosssection planar coils, with their coil planes perpendicular to the local magnetic axis to generate a strong toroidal field, the normal magnetic field on the plasma surface is ∼1/5 of the total magnetic field [7] (β = 0 means there is no plasma current). In an experimental device without a blanket, or with a thin blanket, if the magnetic field is not too strong, commercially available magnets are sufficient. The widely used Nd-Fe-B magnets have a remanence B r of up to 1.55 T and a coercivity H ci of up to 5.0 T at liquid nitrogen temperatures [8]. In addition, it is likely that more advanced magnets will become available due to developments in materials science and magnet manufacturing and engineering.
Various methods have been proposed for stellarator magnet design. Helander et al [5] put forward a magnet design method that solved the problem of curl-free magnetization distribution in fixed-thickness permanent magnets, in which the effective surface current could be determined by NESCOIL (NEumann Solver for fields produced by externals COILs) [9]. Caoxiang Zhu et al [10,11] proposed two algorithms, including a linear method [10] that correlated the surface magnetization with the surface current, which could be linearly solved by existing coil-design programs [9,12] and another topology optimization method [11] called FAMUS (Flexbile Advanced Magnets Used for Stellarators) that adopted a 'density method' [13] to polarize the distribution of permanent magnets in the entire design space, based on discretized magnetic dipoles with homogeneous magnetization and free orientation. A study of geometric concepts for magnet arrays based on FAMUS has been performed with the aid of the MAGPIE (Magnets, piecewise) program [14]. In addition, we have previously proposed a method that can design magnets with uniform magnetization and their magnetic axes perpendicular to the plasma surface, based on Fourier decomposition and magnetic surface charge [7]. Landreman [15] subsequently developed a linear least-squares method that exploited the Biot-Savart law's exact linearity in magnetization density and approximate linearity in magnet size for magnets far from the target region.
The magnets obtained using the above methods vary in shape, size, and even remanence and magnetization orientations. Most recently, in order to obtain a magnet design with (1) simple standard structures, (2) uniform magnetization, and (3) finite discrete magnetization orientations for Halbach arrays [16], we have proposed a two-step magnet design strategy [17] based on the idea of the divide and conquer strategy [18]. This strategy uses the 'local compensation' method to yield an initial magnet design and the 'global fine-tuning' method to further optimize the initial design. A method for designing stellarators using identical magnet blocks has been proposed for the first time, and the MASTER (Magnet Arrangement for Stellarator Towards Engineering Realization) program has been developed based on the proposed design strategy. Halbach solutions using identical magnet blocks can greatly simplify the fabrication, assembly, and maintenance of the magnets, and thus significantly reduce stellarator construction costs.
In order to demonstrate the universality and effectiveness of the two-step magnet design strategy for developing magnet design methods suitable for other magnet layouts, in this paper, we generalize the strategy to a new method for designing stellarators with variable-thickness perpendicular permanent magnets. In this new method, the magnets have the same remanence B r and are discretized into thin magnet pieces that can be approximated by magnetic dipoles at their centers. The magnet pieces are stacked together layer by layer, starting from the winding surface, with their magnetic axes perpendicular to the plasma surface. The magnetization orientations of the magnet pieces in each radial column are kept the same as each other. In the design process, the magnet thickness, i.e. the number of magnet pieces in each radial column changes dynamically over a number of iterations until the accuracy requirement is satisfied. In addition, a smoothing procedure is applied to avoid excessive local magnet thickness.
The new strategy also makes use of the 'local compensation' method and the 'global fine-tuning' method. In the 'local compensation' method, the decision to add, cancel, or retain a piece of magnet depends on the outcome that obtains the minimum absolute value of the local normal magnetic field, |B n (p)|, on the plasma surface. The location p is the point where the absolute value of the normal magnetic field on the plasma surface generated by the magnet piece under evaluation reaches a maximum. A truncation threshold Bthreshold is applied to prevent the method from focusing on reducing the residual field in locations where it is already sufficiently weak. With an appropriate threshold value, large areas of continuous vacancies appear, which can be used for ports and plasma access. The procedure of the 'global fine-tuning' method is similar to that of the 'local compensation' method, but the magnet design criterion is replaced with a decision about which choice will result in the minimum surface integral of the normal magnetic field square B 2 n on the plasma surface, χ 2 B . The vacancies given by the 'local compensation' method can either be excluded from the 'global fine-tuning' method to reserve ports and plasma access, or they can be included in the 'global fine-tuning' method for better B n compensation. The new method is successfully applied to the magnet design of a N fp = 2 quasi-axisymmetric stellarator with strong toroidal magnetic fields generated by 12 identical circular planar coils [7]. The results indicate that the two-step magnet design strategy has great potential to develop permanent magnet design methods suitable for various magnet layouts. This paper is organized as follows. The method is described in section 2. Section 3 presents the numerical validation and the benchmarks of this method and the previously proposed Fourier decomposition method. Section 4 includes a magnet design with ports and a check of the method's accuracy in generating a target magnetic field using field-line tracing. The summary and a discussion are presented in section 5.

Magnet design target
In a permanent-magnet stellarator, the plasma is generally confined in the magnetic field generated by the plasma current B plasma , the planar coils B coil , and the permanent magnets B magnet . Moreover, the plasma boundary, i.e. the plasma surface, must be a magnetic surface [9]. Therefore, the objective of the magnet design is to find a magnet arrangement that can eliminate the normal magnetic field created by the plasma current and the coils on the plasma surface [5], i.e. one that can make where n is the normal vector of the plasma surface. The surface integral of the normal magnetic field square on the plasma surface, χ 2 B , is usually used to evaluate the global normal magnetic field on the plasma surface: where ds is asurface panel of the plasma surface.

Magnet layout and magnetic field calculation
Magnets with poloidal and toroidal resolutions of M × L are mounted on the winding surface with their magnetic axes perpendicular to the plasma surface. The winding surface encloses the plasma at a constant distance d from the plasma surface. Each radial column of magnets is discretized into thin magnet pieces with equal thicknesses of dh. The vertex coordinates of the magnet pieces are extended, starting from the plasma surface, along the corresponding normal vector of the plasma surface; for example, the vertex coordinates of the inner and outer surfaces of the ith-layer magnet pieces are calculated by equation (3): where C inner piece and C outer piece are the coordinates of the inner and outer surfaces of the magnet pieces and C plasma represents the coordinates of the plasma surface. The magnetization orientations of the multi-layer magnet pieces in each radial column are the same, i.e. the part pointing away from the plasma is defined as positive (+), and the part pointing toward the plasma is defined as negative (−). The number of magnet pieces, i.e. the magnet thickness, in each radial column is flexible and can be different from those of the neighboring columns. The design target is to find the appropriate magnet thickness distribution to generate the required magnetic field. Once the magnet design is determined, the magnet pieces in each radial column are considered to be one single magnet block, due to their identical remanence and magnetization orientations. The magnet assembly formation is the same as that proposed in the Fourier decomposition method [7], and the permanent magnets can easily be inserted into the cells of a gridded frame attached to the winding surface and fixed from the back using springs and caps from the back. A sketch of the magnet layout with a relatively large dh in the toroidal scope 0, π/6 is shown in figure 1. The poloidal and toroidal resolutions are 120 × 240. Here, h ± is the magnet thickness; a positive value (red) means that the magnetization orientations of the magnet pieces point away from the plasma and a negative value (blue) means that they point towards the plasma. There are no vacancies among the multi-layer magnet pieces of each radial column, which can make full use of the design space and reduce magnet consumption. An empty radial column indicates that there are no magnets at this location. The thin magnet pieces can be approximated by magnetic dipoles at their centers. The magnetic field generated by the magnetic pieces is calculated using the Biot-Savart law [11]. If we suppose that N p discrete magnet pieces are used to generate the required magnetic field, the calculation formula is as shown in equation (4): where r i is the positional vector between the evaluated field point and the ith magnet piece, |m i | = B r V i /μ 0 is the corresponding magnetic dipole moment, B r is the remanence, and V i is the volume. The normal magnetic field created by the permanent magnets on the plasma surface is the linear superposition of that generated by each magnet piece: where B i piecen is the normal magnetic field generated by the ith magnet piece on the plasma surface.

Magnet design method
The two-step magnet design strategy uses the 'local compensation' method to yield an initial magnet design and the 'global fine-tuning' method to further optimize the initial design, which decomposes the design process of large number of magnet pieces into independent designs for each magnet piece and then iterates to find the solution [17]. In this paper, when we generalize the strategy to the new design, the 'local compensation' method and the 'global fine-tuning' method are both iteration procedures conducted on the last magnet pieces of each radial column according to the poloidal and toroidal resolutions. When there is no magnet in the corresponding radial column, we take the first vacancy in the radial column as a virtual magnet piece. The number of magnet pieces, i.e. the magnet thickness, of each radial column changes dynamically with the iterations until the accuracy requirements are achieved. When a magnet piece is under evaluation, there are three states for each magnet piece to choose: (1) vacancy, (2) existing state with the magnetization orientation pointing away from the plasma and (3) existing state with the magnetization orientation pointing toward the plasma. The three states are labelled '0,' '+1,' and '−1' in the following description.

Smoothing operation.
A smoothing operation is applied when a magnet piece needs to be added to the end of a radial column to prevent the local magnet thickness from being too large, i.e. the maximum layer difference dl between the magnet pieces in each radial column under evaluation and those of its nearby radial columns is limited. This smoothing operation can also prevent the magnet thickness from fluctuating sharply within iterations. As shown in figure 2, a 3 (poloidal) × 3 (toroidal) smoothing window (labelled with a red box) is given as an example. The mean thickness, excluding the radial column under evaluation and the radial column without magnets in the smoothing window, is (|−3| When one magnet piece needs to be added at the end of the radial column under evaluation (labelled with the yellow box), the criterion (2 + 1) − 1.857 dl needs to be satisfied. In order to use as few magnets as possible, the magnet smoothing operation is not applied when one magnet piece needs to be cancelled. In addition, because the magnets are installed perpendicular to the plasma surface, when we need to design ports, it is difficult to reduce the normal magnetic When the radial column under evaluation (labelled with the yellow box) needs to have one magnet piece added at its end, the criterion (2 + 1) − 1.625 dl needs to be satisfied. field in the port region on the plasma surface. Moreover, the magnet thickness distribution has a sharp change from the port region to the body region; this requires the smoothing operation limitation to be removed near the port region for more flexible adjustment of the magnet thickness and thus for better B n compensation. In other words, dl near the port region needs be defined as larger than that in the body region.

2.3.2.
The local compensation method. The 'local compensation' method locally evaluates the influence that the magnet pieces exert on B n . Because the intensity of the magnetic field generated by the magnet piece decays with the third power of the distance to the field point, it is reasonable to focus on the local influence that each magnet piece has on the normal magnetic field on the plasma surface [17]. What we need to do for the magnet piece under evaluation is to determine which state can obtain the minimum absolute value of the local normal magnetic field, |B n (p)|, on the plasma surface, on the premise of meeting the magnet layout limitation. The 'local' p is the point where the absolute value of the normal magnetic field on the plasma surface generated by the magnet piece under evaluation reaches a maximum. For example, as shown in figure 3, the normal magnetic field on the plasma surface created by magnet pieces with magnetization orientations pointing away from the plasma and toward the plasma, B + piecen and B − piecen , show clear peak and valley patterns that are located at p: (60, 120), where the poloidal and toroidal resolutions of the plasma surface are 120 × 240. The peak and valley pattern has a major influence on B n .
In addition, the truncation threshold B threshold is applied to prevent the method from focusing on reducing an already small |B n (p)|; i.e. when |B n (p)| B threshold , no more magnets are added to the radial column under evaluation. By doing so, the 'local compensation' method avoids falling into the senseless infinite optimization of a sufficiently weak normal magnetic field. Of course, the truncation threshold can also be abandoned to obtain better B n compensation. With an appropriate B threshold , large areas of continuous vacancies that can used for ports and plasma access appear, which are essential for stellarator design. When the local optimization target of each magnet piece has been achieved, the absolute value of the normal magnetic field |B n | all over the plasma surface is reduced to a suitable level, and thus χ 2 B is also automatically reduced . Based on all the definitions and explanations given above, the detailed iteration procedures for the 'local compensation' are organized using the following three steps and the corresponding simplified iterative flow chart is presented in figure 4.
Step 1: scan the last magnet piece of each radial column according to the poloidal and toroidal resolutions. Note that B n is updated each time the magnets change. Two cases govern how the state of each magnetic piece is determined: (a) The initial state of the magnet piece under evaluation is '0', which means that the magnet piece is the first vacancy in the corresponding radial column.
1. If the initial absolute value of the local normal magnetic field is not larger than the truncation threshold, i.e.|B n (p)| B threshold , the state of the magnet piece remains '0'. Then turn to the next magnet piece. 2. If the initial absolute value of the local normal magnetic field is larger than the truncation threshold, i.e. |B n (p)| > B threshold , the other states '+1' and '−1' will be evaluated. Because there is no limitation for the first-layer magnet piece, from a choice of '0,' '+1,' and '−1,' the state that can obtain the minimum absolute value of the corresponding local normal magnetic field |B n (p)| is assigned to the magnetic piece under evaluation. If the adopted state is not '0' and the updated |B n (p)| is still larger than B threshold , one new magnet piece added behind it will be evaluated immediately and only the adopted state and '0' will be taken into consideration. In the process of magnet addition, the smoothing operation must be satisfied. When this step ends, turn to the next magnet piece.
(a) The initial state of the magnet piece under evaluation is '+1' or '−1.' First set its state to '0' and update B n . Two subcases are used to determine its specific state.
1. If the absolute value of the updated local normal magnetic field is not larger than the truncation threshold, i.e.|B n (p)| B threshold or if, from a choice of '0,' '+1,' and '−1,' the state that can obtain the minimum absolute value of the corresponding local normal magnetic field |B n (p)| is '0,' '0' is adopted as the state of the magnet piece under evaluation. Then turn to the next magnet piece. 2. If the absolute value of the updated local normal magnetic field is larger than the truncation threshold, i.e. |B n (p)| > B threshold and from a choice of '0,' '+1,' and '−1,' the state that can obtain the minimum absolute value of the corresponding local normal magnetic field |B n (p)| is not '0': (a) If the magnet piece under evaluation is the first layer in the corresponding radial column, because there is no limitation for the first layer magnet piece, the state that can obtain the minimum absolute value of the corresponding local normal magnetic field |B n (p)| is adopted as the state of the magnet piece under evaluation. If the updated |B n (p)| is still larger than B threshold , one new magnet piece added behind it will be evaluated immediately and only the adopted state and '0' will be taken into consideration.
In the magnet addition, the smoothing operation must be satisfied. Then turn to the next magnet piece. (b) If the magnet piece under evaluation is not the first layer in the corresponding radial column and the state that can obtain the minimum absolute value of the corresponding local normal magnetic field |B n (p)| is the same as that of its adjacent inner magnet piece, this state is adopted as the state of the magnet piece under evaluation. If the updated |B n (p)| is still larger than B threshold , one new magnet piece added behind it will be evaluated immediately and only the adopted state and '0' will be taken into consideration. In the process of magnet addition, the smoothing operation must be satisfied. Then turn to the next magnet piece. (c) If the magnet piece under evaluation is not the first layer in the corresponding radial column and the state that can obtain the minimum absolute value of the corresponding local normal magnetic field |B n (p)| is the inverse of the state of its adjacent inner magnet piece, '0' is adopted as the state of the magnet piece under evaluation. The adjacent inner magnet piece is then evaluated immediately. If the adjacent inner magnet piece is the first layer, all the states '0,' '+1,' and '−1' will be taken into consideration, otherwise, only '0' and its initial state will be taken into consideration. Then turn to the next magnet piece.
Step 2: after scanning the last magnet pieces of each radial column according to the poloidal and toroidal resolution once, calculate the surface integral of the total normal magnetic field on the plasma surface, χ 2 B = B 2 n ds.
Step 3: repeat steps 1 and 2 until χ 2 B reaches convergence. In addition, it is possible to assign a maximum magnet thickness limitation in the above procedure (or not). Generally, the normal magnetic field on the plasma surface can be greatly reduced using the 'local compensation' method.

2.3.3.
The global fine-tuning method. The 'global finetuning' method globally evaluates the influence that the magnet pieces have on B n . Using χ 2 B as the optimization target with which to evaluate the global normal field on the plasma surface, the 'global fine-tuning' method is a good supplement to the 'local compensation' method. When the 'local compensation' method reaches convergence, the results are delivered to the 'global fine-tuning' method as an input to further optimize the magnet distribution. The optimization procedure of the 'global fine-tuning' method is the same as that of the 'local compensation' method, but the magnet design criterion is replaced by a consideration of which state of the magnet piece under evaluation can obtain the minimum χ 2 B . In addition, it is also possible to assign a maximum magnet thickness limitation in the design procedure (or not). Note that the 'global fine-tuning' method is just a supplement, most of the magnet design work is still undertaken by the 'local compensation' method. Therefore, it is possible to assign the magnets involved in the 'global fine-tuning' method. For example, if the truncation threshold in the 'local compensation' method is used to design ports and plasma access, the vacancies produced by the 'local compensation' method will not be involved in the 'global fine-tuning' method. When the 'global fine-tuning' method is used to further optimize the magnet distribution, the whole normal field on the plasma surface evaluated by χ 2 B can be greatly reduced.

Numerical validation
In order to validate the new method for designing stellarators with variable-thickness perpendicular magnets, we introduce a two-period (N fp = 2) and β = 0 quasi-axisymmetric equilibrium for the ESTELL stellarator [6]. Here, β = 0 indicates that no magnetic field is produced by the plasma current, B plasma = 0. The major radius of the magnetic axis R 0 = 1.4 m and the field strength on the magnetic axis B 0 = 1 T. A strong toroidal magnetic field is generated by 12 uniformly distributed identical planar coils whose coil planes are perpendicular to the local magnetic axis. Each coil has 15 × 9 turns, and they all carry the same current I coil = 4.4 kA [7]. The normal magnetic field on the plasma surface generated by the coils is required to be reduced by the magnets. Figure 5 shows the planar coils and the generated normal magnetic field on the plasma surface. One of the most widely commercially available neodymium-iron-boron (Nd-Fe-B) magnets with a remanence of 1.4 T at room temperature [8] is chosen for the following magnet design, and the thickness of each magnet piece dh = 2.0 mm in the following calculation. The winding surface remains at a constant distance d = 10 cm from the plasma surface. The spatial resolutions of the plasma surface and the winding surface for magnet assembly are both 120 (poloidal) × 240 (toroidal).

Magnet design validation
To enable fine adjustment of the magnet thickness distribution, the thickness of each magnet piece is set to dh = 2.0 mm. In order to reduce the normal magnetic field on the plasma surface as much as possible, the truncation threshold in the 'local compensation' method is defined as B threshold = 0 Gs. The vacancies given by the 'local compensation' method are also involved in the 'global fine-tuning' method for further optimization. In addition, in order to obtain a smoothed magnet distribution, the smoothing window chosen for the magnet addition process is 5 (poloidal) × 5 (toroidal) and the maximum layer difference dl = 2. As shown in figure 6 (a), after 108 iterations (requiring ∼1.13 h) of the 'local compensation' method carried out on 90 CPUs at 3.9 GHz, the surface integral of the total normal magnetic field on the plasma surface χ 2 B , the maximum normal magnetic field B max n , and the final flux-surface-averaged residual B n relative to the total field |B n /B| converge to 5.78 × 10 −8 T 2 m 2 , 5.9 Gs, and 3.87 × 10 −5 , respectively. Then, after 28 iterations (requiring ∼0.45 h) of the 'global fine-tuning' method carried out on 86 CPUs at 3.9 GHz, χ 2 B , B max n , and |B n /B| are further optimized to 4.21 × 10 −8 T 2 m 2 , 2.9 Gs and 3.68 × 10 −5 , respectively. The black dashed line indicates the iteration boundary between the 'local compensation' method and the 'global finetuning' method. The final residual normal magnetic field on the plasma surface is shown in figure 6 (b). This new method shows great performance in the design of variable-thickness perpendicular permanent magnets for stellarators.

Benchmark using the Fourier decomposition method
Because they share the same magnet arrangements, it is meaningful to compare the present method with the previously introduced Fourier decomposition method [7]. For the Fourier decomposition method, the final χ 2 B , B max n , and |B n /B| are 4.79 × 10 −9 T 2 m 2 , 2.0 Gs, and 1.31 × 10 −5 , which are slightly smaller than the results of 5.78 × 10 −8 T 2 m 2 , 5.9 Gs, and 3.87 × 10 −5 obtained by the present method. However, the difference is very small, which indicates that the present method is effective and reliable. Moreover, the new method involves a much smaller computational cost than that of the Fourier decomposition method. Contour plots of the magnet thickness h ± on the winding surface obtained by the two methods are shown in figures 7 (a) and (c). It is obvious that the two distributions are almost the same as each other. In addition, Figures 7 (b) and (d) show that the 3D magnet distributions in the toroidal scope 0, π/6 also exhibit the same appearance. The main difference between the magnet distributions obtained by the two methods is the distribution smoothness. The magnet design obtained by the Fourier decomposition method is a smooth and continuous magnet distribution, but that obtained by the new method is a discrete distribution. One of the main reasons for this is that the new method cannot adjust the magnet thickness more delicately because of the fixed magnet thickness dh of each magnet piece. Moreover, in order to reduce magnet consumption as much as possible, we do not use the smoothing operation when cancelling one magnet piece, which also leads to the smoothness difference. The vacancies appearing in figure 7 (d) provide direct evidence of this. In conclusion, although there are some differences, the benchmark indicates that the results obtained by the new method fit those obtained by the Fourier decomposition method quite well, which proves the excellent performance of the new magnet design method and validates the universality of the two-step magnet design strategy.

Smoothing operation
The smoothing operation is very effective at preventing the local magnets from being too thick. As shown in figure 8, the magnet distribution designed with a 3 (poloidal) × 3 (toroidal) smoothing window and a maximum layer difference of dl = 15 appears to be more uneven than that in figures 7 (c) and (d), and the magnet thickness shows local peak patterns. However, in this design, χ 2 B , B max n , and |B n /B| are finally optimized to 3.21 × 10 −8 T 2 m 2 , 2.28 Gs, and 3. 11 × 10 −5 , respectively, which is a better result than those of figures 7 (c) and (d) and closer to that obtained using the Fourier decomposition method. Therefore, it is necessary to balance the magnet distribution and the performance of the final residual normal magnetic field on the plasma surface.

Magnet design with ports
Whatever form the magnet arrangement takes, the availability of ports for plasma heating and diagnostics is a basic requirement for a stellarator. As shown in figures 7 (d) and 8 (b), some vacancies or thin magnets appear in the junctions between magnets with opposite magnetization orientations. These vacancies and thin magnets correspond to normal magnetic fields with opposite directions on the plasma surface, which are also relatively weak normal fields. Therefore, applying an appropriate value of B threshold in the 'local compensation' method is a good way to design ports in these places. When the local normal magnetic field on the plasma surface is smaller than the threshold, i.e. |B n (p)| B threshold , no more magnets are added in the radial column under evaluation. Under this condition, places on the winding surface corresponding to points where the initial normal magnetic fields on the plasma surface are already sufficiently weak will be designed without magnets. These vacancies can then be actively reserved for ports in the 'global fine-tuning' method.   Threshold scan to find the optimal magnet design with ports; the solid lines with triangles are the magnet designs produced by the 'local compensation' method, and the dashed lines with dots are the further optimized magnet designs given by the 'global fine-tuning' method. Figure 10. The maximum and average widths of the ports, which are sampled at toroidal angles of φ = 0, π/12, π/6, π/4, π/3, 5π/12, π/2 .

Threshold scan
In order to design suitable ports, a scan of the truncation threshold B threshold is required. In the scan process, the smoothing window chosen is still 5 (poloidal) × 5 (toroidal) and the maximum layer difference dl body = 2 for the body region and dl port = 5 near the port region. The final χ 2 B , B max n , |B n /B| and the total areas of the vacancies for the ports S port of the magnet designs obtained with different thresholds are given in figure 9. The solid lines with triangles represent the magnet designs produced by the 'local compensation' method and the dashed lines with dots are the further optimized magnet designs provided by the 'global fine-tuning' method. It is easy to find that regardless of the threshold value, the values of χ 2 B , B max n , and |B n /B| further optimized by the 'global fine-tuning' method are much smaller than those given by the 'local compensation' method, but S port is almost the same before and after the further optimization, which means that the 'global fine-tuning' method indeed complements the 'local compensation' method in the design of better magnet arrangements. In addition, χ 2 B , B max n , |B n /B| , and S port increase with an increase of B threshold , indicating that the magnet designs obtained with larger thresholds lead to worse performance in terms of B n compensation, but yield larger port areas. S port has a steep increase near B threshold = 2-4 Gs and tends to converge when B threshold > 4 Gs. To determine the appropriate magnet scheme, it is necessary to evaluate the final residual normal magnetic field on the plasma surface as well as S port .
The ports are generally band-shaped in the magnet design obtained with the proposed method, as shown in figure 12 (d). Therefore, in addition to the total vacancy area of the ports S port , the change trend of the port width with different values of B threshold is a more intuitive parameter with which to evaluate the magnet design. As shown in figure 10, the maximum and average port widths sampled at the toroidal angles φ = 0, π/12, π/6, π/4, π/3, 5π/12, π/2 present similar trends to that of S port . For B threshold = 18 Gs, the maximum port width is up to 0.161 m.
The computation times of the local compensation method and the global fine-tuning method for different values of B threshold are shown in figure 11. The computation time of the local compensation method first decreases with increasing values of B threshold and then converges when B threshold > 4 Gs. The computation time of the global fine-tuning method increases with increasing values of B threshold , because the larger B threshold used in the local compensation method leads to worse performance in terms of B n compensation, and thus the global finetuning method needs more iterations to further optimize the magnet design. Because the global fine-tuning method needs to compute the surface integral of B n for each iteration, the computational cost of each iteration is larger than that of the local compensation method. Therefore, the total computation time first decreases with increasing values of B threshold and then increases when B threshold > 6 Gs.

Port design
According to the threshold scan results given above, in order to obtain a magnet arrangement that yields a relatively larger port area while generating the required magnetic field, B threshold = 14 Gs is applied for the magnet design, since S port and the port width increase slowly when B threshold > 14 Gs. As shown in figure 12 (a), after 77 iterations of the 'local compensation' method, χ 2 B , B max n , and |B n /B| converge to 2.43 × 10 −5 T 2 m 2 , 21.0 Gs, and 1.19 × 10 −3 , respectively. After 71 iterations of the 'global fine-tuning' method, χ 2 B , B max n , and |B n /B| are further optimized to 3.62 × 10 −7 T 2 m 2 , 9.0 Gs,  figure 12 (b); it is obvious and reasonable that the relatively stronger residual normal fields are located in the regions on the plasma surface that correspond to the port regions on the winding surface. In addition, as shown in figure 12 (c), the field-line tracing result obtained for the magnetic field calculated from the predesigned planar coils and the magnets, B = B coil + B magnet , indicates that the flux surface matches the target equilibrium magnetic flux surfaces quite well. The relative rotational transform difference at the plasma surface is Δι/ι = 1.1 %. The field-line tracing is carried out in the cylindrical coordinate system with an improved sixth-order Runge-Kutta algorithm [19]. Most importantly, as shown in figure 12 (d), four band-shaped continuous areas with continuous vacancies for ports automatically appear, distributed in the toroidal direction at the junctions between magnets with opposite magnetization orientations. In figure 10, the maximum port width is up to 0.152 m. The overall design is shown in figure 13; the total area of the vacancies on the winding surface for ports is S port = 3.71 m 2 and the total magnet volume isV magnet = 0.86 m 3 .

Summary and discussion
In this paper, we have generalized the two-step magnet design strategy [17] to a new method for designing variablethickness perpendicular permanent magnets for stellarators, which validates the universality of the two-step magnet design strategy. The field line tracing results indicate that the magnet design obtained with the new method, coupled with the planar coils, can perfectly generate the required magnetic equilibrium. Although we only show the designs for a quasiaxisymmetric stellarator, the new method can be applied to any configuration.
In the new method, the magnets are discretized into thin magnet pieces that can be approximated by magnetic dipoles at their centers and stacked layer by layer, starting from the winding surface and with their magnetic axes perpendicular to the plasma surface. No vacancies appear among the multilayer magnet pieces of each radial column. In addition, all the magnets are defined using the same remanence B r , and all the magnet pieces in each radial column have the same magnetization orientation. Once the magnet design is determined, the magnet pieces in each radial column are considered to be a single magnet block. The two-step magnet design strategy includes the 'local compensation' method, which yields the initial magnet design and the 'global fine-tuning' method used to further optimize the initial design. Therefore, the new method is also a two-step procedure. The initial design given by the 'local compensation' method is delivered to the 'global fine-tuning' method as an input for further optimization. For both steps, the new method iteratively scans the last magnet piece of each radial column according to the poloidal and toroidal resolution, and the number of magnet pieces, i.e. the magnet thickness, of each radial column fluctuates over a number of iterations until the corresponding convergences are achieved. The new method has two features, as follows: (a) The truncation threshold B threshold used in the 'local compensation' method shows excellent performance when used to design ports for the perpendicular magnet layout. The magnet thickness pattern is highly similar to that of the normal magnetic field that needs to be created by magnets on the plasma surface. If the algorithm is supplied with an appropriate threshold value, the regions on the winding surface that correspond to the relatively weak initial normal field on the plasma surface will not be designed with magnets, and thus large areas of continuous vacancies appear that can be used for ports and plasma access. Larger threshold values lead to larger port areas but worse performance in terms of B n compensation. Therefore, a scan of the threshold is required to determine an appropriate value. (b) A smoothing operation is introduced to prevent the local magnet thickness from being too large, i.e. the maximum layer difference dl between the magnet pieces in each radial column under evaluation and those in the nearby radial columns is limited. The smoothing operation can also avoid sharp fluctuations of the magnet thickness within iterations. In addition, the maximum layer difference dl near the port region can be defined to be larger than that in the body region to in order to remove the thickness limitation and obtain better B n compensation in the port region. In order to use as few magnets as possible, the smoothing operation is only applied in the magnet addition process. The magnet designs obtained with different smoothing windows and dl values show different magnet thickness distributions and thus lead to different residual normal magnetic field distributions on the plasma surface. Therefore, it is necessary to balance the magnet thickness distribution and the residual normal field on the plasma surface.
This new method was compared to the previously proposed Fourier decomposition method [7] on a two-period (N fp = 2) and β = 0 quasi-axisymmetric equilibrium for the ESTELL stellarator [6]. The results indicate that this new method can yield almost the same magnet design as that produced by the Fourier decomposition method. Compared to the Fourier decomposition method, the new method has three main advantages, as follows.
(a) The new method is more flexible. In the Fourier decomposition method, once the target magnetic field is determined, the final magnet distribution is then uniquely determined unless the method is modified to be suitable for other fixed magnet layouts. However, for this new method, different smoothing operations, truncation thresholds B threshold used in the 'local compensation' method, and even the thicknesses of individual magnet pieces can give different magnet designs, which may help us to discover more possibilities. In particular, with an appropriate B threshold , the automatic appearance of continuous vacancies is helpful for port and plasma access design. (b) Although the design accuracy of the new method is lower than that of the Fourier decomposition method, the payoff is absolutely worthwhile. The magnet design obtained with the Fourier method has a smooth and continuous magnet distribution. If we take this design into engineering, the reduced accuracy generated by the discretization of magnets is not absolutely controllable. This new method has been developed using discrete magnet pieces, which can avoid this problem to some extent. (c) The new method needs fewer computational resources.
The Fourier decomposition needs to calculate many tiny Fourier components to make fine adjustments in the magnet thicknesses; therefore, it inevitably consumes many computational resources. However, in the new method, with the exceptions of the poloidal and toroidal resolutions, the new method has a radial resolution, i.e. the thickness of individual magnet pieces. In this situation, the new method sacrifices the accuracy of the thickness adjustment and thus reduces the demand for computational resources.
In conclusion, the two-step magnet design strategy has been successfully generalized to the design of variable-thickness perpendicular permanent magnets for stellarators, which demonstrates that the two-step magnet design strategy has great potential to develop permanent magnet design methods that are suitable for various magnet layouts. An efficient, robust, and easy to implement design strategy may promote the development of permanent-magnet stellarators.