International Journal of Antennas and Propagation

Volume 2017 (2017), Article ID 9845050, 8 pages

https://doi.org/10.1155/2017/9845050

## Application of a Sparsity Pattern and Region Clustering for Near Field Sparse Approximate Inverse Preconditioners in Method of Moments Simulations

^{1}Computer Science Department, University of Alcalá, Madrid, Spain^{2}NewFasant S.L., Guadalajara, Spain

Correspondence should be addressed to Carlos Delgado

Received 18 May 2017; Revised 28 September 2017; Accepted 8 October 2017; Published 31 October 2017

Academic Editor: Nicolas Pinel

Copyright © 2017 Carlos Delgado et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

This document presents a technique for the generation of Sparse Inverse Preconditioners based on the near field coupling matrices of Method of Moments simulations where the geometry has been partitioned in terms of regions. A distance parameter is used to determine the sparsity pattern of the preconditioner. The rows of the preconditioner are computed in groups at a time, according to the number of unknowns contained in each region of the geometry. Two filtering thresholds allow considering only the coupling terms with a significant weight for a faster generation of the preconditioner and storing only the most significant preconditioner coefficients in order to decrease the memory required. The generation of the preconditioner involves the computation of as many independent linear least square problems as the number of regions in which the geometry is partitioned, resulting in very good scalability properties regarding its parallelization.

#### 1. Introduction

Many of the modern approaches for electromagnetic analysis based on the Method of Moments (MoM) [1] rely on the idea of only storing the near field coupling terms of the impedance matrix, which typically extends about one-quarter of the wavelength under analysis. The coupling effects between distant parts of the geometry are taken into account in the matrix-vector product computation as a part of the iterative solution process, and the specific manner in which this is carried out depends on each approach. Some popular schemes are based on the use of the Multilevel Fast Multipole Algorithm (MLFMA) [2], which includes the processes of aggregation, translation, and disaggregation of multipole expansions in the computation of such products. Other schemes use matrix compression techniques [3, 4] to compute these multiplications efficiently. All these methods circumvent the burden of storing the full MoM matrix, which would easily surpass the memory capacity and cripple the efficiency for the analysis of even moderately sized problems.

Since the full coupling matrix is no longer to be calculated and due to the size of the system to be computed, iterative solvers play a major role in modern simulation methods, although it is worthwhile to mention that some approaches extend the range in which direct solvers can be applied, such as those based on the use of Macro Basis functions [5], which reduce the number of unknowns of the original problem. The necessity of iteration, however, is unavoidable for electrically large or very large problems and with the use of iterative solvers arises the problem of slow convergence for some cases, due to reasons of different nature, such as geometrical features of the model, multiscale problems, or to the intrinsic electromagnetic behavior, such as the presence of resonant regions within the geometry. In any case the use of preconditioners is very important in order to ease the iterative process, which, in many cases, accounts for most of the total simulation time.

It is noteworthy to mention a family of preconditioners that are based on the physical properties of the problem and more specifically on the idea that inverse matrix that is approximately represented by the preconditioner is essentially an approximate solution of the electromagnetic problem. Previous works have been published for the efficient generation of physical-based preconditioners considering quasiplanar structures [6] or exploiting the quasiplanar electromagnetic behavior of conducting rough surfaces [7]. A method for the generation of a preconditioner for problems involving multiple scatterers, based on splitting the system matrix according to the types of material of the subdomains as well as currents on different parts, is shown in [8].

Another group of preconditioners that are also commonly used in the context mentioned above rely on the numerical manipulation of the system matrix rather than on physical properties of the problem under analysis. The Incomplete LU (ILU) decomposition [9–11] and the Sparse Approximate Inverse (SAI) [12–14] show a good behavior and are extensively used at the present time. The ILU preconditioner is considered slightly more efficient for the same amount of data [15], although the SAI approach is much more scalable, since it can be computed by solving independent linear least square (LLS) systems that generate the rows of the preconditioner. This property renders the SAI preconditioner very well suited to be used in very large problems that require a fair number of processing nodes. It should be noted, however, that modern ILU preconditioners based on multifrontal/multilevel schemes allow better performance than the conventional ILU approach while offering a good compromise between robustness and efficiency [16–18]. The application of multilevel ILU preconditioning techniques applied to electromagnetic problems involving the Electric Field Integral Equation has been documented in [19].

