# A Hybrid Technique Based on the Combination of Multilevel Fast Multipole Algorithm and the Geometrical Theory of Diffraction

Authors: María-Jesús Algar; Javier Moreno; Iván González; Lorena Lozano; Felipe Cátedra
Journal: International Journal of Antennas and Propagation
 (2014)

---

## Abstract

This paper proposes a hybrid technique for treating electromagnetic problems of scattering and radiation in which the source structure is described as an array of antennas. This strategy is based on the combination of the rigorous method multilevel fast multipole algorithm (MLFMA) and the high frequency technique geometrical theory of diffraction (GTD). Thanks to the use of MLFMA, the source can be discretized into several cubic regions considering each of them as a source point in order to reduce the number of times required to compute the ray tracing when GTD is applied to obtain the scatter field. In this analysis, objects with complex shapes are described by using nonuniform rational B-splines (NURBS) which is a very common way to model geometrical bodies. Numerical results that demonstrate the accuracy and efficiency in terms of CPU time are shown.


---

## Body

## 1. Introduction

Traditionally, the method of moments has been used to solve scattering or radiation problems using integral equations [1]. This method typically applies a sampling rate of about ten divisions per wavelength to discretize the object and partition it into subdomains, in order to guarantee the accuracy of the results. The induced currents over each subdomain are expressed in terms of basis functions such as Rao-Wilton-Glisson functions [2] or rooftops [3]. The application of this method to solve realistic problems, that is, electrically large problems, is limited by the size and the density of the impedance matrix due to the dense sampling of objects. Therefore, to reduce the computational bottleneck of this method, in terms of memory requirements and CPU time, different approaches have been developed.On the one hand, several rigorous techniques have been developed. One of the most common approaches is the fast multipole method (FMM) [4, 5]. This method only needs to store the near-field terms of the coupling matrix and compute the far-field interactions. Following this schema, this method can improve the complexity of matrix-vector products from O(N2) to O(N3/2). Its multilevel implementation, MLFMA [6], reduces the computational complexity to O(NlogN) since it allows a fast matrix-vector product computation in the iterative solution process. However, these improvements in terms of CPU time are not enough to avoid the constraint leading to the high sampling rate imposed by the MoM discretization. Additionally, both methods demand memory requirements to store the near-field coupling terms.The main objective of FMM and MLFMA is to increase the computation efficiency. However, an alternative solution to the limitation of MoM is to reduce the number of unknowns replacing the type of subdomains basis functions by a set of macrobasis functions [7–10]. In this sense, it is worthwhile to remark the characteristics basis function method (CBFM) [11, 12]. This alternative involves the decomposition of the geometry into blocks. Over each block a set of macrobasis functions (CBFs) is defined in terms of functions of small subdomains. Consequently, the use of these functions reduces the coupling matrix and the computation of the electrical field can be handled applying direct solvers. Another option is the synthetic function expansion (SFX) [13] approach which reduces the complexity of the MOM analysis of complex geometries dividing the overall structure into portions, recognizing that the degrees of freedom of the field coupling between the solutions on these portions are limited, and building basis functions that exploit this property.However, the improvement achieved by the CBFM is not enough when the electrical size of the structure under analysis is very large since the matrix may become very large and the memory needed to store the reduced matrix can represent a problem as well. Therefore, an iterative method is needed to solve the reduced matrix equation. In [14, 15], the stabilized biconjugate gradient method is used to solve this issue. Authors suggested a numerically efficient approach based on the combination of the CBFM and MLFMA. The reduced matrix is computed containing only the near-field terms and the far-field terms are computed using MLFMA. This strategy allows us to reduce the number of unknowns as well as to optimize the memory storage.On the other hand, several research works propose the application of high frequency techniques such as geometrical optics (GO) [16] or physical optics (PO) [17] or the hybridization of these techniques with the conventional method of moments [18, 19].In this work, we propose an alternative to take advantage of the MLFMA with the aim of computing electrically large problems applying the high frequency technique geometrical theory of diffraction (GTD) [20], in which the radiation structure is composed of an array of multiple antennas. In order to solve this task, the strategy developed is based on the combination of the MLFMA together with GTD. The objective of this approach is to reduce the number of times needed to obtain the ray tracing to compute the electrical field of the array of antennas. In this sense, this kind of problems can be solved with a lower CPU-time compared to the direct application of GTD.The geometrical description of the cases studied in this paper has been modeled using nonuniform rational B-splines (NURBS) [21]. This description is well suited to set up an adaptive mesh to the real object with less information than required with planar meshing elements.This paper is organized as follows. Section2 describes the main basis and the formulation of the MLFMA. This is followed, in Section 3, by an explanation of the combination of MLFMA and GTD to speed up the computation of the electrical field when an array of multiple antennas is considered. Section 4 shows the results obtained with the presented approach in order to validate it. Finally, Section 5 summarizes the contents of the paper and some conclusions about the performance and efficiency of the strategy are developed.


