# P-FFT Accelerated EFIE with a Novel Diagonal Perturbed ILUT Preconditioner for Electromagnetic Scattering by Conducting Objects in Half Space

Authors: Lan-Wei Guo; Yongpin Chen; Jun Hu; Joshua Le-Wei Li
Journal: International Journal of Antennas and Propagation
 (2015)

---

## Abstract

A highly efficient and robust scheme is proposed for analyzing electromagnetic scattering from electrically large arbitrary shaped conductors in a half space. This scheme is based on the electric field integral equation (EFIE) with a half-space Green’s function. The precorrected fast Fourier transform (p-FFT) is first extended to a half space for general three-dimensional scattering problems. A novel enhanced dual threshold incomplete LU factorization (ILUT) is then constructed as an effective preconditioner to improve the convergence of the half-space EFIE. Inspired by the idea of the improved electric field integral operator (IEFIO), the geometrical-optics current/principle value term of the magnetic field integral equation is used as a physical perturbation to stabilize the traditional ILUT perconditioning matrix. The high accuracy of EFIE is maintained, yet good calculating efficiency comparable to the combined field integral equation (CFIE) can be achieved. Furthermore, this approach can be applied to arbitrary geometrical structures including open surfaces and requires no extra types of Sommerfeld integrals needed in the half-space CFIE. Numerical examples are presented to demonstrate the high performance of the proposed solver among several other approaches in typical half-space problems.


---

## Body

## 1. Introduction

The efficient analysis of electromagnetic scattering by electrically large and arbitrary shaped perfect electrically conducting (PEC) objects in a half space plays an important role in lots of applications. One of the most popular approaches is to formulate the surface integral equation (SIE) with the kernel of a half-space Green’s function to make the unknowns only relate to the boundary of the object [1]. The effect of the inhomogeneous background is therefore included in the dyadic half-space Green’s function, which is expressed in terms of the Sommerfeld integrals (SIs) [1–3]. The SIE is then converted into a matrix equation by the method of moments (MoM) and solved algebraically.Electric field integral equation (EFIE) is commonly applied due to its great accuracy and versatility [2]. However, the conditioning of the matrix that resulted from the discretization of EFIE is usually poor, especially for electrically large or densely meshed geometries. To speed up the convergence, the magnetic field integral equation (MFIE) is usually adopted to form a well-conditioned combined field integral equation (CFIE) in the free space, at the expense of losing some accuracy and/or versatility (not easily applicable to open or sharp surfaces for instance). In a half space, however, this expense could be even higher due to the complication of the half-space Green’s function; much more components of the computationally expensive Sommerfeld integrals are required in the magnetic-type integral operators [3, 4]. To overcome this difficulty, a well-conditioned improved EFIE (referred to as IEFIE) developed in [5] can be utilized and extended to the half-space problems. In this method, the geometrical-optics current of the MFIE operator is added into the EFIE operator in an iterative way and thus totally gets rid of the magnetic-type Green’s functions. However, the performance of the IEFIE is not always stable.On the other hand, SIE with global coupling leads to a dense impedance matrix, which requires huge memory and CPU resources, scaling asO(N2). For the solution of such a dense matrix system, an integral equation based fast algorithm is always required. A number of techniques for free space can be extended to a half-space or layered-medium background, after necessary manipulations or approximations [4, 6–8]. By splitting the dyadic Green’s function into terms that are in either convolution or correlation forms, we extend the p-FFT method [9] to model the electromagnetic scattering of general three-dimensional (3D) objects in a half space. Again, unlike the free space, the half-space CFIE with the p-FFT implementation requires more memory to store FFT propagation matrices, compared with its EFIE counterpart.Considering all these factors, EFIE is definitely more efficient to be implemented in the half space, if the conditioning issue can be addressed properly. There are many efforts to improve the conditioning issue by employing preconditioners. Simple preconditioners like the block-diagonal (BD) preconditioners [10] or symmetric successive overrelaxation (SSOR) preconditioner [11] are effective only when the matrix has some degree of diagonal dominance. However, they are ineffective for the EFIE matrix due to the weaker diagonal dominance and indefiniteness of EFIE [12]. By exploiting the self-regularizing property of the EFIE, the Caldern multiplicative preconditioner (CMP) gives rise to matrices that are well conditioned [2]. However, the BC basis functions involved usually make the construction cost of CMP higher than other candidates. Multiresolution (MR) preconditioner [13] with the hierarchic MR basis allows a simple diagonal preconditioner to improve the spectrum of the MoM matrix and also allows the use of different preconditioning schemes on different scales [12, 14]. Both incomplete LU (ILU) preconditioners [15–20] and sparse approximate inverse (SAI) preconditioners [21] have been successfully used in nonsymmetric dense systems. Even though ILU methods are inherently sequential while SAI preconditioners are easy to implement parallel computing, it is well known that SAI is not as successful as ILU with the same amount of memory [22].It is known that the dual threshold incomplete LU factorization (ILUT(τ,p)) scheme proposed by Saad [16] is quite powerful and robust in nonsymmetric dense systems. However, when the object is electrically large or densely meshed, the factors L and U are often ill-conditioned due to the large condition number of the original matrix system, resulting in the ILUT solver to be unstable. Several techniques can help to improve the quality of the preconditioner, such as reordering, scaling, diagonal shifting, pivoting, and condition estimators [17–20]. Unfortunately, none of them is stable or robust all the way, since no physical information is carried in these operations. For example, the shift is difficult to predict, and the pivoting may completely destroy the symmetry of the matrix system [18].In this paper, we will construct a robust diagonal perturbed incomplete ILUT (DP-ILUT) preconditioner based on a robust shifting formula for the efficient solution of EFIE with a half-space kernel. The geometrical-optics current/principle value term of the magnetic field integral equation is implemented as a physical perturbation for constructing a well-conditioned preconditioner matrix. This method makes the LU factors more diagonally dominant, which increases the robustness of LU decomposition. It should be noted that a similar idea was proposed in the ILU(0) for free-space problems [15]; however, we would like to emphasize that the method benefits more when extended to ILUT in a half-space problem. Numerical tests will be presented to validate our method.


