High Order Relaxation Methods for Co-simulation of Finite Element and Circuit Solvers

Coupled problems result in very stiff problems whose characteristic parameters differ with several orders in magnitude. For such complex problems, solving them monolithically becomes prohibitive. Since nowadays there are optimized solvers for particular problems, solving uncoupled problems becomes easy since each can be solved independently with its dedicated optimized tools. Therefore the co-simulation of the sub-problems solvers is encouraged. The design of the transmission coupling conditions between solvers plays a fundamental role. The current paper applies the waveform relaxation methods for co-simulation of the finite element and circuit solvers, investigate the contribution of higher order integration methods and also considers the the theoretical modelling of a boost converter. The method is illustrated on a coupled finite element inductor and a boost converter and focuses on the comparison of the transmission coupling conditions based on the waveform iteration numbers between the two sub-solvers. We demonstrate that for lightly coupled systems the dynamic iterations between the sub-solvers depends much on the internal integrators in individual sub-solvers whereas for tightly coupled systems it depends also to the kind of transmission coupling conditions.


Introduction
Coupled problems are modeled by very complex mathematical models with rate of change of some parameters differing with several orders of magnitude, in both the smallest and largest range of variations. They constitute very stiff problems. Their monolithic solution is computationally expensive. The solution of individual problems by optimized solvers is encouraged.
Waveform relaxation (WR) is domain decomposition in time. After sub-problems are formed, each one is modelled by a system of differential equation indexed by an iteration index which allows then solving the sub-problems independently. This decomposition results in various transmission coupling conditions (TCCs) at interface of the sub-systems. The unkown TCCs are first approximated and then updated until a convergence level is satisfied. The waveform iteration numbers for given TCCs differ with several orders.
The main objective study of the current paper is the comparison of convergence of the WR with respect to the kind of the TCCs and their application to co-simulation more complex electromagnetic systems.
Co-simulation stands for coupling two independent solvers; that is a monolithic problem is subdivided into subproblems so that each of the sub-problem is solved by its own dedicated solver. Whereby an iterative process is created in which an artificial interface is first applied on one of the solver and vice-versa until convergence occurs. The coupling or co-simulation of the finite element solver such as Gmsh-GetDp [5,6,7] to an optimized circuit solver such as LTspice will reveal a drastically reduction of the number of the waveform relaxation numbers.
As was investigated before [1,2,3], the voltageimpedance transmission coupling conditions (VZ-TCCs) converge very fast than the current-voltage conditions (VI-TCCs). The current-voltage transmission coupling conditions can be improved by introducing more sophisticated optimized transmission coupling conditions (OTCCs). The first transmission coupling conditions for circuit coupling were introduced first in [4] and recently for coupling fields and circuits in [1]. Both the technique can be improved if one inspects the integrators in each individual sub-solver. In electronic circuits, some of the elements are static so that optimization becomes difficult, therefore artificial parameters are introduced. Some of the parameters are interpreted as circuit elements and others are dimensionless and regarded as circuit gains. These combinations result in similar operations as the Thevenin equivalent circuit. Each of the parameter represents the impedance of the neighbouring subsystem. The convergence of the Voltage-current and zero-order OTCCs for linear problems are compared. Their convergence test is the foundation explanation of the behaviour of the lucky in convergence of the current-voltage transmission coupling conditions. The current paper is organized into five sections. In addition to the current introduction, Section 2 presents the mathematical modeling of the WR, Section 3 presents the WR applied to coupled field/circuit problem, Section 4 presents application cases, Section 5 presents simulation results. Conclusions and future perspectives are presented in the last section.

