1

Parallel Matrix-Free Implementation of Frequency-Domain
Finite Difference Methods for Cluster Computing
Amir Geranmayeh

arXiv:1705.08849v1 [] 23 May 2017

Continental Automotive GmbH (I ID RD EE EL ED), VDO-Str. 1, 64832 Babenhausen, Germany
E-mail: amir.geranmayeh@continental-corporation.com

Abstract— Full-wave 3D electromagnetic simulations of complex planar devices, multilayer interconnects, and chip packages
are presented for wide-band frequency-domain analysis using the
finite difference integration technique developed in the PETSc
software package. Initial reordering of the index assignment
to the unknowns makes the resulting system matrix diagonally
dominant. The rearrangement also facilitates the decomposition
of large domain into slices for passing the mesh information
to different machines. Matrix-free methods are then exploited
to minimize the number of element-wise multiplications and
memory requirements in the construction of the system of
linear equations. Besides, the recipes provide extreme ease of
modifications in the kernel of the code. The applicability of
different Krylov subspace solvers is investigated. The accuracy is checked through comparisons with CST MICROWAVE
STUDIOr transient solver results. The parallel execution of the
compiled code on specific number of processors in multi-core
distributed-memory architectures demonstrate high scalability of
the computational algorithm.
Keywords— finite difference frequency domain (FDFD), finite
integration technique, matrix-free method, parallel programming, computational electromagnetics.

I. I NTRODUCTION
The integrated circuits (IC) industry has previously practiced the extraction of equivalent lumped-element models as
a simulation tool for hardware design. Advances in semiconductor interconnect technology, such as ever increasing
clock frequency and package density, however, demand fullwave analysis of a product-level model to preserve the signal
integrity in high-speed electronic devices [1], [2]. As a general
multi-grid difference-based electromagnetic field solver tool
[3], [4], the finite integration technique (FIT) provides a direct
discretization of the Maxwell’s equations in their integral form
using voltages and fluxes, respectively, along the edges and
surfaces of a pair of interlaced structured grids [5], [6].
To perform matrix-vector products for the matrices that
do not fit in the cache, CPU requires more time to calling
the matrix entries from the main memory than performing
the floating-point multiplications and additions. Matrix-free
methods, instead of saving the matrix in memory, recompute
the operator representing matrix elements with additional
floating-point operations [7]. On the other hand, the matrixfree methods are not able to use black-box preconditioners.
Modern object-oriented programming paradigms supply
an enormous level of abstraction in mechanisms needed
within parallel computing [8]. The software package PETSc
(Portable, Extensible Toolkit for Scientific Computation) consists of an expanding set of data structures and routines such

as parallel matrix and vector assembly which provide building
blocks for facilitating the development of scientific application
codes on parallel (and serial) computers [9]. PETSc uses
the MPI standard for message-passing communications and
it provides a collection of Krylov subspace (KSP) methods
which can be employed as parallel linear solvers for the
frequency-domain FIT matrix systems. This paper works out
how the PETSc shell matrix concept can create a foundation
for building and solving large-scale FIT systems on parallel
processors. It is shown that the computational speed and
memory demands of 3D finite-difference frequency-domain
(FDFD)-based electromagnetic solvers can be scaled down simultaneously using the advanced parallel programming tools.
The remainder of this report is organized as follows. In
Section II, the FIT discretization scheme is briefly explained.
Section III includes implementation aspects of the matrix-freebased electromagnetic-field solver. In Section IV, the accuracy
of the developed code is checked through eigenmode calculations in a rectangular cavity as well as comparisons with
CST MICROWAVE STUDIOr (MWS) discrete-port results
[10]. The overall performance of the parallelized algorithm
is measured on cluster of computers. Numerical simulations
of complex planar devices, multilayer interconnects, and chip
packages are presented.
II. F INITE I NTEGRATION T ECHNIQUE
The FIT renders a consistent transformation of the integral
form of Maxwell’s equations into a set of matrix equations
with a unique solution [5]. Maxwell’s equations are hereby
discretized on a pair of spatially staggered grids, namely the
e with edge lengths (Ln , L
e n ) and
primal and dual grid (G, G)
e
facet surface areas (An , An ), respectively, where the dual cell
e are defined on the barycenter of the primary cell
gridpoints G
mesh G [11]. The state variables of the FIT are the electric
and magnetic grid voltages
Z
Z
a
a
en =
E(r, t) · dl, h n =
H(r, t) · dl,
Ln

