A topology optimization algorithm for magnetic structures based on a hybrid FEM–BEM method utilizing the adjoint approach

A method to optimize the topology of hard as well as soft magnetic structures is implemented using the density approach for topology optimization. The stray field computation is performed by a hybrid finite element–boundary element method. Utilizing the adjoint approach the gradients necessary to perform the optimization can be calculated very efficiently. We derive the gradients using a “first optimize then discretize” scheme. Within this scheme, the stray field operator is self-adjoint allowing to solve the adjoint equation by the same means as the stray field calculation. The capabilities of the method are showcased by optimizing the topology of hard as well as soft magnetic thin film structures and the results are verified by comparison with an analytical solution.

In many applications like magnetic sensors, magnetic storage, electric machines and so on, magnetic fields with specific spatial profiles are exploited in order to achieve a specific goal. The properties of magnetic fields are defined by the properties of the involved magnetic structures. Most commonly, in order to obtain a certain magnetic field, the geometry of a magnetic structure is optimized. This can e.g. be the geometry of a permanent magnet producing a well defined magnetic field, or a soft magnetic structure shaping the field possibly generated by an electric current. The process of optimizing the geometry of these structures can be very time consuming. A common approach is called geometry parameterization where the structure's geometry is parameterized with relatively few parameters and the parameter region is investigated until an optimal solution is found. Geometry parameterization is one representative of a group of techniques summarized as shape optimization 1 . Another approach is to directly optimize the distribution of material within a given domain. This approach is called topology optimization 2 . It drastically increases the degrees of freedom (dofs) but allows for possibly better, geometric realizations of the considered system, independent of any introduced parameterizations. Successful applications of topology optimization of magnetic structures include the optimization of magnetic recording write heads 3 , of rotors of brushless DC motors 4 and of electromagnetic actuators 5,6 .
In this paper a method to optimize the topology of magnetic structures is presented. The presented method is based on a highly efficient hybrid finite element-boundary element method (FEM-BEM) approach 7 solving the magnetostatic Maxwell's problem. In contrast to already presented finite element 3,8 and finite differences 9 algorithms using a FEM-BEM approach has the advantage that only the regions of interest need to be discretized reducing the dofs dramatically. Furthermore, in order to efficiently calculate gradients during optimization, the adjoint approach is utilized. Our method takes the B/H curve of a hard or soft magnetic material linearized at the working point as input and is able to include external fields generated by electric coils. After optimization, the result is an optimized magnetic structure and its stray field dependent magnetization state.
The algorithm is implemented on the basis of magnum.fe 10 which in turn utilizes the finite element library FEniCS 11 . In particular FEniCS is used for the definition of the objective functionals and to automatically differentiate the involved partial derivatives.
The paper is structured as follows. In the first section the solution of the forward problem-the calculation of the stray field of a magnetic structure-is presented. In the following section the density approach to topology optimization is introduced. After that the gradient necessary to perform the optimization and the adjoint equation are introduced. Thereafter details of the numerical optimization are discussed and numerical experiments www.nature.com/scientificreports/ showcasing the capabilities of the method are presented. A detailed derivation of the gradient using the adjoint approach within a "first optimize then discretize" scheme is given in the appendix.

Forward problem
The stray field of a magnetic body consisting of a magnetic material with linear magnetic susceptibility χ and possibly a remanence magnetization M r within an external magnetic field H ext , can be calculated starting from the magnetostatic Maxwell's equation. Using the reduced scalar potential formulation the total magnetic field H is split into an external part H ext (created via currents outside of m ) and the curl-free induced magnetic field H d . Expressing the induced field as gradient of a scalar potential u, H d = −∇u finally yields In order to use the hybrid FEM-BEM approach, next R 3 is divided into a region m containing all magnetic material and a region R 3 \ m . Outside m Eq. (1) then reduces to and the following boundary and jump conditions apply where n is the outward pointing normal vector and [·] ∂� denotes the jump over the boundary ∂� . The resulting system is solved using a hybrid FEM-BEM algorithm 7 .

