二次Lagrangian有限元方程的几何多重网格法

二次Lagrangian有限元方程的几何多重网格法
李明; 崔向照; 李郴良; 赵金娥
【期刊名称】《《高校应用数学学报A辑》》
【年(卷),期】2014(029)004
【总页数】7页(P412-418)
【关键词】二次Lagrangian有限元; 限制算子; 几何多重网格法
【作者】李明; 崔向照; 李郴良; 赵金娥
【作者单位】红河学院数学学院云南蒙自661199; 桂林电子科技大学数学与计算科学学院广西桂林541004
【正文语种】中文
【中图分类】O241.6
Lagrangian有限元法是求解偏微分问题一类重要的离散方法,通常可分为线性Lagrangian有限元(简称线性元)和高次Lagrangian有限元(简称高次元),其中高次元因对问题具有更好的逼近性而被广泛应用于实际计算.但由于高次元方程的刚度矩阵稠密性较强,代数性质较复杂,这增加了其求解的难度,因此发展有效求解高次元方程的数值算法成为科学工程计算领域的研究热点之一.
2005年,舒适教授针对高次元方程,利用有限元基函数的特殊性质,结合部分几何信息,提出了从高次元空间(细网格层)到线性元空间(粗网格层)的网格转移算子,设计了代数两层网格法(A2MG),然后使用通常的代数多重网格法(AMG)求解A2MG法中粗层网格方程,提出了求解高次元方程的AMG01法[1-2].在文[1]思想的基础上,孙
杜杜,肖映雄等人进一步丰富了求解高次元方程的AMG法[3-9].2011年文[10]讨论了基于新外推公式的所有网格层均为二次元空间的瀑布型多重网格法.
在文[1-2]的基础上,将二次元空间作为最细网格层,选取一系列线性元空间作为粗网格层,构造一种新的限制算子,提出一种求解二次Lagrangian有限元方程的几何多重网格法(GMG01),并将其计算量与AMG01法和通常的GMG法进行了比较,同时给出了数值实验.
本文考虑来源于静电场,地球物理等领域的二维泊松方程
其中,u为足够光滑的待求函数,f为已知光滑函数,Ω是一个凸多边形区域,∂Ω是Ω的边界.
取步长h1,使用三角形网格剖分Ω,采用二次元离散问题(1),可得网格空间Vq,h1和方程
取步长hj=h1×2j-1,j=1,2,···,J,基于三角形网格剖分,采用线性元离散问题(1),可得到网格空间Vl,hj和方程
几何多重网格法(GMG)是求解有限元方程(3)常用的一类算法,对于给定的初始值U01,其算法如下:
算法1 GMG法[11]
预处理.置Ai:=(i=1,2,···,J),F1:=j:=1;
步骤1.前磨光m1次,
步骤2.计算残量
步骤3.限制,Fj+1:=
步骤4.若j+1＜J,置j:=j+1,取为零向量,转步骤1.否则,精确求解AJUJ= FJ得置j:=J−1;
步骤5.校正
步骤6.后磨光m2次,
步骤7.若j＞1,置j:=j−1,转步骤5,否则,输出相关数据.
其中,分别表示第j层与第j+1层之间的限制算子和插值算子;“a:=b”表示将b赋值给表示使用某种迭代算子G,以Uj0为初始值,对方程AjUj=Fj迭代m1次得到的迭代解.
如果调用一次GMG法得到的数值解精度不够高,则可将
取作新的初值再次调用GMG法进行计算,如此反复,直至求出满足精度条件的近似解.
为了快速求解方程(2),文[1]根据二次元和线性元基函数之间的关系,结合网格的部分几何信息,通过代数途径建立二次元空间Vq,h1(细网格层）到线性元空间Vl,h1(粗网格层)的网格转移算子,形成求解方程(2)的A2MG法.对于给定的初始值
U1q,0,A2MG法流程如下:
算法2 A2MG法[1]
预处理.
(1)求限制矩阵
(2)求插值算子
(3)求粗化矩阵
步骤1.前光滑m1次,
步骤2.精确求解粗网格方程
步骤3.校正,
步骤4.后光滑m2次
为了适应求解规模更大的二次元方程,文[1]提出使用求解线性元方程的常用的AMG法,近似求解A2MG法(步骤2)中方程构造了如下算法:
算法3 AMG01法[1]
预处理.同算法2;
步骤1.同算法2;
步骤2.调用v次AMG法近似求解得EH;
步骤3.同算法2;
步骤4.同算法2.
文[1-2]利用高次元与线性元基函数之间的关系,构造了网格转移算子,提出了
AMG01法.为了进一步提高这类算法的计算效率,利用节点的邻点信息,将残量进行组合,可构造一种新的限制算子.
首先,考虑在三角形网格下线性元与二次元节点之间的关系.对于步长h,其对应的二次元的节点可分为type-a型点(即三角形的顶点)和type-b型点(即三角形各边的中点)(见图1左).网格中所有type-a型点的下标组成的集合(简称下标集)记为Sa,h.对于为i的邻点的下标组成的集合(简称邻点下标集).如图2(左)中
步长h线性元网格的节点下标集记为Ch,连接三角元各边的中点,形成加密网格(步长h/2,见图1右).记粗网格(步长h)与加密网格(步长h/2)重叠的节点记为type-c 型点,不重叠的节点记为type-d型点,网格中所有type-c型点的下标集为Sc,h/2.对于为i的邻点下标集.如图2(右)中
其次,讨论网格层之间的关系.使用一系列嵌套的线性元空间形成的多个粗网格层代替文[1-2]中的粗网格层.为便于讨论,记步长h1对应的二次元空间Vq,h1为最细网格层,步长对应的线性元空间Vl,hj(j=1,2,3,···,J)为一系列粗网格层.此时相邻网格层间的关系可分为如下两类:
Class-A:细网格层为二次元空间Vq,h1,粗网格层为线性元空间Vl,h1.易知,Vq,h1上的typea型点的下标集Sa,h1与Vl,h1的下标集Ch1存在一一对应的关系,即对∀w∈Ch1存在唯一的i∈Sa,h1使得w,i对应的节点重合.记这种对应关系为映射w=indA(i).
Class-B:细网格层为线性元空间Vl,hj,相邻粗网层为线性元空间
Vl,hj+1(j=1,2,···,J−1).易知,对于Vl,hj+1上的∀w∈Chj+1,在Vl,hj上存在唯一的i∈Sc,hj与之对应.不妨也记为映射w=indB(i).
记分别表示细网格上的残量Rq,h1,Rl,hj中的第t个元素.本文的限制算子设计如下: 对于Class-A,
对于Class-B,
针对网格层Vq,h1和Vl,hj(j=1,2,···,J),使用二次元空间Vq,h1作为最细网格层,线性元空间Vl,hj作为其它粗网格层.使用前面提出的限制算子(式(4)和式(5))代替GMG法所采用的直接映射,五点加权限制算子或九点加权限制算子(见文[11]),可构造求解方程(2)的改进的几何多重网格法(GMG01).对于给定的初始值U01,算法如下:
算法4 GMG01法
预处理.置由式(4),式(5)构造限制算子
步骤1.前磨光m1次,
步骤2.计算残量
步骤3.限制,
步骤4.若j+1＜J,置j:=j+1,转步骤1,否则,精确求解置j:=J−1;
步骤5.校正,
步骤6.后磨光m2次,
步骤7.若j＞1,置j:=j−1,转步骤5,否则,输出相关数据.
其中,表示第j+1层到第j层之间的插值算子.
最后，讨论通常的GMG法,AMG01法和GMG01法的计算量.假设使用的网格层为1个二次元空间和J个线性元空间,求解过程中它们的套迭代次数记为k,且AMG01法(步骤2)中累计调用AMG法的次数记为Pv,AMG法的网格层数取
J.GMG法,GMG01法,AMG01法及其步骤2中AMG法的前后磨光步数均取
M/2(这里M为正偶数).
由文[11]可知,插值算子和限制算子的计算量约总计算量的10%-20%,为了便于讨论,忽略掉这部分计算量.记最细网格层上磨光一次的计算量为WU(work unit)[11].对于二维问题d=2, GMG01法(或GMG法)的计算量为
AMG01法的计算量包含k次套迭代的计算量和调用v次AMG法的计算量.即
在AMG01法的步骤2中,通常需要调用多次AMG法才能求得满足精度条件的数值解,即有v≥k,对比式(6)和式(7),可知在相同的套迭代次数情况下,AMG01法的计算量要比GMG01法(或GMG法）多
为验证新算法的有效性,取如下算例进行数值实验.
算例1真解区域Ω:(0,1)×(0,1),右端函数
算例2真解区域Ω:(0,1)×(0,1),右端函数
采用均匀网格剖分以上算例.选取AMG01法,GMG法作为对比算法,并在相同网格层数(文中取J=4)的前提下比较它们的计算效率.
AMG01法步骤2中所用的AMG法的网格层数取为J,AMG法的粗化算法和插值算子取自文[12].
GMG法所用网格空间和网格层数均与GMG01法相同,GMG法的限制算子为五点加权限制算子,插值算子为双线性插值.
在GMG01法(步骤5)中,分别使用三次多项式插值,三次样条插值,双线性插值的算法对应记为GMG01(C),GMG01(S)和GMG01(L).
GMG01法,AMG01法,GMG法的初始值均取为零向量,前后磨光步数取3.这些算法的套迭代终止条件为套迭代100次或(这里r0,rk分别表示初始值和第k次套迭代得到的数值解对应的残量,||.||为2范数).P
表1-3中,“Mesh”表示二次元网格的规模;“iter”表示套迭代次数;“v”表示整个求解过程中,AMG01法(步骤2)中累计调用AMG法的次数;“cpu”表示算法从
预处理直至结束所需的时间(单位:秒);”∗∗∗”表示套迭代100次仍未得到满足精度条件的解.
由表1-3可以看出:
1)通常的GMG法的计算效率远低于GMG01法和AMG01法.
2)GMG01法所需套迭代次数,与AMG01法基本相当.
3)GMG01法的计算量约为AMG01法的一半.
4)GMG01法和AMG01法的套迭代次数不随求解规模的增加而增加,均具有较强的稳健性.
5)采用三次样条插值和三次多项式插值的GMG01法,略优于采用双线性插值的GMG01法.
以上结果表明,与AMG01法和GMG法相比,本文提出的GMG01法计算量少,计算时间短,稳健性强.
致谢:感谢评审专家提出的宝贵评审意见!
【相关文献】
[1]舒适.一种部分几何和分析信息的代数多重网格法及其应用[D].湘潭:湘潭大学,2004.
[2]Shu Shi,Sun Dudu,Xu Jinchao.An Algebraic Multigrid Method for Higher-order Finite Element Discretizations[J].Computing,2006,(77):347–377.
[3]孙杜杜,舒适.求解三维高次拉格朗日有限元方程的代数多重网格法[J].计算数学,2005,27(1): 101-112.
[4]刘轩,舒适.非结构四边形二次Lagrangian有限元方程的代数多重网格法及收敛性分析[J].高等学校计算数学学报,2005,27:215-222.
[5]赵钍焱.求解二维弹性力学问题高次有限元方程的代数多重网格法[D].湘潭:湘潭大学,2006.
[6]张红梅.求解三维弹性力学问题高次有限元方程的代数多层网格法[D].湘潭:湘潭大学,2007.
[7]张红梅,肖映雄,舒适.三维椭圆问题三次有限元方程的代数多层网格法[J].数值计算与计算机应用,2008,29(4):277-290.
[8]肖映雄,张红梅,舒适.三维弹性问题高次有限元方程的代数多层网格法[J].计算力学学报, 2010,27(6):995-1000.
[9]李明,李郴良,崔向照.求解二次Lagrangian有限元方程的代数两水平法[J].应用数学学报, 2013,36(5):1-9.
[10]李明,李郴良.求解二次Lagrangian有限元方程的新瀑布型多重网格法[J].山西大学学报,2011, 34(2):185-189.
[11]Briggs W L,Henson V E,McCormick S F.A Multigrid Tutorial(第2版)[M].北京:清华大学出版社,2011.
[12]Ruge J,Stuben K.Algebraic multigrid[A],in:McCormick S F,ed.,Multigrid Methods[C]. Philadelphia:SIAM,1987.
[13]谭敏,肖映雄,舒适.一种各向异性四边形网格下的代数多重网格法[J].湘潭大学自然科学学报, 2005,27(1):78-84.
[14]李明,李郴良.非均匀三角形网格下椭圆问题的代数两网格法[J].高校应用数学学报,2013,
28(4):439-446.。