面向多核CPU/眾核GPU 架構(gòu)的非結(jié)構(gòu)CFD共享內(nèi)存并行計(jì)算技術(shù)

張健，李瑞田，鄧亮，*，代喆，劉杰，徐傳福

1.國(guó)防科技大學(xué) 并行與分布計(jì)算全國(guó)重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410073

2.中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000

計(jì)算流體力學(xué)（Computational Fluid Dynamics，CFD）在現(xiàn)代工業(yè)和科學(xué)中發(fā)揮著重要作用，可用于模擬和預(yù)測(cè)流體運(yùn)動(dòng)的物理和化學(xué)特性，有助于改進(jìn)產(chǎn)品（例如汽車、飛機(jī)和建筑物等）的設(shè)計(jì)和制造。目前主流的工業(yè)CFD 軟件幾乎都采用非結(jié)構(gòu)網(wǎng)格，主要是因?yàn)榉墙Y(jié)構(gòu)網(wǎng)格比結(jié)構(gòu)網(wǎng)格更適合處理復(fù)雜的幾何外形，而且可以自適應(yīng)地調(diào)整網(wǎng)格大小和形狀，更好地適應(yīng)流體流動(dòng)的特征。此外，非結(jié)構(gòu)網(wǎng)格能夠更好地處理移動(dòng)邊界和多相流體等復(fù)雜問(wèn)題，具有更好的可擴(kuò)展性和靈活性。CFD 軟件的發(fā)展一直被高性能計(jì)算（High Performance Computing，HPC）技術(shù)推動(dòng)，高性能計(jì)算能力的飛速發(fā)展可以用更精細(xì)的網(wǎng)格和更高保真的方法計(jì)算更復(fù)雜的飛行構(gòu)型［1］，比如大 渦模擬（Large Eddy Simulation，LES）或直接數(shù)值模擬（Direct Numerical Simulation，DNS）［2］、數(shù)值虛擬飛行［3］、真實(shí)物理/化學(xué)過(guò)程模擬［4］等。當(dāng)前的高性能計(jì)算已邁進(jìn)E 級(jí)時(shí)代［5-6］，多核CPU（Central Processing Unit）/眾核GPU（Graphics Processing Unit）架構(gòu)成為主流。截至2022 年末，世界Top 500 超級(jí)計(jì)算機(jī)的前10 名（除排名第2 的日本Fugaku 外）都是基于該類架構(gòu)。近年來(lái)，中國(guó)在可擴(kuò)展并行算法及并行應(yīng)用方面取得世界領(lǐng)先成果［7］，但是面臨新型高性能體系架構(gòu)，必須轉(zhuǎn)變傳統(tǒng)的并行計(jì)算范式［8］。特別是，在當(dāng)前工業(yè)軟件自主可控備受關(guān)注的背景下，實(shí)現(xiàn)自主工業(yè)軟件與國(guó)產(chǎn)高性能計(jì)算的深度融合是重要發(fā)展趨勢(shì)。

傳統(tǒng)的并行化方法是利用區(qū)域分解思想對(duì)計(jì)算網(wǎng)格進(jìn)行剖分，采用MPI（Message Passage Interface）通信庫(kù)支持并行開(kāi)發(fā)。然而隨著分區(qū)數(shù)的增多，不僅通信開(kāi)銷加大、并行效率降低，而且算法收斂效率也急速下降。新型高性能體系架構(gòu)下，應(yīng)采用MPI+X 多層次并行編程模型，即計(jì)算節(jié)點(diǎn)間基于MPI 并行，而節(jié)點(diǎn)內(nèi)采用統(tǒng)一的一種共享內(nèi)存并行編程模型X，一般為OpenMP（Open Multiprocessing）、OpenACC（Open Accelerators）、CUDA（Compute Unified Device Architecture）或OpenCL（Open Computing Language）等。當(dāng)前高性能計(jì)算系統(tǒng)的單節(jié)點(diǎn)計(jì)算能力已得到大幅提高，因此設(shè)計(jì)一種高效的節(jié)點(diǎn)內(nèi)共享內(nèi)存并行算法對(duì)提升多核CPU/眾核GPU 架構(gòu)下非結(jié)構(gòu)CFD 應(yīng)用浮點(diǎn)效率顯得尤其重要。但是，這類非結(jié)構(gòu)CFD 計(jì)算的共享內(nèi)存并行化存在兩個(gè)方面的挑戰(zhàn)：一是大量的間接索引進(jìn)行數(shù)據(jù)規(guī)約操作會(huì)引起數(shù)據(jù)競(jìng)爭(zhēng)（data racing）問(wèn)題，導(dǎo)致無(wú)法直接并行化；二是非規(guī)則、非連續(xù)數(shù)據(jù)訪存導(dǎo)致數(shù)據(jù)局部性差、浮點(diǎn)計(jì)算性能極低。應(yīng)對(duì)這些挑戰(zhàn)的策略大致可分為兩類：第1 類直接對(duì)程序內(nèi)的循環(huán)體進(jìn)行并行，比如使用OpenMP 編程模型，通過(guò)添加簡(jiǎn)單的pragmas 指令以實(shí)現(xiàn)循環(huán)級(jí)并行，其中數(shù)據(jù)競(jìng)爭(zhēng)問(wèn)題采用原子操作、著色或規(guī)約等方法來(lái)解決［9］；第2 類是基于數(shù)據(jù)剖分加任務(wù)并行的策略，在MPI 進(jìn)程內(nèi)再次劃分任務(wù)子域，為每個(gè)任務(wù)創(chuàng)建線程，不同任務(wù)并行執(zhí)行，并根據(jù)數(shù)據(jù)依賴關(guān)系設(shè)計(jì)合理的任務(wù)調(diào)度，避免數(shù)據(jù)競(jìng)爭(zhēng)［10］。兩類方法在性能和平臺(tái)適應(yīng)性方面各有優(yōu)劣。循環(huán)并行最適合于規(guī)則化計(jì)算中的結(jié)構(gòu)化并行，而任務(wù)并行最適合于不規(guī)則或遞歸算法中的非結(jié)構(gòu)化并行［11］。相較于循環(huán)并行，眾核GPU 架構(gòu)下任務(wù)并行的支持還不夠廣泛或成熟。多核CPU 架構(gòu)下的并行編程模型雖然支持任務(wù)并行，但缺乏對(duì)應(yīng)用的特性支持。需針對(duì)特定問(wèn)題進(jìn)行數(shù)據(jù)剖分。因此，需要結(jié)合非結(jié)構(gòu)CFD 應(yīng)用特點(diǎn)開(kāi)展針對(duì)性研究。