## 2. Multilevel Fast Multipole Method Formulation

The application of MLFMA entails the compartmentalization of the geometry into several cubic regions. Figure1 shows for a given geometry how its volume is divided into certain cubes. This process begins with a small cubic size which is typically a quarter of wavelength, as recommended in [5].Figure 1
Region division of a generic geometry in the MLFMA.This method is implemented in three steps: the aggregation, the translation, and the disaggregation step. In the aggregation step, for the first level of cubic regions, the electrical field of each cube is obtained and associated with its geometrical centre, since this point has been chosen to be the aggregation point. This aggregation term is obtained according to Expression (2). In the following steps, this first level of cubic regions is aggregated into bigger ones. The aggregation terms of the lower cubes contained in each cube of the higher level are taken into account by using shifting and interpolation to reduce the computational resources.When the aggregation step is finished, the translations between cubes which are separated but belong to the same level are carried out. Finally, when all the cubes have received the contributions from the rest of cubes in the same level, these contributions are released to the children cubes by means of shifting andanterpolation (as described in [22]) in the disaggregation step.The computation of the matrix-vector products is achieved following this expression:(1)Aij=∫VmjAGG(k^)τmm′(k^,rmm′⟶)Vm′iDIS(k^)d2k^,
where VmjAGG(k^) is the aggregation term, Vm′iDIS(k^) is the disaggregation term, and τmm′ is the translation term and k^ is the direction vector for every multipole-algorithm angular sample.The aggregation term is calculated as follows:(2)VmjAGG(k^)=∑k=1Nm∫u∫ve-jk^r→jm(I--k^k^)Tk(u,v)dudv,
where Nm is the number of subdomain functions within the cubic region, r→jm is the vector that extends from the sampling point to the aggregation point, and Tk(u,v) is the k-order low-level basis function (in our case we use generalized rooftops) defined over each block.Analogously, the disaggregation term can be calculated as(3)Vm′iDIS(k^)=∑k=1Nm∫u∫vejk^r→im′(I--k^k^)Rk(u,v)dudv,
where Rk(u,v) is the k-order testing function (we use generalized razor-blade) defined over each block.The translation term between pointsm (active group) and m′ (victim group) is given by
(4)τmm′(k^,r→mm′)=jk4π∑l=0Ljl(2l+1)hl(1)(krmm′)Pl(r^mm′·k^),
where hl(1)(x) is a spherical Hankel function of the first kind and Pl(x) is a Legendre polynomial defined as follows:
(5)Pl(x)=12πi∮(1-2tx+t2)-1/2t-l-1dt.


## 3. Use of MLFMA Together with GTD

