快速多极子方法的并行技术

1纲要•FMM•DataStructures•Parallelization2纲要•FMM•DataStructures•Parallelization3FMMinComputationalElectromagnetics•EFIE•MFIE•CFIE•Green函数2()(),,()/()()';iSLSLjkkgdStJtErrJIr,r'Jr'|-'|(,')/4|-'|jkrrgrrerr()/2()(),()()()';iSKKgdStJrtnJtnHrJJr'r,r'(1)CFIEEFIEMFIE4积分方程的离散•RWG矢量基函数•MOM离散1()()NiiijJrfr+-()/2,(),;()()/2,(),.iiiiiiiiiiiilAifTlAifTiρrρrrrrfrρrρrrrr()()(1)[()/2()],()()(1)().mnmnnnSiimmSzLKdSvdSfrffrnffrErnHrn=1,1,2,,.NmnnmzjvmNRao-Wilson-Glisson5MethodofMoments•SurfaceisDiscretizedintoPatches(BasisFunctions)•BasisFunctionsInteractthroughtheGreen’sFunction,Grr•GeneratesaDenseMethodofMomentsMatrix11nnnnZIV,,,''ZrrrjiSjifGfds()rf()rtPulse6线性系统：Mx=sM是N×N矩阵，x、s是N矢量●Directsolution(Gausselimination,LUDecomposition,SVD,…)空间复杂度为O(N2)，需要O(N3)次运算；●Iterativemethods,空间复杂度仍为O(N2)，如果K(kN)步收敛，每步需要的矩阵乘向量的运算为O(N2)；—Moore’Law:processorspeeddoublesevery18months—amillionvariable,16generationsofMoore’LawbeforeaO(N2)algorithmwascomparablewithaO(N)algorithm—1GBRAM=10243=1,073,741,824bytes=largestN=32,768●Finding:快速矩阵乘向量的算法(NlogN)；并行实施。7FastMultipoleMethods(FMM)•IntroducedbyRokhlin&Greengardin1987•Calledoneofthe10mostsignificantadvancesincomputingof20thcentury•Speedsupmatix-vectorproducts(sums)ofaparticulartype•以上求和要求O(MN)运算复杂度•对给定的精度，FMM可以获得O(M+N)运算复杂度•可以加速matix-vectorproducts，使O(N2)变为O(NlogN)•加速线性系统求解，如果用迭代方法，k步收敛，每步用矩阵矢量相乘，使计算复杂度由O(N3)变为O(kNlogN)1()(),{}[]{}Njijijjijisxxxs8FMM:Application•MolecularandStellardynamics—computationofforcefieldsanddynamics•Solutionofacousticalscatteringproblems—HelmholtzEquation•ElectromagneticWaveScattering—Maxwell’sEquations•FluidMechanics:Potentialflow,vertexflow—Laplace/PoissonEquations9FMM:Fundament•格林函数的加法定理•jlpl平面波展开||(2)1e(1)(21)()()()||jkllllljkljkdhkrPr+ddrr+djl_—第一类球面Bessel函数hl(2)—第二类球面Hankel函数Pl—Legendre多项式24()()()e()njknnnSajjkdPPdkddrkrk注意到lz时，函数jl(z)和hl(2)(z)幅值大致保持为常数；lz时，函数jl(z)衰减非常快，而hl(2)(z)递增非常快。当dr时，上式在保证精度的情况下截断。则上式可以写为：||(2)1e(1)(21)()()()||jkLllllljkljkdhkrPr+ddrr+dL=kd+cln(kd+)Kd-源点到观察点的最大半径c-是一个依赖希望精度的常数=1最小的相对误差小于0.1=5相对误差小于10-6=10准确到双精度10FastMultipoleBasics直接求解：分解：复杂度：O(MN)复杂度：O(pN+pM){cm}11•FMM形式的矩阵向量乘积1'()'221'()'[()]16NijjijjnmGnmjGmKpipppjjpmGnmjGmzjzjkwTjDA'2''(')(,),ln(),2;e()(')';mpmmpLmmjkpppjjjTkrLkdckdKLdSkrrTkrAI-kkfr()e[()()(1)()]mjkpppipiiidSkrrDI-kkfrkfrn近场部分远场部分电磁场积分方程的多极子展开12FMM{}{}[][][]{}NearyyDTAx1{}[]{}xAx21{}[]{}xTx——聚集过程，将基函数聚集成平面波函数，其结果表示K个平面波12{}[]{}yDx——转移过程，将得到的{x1}平移到另一个盒子的中心，其结果表示该盒子中心的K个平面波——发散过程，将得到的{x2}发散到盒子中的基函数。M2M,M2L,L2L：聚集–转移–发散M:多极子展开L:局部展开13FMM123()1()()2()()3()()(,),()(,),()(,)iiinixEnnixEnnixEnuuuiiiyyxyyxyyx由于E2(n)和E3(n)是互补的，因此对任意的n,23()()23()()()(,)()()innixEnEnuiyyxyy14FMMAlgorithmStep1M2M15Step2M2L16Step3L2LSummation17MLFMMAlgorithm•UpwardPass:mergescatteringmatrices•DownwardPass:constructsplittingandexchangematrices•FinalSummationBasedOn:◆HierarchicalGrouping◆DataStructure18L2LM2MM2MM2LM2LM2LMultilevelMultipoleOperatorsFinestLevelFinest-1LevelL2LUpTreeAcrossTreeDownTree19MLFMMAlgorithm—UpwardPass•Step1:在最细的层盒子求解远场展开系数，xi∈E1(n,l),得到C(n,l)或，这也可以用于xi∈E3(n,l)•Step2:对于l=L-1,…,2,从step1得到值进行递归得到。同样适合xi∈E3(n,l)结果：得到分层组聚集系数(,)1()nly20MLFMMAlgorithm—DownwardPass•Step1:l=2,..L递归进行E4•Step2:任意一个非空组自身加上其邻接组最多可有3d个，其父组及父组的邻接组最多可形成3d×2d，因此次相邻数目最多为p=3d(2d-1)；三维是189，二维是27，一维是3。局部到局部的转移，E3(n,l)和E4(n,l+1)产生E3(n,l+1)结果：可以得到各分层组的转移系数21KeyWords•空间多层组划分•Morton编号•相邻组的作用•远场组的上聚•次相邻组中心的转移•远场组的下推22Grouping●每个盒子在第l层的指标号数为n,那么任意盒子的指标为(n,l);●需要构建的函数，如Parent(n),ChildrenAll(n),Children(X;n,l),NeighborsAll(n,l),Neighbors(X;n,l)23Grouping每个盒子在l=0,..,L层的指标n=Number=0,1,…,2ld-1.12324(,)(,)(,){(,),((,),)}(,){0,...,21}\(,)(,){(((),1),1))}\{(,),((,),)}ldInlnlInlnlNeighbornllInlInlInlChildrenNeighborParentnllnlNeighbornll12312422(,)(,)(,)(,){((,),)}(,)(0,0)\(,)(,){((,),)}(,)((),1)\(,)((),1)ddEnlRnlEnlRnlNeighbornllEnlEEnlnlNeighbornllEnlEParentnlEnlParentnlNeighb表示盒子的空间区域;表示盒子与其相邻的空间区域;表示盒子及其相邻的外部空间区域；表示父盒子及其相邻{(),1),1}orParentnll(的空间区域；1234(,)(,)12(,)(,)(,)(,)34(,)(,)()(,),()(,),()(,),()(,)iiiinlnliixEnlxEnlnlnliixEnlxEnluuuuiiiiyyxyyxyyxyyx由于E2(n,l)和E3(n,l)是互补的，因此对任意的n,l(,)(,)23()()()nlnlyyy24纲要•FMM•DataStructures•Parallelization25DataStructure•构造树—压缩八叉树•建立连接—morton键—排序26构造树•离散点的空间层次划分•对应的四叉树及其压缩四叉树27•确定层数L—根据读入模型的最大几何尺寸确定计算区域D,根据问题的参数确定最小盒子尺寸d,树结构的层数为L=log2(D/d)•第l-1层立方体等分为八个子立方体,形成第l层更小的立方体,l-1是l层的父层,l层是l-1层的子层.•形成相邻组、次相邻组、远场组•第l层的节点数不超过2dl个构造树—八叉树(1)28构造树—八叉树（2）procedureOctree_BuildOctree={empty}Fori=1ton...对所有的点做循环Octree_Insert(i,root)...将点i插入八叉树相应的位置endFor...八叉树中可能有很多空的叶节点,但它们的兄弟节点非空Traversethetree(via,say,breadthfirstsearch)...宽度周游Eliminatingemptyleaves...去掉空的叶节点CompressOctree...压缩八叉树...如果某中间节点仅包含一个子节点，则被其替换每个节点用两个键值描述:包含相同数据的最小单元和最大单元29构造树—八叉树（3）•包含N个叶节点的压缩八叉树节点总数不超过2N-1–因此可以采用数组而不是链表来存储压缩八叉树•MLFMM基于后序周游的压缩八叉树数据结构–从键值最小的叶节点开始周游–每个节点都存储在其子节点之后，且紧挨其子节点存储•节点排序时，使用的是其所表示的最小单元键值–对于兄弟节点，按键值从小到大排序–各节点所表示的单元指的是它所包