NNW-FlowStar軟件［12］（FlowStar）是由國(guó)家數(shù)值風(fēng)洞工程資助建設(shè)的自主流體仿真工業(yè)軟件，核心求解器基于非結(jié)構(gòu)有限體積方法和MPI 并行技術(shù)開(kāi)發(fā)。為適應(yīng)新型硬件發(fā)展趨勢(shì)，需要進(jìn)一步發(fā)展MPI+X 多層次并行技術(shù)。王年華等基于原子操作開(kāi)展了MPI+OpenMP 混合并行策略研究［13］。張曦等開(kāi)展了面向GPU 的圖著色方法性能優(yōu)化工作［14］。Gomes等在開(kāi)源軟件SU2 上實(shí)現(xiàn)了多種循環(huán)并行策略，且支持在著色和規(guī)約方法之間自動(dòng)切換，以平衡效率和魯棒性［15］。Al Farhan和Keyes 基于剖分復(fù)制方法在NASA（National Aeronautics and Space Administration）的PETSC-FUN3D 上實(shí)現(xiàn)了多核架構(gòu)下的多線程并行［16］。Fournier 等在法 國(guó)電力 集團(tuán)的程序Code_Saturne 上基于著色算法實(shí)現(xiàn)了兩種不同的循環(huán)并行算法［17］。綜上，目前針對(duì)非結(jié)構(gòu)網(wǎng)格CFD 共享內(nèi)存并行的研究主要集中在循環(huán)并行類的方法，任務(wù)并行方法的實(shí)踐還較少，僅在有限元軟件中開(kāi)展了少量研究［10，18］。在多核CPU/眾核GPU 架構(gòu)下，具體哪種方法適用于復(fù)雜的工業(yè)級(jí)CFD 軟件缺乏統(tǒng)一認(rèn)識(shí)，同時(shí)現(xiàn)有方法在訪存優(yōu)化方面仍有提升空間，需要進(jìn)一步深入探索。

本文在FlowStar 軟件上實(shí)現(xiàn)多種循環(huán)并行和任務(wù)并行的方法，全面對(duì)比分析各個(gè)方法的優(yōu)缺點(diǎn)以及平臺(tái)適配性，以期為面向多核/眾核架構(gòu)下的非結(jié)構(gòu)CFD 并行計(jì)算提供重要參考。同時(shí)，針對(duì)非結(jié)構(gòu)網(wǎng)格CFD 計(jì)算訪存效率低下問(wèn)題，集成網(wǎng)格重排序、循環(huán)融合、多級(jí)訪存等多種工程實(shí)用的優(yōu)化策略。最后，通過(guò)M6 機(jī)翼和CHN-T1 飛機(jī)等三維非結(jié)構(gòu)網(wǎng)格算例進(jìn)行性能測(cè)試，驗(yàn)證該方法的效果。

1 問(wèn)題描述

1.1 格心型非結(jié)構(gòu)網(wǎng)格有限體積方法

FlowStar求解的主控方程為雷諾平均N-S方程（Reynolds-Averaged Navier-Stokes，RANS），采用基于格心的有限體積法進(jìn)行離散求解。

式中：t 為時(shí)間變量；Ω 為控制體的體積；?Ω 為控制體封閉面的面積；W 為守恒變量；Fc為無(wú)黏通量；Fv為黏性通量；S 為面積；Q 為源項(xiàng)。對(duì)于格心型有限體積方法，控制體就是網(wǎng)格單元體，流場(chǎng)變量與單元質(zhì)心相關(guān)聯(lián)。對(duì)式（1）進(jìn)行離散，得到對(duì)特定控制單元ΩI的形式為

式中：I 為序號(hào)，對(duì)應(yīng)網(wǎng)格單元編號(hào)；NF為該網(wǎng)格單元擁有的面數(shù)目；ΔSm為第m 個(gè)網(wǎng)格面的面積。

對(duì)于結(jié)構(gòu)網(wǎng)格，可以通過(guò)索引IJK 表示計(jì)算空間，直接對(duì)應(yīng)于流場(chǎng)變量在計(jì)算機(jī)中的存儲(chǔ)方式，這種特性可以很方便地獲取網(wǎng)格實(shí)體（點(diǎn)、面、體）間的鄰接關(guān)系。而非結(jié)構(gòu)網(wǎng)格不具備該屬性，網(wǎng)格實(shí)體通常沒(méi)有規(guī)則的排列順序，需要額外的數(shù)據(jù)結(jié)構(gòu)來(lái)存儲(chǔ)網(wǎng)格拓?fù)湫畔?。非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)結(jié)構(gòu)可以抽象地視為點(diǎn)、面、體等網(wǎng)格實(shí)體的集合，物理量定義在這些集合上，并且需要顯式定義網(wǎng)格實(shí)體之間的連接性信息。對(duì)于格心型求解器，數(shù)值操作主要是在控制單元的表面，通過(guò)獲取相鄰單元中心的值來(lái)進(jìn)行計(jì)算，因此采用基于面的數(shù)據(jù)結(jié)構(gòu)是很自然的。例如，當(dāng)求解離散化的控制方程式（2）時(shí)，必須提供控制體面心位置的通量。如圖1 所示，為計(jì)算面f0的通量，需要提供指向共享該面的兩個(gè)單元體（c0和c1）的索引指針，以訪問(wèn)這兩個(gè)控制體的流場(chǎng)變量。

