Computing Forces by ECSW-Hyperreduction in Nonlinear Magnetodynamic FEM Problems

This contribution proposes a novel strategy to apply the hyperreduction principle for the mechanical force calculation of non-linear magnetodynamic problems in an FEM context. The already established energy-conserving sampling and weighting (ECSW) method is therefore adapted to exploit the structure of mechatronic systems, so the hyperreduction is only applied to the non-linear domain. It is consequently named partial hyperreduction, here the pECSW. Also, the force calculation is formulated such that ECSW is applied as sampling-based integration. The force computation requires a new element sampling and weighting process, here called cECSW, arising from the configurational updates of the geometry corresponding with the possible motions of the mobile part. The advantage is the possibility of avoiding treating the full-order model for the force calculation. Storing the reduced order model with two element sets for field and force calculation is sufficient. The strategy is demonstrated on the TEAM20 problem of the Compumag benchmark systems, which is extended to dynamic excitation. The new methods show the ability to produce accurate results for significantly smaller systems.


I. INTRODUCTION
A NALYZING electromagnetic devices using the finite element method (FEM) is a standard procedure.Using non-linear material laws to consider the saturation of the magnetic field is also well-established in the industry.A recent topic of research is the model order reduction (MOR) of non-linear magnetodynamic field simulations.Promising results were achieved by the hyperreduction approach, namely, the energy-conserving sampling and weighting (ECSW) method [1].

A. Related Work
In the realm of computational engineering, MOR and hyperreduction, a term introduced in [2], play pivotal roles in managing complex systems.The proper orthogonal decomposition (POD) is a common MOR technique [3], successfully applied to quasi-static magnetodynamic systems [4], [5].An overview of projection-based methods is given in [6].POD and Krylov subspace techniques were also investigated for linear material in the frequency domain [7], [8].
However, the straightforward projection into reduced space can be insufficient in addressing non-linear problems, where much computational cost arises from evaluating each element's non-linear contribution to the system's energy.This issue has been the focus of many studies [9], [10], [11], with the application of the discrete empirical interpolation method Manuscript received 16 May 2023; revised 25 September 2023; accepted 6 November 2023.Date of publication 13 November 2023; date of current version 10 January 2024.Corresponding author: J. Maierhofer (email: j.maierhofer@tum.de).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TMAG.2023.3332210.
Digital Object Identifier 10.1109/TMAG.2023.3332210(DEIM) [12].In [13], a method is introduced for reducing non-linear magnetodynamic problems by combining POD with manifold interpolation.This approach has demonstrated superior results compared to classical techniques like direct POD reduction and standard interpolation of pre-computed reduced bases.
In order to improve the accuracy and robustness of reduction for non-linear models, the ECSW method, originating from an idea of computer graphics [14], further developed for structural dynamics [15] has emerged as a promising candidate.The ECSW method might offer significant advantages in electrodynamic systems regarding stability and result quality, even for highly reduced models.Initial steps have been made in recent studies [1], [16], [17].In these works, the ECSW method was successfully utilized to compute the magnetic vector potential, yielding promising results, indicating its ability to significantly reduce computational costs and retain good accuracy.
Applying the ECSW method for electromagnetic problems related to mechatronic applications (actuators and sensors) offers the possibility to exploit the ECSW further than has already been done so far.Indeed, since only part of the domain is non-linear (due to saturation), the implementation of ECSW can be significantly improved.Furthermore, in applications where forces need to be derived from the electromagnetic simulation, the ECSW idea can be exploited to significantly reduce the cost of the post-processing step.In this article, we propose new ideas to exploit ECSW for those two important aspects.

B. Objective
As monitoring mechatronic systems with the use of socalled digital twins in embedded controllers becomes more and more important, demand for small, yet sufficient, models with a small footprint in terms of memory and CPU load rises.Two main objectives, therefore, drive the MOR strategy in this article.
1) Save computational time in the transient, non-linear field computation by investing only minimal offline effort.2) Generate a fully self-contained minimal model capable of computing dynamic forces for any current trajectory applied.In all generality, the system's structure is shown in Fig. 1.The air domain (constant µ, no conductivity) encapsulates the excitation domain (copper coils) 0 , where the current density j 0 is given.The core c and the mobile target t are domains with non-linear magnetic permeability [µ(B)] and constant conductivity σ .Both can, therefore, induce eddy currents for a changing magnetic field.The mobile target domain t is moved by the virtual displacement δs to evaluate the mechanical force.

