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

小长方体多大合适
•
•
这里我们讨论小长方体多大合适的问题？
由上述分析可以知道：转移步骤所需计算量很小。然而，聚集和发散步骤并非如此。实 际上，将原来的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）
快速多极子技术的数学原理
ij j i 离散积分方程系数矩阵的元素可表示成 （2.71） 这里P（G）表示作用在格林函数G上的算子。假设{x}和{y}分别代表相距较远的两个小长方体A、B中的 的未知数。那么{x}对{y}的作用可表示成 {y}=[Z]{x} （2.72） 快速多极子技术将此矩阵和矢量相乘分解成聚集、转移、发散三步骤进行。下边具体介绍此分解过程。 很简单，主要靠下面两个数学恒等式。第一个便是关于格林函数的加法定律 e jk |r d | （2.73） l (2) jk (1) (2l 1) jl (kd )hl (kr ) Pl (d r ) |rd | l 0 (2) 这里jl是第一类球面Bessel函数，hl 是第二类球面Hankel函数，Pl是Legendre多项式，以及d<r.值得 注意的是，在l<z时，函数jl(z)和hl(z)幅度大致保持常数；在l〉z时，jl(z)衰减非常快，而hl2(z) 递增非常快。这样当d<<r时，式（2.73）能在保证精度下截断。这样展开（2.73）便可以写成 L （2.74） e jk|r d | jk (1)l (2l 1) jl (kd )hl (2) (kr ) Pl (d r ) |rd | l 0 通常取L=kd+2ln(kd+pi)就能保证较高精度了。第二个恒等式便是式（2.74）中jlPl的平面波展开
快速多极子技术的数学原理（续）
Zij DipTp Apj
K p 1
（2.80）
其中
jk p rmo Dip p gi e dS S
（2.81）
Apj g j e
S'
jk p rm'o'
dS'
(2.82)
注意Dip、Apj是矢量，Tp是标量。于是式（2.72）的矩阵和矢量相乘便可表示成 (2.83)
•
•
多层快速多极子技术实施办法示意图
m’
m’
m 转移过程
m
m’
聚集过程
m’
m 发散过程
m
图4 多层快速多极子技术实施办法示意图
多层快速多极子技术的基本思路
• 为了给出一个完整的多层快速多极子技术实现步骤，不失一般性，下面以求解域是正方形 为例来具体说明。如图5所示，将正方形分成四个小正方形，这是第一层分组，其中“*” 号表示小组中心。后再将每一小正方形一分为四，此过程反复进行，直到最小正方形边长 在半个波长为止。为了叙述方便，这里只给出三层分组。图5(b)和图5(c)分别给出第二层 和第三层的分组情况，其中符号“。”和“.”分别表示第二层和第三层的小组中心。 和快速多极子技术一样，多层快速多极子技术也是将矩阵与矢量相乘分解成聚集、转移、 发散三步骤。然聚集和发散过程是先在最高层进行，后通过移置中心、以及插值完成低层 的聚集和发散。至于转移过程，与快速多极子技术中的执行方式完全一样，只是多层快速 多极子技术只在同一层的次相邻中心（后面将解释）进行。 先考虑聚集过程。以第二层的聚集为例。这层的聚集是通过将其高一层即第三层平面波的 中心先移置、后插值完成的。假设在第三层中的每一小正方体有N3个未知数，这样便有 M=N/N3个小正方体。因为聚集在小正方体中心“.”的平面波个数通常小于未知数个数N3， 因此从中心“.”到中心“。”的移置所需计算量要小于N3xM=N。同样可以分析得到其他层 的移置过程所需计算量也要小于N。下边会阐述每一层中插值过程所需计算量也大致为O （N）。这样执行每一层的聚集过程计算量也就在O（N）量级。发散过程相似，计算量也 在O（N）量级。
FEKO Suite 5.5 基本算法介绍
快速多极子（FMM） 多层快速多极子技术（MLFMM）
FMA overview
•
我们知道解线性方程组的方法可分为两类：一类是直接法，如高斯消元法等； 一类是迭代法，如：共轭梯度法等 用矩量法（MoM）求解线性方程组，它的系数矩阵是满秩的。如用直接法求解， 则计算机内存需要O（N2），运算量达O（N3）；如用迭代法求解，内存一样 需要O（N2），而每次迭代的运算量达O(N2)).如此之多的内存需要量，如此之 大的运算量，大大限制了矩量法的应用范围，在90年代以前，矩量法仅仅适用 于电小尺寸物体（物理尺寸/工作波长< 10）。 20世纪90年代以后，情况发生了改变，目前矩量法已可以计算相当大的电大尺 寸物体，这主要归归功于Rokhlin提出的快速多极子算法，这是一种减少内存需 求，加快矩阵和矢量相乘的技术。
Байду номын сангаас
Z
g P (G ) g dS ' dS
2 jk d 4 ( j) jl (k0d ) P l (d r ) e Pl (k r )d kˆ l
（2.75）
快速多极子技术的数学原理（续）
ˆ ”表示是积分在整个单位球面上进行。此积分可用高斯面积分方法进行。 这里的积分符号“d k 具体来说,就是在区间[0，pi]上取L点，使得cos( )在区间[-1，1]上满足Gauss-Legendre L点 积分公式。这L点的theta值及积分权因子可以直接调用可以直接调用文献[7]中的子程序“”得 到。对于水平方向phi值得选取，可以在区间[0，2pi]上等间隔选取2L个值。于是式（2.75）右 边的积分便可写成： （2.76）
快速多极子技术的数学原理（续）
m’ O’
m
O
图3：快速多极子技术具体实施示意图
多层快速多极子技术的基本思路
• 由前边的讨论可知，快速多极子技术中的组不能太大，因为那样转移过程虽然能非常有 效地计算，但聚集和发散过程都不能有效进行；组也不能太小，因为那样聚集和发散过 程能够虽能有效进行，然转移过程又不能有效计算。为此，我们是通过组的恰当大小来 获得快速多极子技术的最佳效率。 这里我们将介绍一种新的方式来更有效地实现快速多极子技术。其基本思路就是将未知 数分成不同层次的组，低层组大，高层组小，让聚集和发散过程先在最高层进行，后通 过移置、插值完成底层中的聚集和发散，而转移过程只在每层的部分组之间进行。这种 实现方式被称为多层快速多极子技术。 这里，我们以聚集过程为例来具体阐述这一实现方式。如图4所示，假设大组中有4m个 未知数，这样实现聚集需16m2次计算机操作。如果聚集先在小组进行，需4m2次计算机 操作。后将所得的四类以小组中心为起点的平面波移到以大组中心为起点，并相加得到 m个以大组中心为起点的平面波，这又需4m次计算机操作。接着再将m个以大组中心为 起点的平面波插值，得到4m个大组中心平面波，从而完成大组聚集过程。后边我们会说 明此插值过程需64m次计算机操作。因此这种实现方式总共需4m2+4m+64m次计算机操 作。在m很大时，明显少于原来的16m2次计算机操作。
远相互作用示意图
100 个 未知数
100 x 100次计算机操作
100 个 未知数
图1：远相互作用常规实施办法示意图
远相互作用FMM实施方法
小于100个平面波 小于100个平面波
O’
转移过程 100次计算机操作
O
聚集过程 100x100次计算机操作 100x100次计算机操作
发散过程
图2：远相互作用FMA实施办法示意图
（2.78）
pVp (k p d )Tp (kr, k p r )
p 1
K



（2.77）
(2)
l 0 注意式（2.77）右边中的 Vp (k p d ) 与r无关，而 T (kr , k p r )
Tp (kr , k p r ) ( j ) (2l 1)hl
{ y} [ D][T ][ A]{x}
矩阵[D]和[A]的元素已由式（2.81），（2.82）给出。矩阵[T]的元素可写成Tpq= Tp pq 。 显然，[T]是对角阵。表达式(2.83)便是快速多极子技术将直接相互作用分解成聚集、转 移、发散三步骤的数学表达式。在这表达式中，[A]{x}表示聚集过程。它将A中基函数聚 集得到的平面波函数，所得结果{x1}表示K个平面波；[T]{x1}表示转移过程。它将聚集 得到的平面波从A中心转移到B中心，所得结果{y1}表示在B中心的K个平面波；最后 [D]{y1}就是发散过程。它将K个平面波发散到B中的基函数，从而得到最终结果{y}。