en
L

the electric and magnetic grid fluxes
Z
Z
a
a
a
a
dn =
D(r, t) · dA, b n =
en
A

B(r, t) · dA,

An

passing through surfaces of the pair of interlaced structured
grids, as well as the current
Z
a
a
jn =
J(r, t) · dA.
en
A

2

a

a
a

a
a

a
Storing
the discrete components in vectors e
, h, d, b, and
a
a
j , the Faraday’s and Ampère’s laws can be transformed into
a set of grid equations

a
Ce
=−

a
a
a
a
d a
eh
C
=
d+ j
dt

a
d a
b,
dt

(1)

e are, respectively, the discrete
in which the matrices C and C
curl-operators on the primal and dual grid. Owing to the
e = CT where T denotes the transpose of
topological nature, C
matrix. In Cartesian grid systems with Np = Nx × Ny × Nz
primal grid nodes


0 − Pz Py
C =  Pz 0 − Px 
(2)
−Py Px 0
where [Pw ]Np ×Np is the discrete partial differentiation operator along w : x, y, z directions and it contains only two nonzero entries per line, +1 and −1 . The voltages and fluxes are
coupled by the diagonal material matrices
a
a

a
a

a

h = M−1
µ b.

a
d = Mε e
,

(3)

The material matrices hold the locally averaged (inhomogeneous) permittivity and permeability distribution as well as
the metrics of the non-equidistant grid. The permittivity matrix
en and
elements are obtained by averaging on the dual facet A
the permeability matrix elements are obtained by averaging
e n . They may consist of complex values
along the dual edge L
in frequency domain to represent the dielectric or magnetic
losses [12]. Enforcing the Laplace transform to the timedomain representation of Maxwell’s equations in (1), all the
variables are transferred to the spectral domain. Replacing the
d
with the Laplace variable s = jω lets us
time-derivative dt
to analyse each frequency f = ω/2π in the time-harmonic

x, y, z unknown components are indexed one after the other
cell by cell. In other words, the unknown vector components
are initially arranged starting from the one with the lowest cell
indices in x-, y-, z-directions and proceed to the higher cell
indices sequentially. This corresponds to rearrange the vector
components of the classical FIT in the following manner




a
a
e 1x
e 1x
 a 
 a 
 e 1y 
 e 2x 




 a 

 .
 ..
 e 1z 





 a 
 a 
 eN x 
 e 2x 





 a
 a 
 e 2y 
 e 1y 
 a 


 e 2z 

 a


 e 2y 
(6)
.
 −→  ..
 .

 .

 ..





 ..
 a 
 eN y 

 .





 .
 a 
 e 1z 

 ..




 a 
 a 
 e 2z 
 eN x 



 .
 a 

 .
 eN y 

 .
a

eN z

a

eN z

This standard cache-based architecture is used for the layout
of unknowns in PETSc. Of course, to avoid any matrix
permutation, the above index assignment are followed from
the early stage of the discretization, Fig. 1.

a

regime individualy. Eliminating for example h in (1), the soa
called wave equation w.r.t. e
is obtained
a
a

a
a
2
j
CT M−1
µ C e + s Mε e = −s

(4)

with ne ≈ 3Np unknowns. The curl-curl system in (4) can be
a0

symmetrized by the change of variable e =
−1/2
M−1/2
CT M−1
ε
µ CMε

|

{z
A

}

a0

2

a0

e +s e =