## 2. Background: Integral Equations and Fast Iterative Solver

### 2.1. Integral Equations in Half-Space Problems

Consider a 3D PEC object of arbitrary shape located above a lossy half space as illustrated in Figure1. Define ϵk and μk as the permittivity and permeability of layer m(m=0,1), where possible loss is included. By imposing the PEC boundary condition for the electric field and magnetic field, the EFIE and MFIE can be derived as(1)EFIE:n^×n^×LEr,r′·Jr′=-n^×n^×Einc,MFIE:-0.5·I+n^×KHr,r′·Jr′=-n^×Hinc,where I denotes the unit tensor or identity matrix, which comes from the geometrical-optics current term of the MFIE operator, n^ stands for the outward pointing unit normal vector, r is the observation point, and r′ is the source point. The operators LE and KH are defined as [23](2)LEr,r′=-jωμ0∫dr′G¯AJr,r′+∇k02GVJr,r′,KHr,r′=P.V.∫dr′G¯HJr,r′,where P.V. stands for the Cauchy principal value integration and (3)G¯AJ=x^x^Gxxaj+x^z^Gxzaj+y^y^Gyyaj+y^z^Gyzaj+z^x^Gzxaj+z^y^Gzyaj+z^z^Gzzaj,G¯HJ=∇×G¯AJ,where G¯AJ, G¯HJ, and GVJ are the half-space dyadic and scalar Green’s functions [24] and the vectors r=ρρ^+zz^ and r′=ρ′ρ^+z′z^ denote field and source points in the cylindrical system. The CFIE is then obtained by linearly combining the EFIE and MFIE. In analogy with the CFIE, the IEFIE is constructed by adding the residue term of MFIE to EFIE as [5](4)IEFIE:EFIE+0.5αη0·I,where α is the combination coefficient varying between 0 and 1 and η0 is the wave impedance in free space.Figure 1
The geometry of a 3D PEC object above the half space.By discretizing the surface of the PEC object into triangular patches, the current density on the surface can be expanded by using the RWG basis functions [25]. The implementation of the standard MoM to the above equations renders a matrix system of the following form:(5)Z¯·I=Vinc.The linear algebraic equation of IEFIE is then given in the following iterative way:(6)Z¯IEFIE·Ii+1=Z¯EFIE+αη0Z¯IO·Ii+1≈VEFIE+αη0Z¯IO·Ii,where(7)Z¯mnIO=fmr,fnr′=∫fmr·fnr′dr′.The integral is nonzero only when the supports of the test function fm(r) and basis function fn(r′) are overlapped. Thus, the total element number of Z¯IO is no more than 5N, where N is the total number of RWG basis functions. Due to the introduction of the identity operator (IO) originated from MFIE, the eigenvalues of the IEFIE system are shifted away from the origin. Consequently, the condition number of IEFIE is much lower than that of EFIE, making it well conditioned [5].


### 2.2. The P-FFT Procedure for Half-Space Problems