C. System Equations
The discretized equations for electrodynamic systems are founded on Maxwell's equations.For quasi-static problems, the frequency of time-varying fields is low enough so that one can neglect electromagnetic wave emission and a simplification of Maxwell's equations ( 2) and (1) to the so-called magnetodynamic equations can be considered The isotropic constitutive law for the material (3) and Ohm's law (4) complete the governing equations The concept of the vector potential A is defined as and inserted in Faraday's law (1) the magnetodynamic equations are solely written in terms of the vector potential1 Applying the finite element theory with edge type elements to (7) yields the classical discretized first time-order differential equation set (8), with u being the unknown values of the vector potential, M the magnetic mass matrix, g the internal current vector and f the external load (current) vector A more detailed derivation of the formulation used here can be found in [19] and [20].
The internal current density g(u) is derived via the symbolic derivation of the magnetic energy density with respect to the vector potential.The magnetic energy density is obtained by integrating the B H -curve of the material's datasheet After the derivation, the weak form is built, and the discretization is applied.
As a time integration scheme, a backward Euler scheme in combination with a Newton-Raphson method is considered.With r being the residual

II. HYPERREDUCTION APPROACH
The ECSW Method weights selected elements such that the internal currents g(u) produce the same virtual work as the non-reduced system (analogous to what was proposed in the seminal paper [15] where structural dynamic problems were considered and g(u) were non-linear internal mechanical forces).The ECSW method is a combination of a Galerkin projection on a subspace of the solution to reduce the degrees of freedom and a hyperreduction to reduce the non-linear internal current term to a small evaluation set as shortly summarized next.

A. Galerkin Projection
The full solution vectors of a reference simulation span the space of the expected solution.Here, in a pre-processing step, a solution of the full problem is computed while applying an excitation similar to the one applied in further simulations and a singular value decomposition (SVD) of the solution snapshots U = [u 1 , . . ., u m ] is performed.The reduction basis V contains the most important left singular vectors of that decomposition, and the solution is then approximated by u ≈ V q (11) where q are the generalized degrees of freedom in the reduced space.The reduced problem, denoted with the index r , is then obtained by projection of (8) on V Although this projection reduces the size of the problem, the computational burden to advance the simulation in time is only marginally reduced since the internal currents g(V q) and the tangent matrix (dg/du) still need to be evaluated in every iteration step (9), and (10).Hence an additional so-called hyperreduction step is needed, using, for instance, the ECSW approach.

B. Energy Conserving Sampling and Weighting
The method was developed for non-linear finite element dynamic models in mechanical engineering in [15] and [21] around 2014.The goal is to evaluate only a subset Ẽ of elements E and weigh them with a factor ζ e > 0.
The sum of all virtual work contributions, produced in each element by the internal currents, should be approximated as accurately as possible.Therefore, considering the contribution of each element (localization vector L e ) to the reduced current, one can write where and where ζ is only non-zero for the selected set Ẽ, a subset of the full element set E. The approximation ( 14) should be as good as possible for any q.To find a ζ with as many zero entries as possible which approximates the sum of the virtual work to a defined tolerance, m training sets (index s) q s , obtained during a pre-simulation (resulting u s ) and a projection with the basis V are considered as explained next. 1) Training Sets and sNNLS: Considering the full-order solution vectors, filtered with the basis, at each computed time step t s , one can compute the contribution V T L T e g e (L e V q s ) of each element to the reduced force V T g(V q s ).This information is then used to find a minimum set of non-zero ζ e to approximate the reduced currents V T g(V q s ) up to a given tolerance τ . . .
In addition, to guarantee the positive definiteness of the associated energy, all weightings in ζ must be positive.A quasi-optimal set ζ can be found by a sparse non-negative least square (sNNLS) solver.More details on this algorithm can be found in [21].2) Evaluated Elements: As the unassembled internal currents are used to calculate the virtual work of each element, the weighting factors ζ are also at the level of unassembled elements.In Fig. 2, the active elements are marked in blue and their associated dofs in red.Due to the use of Nédélecelements, the edges of the elements represent the dofs of the system.

C. Hyperreduction
After finding the desired weighting vector ζ , the hyperreduced internal currents vector g r,ECSW and tangential stiffness matrix K r,ECSW are given as: This strategy allows to significantly reduce the cost of the simulation with the reduced (12) [1].

