基于GPU 的LBM 遷移模塊算法優(yōu)化

黃斌，柳安軍，潘景山*，田敏，張煜，朱光慧

（1.齊魯工業(yè)大學(xué)（山東省科學(xué)院）山東省計(jì)算中心（國(guó)家超級(jí)計(jì)算濟(jì)南中心），山東 濟(jì)南 251013；2.濟(jì)南超級(jí)計(jì)算技術(shù)研究院高性能計(jì)算實(shí)驗(yàn)室，山東 濟(jì)南 251013；3.哈爾濱工業(yè)大學(xué)能源科學(xué)與工程學(xué)院，黑龍江 哈爾濱 150001）

0 引言

格子玻爾茲曼方法（LBM）是基于介觀定理的離散流體力學(xué)計(jì)算方法，通過(guò)描述粒子在固定格點(diǎn)中的運(yùn)動(dòng)來(lái)計(jì)算流體的運(yùn)動(dòng)狀態(tài)。相比于宏觀計(jì)算流體力學(xué)（CFD）方法需要求解二階偏微分方程（連續(xù)性方程、動(dòng)量方程），LBM 只需要求解一個(gè)一階偏微分方程（粒子數(shù)密度守恒方程）［1-2］，通過(guò)計(jì)算介觀微團(tuán)的碰撞與遷移過(guò)程來(lái)演變宏觀物理過(guò)程。除此之外，作為一種瞬態(tài)演變算法，LBM 還具有描述簡(jiǎn)單、易于編程、復(fù)雜邊界易于設(shè)置等特性［3］。因此，該算法在提出不久后，已經(jīng)廣泛應(yīng)用于計(jì)算流體力學(xué)的相關(guān)領(lǐng)域［4-5］。另一方面，并行計(jì)算的出現(xiàn)加速了數(shù)值模擬計(jì)算的發(fā)展，高性能計(jì)算（HPC）逐漸在流體研究中扮演了重要角色。經(jīng)過(guò)幾十年的發(fā)展，并行計(jì)算已經(jīng)發(fā)展出了多種成熟的編程模型，其中基于圖形處理器（GPU）進(jìn)行加 速的有OpenCL［6］、OpenAcc［7］、統(tǒng)一計(jì) 算設(shè)備架構(gòu)（CUDA）［8-9］等。CUDA 是基于GPU 的一種并行計(jì)算平臺(tái)與編程模型，適合處理大規(guī)模的密集計(jì)算。

與傳統(tǒng)的CFD 算法不同，LBM 算法本身就有著易于并行的特性，適合通過(guò)GPU 來(lái)對(duì)計(jì)算進(jìn)行加速。由于LBM 按照網(wǎng)格點(diǎn)來(lái)設(shè)計(jì)模型，因此計(jì)算任務(wù)可以按照網(wǎng)格點(diǎn)來(lái)分配給各個(gè)計(jì)算單元［10］?，F(xiàn)如今越來(lái)越多的研究人員選擇GPU 來(lái)對(duì)LBM 算法進(jìn)行并行加速，減少了科研計(jì)算時(shí)間：RAHMAN 等［11］利用GPU 計(jì)算基于LBM 算法的冪律非牛頓納米流體在矩形腔內(nèi)的磁流體動(dòng)力學(xué)熱溶質(zhì)自然對(duì)流流動(dòng)；WATANABE 等［12］基于LBM 研究多潮汐渦輪機(jī)，使用10 臺(tái)P100 GPU 加速，在9 h 內(nèi)完成了8.55×108規(guī)模網(wǎng)格點(diǎn)對(duì)10 臺(tái)潮汐渦輪機(jī)的大規(guī)模模擬；KIANI-OSHTORJANI 等［13］在GPU 上 基于LBM 研 究了流體與單個(gè)顆粒團(tuán)簇混合物中的耦合傳熱問(wèn)題。