The main goal of the strategy proposed in this paper is to solve radiation problems in which the source structure is composed of an array of multiple antennas. If we try to solve this kind of problems applying GTD directly, this technique needs too much CPU time to obtain the ray tracing for each antenna of the array for the computation of the electrical field radiated by the whole structure. Therefore, to reduce the computational requirements, a new approach based on the combination of the MLFMA and GTD has been developed. In terms of efficiency, this new method allows us to decrease the number of times in which it is necessary to obtain the ray tracing. However, in terms of accuracy, it is worthwhile to remark that every part of the structure shall be at least at a distanceD of the closer multipole, D being twice the size of the cube considered to extract each multipole. As indicated in [22], this distance guarantees enough accuracy in the MLFMA computations of fields.If the radiation pattern of an array of several antennas is obtained applying MLFMA, it is possible to compute the scatter field radiated by the whole geometry considering the source structure as a set of radiation diagrams. Each radiation diagram corresponds to each cube generated by MLFMA.In this sense, the first step in the computation of the electrical field is the compartmentalization of the source structure into cubical regions. Given an array of antennas like the one shown in Figure2, MLFMA is applied to obtain the cubical regions that cover the whole structure, as has been described in the aggregation step. According to Expression (2), the aggregation term is computed for each cube and the far field radiated by the cube is associated to its geometrical center. This information is stored in an external file. Therefore, the source structure can be treated as a set of source points, locating each of them at the center of its respective cube. Figure 2 shows a schema of this process.Figure 2
Discretization of a generic array of antennas.Once we have computed the electrical field for each source applying a fast ray tracing algorithm [23], it is easy to obtain the total field due to the whole source structure. It is worthwhile to remark that the computation of the ray tracing can be done taking into account multiple iterations between the surfaces of the geometry without forfeiting reduced computational resources.The number of cubical regions depends on the size of the selected window. A typical value of the region size is2λ or 4λ. As the number of cubical regions is lower than the number of antennas that compose the array, the ray tracing algorithm has to be applied fewer times. In this way, this strategy effectively reduces the CPU time required to compute the scatter field of this kind of structures.


## 4. Computational Results

In order to give an idea of the computational time saving achieved with the combination of MLFMA and GTD, the radiation pattern of two representative examples has been computed.


### 4.1. Hispasat Satellite

The first test case is the Hispasat satellite whose mock-up with a size of 1.6 × 1.4 × 1.5 m is depicted in Figure3. The geometrical model considered in this analysis is composed of 68 NURBS. The source consists of an array of 3 × 3 conical horns built to work at 10 GHz and a distance of 0.05 m between the centers of their apertures.Figure 3
Geometrical model of Hispasat with an array of 3 × 3 conical horns.The radiation pattern at a frequency of 10 GHz and for the set of directions in theφ = 0° cut and sweep in θ from θ = −90° to θ = 90° with 181 directions has been computed. The numerical results shown in Figure 4 compare the values obtained for the strong component applying three different methods: the traditional GTD (considering reflection and diffraction effects), the MoM-MLFMA, and the combination of MLFMA together with GTD (also considering reflection and diffraction effects).Figure 4
Radiation pattern of the Hispasat at 10.0 GHz, cutφ = 0°.Some conclusions are drawn after analyzing these results. The accuracy of the results obtained with the new technique is very good since the shape of the graphic is very similar to the results obtained with the other two methods. In the computation of the radiation pattern applying GTD, obtaining the ray tracing for the whole set of conical horns is necessary. However, considering a size window of4λ, the number of sources obtained with the new technique has been reduced from 9 to 7. Therefore, in this case, the ray tracing has only been computed for 7 source points.The main advantage of the new technique is the reduction in the CPU time, as it can be observed in Table1, due to the reduction in the number of sources. This test case has been simulated in an Intel Xeon at 2.2 GHz and 128 GB of RAM with 8 processors.Table 1
Comparison of the CPU time for the analysis of the Hispasat test case.
GTD
MoM-MLFMA
MLFMA and GTD
CPU time
36 min
1 h 16 min
19 min


### 4.2. Passat Car

The next test case consists of the realistic model of the Passat car, whose dimensions are 4.71 × 1.77 × 1.45 m, shown in Figure5. The radiation pattern of this structure has been computed at a frequency of 10 GHz. The source is composed of an array of 7 × 7 vertical dipoles built at 10 GHz and with a separation of 0.05 m.Figure 5
Geometrical model of the Passat with an array of 7 × 7 vertical dipoles.Figures6 and 7 depict the numerical results for φ = 0° and φ = 90° cuts, respectively, with a sweep in θ from θ = −90° to θ = 90°, considering 181 directions. GTD (considering reflection and diffraction effects), MoM-MLFMA, and the combination of MLFMA together with GTD (also considering reflection and diffraction effects) have been applied to obtain the radiation pattern of the Passat shown in Figure 5.Figure 6
Radiation pattern of the Passat at 10.0 GHz, cutφ = 0°.Figure 7
Radiation pattern of the Passat at 10.0 GHz, cutφ = 90°.In this case, considering a size window of4λ, the number of sources obtained with the new technique has been reduced from 49 to 2. Therefore, this reduction in the number of sources implies a reduction in the execution time, as it can be observed in Table 2. This test case has been simulated in an Intel Xeon at 2.20 GHz and 128 GB of RAM with 8 processors.Table 2
Comparison of the CPU time for the analysis of the Passat test case.
GTD
MoM-MLFMA
MLFMA and GTD
CPU time
2 h 37 min
3 h 59 min
9 min