The preconditioners based on the ILU and SAI approaches are often generated considering only the near field part of the coupling matrix, which includes the strongest interactions between basis and testing functions. This allows a fast generation and reasonably reduced size, although the performance of such preconditioners may be lacking in the analysis of problems in which there are strong interactions between parts of the geometry that are physically distant, like the interaction between the feed of a reflector antenna and the main reflector or between parts of certain cavities.

The most common approach in the application of the SAI preconditioner involves the definition of a sparsity pattern, understood as the pattern of the elements of the preconditioner that are not null, which is the same as that of the near field coupling matrix [12]. In this work we propose an improved strategy based on the clustering of the least squares systems that fall under each near field region as well as a distance threshold parameter that allows a finer control of the sparsity pattern of the SAI preconditioner. We have obtained noticeable improvements on the convergence properties of the conventional use of the SAI preconditioner compared to the approach described here. Illustrative examples are given in the Numerical Results.

#### 2. Description of the Approach

The application of the MoM is based on the solution of a linear system of equations as follows:where and are the current coefficient and the excitation vectors, respectively, and matrix contains the coupling terms between each basis function and each testing function of the problem. The expression of term of can be written aswhere it is computed as the inner product of the th testing function and the field radiated by the th active basis function over the th subdomain. This field is given by the operator . Consequently, each row of matrix represents the coupling terms between all the active basis functions of the geometry and one passive testing function and, in turn, a column of contains the coupling terms between one active basis function and all the passive testing functions. Since is dense and unmanageable for many realistic problems it is very common to compartmentalize the original geometry into regions and only calculate and store the coupling terms of the basis and testing functions contained within the same region, as well as those contained within adjacent regions. The typical region size is about *λ*/4, which makes the matrix containing such terms, denoted as in this work, sparse and with a much more manageable size even for large problems. The rest of the coupling terms, noticeably weaker, are contained in the far field coupling matrix . Many modern efficient numerical methods do not require calculating the far field coupling matrix [2–4] and consider those interactions using alternative approaches.

In order to incorporate the near field preconditioning scheme proposed in the present work, it is useful to write the MoM equation separating the coupling matrix into the near field term () and the far field term () as shown as follows:where is the preconditioner. This matrix should ideally be as close as possible to the inverse of . However, it is very important to keep the generation procedure scalable and efficient both in terms of size and CPU time.

Under the assumption that the near field part of the coupling matrix is stronger than the far field term, it is possible to rewrite (3) in the following fashion:

##### 2.1. Region Clustering in the Generation of the Preconditioner

A commonly used SAI preconditioner relies on the generation of by minimizing the* Frobenius norm* of *,* where is the identity matrix. The scalability properties of this preconditioner arise from the possibility of decomposing such expression into the minimization of the norm of the difference between each row of the identity matrix and the preconditioner multiplied by the near field coupling matrix as follows:where indicates the total number of basis functions of the problem, stands for the th row of the identity matrix, and denotes the th row of the preconditioner .

It should be remarked that the technique described in this work relies on the compartmentalization of the computational space in terms of* regions*, so that only the coupling terms between subdomain functions located in the same or in adjacent regions are to be computed. This is a widely adopted scheme by many modern methods. Under this context, it is possible to take advantage of the reutilization of the same LLS matrix to obtain as many rows of the preconditioner as the number of the unknowns contained in each region. This can be illustrated by rewriting (5) in the following manner:where indicates the total number of regions of the problem. Considering at this point that is the number of subdomain functions contained in region and that represents the absolute index of the th subdomain contained in region , the matrix in (6) consists of rows of the unit matrixand, analogously, contains rows of the preconditioner:

##### 2.2. Parametric Sparsity Pattern

