Solving Large-Scale Inverse Magnetostatic Problems using the Adjoint Method

An efficient algorithm for the reconstruction of the magnetization state within magnetic components is presented. The occurring inverse magnetostatic problem is solved by means of an adjoint approach, based on the Fredkin-Koehler method for the solution of the forward problem. Due to the use of hybrid FEM-BEM coupling combined with matrix compression techniques the resulting algorithm is well suited for large-scale problems. Furthermore the reconstruction of the magnetization state within a permanent magnet as well as an optimal design application are demonstrated.

where the magnetic scalar potential u is given by where m is the magnetization and Ω is the magnetic region.
The forward problem requires the solution of the potential u on the region Ω h generated by the magnetization in a magnetic region Ω m (see Fig. 1). This problem is solved by considering a single region Ω = Ω m ∪ Ω h with m(x) = 0 for x ∈ Ω h . The hybrid FEM/BEM method introduced by Fredkin and Koehler 3 is one of the most accurate methods for the solution of this problem and will be used in the following. Consider the following splitting of the solution u: 1 2 Here u 1 is defined by 1 ∂ ∂ = − ⋅ ∂ Ω n m n u on (8) 1  Ω = . u 0 in \ (9) 1 3 This Neumann problem is solved with the finite-element method. While u 1 solves for the right-hand side m and fulfills the jump condition of the normal derivative − m · n, it is not continuous across ∂ Ω as required. This jump is compensated by u 2 which is defined as

Figure 1. Discretized magnetization region Ω m (blue) and measurement region Ω h (brown).
Since the strayfield problem is scale-invariant, length units are omitted. The magnetization is defined on a unit cube. Measurement boxes of thickness 0.04 are located next to each side of the cube, using an airgap of 0.1 (the 2 boxes in front are not shown in the figure).
Scientific RepoRts | 7:40816 | DOI: 10.1038/srep40816 This system is solved by the following double-layer potential For efficiency reasons the double-layer potential is only computed on the boundary ∂ Ω using a Galerkin boundary-element method. Subsequently these values are used as Dirichlet boundary conditions to solve u 2 within Ω using the finite-element method.
All potentials are calculated using piecewise linear basis function ( 1 ) and the derived strayfield would be constant within each element. Thus, a mass lumping procedure needs to be used to project the field onto piecewise linear basis functions which are defined on each vertex of the mesh.
Inverse Problem. The inverse problem can be understood as a PDE constrained optimization problem. Due to the ill-posedness of the inverse problem, some additional information needs to be provided to allow reasonable results. This can be achieved by using Tikhonov regularization, which uses an additional penalty term which should be considered for the optimization. A possible candidate for the objective function is where h m is the prescribed (measured) field in Ω h and α is the Tikhonov constant corresponding to the regularization functional T(m). This functional should be minimized, constrained by with boundary conditions as given above. We aim to solve this problem using a gradient based iterative minimizer. The constraint F gives an implicit expression for u(m) which allows to directly calculate the desired gradient by The inefficient direct calculation of the term ∂ ∂m u can be avoided using the adjoint approach, which makes use of the derivative of the constraint equation to eliminate the problematic term. Compared with a naive implementation using a dense system matrix ∂ ∂ u m , the computational as well as storage costs can be reduced where λ is given by the adjoint equation Since the Poisson problem is self-adjoint, the adjoint system (19) can be solved along the lines of the forward problem. Computing the variational derivative on the RHS yields where the sources (RHS) live on Ω h and the solution is only computed on Ω m . The same boundary conditions as for the forward problem hold. Thus, the above described hybrid method is applied. The gradient of J is then finally given by Note that ∇u and λ ∇ are projected onto 1  before computing (20) and (21) respectively. The algorithm is implemented using Magnum.Fe 10 , which is based on the finite element library FEniCS 11 . This allows a very comfortable definition of the regularization functional. Furthermore, automatic differentiation can be used for the calculation of the corresponding partial derivatives. The algorithm was verified by comparison with a FEM-only implementation, using the dolfin-adjoint library 12 , which allows to automatically derive the adjoint equation for a given forward problem.

