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

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 geometricaloptics 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.


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][2][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 as ( 2 ).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][7][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][16][17][18][19][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(, )) 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  and  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][18][19][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.

Background: Integral Equations and Fast
Iterative Solver

Integral Equations in
Half-Space Problems.Consider a 3D PEC object of arbitrary shape located above a lossy half space as illustrated in Figure 1.Define   and   as the permittivity and permeability of layer  ( = 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 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 L  and K  are defined as [23] L  (r, r  ) where P.V. stands for the Cauchy principal value integration and where G  , G  , and   are the half-space dyadic and scalar Green's functions [24] and the vectors r = ρ + ẑ and r  =   ρ+  ẑ 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] IEFIE : where  is the combination coefficient varying between 0 and 1 and  0 is the wave impedance in free 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: The linear algebraic equation of IEFIE is then given in the following iterative way: where The integral is nonzero only when the supports of the test function   () and basis function   (  ) are overlapped.Thus, the total element number of Z IO is no more than 5, where  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].

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 halfspace Green's function can be termed as where || is the horizontal distance between field and source points.The primary term, which is a function of ( −   ), is The secondary term, a function of ( +   ), 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 nearfield 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]: 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  .For the secondary term, however, G grid (denoted as G  ) has only Hankel-two-level block-Toeplitz structure because the halfspace 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 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 the G  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.

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 ( 0.25 ) iterations in Krylov subspace methods while EFIE formulation converges International Journal of Antennas and Propagation in ( 0.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 matrix M as a nonsingular approximation of Z, the MVM can be written in the following form: 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(, ) 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 However, ILUT(, ) 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 Z EFIE NFI .The diagonal shifting formula for Z EFIE NFI can be written as [19] where ZEFIE NFI 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 ZEFIE NFI 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 This diagonal perturbation makes the system diagonally dominant and is more physical and stable than other shift parameters used in [18][19][20].It therefore increases the robustness of ILUT(, ) 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.

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  to be the average number of nonzero elements in a row of the near-field matrix.
It is noted that a strong indicator of the instability quality of the ILUT preconditioner is an estimate of ‖([][]) −1 ‖, called "condest" [18].This condition estimate is defined as which can be easily computed by using a forward substitution followed by a backward substitution.

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   = 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.
International Journal of Antennas and Propagation

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   = 6.38 − 0.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.
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.

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   = 6.5 − 0.6.The length, width, and height of the aircraft model are 12.5, 7, and 3 m, respectively, as shown in Figure 6.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.
The tank model is 7 × 2 × 1.9 m 3 .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.
To further demonstrate the advantage of the present method, Table 3 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   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.

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.

Figure 1 :
Figure 1: The geometry of a 3D PEC object above the half space.

Figure 3 :
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 :Figure 5 :
Figure 4: The convergence history of GMRES algorithms with different integral operators and preconditioners in the almond example.

Figure 6 :
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 :
Figure 7: The convergence history of GMRES algorithms with different integral operators and preconditioners in the aircraft model example.

Figure 8 :
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 :
Figure 9: The convergence history of GMRES algorithms with different integral operators and preconditioners in the tank model example.

Table 3 :
Comparison of different integral operators and preconditioners for 3D scattering.

Table 4 :
Comparison of memory requirement of EFIE-DP-ILUT and CFIE for 3D scattering.