基于MRT-LBM方法的大規(guī)?？蓴U(kuò)展并行計算研究

劉智翔　方　勇　宋安平　徐　磊　王曉偉2,　周麗萍2,　張　武2,

1(上海大學(xué)通信與信息工程學(xué)院　上?！?00444)2(上海大學(xué)高性能計算中心　上?！?00444)3(上海大學(xué)計算機(jī)工程與科學(xué)學(xué)院　上?！?00444)(zxliu@shu.edu.cn)

基于MRT-LBM方法的大規(guī)?？蓴U(kuò)展并行計算研究

劉智翔1,2方勇1宋安平3徐磊3王曉偉2,3周麗萍2,3張武2,3

1(上海大學(xué)通信與信息工程學(xué)院上海200444)2(上海大學(xué)高性能計算中心上海200444)3(上海大學(xué)計算機(jī)工程與科學(xué)學(xué)院上海200444)(zxliu@shu.edu.cn)

摘要在大規(guī)模三維復(fù)雜流動的數(shù)值模擬中，針對具有良好數(shù)值穩(wěn)定性的多弛豫時間模型格子Boltzmann方法(MRT-LBM)，并結(jié)合大渦模擬湍流模型和曲面邊界插值格式，分析了在D3Q19離散速度模型下的網(wǎng)格生成、流場信息初始化和迭代計算3部分的可并行性.采用MPI編程模型，從分布式集群的特點(diǎn)和計算量負(fù)載均衡的角度出發(fā)，分別提出了適合于大規(guī)模分布式集群的網(wǎng)格生成、流場信息初始化和迭代計算的并行算法.該并行算法也能有效適用于D3Q15和D3Q27離散速度模型.通過在國產(chǎn)神威藍(lán)光超級計算機(jī)上的測試，分別針對求解問題總體計算規(guī)模固定和保持每個計算核中計算量一致的2種情況的并行性能分析，驗(yàn)證了該并行算法在十萬計算核的量級下仍具有良好的加速比和可擴(kuò)展性.

關(guān)鍵詞大規(guī)模并行計算；可擴(kuò)展；負(fù)載均衡；格子Boltzmann方法；多松弛時間模型；大渦模擬

復(fù)雜流體運(yùn)動的數(shù)值模擬一直是大規(guī)?？茖W(xué)與工程計算中最重要且具有挑戰(zhàn)性的研究領(lǐng)域之一.與基于宏觀Navier-Stokes方程的計算流體動力學(xué)和基于微觀分子動力學(xué)的流體力學(xué)模擬方法不同，格子Boltzmann方法(lattice Boltzmann method, LBM)從介觀尺度來描述流體系統(tǒng)[1-2]，由于Boltzmann方程自身本質(zhì)的運(yùn)動學(xué)特性，以及可以根據(jù)經(jīng)典的Chapman-Enskog展開從LBM得到Navier-Stokes方程，使得LBM比基于連續(xù)介質(zhì)假設(shè)的Navier-Stokes方程包含了更多的物理內(nèi)涵.同時LBM可以看成是Navier-Stokes方程差分法逼近的一種無條件穩(wěn)定的格式[2].在LBM中沒有基于連續(xù)介質(zhì)的假設(shè)，而是把流體看成是由許多只有質(zhì)量沒有體積的微粒所組成,這些微?？梢韵蚩臻g若干個方向任意運(yùn)動.通過其質(zhì)量、動量、能量守恒原理，建立表征質(zhì)點(diǎn)在給定的時刻位于空間某一個位置附近的概率密度函數(shù)，再通過統(tǒng)計的方法來獲得質(zhì)點(diǎn)微粒的概率密度分布函數(shù)與宏觀運(yùn)動參數(shù)間的關(guān)系.利用這種簡化的運(yùn)動學(xué)方程形式，既可以避免求解含有復(fù)雜碰撞項的完整積分-微分Boltzmann方程，又可避免對流場中的每個分子進(jìn)行大規(guī)模計算.由于LBM具有算法簡單、計算效率高、并行性好以及能夠方便處理復(fù)雜邊界條件等優(yōu)點(diǎn)，所以其非常適合處理復(fù)雜幾何形狀和復(fù)雜流動大規(guī)模計算問題[3].

目前，常用的LBM方法中主要采用3種碰撞模型：單弛豫時間(single-relaxation-time, SRT)[4-5]、多弛豫時間(multiple-relaxation-time, MRT)碰撞模型[6-8]和熵模型(entropic)[9].SRT碰撞模型只有一個弛豫時間可調(diào)，而MRT模型是將碰撞過程通過線性變換轉(zhuǎn)換到矩空間中進(jìn)行，且更好地考慮了實(shí)際流動中的各向異性和計算條件的參數(shù)優(yōu)化，具有多個可調(diào)的弛豫時間.文獻(xiàn)[6]針對方腔流問題，從數(shù)值穩(wěn)定性和精度等角度對這3種碰撞模型進(jìn)行了詳細(xì)分析，分析出MRT碰撞模型具有更好的數(shù)值穩(wěn)定性和精度.而在針對湍流問題的數(shù)值模擬中，由于問題本身的復(fù)雜性和具有高雷諾數(shù)，直接使用LBM進(jìn)行數(shù)值模擬無法準(zhǔn)確進(jìn)行計算，常需要將大渦模擬湍流模型引入到計算中[10].Krafczyk等人[11]將Smagorinsky渦粘性模型與D3Q15的MRT-LBM進(jìn)行耦合計算，對槽道中物體繞流進(jìn)行了大渦模擬研究.Menon等人[12]將局部動力亞格子應(yīng)力模型引入SRT-LBM對三維合成射流和自由射流進(jìn)行了數(shù)值模擬.Leila等人[13]對LBM在湍流數(shù)值模擬中的應(yīng)用進(jìn)行了詳細(xì)總結(jié)，證明了結(jié)合大渦模擬湍流模型的MRT-LBM進(jìn)行湍流的模擬具有更好的數(shù)值穩(wěn)定性.