III. PARTIAL ECSW
In general, major parts of magnetodynamic systems can be modeled linearly.Only small parts of the domain need a nonlinear consideration.To reduce the order of the model, the linear domain can be sufficiently reduced by projection, while in the non-linear domain, an additional hyper-reduction step is needed as explained previously.Applying the hyperreduction method ECSW only to a part of the domain will be called partial ECSW, short pECSW in the following.Applying pECSW saves offline computational effort as the sNNLS only iterates over a much smaller subset of elements and also avoids a second approximation of the linear domain, Fig. 3.
Realizing the partial ECSW, two subvariants are imaginable.1) Domain decomposition, i.e., split the domain into a linear and a non-linear subdomain and hyperreduce the non-linear force function of the non-linear subdomain.2) Order decomposition of internal currents, i.e., split non-linear internal currents in the entire domain into linear and higher order terms and hyperreduce only higher order terms.

A. pECSW-Domain Decomposition
The first approach is inspired by substructuring techniques where structural systems are partitioned into subsystems, see Fig. 4. The individual subsystems have different sets of equations that are re-assembled after having reduced the subdomains.The difference to the classic ECSW formulation is in the assembly process of the Y matrix where the contribution of each element to the internal force is collected.The full set of elements E now already is a subset of elements, i.e., the elements in the non-linear domain, denoted by an Asterisk, E ⋆ , to distinguish from the subset stemming from the full elements, denoted by a tilde, Ẽ.The selection process for the elements is the same as in the classic ECSW except that in (14), only the elements from The equations for the backward Euler time stepping are adapted accordingly to consider the linear and hyperreduced non-linear domains separately.The problem (8) can thus be solved at t n+1 as As for the classic ECSW, the Newton algorithm needs to solve the next solution step in reduced coordinates q with the corresponding stepping matrix K step,red (q i ) from the time integration scheme.In contrast to the previous stepping matrix algorithm, the reduced stepping matrix is now a combination of the hyperreduced K r,pECSW and the projected K lin .The linearized problem for a given time step at Newton iteration i + 1 writes where

B. pECSW-Order Decomposition
A second way to separate linear and non-linear behavior is interpreted as a local modification technique.As a baseline, a linear system is assumed with a constant permeability that fits the material curve for small magnetic fluxes.The nonlinear domains are then modified with an additional term, which represents the non-linearity.This additional term now contains all the higher-order terms of the magnetic flux density without the linear order.In Fig. 5, the gap between the assumed linear and the true non-linear material law is visualized.
While the matrix representation of the domain decomposition approach (see Fig. 4) shows a splitting of the equations (vertical splitting), the order decomposition splits the force contributions (horizontal splitting).The finally hyperreduced internal current vector is marked in yellow, as it is a modified function, see Fig. 6.
As the method relies on the unassembled internal currents of the finite elements, two weak forms, respectively, for linear (resulting in a matrix K 1 ) and non-linear material (resulting in g(u)) are formulated and evaluated on the mesh.The difference is consequently then the higher order internal currents g ⋆ and is used for the training of the hyperreduction algorithm Finally, depending on whether an element is in a linear domain or in a non-linear domain, its internal force will be computed by The further procedure is analogous to the previously shown pECSW method, i.e., projection of the internal currents in the chosen basis, then solving the next solution vector in generalized coordinates using the modified stepping matrix.

C. pECSW Conclusion
This extension of the ECSW method is not primarily intended to reduce the number of computed elements further to gain lower computational effort in online calculations.Rather, the objective is to increase the trust in the results as the hyperreduction is only applied where needed (namely, only in the non-linear part, leaving the linear part without hyper-reduction treatment).A nice benefit in terms of lower effort arises in the sNNLS algorithm that selects the elements.By only allowing a subset of elements, the algorithm becomes faster.

