聲子BTE應(yīng)用的并行和優(yōu)化研究*

文敏華，劉永志，鮑 華，胡 躍，沈泳星，韋建文，林新華+

1.上海交通大學(xué) 高性能計算中心，上海 200240

2.上海交通大學(xué) 密西根學(xué)院，上海 200240

1 引言

DGX-2是NVIDIA家族最強大的人工智能系統(tǒng)，基于NVSwitch技術(shù)，整合16塊完全互聯(lián)的TeslaV100GPU構(gòu)建的可擴展框架。整個系統(tǒng)顯存達到512 GB，共計擁有40 960個雙精度CUDA（compute unified device architecture）核心。除了多塊GPU卡帶來的強大計算能力，DGX-2中采用的NVSwitch技術(shù)提供了GPU之間的高速互聯(lián)，NVSwitch 擁有18 個51.5 GB/s 的NVLink 端口?；谄鋸姶蟮挠嬎隳芰σ约癎PU 間的高速通信網(wǎng)絡(luò)，DGX-2 在通用計算領(lǐng)域也擁有著巨大潛力。

聲子玻爾茲曼輸運方程（Boltzmann transport equa-tion，BTE）能夠在亞連續(xù)尺度下模擬平衡或非平衡態(tài)的熱傳導(dǎo)現(xiàn)象，并且由于在大范圍長尺度上的有效性持續(xù)受到歡迎[1]。在最近四十年中，半導(dǎo)體技術(shù)快速發(fā)展不僅使得芯片計算能力大大加強，也使得功耗越來越高；當(dāng)前熱問題正成為決定芯片性能、可靠性和成本的主要因素。因此通過聲子BTE方程模擬計算半導(dǎo)體器件的熱傳導(dǎo)對了解其基本物理機制和發(fā)展散熱策略具有重要意義。但是當(dāng)前最為成熟的有限體積法求解聲子BTE方程來模擬工程實際問題時，仍然存在計算量大、計算時間長的困難。

基于將DGX-2強大的計算能力以及高速的通信速率應(yīng)用于聲子BTE 方程求解的初衷，本文做了如下工作：（1）使用CUDA框架對采用有限體積法的聲子BTE求解器的迭代部分進行了針對GPU移植與并行優(yōu)化。在GPU 上實現(xiàn)了迭代計算的主要步驟，如聲子散射項的計算、線性方程組的求解、晶格溫度平衡態(tài)分布的計算等。（2）采用3 種方式，MPI+CUDA，NCCL（NVIDIA collective communications library）函數(shù)以及CUDA-AwareMPI，實現(xiàn)了對聲子BTE方程的多GPU的并行求解?？偠灾?，本文有以下兩點貢獻：

使用CUDA框架首次在GPU上實現(xiàn)了聲子BTE求解的迭代過程，并在單塊V100 上，較Intel Xeon Gold 6248單核性能提高了5.3倍至31.5倍。

在DGX-2 平臺上實現(xiàn)了對聲子BTE 的多GPU并行求解，測試了3 種并行方式對程序性能的影響，其中性能最優(yōu)的NCCL 庫函數(shù)版本在8 臺DGX-2 共128塊V100上實現(xiàn)了83%的強擴展并行效率。

2 相關(guān)工作

有限體積法等確定性方法求解聲子玻爾茲曼輸運方程較為困難。因為只有對角度做到足夠數(shù)量的離散才能夠滿足網(wǎng)格無關(guān)性的驗證，而這就使得該方法會產(chǎn)生數(shù)量龐大的離散方程，對計算量的需求較為龐大。針對這種情況，Ali 等[2]提出了求解聲子BTE的幾種大規(guī)模并行計算的策略和算法：（1）基于聲子模式；（2）基于角方向；（3）混合聲子模式和網(wǎng)格單元的并行方式。對三維器件類硅結(jié)構(gòu)進行非平衡態(tài)的瞬態(tài)模擬，將計算域離散為604 054 個網(wǎng)格單元，使用400 個角度對角方向進行離散，40 個聲子模式對聲子散射曲線進行離散，使用400個進程計算一個時間步長都需要1 h。Ni等[3]也對有限體積法求解聲子BTE 方程進行了并行計算研究，提出了基于空間域和聲子模式的并行策略。其中空間域較為適用全反射性的聲子BTE模型，能夠較好地擴展到128核并行。以上都是對聲子BTE 進行CPU 并行求解，可以看出并行的核數(shù)為數(shù)百核左右，并行規(guī)模并不大。