圖1 格心型求解策略網(wǎng)格拓?fù)鋽?shù)據(jù)結(jié)構(gòu)Fig.1 Grid topology data structure of cell-centered scheme

在完成式（2）的右端項(xiàng)求解后，計(jì)算得到每個(gè)控制單元的殘差R，然后通過(guò)時(shí)間離散方法推進(jìn)求解，逼近最終解。

式中：A 為系統(tǒng)矩陣；ΔW 為解的更新；Δt 為時(shí)間步長(zhǎng)；n 為迭代步。

1.2 非規(guī)則訪存模式下的數(shù)據(jù)競(jìng)爭(zhēng)問(wèn)題

在非結(jié)構(gòu)網(wǎng)格上的計(jì)算一般都是在網(wǎng)格實(shí)體集合上的循環(huán)，每一個(gè)網(wǎng)格實(shí)體執(zhí)行相應(yīng)的函數(shù)內(nèi)核。在FlowStar中，右端項(xiàng)殘差的計(jì)算可以拆分成多個(gè)基于面或體循環(huán)的計(jì)算內(nèi)核。如果一個(gè)集合上的循環(huán)只在執(zhí)行內(nèi)核期間寫入該集合上定義的數(shù)據(jù)，那么循環(huán)的每次迭代都可以并行運(yùn)行。然而如果存在數(shù)據(jù)間接訪存的情況，可能有從多個(gè)集合的迭代同時(shí)更新相同內(nèi)存位置數(shù)值的情況，從而造成數(shù)據(jù)競(jìng)爭(zhēng)問(wèn)題。這種存在間接訪存的循環(huán)在基于非結(jié)構(gòu)網(wǎng)格的CFD 應(yīng)用中很常見(jiàn)，如無(wú)黏通量的計(jì)算就屬于這類情景（算法1）。

在算法1中，面循環(huán)的計(jì)算遵循收集-分散（gather-scatter）的訪存模式。在計(jì)算面的通量時(shí)，首先需要從面左右兩端的體收集變量數(shù)據(jù)，計(jì)算完成后再完成相同體索引位置上的物理量更新。由于同一個(gè)面可能被多個(gè)體共享，如果不加以控制，不同線程同時(shí)更新可能導(dǎo)致數(shù)據(jù)沖突。以圖1 為例，f0、f1、f2如果分屬不同線程，在分別完成各自的通量計(jì)算后，可能同時(shí)更新c0的物理量，導(dǎo)致錯(cuò)誤計(jì)算。因此，需要找到解決數(shù)據(jù)競(jìng)爭(zhēng)問(wèn)題的方案，實(shí)現(xiàn)非結(jié)構(gòu)CFD 共享內(nèi)存并行化。

2 非結(jié)構(gòu)CFD 共享內(nèi)存并行化方法

詳細(xì)介紹在FlowStar 上實(shí)現(xiàn)的5 種共享內(nèi)存并行方法，其中原子操作、規(guī)約和著色方法屬于循環(huán)并行類，剖分復(fù)制、分治法屬于任務(wù)并行類。

2.1 原子操作

原子操作和鎖是典型的保證線程安全的解決方案，可用于保護(hù)對(duì)潛在沖突內(nèi)存位置的訪問(wèn)，其優(yōu)點(diǎn)是幾乎不需要更改代碼，例如通過(guò)調(diào)用CUDA 內(nèi)核的atomicAdd 實(shí)現(xiàn)聚合累加操作（算法2）。然而，當(dāng)頻繁出現(xiàn)沖突時(shí)，方案的可擴(kuò)展性有限。這是因?yàn)楫?dāng)多個(gè)原子更新同時(shí)應(yīng)用到同一個(gè)內(nèi)存位置時(shí)，每個(gè)單獨(dú)的更新都是串行的，且順序是任意的，這會(huì)降低原子更新的吞吐量［19］。

2.2 規(guī)約策略

規(guī)約策略（Reduction）的基本原理是將面循環(huán)計(jì)算的gather和scatter 過(guò)程分離。如算法3 所示，面循環(huán)計(jì)算得到物理量后不更新到體上而是先存起來(lái)，然后再進(jìn)行體循環(huán)，完成物理量的規(guī)約。拆分出的兩個(gè)循環(huán)均無(wú)數(shù)據(jù)沖突，因此可以并行。規(guī)約策略除了需要耗費(fèi)額外的內(nèi)存外，另一個(gè)缺點(diǎn)是數(shù)據(jù)局部性較差，因?yàn)閷?duì)循環(huán)進(jìn)行了拆分，在讀取數(shù)據(jù)和寫入數(shù)據(jù)時(shí)都不能很好地重用數(shù)據(jù)。

2.3 面著色技術(shù)

面著色對(duì)相鄰面使用不同顏色（color）進(jìn)行標(biāo)記，屬于同一顏色的面之間不會(huì)產(chǎn)生任何數(shù)據(jù)沖突，從而保證同一顏色組內(nèi)的計(jì)算可以并行開(kāi)展，如圖2 所示。

圖2 對(duì)網(wǎng)格面進(jìn)行著色Fig.2 Face coloring of meshes

算法4 給出了基于面著色技術(shù)的OpenMP 并行實(shí)現(xiàn)偽代碼。按照顏色順序進(jìn)行計(jì)算，對(duì)每個(gè)顏色內(nèi)的循環(huán)并行化。面著色技術(shù)比較通用，可以移植到其他硬件平臺(tái)上（比如GPU）［20］，但缺點(diǎn)是在不連續(xù)的面上循環(huán)計(jì)算，數(shù)據(jù)局部性較差，甚至低于規(guī)約策略。