IV. CONFIGURATIONAL ECSW
Due to the so-called skin effect, representing eddy currents correctly usually requires a very fine mesh in the vicinity of the domain boundary.The size depends on the frequency of the time-varying fields and its material properties.This element size is an absolute value, which means it does not scale with the overall model's size.The bigger the modeled system, the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.more elements are needed.However, calculating global forces would not require a high element resolution, which leads to the idea of reducing the number of evaluated elements for the post-processing step of calculating the mechanical forces of magnetodynamic systems.
The approach called configurational ECSW here, in short cECSW, is an adaption of the ECSW idea to the postprocessing.It is, therefore, not involved in the actual solution process and is thus independent of the method used for the computation of the magnetic field.In the presented case, the method is aimed at computing the magnetic force on a mobile ferromagnetic target.The cECSW is trained by different configurations of the geometry, which represent the virtual displacement needed for the force computation according to the virtual work principle (shortly explained in Section IV-A).The underlying assumption is that the weights of the selected elements in the computation of the mechanical force is not dependent on the configuration.Two configurations for a horizontal displacement of the target are shown generically as an example for the discussion in Fig. 7.
1) Remark: A strict distinction must be made between parametric and configurational geometry changes.The parametric geometry changes are real changes in the geometry with no limit to their size, as typically encountered in modification loops during design optimization or when performing the analysis by changing the position of the bolt in a solenoid actuator or the rotor rotation in an electric machine.These changes are not covered in this contribution, as they would not fulfill the assumption that the weights of the selected elements are independent of the configuration.The configurational changes handled here are only virtual changes, which indicate what kinematically admissible motions the body could undergo.The mesh topology must stay constant.It is common to just deform one layer of elements, as these virtual changes can be very small, even for a finite difference evaluation as seen in Section IV-A.

A. Force Calculation
The global mechanical force on a mobile target body is derived in general from the total change of energy W due to the change of the system's configuration (i.e., position) with an amplitude denoted by s [19].For an applied, constant current density j , the force energetically conjugated to the considered configurational change (e.g., resulting horizontal force in Fig. 7) is obtained from the change of magnetic coenergy.Considering several possible configurational changes (e.g., vertical and rotational motion) and calling s the set of their amplitudes, the associate conjugate forces are obtained by derivation of the magnetic energy For quasi-static systems (as claimed during the derivation), the force can be computed for any timestep purely from the conservative system.Note that the current density j is not only the externally applied current but also the eddy current distribution for each timestep.The equivalent static problem including the eddy current, defined by M u, is Fig. 6. pECSW-order decomposition approach.The magnetodynamic equation remains one over the entire domain.The non-linear behavior is split into a linear order term, represented as a matrix, and a higher order term, which undergoes the hyperreduction process.Let us recall (see Section I-C) that the finite element formulation searches for a discrete vector potential u such that it minimizes the energy functional of the system.
Defining the relation for the energy W □ = f T u = W + W ′ , which is the sum of energy W and co-energy W ′ , the solution of the magnetic problem satisfies Hence, around the solution of the problem, the co-energy W ′ does not change with an infinitesimal change of the vector potential.Written in general terms before the discretization Using then the chain rule to evaluate the mechanical force ( 24) The global co-energy is obtained by integrating the co-energy density w ′ ( A), which is the integral under the curve of the constitutive material law, Fig. 5.With the material being assumed to be isotropic, B and H are always collinear which allows to write (28) in terms of magnitudes Performing the partial derivations in (27), note that, by definition, the remaining quantities are constant.The expression for the force becomes The change of position of the part leads to a change of the Jacobian, which describes the coordinate transformation between physical coordinates and the isometric coordinates for the FEM elements.The Jacobian of the isogeometric mapping is denoted by J , and its determinant by |J | = det(J ).In the first term of (29), one needs to evaluate the change of B (at constant A) due to a deformation of the initial coordinate system by the perturbation.This can be obtained using the Jacobian of the transformation J from the initial configuration to the perturbed one [22].In all generality, the force arising from the co-energy is written for all quantities given in the original configuration [23] Change due to change of derivation operator curl An overview of the different ways to use (30) for force calculation on a given mesh is given in Fig. 8 and will be discussed next.After a successful field calculation per timestep, the vector potential is further processed to gain the global mechanical forces.While the classical Jacobian method needs all the virtually deformed elements (marked yellow in Fig. 8), the cECSW approach proposed in the following paragraphs samples and weighs those elements in either a nonlinear or a linear way, which reduces the evaluated sub-subset of elements to the green marked ones in Fig. 8.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 8. cECSW scheme.The field solution provides a time series of magnetic vector potentials.In the post-processing, besides the classical Jacobian Method, the cECSW is proposed.Only a subset of virtually deformed elements is considered for the force integration.

B. Force by Classical Jacobian Method
For elements that are not deformed during the virtual displacement of the mobile target, the derivative of the Jacobian with respect to the displacement is zero.That results in a subset of elements that do not contribute to the global force.Therefore, only a subset of elements needs to be evaluated.The integral over the domain changes into integrals over each element e and a sum over all elements in the subset