關(guān)于使用GPU 進行科學(xué)計算方面，在電動力學(xué)領(lǐng)域，Priimak[4]在GPU 上實現(xiàn)了描述半導(dǎo)體超晶格中電子傳輸?shù)幕谟邢薏罘址ǖ牟柶澛斶\方程，在GTX 680上相較于CPU串行版本獲得了118倍的加速。在流體力學(xué)領(lǐng)域，Calore等[5]在GPU上實現(xiàn)并優(yōu)化了可求解大規(guī)模湍流的格子玻爾茲曼方法，在K80 上，相較Xeon-Phi 獲得了4 倍加速比，并且最多擴展到32個GPU并行計算，并行效率約90%，總的性能接近20 Tflops。在科學(xué)計算領(lǐng)域，Bell 等[6]基于不同稀疏矩陣存儲格式實現(xiàn)其相應(yīng)的矩陣向量乘，并在結(jié)構(gòu)化網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格上進行測試，相較于4 核CPU 并行，GPU 計算達到了10 倍以上的加速比。Anzt 等[7]研究了預(yù)處理的BiCGSTAB（biconjugate gradient stabilized method）、CGS（conjugate gradient squared method）和QMR（quasi-minimal residual）等Krylov求解器在GPU上的效果，證明了在性能方面，相比于ILU預(yù)處理器，Jacobi預(yù)處理通常會有更高的效率。

從以上工作看出，GPU 在通用計算領(lǐng)域應(yīng)用較廣，但是目前還沒有在GPU 上對有限體積法求解聲子BTE方程進行相關(guān)研究。本文不僅嘗試在GPU上對該方程進行求解，且使用角方向的并行策略擴展到128塊GPU上進行計算，并取得83%的并行效率。

3 背景介紹

3.1 聲子玻爾茲曼輸運方程

聲子玻爾茲曼輸運方程是一個7 維非線性積分微分方程。目前求解該方程的方法基本有3種：（1）隨機性方法；（2）格子玻爾茲曼方法；（3）確定性方法。隨機性方法如蒙特卡羅法能有效求解這類高維偏微分方程[8]。但它最大的問題在于求解實際工程問題時需要避免結(jié)果的統(tǒng)計波動，但這樣的代價過高[9]。格子玻爾茲曼方法僅用于在簡單的二維結(jié)構(gòu)中求解聲子BTE方程[10-11]。同時，格子玻爾茲曼方法角度與空間離散不能相互獨立，離散方式不合適會極大影響結(jié)果的準(zhǔn)確性[12-13]?；陔x散的確定性方法，如離散坐標(biāo)法（discrete ordinate method，DOM）[14]、有限體積法（finite volume method，F(xiàn)VM）[15]和有限單元法（finite element method，F(xiàn)EM）[16]等廣泛用于聲子BTE的求解。采用離散方法產(chǎn)生的離散方程數(shù)目龐大，對求解的計算量需求也非常龐大。

聲子是晶格振動的能量量子，其能量與頻率有關(guān)。頻率和波矢之間的關(guān)系為色散曲線，在計算中，連續(xù)色散曲線難以直接處理，因此采用兩種簡化的計算模型：灰體模型和非灰模型。灰體模型將所有聲子歸于一個群體；而非灰模型，將色散曲線離散成若干個聲子模式。以下以非灰模型為例介紹聲子玻爾茲曼輸運方程。通過結(jié)合聲子散射項和基于能量的聲子玻爾茲曼輸運方程，可得到如下所示的在弛豫時間近似下的聲子能量玻爾茲曼輸運方程：

3.2 基于有限體積法的聲子BTE求解

采用有限體積法對聲子BTE 方程進行離散，將其變成一系列可求解的低維方程。對角度、計算域及色散曲線進行離散，可得到如下形式的離散方程：