1/2
Mε

−1/2
−sMε

(a)

(b)

a

e,
a
a

j

(5)

−1/2

where Mε
represents the root of inverse permittivity values
and the matrix A has a block-banded structure with only
13 non-zero entries per line. Enforcing the electric boundary
condition (BC), the tangential electric field Etan as well as the
normal component of the magnetic flux density Bn vanish on
the perfectly electric conducting (PEC) walls.
To accelerate the solution of the final system of equations
by parallel processing, the resulting coefficient matrix should
have a band-diagonal structure to minimize communication
between multi-processors with shared or distributed memory.
In order to obtain a diagonally dominant matrix A, instead of
applying the discrete derivative operator to first x, then y, and
finally z components of all the grid nodes, every three in-cell

Fig. 1. Sparsity pattern of the discrete curl matrix C for a 15 ×
10 × 10 grid. (a) The classical FIT (2). (b) The reordered numbering
arrangement according to (6).

The complete-solution procedure briefly involves the following successive numerical stages: preparation of the (possibly complex) permittivity and permeability matrices, construction of the discrete curl matrix, incorporation of the
boundary conditions to the material matrices, setting up the
excitation vector, assemblage of the system matrix equation for
the frequency of interest, successive solution of the obtained
sparse system of equation, and at last post-processing of
the results for monitoring the field quantities on the desired
probes.
III. M ATRIX -F REE M ETHOD
The memory limitations for the storage of the system
matrix and the associated subspaces prohibit the usage of

3

adequately fine grids for an accurate modeling of electricallysmall geometrical complexities. As a remedy, the matrix-free
methods do not require explicit storage of the product of
matrices in (5). Replacing the Laplace operator s in (5) with
the angular frequency jω, the matrix-free algorithm requires
only the application of the linear matrix operator to a vector,
a0

a0

−1/2

−1/2
M−1/2
CT M−1
e −ω 2 I e = −jωMε
ε
µ C Mε
| {z }

a
a

j

x(1)

|
|

|

y(1)

{z

x(2)

|
|

{z

{z

y(2)

{z

}
}

{z

y(3)

MULTIPLICATIONS PER ITERATION AND MEMORY CONSUMPTIONS
FOR DIFFERENT TYPES OF

FIT IMPLEMENTATIONS WITH ne

DEGREES OF FREEDOM .

FIT
Algorithms
Symmetric
Memory Usage
Multiplications

Classical
(4)
(5)
×
X
164 ne 164 ne
13 ne
13 ne

Matrix-Free
(7)
(8)
X
X
80 ne 72 ne
12 ne 9 ne

(7)
utilize the above-mentioned two shell matrix operations

}

0

}

x(3)

TABLE I
C OMPARISON OF THE REQUIRED ELEMENT- WISE

}
(q)

th

where I is the identity matrix and .
represents the q
reuse of the memory buffers x and y. Thus, properly userdefined shell matrix operations can be applied instead of
explicit assemblage of A, using for example the PETSc
matrix operations MatMult() and MatMultTranspose()
between three VecPointwiseMult() vector operations plus
VecAXPY() for summing up the scalar −ω 2 term. Therefore,
the matrix-free implementation of the square matrix-vector
products in FIT requires only 2 temporary vectors, namely
x and y as (7) implies, and 4 temporary vectors when C is
non-square due to the possible elimination of null rows.
For preallocation of the sparse parallel matrix C via
MatCreateMPIAIJ(), the number of nonzeros per row is