Topology optimization
Topology optimization tries to find, within a predefined spatial domain opt , the optimal material distribution with respect to a design goal. This can be achieved using the density approach for topology optimization 12 . Here, a scalar indicator function is introduced, transforming the remanence magnetization M r → M r (ρ) = ρ p M r as well as the susceptibility χ → χ(ρ) = ρ p χ . Here, p is a parameter introduced in order to penalize intermediate values of 0 ≤ ρ ≤ 1 that are allowed during optimization 12 . However, it turns out that the impact of p is more complex as is investigated for each presented numerical experiment (see "Numerical experiments" section). The forward problem then reads where u is an at least twice and ρ an at least once differentiable function. Note, that after transformation into its weak form and shifting of the divergence operator onto the test function by application of partial integration, these conditions are relaxed and after discretization piecewise linear basis function (P 1 ) can be used for the scalar potential u and the scalar indicator function ρ only has to be constant within each element. The numerical optimization is performed by casting the design goal into an objective functional Ĵ (ρ) = J(H(ρ), ρ) and minimizing it using an iterative, gradient based optimization algorithm.

Adjoint approach
In order to efficiently minimize the reduced objective functional Ĵ (ρ) containing the design goal, its gradient with respect to the design variable ρ , ∇ ρĴ is needed. This is where the adjoint approach comes into play, offering a very efficient method to calculate the gradient. As derived in detail in the appendix ("Appendix" section) within a "first optimize then discretize" scheme, the gradient is given by where is the adjoint variable that is obtained by solving the adjoint equation with boundary condition ρ(x) ∈ [0, 1], where ρ(x) = 0: no material 1: material