In the half-space problems, each component of dyadic Green’s function (DGF) can be separated into the primary term (direct interaction in a homogeneous environment) and the secondary term (reflection and transmission due to the interface) [23, 24]. As a result, the corresponding dyadic half-space Green’s function can be termed as(8)G¯ρ,z,z′=G¯ρ,z-z′+G¯ρ,z+z′,where |ρ| is the horizontal distance between field and source points. The primary term, which is a function of (z-z′), is(9)G¯ρ,z-z′=I¯+∇∇k02e-jk0r-r′4πr-r′.The secondary term, a function of (z+z′), unfortunately, however, consists of many oscillatory and slowly convergent Sommerfeld integrals [1]. In our implementation, even though the singularity subtraction [26], the Weighted Averages Algorithm [27], and interpolation are used to accelerate the evaluation of each Sommerfeld integral, it is still computationally expensive.Similar to the classical p-FFT for structures in free space, the p-FFT for half space introduces an auxiliary regular grid that encloses the structure of interest. The matrix-vector multiplication (MVM)Z¯·I is then split into two parts: the near-field interactions (NFIs) Z¯FFI and far-field interactions (FFIs) Z¯NFI. The latter are accelerated by FFTs, while the former are a “precorrected" matrix which is calculated directly with precorrections [9]:(10)Z¯·I≈Z¯NFI·I+Z¯FFT-FFI·I,withZ¯FFT-FFI=Λ†G¯gridΛ,where Λ and Λ† represent the projection and interpolation operators, respectively. The FFT procedure of the primary term is the same as the classical one in free space, where G¯grid is a three-level block-Toeplitz structure matrix G¯T. For the secondary term, however, G¯grid (denoted as G¯H) has only Hankel-two-level block-Toeplitz structure because the half-space Green’s functions are not translationally invariant in the stratification direction [4]. Since each convolution/Hankel form of DGF could be converted into Toeplitz one by using a simple block antidiagonal permutation matrix P¯, G¯gridΛ can still be efficiently multiplied with vector I using 3D FFTs. Overall, we have (11)G¯gridΛI=IDFTDFTG¯T·DFTΛI+DFTG¯HP·DFTP¯-1ΛI,where DFT is the discrete Fourier transform and IDFT is the inverse transform. Therefore, the MVM of the whole system can be accelerated using the common four-stage procedure of p-FFT [9] in the iterative solvers.After simplifying DGF in half space by recombining the asymptotic terms, CFIE in the p-FFT implementation for half space requires almost two times memory request compared with EFIE. More importantly, CFIE also needs much more time-consuming Sommerfeld integrals for theG¯HJ evaluations in filling the impedance matrix. Considering all these factors, EFIE is more efficient to be implemented, if the conditioning issue can be addressed. In the next section, we will propose an efficient preconditioner to improve the convergence property of the EFIE in a half space.


## 2.1. Integral Equations in Half-Space Problems

Consider a 3D PEC object of arbitrary shape located above a lossy half space as illustrated in Figure1. Define ϵk and μk as the permittivity and permeability of layer m(m=0,1), where possible loss is included. By imposing the PEC boundary condition for the electric field and magnetic field, the EFIE and MFIE can be derived as(1)EFIE:n^×n^×LEr,r′·Jr′=-n^×n^×Einc,MFIE:-0.5·I+n^×KHr,r′·Jr′=-n^×Hinc,where I denotes the unit tensor or identity matrix, which comes from the geometrical-optics current term of the MFIE operator, n^ stands for the outward pointing unit normal vector, r is the observation point, and r′ is the source point. The operators LE and KH are defined as [23](2)LEr,r′=-jωμ0∫dr′G¯AJr,r′+∇k02GVJr,r′,KHr,r′=P.V.∫dr′G¯HJr,r′,where P.V. stands for the Cauchy principal value integration and (3)G¯AJ=x^x^Gxxaj+x^z^Gxzaj+y^y^Gyyaj+y^z^Gyzaj+z^x^Gzxaj+z^y^Gzyaj+z^z^Gzzaj,G¯HJ=∇×G¯AJ,where G¯AJ, G¯HJ, and GVJ are the half-space dyadic and scalar Green’s functions [24] and the vectors r=ρρ^+zz^ and r′=ρ′ρ^+z′z^ denote field and source points in the cylindrical system. The CFIE is then obtained by linearly combining the EFIE and MFIE. In analogy with the CFIE, the IEFIE is constructed by adding the residue term of MFIE to EFIE as [5](4)IEFIE:EFIE+0.5αη0·I,where α is the combination coefficient varying between 0 and 1 and η0 is the wave impedance in free space.Figure 1
The geometry of a 3D PEC object above the half space.By discretizing the surface of the PEC object into triangular patches, the current density on the surface can be expanded by using the RWG basis functions [25]. The implementation of the standard MoM to the above equations renders a matrix system of the following form:(5)Z¯·I=Vinc.The linear algebraic equation of IEFIE is then given in the following iterative way:(6)Z¯IEFIE·Ii+1=Z¯EFIE+αη0Z¯IO·Ii+1≈VEFIE+αη0Z¯IO·Ii,where(7)Z¯mnIO=fmr,fnr′=∫fmr·fnr′dr′.The integral is nonzero only when the supports of the test function fm(r) and basis function fn(r′) are overlapped. Thus, the total element number of Z¯IO is no more than 5N, where N is the total number of RWG basis functions. Due to the introduction of the identity operator (IO) originated from MFIE, the eigenvalues of the IEFIE system are shifted away from the origin. Consequently, the condition number of IEFIE is much lower than that of EFIE, making it well conditioned [5].


## 2.2. The P-FFT Procedure for Half-Space Problems

