Abstract
An efficient and memory saving method based on pseudoskeleton approximation (PSA) is presented for the effective and accurate analysis of finite periodic structures. Different from the macro basis function analysis model, our proposed method uses the formulations derived by the local Rao-Wilton-Glisson basis functions. PSA is not only used to accelerate the matrix-vector product (MVP) inside the single unit but also adopted to decrease the calculation burden of the coupling between the different cells. Moreover, the number of decomposed coupling matrices is minimized due to the displacement invariance of the periodic property. Consequently, even compared with the multilevel fast multipole algorithm (MLFMA), the new method saves much more memory resources and computation time, which is also demonstrated by the numerical examples.
1. Introduction
Periodic structures have recently found wide applications in the electromagnetic engineering such as antenna arrays and metamaterials with negative permittivity and negative permeability. Hence, the accurate and efficient analysis of periodic structures becomes quite essential. If the periodic structure is an infinite array, simple methods can be applied based on Floquet’s theorem [1] or periodic Green’s function [2], where only a single cell of the periodic structure is the domain of interest.
However, all periodic structures have finite size in real-life problems, although the size may be very large. Therefore, the numerical algorithms which accurately consider the mutual coupling between all cells should be used if the accurate results are required. The method of moments (MoM) [3] and its fast algorithms such as fast multipole method (FMM) [4], adaptive cross approximation (ACA) [5], and FFT-based methods [6] are flexible approaches to study the surface problems. However, the efficiency of numerical methods is still limited since the periodic property is not used in the algorithm framework. Recently, a hybrid method combining the accurate MoM and periodic method of moment (PMM) [7] has been proposed which can gain the balance between the two methods. Moreover, some physically based entire-domain basis functions [8] have been developed to reduce the number of unknowns. Further, the FMM and FFT techniques are integrated to accelerate the calculation [9].
Compared to the ACA method, pseudoskeleton approximation (PSA) [10] is also an efficient low-rank-based algebraic fast algorithm which makes it a really competitive alternative. In this paper, we propose an efficient method with low-memory requirement based on PSA to perform the analysis for finite periodic structures effectively and accurately. In consideration of the accuracy of the mutual interactions [8] and the simplicity of implementation, our proposed method uses the formulations derived from the local basis functions instead of macro basis function (MBF) [11, 12]. In this paper, PSA is not only used to accelerate the matrix-vector product (MVP) inside the single unit but also adopted to decrease the calculation burden of the coupling between the different cells. Moreover, the number of decomposed coupling matrices is minimized due to the displacement invariance of the periodic property. With these improvements, an efficient method with low-memory usage of finite periodic objects can be achieved. Several numerical examples are given to show the priority of the proposed method compared to the conventional multilevel fast multipole algorithm (MLFMA) [13] for periodic structures.
2. MoM and PSA Formulation
In this section, the basis principle of MoM and PSA is briefly introduced at first. Then, the choice of arguments in PSA is discussed.
2.1. MoM Equations and Its Fast Algorithms
Consider a time-harmonic electromagnetic wave scattering or radiation problem of an arbitrary perfect electrically conducting (PEC) object. The object is excited by an incident electric field , then the electric field integral equation (EFIE) associated with the surface equivalent currents can be expressed by where is the tangential unit vector on the surface of the object and and stand for the angular frequency and permeability, respectively. is 3D scalar Green’s function. The linear system of MoM is obtained by discretizing the unknown vector with Rao-Wilton-Glisson (RWG) [14] basis functions and applying Galerkin’s testing method. Let represent the above EFIE matrix system. The MVP process can be accelerated by fast algorithms which can be written as follows: where is the matrix of near field interactions which are directly computed and stored and stands for couplings between the far-field interactions which will be accelerated together with in the iterative solving process.
2.2. Basic PSA Frameworks
According to the low-rank decomposition, the far-field interaction matrix with rank-deficient property can be approximated by a product of two much smaller submatrices and : where is the effective rank and satisfies . While in the skeleton approximation (SA) theory [10], there is a nonsingular submatrix in . Denote the rectangular matrices as and which contain the columns and rows of , respectively, then is expressed as where and have dimensions with and , respectively. In the PSA method, is reevaluated as where is the pseudoinverse of with dimensions of , then columns and rows are chosen from to get the and . Three aspects need to be noted here: (i) is the inverse of and the computation of inverse is very expensive; (ii) the determination of the value of is a balance between accuracy and efficiency, where is a number large enough so that most important bases will be embedded; (iii) how to choose those working columns and rows of and .
For the (ii) and (iii) problems, we will discuss them in the next subsection. For the first problem, assuming that can be decomposed via singular value decomposition (SVD) as where and are unitary matrices, is a diagonal matrix with nonnegative real numbers, and is the complex conjugate transpose of . In the actual implementation, the dimension of , , and can be further decreased by a preset threshold [15]. Let , , and represent the reduced submatrix of , , and , respectively. Then, calculation of the pseudoinverse of (i.e., ) is straightforward:
By combining (5) and (7), the original far-field interaction matrix can be decomposed as
2.3. Choice of Arguments in PSA
As mentioned in the previous subsection, the choice of the value of and the specific sampling rows and columns determine the performance of PSA. In the randomized PSA (RPSA) [15], is equal to . Then, the problem of (ii) transforms into how to estimate the rank of the original matrix. In [16], Chai and Jiao give the approximation of rank of the 3D EM problem by , where is the studying wave number. In this paper, the rank is approximated by where is the diameter of the bounding box corresponding to the octree structure and is a preset positive parameter. The larger the is, the more accurate the matrix decomposition is. In this paper, when is set as 3, the satisfactory accuracy can be guaranteed.
For the problem of (iii), instead of using random numbers in RPSA, we use a strategy analogous to ACA in this paper. Firstly, initialize from the 0th row as the first row index. Then, find the largest entry in this row, and the corresponding column value in which this entry is located is chosen as the next column index. Similarly, find the largest entry in the current column and get the next row index which should be different from all previous row indexes. This process is carried out iteratively until rows and columns are found and stored. Please refer to [17] for more information.
3. Proposed Method for Periodic Structures
We consider a case of arbitrarily shaped PEC patch, for example, refer in Figure 1. The total impedance matrix contains two parts: self-coupling blocks and mutual coupling blocks. Therefore, both the blocks are analyzed and decomposed by PSA to gain better efficiency and low-memory usage.
3.1. PSA for Decomposition of Self-Coupling Matrix
In the previous mentioned periodic algorithms, the self-coupling block elements are calculated by a direct MoM method which is not efficient for a large array unit. Therefore, PSA is used to decompose the self-coupling matrix. As the same to all fast algorithms of MoM, the low-rank decompositions based on PSA are performed on the far-field groups while the near-field elements are directly computed and saved. Moreover, some operations are taken out to further improve the efficiency of PSA. Firstly, the pseudoinverse of the (i.e., ) will not be directly calculated since the complexity of SVD is very high. Instead, the LU decomposition will be performed when the dimension grows up. Different from the previously proposed PSA, , , and the LU decomposition of will be stored. Moreover, and in (8) are further decomposed by ACA technique. Finally, the original far-field coupling matrix can be decomposed into six submatrices which is showed by the following formula: where (or ) is the ACA decomposition of (or ) and is the LU decomposition of . Note that the inverse of LU matrix of is not directly computed. Actually, the LU back substitute is performed after the MVP based on the matrix in each MVP process.
3.2. PSA for Decomposition of Mutual Coupling Matrices
In the traditional fast algorithms, the matrix decompositions are only performed on the far-field groups especially for the physically related methods such as MLFMA. However, the mutual interactions between two different cells (off-diagonal blocks) are all computed by PSA in our implementation. This treatment has been also used in [18] and demonstrated to be more efficient but losing little accuracy. Moreover, the displacement invariance property is explored in our implementation. Since the mutual coupling remains the same when the distance between two groups’ centers is not changed, the corresponding matrix will be calculated and stored only once under the index of relative distance. In the MVP process, the stored coupling matrices may be used several times when the relative distances are the same. By using this scheme, the number of mutual coupling matrices can be decreased dramatically. For example, the original number of off-diagonal matrix blocks of Figure 1 is while the reduced number will be 62 if the displacement invariance is used.
3.3. Preconditioner Considerations
The preconditioning stage should be also considered carefully while dealing with complex structures. In this paper, the preconditioning matrix for impedance matrix is built based on the self-coupling matrix .
Moreover, the inner-outer iteration scheme can be also applied when facing ill-conditioned matrix systems such as EFIE for radiation analysis. In our implementation, the inner iteration is the solution of the self-coupling matrices. Certainly, the inner solution area can be also extended to improve the convergence rate of outer iteration.
4. Numerical Results
In this section, several numerical experiments are presented to demonstrate the efficiency and low-memory requirement of the proposed method. For all the simulations, the mesh sizes are no less than 0.1λ (λ represents the wave length). The biconjugate gradient stabilized method (BiCGSTAB) is adopted as the iterative solver for the matrix equations, and the threshold of the iterative stopping criterion for residuum is set to 0.001. The sparse approximate inverse preconditioner is used in all simulations. All the computations were carried out on a workstation with four Xeon E5-4620 CPU and 256 GB of RAM with OpenMP technique, and the digits were stored in double precision.
4.1. Scattering Simulation
In order to verify the accuracy and efficiency of the proposed method, we first consider a periodic structure consisting of element cells (also shown in Figure 1), where the unit cell is a 0.5 m PEC sphere. The working frequency is set to be 1 GHz which leads to 20286 unknowns (degrees of freedom) in each patch. Hence, the total number of unknowns is . The gap between two unit cells is set to be 0.5 m. The radar cross sections (RCS) computed by MLFMA, the proposed method (periodic PSA), and the commercial software (FEKO) are illustrated in Figure 2(a). It could be clearly seen that the numerical results from the proposed method have excellent agreement with both the conventional MLFMA and FEKO. While as in Table 1, both the computation time and peak memory usage by periodic PSA are less than the MLFMA method. What is more, another case is considered when the gap between two unit cells is decreased to be 0.1 m. The RCS results and the use of computer resources are shown in Figure 2(b) and Table 1, respectively. Although the rank of mutual coupling matrix gets bigger due to the tight coupling interactions, the proposed method still shows the good efficiency and low-memory usage with losing little accuracy.
(a)
(b)
Furthermore, the scattering of a large-scale periodic structure consisting of element cells is considered. The working frequency is set to be 1 GHz. Different radii of the sphere are considered which lead to about 0.81, 3.38, and 5.46 million numbers of unknowns corresponding to 0.1 m, 0.2 m, and 0.25 m models, respectively. The distance between two centers of neighboring sphere is set to be a fixed value of 1.0 m. Table 2 shows the information of models and the numerical performance between MLFMA and the periodic PSA. It can be seen that for a large-scale problem, the proposed PSA method still needs less computation time and much less memory requirement. It should be noted that in the 3rd example in Table 2 when the sphere radius is set to 0.25 m, the finest box in MLFMA has to be set to 0.20λ to avoid low-frequency breakdown. Figure 3 shows that the RCS results of 0.25 m sphere models computed from the periodic PSA still agree well with FEKO and MLFMA.
4.2. Radiation Simulation
Lastly, we consider a radiation problem, which contains antenna cells with delta-gap excitations. Each cell has 1145 unknowns, and the total number is . The gap between two cells in this model is 0.005 m which is smaller than 0.1λ. Since EFIE is the only choice for radiation problems, the inner-outer iteration scheme is adopted to accelerate the calculation process. Figure 4 shows the current density distribution of the antenna array. The far-field gain patterns computed by the proposed method and MLFMA are shown in Figure 5.
Obviously, the two methods give nearly the same results. However, as in Table 3, the computation time has been reduced from 419 seconds to 172 seconds in solving matrix equations for the radiation problem. Moreover, the memory usage for MLFMA is 8841 MB while the value is only 358 MB for the proposed method. Again, the proposed method performs much better than MLFMA. The first reason is that only 360 off-diagonal blocks are need to be calculated and stored in our method while the original value is . Secondly, the mesh size is about 0.01λ for the antenna cell (dense mesh). As it is known, the lowest-level box size is limited to be no less than 0.20λ since MLFMA suffers from the low-frequency breakdown, which results in a heavy near-field matrix consumption.
5. Conclusion
In this paper, an efficient and memory saving method based on pseudoskeleton approximation (PSA) for the effective and accurate analysis of finite periodic structures is presented. The appropriate choice of sampling rows and columns as well as sampling number in PSA is discussed and confirmed. Based on the algebraic fast algorithm PSA, we have carefully handled the self-coupling and mutual-coupling blocks of the original impedance matrix. Moreover, the number of decomposed coupling matrices is minimized by the employment of the displacement invariance of the periodic property. The simulation results of large-scale scattering and radiation examples show that the proposed periodic PSA method needs less computation time and much less memory compared to conventional MLFMA. Hence, numerical results demonstrate the efficiency and superiority of the proposed method.
Data Availability
No additional data are available.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This work was supported in part by the National Natural Science Foundation of China under Grant 61171035.