對于復(fù)雜流動問題的計算，LBM往往需要大規(guī)模的笛卡兒網(wǎng)格才能準(zhǔn)確進(jìn)行數(shù)值模擬，這時高性能計算將是一種較好的解決辦法.目前，高性能計算集群大多采用分布式結(jié)構(gòu)，在不同的計算節(jié)點(diǎn)間采用消息傳遞接口(message passing interface, MPI)進(jìn)行通信.而千萬億次的超級計算機(jī)(天河、曙光星云、神威)通常都具有數(shù)十萬、百萬的計算核，充分發(fā)揮大規(guī)模計算核的優(yōu)勢來加速問題的求解，需要研究具有良好可擴(kuò)展性的高效并行程序[14].考慮到LBM適合于大規(guī)模計算的優(yōu)點(diǎn)和高性能計算機(jī)的特點(diǎn)，本文將針對MRT-LBM方法，并結(jié)合大渦模擬、曲面邊界插值格式和笛卡兒網(wǎng)格生成，研究能夠求解實(shí)際流動問題、具有高可擴(kuò)展的MRT-LBM并行算法.

1多弛豫時間的格子Boltzmann方法

本節(jié)主要介紹在三維情況下基于D3Q19離散速度模型的多弛豫時間格子Boltzmann方法、大渦模擬湍流模型和曲面邊界的插值格式.

1.1MRT-LBM

格子Boltzmann方法是一種介觀方法，它包括3部分：粒子分布函數(shù)、離散速度模型和分布函數(shù)的演化模型.LBM的分布函數(shù)演化模型[7]可表示為

fi(x+ξiδt,t+δt)-fi(x,t)=Ωi,

(1)

其中，fi是粒子分布函數(shù)，Ωi是碰撞項，δt是時間步長，ξi是第i個方向的離散速度.由于本文針對的是實(shí)際流動問題的數(shù)值模擬，所以僅考慮三維情況.

文中使用D3Q19離散速度模型進(jìn)行描述，其離散速度(如圖1所示)為

(2)

其中,c是格子聲速，其值為c=δxδt=1.碰撞項Ωi通常有3種：單弛豫時間、多弛豫時間和熵模型.在MRT模型中，碰撞是在矩空間中進(jìn)行，則碰撞項可以表示為

Ωi=-{M-1·S·[m(x,t)-meq(x,t)]}i，

(3)

其中，M是變換矩陣，meq是矩空間的平衡態(tài)函數(shù)，其具體值可見文獻(xiàn)[7].S是對角矩陣，其值為

S=diag(0,s1,s2,0,s4,0,s4,0,s4,s9,

s10,s9,s10,s13,s13,s13,s16,s16,s16).

(4)

Fig. 1　The D3Q19 velocity model.圖1　D3Q19離散速度模型

弛豫參數(shù)s9和s13滿足：

(5)

這里ν是粘性系數(shù).其余的弛豫參數(shù)可根據(jù)文獻(xiàn)[8]取為s1=1.19，s2=s10=1.4，s4=1.2，s16=1.98.根據(jù)上面的描述，MRT-LBM的計算過程可以分為2部分：

1) 碰撞過程:

{M-1·S·[m(x,t)-meq(x,t)]}i；

(6)

2) 遷移過程:

(7)

1.2大渦模擬湍流模型

為了準(zhǔn)確模擬高雷諾數(shù)湍流運(yùn)動中的情況，可以在MRT-LBM中引入大渦模擬亞格子模型，文中采用常用的Smagorinsky渦粘模型[15]，其引入的方式是通過改變原來的弛豫時間來實(shí)現(xiàn)(即修改式(5)中的粘性系數(shù))，引入大渦模擬后的有效粘性變?yōu)?/p>

ν=ν0+νt，

(8)

νt=(CsΔ)2|S|,

(9)

這里，Cs是Smagorinsky參數(shù)，|S|為應(yīng)變率張量Sij的模.因此，在MRT-LBM中引入的大渦模擬湍流模型僅需要在每次碰撞時按照式(8)(9)調(diào)整弛豫參數(shù)s9和s13.

1.3曲面插值格式

Fig. 2　Layout of the regularly spaced lattices and curved wall boundary.圖2　網(wǎng)格點(diǎn)和曲面邊界的布局圖

