算法公式原理:
迭代方程:Ax = b
第一步:在细网格上采用普通的GS、雅可比、ILU等迭代三次获取残差,即:r1 = b - A*x
第二步:将细网格残差插值到粗网格上,得到粗网格残差,即:r2 = R*r1 (R为限定矩阵)
第三步:在粗网格采用普通的GS、雅可比、ILU等迭代三次获取误差值e2,即:RAI * e2 = r2
第四步:将粗网格误差e2插值到细网格上,得到细网格误差e1,即:e1 = I * e2 (I为插值矩阵)
第五步:更新细网格的数值解,即:x = x + e1
返回第一步,循环迭代,直至满足收敛标准。
以上关键点:
如何得到粗网格?
第一步:定义网格间的强弱关系:
首先定义一个阈值:q = 0.5,用于后续强弱影响的计算判断
1. 找出每行方程非对角元素系数的最大值,即:maxAij(取绝对值后的最大值,因为系数均为负数)
2. 统计每个网格受到的强影响网格个数和弱影响网格个数,并存储记录好。
每行方程的对角线元素作为当前统计的网格,该行的其他非对角线元素作为被统计的对象。
判断条件:
非对角线系数 Aij > q * maxAij,则该Aij作为对Aii强影响网格,加入Sij集合中。
非对角线系数 Aij < q * maxAij,则该Aij作为对Aii弱影响网格,加入Fij集合中。
以上则可记录每个网格受到相邻网格的强弱影响。
3. 统计Aii强影响其他网格的集合,即记录被Aii强影响的网格,加入集合Gii中,Gii中的网格个数作为影响数,用于粗细网格划分(注意影响数指的的强影响的网格个数)。
判断依据:网格Aii的相邻网格Aij的强影响网格集合中是否有Aii网格,如果有,则说明该相邻网格被网格Aii强影响,则可以将该相邻网格Aij加入到Aii强影响的网格集合中。
到此,我们就完成了网格的强弱影响关系的划分,注意,这里说的相邻网格,指的就是每个方程的非对角线元素Aij,因为物理模型在组装为方程后,每个网格与其相邻网格必然组成一个方程,相邻网格即代表非对角线元素,当前网格即为该方程的对角线元素,即待求的数值解。
第二步:粗略划分粗细网格
1.遍历每一个网格,获取影响数最大的网格Aii(即第一步,第3点中统计的网格个数)
2. 将Aii加入到粗网格集合Cii中
3. Aii的Sij中的网格,可直接加入细网格集合Fii,因为这些网格被Aii强影响,可以用粗网格Aii插值得到。
4. 将对Aii的Sij中的网格强影响的网格的影响数+1,这些影响数+1的网格一定要是未划分的网格。
如此,在进入第1步,寻找剩下的网格中影响数最大的网格Aii,继续划分,知道所有网格划分完毕。
第三步:精细划分粗细网格
第二步中的网格划分可能存在漏洞,需要再次判断划分。
首先统计以下两个集合:
1. 对Aii强影响的粗网格集合Mii。
2. 对Aii强影响的细网格集合Nii。
因为经过第二步的粗细网格划分后,我们后很容易的统计出上面这两个集合来。
开始精细划分:
1. 遍历每一个网格Aii,判断该网格Aii是粗网格还是细网格,如果是粗网格则加入加集合Vii中。
2. 如果是细网格,则遍历该细网格Aii的Nii中的每个细网格,再找出每个细网格的Mii,判断该细网格的Mii的网格是至少有一个是在粗网格集合中的网格。
3. 如果没有则说明当前遍历的细网格Aii没有强影响的粗网格来插值得到。需要将该细网格Aii的Nii中的这个细网格直接加入到粗网格集合Cii中,并在细网格集合Fii中删除这个细网格,以及该Aii的Mii集合加入该细网格,而Nii集合中要删除这个细网格,因为这个细网格已经转化为了粗网格,以便能插值得到细网格Aii。
4. 将Aii的Mii网格从小到大放入集合Vii中,检查并保证Vii的元素不重复。
以上则完成了粗细网格的划分,粗网格在集合Cii中,细网格在集合Fii中。集合Vii用于计算插值矩阵I
如何得到插值矩阵I?
注意这里:插值矩阵I可采用CSR存储格式(行压缩方式)
1. 遍历每个网格,判断是否为集合Cii中的粗网格,如果是,则I矩阵对应系数为1,该网格在集合Vii中的下标即为在I矩阵的列。
2. 不是上述的粗网格,则必然是集合Fii中的细网格,此时 I 矩阵的值需要按强影响粗网格网格、强影响细网格、弱影响网格的系数权重来转化为全部由粗网格系数插值得到的系数权重分布放入矩阵I中。
粗网格对应的位置,即所在的列,仍是相应粗网格在集合Vii中的下标
给出一个小例子作为一定的参考,如下:
主意:上图中,系数权重的分母是相邻网格的系数之和。强影响的细网格需要转化为粗网格,弱影响的网格系数直接移到方程左边。
这里计算系数的涉及到粗细网格转化,同时又涉及到CSR存储的数据结构方法,需要仔细推敲。
尤其是I矩阵粗网格的排列
如何得到限定矩阵R?
只要得到插值矩阵I,将矩阵I进行转置即可得到R矩阵,那么计算RAI矩阵,再进行AMG算法的迭代就容易了
上面这一套流程我自己已经用C++代码实现并验证过了。
最后感谢博主:梅冠华的博文,我基本是按照他的博文来实现的,其中遇到了好些当时没看懂的地方,现在写这篇博文可以作为他的博文的补充,相互借鉴。