為改善圖著色算法的數(shù)據(jù)局部性，在Flow-Star 中進(jìn)一步實(shí)現(xiàn)了分組著色的并行算法。與面著色處理單個(gè)面單元不同的是，分組著色處理的基本單元是擁有固定面數(shù)的分組。如圖3 所示，著色后同一種顏色的面中含有若干個(gè)groupsize的面分組，組與組之間不存在數(shù)據(jù)競(jìng)爭(zhēng)。算法5給出了基于分組著色技術(shù)的OpenMP 并行實(shí)現(xiàn)偽代碼。為了完成細(xì)粒度并行任務(wù)調(diào)度，同一面分組內(nèi)的面序號(hào)（faces）需要重排以實(shí)現(xiàn)連續(xù)存儲(chǔ)，其中調(diào)度模式schedule（static，groupsize）表明每個(gè)線程會(huì)處理groupsize 大小的面分組單元。事實(shí)上，當(dāng)groupsize=1時(shí)，分組著色就退化成了普通的著色算法。分組著色的一個(gè)缺點(diǎn)就是groupsize 的大小選取經(jīng)驗(yàn)性強(qiáng)，需要根據(jù)具體網(wǎng)格大小和機(jī)器緩存進(jìn)行優(yōu)選，如果太大可能造成分組著色不成功，如果太小會(huì)導(dǎo)致數(shù)據(jù)局部性不夠，影響性能。在后續(xù)的性能測(cè)試中，根據(jù)多次測(cè)試優(yōu)選，選取groupsize=512。

圖3 對(duì)網(wǎng)格面進(jìn)行分組著色Fig.3 Face group coloring of meshes

2.4 剖分復(fù)制策略

剖分復(fù)制（Divide &Replication）策略是在為MPI 并行計(jì)算完成區(qū)域分解后，在進(jìn)程內(nèi)部根據(jù)線程數(shù)再進(jìn)行一次分區(qū)，如圖4 所示?？梢岳^續(xù)使用圖分區(qū)工具M(jìn)ETIS［21］進(jìn)行分區(qū)保證負(fù)載均衡。對(duì)面循環(huán)或體循環(huán)的計(jì)算進(jìn)行任務(wù)劃分，從而為每個(gè)OpenMP 線程分配面或體的所有權(quán)。與分組著色的思想類似，都是盡量讓單個(gè)線程處理一組內(nèi)存位置相近的單元來(lái)提升線程內(nèi)部的數(shù)據(jù)局部性，進(jìn)而提升并行性能。不同的是，剖分復(fù)制策略的核心是通過(guò)工作復(fù)制來(lái)保證線程安全［16］，而不是直接在循環(huán)上進(jìn)行并行，因此更接近任務(wù)并行。如果兩個(gè)線程共享1 個(gè)面，那這兩個(gè)進(jìn)程將復(fù)制這個(gè)面的工作。在寫回操作時(shí)，通過(guò)額外的邏輯判斷保證每個(gè)線程只負(fù)責(zé)自己分區(qū)內(nèi)的體單元數(shù)據(jù)更新，從而避免數(shù)據(jù)沖突，如算法6 所示。剖分復(fù)制策略的優(yōu)點(diǎn)是每個(gè)線程都使用其專用緩沖區(qū)，避免了工作線程之間進(jìn)行全局同步的開(kāi)銷，而且不用改變?cè)嫉淖兞看鎯?chǔ)數(shù)據(jù)結(jié)構(gòu)。但該策略如果應(yīng)用于擁有成千上萬(wàn)輕量級(jí)線程的GPU 平臺(tái)上，計(jì)算負(fù)載較小且存在大量的分支判斷，性能提升效果有限。

圖4 對(duì)計(jì)算網(wǎng)格進(jìn)行多級(jí)分區(qū)Fig.4 Hierarchical partitioning on computational grid

2.5 基于分治法的任務(wù)級(jí)并行

分治法（Divide &Conquer，D&C）的思想和剖分復(fù)制策略類似，需要對(duì)MPI 進(jìn)程內(nèi)的網(wǎng)格進(jìn)一步剖分，區(qū)別是遞歸地構(gòu)建任務(wù)樹，并基于任務(wù)樹遍歷實(shí)現(xiàn)task 并行算法［18］，如圖5 所示。在每個(gè)遞歸級(jí)別上，都會(huì)創(chuàng)建3 個(gè)子分區(qū)，包括兩個(gè)獨(dú)立的左右子分區(qū)和1 個(gè)分隔分區(qū)。其中，分隔分區(qū)由分區(qū)交界面上的網(wǎng)格實(shí)體所構(gòu)成，這樣左右子分區(qū)不具備連通性，從而無(wú)數(shù)據(jù)依賴可并行執(zhí)行。所有子分區(qū)（包括分隔分區(qū)）將繼續(xù)遞歸分下去，終止條件是葉子節(jié)點(diǎn)的網(wǎng)格單元數(shù)小于某個(gè)閾值，閾值大小一般選擇緩存區(qū)（Cache）大小，以發(fā)揮最大性能。

圖5 分治遞歸樹建立過(guò)程Fig.5 Building process of divide and conquer recursive tree

執(zhí)行計(jì)算時(shí)，采用task 并行編程模型對(duì)樹節(jié)點(diǎn)進(jìn)行遞歸遍歷執(zhí)行（算法7）：左分區(qū)和右分區(qū)并行執(zhí)行，一旦它們同步，就執(zhí)行分隔區(qū)，從而實(shí)現(xiàn)基于任務(wù)級(jí)的共享內(nèi)存并行。這種并行方式適合于支持task 的并行編程模型，如OpenMP task、Cilk Plus等。分治法的優(yōu)點(diǎn)是保持了良好的數(shù)據(jù)局部性，并將全局同步轉(zhuǎn)換成多個(gè)局部同步，減少了同步開(kāi)銷。但通用GPU 編程模型（如OpenACC、CUDA、OpenCL 等）暫不支持細(xì)粒度任務(wù)調(diào)度機(jī)制，需要顯式地管理任務(wù)樹的調(diào)度，實(shí)現(xiàn)難度較大且應(yīng)用效果有限。

