電磁散射問題的快速計(jì)算

電磁散射問題的快速計(jì)算1電磁散射問題計(jì)算復(fù)雜電大尺寸導(dǎo)體的雷達(dá)散射截面(RCS)例如:飛機(jī)(VFY218)三角形網(wǎng)格剖分積分方程離散迭代法求解預(yù)條件子多層組劃分MLFMM并行化后處理

2計(jì)算原理Maxwell方程組電磁場(chǎng)積分方程離散后的系數(shù)矩陣格林函數(shù)由多極子展開得到矩陣遠(yuǎn)場(chǎng)部分的近似通過多層組劃分和層間插值，迭代法求解過程中的矩陣-向量乘積的復(fù)雜度降低到O(NlogN)3計(jì)算流程網(wǎng)格剖分與幾何信息讀入最細(xì)層的近場(chǎng)信息聚集(內(nèi)插值)近場(chǎng)矩陣近場(chǎng)作用轉(zhuǎn)移(次相鄰組)矩量法迭代法求解向量運(yùn)算矩陣-向量乘積BLAS,LAPACK構(gòu)造分布式八叉樹發(fā)散(外插值)各層的遠(yuǎn)場(chǎng)信息遠(yuǎn)場(chǎng)作用預(yù)條件子線性方程組電磁場(chǎng)積分方程4并行迭代方法向量運(yùn)算(BLAS-1)向量運(yùn)算的并行矩陣-向量乘積(BLAS-2)結(jié)構(gòu)矩陣對(duì)角化(FFT)稠密矩陣稀疏化(FMM,小波變換)矩陣-向量乘積的并行傳統(tǒng):矩陣分塊、區(qū)域分解MLFMM:樹結(jié)構(gòu)并行劃分提高并行效率高效預(yù)條件子(塊對(duì)角、稀疏近似逆)重排運(yùn)算次序，讓計(jì)算與通信的重疊計(jì)算任務(wù)的劃分盡可能保證負(fù)載平衡ZnearI聚集轉(zhuǎn)移發(fā)散ZnearI聚集轉(zhuǎn)移發(fā)散ZnearI聚集轉(zhuǎn)移發(fā)散[Zij][Ij]MLFMM的并行計(jì)算P0P1Pnp-1數(shù)據(jù)交換數(shù)據(jù)交換5FastAlgorithmyearmethodreferencestorageflops1947GE(banded)VonNeumann,Goldstinen5n71950OptimalSORYoungn3n4logn1971CGReidn3n3.5logn1984FullMGBrandtn3n31987FMMGreengard,Rokhlinn2lognn2lognnnIfn=64,thistableimpliesanoverallreductioninflopsof160million,whichmeetstheMoore’sLaw!(doublingin1.5year)SciDACInitiative,DOE,CSGF,2005n6計(jì)算結(jié)果對(duì)比2001年W.C.Chew等開發(fā)出計(jì)算電磁散射問題的并行程序— ScaleME，該程序在SGIOrigin2000上完成了問題規(guī)模超過

1000萬的VFY218飛機(jī)模型的RCS計(jì)算(128個(gè)CPU，內(nèi)存69G， 計(jì)算時(shí)間為7.5hour，入射波頻率為8GHz)。它顯示了MLFMM