在三維復(fù)雜幾何的數(shù)值模擬中，由于MRT-LBM采用的是統(tǒng)一尺度的笛卡兒網(wǎng)格，網(wǎng)格點(diǎn)無法保證恰好在曲面邊界上，因此需要構(gòu)造一種插值格式來處理曲面邊界上的網(wǎng)格點(diǎn).為了便于描述插值格式，這里給出二維情況下的曲面邊界示意圖，如圖2所示：

在MRT-LBM中，我們先將所有的網(wǎng)格點(diǎn)分為3類：流體網(wǎng)格點(diǎn)xf(fluid node)、曲面邊界網(wǎng)格點(diǎn)xw(boundary wall node)和幾何體內(nèi)部網(wǎng)格點(diǎn)xs(solid node).由于幾何體內(nèi)部的網(wǎng)格點(diǎn)在物體內(nèi)部，所以不需要計算.流體網(wǎng)格點(diǎn)和曲面邊界網(wǎng)格點(diǎn)需要計算.因此，需要在曲面邊界網(wǎng)格點(diǎn)上進(jìn)行插值以獲得D3Q19離散速度方向上的全部未知粒子分布函數(shù).本文采用文獻(xiàn)[16]中的YMS曲面插值格式處理曲面邊界網(wǎng)格點(diǎn)xw，其具體插值格式為

(10)

(11)

2MRT-LBM高可擴(kuò)展并行算法

根據(jù)第1節(jié)中的內(nèi)容，MRT-LBM的計算流程如圖3所示：

Fig. 3　Flowchart of MRT-LBM.圖3　MRT-LBM的計算流程圖

本節(jié)將根據(jù)圖3中的MRT-LBM計算流程圖，從笛卡兒網(wǎng)格生成開始，逐步分析MRT-LBM的并行性，并提出基于MRT-LBM的并行算法.MRT-LBM的求解基本流程主要包括3個部分：網(wǎng)格生成、流場信息初始化和迭代計算.下面我們將分別對這3部分進(jìn)行并行性分析.

2.1網(wǎng)格生成部分的并行

MRT-LBM使用的是統(tǒng)一尺度的笛卡兒網(wǎng)格，如圖4所示.笛卡兒網(wǎng)格的并行生成比較容易實(shí)現(xiàn)，其主要的難點(diǎn)在于從計算負(fù)載均衡的角度出發(fā)如何劃分網(wǎng)格和網(wǎng)格點(diǎn)類型(包括流體網(wǎng)格點(diǎn)、邊界網(wǎng)格點(diǎn)和幾何體內(nèi)部網(wǎng)格點(diǎn))的判斷.

Fig. 4　Schema of three-dimensional grid parallel partition.圖4　三維網(wǎng)格并行劃分示意圖

首先我們給出在三維情況下，網(wǎng)格生成的3個主要步驟：

Step1. 根據(jù)流場范圍和網(wǎng)格尺寸Δx，生成1個初始的三維笛卡兒網(wǎng)格，每個網(wǎng)格點(diǎn)的信息包括：網(wǎng)格點(diǎn)的坐標(biāo)值、圖1中19個速度方向的鄰居編號.

Step2. 根據(jù)復(fù)雜幾何的STL文件格式，判斷笛卡兒網(wǎng)格點(diǎn)在幾何體的內(nèi)外，一般采用射線法進(jìn)行判斷，具體判斷算法可見文獻(xiàn)[17].判斷完成后需要標(biāo)記網(wǎng)格點(diǎn)的類型：流場內(nèi)網(wǎng)格點(diǎn)(fluid)、幾何體內(nèi)部網(wǎng)格點(diǎn)(solid).

Step3. 與幾何體相鄰的流場內(nèi)網(wǎng)格點(diǎn)標(biāo)記為邊界網(wǎng)格點(diǎn)(boundary)，同時根據(jù)式(11)計算邊界網(wǎng)格點(diǎn)沿著19個速度方向到幾何體邊界的距離值q.

根據(jù)Step1～Step3，并結(jié)合MRT-LBM的碰撞遷移過程(式(6)(7))和分布式集群的特點(diǎn)，網(wǎng)格生成部分進(jìn)行并行操作步驟如下：

Step1. 根據(jù)計算進(jìn)程數(shù)np，將流場范圍均等分成np份，記錄每一份的流場范圍[(xi,min,yi,min,zi,min),(xi,max,yi,max,zi,max)](i=1,2,…,np).

Step2. 如圖4所示，為了方便MRT-LBM的計算，在流場范圍均勻劃分時需要預(yù)留1塊緩沖區(qū)，即每個計算進(jìn)程負(fù)責(zé)的流場范圍應(yīng)該為[(xi,min-Δx,yi,min-Δx,zi,min-Δx),(xi,max+Δx,yi,max+Δx,zi,max+Δx)](i=1,2,…,np).

Step3.每個計算進(jìn)程根據(jù)各自負(fù)責(zé)的流場范圍，根據(jù)并行操作Step1～Step3各自獨(dú)立生成MRT-LBM所需的計算網(wǎng)格.

從并行操作Step1～Step3中可以看出，每個計算進(jìn)程在理論上分配的網(wǎng)格總數(shù)是一致的，且每個進(jìn)程各自獨(dú)立進(jìn)行網(wǎng)格生成，無需通信，這樣并行劃分的網(wǎng)格生成的好處是能夠被后續(xù)的MRT-LBM計算直接使用，且每個進(jìn)程的網(wǎng)格數(shù)基本一致，具有較好的負(fù)載均衡.