3 非結(jié)構(gòu)CFD 訪存優(yōu)化技術(shù)

非連續(xù)數(shù)據(jù)訪存導(dǎo)致數(shù)據(jù)局部性差是影響非結(jié)構(gòu)CFD 共享內(nèi)存并行性能的主要障礙，引入3 種訪存優(yōu)化技術(shù)以期進(jìn)一步提升并行計(jì)算效率。

3.1 網(wǎng)格重排序

由于非結(jié)構(gòu)網(wǎng)格單元大小和形狀的不規(guī)則性，相鄰網(wǎng)格單元之間的編號(hào)跨度很大，因此在對(duì)數(shù)據(jù)進(jìn)行訪問(wèn)時(shí)，緩存系統(tǒng)可能無(wú)法有效緩存所需的數(shù)據(jù)，訪存效率低下?？梢钥紤]對(duì)網(wǎng)格單元體進(jìn)行重排序，以改善數(shù)據(jù)局部性，提高緩存命中率。常見(jiàn)的排序算法包括Cuthill-McKee［22］（CMK）、空間填充曲線（Space-filling curve）［23］、圖分區(qū)方法［24］等。分別實(shí)現(xiàn)了上述3 種排序算法，圖6 給出了使用不同算法對(duì)某個(gè)翼形網(wǎng)格重排序后得到的系統(tǒng)矩陣形式，其中矩陣的非零項(xiàng)塊子矩陣用填充矩形表示?？梢钥吹剑嘏判蚝笙噜従W(wǎng)格體之間編號(hào)更接近，矩陣的帶寬顯著降低，結(jié)構(gòu)更規(guī)則。在第4.3 節(jié)還將詳細(xì)地對(duì)這3 種排序方法的應(yīng)用效率進(jìn)行對(duì)比。

圖6 使用不同算法進(jìn)行網(wǎng)格重排序后的系統(tǒng)矩陣Fig.6 System matrix after grid reordering using different algorithms

完成體單元排序后，再根據(jù)面左右的最小體編號(hào)，重新對(duì)面編號(hào)進(jìn)行排序，以實(shí)現(xiàn)面循環(huán)計(jì)算中對(duì)面相鄰體單元近似最優(yōu)的內(nèi)存訪問(wèn)，如圖7所示。所有排序算法的復(fù)雜度為O（N）或最多為O（N lg N），所以重排序的開(kāi)銷是值得付出的［25］。

圖7 面循環(huán)中面對(duì)體的訪問(wèn)關(guān)系Fig.7 Face-to-cell access in face-loop

3.2 循環(huán)融合

為了提高代碼的可讀性和模塊化水平，F(xiàn)lowStar 程序?qū)σ恍┌罅垦h(huán)的子程序進(jìn)行了封裝，比如在無(wú)黏通量計(jì)算中，將加載變量、重構(gòu)、通量計(jì)算、更新右端項(xiàng)等封裝成了單獨(dú)的函數(shù)。這種模式很可能導(dǎo)致數(shù)據(jù)局部性下降，因?yàn)橥粋€(gè)數(shù)組在前面的循環(huán)中被依次寫入，后面的循環(huán)中很可能會(huì)再次依次讀入。當(dāng)數(shù)組長(zhǎng)度很大時(shí)，需要讀取的數(shù)據(jù)已經(jīng)不在緩存中，必須重新從內(nèi)存加載。通常可采用循環(huán)融合的優(yōu)化策略以改善數(shù)據(jù)局部性，如圖8 所示，將子程序的核函數(shù)提取到同一個(gè)循環(huán)內(nèi)。這樣，被寫入的值可以被及時(shí)讀到，降低了緩存缺失（cache miss）的概率。此外，循環(huán)體數(shù)目的減少可消除循環(huán)并行時(shí)冗余的全局同步開(kāi)銷，進(jìn)一步提升計(jì)算效率。

圖8 無(wú)黏通量計(jì)算函數(shù)進(jìn)行循環(huán)融合Fig.8 Loop fusion on inviscid flux subroutine

3.3 GPU 共享存儲(chǔ)優(yōu)化方法

非結(jié)構(gòu)CFD 的非連續(xù)訪存行為使得其在GPU 上的性能表現(xiàn)不及預(yù)期，結(jié)合SMEM（Shared Memory）進(jìn)行多級(jí)訪存是一種潛在的優(yōu)化方法。GPU 的SMEM 是一種在block 內(nèi)能訪問(wèn)的內(nèi)存，存儲(chǔ)硬件位于芯片上（on chip），訪問(wèn)速度較快。SMEM 的控制與生命周期管理與L1緩存不同，SMEM 的使用受用戶控制，L1 受系統(tǒng)控制。因此，可以顯式管理SMEM，緩存一些需要反復(fù)讀寫的數(shù)據(jù)，提升訪存效率。文獻(xiàn)［14，26］嘗試基于著色算法進(jìn)行非結(jié)構(gòu)CFD 的SMEM 優(yōu)化，但因著色本身會(huì)加劇訪存的不連續(xù)性，導(dǎo)致性能提升有限。提出了一種新的基于規(guī)約操作（算法3）的SMEM 多級(jí)訪存優(yōu)化方法。GPU 核函數(shù)預(yù)先將體單元數(shù)據(jù)從全局存儲(chǔ)（global memory）加載到SMEM 單元，面到體的規(guī)約操作在SMEM 上完成，再將更新后的SMEM 數(shù)據(jù)寫回全局存儲(chǔ)，整個(gè)數(shù)據(jù)流如圖9 所示。在該示例中，對(duì)于三維非結(jié)構(gòu)網(wǎng)格，每個(gè)體至少有4 個(gè)面單元，所以對(duì)同一內(nèi)存位置要進(jìn)行4次數(shù)據(jù)寫操作，SMEM 上的多次寫入可以增加其數(shù)據(jù)復(fù)用率，在抵消SMEM 初始加載和最終寫回全局存儲(chǔ)的開(kāi)銷后，帶來(lái)額外的性能提升。一般來(lái)說(shuō)，對(duì)于GPU 并行計(jì)算，流場(chǎng)變量或者幾何信息的存儲(chǔ)方式最好是結(jié)構(gòu)體數(shù)組（Struct of Array，SoA），而不是數(shù)組結(jié)構(gòu)體（Array of Struct，AoS）。因?yàn)樵贕PU 上并行執(zhí)行的過(guò)程，其實(shí)是將調(diào)用同一個(gè)核函數(shù)的多個(gè)線程歸并在一個(gè)塊（block）當(dāng)中，用SIMD（Single Instruction Multiple Data）的方式執(zhí)行。因此按照塊進(jìn)行數(shù)據(jù)組織的效率更高，也就是更傾向于SoA。由于規(guī)約策略保持了原有塊組織數(shù)據(jù)結(jié)構(gòu)的連續(xù)性，且數(shù)據(jù)復(fù)用率更高，所以相比著色策略進(jìn)行SMEM優(yōu)化的效果會(huì)更好。