Numerical experiments
In the following, the topology of thin film magnetic flux guide concentrators is optimized towards maximization of the vertical (y-)component of the magnetic field H y inside the target domain h . The objective functional that is to be minimized reads Ĵ = − h H y dV . The geometry and the mesh used is depicted in Fig. 1 where also the the optimization domain ( opt ) and the domain in which the magnetic field is evaluated ( h ) is specified. Note that m = opt ∪ h . Note, that we use an irregular mesh in order to avoid any solution bias introduced by a given regularity. The presented numerical experiments include the usage of hard as well as soft magnetic material.
The potential u and the adjoint variable are calculated using piecewise linear basis functions (P 1 ) and the scalar indicator function ρ as well as χ , M r , H ext and the derived strayfield −∇u and ∇ are constant within each element.
The optimizations are performed using the Scipy implementation of the limited-memory Broyden-Fletcher-Goldfarb-Shanno minimization algorithm with bounds (L-BFGS-B) 13 . Note, that as shown in 14 , common readyto-use optimization algorithms such as this, are constructed to only receive coefficient vectors as input and internally use the euclidean ( ℓ 2 ) inner product to perform all necessary operations (e.g. approximation of the Hessian). Since here an optimization in a Hilbert spaces with the L 2 inner product is performed, this leads to mesh dependent convergence rates. To prevent this, there exists the python package Moola 15 that internally uses the correct inner product during optimization. However, in order to perform topology optimization bounds on the optimization variable ρ have to be enforced since 0 ≤ ρ ≤ 1 and so far, there is no optimization algorithm offered by Moola enabling bounds on the optimization variable. Furthermore, note that the L 2 Riesz representer of the gradient is independent of the flux concentrator's geometrical dimensions. However, the objective functional is not since it includes and integration over the volume h . For the chosen geometry dimensions the gradient size is relatively large compared to the size of the objective functional. This leads to the L-BFGS-B algorithm terminating prematurely with problems during the line search, simply returning the initial distribution of ρ . This is since the Armijo rule, one of the Wolfe conditions necessary to end the line search is never met. The objective functional value exposed to the optimization algorithm, is therefore scaled to avoid this problem. This scaling does not alter the results since internally, the L-BFGS-B algorithm performs a quasi-Newton method on the gradient and the function value is only used during the line search.
For results of optimizations that terminate with dofs of ρ at intermediate values 0 ≤ ρ ≤ 1 a brute force regularization is applied where all dofs with ρ < 0.5 are set to 0 and all dofs with ρ ≥ 0.5 are set to 1.
Permanent magnet optimization. First, the topology for an ideal permanent magnet with susceptibility χ = 0 , remanence magnetization M r = 280000 A/m and the anisotropy axis chosen as the y-axis is optimized towards maximization of the vertical (y-)component of the magnetic field H y inside the domain h , as introduced above.
As presented in 8 for a point like target domain, the optimal topology can be derived analytically by considering the magnetization inside opt as being comprised of single dipoles µ . A dipole at a vertical (y-)distance produces a positive flux B y (x) at a point x for an horizontal (x-)distance of The surface condition x = √ 2y therefore shows that the analytical solution is a cone with an opening angle of α ≈ 35.26 • . Note that the assumption of a point like target domain is still valid, since the changes due to the finite For the ideal permanent magnet, the optimized topology, as shown in Fig. 2, produces a mean field H y = � h H y dV /V � h of 22.83 mT. This solution is found independently of the initial distribution of the density functional ρ init for a penalty parameter of p ≤ 4 . Furthermore, scaling the gradient has no influence on the obtained solution, but increasing the gradient size reduces the number of function evaluations n fev needed. This is since the size of the gradient determines how fast the dofs of ρ are set to 0 or 1. This shows, that the optimization problem in the case of zero susceptibility is linear in ρ as expected since for χ = 0 the only source term in the forward problem (7) is given by ρ p M r where M r is constant. In Fig. 3 the value of the smallest dof in the gradient and n fev are plotted against the gradient scaling factor γ for p = 1 . As soon as the smallest dof is larger 1 for γ = 10 3 , all elements are set to 0 or 1 during the optimization, leaving no dofs with intermediate values of ρ and n fev reaches a minimum of 2.
Realistic hard magnetic material χ = 0.2. Using a realistic hard magnetic material with a susceptibility of χ = 0.2 , the optimal solution found has a mean field of H y of 22.35 mT. This smaller mean field value compared to the ideal permanent magnet is due the effect of demagnetization. The optimal topology found ( Fig. 4) differs only slightly from the topology obtained from the ideal permanent magnet optimization.
The optimal solution is found independently of the initial distribution of the density functional ρ for p = 1 . For larger p-values the solution gets worse showing, that the optimization is no longer independent of the parameter p. However, for p ≤ 4 the difference in the mean field is below 0.5%. Scaling the gradient again reduces the number of necessary function evaluations to perform the optimization n fev . However, since the optimization problem is no longer linear in ρ , the number of function evaluations n fev needed is larger than for   www.nature.com/scientificreports/ the ideal permanent magnet case. In Fig. 5 n fev in dependence of the gradient scaling factor γ for p = 1 is shown. The optimization is fastest for a constant initial distribution of ρ init = 1 where for γ ≥ 10 n fev < 250 . Also, the optimization process terminates before all dofs of ρ are set to 0 or 1, although the number of dofs with an intermediate value of ρ is negligible for all simulations with p = 1 . It has been verified that no better solution is found by scaling the gradient down to γ = 10 −5 . Note, that the topology obtained from the optimization using an ideal permanent magnetic material, for a realistic hard magnetic material with χ = 0.2 produces a magnetic field that is only negligibly smaller (<0.1‰) than the optimal topology obtained using the realistic hard magnetic material during optimization. This confirms, that χ = 0 is a good assumption for a permanent magnetic structure.
Optimization of a soft magnetic structure.
Next, the topology of a soft magnetic thin film is optimized to act as a magnetic flux guide concentrator maximizing the magnetic flux in vertical (y-)direction inside h . Such a structure can be used to e.g. improve the magnetic field detection limits of spin valve sensors 16 by ehancing the external field at the sensor position. The optimization is performed for three different values of χ of 10 2 , 10 3 and 10 4 with an external field of 10 mT applied. The optimized topologies for the three different χ-values are shown in Fig. 6. The found structures show the same "tanga" shape as reported in 16 as optimal design. While the optimal topology differs only slightly for χ = 10 3 and χ = 10 4 , for χ = 10 2 the topology is significantly broader. In fact, switching the optimal topologies of χ = 10 3 and χ = 10 4 changes the generated mean field by under 0.2% . The generated mean field and the increase with respect to the full plate ( ρ = 1 in opt ) is shown in Table 1.
As opposed to the optimization of hard magnetic materials, the optimization using a soft magnetic material depends much more on the initial distribution of the density functional ρ init as well as on the p-parameter. In Fig. 7 the p-dependence of the mean field H y for the regularized solutions for χ = 10 3 is plotted for different initial distributions of ρ init and for different values of the gradient scaling factor γ.
The solution found for a constant ρ init = 1 and a randomly initialized ρ init are clearly inferior to the solutions found with a constant ρ init = 0.5 . For ρ init = 0.5 the mean field increases with increasing p until it plateaus. While decreasing γ does not lead to a higher mean field, the p-value at which the mean field plateaus decreases with decreasing γ . At γ = 1 the plateau starts at p = 26 while for γ = 10 −4 the plateau value is 14. Further decreasing γ does not alter the plateau value anymore. In fact, as shown in Fig. 8 also for single p-values, the mean field increases with decreasing γ until a maximum is found for 10 −4 . Note, that for higher p-values than the plateau value, the mean field differs well under 1% from the mean field obtained by the optimization performed with the plateau value, except for individual outliers with very small mean field values for large p-values.
In Fig. 9, the amount of dofs with intermediate values of ρ after optimization is plotted in dependence of p. While the solutions with p ≤ 4 show no dofs with intermediate values of ρ , with increasing p, their number increases and reaches over 60% for the optimal solutions. This is in contrast to the original purpose of the penalization factor p. However, since the changes of the mean field H y due to regularization is under 1% for all optimizations, the effect of the regularization is still negligible. This is since the remaining intermediate values  www.nature.com/scientificreports/ are very close to 0 and 1 respectively. In Fig. 6 the optimized topologies are shown before and after regularization for comparison. The dependence of the number of function evaluations necessary to perform the optimizations n fev on p is plotted in Fig. 9 for two different values of the gradient scaling factor. It can be seen n fev is lower with a smaller gradient scaling factor. With a gradient scaling factor of 10 −4 , the optimization with p = 14 needs 325 n fev while with a gradient scaling factor of 1 for the optimization with p = 26 n fev = 14631 . In Fig. 10     www.nature.com/scientificreports/ For χ = 10 4 the p-dependence of the mean field plateaus in the same way as for χ = 10 3 , although with a different plateau value of p = 17 for a gradient scaling factor of 10 −4 . The dependency on the gradient scaling factor is also similar as shown in Fig. 8. For χ = 100 the situation is slightly different. Here, for gradient scaling factors smaller 10 −4 the mean field increases with increasing p-value until a maximum is found at p = 5 and the mean field decreases when further increasing p.