2.2流場信息初始化部分的并行

笛卡兒網(wǎng)格生成后，需要進(jìn)行的是網(wǎng)格點(diǎn)上流場信息的初始化，包括：宏觀量(密度、速度、壓力)的初始化和微觀量(粒子分布函數(shù))的初始化.宏觀量的初始化主要是根據(jù)初始的宏觀條件進(jìn)行賦值，而微觀量需要根據(jù)已賦值的宏觀量通過平衡態(tài)分布函數(shù)進(jìn)行計算獲得.流場信息初始化的過程都是在每個網(wǎng)格點(diǎn)上獨(dú)立進(jìn)行，相鄰網(wǎng)格點(diǎn)間無需通信.這里采用2.1節(jié)中的已經(jīng)劃分后的網(wǎng)格(每個MPI進(jìn)行對應(yīng)1個網(wǎng)格區(qū)域)，每個MPI計算進(jìn)程內(nèi)部都可以執(zhí)行如下步驟：

Step1. 根據(jù)網(wǎng)格點(diǎn)的類型和初始的宏觀條件給每個網(wǎng)格點(diǎn)賦初始宏觀量的值.

Step2. 根據(jù)網(wǎng)格點(diǎn)上的初始宏觀量,計算每個網(wǎng)格點(diǎn)上微觀量的值.

由于2.1節(jié)中預(yù)留的緩沖區(qū)不用于計算，僅用于MPI的通信操作，所以在流場信息初始化時不需要涉及到該緩沖區(qū).此外，根據(jù)2.1節(jié)已經(jīng)完成了并行任務(wù)劃分，本部分的計算操作僅需要對每個對應(yīng)的MPI計算進(jìn)程執(zhí)行初始化，無需MPI通信.由于每個MPI進(jìn)程的網(wǎng)格數(shù)基本一致，因此初始化部分是滿足負(fù)載基本均衡，且具有良好的可擴(kuò)展性和執(zhí)行效率.

2.3迭代計算部分的并行

這部分是MRT-LBM計算的核心部分，迭代的每一次計算都包括碰撞過程(如式(6)所示)和遷移過程(如式(7)所示)，且在碰撞過程中還引入了大渦模擬湍流模型(如式(8)所示)，在遷移過程中需要使用曲面邊界YMS插值格式(如式(10)所示).在給出具體的并行步驟前，我們先對涉及2.1節(jié)中緩沖區(qū)部分的MPI通信方式進(jìn)行研究.由于在三維情況下的通信與二維情況的通信類似，為了便于描述，本文給出在二維情況下的一維網(wǎng)格劃分和二維網(wǎng)格劃分的通信模型，如圖5所示：

Fig. 5　Data transmission between processes.圖5　進(jìn)程間數(shù)據(jù)傳遞圖

在圖5中灰色部分為2.1節(jié)中所提的緩沖區(qū)部分，它只用于接收相鄰網(wǎng)格塊的數(shù)據(jù).從圖5(a)可以看出，一維網(wǎng)格劃分情況的MPI數(shù)據(jù)傳遞相對簡單，即相鄰的MPI進(jìn)程分別進(jìn)行1次數(shù)據(jù)的傳遞與接收，如果采用非阻塞MPI通信僅需要2次MPI通信就能進(jìn)行1次MRT-LBM的迭代計算.在圖5(b)中，二維網(wǎng)格劃分的MPI通信較為復(fù)雜，通信的次數(shù)也更多，但是該劃分下的MPI通信與D2Q9離散速度模型是一致的.圖5中的2種網(wǎng)格劃分的進(jìn)程間數(shù)據(jù)傳遞都反映出MRT-LBM僅需要局部通信，無需全局通信.而二維網(wǎng)格劃分與一維網(wǎng)格劃分相比，在同樣的問題計算規(guī)模下，二維網(wǎng)格劃分能夠劃分更多的網(wǎng)格塊，可以擴(kuò)大同一個計算問題下所能使用的計算核數(shù).

類似地，針對三維問題中的D3Q19離散速度模型的MPI通信模型，也包括一維、二維和三維網(wǎng)格劃分，每種劃分都僅需局部通信.隨著劃分維數(shù)的增加，MPI通信的復(fù)雜程度也增加，但是所能劃分的網(wǎng)格塊越多，能夠使用更多的計算核求解同一個問題，并加快求解速度.

結(jié)合圖5中的MPI進(jìn)程數(shù)據(jù)傳遞圖和2.1節(jié)、2.2節(jié)中的并行內(nèi)容，迭代計算部分的并行步驟(每個計算進(jìn)程)如下：

Step1. 按照類似于圖5中的進(jìn)程間數(shù)據(jù)傳遞圖，建立進(jìn)程間MPI通信的編號數(shù)組，為后續(xù)的傳遞做準(zhǔn)備；

Step2. 在除緩沖區(qū)和solid以外的網(wǎng)格點(diǎn)中，根據(jù)式(8)計算引入大渦模擬湍流模型的松弛參數(shù)s9和s13；