圖9 顯式管理SMEM 的數(shù)據(jù)流Fig.9 Data flow by explicit management of SMEM

4 實(shí)驗(yàn)結(jié)果與分析

4.1 測(cè)試用例及環(huán)境

為了分析不同并行算法的性能及相關(guān)優(yōu)化技術(shù)的效果，分別采用ONERA M6 機(jī)翼（圖10）及CHN-T1 運(yùn)輸機(jī)標(biāo)模［27］（圖11）的三維非結(jié)構(gòu)網(wǎng)格開(kāi)展性能測(cè)試。其中M6 的測(cè)試網(wǎng)格按體單元數(shù)目規(guī)模分4 個(gè)層次：粗（12 萬(wàn)）、中（80 萬(wàn)）、細(xì)（360 萬(wàn)）和極細(xì)（1 000 萬(wàn)）。CHN-T1 的網(wǎng)格規(guī)模分粗（600 萬(wàn)）和細(xì)（1 700 萬(wàn)）兩個(gè)層次。

圖10 ONERA M6 機(jī)翼網(wǎng)格Fig.10 Computational grid of ONERA M6 wing

圖11 CHN-T1 民用客機(jī)標(biāo)模網(wǎng)格Fig.11 Computational grid of civil passenger aircraft model CHN-T1

測(cè)試分別在國(guó)家超級(jí)計(jì)算濟(jì)南中心的x86 節(jié)點(diǎn)和中國(guó)空氣動(dòng)力研究與發(fā)展中心的內(nèi)部GPU 工作站開(kāi)展。其中濟(jì)南中心x86 節(jié)點(diǎn)的CPU 為Intel?Xeon?Gold 6258R @2.70 GHz，每個(gè)節(jié)點(diǎn)兩路CPU，共58 個(gè)計(jì)算核心。GPU 工作站配置為Intel?Xeon?Platinum 8368 CPU @2.40 GHz 以及NVIDIA A100 PCIe 80G 加速卡。

4.2 一致性測(cè)試

理論上，不同的并行優(yōu)化方法只會(huì)改變空間離散過(guò)程中通量計(jì)算和殘差更新的順序，不會(huì)對(duì)最終計(jì)算結(jié)果造成影響。所以在進(jìn)行性能分析前，先進(jìn)行一致性測(cè)試以驗(yàn)證程序?qū)崿F(xiàn)的正確性。選取對(duì)M6 機(jī)翼中等網(wǎng)格算例進(jìn)行計(jì)算，計(jì)算狀態(tài)為來(lái)流馬赫數(shù)0.839 5、攻角3.06°、雷諾數(shù)1.814×107、溫度288.15 K。圖12 給出了利用不同方法進(jìn)行計(jì)算得到的殘差隨迭代步收斂曲線，可以看到所有方法的收斂曲線完全重合，表1 給出了最終計(jì)算得到升阻力系數(shù)，也完全一致。測(cè)試結(jié)果驗(yàn)證了所實(shí)現(xiàn)不同并行優(yōu)化方法的正確性。

表1 不同并行優(yōu)化方法計(jì)算得到的最終氣動(dòng)力系數(shù)Table 1 Final aerodynamic coefficients computed using different parallel optimization methods

圖12 不同并行優(yōu)化方法殘差隨迭代步收斂歷程比較Fig.12 Comparison of residual convergence process with iteration steps of different parallel optimization methods

4.3 訪存優(yōu)化效果

圖13 給出了使用不同排序算法模擬M6 機(jī)翼粗網(wǎng)格算例的殘差收斂過(guò)程。所有計(jì)算都取相同的計(jì)算條件，CFL（Courant-Friedrichs-Lewy）數(shù)均設(shè)為100?？v坐標(biāo)殘差取連續(xù)性方程對(duì)數(shù)刻度，橫坐標(biāo)為程序執(zhí)行墻鐘時(shí)間。可以看到，3 種排序方法中，CMK 方法的效果最好，相比原始沒(méi)有重排序的性能，殘差下降4 個(gè)量級(jí)時(shí)所需要時(shí)間要少10%左右。該結(jié)果符合預(yù)期，從圖6 可以看出，CMK 方法排序后數(shù)據(jù)局部性最好，系統(tǒng)矩陣帶寬最小。

圖13 不同重排序算法的殘差收斂歷程比較Fig.13 Comparison of residual convergence histories using different renumbering algorithms

