Parallel-In-Time Simulation of Eddy Current Problems Using Parareal

TL;DR

This paper presents a novel application of the Parareal method for parallel-in-time simulation of eddy current problems, addressing differential algebraic challenges and demonstrating its effectiveness on a realistic induction machine model.

Contribution

It adapts the Parareal method to eddy current problems with differential algebraic equations, enabling efficient time-parallel simulations.

Findings

Successful adaptation of Parareal for eddy current problems

Automatic incorporation of differential algebraic constraints

First demonstration on a realistic induction machine model

Abstract

In this contribution the usage of the Parareal method is proposed for the time-parallel solution of the eddy current problem. The method is adapted to the particular challenges of the problem that are related to the differential algebraic character due to non-conducting regions. It is shown how the necessary modification can be automatically incorporated by using a suitable time stepping method. The paper closes with a first demonstration of a simulation of a realistic four-pole induction machine model using Parareal.

Figures

Equations

$$\sigma {\partial }_t\vec{A}+\nabla \times (\nu \nabla \times \vec{A})={\vec{J}}_{\mathrm{s}}(t)$$

$$\vec{n}\times \vec{A}{|}_{\Gamma }=0\text{ where }\Gamma =\partial \Omega$$

$${\int }_{\Omega }\vec{w}\cdot \sigma {\partial }_t\vec{A}+\nabla \times \vec{w}\cdot (\nu \nabla \times \vec{A})\mathrm{d}\Omega$$

$$\vec{A}(\vec{x},t)\approx \sum _{i=1}^n{\vec{w}}_i(\vec{x})a_i(t)$$

$$\displaystyle\mathbf{M}_{\sigma}\mathrm{d}_{t}\mathbf{a}(t)+\mathbf{K}_{\mathbf{\nu}}\bigl{(}\mathbf{a}(t),\theta(t)\bigr{)}\mathbf{a}(t)$$

$$\displaystyle\mathrm{d}_{t}\theta(t)=\omega(t)\quad\text{and}\quad I\mathrm{d}_{t}\omega(t)+\kappa\theta(t)=T\bigl{(}\mathbf{a}(t)\bigr{)}$$

$$\displaystyle\mathbf{M}\mathrm{d}_{t}\mathbf{u}(t)+\mathbf{K}\bigl{(}\mathbf{u}(t))\mathbf{u}(t)$$

$${\mathbf{u}}_i=\mathbf{P}{\mathbf{u}}_{i,\sigma }+\mathbf{Q}{\mathbf{u}}_{i,0}.$$

$${\mathbf{Q}}^{\top }{\mathbf{K}}_{\nu }\mathbf{Q}{\mathbf{u}}_{0,0}={\mathbf{Q}}^{\top }{\mathbf{X}}_{\mathrm{s}}{\mathbf{i}}_0-{\mathbf{Q}}^{\top }{\mathbf{K}}_{\nu }\mathbf{P}{\mathbf{u}}_{0,\sigma }.$$

$$\displaystyle\Bigl{(}\frac{1}{\delta t}\mathbf{M}+\mathbf{K}\bigl{(}\mathbf{u}_{i+1}\bigr{)}\Bigr{)}\mathbf{u}_{i+1}=\mathbf{f}_{i+1}+\frac{1}{\delta t}\mathbf{M}\mathbf{u}_{i}$$

$${\mathbf{u}}_i=\mathcal{F}(t_i,t_0,{\mathbf{u}}_0)$$

