馬文鵬,李瑞瑩,袁 武,高凌云,翟環(huán)欣

(1. 信陽師范大學(xué) 計算機與信息技術(shù)學(xué)院, 河南 信陽 464000;2. 中國科學(xué)院 計算機網(wǎng)絡(luò)信息中心, 北京 100083; 3. 中國科學(xué)院大學(xué), 北京 100049)

0 引言

在工程計算領(lǐng)域,空氣動力、航空航天、熱能動力、天文物理等復(fù)雜流動問題大都可以采用計算流體力學(xué)(Computational Fluid Dynamics,CFD)的方法進行數(shù)值模擬,有效解決傳統(tǒng)流體力學(xué)實驗帶來的巨大功耗問題。OpenFOAM(Open Source Field Operation and Manipulation,OpenFOAM)[1]是一個完全由C++編寫的CFD軟件,其前后處理功能強大、富含多種物理模型且版本更新快速,被廣泛應(yīng)用于各個工程領(lǐng)域。尤其近年來高性能計算(High Performance Computing,HPC)的迅速發(fā)展,為解決大規(guī)模工程實際問題提供了硬件支撐和有效的可編程實現(xiàn)手段,利用HPC強大的并行算力來提高OpenFOAM的性能已經(jīng)成為主流趨勢,愈來愈多的研究人員開始使用異構(gòu)眾核架構(gòu)來改進CFD傳統(tǒng)的數(shù)值計算方法。

VAN PHUC等[2]基于OpenFOAM對風(fēng)場和自由面流場的模擬進行了多項研究,并使用日本超級計算機“京”成功生成了高達千億單元網(wǎng)格的大規(guī)模瞬態(tài)CFD模擬,驗證了千億單元網(wǎng)格和數(shù)十萬MPI(Message Passing Interface,MPI)的并行潛在能力。LIN等[3]針對OpenFOAM中的線性求解器提出了一種改進的預(yù)處理共軛梯度算法(Preconditioned Conjugate Gradient,PCG)和基于對角線的不完整Cholesky預(yù)處理器方法,在神威超算平臺的單核組上分別實現(xiàn)了8.9倍和3.1倍的速度提升。REN等[4]提出3種減少通信操作的方案對PCG算法進行優(yōu)化,基于超級計算平臺天河二號驗證了OpenFOAM的通信優(yōu)化性能,使大規(guī)模線性方程組迭代求解的效率顯著提升。ROJEK等[5]提出一種與人工智能模塊集成的穩(wěn)態(tài)CFD數(shù)值模擬方法,基于OpenFOAM工具箱MixIT,實現(xiàn)了卷積神經(jīng)網(wǎng)絡(luò)監(jiān)督學(xué)習(xí)算法在GPU異構(gòu)計算平臺上10倍左右的加速。吳家明等[6]基于PETSc(Portable, Extensible Toolkit for Scientific Computation,PETSc)改進線性方程組求解模塊,并將其集成到OpenFOAM框架中,驗證了基于新PETSc框架的線性方程組在百萬級大網(wǎng)格規(guī)模下的并行效率。孫士麗等[7]基于OpenFOAM針對三體船連接橋自由落體入水砰擊模型進行數(shù)值模擬,驗證了數(shù)值計算方法的正確性以及三體船連接橋入水砰擊壓力特性和范圍,為三體船結(jié)構(gòu)強度評估與結(jié)構(gòu)設(shè)計提供了數(shù)值基礎(chǔ)。HE等[8]提出一種面向?qū)ο蟮妮o助應(yīng)用框架,分別對OpenFOAM的原始求解器和湍流模型實現(xiàn)了高效優(yōu)化輔助設(shè)計。