set to 4. To store C with 4 non-zero per line when using
a CSR-format, (4 + 8) 4 ne bytes of memory is required for
the storage of the integer column indices and the nonzero real
values. Indeed, the CSR-format also requires a third array with
length 4ne (in fact ne + 1 entries with typically 4 bytes each)
pointing to the beginning of each row. This third array does
affect the memory consumption here as the average number of
nonzeros per row is relatively small. Additionally, two floatingpoint arrays, each allocating 8ne bytes, are needed to store
−1/2
the diagonal elements of the material matrices Mε
and
−1
Mµ in double-precision. Hence, the matrix-free matrices in
(7) demands 8 ne + (4 + 4 + 8) 4 ne + 8 ne = 80 ne bytes of
memory to store C and diagonals of the material matrices.
The explicit sparse storage of the FIT-system matrix A
with 13 double-precision entries per line in (5) requires
(4 + 8) 13 ne = 156 ne bytes plus the additional 8 ne bytes
−1/2
for saving Mε
needed for the final voltage scaling. Note
that the classical FIT (4) or (5) entail the multiplication of 13
non-zero elements in ne rows of the system matrix for every
iteration, whereas in (7) only (1 + 4 + 1 + 4 + 1) entries have
to be multiplied per row. Note that the summation of −ω 2
term should be added to the operation counts of the matrix
free methods.
On the other hand, considering A = ATA where A =
−1/2
−1/2
−1/2
Mµ CMε
and Mµ
represents the root of inverse
permeability values, one can construct A and alternatively

0

a
a
(M−1/2
CM−1/2
)T (M−1/2
CM−1/2
)e
−ω 2 I e
.
µ
ε
µ
ε
|
{z
}
t
{z
}
|
x
|
{z
}

(8)

y

This yields a more efficient algorithm than (7) due to the
less vector operations and overwriting A on C by applying
MatDiagonalScale() and deallocating the memory reserved
−1/2
shrinks the memory usage to (4 + 4 + 8) 4 ne +
for Mµ
8 ne = 72 ne bytes. This calls for (4 + 4) ne multiplications
per iteration and 2 ne multiplications for the construction of
A in advance.
The KSP solvers in PETSc support matrix-free methods.
Matrix-free methods efficiently carry out matrix-vector products but the associated iterative solvers do not work well without a preconditioner [13]. Hence, a user-defined preconditioner
should be devised [14]. An additional matrix-vector product,
however, is not necessarily needed to extract for example
the diagonal elements of the shell matrix, as a diagonal
preconditioner can be constructed by taking the sum square
elements of the columns of the sparse matrix A, i.e.,
X
X
X
2
−1
[P]ii =
[M−1
[A]T
[A]2ki .
µ ]kk [C]ki [Mε ]ii =
ik [A]kj =
k

k,j=i

k

Table I reflects the computational costs and storage demands
for the introduced formulations.
When an extremely large mesh has to be generated, in the
pre-processing stage the simulation domain can simply be decomposed into slices along the z-axis to accumulate the information associated with slices of material matrices on different
machines. The PETSc provides distributed arrays which can be
decomposed along all three axes. The slicing along the z-axis
is an intuitive option which facilitates the implementation and
might not be the optimal domain decomposition. Other options
for mesh partitioning are using the METIS or Scotch software
packages. The optimal domain decomposition for the Jacoby
method is the one which minimizes the surface-to-volume ratio
of each submesh owned by an MPI rank [15].
IV. N UMERICAL R ESULTS AND D ISCUSSION
The correctness of the code was first checked through
eigenmode calculations of rectangular and circular cavities for
which analytical solution exists.
In the following examples,
a
a
time-harmonic currents i = j are imposed at the input ports

4

−1/2

0

a
as the excitation signal and the voltages u = Mε
e
along
the segments of desired probes are summed up as the output
values. The relative permeability µr =1 and zero-conductivity
σ = 0 S/m are assumed for the material properties, unless
otherwise mentioned. A relative decrease to less than 10−12 in
the residual norm is set as the stop criterion for iterative solvers
and the internal impedance of the source is assumed Z0 =
50 Ω for s-parameter calculations. The obtained impedance
(Z)-parameters of multi-port network devices are compared
with CST MICROWAVE STUDIOr (MWS) transient solver
results with discrete face-port excitation. The CST MWS
benchmark results are obtained on the same mesh grids used
for the presented frequency-domain methods. The presented
schemes also cover the usage of higher-order FIT for faster
convergence, since only the material matrices would then be
affected.

