Multilevel Fast Multipole Algorithm for Characteristic Mode Analysis
Qi I. Dai, Jun Wei Wu, Ling Ling Meng, Weng Cho Chew, Wei E. I. Sha
11 Multilevel Fast Multipole Algorithm forCharacteristic Mode Analysis
Qi I. Dai, Jun Wei Wu, Ling Ling Meng, Weng Cho Chew,
Fellow, IEEE,
Wei E. I. Sha
Abstract —Characteristic mode (CM) analysis poses challengesin computational electromagnetics (CEM) as it calls for efficientsolutions of dense generalized eigenvalue problems (GEP). Multi-level fast multipole algorithm (MLFMA) can greatly reduce thecomputational complexity and memory cost for matrix-vectorproduct operations, which is powerful in iteratively solving largescattering problems. In this article, we demonstrate that MLFMAcan be easily incorporated into the implicit restarted Arnoldi(IRA) method for the calculation of CMs, where MLFMA withthe sparse approximate inverse (SAI) preconditioning techniqueis employed to accelerate the construction of Arnoldi vectors. Thiswork paves the way of CM analysis for large-scale and compli-cated three-dimensional ( -D) objects with limited computationalresources. Index Terms —Characteristic mode, generalized eigenvalueproblem, multilevel fast multipole algorithm, preconditioner.
I. I
NTRODUCTION A FTER its humble beginnings in the 1970s, characteristicmode (CM) theory have gained a recent resurgenceof interest in the field of antenna design and optimization.Initiated by Garbacz and refined by Harrington and Mautz[1–3], CM theory was popularized in the antenna communityby the work of Cabedo-Fabres as it has been shown promisingfor systematic antenna design [4]. Such a systematic approachbecomes increasingly favored since heuristic approaches basedon experience and intuition of antenna engineers hardly sur-vive the demanding requirement of designing complex antennasystems. CM analysis is capable of characterizing arbitraryobject’s radiation and scattering properties only relying on itsgeometry and material properties rather than source config-uration, which makes it a favorable candidate for systematicdesign methods where complicated functionality tradeoffs canbe handled. Moreover, CM analysis provides useful physicalinsight into antenna operation as well as guidance on theexcitation of desired radiation modes [5–7].Although the use of CM analysis has shown success inantenna design, relatively little effort has been devoted tothe numerical computation of CMs. Development of fastCM solvers becomes indispensable for applications such as
This work was supported in part by the National Science Foundation underAward Number 1218552, in part by the Research Grants Council of HongKong (GRF 711609 and 711508), and in part by the University Grants Councilof Hong Kong (Contract AoE/P-04/08).Q. I. Dai, L. L. Meng and W. C. Chew are with the Department of Electricaland Computer Engineering, University of Illinois at Urbana − Champaign,Urbana, IL 61801 USA.J. W. Wu is with the School of Electronics and Information, WuhanUniversity, Hubei, 430072, China.W. E. I. Sha is with the Department of Electrical and Electronic Engineer-ing, the University of Hong Kong, Hong Kong. antenna synthesis where sources are integrated to electricallylarge platforms, and multiscale modeling that contains finestructures. Current strategy can only deal with small problemsas it explicitly generates the impedance matrix by applyinga Galerkin’s procedure to the electric field integral equation(EFIE). The memory cost is O ( N ) where N is the number ofunknowns. When N is in the order of or , one can findall eigenpairs (eigenvalues and eigenvectors) of the dense gen-eralized eigenvalue problem (GEP) by Schur decompositionwith an overall computational complexity of O ( N ) . However,in most instances, only a small portion of the spectrum isdesirable, which can be efficiently solved for using Krylovsubspace iterations such as Lanczos and Arnoldi algorithms.A number of matrix-vector products are performed in suchmethods, each of which has a computational complexity of O ( N ) . Even so, this O ( N ) complexity hinders the use ofCM theory in large object analysis where N may be in theorder of or even larger.The fast multipole method (FMM) was originally proposedby Rokhlin for fast particle simulations where static integralequations are solved [8]. It was later extended for the calcu-lation of electromagnetic scattering by many research groups[9–12]. A two-level FMM reduces the complexity of a matrix-vector product and memory cost from O ( N ) to O ( N . ) .A nonnested two-level FMM incorporated with ray propa-gation physics further reduces the complexity to O ( N / ) [13]. Further reduction in complexity and memory cost wasachieved by the development of the multilevel fast multipolealgorithm (MLFMA). The complexity becomes O ( N log N ) when MLFMA is implemented using signature function, inter-polation, and filtering [14]. Based on translation, interpolation,anterpolation, and a grid-tree data structure, MLFMA reducesthe complexity and memory cost to O ( N log N ) [10–12].In this paper, MLFMA is successfully adopted in CManalysis. We demonstrated that MLFMA can be easily incor-porated into iterative eigensolvers such as the implicit restartedArnoldi (IRA) method [15] and Jacobi-Davidson QZ (JDQZ)method [16]. The sparse approximate inverse (SAI) of theexplicitly available submatrix for near group interaction isused as a preconditioner [17] when Arnoldi vectors in IRAare constructed . Other preconditioning schemes such as theCalderon multiplicative preconditioner (CMP) [18] is brieflydiscussed, which is promising for CM analysis in even larger-scale applications.II. C HARACTERISTIC M ODE F ORMULATION
The widespread CM formulation relies on the EFIE for anarbitrarily shaped perfectly electric conductor (PEC) object a r X i v : . [ m a t h . NA ] D ec given by ikη ˆ n × (cid:90) S d r (cid:48) G ( r , r (cid:48) ) · J ( r (cid:48) ) = − π ˆ n × E inc ( r ) , r ∈ S (1)where η = (cid:112) µ/(cid:15) is the characteristic impedance of free-space,and the Green’s dyadic is G ( r , r (cid:48) ) = (cid:20) I − ∇∇ (cid:48) k (cid:21) g ( r , r (cid:48) ) (2)with g ( r , r (cid:48) ) = e ikR R , R = | r − r (cid:48) | (3)and ˆ n is the unit normal of the conducting surface.In the method of moments (MoM), the unknown surfacecurrent J ( r ) is approximated with the popular Rao-Wilton-Glisson (RWG) basis functions f j ( r ) as J ( r ) = N (cid:88) j =1 I j f j ( r ) (4)where I j are the expansion coefficients. The impedance opera-tor relates the induced current J ( r ) to the incident electric field ˆ n × E inc ( r ) on S , and is projected to form the dense impedancematrix ¯ Z using a Galerkin’s procedure where testing functionsare chosen to be the same as basis functions as Z ij = ikη (cid:10) f i , G , f j (cid:11) S = ikη (cid:90) s d r f i ( r ) · (cid:90) s G ( r , r (cid:48) ) · f j ( r (cid:48) ) d r (cid:48) (5)As suggested by Harrington and Mautz [3], the CMs ofarbitrary PEC objects can be obtained by solving the GEPgive by ¯ XJ n = λ n ¯ RJ n (6)where ¯ R and ¯ X are the resistance and reactance matrices,respectively. They are the real and imaginary parts of ¯ Z , i.e., ¯ Z = ¯ R + i ¯ X (7)The eigenpairs λ n and J n correspond to the characteristic val-ues and characteristic currents, respectively. The characteristicvalue λ n is important as | λ n | indicates the modal behavior.When λ n = 0 , the corresponding J n is at resonance whichis efficient in radiating energy. When | λ n | is large, J n is aninefficient radiating mode which stores energy in the object.More explicitly, when λ n > λ n < , J n stores magnetic(electric) energy.III. M ULTILEVEL F AST M ULTIPOLE A LGORITHM
A. Addition Theorem
The standard MLFMA is based on the addition theorem.Letting r i and r j be the observation point and source pointlocated in group m and m (cid:48) , respectively, the spatial vectorfrom the source point to the observation point is written as[12] r ij = r i − r j = r i − r m + r m − r m (cid:48) + r m (cid:48) − r j = r im + r mm (cid:48) + r m (cid:48) j (8)where r m and r m (cid:48) represent the centers of groups m and m (cid:48) .If groups m and m (cid:48) are not nearby, i.e., | r im + r m (cid:48) j | < | r mm (cid:48) | , their interaction is considered as far interaction, and thespherical-wave function can be expanded using the additiontheorem as e ikr ij r ij = ik π (cid:90) d ˆ k e i k · ( r im + r m (cid:48) j ) α mm (cid:48) ( k , r mm (cid:48) ) (9)where r ij = | r i − r j | , and the translator α mm (cid:48) is defined as α mm (cid:48) ( k , r mm (cid:48) ) = L (cid:88) l =0 i l (2 l + 1) h (1) l ( kr mm (cid:48) ) P l (cid:16) ˆ k · ˆ r mm (cid:48) (cid:17) (10)with L as the truncation number of an infinite series, h (1) l ( · ) the first kind spherical Hankel function of order l , and P l ( · ) the Legendre polynomial of degree l . B. Two-Level Algorithm
In the two-level algorithm, the matrix-vector product iscalculated as [12] N (cid:88) j =1 Z ij I j = (cid:88) m (cid:48) ∈ B m (cid:88) j ∈ G m (cid:48) Z ij I j + ik π (cid:90) d ˆ k V fim (ˆ k ) · (cid:88) m (cid:48) / ∈ B m α mm (cid:48) ( k , r mm (cid:48) ) · (cid:88) j ∈ G m (cid:48) V sm (cid:48) j (ˆ k ) I j , i ∈ G m (11)where V fim (ˆ k ) = (cid:90) S d r (cid:16) ¯ I − ˆ k ˆ k (cid:17) · f i ( r im ) e i k · r im (12) V sm (cid:48) j (ˆ k ) = (cid:90) S d r (cid:16) ¯ I − ˆ k ˆ k (cid:17) · f j ( r m (cid:48) j ) e i k · r m (cid:48) j (13)and G m denotes all elements in group m , B m denotes allnearby groups of group m (including itself). The first termcorresponds to the contribution of near interaction where theimpedance submatrix is explicitly generated, and the secondterm corresponds to the contribution of far interaction evalu-ated by the fast multipole method (FMM). Moreover, for thefar interaction, the summation (cid:80) j ∈ G m (cid:48) V sm (cid:48) j I j represents theaggregation of all unknowns in group m (cid:48) to the group center,the summation (cid:80) m (cid:48) α mm (cid:48) represents the translation fromgroup m (cid:48) to m , and V fim · · · represents the disaggregationfrom the center of group m to all unknowns in this group. C. Multilevel Algorithm
The two-level FMM is extended to a multilevel algorithm,i.e., MLFMA, to further reduce the cost of matrix-vectorproduct operations. In MLFMA, the object S is first enclosedby a box that contains S , which is referred to as level- . Thisbox is then partitioned into eight smaller ones to form level- .Each level- box is again partitioned into eight smaller onesto form level- , and recursively continued until the finest level L f with a box size of . λ to . λ is formed. This procedureresults in an oct-tree structure, where level- is taken as thecoarsest level for best efficiency. The far interaction at the coarsest level is implemented using MLFMA where interpo-lation and anterpolation are adopted in the aggregation fromthe finest level to the coarsest level and the disaggregationfrom the coarsest level to the finest level, respectively. Thisis because the radiation patterns become increasingly richeras one progresses from the finest level to the coarsest level.The radiation patterns are sampled to perform the integrationover the Ewald sphere in (11). A well-documented approach isto apply the uniform sampling and Gauss-Legendre samplingin the φ and θ coordinates, respectively, and local Langrangeinterpolation between levels [12].IV. MLFMA FOR
CM A
NALYSIS
For large scale CM analysis, the desired eigenpairs of(6) are normally solved iteratively. The ARnoldi PACKage(ARPACK), based on the implicitly restarted Arnoldi (IRA)method or the Lanczos variant for symmetric matrices, isappropriate for calculating a few eigenpairs of large sparseor structured matrices, where a matrix-vector product requiresonly O ( N ) floating point operations [15]. However, whenARPACK is applied to the dense matrix pencil (cid:0) ¯ X , ¯ R (cid:1) in(6), the matrix-vector product is challenged by the O ( N ) complexity and storage requirement.In most applications, the eigenvalues of great interest arethose with small magnitudes which correspond to efficientradiating modes. Since it is easy for the IRA method toconverge to large eigenvalues, the desired eigenpairs can befound by solving the following standard eigenvalue problem ¯ X − ¯ RJ n = 1 λ n J n (14)The Arnoldi process for (14) calls for fast matrix-vectorproducts in the form of ¯ X − ¯ Ru with an arbitrary vector u to construct Arnoldi vectors. Even though matrices ¯ R and ¯ X can be explicitly generated if tremendous memory is available,the matrix vector product ¯ X − v where v = ¯ Ru can onlybe iteratively computed, as the complexity of inverting ¯ X is O ( N ) .It is well known that MLFMA reduces the complexity andmemory cost of performing ¯ Zu to O ( N log N ) . The standardMLFMA can be easily modified to expedite the requiredmatrix vector products. To calculate ¯ Xu with an arbitrary u , the following plane-wave decomposition is used provided | r im + r m (cid:48) j | < | r mm (cid:48) | cos( kr ij ) r ij = ik π (cid:90) d ˆ k e i k · ( r im + r m (cid:48) j ) α S mm (cid:48) ( k , r mm (cid:48) ) (15)where α S mm (cid:48) ( k , r mm (cid:48) ) = i L (cid:88) l =0 i l (2 l + 1) y l ( kr mm (cid:48) ) P l (cid:16) ˆ k · ˆ r mm (cid:48) (cid:17) (16)with the spherical Neumann function y l ( · ) of order l . Further-more, the reactance matrix ¯ X is given by X ij = (cid:61) m [ Z ij ] = kη (cid:28) f i , (cid:20) I − ∇∇ (cid:48) k (cid:21) cos ( kR ) R , f j (cid:29) S (17) The required matrix-vector product in terms of a two-levelalgorithm is obtained as N (cid:88) j =1 X ij u j = (cid:88) m (cid:48) ∈ B m (cid:88) j ∈ G m (cid:48) X ij u j + ik π (cid:90) d ˆ k V fim (ˆ k ) · (cid:88) m (cid:48) / ∈ B m α S mm (cid:48) ( k , r mm (cid:48) ) · (cid:88) j ∈ G m (cid:48) V sm (cid:48) j (ˆ k ) u j , i ∈ G m (18)where G m , B m , V fim and V sm (cid:48) j are defined the same asbefore. The corresponding multilevel algorithm simply followsthe same scheme used in the standard MLFMA.The calculation of ¯ Ru is less involved as there is noviolation of the following addition theorem sin( kr ij ) r ij = ik π (cid:90) d ˆ k e i k · ( r im + r m (cid:48) j ) α R mm (cid:48) ( k , r mm (cid:48) ) (19)where α R mm (cid:48) ( k , r mm (cid:48) ) = − i L (cid:88) l =0 i l (2 l + 1) j l ( kr mm (cid:48) ) P l (cid:16) ˆ k · ˆ r mm (cid:48) (cid:17) (20)with the spherical Bessel function j l ( · ) of order l . Moreover,the resistance matrix ¯ R is given by R ij = (cid:60) e [ Z ij ] = − kη (cid:28) f i , (cid:20) I − ∇∇ (cid:48) k (cid:21) sin ( kR ) R , f j (cid:29) S (21)The required matrix-vector product in terms of a two-levelalgorithm is obtained as N (cid:88) j =1 R ij u j = (cid:88) j ∈ G m R ij u j + ik π (cid:90) d ˆ k V fim (ˆ k ) · (cid:88) m (cid:48) (cid:54) = m α R mm (cid:48) ( k , r mm (cid:48) ) · (cid:88) j ∈ G m (cid:48) V sm (cid:48) j (ˆ k ) u j , i ∈ G m (22)where the FMM calculation involves both the nearby andnonnearby groups of group m (except for itself). With regardto the multilevel algorithm, the standard MLFMA bookkeepingcan still be employed for the least modification of the existingcode. However, a more efficient implementation conductstranslation between nearby groups which belong to the sameparent group at each level. Besides, the coarsest level is takento be level- instead of level- .When the problem scale is large, the reactance matrix ¯ X becomes ill-conditioned. It may take many iterations for thegeneralized minimal residual (GMRES) method or biconjugategradient stabilized (BiCGSTAB) method to converge to aprescribed error tolerance as ¯ X − u is solved. On the otherhand, efficient preconditioning schemes for solving ¯ X − u are not well documented. An alternative is to consider theequivalent eigenvalue problem of (6) given by ¯ ZJ n = (1 + iλ n ) ¯ RJ n (23) Since λ n are real numbers in the CM theory, λ n close to zerocan be easily found with IRA by solving eigenvalues withlarge magnitude of the standard eigenvalue problem ¯ Z − ¯ RJ n = 11 + iλ n J n (24)Although IRA may require more precision in solving ¯ Z − u than ¯ X − u to yield accurate eigenpairs, (24) is still favored inthis study as there exist efficient preconditioners for accelerat-ing the EFIE solutions ¯ Z − u . These preconditioners may notwork well in preconditioning ¯ X − u as the Green’s function cos( kR ) /R is noncausal. An SAI preconditioner is imple-mented in this study following [17] due to its simplicity. InMLFMA, the impedance submatrix ¯ Z near for near interactionis explicitly evaluated and stored. The sparse preconditioningmatrix ¯ P with a prescribed sparsity pattern is constructed byminimizing (cid:13)(cid:13) ¯ I − ¯ Z near ¯ P (cid:13)(cid:13) F where the subscript F denotesthe Frobenius norm.A more involved but even better preconditioner is the CMPthat projects ¯ Zx = u to form [18] (cid:101) ¯ Z ¯ G − ¯ Zx = (cid:101) ¯ Z ¯ G − u (25)where (cid:101) Z ij = ikη (cid:68)(cid:101) f i , ¯ G , (cid:101) f j (cid:69) S (26)and the Gram matrix G ij = (cid:68) ˆ n × f i , (cid:101) f j (cid:69) S (27)where (cid:101) f i and (cid:101) f j are chosen as Buffa-Christiansen (BC) basisfunctions. Due to the Calderon identity, (cid:101) ¯ Z ¯ G − ¯ Z is well-conditioned for a uniform mesh. When the mesh is nonuni-form, diagonal preconditioners are applied to guarantee lowcondition numbers for both ¯ G and (cid:101) ¯ Z ¯ G − ¯ Z . The reason forkeeping ¯ G well-conditioned is that one needs to compute ¯ G − u iteratively provided ¯ G is sparse. Moreover, MLFMAhas been embedded within CMP, which offers promise forsolving large scale EFIEs with good accuracy using only afew iterations [19].The IRA method is less appropriate for solving (14) due tothe lack of efficiency in preconditioning ¯ X − u . However, theGEP (6) can be solved with the JDQZ method [16]. Withoutconverting (6) to (14), the correction vector t is obtained byapproximately solving the Jacobi-Davidson correction equa-tion (cid:0) ¯ I − ¯ Ruu ∗ (cid:1) (cid:0) ¯ X − σ ¯ R (cid:1) (cid:0) ¯ I − uu ∗ ¯ R (cid:1) t = − (cid:0) ¯ X − σ ¯ R (cid:1) u , t ∗ ¯ Ru = 0 (28)where (cid:0) σ, u = ¯ Vs (cid:1) is a solution of ¯ V ∗ ¯ X ¯ Vs = σ ¯ V ∗ ¯ R ¯ Vs = σ s (29)with the search space ¯ V containing several R -orthogonal basis.The key is that even when (28) is solved with low accuracy,JDQZ still gives rise to correct solutions of (6). However,JDQZ becomes fairly slow in this case. Hence, it is alsobeneficial if efficient preconditioners are designed to acceleratethe solutions of (28). j i mm' m ixyz aa a a a Fig. 1. Point positions.
V. N
UMERICAL R ESULTS
Numerical tests are first performed to validate Equations(15) and (19). Two cases are considered as shown in Fig. 1,where the box size is a . In both cases, the source point j andsource group center m (cid:48) are located at (0 . a, . a, . and (0 . a, . a, . a ) , respectively. The observation point i and its group center m in Case 1 (dashed arrows) are locatedat (2 . a, . , . a ) and (2 . a, . a, . a ) , respectively,such that | r im + r m (cid:48) j | < | r mm (cid:48) | . The relative errors betweenleft and right hand sides of (15) and (19) are calculated withrespect to a/λ , where λ is the wavelength. The truncationnumber L is chosen to be L ≈ kd + 1 . d / ( kd ) / (30)where the number of digits of accuracy d is set to , and d = √ a in this test. As shown in Fig. 2, excellent agreementis observed except for (15) or cos( · ) at lower frequencies( a < . λ ), which can be resolved with the low-frequency fastmultipole algorithm (LF-FMA) [20]. We then place the obser-vation point i and its group center m at (1 . a, . , . and (1 . a, . a, . a ) , respectively, in Case 2 (solid arrows)such that | r im + r m (cid:48) j | > | r mm (cid:48) | . As expected, (15) is nolonger valid due to the violation of addition theorem. However,the regular part of the Green’s function, i.e., sin( · ) can becomputed with excellent accuracy [Fig. 2]. This holds evenwhen i and j are very close to each other. Although the resultscomputed by (15) and (19) are not purely real, the imaginaryparts are negligible as they are normally times smallerthan the real parts.In this study, iterations in IRA or JDQZ (GMRES orBiCGSTAB) are referred to as outer (inner) iterations. Tofurther test the validity of the proposed method, we performthe CM analysis to a rectangular PEC plate. The dimensionof the plate is . m × . m, and the frequency varies from MHz to
MHz at a step of MHz. Totally unknowns are generated in this simulation, rendering MLFMAa three-level algorithm. By solving (24) with MLFMA andIRA, four CMs whose current patterns at
MHz are plottedin Fig. 3. IRA performs outer iterations to obtain desiredmodes, where GMRES is used to compute ¯ Z − u with atolerance of − . The threshold parameters (cid:15) , (cid:15) and (cid:15) in the SAI preconditioner [17] are specified to be . , . and . , respectively.The same current patterns are also obtained by solving(6) with MLFMA and JDQZ, or by applying the generalizedSchur decomposition (also known as QZ decomposition) to − − − a/ λ R e l a t i v e E rr o r Case 1, cos( ⋅ )Case 1, sin( ⋅ )Case 2, cos( ⋅ )Case 2, sin( ⋅ ) Fig. 2. Relative error of plane wave decomposition for the real and imaginaryparts of Green’s function with respect to a/λ .(a) (b)(c) (d)Fig. 3. Current patterns of plate characteristic modes: (a) J . (b) J . (c) J . (d) J (6) with ¯ X and ¯ R generated by MoM. The magnitudes of unknowns in J computed by these approaches agree wellwith each other, as shown in Fig. 4. Taking the result obtainedby MoM + QZ as a reference, the magnitude error of the entireunknowns is . when J is computed by MLFMA + IRA,and . when computed by MLFMA + JDQZ. For the otherthree current modes, all relative errors are less than . Whenthe tolerance for GMRES in IRA is set to − , the errors arearound for all four modes.Characteristic values of J to J at different frequenciesare calculated with the above approaches, where good agree-ment is observed from Fig. 5. In JDQZ, one can obtaincorrect characteristic values even with a small number ofinner iterations (dimension of searchspace in GMRES), whichindicates that (28) is only solved with poor accuracy eachtime. Fig. 6 shows the required number of JDQZ iterationsfor certain convergence tolerance, i.e., − , with respectto the dimension of GMRES searchspace. Although a smallsearchspace reduces the computational time of each outeriteration, the tradeoff is that it increases the total number ofouter iterations.We next consider a PEC sphere with a radius of . λ . Fig. 4. Magnitude of unknowns in current mode J (a) (b)(c) (d)Fig. 5. Characteristic values at different frequencies: (a) J . (b) J . (c) J .(d) J Fig. 6. Number of JDQZ iterations with respect to dimension of GMRESsearchspace.
We increase the mesh density such that , and
15 015 unknowns are generated. A few degenerate CMs are computedwith MoM + QZ and MLFMA + IRA, respectively, wherethe corresponding characteristic values are listed in Table I.
TABLE IC
HARACTERISTIC V ALUES λ n OF A F EW C OMPUTED M ODES AT D IFFERENT M ESH D ENSITIES
Current modes MoM + QZ MLFMA + IRA
942 3519 942 3519 15015 J -0.1974 -0.2129 -0.1963 -0.2174 -0.2249 J -0.1960 -0.2131 -0.1977 -0.2170 -0.2246 J -0.1966 -0.2132 -0.1969 -0.2172 -0.2245 J J J J J
15 015 unknowns): (a) J . (b) J . (c) J . The GMRES tolerance is set to − in the IRA method.Good agreement is observed between results computed withthe two schemes. In [21] where a few CMs of a vehicle with unknowns are solved for, the memory usage is . GBand . GB for MoM + QZ and MoM + IRA, respectively.Hence, the MoM + QZ scheme has not been applied to theproblem with
15 015 unknowns which requires more memory.However, such a problem can be easily handled with the pro-posed MLFMA + IRA scheme on a small computer. Certaincharacteristic current patterns (magnitudes) are illustrated inFig. 7, which are also found from the results computed byMoM + QZ with less unknowns.We finally consider large platforms, i.e., aerospace vehicles,where we can only apply the MLFMA + IRA scheme for CManalysis. JDQZ becomes inefficient in such examples as thereare lack of preconditioners for expediting the convergence tosolutions with reasonable accuracy when (28) is solved. Wecompute several CMs of an unmanned aerial vehicle (UAV),of which the length, width and height are . m, . m,and . m, or . λ , . λ , and . λ at the frequencyof MHz, respectively. In this simulation, the number ofunknowns is
21 056 . The GMRES tolerance is set to × − ,and (cid:15) , (cid:15) and (cid:15) in SAI are set to . , . and . ,respectively. The computed current patterns (magnitudes) areillustrated in Fig. 8.The CM analysis is next performed on a Boeing 787Dreamliner. The length, width, and height of the plane is . m, . m, and . m, respectively, which are . λ , . λ , and . λ at MHz. Totally
117 834 unknowns areused in this simulation. The GMRES tolerance is again setto × − , and (cid:15) , (cid:15) and (cid:15) in SAI are set to . , . and . , respectively. Certain characteristic currents of theDreamliner are depicted in Fig 9. The color axis is adjustedto clearly manifest the current patterns on the plane body. (a)(b)Fig. 8. Characteristic modes of UAV computed using MLFMA + IRA(
21 056 unknowns): (a) Mode with characteristic value λ = − . . (b)Mode with characteristic value λ = 0 . .(a)(b)Fig. 9. Characteristic modes of Boeing 787 Dreamliner computed usingMLFMA + IRA (
117 834 unknowns): (a) Mode with characteristic value λ = 0 . . (b) Mode with characteristic value λ = − . . Numerical examples in this study are computed with afundamental, non-parallel MLFMA program and the open-source ARPACK on an Intel Core i5-2400 CPU with . GHz clock rate. Computational detail in different examples
TABLE IIC
OMPUTATIONAL D ETAIL IN D IFFERENT N UMERICAL E XAMPLES U SING
MLFMA + IRAExample nev ncv nlv nii Memory CPU timePlate 5 20 4 >
45 20 MB 15 secsSphere 10 120 5 >
170 300 MB 55 minsUAV 4 140 7 >
436 MB 2.1 hrsDreamliner 4 180 7 > ∗ Note that in the examples of UAV and Dreamliner, one may obtain converged eigenpairs with smaller ncv. In the sphere example,
15 015 unknowns are used. are summarized in Table II, where nlv is the total numberof MLFMA levels (including level- ), nev is the number ofeigenvalues needed, ncv is the maximum number of Arnoldivectors used, and nii is the number of inner iterations forcalculating each ¯ Z − u . It can be inferred that the SAI pre-conditioner becomes less efficient as the problem scale growslarge. For highly ill-conditioned impedance matrix, it calls formore efficient preconditioning schemes such as the Calderonprojector to further enhance the performance of the proposedCM analysis. VI. C ONCLUSION
We develop an MLFMA based CM analysis for largescale applications which severely challenge the conventionalMoM based approach. The standard MLFMA program can beeasily modified to perform the required matrix-vector productoperations in the iterative eigensolvers. Numerical examplesare provided to demonstrate the validity and efficiency of theproposed scheme. In even larger scale simulations, MLFMAcan be implemented in parallel on many processors for betterperformance. The same idea can be easily extended to anincorporation of the mixed-form fast multipole algorithm (MF-FMA) [22] and CM analysis, which is beneficial for modelswith fine structures. A
CKNOWLEDGMENT
The authors would like to thank Q. S. Liu for her helpfuldiscussion. R
EFERENCES [1] R. J. Garbacz and R. H. Turpin, “A generalized expansionfor radiated and scattered fields,”
IEEE Trans. AntennasPropag. , vol. 19, pp. 348–358, 1971.[2] R. F. Harrington and J. R. Mautz, “Theory of characteris-tic modes for conducting bodies,”
IEEE Trans. AntennasPropag. , vol. 19, pp. 622–628, 1971.[3] ——, “Computation of characteristic modes for conduct-ing bodies,”
IEEE Trans. Antennas Propag. , vol. 19, pp.629–639, 1971.[4] M. Cabedo-Fabres,
Systematic design of antennas usingthe theory of characteristic modes [Ph.D. thesis] . Univ.Politecnica de Valencia, Spain, 2007.[5] J. Ethier,
Antenna shape synthesis using characteristicmode concepts [Ph.D. thesis] . U. of Ottawa, Ottawa,ON, Canada, 2012.[6] K. A. Obeidat, B. D. Raines, R. G. Rojas, and B. T.Strojny, “Design of frequency reconfigurable antennas using the theory of network characteristic modes,”
IEEETrans. Antennas Propag. , vol. 58, pp. 3106–3113, 2010.[7] J. J. Adams and J. T. Bernhard, “A modal approachto tuning and bandwidth enhancement of an electricallysmall antenna,”
IEEE Trans. Antennas Propag. , vol. 59,pp. 1085–1092, 2011.[8] L. Greengard and V. Rokhlin, “A fast algorithm forparticle simulations,”
J. Comput. Phys , vol. 73, pp. 325–348, 1987.[9] R. Coifman, V. Rokhlin, and S. Wandzura, “The fastmultipole method for the wave equation: A pedestrianprescription,”
IEEE Antennas Propagat. Mag. , vol. 35,pp. 7–12, 1993.[10] J. M. Song and W. C. Chew, “Fast multipole methodsolution using parametric geometry,”
Microwave Opt.Technol. Lett. , vol. 7, pp. 760–765, 1994.[11] J. M. Song, C. C. Lu, and W. C. Chew, “Multilevelfast multipole algorithm for electromagnetic scattering bylarge complex objects,”
IEEE Trans. Antennas Propag. ,vol. 45, pp. 1488–1493, 1997.[12] W. C. Chew, J. M. Jin, E. Michielssen, and J. Song,
Fastand Efficient Algorithms in Computational Electromag-netics . Artech House, 2000.[13] R. L. Wagner and W. C. Chew, “A ray-propagation fastmultipole algorithms,”
Microwave Opt. Tech. Lett. , vol. 7,pp. 435–438, 1994.[14] M. A. Epton and B. Dembart, “Multipole translationtheory for the threedimensional Laplace and Helmholtzequations,”
SIAM J. Sci. Comput , vol. 16, pp. 865–897,1995.[15] R. Lehoucq and D. Sorensen, “Deflation techniques foran implicitly restarted Arnoldi method,”
SIAM J. MatrixAnal. Appl. , vol. 17, no. 4, pp. 789 –821, 1996.[16] D. R. Fokkema, G. L. G. Sleijpen, and van der Vorst.H. A., “Jacobi-Davidson style QR and QZ algorithms forthe reduction of matrix pencils,”
SIAM Journal ScientificComputing , vol. 20(1), pp. 94–125, 1998.[17] J. Lee, C. Zhang, and C. C. Lu, “Sparse inverse precondi-tioning of multilevel fast multipole algorithm for hybridintegral equations in electromagnetics,”
IEEE Trans. An-tennas Propag. , vol. 52, pp. 2277–2287, 2004.[18] F. P. Andriulli, K. Cools, B. H., F. Olyslager, A. Buffa,S. Christiansen, and E. Michielssen, “A multiplicativeCalderon preconditioner for the electric field integralequation,”
IEEE Trans. Antennas Propag. , vol. 56, pp.2398–2412, 2008.[19] J. Peeters, K. Cools, I. Bogaert, F. Olyslager, andD. De Zutter, “Embedding Calderon multiplicative pre-conditioners in multilevel fast multipole algorithms,”
IEEE Trans. Antennas Propag. , vol. 58, pp. 1236–1250,2010.[20] J. S. Zhao and W. C. Chew, “Three dimensional multi-level fast multipole algorithm from static to electrody-namic,”
Microwave Opt. Tech. Lett. , vol. 26, pp. 43–48,2000.[21] D. J. Ludick, E. Lezar, and U. Jakobus, “Characteristicmode analysis of arbitrary electromagnetic structures us-ing feko,”
Conference on Electromagnetics in Advanced
Applications (ICEAA) , pp. 208–211, 2012.[22] L. J. Jiang and W. C. Chew, “A mixed-form fast multipolealgorithm,”