Step5. 根據(jù)式(7)進(jìn)行fluid網(wǎng)格點(diǎn)的遷移，同時根據(jù)遷移后的值計算fluid網(wǎng)格點(diǎn)宏觀量的值；

Step6. 根據(jù)式(10)進(jìn)行邊界網(wǎng)格點(diǎn)的遷移和計算宏觀量的值；

Step7. 判斷是否滿足退出條件，如不滿足則返回Step2繼續(xù)計算，否則計算結(jié)束，輸出計算結(jié)果.

迭代計算在MRT-LBM的整個求解過程中所占比重最大，也最耗時.上述的并行算法僅需要在相鄰進(jìn)程間進(jìn)行通信，無需全局通信，具有良好的可擴(kuò)展性，能夠提高M(jìn)RT-LBM在分布式集群上的執(zhí)行效率，顯著縮短計算時間.在流場信息初始化和迭代計算部分，直接使用了并行網(wǎng)格生成部分得到的計算網(wǎng)格，每個進(jìn)程負(fù)責(zé)的網(wǎng)格數(shù)量保持一致，網(wǎng)格量負(fù)載均衡，但是由于在迭代計算中存在部分進(jìn)程中的solid網(wǎng)格點(diǎn)不參與計算，可能會對迭代計算網(wǎng)格量的負(fù)載均衡產(chǎn)生一些影響，但是由于該部分網(wǎng)格點(diǎn)在總體網(wǎng)格量上所占比例不高，所以對整個MRT-LBM并行算法的可擴(kuò)展性影響不明顯.

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

本節(jié)將針對第2節(jié)中提出的MRT-LBM高可擴(kuò)展并行算法，在國產(chǎn)神威藍(lán)光超級計算機(jī)上進(jìn)行詳細(xì)的并行測試分析.

3.1測試環(huán)境

神威藍(lán)光超級計算機(jī)是國內(nèi)首個全部采用國產(chǎn)中央處理器(CPU)和系統(tǒng)軟件構(gòu)建的千萬億次計算機(jī)系統(tǒng).其中計算節(jié)點(diǎn)共8 704個，每個節(jié)點(diǎn)采用的是具有自主知識產(chǎn)權(quán)的申威1 600(SW1600)處理器，每個處理器有16核、主頻1.0GHz，神威總計算核達(dá)到130 000核以上，系統(tǒng)峰值為1070.16TFlops，測試性能為795.9TFlops.計算節(jié)點(diǎn)的內(nèi)存為16GB，系統(tǒng)內(nèi)存容量為170TB，外存容量為2PB.神威藍(lán)光采用了胖樹結(jié)構(gòu)的互聯(lián)網(wǎng)絡(luò)，節(jié)點(diǎn)間采用了來自Mellanox的QDRInfiniband網(wǎng)絡(luò)，其傳輸速度高達(dá)40Gbps.操作系統(tǒng)為國產(chǎn)“神威睿思”并行操作系統(tǒng).

3.2測試結(jié)果

本文的測試用例選取為三維圓球繞流流動問題.測試主要分為3個部分：MRT-LBM高可擴(kuò)展并行算法的數(shù)值穩(wěn)定性和精度測試、固定問題的計算規(guī)模的并行效率測試和固定單個進(jìn)程計算規(guī)模的并行效率測試.固定問題的計算規(guī)模并行測試的計算網(wǎng)格量分別選取3億、6億、13億、26億、52億和100億這6種固定規(guī)模(分別記為固定規(guī)模1～6).固定單個進(jìn)程計算規(guī)模的并行效率測試保證每個進(jìn)程中的計算網(wǎng)格量為2 000.測試的最大核數(shù)為130 000.

首先對MRT-LBM高可擴(kuò)展并行算法的數(shù)值穩(wěn)定性和精度進(jìn)行測試分析.在三維圓球繞流問題的實(shí)驗(yàn)中，選取來流的馬赫數(shù)Ma=0.2，雷諾數(shù)Re取值范圍是10～1 000，計算的圓球阻力系數(shù)cd隨雷諾數(shù)Re變化曲線見圖6所示.從圖6可以看出，MRT-LBM并行算法的計算結(jié)果與經(jīng)驗(yàn)公式的結(jié)果基本符合，驗(yàn)證了MRT-LBM并行算法的數(shù)值穩(wěn)定性和精度.

Fig. 6　Drag coefficient at variable Re.圖6　阻力系數(shù)隨Re變化曲線

接下來進(jìn)行MRT-LBM高可擴(kuò)展并行算法性能方面測試.根據(jù)圖6中的結(jié)果，在接下來的三維圓球繞流問題的測試計算中，Re取值固定為1 000，Ma取值依然為0.2.

首先我們給出固定問題的計算規(guī)模的并行效率測試結(jié)果和分析.由于固定規(guī)模1～6計算網(wǎng)格量從幾億到百億，問題的計算規(guī)模較大，考慮到神威藍(lán)光超級計算機(jī)單個計算核心所能占用的內(nèi)存資源不到1 GB，且隨著計算問題規(guī)模的不斷擴(kuò)大，所需要使用的內(nèi)存資源也越大，為了能夠更好地進(jìn)行加速比和效率的分析，所以本文中對固定規(guī)模1～6中所采用的最小計算核數(shù)為325核，且隨著網(wǎng)格規(guī)模的增大，最小計算核數(shù)也不斷增大，而本文所采用的最大測試核數(shù)統(tǒng)一為130 000核，測試所采用的MRT-LBM的迭代步驟統(tǒng)一設(shè)為100步.