C. Force Integration by ECSW
Applying the idea of the ECSW method for the force calculation to reduce the computational effort and shrinking the model size means selecting elements and calculating weights.The sum over the element forces (the set of elements which undergo a deformation is denoted with E ′ ) is approximated by a sum over a subset of elements E ⋆ which are then weighted with the factor ζ ⋆ e to compensate for the neglected elements.The localization operator L e picks the solution dofs for the corresponding element e, and thus With the energy being positive, the idea of ECSW can be applied, approximating the energy with positive weighted elements.Assuming that the weights are constant with respect to the virtual deformation, the gradient is shifted behind the weight ζ ⋆ e .This is a reasonable assumption, considering that the virtual configurational change is infinitesimal e (L e u) d e (33) In order to select the elements using the sNNLS algorithm, training snapshots need to be prepared.Considering the m solution snapshots of a pre-simulation, the change of co-energy in each element due to a set of c possible configurational changes is computed.This generates a set of mc training sets from which a minimum set of elements and the corresponding weights can be selected following the ECSW strategy.
In order to avoid the evaluation of the Jacobian derivatives for each training dataset for each element, the gradient is approximated with a finite difference.
1) Non-Linear Material Elements: In the general case, where the deformed elements contain non-linear material, the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
training is performed using the co-energy of every element in E ′ , W ′ e (L e u).The configurational finite steps used for the finite difference are denoted with the superscript c0 and c1.The subscript denotes the element number 1, . . ., n and the snapshot number 1, . . ., mc.Using the same amplitude s for each direction of configurational change allows us to omit the scaling in the following explanations.The adapted Y matrix for the force approximation to be used for the sNNLS algorithm becomes In general, the co-energy for the mc different snapshots need to be integrated as in (28).As explained next, this costly computation can be avoided if the elements are linear.
2) Linear Material Elements: If all the virtually deformed elements are of linear material, a simplification during the preparation of the sNNLS can be conducted.This situation occurs for systems where the non-linear material is completely encapsulated in air and only rigid body motions are allowed.For real application systems, this is a very common assumption.A few transformations are carried out to extract the co-energy directly from the internal currents already known from the training simulations.
The linear energy density is written as the product of the internal current density g(u) and the solution vector u.The solution vector does, by definition, not change with a configuration change, so the derivative is only effective on the internal currents.Written for each element, the force becomes (L e u) T (L e g) (37) With the same justification as for the classic ECSW, the elements can be sampled and weighted with positive values, but now the difference of the internal currents for different configurations is approximated.For each training snapshot i j, i = 1, . . ., m, j = 1, . . ., c, hence each line in (36), an equation for the sNNLS is added This shows that in case the elements in E ′ are linear, the force can be evaluated by cECSW from (39) with a reduced element set and weights obtained from the training with (40) where no co-energy needs to be evaluated, but only known internal currents obtained from the configurational steps are needed.The element forces F e can be evaluated either by finite differences or the local Jacobian method (30).Note that, if using (30), it is sufficient to save the weights for Fig. 9. Selected elements for both stages marked in blue and green.Blue elements serve for the field calculation and green elements for the postprocessing.Some elements are in both subsets.Fig. 10.Geometry TEAM20 [24].
the evaluated elements and their Jacobian relations only.The weights, once computed, stay constant after the training, while the Jacobian values depend on the actual (virtual) configuration change.

D. cECSW Conclusion
The main idea is to find an optimized reduced model that holds for the field calculation and for the post-processing, i.e., calculating forces.For this, the ECSW method is applied as sampling and weighting of elements.This saves in the firstline model storage and, as a nice benefit, computational time for the post-processing.

V. EXAMPLE
The previously presented strategy of force calculation utilizing the ECSW method in both stages is demonstrated on the TEAM20 benchmark problem of the Compumag Society.As the original problem definition is only static, it is extended to a magnetodynamic problem by assigning a homogeneous conductivity to all ferromagnetic materials.The dynamic parameters are inspired by the TEAM10 problem, which handles eddy currents in non-linear systems of a similar structure.
The 3-D geometry is defined as seen in Fig. 10.It consists of three parts.A ferromagnetic yoke structure embraces a copper coil.Inside the coil, a cuboid of the same ferromagnetic material as the yoke is mounted such that it is freely movable in the vertical direction.Applying a current, reluctance forces emerge that tend to pull the free body in.
The solids and the surrounding air are meshed with in total 65 740 tetrahedal elements, see Fig. 11.A finite element order  (fe-order) 1 leads to 76 564 dofs, fe-order 2 results in 338 984 dofs.