下標(biāo)i和向量S分別代表空間單元和角度方向，因此對于每一個聲子模式的每一個角方向及其每一個空間上的單元來說，采用一階迎風(fēng)格式，以方向s1為例，式（2）將會有一個如下形式的離散方程：

在假定當(dāng)前的平衡態(tài)分布函數(shù)后，該離散方程會變成一個線性方程，而對每個角方向以及整體的區(qū)域來說，將會有如下線性方程組：

系數(shù)矩陣中Kij表示第j個單元對第i個單元的影響系數(shù)。對該線性方程的求解可得到該方向下的能量密度分布。關(guān)聯(lián)角度之間及聲子模式的能量密度，可以求解溫度和平衡能量密度分布函數(shù)。

4 聲子BTE方程在DGX-2集群的移植優(yōu)化

在使用有限體積法求解聲子BTE 方程中，根據(jù)對色散曲線的簡化不同分為灰體與非灰模型，本文采用非灰模型來求解聲子BTE方程。為求解非灰模型穩(wěn)態(tài)的能量密度分布和溫度分布，需要將角度間的能量密度進行關(guān)聯(lián)，不僅需要考慮角度間的影響，還要考慮每個聲子模式間的影響。

4.1 聲子BTE方程求解過程在GPU上的實現(xiàn)

聲子BTE方程在GPU求解流程如圖1所示。算法的關(guān)鍵步驟如下：

（1）讀入網(wǎng)格文件、邊界條件以及參數(shù)文件，獲取所有聲子模式的遷移速度、比熱容和弛豫時間。

（2）假定每個聲子模式上所有方向的能量密度分布的值，初始化為0；計算離散后的角度控制方向。

（3）根據(jù)一階迎風(fēng)格式得到式（4）中每個聲子模式及其所有角方向上的系數(shù)矩陣。

Fig.1 Process of solving non-gray phonon BTE model圖1 聲子BTE非灰模型求解過程

（4）進行CPU內(nèi)存與GPU顯存之間數(shù)據(jù)傳輸，主要數(shù)據(jù)為每個單元之間的影響系數(shù)。

（5）計算聲子散射項，使用BiCGSTAB 算法求解線性方程組，求解每個網(wǎng)格單元的能量密度分布。

（6）求解新的溫度分布和平衡態(tài)分布函數(shù)。

（7）計算能量密度分布新舊值之差，驗證是否滿足收斂條件。若滿足收斂條件則停止計算，輸出結(jié)果，否則繼續(xù)（5）、（6）步驟。

由于整個迭代過程均在GPU 上實現(xiàn)，因此迭代過程所需數(shù)據(jù)一次性傳入顯存中，之后在GPU 上完成迭代過程中數(shù)據(jù)的更新，這樣減少了GPU 與CPU間的數(shù)據(jù)傳輸，提高了計算效率。通過上述步驟得到穩(wěn)態(tài)的能量密度的分布和溫度分布，之后可以求得熱流量，從而可以得到如熱導(dǎo)率這樣的宏觀特性，完成對介觀尺度下傳熱問題的模擬。

4.1.1 線程分配策略

在GPU 上實現(xiàn)的迭代計算部分的函數(shù)，如聲子散射項函數(shù)、聲子模式溫度函數(shù)、晶格溫度函數(shù)以及平衡態(tài)分布函數(shù)，不同網(wǎng)格單元間沒有數(shù)據(jù)依賴。因此可使每個CUDA 線程對應(yīng)一個相應(yīng)的網(wǎng)格單元，GPU并行計算過程中，每個線程利用上一次迭代計算的結(jié)果或者是上一階段計算函數(shù)輸出的結(jié)果，完成相應(yīng)網(wǎng)格單元中的數(shù)據(jù)的計算，進而完成整個迭代過程的計算。

4.1.2 數(shù)據(jù)存儲結(jié)構(gòu)