Palabos 是基于LBM 算法設(shè)計(jì)的計(jì)算流體力學(xué)軟件［14］。經(jīng)過(guò)十幾年的發(fā)展，Palabos 的功能逐漸完善，能夠完成更多的流體模型計(jì)算，如血液流動(dòng)模擬［15］、液滴碰撞［16］、多孔介質(zhì)［17］等。這些研究表明了Palabos 已經(jīng)成為一個(gè)熱門(mén)的研究LBM 算法的工具。KOTSALOS 等［18］在研究血液流動(dòng)時(shí)使用了GPU 并行，但是在所提算法中GPU 加速的是有限元法（npFEM），而液體的流動(dòng)計(jì)算是在CPU 上完成計(jì)算，因此，他們的研究中也沒(méi)有將LBM 算法完全GPU 并行。對(duì)于Palabos 來(lái)說(shuō)，使用GPU 加速存在的難點(diǎn)是部分計(jì)算無(wú)法直接按照網(wǎng)格映射到計(jì)算單元，直接并行會(huì)存在數(shù)據(jù)沖突。

LBM 含有兩個(gè)計(jì)算熱點(diǎn)：碰撞（collide）和遷移（streaming）［19］。前者的計(jì)算是對(duì)格點(diǎn)本身的數(shù)據(jù)進(jìn)行計(jì)算，這樣并行后每個(gè)計(jì)算單元直接獲取數(shù)據(jù)就可以完成計(jì)算；后者則需要將自身的數(shù)據(jù)與周邊格點(diǎn)上的數(shù)據(jù)進(jìn)行交換，這部分存在著一定的數(shù)據(jù)依賴，不方便展開(kāi)并行。文獻(xiàn)［20］基于Palabos 對(duì)LBM算法做了并行優(yōu)化，由于遷移計(jì)算存在數(shù)據(jù)依賴，其只將碰撞計(jì)算部分做了并行，最終達(dá)到了1.5 左右的加速比。由此可見(jiàn)，設(shè)計(jì)遷移模塊的并行算法能夠進(jìn)一步提升算法的計(jì)算效率。

本文分析LBM 原算法的實(shí)現(xiàn)邏輯以及并行化的難點(diǎn)，介紹并行算法的設(shè)計(jì)思路，測(cè)試優(yōu)化后算法的計(jì)算效率，并和原算法進(jìn)行對(duì)比。

1 算法的原理與實(shí)現(xiàn)

1.1 LBM

LBM 由玻爾 茲曼方 程演化［21］，在模型 建立時(shí)一次性完成求解非線性偏微分方程組的工作［22］，這樣做使得研究人員在數(shù)值模擬計(jì)算中只需要處理簡(jiǎn)單的線性方程或方程組。通過(guò)離散分布函數(shù)，LBM 使用一個(gè)碰撞算子來(lái)模擬一次迭代內(nèi)分布函數(shù)的演變，由此計(jì)算得到密度、壓強(qiáng)、內(nèi)能等諸多流場(chǎng)信息。LBM 在不同案例中的實(shí)現(xiàn)流程相似，如圖1 所示，核心的2 個(gè)計(jì)算模塊分別是碰撞和遷移。

圖1 LBM 實(shí)現(xiàn)流程Fig.1 Implementation process of LBM

LBM 中的基本模型為DnQm模型（n為空間維度，m為離散格點(diǎn)上的速度分量數(shù)量），本文采用的模型為D3Q19，它的速度分量分布如圖2 所示（彩色效果見(jiàn)《計(jì)算機(jī)工程》官網(wǎng)HTML 版本，下同）。

圖2 D3Q19 速度分量模型Fig.2 D3Q19 velocity component model

LBM 簡(jiǎn)化后的核心公式，計(jì)算平衡態(tài)的分布函數(shù)如下：

其中：ωk為各個(gè)速度分量的權(quán)重；ck為速度分量的方向；cs為無(wú)量綱聲速；u為速度。之后由碰撞算子來(lái)更新分布函數(shù)，如下所示：

其中：τ表示平衡態(tài)的松弛時(shí)間，與宏觀下流體的擴(kuò)散系數(shù)有關(guān)。