近年來,除了研究OpenFOAM的數(shù)值模擬計算方法外,針對其在異構(gòu)架構(gòu)上的性能擴展研究也日益增多,但異構(gòu)架構(gòu)上的移植工作多著眼于多節(jié)點MPI通信集合的性能提升。有些學(xué)者將求解器中線性解法器模塊作為重點,研究其在異構(gòu)架構(gòu)上的大規(guī)模并行效率和優(yōu)化,或是設(shè)計應(yīng)用接口模型為其在異構(gòu)架構(gòu)上的移植擴展提供輔助分析。目前,一個完整的求解器框架在CPU+GPU加速器異構(gòu)計算平臺上的并行方案仍然存在一些挑戰(zhàn)。主要表現(xiàn)在:一是OpenFOAM完全由C++語言開發(fā),且為了此軟件開發(fā)和維護的靈活性,不斷采用C++的最新且高級的語言特性,這會帶來異構(gòu)平臺上編譯器的兼容性問題;二是CPU+GPU的異構(gòu)軟硬件生態(tài)成多樣的發(fā)展趨勢,如NVIDIA的GPU硬件平臺搭配計算統(tǒng)一設(shè)備架構(gòu)(Compute Unified Device Architecture,CUDA)[9]、AMD的GPU搭配ROCm(Radeon Open Compute,ROCm)[10],給同一求解器或算法在不同異構(gòu)平臺上的移植和適配帶來額外的成本;三是C++中封裝層次較深,底層數(shù)據(jù)不可見,給異構(gòu)移植時一些面向底層數(shù)據(jù)的優(yōu)化帶來困難。如應(yīng)用程序使用Field類表示CFD中各種物理場,對應(yīng)于面向過程語言中的一維或高維數(shù)組,這些數(shù)據(jù)結(jié)構(gòu)被封裝在底層,無法輕易進行結(jié)構(gòu)數(shù)組(Array of Structures,AoS)、數(shù)組結(jié)構(gòu)(Structure of Arrays,SoA)[11]方面的訪存優(yōu)化。另外,軟件的迭代和更新,也給某些模塊的移植提供了思路。如在2020年6月發(fā)布的OpenFOAM-v2006版本[12],在求解線性系統(tǒng)方面對接了PETSc軟件,使得PETSc中豐富的基于Krylov子空間[3]的線性解法器可以用于OpenFOAM的求解中,其中矩陣數(shù)據(jù)類型和PETSc數(shù)據(jù)類型的轉(zhuǎn)化接口也為其外接異構(gòu)線性求解器提供了借鑒。

對于上述問題,本文重點針對OpenFOAM的Field類和lduMatrix類設(shè)計了適應(yīng)異構(gòu)平臺的、可靈活操作底層數(shù)據(jù)的結(jié)構(gòu),并通過數(shù)據(jù)轉(zhuǎn)換接口實現(xiàn)兩種數(shù)據(jù)的轉(zhuǎn)換;采用可移植異構(gòu)計算接口HIP(Heterogeneous-Compute Interface for Portability)[10]完整實現(xiàn)了不可壓求解器icoFoam在兩類主流異構(gòu)平臺(NVIDIA和AMD)上的移植和性能分析;針對求解過程中產(chǎn)生的線性方程組,在異構(gòu)平臺上設(shè)計了LDU(Lower Diagonal Upper)格式到CSR(Compressed Sparse Row)格式[13]的轉(zhuǎn)換方案,并耦合hipSPARSE和hipBLAS[14]實現(xiàn)了基于Krylov子空間方法的線性求解器。

1 控制方程和PISO算法介紹

icoFoam是OpenFOAM提供的用來求解不可壓縮層流的求解器,其對應(yīng)的控制方程為:

?·U=0,

(1)

(2)

式中:U是速度矢量,p是壓力,v是黏度系數(shù)。式(1)是連續(xù)性方程,式(2)是動量方程。

icoFoam在非結(jié)構(gòu)網(wǎng)格上采用有限體積方法離散控制方程,使用PISO(Pressure Implicit with Splitting Operator,PISO)算法[15]在每個時間步內(nèi)實現(xiàn)壓力和速度的多次相互修正來獲得下一個時間步的壓力和速度。具體的流程如圖1所示。

圖1 PISO算法流程圖Fig. 1 Flow chart of PISO algorithm