聲子BTE 方程維度較高，迭代計算過程所需的數(shù)據(jù)，如能量密度值，除空間維度，還額外擁有聲子模式以及角方向這兩個維度。其中空間維度使用網(wǎng)格編號表示，因此數(shù)據(jù)一般存儲在三維數(shù)組當(dāng)中。由于GPU 函數(shù)是以網(wǎng)格單元進行的CUDA 線程分配，因此計算過程中相鄰線程間訪問的是數(shù)組中同一聲子模式和角方向下，相鄰網(wǎng)格單元間的數(shù)據(jù)。為了能夠利用GPU合并訪存特性，簡化CPU、GPU端數(shù)據(jù)傳輸，數(shù)據(jù)采用如圖2所示的一維數(shù)組進行存儲。

Fig.2 Data storage format圖2 數(shù)據(jù)存儲結(jié)構(gòu)

4.2 BiCGSTAB算法實現(xiàn)流程

上一節(jié)描述的算法中最關(guān)鍵和耗時的步驟是通過求解離散后的聲子BTE線性方程來確定聲子能量密度分布（4.1節(jié)步驟5）。因為必須針對每個角方向和聲子模式求解聲子BTE 線性方程組，每一次迭代需進行大量線性方程組求解。在本研究中，使用預(yù)處理穩(wěn)定雙共軛梯度法（BiCGSTAB）[17]求解線性方程。該算法如下所示：

Jacobi預(yù)處理方法易于在GPU并行計算，相較于ILU 預(yù)處理方法有更好的性能[7,17]，因此本文采用Jacobi 預(yù)處理方式，該方式為最簡單的預(yù)處理形式，其預(yù)處理矩陣為原矩陣的對角矩陣。

整個求解過程流程如圖3所示。算法的主要部分是稀疏矩陣向量乘（sparse matrix-vector multiplication，SpMV）、點積（Dot）以及標(biāo)量向量乘和向量加（alpha X plus Y，AXPY）。其中CPU 負(fù)責(zé)求解過程的流程控制，核心計算部分在GPU 中進行。上述算法在計算過程中，若出現(xiàn)關(guān)于r0的點積值為零，需要對x0、r0以及其他參數(shù)如ρ、α、ω等重新進行初始化并繼續(xù)進行迭代計算。因此在計算過程中，需要將一些參數(shù)傳輸?shù)絻?nèi)存，使用CPU 判斷是否需要重新初始化以及計算是否收斂。

4.2.1 GPU上稀疏矩陣向量乘的實現(xiàn)

SpMV 在線性方程組求解算法中占有重要的地位，許多學(xué)者對其進行了研究。Bell 等[6,11]在GPU 實現(xiàn)了不同稀疏矩陣格式的SpMV計算，發(fā)現(xiàn)SpMV的性能主要受稀疏矩陣存儲格式以及核函數(shù)線程分配方式的影響。

由于本文研究的應(yīng)用中網(wǎng)格單元的系數(shù)，如3.2節(jié)式（4）所示，僅受相鄰3個網(wǎng)格單元的影響，即每行影響系數(shù)最多為3 個，因此整個系數(shù)矩陣為稀疏矩陣。由于ELL（Ellpack）格式適用于分布較為均衡的稀疏矩陣，且SpMV 實現(xiàn)方式簡單，因此本文采用ELL格式稀疏矩陣。

Fig.3 Implementation of BiCGSTAB on GPU圖3 BiCGSTAB算法在GPU上的實現(xiàn)流程

Fig.4 ELL sparse matrix storage format圖4 ELL稀疏矩陣存儲格式

ELL格式存儲方式如圖4所示，采用兩個二維數(shù)組來存儲一個n×k的矩陣（k為包含非零元素最多行的非零元素數(shù)目）。在實際使用中，使用兩個一維向量按列方向?qū)仃囘M行存儲。使用ELL 格式的SpMV 算法在GPU 內(nèi)易于并行實現(xiàn)，計算過程中每個CUDA 線程計算矩陣的一行，由ELL 存儲方式可知，CUDA線程對矩陣值以及列索引的訪問均是連續(xù)的，能夠充分利用GPU合并訪存特性。

4.2.2 GPU上點積計算的實現(xiàn)

