电磁散射问题的快速计算

电磁散射问题的快速计算1电磁散射问题计算复杂电大尺寸导体的雷达散射截面(RCS)例如:飞机(VFY218)三角形网格剖分积分方程离散迭代法求解预条件子多层组划分MLFMM并行化后处理

2计算原理Maxwell方程组电磁场积分方程离散后的系数矩阵格林函数由多极子展开得到矩阵远场部分的近似通过多层组划分和层间插值，迭代法求解过程中的矩阵-向量乘积的复杂度降低到O(NlogN)3计算流程网格剖分与几何信息读入最细层的近场信息聚集(内插值)近场矩阵近场作用转移(次相邻组)矩量法迭代法求解向量运算矩阵-向量乘积BLAS,LAPACK构造分布式八叉树发散(外插值)各层的远场信息远场作用预条件子线性方程组电磁场积分方程4并行迭代方法向量运算(BLAS-1)向量运算的并行矩阵-向量乘积(BLAS-2)结构矩阵对角化(FFT)稠密矩阵稀疏化(FMM,小波变换)矩阵-向量乘积的并行传统:矩阵分块、区域分解MLFMM:树结构并行划分提高并行效率高效预条件子(块对角、稀疏近似逆)重排运算次序，让计算与通信的重叠计算任务的划分尽可能保证负载平衡ZnearI聚集转移发散ZnearI聚集转移发散ZnearI聚集转移发散[Zij][Ij]MLFMM的并行计算P0P1Pnp-1数据交换数据交换5FastAlgorithmyearmethodreferencestorageflops1947GE(banded)VonNeumann,Goldstinen5n71950OptimalSORYoungn3n4logn1971CGReidn3n3.5logn1984FullMGBrandtn3n31987FMMGreengard,Rokhlinn2lognn2lognnnIfn=64,thistableimpliesanoverallreductioninflopsof160million,whichmeetstheMoore’sLaw!(doublingin1.5year)SciDACInitiative,DOE,CSGF,2005n6计算结果对比2001年W.C.Chew等开发出计算电磁散射问题的并行程序— ScaleME，该程序在SGIOrigin2000上完成了问题规模超过

1000万的VFY218飞机模型的RCS计算(128个CPU，内存69G， 计算时间为7.5hour，入射波频率为8GHz)。它显示了MLFMM

