1

Multiple Right-Hand Side Techniques in
Semi-Explicit Time Integration Methods for
Transient Eddy Current Problems

arXiv:1611.06721v2 [] 23 Sep 2017

Jennifer Dutiné∗ , Markus Clemens∗ , and Sebastian Schöps†
∗ Chair of Electromagnetic Theory, University of Wuppertal, 42119 Wuppertal, Germany
† Graduate School CE, Technische Universität Darmstadt, 64293 Darmstadt, Germany

Abstract—The spatially discretized magnetic vector potential
formulation of magnetoquasistatic field problems is transformed
from an infinitely stiff differential algebraic equation system into
a finitely stiff ordinary differential equation (ODE) system by
application of a generalized Schur complement for nonconducting
parts. The ODE can be integrated in time using explicit time
integration schemes, e.g. the explicit Euler method. This requires
the repeated evaluation of a pseudo-inverse of the discrete curlcurl matrix in nonconducting material by the preconditioned conjugate gradient (PCG) method which forms a multiple right-hand
side problem. The subspace projection extrapolation method
and proper orthogonal decomposition are compared for the
computation of suitable start vectors in each time step for the
PCG method which reduce the number of iterations and the
overall computational costs.

∂Ω
Ωc

Ωs
Ωn

Fig. 1. Computational domain Ω split into three regions: conductive and
nonlinearly permeable (Ωc ), nonconductive with constant permeability (Ωn )
and nonconductive with excitation (Ωs ).

Index Terms—Differential equations, eddy currents.

I. I NTRODUCTION

I

N the design process of transformers, electric machines,
etc., simulations of magnetoquasistatic field problems are
an important tool. In particular in multi-query scenarios,
as needed e.g. in the case of uncertainty quantification or
optimization, using efficient and fast algorithms is important.
The spatial discretization of the magnetic vector potential
formulation of eddy current problems yields an infinitely stiff
differential-algebraic equation system of index 1 (DAE). It
can only be integrated in time using implicit time integration
schemes, as e.g. the implicit Euler method, or singly diagonal
implicit Runge-Kutta schemes [1], [2]. Due to the nonlinear
B-H-characteristic of ferromagnetic materials large nonlinear
equation system have to be linearized, e.g. by the NewtonRaphson method, and resolved in every implicit time step. At
least one Newton-Raphson iteration is required per time step.
The Jacobian and the stiffness matrix have to be updated in
every iteration.
A linearization within each time step is avoided if explicit
time integration methods are used. First approaches for this
were published in [3] and [4], where different methods are
used in the conductive and nonconductive regions respectively.
In [3], the Finite Difference Time Domain (FDTD) method
is applied in the conductive regions, while the solution in
the nonconductive regions is computed using the Boundary
Element Method (BEM) [3]. In [4] an explicit time integration
method and the discontinuous Galerkin finite element method
Corresponding author: J. Dutiné (email: dutine@uni-wuppertal.de)

(DG-FEM) are applied in conductive materials, while the finite
element method based on continuous shape functions and an
implicit time integration scheme are used in nonconductive
domains [4]. In another recent approach, a similar DG-FEM
explicit time stepping approach is used for an H − Φ formulation of the magnetoquasistatic field problem [5].
This work is based on an approach originally presented
in [6], where the magnetoquasistatic DAE based on an
~ ∗ −field formulation is transformed into a finitely stiff orA
dinary differential equation (ODE) system by applying a
generalized Schur complement.
The structure of this paper is as follows: Section II introduces the mathematical formulation of the eddy current
problem and the transformation to an ordinary differential
equation. In Section III the time stepping and the resulting
multiple right-hand side problem are discussed. Here, also the
use of the subspace projection extrapolation method and of
the proper orthogonal decomposition method as multiple righthand side techniques is described. In Section IV the simulation
results for validating the presented approach and the effect of
subspace projection extrapolation method and of the proper
orthogonal decomposition method on a nonlinear test problem
are presented. The main results of this paper are summarized
in Section V.
II. M ATHEMATICAL F ORMULATION
~ ∗ −formulation is given
The eddy current problem in the A
by the partial differential equation
 


