多種任務調度混合的IB-LBM并行優化方法

劉智翔，劉慧超，黃冬梅，2*，周麗萍，蘇 誠

（1.上海海洋大學信息學院，上海201306；2.上海電力大學，上海200090；3.上海大學計算機工程與科學學院，上海200444；4.國家海洋局東海分局東海預報中心，上海200136）

0 引言

格子玻爾茲曼法（Lattice Boltzmann Method，LBM）是一種介于微觀和宏觀尺度上模擬流體運動的方法。LBM的流場演化過程清晰，在處理具有多種流體成分或者具有復雜邊界的流場時更具有優勢；但當流場中包含有復雜的幾何結構物體時，流場和浸入邊界的結合處往往會涉及到繁瑣的網格生成以及復雜的流場求解過程，造成模型的計算復雜性上升，降低了執行效率。

為了簡化流場中浸入邊界的計算，Peskin［1］提出了浸入邊界法（Immersed Boundary Method，IBM）。即用固定的歐拉網格表示流場，將浸入的邊界看作具有高剛度系數且可以發生形變的邊界，用拉格朗日點（Lagrangian Point）表示。當邊界發生形變，會產生一種恢復力，將邊界恢復到原來形狀。將邊界上的力分布到流場中，最后求解整個流場，這就是IBM的基本思想。基于IBM 和LBM 通常應用于歐拉網格這一個共同點，可以將兩種方法進行結合。Feng 等［2］首次成功應用浸入邊界-格子玻爾茲曼方法（Immersed Boundary-Lattice Boltzmann Method，IB-LBM）來模擬剛體粒子的運動。其后經過不斷的突破和發展，IB-LBM 已成為研究流固耦合問題的熱門方法。在研究細胞形變方面，Ma 等［3］將IB-LBM 引入超聲效應的影響，模擬紅細胞（Red Blood Cell，RBC）在超聲場中的聚集和形變，研究了紅細胞的瞬態聚集行為以及形變的詳細過程；對于流體與彈性物體的耦合問題，Afra 等［4］提出了一種基于IB-LBM與魯棒格子彈簧模型（Lattice Spring Model，LSM）相結合的流體與彈性體相互作用模型，研究了穩定流態、旋渦脫落流態和剛體運動極限等不同的流動條件，并提出了一種計算彈性體變形的通用公式；Gong 等［5］將非線性有限元方法（Finite Element Method，FEM）引入到IB-LBM 框架，研究非可壓縮流體中移動的可變形物體的非線性流固耦合問題；在模擬多孔介質流動問題上，Wang 等［6］將IB-LBM 應用于多孔殼體的液體俘獲和形變，以模擬漿環反應器中由多孔殼體包圍剛性固芯組成的膨脹聚乙烯（PolyEthylene，PE）顆粒在流體中的復雜流動行為。

在使用IB-LBM 求解流場時，為了得出比較精確的結果，往往需要規模較大、較密集的流場網格，這就會造成模擬過程時間長的問題。IB-LBM 模擬流場時，計算恢復力、體積力離散，計算平衡分布函數、碰撞和遷移，計算宏觀量，都是局部的，這就給IB-LBM 的并行化提供了良好的條件。Valero-Lara等［7］將IB-LBM 在基于圖形處理器（Graphics Processing Unit，GPU）的異構平臺上進行優化，改進內存管理和降低主機與加速器通信的成本；Lorenz 等［8］提出一種新的加速方法，即在GPU 上執行LBM 模型和流固耦合計算，并通過消息傳遞接口（Message Passing Interface，MPI）對共享內存資源上的計算進行并行化。雖然在GPU 上模擬IB-LBM 可以取得很好的效果，但是由于IBM的部分并行性較低，這一部分的開銷會降低并行計算的總體性能；并且在GPU 上模擬IB-LBM 還會受到帶寬的限制，也會影響計算的性能。另一種并行優化的方式，是采用共享存儲并行編程（Open Multi-Processing，OpenMP）。OpenMP 進行并行計算時，各個線程間共享內存，因此比較適合于共享內存結構的單臺計算機。在多核中央處理器（Central Processing Unit，CPU）結構上，具有效率高、內存開銷少的優勢，且編程語句簡潔直觀，易于編程。