和并行處理技術(shù)在計(jì)算電磁學(xué)中的巨大作用2001年最快的機(jī)器為ASCIWhite(LLNL)，包含8192個(gè)處理機(jī)， 測(cè)試LINPACK時(shí)，峰值為12.3TFLOPS,可求解的稠密線性 方程組規(guī)模達(dá)到518,000(優(yōu)化的高斯消去法,內(nèi)存:6.2TB)計(jì)算直徑為12個(gè)波長(zhǎng)的金屬球，未知量規(guī)模：172,800MLFMM:內(nèi)存474MB,計(jì)算時(shí)間：21分鐘（P4,2GCPU,1GMemory）GE:內(nèi)存2.4TB,計(jì)算時(shí)間：約1個(gè)月(估計(jì)值)7國(guó)外研究現(xiàn)狀1987年，耶魯大學(xué)的L.Greengard和V.Rokhlin發(fā)明了快速多極子方法(FMM)，將求解多點(diǎn)源位勢(shì)的復(fù)雜度降到O(N)；1988年，J.Carrier提出了自適應(yīng)樹結(jié)構(gòu)策略，解決了粒子不均勻分布時(shí)FMM的計(jì)算效率下降的問題；1993年，R.Coifman和V.Rokhlin用FMM求解Helmholtz型積分 方程，將求解積分方程的計(jì)算與存儲(chǔ)的復(fù)雜度降到O(N1.5)；1994年，Illinois大學(xué)計(jì)算電磁中心的W.C.Chew等提出了多層快速多極子方法(MLFMM)，也就是通過多層組劃分和層間插值，將迭代求解電磁場(chǎng)積分方程的復(fù)雜度降低到O(NlogN）。L.Greengard(1999)給出了N體問題的FMM嚴(yán)格誤差估計(jì)；S.Amini(2000)給出了聲波散射問題的FMM的誤差分析；對(duì)于電磁散射問題的FMM嚴(yán)格誤差估計(jì)和復(fù)雜度分析，M.A.Epton和B.Dembart(1995)，J.Rahola(1998)，法國(guó)巴黎大學(xué)的E.Darve(2001)和Q.Carayol，F(xiàn).

Collino(2004)給出了許多重要結(jié)果。

8國(guó)外研究現(xiàn)狀V.Rokhlin（1993）給出了轉(zhuǎn)移算子對(duì)角化的概念；M.A.Epton（1995）綜述了FMM在Laplace方程和HelmHoltz方程中的相關(guān)理論；L.Greengard和V.Rokhlin(1997)提出了新版本的快速多極子方法，并用于求解3維Laplace方程，它增加了指數(shù)展開的概念，并將其用于低頻Helmholtz方程(1998)，克服了傳統(tǒng)FMM方法的不穩(wěn)定性；X.B.Sun(2001)給出了矩陣形式描述的FMM，代替了傳統(tǒng)FMM求解Laplace方程的級(jí)數(shù)和形式。E.Darve和H.Pascal（2004）給出了基于指數(shù)展開，同時(shí)在低頻和高頻下 穩(wěn)定的新版本的FMM在Maxwell方程組求解的詳細(xì)數(shù)學(xué)描述和數(shù)值實(shí)現(xiàn)。類似的方法W.C.Chew也曾提到(2005)。M.Nilsson（2002）提出了適合頻域Maxwell方程的塊QMR迭代FMM方法;法國(guó)CERFACS的G.Alleon等2004年給出了基于MLFMM的有效雙層迭代GMRES(m)算法和稀疏近似逆預(yù)處理的數(shù)值實(shí)現(xiàn)。9國(guó)外研究現(xiàn)狀國(guó)外已開發(fā)的相關(guān)軟件DPMTA(aDistributedParallelMultipoleTreeAlgorithm libraryforN-bodyProblem)是杜克大學(xué)電子工程系的

W.Rankin等于2002年完成的.它提供了FMM并行求解大規(guī)模粒子系統(tǒng)相互作用的N體問題的軟件包EGO

是哥廷根大學(xué)普朗克研究所的Hulmut等開發(fā)的分子動(dòng) 力學(xué)模擬的并行軟件,它采用FMM計(jì)算分子之間的相互作用,

并采用時(shí)間多步法計(jì)算分子的運(yùn)動(dòng)(這種算法稱之為快速多步 結(jié)構(gòu)自適應(yīng)多極子方法,),

可計(jì)算數(shù)萬個(gè)原子的運(yùn)動(dòng)軌道

Q-CHEM

重頭算量子化學(xué)軟件,采用FMM計(jì)算電荷庫(kù)侖作用FMMP(FMMforParticle)由日本的ShujiOgata開發(fā),采用

FMM并行求解三維庫(kù)倫場(chǎng)位勢(shì)、力、勢(shì)能和微觀應(yīng)力張量10國(guó)外研究現(xiàn)狀PEC2D(PerfectlyElectricalConductingin2D)計(jì)算2維理 想導(dǎo)體的雙站散射截面.對(duì)邊界曲線分段,每段中點(diǎn)取為未 知點(diǎn),采用點(diǎn)配法和脈沖基函數(shù)離散電場(chǎng)積分方程(EFIE),

并利用GMRES迭代求解線性方程組,每次迭代采用FMM

加速.(SanjayVelamparambil,Illinois,2003)。ANSYS/FEKO是一款用于3D結(jié)構(gòu)電磁場(chǎng)分析的仿真工具， 基于矩量法對(duì)Maxwell方程組的求解。對(duì)電大散射問題，

FEKO引入了MLFMM，而且有并行版本，它是第一個(gè) 商業(yè)化的MLFMM軟件，2004年ANSYS8.0開始包含F(xiàn)EKO。KIFMA

是加州理工的L.Ying

開發(fā)的無約束積分核FMM算法 的數(shù)值軟件,積分核可以是指數(shù)型函數(shù)或高斯型函數(shù),也可以 是RBF(光線基函數(shù))，在計(jì)算插值和轉(zhuǎn)移時(shí)調(diào)用了FFTW。11國(guó)內(nèi)研究現(xiàn)狀國(guó)防科技大學(xué)的韓明華、彭宇行、李思昆等在973項(xiàng)目(分布式 虛擬環(huán)境技術(shù))支持下，2004年完成了‘基于Linux集群電磁散射 并行計(jì)算實(shí)現(xiàn)’的課題研究；電子科技大學(xué)的盧光輝、孫世新、聶在平等在國(guó)家自然科學(xué)基 金支持下，于2005年在自然完成了面向工程應(yīng)用的復(fù)雜目標(biāo)電 磁散射高效數(shù)值分析軟件A-UEST的開發(fā)；北理工的盛新慶教授一直致力于混合有限元、邊界元和快速多 極子方法（合元極）的研究，相關(guān)工作見2004年出版的《計(jì)算 電磁學(xué)要論》一書；清華大學(xué)工程力學(xué)系的姚振漢等一直致力于將FMM用于彈性力 學(xué)邊界元的研究，比如復(fù)合材料（比如纖維或顆粒夾雜問題） 的大規(guī)模數(shù)值模擬；大連理工大學(xué)的滕斌等做過水流繞射的數(shù)值計(jì)算，將多極子展 開用于求解無源不可壓流的高階邊界元；12電磁場(chǎng)積分方程EFIEMFIECFIEGreen函數(shù)13矩量法（MOM）RWG矢量基函數(shù)MOM離散14

球面的三角網(wǎng)格剖分RWG矢量基函數(shù)球面導(dǎo)體存在解析解，可驗(yàn)證算法和程序的正確性球面的網(wǎng)格剖分相對(duì)簡(jiǎn)單15奇異積分?jǐn)?shù)值積分主值積分項(xiàng)奇異點(diǎn)積分項(xiàng)化為線積分:Gauss公式、面梯度(散度)公式16迭代法Krylov子空間迭代法如:CG,BiCGStab TFQMR,GMRES主要的計(jì)算和存儲(chǔ)來自矩陣-向量的乘積系數(shù)矩陣為復(fù)稠密矩陣FMM加速策略不顯式計(jì)算遠(yuǎn)場(chǎng)矩陣近似計(jì)算遠(yuǎn)場(chǎng)作用預(yù)條件子對(duì)角（塊對(duì)角）部分稀疏近似逆（SPAI）17