圖7和圖8分別給出了固定規(guī)模1～6的加速比和并行效率.從圖7可以看出，本文所提的并行算法加速比將隨著計算規(guī)模的增大而不斷增大，且在13萬核時具有良好的加速比.從圖8可以看出，隨著計算問題的規(guī)模增大，其計算效率越接近理論值.

Fig. 7　The speedup of different fixed scales.圖7　固定規(guī)模1～6的加速比

Fig. 8　The parallel efficiency of different fixed scales.圖8　固定規(guī)模1～6的并行效率

下面進(jìn)一步分析MRT-LBM并行算法的并行性能.考慮MRT-LBM并行算法中核心部分是迭代計算，在迭代計算中分為碰撞過程(見式(6))和遷移過程(見式(7)).根據(jù)MRT-LBM的計算流程如圖3所示，迭代計算在整個計算中所占比例最大，且隨著迭代次數(shù)的不斷增加，整個計算時間基本就是MRT-LBM的迭代計算時間.這里我們針對固定規(guī)模5的情況，分析在使用8 125～130 000計算核時MRT-LBM迭代計算(含流場信息初始化)中MPI通信時間和計算時間的關(guān)系.這里我們將MRT-LBM的迭代計算時間記為TIter，MPI通信時間記為TComm，MPI進(jìn)程的計算時間記為TComp，那么TIter=TComp+TComm.

圖9給出了在固定規(guī)模5的情況下MPI通信和計算2部分分別在MRT-LBM迭代計算時間中所占的比例(計算核數(shù)是按照2的倍數(shù)增加)，其中通信占比為TCommTIter，計算占比為TCompTIter.圖10和圖11分別給出了對固定規(guī)模5計算中的MPI通信時間TComm和計算時間TComp.從圖10和圖11可以看出，隨著計算核數(shù)的增加，迭代計算中的MPI進(jìn)程的計算時間TComp幾乎線性減少;而MPI通信時間TComm并沒有顯著增加，而是在同一個量級上緩慢變化.MPI通信時間TComm的穩(wěn)定主要是由于MRT-LBM中僅需要局部通信，且在同一計算規(guī)模下每個MPI進(jìn)程間的數(shù)據(jù)傳遞與接收量是一致的.因此圖10和圖11可以反映出圖9中的不同計算核下的時間占比，即隨著計算核數(shù)的成倍增加，計算的比重緩慢地不斷減少，而通信的比重緩慢增加.因此MRT-LBM并行算法的整體加速比隨著迭代計算中的MPI進(jìn)程的計算時間TComp占比的提高而不斷增大，但由于存在穩(wěn)定在一定量級上的MPI通信時間TComm，因此并行算法的總體加速比只能無限接近于線性加速比.

Fig. 9　The percentages of communication time and computation time among iterative computation in the fifth fixed-scale.圖9　固定規(guī)模5迭代計算中通信與計算部分的占比

Fig. 10　The MPI communication time in the fifth fixed-scale.圖10　固定規(guī)模5中不同計算核數(shù)下的MPI通信時間

Fig. 11　The iteration computation time in the fifth fixed-scale.圖11　固定規(guī)模5中不同計算核數(shù)下的迭代計算時間

結(jié)合圖7～11可以看出，MRT-LBM并行算法具有良好的可擴(kuò)展性，且在使用具有超大規(guī)模的高性能計算機(jī)時，只有當(dāng)問題的計算規(guī)模足夠大時才能更好地發(fā)揮機(jī)器的高性能.

Fig. 12　Three kinds of computation time of the sixth fixed-scale changed with cores.圖12　固定規(guī)模6中不同部分隨核數(shù)變化的計算時間

為了進(jìn)一步分析MRT-LBM并行算法的性能，接下來我們給出在固定規(guī)模6的計算情況，如圖12所示.由于固定規(guī)模6的計算網(wǎng)格達(dá)到了百億，所以考慮到內(nèi)存的限制因數(shù)，測試的計算核數(shù)從32 500取到130 000，并給出了MRT-LBM并行算法中各部分(包括總計算時間、含流場信息初始化的迭代計算時間和網(wǎng)格生成時間)隨著計算核數(shù)增加而計算時間變化的曲線.

從圖12可以看出，總計算時間和迭代計算時間隨著核數(shù)的增加而近似線性減少.但是網(wǎng)格生成計算時間卻并沒有顯著減少，而是基本穩(wěn)定在同一時間上，其主要原因是由于在進(jìn)行網(wǎng)格類型判斷和邊界距離計算時(2.1節(jié)中的Step2和Step3)，因部分進(jìn)程中不含邊界網(wǎng)格點(diǎn)，導(dǎo)致網(wǎng)格生成時網(wǎng)格類型判斷計算的負(fù)載不均衡，所以當(dāng)計算核數(shù)達(dá)到臨界值時并行網(wǎng)格生成時間幾乎保持不變.但是這樣的網(wǎng)格并行生成算法能夠與MRT-LBM迭代計算相互融合，而無需重新再次劃分計算區(qū)域，節(jié)省了再次劃分網(wǎng)格塊的時間，且MRT-LBM并行算法中網(wǎng)格生成僅是前處理的部分，所花費(fèi)的時間與迭代計算的時間相差巨大，如圖12所示，因此本文采用的并行網(wǎng)格生成算法對程序的整體并行效率并沒有產(chǎn)生明顯影響.