A. Material
The non-linear B H -curve of the steel is given in a tabular form, which is plotted in Fig. 12.A spline interpolation is performed to obtain a continuous function that is integrable in order to find the energy density.For the pECSW-Order-Decomposition, the first section of the curve is fit with a linear permeability of µ r = 1500.To compare the performance, two further linear permeabilities are considered which do not fit the actual material well, µ r ,low = 750 and µ r ,high = 3000.The conductivity is given as σ = 7.505 × 10 6 S m −1 .
1) Framework: The work for this contribution was carried out using NGSolve and Python.NGSolve is an open-source, object-oriented finite element library written in C++ and equipped with a comprehensive Python API.It is closely linked to NETGEN [25], an automatic mesh generator that can mesh 2-D and 3-D geometries created with constructive solid geometry (CSG) or open cascade technology (OCT).Computations were conducted using an Intel Xeon CPU E3-1270 v5 at 3.60 GHz, with NGSolve version 6.2.2301 and Python 3.11.
2) Relative Error: The relative error (RE) compares reduced problem solutions u to their corresponding reference u ref .
In the course of this work, this comparison is limited to scalar values (field quantities evaluated locally or global values).The error at each time step t i is summed to produce Fig. 13.Evaluation points P 1 and P j evaluate the field quantity B and j eddy .An average is built over the cross section area defined by the lines α-β and γ -δ.
a single significant value.The formula is The RE definition is quite harsh, and a small offset over multiple time steps results in a high error number.But it is also a useful one-number criterion to see how well two measurement histories match.
3) Evaluation Points: The simulation values are evaluated using two local points and the average results on two cross section areas.The benchmark is also evaluated both in simulation and experiment at point P 1 , as well as the lines α-β and γ -δ, see [26].Point P j , see Fig. 13, is introduced in this article to measure the eddy current values that counteract the coil's driving current.

B. Reference Simulation
For the static reference, the coil is excited by dc-currents with a resulting applied magnetomotive force between 1000 and 5000 AT.With a cross section of 1738.8 × 10 −6 m 2 this results in the current densities of j = 0, 575(A/mm 2 ) and j = 2, 876(A/mm 2 ), respectively.
The magnetic field at the point P 1 and the averaged magnetic field along two lines α-β and γ -δ as well as the force of the plunger in z-direction (see Fig. 14) coincide well with the measurement from the official results [27].
The dynamic simulations are conducted with an excitation function that forms a smooth step j 0 0.575 A mm −2 to 2.876 A mm −2 τ 0.05 s.The parameters are chosen such that the steel parts can be sufficiently saturated and the eddy current density is not small, and yet the eddy current in the coil can be neglected [24].
As a timestep, two different values (t = 0.001 s and t = 0.005 s) are used to test the convergence of the simulation.Both reveal nearly the exact same behavior, so the timestep t = 0.005 s is chosen.

TABLE I RELATIVE ERROR DUE TO THE PROJECTION (WITHOUT ANY HYPERREDUCTION) IN THE SUBSPACE SPANNED BY POD MODES
No further excitation functions to train the system are presented in the course of this contribution.Further details on different strategies for generating training data can be found in [17].
1) Regularization: Following the theory of modified vector potential, the non-conductive domains lead to a non-unique governing set of equations.This singularity is overcome by a small regularization term that is added to the equation like a very small mass matrix, i.e., a very small conductivity in the order of 1 × 10 6 (S/m).
For the sNNLS, the regularization term is neglected, although this is actually a small discrepancy between the computed vector field and the evaluated energy densities.