Fig. 2. A microstrip line enclosed in a metallic box. For the clarity of
representation, the surrounding box was represented as a transparent
one.

As the first example, a microstrip line enclosed in a PEC
box as depicted in Fig. 2 is analyzed. Fig. 3(a) illustrates
the execution times of different blocks of the program in
modeling the microstrip line at a single frequency f = 0.5
GHz on 96 processors in a distributed memory architecture
using the conjugate gradient method. The x-axis tick labels
represent different stages of the solution process, i.e. respectively, preparation of the discrete curl, permittivity, and
permeability matrices, inclusion of boundary conditions in
the material matrices, stimulation of excitation ports and
assemblage of the frequency-dependant shell system matrix,
iterative solution of the obtained system of equations using
Krylov subspace methods, determination of observation probes
to store the desired part of the solution, and finally the postprocessing impedance calculations. As displayed in Fig. 3(a),
the construction of the shell matrix does not take much time
in comparison with the KSP solver iterations till convergence.
To study the solver speedup under uniform memory access
time on multi-core shared memory processsors, a quad-core
64-bit Intel Xeon X5647 CPU with 12M cache and 24 GB
of RAM is considered. The speedup factor reaches to 3 for
solving the microstrip line problem with hundred thousand
unknowns on 4 cores, that is 75% efficiency, and it falls
down to the same efficiency factor 75% on 2 cores due to
the memory bandwidth limitations for larger problems with

million unknowns.
5

3

-1

0

(a)

(b)

Fig. 3. (a) The timing of different stages of constructing (7) for the
problem sketched in Fig. 2 discretized by 849×608×327 mesh cells
and solving it for ne = 506 384 352 unknowns on 96 nodes using the
conjugate gradient method. (b) The speedup diagram of solving the
matrix equation resulting from (7) for the microstrip line problem
depicted in Fig. 2 discretized by 497 × 354 × 193 mesh cells on
clusters of the workstation using the KSP conjugate gradient method.
The 8 nodes timing is used as reference.

Fig. 3(b) demonstrates the speedup factor in solving the
obtained linear system of equations using the KSP conjugate
gradient method. Fig. 4 compares the elapsed CPU time
for solving (7) on 64 nodes for ne =10 212 276 unknowns
associated with the problem shown in Fig. 2 by different
iterative methods provided in PETSc 3.4 [9]. Although the
conjugate gradient squared (cgs) method apparently converges
faster than the others, its matrix-free implementation does not
deliver accurate results close to the resonant frequencies of the
large-size problems, and hence, it is excluded from the rest
of paper. The generalized minimal residual method (gmres)
family have shown to be relatively slow. It seems that the
explicitly stored Krylov space generated by Gram-Schmidt
orthogonalisation constitutes the overhead for gmres. Note that
the flexible gmres (fgmres) is equivalent to gmres as long as
the preconditionner is not modified during the iterations. It is
also worth noting that the conjugate gradient (cg) method is
only applicable if the system matrix A is symmetric positive
definite (SPD), i.e. if A defined before (8) is invertible.
Fig. 5 demonstrates the balance of workload among the
used 96 processors for the solution of ne =506 384 352 FIT
unknowns using the biconjugate gradient stabilized method.
Less than 0.1 second imbalance is observed in the total run
time of different processors due to I/O and less than 1.6 MB in
the peak memory consumption among the ranks. Note that the
execution time and peak memory usage for the construction
of material matrices on the first rank has also been included in
Fig. 5. Fig. 6(a) displays the strong scaling (i.e. fixed global
problem size) experiment, considering the complete solution
time when the problem size is kept fixed ne =101 868 102 and
solved by the biconjugate gradient stabilized method. Fig. 6(b)
exhibits the weak scaling (i.e. fixed problem size per node)
experiment, cosidering the complete solution procedure with
the biconjugate gradient stabilized method when the number
of unknowns are proportionally scaled with the number of
cluster nodes. The iteration counts for different system sizes
grow almost linearly.
The correctness of the parallelized matrix-free implementation for the four short-listed robust iterative solvers from