在當(dāng)前時間步內(nèi),將t0時刻流場的壓力初值p0和速度初值U0代入公式(2)中,PISO循環(huán)開始。首先數(shù)值求解動量方程得到預(yù)測速度Ur,然后使用壓力來表示速度,將速度代入連續(xù)性方程即公式(1)中,得到壓力泊松方程,由預(yù)測速度Ur求解該方程,得到壓力pλ,繼而根據(jù)新的壓力pλ對速度進行修正,得到修正后的速度Uλ。如果在當(dāng)前時間步內(nèi)未達到預(yù)期的迭代次數(shù),則將Uλ當(dāng)作新的速度值Ur作為下一次迭代的初始速度,當(dāng)滿足迭代次數(shù)時,t0時刻的PISO循環(huán)結(jié)束,跳出當(dāng)前時間步,進入下一時刻繼續(xù)求解壓力和速度的耦合,直至最后時刻,輸出最終壓力pt和速度Ut。

2 并行方案設(shè)計

為了兼容主流NVIDIA和AMD的異構(gòu)計算平臺,圖2展示了本文利用可移植的異構(gòu)并行編程模型HIP設(shè)計icoFoam異構(gòu)求解器的整體方案。

圖2 整體方案設(shè)計Fig. 2 Overall scheme design

從數(shù)據(jù)操作的角度出發(fā),將icoFoam的計算分為兩大類:一類是關(guān)于物理場的操作,對應(yīng)于Field類,如插值、通量計算、邊界條件以及關(guān)于某個場數(shù)據(jù)的線性操作;另外一類是關(guān)于矩陣的操作,對應(yīng)于lduMatrix類,如在矩陣組裝成功后進行線性系統(tǒng)的求解,包括矩陣向量乘等?？紤]到C++類對底層操作數(shù)據(jù)的封裝性,通過設(shè)計對應(yīng)的包含CPU端和GPU端數(shù)據(jù)的C結(jié)構(gòu),用于底層數(shù)據(jù)的優(yōu)化。同時,設(shè)計了數(shù)據(jù)結(jié)構(gòu)的轉(zhuǎn)換接口,以保證異構(gòu)計算之外的其他部分,如網(wǎng)格和流場數(shù)據(jù)的IO、流場數(shù)據(jù)的初始化工作、非結(jié)構(gòu)網(wǎng)格拓撲相關(guān)計算操作均保留OpenFOAM的原始特性。

OpenFOAM-v2006版本首次利用其強大的二次開發(fā)功能和動態(tài)鏈接庫機制對接了PETSc,使得外部豐富高效的Krylov子空間算法能夠提供線性解法器。本文充分利用這一機制以及hipSPARSE支持多個后端異構(gòu)平臺的特性,設(shè)計了基于hipSPARSE的Krylov子空間方法的解法器,通過在數(shù)學(xué)庫函數(shù)上層封裝一層hipSPARSE的方式,使編譯器hipcc在編譯調(diào)用的數(shù)學(xué)庫程序時,能夠針對一套HIP程序在不同異構(gòu)平臺上鏈接對應(yīng)的SPARSE庫,如cuSPARSE或rocSPARSE[14]。

2.1 非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)結(jié)構(gòu)轉(zhuǎn)換

OpenFOAM中包含許多重要的類模板,考慮到類模板對底層操作數(shù)據(jù)的封裝性,針對icoFoam求解器的Field和lduMatrix兩種核心類,將程序的核心代碼利用HIP編程模型重新實現(xiàn),其余部分仍以C++為基礎(chǔ)。

關(guān)于物理場的操作,其核心是對Field類的操作,如梯度計算、速度通量計算、邊界條件以及關(guān)于某個場數(shù)據(jù)的線性操作,均是對Field類的操作。另外一類是關(guān)于矩陣的操作,對應(yīng)于LduMatrix類,如對矩陣的組裝、矩陣存儲格式的轉(zhuǎn)換以及矩陣組裝成功后對線性系統(tǒng)的求解,包括矩陣向量乘等。Field和LduMatrix類都繼承于List類,相當(dāng)于對數(shù)組的操作,經(jīng)過多層次的封裝和繼承,底層數(shù)據(jù)對用戶均不可見,需要提取最底層的核心數(shù)據(jù)。因此基于HIP編程模型設(shè)計了數(shù)組結(jié)構(gòu)的存儲格式和數(shù)據(jù)接口用以底層數(shù)據(jù)的轉(zhuǎn)換。如圖3所示,以體場volFields類為例,設(shè)計同時包含CPU端和GPU端數(shù)據(jù)的結(jié)構(gòu)以及數(shù)據(jù)轉(zhuǎn)換接口,用于底層數(shù)據(jù)的轉(zhuǎn)換,從而提高求解器的訪存優(yōu)化性能。

