Data-Driven Solvers for Strongly Nonlinear Material Response
DData-Driven Solvers for Strongly Nonlinear Material Response
Armin Galetzka , Dimitrios Loukrezis , , Herbert De Gersem , , Technische Universit¨at Darmstadt, Institute for Accelerator Science and Electromagnetic Fields (TEMF)Schlossgartenstrasse 8, 64289 Darmstadt, Germany Technische Universit¨at Darmstadt, Centre for Computational EngineeringDolivostrasse 15, 64293 Darmstadt, Germany Katholieke Universiteit Leuven - Kulak, Wave Propagation and Signal Processing Research GroupEtienne Sabbelaan 53, 8500 Kortrijk, Belgium
Abstract
This work presents a data-driven magnetostatic finite-element solver that is specifically well-suited to copewith strongly nonlinear material responses. The data-driven computing framework is essentially a multiobjectiveoptimization procedure matching the material operation points as closely as possible to given material data whileobeying Maxwell’s equations. Here, the framework is extended with heterogeneous (local) weighting factors - oneper finite element - equilibrating the goal function locally according to the material behavior. This modificationallows the data-driven solver to capture sharp gradients in the constitutive law of strongly nonlinear materials,which constitute problematic cases for standard data-driven solvers with a homogeneous (global) weighting factor,hindering their efficiency and accuracy. The local weighting factors are embedded in the distance-minimizing data-driven algorithm used for noiseless data, likewise for the maximum entropy data-driven algorithm used for noisydata. Numerical experiments based on a quadrupole magnet model with a soft magnetic material show that theproposed modification results in major improvements in terms of solution accuracy and solver efficiency. Forthe case of noiseless data, local weighting factors improve the convergence of the data-driven solver by orders ofmagnitude. When noisy data are considered, the convergence rate of the data-driven solver is doubled. keywords— data-driven computing, data science, electromagnetic field simulation, noisy measurements, non-linear material response, soft magnetic materials
Numerical simulations play an important role in the study of electromagnetic (EM) phenomena and are irreplaceablewhen designing EM devices. Irrespective of the numerical approximation method of choice, e.g. the finite elementmethod (FEM) [20, 29], the boundary element method (BEM) [6], or the finite integration technique (FIT) [39, 40],to name but a few, the underlying initial/boundary value problems (IBVPs) are governed by Maxwell’s equations,which relate EM fields to their sources, and by constitutive laws, which describe material responses. With respectto the latter, in the context of this work, we focus on the relation between the magnetic field strength ~H and themagnetic flux density ~B , which depend on the magnetic permeability tensor µ , respectively, on the reluctivity tensor ν = 1 / µ . With the exception of the often irreducible aleatory uncertainties [37], Maxwell’s equations are acceptedas exact. This is not true for the constitutive laws, which are in nearly all cases too complicated to be exactlydescribed. For some simple materials such as vacuum, the constitutive relation is regarded as exactly known andgiven by ~B = µ ~H , µ being a constant factor or tensor. However, for most magnetic materials the relation between ~H and ~B is strongly nonlinear and not exactly known. In such cases, a constitutive law typically takes the form ofan empirical model, which is most commonly based on fitting an analytical expression to experimental observationdata.The dominant paradigm in numerical EM computations is thus based on material models derived after availabledata. These models are subsequently provided to numerical EM field solvers [5, 16]. For that reason, significanteffort is placed on improving the data-fitting and regression techniques that are used for material modeling purposes.For example, advanced machine learning methods are used in an effort to capture complicated constitutive lawsmore accurately [38]. The modeling procedure can be further enhanced with constraints that the material model1 a r X i v : . [ phy s i c s . c o m p - ph ] A ug ust conform to, based on physical considerations or domain expertise, [15, 31]. If large data sets are available,other approaches advocate to deduce material behavior purely based on the data [32, 33]. Nevertheless, regardlessof implementation specifics and level of sophistication, data-based constitutive laws are uncertain in terms of theactual model describing them, a type of uncertainty called epistemic [37]. This model-form uncertainty unavoidablyaffects the solution of the numerical solver, thus rendering its predictions possibly questionable. Moreover, thequality of these predictions cannot be improved with a finer spatial or temporal resolution, as the material model isa “hard-coded” element of the solver.An alternative and fundamentally different paradigm regarding the incorporation of material data into numericalcomputations, bearing the name “data-driven computing”, emerged recently from the field of computational mech-anics. This approach was introduced in the seminal work of Kirchdoerfer and Ortiz [23] and then further developedand analyzed for a range of problems arising in computational mechanics [7, 8, 13, 18, 25, 27, 36]. For the case ofEM simulations, a comparable formulation based on the variational principle was already proposed by Rikabi et al.[35] in 1988. The main idea behind the data-driven computing paradigm is the reformulation of the IBVP describingthe physical phenomenon under investigation, such that the material modeling step is bypassed altogether, thuseliminating the attached epistemic uncertainty. In this framework, the corresponding data-driven numerical solveroperates directly on raw material data, essentially being material-model-free. In particular, the data-driven solverseeks to minimize the distance between the field states in phase space fulfilling exact laws and the field states inphase space representing the constitutive law, where the latter is described in the form of a measurement data set.In the present work, the term “field state” refers to a ~B - ~H pair.By relying on raw material data, the distance-minimization scheme is naturally sensitive to outliers in the dataset. Such outliers appear, for example, when we consider the case of noisy material data [2]. The possibly dominantinfluence of the outliers has a significant impact on the robustness of data-driven solutions. To address this problem,the distance-minimization scheme can be modified, such that the data are first organized in clusters based on themaximum entropy principle and then the solver minimizes the free energy in phase space [24]. The data-drivencomputing paradigm has been further extended for the purpose of material response identification [27], in whichcase the data-driven solver is effectively “inverted” as to locate states that sample the mechanical response of elasticmaterials. A data-driven solver integrating the material response identification approach of [27] has been put forth in[36], showing significant accuracy and efficiency improvements. Outside the field of computational mechanics, data-driven magnetostatic FEM solvers were developed by the authors and co-workers in [10], for the case of noise-freedata and heterogeneous domains where exact and data-based material models co-exist.The data-driven solvers proposed in the literature so far employ a global weighting factor, thus treating all finiteelements in the computational domain equally, regardless their material filling. However, a global, homogeneousfactor may actually hinder the efficiency and accuracy of the data-driven solver in the case of strongly nonlinearmaterials with constitutive laws featuring sharp gradients, due to their innate inability to capture larger differencesin the material’s response. The present work extends the data-driven computing paradigm in order to accommodatesuch strongly nonlinear material behavior as well. Our contribution is motivated by the response of soft magneticmaterials, the BH -curves of which typically consist of a very steep linear part, followed by the saturation part aftera sharp transition. To address the problem, we modify the data-driven solver’s algorithm as to assign local weightingfactors, i.e. one per finite element. Since no a priori knowledge about the working points along the BH -curve withrespect to each finite element is available, the assignment of local weighting factors is performed adaptively during thedata-driven solver’s iterations. We show that the proposed modification to the data-driven solver results in significantcomputational gains for noiseless and noisy material data alike. In the noiseless case, the data-driven solver utilizinglocal weighting factors converges orders of magnitude faster than the data-driven solver based on global factors. Forthe case of noisy data, the convergence rate is improved by a factor greater than 2. Since data-driven methods basedon both global and local weighting factors are applied and compared against each other, a side contribution of thiswork is the extension of [10] to the case of noisy material data. Moreover, with the introduction of local weightingfactors, the two intrusive data-driven approaches suggested in [10] are rendered unnecessary.The rest of this paper is structured as follows. In Section 2 the main idea behind the data-driven paradigm isdescribed. Section 2.1 presents an illustrative example showing the idea behind data-driven computing. Afterwards,in Section 2.2 we specialize the data-driven formulation to the case of magnetostatics, which is then in Section 2.3converted into the weak formulation and finally discretized in Section 2.4. Furthermore, in Section 2.5 we recall thecase of problems with heterogeneous domains, where both exactly known materials and data-based materials coincide.Lastly, Section 2.6 concludes with the introduction of the 2D data-driven weak formulation. Section 3 is dedicatedto the management of noiseless data. In Section 3.1 we recall the usage of a global weighting factor over the entire2 N I S Fe ‘ Fe R Fe R air F = N I ⇐⇒ (a) HB D b DM HB D b DM M ∩ ˆ D (b) Figure 1: (a) C-yoke and equivalent magnetic circuit. (b) The black line shows the constrained set M of the statesthat fulfill the circuit law. The classical solution is given as the intersection of M and the set ˆ D which correspondsto a regression-based BH -curve (shown in red). The data-driven solution is given by the state in the measurementdata set closest to M (data given as blue dots, solution marked with a red circle).domain, whereas in Section 3.2 we present the novel approach with local weighting factors. In Section 3.3 numericalexperiments are carried out for both approaches on a 2D quadrupole magnet and show considerable improvementswhen local weighting factors are employed. The handling of noisy measurement data is described in Section 4. First,we recall the basic idea of noisy measurement treatment introduced by Kirchdoerfer and Ortiz [24] in Section 4.1.Afterwards, we embed the local weighting factors into the framework of noisy measurement data. Again, numericalexperiments show the enhancement brought by local weighting factors in Section 4.3. Concluding remarks are givenin Section 5. Before deriving the general data-driven magnetostatic formulation, we first present an illustrative example of thedata-driven computing paradigm within the context of magnetic field computations. In particular, we consider theC-shaped iron yoke depicted in Figure 1a, where the magnetic flux Φ is generated by the current I flowing througha coil with N turns. Assuming no fringe fields on the outside of the iron core and a magnetic path with constantcross-section S Fe , the field problem can be approximated with the equivalent magnetic circuit shown in Figure 1a.Furthermore, on circuit level, it holds that B = (cid:12)(cid:12)(cid:12) ~B (cid:12)(cid:12)(cid:12) and H = (cid:12)(cid:12)(cid:12) ~H (cid:12)(cid:12)(cid:12) , thereby assuming the magnetic fields ( H air , B )and ( H Fe , B ) which are homogeneous in the air and iron regions respectively. The magnetomotive force F is thengiven as F = I ~H · d ~s = N I = F Fe + F air , (1)where F Fe , F air refer to the magnetomotive forces in the iron yoke and in the air gap, respectively. Under the sameassumptions, the magnetic field strength is constant, which allows us to express the magnetic flux as Φ = BS Fe .Hence, the magnetic flux can be expressed asΦ = FR Fe + R air = F air R air = F Fe R Fe , (2)where R Fe and R air refer to the magnetic reluctance of the yoke and the air gap, respectively. Combining (1) and(2) results in BS Fe = − H Fe ‘ Fe R air + FR air , (3)where ‘ Fe refers to the average length in the iron yoke, see Figure 1a and H Fe to the - currently unknown - magneticfield strength in the yoke. 3he goal is to find the state, i.e. the pair ( H Fe , B ), that fulfills the circuit law (3). All conforming states arecollected in the set M , such that M = (cid:26) ζ : ζ = ( H Fe , B ) ∈ Z : BS Fe = − H Fe ‘ Fe R air + FR air (cid:27) , (4)where the set Z refers to all possible ( H, B ) states. All states that conform to the circuit law (3) form a line in the( H Fe , B ) phase space (black line in Figure 1b). In order to find a solution to (3), additional information regardingthe relation between H Fe and B is necessary. This information is usually given in the form of measurement data, i.e.as discrete points ( H Fe ,m , B m ) , m = 1 , . . . , M, collected in a set D (blue points in Figure 1b). The traditional wayof solving this problem involves fitting a continuous constitutive law to the given data, e.g. by means of regression(red line in Figure 1b). The corresponding set of phase space states is then given byˆ D = { ζ : ζ = ( H Fe , B ) ∈ Z : B = f ( H Fe ) } , (5)where f ( H Fe ) : R +0 → R +0 defines the nonlinear BH -curve given by the red graph in Figure 1b. The solution is thengiven as the intersection between the states fulfilling (4) and the states compatible with the constitutive law, i.e.( H Fe , B ) = M ∩ ˆ D . Diversely, in the data-driven paradigm, a solution in the intersection M ∩ D is sought.Note that the intersection between M and D is very likely to be an empty set, due to the fact that only a finitenumber of measurement data points are available. Therefore, instead of searching for a state in D , we accept assolution the state that fulfills the governing equation (3), while at the same time being the closest to a state in D .Accordingly, we search for the state that minimizes the distance between the governing equations and the availabledata set, given by ζ opt = argmin ζ ∈M F ( ζ, D ) , (6)where F ( ζ , ζ ) is a distance function between states ζ , ζ in phase space. This procedure is the core idea behind adata-driven solver, and is upscaled to a field formulation in the subsequent section. Next, we present the data-driven formulation for the magnetostatic case. The governing equations are Amp`ere’s lawand Gauss’s law for magnetism, which respectively readcurl ~H = ~J, in Ω , (7a)div ~B = 0 , in Ω , (7b)where Ω denotes the computational domain. These equations are complemented with appropriate boundary condi-tions (BCs) on the boundary of Ω, denoted by Γ, thus leading to the ubiquitous setting of a boundary value problem.We also assume that the two equations are free of any aleatory uncertainty, e.g. with respect to the geometry of thedomain or the source terms. Aleatory uncertainty is of course possible and often inevitable, however, this case willnot be considered here.To solve the boundary value problem, a constitutive relation between ~B and ~H is necessary. In the traditionalframework, a material model is derived after a measurement data set e D containing ( ~H, ~B )-pairs for different operatingpoints, such that e D = n(cid:16) ~H ? , ~B ? (cid:17) , (cid:16) ~H ? , ~B ? (cid:17) , . . . , (cid:16) ~H ?M , ~B ?M (cid:17)o . (8)In the context of data-driven solvers, the material modeling step is omitted. A formulation for both phase spacesis needed instead, i.e. regarding the phase space that covers the governing (exact) equations and the phase spaceaddressing the material relation in view of the measurement set e D . Similar to Section 2.1, the states that conformto the magnetostatic equations (7) can be found in the space M = n ζ ( ~x ) : ζ = (cid:16) ~H ( ~x ) , ~B ( ~x ) (cid:17) ∈ H(curl) × H(div) , ∇ × ~H = ~J, ∇ · ~B = 0 and BC , ~x ∈ Ω o , (9)where H(curl) denotes the space of square-integrable functions with square-integrable curl and H(div) the space ofsquare-integrable functions with square-integrable divergence. In order to derive the minimization problem as in (6),4he local measurement set (8) has to be reformulated such that the measurement set is valid in the entire domainof the underlying material, i.e. a spatial representation is necessary. Thus, we form a global set of discrete materialstates D = n ζ ( ~x ) : ζ ∈ L × L , ζ ( ~x ) ∈ e D , ∀ ~x ∈ Ω o . (10)The data-driven solution is again given by the minimization problem (6). We thus accept as solution the state thatconforms with the magnetostatic formulation and is best compatible with the available measurement data.To properly calculate the distance between two states, a suitable norm in the phase space must be defined. Thatnorm is given by || ζ || µ , ˜ ν = Z Ω (cid:20) e µ ~H · ~H + 12 e ν ~B · ~B (cid:21) dΩ , (11)where ζ = ( ~H, ~B ) ∈ Z = L × L and e µ ( ~x ), e ν ( ~x ) are tensorial and space-dependent weighting factors, respectively,and share units with the magnetic permeability and its inverse, the magnetic reluctivity, which explains the notation.We emphasize that the weighting factors e µ and e ν are of computational nature, i.e. they are not required to fulfillany physical material properties. Nevertheless, if appropriately chosen, they can improve the convergence of thedata-driven solver significantly. For the rest of this work, we consider only diagonal material tensors. However, fulltensor treatment is possible by decomposing the full tensor and applying rotation matrices, such that coordinatesystem and material tensor axes coincide. The distance function F ( ~H, ~B ) can now be written as F ( ~H, ~B ) = Z Ω f (cid:16) ~H, ~B (cid:17) dΩ = Z Ω (cid:20) (cid:16) ~H − ~H ? (cid:17) · e µ (cid:16) ~H − ~H ? (cid:17) + 12 (cid:16) ~B − ~B ? (cid:17) · e ν (cid:16) ~B − ~B ? (cid:17)(cid:21) dΩ , (12)where ~H ? and ~B ? refer to the measurement data from the global set of discrete states defined in (10). With thechosen weighting factors, the distance function (12) returns the magnetic energy mismatch between two states. Theminimization problem (6) can then be rewritten asminimize F (cid:16) ~H, ~B (cid:17) , (13a)subject to (cid:26) curl ~H = ~J, div ~B = 0 . (13b)To fulfill Gauss’s law (7b), we introduce the magnetic vector potential ~A , defined as ~B = curl ~A . Now we canformulate a suitable ansatz to solve (13b). The physical requirement for a divergence-free magnetic flux density ~B is fulfilled by incorporating the magnetic vector potential ansatz directly in (13a) and (13b), whereas Amp`ere’s law(7a) is enforced by means of Lagrange multipliers. The stationary problem then reads L (cid:16) ~H, ~A, ~η (cid:17) = F (cid:16) ~H, curl ~A (cid:17) − G (cid:16) ~η, ~H (cid:17) , (14)where the term G ( ~η, ~H ), which includes the Lagrange multiplier ~η ( ~x ), is given by G (cid:16) ~η, ~H (cid:17) = Z Ω ~η · (cid:16) ~J − curl ~H (cid:17) dΩ . (15)Our goal is to find the stationary points of (14). To that end, we compute the functional derivative of (14) withrespect to ~H, ~A , and ~η . Detailed information about functional derivatives can be found in [14, 41]. Since thefunctional derivative is linear, it is applied separately to (12) and (15). Writing the distance function element-wise,we obtain F (cid:16) ~H, curl ~A (cid:17) = Z Ω he µ x ( H x − H ?x ) + e µ y (cid:0) H y − H ?y (cid:1) + e µ z ( H z − H ?z ) i dΩ+ Z Ω he ν x ( ∂ y A z − ∂ z A y − B ?x ) + e ν y (cid:0) ∂ z A x − ∂ x A z − B ?y (cid:1) + e ν z ( ∂ x A y − ∂ y A x − B ?z ) i dΩ , (16)5here ∂ i ≡ ∂∂ i denotes the derivative with respect to the i -th component. The functional derivative of (15) withrespect to ~A is obviously vanishing. Then we must compute only the derivative of (16), which is given by δFδA x = ∂f∂A x − ∂ x ∂f∂ { ∂ x A x } − ∂ y ∂f∂ { ∂ y A x } − ∂ z ∂f∂ { ∂ z A x } = ∂ y { e ν z ( ∂ x A y − ∂ y A x − B ?z ) } − ∂ z (cid:8)e ν y (cid:0) ∂ z A x − ∂ x A z − B ?y (cid:1)(cid:9) (17) δFδA y = ∂f∂A y − ∂ x ∂f∂ { ∂ x A y } − ∂ y ∂f∂ { ∂ y A y } − ∂ z ∂f∂ { ∂ z A y } = ∂ z { e ν x ( ∂ y A z − ∂ z A y − B ?x ) } − ∂ x { e ν z ( ∂ x A y − ∂ y A x − B ?z ) } (18) δFδA z = ∂f∂A z − ∂ x ∂f∂ { ∂ x A z } − ∂ y ∂f∂ { ∂ y A z } − ∂ z ∂f∂ { ∂ z A z } = ∂ x (cid:8)e ν y (cid:0) ∂ z A x − ∂ x A z − B ?y (cid:1)(cid:9) − ∂ y { e ν x ( ∂ y A z − ∂ z A y − B ?x ) } , (19)where the differential operator δ refers to the derivative with respect to a function [14, 41]. Collecting all derivativesleads to ∇ A F = curl (cid:16)e ν curl ~A − e ν ~B ? (cid:17) = curl (cid:16)e ν curl ~A (cid:17) − curl (cid:16)e ν ~B ? (cid:17) , (20)where ∇ A denotes the functional gradient with respect to the components of ~A . In the same manner, the functionalderivative of (14) with respect to the Lagrange multiplier ~η is given by ∇ η L = ~J − curl ~H. (21)Finally, the derivative of (14) with respect to ~H is given as ∇ H L = ∇ H F − ∇ G L = e µ ( ~H − ~H ? ) + curl ~η. (22)Setting the derivatives to zero in order to solve for the stationary points of (14), we arrive at ∇ A L : curl (cid:16)e ν curl ~A (cid:17) − curl (cid:16)e ν ~B ? (cid:17) = 0 , (23a) ∇ H L : e µ ( ~H − ~H ? ) + curl ~η = 0 , (23b) ∇ η L : ~J − curl ~H = 0 . (23c) To solve the stationary problem (23), we employ the finite element method (FEM). Let us for now assume that aparticular state from D has been identified as the closest which conforms with the magnetostatic equations. Wedenote that state with ( ~H × , ~B × ).As commonly done in the setting of the FEM, to obtain the weak formulations, we first multiply equation (23a)with the test functions ~w ∈ V , where V is to be determined. Integration over the domain Ω yields Z Ω curl (cid:16)e ν curl ~A (cid:17) · ~w dΩ = Z Ω curl (cid:16)e ν ~B × (cid:17) · ~w dΩ . (24)Applying Green’s formula results in Z Ω e ν curl ~A · curl ~w dΩ = Z Ω e ν ~B × · curl ~w dΩ + Z Γ ( e ν curl ~A × ~n ) · ~w d S − Z Γ ( e ν ~B × × ~n ) · ~w d S, (25)6here ~n is the unit normal vector at Γ. The boundary is then split into a Neumann part Γ N and a Dirichlet partΓ D . For the rest of this paper, we consider only homogeneous boundary conditions, which means the Neumannconditions are naturally fulfilled, whereas the Dirichlet conditions are enforced strongly in the chosen function space.Hence, the boundary integral term vanishes. Next, we apply the curl operator on (23b) and insert equation (23c).Multiplying again with the test functions ~w and integrating over the domain Ω leads to Z Ω curl ( e ν curl ~η ) · ~w dΩ = Z Ω h curl ~H × · ~w − ~J · ~w i dΩ , (26)which, after applying Green’s formula, results in Z Ω e ν curl ~η · curl ~w dΩ = Z Ω h ~H × · curl ~w − ~J · ~w i dΩ + Z Γ (curl ~η × ~n ) · ~w d S − Z Γ (cid:16) ~H × × ~n (cid:17) · ~w d S. (27)We then apply the same boundary conditions as in (25). For notation purposes, we introduce the inner product( ~u, ~v ) Ω = Z Ω ~u · ~v dΩ , (28)where ~u, ~v ∈ L (Ω) . Now, we can express (25) and (27) using the bilinear forms a ( ~A, ~w ) := ( e ν curl ~A, curl ~w ) Ω , (29a) a ( ~η, ~w ) := ( e ν curl ~η, curl ~w ) Ω , (29b)where the right-hand sides (RHSs) are given by l × ( ~w ) = ( e ν ~B × , curl ~w ) Ω = Z Ω e ν ~B × · curl ~w dΩ , (30a) l × ( ~w ) = ( ~H × , curl ~w ) Ω − ( ~J, ~w ) Ω = Z Ω h ~H × · curl ~w − ~J · ~w i dΩ . (30b)Note that the superscript × on the functionals in (30a) and (30b) signifies that these functionals incorporate statesgiven by the material data set. Note that the RHSs (30) are of non-standard type as they feature a singular excitationterm through the measurement data.In contrast to the traditional approach, the data-driven formulation requires to solve two linear systems. However,the bilinear forms corresponding to (25) and (27) are identical, therefore, the system matrix must be assembled onlyonce. The weak formulations then read: Find ~A, ~η ∈ V = { ~v ∈ H(curl) : ~v × ~n = 0 on Γ D } such that a ( ~A, ~w ) = l × ( ~w ) , ∀ ~w ∈ V, (31a) a ( ~η, ~w ) = l × ( ~w ) , ∀ ~w ∈ V. (31b)Note that, in their current format, the problems (31) are not well-posed. Both vector fields ~A and ~η are only definedup to a gradient, i.e. they feature a large kernel given by the gradients of H(grad) functions. It follows, that thebilinear forms are not coercive, i.e. a ( ~u, ~u ) (cid:3) || u || . For instance, if u = gradΦ, then it directly follows that a ( ~u, ~u ) = 0, but also || u || = || gradΦ || L . There are several strategies to ensure that (31) is satisfied for all testfunctions, e.g. co-tree gauging [12, 9], current vector potential formulations [3, 42] to ensure a compatible RHS, orimposing the Coulomb gauge, to name but a few.After solving (31), we obtain an updated field solution. The magnetic flux density is directly computed after themagnetic vector potential ~A , whereas the magnetic field strength ~H is obtained by employing (23b). The updatedterms read ~B = curl ~A, (32a) ~H = ~H × + e ν curl ~η. (32b)7ext, we have to find the optimal data points in D , i.e. pairs ( ~H ?m , ~B ?m ) ∈ D that are closest to ( ~H ( ~x ) , ~B ( ~x )) ∀ ~x ∈ Ω.In this manner, we find the states in the space of ( ~H, ~B )-pairs that fulfill the constitutive relation, while at the sametime are closest to the magnetostatic formulation. The optimal pairs are determined with (cid:16) ~H × , ~B × (cid:17) = argmin ( ~H ? , ~B ? ) ∈D n F ( ~H, ~B ) o . (33)These new optimal states are taken into account to solve (31), followed by updating the field solution (32) and min-imizing (33). These computations are performed iteratively, until there is no change in ( ~H × , ~B × ) in two consecutiveiterations. At the beginning, the “optimal” state is randomly initialized within the material data set D .The data-driven algorithm in continuous (weak) formulation is given in Algorithm 1. initialize ( ~H × , ~B × ) randomly on D while convergence not reached do Find ~A, ~η ∈ V such that a ( ~A, ~w ) = l × ( ~w ) ∀ ~w ∈ Va ( ~η, ~w ) = l × ( ~w ) ∀ ~w ∈ V Update ~B = curl ~A~H = ~H × + e ν curl ~η Find ( ~H × , ~B × ) in D adjacent to ( ~H, ~B ) with (33). endAlgorithm 1: Iterative scheme for the data-driven magnetostatic field solver in continuous (weak) formulation.
For the rest of this work, lowest order finite elements (FEs) are employed. When discretizing the weak forms in(31), curl-conforming ansatz and test functions must be used for ~A and ~η . The magnetic vector potential and theLagrange multiplier are thus discretized as ~A ≈ ~a h = N e X i =1 a i ~N i ( ~x ) , (34a) ~η ≈ ~η h = N e X i =1 η i ~N i ( ~x ) , (34b)where ~N i are N´ed´elec basis functions of first kind defined on a triangulation of Ω [30], the subscript h denotes themaximum edge length of the triangulation, N e the number of degrees of freedom (DoFs) allocated on edges, and a i , η i the degrees of freedom. The Galerkin approach is used, such that ansatz and test functions are identical. Thealgebraic representation of the weak forms in (31) read K ˜ ν a = b × (35a) K ˜ ν η = b × (35b)where K ˜ ν ,i,j = ( e ν curl ~N j , curl ~N i ) Ω , (36a) b × ,i = ( e ν ~B × , curl ~N i ) Ω , (36b) b × ,i = ( ~H × , curl ~N i ) Ω − ( ~J, ~N i ) Ω . (36c)and a ∈ R N e , respectively η ∈ R N e are the DoFs of the corresponding vector fields in (34). Again, × indicatesthat material data was used to create the RHS. Furthermore, we introduce the diagonal matrices D Ω as well as the8atrices containing the weighting factors D ˜ ν , respectively D ˜ µ . The matrices are defined by D Ω ,d,e,e = Z T e d Ω , = V e (37) D ˜ ν ,d,e,e = e ν ( ~x e ) , (38) D ˜ µ ,d,e,e = e µ ( ~x e ) , (39)where d ∈ { x, y, z } , e ∈ { , . . . , N elem } , T e denotes the e -th element of the triangulation of Ω and ~x e refers to thecenter point of the corresponding elements. The discrete update terms are given by ~B h = curl ~a h (40a) ~H h = ~H × h + e ν curl ~η h . (40b)We will denote the values of the fields in (40) at the center points ~x e by B ◦ respectively H ◦ . Note that for lowest-orderFEs, one point per element is considered. In the case of higher order FEs, the fields are evaluated at the quadraturepoints of the FEs procedure [10]. To obtain the discrete counterpart of the distance function (12) we introduce thematrix induced scalar products h H , H i D ˜ ν = H > D Ω D ˜ ν H , (41) h B , B i D ˜ µ = B > D Ω D ˜ µ B , (42)for H , B ∈ R N elem . Now, the discrete minimization problem reads (cid:0) H × , B × (cid:1) = argmin ( H ? , B ? ) ∈D F ( H ◦ , B ◦ )= argmin ( H ? , B ? ) ∈D h H ◦ − H ? , H ◦ − H ? i D ˜ ν + 12 h B ◦ − B ? , B ◦ − B ? i D ˜ µ . (43)The distance minimization is carried out for each dimension separately. In practical, real-world applications, it is not uncommon that some regions of the computational domain Ω consistof materials with exactly known properties, such that there exists a closed-form relation between ~H and ~B . Forexample, this situation is encountered in models where iron parts and air-gaps coexist. A straightforward strategy isto provide a measurement data set D containing all possible points obeying the closed-form relation, being an infinitenumber. The data-driven computing framework can then be modified accordingly. In that case, the weighting factorin the domain with known constitutive law should be chosen to be equal to the now known material coefficients.A previous work of the authors proposes three modifications of increasing intrusiveness to the data-driven com-puting framework, in order to deal with mixed models, i.e. with exact and data-based materials co-existing in thecomputational domain [10]. In this work, we shall stick to the first, least intrusive approach, where the exact materialrelation is treated by the data-driven solver itself. Under this approach, the distance function (43) is updated withthe known material coefficients µ ex and ν ex , and the solution of the minimization problem is given by B × = B ◦ + D µ ex H ◦ , (44a) H × = D ν ex B × , (44b)where D ν ex , D µ ex are matrices following the definition in (38), respectively (39). The field solution is then updatedwith (44) at the same spot as the field solution for the parts with unknown material relation. Note that, during theassembly of K ˜ ν and in order to update the terms for the magnetic field strength, a mixed weighting factor has tobe considered. Thus, the weighting factor e ν now contains exactly known material coefficients for the domains withknown material law and weighting factors for the purely data-based domains.Other known material laws can be treated in the same manner. Nevertheless, the relations (44) only hold for thecase of linear materials. For the nonlinear case, the corresponding relations have to be derived by minimizing (43)in view of the nonlinear material law. 9 .6 2D magnetostatic weak formulation In many engineering applications, it is possible to reduce the more accurate but computationally challenging three-dimensional formulation to a less demanding two-dimensional one. For example, such a reduction is possible if thegeometry of the considered device remains almost unchanged along one certain direction. Consequently, the fieldcomponent in that particular direction can be regarded as negligible. In the 2D case, the magnetic vector potential isreduced to one spatial component. Without loss of generality, we choose the z -component, such that ~A = (0 , , A z ),respectively ~η = (0 , , η z ). The edge functions can then be chosen to be ~N = Φ( x, y ) ~e z , (45)where Φ are standard nodal functions defined on a triangulation of the 2D cross-section of Ω and ~e z denotes the unitvector in z -direction. Employing the 2D basis function (45) as test and trial functions in (31) leads to two linearsystems to be solved. Note that the constructed 2D basis function (45) have zero divergence, meaning an additionalgauging can be omitted. For more details on solving the 2D magnetostatic problem, we refer the reader to [10]. Itshould also be pointed out that the structure of Algorithm 1 and its discrete counterpart remain unchanged. Until now, the choice of the weighting factors e ν within the context of data-driven solutions has not been addressed.In contrast to conventional FEM solvers, the factor e ν does not necessarily represent the physical behavior of thematerial, but is rather an element of computational nature. Nevertheless, choosing the weighting factors in anappropriate way can increase the convergence rate of the data-driven algorithm substantially. In the following, wediscuss two different possibilities, i.e. global and local weighting factors. Further, we assume that the measurementpairs ( ~H ?m , ~B ?m ) , m = 1 , . . . , M are noise-free, such that the observed values are unaffected by measurement errorsor other uncertainty sources. The case of noisy data is discussed in Section 4.The most naive choice for the weighting factors would be e ν = 1, respectively e µ = 1. However, calculating adistance in the ( H, B ) phase space already requires stretching and squeezing of the axes because the values for H and B differ by orders of magnitude. Moreover, the ratio between H and B strongly depends on the material, whichwas the issue addressed in [10] and necessitates heterogeneous weighting factors. Furthermore, the ratio between H and B is strongly affected by the nonlinearity, which is our major concern here, and is counteracted by adaptivelyupdating e ν and e µ .In the case of a nonlinear material, the constitutive law reads H ( B ) = ν ( B ) B . At a most recently obtainedstate ( H ◦ , B ◦ ), we can distinguish between two possible linearizations, corresponding to the successive-substitutionmethod and the Newton method, respectively. The corresponding operations points read H − H ◦ = ν ◦ c ( B − B ◦ ) , (46) H − H ◦ = ν ◦ d ( B − B ◦ ) , (47)where ν ◦ c = H ◦ B ◦ , (48) ν ◦ d = d H ◦ d B ◦ , (49)are the chord and differential reluctivity, respectively [11] (Figure 1b). The linearization is carried out for eachFE quadrature point and forces a reassembly the the FE matrices in each iteration step. For the problems un-der consideration, up to 100 (successive-substitution) or up to 10 (Newton) iteration steps are needed to obtainconvergence.The data-driven solver does not dispose of a material law and, thus, inherently does not incorporate an iterationfor treating the nonlinearity. One could consider the data-driven iteration as a replacement thereof. The updates(40) and (43) then correspond to the action of linearization in a traditional solver. Hence, it may not wonder thatparticular choices for e ν and e µ proposed below, show similarities to the standard successive-substitution and Newtonmethods. 10or the rest of this work, we determine the weighting factors according to the differential reluctivity. Note thatone could also consider the chord reluctivity, which still results in an improvement in convergence compared tothe original method. However, the results cannot compete with the ones obtained when employing the differentialreluctivity. When restricting to a global weighting factor the results differ only marginally when using the chord orthe differential reluctivity. A straightforward method to get an estimation of the global weighting factor e ν , is to apply finite differences on themeasurement data D to obtain the differential reluctivity. We consider here the case of anisotropic materials anddecoupled fields, such that H d depends only on B d , where d ∈ { x, y, z } refers to the spatial dimension. Then, thereluctivity tensor has only diagonal elements and reads e ν = e ν x e ν y
00 0 e ν z , (50)where e ν x , e ν y , e ν z are considered to be constant in the case of a global weighting factor. We assume the measurementdata to be sorted in ascending order such that ( H ?d,m , B ?d,m ) < ( H ?d,m +1 , B ?d,m +1 ) , m = 1 , . . . , M −
1, holds element-wise, i.e. H ?d,m < H ?d,m +1 and B ?d,m < B ?d,m +1 . Assuming that the data are noiseless and sampled equidistantly, wecan apply the centered finite differences scheme to obtain ν ( B ), such that ν d ( B d ) = d H d d B d ≈ ν d,m ( B ?d,m ) = H ?d ( B ?d,m + ∆ B ?d ) − H ?d ( B ?d,m − ∆ B ?d )2∆ B ?d , for m = 2 . . . M − , (51)with ∆ B ?d = B ?d,m +1 − B ?d,m = const. , ∀ m = 1 . . . M −
1. To obtain a weighting factor on the boundaries, i.e. m = 1,respectively m = M , forward or backward differences can be employed. In the case of non-equidistant data, we canapply standard forward differences, thus obtaining ν d,m ( B ?d,m ) = H ?d ( B ?d,m + ∆ B ?d ) − H ?d ( B ?d,m )∆ B ?d , with ∆ B ?d = B ?d,m +1 − B ?d,m , for m = 1 . . . M − . (52)The stiffness matrix (36a), the RHS (36b), and the distance function (43) are then initialized with the mean overthe discrete reluctivity curve. Thus, when employing (52), the elements of the reluctivity tensor are given as e ν d = E [ ν d ] , with ν d = [ ν d, , . . . , ν d,M − ] > (53)where E [ X ] = N P Nn =1 x n . Averaging over all differential reluctivities is of course only reasonable if the data pointsare almost equidistant. Particularly for data sets with clusters, a weighted average is necessary. Note that, when usinga global weighting factor, we assign each element the same computational constant, i.e. one constant represents theentire BH -curve. The algorithm employing one global weighting factor is shown in Algorithm 2, with the exceptionof the gray shaded part which is reserved for the local weighting factors assignment, discussed next. The global weighting factor introduced in Section 3.1 faces major difficulties to represent the material behavior incases of material curves with large variations. We consider here the case of soft magnetic materials whose BH -curvefeatures a steep linear part and a sharp transition into saturation. However, we expect that similar material behaviorcan be observed for non-EM material properties as well. For ferromagnetic materials, the relative reluctivity in thelinear part is typically in the region of ≈ − , whereas in the saturation part its value lies in the region of ≈ − until it furthers increases to 1 for full saturation. It is therefore obvious that an averaged differential reluctivity is notcapable to cover the BH -curve in its entirety, as also illustrated in Figure 2, due to the fact that a global weightingfactor covers only one operation point in the BH -curve. For instance, consider a region in the computational domainwith sharp edges and thus high field densities. Consequently, the working point on the BH -curve is expected to belocated in the saturation part, equivalently, the local reluctivity is expected to be large compared to the reluctivity inthe linear part. However, if the averaged reluctivity is predominantly attributed to the linear part of the BH -curve,the field states ζ ◦ = ( H ◦ , B ◦ ) ∈ M which should otherwise be located in the saturation part, are still matched tothe linear part. 11 . . . . . . . , , B [T] H (cid:2) A m (cid:3) . . . . . . · − d H d B ν , H B ν Figure 2: Nonlinear BH -curve, relative differential reluctivity d H d B ν and relative chord reluctivity HB ν .As a remedy to this problem, we propose to allocate a local weighting factor per finite element. We denote theheterogeneous weighting tensor with b µ , respectively b ν for reluctivity. Since we have no a priori knowledge aboutthe working points of each element, the assignment of local weighting factors must be performed adaptively duringthe data-driven iteration. The factors are initialized in the same manner as described in Section 3.1. During theiterative process, we monitor the convergence of the distance between states ζ ◦ ∈ M fulfilling Maxwell’s equationsand the states ζ × = ( H × , B × ) ∈ D in the measurement data set, when minimizing the energy mismatch (cid:15) em ,d = N elem X e =1 V e (cid:20) e µ d,e (cid:16) H ◦ d,e − H × d,e (cid:17) + 12 e ν d,e (cid:16) B ◦ d,e − B × d,e (cid:17) (cid:21) , (54)for each dimension d ∈ { x, y, z } , where V e refers to the volume of element e .A stagnation of the energy mismatch (54) during iteration can be attributed to two reasons. First, for elementswhose BH -operating point is matched to the global weighting factor, the states ζ ◦ are likely to be close to somemeasurement point in D . Second, for elements whose BH -operating point is far off the global weighting factor,the states ζ ◦ are most probably distant to the measurement points and further iterations cannot contribute to animprovement, which is illustrated for an example case in Figure 8a. Once detected, this stagnation triggers a switchfrom global to local weighting factors, where the local reluctivity is estimated using the local material solution, i.e.with the states ( H × , B × ) selected by the energy mismatch minimization. Note that it is also possible to switchfrom global to local weighting factors after a certain number of iterations. However, some iterations with a globalaveraged weighting factor is preferable, since the states ζ × are randomly initialized on D and consequently an earlylocal treatment could estimate factors that are not optimal. The local weighting factors are given by b ν d ( ~x ) = N elem X e =1 ν d (cid:16) B × d,e (cid:17) e ( ~x ) with d ∈ { x, y, z } , (55) e ( ~x ) = (cid:26) , ~x ∈ T e , ~x T e , (56)where T e refers to the e -th element of the triangulation of Ω. Keeping these considerations in mind, we adapt thedata-driven algorithm as shown in Algorithm 2. Now, each element has a specific local weighting factor b ν d ( ~x e ),estimated with (51) or (52). This process seems similar to the material constant assigned with a standard Newtonsolver. However, we point out that, in contrast to a Newton solver, we do not explicitly model the material and,more importantly, we do not rely on a physical representation of the material. Since the weighting factors are nowchosen optimally with respect to the current local field solution for each element, we can expect an improvement inthe convergence of the data-driven solver. 12igure 3: Cross-section of a quadrupole magnet. Thegray area is iron, the red areas are the coils, and thewhite area is vacuum. The shaded polygon with theblue boundary shows the computational domain. B n | Γ D = 0 H t | Γ N = 0Figure 4: Computational domain and mesh accordingto one eighth of the quadrupole magnet. On the blueboundary, homogeneous Dirichlet boundary conditionsare imposed, whereas on the black boundary homogen-eous Neumann conditions are imposed. initialize ( H × , B × ) randomly on D estimate e ν on D while convergence not reached do Solve variational problem (35) for a and η Update ~B h = curl ~a h ~H h = ~H × h + e ν curl ~η h Find ( H × , B × ) in D adjacent to ( H ◦ , B ◦ ) with distance function (11) only for local weightingfactor assignment if stagnation detected then assign local reluctivity b ν to elements (55)assemble stiffness matrix K ˆ ν update distance function (11) with local reluctivity b ν endendAlgorithm 2: Data-driven magnetostatic field solver. The gray part shows the extension according to the localweighting factors.
For our numerical investigations, we employ the model of a quadrupole magnet. Such magnets are used in acceleratorfacilities to focus particle beams. Due to translational invariance along the z -direction, the magnet can be modeledin 2D without significant compromises in terms of modeling accuracy. The 2D magnet model is depicted in Figure 3and consists of three different domains, namely, the iron yoke (gray), the coils (red), and vacuum (white). Exploitingthe symmetry of the model and assuming that the permeability of the yoke is sufficiently high, such that the magneticflux leaving the yoke can be neglected, it is sufficient to consider only one eighth of the quadrupole’s geometry, shownin Figures 3 (marked area) and 4. The iron yoke has an anisotropic permeability, which is linear in the y -directionand nonlinear in the x -direction, following the Brauer model [5].The material responses of the coils and vacuum are considered to be exactly known, thus, those regions aretreated as described in Section 2.5. The magnetic material is fully treated by the data-driven solver as discussed inSection 2.4. The permeability of the yoke follows the behavior of a soft magnetic material in the x -direction. Dueto their small hysteresis, soft magnetic materials are widely used in inductors, electrical machines, and other devicesto reduce losses. Moreover, soft magnetic materials feature strongly nonlinear magnetization curves. As illustratednext, the original data-driven algorithm proposed in [23], equivalently, Algorithm 2 in this paper, has difficulties inproviding accurate solutions, due to the fact that the global weighting factor e ν which applies to the entire domain13 must fit not only to the linear part of the BH -curve, but also to the saturation part as well. On the contrary,the data-driven solver with local weighting factors is able to resolve the nonlinearity of the material law, providingsignificant efficiency gains.Both data-driven solvers, i.e. respectively utilizing global or local weighting factors, are employed to providea data-driven solution of the quadrupole’s magnetic field. The data points in D are generated by sampling thematerial curves in x - and y -direction equidistantly. To verify the results obtained with either data-driven field solver,we introduce three error metrics. The first error quantifies the energy mismatch between the most recently computedstates ζ ◦ ∈ M fulfilling the governing equations and the reference solution, and is defined as (cid:15) em = h H ◦ − H ref , H ◦ − H ref i D ν ref + h B ◦ − B ref , B ◦ − B ref i D µ ref h H ref , H ref i D ν ref + h B ref , B ref i D µ ref ! (57)The other two metrics are the relative field errors (cid:15) H = h H ◦ − H ref , H ◦ − H ref i D ν ref h H ref , H ref i D ν ref ! , (58) (cid:15) B = h B ◦ − B ref , B ◦ − B ref i D µ ref h B ref , B ref i D µ ref ! , (59)with respect to the magnetic field strength ~H and the magnetic flux density ~B . Note that D µ ref and D ν ref refer to thepermeability and reluctivity coefficients obtained from the reference solution, which are computed with a standardNewton solver.To monitor the convergence of the errors defined in (57), (58), and (59), the data-driven solvers are employedwith measurement data sets D of increasing size. The convergence behavior of both data-driven solvers with respectto all three errors is shown in Figures 5a and 5b. It can easily be observed that, for the same data set, the data-drivensolver utilizing local weighting factors is orders of magnitude more accurate than the original, with respect to allthree error metrics. In particular, local weighting factors lead to a linear convergence rate. Moreover, the accuracygain due to the use of local weighting factors becomes even more pronounced as the size of the data set increases.10 − − − − − − O ( N − ) Number of data points E n e r g y m i s m a t c h (cid:15) e m globallocal (a) Convergence of the energy mismatch for increasing dataset size. − − − − − − − Number of data points R e l a t i v ee rr o r s (cid:15) H , (cid:15) B (cid:15) H global (cid:15) B global (cid:15) H local (cid:15) B local (b) Convergence of the relative errors in the ~H - and ~B -fieldfor increasing data set size. Figure 5: Convergence of data-driven solutions with respect to data set size. The solid lines refer to the data-drivensolver utilizing local weighting factors, while the dashed lines refer to the standard data-driven solver based on aglobal weighting factor.The efficiency gains due to the local weighting factors become even more obvious if we additionally observe thenumber of outer-loop iterations performed by each data-driven solver, again considering data sets of increasing size.14he convergence of the standard data-driven solver with a global weighting factor in terms of outer-loop iterationsis shown in Figure 6. The convergence behavior of the modified solver utilizing local weighting factors is shown inFigure 7. As can be observed, switching from global to local material assignment, as suggested in Algorithm 2, leadsto a substantial improvement in terms of computational demand, now quantified by the number of solver iterations.This is especially true in the case of large data sets, i.e. with 10 or more points. − − − N = 10 , N = 10 N = 10 N = 10 Number of iterations E n e r g y m i s m a t c h (cid:15) e m (a) Convergence of the energy mismatch for increasingsolver iterations over data sets of increasing size N . − − N = 10 , N = 10 N = 10 N = 10 Number of iterations R e l a t i v ee rr o r (cid:15) H (b) Convergence of the relative error in the ~H -field for in-creasing solver iterations over data sets of increasing size N . Figure 6: Convergence of the data-driven solution with respect to the number of solver iterations for the case of aglobal weighting factor. The dots indicate that convergence is reached. − − − − − − N = 10 N = 10 N = 10 N = 10 N = 10 global localNumber of iterations E n e r g y m i s m a t c h (cid:15) e m (a) Convergence of the energy mismatch for increasingsolver iterations over data sets of increasing size N . − − − − − − N = 10 N = 10 N = 10 N = 10 N = 10 global localNumber of iterations R e l a t i v ee rr o r (cid:15) H (b) Convergence of the relative error in the ~H -field for in-creasing solver iterations over data sets of increasing size N . Figure 7: Convergence of the data-driven solution with respect to the number of solver iterations for the case oflocal weighting factors. The dots indicate that convergence is reached. The algorithm switches from global to localmaterial assignment after 20 iterations. Note that the horizontal scale differs from the one in Figure 6.This improvement in accuracy and efficiency is to be expected due to the fact that the x -component of the BH -15urve of the soft magnetic material comprising the yoke features a strongly nonlinear behavior. To better illustratethe deficiency of the standard data-driven solver and the advantages of local material assignment, Figure 8 showsthe x -component of the field solution (black points), i.e. the states ( H ◦ x , B ◦ x ) ∈ M identified by the solver, in theiron part of the quadrupole. The available material data are depicted as well (blue points). Both figures refer to thefield solution after each algorithm has converged. As predicted, due to the large differences in the reluctivity andits wide value range, the standard data-driven solver based on a global weighting factor fails to capture the materialresponse, in particular regarding the saturation region, as shown in Figure 8a. On the contrary, the data-drivensolver based on local weighting factors yields states that accurately match the material data and the underlyingconstitutive relation, thus leading to the accuracy and efficiency improvements shown in Figures 5 and 7. − − . . · − − H x (cid:2) Am (cid:3) B x [ T ] (a) Global weighting factor e ν . − − . . · − − H x (cid:2) Am (cid:3) B x [ T ] (b) Local weighting factors b ν . Figure 8: ( H x , B x )-pairs corresponding to measurement data (blue points) and to the data-driven field solution aftersolver convergence (black points). Both data-driven simulations have been carried out with N = 10 data points. In real-world applications, it is unreasonable to assume that the measurement data are free of noise. Measurementnoise is introduced by varying ambient conditions, tolerances of the measurement device, or other factors thatcan influence the quality of the measurement. A common assumption is that the added noise follows a Gaussiandistribution. Here we assume independent and identically distributed Gaussian noise with zero mean and variances σ H , respectively σ B , added on the measurement data. A noisy measurement sample m is then defined as( H ?m , B ?m ) = ( ˆ H ?m + ε H , ˆ B ?m + ε B ) , (60)where ˆ H ? , ˆ B ? are the noiseless measurements and ε H ∼ N (0 , σ H ), ε B ∼ N (0 , σ B ) denote the Gaussian noise terms.The data-driven methods presented and employed in Section 3 are not suitable for the solution of problems withnoisy data due to the fact that the distance-minimization approach is highly sensitive to outliers [24]. That is, sinceall data points are considered of equal importance, the solver could choose an outlying point, i.e. one far awayfrom material data clusters, which however conforms better to a solution of the Maxwell equations. Nevertheless,data clusters signify that the corresponding data points are more probable to represent the underlying material law.Clustered data points should then be considered of increased importance compared to outliers in the data set. It istherefore of relevance to attach a measure of importance to the single data points. This is done by first performinga clustering analysis on the measurement set D . To identify the clusters, a scheme based on maximum entropyestimation [19, 22] and simulated annealing is proposed in [24]. The same idea is used in the present work also,albeit combined with the local weighting factors introduced in Section 3.2. We briefly recall the basic idea behindthe maximum entropy estimation and refer the reader to [24] for more details.16 .1 Global weighting factor We assume that a local state ζ e ∈ M fulfilling the magnetostatic equations in an element e of the triangulation ofΩ is available. Prior to performing the distance minimization as in Section 3, the measurement data are weightedaccording to their importance. A weighting factor p e,m is assigned to each data point ζ ?m = ( H ?m , B ?m ), where it holdsthat M X m =1 p e,m = 1 , (61a) p e,m ∈ [0 , . (61b)The weights should be unbiased, while points distant from the local state ζ e should have small weights compared tocloser ones. A distribution satisfying these constraints is given by p e,m ( ζ e , β ) = Z e ( ζ e , β ) − e − β || ζ e − ζ ?m || µ , ˜ ν , (62a) Z e ( ζ e , β ) = M X m =1 e − β || ζ e − ζ ?m || µ , ˜ ν , (62b)where the constant β controls the influence of data points distant to the local state ζ e . The distance weighted center¯ ζ × e per element is then given by ¯ ζ × e = (cid:0) ¯ H × e , ¯ B × e (cid:1) = M X m =1 p e,m ( ζ e , β ) ζ ?m . (63)In the context of machine learning, (62a) is known as the softmax function [34], whereas in statistical mechanics itis connected to the Boltzmann distribution. The normalization constant (62b) is known as the partition function.Applying this weighting, the data-driven solver selects clusters in the β -controlled neighborhood of ζ e , instead ofoutliers. Effectively, the weight p e,m represents the probability that the data point ζ ?m is close to the local state ζ e .For instance, if the measurement point m is matched to the local state, the weight for that data point is p e,m = 1and all other weights vanish due to (61a), i.e. p e,j = 0 , ∀ j = 1 , · · · , M ∧ m = j . Note that in the case of β → ∞ ,the classical distance minimization scheme is recovered [24].The extension of the data-driven algorithm is then straightforward. Instead of employing the distance minimiz-ation function (43), we seek for the states ζ e ∈ M that minimize the free energy in phase space, given by ζ opt ,e ∈ argmin ζ e ∈M { P e ( ζ e , β ) } , with P e ( ζ e , β ) = − β − log ( Z e ( ζ e , β )) . (64)This procedure is then applied to all local states, i.e. all ( H , B )-pairs. Consequently the total free energy is given by P ( ζ, β ) = N elem X e =1 P e ( ζ e , β ) = N elem X e =1 − β − log ( Z e ( ζ e , β )) . (65)A bottleneck of performing the minimization with respect to the total free energy, is that the function (65) can bestrongly non-convex [24]. Hence, an iterative solver could converge to a local minimum. However, P ( ζ, β ) is convexfor a sufficiently small value of the constant β . The issue can then be resolved by initializing β with a suitablevalue such that P ( ζ, β ) is convex, and then increasing that value based on simulated annealing [26], until the globalminimum is reached, as suggested in [24]. In contrast to the noiseless case, the data in the noisy case are disordered, making it less trivial to assign thedifferential reluctivity to a local state ζ e . In the case of ordered data, it is straightforward to apply the tools fromdifferential calculus. Otherwise, a beforehand data analysis is essential.Depending on the size of our measurement data set D , we create K ∈ N clusters, such that the ( H, B )-pair ineach cluster share a “similar” differential reluctivity. To obtain the clusters, well-known K-means algorithms can be17 1 ,
000 2 ,
000 3 ,
000 4 ,
000 5 , . . H x (cid:2) Am (cid:3) B x [ T ] Figure 9: Noisy BH -data set in blue. K-means clustersin orange. 0 1 ,
000 2 ,
000 3 ,
000 4 ,
000 5 , , , , H x (cid:2) Am (cid:3) µ x / µ Figure 10: Estimated µ x , allocated at K-means clustercenters.employed [1, 4, 21, 28]. The reluctivity of each cluster is then approximated with the linear coefficient of a standardlinear regression model. Figure 9 shows exemplarily a measurement set with noisy data in blue and the clustersin orange. The corresponding estimation of the permeability µ is shown in Figure 10. To estimate the differentialreluctivity or permeability of each cluster, a regression technique that is robust to outliers is needed. To that end,we employ Huber regression [17]. In contrast to conventional regression, points are classified as inliers and outliersand outlier points receive a linear weight that reduces their importance. We note that, depending on the size andquality of the measurement data set, nonphysical solutions for the local permeability can occur, which will then haveto be corrected.Particularly for noisy data, we can take advantage of the fact that the weighting factors e µ and e ν do not necessarilyhave physical values in the data-driven context. Instead, the factors serve as indicators for the operating regionsof each element. Hence, in the noisy case, the data-driven algorithm utilizing local weighting factors is identicalto the structure of Algorithm 2, meaning that the manipulation of e µ and e ν is embedded in the maximum-entropydata-driven algorithm suggested in [24]. Algorithm 3 shows the embedding of the local weighting factors assignmentinto the maximum entropy data-driven solver. 18 nitialize cluster centers ( ¯ H × , ¯ B × ) and β [24] estimate e ν on the clustered measurement set D while convergence not reached do Calculate weighted centers ( ¯ H × , ¯ B × ) [24]Find ( H × , B × ) in D adjacent to ( ¯ H × , ¯ B × ) with distance function (11)Solve variational problem (35) for a and η Update ~B h = curl ~a h ~H h = ~H × h + e ν curl ~η h Update β according to [24] only for local weightingfactor assignment if stagnation detected then assign local reluctivity b ν to elementsassemble stiffness matrix K ˆ ν update distance function (11) with local reluctivity b ν endendAlgorithm 3: Maximum entropy data-driven magnetostatic field solver. The gray part shows the extensionaccording to the local weighting factors.
The model under consideration is again the quadrupole magnet presented in Section 3.3. The maximum-entropydata-driven solver is applied both with global and local weighting factors. The noisy measurement data sets aregenerated by sampling the material curves equidistantly. The measurement data are then polluted with additiveGaussian noise according to formula (60), where the chosen variances are σ H = 10 A / m and σ B = 0 .
04 T. Figure 9shows exemplarily the noisy measurement set in x -direction for N = 10 data points. Similar to the numericalinvestigations of Section 3, we monitor the convergence of the data-driven solvers for data sets of increasing sizeusing the errors (57), (58), and (59). Additionally, to estimate the variation of the data-driven solutions due to noiseand assess the robustness of each data-driven solver, 100 different noisy data sets are generated for each data setsize.Using data sets of increasing size, the convergence of the data-driven solutions with respect to the average energymismatch, computed over the 100 different data sets, is shown in Figure 11a. The same figure also depicts the interval ± σ around the average energy mismatches (in gray). The convergence with respect to the average relative fielderrors (cid:15) H and (cid:15) B , respectively defined in (58) and (59), is presented in Figure 11b. It is obvious from both figuresthat the maximum-entropy data-driven solver modified with local weighting factors results in accuracy and efficiencygains, improving the convergence rate of the method by a factor slightly greater than 2. Compared to the noiselesscase, i.e. looking at the results of Section 3.3, the convergence of both solvers is affected significantly by the additionof noise in the data. The impact of noise on the convergence of data-driven solvers has already been documented in[24] for the case of a global weighting factor. This impact is also present when local weighting factors are utilized.Note that Kirchdoerfer and Ortiz presented a convergence order of O ( N − . ) for the case of noisy data. However,as already discussed in Section 3, the convergence rate is here additionally affected by the strong nonlinearity in theresponse of the soft magnetic material. Hence, the slower convergence rate is to be expected.The convergence analysis with respect to the number of solver iterations is shown in Figures 12a and 12b. Bothfigures show the average energy mismatch when global (dashed lines) and local (solid lines) weighting factors areutilized. Figure 12b shows additionally the ± σ intervals around the average energy mismatch for measurementssets of size N = 10 and N = 10 . Compared to the noiseless case, see Section 3.3, a substantial decrease in thenumber of iterations can be observed when the data-driven solver with a global weighting factor is employed. On thecontrary, when local weighting factors are used, the number of iterations is comparable to the noiseless case. Thus,contrary to the noiseless case, assigning local weighting factors does not result in a significant improvement in termsof solver iterations. This difference can be attributed to the simulated annealing schedule, due to the fact that thesolver iterations are now predominantly determined by controlling the constant β , see also [24].The standard deviation σ of the energy mismatch in dependence to the size of the data sets is shown in Figure 13aand in dependence to the number of solver iterations in Figure 13b. Two main conclusions can be drawn from the190 − . − − . − . − . O ( N − . ) O ( N − . ) Number of data points E n e r g y m i s m a t c h globallocal (a) Convergence of the average energy mismatch for increasingdata set size. The shaded areas show ± σ around the averageerrors. − − Number of data points R e l a t i v ee rr o r s (cid:15) H , (cid:15) B (cid:15) H global (cid:15) B global (cid:15) H local (cid:15) B local (b) Convergence of the average relative errors in the ~H - and ~B -field for increasing data set size. Figure 11: Convergence of data-driven solutions with respect to data set size for noisy measurement data. The solidlines refer to the data-driven solver utilizing local weighting factors, while the dashed lines refer to the standarddata-driven solver based on a global weighting factor.0 20 40 60 80 10010 − − global localNumber of iterations E n e r g y m i s m a t c h N = 10 , global N = 10 , global N = 10 , global N = 10 , local N = 10 , local N = 10 , local(a) Convergence of the energy mismatch for increasing solveriterations over noisy data sets of increasing size.
20 40 60 80 10010 − − Number of iterations E n e r g y m i s m a t c h N = 10 , global N = 10 , global N = 10 , local N = 10 , local(b) Convergence of the energy mismatch for increasing solveriterations over noisy data sets of increasing size. The shadedareas show ± σ around the average errors. Figure 12: Convergence of data-driven solutions with respect to the number of solver iterations for noisy measurementdata. The solid lines refer to the data-driven solver utilizing local weighting factors, while the dashed lines refer tothe standard data-driven solver based on a global weighting factor.presented results. First, both figures show that the standard deviation decreases with larger data sets when localweighting factors are used. Thus, as should be expected, the robustness of the method is improved when larger datasets are available. Contrarily, when a global weighting factor is used, no improvement in the variance is observedwhen increasing the data set size from N = 10 to N = 10 . Second, for small data sets of size N = 10 , the standarddeviation for the modified data-driven solver is slightly worse than the one of the standard solver. Local weighting200 − . − − . − . Number of data points S t a nd a r dd e v i a t i o n σ globallocal (a) Standard deviation of the energy mismatch for noisy datasets of increasing size. − − global localNumber of iterations S t a nd a r dd e v i a t i o n σ N = 10 , global N = 10 , global N = 10 , global N = 10 , local N = 10 , local N = 10 , local(b) Standard deviation of the energy mismatch for increas-ing solver iterations over noisy data sets of increasing size. Figure 13: Standard deviation of the data-driven solutions for noisy measurement data. The solid lines refer to thedata-driven solver utilizing local weighting factors, while the dashed lines refer to the standard data-driven solverbased on a global weighting factor.factors still offer an advantage for larger data sets and should be the method of choice in that case. Nevertheless, inboth cases, the differences are only marginal.Overall, based on the variance improvement for larger data sets, as well as because of the convergence gains, localweighting factors are advantageous in the noisy case also. However, the advantages are not as pronounced as in thenoiseless case.
This work presented a data-driven, material-model-free, magnetostatic FEM solver, which addresses the challengingcase of highly nonlinear material responses. The data-driven computing framework was extended with the introduc-tion of local weighting factors, which indicate local operation points in the nonlinear domain. The novel approachwas applied for the case of noiseless data by modifying an existing data-driven magnetostatics solver presented in[10], based on the distance-minimization scheme first suggested in [23]. Likewise, for the case of noisy data, localweighting factors are embedded in the maximum-entropy data-driven algorithm proposed in [24]. The modifieddata-driven solvers were then compared to the standard data-driven solvers based on a global weighting factor.The numerical results presented in Section 3.3 for the noise-free case illustrate that the proposed modification tothe data-driven algorithm leads to an impressive improvement in terms of convergence, accuracy, and computationalefficiency. Therefore, in the case of noiseless data, local weighting factors should be preferred over the standardapproach. The paper shows that, for models with a pronounced nonlinearity in the material behavior, the useof local weighting factors is essential. The results also suggest that the two data-driven approaches of increasedintrusiveness suggested in [10] can be discarded altogether and replaced seamlessly by the approach presented in thiswork.In the case of noisy measurement data, the proposed modification of the maximum-entropy data-driven solverwith local weighting factors improves the convergence rate of the solver by a factor of 2. The improvements inaccuracy and efficiency become more pronounced for larger data sets, i.e. when more than 10 data points areavailable. Nevertheless, the addition of noise affects significantly the convergence behavior of both standard andmodified data-driven solvers, therefore, the advantages due to local weighting factors are not as clear as in thenoiseless case. Similar conclusions are drawn when monitoring the variance of the data-driven solution, which isan indicator of the robustness of each data-driven method. In particular, the data-driven solutions based on local21eighting factors show a reduced variance for increasing data set size, whereas a global weighting factor-based solversseem to stagnate. However, the variance differences are marginal, thus, the robustness of both methods can be seenas comparable. Overall, the modification of the data-driven solver with local weighting factors is advantageous inthe noisy case also, particularly if data sets with N = 10 or more data points are available. Acknowledgment
D. Loukrezis and H. De Gersem would like to acknowledge the support of the Graduate School of Excellence forComputational Engineering at the Technische Universit¨at Darmstadt. D. Loukrezis further acknowledges the supportof the BMBF via the research contract 05K19RDB. A. Galetzka’s work is supported by the DFG through theGraduiertenkolleg 2128 “Accelerator Science and Technology for Energy Recovery Linacs”.
References [1] D. Arthur and S. Vassilvitskii. K-means++: the advantages of careful seeding. In
In Proceedings of the 18thAnnual ACM-SIAM Symposium on Discrete Algorithms , 2007.[2] J. Ayensa-Jim´enez, M. H. Doweidar, J. A. Sanz-Herrera, and M. Doblar´e. A new reliability-based data-drivenapproach for noisy experimental data with physical constraints.
Computer Methods in Applied Mechanics andEngineering , 328:752–774, 2018.[3] O. Biro, K. Preis, W. Renhart, G. Vrisk, and K. R. Richter. Computation of 3-d current driven skin effectproblems using a current vector potential.
IEEE Transactions on Magnetics , 29(2):1325–1328, March 1993.[4] J. Bl¨omer, C. Lammersen, M. Schmidt, and C. Sohler.
Theoretical Analysis of the k-Means Algorithm – ASurvey , pages 81–116. Springer International Publishing, Cham, 2016.[5] J. Brauer. Simple equations for the magnetization and reluctivity curves of steel.
IEEE Transactions onMagnetics , 11(1):81–81, 1975.[6] A. Buffa and R. Hiptmair. Galerkin boundary element methods for electromagnetic scattering. In
Topics incomputational wave propagation , pages 83–124. Springer, 2003.[7] S. Conti, S. M¨uller, and M. Ortiz. Data-driven problems in elasticity.
Archive for Rational Mechanics andAnalysis , 229(1):79–123, 2018.[8] S. Conti, S. M¨uller, and M. Ortiz. Data-driven finite elasticity.
Archive for Rational Mechanics and Analysis ,pages 1–33, 2020.[9] E. Creus´e, P. Dular, and S. Nicaise. About the gauge conditions arising in finite element magnetostatic problems.
Computers & Mathematics with Applications , 77(6):1563 – 1582, 2019. 7th International Conference on AdvancedComputational Methods in Engineering (ACOMEN 2017).[10] H. De Gersem, A. Galetzka, I. G. Ion, D. Loukrezis, and U. R¨omer. Magnetic field simulation with data-drivenmaterial modeling.
IEEE Transactions on Magnetics , 56(8):1–6, 2020.[11] H. De Gersem, I. Munteanu, and T. Weiland. Construction of differential material matrices for the orthogonalfinite-integration technique with nonlinear materials.
IEEE Transactions on Magnetics , 44(6):710–713, 2008.[12] P. Dular, A. Nicolet, A. Genon, and W. Legros. A discrete sequence associated with mixed finite elements andits gauge condition for vector potentials.
IEEE Transactions on Magnetics , 31(3):1356–1359, May 1995.[13] R. Eggersmann, T. Kirchdoerfer, S. Reese, L. Stainier, and M. Ortiz. Model-free data-driven inelasticity.
Computer Methods in Applied Mechanics and Engineering , 350:81–99, 2019.[14] I. Gelfand and S. Fomin.
Calculus of Variations . Dover Publications Inc., 1963.[15] A. Gupta, A. Cecen, S. Goyal, A. K. Singh, and S. R. Kalidindi. Structure-property linkages using a data scienceapproach: application to a non-metallic inclusion/steel composite system.
Acta Materialia , 91:239–254, 2015.2216] B. Heise. Analysis of a fully discrete finite element method for a nonlinear magnetic field problem.
SIAM Journalon Numerical Analysis , 31(3):745–759, 1994.[17] P. Huber and Ronchetti.
Regression , chapter 7, pages 149–198. John Wiley & Sons, Ltd, 2009.[18] R. Ibanez, E. Abisset-Chavanne, J. V. Aguado, D. Gonzalez, E. Cueto, and F. Chinesta. A manifold learn-ing approach to data-driven computational elasticity and inelasticity.
Archives of Computational Methods inEngineering , 25(1):47–57, 2018.[19] E. T. Jaynes. Information theory and statistical mechanics.
Phys. Rev. , 106:620–630, May 1957.[20] J.-M. Jin.
The Finite Element Method in Electromagnetics . Wiley-IEEE Press, 3rd edition, 2014.[21] L. Kaufmann and P. Rousseeuw. Clustering by means of medoids.
Data Analysis based on the L1-Norm andRelated Methods , pages 405–416, 1987.[22] H. K. Kesavan.
Jaynes’ Maximum Entropy Principle , pages 1779–1782. Springer US, Boston, MA, 2009.[23] T. Kirchdoerfer and M. Ortiz. Data-driven computational mechanics.
Computer Methods in Applied Mechanicsand Engineering , 304:81 – 101, 2016.[24] T. Kirchdoerfer and M. Ortiz. Data driven computing with noisy material data sets.
Computer Methods inApplied Mechanics and Engineering , 326:622 – 641, 2017.[25] T. Kirchdoerfer and M. Ortiz. Data-driven computing in dynamics.
International Journal for Numerical Methodsin Engineering , 113(11):1697–1710, 2018.[26] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi. Optimization by simulated annealing.
Science , 220(4598):671–680, 1983.[27] A. Leygue, M. Coret, J. R´ethor´e, L. Stainier, and E. Verron. Data-based derivation of material response.
Computer Methods in Applied Mechanics and Engineering , 331:184–196, 2018.[28] S. P. Lloyd. Least squares quantization in PCM.
IEEE Trans. Inf. Theory , 28:129–136, 1982.[29] P. Monk et al.
Finite element methods for Maxwell’s equations . Oxford University Press, 2003.[30] J. C. N´ed´elec. Mixed finite elements in R . Numerische Mathematik , 35(3):315–341, 1980.[31] C. Pechstein and B. J¨uttler. Monotonicity-preserving interproximation of B-H-curves.
Journal of Computationaland Applied Mathematics , 196(1):45–57, 2006.[32] K. Rajan.
Informatics for materials science and engineering: data-driven discovery for accelerated experiment-ation and application . Butterworth-Heinemann, 2013.[33] K. Rajan. Materials informatics: The materials “gene” and big data.
Annual Review of Materials Research ,45:153–169, 2015.[34] C. Rasmussen and C. Williams.
Gaussian Processes for Machine Learning . Adaptive Computation and MachineLearning. MIT Press, Cambridge, MA, USA, 2006.[35] J. Rikabi, C. F. Bryant, and E. M. Freeman. An error-based approach to complementary formulations of staticfield solutions.
International Journal for Numerical Methods in Engineering , 26(9):1963–1987, 1988.[36] L. Stainier, A. Leygue, and M. Ortiz. Model-free data-driven methods in mechanics: material data identificationand solvers.
Computational Mechanics , 64(2):381–393, 2019.[37] T. J. Sullivan.
Introduction to uncertainty quantification , volume 63. Springer, 2015.[38] A. S. Veeramani, J. H. Crews, and G. D. Buckner. Hysteretic recurrent neural networks: a tool for modelinghysteretic materials and systems.
Smart Materials and Structures , 18:075004 (15pp), 2009.2339] T. Weiland. Time domain electromagnetic field computation with finite difference methods.
InternationalJournal of Numerical Modelling: Electronic Networks, Devices and Fields , 9(4):295–319, 1996.[40] T. Weiland. Finite integration method and discrete electromagnetism. In
Computational Electromagnetics ,pages 183–198. Springer, 2003.[41] R. Weinstock.
Calculus of Variations: With Applications to Physics and Engineering . Dover Publications, 1974.[42] Zhuoxiang Ren. Influence of the rhs on the convergence behaviour of the curl-curl equation.