1.2 熱點(diǎn)分析

針對(duì)算法進(jìn)行熱點(diǎn)分析是并行算法優(yōu)化前的必要步驟，本文選取Palabos 中的經(jīng)典案例——三維方腔流動(dòng)模型。表1 所示為三維方腔流動(dòng)計(jì)算的計(jì)算函數(shù)熱點(diǎn)，該測(cè)試的網(wǎng)格數(shù)為1 283，其中碰撞與遷移部分計(jì)算時(shí)間占比超過(guò)70%。

表1 三維方腔流動(dòng)各個(gè)計(jì)算熱點(diǎn)占比 Table 1 Each computing hotspot proportion of three-dimensional square cavity flow case %

之前的工作大多只將碰撞計(jì)算部分放在GPU上進(jìn)行加速［20］。根據(jù)模型中的網(wǎng)格，將格點(diǎn)映射到GPU 中的計(jì)算單元，每個(gè)計(jì)算單元負(fù)責(zé)一個(gè)格點(diǎn)中的數(shù)據(jù)計(jì)算。因?yàn)榇嬖跀?shù)據(jù)依賴，完成碰撞計(jì)算后無(wú)法在GPU 繼續(xù)完成遷移部分的計(jì)算，需要將計(jì)算結(jié)果傳回CPU 計(jì)算，異構(gòu)系統(tǒng)上的數(shù)據(jù)傳輸產(chǎn)生了大量的時(shí)間消耗，所以將遷移計(jì)算部分也放在GPU上不僅能夠達(dá)到更好的優(yōu)化效果，同時(shí)對(duì)于遷移計(jì)算來(lái)說(shuō)還可以省去數(shù)據(jù)傳輸?shù)牟襟E，提升數(shù)據(jù)在GPU 端的利用率。

1.3 遷移計(jì)算的實(shí)現(xiàn)

在遷移計(jì)算中，每個(gè)格點(diǎn)中的速度分量需要按照一定的規(guī)律與周?chē)顸c(diǎn)中的速度分量進(jìn)行數(shù)據(jù)交換。Palabos 串行計(jì)算下遷移計(jì)算的偽代碼如算法1所示。

上述偽代碼中主要執(zhí)行的是swap 函數(shù)，其運(yùn)算法則如圖3 所示。其中：數(shù)據(jù)a與b代表當(dāng)前格點(diǎn)中位置i與位置i+9 上的速度分量，數(shù)據(jù)c與d代表相鄰格點(diǎn)中位置i與位置i+9 上的速度分量，1≤i≤9。

圖3 swap 函數(shù)的運(yùn)算法則Fig.3 Operation law of swap function

根據(jù)上述偽代碼以及圖3 中sawp 函數(shù)的運(yùn)算規(guī)律可以看出，在遷移計(jì)算中，除了中心點(diǎn)0 之外，剩余的18 個(gè)點(diǎn)都需要參與數(shù)據(jù)交換。每次交換涉及3 個(gè)數(shù)據(jù)：兩個(gè)為格點(diǎn)本身的數(shù)據(jù)，一個(gè)為周?chē)顸c(diǎn)上的數(shù)據(jù)。18 個(gè)格點(diǎn)可分為9 組，每組都有一個(gè)固定的數(shù)據(jù)交換方向，如表2 所示。交換方式為：當(dāng)前格點(diǎn)上位置為i和i+9 的速度分量與對(duì)應(yīng)方向格點(diǎn)中位置i上的速度分量相互交換。

表2 數(shù)據(jù)交換方向與速度分量的對(duì)應(yīng)關(guān)系 Table 2 Correspondence between data exchange direction and velocity component

1.4 并行化存在的難點(diǎn)