圖3 Field類數(shù)據(jù)格式轉(zhuǎn)換Fig. 3 Field class data format conversion

2.2 速度通量并行求解算法

Field類的核心工作是通過通量計算、插值計算和通量修正3個階段實現(xiàn)對物理場phiHbyA的求解,最終返回1個面場,3個函數(shù)的共性是將網(wǎng)格格心數(shù)據(jù)插值為網(wǎng)格面數(shù)據(jù)。

如圖4所示,OpenFOAM中通常將體場變量存儲在網(wǎng)格格心,計算之前采用中心差分法將所有網(wǎng)格面對應(yīng)的兩個控制單元owner(簡稱O)和neighbour(簡稱N)的體場變量插值到面上。數(shù)據(jù)由體場到面場的操作中,每個網(wǎng)格面是獨立的,由于網(wǎng)格的邊界場數(shù)據(jù)不用進行插值操作,所以在計算過程中將內(nèi)場和邊界場分開處理,并組織與內(nèi)面數(shù)量相等的線程數(shù)同時計算內(nèi)場和邊界場。

圖4 網(wǎng)格格心數(shù)據(jù)的插值Fig. 4 Interpolation of grid cell center data

在線程組織方面,以網(wǎng)格內(nèi)部面為單位來分配線程任務(wù),假設(shè)網(wǎng)格內(nèi)面的數(shù)量共有nIFaces個,那么有nIFaces個線程以blockSize的數(shù)量來分塊,線程塊數(shù)為「nIFaces/blockSize?個。在算法1中,由Vol表示體場,Sur表示面場,每個網(wǎng)格面上都有一個線程來實現(xiàn)內(nèi)場和邊界場上的體場計算面場的操作。

算法1:體場計算面場算法和并行實現(xiàn)

輸入:Vol_U, Sur_Sf, 插值系數(shù)w, 所有網(wǎng)格面的O,N

輸出:內(nèi)場Sur_phii、邊界場Sur_phib

1. tid← blockDim.x * blockIdx.x + threadIdx.x

2. if tid

3. Sur_Ui[tid] ← w * (Vol_Ui[O[tid]] + Vol_Ui[N[tid]])

4. Sur_phii[tid] ← Sur_Sfi[tid]Sur_Ui[tid]

5. endif

6. if tid

7. Sur_phib[tid] ← Sur_Sfb[tid] · Vol_Ub[tid]

8. endif

9. return Sur_phii, Sur_phib

如圖5所示,Field類的矢量場在CPU中均采用典型的AoS模式來存儲其空間上相鄰的數(shù)據(jù)分量,例如x、y和z。這種順序訪存模式在GPU中導(dǎo)致訪問單個類型變量時至少產(chǎn)生3倍的時間開銷。因此計算過程中使用SoA結(jié)構(gòu)將矢量場分量每一字段中所有的值存儲到各自的數(shù)組中,使線程能夠順序且連續(xù)地訪問數(shù)據(jù),更高效地利用GPU的全局內(nèi)存。

圖5 AoS和SoA模式的內(nèi)存布局Fig. 5 Memory layout for AoS and SoA patterns

2.3 矩陣組裝并行算法與實現(xiàn)

如圖6所示,OpenFOAM中的矩陣采用LDU的格式來存儲線性方程組的矩陣系數(shù):lower(簡稱l)、diag(簡稱d)和upper(簡稱u)分別對應(yīng)下三角、對角線和上三角。其中,將網(wǎng)格內(nèi)面對應(yīng)的單元O和N代替行列索引upperAddr和lowerAddr作為矩陣上下2個三角元素的索引值共用的2個索引向量[8]。

圖6 OpenFOAM中的矩陣存儲格式Fig. 6 Matrix storage format in OpenFoam