針對IB-LBM 模擬大規模或密集網格時，模擬時間長、效率低下的缺陷，針對IB-LBM 局部計算的特點，通過OpenMP中的編譯制導語句將IB-LBM 并行化。但直接將IB-LBM 作并行化，是一種粗糙的并行化，為此在并行化的基礎上還提出了進一步的優化策略，即將IB-LBM 進行了整體上的結構化分解，對每一結構部分測試不同任務調度方式的優化性能。在不同的線程數下，即2 核、4 核、8 核、16 核下，分別進行大量實驗和分析后，得到了每一線程數下最優的組合策略。

1 浸入邊界-格子玻爾茲曼方法

IB-LBM 在模擬流體流過移動邊界或復雜邊界的情況下，是一種比較靈活且高效的計算，它的計算區域由兩類格點組成——拉格朗日點和歐拉點（Eulerian Point）（見圖1）。其中：歐拉點表示整個流場，拉格朗日點表示浸入物體的邊界，通過邊界的形變計算物體受力，即恢復力，再使用狄拉克函數（Dirac Delta Function）將恢復力離散到流場的歐拉點上，最后求解流場的宏觀密度和速度。

圖1 IB-LBM模型二維流場邊界結構Fig.1 2D flow field boundary structure of IB-LBM model

IB-LBM 模型中求解流場部分應用最廣泛的模型，是由Qian等［9］提出的DnQm（n是空間維數，m是離散速度數）模型，在二維的情況下，最常用的是D2Q9 模型，如圖2 所示，其中ei和fi分別為離散速度方向和分布函數方向。

圖2 D2Q9模型的局部結構Fig.2 Local structure of D2Q9 model

對于流場的碰撞和遷移步驟可以表示為：

其中：fi(x，t)為粒子分布函數；x和t分別為當前粒子所在的格點位置和當前的格子時間；ei為D2Q9模型中的粒子離散的方向；Δt 為格子時間步；τ 為松弛時間，它的取值由流體的運動黏性系數μ（見式（2））決定。

其中：v為當前格點上的宏觀速度，ωi是離散速度的權重系數，ei是離散方向。在D2Q9模型中，各參數的取值為：

在流場的碰撞遷移過程中，邊界的計算要進行單獨的處理。在平直邊界條件的情況下，針對邊界條件，本文采用了非平衡態外推格式［10］，其表達式為：

其中：xb為上下邊界的邊界格點，ρ(xf，t)代表邊界格點xb的相鄰格點xf的密度，vb為邊界格點xb的速度。

在求解IB-LBM 模型的過程中，核心問題是拉格朗日點上力的計算和插值。對于恢復力的計算，本文應用的為懲罰力計算方法［2］。懲罰力法的基本思想是應用胡克定律（Hook’s law），將浸入物體的邊界看作是由拉格朗日點組成的彈性邊界，且該邊界具有合適的剛度系數ε。物體受力產生形變，拉格朗日點的位置移動到參考點的位置（圖3），通過胡克定律來計算因形變產生的恢復力：

圖3 物體形變的邊界點和參考點Fig.3 Boundary point and reference point of object deformation

在計算過程中，拉格朗日點的位置是不變的，在每個格子時間步后，通過狄拉克函數來計算拉格朗日點移動后的參考點的位置。流場中體積力F(x，t)，通過邊界點上的恢復力Fb插值到流場的格點上來得到，其表達式為：

其中：X 為流場中的歐拉點；xi為邊界上第i 個拉格朗日點；Fb為相應的恢復力；n 是邊界上拉格朗日點的總數是一個狄拉克函數，可以通過該函數將邊界點上的恢復力插值到流場的歐拉點上。狄拉克函數有多種表達式，本文使用的為：

