**1. Introduction**

An accurate numerical solution of Electromagnetic scattering problems is critically demanded in the simulation of many real-life applications, such as in the design of industrial processes and in the study of wave propagation phenomena. Electromagnetic (EM) scattering problems address the physical issue of computing the diffraction pattern of the EM radiation that is propagated by a complex body, illuminated by an incident wave. An explicit solution is possible only for simple targets, e.g. for spherical bodies; complicated geometries impose to use approximate numerical techniques.

Until the emergence of high-performance computing in the early eighties, the analysis of scattering problems was afforded by using approximate high frequency techniques such as the shooting and bouncing ray method (SBR) (Lee et al. (1988)). Ray-based asymptotic methods like SBR and the uniform theory of diffraction are based on the idea that EM scattering becomes a localized phenomenon as the size of the scatterer increases with respect to the wavelength. In the last decades, due to impressive advances in computer technology and the introduction of innovative algorithms with limited computational and memory requirement, a more rigorous numerical solution has become possible for many practical applications.

Finite-difference (FD) (Kunz & Luebbers (1993); Taflove (1995)), finite-element (FE) (Silvester & Ferrari (1990); Volakis et al. (1998)) and finite-volume (FV) methods (Bonnet et al. (1998); Botha (2006)) can be used to discretize the Maxwell's equations into a finite volume surrounding the scatterer, giving rise to sparse systems of linear equations. Upon inversion of the system, a solution is computed for all excitations. More recently, alternative approaches based on integral equations are becoming increasingly popular for solving high-frequency EM scattering problems. They reformulate the Maxwell's equations in the frequency domain and solve for the electric and the magnetic currents induced on the surface of the object. Thus integral methods require only a simple description of the surface of the target by means of triangular facets (see an example of discretization in Figure 1). This means that a 3D problem is reduced to solving a 2D surface problem, simplifying considerably the mesh generation especially in the case of moving objects. No artificial boundaries need to be imposed and boundary conditions are automatically satisfied in the case of perfectly conducting objects.

Another interest to use surface discretizations is that they noticeably reduce the effect of grid dispersion errors. Grid dispersion errors occur when a wave has a different phase velocity on

Fast Preconditioned Krylov Methods for Boundary Integral Equations

(2005)); we denote by *<sup>G</sup>*(|*<sup>y</sup>* <sup>−</sup> *<sup>x</sup>*|) = *<sup>e</sup>ik*|*y*−*x*<sup>|</sup>

for Boundary Integral Equations in Electromagnetic Scattering

 Γ (

*Rext j* is defined as

and is evaluated in the domain exterior to the object.

equation is applied to a set of tangential test functions *j*

*<sup>ϕ</sup>i*(*x*) · *<sup>ϕ</sup>j*(*x*) <sup>−</sup> <sup>1</sup>

for each 1 ≤ *i* ≤ *n*. The set of equations (2) can be recast in matrix form as

compute the set of coefficients {*λi*}1≤*i*≤*<sup>n</sup>* such that

*G* (|*y* − *x*|)

*Rext j* ∧*ν*).

*Rext <sup>j</sup>*(*y*) =

divergence operator defined as

equations in electromagnetism.

The operator

∑ 1≤*i*≤*n* *λi* Γ Γ