Mathematical modelling of WR method
Coupled problems are always modeled by a general differential equation of the form [8]: Classically the differentiability condition suggests that x(t) 2 C 1 ([0, 1)) and u(t) be piecewise continuous, and by using the implicit function theorem [9], if the condition @G @ẋ is non-singular, then equation (1) is algebraically solvable forẋ, and hence transformed into an explicit initial value problem (IVP) formulated as follows: If instead the matrix @G @ẋ is singular, then equation (1) can not be solved directly, operations of transformation using projectors have to be performed. If so, then a decoupled system of IVP coupled to an algebraic equation with the initial condition is then found as follows: where y 2 R m is the differential variable, z 2 R n is the algebraic variable, and f : R m ⇥ R n ⇥ R l ⇥ I ! R m ⇥ I; h:R m ⇥ R n ⇥ R l ⇥ I ! R n ⇥ I .

WR method algorithm
Various WR methods can be described. The algorithm 1 describes the Gauss-Seidel WR applied on the general coupling equation (3).
Considering r coupled problems, and introducing the WR iteration index k, and denotinḡ . . , z r }, the Gauss Seidel WR method is shown in the algorithm 1.
where t 0 is the initial simulation time instant, and t f the final simulation time instant. 2: set j = 0, k = 0, and for t 2 [T 0 , T 1 ] guess initial waveforms y 0 i (t) = x 0i (t) such that y 0 i (t 0 ) = x(t 0 ) and then solve 0 = h(y 0 i (t), z 0 i (t), u i (t)) for the algebraic initial waveform z 0 i (t), 3: for i 2 {1, 2, ..., r} , solve: i go 3:. Otherwise go 5 : , t = 0, and then go 3 : . Otherwise go to 6: 6: output y k+1 i (t), z k+1 i (t), 8i 2 {i, . . . , r} end 3. WR method applied to field/circuit problem The details of electromagnetic properties can be modelled by the Maxwell's equations [10,11]. Those accounts for the distribution of electromagnetic field in the magnetic devices. They provide a link between the magnetic field intensity h, the electric field intensity e, the magnetic field induction density b and the electric field density d to their sources which are the magnetization m, electric polarization p, the electric charge density ⇢ and the electric current density j. The Maxwell's equations together with the constitutive equations writes: curl e @b @t = 0, where µ 0 is the permeability of the vacuum, and ✏ 0 is the permittivity of the vacuum. The fields are assumed to be space coordinate and time dependent i.e. for example to say h = h(x, y, z; t).
We consider the quasi-static approximation. When considering the magnetoquasistatic i.e max x2R @d @t << max x2R |j| , and an a V formulation i.e by the transformations b = rot a and the consequence of it e = @a @t grad V , the equation (4) becomes rot(µ 1 curl a) + @a @t where is the conductivity and j s is the exterior fixed current density. The equation (10) is an a V strong form of the magnetoquasistatic approximation. A finite element formulation in terms of the vector magnetic and scalar electric potentials quite similar to the procedure in the work [5] is found by introducing a test function a' such that where ⌦ is the magnetic device domain and its surrounding, ⌦ C is the conducting subdomain of the magnetic device and ⌦ s is the subdomain where the current is imposed. The equation (11) can be descretized for example by using the Whitney elements: . The GetDp-Gmesh finite element solver such as described in [5,7] provide a mathematical like user friendly similar to the mathematical formulation as in equation (11).
The circuits equations are derived using the modified nodal analysis technique. In order to find connections to the circuit, the voltage, current or the impedance including the resistance and the inductance are found from the finite element post-processing and imposed as the interface waveforms between the circuit and the magnetic device. In that case the coupled field-circuit system reads to the following system of equations: where, with x = [u, i L , i V ] gathers the coefficients of the electric potentials as well as the inductor currents and voltage source currents. Matrices , 0 , and F are built classically using the modified nodal analysis technique. In this paper, the modified nodal analysis technique is used to derive equations for a boost-converter (see Figure 3) coupled to a finite element inductor. Various incidence sub-matrices of equation (12) are given as follows: 1 RL , where R D is the diode variable resistance which is defined as follows: where U D is the diode branch voltage. R T is the transistor variable resistance which is defined as follows: where U T is the transistor branch voltage and P s is the transistors triggering switching periodic function, it is a function of the duty cycle. The coupling of systems is determined by the kind of the TCCs between them. In the case of coupled field/circuit problem, the TCCs can be classified into four [4]: The voltage-impedance TCCs (VZ-TCCs), the voltage-current TCCs (VI-TCCs), the zero-order optimized voltage-current TCCs (VI-Opt-0 TCCs), and the first-order optimized voltage current TCCs (VI-Opt-1 TCCs).
In this case we only recall here two kind of the TCCs. Figure 1 shows the VI and VI-Opt-0 TCCs. The VZ and Figure 2 shows VI-Opt-1 TCCs can be found in [1].