將物體邊界點上的恢復力插值到流場的歐拉點上之后，需要根據所采用的DnQm 模型，將歐拉點上的體積力進行離散［11］，其表達式為：

根據上述公式的計算后，流場中每個格點的宏觀密度和速度表達式為：

求得流場的宏觀速度后，通過狄拉克函數，可以求得邊界點上拉格朗日點的速度，其表達式為：

其中：vi為邊界上第i個拉格朗日點的速度，X為流場中的歐拉點，v(X，t)為相應的歐拉點的速度，N 為流場中歐拉點的總個數。求得邊界點的速度后，即可求得受力形變后，拉格朗日點移動后的參考點的位置：

最后，總結IB-LBM的計算步驟為：

1）設置流場及浸入邊界的初始參數，包括流場的大小、流場的初始密度、流體在入口處的初速度、浸入邊界的大小和位置等。

2）通過式（8）計算浸入邊界的恢復力，并應用式（10）將力插值到流場的歐拉點上。

3）將歐拉點上的力用式（11）進行離散，然后執行式（1），進行流場的碰撞遷移；與此同時，在遷移時執行式（7），對流場的邊界進行單獨的處理。

4）執行式（12）和式（13），計算流場的宏觀密度和速度。

5）應用式（14），將歐拉點的速度插值到邊界點上；然后執行式（15），計算拉格朗日點移動后，相應的參考點的位置。

6）回到步驟2），進行下一個格子時間步的演化，直到達到規定的時間步。

2 并行優化策略

在使用IB-LBM 進行問題模擬時，由于要得到穩定、準確的模擬結果，所以流場要進行多次演化，那么程序的負載均衡問題是影響程序性能的首要問題之一。在OpenMP 的并行優化中，提供了三種任務調度方式：靜態調度（static，以下簡稱為s 調度）、動態調度（dynamic，以下簡稱為d 調度）、啟發調度（guided，以下簡稱為g調度）。通過調節各線程的任務調度方式，優化線程的計算負載，以充分利用線程資源，提高程序計算性能。

s調度方式，同樣也是當程序中for或者parallel for 制導指令沒有選擇調度方式時，程序默認的調度方式。s調度在編譯時已經確定某一線程將要執行哪些迭代。假設計算任務為T個線程執行N 次迭代循環，那么每個線程分配到的計算任務約為連續的次迭代。

d 調度方式，線程動態申請計算任務，即當線程的計算任務完成后，繼續申請下一組計算任務。雖然d 調度可以避免平均分配所帶來的負載不均衡問題，但線程的任務動態申請會造成額外的開銷。

g調度方式，在程序開始時會給每個線程分配比較多的計算任務，之后每次分配給線程的計算任務會按指數級下降。

為了便于對IB-LBM 進行并行優化，首先根據計算步驟將IB-LBM 進行結構分解（如圖4），分解后每一部分都是獨立的，在結構塊內部不存在數據依賴，只在結構塊與結構塊之間存在數據依賴。

圖4 IB-LBM結構分解Fig.4 IB-LBM structural decomposition

算法并行最難的問題之一是要保證各個處理器的負載均衡，并減少計算資源浪費的情況。將IB-LBM 結構分解后，可以看出，恢復力的計算和參考點位置的計算這兩個結構塊的計算量大小與邊界點的數量成正比關系。在計算規模較小，即邊界點數量較少的情況下，不宜進行算法的并行，因為并行計算后的總時間，包括線程的開辟與銷毀、任務調度的開銷以及任務進行計算的時間，將大于串行計算的時間，起不到優化效果；而在計算規模較大，即邊界點較多的情況下，由于s 調度和g調度的任務調度的調度開銷要小于d調度，且兩個結構塊的計算任務簡單，各線程的處理速度大致相同，因此恢復力的計算和參考點位置的計算兩個結構塊適合使用s 調度或g調度。對于體積力的計算、碰撞和遷移、邊界處理、宏觀量的計算四個結構塊，其計算量大，且計算任務復雜，導致各個線程的計算任務結束的時間不同，有快有慢。為了避免出現資源浪費，即出現線程空閑的情況，需要根據線程的計算情況，及時地為線程分配計算任務。在三種任務調度方式中，d 調度和g 調度這兩種調度方式更加靈活，可以根據各個線程的完成情況，及時為空閑線程分配計算任務，因此，體積力的計算、碰撞和遷移、邊界處理、宏觀量的計算四個結構塊更適合于應用d調度或g調度。

