錢 皓, 吳 丹, 呂俊良

(吉林大學(xué) 數(shù)學(xué)學(xué)院, 長(zhǎng)春 130012)

為提高計(jì)算效率, 充分利用機(jī)器資源, 設(shè)計(jì)適合并行化的計(jì)算策略已成為計(jì)算數(shù)學(xué)領(lǐng)域內(nèi)的重要課題.在求解偏微分方程的數(shù)值方法上, 目前對(duì)差分法和有限體積法的并行化研究已有許多成果[1-5], 其中區(qū)域分解方法在格式并行化中具有重要作用, 它將計(jì)算區(qū)域劃分為若干個(gè)子區(qū)域, 并對(duì)子區(qū)域進(jìn)行特殊的界面處理, 將復(fù)雜問題轉(zhuǎn)化為若干個(gè)子問題, 使構(gòu)造出的格式有高度并行的特點(diǎn).有限體積元法作為計(jì)算流體力學(xué)的常用數(shù)值方法, 具有高階逼近且能自然保持局部守恒性的優(yōu)點(diǎn).目前該方法的研究主要包括格式的穩(wěn)定性分析[6-8]和收斂性分析[9-11]等.此外, 目前在有限體積元法并行化研究方面也取得了一些成果, 其中文獻(xiàn)[12]在四邊形網(wǎng)格上提出了雙線性元有限體積的并行格式和保正修正后的并行格式.本文采用文獻(xiàn)[12]中的區(qū)域分解思想, 通過預(yù)估校正方法, 在隨機(jī)三角形網(wǎng)格上構(gòu)造一種非條件穩(wěn)定的線性元有限體積并行格式, 并用數(shù)值算例驗(yàn)證該并行格式的收斂性和并行效率.

1 標(biāo)準(zhǔn)有限體積法

1.1 模型問題

考慮擴(kuò)散問題

(1)

其中Ω是一個(gè)帶有邊界?Ω的二維多邊形區(qū)域,X=(x,y)T∈2,f(X,t)為源項(xiàng),g(X,t)為邊界條件,φ(X)為初值條件,κ=(kij(x,y))2×2∈2×2為擴(kuò)散張量.假設(shè)系數(shù)矩陣κ對(duì)稱且有一致的上下界, 即存在常數(shù)0<γ1<γ2, 使得γ1(ξ,ξ)≤(κξ,ξ)≤γ2(ξ,ξ)對(duì)任意實(shí)向量ξ∈2和都成立.并設(shè)函數(shù)f(X,t),g(X,t),φ(X)充分光滑.

1.2 網(wǎng)格剖分與函數(shù)空間

圖1 三角形單元△P1P2P3的對(duì)偶剖分Fig.1 Dual partition of triangular unit △P1P2P3

圖2 頂點(diǎn)P0所在的對(duì)偶單元Fig.2 Dual unit of vertex P0

對(duì)于給定的t∈[0,T], 取試探函數(shù)空間為

其中P1(K)表示在單元K上的一次多項(xiàng)式集合,uh(X,t)|K由K的3個(gè)頂點(diǎn)在t時(shí)刻的函數(shù)值唯一確定.取檢驗(yàn)函數(shù)空間Vh為相應(yīng)于對(duì)偶剖分的分片常數(shù)空間, 即

其中P0(K)表示在單元K*上的常數(shù)多項(xiàng)式集合.事實(shí)上, 對(duì)于每個(gè)vh∈Vh, 均可唯一地表示為

1.3 有限體積元離散格式

(2)

成立, 其中

(3)

成立, 其中

2 線性元有限體積并行格式

2.1 區(qū)域分解

圖3 分成4個(gè)并行子區(qū)域時(shí)的節(jié)點(diǎn)Fig.3 Nodes when domain is divided into 4 parallel sub domains

將區(qū)域Ω分成R個(gè)子區(qū)域Ωi(i=1,2,…,R), 用○和●表示兩個(gè)相鄰子區(qū)域之間界面上的網(wǎng)格節(jié)點(diǎn), □和■表示子區(qū)域內(nèi)部的網(wǎng)格節(jié)點(diǎn),×表示?Ω上的網(wǎng)格節(jié)點(diǎn).由不同子區(qū)域共享的子邊界稱為區(qū)域分解的交界面, 不妨設(shè)有K個(gè), 用Bk表示第k(k=1,2,…,K)個(gè)交界面上節(jié)點(diǎn)的集合.同時(shí)位于多個(gè)交界面上的節(jié)點(diǎn)稱為區(qū)域分解的交點(diǎn), 記作Oj(j=1,2,…,J).圖3為Ω被分成4個(gè)子區(qū)域Ωi(i=1,2,3,4)時(shí)的節(jié)點(diǎn)分布情況.