點積計算需要進行GPU規(guī)約計算。為了規(guī)約的方便，劃分塊時，塊內(nèi)線程數(shù)為2的指數(shù)倍。該方式可使每個線程訪問共享內(nèi)存時不會發(fā)生bankconflict?？紤]到GPU 中warp 內(nèi)線程是同步的，不必進行顯式線程同步。因此當(dāng)活動線程均在同一warp 時，可直接進行warp 內(nèi)規(guī)約計算，進而減少了線程同步的次數(shù)，提高了計算效率。

由于采用多個塊（Block）進行點積計算，因此整個規(guī)約的過程分為兩個階段：第一階段采用多個Block，首先計算出向量中每個值的乘積，之后進行塊內(nèi)規(guī)約計算，并將當(dāng)前塊內(nèi)規(guī)約結(jié)果順序存入中間數(shù)組中；第二階段僅使用1個Block，對中間數(shù)組進行規(guī)約求和，規(guī)約方式同上，最終可求得向量點積的值。

4.2.3 GPU上AXPY實現(xiàn)

在GPU 上并行計算AXPY 較為簡單，為向量中每個數(shù)據(jù)分配一個線程即可。同時可將前后步驟間有數(shù)據(jù)依賴且并行方式相同的函數(shù)合并到一個CUDA 內(nèi)核中，這樣就可避免雙倍的數(shù)據(jù)加載，減少了內(nèi)核數(shù)量，節(jié)約內(nèi)核調(diào)用開銷。例如在求解線性方程組中，會有如下所示的相鄰計算。

其中，K為Jacobi 預(yù)處理矩陣，實際計算過程為向量乘，并行過程與AXPY 一致，因此可以將兩步操作融合成一個CUDA 函數(shù)。該方式不僅簡化了代碼，而且更為重要的是利用了數(shù)據(jù)的局部性，每一個線程計算的結(jié)果可直接應(yīng)用到下一階段計算中。

4.3 基于角方向的大規(guī)模并行策略

聲子玻爾茲曼輸運方程通常有3 種并行策略：（1）計算域分解，該方式能做到較大規(guī)模擴展，但是本文應(yīng)用針對的是二維計算域，網(wǎng)格規(guī)模較小，若將其分解成更小的網(wǎng)格，則難以發(fā)揮GPU 多線程的優(yōu)勢；（2）基于聲子模式的并行，該方式雖然較為容易實現(xiàn)，但是聲子模式的數(shù)量一般為數(shù)十個，不能做到較大規(guī)模的并行；（3）基于角方向的并行，角方向的數(shù)量比聲子模式大一個數(shù)量級，采用角方向并行能夠進行數(shù)百個進程并行計算。因此為了做到較大規(guī)模的并行計算，本研究采用基于角方向的并行策略。本文使用3種不同方法實現(xiàn)了該并行策略，其中第一種實現(xiàn)方式如圖5所示。

Fig.5 Parallel strategy based on angle directions圖5 基于角方向的并行策略

4.3.1 MPI+CUDA

如圖5 所示使用的是MPI+CUDA 的并行策略，即一個MPI 進程控制一個GPU 設(shè)備。在初始化階段，獲取MPI進程號rank，以及節(jié)點上GPU設(shè)備數(shù)量gpu_count，并將MPI進程映射到GPU設(shè)備上gpu_id=rank%gpu_count。在DGX-2 節(jié)點上，一個節(jié)點有16個V100，因此一個節(jié)點分配16 個MPI 進程，每個進程對應(yīng)一個GPU設(shè)備。

求解聲子BTE 方程的迭代過程中，不同角方向的聲子散射所需數(shù)據(jù)為上一次迭代計算結(jié)果，與本次迭代計算無關(guān)。因此不同角方向可以并行執(zhí)行。不同MPI 進程計算不同角方向，并在GPU 上完成該角方向線性方程組求解。由于在計算聲子模式的溫度分布中，每個網(wǎng)格單元需要對其所有角方向的能量密度進行求和。因此每個進程求得該角方向的能量密度后，需將結(jié)果數(shù)據(jù)從設(shè)備端傳輸?shù)街鳈C端，然后通過MPI_Allgather函數(shù)使得每個進程均可獲得所有角方向的能量密度；最后將所有角方向能量密度傳輸?shù)紾PU中完成對聲子模式溫度的計算。然后繼續(xù)計算其他聲子模式下角方向的相關(guān)數(shù)據(jù)，最終求得當(dāng)前迭代過程中的晶格溫度和平衡態(tài)分布。