In the half-space problems, each component of dyadic Green’s function (DGF) can be separated into the primary term (direct interaction in a homogeneous environment) and the secondary term (reflection and transmission due to the interface) [23, 24]. As a result, the corresponding dyadic half-space Green’s function can be termed as(8)G¯ρ,z,z′=G¯ρ,z-z′+G¯ρ,z+z′,where |ρ| is the horizontal distance between field and source points. The primary term, which is a function of (z-z′), is(9)G¯ρ,z-z′=I¯+∇∇k02e-jk0r-r′4πr-r′.The secondary term, a function of (z+z′), unfortunately, however, consists of many oscillatory and slowly convergent Sommerfeld integrals [1]. In our implementation, even though the singularity subtraction [26], the Weighted Averages Algorithm [27], and interpolation are used to accelerate the evaluation of each Sommerfeld integral, it is still computationally expensive.Similar to the classical p-FFT for structures in free space, the p-FFT for half space introduces an auxiliary regular grid that encloses the structure of interest. The matrix-vector multiplication (MVM)Z¯·I is then split into two parts: the near-field interactions (NFIs) Z¯FFI and far-field interactions (FFIs) Z¯NFI. The latter are accelerated by FFTs, while the former are a “precorrected" matrix which is calculated directly with precorrections [9]:(10)Z¯·I≈Z¯NFI·I+Z¯FFT-FFI·I,withZ¯FFT-FFI=Λ†G¯gridΛ,where Λ and Λ† represent the projection and interpolation operators, respectively. The FFT procedure of the primary term is the same as the classical one in free space, where G¯grid is a three-level block-Toeplitz structure matrix G¯T. For the secondary term, however, G¯grid (denoted as G¯H) has only Hankel-two-level block-Toeplitz structure because the half-space Green’s functions are not translationally invariant in the stratification direction [4]. Since each convolution/Hankel form of DGF could be converted into Toeplitz one by using a simple block antidiagonal permutation matrix P¯, G¯gridΛ can still be efficiently multiplied with vector I using 3D FFTs. Overall, we have (11)G¯gridΛI=IDFTDFTG¯T·DFTΛI+DFTG¯HP·DFTP¯-1ΛI,where DFT is the discrete Fourier transform and IDFT is the inverse transform. Therefore, the MVM of the whole system can be accelerated using the common four-stage procedure of p-FFT [9] in the iterative solvers.After simplifying DGF in half space by recombining the asymptotic terms, CFIE in the p-FFT implementation for half space requires almost two times memory request compared with EFIE. More importantly, CFIE also needs much more time-consuming Sommerfeld integrals for theG¯HJ evaluations in filling the impedance matrix. Considering all these factors, EFIE is more efficient to be implemented, if the conditioning issue can be addressed. In the next section, we will propose an efficient preconditioner to improve the convergence property of the EFIE in a half space.


## 3. Diagonal Perturbed ILUT Preconditioner

In addition to being indefinite and non-Hermitian, some of the nonstored far-field interactions of EFIE matrices may be stronger than the NFIs [18]. This causes EFIE system much more difficult to solve than CFIE. For closed targets, the CFIE formulation typically converges in O(N0.25) iterations in Krylov subspace methods while EFIE formulation converges in O(N0.5) iterations [20]. Due to the large condition number of EFIE system, the iterative solver may even diverge in many cases. Therefore, a robust preconditioner should be constructed for EFIE.The ILU-class preconditioners are an effective preconditioning technique for solving integral equation systems. Defining the right preconditioner matrixM¯ as a nonsingular approximation of Z¯, the MVM can be written in the following form:(12)Z¯·M¯-1·u=V,I=M¯-1·u.In ILU-class preconditioners, the matrix M¯ is written in the form of M¯=L¯·U¯, where L¯ and U¯ are approximations of the L and U factors of the standard triangular LU decomposition of matrix Z¯ [28]. Among such ILU-class preconditioners, ILUT(τ,p) proposed by Saad is known to yield more accurate factorizations than the other level-of-fill methods with the same amount of fill-in [16]. Since Z¯FFT-FFI is not numerically available, the preconditioner for p-FFT is always constructed based on Z¯NFI. Hence, the preconditioning matrix is constructed as(13)M¯-1=ILUTZNFIEFIE.However, ILUT(τ,p) preconditioner may encounter difficulties in solving EFIE matrix. This is because small diagonal entries in EFIE matrix may increase the possibility of encountering zero pivots during the LU factorization, decreasing the stability and accuracy of the factorization results [19]. Lots of studies show that the small diagonal entries are responsible for increasing the instability of the LU factorization of matrix ZNFIEFIE. The diagonal shifting formula for ZNFIEFIE can be written as [19](14)M¯≈Z~NFIEFIE=Z¯NFIEFIE+C¯,where Z~NFIEFIE is the matrix after shifting and C¯ is the amount of shifts for the diagonal entries. Unfortunately, however, it is difficult to select suitable shift parameters C¯. A small shifting value may result in a seemingly more accurate but unstable LU factorization; on the other hand, large perturbations may introduce too much inaccuracy in the preconditioner [19].Based on the investigation of the IEFIE structure in (6) and the eigenvalue distributions shown in [5], it seems that Z¯IO in (7) could filter out those small eigenvalues in EFIE. Thus, Z¯IO could be naturally used as the diagonal entries and result in a Z~NFIEFIE matrix with a better condition number. The diagonal shifting matrix C¯=Z¯IO is diagonally dominant and very sparse, whose nonzero elements are near-field interactions. More importantly, the IO term is the residue term of MFIE, which corresponds to the geometrical-optics current and carries the information of basis functions. The magnitudes of the elements of Z¯IO are comparable with those in the EFIE matrix system Z¯EFIE. Hence, such preconditioner system is stable and physical system. By substituting (14) into (13), the diagonal perturbed ILUT preconditioner can be written as(15)M¯-1≈ILUTZ~NFIEFIE=ILUTZ¯NFIEFIE+αη0Z¯IO.This diagonal perturbation makes the system diagonally dominant and is more physical and stable than other shift parameters used in [18–20]. It therefore increases the robustness of ILUT(τ,p) for EFIE systems. A similar idea has also been proposed in [15] for free-space case; however, ILUT used in our paper is better than the ILU(0) preconditioner used in [15]. The reason is that ILUT contains more information of near-field interactions than ILU(0) [28]. In the next section, robustness and efficiency of the employed preconditioner for half-space applications are demonstrated over a bunch of numerical examples.