由以上分析可知，程序在串行執(zhí)行時(shí)，按照網(wǎng)格的坐標(biāo)進(jìn)行遍歷，一個(gè)格點(diǎn)與周?chē)顸c(diǎn)的數(shù)據(jù)交換會(huì)對(duì)下一個(gè)格點(diǎn)獲取的數(shù)據(jù)產(chǎn)生影響，大部分的格點(diǎn)需要在之前的網(wǎng)格完成計(jì)算后才能進(jìn)行遷移計(jì)算。這樣就使得遷移計(jì)算不能像碰撞計(jì)算一樣，將每個(gè)格點(diǎn)的遷移直接映射到一個(gè)計(jì)算單元上，使每個(gè)格點(diǎn)同時(shí)執(zhí)行遷移計(jì)算。

LATT 等［23］在設(shè)計(jì)Palabos 的GPU 加速算法時(shí)也指出：遷移計(jì)算不能直接映射到一個(gè)格點(diǎn)中，因?yàn)榇嬖诜潜镜氐臄?shù)據(jù)訪問(wèn)。他們重新設(shè)計(jì)了算法，并在迭代中增加了用來(lái)推導(dǎo)數(shù)據(jù)來(lái)源的全局索引。這樣雖然解決了非本地?cái)?shù)據(jù)的訪問(wèn)問(wèn)題，但是增加了運(yùn)算成本。因此，設(shè)計(jì)合理的并行算法使得遷移計(jì)算能夠映射到每個(gè)計(jì)算單元上，能夠在不增加運(yùn)算成本的情況下解決遷移計(jì)算的數(shù)據(jù)依賴問(wèn)題。

2 并行算法設(shè)計(jì)

2.1 模型降維

D3Q19 模型為三維模型，每個(gè)格點(diǎn)有x、y、z3 個(gè)下標(biāo)。通過(guò)分析算法1 可知，swap 函數(shù)獲取相鄰格點(diǎn)和當(dāng)前格點(diǎn)的坐標(biāo)至少有一個(gè)下標(biāo)是一致的，同時(shí)每個(gè)位置速度分量進(jìn)行交換的相鄰格點(diǎn)方向固定。因此，可將三維模型轉(zhuǎn)換為二維模型來(lái)進(jìn)行分析。如圖4 所示，根據(jù)速度分量的方向，三維的數(shù)據(jù)交換轉(zhuǎn)換為9 種二維數(shù)據(jù)交換。降低維度一方面使得數(shù)據(jù)依賴性降低，另一方面也使離散的算法更方便并行。將模型展開(kāi)成二維后，可以根據(jù)速度分量方向來(lái)分析遷移計(jì)算在各個(gè)方向的運(yùn)算規(guī)律。

圖4 三維模型轉(zhuǎn)換為二維模型的示意圖Fig.4 Schematic diagram of converting a three-dimensional model to two-dimensional models

網(wǎng)格中所有格點(diǎn)上固定位置的速度分量都和同方向的相鄰格點(diǎn)產(chǎn)生數(shù)據(jù)交換，例如坐標(biāo)為x、y、z的格點(diǎn)，與位置1 上的速度分量進(jìn)行交換的格點(diǎn)坐標(biāo)為x-1、y、z。因此，相同位置速度分量進(jìn)行的數(shù)據(jù)交換都在同一條射線中完成，如圖5 所示。圖5（a）中射線與軸線平行，對(duì)應(yīng)的速度分量位置為1、2、3；圖5（b）中射線與正對(duì)角線平行，對(duì)應(yīng)的速度分量位置為4、6、8；圖5（c）中射線與正對(duì)角線垂直，對(duì)應(yīng)的速度分量位置為5、7、9。在二維模型的情況下，每個(gè)格點(diǎn)需要交換的數(shù)據(jù)只有固定方向上的格點(diǎn)。因此，二維模型數(shù)據(jù)交換可分解為多個(gè)方向一致、平行排列的一維射線。LBM 中流體模型邊界格點(diǎn)的碰撞與遷移計(jì)算與內(nèi)部不同，因此，進(jìn)行優(yōu)化的遷移計(jì)算不包括邊界上格點(diǎn)。