下面我們進(jìn)一步分析MRT-LBM并行算法在固定單個進(jìn)程計算規(guī)模的并行效率.在固定每個MPI進(jìn)程的計算網(wǎng)格量為2500的前提下，計算核數(shù)從10 000核線性增加到130 000核，給出隨著核數(shù)增加和網(wǎng)格量增大的情況下MRT-LBM并行算法的總計算時間、迭代計算時間和網(wǎng)格生成計算時間.如圖13、圖14所示，給出了在MRT-LBM的迭代計算中MPI通信和MPI進(jìn)程的計算時間占比.

Fig. 13　Three kinds of computation time of computation scale increased with cores.圖13　計算規(guī)模隨核數(shù)增大的各部分計算時間

Fig. 14　The percentages of communication time and computation time among iterative computation.圖14　迭代計算中通信與計算部分的占比

從圖13可以看出，隨著核數(shù)的增加和計算規(guī)模的增大，網(wǎng)格生成的時間、迭代計算時間和總時間幾乎穩(wěn)定在一個量級上.而從圖14可以看出，由于每個MPI進(jìn)程間的數(shù)據(jù)傳遞與接收量是一致的，且是局部網(wǎng)格塊間進(jìn)行通信，所以MPI通信的時間所在MRT-LBM迭代計算中的占比是幾乎保持不變的，同時由于計算網(wǎng)格量的一致性，所以計算時間的占比也是幾乎不變的.因此MRT-LBM并行算法在固定單個進(jìn)程計算規(guī)模時仍然具有良好的可擴(kuò)展性.

綜合以上測試結(jié)果及分析，MRT-LBM高可擴(kuò)展并行算法在13萬核的量級上依然能夠保持良好的數(shù)值穩(wěn)定性、精度和加速比，且通信時間無明顯增加，具有良好的可擴(kuò)展性.

4結(jié)束語

在三維復(fù)雜流動的大規(guī)模數(shù)值模擬中，本文針對MRT-LBM具有邊界處理簡單、并行性好等特點(diǎn)，結(jié)合大渦模擬湍流模型、曲面邊界插值格式、笛卡兒網(wǎng)格生成，詳細(xì)分析了MRT-LBM算法的并行性，并給出了基于D3Q19離散速度模型的適合于大規(guī)模分布式集群的MRT-LBM高可擴(kuò)展并行算法.該并行算法也同樣適用于基于D3Q15和D3Q27離散速度模型的MRT-LBM.通過在國產(chǎn)神威藍(lán)光超級計算機(jī)上的并行測試，驗(yàn)證了本文所提的MRT-LBM并行算法具有良好的數(shù)值穩(wěn)定性和精度、并具有良好的可擴(kuò)展性，在13萬計算核上仍能夠獲得80%以上的并行效率.

MRT-LBM并行算法有望在千萬億次計算機(jī)上精確模擬三維實(shí)際流動問題中的湍流現(xiàn)象，為湍流的理論研究和發(fā)展提供相關(guān)幫助.目前，我們還將從網(wǎng)格技術(shù)角度出發(fā)，研究多塊多重笛卡兒網(wǎng)格的自適應(yīng)MRT-LBM并行算法.

參考文獻(xiàn)

[1]Lallemand P, Luo L. Theory of the lattice Boltzmann method: Dispersion, disspation, isotropy, Galilean invariance and stability[J]. Physical Review E, 2000, 61(6): 6546-6562

[2]Succi S. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[M]. Oxford, UK: Oxford University Press, 2001

[3]Aidum C K, Clausen J R. Lattice-Boltzmann method for complex flows[J]. Annual Review of Fluid Mechanics, 2010, 42(1): 439-472

[4]Freitas R K, Henze A, Meinke M, et al. Analysis of lattice-Boltzmann methods for internal flows[J]. Computers & Fluids, 2011, 47(1): 115-21

[5]Qian Yuehong, D’Humières D, Lallemand P. Lattice BGK models for Navier-Stokes equation[J]. Europhysics Letters, 1992, 17(6): 479-484

[6]Luo Lishi, Liao Wei, Chen Xingwang, et al. Numerics of the lattice Boltzmann method: Effects of collision models on the lattice Boltzmann simulations[J]. Physical Review E, 2011, 83(5): 056710

[7]D’Humières D, Ginzburg I, Krafczyk M, et al. Multiple-relaxation-time lattice Boltzmann models in three dimensions[J]. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 2002, 360(1792): 437-451

[8]Shan Xiaowen, Chen Hudong. A general multiple-relaxation-time Boltzmann collision model[J]. International Journal of Modern Physics C, 2011, 18(4): 635-643

[9]Tosi F, Ubertini S, Succi S, et al. Optimization strategies for the entropic lattice Boltzmann method[J]. Journal of Scientific Computing, 2006, 30(3): 368-387

