06_2多层快速多极子技术MLFMM

FEKOSuite5.5基本算法介绍快速多极子（FMM）多层快速多极子技术（MLFMM）FMAoverview•我们知道解线性方程组的方法可分为两类：一类是直接法，如高斯消元法等；一类是迭代法，如：共轭梯度法等•用矩量法（MoM）求解线性方程组，它的系数矩阵是满秩的。如用直接法求解，则计算机内存需要O（N2），运算量达O（N3）；如用迭代法求解，内存一样需要O（N2），而每次迭代的运算量达O(N2)).如此之多的内存需要量，如此之大的运算量，大大限制了矩量法的应用范围，在90年代以前，矩量法仅仅适用于电小尺寸物体（物理尺寸/工作波长10）。•20世纪90年代以后，情况发生了改变，目前矩量法已可以计算相当大的电大尺寸物体，这主要归归功于Rokhlin提出的快速多极子算法，这是一种减少内存需求，加快矩阵和矢量相乘的技术。快速多极子技术的基本思路•大家都知道，迭代法的运算量主要取决于矩阵与矢量相乘的运算量。从物理意义上看，矩阵与矢量相乘虽实际上是源点对场点的作用，然由于整个考虑的源点和场点是重合的，因此也可以形式地认为是等效电流之间的相互作用。•快速多极子技术的基本思路是首先将未知等效电流分成小组。分小组可以按如下方式进行：首先用一个适当大小的长方体将物体刚好包住，然后将此长方体分成小长方体（小长方体究竟多大合适，下面再作具体分析），并将非空小长方体标出储存。此处非空小长方体是指其内有未知等效电流的小长方体，也就是被物体边界相割的小长方体。对任何一个非空小长方体，其他的非空小长方体可以分成两类：一类为近相互作用，一类是远相互作用。通常，两小长方体中心之间的距离小于半个波长的为近相互作用，否则为远相互作用。近远相互作用介绍•下边来分析两小长方体A和B的远相互作用。设A和B内分别都有100个未知数，如图1所示。如果用通常方式来执行他们之间的相互作用，则需要100*100次计算机操作。而快速多极子技术是用一种新的方式来执行A和B之间的远相互作用。其基本思路是将整个相互作用过程分解成三步：聚集、转移、发散。聚集就是将分布在A内的100个未知数所对应的等效电流聚集在A的中心。其目的是获得一组具有下列转移特性的新函数：A内所有等效电流对远处的作用可以由执行这组函数的转移完成；转移就是将聚集过程中得到的一组函数由A的中心转移到B的中心；发散就是将转移到B中心的那组函数发散到B内所有100个未知数所对应的等效电流上，从而完成A和B的远相互作用。此种作用方式由图2表示。下边会阐述平面波函数具有上述转移特性，而且在能够保证高精度情况下，所需平面波个数少于原未知数个数。这就是说，完成新函数从A中心到B中心的转移，只需要少于100次的计算机操作。这就是快速多极子技术能够加快完成A和B远相互作用的原因。作用过程的分解来源于积分方程中格林函数的多极子展开，故此项技术称为快速多极子技术。由于格林函数的多极子展开在近相互作用时很难达到满意精度，则这种新作用方式只适用于远相互作用。这也就是我们将相互作用分成近相互作用和远相互作用的原因。远相互作用示意图图1：远相互作用常规实施办法示意图100个未知数100个未知数100x100次计算机操作远相互作用FMM实施方法图2：远相互作用FMA实施办法示意图小于100个平面波小于100个平面波聚集过程发散过程转移过程100次计算机操作100x100次计算机操作100x100次计算机操作O’O小长方体多大合适•这里我们讨论小长方体多大合适的问题？•由上述分析可以知道：转移步骤所需计算量很小。然而，聚集和发散步骤并非如此。实际上，将原来的100个未知数所对应的100个基函数聚集成大致100个平面波函数，需要100x100次计算机操作；将100个平面波函数发散给100个未知数所对应的100个基函数，也需要100x100次计算机操作。因此，如果只考虑两个小长方体的远相互作用，快速多极子技术所需计算量是超过通常方式的计算量。然而，当我们考虑100个小长方体的远相互作用时，情况就不同了。此时有100x100=10000个未知数，用通常方式完成它们的相互作用需10000x10000次计算机操作。如用快速多极子技术，每个小长方体中的聚集需100x100次计算机操作，现有100个小长方体，因此整个聚集需100x100x100次计算机操作。同样道理，整个发散需100x100x100次计算机操作。至于转移步骤，因为完成一次转移需100次计算机操作，现有100个小长方体，需100x100次转移，因此整个转移也需100x100x100次计算机操作，故用快速多极子技术完成整个相互作用需3*100x100x100次计算机操作，大大少于通常方式的计算量。由此可见，小长方体尺寸不能过大，因为过大会导致聚集和发散两步骤地计算量过大；小长方体尺寸不能过小，因为过小会导致转移步骤的计算量过大。严格来说，假如有N个未知量，分成M组，这样每组大致有N/M个未知数。根据上面分析，聚集和发散步骤所需计算机操作都是O（N2/M），转移步骤需O（NM）次计算机操作。因此完成整个相互作用需要O（N2/M+MN）次计算机操作。不难知道，在M=N1/2时，完成整个相互作用需要的计算机操作次数最少为O（N3/2）。不难分析，此时的内存需要量也为O（N3/2）快速多极子技术的数学原理离散积分方程系数矩阵的元素可表示成（2.71）这里P（G）表示作用在格林函数G上的算子。假设{x}和{y}分别代表相距较远的两个小长方体A、B中的的未知数。那么{x}对{y}的作用可表示成{y}=[Z]{x}（2.72）快速多极子技术将此矩阵和矢量相乘分解成聚集、转移、发散三步骤进行。下边具体介绍此分解过程。很简单，主要靠下面两个数学恒等式。第一个便是关于格林函数的加法定律（2.73）这里jl是第一类球面Bessel函数，hl(2)是第二类球面Hankel函数，Pl是Legendre多项式，以及dr.值得注意的是，在lz时，函数jl(z)和hl(z)幅度大致保持常数；在l〉z时，jl(z)衰减非常快，而hl2(z)递增非常快。这样当dr时，式（2.73）能在保证精度下截断。这样展开（2.73）便可以写成（2.74）通常取L=kd+2ln(kd+pi)就能保证较高精度了。第二个恒等式便是式（2.74）中jlPl的平面波展开（2.75）||(2)0(1)(21)()()()||jkrdlllllejkljkdhkrPdrrd()'ijijZgPGgdSdS||(2)0(1)(21)()()()||jkrdLlllllejkljkdhkrPdrrd20ˆ4()()()()ljkdllljjkdPdrePkrdk这里的积分符号“”表示是积分在整个单位球面上进行。此积分可用高斯面积分方法进行。具体来说,就是在区间[0，pi]上取L点，使得cos()在区间[-1，1]上满足Gauss-LegendreL点积分公式。这L点的theta值及积分权因子可以直接调用可以直接调用文献[7]中的子程序“”得到。对于水平方向phi值得选取，可以在区间[0，2pi]上等间隔选取2L个值。于是式（2.75）右边的积分便可写成：（2.76）这里K=2L2，wp是权因子，,,表示单位球面上取样点的球坐标。将式（2.76）代入式（2.74），并将求和次序交换可得其中（2.77）（2.78）（2.79）注意式（2.77）右边中的与r无关，而与d无关。这表明式(2.77)已将格林函数表示的直接相互作用分解成远距离的转移和近距离的聚集或发散。为了更简明的阐述，不失一般性，以P(G)=G为例说明。如图3所示，取，利用（2.77）便可得矩阵元素（2.71）的FMM表达式快速多极子技术的数学原理（续）2ˆdk(2)0(,)()(21)()()LlppllplTkrkrjlhkrPkr()ppVkd00''',momorrrdrr(,)pTkrkr(sincos,sinsin,cos)ppppppkppkkk(,)pp||1()(,)||4jrdKppppppejkVkdTkrkrrd快速多极子技术的数学原理（续）（2.80）其中（2.81）(2.82)注意Dip、Apj是矢量，Tp是标量。于是式（2.72）的矩阵和矢量相乘便可表示成(2.83)矩阵[D]和[A]的元素已由式（2.81），（2.82）给出。矩阵[T]的元素可写成Tpq=。显然，[T]是对角阵。表达式(2.83)便是快速多极子技术将直接相互作用分解成聚集、转移、发散三步骤的数学表达式。在这表达式中，[A]{x}表示聚集过程。它将A中基函数聚集得到的平面波函数，所得结果{x1}表示K个平面波；[T]{x1}表示转移过程。它将聚集得到的平面波从A中心转移到B中心，所得结果{y1}表示在B中心的K个平面波；最后[D]{y1}就是发散过程。它将K个平面波发散到B中的基函数，从而得到最终结果{y}。m'o''d'pjkrpjjSAgeSmodpjkrippiSDgeS1KijipppjpZDTA{}[][][]{}yDTAxppqT快速多极子技术的数学原理（续）O’Om’m图3：快速多极子技术具体实施示意图多层快速多极子技术的基本思路•由前边的讨论可知，快速多极子技术中的组不能太大，因为那样转移过程虽然能非常有效地计算，但聚集和发散过程都不能有效进行；组也不能太小，因为那样聚集和发散过程能够虽能有效进行，然转移过程又不能有效计算。为此，我们是通过组的恰当大小来获得快速多极子技术的最佳效率。•这里我们将介绍一种新的方式来更有效地实现快速多极子技术。其基本思路就是将未知数分成不同层次的组，低层组大，高层组小，让聚集和发散过程先在最高层进行，后通过移置、插值完成底层中的聚集和发散，而转移过程只在每层的部分组之间进行。这种实现方式被称为多层快速多极子技术。•这里，我们以聚集过程为例来具体阐述这一实现方式。如图4所示，假设大组中有4m个未知数，这样实现聚集需16m2次计算机操作。如果聚集先在小组进行，需4m2次计算机操作。后将所得的四类以小组中心为起点的平面波移到以大组中心为起点，并相加得到m个以大组中心为起点的平面波，这又需4m次计算机操作。接着再将m个以大组中心为起点的平面波插值，得到4m个大组中心平面波，从而完成大组聚集过程。后边我们会说明此插值过程需64m次计算机操作。因此这种实现方式总共需4m2+4m+64m次计算机操作。在m很大时，明显少于原来的16m2次计算机操作。多层快速多极子技术实施办法示意图m’m’m’m’mmmm聚集过程发散过程转移过程图4多层快速多极子技术实施办法示意图多层快速多极子技术的基本思路•为了给出一个完整的多层快速多极子技术实现步骤，不失一般性，下面以求解域是正方形为例来具体说明。如图5所示，将正方形分成四个小正方形，这是第一层分组，其中“*”号表示小组中心。后再将每一小正方形一分为四，此过程反复进行，直到最小正方形边长在半个波长为止。为了叙述方便，这里只给出三层分组。图5(b)和图5(c)分别给出第二层和第三层的分组情况，其中符号“。”和“.”分别表示第二层和第三层的小组中心。•和快速多极子技术一样，多层快速多极子技术也是将矩阵与矢量相乘分解成聚集、转移、发散三步骤。然聚集和发散过程是先在最高层进行，后通过移置中心、以及插值完成低层的聚集和发散。至于转移过程，与快速多极子技术中的执行方式完全一样，只是多层快速多极子技术只在同一层的次相邻中心（后面将解释）进行。•先考虑聚集过程。以第二层的聚集为例。这层的聚集是通过将其高一层即第三层平面波的中心先移置、后插值完成的。假设在第三层中的每一小正方体有N3个未知数，这样便有M=N/N3个小正方体。因为聚集在小正方体中心“.”的平面波个数通常小于未知数个数N3，因此从中心“.”到中心“。”的移置所需计算量要小于N3xM=N。同样可以分析得到其他层的移置过程所需