Case of the voltage-current TCCs
In this case, the artificial signals at interface between the field and circuit models are the current on the side of the circuit, and the voltage on the side of the field (see Figure 1 b). These TCCs can be presented as follows: We consider the case of linear coupling functions, the coupling model is given as follows (18) The convergence rate of the VI-TCCs is determined by the spectral radius of the iteration matrix. In order to find the spectral radius we use an integral transform such as the L,aplace transform. This will only give ways of finding the coupling parameters before the system resolution. The whole solution of the system is done only in the time domain. Hence using the Laplace transform on the system of equations (17) and (18), the iteration equation is found as follows:ŷ The spectral radius of ⇧ V I is found as follows: The iteration equation shows that the case of VI-TCCs has possibility to converges at least two iterations: 2  k  M.

Case of zero-order optimized voltage-current TCCs
This section attempts to find a technique of optimizing the TCCs, a technique that brings more freedom of choices on the TCCs. Though the previous case converges in few iterations, there is a path of optimizing for more reduction of the iteration numbers. In this case, the artificial signals at interface between the field and circuit models are provided in such way that the draw back between the decoupled models is considered by coupling parameters that play the role of impedances of sub-problems. The convergence can be fast if the impedance involves artificial signals derivatives, which contribute on the rank deficiency annihilation. By introducing the coupling parameters ↵ and , the zero-order optimized TCCs are provided as follows: and then the WR equations for field/circuit problem are presented as follows: The convergence of the VI-Opt-0-TCCs is determined by the spectral radius of the iteration matrix. Applying the Laplace transform on equations (21) and (22) yields: Solving the systems (23-24) for X 01 and X 02 yields: and putting equation (25) into equations (24), and after rearranging the similar terms yields the following iteration scheme,ŷ where, and the optimal spectral radius is found by a mini-max optimization technique [4] as follows: where s = i! and ! is the circular frequency. After finding the optimization parameters ↵ and they are after are inserted in the time domain model for complete resolution.