表2 給出了循環(huán)融合的優(yōu)化效果，對(duì)M6 機(jī)翼兩套網(wǎng)格分別使用不同線程進(jìn)行計(jì)算，統(tǒng)計(jì)了前100 步迭代的平均耗時(shí)。從結(jié)果來(lái)看，循環(huán)融合后的性能提升效果明顯，運(yùn)行時(shí)間能夠縮短10% 左右。所以，在之后的測(cè)試中，默認(rèn)采用CMK 網(wǎng)格重排序和循環(huán)融合的性能優(yōu)化策略。

表2 循環(huán)融合前后計(jì)算時(shí)間（100 步）對(duì)比Table 2 Comparison of executing time（100 iterations）before and after loop merge

4.4 CPU 平臺(tái)并行計(jì)算性能分析

第2 節(jié)中描述了不同的共享內(nèi)存并行方法，采用M6 機(jī)翼和CHN-T1 標(biāo)模完成不同方法在CPU 平臺(tái)的性能對(duì)比。首先在單節(jié)點(diǎn)內(nèi)使用不同OpenMP 線程對(duì)M6 機(jī)翼網(wǎng)格進(jìn)行測(cè)試。由于單節(jié)點(diǎn)內(nèi)存有限，選取粗、中兩個(gè)層次的網(wǎng)格進(jìn)行了計(jì)算，得到并行效率數(shù)據(jù)，如圖14 所示。從結(jié)果來(lái)看，剖分復(fù)制策略的并行效率最高，可擴(kuò)展性最好，規(guī)約策略次之，原子操作的效果最差。兩種著色策略對(duì)比，基于簡(jiǎn)單著色的方案在線程數(shù)較少時(shí)要優(yōu)于分組著色策略，然而隨著線程數(shù)增多，分組著色的并行效率保持更高，驗(yàn)證了分組著色能夠通過(guò)改善數(shù)據(jù)局部性提升性能?；贒&C 樹的任務(wù)級(jí)并行效果不如剖分復(fù)制策略，可能是在遞歸創(chuàng)建任務(wù)樹的過(guò)程中，每個(gè)子分區(qū)的網(wǎng)格編號(hào)無(wú)法保持連續(xù)，導(dǎo)致訪存效率低下，因此還有性能提升空間。

圖14 不同OpenMP 并行策略的并行效率對(duì)比Fig.14 Comparison of parallel efficiency of different OpenMP parallel strategies

圖15 展示了幾個(gè)典型狀態(tài)下使用不同并行策略模擬M6 中等網(wǎng)格計(jì)算前100 步迭代的平均耗時(shí)，結(jié)果規(guī)律和并行效率一致。

圖15 不同OpenMP 并行策略的運(yùn)行時(shí)間對(duì)比Fig.15 Comparison of running time of different OpenMP parallel strategies

圖14和圖15 是純OpenMP 并行計(jì)算性能的測(cè)試，接下來(lái)利用CHN-T1 細(xì)網(wǎng)格算例開(kāi)展跨節(jié)點(diǎn)性能測(cè)試。圖16 給出了純MPI 并行和MPI+OpenMP 混合并行模式下以單節(jié)點(diǎn)計(jì)算時(shí)間為基數(shù)的加速比曲線。測(cè)試所用的線程數(shù)和CPU核心數(shù)一致，混合并行每個(gè)節(jié)點(diǎn)使用2 個(gè)進(jìn)程，OpenMP 并行采用性能最佳的剖分復(fù)制策略。如采用2 個(gè)節(jié)點(diǎn)（每個(gè)節(jié)點(diǎn)56 個(gè)核心）計(jì)算時(shí)，純MPI 模式就采用112 個(gè)MPI 進(jìn)程，混合模式采用4 個(gè)MPI 進(jìn)程以及每個(gè)進(jìn)程28 個(gè)OpenMP 線程的組合。從結(jié)果來(lái)看，在小規(guī)模下，純MPI和MPI+OpenMP 混合模式的并行效率均非常高。超過(guò)16 個(gè)節(jié)點(diǎn)后加速比開(kāi)始下降，這是因?yàn)殡S著分區(qū)數(shù)增多，每個(gè)分區(qū)的網(wǎng)格的數(shù)量急速下降，通信開(kāi)銷比例增大，降低了并行計(jì)算的效率。而MPI+OpenMP 混合模式的并行效率要高于純MPI 模式。分析其原因，同一節(jié)點(diǎn)內(nèi)的兩個(gè)核心之間的通信比節(jié)點(diǎn)之間的通信帶寬更高、延遲更低，即使通過(guò)高速網(wǎng)絡(luò)互聯(lián)也是如此，采用混合模式的分區(qū)數(shù)要遠(yuǎn)小于純MPI 模式，而且采用共享內(nèi)存的并行減少了了大量跨節(jié)點(diǎn)通信的開(kāi)銷，所以MPI+OpenMP 混合模式擴(kuò)展性更好，更加適用于大規(guī)模算例的計(jì)算。

圖16 并行計(jì)算CHN-T1 的加速比Fig.16 Speed-up ratios of parallel computation of CHN-T1

4.5 GPU 平臺(tái)并行計(jì)算性能分析

如2.5 節(jié)中所述，由于通用GPU 并行編程模型對(duì)任務(wù)并行模式的支持能力較弱，通常需采用或者設(shè)計(jì)針對(duì)任務(wù)并行的GPU 運(yùn)行時(shí)調(diào)度機(jī)制，開(kāi)發(fā)工作難度大，且平臺(tái)/應(yīng)用適用性較弱。因此，在FlowStar 上實(shí)現(xiàn)了原子操作、面著色、規(guī)約3 種基于循環(huán)并行的方法。首先利用M6 所有規(guī)模的網(wǎng)格算例進(jìn)行了測(cè)試。對(duì)于存在寫沖突的子程序是GPU 并行的最大性能障礙，對(duì)這些子程序進(jìn)行并行的加速比最能體現(xiàn)不同方法的性能優(yōu)劣。圖17 給出了FlowStar 中UpdateRHS函數(shù)GPU 并行化后相比CPU 串行計(jì)算取得的加速比?？梢钥闯?，規(guī)約和面著色策略加速效果相差無(wú)幾，且均優(yōu)于原子操作。此外，網(wǎng)格規(guī)模越大，加速效果越明顯，這是因?yàn)殡S著問(wèn)題規(guī)模的增加，GPU 計(jì)算負(fù)載提高，且CPU和GPU之間的數(shù)據(jù)移動(dòng)的開(kāi)銷比重降低，能夠充分發(fā)揮GPU眾核并行優(yōu)勢(shì)。