C. Field Computation-pECSW Results
To reduce the system with the proposed two-stage method, first, a good field computation is aimed.Therefore, the entire system needs to be computed, and the resulting magnetic vector snapshots are fed into the SVD to find a solution basis.In this case, ten modes are sufficient.
The POD is computed from the first 20 timesteps.Projecting the solution on a subspace of {5, 10, 20} modes shows that the results at the evaluated points can be approximated within 1 % by the use of ten modes.A beginning of overfitting can be detected for 20 modes trained from 20 snapshots, see Table I.Applying the presented variations of the hyperreduction method ECSW leads to model sizes, which are a fraction of the full-size model.As the speedup is directly related to the number of evaluated elements, it is not explicitly discussed here.Of course, some data-handling overhead remains constant through all the simulations.The timings given in Table II, are the wall-times of the computation of one timestep, averaged over all timesteps.The system's element count decreases from 65 740 to a range between 10 and 100.And nearly with the same ratio do the computation times reduce for the evaluation of the elements.Table II also shows the offline times to select the relevant elements and their corresponding weights.The offline time reduces significantly due to the pECSW approach.The finite elementorder influences the computational effort, but the rates of reduction itself are quite independent.
1) Regularization: Due to the projection into the subspace, the regularization of the air domain is not necessary anymore.The standard ECSW method does not violate the well-posed property of the problem.The final stepping matrix unveils condition numbers around 5-25, which leads to easy and robust solutions in the linear solver.The pECSW-DD can be upgraded with a regularization term for the linear domain, which lowers the condition numbers but is not necessary.The pECSW-OD does not promise a bounded condition number.There were cases where the stepping matrix became singular after a few time steps.This behavior could be hindered but not consequently inhibited with a regularization term in the linear domain.
2) Constant Number of Elements: The plot in Fig. 15 shows the magnetic flux density B z at P 1 and the eddy current density j x at P j for a ramp-like excitation.The results of the standard ECSW, pECSW domain decomposition (pECSW-DD), and pECSW order decomposition (pECSW-OD) are plotted for two reduction levels (11 and 55 elements).These element numbers correspond to the number of elements the standard ECSW chose for τ = 0.1 and τ = 0.001.All three methods now have the same computational effort, and a fair comparison of the capability of the reduced model is possible.
The takeaway message from Fig. 15 indicates that both pECSW methods perform better than the standard ECSW.This was expected as for the pECSW all evaluated elements contribute to the approximation of the non-linearity.The order decomposition method works the best.For higher element numbers considered in the ECSW hyperreduced model, the differences between the methods vanish, and all the results match the reference.However, the pECSW-OD has the major drawback of becoming potentially unstable due to a bad matrix condition.
3) Conclusion: Both pECSW approaches fulfill the expectations in terms of building a reduced-order model to compute the magnetic vector potential.Discussing whether it is necessary or even worth applying the pECSW approach instead of the standard ECSW, which considers the full domain despite its non-linearity, two points are to be highlighted as follows.
a) Implementation effort: If an intrusive implementation is possible, the extra effort for pECSW is very small as all the techniques are already there.
b) Offline costs: The offline costs are remarkably smaller (10-20 times smaller) for the pECSW approaches.This leads to the strong recommendation to use a pECSW method.c) Convergence: The convergence of the pECSW-OD is not always guaranteed, and therefore, the method is not recommended, although it would provide the best results when stable.
In summary, the most promising implementation of the ECSW method to compute the field quantities is the classic ECSW, which performs very well and does not need any further implementation complexity but needs some more elements to exhibit similar accuracy.Adding one step of complexity, i.e., the partitioning of the system, the pECSW-DD is the most intuitive and robust hyperreduction approach.

D. Force Calculation-cECSW Results
After investigating the field calculation, the application of the ECSW approach to post-processing in terms of force calculation is discussed in the example.The resulting element selections are visualized in Fig. 16.
Due to the ramp excitation, the final force value converges to the static forces.For a single rigid body mode of the Fig. 16.TEAM20 selected elements.Blue elements for the pECSWalgorithm, green elements for the cECSW-algorithm.Non-selected air elements and coil elements are invisible.Fig. 17.Comparison the global force from standard ECSW and cECSW.After t = 0.2 s, the system is not in a steady state and, therefore, is not fully to the static value.
mobile target in the z-direction, 14 elements were selected for a sNNLS tolerance of τ = 1 × 10 −3 in the ECSW for the forces.The offline costs were around t = 1 s.The online costs for a single force integration decrease from 0.7 to 0.002 s.The force results are very accurate; refer to Fig. 17.The final excitation current at t = 0.2 s is 98 % of the static reference of 5000 AT.