Application cases 4.1. Open loop boost-converter
A boost-converter is a step up converter; it rises up a DC voltage to a higher level. Its circuit is composed of a very fast switching transistor and a diode which commutate respectively at proportion fraction D called a duty cycle and 1 D of a constant time period T s . These commutations generate a square voltage waveform which is regulated by an output capacitor so that the load gets a DC voltage with possibly a small variation or a voltage ripple V . Figure  3 shows the topology of the open loop boost-converter. As far as the transistor closes, the diode is reverse biased, and then the capacitor supplies the load while the input inductor is storing energy. Also as the transistor opens, the diode conducts and then, the source and the inductor supply the load which then results in a boosted voltage. A boostconverter can work either in a continuous current mode (CCM), which happens as the inductor current stays only positive; or in a discrete current mode (DCM), where the inductor current can reach very small even null values at some instants in time.
Assuming the inductor current always positive, and considering the volt-second balance [13], the inductor average voltage can be given as follows: , = 1 D, U in the applied power voltage, I in the input current, U O the output voltage, R ON the conducting diode resistance, U D the diode junction threshold voltage. Like wise, considering a capacitor charge balance [13], the capacitor average current is found as follows: , From solving equation (31) and (32), the output voltage of a boost-converter is found as follows: For ideal switching elements, a more simple expression is found as follows: From this expression, it is clear that the output voltage for a boost-converter is a function of the duty cycle and the load resistance. Assuming ideal elements, the difference equation on the inductor is found from and the difference equation on the capacitor is found from The inductor current and output voltage ripples are calculated from the inductor current and output voltage difference equations (36) and (37), and in the case of the transistor is on, those can be found as follows: Assuming ideal elements, the equations (34)-(35) make it possible to fulfil the output required voltage. By predefining the voltage and current ripples, the boost circuit designing parameters are found from equations (40) and (41): In the case of the DCM, the average inductor current is By imposing a constant volt-second requirement on the inductor, one has: where the factor D 2 is defined by By putting equation (44) into equation (44) one finds: hence by solving the quadratic equation (44), the outputinput voltage relation is found as follows: where M = 2L RTs . In that case, the designing parameters for the boost-converter in a DCM are chosen as follows where I O(crit) is the required minimum output critical current. Further details on the theory of converters can be found in [12,13,14,15]. WR method and associated TCCs are applied on a boost converter with a magnetic inductor as its load. Figure 4 presents a coupled system of a magnetic inductor and a boost converter [1]. Several reasons to consider this case include the need of a case much closer to real life applications where particular features are considered like the multiscale in time, the stiffness due to hard switching which causes abrupt change in the impedance of the subsystems. That helps to differentiate among the TCCs those wich withstand for co-simulation of the more tightly coupled subsysytems.

Coupled open loop boost-converter and FE airgapped inductor
A coupled boost-converter and an FE air-gaped inductor are considered. The transistor of the boost-converter is driven by a pulse width modulation (PWM) signal. The air-gaped FE inductor is powered by a quasi-constant voltage, so that the time constants of the two sub-systems are not the same. Obviously the smallest time step should be a small fraction of the time constant of the boost circuit. This scenario constitute a typical multiscale problem. An efficient computation of such problem should be that one which is able to localize all local events of the multi-scale structure and adapt system integration accordingly. Special solvers such as LTspice are dedicated only for circuit simulations and others such as Gmsh-GetDp are dedicated the finite element. As it was stated before, not all events can be resolved in a single environment. The solution of the FE at the same stepping as the circuit can be time consuming even impossible in some cases. Hence a procedure of domain decomposition and cosimulation of circuit and field solvers is here encouraged.
As shown on Figure 3, the boost-converter load is an FE air-gapped inductor. The geometry of the air-gapped inductor is shown on the Figure 4. Table 1 shows the material and geometric dimension of the FE air-gapped inductor considered, where µ 0 = 4⇡e 7 is the permeability of the air.