和并行处理技术在计算电磁学中的巨大作用2001年最快的机器为ASCIWhite(LLNL)，包含8192个处理机， 测试LINPACK时，峰值为12.3TFLOPS,可求解的稠密线性 方程组规模达到518,000(优化的高斯消去法,内存:6.2TB)计算直径为12个波长的金属球，未知量规模：172,800MLFMM:内存474MB,计算时间：21分钟（P4,2GCPU,1GMemory）GE:内存2.4TB,计算时间：约1个月(估计值)7国外研究现状1987年，耶鲁大学的L.Greengard和V.Rokhlin发明了快速多极子方法(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.

Collino(2004)给出了许多重要结果。

8国外研究现状V.Rokhlin（1993）给出了转移算子对角化的概念；M.A.Epton（1995）综述了FMM在Laplace方程和HelmHoltz方程中的相关理论；L.Greengard和V.Rokhlin(1997)提出了新版本的快速多极子方法，并用于求解3维Laplace方程，它增加了指数展开的概念，并将其用于低频Helmholtz方程(1998)，克服了传统FMM方法的不稳定性；X.B.Sun(2001)给出了矩阵形式描述的FMM，代替了传统FMM求解Laplace方程的级数和形式。E.Darve和H.Pascal（2004）给出了基于指数展开，同时在低频和高频下 稳定的新版本的FMM在Maxwell方程组求解的详细数学描述和数值实现。类似的方法W.C.Chew也曾提到(2005)。M.Nilsson（2002）提出了适合频域Maxwell方程的块QMR迭代FMM方法;法国CERFACS的G.Alleon等2004年给出了基于MLFMM的有效双层迭代GMRES(m)算法和稀疏近似逆预处理的数值实现。9国外研究现状国外已开发的相关软件DPMTA(aDistributedParallelMultipoleTreeAlgorithm libraryforN-bodyProblem)是杜克大学电子工程系的

W.Rankin等于2002年完成的.它提供了FMM并行求解大规模粒子系统相互作用的N体问题的软件包EGO

是哥廷根大学普朗克研究所的Hulmut等开发的分子动 力学模拟的并行软件,它采用FMM计算分子之间的相互作用,

并采用时间多步法计算分子的运动(这种算法称之为快速多步 结构自适应多极子方法,),

可计算数万个原子的运动轨道

Q-CHEM

重头算量子化学软件,采用FMM计算电荷库仑作用FMMP(FMMforParticle)由日本的ShujiOgata开发,采用

FMM并行求解三维库伦场位势、力、势能和微观应力张量10国外研究现状PEC2D(PerfectlyElectricalConductingin2D)计算2维理 想导体的双站散射截面.对边界曲线分段,每段中点取为未 知点,采用点配法和脉冲基函数离散电场积分方程(EFIE),

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

加速.(SanjayVelamparambil,Illinois,2003)。ANSYS/FEKO是一款用于3D结构电磁场分析的仿真工具， 基于矩量法对Maxwell方程组的求解。对电大散射问题，

FEKO引入了MLFMM，而且有并行版本，它是第一个 商业化的MLFMM软件，2004年ANSYS8.0开始包含FEKO。KIFMA

是加州理工的L.Ying

开发的无约束积分核FMM算法 的数值软件,积分核可以是指数型函数或高斯型函数,也可以 是RBF(光线基函数)，在计算插值和转移时调用了FFTW。11国内研究现状国防科技大学的韩明华、彭宇行、李思昆等在973项目(分布式 虚拟环境技术)支持下，2004年完成了‘基于Linux集群电磁散射 并行计算实现’的课题研究；电子科技大学的卢光辉、孙世新、聂在平等在国家自然科学基 金支持下，于2005年在自然完成了面向工程应用的复杂目标电 磁散射高效数值分析软件A-UEST的开发；北理工的盛新庆教授一直致力于混合有限元、边界元和快速多 极子方法（合元极）的研究，相关工作见2004年出版的《计算 电磁学要论》一书；清华大学工程力学系的姚振汉等一直致力于将FMM用于弹性力 学边界元的研究，比如复合材料（比如纤维或颗粒夹杂问题） 的大规模数值模拟；大连理工大学的滕斌等做过水流绕射的数值计算，将多极子展 开用于求解无源不可压流的高阶边界元；12电磁场积分方程EFIEMFIECFIEGreen函数13矩量法（MOM）RWG矢量基函数MOM离散14

球面的三角网格剖分RWG矢量基函数球面导体存在解析解，可验证算法和程序的正确性球面的网格剖分相对简单15奇异积分数值积分主值积分项奇异点积分项化为线积分:Gauss公式、面梯度(散度)公式16迭代法Krylov子空间迭代法如:CG,BiCGStab TFQMR,GMRES主要的计算和存储来自矩阵-向量的乘积系数矩阵为复稠密矩阵FMM加速策略不显式计算远场矩阵近似计算远场作用预条件子对角（块对角）部分稀疏近似逆（SPAI）17

FMM的思想转换算子:M2M,M2L,L2L

聚集–转移–发散

M:多极子展开L:局部展开

18级数展开球面波加法定理Gauss型数值积分平面波展开19多极子展开转移项矢量加法原理FMM形式的Green函数根据矢量恒等式可以得到20FMMaccelerateFMM形式的矩阵向量乘积近场部分远场部分

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

21

MLFMM的思想空间多层组划分Morton编号相邻组的作用远场组的上聚次相邻组中心 的转移远场组的下推22积分区域的多层组划分23树结构代码上聚A—M2M内插值下推C—M2LB—L2L外插值二叉树的例子24Morton次序多层组划分、编号、构建树二维计算区域对应的分布式四叉树25

树结构的并行划分MLFMM并行实现分布式八叉树负载平衡相互作用表列相邻结点的通信次相邻点的通信MPI_Alltoall26MLFMMUpwardPassDownwardPassSummationInterpolationAnterpolationInteractionlistComplexity:O(NlogN)27MLFMM:aninsight近场矩阵-向量乘(相邻组的作用):直接计算通过转移计算次相邻组的作用通过插值实现上聚和下推(比次相邻更远的作用)矩阵-向量乘积(源点对场点的作用)的复杂度:O(NlogN)28多层快速多极子方法的两层示意图29层间插值Lagrange多项式插值球谐变换的全局插值插值系数的计算可分裂为如下两步：采用FFT，复杂度降低为O(qlogq)30单位球面数值积分取样点插值满足Gauss-Legendre数值积分,

31

迭代收敛曲线TFQMR迭代求解EFIE,MFIE,CFIE(球面,3600个未知量)32

雷达散射截面(RCS)单位球面(未知量个数3600,入射频率300MHz)33MLFMM-TFQMR,EFIE,球面问题规模计算近场上聚转移下推其它17424123.359611.493832.665716.54336.218541616340.558453.8093189.494579.710020.3578MLFMM-TFQMR,CFIE(alpha=0.2),球面未知量个数129623043600518470569216入射频率(MHz)190250300360420480内存大小(MB)141925324050CPU计算时间(s)44.10167.23176.33316.87679.891842.98迭代次数141249164220315624问题规模计算近场上聚转移下推其它总时间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结果分析结论基本上能体现出MLFMM的O(NlogN)复杂度能保证较高的计算精度（基于严格误差估计）问题每组的相邻组有33-1=26个，次相邻组却有63-33=189个，而且在最细层计算复杂度为O(MLN),近似于O(N1.5),这导致转移项的计算时间太长.于是Rokhlin、Velamparambil、Alpert等提出了一些新的策略(如预先保存转移模式，采用快速Legendre展开）低频时不稳定(Greengard

提出了指数展开)收敛可能很慢，比如EFIE迭代次数远多于MFIE和CFIE.

如何构造出好的预条件子改进数值性能，使其收敛更快？39

今后的研究方向软件上：（1）MLFMM的并行实现有许多技术细节要解决；（2）高效快速的并行迭代方法和预条件子；（3）计算性能的优化，向量运算调用MKL的BLAS和Lapack， 采用内存循环优化存储，采用动静态调度提高并行效率；（4）自适应技术与FMM的耦合(自适应八叉树)理论上：（1）级数展开的理论误差估计、指数展开等；（2）算法的优化，如转移项的修正内插值(快速Legendre展开、 修正Lagrange插值)或全局频谱插值、图排序的并行划分；（3）多极子在其它领域的应用。40误差估计多极子展开的截断误差层间插值的误差数值积分的误差迭代终止误差前后处理的误差Rokhlin的一个经验公式:L=kd+cln(kd+π)41指数展开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边界条件（辐射条件）边界积分方程格林函数对Green函数的级数展开方式与电磁波相同,FMM可将边界元方法（BEM）得到的线性方程组远场部分稀疏化，从而将积分方程数值求解的复杂度降低到O(NlogN)，因此FMM可以作为BEM的一种加速技术(FM-BEM)

。47量子力学Shroedinger方程对于定态问题它是Helmholtz型的,其中V为势能,E为能量.

氢原子中的库伦势:

汤川(Yukawa)势:积分形式散射问题Klein-Gordan

方程48图形处理(快速插值)光线基函数(RBF)

级数展开(如Laurent展开,Taylor展开……要求级数收敛)

插值函数线性方程组迭代法求解上面的线性方程组时，可以通过级数展开减小矩阵向量乘积的复杂度P—多项式插值(如样条插值)；a—

待定插值系数49快速高斯变换(FGT)Gauss型求和Gauss函数Hermite

展开

(在远场收敛)Hermite

函数热方程初边值问题 温度场对流-反应-扩散方程核密度估计(数理统计中的Gauss型分布)50

参考资料盛新庆,计算电磁学要论,科学出版社,2004.王长清,现代计算电磁学基础,北京大学出版社,2005.梁昆淼,数学物理方法,高等教育出版社,1995.胡俊,聂在平等,三维电大目标散射求解的多层快速多极子 方法,电波科学学报,2004,19-5.孙世新,卢光辉等,并行算法及其应用,机械工业出版社,2005.卢光辉,孙世新,并行处理技术在电大尺寸复杂目标电磁散射中的应用.电子学报,2003,31(6).谭云华,周乐柱,三维电大尺寸复杂群目标的单站RCS的快速多极子分析,北京大学学报(自然科学版),2004,40:5.余德浩,自然边界元方法的数学理论,科学出版社,北京,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