FMM的思想轉(zhuǎn)換算子:M2M,M2L,L2L

聚集–轉(zhuǎn)移–發(fā)散

M:多極子展開L:局部展開

18級(jí)數(shù)展開球面波加法定理Gauss型數(shù)值積分平面波展開19多極子展開轉(zhuǎn)移項(xiàng)矢量加法原理FMM形式的Green函數(shù)根據(jù)矢量恒等式可以得到20FMMaccelerateFMM形式的矩陣向量乘積近場(chǎng)部分遠(yuǎn)場(chǎng)部分

Gn:neargroupComplexity=N(C1M+C2N/M)Optimal:O(N1.5)

21

MLFMM的思想空間多層組劃分Morton編號(hào)相鄰組的作用遠(yuǎn)場(chǎng)組的上聚次相鄰組中心 的轉(zhuǎn)移遠(yuǎn)場(chǎng)組的下推22積分區(qū)域的多層組劃分23樹結(jié)構(gòu)代碼上聚A—M2M內(nèi)插值下推C—M2LB—L2L外插值二叉樹的例子24Morton次序多層組劃分、編號(hào)、構(gòu)建樹二維計(jì)算區(qū)域?qū)?yīng)的分布式四叉樹25

樹結(jié)構(gòu)的并行劃分MLFMM并行實(shí)現(xiàn)分布式八叉樹負(fù)載平衡相互作用表列相鄰結(jié)點(diǎn)的通信次相鄰點(diǎn)的通信MPI_Alltoall26MLFMMUpwardPassDownwardPassSummationInterpolationAnterpolationInteractionlistComplexity:O(NlogN)27MLFMM:aninsight近場(chǎng)矩陣-向量乘(相鄰組的作用):直接計(jì)算通過轉(zhuǎn)移計(jì)算次相鄰組的作用通過插值實(shí)現(xiàn)上聚和下推(比次相鄰更遠(yuǎn)的作用)矩陣-向量乘積(源點(diǎn)對(duì)場(chǎng)點(diǎn)的作用)的復(fù)雜度:O(NlogN)28多層快速多極子方法的兩層示意圖29層間插值Lagrange多項(xiàng)式插值球諧變換的全局插值插值系數(shù)的計(jì)算可分裂為如下兩步：采用FFT，復(fù)雜度降低為O(qlogq)30單位球面數(shù)值積分取樣點(diǎn)插值滿足Gauss-Legendre數(shù)值積分,