4.3.2 CUDA-AwareMPI

如4.3.1 小節(jié)所述的常規(guī)的MPI 實現(xiàn)，需首先將數(shù)據(jù)傳輸至CPU內(nèi)存中，之后使用MPI函數(shù)；而使用CUA-AwareMPI 實現(xiàn)，可直接傳輸GPU 中的數(shù)據(jù)。從圖6、圖7 兩種實現(xiàn)方式可以看出，使用CUDA-AwareMPI的實現(xiàn)，其代碼更簡潔，編程更容易。

Fig.6 Traditional MPI send and receive scheme圖6 傳統(tǒng)MPI send receive模式

Fig.7 CUDA-AwareMPI implementation圖7 CUDA-AwareMPI實現(xiàn)方式

除此之外，使用CUDA-AwareMPI 能夠使數(shù)據(jù)傳輸?shù)牟僮鞫紩涣魉夷軌蛲该骰厥褂肎PUDirect-RDMA的加速技術(shù)。該技術(shù)使得GPU中的數(shù)據(jù)可以直接被送到網(wǎng)卡進行傳輸，從而消除了GPU到CPU設(shè)備的時間消耗。這樣也就顯著增大了GPU和其他節(jié)點的通信效率。

使用CUDA-AwareMPI 與4.3.1 小節(jié)所述的并行方式類似，僅需進行兩點改動即可，一是省略CPU和GPU 間數(shù)據(jù)通信的過程，二是將MPI 函數(shù)傳輸?shù)臄?shù)據(jù)改為GPU緩沖區(qū)數(shù)據(jù)。

4.3.3 NCCL函數(shù)

NCCL[18]是NVIDIA提供的GPU間的集合通信函數(shù)庫，實現(xiàn)了all-reduce、all-gather 和reduce-scatter 等集合通信函數(shù)?？梢栽谑褂肞CIe、NVLink 和NVSwitch的平臺上實現(xiàn)高通信速率。其中2.0 版本可以進行跨節(jié)點間的GPU 通信，并可進行網(wǎng)絡(luò)拓?fù)錂z測以及自動使用GPU DirectRDMA。為了在原并行程序的基礎(chǔ)上使用NCCL函數(shù)，采用一個進程分配一個GPU的方式。

NCCL庫的使用方式為，首先完成對MPI的初始化，之后在0 號進程上創(chuàng)建unique ID 并廣播到其他進程，最后創(chuàng)建NCCL 的通信子。上述步驟完成后，即可使用NCCL中的集合通信函數(shù)。

使用NCCL 函數(shù)對聲子BTE 進行并行求解，僅需在4.3.1小節(jié)所述方式上進行如下3點修改：（1）增加NCCL相關(guān)初始化函數(shù)；（2）略去CPU與GPU數(shù)據(jù)傳輸過程；（3）使用ncclAllGather函數(shù)進行集合通信。

5 實驗結(jié)果及分析

5.1 實驗環(huán)境

本研究實驗所采用硬件配置如下，CPU 為Intel Xeon Gold 6248,架構(gòu)為Cascade Lake，主頻為2.5 GHz，核心數(shù)為20，雙精度浮點性能為1.6 Tflops。節(jié)點內(nèi)存大小為192 GB；GPU版本測試使用的是NVIDIA的DGX-2平臺，其中該平臺采用的GPU為V100，5 120個單精度CUDA核心，2 560個雙精度CUDA核心，雙精度浮點性能達到7.5 Tflops，顯存帶寬可達900 GB/s；DGX-2 內(nèi)部通過NVSwtich 橋接器支持16 個全互聯(lián)的V100 GPU卡。