圖5 二維模型下交換射線在網(wǎng)格中的排布方式Fig.5 Arrangement mode of exchange rays in grid in a two-dimensional model

通過(guò)分析一維模型上的數(shù)據(jù)交換規(guī)律，對(duì)整體遷移計(jì)算完成前后數(shù)據(jù)的區(qū)別進(jìn)行對(duì)比，可以得到數(shù)據(jù)的交換規(guī)律。

2.2 數(shù)據(jù)定位

由圖5 可知，數(shù)據(jù)交換的方向有3 種：當(dāng)方向與軸線平行時(shí)，每條射線上參與數(shù)據(jù)交換的格點(diǎn)數(shù)相同；當(dāng)方向與對(duì)角線平行或垂直時(shí)，距離對(duì)角線越遠(yuǎn)，參與數(shù)據(jù)交換的格點(diǎn)數(shù)越少。下文分別進(jìn)行討論。

2.2.1 方向與軸線平行情況討論

當(dāng)數(shù)據(jù)交換方向與軸線平行，同方向上參與交換的格點(diǎn)數(shù)量相同。按照?qǐng)D3 中swap 函數(shù)的交換規(guī)則，對(duì)同一條射線上的多個(gè)格點(diǎn)數(shù)據(jù)進(jìn)行交換處理，交換方式如圖6 所示。從左至右為串行情況下遷移計(jì)算的方向，第1 個(gè)格點(diǎn)為邊界格點(diǎn)，不參與計(jì)算，但內(nèi)部的遷移計(jì)算會(huì)對(duì)邊界格點(diǎn)數(shù)據(jù)產(chǎn)生影響。當(dāng)參與計(jì)算的格點(diǎn)數(shù)大于3 時(shí)（正常計(jì)算時(shí)遠(yuǎn)大于這個(gè)值），除了第1、2 個(gè)以及最后一個(gè)格點(diǎn)，中間部分格點(diǎn)的數(shù)據(jù)交換方式一致。

圖6 同方向上的格點(diǎn)數(shù)大于等于3 時(shí)數(shù)據(jù)的交換方式Fig.6 Data exchange mode when the number of lattice points in the same direction is greater than or equal to 3

按照swap 函數(shù)的計(jì)算規(guī)律，遷移計(jì)算完成后同一條射線上所有格點(diǎn)上的數(shù)據(jù)來(lái)源存在4 種不同的類型，如表3 中類型A、B、C、D 所示。

表3 格點(diǎn)數(shù)據(jù)交換完成后獲取到的數(shù)據(jù)來(lái)源 Table 3 Data source obtained after the completion of the lattice point data exchange

2.2.2 方向與對(duì)角線平行或垂直情況討論

當(dāng)方向與對(duì)角線平行或垂直時(shí)，大部分的數(shù)據(jù)交換情況與2.2.1 節(jié)中相同，但是存在射線上參與交換的格點(diǎn)數(shù)小于3 的情況。針對(duì)格點(diǎn)數(shù)為1 或2 的情況進(jìn)行數(shù)據(jù)交換處理，最終結(jié)果如圖7 所示。

圖7 同方向上的格點(diǎn)數(shù)為1 或2 時(shí)數(shù)據(jù)的交換方式Fig.7 Data exchange mode when the number of lattice points in the same direction is 1 or 2

按照swap 函數(shù)的交換規(guī)律，存在一種方式與2.2.1 節(jié)中不同，如 表3 中類 型E 所示。

2.3 區(qū)域劃分

完成數(shù)據(jù)定位后可以按照?qǐng)D5 中一條射線上的格點(diǎn)數(shù)量來(lái)對(duì)整體網(wǎng)格區(qū)域進(jìn)行劃分。例如：當(dāng)射線上的點(diǎn)大于等于3 個(gè)時(shí)，射線上第1 個(gè)點(diǎn)（包括不在遷移計(jì)算區(qū)域的邊界格點(diǎn)）的交換類型為A，第2 個(gè)點(diǎn)為類型B，中間部分的點(diǎn)為類型C，最后一個(gè)點(diǎn)為類型D；當(dāng)射線上的點(diǎn)為2 個(gè)時(shí)，第1 個(gè)點(diǎn)為類型A，第2 個(gè)點(diǎn)為類型B，第3 個(gè)點(diǎn)為類型D；當(dāng)射線上的點(diǎn)只有1 個(gè)時(shí)，第1 個(gè)點(diǎn)為類型A，第2 個(gè)點(diǎn)為類型E。