E. Reduced Order Model-Results
In the following, the reduced model created beforehand (pECSW-DD for the field calculation and cECSW for the force, trained once with the ramp seen above) will be tested for three different excitations against their full-order reference solution.The results are visualized in the radar chart in Fig. 19.Each ray represents a relative error measure (local and global).The shorter the circular path, the better.
The first excitation is the same ramp as for the training, see (42).This leads to the so-called onboard validation of the reduced model.The two following trajectories are unknown to the system and consequently represent a cross-validation of the reduced order model.The second excitation function Cross validation of the reduced order model with an impulse excitation for two excitation amplitudes.forms a sine The third excitation function forms an impulse All experiments are conducted with two current density peak amplitudes, namely, 3000 and AT.
The results for the impulse validation [excitation (44)] are plotted in Fig. 18.The two amplitudes nicely show the nonlinear saturation effect for the magnetic flux density.The reduced model leads to a small deviation in the relaxation phase right after impulse.Also, the effect of the eddy currents is intuitive and well-represented by the reduced model.In all categories, the pECSW method performs better.The global quantities, like the magnetic energy and the force, match nearly perfectly with the reference.In order to generate a weighted single error measure over multiple evaluation points, the area and the circumference length in the radar plot of the reduction error Fig. 19 are computed and displayed in the corresponding legend.
Two statements can be drawn from Fig. 19.First, the reduced order model represents its own training simulation the best.This was to be expected.This is often called onboard validation.Second, the estimation of the force is relatively robust.In general, the simulation results show a changing accuracy when the excitation amplitude is changed, but no clear correlation was found between the excitation level and the accuracy of the reduction strategies.
Conclusion: The cECSW method brings major advantages to the computation of mechanical forces in electromechanical systems.With relatively small offline effort, the model can be reduced significantly.This allows the model to be stored in very limited memory environments and the force to be computed with only very small computational power.

A. Summary
Model order reduction not only reduces the wall time for the engineer but also provides a model with a reduced memory footprint.The family of ECSW methods is helpful for both achievements.After a short recap of the classic ECSW method and its application to magnetodynamic systems, two adaptions were presented.One for the objective of a lean and reliable computation of the magnetic field, including eddy-currents for any unknown excitation.The magnetic field computation for a modified TEAM20 example, once trained, shows good behavior for different excitations and timesteps and moreover is free of the empirical regularization factor.On the other hand, the application of the ECSW idea to speed up the evaluation of the mechanical force allowed us to develop a strategy that can be easily interpreted physically, incurs only a little modification of the analysis while still evaluating the force with excellent accuracy.

B. Future Work
Finding a sound reduction basis is identified as key to a successful reduction.The hyperreduction part is found to be robust and does not degrade the result markedly.Expanding the idea of pECSW could also lead to different ways than POD for a reduction basis.A kind of non-linear condensation is thinkable as pre-processing.This could help in reducing and, therefore, accelerating the full system in advance before calculating training snapshots for the nonlinear domain.

C. Closing
In terms of the ECSW method has been shown to produce accurate results for a wide range of systems.It has the advantage of being relatively simple to implement compared to other model order reduction techniques and is well-suited for large-scale systems where computational resources are limited.Due to the method's general origin, it is adaptable to different objective functions.Here, we have proposed two efficient improvements that were shown for magnetodynamic systems.The first improvement consists of applying hyperreduction only to domains that are considered non-linear.A second improvement is proposed, where hyperreduction ideas are extended to the computation of mechanical forces during the post-processing of electromagnetic analysis.Concluding, both adaptions perform well in the given context.

Fig. 2 .
Fig.2.Generic example of selected ECSW elements and the associated dofs.

Fig. 3 .
Fig. 3. Generic example of selected elements in the non-linear domain.(a) Elements that contribute to edge values of the non-linear material.(b) ECSW selected elements only in the non-linear material.

Fig. 4 .
Fig.4.pECSW-domain decomposition approach.The magnetodynamic equation is different for the linear and the non-linear subdomain.Only the non-linear domain is considered for the hyperreduction process.

Fig. 7 .
Fig. 7. Generic example of configurational training.(a) First training configuration.The rigid, moveable body is marked in red.(b) Second training configuration.Virtually displaced body (exaggerated for visualization).

Fig. 14 .
Fig. 14.TEAM20 solution for the B-field and forces compared to the official measurements.

Fig. 15 .
Fig. 15.Result curves for the three ECSW methods for two numbers of elements considered.The reference is solid blue.

Fig. 18 .
Fig.18.Cross validation of the reduced order model with an impulse excitation for two excitation amplitudes.

Fig. 19 .
Fig.19.Relative error RE (41), evaluated for different quantities (global energy W mag and force F z , local field quantities B and j).The graph contains the three test trajectories.As reference serves their full simulation, respectively.A is the area enclosed by the error points.L is the circumference length.

TABLE II COMPARISON
OF ECSW TIME CONSUMPTION RESULTS.THE OFFLINE TIME DENOTES THE TIME FOR THE SNNLS TO FIND THE SELECTED ELEMENTS AND CORRESPONDING WEIGHTS.THE DURATION TO FIND THE NEXT TIME INTEGRATION STEP IS AVERAGED