Noting that will be generally a full matrix it is necessary to define a sparsity pattern that enforces null values for some of its coefficients. It is common to see approaches in the literature that consider such a sparsity pattern to be identical as that of the original near field coupling matrix [12]. We introduce an alternative approach by defining in this context the* sparsity distance parameter * as the distance threshold that determines whether the subdomains contained in a region are to be included in the sparsity pattern. Let us denote for this purpose as the region in which subdomain is contained and as the distance between the centres of regions and , where an are two arbitrary subdomains. The sparsity mask matrix for region is defined asThe matrix, applied for all the regions, defines the sparsity pattern used to generate the preconditioner. With this consideration, the right hand side of expression (6) can be approximated aswhere is an approximation of obtained as the solution of the LLS where the matrices marked with the tilde overscore denote the application of the sparsity pattern by computing the Hadamard product with the* sparsity mask* and eliminating the null rows and columns of the result:In this expression denotes the element-wise matrix multiplication operator and represents the matrix that results from the elimination of the null rows and columns of . The term , in turn, is an operator that returns a filtered version of discarding the coefficients with a module lower than *τ* times the largest self-impedance in the region where they belong. The *τ* parameter, named* preprocessing threshold parameter* through the rest of this document, allows a reduction of the size of the LLS problems to be solved for the generation of the preconditioner by eliminating low-magnitude rows and columns, resulting in faster computation. Considering these comments can be written as follows:The preconditioner, therefore, can be computed by calculating the result of independent LLS problems, which can be distributed among an arbitrary number of processing nodes, and assembling the preconditioner using the partial matrices. It is convenient, however, to filter these results to retain only the elements that pose a significant weight for each row of the preconditioner. Experimental results have shown that for typical applications the size of the preconditioner can be reduced at least by an order of magnitude without a compromise in the convergence performance. Expression (12) illustrates the row filtering of introducing the* postprocessing threshold parameter *:where retains the elements with an absolute value over times the largest element in each row of the preconditioner:The sparsity distance parameter *ξ* introduced in this section allows a fine control of the size of the preconditioning matrix. Large values (similar to the size of the geometry) imply no restrictions in the sparsity pattern of the resulting preconditioner, which can lead to the computation of the exact inverse of with a prohibitive computational cost. Low values, around the size of the regions in which the geometry has been previously partitioned, render an approximate inverse with the same size as that can be computed fast. These sizes can then be modified by the filtering threshold that discards the least significant terms. Our experience, as shown in the next section, shows that sparsity distance parameter values ranging from 0.5 *λ* to present a good tradeoff between size and performance.

#### 3. Numerical Results

Some test cases have been selected in this section in order to illustrate the performance of the approach described in the present document. All the simulations have been run on an HP Z820 workstation using 16 Intel Xeon E5-2660 2.2 GHz processing cores and 128 GB of RAM. In all the cases we have applied the Method of Moments combined with the Multilevel Fast Multipole Algorithm using 0.25 *λ* as the region size unless otherwise indicated. The SAI preconditioners have been obtained distributing the computation of shown in expression (13) for all the regions of the problem between the computing nodes using the OpenMP paradigm [20]. The Generalized Minimal Residual (GMRES) iterative solver [21] with a restart parameter of 300 has been used for all the simulations presented in this section. The parallelization scheme of the SAI preconditioner followed in the simulations shown in this section consists of a loop that iterates over the regions in which the geometry has been previously partitioned. For each iteration the corresponding LLS matrix is assembled and the LAPACK numerical library [22] is used to obtain the solution. Since we use OpenMP to merely distribute the iterations of this loop over all the processing nodes, it is extremely important that the iterations are independent of each other, as is the case for the proposed preconditioner.

In the results presented in this section *θ* represents the polar angle measured from the -axis and *ϕ* denotes the azimuthal angle measured from the -axis, with 0° ≤ *θ* ≤ 180°, 0° ≤* φ* < 360°. These angular coordinates have been illustrated in Figure 1 for the sake of clarification.