2.2 并行算法

算法1線性元有限體積并行算法.

注1本文利用前兩個(gè)時(shí)間層上子區(qū)域交界面處的值, 采用算術(shù)平均意義下的線性外推方法, 對(duì)子區(qū)域交界面處的值進(jìn)行預(yù)估, 將復(fù)雜問題轉(zhuǎn)化成了若干個(gè)Dirichlet邊值子問題進(jìn)行并行求解；再利用求出的值對(duì)交界面進(jìn)行修正.故對(duì)于每個(gè)子進(jìn)程, 其所需保存的數(shù)據(jù)為前兩個(gè)時(shí)間層上子區(qū)域交界面處的值與前一個(gè)時(shí)間層上子區(qū)域內(nèi)部的值.

3 數(shù)值實(shí)驗(yàn)

下面用一些數(shù)值算例檢驗(yàn)線性元有限體積并行格式的收斂性與并行效率.引入L2,H1和L∞的范數(shù):

定義格式的收斂階為

(4)

其中h1,h2表示網(wǎng)格的尺度,δ(h1),δ(h2)表示該網(wǎng)格單元數(shù)下的誤差.為測(cè)試格式的并行效率, 定義計(jì)算公式:

(5)

其中Np為并行程序所用的核數(shù),T1為算法串行花費(fèi)的時(shí)間,Tp為使用p核并行算法花費(fèi)的時(shí)間.令τ=T/N, 則網(wǎng)比μ=τ/h2.

本文采用Fortran編譯語言, 并使用MPICH2的函數(shù)庫(kù), 通過MPI通訊方式實(shí)現(xiàn)程序的并行.數(shù)值實(shí)驗(yàn)在北京超級(jí)云計(jì)算中心M分區(qū)(BSCC-M)上進(jìn)行, 本地環(huán)境為Windows10, 運(yùn)行環(huán)境為L(zhǎng)inux-CentOS7, 單個(gè)計(jì)算節(jié)點(diǎn)的處理器為512 GB內(nèi)存, 128核的AMD EPYC CPU.為達(dá)到并行的效果, 本文對(duì)每個(gè)子區(qū)域都采用單獨(dú)一個(gè)核進(jìn)行運(yùn)算, 因此程序所需的核數(shù)等于子區(qū)域的數(shù)量.在求解程序中的線性方程組時(shí), 本文采用GMRES(m)迭代算法, 迭代終止值為10-16.

例1各向同性擴(kuò)散問題.

將計(jì)算區(qū)域分為16個(gè)子區(qū)域, 設(shè)起始時(shí)刻為0, 終止時(shí)刻為0.1 s, 網(wǎng)比為5, 分別計(jì)算如圖4所示的攝動(dòng)網(wǎng)格Ⅰ和如圖5所示的Sine網(wǎng)格上解的誤差, 得到其計(jì)算精度, 分別如圖6和圖7所示.設(shè)起始時(shí)刻為0, 終止時(shí)刻為0.001 s, 時(shí)間剖分?jǐn)?shù)為1 000, 網(wǎng)格尺寸為1/1 000, 固定一個(gè)攝動(dòng)網(wǎng)格Ⅰ, 記錄不同并行區(qū)域數(shù)時(shí)并行格式所需的計(jì)算時(shí)間, 并計(jì)算并行效率, 結(jié)果列于表1.實(shí)驗(yàn)結(jié)果表明, 對(duì)于各向同性擴(kuò)散問題, 并行格式的數(shù)值解按L2模2階收斂, 按H1模1階收斂, 按L∞模2階收斂, 且具有良好的并行效率.

圖4 攝動(dòng)三角形網(wǎng)格ⅠFig.4 Distorted trianglar mesh Ⅰ

圖5 Sine三角形網(wǎng)格Fig.5 Sine trianglar mesh

圖6 攝動(dòng)網(wǎng)格Ⅰ上的計(jì)算精度Fig.6 Computational accuracy on distorted mesh Ⅰ

圖7 Sine三角形網(wǎng)格上的計(jì)算精度Fig.7 Computational accuracy on Sine trianglar mesh

