A review of some recent methods for the efficient computation of two-dimensional (2-D) periodic Green`s functions in layered dielectric media is presented in this contribution. Such work is the result of a cooperation between "La Sapienza" University of Rome (SAPIENZA) and the Katholieke Universiteit Leuven (KUL) within a software integration activity of the Antenna Center of Excellence (ACE).

Frequency selective surfaces (FSS), electromagnetic band-gap structures (EBG), metamaterials, 2-D phased planar arrays, and 2-D planar leaky-wave radiators are common applications based on 2-D periodic printed structures. In all these cases, a unit cell characterized by a complex geometry of the metallization requires the availability of a versatile full-wave method for the analysis and the design process. The method of moments (MoM) in the spatial domain, in conjunction with a mixed potential integral equation (MPIE) formulation, can satisfy this requirement. Nevertheless, a crucial point in the analysis is related to the efficient derivation of the 2-D periodic vector and scalar mixed-potential Green`s functions in the space domain. In fact, by applying Floquet theorem, the generic component of the 2-D periodic Green`s function in the space domain can be expressed in terms of a spectral sum of an infinite number of spatial harmonics that is slowly convergent.

With the aim of obtaining a software for the efficient computation of 2-D periodic vector and scalar mixed-potential Green’s functions in multilayer media, this series has been accelerated. The asymptotic part of the addends is determined, when source and observation points are in the same layer, as the spectral Green`s function of an asymptotically equivalent homogeneous medium and it is subtracted from the original spectral series. The result is a spectral series that is rapidly convergent, and an asymptotic series that corresponds to the free-space Green’s function with 2-D periodicity of a homogeneous medium. The spectral series is then evaluated by using Shanks’ transformation, while the free-space Green`s function is accelerated either by using Kummer`s decomposition, Poisson`s formula, and Shanks’ transformation or by applying Ewald`s transformation. The proposed approach is then extended to the case of source and observation points that lie in different layers. This acceleration is coded in a routine by SAPIENZA, while KUL provides the non-periodic spectral-domain mixed-potential Green`s functions.

Integral equations applied to scattering problems for arbitrarily shaped conducting and inhomogeneous dielectric objects are usually solved with the Method of Moments (MoM). This procedure leads to a dense system of linear equations with the memory requirement O(N² ), with N being the number of unknowns. It is evident that computational demands increase drastically as the problem size grows. In the case of compact penetrable volumetric scatterers, practically available computer resources limit the scatterer size to an order of one wavelength.

Fast integral equation solvers, such as the Multilevel Fast Multiple Method or the Adaptive Integral Method (AIM), are able to reduce the memory demands for volumetric problems to the order of O(N logN) and O(N), respectively. Alternatively, higher-order basis functions can be employed in the conventional MoM to significantly reduce the number of unknowns, which in many practical cases is more memory and computationally efficient as compared to the fast solvers based on low-order basis functions (RWG, rooftop, or pulse). The obvious step further is to use the fast solvers with higher-order basis functions.

In this paper, AIM is employed to accelerate the solution of the volume integral equation with higher-order hierarchical basis functions. The main idea behind AIM is to split the MoM matrix in two parts, responsible for nearby and far interactions, as Z=Z^near+Z^far. Only the elements of Z^near are computed and stored explicitly. To account for the far interactions, the basis functions are expanded in terms of Dirac delta functions defined at nodes of a regular rectangular grid. Then, the computation of the matrix-vector product Z^farX, with X denoting the unknown solution vector, at each iteration of the iterative solution process is sped up by application of the FFT. The accuracy of the expansion is controlled by the expansion order and the support size b of the basis function. To assure a reasonable expansion order, kb should be less than 1, where k is the wavenumber in the background medium. In the case of low-order basis functions (RWG or rooftop), for which AIM was first developed, they are assumed to span two neighbor cells. Hence, the support b is the sum of the corresponding cell sizes. To retain the accuracy of the expansion for the higher-order basis functions, which become advantageous when they are defined in relatively large cells, we modify the AIM expansion procedure, so that the two parts of the first- and higher-order rooftop basis functions are expanded separately in each cell. Thus, we effectively reduce the support b in the expansion to the size of a single cell, which implies that larger cells with higher expansion orders can be utilized.

Numerical examples will be given to validate the proposed technique as well as to compare its memory requirements and computational costs with AIM based on low-order basis functions. It will be shown that, due to the higher-order convergence, the proposed technique requires less memory and computational time.

In recent years, a wide range of methods have been developed for accelerating the solution of electromagnetic scattering and radiation problems with the Method of Moments (MoM). Many of these methods rely on the subdivision of the problem geometry into sub-scatterers, often recursively (multilevel). The interactions between these sub-scatterers correspond to sub-matrices of the MoM impedance matrix. The acceleration methods exploit the well known fact that for well separated sub-scatterers, the corresponding sub-matrices are low-rank and can be effectively compressed, using either techniques based on the problem physics, like the Fast Multipole Method [1] and the Matrix Decomposition Method [2] or purely mathematical methods, for example [3]. The compressed impedance matrix is then used in the repeated matrix-vector products at the core of all iterative solution methods for matrix equations.
However, the efficiency of the iterative solution of a linear system depends on the condition of the coefficient matrix. Certain problems can yield badly conditioned matrices, leading to very slow convergence, or even failure to converge. Furthermore, the behaviour is generally difficult to predict. Therefore, direct solution of the linear system remains a viable alternative for problems small enough to allow for it.
In this paper we propose a way to accelerate the direct solution of the MoM linear system using any one of the above mentioned techniques. Our method is based on an algorithm we previously published in [4] and [5] in the framework of building an out-of-core preconditioner for iterative methods. It is a block LU decomposition algorithm, the numerically dominant operations of which are products of sub-matrices. When the latter are replaced by their compressed counterparts, the compression is at least partly retained in the resulting block LU factors. A reduction of a factor five to ten, both in computation time and storage requirements can be achieved for some typical scattering problems.

[1] R. Coifman, V. Rokhlin and S. Wandzura, 'The Fast Multipole Method for the Wave Equation: A Pedestrian Prescription,' IEEE A&P Magazine, Vol 35, No 3 June 1993
[2] E. Michielsen and A. Boag, 'A Multilevel Matrix Decomposition Algorithm for Analyzing Scattering from Large Structures,' IEEE Trans A&P, Vol 44, No 8 Aug. 1996
[3] S. M. Seo and J.F. Lee, 'A Single Level Low Rank IE-QR Algorithm for PEC Scattering Problems Using EFIE Formulation,' IEEE Trans A&P, Vol 52, No. 8, Aug. 2004
[4] A. Heldring, J.M. Rius, L.P. Ligthart and A. Cardama, 'Accurate Numerical Modeling of the Tara Reflector System,' IEEE Trans A&P, Vol. 52, No. 7, July 2004
[5] A. Heldring, J.M. Rius and L.P. Ligthart, 'New Block Ilu Preconditioner Scheme for Numerical Analysis of Very Large Electromagnetic Problems,' IEEE Trans Magn. Vol 38, No. 2, March 2002.

Recently, a fast physical optics algorithm (FPO) has been introduced. (Boag, IEEE Trans. Antennas Propagat., vol.52, pp. 197-204, Jan. 2004). This algorithm provides a speed-up for computing the physical optics (PO) integral over complex bodies for a range of aspect angles and frequencies. In this paper, this algorithm is further developed in order to compute the scattered field from non-uniform triangulations of complex bodies. In the original "uniform" FPO algorithm, only the radiation patterns of the smallest subdomains in the bottom level are directly computed and the radiation patterns of the larger subdomains in the upper levels are computed via interpolation and aggregation. In this modified "non-uniform" algorithm, the radiation patterns of the larger triangles, which are too large to fit in the bottom-level subdomains, are directly computed and incorporated in the appropriate aggregation levels.

Real-life radar applications usually require the solution of electromagnetic scattering from electrically large bodies with complex geometries. These types of problems are usually solved using high-frequency techniques, such as PO, physical theory of diffraction (PTD), etc. Since obtaining an explicit solution of the PO integral for such complex bodies is difficult, dividing the surface of such bodies into triangles and evaluating the PO integral on each triangle is commonly used. Although the choice of the triangle size is an open question, the requirement is the accurate representation of the surface curvature by the triangulation. To satisfy this requirement with a minimum number of triangles, non-uniform triangulation, where the triangle size is determined by the surface curvature, can be used.

Following the triangulation of the target surface, the collection of triangles is recursively divided into cubic subdomains until the cube size is in the order of a wavelength. Radiation patterns of each triangle are combined at the appropriate aggregation level to obtain the overall pattern of the whole target.

(a) (b) (c) (d) Figure 1. (a) Example of a non-uniform mesh. (b) Cubic subdomains in the 1st level. (c) Cubic subdomains in the 2nd level. (d) Cubic subdomains in the 3rd level.

We consider the solution of very large three-dimensional problems of computational electromagnetics with the multilevel fast multipole algorithm (MLFMA). Accurate solutions of quite large problems can be performed iteratively with the accelerated matrix-vector products by the MLFMA since this algorithm reduces the complexity of the processing time and the memory requirement to O(NlogN). However, most of the real-life applications require the modeling of the geometries with millions of unknowns. Therefore, in addition to the advances in the solution algorithms, it becomes essential to take advantage of the advances in the computer hardware in order to solve these very large problems. For this purpose, we construct parallel clusters of relatively inexpensive processors connected via special fast networks.

Since the MLFMA requires the evaluation of all electromagnetic interactions, it is not trivial to distribute the problem among the processors. Communications between the processors and duplications of the computations reduce the efficiency of the parallelization. For this reason, we investigate the separate parts of the matrix-vector products to understand and solve the inefficiencies. Specifically, we consider balancing the distribution of the tree structure among the processors for the aggregation and disaggregation steps by using optimizations based on actual measurements. A trade-off between data replication and communication is also investigated. For the translation step, we consider the overlapping of the communications and computations, and also investigate various communication schemes to synchronize the processors.

In the parallelization of the MLFMA, we note that the efficiency is strongly related to the geometry of the problem. Therefore, it is natural to choose different strategies for the solution of different problems. As an example, when the geometry of the problem has a large dimension in one direction so that its modeling requires a relatively low number of unknowns, the setup time becomes as important as the iterative solution. A good parallelization strategy should take into consideration every potential bottleneck of the algorithm.

In the presentation, we will provide a detailed report of our efforts for the solution of very large electromagnetic problems accurately and efficiently with parallel MLFMA.

In the numerical analysis of large-scale problems as those involving antennas and arrays, the use of iterative solvers is mandatory. A number of techniques, commonly named "fast methods", have been introduced in the past years with the aim of drastically reducing the numerical complexity, given in terms of computation time (MoM matrix filling time plus MoM linear system solution time), and required memory occupation. These methods generally act on the cost and memory occupation needed at each step of the iterative solver algorithm, but they do not address the issue of the number of iterations, required to achieve a given accuracy, that could instead reduced by the use of a suitable preconditioner.

Standard preconditioning methods of linear algebra have smooth advantages when applied to fast methods, since they require the storage of a (very large) system matrix, or at least of the inversion of the (sparse) matrix that account for the near field interactions. Here, the preconditioning issue will be addressed by having recourse to the Multiresolution (MR) scheme (Vipiana et al., IEEE Trans. AP, July 2005) that manages to generate "wavelet-like" hierarchical multiscale vector functions for any meshed geometry and well controls the MoM condition number, by the simply application of a diagonal preconditioner. The MR functions are expressed in terms of sub-domain functions (in the present case, the RWG functions defined on a triangular mesh), therefore their generation leads to the construction of a changing basis matrix T that is very sparse (of the order of O(Nlog2N)) and frequency independent.

When used in conjunction with a fast method, the MR preconditioner is S_S = TD_S^-1/2 where D_S = diag(TTZ_ST) and Z_S is the (sparse) matrix that represents the near field interactions. The generation of the diagonal matrix D_S requires a computational effort of the order of O(Nlog2N) and an additional storage of O(N). Since S_S shares the same memory space of the matrix T, the overall storage needed by S_S is O(Nlog2N), while its computational cost is O(N[log2N]2). Finally, the additional complexity introduced by the MR preconditioning at each step of the iterative method is of the order of O(3Nlog2N). This complexity is comparable or smaller than the complexity of the majority of the fast MoM methods and is well compensated by the increased convergence of the iterative solver .

The MR preconditioner have been first applied to the Sparse Matrix/Adaptive Integral Method (SM/AIM) (De Vita et al., MMET02) for the analysis of large arrays (De Vita et al., AP-S 2004), obtaining a total reduction of the numerical complexity with respect to the standard MoM of almost three orders of magnitude. We are now extending the use of the MR preconditioner to the FMM, for which the MR scheme is particularly suited: preliminary results show a behavior similar to that obtained with the SM/AIM.

The multilevel fast multipole algorithm (MLFMA) is widely used to solve large linear systems that are produced by the integral-equation formulation of the computational electromagnetics (CEM). MLFMA reduces the complexity of the matrix-vector multiplication to O(nlogn) time, allowing large systems to be solved iteratively. Among the various integral equations of CEM, the electric-field integral equation (EFIE) is the only formulation applicable to open geometries. The systems generated from EFIE are ill-conditioned; hence both the efficiency and the reliability of the iterative solver depend on the preconditioning operation.

MLFMA keeps only the near-field part of the dense coefficient matrix. One can factorize this sparse near-field matrix and use it for preconditioning. However, due to the fill-in phenomenon, this approach becomes infeasible. Hence, instead of using the exact factors, some of the entries are discarded and the arising incomplete factors are used as a preconditioner. This approach, known as incomplete LU (ILU) preconditioning, has proven to be successful in many applications. (Yousef Saad, Iterative Methods for Sparse Linear Systems, SIAM, Philadelphia.)

In this paper, we show that a thresholding-based ILU preconditioner (ILUT), is effective in most cases. However, for some geometries, no progress is observed towards convergence. We show that this occurs because of the instable incomplete factors. To overcome the problem, we apply column pivoting as suggested in a previous work. (E. Chow and Y. Saad, Experimental study of ILU preconditioners for indefinite matrices, J. Comput. Appl. Math., 86 (1997), pp. 387-414.) With this remedy, we obtain convergence with iteration counts being close to those using exact factorization of the near-field matrix as a preconditioner.

To further increase the effectiveness of the preconditioner, we can use a flexible solver, which allows the preconditioner to vary from step to step. Then, preconditioning of the original system can be accomplished using another iterative solver. This approach results in a two-level iterative solver: The outer solver is used for the solution of the original system, and the inner solver is used for preconditioning the original system. We propose to solve the near-field system with the inner solver, and use ILUT as the preconditioner for the near-field system. Our experiments reveal that, by using only a few iterations for the inner system, the number of outer iterations further drops compared to using ILUT alone.

Multilevel acceleration methods have been intensely studied over the last two decades with the purpose to combine an iterative solution with an efficient performance of matrix-vector products, providing considerably reduced computational requirements [1]

These methods typically require a number of operations per iteration, for very large objects, of the order of O(NlogN) for a problem with N degrees of freedom, which is much less than the O(N²) operation count required by direct matrix-vector multiplication.

In this paper, we develop an NlogN algorithm using a numerically efficient evaluation of fields produced by given current distributions. It is applied to the fast solution of 2-D large PEC objects discretized by the Method of Moments.

The algorithm is based on a recursive procedure where the 2-D scatterer surface is decomposed into a set of subdomain cells. A cell enclosing the object is subdivided into smaller cells at multiple levels, in the form of an octal tree. The fields produced by each of them are computed separately by means of the interaction between source and observation cells.

For the near cell interaction (touching cells) at the finest level, the field is evaluated through direct evaluation of the Electric Field Integral Equation, while for the far cell interaction (non-touching cells) at multiple levels, the field is evaluated at one or a few evaluation points inside it, and then interpolated over the field basis functions.

Before the interpolation, the phase common to all source points in the subdomain cell is removed to cancel rapid oscillations due to the Green’s function of the EFIE. After the interpolation the phase is redistributed and finally the fields evaluated in every subdomain cell are aggregated into the total field. The procedure is similar to multilevel computation of physical optics integral proposed by A. Boag [2]

In the paper, we outline the proposed multilevel algorithm, demonstrate its NlogN operation count both theoretically and numerically, and assess its accuracy as a function of the governing parameters.

[1] W. C. Chew, J-M. Jin, C.-C. Lu, E. Michielssen, and J. M. Song, "Fast Solution Methods in Electromagnetics," IEEE Transactions on Antennas and Propagation, Vol. 45, No. 3, pp. 533-543, March 1997.

[2] A. Boag, and C. Letrou, "Multilevel Fast Physical Optics Algorithm for Radiation from Non-Planar Apertures," IEEE Transactions on Antennas and Propagation, Vol. 53, No. 6, pp. 2064-2072, June 2005.

Performance optimization of the Method of Auxiliary Sources (MAS) remains an important and challenging issue. It is strongly related to the location of the auxiliary sources (AS's) and has been resolved so far only for a few canonical geometries (i.e. [1]). In general, the auxiliary surface should enclose the singularities of the scattered field as tightly as possible. Although the exact position of these singularities is unknown for arbitrary geometries, the radius of the smallest circle, or sphere (in 2D and 3D respectively), encompassing them all, can be determined analytically [2]. An auxiliary surface, conformal to the scatterer surface, circumscribing this circle/sphere is expected to yield optimal results, i.e. minimize the boundary condition error. In [3] it is inferred that this should indeed be the case for several geometries, including an elliptic cylinder, where singularities coincide with the foci. However, the actual error behavior in such a situation is usually unclear, due to additive numerical noise. As a matter of fact, when the AS's lie close to singularities, the condition number of the MAS linear system rises dramatically, rendering numerical inversion inaccurate, and effectively increasing the computational error. To alleviate this problem, in this paper matrix preconditioning is employed for 2D configurations, which can be considered as perturbations of a circle. Given that the geometry is close enough to an exact circle, the corresponding MAS matrix is also a perturbation of the matrix pertaining to the circle. Since the inverse of the latter is known analytically [1], its multiplication with the original, perturbed matrix can be performed exactly. The product yields a perturbation of the unity matrix, which is obviously well conditioned, and easily invertible numerically. Alternatively, the original matrix can be multiplied from the right with the eigenvector matrix of the circular case and from the left with its transpose. This bilateral multiplication yields an almost diagonal matrix, whose most significant elements are essentially perturbations of the eigenvalues of the circular case. Again, this resulting quasi-diagonal matrix is well-conditioned, and hence readily invertible via a numerical scheme. Results are shown for several, circle-like geometries, validating the preconditioning concept, and leading to high precision MAS solutions of complex scattering problems.

Acknowledgments The Project is co-funded by the European Social Fund (75%) and National Resources (25%) (Herakletus).

[1] H. T. Anastassiu, D. G. Lymperopoulos and D. I. Kaklamani, "Accuracy Analysis and Optimization of the Method of Auxiliary Sources (MAS) for Scattering by a Circular Cylinder", IEEE Trans. on Antennas and Propagation vol. 52, no. 6, June 2004, pp. 1541-1547. [2] A. G. Kyurkchan, B. Y. Sternin and V. E. Shatalov, "Singularities of Continuation of Wave Fields", Physics Uspechi, 39 (12), 1996, pp. 1221-1242. [3] H. T. Anastassiu, A. T. Vouldis and G. K. Avdikos, "Optimization Schemes for the Method of Auxiliary Sources Applied to the Scattering Problem from Circular-like Metallic and Dielectric Objects", WSEAS Transactions on Communications, no 10, vol. 4, 2005., pp. 1138-1145.

Computation of radar cross section is a vital step in the design process of antenna arrays or shape optimization of electromagnetically invisible bodies. Though integral equations based methods are often more efficient for such problems, they require special treatment when the geometry involved is inhomogeneous, which is often the case in practical design problems. The finite element method (FEM) with its inherent capability to model inhomogeneous media is therefore a potential alternative [1].
In FEM with total field formulation for scattering problems, plane wave excitation is imposed on the boundary of the computational domain enclosing the body whose scattered fields are to be computed. Such an excitation is represented on the right hand side (excitation vector) of the resultant matrix system. For excitation at different angles of incidence, the FEM system matrix remains the same so long as the finite element meshes are unchanged. Based on this fact, we demonstrate a simple and efficient technique to reduce the CPU time involved in the computation of monostatic RCS. Using the excitation vectors for N_p angles of incidence, an excitation matrix Y can be assembled. This excitation matrix is often highly rank deficient. We perform a QR-decomposition based on the modified Gram-Schmidt orthogonalization[2]. The number of columns in Q is equal to the rank r of Y, and typically r<< N_p. Thus, only r matrix solutions are needed to compute the monostatic RCS over the range of incident angles instead of N_p. The matrix solution time is therefore reduced by a factor of r/N_p. It is noted that for the total field formulation, each column of Y has non-zero entries only for the degrees of freedom on the domain boundary over which the incident plane wave is injected. The memory requirement for the Y matrix is O(n^2/3). For the matrix solution, we use a p-type multiplicative Schwarz preconditioned iterative solver[3]. In the case of FEM codes with adaptive mesh refinement, an optimal mesh for all N_p angles of incidence can simply be generated by performing error estimation on the solution obtained with the row sum of Y as the excitation vector. It is fairly easy to extend the proposed technique to existing FEM codes. Details of the QR-decomposition and the speed up achieved using the technique for monostatic RCS computation from antenna arrays will be shown in the presentation.

REFERENCES
[1] Sylvester, P.P. and Ferrari, R.L., Finite Elements For Electrical Engineers. 3rd Ed., Cambridge University Press.
[2] Golub, G. H and Van Loan, C. F, Matrix Computations. 3rd Ed., The Johns Hopkins University Press
[3] Lee, J.-F. and D. K. Sun, IEEE Trans. MTT, 864, 2004