根據上述內容，逐一測試三種調度方式，在IB-LBM 的每一結構步驟中，并行性能最佳的調度方式被挑選出。經過大量的實驗和對比分析后，選擇出整個IB-LBM 并行性能最優的混合調度方式。之后，根據上面的實驗方式，測試不同線程數下的并行性能最優的混合調度方式。實驗結果如表1。由表1 可以看出，使用較多的是d 調度和g 調度。雖然d 調度動態申請計算任務的方式會造成額外的時間開銷，且g 調度缺乏一定的靈活性，但在三種調度方式中，這兩種調度方式的優化性能更好一點。因此在使用不同數量的線程模擬時，這兩種調度方式的應用比較多；而s 調度相對于其他兩種調度方式，靈活性最差，導致其計算效率較低。

表2 詳細記錄了最優策略下的加速比，同時與理想狀態下的加速比和應用單一任務的加速比進行了比較。

從表2 可以看出，在線程數量較少的情況下，實際的加速比接近理想狀態的加速比。隨著線程數量的增多，線程的開辟和銷毀帶來的額外時間開銷，以及任務調度的時間開銷，對模型執行的整體效率的影響也逐漸加大，但整體的優化效果仍是比較好的，并行計算的性能有了大幅度的提高。同時，在線程數量較少的情況下，多種任務調度方式混用的優化效果與單一任務調度方式優化效果差距很小，可以得出結論：在線程較少的情況下，任務調度方式對優化效率的影響不大。在線程數量較多的情況下，應用單一任務調度與多種任務調度混合的優化效果差距比較明顯；在整體上，d 調度方式與多種調度方式混用的優化效果差距很小。

表1 多種任務調度方式混用的最佳組合Tab.1 The best combination of multiple task scheduling modes

表2 理想狀態、多種任務調度與單一調度的加速比Tab.2 Speedup of ideal state，multiple task scheduling and single scheduling

表3 中為多種任務調度最優組合下，不同的流場規模在不同線程數量時，具體的加速比數值。從表3 可以得出結論：在流場規模較小的情況下，線程的開辟和銷毀，以及任務調度所造成的額外時間開銷，在總的運行時間中影響較大，降低了并行效率；而在流場規模較大的情況下，額外的時間開銷在總時間中的占比較小，因而對并行效率的影響較小。

表3 理想狀態和不同流場規模的加速比Tab.3 Speedup of ideal state and different flow field scales

3 實驗模擬與分析

本章將模擬一個圓柱繞流問題。在流場的幾何結構上，設置圓柱的直徑大小為D，流場的大小為40D×15D，幾何結構如圖5。

圖5 流場幾何結構Fig.5 Flow field geometrical structure

實驗環境為Red Hat Enterprise Linux 6.3 操作系統，應用GCC 4.4.6編譯器進行編譯。硬件設施是Intel E5-4650 CPU，主頻為2.7 GHz，內存為512 GB RDIMM DDR3 1 600 MHz。

求解流場過程中，進行碰撞時為全部格點進行碰撞，流場格點進行遷移，對邊界格點進行額外的處理。在邊界處，采用非平衡態外推格式的邊界處理條件。在整個流場的入口處，在每個格點X 方向上添加一個初速度U。流場中，流體的流動情況可以用雷諾數（Re）來表示，計算雷諾數的表達式為其中：ρ 為流體的密度，設置流場的初始參數時，一般將ρ 設置為1.0，之后ρ 會隨著流場的演變而變化；U為流場中在入口處設置的流體的初速度，本文中U=0.1；D為圓柱的直徑；μ為流體的黏度系數。