性能測試對比采用CPU-GPU 并行計算和僅采用CPU 串行計算的加速效果。其中GPU 計算采用的CUDA 框架版本為9.2，串行CPU版本中進行線性方程求解采用的是PETSc（portable，extensible toolkit for scientific computation）[19]中函數(shù)，其中PETSc版本為3.10。對4個算例進行了GPU加速性能的測試，其中每個測試算例均為二維計算域，且計算域大小、邊界條件、遷移速度和弛豫時間等參數(shù)均相同，網(wǎng)格及網(wǎng)格數(shù)量不同。算例網(wǎng)格數(shù)量分別為6 282、11 974、18 060及24 912。其中網(wǎng)格數(shù)量為24 912的算例，在實際二維問題的研究中屬于規(guī)模較大的算例。串行版本中不同網(wǎng)格算例內(nèi)存使用量分別為0.90 GB、1.85 GB、3.68GB、6.70GB，GPU顯存使用量約為0.30GB、0.81GB、1.00 GB、1.24 GB。內(nèi)存數(shù)據(jù)量較大主要是因為CPU計算過程中需要使用二維矩陣，并將其壓縮為ELL格式稀疏矩陣，而GPU顯存中僅存儲壓縮后的ELL稀疏矩陣。聲子BTE方程求解過程中，離散聲子色散曲線的聲子模式個數(shù)為1，離散的角方向的個數(shù)為256。

5.2 GPU加速性能測試

圖8 所示為使用單卡V100 相較于Intel Xeon Gold 6248串行版本的迭代過程的加速效果。縱坐標(biāo)為單次迭代計算的時間，橫坐標(biāo)為不同網(wǎng)格大小的算例。由于本文旨在在GPU上實現(xiàn)BTE迭代求解過程，因此僅比較了GPU 版本與CPU 串行版本間的加速性能，未與CPU并行進行對比。

Fig.8 Speedup of GPU for different mesh counts圖8 不同網(wǎng)格數(shù)目的GPU加速比

圖8中顯示的加速比針對的是迭代部分的時間，這是因為整個迭代的步數(shù)在103至104量級之間，迭代過程之外部分的耗時相對較少，對加速比的影響較小，因此僅比較迭代過程的加速比。迭代過程的CPU-GPU 混合計算需考慮數(shù)據(jù)傳輸?shù)拈_銷，本文的實現(xiàn)方式是在迭代計算前將GPU迭代所需數(shù)據(jù)全部傳入顯存中，迭代過程中的數(shù)據(jù)在GPU 顯存中獲取并計算更新。因此迭代過程中，GPU 僅需將與計算流程有關(guān)的系數(shù)傳入內(nèi)存。迭代過程中主機端與設(shè)備端數(shù)據(jù)傳輸為迭代時間的1.1%～1.2%。這是由于迭代過程中傳輸?shù)臄?shù)據(jù)多為單個參數(shù)，即單個浮點數(shù)，數(shù)據(jù)量小。

從圖8 中可以看出，隨著測試算例網(wǎng)格數(shù)量增加，GPU加速效果越來越明顯。算例1中網(wǎng)格數(shù)量為6 282，算例4為24 912，算例4網(wǎng)格數(shù)量約為算例1的4 倍，其串行CPU 版本迭代時間為算例1 的15.6 倍，而其GPU 版本迭代時間僅為算例1 的2.6 倍。這是由于GPU 并行過程中，多數(shù)函數(shù)中CUDA 線程數(shù)與網(wǎng)格數(shù)有關(guān)，而當(dāng)網(wǎng)格數(shù)量較小時，其使用的CUDA線程（thread）數(shù)也較少，未能充分利用GPU 的性能。同時GPU可以通過線程塊的切換掩蓋數(shù)據(jù)的訪存延遲，因此當(dāng)網(wǎng)格數(shù)量較大，線程塊數(shù)量較多時，可以充分利用該特性，提高并行效率，因而在本例中GPU對越大規(guī)模數(shù)據(jù)進行并行加速，其效果越明顯。

5.3 多GPU并行性能測試