5

(a)

(b)

Fig. 6.
(a) The strong scaling experiment with the simulation
of the microstrip line at f = 0.5 GHz for the fixed problem
size of ne =101 868 102 on a distributed-memory machine using
the biconjugate gradient stabilized method. (b) The weak scaling
experiment on the simulation of the microstrip line with the help
of biconjugate gradient stabilized method for a constant ratio of
workload to the used nodes in the cluster.
Fig. 4. The solution time for different PETSc iterative methods:
the biconjugate gradient stabilized (bcgs), conjugate gradient (cg),
conjugate gradient squared (cgs), conjugate residual (cr), generalized
minimal residual (gmres), flexible gmres (fgmres), two variants of
transpose-free quasi-minimal residual (tfqmr and tcqmr) methods.

(a)

(b)

Fig. 7. (a) A four-port coplanar coupler without ground plane.
The exact geometrical details can be found in CST-MWS examples
[10]. (b) The input impedance Z11 and the forward gain Z21 for the
coplanar coupler shown in Fig. 7(a).

(a)

(b)

Fig. 5. The load-balance in simulation of the microstrip line at
f = 0.5 GHz discretized by 849 × 608 × 327 mesh cells resulting
in ne = 506 384 352 unknowns distributed on 96 nodes solved by the
biconjugate gradient stabilized method.

Fig. 4, namely the bcgs, cr, cg, and gmres methods, are
now investigated on four complex case studies. As another
example of planar devices, a four-port coupler consisting of
coplanar lines without ground plane planted on a dielectric
substrate with the relative permittivity r =2.7 is considered.
The structure is placed up on a 500 µm thick aluminium layer
as shown in Fig. 7(a). The thickness of the thin substrate is
13 µm. Ports are placed vertically between the striplines and
the underlay ground to let the two (even and odd) quasi-TEM
modes be able to propagate along the lines. Perfect magnetic
conductors (PMC) are considered on all sides of the structure
to avoid the excitation of auxiliary box modes. The structured
mesh Nx =46, Ny =73, Nz =16 results in ne =161 184 degrees
of freedom for which (7) is solved at Nf =200 equidistant
frequencies within the range [fmin , fmax ] = [0.26, 26] GHz
using the generalized minimal residual method. As shown in
Fig. 7(b), the matrix-free method (MFM)-FIT results agree
well with the CST-MWS time-domain solver results on the
same hexahedral mesh. Note that, the present implementation
is equipped with the perfect boundary approximation (PBA)r
technique for accurate modeling of curved objects.
The next structure to be analyzed contains two vias for
transmitting the signal from one side of the printed circuit
board (PCB) to the other side while bridging over a metal
barrier, as depicted in Fig. 8(a). The thickness of the substrate

6

is 1.25 mm. The PEC boundary condition is considered for
all sides of the structure, except the top and bottom sides
which are terminated by PMC. The hexahedral mesh Nx = 31,
Ny = 77, Nz = 22 gives ne =157 542 degrees of freedom
in (8). Fig. 8(b) shows the associated impedances for the
frequency range [fmin , fmax ] = [1, 7.5] GHz obtained using
the conjugate gradient method when the ports are placed 2.5
mm distant from the PEC walls.

(a)

(b)

Fig. 8. (a) A transmission line using strip lines in three dimensions
placed on a dielectric substrate. The strip lines on the top surface of
the substrate are connected to a short strip line on the bottom surface
by small cylinders enabling the bypass of a possible barrier. (b) The
input impedance Z11 and the forward gain Z21 for the microstrip
transmission line in Fig. 8(a).