從流場的模擬效果圖（圖6）中可以看出，隨著流場的演化，圓柱繞流的流場演化過程達到穩定狀態時，圓柱靠近入口處會形成高壓區域，在圓柱的兩側會交替形成低壓區域，即圓柱上下兩側交替出現脫落渦。

圖6 流場的模擬效果圖Fig.6 Simulation of flow field

圖7 為流場在不同時刻的流線圖。其中T 為無量綱時間尺度，其計算公式為：

其中：U 為流場的入口速度，t 為當前的格子時間，D 為圓柱直徑。

圖7 流場在不同時刻的流線圖Fig.7 Streamlines of flow field at different times

從圖7 中可以觀察得出，在圓柱周圍，沒有流線滲透進入到圓柱中，與實際情況相符，說明圓柱邊界與流場實現了無滑移邊界，因此，可以說明IB-LBM 對于模擬圓柱繞流問題是比較準確的。

為了進一步驗證本文所提方法的準確性，將給出圓柱繞流問題在Re=20 和40 情況下的阻力系數、回流區域長度（L，L的度量尺度為圓柱的半徑），并與基于平滑輪廓-格子玻爾茲曼 方 法（Smoothed Profile-Lattice Boltzmann Method，SPLBM）［12-14］的計算結果進行對比，結果如表4 所示。從表4 中可以看出，本文優化方法所計算出的結果與對比文獻中的計算結果一致。針對圓柱繞流問題在Re=100 和200 時，即流場不穩定流動的情況下，給出了阻力系數與斯特勞哈爾數(St)的計算結果，并與采用其他數值模擬方法［15-19］來計算圓柱繞流問題的文獻方法的計算結果進行對比，結果如表5。從表5 中可以看出，本文的阻力系數和St 與文獻中的計算結果保持一致。圖8展示了雷諾數為200時，阻力系數和升力系數隨時間變化的曲線，曲線呈周期性變化，驗證了流場模擬結果的準確性。。

表4 Re=20，40時阻力系數和回流區域長度（L）的對比Tab.4 Comparison of drag coefficient and recirculation zone length（L）when Re=20，40

表5 Re=100，200時阻力系數和斯特勞哈爾數（St）對比Tab.5 Comparison of drag coefficient and Strouhal number（St）when Re=100，200

圖8 Re=200時阻力和升力系數變化曲線Fig.8 Drag and lift coefficient change curves when Re=200

4 結語

本文通過OpenMP 編譯處理方案對IB-LBM 并行計算的性能問題進行了并行優化，并針對OpenMP 中不同類型的任務調度方式：s調度方式、d 調度方式、g 調度方式，進行了實驗測試。使用單一的任務調度時，只能對IB-LBM 進行粗粒度的并行優化，并行性能較低，因此進一步采用多種任務調度混合使用的方式對IB-LBM 進行了新的細粒度的并行優化。通過將IB-LBM 結構化分解，并對各結構步驟進行并行優化，測試三種任務調度的優化性能。通過上述實驗方法，在進行大量實驗和分析后，得出了在不同數量的線程參與模型的并行計算，即2 核、4 核、8 核、16 核的情況下，三種任務調度方式混合使用的最優組合策略。優化結果通過加速比來直觀地展示，可以得出，在使用的線程數量不同的情況下，多種任務調度混合使用的最優組合方式是不同的：在2 核和4 核的情況下，混合任務調度的優化效果差距不大，加速比接近理想狀態；在8核、16 核的情況下，由于開辟線程和任務調度的額外開銷，對優化后的計算性能產生了影響，但并行性能仍是較好的。同時實驗得出，線程較多的情況下，在放大規模后，由于額外的時間開銷在總的執行時間中產生的影響降低，加速比也有較大的提升。雖然本文對IB-LBM 進行了效果比較好的優化，但IBM 部分的優化效果仍然不高，有待進一步研究。