[10]Yu Huidan, Luo Lishi, Girimaji S S. LES of turbulent square jet flow using an MRT lattice Boltzmann model[J]. Computers & Fluids, 2006, 35(1): 957-965

[11]Krafczyk M, Tolke J, Luo Lishi. Large-eddy simulations with a multiple-relaxation-time LBE model[J]. International Journal Modern Physics B, 2003, 17(1): 33-39

[12]Menon S, Soo J H. Simulation of vortex dynamics in three-dimensional synthetic and free jets using the large-eddy lattice Boltzmann method[J]. Journal of Turbulence, 2004, 5(5): 032

[13]Jahanshaloo L, Pouryazdanpanah E, Che Sidik N A. A Review on the Application of the Lattice Boltzmann Method for Turbulent Flow Simulation[J]. Numerical Heat Transfer Applications, 2013, 64(11): 938-953

[14]Li Leisheng, Wang Chaowei, Ma Zhitao, et al. PetaPar: A scalable and fault tolerant petascal free mesh simulation system[J]. Journal of Computer Research and Development, 2015, 52(4): 823-832 (in Chinese)(黎雷生, 王朝尉, 馬志濤, 等. 千萬億次可擴(kuò)展可容錯自由網(wǎng)格數(shù)值模擬系統(tǒng)[J]. 計算機(jī)研究與發(fā)展, 2015, 52(4): 823-832)

[15]Hou S, Sterling J, Chen S, et al. A lattice Boltzmann subgrid model for high Reynolds number flows[J]. Pattern Formation and Lattice Gas Automata, 1996, 6(1): 151-166

[16]Yu Dazhi, Mei Renwei, Shyy Wei. A unified boundary treatment in lattice Boltzmann method[C]Proc of the 41st Aerospace Sciences Meeting and Exhibit. Reston, VA: AIAA, 2003: 0953

[17]Uzgoren E, Sim J, Shyy W. Marker-based, 3-D adaptive Cartesian grid method for multiphase flow around irregular geometries[J]. Communations in Computational Physics, 2009, 5(1): 1-41

Liu Zhixiang, born in 1986. PhD. His research interests include parallel algorithm and high performance computing.

Fang Yong, born in 1964. Professor and PhD supervisor. His research interests include blind signal processing, communication signal processing, and adaptive information system (yfang@staff.shu.edu.cn).

Song Anping, born in 1966. PhD and associate professor. His research interests include parallel computing, medical image processing and big data processing (apsong@shu.edu.cn).

Xu Lei, born in 1987. PhD candidate. Student member of China Computer Federation. His research interests include parallel algorithm and high performance computing (leixushu@shu.edu.cn).

Wang Xiaowei, born in 1970. Master. His research interests include high performance computing and computer architecture (wxw@shu.edu.cn).

Zhou Liping, born in 1984. Master. Member of China Computer Federation. His research interests include parallel algorithm and high performance computing (lpzhou@shu.edu.cn).

Zhang Wu, born in 1957. Professor and PhD supervisor. Senior member of China Computer Federation. His research interests include parallel algorithm and high performance computing.

Large-Scale Scalable Parallel Computing Based on LBM with Multiple-Relaxation-Time Model

Liu Zhixiang1,2, Fang Yong1, Song Anping3, Xu Lei3, Wang Xiaowei2,3, Zhou Liping2,3, and Zhang Wu2,3

1(SchoolofCommunicationandInformationEngineering,ShanghaiUniversity,Shanghai200444)2(HighPerformanceComputingCenter,ShanghaiUniversity,Shanghai200444)3(SchoolofComputerEngineeringandScience,ShanghaiUniversity,Shanghai200444)

AbstractIn the large-scale numerical simulation of three-dimensional complex flows, the multiple-relaxation-time model (MRT) of lattice Boltzmann method (LBM) has better property of numerical stability than single-relaxation-time model. Based on the turbulence model of large eddy simulation (LES) and the interpolation scheme of surface boundary, three iteration calculations of grid generation, initialization of flow information and parallelism property are analyzed respectively under the discrete velocity model D3Q19. Distributed architecture and the communication between different compute nodes using message passing interface (MPI) are often used by current high performance computing clusters. By considering both the features of distributed clusters and the load balance of calculation and using MPI programming model, the grid generation, initialization of flow information and the parallel algorithm of iteration calculation suitable for large-scale distributed cluster are studied, respectively. The proposed parallel algorithm also can be suitable for D3Q15 discrete velocity model and D3Q27 discrete velocity model. Two different cases, solving problem with fixed total calculation and solving problem with fixed calculate amount in every computing cores, are considered in the process of numerical simulation. The performances of parallelism are analyzed for these two cases, respectively. Experimental results on Sunway Blue Light supercomputer show that the proposed parallel algorithm still has good speedup and scalability on the order of hundreds of thousands of computing cores.

Key wordslarge-scale parallel computing; scalability; load balance; lattice Boltzmann method (LBM); multiple-relaxation-time model; large eddy simulation

收稿日期：2014-12-30；修回日期：2015-05-05

基金項目：國家自然科學(xué)基金重大研究計劃培育項目(91330116)

通信作者：張武(wzhang@shu.edu.cn)

中圖法分類號TP301.6

This work was supported by the Major Program of the National Natural Science Foundation of China (91330116).