Conclusion
An algorithm to optimize the topology of magnetic structures based on a hybrid FEM-BEM method and the density approach for topology optimization has been presented. The gradients necessary to efficiently solve the optimization are obtained very efficiently using the adjoint approach. It has been shown that within a "first optimize then discretize" scheme the continuous forward operator is self-adjoint and therefore, the same method can be used to solve the forward as well as the adjoint problem.
The dependence on the p-parameter, originally introduced to penalize intermediate values of the density indicator function 0 ≤ ρ ≤ 1 has been investigated for the presented numerical experiments. For ideal permanent magnetic structures ( χ = 0 ) the linearity of the optimization problem can be utilized by enlarging the gradient. Scaling the gradient by a factor > 10 3 , with p = 1 the optimization needs only two function evaluations and the found solution is in accordance with analytical expectations and shows no intermediate values of the indicator function ρ . For a realistic permanent magnetic material with χ = 0.2 the optimized topology produces a magnetic field that differs by < 0.1 ‰ from the topology found for the ideal permanent magnetic material. This shows that for thin films with negligible demagnetization, χ = 0 is good approximation for a realistic permanent magnetic material with χ = 0.2 . Furthermore, the number of function evaluation needed using a realistic hard magnetic material is in the hundreds and hence much larger than the number of evaluations needed for the assumption of an ideal permanent magnetic material.
For the optimization of soft magnetic materials the optimal solutions are found for a gradient scaling factor ≤ 10 −4 . The solutions found get better with increasing value of the p-parameter. For χ = 10 2 the best solution is found for p = 5 and gets worse for larger p-values. For χ = 10 3 and χ = 10 4 the solutions improve with p until it plateaus for p ≥ 14 and p ≥ 17 respectively. Note, that in contrast to its original purpose, the p-parameter does not decrease the number of dofs with intermediate values of ρ . On the contrary, the percentage of dofs with intermediate values of ρ increases with increasing p. However, using a brute force regularization forcing all dofs with ρ < 0.5 to 0 and all dofs with ρ ≥ 0.5 to 1 only negligibly impacts the field value created by the structures. The optimized topologies have a shape similar to the optimal shape found experimentally in 16 .