表1 攝動(dòng)網(wǎng)格Ⅰ上的并行效率(T=0.001, N=1 000)

例2各向異性擴(kuò)散問題.

考慮各向異性的擴(kuò)散方程(1), 其擴(kuò)散系數(shù)為

其中α為可變參數(shù), 方程的精確解為u(x,y,t)=e-π2tsin(πx)sin(πy).

圖8 當(dāng)α=10-6時(shí)的計(jì)算精度Fig.8 Computational accuracy when α=10-6

圖9 當(dāng)α=102時(shí)的計(jì)算精度Fig.9 Computational accuracy when α=102

將區(qū)域分為16個(gè)子區(qū)域, 設(shè)起始時(shí)刻為0, 終止時(shí)刻為0.2 s, 網(wǎng)比μ=1, 分別取α為10-6和102, 采用線性元有限體積并行格式, 計(jì)算攝動(dòng)網(wǎng)格Ⅰ上數(shù)值解的誤差, 結(jié)果分別如圖8和圖9所示.設(shè)起始時(shí)刻為0, 終止時(shí)刻為0.001 s, 時(shí)間剖分?jǐn)?shù)為1 000, 網(wǎng)格尺寸為1/1 000,α=10-2, 計(jì)算在同一個(gè)攝動(dòng)網(wǎng)格Ⅰ上使用不同核數(shù)時(shí), 格式的并行效率, 結(jié)果列于表2.實(shí)驗(yàn)結(jié)果表明, 對(duì)于各向異性的擴(kuò)散問題, 該并行格式能達(dá)到最佳收斂, 且并行效率達(dá)到0.7以上.

表2 α=10-2時(shí)的并行效率(T=0.001, N=1 000)

例3間斷的擴(kuò)散系數(shù)問題.

考慮帶有非連續(xù)擴(kuò)散系數(shù)的拋物方程(1), 其擴(kuò)散系數(shù)為

方程的精確解為

將區(qū)域分為9個(gè)子區(qū)域, 設(shè)起始時(shí)刻為0, 終止時(shí)刻為0.01 s, 網(wǎng)比為1, 分別計(jì)算如圖10所示的攝動(dòng)網(wǎng)格Ⅱ和如圖11所示的梯形網(wǎng)格上數(shù)值解的誤差, 得到其計(jì)算精度，分別如圖12和圖13所示.設(shè)起始時(shí)刻為0, 終止時(shí)刻為0.001 s, 時(shí)間剖分?jǐn)?shù)為N=1 000, 網(wǎng)格尺寸為1/900, 用線性元有限體積并行格式, 在梯形網(wǎng)格上計(jì)算不同核數(shù)時(shí)的并行效率, 結(jié)果列于表3.實(shí)驗(yàn)結(jié)果表明, 對(duì)于間斷的擴(kuò)散問題, 該并行格式的誤差也能達(dá)到最佳收斂階, 且并行效率良好.

圖10 攝動(dòng)網(wǎng)格ⅡFig.10 Distorted mesh Ⅱ

圖11 梯形網(wǎng)格Fig.11 Trapezoidal mesh

圖12 攝動(dòng)網(wǎng)格Ⅱ上的計(jì)算精度Fig.12 Computational accuracy on distorted mesh Ⅱ

圖13 梯形網(wǎng)格上的計(jì)算精度Fig.13 Computational accuracy on trapezoidal mesh

表3 梯形網(wǎng)格上的并行效率(T=0.001, N=1 000)

綜上所述, 本文從微分方程的數(shù)值格式出發(fā), 在扭曲的三角形網(wǎng)格上, 提出了一種適用于在多核機(jī)器上運(yùn)算的線性元有限體積并行格式.該格式采用預(yù)估校正技術(shù)實(shí)現(xiàn)了其可并行性.為不失去原有格式的收斂性, 本文在子區(qū)域交界面上采用了前兩個(gè)時(shí)間層上值的線性外推方法進(jìn)行預(yù)估.該并行格式只需在相鄰的子區(qū)域之間進(jìn)行數(shù)據(jù)的通訊, 且只有區(qū)域交界面周圍的點(diǎn)參加了通訊, 故不會(huì)產(chǎn)生大量的時(shí)間消耗.數(shù)值實(shí)例結(jié)果表明, 本文格式不僅可極大減少計(jì)算時(shí)間, 而且能保持原有的收斂速度.