(PHD) Schops
(PHD) Schops
(PHD) Schops
Multirate Time-Integration of
Field/Circuit Coupled Problems
Dissertation
und/and
von/by
May, 2011
Die Dissertation kann wie folgt zitiert werden:
urn:nbn:de:hbz:468-20110718-110524-8
[http://nbn-resolving.de/urn/resolver.pl?urn=urn%3Anbn%3Ade%3Ahbz%3A468-20110718-110524-8]
Preface
This treatise was written during my employment at the Chair of Angewandte Mathematik
und Numerische Analysis of the Bergische Universität Wuppertal and my post-graduate
research at the Wave Propagation and Signal Processing Research Group of the Katholieke
Universiteit Leuven funded by a Jahresstipendium für Doktoranden from the DAAD (Ger-
man Academic Exchange Service). At the Bergische Universität Wuppertal I have been a
member of the Institute of Mathematical Modelling, Analysis and Computational Mathe-
matics (IMACM). Correspondingly this thesis consists of parts devoted to modeling, anal-
ysis and computational mathematics.
In concluding this treatise, I would like to thank all those without whose help this thesis
would not have been possible.
First of all, I want to thank Prof. Dr. Michael Günther for his supervision and the op-
portunity to work at the Chair of Applied Mathematics in Wuppertal. I thank Prof. Dr.-
Ir. Herbert De Gersem for his supervision, scientific guidance who was a source of moti-
vation and many fruitful discussions. Herbert De Gersem introduced me to the electrical
engineering community and it was only due to his support that I was able to understand
and tackle the engineering problems of this treatise. Finally I thank Dr. Andreas Bartel for
his scientific supervision, his support of my research and his rigorous mathematical review
of all results.
I would like to thank Prof. Dr. Markus Clemens, Prof. Dr. Karl Meerbergen and
Prof. Dr. Caren Tischendorf who agreed to be members of my doctoral committee.
During my doctoral studies I collaborated with external researchers. The topics ranged
from general discussions of mathematical ideas to the solution of concrete problems. In this
context I am particularly grateful for the cooperations with Sascha Baumanns, Dr. Markus
Brunk and Prof. Dr. Markus Clemens.
I thank my colleagues from Wuppertal and Leuven, present and former, for many math-
ematical discussions but also for being friends who it was always a pleasure to work with.
I would like to name just a few. In particular I want to thank Dr. Markus Brunk and
Dr. Michael Striebel from Wuppertal and Bart Van Damme and Stijn Rebry from Leu-
ven/Kortrijk. Further thanks go to my Italian colleagues Dr. Giuseppe Alı́, Dr. Massimil-
iano Cuplo and Dr. Carmelo Scordia.
I would also like to thank Alun Davies for proof reading this thesis.
Finally, I would like to thank my parents Brigitte and Fred Schöps, my girlfriend Julia
Heuer and my friends for giving me all their love and support.
III
Contents
List of Symbols IX
1 Introduction 1
1.1 Multiscale and Multirate . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Related Works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 Electromagnetism 5
2.1 Maxwell’s Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.1 Boundary and Initial Conditions . . . . . . . . . . . . . . . . . . . . 6
2.1.2 Partitioning into Regions and Materials . . . . . . . . . . . . . . . . 8
2.2 Space Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2.1 Maxwell’s Grid Equations . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.2 Semi-discrete Problem Formulation . . . . . . . . . . . . . . . . . . 12
2.2.3 Gauging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
4 DAE-Index Analysis 29
4.1 Tractability Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.2 Index Analysis for the Field/Circuit Problem . . . . . . . . . . . . . . . . . 30
4.2.1 Index-1 Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.2.2 Index-2 Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5 Multirate Methods 39
5.1 Single Rate DAE Time-integration . . . . . . . . . . . . . . . . . . . . . . 40
5.2 Multirate Bypassing of MQS Schur Complements . . . . . . . . . . . . . . 41
5.2.1 Bypassing as Multirate Time-Integration . . . . . . . . . . . . . . . 44
5.2.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.3 Multirate Cosimulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.3.1 Abstract DAE-DAE Coupling . . . . . . . . . . . . . . . . . . . . . 47
V
Contents
6 Numerical Examples 84
6.1 DAE-Index in Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
6.2 Multirate Bypassing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.2.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
6.3 Cosimulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
6.3.1 Field/Circuit Problem . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.3.2 Semiconductor/Circuit Problem . . . . . . . . . . . . . . . . . . . . 95
6.4 Domain Substructuring of a Transformer . . . . . . . . . . . . . . . . . . . 99
6.4.1 2D Test Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
6.4.2 3D Test Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
6.4.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
7 Conclusions 103
7.1 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
Bibliography 116
VI
List of Figures
VII
List of Figures
VIII
List of Symbols
Electromagnetism
A~ magnetic vector potential . . . . . . . . . . . . . . . . . . . . . . . . . 5
B~ magnetic flux density . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
D~ electric flux density . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
E~ electric field strength . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
ε electric permittivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
Γ boundary of the computational domain Ω . . . . . . . . . . . . . . . . . 8
H~ magnetic field strength . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
J~ electric current density . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
~n outward normal unit vector (at boundary) . . . . . . . . . . . . . . . . . 7
ν magnetic reluctivity (inverse permeability) . . . . . . . . . . . . . . . . . 5
Ω computational domain . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
(r)
Ω region, i.e., r-th subdomain of Ω . . . . . . . . . . . . . . . . . . . . . . 8
ϕ electric scalar potential . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
ρ electric charge density . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
σ electric conductivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
Discrete Electromagnetism
⌢
a line-integrated magnetic vector potential . . . . . . . . . . . . . . . . . 11
⌢
⌢
b area-integrated magnetic flux density . . . . . . . . . . . . . . . . . . 10
C curl matrix (C̃ = C⊤ is the dual curl matrix) . . . . . . . . . . . . . . 10
⌢
⌢
d area-integrated electric flux density . . . . . . . . . . . . . . . . . . . 10
S divergence matrix (S̃ is the dual divergence matrix) . . . . . . . . . . . 10
⌢
e line-integrated electric field strength . . . . . . . . . . . . . . . . . . . 10
⌢
h line-integrated magnetic field strength . . . . . . . . . . . . . . . . . . 10
⌢
⌢
j area-integrated electric current density . . . . . . . . . . . . . . . . . . 10
kν curl-curl matrix (chord) . . . . . . . . . . . . . . . . . . . . . . . . . 14
Kν differential curl-curl matrix . . . . . . . . . . . . . . . . . . . . . . . 14
Mε electric permittivity matrix . . . . . . . . . . . . . . . . . . . . . . . 11
Mν (nonlinear) magnetic reluctivity matrix . . . . . . . . . . . . . . . . . 11
Mσ electric conductivity matrix . . . . . . . . . . . . . . . . . . . . . . . 11
Ωe element for discretization (subdomain of Ω) . . . . . . . . . . . . . . . 10
Φ discrete electric scalar potential . . . . . . . . . . . . . . . . . . . . . 11
Qσ constant projector onto Ker Mσ . . . . . . . . . . . . . . . . . . . . . 15
Pσ complementary projector Pσ = I − Qσ . . . . . . . . . . . . . . . . . . 15
q volume-integrated electric charge density . . . . . . . . . . . . . . . . 10
IX
List of Figures
Electric Network
A⋆ (reduced) circuit incidence matrix . . . . . . . . . . . . . . . . . . . . 18
C Jacobian of the capacitance relation (matrix of lumped capacitances) . . 19
gR constitutive relation for resistances . . . . . . . . . . . . . . . . . . . . 18
iD currents through semiconductor devices . . . . . . . . . . . . . . . . . 19
iL current through inductances . . . . . . . . . . . . . . . . . . . . . . . 18
iM currents through MQS field devices . . . . . . . . . . . . . . . . . . . 19
is constitutive relation for current sources . . . . . . . . . . . . . . . . . 18
iV current through voltage sources . . . . . . . . . . . . . . . . . . . . . 18
L Jacobian of the inductance relation (matrix of lumped inductances) . . . 19
φ magnetic flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
φL constitutive relation for inductances . . . . . . . . . . . . . . . . . . . 18
q electric charge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
qC constitutive relation for capacitances . . . . . . . . . . . . . . . . . . . 18
G Jacobian of the resistance relation (matrix of lumped conductances) . . 19
u nodal potentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
vs constitutive relation for voltage sources . . . . . . . . . . . . . . . . . 18
MQS device
γ arbitrary path through the conductor part of the MQS device . . . . . 21
iM currents through the MQS devices . . . . . . . . . . . . . . . . . . . 24
Lstr nonlinear inductance matrix (for stranded conductors) . . . . . . . . . 25
Mσ̃ modified electric conductivity matrix (for stranded conductors) . . . . . 25
(m)
Ω conductor region Ω(m) ⊂ Ω(M ) . . . . . . . . . . . . . . . . . . . . . . 20
XM full coupling matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
X̄M sparse coupling matrix . . . . . . . . . . . . . . . . . . . . . . . . . . 22
Rsol diagonal matrix of dc resistances (for solid conductors) . . . . . . . . . 25
Rstr diagonal matrix of dc resistances (for stranded conductors) . . . . . . . 25
vM voltage drops at the MQS devices . . . . . . . . . . . . . . . . . . . . 24
Semiconductor
C doping concentration . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
CD extracted capacitance modeling the displacement current . . . . . . . . 27
iD total current through semiconductor . . . . . . . . . . . . . . . . . . . 27
iSD current through semiconductor without displacement current . . . . . . 27
M⋆ drift diffusion mass matrices (where ⋆ ∈ {n, p}) . . . . . . . . . . . . . 27
K⋆ drift diffusion stiffness matrices (where ⋆ ∈ {n, p}) . . . . . . . . . . . . 27
µ⋆ mobility parameters (where ⋆ ∈ {n, p}) . . . . . . . . . . . . . . . . . 27
n electron density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
n discrete electron density . . . . . . . . . . . . . . . . . . . . . . . . . 27
(D)
Ω semiconductor domain . . . . . . . . . . . . . . . . . . . . . . . . . . 26
X
List of Figures
p hole density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
p discrete hole density . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
q elementary charge . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
r⋆ discrete recombination term (where ⋆ ∈ {n, p}) . . . . . . . . . . . . . 27
R generation-recombination term . . . . . . . . . . . . . . . . . . . . . . 27
UT thermal voltage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
Multirate Applications
En estimated energy level . . . . . . . . . . . . . . . . . . . . . . . . . . 44
gD nonlinear semiconductor resistance . . . . . . . . . . . . . . . . . . . . 73
(i)
Lh generalized inductance matrix (step size h, iteration i) . . . . . . . . . 43
LC maximum Lipschitz constant (w.r.t. zD ) . . . . . . . . . . . . . . . . . 76
LD maximum Lipschitz constant (w.r.t. yC ) . . . . . . . . . . . . . . . . . 76
LM extracted nonlinear inductance from field model . . . . . . . . . . . . . 65
φM magnetic flux through field model . . . . . . . . . . . . . . . . . . . . 65
ψ time integrated voltage (energy estimation) . . . . . . . . . . . . . . . 69
gM dc resistance of field model . . . . . . . . . . . . . . . . . . . . . . . . 65
Index Analysis
A mass matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
XI
List of Figures
XII
1 Introduction
Faster innovation cycles and the increasing complexity of electrical devices necessitate their
simulation. Nowadays engineers design new devices on desktop computers thus reducing
the construction of prototypes to a minimum.
Such virtual prototyping consists of several steps: the engineer draws the geometry of
a particular device (e.g. a transformer) using a computer aided design tool. He assigns
material data and boundary conditions. Then he carries out a simulation to verify that
his device is adequate for its purpose. For example the engineer might be interested in the
behavior of a transformer after connecting it to the electric power network. The power
network supplies energy and this excites the device. Now, to predict the performance
of it, the physical quantities must be computed by solving differential equations for the
given excitation, e.g., the time-varying currents and voltages. Typically, the corresponding
simulation in time domain (‘time-integration’) is computationally very expensive because
millions of equations must be solved repeatedly at various time points.
The main goal of this treatise is to develop methods speeding up those simulations.
The idea behind it is trivial: different parts of a given problem are not equally important.
For example the simulation of a power plant is not necessary to predict the behavior
of a transformer plugged into an electric socket. Thus some effects can be neglected or
simulated by simpler models, while other problem parts are crucial and thus need models
of higher quality. Obviously, higher quality automatically implies higher computational
costs.
This treatise deals with the coupling and simulation of device models of different quality
and scales, in other words multirate methods are applied to a hierarchy of multiscale
models.
1
1 Introduction
L R Rstr,1 Rstr,2
voltage [V]
200
u1 u3 -400
0 0.005 0.01 0.015 0.02 0.025 0.03
time [s]
Figure 1.1: Introduction to multirate. Example configuration (a) and multirate behavior in
example (b), see [118].
ular allows for multirate and multimethod, but the presence of algebraic constraints can
handicap higher-order integration and may cause divergence. In this treatise it will be
shown that the stability and convergence of a cosimulation can be verified by analyzing
the coupling interface. The design of the coupling interface is crucial: it determines the
convergence speed and may or may not support multirate time-integration.
In problems with different time-scales the representation of the slow subsystem by a re-
duced order model for as many time steps as possible is advantageous, because the weaker
the coupling, the more efficient the cosimulation. For example in most magnetoquasistatic
applications a lumped linear inductance model extracted from the field equations is suf-
ficient to reflect the device characteristics for a reasonable time span. The nonlinearity
of the inductance due to magnetic saturation can be slowly resolved in time because of
the inertia of energy transport. Thus this treatise proposes to solve a network problem
consisting of entirely lumped models and fit the parameters on the fly by an outer iteration
with PDE device models, Fig. 1.2b.
For that the cosimulation approach is adapted and sufficient conditions for stability
and convergence are derived and the convergence rate is well understood. The parameter
coupling approach is rather general and can be applied to any DAE-DAE problem and
thus a second example from electrical engineering is given: an electric circuit coupled to a
semiconductor model.
2
1.3 Overview
102 iterations 1
iterations 2
iterations 3
100 iterations 4
relative error in current
10-2
PDE network
10-4
10-6
model model
10-8
Figure 1.2: Introduction to cosimulation. Order in dependence of window size (a) and cosim-
ulation (b).
1.3 Overview
The focus of this thesis is the efficient time-transient simulation of device/circuit coupled
problems by exploiting properties of the multiscale models (time scales, symmetries etc).
Several methods are proposed in this thesis, but the cosimulation is the most natural
choice. One has the opportunity to freely choose different methods for each subproblem.
On the other hand cosimulation requires more analysis of the problem formulation and the
interfaces.
The structure of the work is as follows: the next chapter introduces the electromagnetic
field problem. The quasistatic approximations to Maxwell’s equations are recapitulated
and the space discretization is briefly discussed. Maxwell’s equations are the foundation
3
1 Introduction
of the multiscale models that are derived in Chapter 3. It is subdivided into Section 3.1
on electric network modeling, Section 3.2 on magnetoquasistatic devices and finally Sec-
tion 3.3 that recapitulates the drift-diffusion model of a semiconductor. The results in
those Sections are standard, although Section 3.2 contains some new aspects from [118, 8].
Chapter 5 is the core of this thesis. It starts with a brief summary of numerical time-
integration (Section 5.1). In the following Section 4 the time-discretization properties of the
field/circuit coupled problem are analyzed using the differential-algebraic index concept.
The section is based on a recent publication compiled in cooperation with Sascha Bau-
manns and Andreas Bartel, [8], but also includes results from material that was presented
at SCEE 2008, [117]. Section 5.2 discusses a new multirate technique based on standard
time-integration: the field subsystem is reduced by Schur complements and updated ac-
cording to its energy-level. This approach was presented at SCEE 2010 and was devised in
collaboration with Andreas Bartel and Herbert De Gersem, [116]. The cosimulation stabil-
ity and convergence theory for general index-1 DAEs is presented in Section 5.3. It utilizes
results from [7], but features the multi-window propagation of splitting errors in dynamic
iteration schemes, similarly to [5]. Furthermore important special cases are studied with
more rigor. The dynamic iteration schemes are applied to the field/circuit (Section 5.3.5)
and semiconductor/circuit coupled problems (Section 5.3.6). The convergence theory and
its application are the result of a collaboration with Andreas Bartel, Markus Brunk and
Michael Günther, [9]. The multirate coupling in Section 5.3.5 is developed together with
Herbert De Gersem and was presented at Compumag 2009 and EPNC 2010, [118]. Fur-
thermore some of the convergence studies for the semiconductor case in Section 5.3.6 were
presented at SCEE 2010, [1]. In Section 5.4 a domain substructuring approach is discussed
that exploits the particular structure, i.e., linear/nonlinear and conductive/nonconductive
domains, of the eddy current problem. The method was presented at IGTE 2010 in col-
laboration with Markus Clemens, Andreas Bartel and Herbert De Gersem, [36].
Chapter 6 shows the significance and practicability of the various results and methods:
the index-2 time-integration of a magnetostatic device model (Section 6.1), the multirate
bypassing of a magnetoquasistatic transformer model (Section 6.2), the cosimulation of
both MQS and semiconductor devices coupled to circuits (Section 6.3.1 and 6.3.2), and
finally the domain substructuring approach is applied to 2D and 3D transformer models
(Section 6.4). With the exception of the semiconductor/circuit cosimulation all the ex-
periments were performed using the software FIDES (within Octave/CoMSON DP). The
numerical results for the semiconductor problem were obtained by Markus Brunk using
Matlab, [9]. The thesis closes with conclusions in Chapter 7.
4
2 Electromagnetism
2.1 Maxwell’s Equations
Electromagnetic phenomena are described on the macroscopic level by Maxwell’s partial
differential equations, [93]
~ ~
~ = − ∂B ,
∇×E ~ = ∂ D + J~ ,
∇×H ~ =ρ,
∇·D ~ =0,
∇·B (2.1)
∂t ∂t
where E~ = E(~ ~ r, t) is the electric field strength, B ~ = B(~ ~ r, t) the magnetic flux density,
~ = H(~
H ~ r , t) the magnetic field strength, D ~ = D(~~ r , t) the electric flux density, ρ = ρ(~r, t)
the electric charge density and J~ = J(~ ~ r , t) is the electric current density. All quantities
depend on space ~r ∈ Ω and time t ∈ I. These equations are related to each other by the
additional material relations
~ = ~~εE,
D ~ J~ = ~~σ E,
~ ~ = ~~ν B,
H ~ (2.2)
where the permittivity ~~ε = ~~ε(~r), conductivity ~~σ = ~~σ (~r) and reluctivity (inverse perme-
ability) ~~ν = ~~ν (~r, ||B||
~ 2 ) are rank-2 tensors. All field quantities can be expressed by the
magnetic vector potential (MVP) A ~ : I × Ω → R3 and its integration constant, the electric
scalar potential ϕ : I × Ω → R
~
~ = ∇×A
B ~ and ~ = − ∂ A − ∇ϕ.
E (2.3)
∂t
Assumption 2.1 (Quasistatics). We assume that the problems considered in this treatise
are either
~
∂D
i) magnetoquasistatic ( MQS), i.e., ∂t
=0
~
∂B
ii) electroquasistatic ( EQS), i.e., ∂t
=0
5
2 Electromagnetism
Assumption 2.1 i) and inserting Maxwell’s equations into each other, starting from
Ampère’s law, i.e., the second equation of (2.1), gives the magnetoquasistatic approxi-
mation, [71]
~
~~σ ∂ A + ∇ × ~~ν ∇ × A
~ = −~~σ ∇ϕ. (2.4)
∂t
This curl-curl equation is of parabolic type, but if we further disregard changes in the
magnetic field (Ass. 2.1 iii), we end up with the magnetostatic approximation
∇ × ~ν ∇ × A = −~~σ ∇ϕ
~ ~ (2.5)
which is an elliptic problem. From the electric permittivity material law and Ohm’s Law,
i.e., the first two material relations of (2.2) and using the assumption of statics (Ass. 2.1
iii) one deduces the electrostatic approximations
−∇ · ~~ε∇ϕ = ρ and −∇ · ~~σ ∇ϕ = 0 , (2.6)
that are elliptic equations of Poisson’s type, [71]. The formulations above are not the only
possible choice. There are other formulations (e.g. [18]) based on various quantities. Those
are not used in the following.
The choice of the BCs depends on several restrictions: on one hand at least one boundary
must be defined as a Dirichlet type to yield a uniquely solvable problem, [105]. On the
6
2.1 Maxwell’s Equations
Figure 2.1: Boundary conditions and flux lines. (a) Dirichlet condition imposed on left and
bottom boundary, (b) Neumann condition on left and Dirichlet on bottom boundary, (c) anti-
periodic conditions on both boundaries.
other hand the circuit coupling, as discussed in Section 3.2, adds further constraints. In
the case of excitation by source terms, the boundaries are commonly of Dirichlet type
(unless symmetries are exploited), while a coupling via boundary conditions requires more
sophisticated assumptions, [75]. In the magnetic vector potential formulation (2.4) the
conditions read thus:
Definition 2.1 (Dirichlet boundary). Let Γdir denote the Dirichlet boundary and ~n its
outward normal unit vector then
~ r) × ~n = A
A(~ ~ dir with ~r ∈ Γdir
sets the tangent components of the magnetic vector potential at the boundary to the
~ dir . In the case of A
prescribed value A ~ dir = 0 the condition is called homogeneous Dirichlet
BC.
The condition is called electric BC, often also described as “flux wall” or “current gate”.
It does not allow the magnetic flux to pass through the border (the flux lines stay parallel
to the boundary). On the other hand it corresponds to perfectly conducting borders, i.e.,
with infinite electric conductance, see Fig. 2.1a.
Definition 2.2 (Neumann boundary). Let Γneu denote the Neumann boundary and ~n its
outward normal unit vector then
~~ν ∇ × A(~
~ r ) × ~n = H
~ neu with ~r ∈ Γneu
sets the tangent components of the magnetic field strength at the boundary to the pre-
~ neu . In the case of H
scribed value H ~ neu = 0 the condition is called a homogeneous Neumann
BC (for the magnetic vector potential).
The condition is called magnetic BC, often also described as “flux gate” or “current
wall”. It forces the magnetic flux to leave the computational domain perpendicular to the
border and corresponds on the other hand to perfectly resistive borders, see Fig. 2.1b.
7
2 Electromagnetism
boundaries connected by a mapping ~s(~r) : Γanti,+ → Γanti,- and let ~n+ , ~n− describe their
outer normal unit vectors, then
~ r ) · ~n+ = −∇ × A(~
∇ × A(~ ~ s) · ~n−
forces the same magnetic flux entering/leaving through one boundary Γanti,+ to equally
enter/leave through the other one Γanti,- (for an equal spatial distribution), see Fig. 2.1c.
Finally the problem description is complete with the specification of the time interval of
interest, i.e., I := [t0 , te ] and an initial value for the parabolic setting.
Assumption 2.4 (Initial value). The initial conditions are given for t0
~ 0 , ~r) = A
A(t ~ 0 (~r) for all ~r ∈ Ω. (2.7)
with interfaces Γ(r,s) = Ω̄(r) ∩ Ω̄(s) that are piecewise C 2 and Lipschitz-continuous, [72].
We assume that the materials properties are constant within each region.
Assumption 2.6 (Permittivity). We assume an isotropic permittivity that is given by a
positive constant in each region, i.e.,
~
∂D
= ~~ε(~r) = ε(r) I with ε(r) > 0 for ~r ∈ Ω(r) . (2.9)
~
∂E
Assumption 2.7 (Conductivity). We assume a region-wise constant non-negative con-
ductivity, whose axes of anisotropy are aligned with the coordinate system, i.e.,
∂ J~ ~
= ~σ (~r) = ~~σ (r) = diag(σx(r) , σy(r) , σz(r) )
(r)
with σ{x,y,z} ≥ 0 for ~r ∈ Ω(r)
∂E ~
(2.10)
Assumption 2.8 (Reluctivity). We assume that the reluctivity tensor is given by Brauer’s
model, [24]. It is defined as the superposition of an isotropic and anisotropic tensor, i.e.,
8
2.2 Space Discretization
div S̃
ϕ [V] ρ [C/m3 ] Φ [V] q [C]
⌢ [Wb] ⌢
⌢
~ [Wb/m]
A ~ [A/m2 ]
J a j [A]
-grad ]
⊤
S̃
/m [ S]
dd
dd
σ [S
dd
dd
t
t
Mσ
t
−
−
−
−
curl C̃
ε [F/m] ⌢ Mε [F]
~ [V/m]
E ~ [C/m2 ]
D e [V] ⌢
⌢
d [C]
µ [H/m] ⌢
⌢ Mµ [H] ⌢
~ [T]
B ~ [A/m]
H b [Wb] h [A]
curl C
div -grad S S⊤
Figure 2.2: Maxwell’s house. Maxwell’s equations and differential forms, i.e., Tonti’s diagram
for the discrete and continuous formulations, [129, 34].
Only the isotropic material scalar νiso depends nonlinearly on the magnetic flux density
~ 2 and it is assumed that this tensor and all derivatives are positive definite and bounded
||B||
~ iso (r) ~ 2 )B
~ ~ aniso
∂H ∂νiso (||B|| ∂H
and ~~νaniso :=
(r) ~ 2) (r) (r)
ν2 > νiso (||B|| > ν1 , νd,iso := =
~
∂B ∂B~ ∂B ~
9
2 Electromagnetism
primary
primary cell dual
point point
dual dual
cell cell
dual
primary point
cell primary point
Assumption 2.9 (Discretization). The primary grid discretizes the regions, i.e.,
[
Ω̄(r) = Ω̄e with Ωe ∩ Ωj = ∅ for e 6= j , (2.12)
Ωe ⊂Ω(r)
where Ω̄e denotes (the closure of ) elements (triangle, hexahedron etc.) of the primary grid.
This treatises deals with polytopal elements only and thus Ass. 2.9 implies that the
computational domain Ω̄ and all regions Ω̄(r) are polytopes. Again, this assumption is
only made for simplicity of notation and can be overcome, e.g. by curved elements, [33].
Definition 2.4 (Primary objects). The pairwise intersection of all elements of the primary
grid yields vertices (0D polytope) denoted by Pi with i = 1, . . . , n0 , edges (1D polytope)
denoted by Li with i = 1, . . . , n1 , facets (2D polytope) denoted by Ai with i = 1, . . . , n2 ,
and cells (3D polytope) denoted by Vi with i = 1, . . . , n3 .
Similarly, one defines objects for the dual grid, see Section A.1.1.
with discrete curl operators C and C̃, divergence operators S and S̃ on the primal and
dual grid, respectively. The variables are time dependent quantities, whose diacritical
bows correspond to the dimension of the underlying object, i.e., edges, facets and cells,
[21]. On the one hand there are the line-integrals of electric and magnetic field strength ⌢
e(t)
⌢
and h(t). They are located on primary and dual edges. For example the line-integrated
electric field strength and magnetic vector potential are time dependent functions I → Rn1
10
2.2 Space Discretization
measured in Volt (V) and Weber (Wb), respectively, see Fig. 2.2b. They read
Z Z
⌢
ei := ~
E · ds and ⌢
ai := ~ · ds
A i = 1, . . . , n1 , (2.14)
Li Li
and are defined on the primary edges Li ⊂ Ω. The scalar n1 denotes the number of all
primary edges. The other quantities, i.e., the
⌢
discrete
⌢
magnetic
⌢
flux density, the discrete
⌢ ⌢ ⌢
current density and the displacement field b(t), j (t) and d(t), respectively, are defined as
surface integrals. The first is located on primary facets and the latter ones are defined on
dual facets. Finally the electric flux density q is defined on dual volumes.1
The constitutive material laws relate the primary and dual quantities to each other, see
Fig. 2.2b, such that Maxwell’s Grid Equations are completed by
⌢
⌢ ⌢
⌢ ⌢ ⌢
⌢
d = Mε ⌢
e, j = Mσ ⌢
e, h = Mν b , (2.15)
⌢
⌢
where the matrices Mε , Mσ and Mν = Mν (b) represent the permittivities, conductivities
and (nonlinear) reluctivities, respectively. The matrices of permittivities and reluctivities
are always positive definite, while the matrix of conductivities is typically positive semi-
definite due to non-conducting parts in the computational domain, see Section 2.2.3.
In the FEM context the material matrices are typically constructed by assembling local
material matrices defined per element. On the other hand in the classical FIT there is no
need for the construction of local matrices due to the sophisticated index mappings on the
underlying (structured) grid. Nonetheless for unity of notation the element-wise approach
can be used for both methods, [42, 41]. This approach is generalized in Section A.2.2 to
apply to anisotropic materials.
For both discretization the material properties were assumed region-wise constant,
Ass. 2.6-2.8 and this property is inherited by the elements. It holds for the element in
the r-th region (~re ∈ Ωe ⊂ Ω(r) ):
~~ε(~re ) = ~~ε(r), ~~σ (~re ) = ~~σ (r) and ~~ν (~re , ||Be||2 ) = ~~ν (r) (||Be ||2)
with constant permittivity ~~ε(r) , conductivity ~~σ (r) and reluctivity ~~ν (r) . The nonlinear reluc-
tivity is evaluated using an element-wise averaged magnetic flux density ||Be ||2 .
The rough material approximation above implies that curved material boundaries are
discretized by staircase approximations. This simplifies the notation because no material
parameter averaging is necessary, but in practice this limitation is overcome by subgridding
or other more elaborate schemes, [126, 83]. In either case the global material matrices
X X ⌢
⌢
X ⌢
⌢
Mε := Mε,r , Mσ := Mσ,r , Mν (b) := Mν,r (b) , (2.16)
r r r
can be written as the sum of region-wise matrices, assembled from the element contribu-
tions, see Section A.2.2.
Definition 2.5 (Regional index sets). The indices of all primary objects belonging to a
region can be expressed by index sets. For example the edges for the (closed) region r are
1
For simplicity of notation the electric flux density is denoted without bows, although located on volumes.
11
2 Electromagnetism
given by
L Ω̄(r) := i ∈ N | Li ∈ Ω̄(r) .
Note that an object may be included in one or more index sets because of intersections of
closures (for example the region-wise matrices overlap only at region boundaries).
Lemma 2.10 (Conductive regions). On a (conductive) region r with ~~σ (r) > 0 the region-
wise material matrix is positive definite, i.e.,
⌢⊤
e Mσ,r ⌢e>0 if i ∈ L Ω̄(r) such that ⌢
e i 6= 0.
Proof. Direct consequence of the element assembly, see for example the FIT case in Sec-
tion A.2.2.
⌢
⌢ d⌢
b = C⌢
a and ⌢
e=− a + S̃⊤ Φ . (2.17)
dt
For 3D discretizations the matrix C is singular, such that the vector potential is not
uniquely defined (as in in the continuous case). Section 2.2.3 discusses possible regular-
ization techniques. In further correspondence with the continuous formulation, i.e., (2.4)
and (2.5), one obtains the semi-discrete curl-curl equations for the magnetoquasistatic
(Ass. 2.1 i) and magnetostatic regimes (Ass. 2.1 iii) as
d⌢ ⌢
⌢ ⌢
⌢
Mσ a + C̃Mν (C⌢
a)C⌢
a = js and C̃Mν (C⌢
a)C⌢
a = js (2.18)
dt
⌢
⌢
with a given (facet-integrated) source current density j s . The degrees of freedom (DoF)
of this problem are the line-integrated vector potentials ⌢ a. The mass matrix Mσ is the
algebraic origin for eddy currents and thus magnetoquasistatic problems are called eddy-
current problems. The stiffness matrix is the derivative of the curl-curl term with respect
to ⌢
a. It is called the (differential) curl-curl matrix
∂ ⌢ ⌢
C̃Mν (C a)C a = C̃Mν,d (C⌢
a)C , (2.19)
∂⌢a
whose derivation is given in more detail for the FIT case in Section A.3. Finally the
semi-discrete Poisson problems in the electrostatic approximation (Ass. 2.1 iii) read
12
2.2 Space Discretization
(a) spy plot of (regularized) curl-curl matrix (b) curl-curl discretization stencil
Figure 2.4: The regularized curl-curl matrix consists of entries from the original matrix and
additions from the regularization, see [36].
similarly to (2.6). The equations (2.18) and (2.20) yield well-posed problems when bound-
ary conditions and for the initial value problems the corresponding initial values are pre-
scribed. For the magneto(quasi)static problems the discretization of the boundary condi-
tions, Section 2.1.1 yields restrictions for MVP; for example in the homogeneous Dirichlet
case the tangential components vanish
⌢
ai = 0 if i ∈ L (Γdir ) ,
⌢
ar = ϑ ⌢
as if r ∈ L (Γanti,+ ) and Ls = ~s(Lr )
where ϑ = ±1 depends on the orientation of the edges and the type of periodicity. The
homogeneous Neumann condition is naturally fulfilled by the discretization. In any case
the conditions are built into the system matrices, such that superfluous components are
removed, e.g. [103, Section 8.4.5].
2.2.3 Gauging
On a 3D domain the curl-operator inherits the non-uniqueness from its continuous coun-
terpart, i.e.,
whereby the gradient operator −S̃⊤ has full column rank, Section A.1.1. Thus the curl-
curl matrix (2.19) has zero eigenvalues and it is not invertible. This is inconvenient for
the structural analysis and renders direct matrix factorizations impossible. This can be
overcome by regularizations in a similar way to the continuous Coulomb gauge, [20, 37,
81]. Those regularizations remove the nullspace by addition
kν (⌢
a) := C̃Mν (C⌢
a)C + Zσ , where Φ⊤ S̃Zσ S̃⊤ Φ > 0 for all Φ 6= 0, (2.22)
such that kν (⌢
a) is positive definite (and Zσ is positive semi-definite).
13
2 Electromagnetism
eig(C̃Mν C)
eig(Kν ) regularization
eig(λMσ + C̃Mν C)
eig(λMσ + Kν ) regularization
Figure 2.5: Shift of eigenvalues due to grad-div regularization, parameter λ = 103 , cf. [36].
is called the Grad-Div regularization. It utilizes the gradient and divergence operators S̃
from above and suitable (artificial) material matrices: M1 maps primary edges to dual
facets and the norm matrix M2 maps dual points to primary volumes, [37].
For a homogenous material distribution and an equidistant grid, the Grad-Div regular-
ized curl-curl matrix (see Fig. 2.4) corresponds to the discrete vector Laplacian (−∇2 ),
[37]. Suitable choices for the material matrices M1 and M2 are discussed in [36].
When considering a magnetoquasistatic problem (2.18), we are interested in the regu-
larity of the matrix pencil (c.f. Fig. 2.5)
[Mσ , Kν (⌢
a)] := λMσ + Kν (⌢
a), with λ > 0 (2.25)
which occurs naturally when solving problems in the time domain, e.g. λ = 1/h for the
implicit Euler, see Section 5.1. The conductivity matrix Mσ already shifts some of the
zero eigenvalues of the curl-curl matrix onto the positive real axis, see Section 5.2. The
regularization (2.22) would affect those eigenvalues superfluously. Hence a gauging in the
nonconductive domain is sufficient, [37, 36]. This domain can be addressed formally by the
following projector:
14
2.2 Space Discretization
Pσ = M+σ Mσ (2.26)
In the MQS case it is sufficient to regularize only in the non-conductive parts, e.g. by
restricting the regularization to Zσ = Q⊤σ Zσ Qσ . The restriction can be incorporated into
the artificial material matrices in Definition 2.7. This is the ‘Grad-Div Gauging Type II’
in [36]. In either case the following result holds (in particular for the Grad-Div Gauging)
2. else Mσ S̃⊤ Φ 6= 0 and thus the first summand is positive and the second summand is
non-negative.
Independently of the particular choice for the regularization, the following assumption is
made:
Assumption 2.13 (Gauge). In the magnetostatic case the curl-curl matrix is fully regu-
larized (Definition 2.6) with a positive semi-definite matrix Zσ such that kν and Kν are
positive definite. In the MQS case the same regularization is applied but only for elements
in Ker Mσ .
Remark 2.1 (Weak Gauging). No explicit regularization is necessary when solving linear
systems with system matrices of the form (2.25) with Krylov subspace methods. They
exhibit a weak gauging property, [37], i.e., they will not alter the initial guess in the
nullspace of the system matrix, [78].
15
2 Electromagnetism
2.3 Conclusions
This chapter has introduced briefly the continuous and discrete formulation of Maxwell’s
equations and their common quasistatic approximations. General properties of the discrete
material and differential operators for spatial discretizations based on edge-elements have
been presented. The notation was based on FIT, but without loss of generality.
Special emphasis was put on a uniquely solvable discrete problem formulation, because
the following chapters will discuss time-integration and the corresponding theory and meth-
ods rely on uniqueness.
This chapter has established the electromagnetic framework into which the following
chapter embeds (multiscale) models for electromagnetic networks and distributed device
models, i.e., the magnetoquasistatic and semiconductor devices.
16
3 Multiscale Device Models
In the previous chapter the electromagnetic phenomena of a single problem (‘device’) were
described on the macroscopic scale. In electric circuit simulation, e.g. [32], the behavior
of the combination of a large number of such devices is in focus. This chapter focusses on
effects in circuit simulation that must be resolved on different spatial scales. This will be
clarified by the description of a transformer and a semiconductor, whose spatial dimensions
differ by many orders of magnitude.
A multiscale simulation with PDE models of all devices is obviously computationally
’
inappropriate. Instead the complexity of some devices is reduced, e.g., by disregarding the
spatial distribution (lumping). Then a network of idealized basic elements is considered,
where each element describes only one effect, e.g. resistances, inductances, capacitances
and sources. A single physical device is approximated by several basic elements. They are
called equivalent, compact or companion models. Most often, the network equations for
all devices are set up element-wise according to the modified nodal analysis (MNA), which
is introduced in Section 3.1. Mathematically speaking the result is a set of differential
algebraic equations (DAEs).
However, many devices cannot be described sufficiently accurately in terms of a few
idealized lumped elements. The corresponding equivalent circuits become too complex
and contain hundreds of parameters, most of them with no physical interpretation, [44].
Especially if one is interested in the effects on all scales at the same time, a hierarchical
modeling of PDE devices and the electric network DAEs is the only viable strategy (mixed
mode simulation, [62]). The result is a system of partial differential algebraic equations
(PDAEs), which models the distributed effects by computationally expensive PDE models
only where necessary.
R1 R2
1 2 3 4 5
v(t) L1 L2 Rload
0
6
PDE
17
3 Multiscale Device Models
In the next sections we present different models (see Fig. 3.1) and couplings. The first
model in Section 3.1 is a recapitulation of the electric network description in terms of the
modified nodal analysis (MNA). This formulation is the standard approach in commer-
cial tools for circuit analysis (e.g. SPICE-like simulators as Titan by Infineon), [108]. In
Section 3.2 the magnetoquasistatic field device and the corresponding field/circuit cou-
pling is established. The device model yields a current/voltage relation and this makes
the circuit coupling straightforward, independently of the type of circuit analysis. This
treatise follows [113, 8] and features the MNA. Other approaches use the standard nodal
analysis [46, 137] or study the coupling and numerical treatment for a coupled system
using loop/branch analysis, e.g. [15]. Finally in Section 3.3 the standard drift-diffusion
model of a semiconductor is reproduced, [123, 119, 62]. This well-known model features
a microscopic correction of the electrostatic Maxwell’s equations in terms of holes and
electrons.
where each row of A⋆ refers to a network node. One node is identified as the mass node
(‘ground’) and its row is skipped in the reduced matrix A⋆ . The flux/charge oriented
modified nodal analysis (MNA) yields equations of the form, [64, 57, 54]
d
AC q + AR gR (A⊤R u, t) + AL iL + AV iV + AI is (t) = 0, (3.1a)
dt
q − qC (A⊤C u, t) = 0, (3.1b)
d
φ − A⊤L u = 0, (3.1c)
dt
φ − φL (iL , t) = 0, (3.1d)
A⊤V u − vs (t) = 0, (3.1e)
18
3.1 Electric Network
to classes of controlled sources is straightforward, [54]. The unknowns of the system are
most importantly the node potentials u : I → Rnu (without ground), they correspond
to the voltage drop between the node and ground. Further unknowns are the currents
iL : I → RnL , iV : I → RnV through inductors and voltage sources and the charges
and fluxes q : I → RnC and φ : I → RnL , respectively (where nL , nV and nC denote
their respective number), see [57, 65]. The problem is completed with (consistent) initial
conditions at time t = t0 , [8].
The flux/charge-oriented MNA above is reduced to the traditional MNA, when the
unknowns for the fluxes and charges are eliminated. This yields a smaller system, but does
not guarantee the conservation of energy, [65]. Structurally both approaches are nearly
equivalent, [54] and in either case a mathematically consistent description must fulfill
some topological conditions. The constitutive relations should be passive, i.e.,
Assumption 3.1 (Local passivity). The functions qC (v, t) , φL (i, t) and gR (v, t) are con-
tinuous differentiable with positive definite Jacobians:
19
3 Multiscale Device Models
im im
Ω(m) vm Ω(m) vm
Figure 3.2: MQS device types. Solid conductor connected at boundary and forming a loop
within the domain.
The network is now soundly defined. It is the coupling framework for all devices. The
elements (multiscale models from 0D to 3D) may communicate only via lumped network
quantities. Thus the following PDE models must have corresponding boundary conditions
and source terms, such that physical correctness of the overall model is ensured. This is
the topic of the following sections.
20
3.2 Magnetoquasistatic Conductor Model
im contact at
the boundary
edges
Ω(m)
Ω S reference plane
Γref vm
contact at
the boundary
(a) Solid conductor coupling, [71] (b) Cartesian reference plane, see [43]
Figure 3.3: MQS device coupling. Source term coupling at the reference plane.
(albeit different from the engineering point of view), [117]. The important differences are
discussed in Section 3.2.1.
Let us consider in the following a single (solid) conductor (see Fig. 3.3a). Its region shall
be denoted by Ω(m) . It corresponds to a circuit branch, which is connected by two perfect
conducting contacts. The 0D-voltage drop vm must be distributed onto the 3D-domain;
this defines an electric field.
For that we denote by γ an arbitrary path through the conductor domain and S is the
area inside the loop described by the path γ, see Fig. 3.3a. The wires connecting the
conductor to the circuit are regarded as 0D objects and therefore they are not considered
in the field model, see [71, Chapter 8]. The voltage drop is given by Faraday’s Law
Z Z
~ · d~s − vm = − d ~ · dA
~.
E B
dt
γ S
d ~ ~ =∇×A
~ and applying Stokes theorem yields
Replacing E = − dt A − ∇ϕ, B
Z Z Z
d ~ d ~ · d~s .
− A · d~s − ∇ϕ · d~s − vm = − A
dt dt
γ γ ∂S
Finally the electric field is only defined via its integral. This can be exploited in the
construction of a discrete distribution. The discrete electric field ⌢
e M is located on primary
edges. Since we are only interested in line integrals on those edges, we look at an arbitrary
21
3 Multiscale Device Models
Figure 3.4: Full and sparse coupling. Stranded conductor coupling vectors and for the first
coil of a transformer model discretized by FIT. This corresponds to a 1A excitation.
discrete path ⌢ γ ∈ {−1, 0, 1}n1 from one contact to the other (within Ω(m) ). Due to the
linearity of Ohm’s law, it is sufficient to consider an applied voltage vm = 1V and define
a corresponding distribution matrix X̄m ∈ Rn1 , such that X̄⊤m ⌢γ = 1 for all ⌢
γ.
In [43, 8] a construction is proposed that imposes the voltages only onto edges that cross
a reference plane Γref (see Fig. 3.3b). This implies that the reference plane is a set of dual
facets. The definition of the coupling reads for a Cartesian grid with an orthogonal aligned
reference plane
(
±1 if edge Li ⊂ Ω(m) crosses the reference plane,
X̄m i = (3.4)
0 else,
where the sign depends on the directions of the edges. This coupling vector X̄m is sparse;
it features only nonzero entries at the reference plain (,,2D coupling”). This is computa-
tionally less costly than a full coupling that exhibits entries in the whole 3D domain Ω(m) ,
[43]. Fig. 3.4 shows the coupling pattern of both approaches.
Nonetheless the full coupling is favored in the following, because it can be constructed
in a way that the source current density is divergence free. The divergence property of the
sparse coupling X̄m is corrected by previously solving the following Poisson problem
with unit excitation vm = 1V and projector Qσ,m onto Ker Mσ,m , see Definition 2.8.
Lemma 3.4 (Excitation). Let a sparse coupling vector X̄m be given as defined in (3.4)
and let the potential Φm denote the solution of (3.5). Then the full coupling vector reads
Xm := X̄m − S̃⊤ Φm .
22
3.2 Magnetoquasistatic Conductor Model
The current excitation is divergence free and the voltage drops are only applied in conductive
regions, i.e., Xm = Pσ Xm .
Proof. From Lemma 2.10 follows X̄m ⊂ im Mσ,m and thus Qσ X̄m = 0. The rest is clear
from equation (3.5).
By linearity, the coupling matrix Xm allows the application of an arbitrary voltage
vm = A⊤m u by multiplication. The applied source current density is given by
⌢
⌢
j m = Mσ,m S̃Φ = Mσ,m Xm vm . (3.6)
In the sparse coupling the total current through the conductor is given by integrating over
the reference cross section. In the full coupling approach, the coupling vector Xm averages
over all cross section integrals. We find by using Ohm’s Law in the solid conductor region
⌢
⌢ d⌢
im = X⊤m j = X⊤m Mσ,m ⌢
e = X⊤m Mσ,m Xm vm − X⊤m Mσ,m a. (3.7)
dt
The excitation is easily generalized to the case of an arbitrary number of conductors
⌢
⌢
X
jM = Mσ,m Xm vm . (3.8)
m∈M
where M is the index set of all conductor regions (their total number is denoted by nM =
|M|). Thus the field device has nM terminals. In the case of the transformer shown in
Fig. 3.4, the set consists of two regions that correspond to the primary and secondary coil.
In practice the region-wise construction of multiple neighbored conductors can suffer
from smearing effects, [117]. The conductivities and thus the coupling vectors are not
clearly separated at material boundaries. For a nonzero voltage drop vm the following
currents are not the same
Mσ Xm vm 6= Mσ,m Xm vm ,
if a conductive region is adjacent to the coupling region Ω(m) . This may cause the model to
behave unexpectedly. Let us consider a transformer discretized by 2D FEM, see Fig. 3.5.
A voltage excitation in the (stranded) conductor region Ω(m) causes no current flow in
the region itself because stranded conductors are commonly modeled as nonconductive
material to disable eddy currents, see Section 3.2.1. On the other hand there is a current
within the neighboring iron core due to smearing at the boundary.
This smearing is a typical discretization error. It can be prevented by using adequate
material matrices, a sufficiently large insulation region around each conductor or the error
can be made arbitrarily small by mesh refinements. The problem with the latter approach
is that the insulation layers may be very thin but have to be resolved by the mesh. In the
following it is assumed that there is no smearing.
Assumption 3.5 (No smearing). There is no smearing, i.e., the images of the region-wise
conductivity matrices are distinct
23
3 Multiscale Device Models
(a) transformer (b) dots denote dislocated voltages/currents in the iron core
Figure 3.5: Smearing. Current in the iron core while excitation only applied to the coil Ω(m) .
This allows for the following simplified notation of the device equations using the curl-curl
equation (2.28) and the voltage drop vM = A⊤M u
d⌢
Mσ a + kν (⌢ a = Mσ XM A⊤M u ,
a)⌢ with ⌢
a(t0 ) = ⌢
a0 , (3.10a)
dt
where the total current (3.7) is equivalently given by
iM = X⊤M kν (⌢
a)⌢
a. (3.10b)
Each column of the matrix XM is a coupling vector Xm that corresponds to the region
Ω(m) . The branch currents are gathered in the MQS device’s current vector iM .
with a divergence-free stranded conductor coupling vector Xstr . The construction is similar
to the solid conductor coupling vector, starting from a sparse coupling vector X̄str . The
only difference is the scaling to account for the area of the reference plane and the number
24
3.2 Magnetoquasistatic Conductor Model
of windings, [49, 48, 43]. Using the no-smearing Ass. 3.5, the coupled system for solid and
stranded conductors reads
d⌢
Mσ̃ a + kν ⌢
a = Mσ Xsol vsol + Xstr istr (3.12a)
dt
with gauging, initial and boundary conditions for ⌢
a and the coupling equations
d ⊤
− a + R−1
Xsol Mσ ⌢ sol vsol = isol , R−1 ⊤
sol := Xsol Mσ Xsol , (3.12b)
dt
d
− X⊤str ⌢a + Rstr istr = vstr , Rstr := Xstr M+σ Xstr , (3.12c)
dt
where the superscript ‘+’ denotes a Moore-Penrose pseudo-inverse [61]. The matrices R⋆
denote extracted dc resistances, cf. (3.7). The modified conductivity matrix
X
Mσ̃ := Mσ,r (3.13)
Ωr 6=Ωstr
corresponds to the global conductivity matrix Mσ defined in (2.16), but without the con-
tributions from stranded conductors. This has consequences for the gauging, Section 2.2.3
because modified conductivity matrix Mσ̃ has a larger nullspace
Similarly to the definition of the lumped dc resistances in (3.12b) and (3.12c), an in-
ductance
⌢
matrix can be extracted. The approach corresponds to an 1 Ampere excitation
⌢
j str = Xstr of the curl-curl matrix, [111].
Lemma 3.7 (Inductance extraction). If the stranded conductor coupling vector is diver-
gence free, i.e., SXstr = 0 and the curl-curl matrix kν is symmetric positive semi-definite,
then the extracted inductance
a) := X⊤str k+
Lstr (⌢ ⌢
ν ( a)Xstr (3.14)
Proof. We can add a grad-div regularization, without altering the inductance matrix, i.e.,
Lstr (·) = X⊤str k+
ν (·)X str = X⊤
k+
str ν (·)X str + X⊤ ⊤
str S SX str = X⊤
str k+
ν (·) + S⊤
S Xstr .
The inner matrix is symmetric positive definite, see Remark 4.1 and thus the nullspace of
this matrix is given by Xstr . This matrix has full column rank, thus we conclude symmetric
positive definiteness.
This reveals that a stranded conductor model corresponds to the series connection of a
lumped (nonlinear) inductance Lstr and a resistance Rstr .
25
3 Multiscale Device Models
Lemma 3.8 (Structural equivalence of conductor models, [117]). The solid and stranded
conductor models are structural equivalent, i.e., the stranded conductor can be written as
a solid conductor with a particular conductivity matrix.
Proof. To show the structural equivalence of the solid and stranded model, [117], we left-
multiply (3.12c) by Xstr R−1
str and add the result to (3.12a)
d⌢
Mσ̃ +Xstr R−1 ⊤
str Xstr a) = Mσ Xsol vsol + Xstr R−1
a + kν (⌢ ⊤ + −1
str Xstr Mσ,str Xstr Rstr vstr . (3.15a)
| {z } dt | {z }| {z }
=:Mequiv =:Mequiv =:Xequiv
R−1 ⊤ + ⌢
str Xstr Mσ,str kν ( a) = istr . (3.15b)
| {z }
=X⊤
equiv
This is the same structure as for solid conductors only, i.e., (3.10a)-(3.10b). Note that the
new (artificial) conductivity matrix Mequiv contributes with only one (positive) eigenvalue
for each stranded conductor to the global conductivity matrix.
Similarly other conductor models can be rewritten in the form of a solid conductor with
special conductivity matrix, e.g. the foil model, [39]. In the following chapters we will
chose the formulation that is more convenient for the corresponding section. Due to the
equivalence above the analysis and methods remain applicable for all conductor models.
For example the following field/circuit coupled system Section 3.4 is given in the solid
conductor formulation.
26
3.3 Semiconductor Device Model
where the refined electric flux density ρ is given in dependence of the elementary charge q,
the electron and hole densities n, p and the doping concentration C(~r). The densities are
determined by the additional conservation laws
∂n
− q −1 ∇ · J~n = −R, J~n = µn (UT ∇n − n∇ϕ), (3.17a)
∂t
∂p
+ q −1 ∇ · J~p = −R, J~p = −µp (UT ∇p + p∇ϕ), (3.17b)
∂t
where R = R(n, p) denotes the generation-recombination term, µn , µp are the mobility
parameters and UT is the thermal voltage. The total current leaving the device at terminal
(D) (D)
Γk ⊂ Γdir (k = 1, 2) is then given by
Z ∂
iD = ~ ~ ~ ~
Jdisp − (Jn + Jp ) · dA where J~disp = ~~ε ∇ϕ (3.17c)
(D)
Γk ∂t
d
Mn (Φ) n + Kn (Φ)n = rn (p, Φ), S̃Mε S̃⊤ Φ = q(n, p, vD ), (3.18a)
dt
d d
Mp (Φ) p + Kp (Φ)p = rp (n, Φ), iD = jD (n, p, Φ, Φ), (3.18b)
dt dt
with (regular) matrix functions Mn , Mp , Kn , Kp and the Laplacian S̃Mε S̃⊤ . The bold
symbols represent the vectors containing the discrete approximations of the corresponding
continuous quantities in (3.17). The lumped current iD is the discrete approximation of
the total current and it can be obtained in a post-processing step, i.e., it is not necessary
to solve a differential equation for Φ. Finally vD denotes the the voltage drop applied
to the device. It is determined by the surrounding circuit. The boundary conditions are
incorporated in the functions rn , rp and q.
The displacement current in (3.18b) can be expressed in terms of the time derivative of
the applied voltage drop vD , [2]. This capacitance extraction gives the following equivalent
definition of the coupling current
d
iD = CD vD − iSD with iSD := jSD (n, p, Φ). (3.18c)
dt
The capacitance CD may either by extracted from the discrete model or computed ana-
lytically, e.g. CD = εAD /lD for a cubic diode with isotropic permittivity ε, length lD and
cross-section AD .
This reveals that the semiconductor model corresponds to the parallel connection of a
lumped capacitance and a lumped (nonlinear) resistance.
27
3 Multiscale Device Models
Remark 3.3 (Singular mass matrices). The spatial discretization (e.g. using mixed finite
elements) of the semiconductor (3.18) might yield singular mass matrices Mn , Mp . This
turns n, p in some discretization nodes into algebraic variables, [26]. However, this does
not affect the following analysis; projectors are defined similarly to Definition 2.8. This
particular case is disregarded in the following.
3.4 Conclusions
At the conclusion of this chapter we assemble the equations of the spatial discrete models
into one system of equations: the extended circuit model (3.2), (3.1b)-(3.1e), the MQS
device model (3.10a)-(3.10b) and the drift-diffusion model (3.18). The coupled PDAE sys-
tem is given by the extended network equations (for simplicity of notation in the notation
of traditional MNA)
d
AC qC (A⊤C u, t) + AR gR (A⊤R u, t) + AL iL + AV iV + AI is (t)
dt (3.19a)
d
+ AM X⊤M kν (⌢
a)⌢
a + AD jD (n, p, Φ, Φ) = 0
dt
d
φ (iL , t) − A⊤L u = 0 (3.19b)
dt L
A⊤V u − vs (t) = 0 (3.19c)
d⌢
Mσ a + kν (⌢ a − Mσ XM A⊤M u = 0 (3.19d)
a)⌢
dt
d
Mn (Φ) n + Kn (Φ)n − rn (p, Φ) = 0 (3.19e)
dt
d
Mp (Φ) p + Kp (Φ)p − rp (n, Φ) = 0 (3.19f)
dt
S̃Mε S̃⊤ Φ − q(n, p, A⊤D u) = 0 (3.19g)
with the unknown node potentials u : I → Rn without ground and the currents iL :
I → RnL , iV : I → RnV through inductors and voltage sources (where nL and nV denote
their number), the magnetic vector potential ⌢ a : I → Rn1 for the MQS device and the
electric scalar potential Φ : I → Rn0 , the electron and hole densities n, p : I → Rn0
for semiconductors. In the spirit of traditional MNA, the device currents iM and iD in
Kirchhoff’s current law (3.2) are replaced by the corresponding current assignments (3.10b)
and (3.18b).
The following chapter discusses the expected numerical difficulties during the time-
integration of problem (3.19a)-(3.19d) with AD = [] in terms of the DAE-index. Later
on, in Section 5 efficient multirate strategies for the solution in the time domain are dis-
cussed.
28
4 DAE-Index Analysis
Differential/algebraic equations are structurally different from ordinary differential equa-
tions. The additional algebraic constraints make those problems ‘infinitely’ stiff, so only
implicit time-integration methods can be applied. Furthermore the errors of numerical ap-
proximations are more critical. To classify this criticality, the DAE-index was introduced,
e.g. [25, 68]. There are different concepts, but roughly speaking the idea is that the index
corresponds to the highest derivative of an input function that enters the problem. This
highest derivative affects the numerical approximation most severely. In the index-0 case,
which corresponds to an ordinary differential equation, there is no such derivative. On the
other hand systems of index larger than 1 are called higher index problems and suffer from
higher derivatives. They may require special numerical treatment.
In this chapter we will analyze the DAE-index of the coupled system (3.19). It consists
of contributions from three different subproblems. The index results for circuits containing
semiconductor devices based on the drift-diffusion equation are well-known, e.g. [13, 120,
127]. Thus this section focuses on the index-analysis of the field/circuit coupled problem.
The extension to all three problems is straightforward when using topological assumptions
that keep the circuit’s branches of both PDE problems separate.
In the following we employ the tractability index. This concept is often used in the
circuit analysis community. It gives a detailed view of the structure of the equations by
using projectors and especially for the MNA the tractability index is only determined
by the circuit’s topology, [54]. However, we expect for the field/circuit problem the same
results when employing other index concepts.
29
4 DAE-Index Analysis
assumed continuous in their arguments and smooth. Their partial derivatives are
∂ ∂
d′x (x, t) := d (x, t) and b′x (x, t) := b (x, t) .
∂x ∂x
The unknown solution is x = x(t) ∈ D ⊂ Rl , t ∈ I = [t0 , te ].
Definition 4.1 (Properly stated leading term, [92]). The DAE (4.1) has a properly stated
leading term if and only if
and if there is a representing projector R ∈ C 1 (I, Rn ), Ker A = Ker R (t), Im d′x (x, t) =
Im R (t) and d (x, t) = R (t) d (x, t) for all x ∈ D and t ∈ I.
The following definitions are used to discuss the index, see [91].
Definition 4.2 (Matrix chain and subspaces). Given the DAE (4.1), we define recursively
the following objects:
Definition 4.3 (Tractability index, [92]). The DAE (4.1) with a properly stated leading
term is called DAE of (tractability) index-0 if and only if
30
4.2 Index Analysis for the Field/Circuit Problem
linear 3D model with an attached network. Now following [8], the nonlinear field/circuit
problem in 3D is studied using the tractability index concept, [63]. This allows for a deeper
understanding of the index-2 case, for example compared to [117]. For the tractability
analysis the coupled system without semiconductors (AD = []), i.e., equations (3.19a-
3.19d), is given in the notation of traditional MNA with a properly stated leading term,
i.e., in the form of (4.1)
AC 0 0 +
0 I 0 d AC AC qC (·)
φL (·)
0 0 0 dt
M+σ Mσ ⌢a
0 0 Mσ
(4.2)
AR gR (·) + AL iL + AV iV + AM X⊤ kν (·)⌢ a + AI is (t)
−A⊤L u
+ = 0.
A⊤V u − vs (t)
⌢ ⊤
kν (·)a − Mσ XAM u
where the unknown vector is given by x⊤ = u⊤ i⊤L i⊤V ⌢ a⊤ . Let Qσ and Pσ = I − Qσ
be defined as in Definition 2.8, i.e., Qσ is the constant projector onto Ker Mσ , such that
Ass. 2.13 ensures
Ker Mσ + Q⊤σ Kν (·)Qσ = {0}
Compared with the circuit only case, R features an additional row and column for Pσ ,
because of the additional curl-curl equation. The curl-curl equation contributes also with
the conductivity matrix Mσ to the overall mass matrix
AC C (·) A⊤C 0 0 0
0 L (·) 0 0
G0 (x, t) := Ad′x (x, t) = . (4.3)
0 0 0 0
0 0 0 Mσ
If the mass matrix above is regular, all equations are differential equations, such that the
problem is an ODE. This is the case for the following class of circuits, [8].
Theorem 4.1 (Index-0). Let Ass. 2.13, 3.1 and 3.2 be fulfilled. Then the DAE (4.2) has
index-0 if and only if there is no voltage source, a tree containing capacitors only and
31
4 DAE-Index Analysis
Proof. The matrix G0 is nonsingular iff all blocks on the diagonal have full rank. The
matrices C and L are positive definite by assumption and Ker A⊤C = {0} is trivial iff
the circuit has a tree containing capacitors only, see Remark 3.1. The third row/column
vanishes iff there are no voltage sources. Finally, Ker Mσ = {0} iff the circuit does not
include MQS devices (this is obviously the classical MNA case, [54]) or the domains of all
MQS device are conductive (Mσ has full rank).
and
AR G (·) A⊤R QC 0 AV AM X⊤ Kν (·)Qσ
−A⊤L QC 0 0 0
b′x (x, t) Q0 = ,
A⊤V QC 0 0 0
−Mσ XA⊤M QC 0 0 Kν (·)Qσ
where Q0 is a constant projector onto Ker G0 (x, t). Again, the matrices Q0 and b′x (x, t) Q0
are similar to the well-known 3×3 block structure in MNA index analysis. Thus the index-1
proof as given in [8] is a straightforward extension of [54].
Theorem 4.2 (Index-1). Let Ass. 2.13, 3.1 and 3.2 be fulfilled and the circuit contains at
least an MQS device, a voltage source or there is no tree containing capacitors only. Then
the DAE (4.2) has index-1 if and only if there is neither
(a) a LIM-cutset, i.e., a cutset consisting of inductances, current sources and MQS
devices only, nor
(b) a CV -loop, i.e., a loop consisting of capacitances and at least a voltage source only.
Proof. The intersection (N0 ∩ S0 ) (x, t) is analyzed in the following, see Definition 4.2.
Let W0 (x, t) denote a projector along Im G0 (x, t). Consequently W⊤0 (x, t) is a projector
onto Ker G⊤0 (x, t). Symmetry gives Ker G⊤0 (x, t) = Ker G0 (x, t) and thus we can choose
W⊤0 (x, t) = Q0 . We find starting from Definition 4.2
with
32
4.2 Index Analysis for the Field/Circuit Problem
and consequently
We show that z = 0 when z ∈ (N0 ∩ S0 ) (x, t). Let z⊤ = [z⊤1 z⊤2 z⊤3 z⊤4 ], then follows from
Q0 z = z
QC z1 = z1 , (4.5)
Qσ z4 = z4 , (4.6)
z2 = 0, (4.7)
where Mσ + Kν (·) is positive definite by Ass. 2.13 and thus we achieve z4 = 0 using (4.6).
Now we are back to the ‘classical’ case, i.e., without MQS device. From (4.5-4.10) it follows
Left-multiplication of (4.11) by z⊤1 and using (4.12) gives z1 ∈ Ker A⊤R QC and thus
Q⊤C AV z3 = 0. Now, from (4.13) follows
Then (N0 ∩ S0 ) (x, t) = {0} holds iff there are neither LIM-cutsets nor CV -loops with at
33
4 DAE-Index Analysis
The intersection (4.4) can be given elegantly in terms of additional projectors from the
classical MNA index-analysis, [54]: the constant projector QCRV = QC QV-C QR-CV onto
Ker[AC AR AV ]⊤ , where QC-V , QV-C and QR-CV are constant projectors onto Ker Q⊤C AV ,
Ker A⊤V QC and Ker A⊤R QC Q⊤C-V From those projectors this result follows immediately:
is constant.
This result will be exploited in the next section in Theorem 4.5, where we prove that the
DAE-index is at most 2.
Remark 4.1. Let K, K+ and Q denote an arbitrary matrix, its Moore-Penrose pseudo-
inverse and a projector, respectively. If K = Q⊤ KQ then it follows that K+ = QK+ Q⊤ .
With Remark 4.1 the following lemma from [8] can be proved, which is structurally
similar to the extraction of a lumped inductance from the PDE model, see Lemma 3.7.
Lemma 4.4 (Consistent excitation). Let Ass. 2.13 and 3.4 be fulfilled, then
+
⊤ + ⊤
X Hk (·) X with Hk (·) = Kν (·) Kν (·) − Qσ Kν (·)Qσ Kν (·)
is positive definite.
with a curl-curl matrix Zk (·) := Kν (·) + M⊤ S̃⊤ S̃M that is fully regularized and the block-
elimination T (·) := I − (P⊤σ Zk (·) Qσ )(Q⊤σ Zk (·) Qσ )+ .
34
4.2 Index Analysis for the Field/Circuit Problem
The matrix T (·) is regular with T−1 (·) = I + (P⊤σ Zk (·) Qσ )(Q⊤σ Zk (·) Qσ )+ and X has full
column rank by construction. Thus the definiteness of Zk (·) must be shown, where the
only interesting elements are from Ker C = im S̃⊤ , i.e.,
x⊤ S̃ Zσ + M⊤ S̃⊤ S̃M S̃⊤ x > 0 for all x 6= 0
1. if S̃⊤ x ∈ Ker M then there is an y such that S̃⊤ x = Qσ y. Thus the second summand
vanishes and the first summand is positive because of Ass. 2.13.
2. else M S̃⊤ x 6= 0 and thus the second summand is positive (because its kernel is
Ker M S̃⊤ ) and the first summand is non-negative.
Remark 4.2 (Schur complement). The matrix Hk in Lemma 4.4 corresponds to a Schur-
complement. Let us assume a convenient partitioning of the curl-curl-equations into equa-
tions for conducting and nonconductive domains Mσ = diag(M11 , 0). Then the projector
onto Ker Mσ can be given as Qσ = diag(0, I). With a corresponding block partitioning of
the curl-curl matrix
K11 (·) K12 (·) H11 0
Kν (·) = follows Hk (·) = . (4.16)
K⊤12 (·) K22 (·) 0 0
with H11 := K12 (·)K+22 (·)K⊤12 (·). The Schur complement is the simplest form of domain
substructuring, [104]. It yields a reduced system because only the first block must be solved.
This can be exploited when solving linear systems, see Section 5.4.
Now, having computed the important ingredient Lemma 4.4, the next element of the
matrix chain Definition 4.2, i.e., G1 (x, t) = G0 (x, t) + b′x (x, t) Q0 is analyzed to extract
the topological conditions for the index-2 case.
AC C (·) A⊤C + AR G (·) A⊤R QC 0 AV AM X⊤ Kν (·)Qσ
−A⊤L QC L (·) 0 0
G1 (x, t) = ⊤
.
AV QC 0 0 0
⊤
−Mσ XAM QC 0 0 Mσ + Kν (·)Qσ
Following [8], the next proof will not be based on the classical procedure as given in [54],
where the projector Q1 onto Ker G1 (x, t) is computed. Instead it will be shown that the
intersection N1 ∩ S1 is trivial, see Definition 4.2 in Section 4.1.
Theorem 4.5 (Index-2). Let Ass. 2.13, 3.1, 3.2, and 3.4 be fulfilled and let the circuit con-
tain at least one MQS device or one voltage source or there is no tree containing capacitors
only. Then the DAE (4.2) has index-2 if and only if there is either
(a) a LIM-cutset, i.e., a cutset consisting of inductances, current sources and MQS
devices only, or
(b) a CV -loop, i.e., a loop consisting of capacitances and at least a voltage source only.
35
4 DAE-Index Analysis
is defined and it holds true Im G1 (x, t) ⊂ Ker W1 (x, t). Now, we can reformulate the
object
with Hk as defined in Lemma 4.4. The matrix imposes the following conditions on z (this
is now almost the same procedure as in the proof of Theorem 4.2):
The left-multiplication of equation (4.22) by (Qσ z4 )⊤ together with the Gauging Ass. 2.13
(Mσ + Kν (·) is positive definite) gives
Qσ z4 = 0, i.e., z4 = P σ z4 . (4.23)
z4 = XA⊤M QC z1 . (4.24)
36
4.2 Index Analysis for the Field/Circuit Problem
Still using standard techniques, equations (4.20), (4.25) and the definition of QC yield
and thus QC z1 = QCRV QC z1 . Now, substituting (4.20) and (4.24) in (4.17) gives
because HC (·) = AC C (·) A⊤C + Q⊤C QC is positive definite. Multiplying of (4.19) from
the left by Q⊤C leads to Q⊤C AV z3 = 0 and z3 ∈ Im QC-V respectively. Finally (4.18) and
z3 ∈ Im QC-V give
37
4 DAE-Index Analysis
onto (N0 ∩ S0 ) (x, t) and the complimentary projector U := I − T. The stiffness term can
be split accordingly b (x, t) = b (Ux, t) + BTx using a problem-specific matrix B and
d (x, t) = d (Ux, t). Thus it can be concluded that the index-2 variables enter our system
linearly.
4.3 Conclusions
In this section the structural properties of magnetoquasistatic devices in electrical circuits
modeled by MNA were discussed. The coupled multiscale system of lumped devices (re-
sistors, inductors, capacitors, independent current and voltage sources) and MQS devices
was formulated with a proper leading term and analyzed by the tractability index concept.
The field/circuit coupled problem was proved to be numerically harmless, i.e., it is
index-2 (with linear index-2 components) at most and index-1 under rather mild conditions
(analogously to the classical network case). The MQS devices were plugged into the circuit
as controlled current sources, but the analysis shows that they behave topologically as
inductances. This corresponds to the physical effects covered by the eddy current problem.
The index results are numerically verified in Section 6.1 by applying Euler’s method as
a time-integrator to a PDE inductance model. More efficient (multirate) time-integration
methods are discussed in the following section.
38
5 Multirate Methods
As we mentioned in the introduction, the coupled multiscale problems of Section 3 are
a challenging task for time-integrators: the intrinsic time rates of subsystems differ by
several orders of magnitude and the type of the subsystems is different, i.e., concerning
their symmetry and definiteness. Unfortunately the structural properties can also change.
The coupling of DAE index-1 subsystems can yield an arbitrarily high index problem, e.g.
[51].
The introductory example Figure 5.1a shows an electric circuit configuration, where
a part is refined by a magnetoquasistatic field model using partial differential equations
(PDEs). Due to the different modeling techniques, i.e., network and space discretization,
the underlying equations are different in shape, the circuit system is non-symmetric and
typically tackled by direct solvers, while the field is symmetric and solved by iterative
methods. Furthermore the voltages, Figure 5.1b, are only pulsed in the circuit part but
slow and sinusoidal in the field part. The discretization has to resolve the dynamics of the
coupled system as a whole and thus it produces a series of time steps that matches the
dynamics of the most active component (i.e., the one working at the highest frequency).
Due to switches, filters or high integration there may only be a small number of devices
active at any given moment, while the others remain latent. The time-integrator will
resolve those parts with an unnecessarily high time resolution causing an avoidable high
computational cost.
Standard single-rate time-integration (described in Section 5.1) is inefficient for those
problems. A possible work-around is presented in Section 5.2: a Schur complement ap-
proach for the MQS device is introduced, cf. [134, 58]. It allows us to use different linear
solvers for the subsystems, i.e. a direct solver for the circuit and an iterative one for the
field. Furthermore the Schur complement is used for several time steps (bypassing), such
that multiple time scales are exploited to some extend.
A more general way is to treat every subsystem independently within a cosimulation
scheme, Section 5.3. This permits the use of multirate and multimethod techniques and
600
field model u1
u3
Lstr 400
voltage [V]
u1 u3 -400
0 0.005 0.01 0.015 0.02 0.025 0.03
time [s]
39
5 Multirate Methods
simplifies the coupling of simulator packages. The drawback is the decline in stability:
the presence of algebraic constraints can handicap higher-order integration, [136] and may
cause divergence, [1, 5].
In this thesis stability and convergence are proved for an iterative cosimulation scheme
(‘dynamic iteration’) by analyzing the coupling interface, see Section 5.3. The analysis in-
cludes particularly the application in field/circuit and semiconductor/circuit problems as
discussed in Section 5.3.5 and Section 5.3.6, respectively. In Section 5.3.5 the cosimulation
approach is adapted for efficient multirate time-integration of field/circuit coupled prob-
lems by using reduced order models: ‘fitting parameters on the fly’, [118]. The conditions
for stability and convergence are obtained and the convergence order is well understood.
using coefficients αi (k-th order BDF). For the special choice k = 1 with α0 = −α1 = 1,
the method is the implicit Euler scheme. The nonlinear time-discrete problems are usually
40
5.2 Multirate Bypassing of MQS Schur Complements
solved by Newton-Raphson
α0
J(i)
n xn
(i+1)
= −F(i) + J(i) x(i) with J(i)
n := G0 x(i) ′ (i)
n , tn + bx xn , tn
| n {z n n} h (5.2)
=:r
F(i)
n := F(xn−k , . . . , xn−1 , x(i)
n , tn−k , . . . , tn )
where G0 := Ad′x (x, t) and b′x are the (differential) mass and stiffness matrices, respec-
tively (see Definition 4.2). The problems in this treatise exhibit only constant mass ma-
trices, such that G0 = const and thus the time-derivative is discretized directly.
In Section 2.2.3 the regularization of the curl-curl operator was discussed, such that the
(i)
matrix pencil Jn above is invertible, cf. equation (2.25). This guarantees well-posed linear
problems for direct solvers in the Newton-Raphson iteration (5.2). Alternatively iterative
solvers could be used, which would benefit from their weak gauging property as noted in
Remark 2.1, [37].
The single rate approach is not efficient for problems with different time-scales. Some
components of the unknown vector may behave slowly and almost linearly in time, while
other are fast and nonlinear. This limitation can be mitigated by using bypassing tech-
niques, Section 5.2 or completely overcome by cosimulation Section 5.3.
with element-wise assembly matrices Qe that organize the contributions such that the
1
Element refers here to circuit elements and they are not related to elements used for space discretization.
41
5 Multirate Methods
MNA equations (3.19) are obtained, [116]. In particular they incorporate the incidence
(i) (i)
matrices Ae . One calls the tuple Je , Fe the element stamp. It consists of internal and
external variables, i.e., variables used only inside the particular element and variables that
are related to other elements by the simulator, [56, 58].
Following [116, 134] we want to speed up the solution of the Newton system in
field/circuit applications by eliminating the magnetic vector potential ⌢ a. We focus here on
the MQS device’s contribution to the time-discretization of system (3.19). This is in the
notation above a particular contribution to the function F in (5.3). It reads
I 0 0 0 0 0
(i) (i) 1 (i)
FM := RM −I 0 xM + 0 0 X⊤M ρxM ,
h
−XM 0 kν (⌢ a (i) ) 0 0 Mσ̃
The first row contains the contribution to the current balance equation, the second line
is the coupling equation and finally the last row represents the MQS curl-curl equation.
Here the model is excited by stranded conductors (the solid conductor case is analogous).
We obtain the following Jacobian contribution:
I 0 0
(i) (i) α0
JM := RM −I αh0 X⊤M with Kh := Kν (⌢ a (i) ) + Mσ̃ (5.4)
(i) h
−XM 0 Kh
where the only unknown ⌢ a is the magnetic vector potential, which is internal, i.e., it is not
used outside the MQS stamp. Only the current/voltage relation of the series connection
of a (nonlinear) inductor and a resistor needs to be revealed to the host circuit simulator.
This is the case if ⌢
a is eliminated from the Newton system by the Schur complement.
This is beneficial for all kinds of large elements, e.g. it allows cache-optimized stamping
of semiconductor models, [58], or field devices, [134, 116]. In either case more compact
stamps are obtained, which fit better into the overall MNA framework.
The unknown ⌢ a is removed and one ends up with a reduced stamp in terms of x̃⊤M =
⊤ ⊤
(iM , vM ). The corresponding reduced Jacobian reads
" #
I 0 0 I 0
(i) (i) I 0
J̃M :=
(i)
−1 JM 0 I =
(i) (5.6)
0 I − αh0 X⊤M Kh RM + αh0 Lh −I
0 0
42
5.2 Multirate Bypassing of MQS Schur Complements
(i)
using Kh from (5.4). For static models, equation (5.7) extracts an inductance by one
(i)
Ampère excitation, [118], but here Lh also takes eddy current effects into account (due
to the presence of the conductance matrix Mσ̃ ). Thus the inductance depends on the
frequency (via the step size h) and therefore the matrix must be recomputed or interpolated
for any change of h.
In the following the computational costs are discussed in terms of direct solvers (matrix
factorizations). Their application is beneficial for our approach, because the factorization
can be stored. For the Schur complement in the Newton iteration i+1, we need to compute
(i) (i)
Lh . When applying a direct solver, the matrix Kh is factorized (one LU decomposition)
and afterwards forward/backward substitutions are carried out for each circuit branch
(m = 1, . . . , nM ):
(i) (i) (i) ⌢(i)
a(i)
Kh ⌢m = Xm s.t. Lh = X⊤M ⌢
aM with aM := [⌢
a1, . . . , ⌢
a nM ] (5.8)
by sparse inner products. Also, the magnetic vector potential for the right-hand-side
voltage must be computed. To this end, we solve the following equation inside each Newton
iteration (derived from Jacobian (5.4) and right-hand side (5.5)):
43
5 Multirate Methods
200 0.06
pwm
150 sine
0.05
100
0.04
inductance [H]
50
voltage [V]
0 0.03
−50
0.02
−100
0.01
−150
−200 0
0 0.02 0.04 0.06 0.08 0 0.002 0.004 0.006 0.008 0.01
time [s] time [s]
(a) fast signals applied to MQS device (b) slow change of the inductance
Figure 5.2: Time scales in a field/circuit coupled problem: the inductivity changes according
to the slow sine wave, which is the low frequency part of the pulsed signal.
This improves the result of [134] for a similar setting. It has been shown that only one
factorization and nM + 1 forward/backward substitutions, i.e., equations (5.8) and (5.9),
are necessary for the Schur complement in each Newton iteration. The choice of the
linear solver for the Schur complement is independent of the solver used in the circuit
host simulator (typically a direct solver). So, for example an iterative method such as
(block) PCG could be used for the Schur complement, [110]. This allows us to solve 3D
problems within the circuit simulator without additional gauging, because of the weak
gauging property, Remark 2.1. This highly improves the efficiency of the linear solver,
especially if multiple right-hand-sides are supported to solve (5.8-5.9), e.g., [142].
with the initial energy level of the device E0 . In the i-th Newton iteration of time step n
the energy is approximated by
n−1
!
X (i) (i)
En ≈ E0 + h iM,j vM,j + iM,j vM,j
j=0
and this can easily be compared to the initial energy level E0 using a (relative) norm.
Updates of the nonlinearity are only necessary if the energy level changes significantly.
Consequently, the computation of the material matrices can often be bypassed, Fig. 5.3.
Then the model behaves (nearly) linearly and only one forward/backward substitution
44
5.3 Multirate Cosimulation
circuit
t0 h te
Figure 5.3: Bypassing. The field and circuit problems are decoupled: the solid arrows denote
time steps where the system is solved by an ordinary Newton-Raphson iteration. Time steps
with the dashed arrows indicate the usage of the Schur complement.
for the right-hand-side per iteration is necessary (5.9). This defines a simplified Newton
algorithm, where the Jacobian (5.6) is frozen for several iterations and possibly several
time steps if the (relative) change of energy does not exceed a threshold and thus the
reluctivity is (nearly) constant.
Furthermore, if the problem is rather latent the right-hand-side evaluation can be by-
passed as well. The vector potential needs no update and thus the field problem is decou-
pled from the circuit, where it is represented by an inductance matrix. This is obviously
not free of risks, because the bypassing of the right-hand-side changes the fixed point of
the Newton scheme. The algorithm is given in Listing 1 where the energy level and the
reluctivity ν are monitored by relative norms. This algorithm relieves the host simulator of
the burden of solving unnecessarily large system of equations, especially if the nonlinearity
is weak. On the other hand the host-simulator still can use a Jacobian for its Newton
iteration. The drawback are the additional iterations due to the inferior convergence of
simplified Newton, albeit solving a sequence of reduced systems.
5.2.2 Conclusions
The Schur complement approach yields significantly smaller element stamps that are equiv-
alent to a series connection of an inductance and a resistance (this is in line with DAE-index
result, Section 4). The additional costs of the complement computation can be disregarded
if linear solvers with multiple right-hand side techniques are available. In particular iter-
ative solvers will improve the efficiency of this solution process. Due to bypassing, the
field/circuit coupling is weakened and thus the time-integration of the circuit is cheapened
because only basic elements are evaluated. This decoupling exploits the multirate time
behavior of the coupled system if present, see Section 6.2.
Although the bypassing allows for a decoupling, the circuit simulator still controls the
simulation: the same time-integrator is used for all subproblems, the circuit simulator
decides about step sizes and Newton iterations etc., see Fig. 5.3. If the level of decoupling
is to be further increased the use of waveform relaxation approaches is appropriate. This
is the topic of the next section.
45
5 Multirate Methods
that is able to mitigate some of the most severe problems (multirate and the efficency of
linear solvers). On the other hand, problems often come along with their own specific
simulation packages and thus the coupled system cannot be solved monolithically (as a
whole system) in the time domain. This renders standard time-integration impossible
and even the bypassing approach is only feasible if the simulation packages supply the
corresponding interfaces, e.g., such that we can access the magnetic vector potential to
ensure correct initial values.
One way to overcome this impasse is to apply cosimulation methods, i.e. methods for
the coupling of simulator packages. Each model is simulated separately by its package,
e.g. network and PDE models. Thus all problems may be solved on their own time scale,
with tailor-made methods (multimethod). Information on how the models interact are
exchanged only at synchronization points (on ‘time windows’).
In the following we will revisit the waveform relaxation or dynamic iteration schemes,
which solve the subproblems iteratively and exchange coupling information in each sweep,
[88, 141]. Fig. 5.4 depicts the general idea of separate time-stepping, windowing and it-
eration. Those schemes are known to be unconditionally stable for coupled ODEs, [97,
27] and their convergence is well understood for ODEs that stem from circuit simulation
and space discretized PDEs, [135, 77]. On the other hand instabilities are known when
those schemes are applied to DAEs. The introduction of a windowing technique has been
proposed to speed-up convergence and avoid instabilities. Nonetheless the contraction of
the underlying fixed point iteration is only guaranteed if a stability constraint is fulfilled.
This constraint forces the algebraic coupling to be weak, but then convergence on a single
window is enforced, [76, 7].
The error transport for multiple windows has been analyzed so far for a special class
of DAEs, where the coupling is established by Lagrangian multipliers, [5]. In [51] the
error transport of the general case was considered using a differential equation for the
error propagation, but the DAE system was reduced to its underlying ODE. This has
simplified several aspects, that are discussed in [9], which is the basis of the following
analysis. It generalizes the approach of [5] to the most general form of index-1 DAEs,
for which a similar stability constraint is derived. Global convergence and stability (with
error propagation) will be guaranteed if the splitting error remains in a neighborhood of the
analytical solution. Sections 5.3.5 and 5.3.6 will apply the theory to the electromagnetic
46
5.3 Multirate Cosimulation
field
t0 H1 t
circuit
Figure 5.4: Dynamic iteration. The subsystems (e.g. field and circuit) are discretized on their
own time scales using different time steps (solid arrows). Several time steps are a time window
and they are computed iteratively (the ‘sweeps’ are denoted by the dashed arrows).
models in electric networks of the previous sections and show that dynamic iteration is
indeed a strategy for exploiting multirate behavior using reduced order models similar to
those illustrated in [106].
where the variables are renamed to x⊤ = [y⊤ , z⊤ ] with the vector functions f and g. This
formulation addresses the whole problem abstractly and thus allows us to extract the un-
derlying principles without the dispensable details. The following analysis of system (5.11)
is limited to the index-1 case. In the field/circuit case this can be assured by Theorem 4.2
that demands loop/cutset conditions for the circuit. This prerequisite is mathematically
formalized by the following index-1 assumption:
Assumption 5.1 (Monolithic index-1). The differential algebraic initial value prob-
lem (5.11)
and the initial values are particularly consistent, i.e., y0 and z0 solve the algebraic
equation (5.11b),
(b) the right-hand-side functions f and g are supposed to be sufficiently often differen-
tiable in the neighborhood of the solution,
(c) the Jacobian ∂g/∂z is non-singular in the neighborhood of the solution (the coupled
problem is index-1).
47
5 Multirate Methods
Remark 5.1 (Index concepts for Hessenberg systems). Please note that we do not dis-
tinguish between different DAE index-concepts for semi-explicit systems. Equation (5.11)
describes the special case of an index-1 Hessenberg system and for Hessenberg systems
of index-1, index-2 and index-3 it has been shown in [67, Page 13] that the important
index-concepts coincide, [29].
In the index-1 case above the vector y contains only the differential variables, i.e., the
variables that are defined by first derivatives with respect to time, while z contains the
algebraic variables that are not described by any derivatives.
Let the coupled problem (5.11) consist of r subsystems. For example a network con-
taining an MQS and a semiconductor device consists of r = 3 subsystems. Let us further
assume that the equations are partitioned accordingly, then the i-th subsystem is given by
(for i = 1, . . . , r)
where the global right-hand-sides are assembled in the obvious way f⊤ = [f⊤ ⊤
1 , . . . , fr ] and
g⊤ = [g⊤1 , . . . , g⊤r ]. Similar to the monolithic Ass. 5.1 we require for each subsystem
Assumption 5.2 (Subsystem index-1). The Jacobian
where w1 (t) and w2 (t) are given inputs and each subsystem is index-1
I I
∂g1 /∂z1 = ∂g2 /∂z2 = I, but w1 = z2 , w2 = z1 yields ∂g/∂z =
I I
48
5.3 Multirate Cosimulation
where w1 (t) and w2 (t) are given inputs and each subsystem is index-2
0 I
∂g1 /∂z1 = ∂g2 /∂z2 = 0, but w1 = z2 , w2 = z1 yields ∂g/∂z =
I 0
being close to the exact solution (5.12). As indicated before, the dynamic iteration scheme
operates on time windows [Tn , Tn+1 ], such that
with window size Hn := Tn+1 − Tn . In the multirate context those windows are called
macro steps in contrast to the micro steps h of the numerical time-integration, see Fig. 5.4.
Assuming a numerical approximation is computed on the window [Tn−1 , Tn ], the dynamic
iteration defines a new approximation on the consecutive window
x̃|[Tn ,Tn+1 ] ∈ Cn1,0 with Cn1,0 := C 1 ([Tn , Tn+1 ], Rny ) × C([Tn , Tn+1 ], Rnz )
using an extrapolation step followed by (one or more) iteration steps. The corresponding
steps are defined by operators as in [5].
1,0
Extrapolation step. Let the operator Φn : Cn−1 → Cn1,0 denote a continuous extrapola-
tion from the old window [Tn−1 , Tn ] to the new window [Tn , Tn+1 ]. This defines an initial
guess of the new approximation
" # " # " #
(0)
ỹn ỹ|[Tn−1 ,Tn ] Φy,n
(0) := Φn with Φn = . (5.16)
z̃n z̃|[Tn−1 ,Tn ] Φz,n
(0)
Actually the initial value x̃n (Tn ) is fixed from the previous window. Its constant extrap-
olation on the new window is the most common choice for an initial guess:
" #
(0)
ỹn (t) ỹn (Tn )
(0) = for all t ∈ [Tn , Tn+1 ].
z̃n (t) z̃n (Tn )
49
5 Multirate Methods
This operator introduces an error in O(Hn ), which can be improved by linear or higher
order polynomial extrapolation. It is uniformly Lipschitz-continuous, independently of Hn ,
see [5].
Iteration step The extrapolation has defined an initial guess. It is followed by an iteration
step defined by the mapping Ψn : Cn1,0 → Cn1,0
" # " # " # " #
(k−1) (k) (k−1)
ỹn ỹn ỹn Ψy,n
(k−1) → (k) := Ψn (k−1) with Ψn = (5.17)
z̃n z̃n z̃n Ψz,n
where the differential and algebraic splitting functions F and G, respectively, are consis-
tent, i.e., they fulfill the following compatibility condition:
Definition 5.1 (Consistent splitting functions). Splitting functions F and G are called
consistent if they are sufficiently differentiable and fulfill the compatibility condition
F y, y, z, z = f y, z and G y, y, z, z = g y, z . (5.19)
Remark 5.2 (Fixed-point). The following statements follow immediately from the defini-
tions above:
• the compatibility condition implies the exact solution x is a fixed-point of the itera-
tion operator Ψn .
• for the partitioned system (5.13) with corresponding unknowns ỹ⊤n = [ỹ⊤1,n . . . ỹ⊤r,n ],
z̃⊤n = [z̃⊤1,n . . . z̃⊤r,n ], the iteration operator Ψn is defined by r initial-value problems:
(k) (k) (k−1)
ỹ˙ i,n = Fi (ỹn(k) , ỹn(k−1) , z̃(k) (k−1)
n , z̃n ), with ỹi,n (Tn ) = ỹi,n (Tn ),
(5.21)
0 = Gi (ỹn(k) , ỹn(k−1) , z̃(k) (k−1)
n , z̃n )
All common iteration schemes can be encoded by the splitting functions above, i.e., Pi-
card, Jacobi or Gauß-Seidel-type schemes, [9]. For example the Gauß-Seidel-type scheme
is visualized in Fig. 5.5 for the case of two subsystems. It computes sequentially the so-
lutions of all r subsystems. It iterates such that the i-th subsystem in the k-th iteration
(k)
solves for the new solution xi,n where it utilizes the latest available data for the variables
50
5.3 Multirate Cosimulation
PSfrag
(k−1) (k)
x̃1,n x̃1,n
Subsystem 1 Subsystem 2
(k−1) (k)
x̃2,n x̃2,n
Figure 5.5: Schematic representation of the k-th Gauß-Seidel iteration on the n-th time window
for r = 2 subsystems; schematic as given in [51].
of the previous subsystems (1, . . . , i − 1) and old data (from the previous iteration) for the
variables of following subsystems (i + 1, . . . , r).
Remark 5.3. The Gauß-Seidel-type iteration scheme of Definition 5.2 defines consistent
splitting functions F and G, i.e., they fulfill the compatibility condition and inherit their
smoothness and differentiability from f and g, see Definition 5.1. Please note that the split-
ting functions reflect changes in the computational sequence. In other words: it matters
in which sequence the subsystems (5.21) are solved.
For coupled ODEs (short ‘ODE-ODE’) various splitting schemes can be shown to be
convergent, [27] and similar results are known for the special case of an ODE coupled to
an algebraic equation (short ‘ODE-AE’), [136].
(i) The Gauß-Seidel type scheme for r = 2 ODE subsystems has the following form
(k) (k) (k−1)
" #
ỹ˙ 1 = f1 (ỹ1 , ỹ2 ), (k) (k−1)
f1 (ỹ1 , ỹ2 )
(k) (k) (k)
with F= (k) (k) .
ỹ˙ 2 = f2 (ỹ1 , ỹ2 ), f2 (ỹ1 , ỹ2 )
(ii) The Gauß-Seidel type scheme for an ODE and AE subsystem, i.e., the ‘fractional
51
5 Multirate Methods
(k−1)
Please note that the old waveforms ỹ2 and z̃(k−1) enter in both examples only in the
first differential equation. The importance of that fact becomes clear by Cor. 5.6.
In contrast to the special cases above, the general DAE case is more involved: the
algebraic constraints (5.22b) may depend on old algebraic variables, i.e., those of the
(k−1)
previous iteration z̃n . They can cause divergence and thus the scheme is carefully
analyzed in the following: Section 5.3.3 gives a fixed point analysis in function space of
the iteration scheme and Section 5.3.4 carries out an error analysis (using the fixed point
argument) to prove stability and convergence.
The reader who is not interested in mathematical details and proofs may skip to Sec-
tion 5.3.5 or Section 5.3.6 for applications in field/circuit and semiconductor/circuit cou-
pling, respectively.
where kv(t)k2,∞ := maxt kv(t)k2 utilizes the maximum-norm in time and the Euclidean
norm in space.
The splitting functions must fulfill the following smoothness properties on that neighbor-
hood, which are typically inherited from the original right-hand-side functions f and g:
Assumption 5.3 (Smoothness in the neighborhood of the solution). Let the problem (5.11)
with consistent splitting functions F, G be given, then it is assumed that there is a d0 > 0
with
(a) the (differential) splitting Function F is Lipschitz-continuous on Ud0 ,n with constant
LF > 0
(b) the (algebraic) splitting Function G is totally differentiable with Lipschitz-continuous
derivatives on Ud0 ,n
52
5.3 Multirate Cosimulation
Ass. 5.3 ensures that the split problems (5.18) are index-1 and have a well-defined so-
lution. For the Gauß-Seidel iteration the smoothness of the neighborhood in Ass. 5.3 is
implied by the smoothness of the right-hand-sides, see Ass. 5.2 and Definition 5.2.
Now, let the functions X, X̃ ∈ Ud0 ,n be given and further k denote the number of
iterations on the n−th time window. This allows for the following abbreviations, [9]
Ynk := Ψky,n X, Zkn := Ψkz,n X, and Ỹnk := Ψky,n X̃, Z̃kn := Ψkz,n X̃, (5.24)
that are introduced to measure the distance of a approximation from the exact solution
after k iterations, i.e., the following abbreviation for the differences
Based on the result in [7] for particular waveforms, the following general result yields an
estimate for the dynamic iteration of two arbitrary waveforms X and X̃ on Ud,n
Lemma 5.4 (Recursion estimate). Let Ass. 5.1 and Ass. 5.3 be fulfilled for the consistent
splitting functions F, G. Then there is a constant C > 1, such that for a distance
respectively.
Proof. The proof follows [9] and generalizes [5, 7]. It is split into two parts: in the first
part the recursion estimate (5.27) for the differential unknowns ∆ky,n is shown (similarly to
the classical proof of the Picard-Lindelöf iteration). The second part proves the estimate
for the algebraic unknowns ∆kz,n using a homotopy between the two waveforms (5.26), see
Fig. 5.6.
Estimate for the differential components. Inserting the two waveforms (5.26) into (5.18a)
yields two differential equations. Subtracting one from the other and then integrating over
53
5 Multirate Methods
y(t) Ud0 ,n
Ỹnk d0
Ỹnk−1
Ynk−1 Ynk
∆0y,n (Tn )
y|[Tn,Tn+1 ]
Ud,n
Tn t1 t2 Tn+1 t
Figure 5.6: Contraction of waveforms. A differential component is shown for the k-th iteration
of the time window Tn , see the proof of Lemma 5.4.
using the consistency of F and its Lipschitz-continuity on Ud0 ,n (Ass. 5.3). To obtain the
new waveforms at iteration k the solvability of the ODE (5.18a) is required, which can
be assured for a sufficiently small time window [Tn , τ ] by standard ODE theory. The
smallness is concretized in (5.33) by the constant C. From the fact that the initial offset
is given by error propagation and thus cannot be improved by iterations, follows
Zτ
0
≤ ∆y,n (Tn )k2 + LF ∆ky,n 2 + ∆y,n
k−1
2
+ ∆kz,n 2 + ∆z,n k−1
2
dt. (5.29)
Tn
The Implicit Function Theorem (with Ass. 5.3) allows us to solve the algebraic equa-
tion (5.18b) for Z(k) = ζ̂(Y (k) , Y (k−1) , Z(k−1) ) and analogously for Z̃(k) . This yields
∆kz,n 2
= ζ̂(Yn(k) , Yn(k−1) , Z(k−1)
n ) − ζ̂(Ỹn(k) , Ỹn(k−1) , Z̃(k−1)
n ) 2
≤ Lζ̂ ∆ky,n 2 + ∆y,n k−1
2
k−1
+ ∆z,n 2
(5.30)
with the Lipschitz constant Lζ̂ > 0. Inserting this result into (5.29) yields the maximum
∆ky,n ≤ ∆0y,n (Tn ) 2
+ L0 Hn ∆ky,n + ∆y,n
k−1 k−1
+ ∆z,n
54
5.3 Multirate Cosimulation
This proves the differential part of the estimate (5.27). To verify that the new waveforms
are in the neighborhood, the distance from the exact solution is measured. This can be
done using the estimates above but for the special cases, that one of the waveforms is
identified as the exact solution, i.e., the fixed point. We set
Then the following a priori estimates are found from the estimate (5.32)
Estimate for the algebraic components. In the second part of the proof, the inequality for
the algebraic component is shown. The key is the following homotopy for θ ∈ [0, 1]:
Y (k),θ (t) := θỸnk (t) + (1 − θ)Ynk (t), and Z(k),θ (t) := θZ̃kn (t) + (1 − θ)Zkn (t).
Insertion of the homotopies into the splitting function G defines implicitly an ‘overloaded’
version of G in the only parameter θ
∂G
G(θ) := G Y (k),θ , Y (k−1),θ , Z(k),θ , Z(k−1),θ and Gu (θ) := (θ).
∂u
where u denotes an arbitrary argument of the splitting function. The waveforms fulfill the
algebraic constraint and thus it holds that G(0) = G(1) = 0. Therefore the integral
0 = G(1) − G(0)
Z1
k k−1 k k−1
= Gy(k) (θ)∆y,n + Gy(k−1) (θ)∆y,n + Gz(k) (θ) ∆z,n + Gz(k−1) (θ) ∆z,n dθ (5.34)
0
is obtained, where all partial derivatives are identified as the error abbreviations (5.25),
∂
e.g., ∂θ Y (k),θ = Ỹnk − Ynk = ∆ky,n . The Estimate (5.33) and Ass. 5.3 guarantee that all
arguments of G are in the neighborhood Ud,n . Thus the Lipschitz continuity of G on Ud,n
(Cd ≤ d0 ) with the constant LG′ yields
Gu (θ) − Gu (0) 2 ≤ LG′ θỸnk + (1 − θ)Ynk − Ynk 2 + . . .
k−1 k−1 k−1
. . . + θZ̃n + (1 − θ)Zn − Zn 2
55
5 Multirate Methods
= LG′ θ ∆ky,n 2
k−1
+ ∆y,n 2
+ ∆kz,n 2
k−1
+ ∆z,n 2
≤ 18(1 + Lζ̂ )LG′ d. (5.35)
Now, Ass. 5.3 guarantees the regularity of Gz(k) (0) and thus left-multiplication of (5.34)
yields Z1
0 = G−1 z(k)
(0) Gz(k) (0) + Gz(k) (θ) − Gz(k) (0) ∆kz,n
k−1
0 + Gz(k−1) (0) + Gz(k−1) (θ) − Gz(k−1) (0) ∆z,n
(5.36)
+ Gy(k) (θ)∆ky,n
k−1
+ Gy(k−1) (θ)∆y,n dθ.
Furthermore Ass. 5.3 (smoothness) guarantees that the Jacobians G−1 z(k)
, Gz(k−1) , Gy(k) and
Gy(k−1) are uniformly bounded on Ud0 ,n . Let the corresponding constant be denoted by
cg . Then solving equation (5.36) for G−1
z(k)
(0)Gz(k) (0)∆kz,n = ∆kz,n and application of the
maximum norm in conjunction with estimate (5.35) yields
c̃ k−1 c̃ k
k
δz,n ≤ ||G−1
z(k)
Gz(k−1) || 2,∞ + d δz,n + d δz,n + c2g δy,n
k k−1
+ δy,n ,
2 2
with the constant c̃ := 36(1 + Lζ̂ )LG′ cg and
G−1 G (k−1)
z(k) z 2,∞
= max G−1 G (k−1)
z(k) z 2
Ynk (t), Ynk (t), Zkn (t), Zkn (t) .
t∈[Tn ,Tn+1 ]
1
Finally the estimate (5.32) with Hn < Hmax and d < min d0 /C, 2c̃ give
k
δz,n ≤ 3(1 + c̃d)c2g ∆y,n
k−1 k−1
(tn ) 2 + δy,n
c̃ k−1
+ (1 + c̃d) 2c2g L0 H + ||G−1 z(k)
G z(k−1) || 2,∞ + d δz,n
2
k−1 k−1 k−1
≤ C ∆y,n (Tn ) 2 + δy,n + CHn + αn δz,n . (5.37)
and this proves the algebraic part of the estimate (5.27). Using d < d0 the global constant
2 2 c̃
C > max 2L0 , 7, 9Lζ̂ , 3(1 + c̃d0 )cg , (1 + c̃d0 )cg L0 ,
2
is large enough to deduce all estimates, i.e., (5.32), (5.33) and (5.37).
Lemma 5.4 and especially the proof above distinguish between the different (Lipschitz)
constants and their respective origin. On the other hand the estimates in [5] are rougher.
Naturally the recursion estimate (5.27) can be brought to the form of Lemma 3.1 as it is
given in [5].
Remark 5.4 (Rougher recursion estimate). In the infinity norm the initial offset is
k−1 k−1
bounded by the maximum error on time window, i.e., |∆y,n (Tn )| ≤ δy,n . Thus (5.27)
implies the following (rougher) estimate
k 0
δy,n δy,n 1 k−1
k
≤ K k−1 + ∆y,n (Tn ) 2 , (5.38)
δz,n δz,n 0
56
5.3 Multirate Cosimulation
with a (possibly) larger constant C. The rougher estimate matches the structure in [5, 7].
The consequences of the iterative application of estimates of the structure (5.38) are
discussed in [5]. The following result is easily transferred to the general index-1 setting, it
holds that
Proposition 5.5 (Iteration estimate). Let the same assumptions as for Lemma 5.4 be
fulfilled with a constant C > 1 > αn . Then there is a new constant C0 > C such that for
all k ≥ 1 and Hn ≤ Hmax it holds:
k !
max(0,k−2) k−2 0
δy,n C(4C + 1)Hn µn 4CHn µn δy,n 1+C0Hn 0
k ≤ 0 + δy,n (Tn )
δz,n 4Cµnk−1 µkn + (µn − αn )k δz,n C0
(5.39)
where
2CHn
µn = µ(αn , Hn ) := αn + αn
√ . (5.40)
2C
+ Hn
Proof. The total error after k iterations corresponds to the iterative application of estimate
(5.38). This is the multiplication by the k-th power of the matrix K, see (5.28). Thus the
claim is deduced using the eigenvalues of K, i.e.,
1 p
2 2
λ1,2 (K) = αn + 2 CHn ± αn + 4 C Hn , (5.41)
2
the details of the proof can be found in [5].
Without discussing the contraction of the recursion above (i.e., spectral radius ρ(K) < 1)
and shifting the analysis of the contraction factor αn to Theorem 5.7, the following corollary
is immanent
Corollary 5.6 (Simple coupling and convergence Rate). The eigenvalues (5.41) determine
the rate of contraction of the recursion above for the limit Hn → 0 and α < 1.
(ii) Given a splitting (‘simple coupling’), where no algebraic constraint depends on old
algebraic iterates, i.e.,
Gz(k−1) = 0,
√
then αn = 0 and the convergence rate is O( Hn ).
(iii) Given a splitting, where no algebraic constraint depends on old iterates, i.e.,
57
5 Multirate Methods
Proof. The first claim (i) is shown by Taylor expansion of (5.41). It holds
p
αn2 + 4 C 2Hn = αn (1 + 2 C 2Hn /αn2 ) + O(Hn2 ).
and this, together with the assumptions, concludes the proof. The other claims (i) and
(ii) are shown in the same way as Lemma 5.4. One exploits that some Jacobians vanish,
namely Gz(k−1) and Gy(k−1) .
The different convergence rates imply that it is beneficial to design the coupling interface
in such a way that the contraction factor αn vanishes. This is for example obtained if
algebraic couplings are avoided. This guarantees immediately a higher order convergence
(cf., ‘simple coupling’ in [7]). Whether a differential coupling can be achieved or not,
depends obviously on the particular problem, the splitting scheme and especially on the
computational sequence of the subsystems, [5]. This will be discussed in more detail in
Section 5.3.5 and Section 5.3.6 for examples in field/circuit and semiconductor/circuit
applications, respectively.
Remark 5.5 (Convergence rate of the fractional step method). Cor. 5.6 includes two im-
portant special cases: the fractional step method, see Example 5.2, exhibits a convergence
rate in the order of O(Hn ). This special case is discussed in [136]. The same rate applies
to the ODE-ODE Gauß-Seidel splitting, but the theory in this treatise is obviously not
tailored for a deeper understanding of the ODE case, see [27].
Now, having established a recursion estimate for a finite number of iterations, it will be
used in the next section to obtain convergence and stability of the scheme.
Definition 5.4 (Global splitting error). The global error for k iterations on the n-th time
58
5.3 Multirate Cosimulation
y(t)
dy,n
ey,n ǫy,n
dy,1
dy,0
T0 T1 T2 ... Tn+1 T
Figure 5.7: Lady Windermere’s Fan. Error propagation for the differential component. The
solid line depicts the exact solution, the dotted lines are approximations, [69].
where ỹ(t) and z̃(t) denote the numerical approximations of the exact solution given by
y(t) and z(t), see (5.20).
The global error consists of the local splitting error of the current window and of the
splitting errors of previous windows, that are propagated by the initial values of each
window.
describes the difference between the exact solution x|[Tn ,Tn+1 ] and the approximation after
k iterations that is obtained when starting from the exact data x|[Tn−1 ,Tn ] on the n-th
window.
The definition above exploits the fact that the exact solution x is a fixed point of Ψn , i.e.,
x = Ψkn x and thus the definition above gives indeed the local error.
describes the difference between the solutions that are obtained by k iterations on window
n when starting from the exact solution x|[Tn−1 ,Tn ] and the approximation x̃|[Tn−1 ,Tn ] .
59
5 Multirate Methods
The splitting error definitions imply the additive error decomposition in the style of Lady
Windermere’s Fan, Fig. 5.7,
ǫy,n dy,n ey,n
= + (5.44)
ǫz,n dz,n ez,n
where the superscript for the iteration number of the error on window n is disregarded if
a fixed number kn of iteration is applied, e.g., ǫy,n := ǫky,n
n
.
Theorem 5.7 (Contraction). Let an index-1 DAE (5.11) with Ass. 5.1, the constant
extrapolation operator and consistent splitting functions (Ass. 5.3) be given. Then for d
and H < H0 small enough, the hypothesis
G−1
z(k)
Gz(k−1) 2,∞
<1 (5.45)
implies that the local error mapping decreases strictly for all k > 1, i.e.,
k−1
dx,n 2,∞
> dkx,n 2,∞
, (5.46)
Proof. The proof assumes a constant extrapolation with accuracy O(Hn ), see (5.16). The
generalization to higher order schemes is straightforward. Contraction is shown by induc-
tion on k.
Induction Basis. For the differential equation at iteration k = 0 follows
Zτ
0
y(τ ) − y (τ ) 2 = f(y, z) dt 2
≤ cf H n , τ ∈ [Tn , Tn+1 ], (5.47)
Tn
with cf := ||f(y|[Tn ,Tn+1 ] , z|[Tn ,Tn+1 ] )||2,∞ . The index-1 Ass. 5.1 together with the implicit
function theorem guarantees the solvability for the algebraic component, i.e., z = φ(y)
and thus
z(τ ) − z0 (τ ) 2
= φ(y(τ )) − φ(y(0)) 2 ≤ Lφ cf Hn (5.48)
where Lφ denotes the corresponding Lipschitz constant. Thus for Hn sufficiently small,
i.e., Hn < H0 := cf (L1φ +1) , the extrapolated waveform is in the neighborhood, i.e.,
Induction step. The error decreases if the matrix K in the recursion estimate (5.27) has a
spectral radius ρ(K) < 1, see Prop. 5.5. From the eigenvalues (5.41) it follows that αn < 1
60
5.3 Multirate Cosimulation
Proposition 5.8 (Local error estimate). Let Ass. 5.1 and Ass. 5.3 be fulfilled for the
consistent splitting functions F, G. Then for hypothesis (5.45) and H sufficiently small
(H < H0 ), there is a constant Cd⋆ independent of window size Hn , contraction factor αn
and iteration number kn , such that the local error is bounded in terms of the step size
dy,n 2
+ Hn dz,n 2
≤ Cd⋆ Hn δn0 , (5.50)
Proof. Theorem 5.7 states Φn x|[Tn−1 ,Tn ] ∈ Ud,n given H small enough, see (5.49). Thus
the recursion estimate (5.39) in Prop. 5.5 is applicable (αn < 1) for the two particular
waveforms X := x|[Tn ,Tn+1 ] (exact solution) and X̃ := Φn x|(Tn−1 ,Tn ] (extrapolation of exact
0
data). There is no initial offset: δy,n (Tn ) = 0. Finally summation of both the differential
and algebraic estimates conclude the proof using a sufficiently large constant Cd⋆ .
For a sufficiently small window size Hn and assuming hypothesis (5.45), it follows µn < 1.
Thus the dynamic iteration on a single window convergences to a fixed point as k → ∞.
After stopping this fixed point iteration after a finite number of iterations kn a local
splitting error remains. The propagation of this error to subsequent windows is analyzed
in the following section.
Error Propagation
Due to the windowing technique errors on previous windows accumulate and propagate to
the current window, Fig. 5.7. Stability and convergence require that this error is control-
lable. For this the differential and algebraic propagation errors ey,n and ez,n are analyzed.
Following [9] the Prop. 5.5 is utilized to obtain (cf. [5]):
Proposition 5.9 (Propagation error). Let an index-1 DAE (5.11) with Ass. 5.1, the
constant extrapolation operator and consistent splitting functions (Ass. 5.3) be given. If
µn < 1, then there is a constant Ce⋆ > 0, such that the propagation error is bounded
" # " #
ey,n 2 1 + Ce⋆ Hn Ce⋆ Hn ǫy,n−1 2
≤ · (5.52)
ez,n 2 Ce⋆ αn⋆ ǫz,n−1 2
61
5 Multirate Methods
Proof. The application of Prop. 5.5 with the waveforms X := Φn x|(Tn−1 ,Tn ] (extrapolation
of exact data) and X̃ := Φn x̃|(Tn−1 ,Tn ] (extrapolation of erroneous data) yields an initial
offset at Tn , which is bounded by the total error on the previous time window, i.e.,
∆0y,n (Tn ) 2
≤ y|(Tn−1 ,Tn ] − ỹ|(Tn−1 ,Tn ] 2 .
The extrapolation operator is a uniformly Lipschitz continuous mapping, see (5.16). Let
LΦ denote the corresponding constant. It follows
" 0 # " # " #
δy,n y|[Tn−1 ,Tn ] − ỹ|[Tn−1 ,Tn ] 2 ey,n−1 2
0
≤ LΦ = LΦ ,
δz,n z|[Tn−1 ,Tn ] − z̃|[Tn−1 ,Tn ] 2 ez,n−1 2
such that the proof is completed by the application of Prop. 5.5 for the particular choice
of waveforms above.
Finally, combining all the previous results, the following theorem guarantees stability and
from that the global convergence result is deduced by iterative application, [5].
Theorem 5.10 (Stability). Let an index-1 DAE (5.11) with Ass. 5.1, the constant extra-
polation operator and consistent splitting functions (Ass. 5.3) be given. If the contractivity
constant is bounded
km
αm ≤ ᾱ < 1 and LΦ αm ≤ ᾱ for 0 ≤ m ≤ n,
ǫy,m 2
+ ǫz,m 2
≤ d for 0 ≤ m < n,
then there is a constant C ⋆ > 0 (independent of the window number n and of the window
sizes Hm ) such that the global error on the n-th time window satisfies
0
ǫy,n 2
+ ǫz,n 2
≤ C ⋆ max δm ≤ d (5.54)
0≤m<n
2CHm p
µ m = αm + αm
√ < α m + C α,m Hm (5.55)
2C
+ Hm
km
where Cα,m is a sufficiently large constant. Then the assumption LΦ αm ≤ ᾱ < 1 yields
p km p
⋆
αm = LΦ µkmm + (µm − αm )km < LΦ αm + Cα,m Hm ) + (Cα,m Hm )km ) < 1,
62
5.3 Multirate Cosimulation
! !
ǫy,n 2
1 + Ce⋆ H Ce⋆ H ǫy,n−1 2
Cd⋆ Hδn0
≤ · + .
ǫz,n 2
Ce⋆ α⋆ ǫz,n−1 2
Cd⋆ δn0
This proves the first half of the inequality (5.54), while the second half is clear by the
0
definition of the extrapolation operator, i.e., for the constant extrapolation δm = O(Hm ).
Thus the extrapolation errors can be made arbitrarily small by the window size Hm .
Finally Theorem 5.10 above is applied iteratively, [9, 5]. This guarantees that the ap-
proximating waveform remains in the neighborhood of the exact solution, which depends
only on the sizes of the time windows Hn . This is the desired global convergence and
stability result.
Corollary 5.11 (Convergence and stability). Let the same assumptions as for Theo-
rem 5.10 be fulfilled. Then there is a constant C ⋆ , such that
0
ỹ|[0,te ] − y|[0,te ] 2
+ z̃|[0,te ] − z|[0,te ] 2
≤ C ⋆ · max δm ,
0≤m<N
0
where δm is the extrapolation error on the m-th window.
Let us conclude that the only important constraint for stability and convergence is the
hypothesis (5.45), i.e., the algebraic-to-algebraic coupling, [5, 7]
G−1
z(k)
Gz(k−1) 2,∞
< 1,
while the dependence on the window size is a natural condition that cannot be circum-
vented. In the following Sections 5.3.5 and 5.3.6 cosimulations of applications from elec-
trical engineering are discussed; the focus is especially on the hypothesis above and how
this coupling can be avoided in practice.
63
5 Multirate Methods
there is no error control and convergence analysis. Thus this treatise features the dynamic
iteration approach as introduced in Section 5.3.2 and analyzed in Sections 5.3.3-5.3.4.
The DAE-index of the subsystems and the coupling interface are the crucial points
for the convergence of a dynamic iteration scheme. The hypothesis (5.45) identified the
exchange of algebraic variables as the mathematical reason for divergence of the iteration
scheme. Thus an interface is required that circumvents those problems. Let us recapitulate
the coupled field/circuit system in order to derive an adequate interface. It reads for the
flux/charge oriented MNA and excited by stranded conductors, see Section 3.2
d
AC q + AR gR (A⊤R u, t) + AL iL + AV iV + AI is (t) + AM iM = 0,
dt
q − qC (A⊤C u, t) = 0,
d (5.56a)
φ − A⊤L u = 0,
dt
φ − φL (iL , t) = 0,
A⊤V u − vs (t) = 0
64
5.3 Multirate Cosimulation
replacemen
iM LM
Figure 5.8: Field/circuit coupling interfaces. In interface (a) the coupling to the circuit is
given by the current through the device and in model (b) by an extracted inductance, [9].
(a) Source coupling. For the given voltage vM the system (5.56) is used to compute
(k)
a new waveform for the current iM . Then it is inserted as a time-dependent current
source into the network equations (5.56a), see Fig. 5.8a. This interface corresponds
to the model in [14]. The circuit input on the n-th time window is given by
(k)
is (t) := iM (t), for t ∈ [Tn , Tn+1 ]
and AI = AM , where it is assumed for simplicity of notation that the current through
the MQS device is the only current source in the circuit.
(b) Parameter coupling. For the given voltage vM the system (5.56) is used to compute
a(k) . Then a lumped parameter model,
the saturation level, i.e., the vector potential ⌢
e.g., a time-dependent inductance matrix and a constant resistance are extracted,
see, Fig. 5.8b. For the n-th time window, i.e., t ∈ [Tn , Tn+1 ], it holds
gM := R−1
M AM u with RM := X⊤M M−1
σ XM ,
(k) (k) (k)
φM := LM (t)iL + φeddy (t) with LM (t) := X⊤M k−1
ν
⌢(k)
a (t) XM (5.57)
(k) (k) (k)
and φeddy (t) := LM (t)iM (t) − X⊤M ⌢
a(k) (t) ,
(c) MOR coupling. This is a generalization of (b), where LM may describe an arbitrary
(linear) system obtained by Model Order Reduction (MOR): for the given voltage vM
the system (5.56) is used to compute the MVP ⌢ a(k) (t). Then a standard model order
reduction can be applied to the system (5.56), e.g. proper orthogonal decomposition,
[106]. It may utilize the waveform of the MVP on t ∈ [Tn , Tn+1 ], for example such
that the curl-curl term loses its nonlinear character, i.e., k−1
ν
⌢(k)
a (t) depends only
on time.
In practice the subproblems (field and circuit) are solved numerically by a single rate
time-integrator, see Section 5.1. Hence the solution is not a continuous waveform, but a
series of discrete solutions at time steps ti . A (continuous) waveform can be reconstructed
from them for example using dense-output or spline interpolation, [69]. Obviously, the
65
5 Multirate Methods
0.8
0.4 H
0.2
0 source coupling
exact solution
-0.2
0 0.2 0.4 0.6 0.8 1
time [s]
Figure 5.9: Qualitative behavior (‘festoon-like’) of the waveforms when using the source cou-
pling (a) in Definition 5.7. Window size H = 10−2 .
interpolation must be of sufficient quality (order), so that after exchanging the waveforms
the other time-integrator can benefit from it.
In a monolithic simulation, the interfaces (a) and (b) in Definition 5.7 are equivalent:
one can interpret (5.57) as the Schur complement of the system (5.56). In fact, this is
another way of looking at the bypassing approach in Section 5.2. On the contrary approach
(c) will introduce a modeling error depending on the reduction technique, [106]. For the
weak coupling by a dynamic iteration scheme, the models behave differently: the source
approach (a) is only a black-box coupling. Additional evaluations of the source model at
time points ti ∈ (Tn , Tn+1 ) do not reflect the physical behavior, see Fig. 5.9. On the other
hand the lumped parameter approach (b) still models the inductive effect (Faraday’s Law)
correctly. The drawback is the additional computational cost: the inductance matrices LM
must be computed for each time step. This drawback applies even more to general MOR
approaches. Nonetheless these approaches pay off in practice due to the better decoupling.
Remark 5.6 (Interface and DAE-index of the MQS device, [117]). It has been shown
in Example 5.1 that the design of the coupling interface may change the DAE-index.
For the field/circuit case the models in Definition 5.7 guarantee that the MQS device-
subsystem remains an index-1 problem. On the other hand the current-driven case, e.g,
where equation (5.56c) is fed by a given current and the voltage vM is computed and
supplied to the network equations (3.1e) corresponds to an index-2 system (Theorem 4.5)
and thus the dynamic iteration theory of Section 5.3.2 is not applicable, [117].
For the circuit system the standard loop and cutset conditions from Theorem 4.2 are
assumed and thus the circuit system (5.56a) is index-1. Moreover, the coupled system is
index-1 (another application of Theorem 4.2). Then the following problem description is
derived, where the MQS device subsystem is abstractly addressed by the functions with
subscript M and the circuit subsystem by subscript C (the structure below is the same for
all interfaces in Definition 5.7)
66
5.3 Multirate Cosimulation
with regular ∂gM /∂zM and ∂gC /∂zC . The variables of the field and circuit equations are
⌢
Qσ a
⌢ u
Pσ a iM
and yC := q , zC := iL
yM := , zM := vM
φM φ
iV
LM
where y denotes differential (i.e., defined by differential equations) and z algebraic compo-
nents (i.e., only defined by algebraic constraints). Similar to flux/charge oriented MNA,
the magnetic φM := XM ⌢ a is typically not computed as an explicit unknown of the sys-
tem. Projector Pσ picks out the differential part of the magnetic vector potential, i.e., the
components defined in conductive materials, see Definition 2.8. A detailed derivation of
structure (5.58) is given in [117].
Now, having defined the subproblems and identified the coupling variables, a splitting
scheme must be chosen that defines the computational sequence of the waveforms. From
its structure a convergence guarantee can be deduced (depending on hypothesis (5.45)).
Following [118, 9] a Gauß-Seidel-type dynamic iteration scheme is applied to (5.58), so the
newest data available is always exploited. Let us start with the computation of the field
subproblem. The splitting scheme reads
" #
(k) (k−1)
f (z , z )
F y(k) , y(k−1) , z(k) , z(k−1) := M M(k) C(k) , (5.59a)
fC (yC , zC )
" #
(k) (k)
gM (yM , zM )
G y(k) , y(k−1) , z(k) , z(k−1) := (k) (k) (k) , (5.59b)
gC (yC , zC , zM )
(k−1)
with y⊤ := [y⊤M y⊤C ] and z⊤ := [z⊤M z⊤C ]. The only old iterate zC , i.e., the voltage
drop defined by the circuit, enters a differential equation via the function fM . Thus the
contraction factor α vanishes, see Cor. 5.6 (iii).
Theorem 5.12 (Convergence of field/circuit cosimulation). Let the assumptions of The-
orem 5.10 be given. Then the Gauß-Seidel-type dynamic iteration of the field (5.56) and
circuit subsystems (3.1) coupled by one of the interfaces from Definition 5.7 is uncondi-
tionally stable and convergent with a window-wise convergence rate O(Hn ), if the iteration
starts with the computation of the MQS device.
Proof. Application of Cor. 5.6 to the splitting functions (5.59).
A reordering of the computational sequence creates a mutual algebraic dependence. In
contrast to the previous scheme (5.59) the Gauß-Seidel scheme
" #
(k) (k)
f (zM , zC )
F̃ y(k) , y(k−1) , z(k) , z(k−1) := M (k) (k) , (5.60a)
fC (yC , zC )
" #
(k) (k)
g M (y , z )
G̃ y(k) , y(k−1) , z(k) , z(k−1) := M M
(k) (k) (k−1) , (5.60b)
gC (yC , zC , zM )
G−1
z(k)
Gz(k−1) 2,∞
> 0.
67
5 Multirate Methods
200 0.06
pwm
150 sine
0.05 150 maximum at 45.7764 Hz
100
0.04
inductance [H]
50
voltage [V]
Signal
100
0 0.03
−50
0.02
−100 50
0.01
−150
−200 0 0
0 0.02 0.04 0.06 0.08 0 0.002 0.004 0.006 0.008 0.01 0 200 400 600 800 1000
time [s] time [s] Frequency [Hz]
(a) PWM and sine waveforms (b) self-inductance (c) PWM frequency spectrum
Thus the smallness of the contraction factor α is not automatically fulfilled, see equation
(5.28). Consequently convergence is not guaranteed by the theory described above and
divergence might occur. Even if 0 < α < 1 holds true in (5.60) the convergence properties
of (5.58) would be better, see Cor. 5.6. Thus scheme (5.59) is employed in the following,
see Listing 2 on Page 71.
As previously mentioned, cosimulation is especially efficient if the adaptive time-
integration of the subproblems can exploit different time scales.
The drawback of single rate time-integration is that it resolves the dynamics of a system
as a whole, cf., Sections 5.1 and 5.2.1. Thus it yields a series of time steps that matches
the dynamics of the most active component, i.e., the one working at the highest frequency.
In coupled, multiphysical systems (e.g. electromagnetic problems with heating effects)
one can split the equations corresponding to their time constants on the basis of physical
reasoning. Let us consider a 2D model of an induction motor, [40]. Its rotor position
can be updated using an additional ODE and, furthermore, changes in the conductivity
due to the MQS device heating up can be modeled by another PDE. Here the separation
into subproblems allows immediately for a efficient cosimulation exploiting the inherent
time-constant of the subproblems. This approach is particularly convenient when different
physical effects are simulated by different software tools.
In contrast to this, the field and circuit subproblems describe the same physical phenom-
ena and one speaks of ‘refined modeling’, [11]. Hence the subproblems feature similar time
constants. Nevertheless, due to switches or filters in the circuit there may only a subset
of the devices active at any given time, while the others remain latent. Then cosimulation
is more efficient than the single-rate approach, see the example in Section 6.3.1 where a
low-pass filter causes dynamics at different rates.
If the circuit topology does not provide such a splitting, the time rates are not well
separated in the coupling variables. Let us consider another example, where the voltage
that is applied to the MQS device is a pulse-width-modulated (PWM) sine-wave switching
at 20 kHz, see Fig. 5.10. The saturation, i.e., the nonlinearity of the inductance matrix,
is characterized by the underlying, much slower sine wave, see Fig. 5.10b. Mathematically
speaking the saturation of nonlinear materials in the PDE model depends on the energy
supplied, see (5.10) in the section on bypassing Section 5.2. This is related to the time
68
5.3 Multirate Cosimulation
1.2 200
pwm sine
1 sine 150 spline
integrated voltage [Vs]
0.8 100
50
voltage [V]
0.6
0
0.4
−50
0.2
−100
0 −150
−0.2 −200
0 0.02 0.04 0.06 0.08 0 0.02 0.04 0.06 0.08
time [s] time [s]
(a) integral of PWM and sine wave (b) original waveform and spline
Figure 5.11: Spline approximation of the integral waveform and its derivative.
Zt
ψ(t) = vM (s) ds with t ∈ [Tn , Tn+1 ] .
Tn
Consequently the relevant time rates of the nonlinear behavior is given by the time rates
of the integral above, even if the voltage applied is a much faster switching signal. The
waveform in Fig. 5.10a shows the PWM voltage while Fig. 5.11a depicts its integral with
respect to time. The integral is a step function approximating the cosine at a frequency of
50Hz, see Fig. 5.10c. For high frequencies the approximation is very accurate (below the
accuracy of the nonlinear curve), so one could use the smooth signal instead of the step
function.
The same approach is valid if the signal (voltage) is composed of a fast and a slow part
Problems often exhibit a time rate of interest that is given by vslow , and either the amplitude
of the fast signal max |ψfast | ≪ max |ψslow | is negligible or the energy of the fast switching
voltage vfast (t) evolves at a slower time rate. The impact of the fast signal on the nonlinear
effects can be disregarded, similar to thermal coupling that takes effect only at a slow time
rate due to energy transport, [47, 7].
69
5 Multirate Methods
The known waveform vfast (t) is integrated in time on the interval In := [Tn , Tn+1 ], which
corresponds to a summation of the discrete solutions vM (ti ) (ti ∈ In ) multiplied by the
corresponding step sizes hi
n
X
ψn := vM (ti )hi . (5.61)
i=1
Now we shall define a cubic spline interpolation ψ̃(t) of ψ(t) that gives a smooth time-
integrated voltage. The interpolation knots should be chosen accordingly to the dynamics
of the low-frequency part of the waveform. For example the number of knots can eas-
ily be estimated by a Fourier analysis or might be known beforehand. In the example
Fig. 5.10b) 20 uniformly distributed knots per period were chosen, and these yield a sat-
isfactory approximation of the sinusoidal waveform. It is crucial that the spline is a good
approximation of (5.61), because otherwise the energy balance will be violated.
Finally the cubic polynomials of the spline interpolation are differentiated piecewise with
respect to time yielding a slowly varying spline approximation ṽ slow (t) to the low-frequency
part of the waveform. The adaptive time-integrator will require only a few time steps for
the smooth slowly-varying signal compared to the hundreds of steps that are necessary to
sample the original, non-smooth and fast switching PWM signal.
This smoothing requires only a small change in Listing 2: the integration, resampling and
derivation of the spline for the voltage excitation mast be done in Step 2d. This approach
can reduce the computational effort significantly, but it introduces a model error that
relates to the frequency of the pulsed signal and the accuracy of the spline interpolation.
Conclusions
In this section the stability and convergence of the field/circuit cosimulation were mathe-
matically analyzed, following [118, 9]. It was shown that different coupling interfaces are
feasible, but that the parameter coupling is superior because it reflects the underlying laws
of physics and exploits multirate behavior (located in the nonlinearity of the MQS device).
This is documented by numerical examples in Sections 6.3 and 6.3.1.
70
5.3 Multirate Cosimulation
where u(k) is the waveform of the node potentials from the last iteration on the
previous window. Set k := 1 and start dynamic iteration, i.e., go to Step 2).
4) Next window. If Tn+1 ≥ te then go to Step 5), else set new initial values
Determine the new window size Hn+1 , e.g. from the step size predictor of the time-
integrators in Step 2a) and 2b). Go to Step 1) with n := n + 1.
5) Stop.
71
5 Multirate Methods
d
AC q + AR gR (A⊤R u, t) + AL iL + AV iV + AI is (t) + AD iD = 0,
dt
q − qC (A⊤C u, t) = 0,
d (5.62a)
φ − A⊤L u = 0,
dt
φ − φL (iL , t) = 0,
A⊤V u − vs (t) = 0,
72
5.3 Multirate Cosimulation
iSD iD
circuit device circuit device
vD vD
(a) source coupling (b) parameter coupling
Figure 5.12: Semiconductor/circuit interfaces. (a) the displacement current is modeled in the
PDE device (semiconductor) and (b) given by a parallel lumped capacitance in the network
model (circuit), see [1, 9].
to the model by [7]. The waveform on the n-th time window is given by
(k)
is (t) := iD (t), for t ∈ [Tn , Tn+1 ]
and AI = AD , where it is assumed for simplicity of notation that the current through
the semiconductor is the only current source in the circuit.
(b) Parameter coupling. The capacitance CD can be extracted from the semiconduc-
tor model beforehand. From the given voltage vD the current iSD is computed, i.e.,
the one without the displacement current, see equation (5.62b) and Fig. 5.12b. For
the n-th time window, i.e., t ∈ [Tn , Tn+1 ], it holds
d d (k)
AC q := AD CD A⊤D u and is (t) := iSD (t) (5.63)
dt dt
with the incidence matrix AI = AC = AD , where it is assumed for simplicity of
notation that the semiconductor is the only capacitance and current source in the
circuit. The new characteristic equations above replace the corresponding definitions
in the network equations (5.62a). Please note that (5.63) uses the traditional MNA
for the constant capacitance CD . This does not compromise the conservation laws
because the charge is given here by a linear relation, [57].
Remark 5.7 (Parameter coupling for semiconductors). The parameter coupling for semi-
conductors in (5.63) uses additional information on the model for the capacitive effect, but
it is still a source coupling (it directly utilizes the current iSD ). The parameter approach
can be further improved by using the nonlinear lumped conductance
(k)
iSD (t)
gD (A⊤D u, t) := (k)
A⊤D u for t ∈ [Tn , Tn+1 ]
vD (t)
(k)
instead of the plain current iSD . The equation above is given for the special case of a diode,
i.e., scalar voltage and current.
In the case of a monolithic coupling, the interfaces (a)-(b) are equivalent, [9]. On the
other hand for a weak coupling by a dynamic iteration scheme, they behave differently,
cf. Definition 5.7. The additional lumped device in the interface (b), i.e., the capacitance
CD , is a simple compact model for the displacement current of the diode. It reflects the
physical law correctly, even though in lumped form. This idea is similarly used in [51],
73
5 Multirate Methods
In the following section the two coupling interfaces are analyzed in the framework of
dynamic iteration, i.e., Section 5.3.4.
It is assumed that the circuit system fulfills the standard loop and cutset conditions from
Theorem 4.2 for both coupling interfaces, i.e., neither the current sources iD and iSD intro-
duce LI-cutsets nor the capacitance CD creates a CV -loop, see Section 4. For the semicon-
ductor subsystem (5.62c) a prescribed voltage drop vD = A⊤D u also yields an subsystem
index-1 DAE , [26]. Moreover, it can be assured that the coupled system is monolithic
index-1, [120]. This allows the application of the dynamic iteration theory above and the
following analysis is eligible, [1, 9].
where the circuit subsystem is denoted by subscript C and the semiconductor subsystem
by D. The partial derivatives ∂gD /∂zD and ∂gC /∂zC are regular due to the index-1
assumptions above. The differential and algebraic components of the subsystems are given
by
u
n Φ q
yD := , zD := and yC := , zC := iL ,
p iD φ
iV
where the space discrete electric scalar potential is denoted by the vector Φ and the device
current by iD , see the definitions in (5.62c).
The circuit is given in flux/charge oriented MNA and thus the node potentials u and
consequently the voltage drop vD = A⊤D u are algebraic variables of the circuit. The device
subsystem depends only on the algebraic circuit variables zC in the algebraic equation
gD . In turn the algebraic variable zC (that contains the device current) may enter the
differential fC and/or the algebraic equations gC of the circuit subsystem (5.64) (that
depends on the circuit topology), [1, 9]. In a nutshell: the algebraic equations of both
subsystems depend on the algebraic variable of the other subsystem. As a consequence
the contraction factor does not vanish, i.e., α > 0, independently of the computational
sequence of the subsystems in the dynamic iteration scheme. For example the splitting
74
5.3 Multirate Cosimulation
with y⊤ := [y⊤D y⊤C ] and z⊤ := [z⊤D z⊤C ]. The superscript (k) denotes the iteration number. It
encodes the computational order. Due to the dependence of G on an old algebraic iterate
convergence cannot be guaranteed by structural analysis. The parameters of all devices
and the circuit topology will have a serious influence, see Section 6.3.2. Thus the following
result is obtained, [1]
Proof. From the splitting functions (5.65) it follows that the hypothesis (5.45) is not triv-
ially fulfilled, i.e., convergence and stability are only guaranteed for cases where α < 1.
Parameter coupling. The interface (b) of Definition 5.8 replaces the current iD in the
balance equation (3.1a) by the current iSD with a parallel capacitance CD in traditional
MNA notation, i.e., by introducing some node potentials as differential unknowns:
where QD is a projector onto the kernel of A⊤D and PD = I − QD its complement. The
projectors separate the differential from the algebraic components, i.e., PD picks out the
difference of the node potentials at the capacitance CD and QD adresses the other po-
tentials. Consequently the current iSD enters only in fC via zD , because it is related to a
differential equation.
When starting, as before, with the semiconductor device computation, the splitting
functions read (device first)
" #
(k) (k)
fD (y , z )
F(y(k) , y(k−1) , z(k) , z(k−1) ) := D D
(k) (k) (k) ,
fC (yC , zC , zD )
"
(k) (k) (k−1)
# (5.67)
(k) (k−1) (k) (k−1) gD (yD , zD , yC )
G(y , y ,z ,z ) := (k) (k) .
gC (yC , zC )
75
5 Multirate Methods
The capacitive path between the coupling nodes ensures that the voltage drop vD is part
of the differential variables yC and thus the only old iterate used is differential. Conse-
quently the contraction factor α vanishes
√ for the splitting functions (5.67) and Cor. 5.6 (ii)
guarantees a convergence rate of O( H).
On the other hand for the reversed computational order (circuit first)
" #
(k) (k)
fD (y , z )
F̂(y(k) , y(k−1) , z(k) , z(k−1) ) := D D
(k) (k) (k−1) ,
fC (yC , zC , zD )
"
(k) (k) (k)
# (5.68)
g (y , z , y )
Ĝ(y(k) , y(k−1) , z(k) , z(k−1) ) := D D (k)D (k) C ,
gC (yC , zC )
there is no dependence on old iterates in any algebraic equation, and thus Cor. 5.6 (iii)
promises a higher convergence rate, i.e., O(H), [9]. The alleged difference in the convergence
rates is analyzed in the next section.
Definition 5.9 (Lipschitz constants). Let Assumptions 5.1 and 5.2 be fulfilled and let
the Lipschitz-continuous functions ζ C and ζ D define zC and zD by the implicit function
theorem applied to gC and gD from (5.66), respectively. Then the following Lipschitz
constants are defined
• let L denote the maximum of the Lipschitz constants of f⋆ and ζ ⋆ w.r.t. y⋆ and z⋆
for ⋆ ∈ {D, C}
The latter Lipschitz constants LC and LD are a measure of the strength of the mutual
coupling between the semiconductor and circuit subsystems. Following [9], a problem-
specific version of the recursion estimate Lemma 5.4 is found
Lemma 5.14 (Refined recursion estimate). Let the assumptions of Lemma 5.4 be fulfilled.
Then for the Gauß-Seidel-type dynamic iteration scheme (5.67) the recursion estimate
k k−1
δy,n CD Hn 0 δy,n 1 + CHn
≤ k−1 + ∆0y,n (Tn ) 2 (5.69)
k
δz,n C 0 δz,n C
| {z }
=:KD
76
5.3 Multirate Cosimulation
holds true and similarly for the reversed order scheme (5.68)
k k−1
δy,n 0 CHn δy,n 1 + CHn
≤ + ∆0y,n (Tn ) (5.70)
k
δz,n 0 CC Hn δz,n k−1
C 2
| {z }
=:KC
where Hmax > Hn denotes the maximum time window size and C > 1 is a sufficiently large
constant.
Proof. The proof is basically the same as for Lemma 5.4, but here the Lipschitz constants
from Definition 5.9 are analyzed separately. Instead of equation (5.30) one finds for the
particular splitting scheme (5.67)
∆kz,n 2
≤ L ∆ky,n 2
k−1
+ LD ∆y,n 2
(5.72)
with Lipschitz constants L, LC and LD as defined in Definition 5.9. Now the insertion of
k
(5.72) into (5.73) and solving for the differential difference δy,n yields
k
δy,n ≤ (1 + CHn ) ∆0y,n (Tn ) 2
k−1
+ CD Hn δy,n , (5.74)
with the constant CD as defined above, C as defined in (5.76) and using the upper bound
of the time window size H < Hmax . Then (5.72) and (5.74) imply
k
δz,n ≤ C ∆0y,n (Tn ) 2
k−1
+ Cδy,n . (5.75)
This concludes the proof of (5.69); the other estimate (5.70) is shown analogously.
77
5 Multirate Methods
Proof. Stability and convergence follow from Cor. 5.6. Thus only the convergence rate
O(H) must be shown. The spectral radii of the iteration matrices KC and KD are given
by ρ(KC ) = CC Hn and ρ(KD ) = CD Hn , respectively. Thus for H ≥ max Hn follows
convergence rate O(H), where CC and CD are estimates for the leading coefficients.
Although both computational sequences have the same order, the proof shows that the
speed of the iteration scheme depends on different leading order coefficients, i.e., CC and
CD . Again, their estimates differ only by the Lipschitz constants LC and LD , reflecting the
strength of the coupling via differential and algebraic equations, see (5.71). This difference
can be observed in numerical simulations, see Section 6.3.2.
Conclusions
In this section the stability and convergence of the semiconductor/circuit cosimulation was
mathematically analyzed, following [1, 9]. It was shown that only parameter coupling guar-
antees convergence for both computational sequences. This is documented by a numerical
example using the different coupling interfaces and sequences in Section 6.3.2.
78
5.4 Domain Substructuring in MQS Devices
(a) nonlinear domain: iron core (b) linear domain: copper coils (conductors)
Figure 5.13: Domain substructuring of a transformer. Iron core exhibits eddy currents and
a nonlinear permeability. The surrounding air and the coils are modeled by static and linear
equations (no eddy currents, strands below skin-depth), cf. Fig. 3.4.
J(⌢
a) := λMσ̃ + Kν (⌢
a) (5.77)
where λ is related to the time step size. For the implicit Euler method with constant time
step size h, it holds: λ = 1/h, Section 5.1. The differential curl-curl matrix Kν (⌢ a), see
Definition 2.6 and the conductivity matrix Mσ̃ have typically a common nullspace in 3D
formulations. This is removed by a regularization technique, e.g. Grad-Div, as given in
Section 2.2.3.
In relevant technical applications, typically only few materials are modeled in such a
way that (non-)linear or even hysteretic behavior is taken into account. Furthermore
many regions are non-conductive, either because the material is nonconductive (air) or
because it is modeled non-conductively to prevent eddy currents, e.g. without resolving
laminations or windings in the discretization (e.g. as for copper coils, see Section 3.2).
This turns the corresponding degrees of freedom into algebraic variables and increases the
nullspace of the conductivity matrix Mσ̃ .
Consequently the models contain large regions described only by linear algebraic rela-
tions. The entries of the system matrix J related to edges that are strictly linear model
parts do not change during the overall simulation. This is commonly exploited in profes-
sional implementations, e.g., the assembly of the corresponding elements is bypassed. But
this does not fully exploit the structure: the linear model parts are still solved unnecessarily
in every step of the nonlinear iteration in the time-integration, [36].
The degrees of freedom of the eddy-current problem (2.18) can be separated by projec-
79
5 Multirate Methods
where K11 and K22 are the curl-curl matrices of the respective subdomains and K12 is a
interface matrix. The regularity of K22 follows from the Gauging Assumption 2.13. Thus
(5.78) fulfills the standard criterion for a differential-index-1 equation, Section 4, [101].
The first equation (5.78a) is an ordinary differential equation with a positive definite
conductivity matrix M11 and the second equation (5.78b) is an algebraic constraint. How-
ever time-discretization turns (5.78) into a nonlinear algebraic problem that has to be
solved by Newton-Raphson.
For a given ⌢ a1 we can solve the algebraic equation for the non-conductive domain and
reinsert the resulting ⌢a2 into the differential equation (5.78a). This is in terms of discrete
domain substructuring the Schur complement
Σ11 (⌢
a 1 ) := K11 (⌢
a1 ) − H11 with H11 := K12 (K22 )−1 K⊤12
that was already used in the index-analysis, see equation (4.16). This approach is known
in DAE theory as index-reduction, because we have transformed the index-1 DAE (5.78)
into an ordinary differential equation (index-0)
⌢
⌢
M11 dt ⌢
a1 + Σ11 (⌢
a 1 )⌢
a1 = j 11 (5.79)
⌢
⌢ ⌢
⌢ ⌢
⌢
with the reduced right-hand-side j 11 := j s,1 − K12 (K22 )−1 j s,2 . If the application of the
inverse K22 is feasible, the Schur complement system is iteratively solved by Listing 3, that
is basically the preconditioned conjugate gradients method with some additional linear
solves for the complement, [110].
In Listing 3, the matrix J11 := λM11 + K11 (·) denotes the upper left block of the matrix
pencil (5.77) and ⌢ a1 (t0 ) is the initial value for time-integration and P is a matrix for
preconditioning. So far, this method only differs from the Schur complement equipped by
standard CG by the fact that, here, the Schur complement is not explicitly formed, [36].
80
5.4 Domain Substructuring in MQS Devices
end
λ. This shifts the corresponding eigenvalues of the curl-curl matrix (further) to the right
on the positive axis, while the eigenvalues of the non-conductive domain remain unaltered,
Fig. 5.14a.
By construction, all eigenvalues of the Schur system (5.79) are affected by conductivities
(M11 has full rank). This diminishes the influence of the material jump from the system’s
spectrum and improves the speed of convergence of the CG method. Additionally, the
(maximum) number of iterations decreases because the system has less DoFs than before.2
Fig. 5.14b shows the improvement in the convergence rate using the 2D example of Sec-
tion 6.4. This plot does not take the additional computational costs into account: the
higher convergence speed comes at the price of a more expensive iteration since in each
step an inner solver computes the Schur system, see Listing 3
⌢
⌢
(i+1) (i)
K22 p2 = j s,2 − K⊤12 p1 .
2
In exact arithmetics CG needs fewer iterations than the number of DoF, [110].
81
5 Multirate Methods
101
full system
eig(C̃Mν C)
relative residual
100 schur
eig(Mσ̃ ) 10−1
10−2
eig(λMσ̃ + C̃Mν C) 10−3
10−4
0 4 6 8 10 12 0 100 200 30 400 500
10 10 10 10 10 iteration number
(a) eigenvalue spectrum for λ = 103 (b) CG convergence
Figure 5.14: Eigenvalues and PCG convergence. The eigenvalues in (a) are given for the FIT
discretization of a 3D transformer model: 882 of 4752 eigenvalues of the curl-curl matrix are
zero, 2457 eigenvalues of the conductivity matrix are zero and finally their sum has 673 zeros.
The plot (b) shows the convergence for the full and Schur system (5.79) for a 2D model (the
models are discussed in Section 6.4, see [36]).
• an iterative method, preferably equipped with deflation, in order to benefit from the
repeated solving with the same system matrix, e.g. [52, 35]
• MOR (e.g. using proper orthogonal decomposition, Krylov techniques, etc.) con-
structed from the first solutions of the time-stepping process;
The first three ‘solvers’ compute (approximations of) the solution that belong to the origi-
nal problem, whereas the last two alternatives solve a modified problem (in the static and
linear regions).
The second approach is discussed in Section 6.2 for the eddy-current problem of a trans-
former model in 2D and 3D formulations. It features a sparse Cholesky factorization of
the matrix K22 . In 2D this approach works very well, while in 3D the additional burden
of the forward/backward substitutions in each iterations may become dominant, [36].
As we have said before, in 3D the block K22 has a non-trivial nullspace (the gradient
fields, see (2.21)). Thus standard factorizations fail and a regularization must be applied,
Section 2.2.3. For the Schur complement method an additional constraint for the regular-
ization is crucial to preserve: the separation between conductive and non-conductive model
82
5.4 Domain Substructuring in MQS Devices
regions. This ensures that the Schur complement removes the (conductivity) material jump
completely and only then can a significant gain in PCG’s convergence be expected, [36].
5.4.4 Conclusions
In this section we have proposed a variant of the Schur complement method for the eddy-
current problem. This adapted version has been shown to exploit the material structure,
i.e., the static and linear part. It reduced the differential-algebraic problem to an ordinary
differential equation and the resulting system matrix has an advantageous eigenvalue spec-
trum. This speeds up the convergence of the preconditioned conjugate gradient algorithm.
Computational examples are given in Section 6.2.
83
6 Numerical Examples
In this chapter the multirate methods and analysis given in Chapter 5 are numerically
verified by examples from electrical engineering.
The software is written in Octave within the framework of the demonstrator platform
of the CoMSON project (Coupled Multiscale Simulation and Optimization in Nanoelec-
tronics). For the circuit simulation part the OCS package (Octave’s Circuit Simulator) is
used. It is coupled to the new package FIDES (Field Device Simulator) for magnetoqua-
sistatic field device simulations, [113]. The 3D examples are obtained from handmade FIT
discretizations, while 2D examples are designed in FEMM and discretized by Triangle [94,
121]. Visualizations are obtained by Paraview, [98]. The workflow is shown in Fig. 6.1.
The MQS devices may consist of several conductor models (stranded or solid), connected
to the electric circuit as one multiport device, Section 3.2. The connection to OCS is
established by calling a corresponding device file in the circuit netlist, that is an IFF-file,
[55]. For the strong field/circuit coupling (monolithic coupling) the device file Mfidesmono
defines the (full or reduced) element stamp, such that the field equations are solved along
with the circuit equations by the same time stepping scheme, Section 5.2. Consequently
the device file adds additional unknowns (external and internal variables) to the circuit
problem: each conductor inside the model is excited by a voltage drop and hence it is
represented in the circuit by two pins (2 external variables), Section 3.2.
The internal variables are the magnetic vector potential ⌢
a and the currents through each
conductor model iM . They require several input parameters: the filename of the model,
followed by the specification of external variables. These specifications must meet the
topology of the field model, e.g., the number of unknown currents must match the number
of columns of the coupling matrix XM , [114].
FEMM design
Paraview visualization
Figure 6.1: Flow chart of software packages. The focus in this treatise is on the efficient
coupling of simulator packages, i.e., the dashed box. The analysis and methods from Section 5
belong there.
84
6.1 DAE-Index in Applications
RM RM
vM (t) iM LM iM (t) vM LM
Figure 6.2: DAE-index example circuits. (a) voltage-driven and (b) current-driven devices
have a different DAE-index, [8].
The procedure for the cosimulation (dynamic iteration) follows a different approach: the
resistances are defined separately and the fluxes are given by a special nonlinear inductance
device Mfidesinduct. This device file receives the extracted inductances from the field
subproblem using an outer iteration, see interface (5.57) and interpolates if necessary, [114].
The following examples are discussed in their corresponding sections
5.1 DAE-index: the increased error when solving a DAE-index-2 system monolithically
for the MQS device is analyzed using a linear axisymmetric inductance example,
5.2 bypassing: the multirate bypassing for reduced field stamps is applied to a pulsed
circuit coupled to a nonlinear 2D model of a transformer,
5.3 cosimulation: examples of the Gauß-Seidel-like dynamic iteration of MQS and
semiconductor devices with circuits; the results for the semiconductor example are
taken from [9],
5.4 substructuring: the domain substructuring method is applied to a transformer in
2D and 3D formulations.
Kν ⌢
a = XM iM
d ⊤⌢
− X a + RM iM = vM
dt M
either for given voltages v(t) or currents i(t). The magnetostatic system only exhibits
inductive effects (no eddy currents) and the coil is modeled by a stranded conductor, see
Section 3.2.1. The corresponding axisymmetric PDE model was discretized by FEMM,
85
6 Numerical Examples
iron core
air
Figure 6.3: Inductor example. Axisymmetric inductor model from FEMM, discretized by
Triangle, [94]. The coupling is established via stranded conductor models, see [9].
[94]. The device model, Fig. 6.3, is taken from the FEMM examples section in the online
tutorial (‘Inductance Calculation Example’, file: induct1a.fem1 ). In this simple case
the tractability concept matches the Kronecker and differentiation indices and thus it
corresponds to the special cases in [117, 131].
Solely the circuit topology determines the DAE-index of the coupled problem, see The-
orem 4.2 and Theorem 4.5. The index-1 and index-2 cases are obtained for the different
choices (see Fig. 6.2)
(a) voltage source connected to the device, i.e., incidence matrices AV = [1] and AM =
[−1], states an index-1 problem,
(b) current source connected to the device, i.e., incidence matrices AI = [1] and AM =
[−1], states an index-2 problem (LIM-cutset).
and thus the following lumped quantities are easily extracted from the PDE model
For the time-discretization the implicit Euler scheme was used with fixed step sizes
h = 10−11 s, 10−10 s, . . ., 10−6 s, Section 5.1. The application of higher-order methods is
1
see http://www.femm.info/wiki/InductanceExample
86
6.2 Multirate Bypassing
−6
10 −6
10
step size 1e−11
step size 1e−10
step size 1e−09
−8
10 step size 1e−08 −8
10
step size 1e−07
step size 1e−06
relative error
−10
relative error
10 −10
10
−12 −12
10 10
step size 1e−11
step size 1e−10
−14
10 −14 step size 1e−09
10
step size 1e−08
step size 1e−07
step size 1e−06
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
time [t] x 10
−5 time [t] x 10
−5
(a) relative error in current (index-1) (b) relative error in voltage (index-2)
Figure 6.4: Index-1 vs. index-2 errors. Errors for (a) the index-1 setting (driven by voltage
source) and (b) the index-2 setting (driven by current source), [9].
straightforward but it will not give further insight (e.g. using FIDES and a RADAU
scheme supplied by odepkg, see Fig. 6.1).
The numerical solutions of the index-1 and index-2 case are compared to the analytical
reference solution as given above (either the current (6.1) or voltage (6.2)). Fig. 6.4 shows
the numerical errors due to time-integration for both cases. The relative errors of both
problems behave very differently: in the index-1 setting, Fig. 6.4a, the relative error of
the current decreases with the step size. Small oscillations (near the machine precision)
occur at the smallest step size h = 10−11 s. On the other hand, in the index-2 setting
Fig. 6.4b, the error oscillates at a high amplitude for step sizes below h = 10−8 s. This
is a typical index-2 phenomenon: the error increases while the step size decreases. These
numerical results underline the difference that is mathematically described by Theorem 4.2
and Theorem 4.5.
Although index-2 problems are rather ill-conditioned, Fig. 6.4b shows clearly that the
index-2 errors are not propagated in time. Remark 4.4 explains this behavior: the index-
2 components enter the system only linearly and cannot affect subsequent time-steps,
[6, 128]. Nonetheless one must not use the index-2 variables (the voltages) for step size
control, because the (numerical) oscillations might be detected by the control and hence
the predictor would suggest unreasonably small step sizes.
87
6 Numerical Examples
24 cm PDE model
6
10
8 cm 3 cm
LM 5
10
3
10
22 cm
2
10
iL1 iL2 Rload
VSRC 1
10
0
9 cm 10
8 cm 0 0.5 1 1.5 2
magnetic flux B [T]
LU decompositions
-4 Newton
10 tolerance
current [A]
400 3
-6
10
10
-8
200 10 10
2
-10
10
0 1 standard Newton
-12 10
10 simplified Newton
-14 0
bypassed Newton
-200 10 10
0 0.002 0.004 0.006 0.008 0.01 0 0.002 0.004 0.006 0.008 0.01 0 0.002 0.004 0.006 0.008 0.01
time [s] time [s] time [s]
(d) currents through coils (e) relative errors (f) number of factorizations
Figure 6.5: Bypassing example. Device and circuit problem description, reference solution,
errors and decompositions, cf. [116].
The model has been drawn and discretized by FEMM and simulated by the software
packages OCS, FIDES and odepkg, see Fig. 6.1. The simulation is carried out only for
the startup phase, i.e., until the saturation phase is reached, Fig. 6.5d. Afterwards the
computation can be continued by a linear model without difficulty.
Fig. 6.5b shows the field/circuit coupling, where a pulse width modulated (PWM) voltage
source is connected to the primary side of a transformer. The PWM voltage is switching
at a frequency of 20kHz. The secondary side is connected to a load resistance Rload = 10Ω.
For the simulation the different strategies of Section 5.2 were implemented in
Mfidesschur for OCS, [114]: this is (a) the standard Newton-Raphson scheme that serves
as a reference, (b) a simplified Newton iteration with Schur complements and (c) a sim-
plified Newton with bypassing of right-hand-side evaluations.
The circuit device element Mfidesschur uses the bypassing heuristics defined in List-
ing 1. This algorithm allows the element to decide independently whether a new factoriza-
tion of (5.8) or a new right hand-side evaluation (5.9) is necessary or not. The advantage
is that this implementation only requires changes to the MQS device element. No other
parts of the host circuit simulator need adjustment. The other (basic) elements remain
responsible for their respective contributions and the outer (standard) Newton-Raphson
scheme is still available for them.
On the other hand, if the step size control of the host simulator is accessible, it should
be configured so as to be as conservative as possible, i.e., the step size h should be kept
constant as long as possible. In fact, step size changes require a recomputation of the
88
6.2 Multirate Bypassing
Table 6.1: Computational costs for the different Newton strategies. The ‘full Newton’ solves
the full system of equations without bypassing. The ‘simplified Newton’ solves the reduced
system and bypasses some Jacobian updates. The ‘bypassed Newton’ solves the reduced system
and bypasses Jacobian and right-hand-side updates.
Schur complement Lh but alternatively one may continue with a ‘wrong’ Jacobian, cf.
equation (5.7). For the following example from [116] this is not implemented: only back-
ward Euler with fixed step sizes is used, because the examples are determined by pulsed
inputs and step size prediction would yield high amounts of rejected steps. In principle
adaptive higher-order methods are available by odepkg.
Fig. 6.5f shows the different bypassing strategies in comparison with the reference ap-
proach. The reference method is a non-optimized Newton algorithm; it evaluates the
material curve (unnecessarily) in every iteration, see Table 6.1. On the other hand both
bypassing strategies (simplified and bypassed Newton) detect the linearity in the material
curve (for t ≤ 0.003s) and skip the superfluous evaluations and matrix factorizations. De-
pending on the error, the nonlinear effect is not important here and the reduced (lumped)
models are sufficient, see Fig. 6.5e. Those models preserve correctly the characteristics on
the fast scale, e.g., they resemble the switching in the current due to PWM, Fig. 6.5d,
while other multirate techniques would probably fail here, especially when waveforms are
interpolated and not models (‘the multirate behavior is in the nonlinearity’).
In the highly nonlinear saturation phase, 0.003s < t ≤ 0.007s, all approaches require
approximately the same (high) number of Jacobian updates. Without those updates the
Newton iteration might diverge or even converge to a wrong solution if the bypassing
of right-hand-side evaluations is not controlled. Each bypassing of the right-hand-side
assumes linearity and as a consequence the Newton iteration requires less Jacobian updates
but the error increases, see Fig. 6.5. Finally after the saturation level is reached, t > 0.007,
the field problem behaves again linearly and the updates of the simplified and bypassing
Newton are clearly reduced, [116].
6.2.1 Conclusion
The bypassing approach was shown to exploit reliably multirate potential in the nonlinear
field part. Of course the efficiency depends on the particular choice of error norms, toler-
ances and the device’s characteristics. Nonetheless, the numerical experiments indicated
that the heuristic is rather insensitive to changes in those parameters and that the compu-
tational costs can be significantly reduced even when using conservative parameters (i.e.,
small tolerances).
Obviously, the changes in the saturation cause the high computational costs. For a
transformer example this typically occurs only during the start-up phase, but in an in-
duction machine, where the saturation follows the rotation, one is forced to recompute
the Schur complement in every turn. Nonetheless the rotation is still determined by the
89
6 Numerical Examples
v(t) Rload
0
6
PDE model
(a) 2D transformer model given (b) lumped rectifier circuit, with the MQS device, four diodes,
by FEMM, [94] two winding resistances and one load resistance
200 80 4
primary
150 3.5
60 secondary
100 3
inductance [H]
voltages [V]
50 2.5
current [A]
40
0 2
-50 20 1.5
-100 1
0
source primary
-150 0.5
rectified secondary
-200 -20 0
0 0.005 0.01 0.015 0.02 0 0.005 0.01 0.015 0.02 0 0.005 0.01 0.015 0.02
time [s] time [s] time [s]
(c) voltage drop at voltage (d) currents through the trans- (e) self-inductances of trans-
source (dashed), and load resis- former’s primary and secondary former’s primary and secondary
tance (solid) coil coil
energy (5.10) and not by the fast frequency of the pulsed inputs and thus solving on the
slow scale should be sufficient for most applications. Furthermore an optimized version of
the bypassing could save previously computed inductances, e.g. in dependence of the rotor
angle, and reuse them in the subsequent simulation.
6.3 Cosimulation
In this section cosimulations of exemplary field/circuit and semiconductor/circuit problems
are performed using Gauß-Seidel-like dynamic iteration schemes, see Sections 5.3.5 and
5.3.6. For both settings the different coupling interfaces
are implemented and numerically compared. Furthermore it will shown that the conver-
gence speed predicted in Section 5.3 is verified by the simulations.
The focus in the field/circuit case is primarily on the multirate potential that can be
exploited in cosimulation, while the main topic in semiconductor/coupling is the compu-
90
6.3 Cosimulation
101 101
100 100
time-integrated error [As]
10-2 10-2
10-3 10-3
10-5 5.10-5 10-4 10-4 10-3 10-2
window size H [s] window size H [s]
Figure 6.7: Convergence of field/circuit cosimulation. Figures show the splitting error in the
primary current on the full time interval in dependence of the window size H.
tational sequence and the convergence properties. Due to similar time scales for semi-
conductors and integrated circuits there is typically no multirate potential in this kind of
simulations (with the obvious exception of devices in latent branches, [124]).
91
6 Numerical Examples
102 102
101 101
O(H)
time-integrated error [As]
10-1 10-1
10-2 10-2
10-3 10-3
10-4 10-4
1 5 10 15 1 5 10 15
number of window iterations k number of window iterations k
(a) source coupling for H = 10−4 (b) parameter coupling for H = 5 · 10−3
Figure 6.8: Contraction of field/circuit cosimulation. Splitting errors in the primary current
on the first time window [0, H] in dependence of the number of iterations.
For the numerical experiments, the time-integrator has been the implicit Euler scheme
with fixed time step h = 10−5 s, see Section 5.1. The dynamic-iteration scheme (5.58) has
been performed using various window sizes and the two interfaces (a) source coupling
and (b) parameter coupling, see Definition 5.7. The interfaces are implemented in
the functions Mfidescurrent and Mfidescurrent within FIDES, [114]. The cosimulation
results have been compared with a monolithic simulation (implicit Euler with h = 10−5 s).
Thus the comparison neglects the time-discretization error (due to the implicit Euler) and
visualizes only on the splitting error (due to the iteration scheme).
Contraction
The splitting errors of the source coupling interface are shown in Fig. 6.7a. The simula-
tions were performed with constant extrapolation and up to 5 iterations per window. No
interpolation was necessary to recover the waveforms from the discrete currents, due to the
same step size h in both subproblems. Convergence was only obtained for very small time
window sizes near the time step size, i.e., H = 2 · 10−5 s and 5 · 10−5 s. The parameter cou-
pling interface comes with additional costs: two additional linear systems must be solved
afterwards to extract the 2 × 2 inductance matrix LD (t), but the effort pays out as shown
in Fig. 6.7b. The dynamic iteration with the parameter interface converges faster and for
the much larger window sizes H = 5 · 10−5 s, . . . , 10−2 s as the source coupling (up to 5 iter-
ations per window, constant extrapolation and spline interpolation). Both results confirm
neatly the theoretic results of Section 5.3: the error decreases with the time window size if
the the window size is sufficiently small. The second convergence study, Fig. 6.8 shows the
time-integrated error on the first time window versus the number of iterations. The source
coupling converges badly on the first time window (H = 10−4 s) and it divergences on a
subsequent window (i.e., the corresponding window size Hn is not small enough). For the
parameter coupling interface, the convergence is much better even for the larger window
size H = 5 · 10−3 s. The iteration order is here approximately linear in H, which matches
the expected order in Theorem 5.12. For more than 4 iterations the splitting error is in
the order of the time-integration error.
92
6.3 Cosimulation
600
u1
PDE model u3
400
LM
L R RM,1 RM,2 200
voltage [V]
0
u1 u3
-400
0 0.005 0.01 0.015 0.02 0.025 0.03
time [s]
200
adaptive steps
voltage [V]
5 h
0
0 sweep
-200 T1 T2
-5 -400
0 0.005 0.01 0.015 0.02 0.025 0.03 0 0.005 0.01 0.015 0.02 0.025 0.03
time [s] time [s]
Figure 6.9: Multirate cosimulation example. Nonlinear field/circuit configuration (a) exhibit-
ing different time constants in the voltages u1 and u3 due to a fast switching PWM voltage
source (b); currents (c) and partitioning into time windows H and time steps h (d), [118].
For small window sizes near the step size of the time discretization (H ≈ h) the dis-
cretization error may dominate the asymptotic behavior; then additional iterations are
superfluous. In the present case of the implicit Euler method, its accuracy is first or-
der, thus depending on the quality of the initial guess, one or two iterations are typically
sufficient. This changes for higher order methods.
The rectifier example could not benefit from multirate time-integration: the time step-
ping was fixed to simplify the convergence analysis. This will be generalized in the following
multirate example.
93
6 Numerical Examples
300
reference
guess no iterations
200 iteration 1
sweep control
iteration 2
iteration 3 10.5%
100
-100 2.6%
-200
-300 0.8%
-400
0.028 0.0285 0.029 10-5 5.10-5 10-4 2.10-4
time [s] window size H [s]
(a) waveforms are improved by additional itera- (b) relative error in the current in dependence of
tions (H = 10−2 s) the window size
hL = 10−4 s would be sufficient to render the dynamics of the field model (Fig. 6.9b).
The interface couples via inductance parameters as defined in equation (5.57). All time-
integrations have been performed by the implicit Euler method, Section 5.1. The circuit
subproblem is discretized by a fixed step size of h = 10−6 s, which is reasonable for fast
switching PWM signal, while the field subsystem is solved adaptively. The low order
method was chosen to allow easy comparison, although high order adaptive multi-method
time-integration is straightforward in FIDES.
The dynamic iteration is performed with linear extrapolation, spline interpolation, fixed
time window sizes from H = 10−5 s to 2·10−4 s, either with (k ≤ 3 iterations) or without
sweep control (1 sweep), [118]. Table 6.10c shows the computational effort expressed in
linear systems solved. The first summand relates to the time-integration and the second
relates to the inductance extraction. The costs of the monolithic integration (fixed step
size href = 10−6 s) are included for comparison. The relative errors shown in Table 6.10c
and Fig. 6.10b are always given with respect to the monolithic reference solution (scaled
by the maximal current 15.3 A). Due different time-steppings for the field subsystem, the
depicted errors consists of both splitting and time-discretization errors.
The dynamic iteration (‘weak’ coupling) gives for the window size H = 10−4s the same
level of accuracy as the monolithic approach (‘strong’), but requires less than half of the
94
6.3 Cosimulation
Figure 6.11: Semiconductor/circuit cosimulation example. Model and parameters from [9].
computational effort. The additional iterations improve the accuracy but also increase the
computational costs, Fig. 6.10a. The iterations cannot improve the order of the method
Fig. 5.14b. This is mainly because the order of the time-discretization dominates the
plot. If one uses higher order methods additional iterations are important to conserve
the quality of approximation. Fig. 6.10d shows the convergence of the cosimulation when
using a higher order Runge-Kutta method (RADAU5). The convergence order in terms of
the window size increases clearly with the number of dynamic iterations (for this plot the
same example was used, but the applied input signal was smoothened to benefit from the
higher order of RADAU5).
Larger time windows (H = 2·10−4s) cause larger errors (> 12%) but at reduced com-
putational costs. On the other hand simulations with H ≤ 5 · 10−5 s require less than 70%
of the computational effort of the monolithic coupling (h = 10−5 s) while being significantly
more accurate, see Fig. 5.14b. For small window sizes the sweep control does not require
iterations and thus the curves of both methods coincide in Fig. 5.14b, [118].
Conclusions
Two examples underlined the convergence and stability of the dynamic iteration approach
for the field/circuit coupled problem. The convergence rate was discussed and it has been
shown that the scheme automatically exploits multirate behavior due to the decoupling
of subproblems. The source and parameter interface were implemented. The parameter
approach together with outer iterations allows for enlarged window sizes. For optimal
results higher order time-integration is needed and a combined window size and sweep
control will further improve the efficiency.
95
6 Numerical Examples
0 −2
10 10
circuit first circuit first
device first −3
device first
10
relative error
relative error
−1 −4
10 10
−5
10
−2 −6
10 10
0 2 4 6 8 10 0 2 4 6 8 10
number of window iterations k number of window iterations k
(a) source coupling (b) parameter coupling
Contraction
For the convergence study we analyzed the academic test example shown in Fig. 6.11a. The
diode-block consists of ND = 1500 diodes. Due to their parallel connection it is sufficient
to simulate a single PDE device and multiplying the output current by the number of
devices.
The dynamic iteration is applied with constant extrapolation and 10 iterations per win-
dow on the time interval I = [0, 10] · 10−12 s. The underlying time-discretization was
performed by the implicit Euler method. The time windows and time steps are chosen
to be the same H = h = 0.1 · 10−12 s, i.e., after each time step both subproblems are
synchronized. The detailed algorithm is similar to Listing 2; it is described in more detail
in [1].
Now the contraction analysis follows. The solutions of the cosimulation after 10 itera-
tions is compared with the monolithic reference. The reference allows the verification of
convergence of the dynamic iteration scheme. It is performed with same method and step
size, such that the comparison shows the splitting error only (time-integration errors are
neglected).
(a) Source Coupling. As predicted in Lemma 5.13 the dynamic iteration scheme does
not converge reliably for the source coupling. The convergence depends on the contraction
parameter α in (5.28). In the present numerical example the amplification of the diode’s
current by the factor 1500 directly affects the corresponding Lipschitz constants and thus
causes divergence. This is shown in Fig. 6.12a where the relative error of the network com-
ponents (with respect to the monolithic reference solution) is plotted in dependence of the
number of iterations for the time window [2.2, 2.3] · 10−12s. The dynamic iteration schemes
converges (slowly) on the previous windows until it diverges on the window depicted. The
same problem occurs independently of the computational sequence, i.e., the iteration for
both device-first and circuit-first do not converge. The different starting values for the
sequences are due to different errors on previous time windows (error propagation).
(b) Parameter Coupling. In the second approach the displacement current of the
diode is extracted and modeled by a lumped parallel capacitance, see Definition 5.8. It
96
6.3 Cosimulation
0
10
0
circuit first 10 circuit first
device first device first
−2
10 −2
relative error
10
relative error
−4
10 −4
10
−6 −6
10 10
0 2 4 6 8 10 0 2 4 6 8 10
number of window iterations k number of window iterations k
relative error
−2
10
−2
10
−3
10
−3
−4 10
10
−5 −4
10 10
0 2 4 6 8 10 0 2 4 6 8 10
number of window iterations k number of window iterations k
Figure 6.13: Contraction and Lipschitz constants. Splitting errors on the time window
[0.4, 0.5] · 10−12 s for different numbers of diodes, see [9].
is for the given example CD = 10−17 F for each diode and this accumulates to a single
capacitance of 1.5 · 10−14 F. Consequently the interface couples via the current iSD . In
contrast to the source coupling approach above the parallel capacitance aids the dynamic
iteration scheme. One obtains a robust algorithm that yields a sufficiently accurate solution
after a few iterations, see Fig. 6.12b. The convergence plot depicts the relative error (w.r.t.
the monolithic reference solution) of the network components in dependence of the number
of iterations for the same interval as above, i.e., [2.2, 2.3] · 10−12 s. Moreover, due to
significantly better convergence on the previous time windows the initial error is reduced.
Computational Sequence
The convergence study above revealed that for the parameter coupling, i.e., interface (b),
both computational sequences yield a convergent cosimulation, but at slightly different
speeds, see Fig. 6.12. According to Theorem 5.15 the speed is determined by the Lipschitz
constants and they can be influenced by the number of diodes ND . More precisely the
number of diodes affects the Lipschitz constant LF from Definition 5.9. Figure 6.13 shows
the corresponding convergence plots of the dynamic iteration scheme as above applied to
the same problems but for varying numbers ND = 1, 10, 100 and 1000. As before the
splitting error is computed with respect to the reference solution and again the time step
size h = 0.1 · 10−12 s was applied.
97
6 Numerical Examples
−1 −1
10 10
circuit first circuit first
device first device first
−2 −2
10 10
Error
Error
O(H)
−3 O(H) −3
10 10
−4 −4
10 −14 −13 −12
10 −14 −13 −12
10 10 10 10 10 10
window size H [s] window size H [s]
Figures 6.13a and 6.13b show that for a small number of diodes (ND ≤ 10), i.e., with a
small Lipschitz constant LC , the speed of convergence is nearly the same for both computa-
tional sequences. However, for an increasing number of diodes (ND > 10) the convergence
speed is clearly superior when solving the circuit subproblem first: for ND = 100 diodes
the same level of accuracy is obtained with one iteration less and for ND = 1000 diodes one
saves nearly two iterations. Those results reflect the effect of the leading order coefficients
CD and CC as predicted in (5.71), cf. [9].
For further increasing values of ND (and thus LC ) the advantage of the device-first approach
is expected to improve further on and it will require fewer iterations than the reversed order
approach. This shows that a deeper knowledge on the strength of the coupling, i.e., a good
estimation of the Lipschitz constants (LD , LC ) helps to determine an optimal sequence for
solving the subsystems. Similarly it is known that reordering the computational sequence
can turn a splitting scheme with contraction factor α > 1 into a convergent setting, see
[5].
Conclusions
Finally global convergence and stability of the parameter coupling scheme is numerically
verified on the time interval I = [0, 2]·10−12s by performing simulations for decreasing time
window sizes, i.e., H = 10−12 , 10−13 , 10−14 . Simulation are computed for both sequences
(circuit-first and device-first), where only one Gauß-Seidel iteration (k = 1) is applied.
Time-integration is still obtained by the implicit Euler method. On one hand the additional
iterations decrease the splitting error by O(H), but on the other hand the implicit Euler
method is globally first order accurate O(h) = O(H) (step size equals window size) and
thus more iterations will not improve the total order of the scheme (when considering
splitting and time-discretization errors).
The relative errors (with respect to the monolithic reference solution using the same
time-stepping) are depicted in Fig. 6.14 for different window sizes. This is an order plot
for the convergence of the splitting scheme. For both sequences the convergence order is
approximately linear in the time window size H. Again, the larger Lipschitz-constants LC
98
6.4 Domain Substructuring of a Transformer
air
coil 1 iron core
(nonlinear) coil 2
(a) nonlinear, conductive region (iron), linear, non- (b) transformer model in 3D, eddy cur-
conductive regions: coils and air rents at the surface of the solid iron core
Figure 6.15: Domain substructuring example: (a) material regions, (b) field distribution in
3D, (c) region-wise degrees of freedom, see [36].
(for ND = 1000) causes the circuit-first approach to perform better than the device-first
approach (although both sequences are still of the same order). This underlines the results
of Theorem 5.15.
99
6 Numerical Examples
eig(C̃Mν C)
eig(λMσ )
Figure 6.16: Eigenvalues in substructuring example. The eigenvalues of the scaled conduc-
tivity matrix dominate the eigenvalues of the curl-curl matrix, time step h = 10−5 s, see also
Fig. 5.14a and [36].
is straightforward and it is available within FIDES (using odepkg, see Fig. 6.1). The
embedded Newton-Raphson scheme solves the linear problems directly by SuiteSparse
(Cholmod, [31]) or by PCG with a Jacobi preconditioner for the original, the regularized
and the Schur complement system (5.79).
The 2D and 3D problem formulations are structurally different: the block K22 is only
invertible in the 2D case where it corresponds to the discretization of the Laplacian. In
the 3D FIT discretization the same block must be regularized to allow for direct solvers. A
grad-div term was added in the non-conductive region as described in [36]. Furthermore it
will be shown that the sparsity patterns in 2D and 3D are different and thus direct solvers
behave differently. Let us start with the 2D problem.
100
6.4 Domain Substructuring of a Transformer
Table 6.2: Substructuring of the 2D model. Time transient simulation, 10 time steps, 48
Newton iterations, fixed step-size h = 10−3 . ‘iterations’ denotes the total number of all CG
iterations and ‘condition’ refers to the averaged condition number as approximated by PCG.
‘Cholmod’ and ‘PCG’ refer to direct/iterative solves of the full system, while ‘Schur PCG’
denotes the proposed algorithm from Listing 3 with Jacobi preconditioning, [36].
time iterations condition time iterations condition
Cholmod 37.0min - - Cholmod 55.3min - -
PCG 18.7min 13186 9.4 · 103 PCG 21.3min 14057 1.0 · 104
Schur PCG 6.5min 4083 7.4 · 102 Schur PCG 7.1min 4163 7.7 · 102
(a) 2D transformer with 333732 DoFs. (b) 2D transformer with 352990 DoFs (refined).
Table 6.3: Substructuring of the 3D model. Time transient simulation, 48 steps, step-sizes
h = 10−5 and 10−4 ; ‘iterations’ denotes the total number of all PCG iterations and ‘condition’
refers to the averaged condition number as approximated by PCG. PCG was used to solve
the curl-curl equation (‘PCG original’), its regularized version (‘PCG gauged’) and the Schur
system (5.79) (‘Schur method’); regularization is carried out by an edge-wise grad-div term,
[36].
time iterations condition time iterations condition
PCG original 12.3min 13423 1.4 · 104 PCG original 17.8min 19019 3.0 · 104
PCG gauged 49.8min 38425 1.2 · 105 PCG gauged 76.1min 56935 1.5 · 105
Schur method 7.2min 815 1.9 · 101 Schur method 20.6min 2456 1.4 · 102
(a) 3D example with step size h = 10−5 (b) 3D example with step size h = 10−4
The time step pushes the eigenvalues of the conductivity matrix further along the positive
axis, while the non-conductive ones remain unaltered. Thus especially in the case of tiny
step sizes, the eigenvalues of the conductivity matrix dominate the spectrum.
Table 6.3 shows the total simulation time, the number of PCG iterations and the es-
timated condition number for (i) the original system, (ii) the system with applied grad-
div regularization and finally (iii) the Schur system Listing 3. Again, the Schur com-
plement method reduces the number of PCG iterations substantially, but for large step
sizes this does not compensate for the higher computational costs due to the additional
forward/backward substitutions. The sparsity pattern is not well preserved by Cholmod
and thus the matrix/vector operations increase the averaged time for a single PCG itera-
tion from 0.05s to 0.5s (other direct solvers using other ordering strategies might perform
better3 ).
Only for small step sizes h ≤ 10−5 the gap between the eigenvalues becomes dominant
and the fewer number of iterations compensates for the additional costs per iteration.
In 3D, the direct solver fails for the singular matrix pencil and when applied to the
regularized system, the memory requirements are so huge that the simulation takes an
unreasonable time of ten hours.
3
Cholmod uses METIS for ordering, http://www-users.cs.umn.edu/~karypis/metis
101
6 Numerical Examples
6.4.3 Conclusions
The 2D and 3D computational examples have demonstrated the general feasibility of the
Schur complement method. The direct linear system solvers perform very well for (small
scale) 2D problems, but they become too expensive for large 3D problems. The usage
of low-rank approximations, model order reduction techniques and inner iterative solvers
with deflation is expected to further improve the proposed method. Combinations of the
method with other more sophisticated preconditioners or multigrid solvers is the subject
of current research, [115].
102
7 Conclusions
In this thesis multiscale models from electrical engineering, i.e., the lumped electric net-
work, the magnetoquasistatic field device and the drift-diffusion model of a semiconductor
device were derived, analyzed and mutually coupled. This culminates in a system of partial
differential algebraic equations (PDAEs) that may contain device models of any dimension
(from 0D to 3D) interconnected by the surrounding electric network.
For the spatial discretized PDAE system new differential-algebraic-index results in terms
of the tractability concept were obtained: in particular it was shown that the field/circuit
coupled problem is at most index-2 and that the index-2 variables enter the system only
linearly. Thus the error propagation during time-integration is not seriously affected.
Numerical test examples have underlined this result and show that the index-2 case is
rather harmless but nonetheless the time-integration error can be decreased when modeling
the circuit topology according to the index-1 results.
Two different classes of multirate time-integration methods for the given PDAE problem
were discussed. They both follow a hierarchical multiscale approach where lumped models
(for example generalized inductances) are fitted on the fly by spatially distributed PDE
models (for example magnetoquasistatic Maxwell’s Equations).
This is in particular accomplished by the multirate Schur complement approach in which
the PDE device equations are eliminated from the circuit equations such that only a lumped
model remains. The lumped model is updated according to an energy-based scheme such
that the problem-specific multirate potential is exploited. Furthermore the complement
makes multimethod approaches feasible, i.e., different linear solvers can be used for the
subproblems. The acceleration of the time-integration is significant: a numerical example
shows a 40× speed-up compared to a standard approach.
Secondly the multirate cosimulation (using dynamic iteration) was adapted for this hier-
archical approach. For that purpose we introduced a parameter coupling interface that uses
lumped surrogate models. This interface was shown to be superior to the common source
coupling approach. Again, this unlocks multirate potential and allows even the application
of different time-integrators for each subproblem – with the drawback of an increased com-
putational overhead due to additional iterations. Nonetheless numerical examples from
field/circuit coupling show a clear reduction in the computational costs compared to a
monolithic single rate approach.
Furthermore, it was mathematically proved that the dynamic iteration applied to both
field/circuit and semiconductor/circuit coupled problems is always convergent if the cou-
pling interface is modeled by the parameter approach. To this end we carried out a fixed
point analysis in function space and a propagation error analysis. This allowed us to com-
pare the different interface models and we discussed their advantages and disadvantages.
The underlying principles were given within an abstract framework such that this approach
can be easily applied to other fields of application.
103
7 Conclusions
7.1 Outlook
This treatise raises new interesting questions and problems: the DAE-index analysis and
multirate cosimulation of a coupled system consisting of all three subproblems, i.e., field,
circuit and semiconductor seems to be the inevitable next step. A tentative analysis reveals
that the interaction of coupling interface, number of subsystems and convergence rate offers
interesting new aspects. Along with the analysis, numerical experiments should be carried
out for this enlarged system.
Another interesting aspect is the relation of the order of the time-integrator and the
overall cosimulation. The results in Section 6.3.1 imply that higher order cosimulation
requires additional iteration. This is clear for simulations where time window and step
size coincide, but what happens in the general (multirate) case? If this is well understood
and one knows the approximation order of each iterate, this will allow us to predict the
splitting error and consequently to derive a very efficient control for the window size and
the number of iterations.
104
A Discretization Properties
The three indices ix , iy and iz , given nx , ny and nz ∈ N, are combined into one space
index, which allows us to number the objects consecutively:
105
A Discretization Properties
primary
point
-1
-1 1
dual
cell
dual
primary point
cell 1
(a) primary and dual cell (b) primary curl
where the partial differential operators Pw ∈ {−1, 0, 1}n0 ×n0 for each spatial direction
w ∈ {x, y, z} are defined as
−1 for q = p,
(Pw )p,q := δp+kw ,q − δp,q = 1 for q = p + kw , (A.2)
0 else;
where I denotes the identity matrix and kw is defined as in (A.1). For inner edges each par-
tial differential operator exhibits exactly two non-zeros (1 and −1) per row. Consequently
the primary curl matrix has four entries (−1, −1, 1, 1) in each row which pick out the
corresponding line-integrals from a vector that lives on the primary edge, see Fig. A.1b.
Lemma A.1. Two discrete partial differential operators Pv and Pw with v, w ∈ {x, y, z}
commute
Pv Pw = Pw Pv . (A.3)
Proof. Let Dp,q := δp+1,q denote a shift matrix then Pw = Dkw − I and
106
A.1 Discrete Operators
z
y
x
ex(7)
)
)
(6
(5
ey
ey
ex(5)
ez(4)
ez(3)
ez(2)
ez(1)
ex(3)
)
(1
(2
ey
ey
ex(1)
The result of Lemma A.1 corresponds to the interchange of partial derivatives in the
continuous case. A similar proof is given for example in [139, 12].
Theorem A.2. The product of the discrete divergence and the discrete curl matrices is
zero
S∞ C∞ = 0 (A.4)
S̃∞ C̃∞ = 0. (A.5)
Proof. Both relations are simple applications of Lemma A.1. For example we find for the
operators on the primary grid
S∞ C∞ = Py Pz − Pz Py Pz Px − Px Pz Px Py − Px Py = 0
which is zero because the Pw commute. The dual case is analogous. This proof and its
consequences are discussed in several publications, e.g. [12].
Example A.1. Let us examine the (primary) FIT grid of dimensions 2 × 2 × 2 as shown
in Fig. A.2. It describes only one complete volume with 6 facets and 12 edges, but another
7 volumes, 18 facets and 12 edges are superfluously included in the numbering scheme.
Hw = {1 ≤ n(ix , iy , iz ) ≤ n0 |iw = nw }
107
A Discretization Properties
We denote the set of all points connected to at least one phantom edge by Hxyz := Hx ∪
Hy ∪ Hz and define Hxy , Hxz , Hyz accordingly. Each index set Hw with w ∈ {x, y, z} gives
raise to a diagonal matrix Lw ∈ {0, 1}n0×n0 with
(
1 for p = q and p ∈
/ Hw ,
(Lw )pq = (A.6)
0 else.
Lemma A.3. The matrices Lw are projectors, i.e., L2w = Lw and have the following
properties for two directions v, w ∈ {x, y, z} with v 6= w
Lw Lv = Lv Lw , (A.7)
Lw Lv Pv = Lv Pv Lw , (A.8)
Proof. The idem-potency of the projectors and their commuting is trivial, since they are
diagonal matrices containing only ones and zeros. The left-hand-side of the third equa-
tion (A.8) reads
−1 for p = q ∧ p ∈
/ Hw ∪ Hv ,
(Lw Lv Pv )pq = 1 for q = p + kv ∧ p ∈
/ Hw ∪ Hv ,
0 else.
because p ∈
/ Hv implies that
p = ix + iy ky + iz kz with iv < nv
Only the edges in the set Hw exist in the grid, i.e., are degrees of freedom. We define
analogously that facets and volumes exist if none of their edges are phantoms. We assemble
the projectors corresponding to the index sets in x, y z-order, such that
IP = IṼ := I
IL = IÃ := blkdiag(Lx , Ly , Lz )
IA = IL̃ := blkdiag(Ly Lz , Lx Lz , Lx Ly )
IV = IP̃ := Lx Ly Lz
108
A.1 Discrete Operators
Figure A.3: The number of total and non-phantom objects in Maxwell’s house
where IP , IL , IA and IV denote the projectors for all points, edges, facets and volumes
in the primary grid; the ones with a tilde are the corresponding counterparts on the dual
grid, see Fig. A.3. These definitions match the matrices in [45, Section 2.2.2].
Corollary A.4. The matrices IP , IL , IA and IV are projectors.
Projected Operators
The curl operator relates the edges to the flux through their facet and therefore we ignore
contributions from phantom edges and facets. Using the matrices from above we can apply
the curl-operator on finite domains by defining (cf. [45, Section 2.2.2])
The divergence operators are a mapping between facets and volumes; the primary and
dual operator read using the projectors
Remark A.1. The curl operator C∞ was constructed from only three partial differential
matrices Pw with w ∈ {x, y, z}, but the finite operator C in equation (A.9) is constructed
from 6 pairwise distinct blocks, e.g.
Ly Lz Py Lz 6= Lx Ly Py Lx (A.11)
and those blocks differ again from the blocks in both, the primary and dual divergence
operators S and S̃
Lx Ly Lz Py Lx Lz 6= Py Lx Ly Lz (A.12)
Nonethless the blocks have a redundancy in the projectors, which is revealed by the
following corollary:
Corollary A.5. Let S denote the primary divergence operator as defined in equa-
tion (A.10) and C the primary curl operator from equation (A.9), then
C := IA C∞ IL = IA C∞
109
A Discretization Properties
S := IV S∞ IA = IV S∞
Proof. This is just a consequence of Lemma A.3. We show here only the second equation;
it holds
S = IV S∞ IA = Lx Ly Lz Px Ly Lz Py Lx Lz Pz Lx Ly = Lx Ly Lz Px Py Pz = IV S∞
because the projectors on the left commute, so they can be reordered. Then the operators
commute with the projectors on their right and finally the projectors are idem-potent.
Remark A.2. Sometimes in the literature, the partial differential matrices are directly
constructed for finite domains, for example
−1 for q = p and j ≤ J − 1,
(P̄y )pq = +1 for q = p + ky and j ≤ J − 1,
0 else;
P̄y = Ly Py
in our notation and it unveils that the boundary constraints of P̄y are not sufficient,
compare equations (A.11) and (A.12).
SC = 0,
S̃C̃ = 0.
Proof. We show the first equation using the results from Corollary A.5 for S and C, we
obtain directly
SC = IV S∞ C∞ IL = 0
because S∞ C∞ = 0 holds as a stated in Theorem A.2. The dual case can be proved
analogously.
Curl-Curl Matrix
The projectors apply also naturally to material matrices. In this case the consistency of
the discrete problem is ensured, since material properties might be assigned to phantom
110
A.2 Material Matrices
objects, e.g. if all diagonal entries of the reluctivity matrix are positive. Let us discuss in
the following the curl-curl matrix as it is used in the vector-potential formulations.
Corollary A.7. For the curl-curl matrix it is equivalent to either impose the boundary
conditions on the material matrix only or on the curl matrices only
Proof. Starting from the curl-curl matrix expressed in the infinite operators with the ma-
terial matrix Mν = IL̃ Mν,∞ IA mapping from primary facets to dual edges
which shows the equivalence in imposing the boundary condition in both the curl matrices
and material matrix Mν , or just in one of them.
Conclusion
The phantom edges do not belong to the degrees of freedom. It was shown that they can
be disregarded in the construction of the differential operators, if the material matrices are
assembled correctly.
where p refers to the dimension of the object, e denotes the element and 1i a column vector
with a 1 at the i-th position and 0 otherwise (the non-negativity follows from Ass. A.8).
The length of the vectors is given by the overall number of either points (n0 ), edges (n1 )
or facets (n2 ). The local numbers are denoted by mp , see Table A.4a. The index mappings
πp,e are injective embeddings identifying the local to its global number
i.e., full column rank. For example in the 3D FIT case the incidence matrix Q2,e for the
facets (Fig. A.4b), is constructed by the index mappings
111
A Discretization Properties
where kx , ky and kz are defined as in (A.1) and n1 is the global number of edges.
Assumption A.9 (No phantom objects). It is assumed that degrees of freedom are not
located on phantom objects, i.e., that each p-dimensional object belongs to one or more grid
elements, see Section A.1.3
∀ip ∃(e, jp ) πp,e (jp ) = ip , (A.13)
where e refers to an element, ip and jp are the global and local number of the same p-
dimensional object.
where the diagonal 3 × 3-matrices D1,e and D2,e contain the primary edge lengths and
primary facets areas of the element ⌢
e. The ⌢local area-integrated magnetic flux densities
⌢ ⌢
(i.e., magnetic fluxes) are given by be := Q⊤2,e b. Allocation of the fluxes at the facet centers
and averaging the opposing ones yields the local flux density Be in the element center
1 ⌢
⌢
Be := D−1
2,e ⊗ 1 1 be and ||Be||2 = (Be )21 + (Be )22 + (Be )23
2
112
A.2 Material Matrices
is the corresponding magnitude of the magnetic flux. The definition (A.14) yields for
(r)
diagonal material tensors, i.e., νaniso = 0, diagonal and positive definite matrices. In this
case it coincides with the classical FIT approach. On the other hand the definition above
preserves symmetry even for anisotropic materials, which is in contrast to other approaches,
[99]:
Lemma A.10. The matrix Mν,e is symmetric positive definite if the material tensor ν is
so.
Proof. For the hexahedral FIT-discretization the volume matrix of each element is given
by D3,e = D1,e D2,e = |Ωe | · I; rewriting Mν,e in the form T⊤ ν (r) T yields symmetry by the
assumption above (with a suitable matrix T). On the other hand positive semi-definiteness
follows from the positive semi-definiteness of each Kronecker product and finally the full
rank of the matrix sum is proved by Sylvester’s inequality (the product of both summands
is zero).
and finally the global matrices are defined in the obvious way
X X ⌢
⌢
X ⌢
⌢
Mε := M(r)
ε , Mσ := M(r)
σ , Mν (b) := M(r)
ν (b) , (A.15)
i i i
with symmetric positive definite local matrices Mε,e and Mν,e (as shown above for the
reluctivity matrix). The matrix Mσ,e is in general only symmetric positive semi-definite
due to non-conducting areas (e.g. due to air). In the special case of FIT the construction
yields diagonal global material matrices if only isotropic materials or anisotropic materials
with principal directions coinciding with the grid directions are present.
In all cases the symmetry (diagonality) of the local matrices obviously carries over to
the global material matrix. Even the definiteness is ensured since each object (edge, facet)
is assigned at least one material parameter, mathematically speaking:
Lemma A.11. The global matrix Mξ with ξ ∈ {ε, σ, ν} is positive definite if all element
contributions Mξ,e are so; it is positive semi-definite if Mξ,e are only semi-definite.
Proof. Let p denote edges or facets depending on the material property ξ and let there
be a vector x 6= 0 ∈ RNp and ip denote the index of a non-zero degree of freedom. Then
follows from (A.13) that there is an element e such that y := Qp,e x 6= 0 and hence the
corresponding summand is positive y⊤ Mξ,e y > 0. Therefore the whole sum is positive,
since all other summands are non-negative. The semi-definite case is trivial.
A.2.3 Nonlinearity
In real world applications the material properties depend on additional conditions for
example due to nonlinear material relations as in the case of the reluctivity matrix in
113
A Discretization Properties
(A.15). In this case the additional dependence has to be introduced into the global material
and local matrices. This yields matrix-valued functions instead of the constant matrices
above. The following analysis applies to nonlinear reluctivities, but also to temperature-
dependent conductivities: let be given a positive scalar valued function f (e.g., f = νiso )
f : Rnp → R>0
whose argument αe can be any quantity located at the points, edges, facets or centers of
primary element (for example the averaged flux magnitude ||Be ||2 is located at the center).
Although the function itself depends nonlinearly on the parameter, its value affects the
local element matrix only affine linearly
X
Mξ (αe ) := Qp,e f (αe )Mξ,e Q⊤p,e . (A.16)
e
where p = {1, 2} is chosen accordingly to the material, as defined above. The parameter-
dependent material matrix has the same properties as before: it is still symmetric (diagonal
if applicable) and positive (semi-)definite, with a kernel that does not depend on the
parameter (the image is R>0 ).
Lemma A.12. The kernel of Mξ (αe ) as defined in (A.16) is constant, i.e., independently
of the parameter αe .
Proof. Let αe and αe∗ denote two distinct parameters and x ∈ Ker Mξ (αe ), then follows
Ker Mξ (αe ) ⊂ Ker Mξ (αe∗ ) from the scalar nature of fe
X
0 = x⊤ Mξ (αe )x = fe (αe )x⊤ Qp,e Mξ,e Q⊤p,e x
X e
= ∗ ⊤
fe (αe )x Qp,e Mξ,e Q⊤p,e x = x⊤ Mξ (αe∗ )x,
e
because all summands are necessarily zero, independently of fe > 0. Now using the same
arguments for x∗ ∈ Ker Mξ (αe∗ ) yields Ker Mξ (αe ) = Ker Mξ (αe∗ ).
The arguments above are only valid for effects that have an isotropic (scalar valued)
impact on the material and do not force the material property to vanish (e.g. switches).
Typical applications are in the modeling of temperature dependent conductivities or sat-
uration effects of the reluctivities.
dkν (⌢
a) d ⌢ ⌢
Kν = = C̃Mν (C a)C a + Zσ
d⌢a d⌢a
d ⌢ ⌢
⌢ ⌢
d ⌢
⌢ db
⌢ ⌢
⌢
⌢
114
A.3 Differential Curl-Curl Matrix
⌢
⌢
and due to the linearity the derivative with respect to b, it can be passed through the
sums in equation (A.15) right to the nonlinear isotropic material tensor and thus
(r) ⌢
dν (||Be ||2 )
(r) ⌢
⌢ (r) ⌢⌢ 1 0 1 0 ⌢⊤
⌢ ⌢
νFIT,d (be ) := νFIT (be ) + iso D(r)
ν ⊗ 0 1
D−1
2,e ⊗ 0 1
be be,av ,
d||Be ||2
In conclusion, the differential reluctivity (and thus the differential curl-curl matrix) are
assembled by using the differential reluctivity tensor νFIT,d instead of the chord reluctivity
tensor νFIT in (A.14).
⌢
⌢ 1 1 0 (r) ⌢
⌢
−1 1 0
Mν,d,e (be ) := D1,e ⊗ 0 1 νFIT,d (be ) D2,e ⊗ 0 1 (A.17)
4
The resulting matrix can be nonsymmetric depending on the construction, [41]. However,
for physical correct material curves positive definiteness can be shown.
Corollary A.13. For Brauer’s model the differential reluctivity matrix is positive definite.
This proves that the differential reluctivity material tensor is positive definite and thus the
material matrix is positive definite by Lemma A.11.
115
Bibliography
[1] Giuseppe Alı̀, Andreas Bartel, Markus Brunk, and Sebastian Schöps. A convergent
”
iteration scheme for semiconductor/circuit coupled problems“. In: Scientific Com-
puting in Electrical Engineering SCEE 2010. Ed. by Bas Michielsen. Mathematics
in Industry. To be published. Berlin: Springer, 2011. isbn: 978-3-540-71979-3.
[2] Giuseppe Alı̀, Andreas Bartel, and Michael Günther. Parabolic differential-
”
algebraic models in electrical network design“. In: SIAM Journal on Multiscale
Modeling and Simulation 4.3 (2005), pp. 813 –838. doi: 10.1137/040610696.
[3] Ana Alonso and Alberto Valli. A domain decomposition approach for hetero-
”
geneous time-harmonic Maxwell equations“. In: Computer Methods in Applied
Mechanics and Engineering 143.1-2 (1997), pp. 97 –112. issn: 0045-7825. doi:
10.1016/S0045-7825(96)01144-9.
[4] Angelo Marcello Anile, Giuseppe Alı̀, and Giovanni Mascali, eds. Scientific Com-
puting in Electrical Engineering SCEE 2004. Mathematics in Industry 9. Berlin:
Springer, 2006.
[5] Martin Arnold and Michael Günther. Preconditioned Dynamic Iteration for Cou-
”
pled Differential-Algebraic Systems“. In: BIT Numerical Mathematics 41.1 (2001),
pp. 1 –25. doi: 10.1023/A:1021909032551.
[6] Martin Arnold, Karl Strehmel, and Rüdiger Weiner. Errors in the numerical solu-
tion of nonlinear differential-algebraic systems of index 2. Tech. rep. Martin-Luther-
University Halle, 1995.
[7] Andreas Bartel. Partial Differential-Algebraic Models in Chip Design - Thermal and
Semiconductor Problems. Fortschritt-Berichte VDI, Reihe 20. Dissertation. Düssel-
dorf: VDI Verlag, 2004.
[8] Andreas Bartel, Sascha Baumanns, and Sebastian Schöps. Structural analysis
”
of electrical circuits including magnetoquasistatic devices“. In: Applied Numerical
Mathematics (2010). Submitted.
[9] Andreas Bartel, Markus Brunk, Michael Günther, and Sebastian Schöps. Dynamic
”
iteration schemes for coupled problems in circuit/semiconductor simulation“. In
preparation. 2011.
[10] Andreas Bartel, Michael Günther, and Anne Kværnø. Multirate Methods in Elec-
”
trical Circuit Simulation“. In: Progress in Industrial Mathematics at ECMI 2000.
Ed. by Angelo Marcello Anile, Vincenzo Capasso, and Antonio M. Greco. Berlin:
Springer, 2002, pp. 258 –265.
116
Bibliography
[11] Andreas Bartel and Roland Pulch. A Concept for Classification of Partial
”
Differential Algebraic Equations in Nanoelectronics“. In: Progress in Industrial
Mathematics at ECMI 2006. Ed. by Luis L. Bonilla, Jose M. Vega, Miguel
Moscoso, and Gloria Platero. Berlin: Springer, 2006. isbn: 978-3-540-71991-5. doi:
10.1007/978-3-540-71992-2_79.
[12] Michael Bartsch et al. Solution of Maxwell’s Equations“. In: Computer Physics
”
Communications 73.1-3 (1992), pp. 22 –39. doi: 10.1016/0010-4655(92)90026-U.
[13] Sascha Baumanns, Monica Selva Soto, and Caren Tischendorf. Consistent initial-
”
ization for coupled circuit-device simulation“. In: Scientific Computing in Electrical
Engineering SCEE 2008. Ed. by Janne Roos and Luis R. J. Costa. Mathematics in
Industry 14. Springer, 2010, pp. 297 –304.
[14] Gary Bedrosian. A New Method for Coupling Finite Element Field Solutions
”
with External Circuits and Kinematics“. In: IEEE Transactions on Magnetics 29.2
(1993), pp. 1664 –1668. doi: 10.1109/20.250726.
[15] Galina Benderskaya. Numerical Methods for Transient Field-Circuit Coupled Sim-
”
ulations Based on the Finite Integration Technique and a Mixed Circuit Formula-
tion“. PhD thesis. TU Darmstadt, 2006.
[16] Michele Benzi, Gene H. Golub, and Jörg Liesen. Numerical solution of
”
saddle point problems“. In: Acta Numerica 14.1 (2005), pp. 1–137. doi:
10.1017/S0962492904000212.
[17] Oszkár Bı́ró and Kurt Preis. On the use of the magnetic vector potential in the
”
finite-element analysis of three-dimensional eddy currents“. In: IEEE Transactions
on Magnetics 25.4 (July 1989), pp. 3145 –3159. doi: 10.1109/20.34388.
[18] Oszkár Bı́ró, Kurt Preis, and Kurt R. Richter. Various FEM formulations for the
”
calculation of transient 3D eddy currents in nonlinear media“. In: IEEE Transac-
tions on Magnetics 31.3 (May 1995), pp. 1307–1312. doi: 10.1109/20.376269.
[19] Alain Bossavit. Computational Electromagnetism: Variational Formulations, Com-
plementarity, Edge Elements. San Diego: Academic Press, 1998. isbn: 0121187101.
[20] Alain Bossavit. Stiff problems in eddy-current theory and the regularization of
”
Maxwell’s equations“. In: IEEE Transactions on Magnetics 37.5 (Sept. 2001),
pp. 3542–3545. issn: 0018-9464. doi: 00189464/01$10.00.
[21] Alain Bossavit. Whitney forms: a class of finite elements for three-dimensional
”
computations in electromagnetism“. In: IEE Proceedings 135.8 (1988), pp. 493–
500. doi: 10.1049/ip-a-1:19880077.
[22] Alain Bossavit and Lauri Kettunen. Yee-like schemes on staggered cellular grids: a
”
synthesis between FIT and FEM approaches“. In: IEEE Transactions on Magnetics
36.4 (July 2000), pp. 861–867. doi: 10.1109/20.877580.
[23] Hans Georg Brachtendorf, G. Welsch, R. Laur, and A. Bunse-Gerstner. Numer-
”
ical steady state analysis of electronic circuits driven by multi-tone signals“. In:
Electrical Engineering (Archiv für Elektrotechnik) 79.2 (1996), pp. 103–112. doi:
10.1007/BF01232919.
117
Bibliography
118
Bibliography
[38] Herbert De Gersem, Markus Clemens, and Thomas Weiland. Coupled finite-
”
element, spectral-element discretisation for models with circular inclusions and far-
field domains“. In: IET Science, Measurement & Technology 149.5 (2002), pp. 237
–241. doi: 10.1049/ip-smt:20020620.
[39] Herbert De Gersem and Kay Hameyer. A Finite Element Model for Foil Winding
”
Simulation“. In: IEEE Transactions on Magnetics 37.5 (Sept. 2001), pp. 3472 –
3432. doi: 10.1109/20.952629.
[40] Herbert De Gersem, Kay Hameyer, and Thomas Weiland. Field-Circuit Coupled
”
Models in Electromagnetic Simulation“. In: Journal of Computational and Applied
Mathematics 168.1-2 (2004), pp. 125 –133. doi: 10.1016/j.cam.2003.05.008.
[41] Herbert De Gersem, Irina Munteanu, and Thomas Weiland. Construction of Differ-
”
ential Material Matrices for the Orthogonal Finite-Integration Technique With Non-
linear Materials“. In: IEEE Transactions on Magnetics 44.6 (June 2008), pp. 710–
713. doi: 10.1109/TMAG.2007.915819.
[42] Herbert De Gersem, Irina Munteanu, and Thomas Weiland. Differential material
”
matrices for the Finite Integration Technique“. In: The European Physical Journal-
Applied Physics 39 (2007), pp. 165–169. doi: 10.1051/epjap:2007119.
[43] Herbert De Gersem and Thomas Weiland. Field-Circuit Coupling for Time-
”
Harmonic Models Discretized by the Finite Integration Technique“. In:
IEEE Transactions on Magnetics 40.2 (Mar. 2004), pp. 1334 –1337. doi:
10.1109/TMAG.2004.824536.
[44] Georg Denk. Circuit Simulation for Nanoelectronics“. In: Scientific Computing in
”
Electrical Engineering SCEE 2004. Ed. by Angelo Marcello Anile, Giuseppe Alı̀,
and Giovanni Mascali. Mathematics in Industry 9. Berlin: Springer, 2006. doi:
10.1007/978-3-540-32862-9_2.
[45] Martin Dohlus. Ein Beitrag zur numerischen Berechnung elektromagnetischer
”
Felder im Zeitbereich“. German. Dissertation. TU Darmstadt, 1992.
[46] Thomas Dreher and Gérard Meunier. 3D Line Current Model of Coils and External
”
Circuits“. In: IEEE Transactions on Magnetics 31.3 (May 1995), pp. 1853 –1856.
doi: 10.1109/20.376398.
[47] Johan Driesen, Ronnie J. M. Belmans, and Kay Hameyer. Methodologies for cou-
”
pled transient electromagnetic-thermal finite-element modeling of electrical energy
transducers“. In: IEEE Transactions on Industry Applications 38.5 (Sept. 2002),
pp. 1244–1250. doi: 10.1109/TIA.2002.803024.
[48] Patrick Dular and Johan Gyselinck. Modeling of 3-D stranded inductors with the
”
magnetic vector potential formulation and spatially dependent turn voltages of re-
duced support“. In: IEEE Transactions on Magnetics 40.2 (Mar. 2004), pp. 1298
–1301. doi: 10.1109/TMAG.2004.825153.
[49] Patrick Dular, François Henrotte, and Willy Legros. A general and natural
”
method to define circuit relations associated with magnetic vector potential for-
mulations“. In: IEEE Transactions on Magnetics 35.3 (May 1999), pp. 1630–1633.
doi: 10.1109/20.767310.
119
Bibliography
[50] Falk Ebert. Convergence of relaxation methods for coupled systems of ODEs and
DAEs. Tech. rep. 200. MATHEON, 2004.
[51] Falk Ebert. On Partitioned Simulation of Electrical Circuits using Dynamic Iter-
”
ation Methods“. PhD thesis. Berlin: TU Berlin, 2008.
[52] Jocelyne Erhel and Frederic Guyomarc’h. An Augmented Conjugate Gradient
”
Method for Solving Consecutive Symmetric Positive Definite Linear Systems“. In:
SIAM Journal on Matrix Analysis and Applications 21.4 (2000), pp. 1279–1299.
doi: 10.1137/S0895479897330194.
[53] Diana Estévez Schwarz. Consistent initialization for index-2 differential algebraic
”
equations and its application to circuit simulation“. PhD thesis. Berlin: Humboldt-
Universität, 2000.
[54] Diana Estévez Schwarz and Caren Tischendorf. Structural analysis of
”
electric circuits and consequences for MNA“. In: International Journal
of Circuit Theory and Applications 28.2 (2000), pp. 131 –162. doi:
10.1002/(SICI)1097-007X(200003/04)28:2<131::AID-CTA100>3.0.CO;2-W.
[55] Carlo de Falco. The IFF File Format Specification. CoMSON Project. Wuppertal,
2008.
[56] Carlo de Falco, Georg Denk, and Reinhart Schultz. A Demonstrator Platform
”
for Coupled Multiscale Simulation“. In: Scientific Computing in Electrical En-
gineering SCEE 2006. Ed. by Gabriela Ciuprina and Daniel Ioan. Mathematics
in Industry 11. Berlin: Springer, 2007, pp. 63–71. isbn: 978-3-540-71979-3. doi:
10.1007/978-3-540-71980-9_5.
[57] Uwe Feldmann and Michael Günther. CAD-based electric-circuit modeling in in-
”
dustry I: mathematical structure and index of network equations“. In: Surveys on
Mathematics for Industry 8.2 (1999), pp. 97 –129.
[58] Uwe Feldmann, Masataka Miyake, Takahiro Kajiwara, and Mitiko Miura-
Mattausch. On Local Handling of Inner Equations in Compact Models“. In: Sci-
”
entific Computing in Electrical Engineering SCEE 2008. Ed. by Janne Roos and
Luis R. J. Costa. Mathematics in Industry 14. Springer, 2010, pp. 143–150.
[59] Martin J. Gander and Albert E. Ruehli. Optimized waveform relaxation methods
”
for RC type circuits“. In: IEEE Transactions on Circuits and Systems 51.4 (Apr.
2004), pp. 755 –768. issn: 1549-8328. doi: 10.1109/TCSI.2004.826193.
[60] Charles William Gear and D. R. Wells. Multirate Linear Multistep Methods“. In:
”
BIT Numerical Mathematics 24.4 (1984), pp. 484 –502. doi: 10.1007/BF01934907.
[61] Gene H. Golub and Charles F. van Loan. Matrix Computations. 3rd ed. The Johns
Hopkins University Press, 1996.
[62] Klaus-Tibor Grasser. Mixed-Mode Device Simulation“. PhD thesis. TU Wien,
”
1999.
[63] Eberhard Griepentrog and Roswitha März. Differential-Algebraic Equations and
Their Numerical Treatment. Teubner, Leipzig, 1986.
[64] Michael Günther and Uwe Feldmann. The DAE-Index in Electric Circuit Simu-
”
lation“. In: Mathematics and Computers in Simulation 39 (Nov. 1995), pp. 573 –
582.
120
Bibliography
[65] Michael Günther, Uwe Feldmann, and Jan ter Maten. In: Numerical Methods in
Electromagnetics. Ed. by Wil Schilders, Philippe Ciarlet, and Jan ter Maten. Vol. 13.
Handbook of Numerical Analysis. North-Holland, 2005. Chap. Modelling and Dis-
cretization of Circuit Problems. isbn: 978-0-444-51375-5.
[66] Wolfgang Hackbusch. Direct Domain Decomposition using the Hierarchical Ma-
”
trix Technique“. In: Fourteenth International Conference on Domain Decomposi-
tion Methods. Ed. by Ismael Herrera, David E. Keyes, Olof B. Widlund, and Robert
Yates. 2003.
[67] Ernst Hairer, Christian Lubich, and Michel Roche. The Numerical Solu-
tion of Differential-Algebraic Systems by Runge-Kutta Methods. Lecture Notes
in Mathematics. Berlin: Springer, 1989. isbn: ISBN 3-540-51860-6. doi:
10.1007/BFb0093947.
[68] Ernst Hairer, Syvert P. Nørsett, and Gerhard Wanner. Solving Ordinary Differential
Equations II: Stiff and Differential-Algebraic Problems. 2nd ed. Springer Series in
Computational Mathematics. Berlin: Springer, 2002. isbn: 978-3-540-604525-0.
[69] Ernst Hairer, Syvert P. Nørsett, and Gerhard Wanner. Solving Ordinary Differ-
ential Equations I: Nonstiff Problems. 2nd ed. Springer Series in Computational
Mathematics. Berlin: Springer, 2000. isbn: 978-3-540-56670-0.
[70] A.Y. Hannalla and D.C. MacDonald. Numerical analysis of transient field problems
”
in electrical machines“. In: Proceedings of the Institution of Electrical Engineers
123.9 (Sept. 1976), pp. 893–898. issn: 0020-3270. doi: 10.1049/piee.1976.0190.
[71] Hermann A. Haus and James R. Melcher. Electromagnetic Fields and Energy. Pren-
tice Hall, 1989. isbn: 013249020X.
[72] Bodo Heise. Analysis of a Fully Discrete Finite Element Method for a Nonlinear
”
Magnetic Field Problem“. In: SIAM Journal on Numerical Analysis 31.3 (1994),
pp. 745–759.
[73] Ralf Hiptmair and Oliver Sterz. Current and voltage excitations for the eddy cur-
”
rent model“. In: International Journal of Numerical Modelling: Electronic Networks,
Devices and Fields 18.1 (2005), pp. 1–21. doi: 10.1002/jnm.555.
[74] Chung-Wen Ho, Albert E. Ruehli, and Pierce A. Brennan. The modified nodal
”
approach to network analysis“. In: IEEE Transactions on Circuits and Systems
22.6 (June 1975), pp. 504 –509. issn: 0098-4094. doi: 10.1109/TCS.1975.1084079.
[75] Daniel Ioan and Irina Munteanu. Missing link rediscovered: the electromagnetic
”
circuit element concept“. In: The Japan Society of Applied Electromagnetics and
Mechanics 8 (1998), pp. 302–320.
[76] Zdzislaw Jackiewicz and Marian Kwapisz. Convergence of Waveform Relaxation
”
Methods for Differential-Algebraic Systems“. In: SIAM Journal on Numerical Anal-
ysis 33.6 (Dec. 1996), pp. 2303–2317. doi: 10.1137/S0036142992233098.
[77] Jan Janssen and Stefan Vandewalle. On SOR Waveform Relaxation Methods“. In:
”
SIAM Journal on Numerical Analysis 34 (1997), pp. 2456–2481.
[78] Enrique Francisco Kaasschieter. Preconditioned conjugate gradients for solving
”
singular systems“. In: Journal of Computational and Applied Mathematics 24.1-2
(1988), pp. 265 –275. doi: 10.1016/0377-0427(88)90358-5.
121
Bibliography
122
Bibliography
[92] Roswitha März. Solvability of linear differential algebraic equations with properly
”
stated leading terms“. In: Results in Mathematics 45.1-2 (2004). Preprint, pp. 88–
105.
[93] James Clerk Maxwell. A Dynamical Theory of the Electromagnetic Field“. In:
”
Royal Society Transactions CLV (1864), pp. 459 –512.
[94] David Meeker. Finite Element Method Magnetics User’s Manual. Version 4.2
(09Nov2010 Build). 2010.
[95] Ronny Mertens et al. An algebraic multigrid method for solving very large electro-
”
magnetic systems“. In: IEEE Transactions on Magnetics 34.5 (Sept. 1998), pp. 3327
–3330. doi: 10.1109/20.717782.
[96] Bas Michielsen, ed. Scientific Computing in Electrical Engineering SCEE 2010.
Mathematics in Industry. To be published. Berlin: Springer, 2011. isbn: 978-3-540-
71979-3.
[97] Ulla Miekkala and Olavi Nevanlinna. Convergence of dynamic iteration methods
”
for initial value problems“. In: SIAM Journal on Scientific and Statistical Comput-
ing 8 (July 1987), pp. 459–482. doi: 10.1137/0908046.
[98] Kenneth Moreland. The ParaView Tutorial. 3.8. Sandia National Laboratories.
2010.
[99] Victor C. Motrescu and Ursula van Rienen. Computation of electrostatic fields
”
in anisotropic human tissues using the Finite Integration Technique (FIT)“. In:
Advances in Radio Science 2 (2004), pp. 309–313. doi: 10.5194/ars-2-309-2004.
[100] Takayoshi Nakata et al. Improvements of convergence characteristics of Newton-
”
Raphson method for nonlinear magnetic field analysis“. In: IEEE Transac-
tions on Magnetics 28.2 (Mar. 1992), pp. 1048 –1051. issn: 0018-9464. doi:
10.1109/20.123861.
[101] André Nicolet and François Delincé. Implicit Runge-Kutta Methods for Transient
”
Magnetic Field Computation“. In: IEEE Transactions on Magnetics 32.3 (May
1996), pp. 1405 –1408. doi: 0.1109/20.497510.
[102] Kurt Preis, H. Stogner, and Kurt R. Richter. Finite element analysis of anisotropic
”
and nonlinear magnetic circuits“. In: IEEE Transactions on Magnetics 17.6 (Nov.
1981), pp. 3396–3398. doi: 10.1109/TMAG.1981.1061469.
[103] Alfio Quarteroni. Numerical Models for Differential Problems. Modeling, Simula-
tion and Applications. Heidelberg: Springer, 2009. isbn: 978-88-470-1070-3. doi:
10.1007/978-88-470-1071-0.
[104] Alfio Quarteroni and Alberto Valli. Domain Decomposition Methods for Partial Dif-
ferential Equations. Numerical Mathematics and Scientific Computation. Oxford:
Oxford University Press, 1999.
[105] Alfio Quarteroni and Alberto Valli. Numerical Approximation of Partial Differential
Equations. Vol. 23. Springer Series in Computational Mathematics. Berlin: Springer,
2008.
[106] Muruhan Rathinam, Linda, and Linda R. Petzold. Dynamic Iteration Using Re-
”
duced Order Models: A Method For Simulation Of Large Scale Modular Systems“.
In: SIAM Journal on Numerical Analysis 40 (2002), pp. 1446–1474.
123
Bibliography
[107] Stefan Reitzinger and Joachim Schöberl. An algebraic multigrid method for fi-
”
nite element discretizations with edge elements“. In: Numerical Linear Algebra with
Applications 9.3 (2002), pp. 223 –238. doi: 10.1002/nla.271.
[108] Ricardo Riaza. Differential-Algebraic Systems: Analytical Aspects and Circuit Ap-
plications. World Scientifc, 2008.
[109] Janne Roos and Luis R. J. Costa, eds. Scientific Computing in Electrical Engineer-
ing SCEE 2008. Mathematics in Industry 14. Springer, 2010.
[110] Yousef Saad. Iterative methods for sparse linear systems. 2nd ed. Boston: SIAM,
2003. isbn: 978-0898715347.
[111] Sheppard J. Salon. Finite Element Analysis of Electrical Machines. Kluwer, 1995.
[112] Valeriu Savcenco, Willem Hundsdorfer, and Jan G. Verwer. A multirate time step-
ping strategy for parabolic PDE. Tech. rep. MAS-E0516. Amsterdam: Centrum voor
Wiskunde en Informatica, 2005.
[113] Sebastian Schöps. Coupling and Simulation of Lumped Electric Circuits Refined
”
by 3-D Magnetoquasistatic Conductor Models Using MNA and FIT“. MA thesis.
Wuppertal: Bergische Universität, 2008.
[114] Sebastian Schöps. Field Device Simulator User’s Manual. 2009-08-31. 2009.
[115] Sebastian Schöps, Andreas Bartel, and Markus Clemens. Efficient Solution of Tran-
”
sient Eddy Current Problems Using Linear/Nonlinear Domain Substructuring“. In
preparation. 2011.
[116] Sebastian Schöps, Andreas Bartel, and Herbert De Gersem. Multirate Time-
”
Integration of Field/Circuit Coupled Problems by Schur Complements“. In: Sci-
entific Computing in Electrical Engineering SCEE 2010. Ed. by Bas Michielsen.
Mathematics in Industry. To be published. Berlin: Springer, 2011. isbn: 978-3-540-
71979-3.
[117] Sebastian Schöps, Andreas Bartel, Herbert De Gersem, and Michael Günther.
DAE-Index and Convergence Analysis of Lumped Electric Circuits Refined by
”
3-D MQS Conductor Models“. In: Scientific Computing in Electrical Engineering
SCEE 2008. Ed. by Janne Roos and Luis R. J. Costa. Mathematics in Industry 14.
Springer, 2010, pp. 341 –350. doi: 10.1007/978-3-642-12294-1.
[118] Sebastian Schöps, Herbert De Gersem, and Andreas Bartel. A Cosimula-
”
tion Framework for Multirate Time-Integration of Field/Circuit Coupled Prob-
lems“. In: IEEE Transactions on Magnetics 46.8 (2010), pp. 3233 –3236. doi:
10.1109/TMAG.2010.2045156.
[119] Siegfried Selberherr. Analysis and Simulation of Semiconductor Devices. Springer,
1984.
[120] Monica Selva Soto and Caren Tischendorf. Numerical Analysis of DAEs from Cou-
”
pled Circuit and Semiconductor Simulation“. In: Applied Numerical Mathematics
53.2–4 (2005), pp. 471–488. doi: 10.1016/j.apnum.2004.08.009.
124
Bibliography
125
Bibliography
126
Index
127
Index
128
Index
regularization, 13
Coulomb’s gauge, 5, 13
grad-div, 13
reluctivity, 5, 11
nonlinear, 9, 11, 43, 44, 113
tensor, 9
Runge-Kutta methods, 40, 95
129