Session:  Session 3A10A  Fast and Multilevel Techniques (08j) 
Type:  Oral Antenna 
Date:  Wednesday, November 08, 2006 
Time:  08:30  12:20 
Room:  Gallieni B 
Chair:  
Cochair:  
Remarks: 
Seq  Time  Title  Abs No  
1  08:30 
Methods for the Accelerated Computation of Green`s Functions with 2D Periodicity in Layered Media
Baccarelli, P.^{1}; Galli, A.^{2}; Paulotto, S.^{2}; Valerio, G.^{2}; Volski, V.^{3}; Vandenbosch, G. A. E.^{3} ^{1}ITALY; ^{2}"La Sapienza" University of Rome, ITALY; ^{3}Katholieke Universiteit Leuven, BELGIUM A review of some recent methods for the efficient computation of twodimensional (2D) 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 bandgap structures (EBG), metamaterials, 2D phased planar arrays, and 2D planar leakywave radiators are common applications based on 2D periodic printed structures. In all these cases, a unit cell characterized by a complex geometry of the metallization requires the availability of a versatile fullwave 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 2D periodic vector and scalar mixedpotential Green`s functions in the space domain. In fact, by applying Floquet theorem, the generic component of the 2D 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 2D periodic vector and scalar mixedpotential 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 freespace Green’s function with 2D periodicity of a homogeneous medium. The spectral series is then evaluated by using Shanks’ transformation, while the freespace 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 nonperiodic spectraldomain mixedpotential Green`s functions. 

2  08:50 
Adaptive Integral Method for HigherOrder Hierarchical Method of Moments
Kim, O.S.; Meincke, P. Oersted.DTU, Electromagnetic Systems, Technical University of Denmark, DENMARK 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^{2 }), 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, higherorder 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 loworder basis functions (RWG, rooftop, or pulse). The obvious step further is to use the fast solvers with higherorder basis functions. In this paper, AIM is employed to accelerate the solution of the volume integral equation with higherorder 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 matrixvector product Z^{far}X, 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 loworder 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 higherorder 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 higherorder 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 loworder basis functions. It will be shown that, due to the higherorder convergence, the proposed technique requires less memory and computational time. 

3  09:10 
Fast Direct Solution of the Method of Moments Linear System
Heldring, A.; Tamayo, J. M.; Rius, J. M. Universitat Politecnica de Catalunya, SPAIN
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 subscatterers, often recursively (multilevel). The interactions between these subscatterers correspond to submatrices of the MoM impedance matrix. The acceleration methods exploit the well known fact that for well separated subscatterers, the corresponding submatrices are lowrank 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 matrixvector products at the core of all iterative solution methods for matrix equations. 

4  09:30 
Multilevel PO Algorithm for NonUniform Triangulations
Gurel, L. Bilkent University, TURKEY Recently, a fast physical optics algorithm (FPO) has been introduced. (Boag, IEEE Trans. Antennas Propagat., vol.52, pp. 197204, Jan. 2004). This algorithm provides a speedup 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 nonuniform 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 "nonuniform" algorithm, the radiation patterns of the larger triangles, which are too large to fit in the bottomlevel subdomains, are directly computed and incorporated in the appropriate aggregation levels. Reallife radar applications usually require the solution of electromagnetic scattering from electrically large bodies with complex geometries. These types of problems are usually solved using highfrequency 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, nonuniform 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.


5  09:50 
Efficient Parallelization of Multilevel Fast Multipole Algorithm
Ergul, O.; Gurel, L. Bilkent University, TURKEY We consider the solution of very large threedimensional problems of computational electromagnetics with the multilevel fast multipole algorithm (MLFMA). Accurate solutions of quite large problems can be performed iteratively with the accelerated matrixvector 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 reallife 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 matrixvector 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 tradeoff 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. 

6  10:40 
MultiResolution Enhancement of Fast Method of Moments Analysis of Antennas
Pirinoli, P.^{1}; Freni, A.^{2}; De Vita, P.^{2}; Vipiana, F.^{1}; Vecchi, G.^{1} ^{1}Politecnico di Torino, ITALY; ^{2}University of Florence, ITALY In the numerical analysis of largescale 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 "waveletlike" 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 subdomain 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_{S}T) 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., APS 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. 

7  11:00 
Incomplete LU Preconditioning for the ElectricField Integral Equation
Malas, T.; Gürel, L. Bilkent University, TURKEY The multilevel fast multipole algorithm (MLFMA) is widely used to solve large linear systems that are produced by the integralequation formulation of the computational electromagnetics (CEM). MLFMA reduces the complexity of the matrixvector multiplication to O(nlogn) time, allowing large systems to be solved iteratively. Among the various integral equations of CEM, the electricfield integral equation (EFIE) is the only formulation applicable to open geometries. The systems generated from EFIE are illconditioned; hence both the efficiency and the reliability of the iterative solver depend on the preconditioning operation. MLFMA keeps only the nearfield part of the dense coefficient matrix. One can factorize this sparse nearfield matrix and use it for preconditioning. However, due to the fillin 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 thresholdingbased 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. 387414.) With this remedy, we obtain convergence with iteration counts being close to those using exact factorization of the nearfield 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 twolevel 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 nearfield system with the inner solver, and use ILUT as the preconditioner for the nearfield 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. 

8  11:20 
Multilevel Field Interpolation Algorithm for Large PEC Objects
Espinosa, H. G.^{1}; Heldring, A.^{1}; Rius, J. M.^{1}; Mosig, J. R.^{2} ^{1}Universitat Politecnica de Catalunya, SPAIN; ^{2}Ecole Polytechnique Federale de Lausanne, SWITZERLAND
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 matrixvector 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^{2}) operation count required by direct matrixvector 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 2D large PEC objects discretized by the Method of Moments.
The algorithm is based on a recursive procedure where the 2D 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 (nontouching 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, JM. 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. 533543, March 1997.
[2] A. Boag, and C. Letrou, "Multilevel Fast Physical Optics Algorithm for Radiation from NonPlanar Apertures," IEEE Transactions on Antennas and Propagation, Vol. 53, No. 6, pp. 20642072, June 2005. 

9  11:40 
Preconditioning Techniques for the Method of Auxiliary Sources Applied to Geometries Characterized as Perturbations of a Circle
Vouldis, A. T.^{1}; Avdikos, G. K.^{1}; Anastassiu, H. T.^{2} ^{1}National Technical University of Athens, GREECE; ^{2}Hellenic Aerospace Industry, GREECE 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 quasidiagonal matrix is wellconditioned, and hence readily invertible via a numerical scheme. Results are shown for several, circlelike geometries, validating the preconditioning concept, and leading to high precision MAS solutions of complex scattering problems. Acknowledgments The Project is cofunded 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. 15411547. [2] A. G. Kyurkchan, B. Y. Sternin and V. E. Shatalov, "Singularities of Continuation of Wave Fields", Physics Uspechi, 39 (12), 1996, pp. 12211242. [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 Circularlike Metallic and Dielectric Objects", WSEAS Transactions on Communications, no 10, vol. 4, 2005., pp. 11381145. 

10  12:00 
Fast Monostatic RCS Computation in FEM Based Solvers Using QR Decomposition
Venkatarayalu, N.^{1}; Gan, Y. B.^{1}; Zhao, K.^{2}; Lee, J. F.^{2} ^{1}Temasek Laboratories, National University of Singapore, SINGAPORE; ^{2}ElectroScience Laboratory, The Ohio State University, UNITED STATES
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]. 