線性元有限體積元法的一種并行格式

錢 皓, 吳 丹, 呂俊良

(吉林大學 數學學院, 長春 130012)

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

1 標準有限體積法

1.1 模型問題

考慮擴散問題

(1)

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

1.2 網格剖分與函數空間

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

圖2 頂點P0所在的對偶單元Fig.2 Dual unit of vertex P0

對于給定的t∈[0,T], 取試探函數空間為

其中P1(K)表示在單元K上的一次多項式集合,uh(X,t)|K由K的3個頂點在t時刻的函數值唯一確定.取檢驗函數空間Vh為相應于對偶剖分的分片常數空間, 即

其中P0(K)表示在單元K*上的常數多項式集合.事實上, 對于每個vh∈Vh, 均可唯一地表示為

1.3 有限體積元離散格式

(2)

成立, 其中

(3)

成立, 其中

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

2.1 區域分解

圖3 分成4個并行子區域時的節點Fig.3 Nodes when domain is divided into 4 parallel sub domains

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

2.2 并行算法

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

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

3 數值實驗

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

定義格式的收斂階為

(4)

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

(5)

其中Np為并行程序所用的核數,T1為算法串行花費的時間,Tp為使用p核并行算法花費的時間.令τ=T/N, 則網比μ=τ/h2.

本文采用Fortran編譯語言, 并使用MPICH2的函數庫, 通過MPI通訊方式實現程序的并行.數值實驗在北京超級云計算中心M分區(BSCC-M)上進行, 本地環境為Windows10, 運行環境為Linux-CentOS7, 單個計算節點的處理器為512 GB內存, 128核的AMD EPYC CPU.為達到并行的效果, 本文對每個子區域都采用單獨一個核進行運算, 因此程序所需的核數等于子區域的數量.在求解程序中的線性方程組時, 本文采用GMRES(m)迭代算法, 迭代終止值為10-16.

例1各向同性擴散問題.

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

圖4 攝動三角形網格ⅠFig.4 Distorted trianglar mesh Ⅰ

圖5 Sine三角形網格Fig.5 Sine trianglar mesh

圖6 攝動網格Ⅰ上的計算精度Fig.6 Computational accuracy on distorted mesh Ⅰ

圖7 Sine三角形網格上的計算精度Fig.7 Computational accuracy on Sine trianglar mesh

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

例2各向異性擴散問題.

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

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

圖8 當α=10-6時的計算精度Fig.8 Computational accuracy when α=10-6

圖9 當α=102時的計算精度Fig.9 Computational accuracy when α=102

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

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

例3間斷的擴散系數問題.

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

方程的精確解為

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

圖10 攝動網格ⅡFig.10 Distorted mesh Ⅱ

圖11 梯形網格Fig.11 Trapezoidal mesh

圖12 攝動網格Ⅱ上的計算精度Fig.12 Computational accuracy on distorted mesh Ⅱ

圖13 梯形網格上的計算精度Fig.13 Computational accuracy on trapezoidal mesh

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

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