~
∂A
~ ∇×A
~ = ~JS ,
κ
+∇× ν ∇×A
(1)
∂t

c 2017 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of
any copyrighted component of this work in other works.

2

~ is the time-dependent
where κ is the electrical conductivity, A
magnetic vector potential, ν is the reluctivity that can be
~ S iS (t), where
nonlinear in ferromagnetic materials and ~JS = X
~ S distributes
iS (t) is the time-dependent source current and X
the current density spatially. Furthermore, initial values and
boundary conditions are needed. The weak formulation of (1)
~
leads to the variational problem: find A
Z
~
∂A
w
~ ·κ
dΩ +
∂t
Ω
Z
Z


~ ∇×A
~ dΩ =
∇×w
~ ·ν ∇×A
w
~ · J~s dΩ
Ω

Ω

for all w
~ ∈ H0 (curl, Ω) where we denote the spatial domain
with Ω and assume Dirichlet conditions at the boundary
∂Ω, see Fig. 1. Discretization and choosing test and ansatz
functions from the same space according to Ritz-Galerkin, i.e.,
~ r, t) ≈
A(~

N
dof
X

w
~ i (~r)ai (t)

where Mc is the conductivity matrix, Kc is the nonlinear
curl-curl reluctivity matrix in conducting regions, Kn is the
typically constant curl-curl matrix in nonconducting regions,
Kcn is a coupling matrix, and jSn is the source current
typically defined in the nonconducting domain only. The
conductivity matrix in (3) is not invertible and therefore the
problems consists of differential-algebraic equations (DAEs).
The numerical solution of these systems is more difficult than
in the case of ordinary differential equations (ODEs). The level
of difficulty is measured by the DAE index, which can be
roughly interpreted as the number of differentiations needed to
obtain an ODE from the DAE [1]. System (3) is essentially an
index-1 DAE with the speciality that the algebraic constraint,
i.e., the second equation in (3), is formally not uniquely
solvable for an without defining a gauge condition due to the
nullspace of the discrete curl-curl operator Kn . However, it is
well known that many iterative solvers have a weak gauging
property, e.g. [7], such that a formal regularization can be
avoided.
Relying on this weak gauging property, the generalized
Schur complement
(4)

where K+
n represents a pseudo-inverse of Kn in matrix form,
is applied to (3) and transforms the DAE into
d
ac + KS (ac )ac
dt
an

= −Kcn K+
n js,n ,

(7)

are computed in the m-th time step, where ∆t is the time step
size. The Courant-Friedrich-Levy (CFL) criterion determines
the maximum stable time step size of explicit time integration
methods [1]. For the explicit Euler method
∆t ≤

2

λmax M−1
c KS (ac )

(8)

is an estimation for the maximum stable time step size,
where λmax is the maximum eigenvalue [8]. The maximum
eigenvalue can be estimated using the power method [10].
III. M ULTIPLE R IGHT-H AND S IDE P ROBLEM

leads to a spatially discretized symmetric equation system in
time domain. Separation of the degrees of freedom (dofs) into
two vectors ac storing the dofs allocated in conducting regions
(if ~r ∈ Ωc ) and an holding the dofs allocated in nonconducting
regions (if ~r ∈ Ωn ∪ Ωs ) yields the DAE system
   

   
Kc (ac ) Kcn ac
0
Mc 0 d ac
+
=
, (3)
jSn
0 0 dt an
Kn an
K>
cn

Mc

+ m
+ > m
am
n = Kn js,n − Kn Kcn ac

(2)

i=0

>
KS (ac ) := Kc (ac ) − Kcn K+
n Kcn ,

evaluated using the preconditioned conjugate gradient method
(PCG) [9]. The finitely stiff ODE (5a) can be integrated
explicitly in time, e.g. by using the explicit Euler method.
Using this time integration method, the expressions


m−1
m
m−1 m−1
am
+∆tM−1
Kcn K+
)ac
, (6)
c = ac
c
n js,n −KS (ac

(5a)

+ >
= K+
n js,n − Kn Kcn ac . (5b)

A regularization of Kn by a grad-div or tree/cotree gauging
can be used alternatively [6], [8]. Here, the pseudo-inverse is

As the matrix Kn remains constant within each explicit time
step, the repeated evaluation of a pseudo-inverse K+
n in (6),
(7) forms a multiple right-hand side (mrhs) problem of the
form
Kn ap = jp ⇔ ap = K+
(9)
n jp .
Here, jp represents one of the right-hand side vectors
> m−1
> m
. Instead of computing the matrix
jm
s,n , Kcn ac , and Kcn ac
+
Kn explicitly, a vector ap is computed according to (9) using
the preconditioned conjugate gradient (PCG) method [9].
Improved start vectors for the PCG method can be obtained
by the subspace projection extrapolation (SPE) method or
the proper orthogonal decomposition (POD) method. In the
SPE method, the linearly independent column vectors of a
matrix USPE are formed by a linear combination of an
orthonormalized basis of the subspace spanned by solutions
ap from previous time steps. The modified Gram-Schmidt
method is used for this orthonormalization procedure [11]. The
improved start vector x0,SPE is then computed by [12]
−1 >
x0,SPE := USPE U>
USPE K>
(10)
SPE Kn USPE
cn jp .
As only the last column vector in the matrix USPE changes
in every time step, all other matrix-column-vector products in
computing Kn USPE in (10) are reused from previous time
steps in a modification of the procedure in [12] referred to as
the ”Cascaded SPE” (CSPE) [9].
When using the POD method for the PCG start vector
generation, NPOD solution vectors from previous time steps
form the column vectors of a snapshot matrix X which is
decomposed into
X = UΣV>
(11)
using the singular value decomposition (SVD) [13], [14],
[15]. Here, U and V are orthonormal matrices and Σ is a
diagonal matrix of the singular values ordered by magnitude
(σi ≥ σj for i < j). The index k is chosen such that the information of the largest singular values is kept
σk
≤ εPOD .
(12)
σ1

3

The threshold value εPOD is here chosen as εPOD := 10−4 .
A measure how much information is kept can be computed
by the relative information criterion [16]
k
P

σi

i=1
NP
POD

!

≈ 1.

(13)

σi

i=1

Defining UPOD = [U:,1 , ... , U:,k ] as the first k columns of U
allows to compute an improved start vector x0,POD by

−1 >
x0,POD := UPOD U>
UPOD K>
POD Kn UPOD
cn jp . (14)
The repeated evaluation of M−1
in (6) also forms a mrhs
c
problem, and both the POD and the CSPE method can be used
for computing improved start vectors for the PCG method.
In the case of small matrix dimensions of the regular matrix
Mc , the inverse can also be computed directly using GPUacceleration.
IV. N UMERICAL VALIDATION
The ferromagnetic TEAM 10 benchmark problem is used
for numerical validation of the presented explicit time integration scheme for magnetoquasistatic fields [17]. The domain
consists of two square-bracket-shaped steel plates opposite
of each other and a rectangular steel plate between them,
resulting in two 0.5 mm wide air gaps. The model geometry is shown in Fig. 2. The position where the magnetic
field is evaluated is marked as S1. The excitation current
iS = (1 − exp(−t/τ )), where τ = 0.5 s, is applied for a
time interval of 120 ms starting at t = 0 s [17]. The resulting
magnetic flux density is computed for this time interval.
The finite element method (FEM) using 1st order edge
elements is used for the spatial discretization [18]. All simulations are computed on a workstation with an Intel Xeon E5
processor and an NVIDIA TESLA K80 GPU. The conjugate
gradient method is preconditioned by an algebraic multigrid
method [19]. The matrix Mc is inverted using the Magmalibrary and GPU-acceleration [20].
A fine mesh resulting in about 700,000 dofs and the implicit
Euler method are used to validate the simulation code. A good
agreement between the measured results published in [17] and
the simulation of this fine spatial discretization is shown in Fig.
2. The required simulation time of this simulation using the
implicit Euler method is 5.38 days using an in-house implicit
time integration magnetoquasistatic code.
For benchmarking the proposed mrhs techniques for the
(semi-)explicit time integration scheme, a model with a coarse
spatial discretization yielding about 30,000 dofs and the explicit Euler method is used. For this spatial discretization, the
resulting maximum stable time step size according to (8) is
∆tCFL = 1.2 µs. Both meshes are presented in Fig. 3 The
results for the average magnetic flux density are compared with
the results obtained using the same discretization in space and
the implicit Euler method for time integration and show good
agreement, depicted in Fig. 2. The resulting field plots for both
spatial discretizations are shown in Fig. 4. The simulation time
for the implicit time integration method is still 2.58 h.

Fig. 2. Comparison of results for the average magnetic flux density evaluated
at position S1 and model geometry as inset.

TABLE I
S IMULATION T IME AND AVERAGE N UMBER OF PCG I TERATIONS U SING
D IFFERENT S TART V ECTORS x0
Start vector

Avg. Number of PCG Iterations

x0 := am−1
p

3.16

Simulation Time
2.35 h

x0 := x0,POD

2.18

17.35 h

x0 := x0,CSPE

1.02

1.62 h

The effect of computing improved start vectors using POD
or CSPE on the average number of PCG iterations and on
the simulation time is compared to using the solution from
as start vector for the PCG
the previous time step am−1
p
method. An overview is presented in Table I and shows that
both the CSPE and the POD start vector generation methods
significantly reduce the number of PCG iterations. When using
CSPE the number of column vectors in the operator USPE in
(10) is increased during the simulation to improve the spectral
information content of USPE . This number remains below 20.
Thus, only small systems have to be solved for the inversion
in (10) and the effort to perform all computations of the
CSPE method is low. This is also confirmed by the simulation
time which is shortest when using CSPE. The simulation time
resulting from using explicit time integration and CSPE for
start vector generation is 63 % of the simulation time of the
implicit reference simulation. A bar plot showing the reduced
simulation time by using the explicit Euler scheme and CSPE
compared to using the standard formulation and the implicit
Euler method for time integration is depicted in Fig. 5.
In case of the POD, the amount of information kept
according to (13) is > 0.99 during the entire simulation.
However, the computational effort for performing the SVD and
the computations in (14) is higher than the effort for CSPE.
Although the number of PCG iterations is further decreased,
the simulation time resulting from using POD for start vector
generation is higher than when using apm−1 as start vector for
the PCG method due to the costs of the repeated SVD.

4

using the POD and the CSPE method. Although both reduce
the number of PCG iterations needed, the computational
effort of the CSPE is significantly lower than for the POD
method. Reducing the computational effort of the POD, e.g. by
accelerating the computation of the SVD is subject to further
investigations. Using the CSPE method, the overall simulation
time was reduced by 37 % compared to the simulation time of
the implicit reference simulation.

Fig. 3. Meshes resulting in about 700,000 dofs (left) and in about 30,000
dofs (right).

ACKNOWLEDGMENT
This work was supported by the Deutsche Forschungsgemeinschaft (DFG) under grant numbers CL143/11-1 and
SCHO1562/1-1. The third author is supported by the Excellence Initiative of the German Federal and State Governments
and The Graduate School of CE at TU Darmstadt.
R EFERENCES

| B |/ T
0.00

0.35

0.70

1.05

1.50

~ for the spatial discretization
Fig. 4. Field plots of the magnetic flux density B
with about 700,000 dofs (left) and with about 30,000 dofs (right).

Fig. 5. Comparison of simulation times.

V. C ONCLUSION
The magnetic vector potential formulation of eddy current
problems was transformed into an ODE system of finite
stiffness using a generalized Schur complement. The resulting
ODE system was integrated in time using the explicit Euler
method. A pseudo-inverse of the singular curl-curl matrix in
nonconducting material was evaluated using the PCG method.
Improved start vectors for the PCG method were calculated

[1] E. Hairer and G. Wanner, Solving ordinary differential equations II: stiff
and differential-algebraic problems, 2nd rev. ed., Springer, Berlin, 1996.
[2] A. Nicolet, ”Implicit Runge-Kutta methods for transient magnetic field
computations,” IEEE Trans. Magn., vol. 32, pp. 1405–1408, May 1996
[3] T. V. Yioutlsis, K. S. Charitou, C. S. Antonopoulos, and T. D. Tsiboukis,
”A finite difference time domain scheme for transient eddy current
problems,” IEEE Trans. Mag., vol. 37, pp. 3145–3149, Sep. 2001.
[4] S. Außerhofer, O. Bı́ro, and K. Preis, ”Discontinuous Galerkin finite
elements in time domain eddy-current problems,” IEEE Trans. Magn.,
vol. 45, pp. 1300–1303, Feb. 2009.
[5] J. Smajic, ”DG-FEM for time domain H − Φ eddy-current analysis,”
presented at the 17th biennial Conference on Electromagnetic Field
Computation (CEFC) 2016, Nov. 2016.
[6] M. Clemens, S. Schöps, H. De Gersem, and A. Bartel, ”Decomposition
and regularization of nonlinear anisotropic curl-curl DAEs,” Compel, vol.
30, pp. 1701–1714, 2011.
[7] M. Clemens and T. Weiland, ”Transient eddy-current calculation with the
FI-method,” IEEE Trans. Magn., vol. 35, pp. 1163-1166, May 1999.
[8] S. Schöps, A. Bartel, and M. Clemens, ”Higher order half-explicit time
integration of eddy current problems,” IEEE Trans. Magn., vol. 48, pp.
623-626, Feb. 2012.
[9] J. Dutiné, M. Clemens, S. Schöps, and G. Wimmer, ”Explicit time
integration of transient eddy current problems,” presented at the 10th
International Symposium on Electric and Magnetic Fields (EMF) 2016.
Full paper to appear in J. Num. Mod.: ENDF.
[10] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., The
Johns Hopkins University Press, Baltimore, 2013.
[11] L. N. Trefethen and D. Bau, Numerical linear algebra, Society for
Industrial and Applied Mathematics, Philadelphia, 1997.
[12] M. Clemens, M. Wilke, R. Schuhmann, and T. Weiland, ”Subspace
projection extrapolation scheme for transient field simulations,” IEEE
Trans. Magn., vol. 40, pp. 934–937, Apr. 2004.
[13] T. Henneron and S. Clénet, ”Model order reduction of non-linear
magnetostatic problems based on POD and DEI methods,” IEEE Trans.
Magn., vol. 50, pp. 33–36, Feb. 2014.
[14] T. Henneron and S. Clénet, ”Model-order reduction of multiple-input
non-linear systems based on POD and DEI methods,” IEEE Trans. Magn.,
vol. 51, pp. 1–4, Mar. 2015.
[15] Y. Sato, M. Clemens, and H. Igarashi, ”Adaptive subdomain model
order reduction with discrete empirical interpolation method for nonlinear
magneto-quasi-static problems,” IEEE Trans. Magn., vol. 52, pp. 1–4,
Mar. 2016.
[16] D. Schmidthäusler, S. Schöps, and M. Clemens, ”Linear subspace reduction for quasistatic field simulations to accelerate repeated computations,”
IEEE Trans. Magn., vol. 50, pp. 421–424, Feb. 2014.
[17] T. Nakata and K. Fujiwara, ”Results for benchmark problem 10 (steel
plates around a coil),” Compel, vol. 9, pp. 181–192, 1990.
[18] A. Kameari, ”Calculation of transient 3d eddy current using edgeelements,” IEEE Trans. Magn., vol. 26, pp. 466–469, Mar. 1990.
[19] K. Stüben, Algebraic multigrid (AMG): an introduction with applications, 1999.
[20] W. Bosma, J. Cannon, and C. Playoust, ”The Magma algebra system
I: the user language,” J. Symbolic Comput., vol. 24, pp. 235–265, Sep.
1997.

