基于GPU 的LBM 遷移模塊算法優化

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

（1.齊魯工業大學（山東省科學院）山東省計算中心（國家超級計算濟南中心），山東 濟南 251013；2.濟南超級計算技術研究院高性能計算實驗室，山東 濟南 251013；3.哈爾濱工業大學能源科學與工程學院，黑龍江 哈爾濱 150001）

0 引言

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

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

Palabos 是基于LBM 算法設計的計算流體力學軟件［14］。經過十幾年的發展，Palabos 的功能逐漸完善，能夠完成更多的流體模型計算，如血液流動模擬［15］、液滴碰撞［16］、多孔介質［17］等。這些研究表明了Palabos 已經成為一個熱門的研究LBM 算法的工具。KOTSALOS 等［18］在研究血液流動時使用了GPU 并行，但是在所提算法中GPU 加速的是有限元法（npFEM），而液體的流動計算是在CPU 上完成計算，因此，他們的研究中也沒有將LBM 算法完全GPU 并行。對于Palabos 來說，使用GPU 加速存在的難點是部分計算無法直接按照網格映射到計算單元，直接并行會存在數據沖突。

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

本文分析LBM 原算法的實現邏輯以及并行化的難點，介紹并行算法的設計思路，測試優化后算法的計算效率，并和原算法進行對比。

1 算法的原理與實現

1.1 LBM

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

圖1 LBM 實現流程Fig.1 Implementation process of LBM

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

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

LBM 簡化后的核心公式，計算平衡態的分布函數如下：

其中：ωk為各個速度分量的權重；ck為速度分量的方向；cs為無量綱聲速；u為速度。之后由碰撞算子來更新分布函數，如下所示：

其中：τ表示平衡態的松弛時間，與宏觀下流體的擴散系數有關。

1.2 熱點分析

針對算法進行熱點分析是并行算法優化前的必要步驟，本文選取Palabos 中的經典案例——三維方腔流動模型。表1 所示為三維方腔流動計算的計算函數熱點，該測試的網格數為1 283，其中碰撞與遷移部分計算時間占比超過70%。

表1 三維方腔流動各個計算熱點占比 Table 1 Each computing hotspot proportion of three-dimensional square cavity flow case %

之前的工作大多只將碰撞計算部分放在GPU上進行加速［20］。根據模型中的網格，將格點映射到GPU 中的計算單元，每個計算單元負責一個格點中的數據計算。因為存在數據依賴，完成碰撞計算后無法在GPU 繼續完成遷移部分的計算，需要將計算結果傳回CPU 計算，異構系統上的數據傳輸產生了大量的時間消耗，所以將遷移計算部分也放在GPU上不僅能夠達到更好的優化效果，同時對于遷移計算來說還可以省去數據傳輸的步驟，提升數據在GPU 端的利用率。

1.3 遷移計算的實現

在遷移計算中，每個格點中的速度分量需要按照一定的規律與周圍格點中的速度分量進行數據交換。Palabos 串行計算下遷移計算的偽代碼如算法1所示。

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

圖3 swap 函數的運算法則Fig.3 Operation law of swap function

根據上述偽代碼以及圖3 中sawp 函數的運算規律可以看出，在遷移計算中，除了中心點0 之外，剩余的18 個點都需要參與數據交換。每次交換涉及3 個數據：兩個為格點本身的數據，一個為周圍格點上的數據。18 個格點可分為9 組，每組都有一個固定的數據交換方向，如表2 所示。交換方式為：當前格點上位置為i和i+9 的速度分量與對應方向格點中位置i上的速度分量相互交換。

表2 數據交換方向與速度分量的對應關系 Table 2 Correspondence between data exchange direction and velocity component

1.4 并行化存在的難點

由以上分析可知，程序在串行執行時，按照網格的坐標進行遍歷，一個格點與周圍格點的數據交換會對下一個格點獲取的數據產生影響，大部分的格點需要在之前的網格完成計算后才能進行遷移計算。這樣就使得遷移計算不能像碰撞計算一樣，將每個格點的遷移直接映射到一個計算單元上，使每個格點同時執行遷移計算。

LATT 等［23］在設計Palabos 的GPU 加速算法時也指出：遷移計算不能直接映射到一個格點中，因為存在非本地的數據訪問。他們重新設計了算法，并在迭代中增加了用來推導數據來源的全局索引。這樣雖然解決了非本地數據的訪問問題，但是增加了運算成本。因此，設計合理的并行算法使得遷移計算能夠映射到每個計算單元上，能夠在不增加運算成本的情況下解決遷移計算的數據依賴問題。

2 并行算法設計