31

迭代收斂曲線TFQMR迭代求解EFIE,MFIE,CFIE(球面,3600個(gè)未知量)32

雷達(dá)散射截面(RCS)單位球面(未知量個(gè)數(shù)3600,入射頻率300MHz)33MLFMM-TFQMR,EFIE,球面問題規(guī)模計(jì)算近場(chǎng)上聚轉(zhuǎn)移下推其它17424123.359611.493832.665716.54336.218541616340.558453.8093189.494579.710020.3578MLFMM-TFQMR,CFIE(alpha=0.2),球面未知量個(gè)數(shù)129623043600518470569216入射頻率(MHz)190250300360420480內(nèi)存大小(MB)141925324050CPU計(jì)算時(shí)間(s)44.10167.23176.33316.87679.891842.98迭代次數(shù)141249164220315624問題規(guī)模計(jì)算近場(chǎng)上聚轉(zhuǎn)移下推其它總時(shí)間129621.68385.79686.82658.87380.919244.109216498.1216240.1441711.0337351.878341.80231842.9834每次迭代的計(jì)算時(shí)間MLFMM-TFQMR，CFIE35迭代所需內(nèi)存MLFMM–TFQMR,CFIE36BeforeinstallNETGENinDeepComp6800,youwillneed:OpenGl,Tcl/Tk,TIX.ItcanbeusedinFEM/FDM/BEM,andsolverslikeUG37RCSofaunitcube(perfectelectricallyconductor),MLFMM

6146points,12288patches,18432edges;900MHzincidentwave

Numberofnodesatlevel0、1、2、3、4oftheoctreeare：

1(root)，7，26，146，596(leafboxes)

NumberofiterationsofTFQMRare:

EFIE,278(835mults,1838s);MFIE,35;CFIE,1438結(jié)果分析結(jié)論基本上能體現(xiàn)出MLFMM的O(NlogN)復(fù)雜度能保證較高的計(jì)算精度（基于嚴(yán)格誤差估計(jì)）問題每組的相鄰組有33-1=26個(gè)，次相鄰組卻有63-33=189個(gè)，而且在最細(xì)層計(jì)算復(fù)雜度為O(MLN),近似于O(N1.5),這導(dǎo)致轉(zhuǎn)移項(xiàng)的計(jì)算時(shí)間太長(zhǎng).于是Rokhlin、Velamparambil、Alpert等提出了一些新的策略(如預(yù)先保存轉(zhuǎn)移模式，采用快速Legendre展開）低頻時(shí)不穩(wěn)定(Greengard

提出了指數(shù)展開)收斂可能很慢，比如EFIE迭代次數(shù)遠(yuǎn)多于MFIE和CFIE.

如何構(gòu)造出好的預(yù)條件子改進(jìn)數(shù)值性能，使其收斂更快？39

今后的研究方向軟件上：（1）MLFMM的并行實(shí)現(xiàn)有許多技術(shù)細(xì)節(jié)要解決；（2）高效快速的并行迭代方法和預(yù)條件子；（3）計(jì)算性能的優(yōu)化，向量運(yùn)算調(diào)用MKL的BLAS和Lapack， 采用內(nèi)存循環(huán)優(yōu)化存儲(chǔ)，采用動(dòng)靜態(tài)調(diào)度提高并行效率；（4）自適應(yīng)技術(shù)與FMM的耦合(自適應(yīng)八叉樹)理論上：（1）級(jí)數(shù)展開的理論誤差估計(jì)、指數(shù)展開等；（2）算法的優(yōu)化，如轉(zhuǎn)移項(xiàng)的修正內(nèi)插值(快速Legendre展開、 修正Lagrange插值)或全局頻譜插值、圖排序的并行劃分；（3）多極子在其它領(lǐng)域的應(yīng)用。40誤差估計(jì)多極子展開的截?cái)嗾`差層間插值的誤差數(shù)值積分的誤差迭代終止誤差前后處理的誤差Rokhlin的一個(gè)經(jīng)驗(yàn)公式:L=kd+cln(kd+π)41指數(shù)展開M2L增加為M2exp,exp2exp,exp2LExponentialexpansionsPropagativetermEvanescentterm其中M2e(將多極子開開轉(zhuǎn)化為平面波展開)；e2e:(平面波展開的轉(zhuǎn)移)；e2L:(將平面波展開轉(zhuǎn)化為局部展開)坐標(biāo)旋轉(zhuǎn)Gauss-Laguerre

Quadrature42轉(zhuǎn)移項(xiàng)快速計(jì)算轉(zhuǎn)移模式FMM-1D43分子動(dòng)力學(xué)和恒星動(dòng)力學(xué)靜電場(chǎng)和引力場(chǎng)位勢(shì)的多極子展開矩陣向量乘積形式44彈性力學(xué)Navier方程Kelvin基本解(3D)邊界積分方程FM-BEM——(Yamada,Yoshida,Nishimura,etc.)45流體力學(xué)

Stokes方程N(yùn)umman

邊界條件單層位勢(shì)雙層位勢(shì)速度場(chǎng)邊界積分方程46聲 學(xué)聲波滿足的Helmholtz方程

Robin邊界條件,Sommerfield邊界條件（輻射條件）邊界積分方程格林函數(shù)對(duì)Green函數(shù)的級(jí)數(shù)展開方式與電磁波相同,FMM可將邊界元方法（BEM）得到的線性方程組遠(yuǎn)場(chǎng)部分稀疏化，從而將積分方程數(shù)值求解的復(fù)雜度降低到O(NlogN)，因此FMM可以作為BEM的一種加速技術(shù)(FM-BEM)

。47量子力學(xué)Shroedinger方程對(duì)于定態(tài)問題它是Helmholtz型的,其中V為勢(shì)能,E為能量.

氫原子中的庫(kù)倫勢(shì):

湯川(Yukawa)勢(shì):積分形式散射問題Klein-Gordan

方程48圖形處理(快速插值)光線基函數(shù)(RBF)

級(jí)數(shù)展開(如Laurent展開,Taylor展開……要求級(jí)數(shù)收斂)

插值函數(shù)線性方程組迭代法求解上面的線性方程組時(shí)，可以通過級(jí)數(shù)展開減小矩陣向量乘積的復(fù)雜度P—多項(xiàng)式插值(如樣條插值)；a—

待定插值系數(shù)49快速高斯變換(FGT)Gauss型求和Gauss函數(shù)Hermite

展開

(在遠(yuǎn)場(chǎng)收斂)Hermite

函數(shù)熱方程初邊值問題 溫度場(chǎng)對(duì)流-反應(yīng)-擴(kuò)散方程核密度估計(jì)(數(shù)理統(tǒng)計(jì)中的Gauss型分布)50

參考資料盛新慶,計(jì)算電磁學(xué)要論,科學(xué)出版社,2004.王長(zhǎng)清,現(xiàn)代計(jì)算電磁學(xué)基礎(chǔ),北京大學(xué)出版社,2005.梁昆淼,數(shù)學(xué)物理方法,高等教育出版社,1995.胡俊,聶在平等,三維電大目標(biāo)散射求解的多層快速多極子 方法,電波科學(xué)學(xué)報(bào),2004,19-5.孫世新,盧光輝等,并行算法及其應(yīng)用,機(jī)械工業(yè)出版社,2005.盧光輝,孫世新,并行處理技術(shù)在電大尺寸復(fù)雜目標(biāo)電磁散射中的應(yīng)用.電子學(xué)報(bào),2003,31(6).譚云華,周樂柱,三維電大尺寸復(fù)雜群目標(biāo)的單站RCS的快速多極子分析,北京大學(xué)學(xué)報(bào)(自然科學(xué)版),2004,40:5.余德浩,自然邊界元方法的數(shù)學(xué)理論,科學(xué)出版社,北京,1993.

楊德全,趙忠生,邊界元理論及應(yīng)用,北京理工大學(xué)出版社,2002

51ReferenceL.GreengardandV.Rokhlin,“AFastAlgorithmForParticleSimulations",J.Comput.Phys.,vol.73,pp.325-348,1987L.GreengardandV.Rokhlin.AnewversionofthefastmultipolemethodfortheLaplaceequationinthreedimensions.Acta

Numerica6,pages229–269,1997.L.Greengard.TheRapidEvaluationOfPotentialFieldsinParticleSystems.MITPress,Cambridge,Massachusettes,1998.V.Rokhlin,Rapidsolutionofintegralequationsofscatteringtheoryintwodimensions,JournalofComp.Phys.,vol.86,414-439,1990.V.Rokhlin,R.Coifman,S.Wandzura,Thefastmultipolemethodforthewaveequation:apedestrianprescription,IEEEAntennasandPropagationMagazine,vol.35,pp.7-12,1993.52ReferenceW.C.ChewandJ.M.Song.Multilevelfastmultipolealgorithmforsolvingcombinedfieldintegralequationofelectromagneticscattering.MicrowaveandOptical