## 4.1. Hispasat Satellite

The first test case is the Hispasat satellite whose mock-up with a size of 1.6 × 1.4 × 1.5 m is depicted in Figure3. The geometrical model considered in this analysis is composed of 68 NURBS. The source consists of an array of 3 × 3 conical horns built to work at 10 GHz and a distance of 0.05 m between the centers of their apertures.Figure 3
Geometrical model of Hispasat with an array of 3 × 3 conical horns.The radiation pattern at a frequency of 10 GHz and for the set of directions in theφ = 0° cut and sweep in θ from θ = −90° to θ = 90° with 181 directions has been computed. The numerical results shown in Figure 4 compare the values obtained for the strong component applying three different methods: the traditional GTD (considering reflection and diffraction effects), the MoM-MLFMA, and the combination of MLFMA together with GTD (also considering reflection and diffraction effects).Figure 4
Radiation pattern of the Hispasat at 10.0 GHz, cutφ = 0°.Some conclusions are drawn after analyzing these results. The accuracy of the results obtained with the new technique is very good since the shape of the graphic is very similar to the results obtained with the other two methods. In the computation of the radiation pattern applying GTD, obtaining the ray tracing for the whole set of conical horns is necessary. However, considering a size window of4λ, the number of sources obtained with the new technique has been reduced from 9 to 7. Therefore, in this case, the ray tracing has only been computed for 7 source points.The main advantage of the new technique is the reduction in the CPU time, as it can be observed in Table1, due to the reduction in the number of sources. This test case has been simulated in an Intel Xeon at 2.2 GHz and 128 GB of RAM with 8 processors.Table 1
Comparison of the CPU time for the analysis of the Hispasat test case.
GTD
MoM-MLFMA
MLFMA and GTD
CPU time
36 min
1 h 16 min
19 min


## 4.2. Passat Car

The next test case consists of the realistic model of the Passat car, whose dimensions are 4.71 × 1.77 × 1.45 m, shown in Figure5. The radiation pattern of this structure has been computed at a frequency of 10 GHz. The source is composed of an array of 7 × 7 vertical dipoles built at 10 GHz and with a separation of 0.05 m.Figure 5
Geometrical model of the Passat with an array of 7 × 7 vertical dipoles.Figures6 and 7 depict the numerical results for φ = 0° and φ = 90° cuts, respectively, with a sweep in θ from θ = −90° to θ = 90°, considering 181 directions. GTD (considering reflection and diffraction effects), MoM-MLFMA, and the combination of MLFMA together with GTD (also considering reflection and diffraction effects) have been applied to obtain the radiation pattern of the Passat shown in Figure 5.Figure 6
Radiation pattern of the Passat at 10.0 GHz, cutφ = 0°.Figure 7
Radiation pattern of the Passat at 10.0 GHz, cutφ = 90°.In this case, considering a size window of4λ, the number of sources obtained with the new technique has been reduced from 49 to 2. Therefore, this reduction in the number of sources implies a reduction in the execution time, as it can be observed in Table 2. This test case has been simulated in an Intel Xeon at 2.20 GHz and 128 GB of RAM with 8 processors.Table 2
Comparison of the CPU time for the analysis of the Passat test case.
GTD
MoM-MLFMA
MLFMA and GTD
CPU time
2 h 37 min
3 h 59 min
9 min


## 5. Conclusions

In this paper, we present the development of a hybrid technique based on the combination of the MLFMA and GTD to obtain the radiation pattern when an array of several antennas is considered. The main purpose of this technique is to reduce the CPU time required to obtain the ray tracing for each antenna of the array without losing accuracy in the results. In this sense, the aggregation step of the MLFMA reduces the number of source points, which implies a reduction in the number of times that is required to compute the ray tracing. The numerical results obtained applying the new approach show good agreement with the traditional high frequency technique GTD or MoM-MLFMA. In spite of these methods, this simple and efficient technique reduces the CPU time required to solve this kind of problems obtaining a high level of accuracy in the results.



---

*Source: 972584-2014-07-02.xml*