The next case study is composed of several different materials. An RJ45 connector, consisting of a socket jack and
4 differential pairs of wires for the signal transmission, is
analyzed. The wires of the male socket are fixed to a substrate
plate, the other wires are connected to a metallic ground plane
for shielding purposes, Fig. 9(a). The thickness of the metal
connectors are 0.5 mm, the substrate thickness is 1 mm, and
the insulator grip thickness is 1.5 mm. The simulation domain
is confined to a 34 × 14 × 11.75 mm3 metallic box. The
generated mesh Nx =63, Ny =45, Nz =28 yields to ne =238 140
degrees of freedom in (8) which are solved using the KSP
conjugate residual method. The discrete port 1 imposes a unit
current as the input at the plug end of the connector and the
crosstalk voltage along the coupled port 2 at the same side of
the connector is shown in Fig. 9(b).

(a)

(b)

Fig. 9. (a) An RJ45 connector containing 4 pairs of wire conductors
and the plastic socket. The model dimensions can be found in CST
STUDIO SUITE TM examples. (b) The input impedance Z11 and the
crosstalk impedance Z21 for the connector RJ45 in Fig. 9(a).

As the final example, a quad flat IC package consisting of a
silicon chip with the relative permittivity r =12.3 encapsulated
within a nonconductive plastic compound is simulated. The

surface mounted package is concentrically soldered to the PCB
with a PEC ground layer, Fig. 10(a). The thickneParallelss
of the metal strips is 0.1 mm, the substrate thickness is 0.5
mm, and the plastic thickness is 2.1 mm. The electric BC is
considered for all sides of the 13.4 × 13.4 × 3.075968 mm3
surrounding box. The structured mesh Nx =65, Ny =65, Nz =14
results in ne =177 450 degrees of freedom. The discrete port
1 imposes a specific current as the input and (7) is solved
for Nf =200 samples in the frequency range [fmin , fmax ] =
[0.06, 6] GHz using the biconjugate gradient stabilized method.
The induced voltage along the port 1 and the coupled voltage
along the neighboring pad are read as the outputs. The
broadband results shown in Fig. 10(b) comply with the CSTMWS transient solver results on the same hexahedral mesh.
The enormous abstraction of the code in PETSc environment
eases the solver modification, for example when mutual field
coupling between many ports in Fig. 10(a) are to be calculated.
In fact, one only needs to update the excitation vector on the
right hand side of (4) or (5), likewise (7) or (8), by activating
all pins of the chip at single run and the solver has not to be
run for each single port separately.

(a)

(b)

Fig. 10.
(a) A 52-pin IC package consisting of a silicon chip
encapsulated within a nonconductive plastic compound. The model
size can be found in CST-MWS examples [10]. (b) The input
impedance Z11 and the coupling impedance Z21 for the IC package
in Fig. 10(a).

V. C ONCLUSIONS
The matrix-free methods were adopted to the FIT scheme
to minimize the number of multiplications and memory requirements in the construction of the final system of linear
equations. The matrix-free FIT paves the path for realistic
modeling of electrically large-scale radiation problems. The
introduced recipe provides extreme ease of modifications in
the kernel of the applied algorithm, when material tensors
are diagonal which is the case in most practical applications with isotropic materials. Initial index reordering of
the unknowns was applied to make the FIT system matrix
diagonally dominant. The rearrangement also facilitated the
decomposition of large domain into slices for passing the mesh
information to different machines in the pre-processing stage.
Using the PETSc framework for high-performance computing, the accuracy and efficiency of different KSP solvers on

7