$$\mathbf{H}(\mathbf{U}):=\{\begin{array}{ll}{\mathbf{U}}_0-{\mathbf{u}}_0& =0,\\ {\mathbf{U}}_1-\mathcal{F}(T_1,T_0,{\mathbf{U}}_0)& =0,\\ & ⋮\\ {\mathbf{U}}_{N-1}-\mathcal{F}(T_{N-2},T_{N-1},{\mathbf{U}}_{N-2})& =0.\end{array}$$

$$\frac{\partial\mathbf{H}}{\partial\mathbf{U}}\bigl{(}\mathbf{U}^{(k-1)}\bigr{)}\;(\mathbf{U}^{(k)}-\mathbf{U}^{(k-1)})=-\mathbf{H}(\mathbf{U}^{(k-1)}).$$

$$\frac{\partial H_{j}}{\partial\mathbf{U}}\bigl{(}\cdot\bigr{)}=\Bigl{[}\mathbf{0},\ldots\mathbf{0},-\frac{\partial\mathcal{F}}{\partial\mathbf{U}}\bigl{(}T_{j},T_{j-1},\cdot\bigr{)},\mathbf{I},\mathbf{0}\ldots,\mathbf{0}\Bigr{]}$$

$${\mathbf{U}}_j^{(k)}$$

$$\displaystyle\qquad\quad+\frac{\partial\mathcal{F}}{\partial\mathbf{U}}\bigl{(}T_{j},T_{j-1},\mathbf{U}_{j-1}^{(k-1)}\bigr{)}\;({\mathbf{U}^{(k)}_{j-1}}-{\mathbf{U}^{(k-1)}_{j-1}}),$$

$$\approx \mathcal{F}(T_j,T_{j-1},{\mathbf{U}}_{j-1}^{(k-1)})$$

$$+\mathcal{G}(T_j,T_{j-1},{\mathbf{U}}_{j-1}^{(k)})-\mathcal{G}(T_j,T_{j-1},{\mathbf{U}}_{j-1}^{(k-1)})$$

$${\bar{\mathbf{u}}}_j=\mathcal{G}(T_j,T_{j-1},{\mathbf{U}}_{j-1})$$

$${\tilde{\mathbf{u}}}_j=\mathcal{F}(T_j,T_{j-1},{\mathbf{U}}_{j-1})$$

$${\mathbf{U}}_j^{(k)}={\tilde{\mathbf{u}}}_j^{(k-1)}+{\bar{\mathbf{u}}}_j^{(k)}-{\bar{\mathbf{u}}}_j^{(k-1)}.$$

$$\mathbf{P}{\mathbf{U}}_j^{(k)}=\mathbf{P}{\tilde{\mathbf{u}}}_j^{(k-1)}+\mathbf{P}{\bar{\mathbf{u}}}_j^{(k)}-\mathbf{P}{\bar{\mathbf{u}}}_j^{(k-1)}$$

Full text

I Introduction

The numerical simulation of the eddy current problem in time domain is computationally expensive due to implicit time stepping. This is particularly challenging if long time periods have to be considered as for example when the start-up of an electrical machine is simulated, possibly with surrounding circuitry [1]. Most implementations, as in GetDP [2], use low-order time stepping schemes as for example the $\theta$-method, [3]. On the other hand, higher order time integration methods, e.g. [4, 5], have been proposed but are rarely used in practice. Recently, explicit methods gained interest as computational hardware architectures seem to favor those algorithms [6, 7, 8]. Another approach is time domain parallelization [9]. For example the reformulation of the time stepping process as one big system of equations has been proposed in [10]. However, it requires to rewrite the time-stepping code.

In this contribution the usage of the (non-intrusive) Parareal method [11, 12] is proposed and its application to a real world electrical engineering problem, i.e., an electrical machine, is shown. Furthermore, the method is adapted to the particular challenges of space discretized eddy current problems related to the differential algebraic character of the equation. The Parareal method has already been applied to wide range of problems in mathematics and in physics. These problems include for example linear and nonlinear parabolic problems [13], molecular dynamics [14], stochastic differential equations [15], Navier Stokes [16], quantum control [17].

This paper is structured as follows: after this introduction, Section II discusses the modeling and the differential algebraic character of the system. Section III introduces the Parareal algorithm and its adaption. Finally, the numerical results are presented in Section IV. Conclusions are given in Section V.

II Modeling and Discretization

When disregarding displacement currents, and introducing the magnetic vector potential $\vec{A}$ as unknown, one obtains the eddy current problem in A⋆-formulation

[TABLE]

on the domain $\Omega\times\mathcal{I}$ and $\mathcal{I}:=(t_{0},t_{\text{end}}]$, see Fig. 1. The problem is well posed when supplying a gauge condition, suitable boundary conditions, e.g. Dirichlet

[TABLE]

and an initial value $\vec{A}(\vec{r},t_{0})=\vec{A}_{0}(\vec{r})$ with $\vec{r}\in\Omega$. The material is described by the conductivity $\sigma$ and nonlinear reluctivity $\nu$; the current density $\vec{J}_{\mathrm{s}}=\sum_{k}\vec{\chi}_{\mathrm{s},k}\,i_{k}$ is given by stranded-conductor winding functions $\vec{\chi}_{\mathrm{s},k}$, which homogeneously distribute the currents $i_{k}$.

II-A Discretization by Finite Elements

Equation (1) is reformulated in the following weak form: find $\vec{A}\in H_{0}(\mathrm{curl},\Omega)$ such that

[TABLE]

for all $\vec{w}\in H_{0}(\mathrm{curl},\Omega)$. Discretization by a finite set of edge elements [18]

[TABLE]

yields for the induction machine ”im_3kw”, cf. Fig. 2, the following system of differential algebraic equations (DAEs)

[TABLE]

where $\mathbf{a}(t)\in\mathbb{R}^{n}$ is the vector of (line-integrated) magnetic vector potentials, $\mathbf{M}_{\sigma}\in\mathbb{R}^{n\times n}$ denotes the (singular) mass matrix representing the conductivities and $\mathbf{j}_{\mathrm{s}}(t)\in\mathbb{R}^{n}$ describes the discretized source current density. Finally, $\mathbf{K}_{\nu}(\mathbf{a},\theta)\in\mathbb{R}^{n\times n}$ is the curl-curl matrix which depends on rotor angle and flux. Movements are considered by the moving band approach [19] determined by the mechanical equation

[TABLE]

with initial values $\theta(t_{0})=\theta_{0}$ and $\omega(t_{0})=\omega_{0}$, where $\omega$ is the angular velocity, $I$ the inertia, $\kappa$ the torsion coefficient and $T$ defines the mechanical excitation given by the magnetic field.

Let us address in the following the current driven coupled problem (2) and (3) as

[TABLE]

with unknown $\mathbf{u}^{\!\top}=[\mathbf{a}^{\!\top},\theta,\omega]$ and the obvious definitions for $\mathbf{M}$, $\mathbf{K}$ and $\mathbf{f}$.

II-B Differential algebraic equations

The solution of the DAE (1) is straightforward since the system is an index-1 DAE [21, 4]. It can be treated with standard techniques, while higher index problems are increasingly more difficult to solve [3].

Let $\mathbf{M}^{+}$ be the Moore-Penrose pseudo inverse of $\mathbf{M}$ such that $\mathbf{P}=\mathbf{M}^{+}\mathbf{M}$ and $\mathbf{Q}=\mathbf{I}-\mathbf{P}$ denote projectors decomposing the vector potential $\mathbf{u}_{i}=\mathbf{u}(t_{i})$ at each time instance into its differential and algebraic components, respectively

[TABLE]

When solving (2) for given currents $\mathbf{i}_{i}=\mathbf{i}(t_{i})$, only initial conditions for the differential components $\mathbf{u}_{0,\sigma}$ may be prescribed. The algebraic components $\mathbf{u}_{0,0}$ must be consistently determined by solving the constraint

[TABLE]

However, when using the implicit Euler method to solve an initial value problem with inconsistent data, i.e, $\mathbf{u}_{0,\sigma}$ and $\mathbf{u}_{0,0}$ do not fulfill (5), a projection is automatically carried out: the time stepping instruction for $t_{i}$ to $t_{i+1}=t_{i}+\delta t$

[TABLE]

ignores inconsistent algebraic components after the first step due to the term $\mathbf{M}\mathbf{u}_{i}=\mathbf{M}\mathbf{P}\mathbf{u}_{i}=\mathbf{M}\mathbf{u}_{i,\sigma}$. This is generally not the case for higher index DAEs and other time-stepping schemes, see e.g. [21].

The implicit Euler method is a numerical implementation of the solution operator $\mathcal{F}:\mathcal{I}\times\mathcal{I}\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ such that

[TABLE]

which propagates $\mathbf{u}_{0}$ through time. Let us define another coarse propagator denoted by $\mathcal{G}:\mathcal{I}\times\mathcal{I}\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ of lower precision, e.g., implicit Euler with a time step $\Delta t\gg\delta t$.

III The Parareal Method

Let us split the total time interval into smaller intervals $\mathcal{I}_{j}:=(T_{j-1},T_{j}]$ with $t_{0}=T_{0}<T_{1}<\ldots<T_{N}=t_{\text{end}}$ according to the number of CPUs $N$ available. On each interval, one defines the equation (4) with initial value $\mathbf{U}_{j-1}:=\mathbf{u}(T_{j-1})$, final value $\mathbf{U}_{j}:=\mathbf{u}(T_{j})$. Continuity at the interfaces $T_{j}$ is established by matching conditions, see Fig. 3

[TABLE]

In other words, the problem of matching can be considered as the unknown of a nonlinear equation $\mathbf{H}:\mathbb{R}^{N\cdot n}\to\mathbb{R}^{N\cdot n}$ in the variable ${\mathbf{U}}^{\!\top}=[\mathbf{U}_{0}^{\!\top},...,\mathbf{U}_{j}^{\!\top},...,\mathbf{U}_{N-1}^{\!\top}]$.

III-A Interpretation as Newton’s Method

The system (7) can be solved by Newton’s method. It reads using the superscript $(k)$ to account for the iterations

[TABLE]

The $j$-th row of the Jacobian matrix $\partial\mathbf{H}/\partial\mathbf{U}$ has only two entries in columns $j-1$ and $j$. It is given by

[TABLE]

where $\mathbf{I}$ denotes the identity of dimension $n$. After rearranging the terms, equation (8) is equivalent to the explicit update formula

[TABLE]

which is a Quasi-Newton method as proposed e.g. in [12].

Due to the splitting of the time interval, one may take advantage of the parallel architecture of modern computers to speed up the time integration similar to multiple shooting methods [11]. The pseudo code of the resulting Parareal algorithm is shown in Alg. 1.

III-B Discussion of the Algorithm

The algorithm solves two kinds of problems in a nested loop until convergence is reached: a cheap problem defined on a coarse time and possibly spatial grid is solved sequentially (line 5, Alg. 1) to propagate missing initial conditions and high-fidelity problems are solved in parallel on the intervals $\mathcal{I}_{j}$ (line 9, Alg. 1). For both problems the implicit Euler method (6) can chosen, or alternatively a higher order method.

The solution of the cheap problem at $T_{j}$ is denoted by

[TABLE]

which is computed by propagating the initial value ${\mathbf{U}}_{j-1}$ from $T_{j-1}$ to $T_{j}$ by coarsely discretizing (2) in time, i.e., using large time steps $\Delta t$. The solution of the high-fidelity problem is

[TABLE]

obtained by solving (2) with initial condition ${\mathbf{U}}_{j-1}$ using fine discretizations, i.e., small $\delta t$. This allows to rewrite the update equation (11) in Alg. 1 (line 6) as

[TABLE]

It can be shown, [12, Theorem 1], that Alg. 1 yields the correct solution until time $T_{k}$ after $k$ iterations, so the correct solution is obtained after at most $N$ iterations. This implies that Parareal does not take longer than the sequential time stepping procedure using the fine solver $\mathcal{F}$ from (13) if neglecting the computational costs of the coarse solution operator $\mathcal{G}$. However, in this case Parareal requires up to $N$-times more CPU time due to its parallel processing.

III-C Discussion of the Algebraic Equations

For the eddy current problem line 6 of Alg. 1 must be adapted to reflect the differential algebraic character, i.e.

[TABLE]

with a subsequent solve of (5) to obtain a consistent $\mathbf{Q}\mathbf{u}^{(k)}_{j}$. However, when using Implicit Euler as shown in (6), this step is automatically taken care of. Similarly, the norm in line 3 should be adapted to only account for differential components, e.g. by considering a projection or the eddy current losses.

IV Numerical Example

Parareal is particularly interesting for problems with multi-tone solutions, e.g. due to slotting as present in the machine ”im_3kw” as visible in Fig. 4. This was also observed in molecular-dynamics [14]. However, speed-ups have also been observed for less favorable problems [12].

The algorithm was implemented in GNU Octave [22]. GetDP is used for the simulation of the 2D model with 8308 degrees of freedom [2, 20]. They are executed on an Intel Xeon cluster with $80\times 2.00$GHz cores, i.e., $8\times$E7-8850 and $1$TB DDR3 memory.

IV-A Sequential Time Steppers

As sequential time stepper the implicit Euler scheme is applied for 10 electrical periods, i.e., $\mathcal{I}=(0,0.2]$ s. The time grid of this simulation is refined from $\delta t=10^{-3}$ s and $\delta t=10^{-4}$ s to $\delta t=10^{-5}$ s. The coarsest sequential simulation takes approx. 15min, while the ones with finer grids correspondingly more time, i.e., 2h and 20h, respectively. Each time step requires ca. $0.3$s for matrix reassembly, Newton and solving the systems.

IV-B Parareal

The Parareal implementation uses OpenMP parallelized calls of GetDP. The implicit Euler method is used with time step sizes $\delta t=10^{-5}$ s and $\Delta t=10^{-3}$ s for the fine and coarse problem, respectively. Fig. 5 shows the errors in comparison with the sequential reference simulation at $\delta t=10^{-5}$ s. For the Parareal simulation $N=40$ cores have been used. After $k=4$ iterations, a relative $l_{2}$ accuracy of $10^{-2}$ has been obtained. The sudden increase of the error after approx. $T_{4}$ can be explained by the convergence of the Parareal algorithm: it is expected to reproduce the sequential time-stepper’s solution for $t<T_{k}$ ($k=4$). The computation has a potential speed-up of $N/k=10$ with respect to the reference simulation when neglecting communication costs and the coarse grid solution. The speed-up is obtained since the effective length of the time interval was reduced by a factor of $N=40$, but iterated $k=4$ times. The parallelization allows to obtain errors below $1$% within the same time that a sequential simulation needs to get errors of $100$%. However, due to suboptimal implementation the actual speed-up was only ca. 3 times.

V Conclusion

This paper proposes the usage of Parareal for eddy current problems and discussed necessary modification due to non-conducting regions in the computational domain. A first implementation shows quick convergence of the iterative algorithm and promises a high speed-up. Future research will investigate the optimal choice of the coarse propagator.

Acknowledgement

This work was supported by DFG grants SCHO 1562/1-1 and CL 143/11-1, the Excellence Initiative of German Federal and State Governments and the Graduate School for CE at TU Darmstadt. The authors would like to thank R.V. Sabariego and N. Marsic for the fruitful discussions on GetDP.