icoFoam求解器中,速度場和壓力場對應(yīng)的線性方程組信息分別存儲在速度方程UEqn和壓力方程PEqn中,離散后都會生成一個LDU格式的矩陣用于求解線性方程組。針對與速度方程UEqn相關(guān)的矩陣和源項組裝,將其分為一階時間導(dǎo)數(shù)項、對流項以及擴散項3個階段。

一階時間導(dǎo)數(shù)項中,基于網(wǎng)格單元求解對角線系數(shù)和源項;對流項和擴散項中,基于網(wǎng)格內(nèi)面計算矩陣系數(shù),基于網(wǎng)格單元計算邊界貢獻。因此依據(jù)網(wǎng)格特性來組織線程,例如以nIFaces個線程、「nIFaces/blockSize?個線程塊同時計算矩陣系數(shù)和邊界貢獻。在計算對角線系數(shù)時,為了避免上、下三角系數(shù)取到同一個對角線的值而引起訪問沖突,利用原子操作使讀—計算—寫能夠連貫執(zhí)行,確保只有當(dāng)前線程訪問計算對角線系數(shù)所需資源,不受其他線程干擾。

針對UEqn,設(shè)計算法2來組裝矩陣系數(shù)和源項,pEqn中僅包含擴散項,UEqn中擴散項的組裝方法同樣適用于PEqn。

算法2:組裝矩陣系數(shù)算法和并行實現(xiàn)

輸入:1.網(wǎng)格體積C, 系數(shù)r, w, 速度Ut-1, 通量phi

輸出:矩陣系數(shù)l, d, u、源項S、邊界貢獻Ic、Bc

1. tid← blockDim.x * blockIdx.x + threadIdx.x

2. if tid

3. d[tid] ← r * C[tid]

4. S[tid] ← r * Ut-1[tid] * C[tid]

5. endif

6. if tid

7. l[tid] ← w * phi[tid]

8. u[tid] ← l[tid] + phi[tid]