Simulation Results
For a fixed time window of 24 times the switching period of the transistor, simulations of the circuit are done by the LTspice software whereas the field of finite element inductor is simulated using the GetDp-Gmsh. For transferring data from one software to another an interpolation process is first conducted to adjust the stepsize of the waveforms at interface. The built-in piece wise linear (pwl) function  Table 1: FE numerical data for inductor geometric dimension in a planar symmetry, and the magnetic parameters.   helps for interpolation in LTspice and the built-in Interpola-tionLinear function do the same in GetDp-Gmsh. As stated before, the previous research on the coupling transmission conditions pushed us also on the consideration to provide meaning why the number of iterations were increasing on the case of VI-TCCs. Co-simulation starts with interface artificial physics values or wave forms being interpolated and transferred to the field part or vice versa to circuit and that requires several iterations before convergence occur.
The Figure 5 shows the voltage waveforms at interface between the circuit and the field. This shows that the voltage has been boosted from the input voltage of 100V to output voltage of 200V. The discrepancies in the two waveform are due to that the LTspice is more appropriate for simulation of the circuit since all the phenomena in the circuit such as the transistor and the diode are treated accordingly than in the software with very low order integration methods. Figure 6 shows the step size taken by the two sub-solvers. It is shown that the circuit is adaptive whereas the field only needs a very large step size. The multiscale is clearly shown here and what is more important is that such kind of tiny time steps is not required for the finite element sub-problem which is not concerned by the status of the circuit. It is observed that the time steps used by the circuit solver depicts very well the stiffness of the problem where the circuit takes two extremes of the time constants. Noting also that the time integrator just update the time step from small and takes a possible large time step accordingly. Figure 7 shows the voltage of the diode branch and Figure  8 presents a snapshot of the diode voltage and highlights different stage the time integration of the circuit is taking. There are several occasions that occur; the case when the diode and transistor commutates individually or together.
The Figure 9 shows the convergence errors for various TCCs considered. It is revealed here that the number of WR iterations depends much on the kind of the time integrator of the sub-solver. For this moderately coupled problem it takes 8 iterations for the VI-TCCs with an adaptive integrator, 9 iterations for the VI-Opt-0-TCCs and 14 for the    VI-TCCs with non-adaptive integrator. Though for lightly coupled problems, the difference in terms of WR iterations is not relatively high for the coupling with adaptive integrators and the optimized conditions with non adaptive integrators, the difference is much higher for tightly coupled problems as it is going to be shown in the next paragraph.
In contrast to the case for which the circuit and the field are tightly coupled things will behave in another way; The figure 10 shows the voltage waveform at interface between the circuit and the field. The features in it shows how much the circuit is interacting with the field. Actuary as the transistor switches on and off, the inductance also follows by trying to opposite the change in the voltage.
The Figure 11 shows the convergence errors for various TCCs considered in this case of the more tightly coupled circuit and the field. It is revealed here that the number of WR iterations depends much on the kind of the time integrator of the sub-solver. For the two kind of waveform relaxation method which include the Gauss-Jacobi WR and the Gauss-Seidel, it is shown that for this tightly coupled problem, it takes 54 WR iterations for the classical Guss-Jacobi VI-TCCs i.e that one using low order time integrator whereas it takes 42 WR iterations for the classical Guss-Jacobi VI-TCCs using adaptive time integrator, Also it takes 48 WR iterations for the classical Guss-Jacobi using optimized VI-TCCs and finally 30 WR iterations for the classical Guss-Seidel VI-TCCs with non-adaptive time integrator. It is now clear that the number of the waveform relaxation is drastically reduced when each of the problem is solved according to its time constants characteristics and also when the optimized transmission coupling conditions are considered..

Conclusion and future perspective
The paper has focused on the application of the WR methods for co-simulation of field and circuit solvers with emphasis on two kind of the TCCs. In addition to the VI-TCCs and the VI-Opt-0-TCCs, the VI-TCCs with the circuit solved by more optimized solver such as LTspice proves to be more efficient than other types of the VI-TCCs. WR applied on real cases face convergence problems; we enumerated those difficulties and suggest how to handle convergence problems. Using the more optimized solvers revealed the more reduction of iterations compared to less optimized solvers.
The subtle problem was an increase in the number of the WR relaxation when the VI-TCCs are applied. The fact that those depend on the stiffness of the problem is tackled here by using high order adaptive time integrators and the more accurate is found by using the commercial solver LTspice. The optimized TCCs as introduced in [1,4] were shown to be improved by including the adaptive time integrators.
The major contribution of this paper is a clear explanation on how the convergence of the classical VI-TCCs and the optimized VI-TCCs can be improved. The observations of the results as in the Figure 11 suggest that the co-simulation using specific solvers for each sub-problem is the more advisable solution to coupled problems.
The future continuation of the work will focus on comparison of the TCCs on more nonlinear problems such as considering the saturation and the hysteresis on the field sub-problem and also including feedback in the control circuit will increase more complex features and decision on the use of the TCCs.