采用4.3節(jié)基于角方向的并行策略，對多GPU并行強擴展性能進行測試。由于本文旨在實現(xiàn)聲子BTE 程序的GPU 加速以及實現(xiàn)多GPU 版本的聲子BTE 程序，因此多GPU 并行的加速比是以單GPU 迭代時間為基準(zhǔn)，并未與CPU 并行進行對比。測試使用算例為上述網(wǎng)格大小為24 912 的算例。每個DGX-2 節(jié)點有16 塊V100GPU，共有8 個節(jié)點。最多擴展到128個進程，其強擴展性能如圖9所示。

Fig.9 Strong scalability for different parallel methods圖9 不同并行方式的強可擴展性

圖9 中橫坐標(biāo)為GPU 個數(shù)，一個MPI 進程對應(yīng)一個GPU設(shè)備，圖中標(biāo)注出了并行效率最好的NCCL庫函數(shù)方式，以及并行效率最差MPI+CUDA 方式的加速比。其中使用MPI+CUDA 的版本，擴展到128卡上，相較于1 塊GPU 卡，其加速比為68，其并行效率為53%；使用CUDA-AwareMPI 技術(shù)的并行版本，128卡上加速比為97，并行效率為76%；使用NVIDIA的NCCL庫實現(xiàn)并行的版本，128卡上加速比為107，其并行效率為83%。采用NCCL 庫以及CUDA-AwareMPI的方式，在同一節(jié)點可利用NVSwitch實現(xiàn)GPU間直接高效通信，節(jié)點間通信則利用了GPUDirect-RDMA技術(shù)。而MPI+CUDA方式，首先需要通過PCIe（peripheral component interconnect express）將GPU 中的數(shù)據(jù)傳輸至CPU，之后進行節(jié)點內(nèi)通信以及將數(shù)據(jù)傳入IB（infiniBand）網(wǎng)卡進行節(jié)點間通信。由于通信方式的不同，使得使用NCCL函數(shù)相比MPI+CUDA實現(xiàn)方式性能提升57%。由此可以看出，基于NVLink和NVSwitch 技術(shù)的GPU 間通信相比傳統(tǒng)方式更快速和高效。

6 結(jié)束語

本文提出了采用非結(jié)構(gòu)化網(wǎng)格的有限體積法求解聲子玻爾茲曼輸運方程（BTE）的GPU并行加速方法，并在8臺DGX-2集群上進行了移植和并行優(yōu)化。在GPU上實現(xiàn)了整個迭代過程，減少了主機端與設(shè)備端的數(shù)據(jù)傳輸，同時在求解線性方程組時，通過規(guī)約過程的循環(huán)展開以及內(nèi)核融合等優(yōu)化方法，在單塊V100GPU上，網(wǎng)格單元為2.4萬的算例，相較于Intel Xeon Gold 6248上的串行版本獲得了31.5倍的性能提升。

同時，本文實現(xiàn)和評估了3 種多節(jié)點多GPU 并行的方法。采用MPI+CUDA的并行方式并行效率最低，而采用CUDA-AwareMPI以及使用NCCL函數(shù)均能充分利用NVSwitch和GPUDirectRDMA帶來的通信性能增益。其中使用NCCL 庫函數(shù)的并行方法最為高效，在8 臺DGX-2（128 塊NVIDIAV100）上獲得83%的并行效率，比純MPI版本提升達57%。

目前僅實現(xiàn)了基于角方向的并行策略，該方法需對網(wǎng)格能量密度進行集合通信，由于當(dāng)前網(wǎng)格規(guī)模并不大，因此對性能影響較小。若對三維問題并行求解，網(wǎng)格規(guī)模會提升一個數(shù)量級，則其對性能的影響會較為明顯。因此未來會考慮實現(xiàn)基于聲子模式的并行，該方式通信頻率較少，但是會有負(fù)載不均衡問題，因此需采取相應(yīng)策略進行優(yōu)化。此外本研究中僅采用Jacobi 預(yù)處理方式，雖然其較容易并行，但是收斂較差，因此接下來考慮測試不同的預(yù)處理矩陣的性能。同時將會測試更多的算例并且嘗試使用不同的矩陣格式與線性方程求解器來優(yōu)化單塊GPU上的性能。

致謝感謝上海交通大學(xué)高性能計算中心程盛淦老師在DGX-2設(shè)備的使用上提供的幫助。