## 4. Numerical Results

In our experiments, starting with the zero initial guess, the restarted version of the generalized minimal residual (GMRES) is used as the iterative method. The CPU times reported in this section are obtained on a 64-bit server configured by Intel Xeon CPU with 3.47 GHz and 128 GB of memory. For ILUT, we set the drop toleranceτ as 10-4 and set p to be the average number of nonzero elements in a row of the near-field matrix.Six methods are studied and compared in the context of a half space and accelerated by the half-space p-FFT; they are (1) EFIE-NO, EFIE with no preconditioner; (2) EFIE-MILU0, EFIE with modified ILU0 proposed in [15]; (3) EFIE-ILUT, EFIE with conventional ILUT(τ,p); (4) EFIE-ILUTP5, EFIE with the common used ILUTP, where 0.5 is used as diagonal entries [16, 18]; (5) EFIE-DP-ILUT, EFIE with the proposed diagonal perturbations ILUT(τ,p); and (6) CFIE with no preconditioner as the reference. In all cases, the convergence threshold of GMRES is set at 0.001, and the combination factor α of the diagonal perturbations ZIO is set as 0.3.It is noted that a strong indicator of the instability quality of the ILUT preconditioner is an estimate of([L][U])-1, called “condest” [18]. This condition estimate is defined as(16)condest=LU-1·e∞,e=1,…,1Twhich can be easily computed by using a forward substitution followed by a backward substitution.


### 4.1. Eigenspectrum Distribution

This section analyzes the eigenspectrum distribution and the condition number of the preconditioned linear system for a simple case. A 9-inch almond with 2,883 unknowns is located 0.02 m above the half space with the bottom layer of relative permittivityεr=4. The plane wave is incident from θinc=80∘ and ηinc=90∘ at 3 GHz. In Figure 2, we compare the eigenspectrum distribution of the preconditioned linear system equation (12) obtained with different ILUT preconditioners. Table 1 shows the condition number which is defined as the ratio between the highest norm of the eigenvalue λmax and the lowest norm of the eigenvalue λmin; that is, σ=λmax/λmin. A spectrum clustered away from the origin always implies a rapid convergence of Krylov subspace methods. As shown in Figure 2, both ILUT and DP-ILUT (that we proposed in this paper) move the eigenvalues of the preconditioned system towards (1,0) and thus avoid a possible eigenvalue cluster close to zero. However, ILUTP5 has less effect in shifting the eigenvalues. It is also shown that the eigenvalues are more concentrated in the shifted-unit-circle for our DP-ILUT preconditioner. It should be emphasized that the reason for the failure of the ILUTP5 preconditioner is that the large perturbation 0.5 introduces too much inaccuracy in the preconditioner. This graphical observation is supported by the condition number of the linear system equation (12). Indeed, when DP-ILUT preconditioner is used, the condition number and the required number of matrix-vector products to reach convergence are smaller than those of other preconditioners. Table 2 lists the number of “condests” for different ILU preconditioners. It can be seen that the small diagonal entries sre responsible for increasing the stability of the LU factorization of matrix.Table 1
Condest, condition number, and iteration number of different preconditioned linear systems equation (12).
EFIE-NO
EFIE-ILUTP5
EFIE-ILUT
EFIE-DP-ILUT
Condest
—
2
109,938
4,105
Condition number
4,914
4,885
121
36
Iterative number
229
224
53
31Table 2
Condest, of different preconditioned linear systems equation (12) for 3D scattering example.
ILUT
ILUTP5
MILU0
DP-ILUT
Almond
1.14
E
+
11
2.5
693,098
1,956,402
Open cavity
133
3.4
3.5
2.7
Aircraft model
1050
3.5
21.3
23.6
Tank model
1467
5.6
125
163Figure 2
Eigenvalue distribution of the linear systems resulting from the preconditioned linear system equation (12) obtained with (a) EFIE-NO; (b) EFIE-ILUTP5; (c) EFIE-ILUT; and (d) EFIE-DP-ILUT. The red solid circle is the shifted-unit-circle.
(a)
(b)
(c)
(d)


### 4.2. Accuracy and Convergence Test