Numerical Experiments
The presented algorithm is validated by the reconstruction of the flower-state within a magnetic unit cube. The strayfield is calculated within measurement planes next to each side of the cube (see Fig. 1). The magnetic state of the cube is parametrized by For this specific flower-state the absolute value of the magnetization is known to be constant. Thus, the solution of the inverse problem could be simplified by using the following penaltization functional The assumption of a constant magnetization may be a good approximation for (isotropic) permanent magnetic materials. Due to the large magnetic remanence and the relatively small susceptibility the induced magnetization may be negligible (see ref. 13 for a simple model of isotropic permanent magnetic materials).
A Gaussian noise with zero mean and a standard deviation σ = 10 −4 has been added to the field, calculated by the forward problem, which should simulate unavoidable measurement errors. The minimization problem is solved by a gradient descent method combined with a line-search strategy. As expected, reconstruction without using a proper regularization leads to large magnetization vectors near to the edges of the unit cube. Increasing the regularization parameter α, first avoids the over-fitting of the noisy measurement data, but finally leads to blurring of the reconstruction results if α gets too large. Determining the optimal alpha is a crucial step for the solution of an inverse problem. The results for the reconstruction of a flower state with c tilt = 2 using α = 10 −3 is visualized in Fig. 2. Although there are some deviations of the reconstructed magnetization from the reference state. A cut through the y = 0 plane is visualized. Starting from the initial parametrized flower state in Ω m the magnetic strayfield is calculated within the fieldboxes Ω h (green arrows). In order to simulate measurement errors a Gaussian noise with σ = 10 −4 has been added to the forward strayfield. The reconstructed magnetization as well as strayfield are computed using α = 10 −3 (red arrows). The relative differences of initial and reconstructed states are indicated by the gray-scale. Maximal relative errors of the x-components amount to 0.25 for the magnetization, and 5 · 10 −3 for the induced magnetic field, respectively.
It can be seen that the created strayfield is nearly identical. As stated above this is a clear indication of the bad condition of the inverse problem.
Using the so-called L-curve method 14 , allows to visualize the trade-off between the reconstruction of the strayfield and the fulfillment of the regularization constraint. Plotting the regularization norm (also called solution norm) T m ( ) over the residual norm ∇ − − h u m for different regularization parameters α, shows an L-shaped curve. The optimal α opt can be selected at the corner of the L-curve which means that α is large enough to reduce the regularization norm significantly, but it does not change the residual norm too much. The resulting L-curves for the reconstruction of the flower state for different noise levels are summarized in Fig. 3. An optimal value α opt ≈ 5 · 10 −5 can be found.
The application of the presented method to optimal design problems should now be demonstrated by means of a Halbach cylinder configuration. The goal is to find a magnetization configuration within a cylindrical domain, which creates a homogeneous strayfield inside of the cylinder. The magnetization domain Ω m has an  Starting from a homogeneous strayfield the presented algorithm reproduces a Halbach like magnetization configuration within Ω m (red arrows). A constant-norm regularization with α = 10 4 is used and shows a nearly perfect match with the analytical solution (green arrows). The resulting strayfield is calculated inside Ω h and shows a nearly homogeneous distribution (red arrows). The relative errors of the magnetization magnitude and the reconstruced strayfield are indicated by the gray-scale, and their maximum amount to 2% and 6%, respectively. outer radius r o = 1.0, an inner radius r i = 0.6, and a height of h = 2.0, while a cylindrical measurement domain Ω h with radius r m = 0.5 and the same height is used. The magnetization vectors should have constant norm, which suggests using the constant-norm penaltization functional (24). The analytic solution for a cylindical Halbach array 15 can be expressed in cylindrical coordinates as where ρ, φ are the cylindrical coordinates, with the corresponding unit vectors e ρ , e φ . As visualized in Fig. 4 there is a nearly perfect match of the reconstructed and the analytical Halbach configuration.

Conclusion
An efficient algorithm for the solution of inverse problems has been introduced. The use of the Finite Element library FEniCS allows to easily define application specific regularization functionals in a very flexible way. Thus, the implemented algorithm is suitable for a wide range of applications including reverse engineering of magnetic components, design and optimization of magnetic circuits and topology optimization, respectively. Using the hybrid FEM-BEM method proposed by Fredkin-Koehler allows to handle the open-boundary problem accurately and without the need for global mesh including a large airbox. Source identification has been validated by the successful reconstruction of the magnetic flower-state within a unit cube by means of Tikhonov regularization. The selection of a suitable regularization parameter has been demonstrated using the L-curve method. Finally, the application of the method to an optimal design problem has been demonstrated by means of an Halbach cylinder, which is nearly perfectly reproduced.