Appendix
Deriving the gradient ∇ ρĴ using the adjoint approach. In the following the adjoint approach is applied to obtain the gradient ∇ ρĴ of the reduced objective functional Ĵ (ρ) with respect to the design variable ρ . A "first optimize then discretize" scheme 17 is used, within which the gradient is derived using continuous operators. These operators are, as will turn out, self-adjoint where the discrete forward operator of the hybrid FEM-BEM algorithm is not. This simplifies the numerics involved in solving the adjoint equation. Following 17 , in order to derive the gradient ∇ ρĴ we consider the derivative Here, du dρ is the numerically challenging term since u depends on ρ via the forward problem (7). In the discretized form, naive evaluation of du dρ by means of finite difference would require n ρ solutions of the forward problem, where n ρ is the number of dofs of ρ . Looking at the derivative of the forward problem the critical term can be expressed as Inserting this into Eq. (12) gives an expression for the derivative where du dρ does not appear explicitly anymore. Consider now that the goal is to perform an optimization on ρ , an at least once differentiable element of a suitable Hilbert space R. This means, that an initial guess ρ 0 is updated iteratively using a suitable update candidate to obtain a new ρ i that further reduces the objective functional. Since the derivative dĴ dρ : R → R is a functional and therefore an element of the dual space R * of R, it is not a suitable candidate to be used as update during optimization. What is needed is an element of the primal space similar to ρ , the gradient ∇ ρĴ ∈ R 14 .
Since within a finite element setting all function spaces are Hilbert spaces, the Riesz representation theorem can be used to obtain the gradient. The Riesz representation theorem states, that the dual space H * of a Hilbert space H is isometric to itself. This means that the identification H = H * can be done and for a functional I ∈ H * applied to u ∈ H a unique element w of the primal space H can be found such that where (·, ·) H is the inner product defined on H and w ∈ H is the Riesz representer of I with respect to the inner product on H, also denoted as w = R H (I) . Note that, for the presented example the appropriate inner product is the L 2 inner product that arises naturally when the forward problem is brought into its weak form by multiplication with a suitable test function. The specific reference is therefore dropped from this point on an each inner product as well as the Riesz representers are assumed to be with respect to L 2 if not stated otherwise.
The gradient ∇ ρĴ ∈ R can now be identified as the  where partial integration is applied twice and the boundary integrals vanish if the same boundary condition as valid for the scalar potential u is also also demanded from w and z 18 .
The adjoint equation becomes with boundary condition Note that this equation has the same form as the forward problem F itself. This is a consequence of the selfadjointness of the operator ∂F ∂u that is given since the above calculation shows that ∂F ∂u = ∂F ∂u * . Therefore, also has to be twice differentiable and the adjoint equation can be solved by the same method as described in "Forward problem" section, as will be seen below in "Solving the adjoint equation" section. The operator ∂F ∂ρ * is obtained by again first deriving the operator ∂F ∂ρ using the Fréchet derivative Using Eq. (18) gives Here again partial integration was applied and the boundary integral vanished for s satisfying the boundary condition (26), which is fulfilled since in the expression for the gradient ∂F ∂ρ * is applied to that also satisfies this boundary condition.
Finally, the gradient can be written as Solving the adjoint equation. In order to solve the adjoint Eq. (25) we first have a closer look at it. On the right hand side is the partial derivative of the objective functional J with respect to the scalar potential u. However, we will not directly impose conditions on the potential u but on the magnetic field H that is defined as H = −∇u . Consider the weak form of the adjoint equation   where the boundary term vanishes due to the boundary condition on . Assuming that the objective functional J is defined locally on an area opt , this is the weak form of that is identical to the forward problem in Eq. (7) with instead of u, R ∂J ∂H replacing ρ p M r and H ext = 0 . We can therefore use the same algorithm to solve the adjoint equation as is used to solve the forward problem. Note, that it is perfectly fine to assume the same jump conditions on as on u since they are contained in the weak form of the forward problem anyways as shown in 19 .
Obtaining the Riesz representers. To obtain the Riesz representers R ∂J ∂H and R ∂J ∂ρ appearing inside the adjoint equation and the expression for the gradient, the coefficient vectors J H and J ρ belonging to the functionals ∂J ∂H and ∂J ∂ρ respecively, are calculated using FEniCS 11 . These coefficient vectors can be used to calculate the action of the functionals onto an element of the primal space r 14 . For the functional ∂J ∂ρ this reads where (·, ·) ℓ 2 represents the euclidiean inner product. On the other hand, using the Riesz representer R ∂J ∂ρ the action can be written as with ψ i (x) being the discrete basis functions and M = � ψ i ψ j dV being the mass matrix. Since the two equations must deliver the same result and furthermore M T = M , the Riesz representer with respect to the L 2 inner product can be obtained from the coefficient vector as (30)