in Electromagnetic Scattering 3 Eqn. (1) is called Electric Field Integral Equation (EFIE) (see Bilotti & Vegni (2003); Li et al.

<sup>157</sup> Fast Preconditioned Krylov Methods

Γ is the boundary of the object, *k* the wave number and *Z*<sup>0</sup> = *μ*0/*ε*<sup>0</sup> the characteristic impedance of vacuum (*�* is the electric permittivity and *μ* the magnetic permeability). Given a continuously differentiable vector field *j*(*x*) represented in Cartesian coordinates on a 3D Euclidean space as *j*(*x*1, *x*2, *x*3) = *jx*<sup>1</sup> (*x*1, *x*2, *x*3)*ex*<sup>1</sup> + *jx*<sup>2</sup> (*x*1, *x*2, *x*3)*ex*<sup>2</sup> + *jx*<sup>3</sup> (*x*1, *x*2, *x*3)*ex*<sup>3</sup> , where *ex*<sup>1</sup> ,*ex*<sup>2</sup> ,*ex*<sup>3</sup> are the unit basis vectors of the Euclidean space, we denote by div*j*(*x*) the

div*<sup>j</sup>*(*x*) = *<sup>∂</sup>jx*<sup>1</sup>

disconnected parts, breaks on the surface; hence, it is very popular in industry.

*j <sup>t</sup>* + 1 2 Γ *j*. *j <sup>t</sup>* <sup>=</sup> <sup>−</sup> Γ

Γ

*∂x*<sup>1</sup> + *∂jx*<sup>2</sup> *∂x*<sup>2</sup> + *∂jx*<sup>3</sup> *∂x*<sup>3</sup> .

For closed targets, the Magnetic Field Integral Equation (MFIE) can be used, which reads

The EFIE formulation can be applied to arbitrary geometries such as those with cavities,

Both formulations suffer from interior resonances which make the numerical solution more problematic at some frequencies known as resonant frequencies, especially for large objects. The problem can be solved by combining linearly EFIE and MFIE. The resulting integral equation, known as Combined Field Integral Equation (CFIE), is the formulation of choice for closed targets. We point the reader to Gibson (2008) for a thorough presentation of integral

On discretizing Eqn. (1) in space by the Method of Moments (MoM) over a mesh containing *n* edges, the surface current *<sup>j</sup>* is expanded into a set of basis functions {*<sup>ϕ</sup>i*}1≤*i*≤*<sup>n</sup>* with compact support (the Rao-Wilton-Glisson basis, Rao et al. (1982), is a popular choice), then the integral

<sup>4</sup>*π*|*<sup>y</sup>* <sup>−</sup> *<sup>x</sup>*<sup>|</sup> the Green's function of Helmholtz equation,

(*<sup>H</sup> inc* <sup>∧</sup>*<sup>ν</sup>*).

*grad yG*(|*<sup>y</sup>* <sup>−</sup> *<sup>x</sup>*|) <sup>∧</sup>*<sup>j</sup>*(*x*)*dx*,

*t*

*<sup>k</sup>*<sup>2</sup> *div*<sup>Γ</sup>*<sup>ϕ</sup>i*(*x*) · *div*<sup>Γ</sup>*<sup>ϕ</sup>j*(*y*)

. Selecting *j*

<sup>=</sup> *<sup>i</sup> kZ*<sup>0</sup> Γ 

 *dxdy* =

*Aλ* = *b*, (3)

*<sup>t</sup>* = *ϕj*, we are led to

*Einc*(*x*) · *ϕj*(*x*)*dx*, (2)

*j t* .

Fig. 1. Example of surface discretization in an integral equation context. Each unknown of the problem is associated to an edge in the mesh. Courtesy of the EMC-CERFACS Group in Toulouse.

the grid compared to the exact solution; they tend to accumulate in space and may introduce spurious solutions over large 3D simulation regions (Bayliss et al. (1985); Jr. (1994); Lee & Cangellaris (1992)). For second-order accurate differential schemes, to alleviate this problem the grid density may grow up to <sup>O</sup>((*kd*)3) unknowns in 2D and of <sup>O</sup>((*kd*)4.5) in 3D, where *<sup>k</sup>* is the wavenumber and *d* is the approximate diameter of the simulation region. Therefore, the overall solution cost may increase considerably also for practical (*i.e.* finite) values of wavenumber (Chew et al. (1997)).

Boundary element discretizations are applied in many scientific and engineering areas beside electromagnetics and acoustics, e.g. in biomagnetic and bioelectric inverse modeling, magnetostatic and biomolecular problems, and many other applications (Forsman, Gropp, Kettunen & Levine (1995); Yokota, Bardhan, Knepley, Barba & Hamada (2011)). The potential drawback is that they lead, upon discretization, to large and dense linear systems to invert. Hence fast numerical linear algebra methods and efficient parallelization techniques are urged for solving large-scale boundary element equations efficiently on modern computers. In this chapter we overview some relevant techniques. In Section 2 we introduce the boundary integral formulation for EM scattering from perfectly conducting objects. In Section 4 we discuss fast iterative solution strategies based on preconditioned Krylov methods for solving the dense linear system arising from the discretization. In Section 5 we focus our attention on the design of the preconditioner, that is a critical component of Krylov methods in this context. We conclude our study in Section 5 with some final remarks.