將數(shù)據(jù)的交換類型映射到網(wǎng)格中，每個(gè)格點(diǎn)可以通過(guò)自身坐標(biāo)來(lái)獲取完成遷移計(jì)算需要的數(shù)據(jù)。根據(jù)格點(diǎn)的坐標(biāo)對(duì)網(wǎng)格進(jìn)行區(qū)域劃分，可以得到速度分量分組下不同數(shù)據(jù)交換類型的區(qū)域。由于在相同射線的排布方式下區(qū)域劃分的方式一致，因此只列出第一個(gè)位置上的映射區(qū)域，剩余位置區(qū)域劃分通過(guò)改變坐標(biāo)可得，結(jié)果如表4 所示，其中，N表示網(wǎng)格的長(zhǎng)度。

表4 5 種交換類型在整體網(wǎng)格下的映射區(qū)域 Table 4 Mapping areas of five exchange types under the global grid

通過(guò)模型降維、數(shù)據(jù)定位和區(qū)域劃分3 個(gè)步驟，可以解決串行計(jì)算中存在的數(shù)據(jù)依賴問(wèn)題。程序在CPU 端處理完數(shù)據(jù)后，將數(shù)據(jù)傳輸?shù)紾PU 端并存儲(chǔ)在全局內(nèi)存中，每個(gè)計(jì)算核心通過(guò)CUDA 模型中的線程和塊索引獲取格點(diǎn)的坐標(biāo)，可直接完成碰撞部分的計(jì)算。為避免產(chǎn)生訪存沖突，需要等待所有計(jì)算核心完成碰撞計(jì)算后才能開(kāi)始遷移計(jì)算，因此，碰撞計(jì)算和遷移計(jì)算分為2 個(gè)核函數(shù)來(lái)完成。之后計(jì)算核心單元通過(guò)2.3 節(jié)中區(qū)域劃分的方式來(lái)獲取當(dāng)前格點(diǎn)計(jì)算所需要的數(shù)據(jù)，完成遷移計(jì)算任務(wù)。

3 并行算法測(cè)試與結(jié)果分析

3.1 實(shí)驗(yàn)環(huán)境

本文并行算法的測(cè)試在山河超級(jí)計(jì)算機(jī)單節(jié)點(diǎn)上完成，節(jié)點(diǎn)的相關(guān)信息如表5 所示。

表5 測(cè)試環(huán)境 Table 5 Testing environment

3.2 并行效果對(duì)比

對(duì)于LBM 算法來(lái)說(shuō)，由于需要大量的迭代計(jì)算，因此優(yōu)化每次迭代的時(shí)間消耗能夠節(jié)省大量的時(shí)間成本。根據(jù)上文設(shè)計(jì)的并行算法，測(cè)試在128×128×128 規(guī)模網(wǎng)格上并行優(yōu)化的效果，結(jié)果顯示：原程序的每次迭代平均運(yùn)行時(shí)間在0.36 s 左右，經(jīng)過(guò)優(yōu)化后時(shí)間減少到0.21 s 左右。計(jì)算的加速比在1.7 左右，表明了本文設(shè)計(jì)的并行算法具有可行性。