圖17 不同規(guī)模M6 算例UpdateRHS 函數(shù)加速比Fig.17 Speed-up ratios of UpdateRHS subroutine for different grid levels of M6 wing

3.3 節(jié)講述了通過(guò)顯式管理SMEM，緩存一些需要反復(fù)讀寫的數(shù)據(jù)可以進(jìn)一步提升訪存效率。圖18 給出了SMEM 優(yōu)化前后UpdateRHS函數(shù)計(jì)算加速比對(duì)比。可以看到基于面著色方法進(jìn)行SMEM 優(yōu)化的效果不明顯（圖中紅色曲線），結(jié)論和文獻(xiàn)［14］一致。進(jìn)一步參考文獻(xiàn)［26］，利用分組多級(jí)著色策略（黑色曲線）的加速效果同樣不好，文獻(xiàn)中取得較好性能加速的原因主要在于將AoS 的變量存儲(chǔ)方式改成SoA，而FlowStar 的數(shù)據(jù)結(jié)構(gòu)已經(jīng)是SoA。相比之下，基于規(guī)約策略（藍(lán)色曲線）的SMEM 優(yōu)化效果非常好，優(yōu)化過(guò)后可以進(jìn)一步將性能提升3倍，驗(yàn)證了所提方法的有效性。

圖18 共享存儲(chǔ)優(yōu)化前后UpdateRHS 計(jì)算加速比對(duì)比Fig.18 Comparison of speed-up ratios of UpdateRHS before and after shared memory optimization

最后，利用優(yōu)化前后的GPU 并行程序?qū)HN-T1 粗細(xì)兩套網(wǎng)格進(jìn)行了計(jì)算，統(tǒng)計(jì)了相比CPU 串行計(jì)算，各個(gè)子模塊以及程序整體運(yùn)行時(shí)間的性能加速比，如圖19 所示?？梢钥吹郊铀傩Ч浅Ｃ黠@，程序整體運(yùn)行的加速比能達(dá)到127。特別是通量計(jì)算部分，子程序的加速比可以達(dá)到400 以上，因?yàn)橥坑?jì)算部分大部分都是面循環(huán)，通過(guò)規(guī)約加SMEM 優(yōu)化的策略在解決數(shù)據(jù)競(jìng)爭(zhēng)問(wèn)題的基礎(chǔ)上緩解了間接訪存情況，并且提升了數(shù)據(jù)局部性。LU-SGS 算法部分加速比略低，是因?yàn)殡[式求解算法固有的串行屬性，這里采用了簡(jiǎn)單的體著色并行策略，性能依然有提升空間，后續(xù)還將繼續(xù)深入研究。

圖19 CHN-T1 算例GPU 并行計(jì)算整體加速比Fig.19 Overall speed-up ratios of parallel computation of CHN-T1 on GPU

5 結(jié)論

針對(duì)非結(jié)構(gòu)CFD 應(yīng)用，在多核CPU 平臺(tái)和眾核GPU 平臺(tái)上分別實(shí)現(xiàn)了基于循環(huán)級(jí)與任務(wù)級(jí)的多種共享內(nèi)存并行方法，設(shè)計(jì)了網(wǎng)格重排序、循環(huán)融合、多級(jí)訪存等性能優(yōu)化技術(shù)。具體地，針對(duì)多核CPU 平臺(tái)，全面覆蓋了已有的5 種共享內(nèi)存并行方法；針對(duì)眾核GPU 平臺(tái)，提出了基于SMEM 優(yōu)化的規(guī)約策略，提升了GPU 并行計(jì)算性能。所有方法均在工業(yè)級(jí)CFD 軟件中進(jìn)行了系統(tǒng)驗(yàn)證評(píng)估，得到以下主要結(jié)論：

1）對(duì)于多核CPU 平臺(tái)，基于任務(wù)級(jí)并行策略（特別是剖分復(fù)制方法）的數(shù)據(jù)局部性更好，線程的全局同步開(kāi)銷更小，性能更高；相比純MPI并行模式，MPI+OpenMP 的混合并行模式的規(guī)模擴(kuò)展性更好，更加適用于大規(guī)模并行計(jì)算。

2）對(duì)于眾核GPU 平臺(tái)，實(shí)現(xiàn)了3 種基于循環(huán)的并行策略，其中面著色和規(guī)約策略均優(yōu)于原子操作，基于SMEM 優(yōu)化的規(guī)約策略加速效果明顯；相比串行CPU 計(jì)算，單GPU 并行計(jì)算能取得127 倍的加速比。

3）提高數(shù)據(jù)局部性在提升非結(jié)構(gòu)CFD 計(jì)算性能中起著重要性的作用，采用Cuthill-McKee網(wǎng)格重排序以及循環(huán)融合技術(shù)能夠改善數(shù)據(jù)局部性，可以分別使整體性能提升10%；基于SMEM 優(yōu)化的規(guī)約策略對(duì)于有數(shù)據(jù)競(jìng)爭(zhēng)的熱點(diǎn)函數(shù)，性能可提升3倍。

在下一步工作中，將圍繞隱式求解方法開(kāi)展更加深入的共享內(nèi)存并行算法研究，進(jìn)一步提升非結(jié)構(gòu)CFD 軟件的并行計(jì)算效率和擴(kuò)展性。