To demonstrate the accuracy and efficiency of the proposed scheme, the bistatic RCS of two benchmark targets is calculated. They consist of an almond with 43,743 unknowns at 9 GHz and an open elliptical cavity with 119,600 unknowns at 0.3 GHz. Both are situated above the half space with the bottom layer of relative permittivityεr=6.38-j0.663.A 9-inch almond is located 0.02 m above the interface. The plane wave is incident fromθinc=80∘ and ηinc=90∘ at 9 GHz. Figure 3 shows the bistatic RCS patterns of the almond by different methods. The results are further compared with the one calculated by FEKO. All of them agree well, except for CFIE, which is less accurate. This is because EFIE can produce numerical results with much higher accuracy than CFIE for irregular structures. Figure 4 shows the convergence history of the GMRES algorithm. It can be observed that EFIE-ILUT is counteracted in this example due to the ill-conditioned factors of the ILUT preconditioning matrix. Thanks to the diagonal perturbations, both EFIE-MILU0 and EFIE-DP-ILUT have better convergence than EFIE-NO. However, due to the use of 0.5 as the diagonal entries, ILUTP5 has no positive effect. It shows that the stability of LU factorization is an important factor for ILU-class preconditioners; at the same time, the accuracy of the diagonal shifting formula is also very important.Figure 3
The H-polarized bistatic RCS patterns of an almond shape above the half space in theθ=80∘ cut at 9 GHz.Figure 4
The convergence history of GMRES algorithms with different integral operators and preconditioners in the almond example.The open elliptical cavity, which is located 0.33λ above the interface, has a side length, major axis, and minor axis of 10 m, 10 m, and 0.5 m, respectively. The plane wave is incident from θinc=80∘, ϕinc=0∘ at 300 MHz. Results from different methods are investigated, where the FEKO simulation is also taken as the reference. As shown in Figure 5 and Table 3, the EFIE-DP-ILUT is as efficient as the conventional EFIE-ILUTP5 preconditioner.Table 3
Comparison of different integral operators and preconditioners for 3D scattering.
Unknowns
p (ILUT)
EFIE-DP-ILUT
EFIE-ILUT
EFIE-ILUTP5
CFIE
N
iter
T
sol
(sec)
N
iter
T
sol (sec)
N
iter
T
sol (sec)
N
iter
T
sol (sec)
Almond
43,743
200
53
2,034
800+
5000+
742
5,028
49
2,239
Open cavity
120,105
55
84
12,126
180
13,859
80
11,926
—
—
Aircraft model
69,534
100
276
8,752
800+
37,000+
800+
37,000+
134
8,372
Tank model
701,262
80
351
60,392
800+
130,000+
773
121,857
123
74,562Figure 5
The convergence history of GMRES algorithms with different integral operators and preconditioners in the open elliptical cavity example. Both the top view and the front view of the open elliptical cavity are illustrated.


### 4.3. Large-Scale Application

Next, the validation of the proposed scheme for large-scale computation is demonstrated by analyzing electromagnetic scattering from two complex 3D structures: aircraft and a tank model. Both are located above the half space with the bottom layer of relative permittivityεr=6.5-j0.6.The length, width, and height of the aircraft model are 12.5, 7, and 3 m, respectively, as shown in Figure6. It is located 0.3 m above the interface, and the plane wave is incident from θinc=60∘ and ϕinc=0∘. The frequency is 400 MHz and the corresponding electrical size is 16.7λ, resulting in 69,534 unknowns. Figure 6 shows the bistatic RCS in the vertical polarization. It is shown that the proposed EFIE-DP-ILUT agrees well with the result from CFIE. The convergence is shown in Figure 7, where all the other methods behave worse than the proposed ILUT-DP. This example also shows the instability of EFIE-ILUTP5 preconditioner.Figure 6
The V-polarized bistatic RCS patterns of the aircraft model above the half space in theθ=80∘ cut at 400 MHz with EFIE-DP-ILUT and CFIE.Figure 7
The convergence history of GMRES algorithms with different integral operators and preconditioners in the aircraft model example.The tank model is7×2×1.9m3. It is located above the interface, and the incidence plane wave in 1.5 GHz is fixed at θinc=60∘ and ϕinc=0∘. The surface of the model is discretized with ~λ0/10 average edge length resulting in 702,483 unknowns. The bistatic RCS in the vertical polarization is computed using proposed EFIE-DP-ILUT and compared with the one from conventional CFIE code. As is shown in Figure 8, EFIE-DP-ILUT agrees well with CFIE. Figure 9 shows that all the other methods EFIE-MILU0, EFIE-NO, EFIE-ILUT, and EFIE-ILUTP5 cannot converge to the threshold after 700 iterations.Figure 8
The V-polarized bistatic RCS patterns of the tank model above the half space in theθ=80∘ cut at 1.5 GHz with EFIE-DP-ILUT and CFIE.Figure 9
The convergence history of GMRES algorithms with different integral operators and preconditioners in the tank model example.To further demonstrate the advantage of the present method, Table3 summarizes the number of iterations and the solving CPU time with different integral operators and preconditioners for the aforementioned 3D structures. All the four examples show that the proposed DP-ILUT is definitely more robust than other preconditioners. Although CFIE may require fewer iterations, it needs more CPU time for the evaluation of GHJ in filling Green function interpolation tables and also for the MFIE part of the p-FFT four-stage procedure. As shown in Table 4, the total CPU time of EFIE-DP-ILUT is comparable with CFIE, but the memory requirement is less than CFIE. In the large-scale application, the traditional EFIE-ILUT, EFIE-ILUTP5, and EFIE-MILU0 proposed by [15] do not converge in 700 iterations, while EFIE-DP-ILUT converges rapidly. It is worth pointing out that even though EFIE-ILUTP5 yields a more smaller “condest” value compared to EFIE-DP-ILUT, the ILUTP5 preconditioner is not always stable. More specifically, when sharp surfaces are involved in the geometrical model, our DP-ILUT method is more stable than the common used ILUTP5.Table 4
Comparison of memory requirement of EFIE-DP-ILUT and CFIE for 3D scattering.
Almond
Open cavity
Aircraft
Tank
EFIE-DP-ILUT (MB)
593
760
1,434
15,100
CFIE (MB)
623
—
2,106
21,585