the shell matrix were investigated. The biconjugate gradient
stabilized and the conjugate residual methods are shown to
be optimal with respect to the solution time among the other
alternative iterative solvers for the presented problems with
different sizes. Performing large-scale simulations with 1.5
billion unknowns on 64 processors takes about several hours
per frequency sample with 6 GB peak memory usage. The
method also permitted running 102 millions of unknowns on
a shared-memory multiprocessor system in less than half a
day.
ACKNOWLEDGEMENT
This work was initiated under the German Federal Ministry of Education and Research (BMBF) funds to the
MoreSim4Nano partner in TU Darmstadt with contract
03MS613G.
R EFERENCES
[1] Y. Shao, Z. Peng, and J. F. Lee, “Full-wave real-life 3-D package signal
integrity analysis using nonconformal domain decomposition method,”
IEEE Trans. Microw. Theory Tech., vol. 59, no. 2, pp. 230–241, Feb.
2011.
[2] W. Yu, X. Yang, Y. Liu, and R. Mittra et al., “New development of
parallel conformal FDTD method in computational electromagnetics
engineering,” IEEE Antennas Propag. Mag., vol. 53, no. 3, pp. 15–41,
2011.
[3] C. S. Lavranos and G. A. Kyriacou, “Eigenvalue analysis of curved
waveguides employing an orthogonal curvilinear frequency-domain
finite-difference method,” IEEE Trans. Microw. Theory Tech., vol. 57,
no. 3, pp. 594–611, Mar. 2009.
[4] P. Przybyszewski, “Fast finite difference numerical techniques for the
time and frequency domain solution of electromagnetic problems,”
Supervised by M. Mrozowski, Technical University of Gdansk, 2001.
[5] R. Schuhmann and T. Weiland, “Conservation of discrete energy and
related laws in the finite integration technique,” Progress In Electromagnetics Research, vol. 32, pp. 301–316, 2001.
[6] E. Gjonaj, T. Weiland, I. Munteanu, and P. Thoma, “A parallel electromagentic simulation approach for the signal integrity analysis of IC
packages,” in Proc. IEEE Int. Electromagn. Compat. Symp. (EMC’07),
Honolulu, HI, 2007, pp. 1–5.
[7] T. Ren, T. Kalscheuer, S. Greenhalgh, and H. Maurer, “Boundary
element solutions for broadband 3D geo-electromagentic problems accelerated by multi-level fast mutlipole method,” Int. J. Geophys., vol.
192, no. 2, pp. 473–499, Feb. 2013.
[8] T. Iwashita, Y. Hirotani, T. Mifune, T. Murayama, and H. Ohtani, “Largescale time-harmonic electromagnetic field analysis using a multigrid
solver on a distributed memory parallel computer,” Parallel Computing,
vol. 38, no. 9, pp. 485–500, Sep. 2012.
[9] S. Balay, J. Brown, K. Buschelman, V. Eijkhout, W. Gropp, D. Kaushik,
M. Knepley, L. C. McInnes, B. Smith, and H. Zhang, PETSc Users
Manual Revision 3.4, Argonne, IL, May 2013.
[10] “CST STUDIO SUITE TM 2013,” CST - Computer Simulation Technology AG, Bad Nauheimer Str. 19, 64289 Darmstadt, Germany, Tech.
Rep., 2013.
[11] M. Clemens and T. Weiland, “Discrete electromagnetism with the finite
integration technique,” Progress In Electromagnetics Research, vol. 32,
pp. 65–87, 2001.
[12] R. C. Rumpf, C. R. Garcia, E. A. Berry, and J. H. Barton, “Finitedifference frequency-domain algorithm for modeling electromagnetic
scattering from general anisotropic objects,” Progress In Electromagnetics Research B, vol. 61, pp. 55–67, 2014.
[13] W. Shin, “3D finite-difference frequency-domain method for plasmonics
and nanophotonics,” Ph.D. dissertation, Stanford University, 2013.
[14] A. Chabory, B. de Hon, W. Schilders, and A. Tijhuis, “Preconditioned finitedifference frequency-domain for modelling periodic dielectric structures - comparisons with fdtd,” in Proc. 38th European
Microwave Conf. (EuMC’08), Amsterdam, the Netherlands, Oct 2008.
[15] G. Hager and G. Wellein, Introduction to High Performance Computing
for Scientists and Engineers, 1st ed. CRC Press, Taylor & Francis
Group, 2011.

