聲子BTE應用的并行和優化研究*

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

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

2.上海交通大學 密西根學院，上海 200240

1 引言

DGX-2是NVIDIA家族最強大的人工智能系統，基于NVSwitch技術，整合16塊完全互聯的TeslaV100GPU構建的可擴展框架。整個系統顯存達到512 GB，共計擁有40 960個雙精度CUDA（compute unified device architecture）核心。除了多塊GPU卡帶來的強大計算能力，DGX-2中采用的NVSwitch技術提供了GPU之間的高速互聯，NVSwitch 擁有18 個51.5 GB/s 的NVLink 端口。基于其強大的計算能力以及GPU 間的高速通信網絡，DGX-2 在通用計算領域也擁有著巨大潛力。

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

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

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

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

2 相關工作

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

關于使用GPU 進行科學計算方面，在電動力學領域，Priimak[4]在GPU 上實現了描述半導體超晶格中電子傳輸的基于有限差分法的玻爾茲曼輸運方程，在GTX 680上相較于CPU串行版本獲得了118倍的加速。在流體力學領域，Calore等[5]在GPU上實現并優化了可求解大規模湍流的格子玻爾茲曼方法，在K80 上，相較Xeon-Phi 獲得了4 倍加速比，并且最多擴展到32個GPU并行計算，并行效率約90%，總的性能接近20 Tflops。在科學計算領域，Bell 等[6]基于不同稀疏矩陣存儲格式實現其相應的矩陣向量乘，并在結構化網格和非結構化網格上進行測試，相較于4 核CPU 并行，GPU 計算達到了10 倍以上的加速比。Anzt 等[7]研究了預處理的BiCGSTAB（biconjugate gradient stabilized method）、CGS（conjugate gradient squared method）和QMR（quasi-minimal residual）等Krylov求解器在GPU上的效果，證明了在性能方面，相比于ILU預處理器，Jacobi預處理通常會有更高的效率。

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

3 背景介紹

3.1 聲子玻爾茲曼輸運方程

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

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

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

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

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

在假定當前的平衡態分布函數后，該離散方程會變成一個線性方程，而對每個角方向以及整體的區域來說，將會有如下線性方程組：

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

4 聲子BTE方程在DGX-2集群的移植優化

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

4.1 聲子BTE方程求解過程在GPU上的實現

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

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

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

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

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

（4）進行CPU內存與GPU顯存之間數據傳輸，主要數據為每個單元之間的影響系數。

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

（6）求解新的溫度分布和平衡態分布函數。

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

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

4.1.1 線程分配策略

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

4.1.2 數據存儲結構

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

Fig.2 Data storage format圖2 數據存儲結構

4.2 BiCGSTAB算法實現流程

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

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

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

4.2.1 GPU上稀疏矩陣向量乘的實現

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

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

Fig.3 Implementation of BiCGSTAB on GPU圖3 BiCGSTAB算法在GPU上的實現流程

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

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

4.2.2 GPU上點積計算的實現

點積計算需要進行GPU規約計算。為了規約的方便，劃分塊時，塊內線程數為2的指數倍。該方式可使每個線程訪問共享內存時不會發生bankconflict。考慮到GPU 中warp 內線程是同步的，不必進行顯式線程同步。因此當活動線程均在同一warp 時，可直接進行warp 內規約計算，進而減少了線程同步的次數，提高了計算效率。

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

4.2.3 GPU上AXPY實現

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

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

4.3 基于角方向的大規模并行策略

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

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

4.3.1 MPI+CUDA

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

求解聲子BTE 方程的迭代過程中，不同角方向的聲子散射所需數據為上一次迭代計算結果，與本次迭代計算無關。因此不同角方向可以并行執行。不同MPI 進程計算不同角方向，并在GPU 上完成該角方向線性方程組求解。由于在計算聲子模式的溫度分布中，每個網格單元需要對其所有角方向的能量密度進行求和。因此每個進程求得該角方向的能量密度后，需將結果數據從設備端傳輸到主機端，然后通過MPI_Allgather函數使得每個進程均可獲得所有角方向的能量密度；最后將所有角方向能量密度傳輸到GPU中完成對聲子模式溫度的計算。然后繼續計算其他聲子模式下角方向的相關數據，最終求得當前迭代過程中的晶格溫度和平衡態分布。

4.3.2 CUDA-AwareMPI

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

Fig.6 Traditional MPI send and receive scheme圖6 傳統MPI send receive模式

Fig.7 CUDA-AwareMPI implementation圖7 CUDA-AwareMPI實現方式

除此之外，使用CUDA-AwareMPI 能夠使數據傳輸的操作都會被流水化且能夠透明化地使用GPUDirect-RDMA的加速技術。該技術使得GPU中的數據可以直接被送到網卡進行傳輸，從而消除了GPU到CPU設備的時間消耗。這樣也就顯著增大了GPU和其他節點的通信效率。

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

4.3.3 NCCL函數

NCCL[18]是NVIDIA提供的GPU間的集合通信函數庫，實現了all-reduce、all-gather 和reduce-scatter 等集合通信函數。可以在使用PCIe、NVLink 和NVSwitch的平臺上實現高通信速率。其中2.0 版本可以進行跨節點間的GPU 通信，并可進行網絡拓撲檢測以及自動使用GPU DirectRDMA。為了在原并行程序的基礎上使用NCCL函數，采用一個進程分配一個GPU的方式。

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

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

5 實驗結果及分析

5.1 實驗環境

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

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

5.2 GPU加速性能測試

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

Fig.8 Speedup of GPU for different mesh counts圖8 不同網格數目的GPU加速比

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

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

5.3 多GPU并行性能測試

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

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

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

6 結束語

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

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

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

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