電磁散射問題的快速計算

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

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

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

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

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

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

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

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

并采用時間多步法計算分子的運動(這種算法稱之為快速多步 結構自適應多極子方法,),

可計算數(shù)萬個原子的運動軌道

Q-CHEM

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

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

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

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

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

是加州理工的L.Ying

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

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

FMM的思想轉換算子:M2M,M2L,L2L

聚集–轉移–發(fā)散

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

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

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

21

MLFMM的思想空間多層組劃分Morton編號相鄰組的作用遠場組的上聚次相鄰組中心 的轉移遠場組的下推22積分區(qū)域的多層組劃分23樹結構代碼上聚A—M2M內插值下推C—M2LB—L2L外插值二叉樹的例子24Morton次序多層組劃分、編號、構建樹二維計算區(qū)域對應的分布式四叉樹25

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

31

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

雷達散射截面(RCS)單位球面(未知量個數(shù)3600,入射頻率300MHz)33MLFMM-TFQMR,EFIE,球面問題規(guī)模計算近場上聚轉移下推其它17424123.359611.493832.665716.54336.218541616340.558453.8093189.494579.710020.3578MLFMM-TFQMR,CFIE(alpha=0.2),球面未知量個數(shù)129623043600518470569216入射頻率(MHz)190250300360420480內存大小(MB)141925324050CPU計算時間(s)44.10167.23176.33316.87679.891842.98迭代次數(shù)141249164220315624問題規(guī)模計算近場上聚轉移下推其它總時間129621.68385.79686.82658.87380.919244.109216498.1216240.1441711.0337351.878341.80231842.9834每次迭代的計算時間MLFMM-TFQMR，CFIE35迭代所需內存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結果分析結論基本上能體現(xiàn)出MLFMM的O(NlogN)復雜度能保證較高的計算精度（基于嚴格誤差估計）問題每組的相鄰組有33-1=26個，次相鄰組卻有63-33=189個，而且在最細層計算復雜度為O(MLN),近似于O(N1.5),這導致轉移項的計算時間太長.于是Rokhlin、Velamparambil、Alpert等提出了一些新的策略(如預先保存轉移模式，采用快速Legendre展開）低頻時不穩(wěn)定(Greengard

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

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

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

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

Stokes方程Numman

邊界條件單層位勢雙層位勢速度場邊界積分方程46聲 學聲波滿足的Helmholtz方程

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

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

氫原子中的庫倫勢:

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

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

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

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

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

展開

(在遠場收斂)Hermite

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

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

楊德全,趙忠生,邊界元理論及應用,北京理工大學出版社,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