## 4.1. Eigenspectrum Distribution

This section analyzes the eigenspectrum distribution and the condition number of the preconditioned linear system for a simple case. A 9-inch almond with 2,883 unknowns is located 0.02 m above the half space with the bottom layer of relative permittivityεr=4. The plane wave is incident from θinc=80∘ and ηinc=90∘ at 3 GHz. In Figure 2, we compare the eigenspectrum distribution of the preconditioned linear system equation (12) obtained with different ILUT preconditioners. Table 1 shows the condition number which is defined as the ratio between the highest norm of the eigenvalue λmax and the lowest norm of the eigenvalue λmin; that is, σ=λmax/λmin. A spectrum clustered away from the origin always implies a rapid convergence of Krylov subspace methods. As shown in Figure 2, both ILUT and DP-ILUT (that we proposed in this paper) move the eigenvalues of the preconditioned system towards (1,0) and thus avoid a possible eigenvalue cluster close to zero. However, ILUTP5 has less effect in shifting the eigenvalues. It is also shown that the eigenvalues are more concentrated in the shifted-unit-circle for our DP-ILUT preconditioner. It should be emphasized that the reason for the failure of the ILUTP5 preconditioner is that the large perturbation 0.5 introduces too much inaccuracy in the preconditioner. This graphical observation is supported by the condition number of the linear system equation (12). Indeed, when DP-ILUT preconditioner is used, the condition number and the required number of matrix-vector products to reach convergence are smaller than those of other preconditioners. Table 2 lists the number of “condests” for different ILU preconditioners. It can be seen that the small diagonal entries sre responsible for increasing the stability of the LU factorization of matrix.Table 1
Condest, condition number, and iteration number of different preconditioned linear systems equation (12).
EFIE-NO
EFIE-ILUTP5
EFIE-ILUT
EFIE-DP-ILUT
Condest
—
2
109,938
4,105
Condition number
4,914
4,885
121
36
Iterative number
229
224
53
31Table 2
Condest, of different preconditioned linear systems equation (12) for 3D scattering example.
ILUT
ILUTP5
MILU0
DP-ILUT
Almond
1.14
E
+
11
2.5
693,098
1,956,402
Open cavity
133
3.4
3.5
2.7
Aircraft model
1050
3.5
21.3
23.6
Tank model
1467
5.6
125
163Figure 2
Eigenvalue distribution of the linear systems resulting from the preconditioned linear system equation (12) obtained with (a) EFIE-NO; (b) EFIE-ILUTP5; (c) EFIE-ILUT; and (d) EFIE-DP-ILUT. The red solid circle is the shifted-unit-circle.
(a)
(b)
(c)
(d)


## 4.2. Accuracy and Convergence Test

To demonstrate the accuracy and efficiency of the proposed scheme, the bistatic RCS of two benchmark targets is calculated. They consist of an almond with 43,743 unknowns at 9 GHz and an open elliptical cavity with 119,600 unknowns at 0.3 GHz. Both are situated above the half space with the bottom layer of relative permittivityεr=6.38-j0.663.A 9-inch almond is located 0.02 m above the interface. The plane wave is incident fromθinc=80∘ and ηinc=90∘ at 9 GHz. Figure 3 shows the bistatic RCS patterns of the almond by different methods. The results are further compared with the one calculated by FEKO. All of them agree well, except for CFIE, which is less accurate. This is because EFIE can produce numerical results with much higher accuracy than CFIE for irregular structures. Figure 4 shows the convergence history of the GMRES algorithm. It can be observed that EFIE-ILUT is counteracted in this example due to the ill-conditioned factors of the ILUT preconditioning matrix. Thanks to the diagonal perturbations, both EFIE-MILU0 and EFIE-DP-ILUT have better convergence than EFIE-NO. However, due to the use of 0.5 as the diagonal entries, ILUTP5 has no positive effect. It shows that the stability of LU factorization is an important factor for ILU-class preconditioners; at the same time, the accuracy of the diagonal shifting formula is also very important.Figure 3
The H-polarized bistatic RCS patterns of an almond shape above the half space in theθ=80∘ cut at 9 GHz.Figure 4
The convergence history of GMRES algorithms with different integral operators and preconditioners in the almond example.The open elliptical cavity, which is located 0.33λ above the interface, has a side length, major axis, and minor axis of 10 m, 10 m, and 0.5 m, respectively. The plane wave is incident from θinc=80∘, ϕinc=0∘ at 300 MHz. Results from different methods are investigated, where the FEKO simulation is also taken as the reference. As shown in Figure 5 and Table 3, the EFIE-DP-ILUT is as efficient as the conventional EFIE-ILUTP5 preconditioner.Table 3
Comparison of different integral operators and preconditioners for 3D scattering.
Unknowns
p (ILUT)
EFIE-DP-ILUT
EFIE-ILUT
EFIE-ILUTP5
CFIE
N
iter
T
sol
(sec)
N
iter
T
sol (sec)
N
iter
T
sol (sec)
N
iter
T
sol (sec)
Almond
43,743
200
53
2,034
800+
5000+
742
5,028
49
2,239
Open cavity
120,105
55
84
12,126
180
13,859
80
11,926
—
—
Aircraft model
69,534
100
276
8,752
800+
37,000+
800+
37,000+
134
8,372
Tank model
701,262
80
351
60,392
800+
130,000+
773
121,857
123
74,562Figure 5
The convergence history of GMRES algorithms with different integral operators and preconditioners in the open elliptical cavity example. Both the top view and the front view of the open elliptical cavity are illustrated.