通過(guò)程序在執(zhí)行時(shí)各個(gè)部分的計(jì)算時(shí)間可以了解算法的主要時(shí)間消耗，有利于優(yōu)化的進(jìn)行。中央處理器（CPU）端可以使用gprof 工具來(lái)計(jì)算函數(shù)的執(zhí)行時(shí)間，而GPU 端可以使用nvprof［24］工具來(lái)統(tǒng)計(jì)數(shù)據(jù)傳輸以及計(jì)算時(shí)間。在串行版本中，程序完成1×105步的迭代時(shí)間消耗在10 h 左右，經(jīng)過(guò)并行優(yōu)化，CPU-GPU 混合版本的計(jì)算時(shí)間消耗在5.8 h 左右。圖8 展示了CPU 版本和混合版本在128×128×128 網(wǎng)格上單次迭代中各個(gè)計(jì)算部分的時(shí)間消耗對(duì)比。從圖中可以看出，碰撞和遷移部分在并行后計(jì)算時(shí)間占比減少，由原本的71%減少到59%。

圖8 串行算法和并行算法迭代中各個(gè)部分的執(zhí)行時(shí)間對(duì)比Fig.8 Comparison of execution time of each part in iteration by serial algorithm and parallel algorithm

每秒百萬(wàn)網(wǎng)格更新數(shù)（MLUPS）是衡量LBM 算法性能的指標(biāo)，計(jì)算方法如下：

其中：Nx、Ny、Nz分 別為模型中x、y、z軸的長(zhǎng)度；I為迭代次數(shù)；T為計(jì)算時(shí)間。MLUPS 指標(biāo)可反映流體算法的計(jì)算效率。如圖9 所示，通過(guò)對(duì)比不同規(guī)模下串行計(jì)算、只有碰撞部分并行以及本文算法的MLUPS 指標(biāo)，表明設(shè)計(jì)的算法提升了流體計(jì)算的效率。

圖9 不同網(wǎng)格規(guī)模下3 種算法的效率對(duì)比Fig.9 Efficiency comparison of three algorithms at different grid scales

3.3 加速比與弱可擴(kuò)展性

并行算法需要良好的擴(kuò)展性［25-26］，當(dāng)計(jì)算的規(guī)模增大時(shí)，并行程序仍然需要保持良好的優(yōu)化效果。針對(duì)這一問(wèn)題，對(duì)并行算法進(jìn)行其他網(wǎng)格維度下的計(jì)算測(cè)試。同時(shí)，對(duì)計(jì)算的加速比與僅優(yōu)化碰撞部分的算法進(jìn)行比較，如圖10 所示。

圖10 不同網(wǎng)格規(guī)模下加速比對(duì)比Fig.10 Comparison of acceleration ratios at different grid scales

由圖10 可以看出，隨著網(wǎng)格規(guī)模的擴(kuò)大，算法的并行效果能保持一定的弱可擴(kuò)展性。在網(wǎng)格規(guī)模為512×512×512 時(shí)，具有1.94 的加速比，同時(shí)并行算法的計(jì)算效率相比于只并行碰撞計(jì)算部分的情況提高了30%左右。

4 結(jié)束語(yǔ)

隨著高速計(jì)算機(jī)的發(fā)展，計(jì)算流體力學(xué)逐漸成為與理論流體力學(xué)、實(shí)驗(yàn)流體力學(xué)同樣重要的研究方向。充分發(fā)揮計(jì)算機(jī)的優(yōu)勢(shì)是計(jì)算流體力學(xué)發(fā)展過(guò)程中不可或缺的一環(huán)，而這正是設(shè)計(jì)并行算法的意義所在。本文詳細(xì)分析了LBM 算法中的遷移計(jì)算部分在串行計(jì)算中的實(shí)現(xiàn)邏輯，證明了該部分并行的可行性，然后通過(guò)模型降維、數(shù)據(jù)定位、區(qū)域劃分等方法，基于CUDA 設(shè)計(jì)了并行算法。該算法成功解決了遷移計(jì)算中存在的數(shù)據(jù)依賴問(wèn)題。本文算法是基于單GPU 設(shè)計(jì)的，因此，下一階段工作將針對(duì)多GPU 的情況對(duì)算法進(jìn)行優(yōu)化，進(jìn)一步提高LBM算法的計(jì)算效率。