9. atomicAdd(&d[N[tid]],(-l[N[tid]])

10. atomicAdd(&d[O[tid]],(-u[O[tid]])

11. endif

12. if tid

13. Ic[tid]← 0.0

14. Bc[tid] ← phi[tid] * Ut-1[tid]

15. endif

16. return l, d, u, S, Ic, Bc

2.4 線性求解器并行算法與實現(xiàn)

OpenFOAM中關(guān)于流場的核心運算是通過網(wǎng)格面來組織的,如圖6所示,其計算流程中線性系統(tǒng)的系數(shù)矩陣的存儲格式也與面單元的拓撲結(jié)構(gòu)一一對應(yīng)。LDU格式的矩陣表示與以面為單位組織計算的特點非常契合,但是由于矩陣數(shù)據(jù)被拆分存儲,也給在加速平臺上矩陣按照整行(整列)為單位的操作(如稀疏矩陣向量乘操作)帶來了不便,同時整行數(shù)據(jù)的分散存儲阻礙了數(shù)據(jù)局部性、訪存融合的優(yōu)化。在加速平臺上,提出了一種LDU到CSR格式的轉(zhuǎn)換算法,且基于hipSPARSE庫,針對速度方程和壓力方程,耦合BiCGStab和PCG兩種Krylov子空間方法作為線性系統(tǒng)的解法器。

算法3描述了該轉(zhuǎn)換算法的過程,所有操作均在加速器上完成。

算法3：LDU格式到CSR格式的轉(zhuǎn)換算法和并行實現(xiàn)

輸入:LDU格式的矩陣上、下三角、對角數(shù)據(jù):

1. u, l([0:nIFaces-1]), d([0:nCells-1])

2. 所有網(wǎng)格面的O和N, ([0:nIFaces-1])

輸出:CSR格式的矩陣:csr_r, csr_c, csr_v

1. 將LDU格式合并成非排序的coo格式:

tid← blockDim.x * blockIdx.x + threadIdx.x

if tid

(coo_r, coo_c, coo_v)[tid] ← (N, O, l)[tid]

endif

if tid

(coo_r, coo_c, coo_v)[nIFaces+tid] ← (tid, tid, d[tid])

endif

if tid

nIFC ← nIFaces + nCells

(coo_r, coo_c, coo_v)[tid + nIFC] ← (O, N ,u)[tid]

endif

2.將Coo格式按照行壓縮排序,得到置換P:

P ← hipSparseXcooSortByRow(coo_r, coo_c)

3.將置換應(yīng)用到coo_v得到排序后的矩陣:

sorted_coo_v ← hipSparseDgthr(P, coo_v)

4.將coo_r按照行壓縮存儲得到csr_r:

csr_r ← hipsparseXcoo2csr(coo_r)

5.csr_c ← coo_c; csr_v ← sorted_coo_v

6.return csr_r, csr_c, csr_v

算法3中,首先由LDU格式的存儲規(guī)則,可以獲得矩陣每個元素的行號和列號,發(fā)起一個kernel函數(shù),將矩陣的下三角、上三角和對角元素的行號、列號及元素的值合并到一個未排序的COO格式存儲的矩陣中;隨后的任務(wù)就轉(zhuǎn)化為將COO格式的矩陣按照行壓縮進行排序,基于hipSPARSE庫函數(shù),此部分分為兩個步驟:第一步是通過hipsparseXcooSortByRow獲取一個從未排序到按行壓縮排序的置換;第二步將置換應(yīng)用到原始的COO格式的矩陣的元素中,得到轉(zhuǎn)換過后的元素;最后,將COO格式表示的行信息coo_r壓縮得到CSR格式表示的行信息csr_r。該算法最耗時的部分是得到置換數(shù)組P的過程,屬于訪存密集型操作,但是該步驟在針對由速度或者壓力方程得到的系數(shù)矩陣時,只需要執(zhí)行一次,因為在整個時間迭代過程中,每個時間步對應(yīng)的系數(shù)矩陣非零元的位置是保持不變的。在得到CSR格式存儲的矩陣之后,基于hipSPARSE和hipBLAS庫,實現(xiàn)了Krylov子空間的PCG和BiCGStab兩個線性解法器,用于求解壓力和速度對應(yīng)的線性系統(tǒng)。

3 測試環(huán)境及算例

3.1 測試環(huán)境

實驗是在分別包含AMD gfx906 GPU和NVIDIA TeslaV100 GPU的兩個異構(gòu)計算節(jié)點上進行的。其中,節(jié)點1采用的是2.5 GHz的CPU主頻處理器,共有128 GB內(nèi)存,并配置4塊16 GB gfx906 GPU加速器,rocm版本為4.3(dtk 21.10)。節(jié)點2采用的是2.4 GHz的Intel(R) Xeon(R) E5-2640 v4 CPU主頻處理器,共有128 GB內(nèi)存,并配置4塊16 GB的TeslaV100顯卡,CUDA版本為10.0。

3.2 算例設(shè)置

實驗使用OpenFOAM-v2006版本的不可壓縮層流頂蓋驅(qū)動方腔流(Lid-driven Cavity Flow,Cavity)[1]為測試算例,如圖7所示。該算例通過將z軸方向的邊界條件設(shè)置為empty類型來實現(xiàn)二維的計算。

圖7 Cavity算例模型Fig. 7 Example model of Cavity

針對Cavity算例在z軸方向進行擴展,將front和back的邊界條件從原始的empty類型設(shè)置為symmetryPlane類型,movingWall與fixedWalls的邊界條件保持不變,其網(wǎng)格控制文件blockMesh的設(shè)置如表1所示。

表1 Cavity中blockMesh設(shè)置Tab. 1 BlockMesh setting of Cavity

在使用Cavity算例測試時,同時將x、y和z軸方向擴展為20至40,60,…,100,設(shè)置算例運行時間為1×10-4s,時間步長為1×10-6s,共迭代100個時間步。在測試過程中關(guān)閉I/O操作,避免流場信息的頻繁寫入、寫出影響測試結(jié)果的穩(wěn)定性。

4 實驗結(jié)果與分析

采用HIP實現(xiàn)的icoFoam的同一套不可壓異構(gòu)求解器,分別在ROCm和CUDA兩個平臺上進行測試。與2.2、2.3、2.4小節(jié)對應(yīng),在測試過程中分別統(tǒng)計這三部分在加速器上的時間,并與icoFoam在CPU上的時間進行對比分析。CPU版本的串行程序在編譯過程中使用了-O3優(yōu)化選項。為了保證對比實驗的有效性,使用CPU最優(yōu)的測試時間和GPU測試時間計算加速比,且為了避免測試的偶然性,統(tǒng)計測試時間使用連續(xù)3次在GPU上測試時間的平均值。

表2給出了不同網(wǎng)格規(guī)模下分別在ROCm和CUDA平臺上面場相關(guān)操作的時間(tROCm、tCUDA)和對應(yīng)的CPU上時間(tCPU)的對比。由表2的實驗數(shù)據(jù)可以看出,經(jīng)過100個時間步的迭代后,不同網(wǎng)格規(guī)模下的ROCm平臺計算面場的算法加速效果在3.7～66.0倍。CUDA平臺加速效果在5.4～87.2倍,隨著網(wǎng)格規(guī)模的增大,加速比呈現(xiàn)遞增趨勢。說明與icoFoam串行計算相比,GPU的并行計算方法具有良好的加速效果,時間計算成本顯著降低。

表2 面場串行和并行計算時間對比Tab. 2 Comparison of serial and parallel computation time of surface field

表3給出了不同網(wǎng)格規(guī)模下組裝系數(shù)矩陣在各平臺上的對比。根據(jù)表3的實驗結(jié)果,發(fā)現(xiàn)GPU的矩陣組裝并行方案加速效果明顯,在ROCm計算平臺上的加速比為2.6～43.3倍,CUDA計算平臺上的加速比為5.4～66.4倍。值得注意的是,算法2在進行LDU格式矩陣系數(shù)組裝時,使用了雙精度原子操作,盡管如此,相比于CPU的時間,原子操作的實現(xiàn)依然有較大的并行優(yōu)勢,這主要得益于原子操作在兩個平臺上的高效實現(xiàn)。

表3 矩陣組裝串行和并行計算時間對比Tab. 3 Comparison of serial and parallel computation time of matrix assembly

表4給出了算法3中從LDU格式到CSR格式的轉(zhuǎn)換時間,由于在時間迭代過程中,兩種格式在非零元位置的對應(yīng)關(guān)系上保持不變,因此,算法3中得到置換的操作僅需要執(zhí)行一次即可。表5和表6給出了分別使用PCG和BiCGStab求解一次壓力和某個速度分量相關(guān)的線性系統(tǒng)的時間,雖然矩陣格式轉(zhuǎn)換時間(20 ms左右)相比于求解一次線性系統(tǒng)的時間并不占優(yōu)勢,但是由于轉(zhuǎn)換在整個求解流程中僅需一次,而每個時間步都需要求解一個線性系統(tǒng),因此,矩陣轉(zhuǎn)換時間相對于整個計算流程時間,如計算100個時間步,在1%左右,并不會影響整體計算效率。

表4 LDU格式到CSR格式的轉(zhuǎn)換時間統(tǒng)計Tab. 4 Conversion time statistics of PCG algorithm LDU format to CSR format

表5和表6的實驗結(jié)果表明,基于hipSPARSE和hipBLAS實現(xiàn)的兩個Krylov子空間的線性解法器能夠同時支持ROCm和CUDA平臺,并能獲得18倍和37倍的加速效果。由于求解線性系統(tǒng)占整個計算的大部分時間,因此對于icoFoam,整體在異構(gòu)平臺上的加速分別為20.5倍和38.4倍。

5 結(jié)論

本文面向CPU+GPU眾核處理器使用可移植性跨平臺接口(HIP)實現(xiàn)了不可壓縮求解器icoFoam在兩類主流異構(gòu)平臺NVIDIA和AMD上的移植,基于異構(gòu)眾核加速器測試了不同算例規(guī)模下的求解器并行算法的計算時間。經(jīng)過與OpenFOAM串行程序的對比分析,基于兩類異構(gòu)計算平臺的整個求解器并行算法的設(shè)計,能夠大幅度提升OpenFOAM的計算效率,大規(guī)模并行計算能力的提升分別高達20.5倍和38.4倍。