## 4.3. Large-Scale Application

Next, the validation of the proposed scheme for large-scale computation is demonstrated by analyzing electromagnetic scattering from two complex 3D structures: aircraft and a tank model. Both are located above the half space with the bottom layer of relative permittivityεr=6.5-j0.6.The length, width, and height of the aircraft model are 12.5, 7, and 3 m, respectively, as shown in Figure6. It is located 0.3 m above the interface, and the plane wave is incident from θinc=60∘ and ϕinc=0∘. The frequency is 400 MHz and the corresponding electrical size is 16.7λ, resulting in 69,534 unknowns. Figure 6 shows the bistatic RCS in the vertical polarization. It is shown that the proposed EFIE-DP-ILUT agrees well with the result from CFIE. The convergence is shown in Figure 7, where all the other methods behave worse than the proposed ILUT-DP. This example also shows the instability of EFIE-ILUTP5 preconditioner.Figure 6
The V-polarized bistatic RCS patterns of the aircraft model above the half space in theθ=80∘ cut at 400 MHz with EFIE-DP-ILUT and CFIE.Figure 7
The convergence history of GMRES algorithms with different integral operators and preconditioners in the aircraft model example.The tank model is7×2×1.9m3. It is located above the interface, and the incidence plane wave in 1.5 GHz is fixed at θinc=60∘ and ϕinc=0∘. The surface of the model is discretized with ~λ0/10 average edge length resulting in 702,483 unknowns. The bistatic RCS in the vertical polarization is computed using proposed EFIE-DP-ILUT and compared with the one from conventional CFIE code. As is shown in Figure 8, EFIE-DP-ILUT agrees well with CFIE. Figure 9 shows that all the other methods EFIE-MILU0, EFIE-NO, EFIE-ILUT, and EFIE-ILUTP5 cannot converge to the threshold after 700 iterations.Figure 8
The V-polarized bistatic RCS patterns of the tank model above the half space in theθ=80∘ cut at 1.5 GHz with EFIE-DP-ILUT and CFIE.Figure 9
The convergence history of GMRES algorithms with different integral operators and preconditioners in the tank model example.To further demonstrate the advantage of the present method, Table3 summarizes the number of iterations and the solving CPU time with different integral operators and preconditioners for the aforementioned 3D structures. All the four examples show that the proposed DP-ILUT is definitely more robust than other preconditioners. Although CFIE may require fewer iterations, it needs more CPU time for the evaluation of GHJ in filling Green function interpolation tables and also for the MFIE part of the p-FFT four-stage procedure. As shown in Table 4, the total CPU time of EFIE-DP-ILUT is comparable with CFIE, but the memory requirement is less than CFIE. In the large-scale application, the traditional EFIE-ILUT, EFIE-ILUTP5, and EFIE-MILU0 proposed by [15] do not converge in 700 iterations, while EFIE-DP-ILUT converges rapidly. It is worth pointing out that even though EFIE-ILUTP5 yields a more smaller “condest” value compared to EFIE-DP-ILUT, the ILUTP5 preconditioner is not always stable. More specifically, when sharp surfaces are involved in the geometrical model, our DP-ILUT method is more stable than the common used ILUTP5.Table 4
Comparison of memory requirement of EFIE-DP-ILUT and CFIE for 3D scattering.
Almond
Open cavity
Aircraft
Tank
EFIE-DP-ILUT (MB)
593
760
1,434
15,100
CFIE (MB)
623
—
2,106
21,585


## 5. Conclusion

In this paper, the EFIE with half-space dyadic Green’s function is presented to analyze the scattering from PEC targets of arbitrarily shape in a half space. Inspired by the idea of IEFIE, a novel diagonal perturbation of ILUT preconditioner is presented and investigated. It increases the robustness of LU decomposition and is more stable and physical than other existing diagonal shifting formulae. The proposed robust preconditioner is simple to be combined with the extended 3D p-FFT, making the EFIE formulation applicable in solving complex and electrically large problems in a half space.



---

*Source: 813273-2015-12-22.xml*