2.1 模型降維

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

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

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

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

通過分析一維模型上的數據交換規律，對整體遷移計算完成前后數據的區別進行對比，可以得到數據的交換規律。

2.2 數據定位

由圖5 可知，數據交換的方向有3 種：當方向與軸線平行時，每條射線上參與數據交換的格點數相同；當方向與對角線平行或垂直時，距離對角線越遠，參與數據交換的格點數越少。下文分別進行討論。

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

當數據交換方向與軸線平行，同方向上參與交換的格點數量相同。按照圖3 中swap 函數的交換規則，對同一條射線上的多個格點數據進行交換處理，交換方式如圖6 所示。從左至右為串行情況下遷移計算的方向，第1 個格點為邊界格點，不參與計算，但內部的遷移計算會對邊界格點數據產生影響。當參與計算的格點數大于3 時（正常計算時遠大于這個值），除了第1、2 個以及最后一個格點，中間部分格點的數據交換方式一致。

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

按照swap 函數的計算規律，遷移計算完成后同一條射線上所有格點上的數據來源存在4 種不同的類型，如表3 中類型A、B、C、D 所示。

表3 格點數據交換完成后獲取到的數據來源 Table 3 Data source obtained after the completion of the lattice point data exchange

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

當方向與對角線平行或垂直時，大部分的數據交換情況與2.2.1 節中相同，但是存在射線上參與交換的格點數小于3 的情況。針對格點數為1 或2 的情況進行數據交換處理，最終結果如圖7 所示。

圖7 同方向上的格點數為1 或2 時數據的交換方式Fig.7 Data exchange mode when the number of lattice points in the same direction is 1 or 2

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

2.3 區域劃分

完成數據定位后可以按照圖5 中一條射線上的格點數量來對整體網格區域進行劃分。例如：當射線上的點大于等于3 個時，射線上第1 個點（包括不在遷移計算區域的邊界格點）的交換類型為A，第2 個點為類型B，中間部分的點為類型C，最后一個點為類型D；當射線上的點為2 個時，第1 個點為類型A，第2 個點為類型B，第3 個點為類型D；當射線上的點只有1 個時，第1 個點為類型A，第2 個點為類型E。

將數據的交換類型映射到網格中，每個格點可以通過自身坐標來獲取完成遷移計算需要的數據。根據格點的坐標對網格進行區域劃分，可以得到速度分量分組下不同數據交換類型的區域。由于在相同射線的排布方式下區域劃分的方式一致，因此只列出第一個位置上的映射區域，剩余位置區域劃分通過改變坐標可得，結果如表4 所示，其中，N表示網格的長度。

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

通過模型降維、數據定位和區域劃分3 個步驟，可以解決串行計算中存在的數據依賴問題。程序在CPU 端處理完數據后，將數據傳輸到GPU 端并存儲在全局內存中，每個計算核心通過CUDA 模型中的線程和塊索引獲取格點的坐標，可直接完成碰撞部分的計算。為避免產生訪存沖突，需要等待所有計算核心完成碰撞計算后才能開始遷移計算，因此，碰撞計算和遷移計算分為2 個核函數來完成。之后計算核心單元通過2.3 節中區域劃分的方式來獲取當前格點計算所需要的數據，完成遷移計算任務。

3 并行算法測試與結果分析

3.1 實驗環境

本文并行算法的測試在山河超級計算機單節點上完成，節點的相關信息如表5 所示。

表5 測試環境 Table 5 Testing environment

3.2 并行效果對比

對于LBM 算法來說，由于需要大量的迭代計算，因此優化每次迭代的時間消耗能夠節省大量的時間成本。根據上文設計的并行算法，測試在128×128×128 規模網格上并行優化的效果，結果顯示：原程序的每次迭代平均運行時間在0.36 s 左右，經過優化后時間減少到0.21 s 左右。計算的加速比在1.7 左右，表明了本文設計的并行算法具有可行性。

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

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

每秒百萬網格更新數（MLUPS）是衡量LBM 算法性能的指標，計算方法如下：

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

圖9 不同網格規模下3 種算法的效率對比Fig.9 Efficiency comparison of three algorithms at different grid scales

3.3 加速比與弱可擴展性

并行算法需要良好的擴展性［25-26］，當計算的規模增大時，并行程序仍然需要保持良好的優化效果。針對這一問題，對并行算法進行其他網格維度下的計算測試。同時，對計算的加速比與僅優化碰撞部分的算法進行比較，如圖10 所示。

圖10 不同網格規模下加速比對比Fig.10 Comparison of acceleration ratios at different grid scales

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

4 結束語

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