易培淮，李衛(wèi)東，林 韜，鄒佳恒，鄧子艷，劉 言

（1.中國科學(xué)院大學(xué) 核科學(xué)與技術(shù)學(xué)院，北京 100049；2.中國科學(xué)院高能物理研究所，北京 100049）

0 概述

江門中微子實(shí)驗(yàn)（JUNO）［1-2］通過探測核電站反應(yīng)堆產(chǎn)生的電子反中微子來測量中微子質(zhì)量等級［2］和精確測量中微子振蕩參數(shù)，具有豐富的物理目標(biāo)。用于探測中微子的中心探測器是世界上能量精度最高、規(guī)模最大的液體閃爍體探測器［3］。宇宙線繆子作為JUNO 實(shí)驗(yàn)重要的本底之一［4-5］，需要經(jīng)過精細(xì)重建后才能夠排除，而精細(xì)的重建算法需要大量的模擬事例作為輸入［6-7］。但繆子穿過探測器時(shí)會產(chǎn)生百萬量級的閃爍光和契侖科夫光，導(dǎo)致基于Geant4［8］的全模擬計(jì)算量巨大。為此，文獻(xiàn)［9］提出一種快速模擬方法，首先用Geant4 模擬預(yù)先產(chǎn)生PMT 的擊中光子數(shù)和擊中時(shí)間的分布，然后在快速模擬時(shí)讀入這些分布直方圖并對其進(jìn)行抽樣，用抽樣的方法代替耗時(shí)的光子傳輸模擬。盡管該方法可使單個(gè)事例模擬時(shí)間從約1 h 縮減至10 min 左右，但傳統(tǒng)的串行計(jì)算方式仍然需要消耗較長的時(shí)間。

基于異構(gòu)的GPU 應(yīng)用［10-12］以及跨節(jié)點(diǎn)的信息傳遞接口（Message Passing Interface，MPI）［13］被廣泛地應(yīng)用于高性能計(jì)算，并獲得數(shù)倍至數(shù)萬倍的加速比［14-16］。在國際高能物理模擬領(lǐng)域，目前出現(xiàn)了利用多線程加速模擬計(jì)算，比較典型的是歐洲核子中心大型強(qiáng)子對撞機(jī)（LHC）上的ATLAS 實(shí)驗(yàn)?zāi)M軟件［17］，以及利用多節(jié)點(diǎn)進(jìn)行模擬的Alice 實(shí)驗(yàn)?zāi)M軟件［18］。此外，各個(gè)高能物理實(shí)驗(yàn)開發(fā)了快速模擬軟件［19-20］來提升模擬速度，這些快速模擬軟件都只能在CPU 上運(yùn)行。如何運(yùn)用GPU 并行Geant4 模擬［21］以及將快速模擬軟件移植到GPU 上運(yùn)行仍處于初步研究階段，到目前為止尚無發(fā)現(xiàn)成功應(yīng)用案例。

繆子在探測器中產(chǎn)生的閃爍光具有各向同性的屬性，并且模擬過程相互獨(dú)立，適合在GPU 上并行。本文提出一種基于GPU 的繆子快速模擬方法，將閃爍光的產(chǎn)生和模擬置于GPU 中進(jìn)行并行計(jì)算，以提高運(yùn)算速度。對于多節(jié)點(diǎn)GPU 集群，該方法采用MPI 分布式通信在GPU 集群中實(shí)現(xiàn)光子模擬分布式計(jì)算。

1 繆子快速模擬算法

JUNO 模擬軟件由物理產(chǎn)生子、探測器模擬、電子學(xué)模擬3 個(gè)部分構(gòu)成，如圖1 所示。物理產(chǎn)生子用于產(chǎn)生原初粒子信息，之后再由探測器模擬對原初粒子進(jìn)行運(yùn)輸和跟蹤。模擬軟件基于Geant4 和SNiPER［22］開發(fā)，用于模擬粒子和探測器之間的相互作用，并記錄探測器中的響應(yīng)信息，最后由電子學(xué)模擬對擊中信息數(shù)字化，生成與真實(shí)實(shí)驗(yàn)相同的波形數(shù)據(jù)。JUNO 實(shí)驗(yàn)采用1 GHz 的采樣率采集波形，電子學(xué)模擬中按照1 ns 間隔讀出最終的波形數(shù)據(jù)。

圖1 JUNO 模擬軟件示意圖Fig.1 Schematic diagram of JUNO simulation software

快速模擬方法將繆子徑跡分割成很多份短的步長（Step），如圖2 所示。對于每個(gè)Step，需要根據(jù)位置和沉積能量產(chǎn)生探測器的響應(yīng)。Step 所在位置可以確定徑向位置R以及該位置與每個(gè)PMT 之間的夾角θ，可以唯一確定需要加載PMT 的光電子（nPE）直方圖和擊中時(shí)間（Hit time）直方圖數(shù)據(jù)?？焖倌M的第1 步是根據(jù)每個(gè)Step 的可見能量對nPE 直方圖進(jìn)行抽樣，產(chǎn)生該P(yáng)MT 上的nPE 數(shù)目。第2 步是為每個(gè)PE 產(chǎn)生一個(gè)擊中時(shí)間。Hit time 直方圖是使用初始時(shí)間為0 的事例構(gòu)建的，包括光子產(chǎn)生時(shí)間和光子飛行到達(dá)PMT 的時(shí)間，最終的擊中時(shí)間等于Step 當(dāng)前的時(shí)間加上抽樣得出時(shí)間。

圖2 繆子快速模擬幾何關(guān)系示意圖Fig.2 Schematic diagram of Muon fast simulation geometrical relationship

快速模擬的抽樣計(jì)算需要預(yù)先加載36 000 份nPE 直方圖和Hit time 直方圖，才能夠?qū)崿F(xiàn)精確的快速模擬。快速模擬也需要加載約18 000 支PMT 的位置坐標(biāo)，用于計(jì)算夾角θ。以上的直方圖數(shù)據(jù)和PMT 坐標(biāo)數(shù)據(jù)都可以在模擬初始化階段加載完成。

2 GPU 繆子快速模擬

模擬繆子在穿過探測器時(shí)，Geant4 會產(chǎn)生繆子的步長信息，快速模擬方法對每個(gè)步長進(jìn)行計(jì)算得出每個(gè)PMT 上的擊中數(shù)目和時(shí)間。針對模擬過程中步長和PMT 數(shù)目十分龐大的挑戰(zhàn)，提出一種GPU繆子快速模擬方法，并完成了基于CUDA［23］的快速模擬并行軟件。該方法的核心思想是：首先將繆子的步長集合拆分成子集，分配給不同的GPU 卡處理，實(shí)現(xiàn)GPU 卡間的并行；接著在每一個(gè)GPU 卡上，把原來CPU 串行遍歷PMT 改成按照GPU 線程并行模擬PMT。GPU 快速模擬方法通過上述2 個(gè)層次的并行加快了運(yùn)行速度。

GPU 繆子快速模擬方法在初始化階段將預(yù)生成的直方圖數(shù)據(jù)經(jīng)過序列化存儲于GPU 全局內(nèi)存中，供后續(xù)內(nèi)核算法訪問。在事例循環(huán)階段，該方法首先從Geant4 收集步長信息，并進(jìn)行實(shí)時(shí)序列化，用于后續(xù)傳輸至GPU 內(nèi)存。該方法使用多GPU 及MPI進(jìn)行并行模擬，對步長信息自動劃分，并將劃分后的數(shù)據(jù)傳輸至對應(yīng)的GPU 全局內(nèi)存，解決了因單塊GPU 顯存的限制無法處理所有步長信息的問題，從而實(shí)現(xiàn)大規(guī)模并行。該方法實(shí)現(xiàn)的CUDA Kernel 函數(shù)將一個(gè)PMT 的模擬映射到一個(gè)CUDA 線程，通過啟動大量的線程完成大規(guī)模的模擬計(jì)算，并且將中間結(jié)果存儲于GPU 全局內(nèi)存。在每個(gè)事例終止時(shí)，合并所有GPU 的數(shù)據(jù)，并將結(jié)果傳回CPU，從而完成一個(gè)事例的模擬。GPU 繆子快速模擬的Kernel 函數(shù)主要負(fù)責(zé)光子的模擬，用于替換串行的CPU 算法。CPU 算法中對于每個(gè)Step 都需要遍歷20 000 個(gè)PMT 的計(jì)算，GPU 算法則將20 000 個(gè)PMT 的計(jì)算遷移至GPU 并行，極大地縮短了運(yùn)行時(shí)間。

GPU 繆子快速模擬算法的具體實(shí)現(xiàn)包括1 個(gè)運(yùn)行于CPU 的函數(shù)和2 個(gè)運(yùn)行于GPU 的CUDA Kernel函數(shù)，如算法1 所示。

在算法1 中，第1 行～第8 行是CPU 部分的代碼實(shí)現(xiàn)，主要負(fù)責(zé)CPU 與GPU 之間的數(shù)據(jù)拷貝以及CUDA Kernel 函數(shù)的調(diào)用。對于每個(gè)物理事例，收集Geant4產(chǎn)生的步長信息，然后拷貝至GPU，對應(yīng)于第3 行、第4行。隨后，調(diào)用CUDA Kernel函數(shù)，通過block和thread 2 個(gè)參數(shù)設(shè)置并發(fā)線程數(shù)，且總的并發(fā)數(shù)不少于PMT數(shù)目，對應(yīng)于第5 行。模擬完成后，再將結(jié)果合并拷貝回CPU，對應(yīng)于第6 行、第7 行。

第9 行～第18 行是CUDA Kernel 函數(shù)的實(shí)現(xiàn)，在GPU 上并行執(zhí)行，主要負(fù)責(zé)PMT 的模擬。由于總的并發(fā)線程數(shù)將大于PMT 數(shù)目，對于超過PMT 數(shù)目的線程，設(shè)置不執(zhí)行后續(xù)的計(jì)算，對應(yīng)于第10 行～第12 行。通過CUDA 提供的線程索引，計(jì)算出當(dāng)前線程將要模擬的PMT ID，對應(yīng)于第13 行。第14 行～第17 行對傳入的步長信息進(jìn)行循環(huán)，并調(diào)用另一個(gè)Kernel 函數(shù)SimulationHits 進(jìn)行光子的模擬。模擬產(chǎn)生的結(jié)果被存儲于GPU 全局內(nèi)存，用于后續(xù)返回。

第19 行～第37 行對應(yīng)于SimulationHits 的實(shí)現(xiàn)，該函數(shù)的輸入?yún)?shù)為PMT ID 以及步長信息，返回的結(jié)果為當(dāng)前步長生成的光子數(shù)目以及對應(yīng)的擊中時(shí)間。該函數(shù)從步長的位置計(jì)算得出徑向距離R；再從PMT 與當(dāng)前步長的位置關(guān)系得出角度θ?；诰嚯xR和角度θ2 個(gè)變量，從全局內(nèi)存中分別定位到nPE 和Hit Time 2 套直方圖數(shù)據(jù)，對應(yīng)于第21 行。第22 行～第35 行實(shí)現(xiàn)了從2 套直方圖分別抽樣nPE和hit time。由于直方圖是基于1MeV 能量產(chǎn)生的，算法中將當(dāng)前步長的能量分成整數(shù)energyInt 和小數(shù)energyFraction 2 個(gè)部分：對于整數(shù)部分，重復(fù)進(jìn)行抽樣產(chǎn)生hit time；對于小數(shù)部分，若產(chǎn)生的隨機(jī)數(shù)小于小數(shù)部分，再產(chǎn)生hit time，對應(yīng)于第31 行～第34 行。

對于偽代碼中SimulationHits 算法所描述的直方圖抽樣的過程，即代碼第26 行和第32 行中的getRandom，實(shí)現(xiàn)過程如下：首先在初始化加載直方圖數(shù)據(jù)時(shí)，將概率分布直方圖轉(zhuǎn)化為累計(jì)分布函數(shù)直方圖，直方圖在全局內(nèi)存中以定長數(shù)組的形式保存。GPU 線程通過Step 信息中的R和θ找到對應(yīng)的累積分布直方圖后，調(diào)用CUDA 隨機(jī)數(shù)算法cuRAND 生成在范圍為0～1 的概率值，遍歷直方圖數(shù)組中找到剛好不小于概率值的元素，元素的下標(biāo)即為抽樣的結(jié)果。

在并行模擬時(shí)，線程束（Warp）內(nèi)的32 個(gè)線程從GPU 全局內(nèi)存讀取直方圖數(shù)據(jù)用于抽樣計(jì)算，如圖3 所示。直方圖在內(nèi)存的放置策略按照距離R和角度θ2 級索引的方式，保證Warp 中的線程能夠訪問連續(xù)內(nèi)存空間的直方圖：第1 級索引為距離R，長度為200；第2 級索引為角度θ，長度為180。在圖3中，32 個(gè)線程訪問的直方圖具有相同的距離R和相近的角度θ，能夠合并訪問GPU 全局內(nèi)存空間，從而提高數(shù)據(jù)訪問效率。

圖3 GPU 內(nèi)存優(yōu)化策略Fig.3 GPU memory optimization strategy

除訪問直方圖數(shù)據(jù)外，GPU 線程還要訪問步長信息。預(yù)先產(chǎn)生的直方圖數(shù)據(jù)數(shù)目固定，而步長的數(shù)目則與模擬事例的能量相關(guān)，能量越高產(chǎn)生的步長數(shù)目越多，所以對步長信息的實(shí)現(xiàn)與直方圖數(shù)據(jù)不同。針對CUDA 多次啟動核函數(shù)開銷高的問題，GPU 算法并不對每個(gè)Step 單獨(dú)調(diào)用快速模擬的方法，而是創(chuàng)建步長信息的集合，將Step 信息集中收集，待Step 收集完成后，再一次性或分批調(diào)用并行化快速模擬算法，從而提高算法效率。通過引入步長集合的設(shè)計(jì)，將一個(gè)高能的繆子事例拆分成多批進(jìn)行模擬，也可將多個(gè)低能的繆子事例合并成一批進(jìn)行模擬。

繆子快速模擬并行軟件以面向?qū)ο蠓椒ㄩ_發(fā)，并與JUNO 離線軟件系統(tǒng)整合，如圖4 所示。以SNiPER和Geant4 作為底層軟件，繆子快速模擬軟件包括CPUFastSim、Parallel-FastSim 和Helper 3 個(gè)組件，分別實(shí)現(xiàn)CPU串行模擬、GPU并行模擬和輔助工具。Helper主要提供從文件加載直方圖等通用功能。在Parallel-FastSim組件中，除了運(yùn)行于GPU的Kernel函數(shù)和Buffer以CUDA 實(shí)現(xiàn)，其他的類都派生自SNiPER 框架，并以C++實(shí)現(xiàn)運(yùn)行于CPU 中。負(fù)責(zé)GPU 快速模擬的GPUFastSimTool 派生自ToolBase，與串行版本CPUFastSimTool 具有相同的接口，能夠在Geant4 的EventAction 與SteppingAction 中被調(diào)用，軟件能夠支持串行和并行模式，并且復(fù)用FastSimHelper 中提供的功能。負(fù)責(zé)MPI 通信的MPISvc 派生自SvcBase，能夠在配置跨節(jié)點(diǎn)運(yùn)行時(shí)按需被GPUFastSimTool 調(diào)用。負(fù)責(zé)GPU Context 的GPUSvc 類派生自SvcBase，這樣既可以被GPUFastSimTool調(diào)用，也可以被MPISvc調(diào)用。通過配置參數(shù)，軟件在運(yùn)行時(shí)將動態(tài)創(chuàng)建相應(yīng)的對象，從而支持單GPU、多GPU、MPI-GPU3 種加速模式。

圖4 GPU 繆子快速模擬的UML 類圖Fig.4 UML class diagram of GPU Muon fast simulation

多GPU模式的繆子快速模擬軟件復(fù)用了和單GPU一致的流程，如圖5 所示。在初始化階段，GPUSvc 按照設(shè)定的GPU 數(shù)量初始化GPU 設(shè)備列表，并檢測是否支持NVLink。該功能可在任意2 塊GPU 之間互傳數(shù)據(jù)，速度遠(yuǎn)高于PCIe 總線傳輸。在多GPU 模式下，初始化和結(jié)果合并階段使用GPU 間點(diǎn)對點(diǎn)傳輸?shù)姆绞剑蕴岣邤?shù)據(jù)傳輸效率。CPU 將直方圖數(shù)據(jù)傳遞給編號為0 的GPU，然后其他GPU 通過GPUDirect 拷貝編號為0 的GPU 中的數(shù)據(jù)。在事例循環(huán)階段，采取CPU 多線程來控制和調(diào)度多GPU：每個(gè)CPU 線程控制各自的GPU 設(shè)備，將Step 信息均勻分配至GPU，并異步啟動Kernel函數(shù)。CPU 主線程負(fù)責(zé)等待全部控制線程執(zhí)行結(jié)束，實(shí)現(xiàn)多GPU 設(shè)備同步。GPU Kernel 函數(shù)的結(jié)果由線程保存在各自開辟的定長數(shù)組中，以避免線程間資源競爭。對結(jié)果進(jìn)行合并的方式為：單GPU 模式的結(jié)果直接由GPU 傳回CPU；多GPU 模式的結(jié)果先匯總至編號0 的GPU 內(nèi)，再由編號0 的GPU 傳輸給CPU。GPU 計(jì)算完成后全局內(nèi)存不釋放，僅將結(jié)果置零，直到所有事例計(jì)算完成才釋放全局內(nèi)存。

圖5 GPU 繆子快速模擬時(shí)序圖Fig.5 Sequence diagram of GPU Muon fast simulation

3 MPI-GPU 并行計(jì)算

盡管單節(jié)點(diǎn)GPU 并行計(jì)算可以加速光子模擬，但同一節(jié)點(diǎn)可搭載的GPU 數(shù)量有限。計(jì)算集群GPU 資源豐富，若聚合集群中多個(gè)節(jié)點(diǎn)的GPU 算力，則光子模擬速度能進(jìn)一步提升。將算法移植到集群上并行運(yùn)行的關(guān)鍵是節(jié)點(diǎn)間信息的分發(fā)與收集，而MPI 能提供跨節(jié)點(diǎn)高性能數(shù)據(jù)通信。MPI-GPU 并行計(jì)算的核心思想是將步長集合進(jìn)行切割，通過MPI 消息將步長子集發(fā)送給GPU 節(jié)點(diǎn)進(jìn)行模擬計(jì)算，最后合并各個(gè)節(jié)點(diǎn)模擬的結(jié)果，從而完成跨節(jié)點(diǎn)并行計(jì)算。

如圖6 所示，MPI-GPU 快速模擬包含主節(jié)點(diǎn)和多個(gè)從節(jié)點(diǎn)，主節(jié)點(diǎn)負(fù)責(zé)跨節(jié)點(diǎn)部分的流程控制、步長子集的分配以及模擬結(jié)果的合并，從節(jié)點(diǎn)完成相應(yīng)步長子集的快速模擬。由于網(wǎng)絡(luò)存在延遲，上述節(jié)點(diǎn)都以多線程的方式啟動，并由專用線程負(fù)責(zé)節(jié)點(diǎn)之間的通信。在圖6 中，左右2 個(gè)框圖中的計(jì)算運(yùn)行于主節(jié)點(diǎn)。主節(jié)點(diǎn)收集步長信息并暫存于數(shù)組中，收集完畢后將其拆分形成任務(wù)，并放入隊(duì)列中。主節(jié)點(diǎn)中的專用線程會監(jiān)聽從節(jié)點(diǎn)的請求，并將任務(wù)分配給相應(yīng)的從節(jié)點(diǎn)。在任務(wù)調(diào)度時(shí)，會出現(xiàn)拆分后的任務(wù)數(shù)目多于子節(jié)點(diǎn)數(shù)目或者子節(jié)點(diǎn)運(yùn)行失敗的情況，主節(jié)點(diǎn)會將未執(zhí)行或執(zhí)行失敗的任務(wù)分發(fā)給空閑的從節(jié)點(diǎn)。完成相應(yīng)任務(wù)后的從節(jié)點(diǎn)會發(fā)送消息通知主節(jié)點(diǎn)，由主節(jié)點(diǎn)對從節(jié)點(diǎn)的運(yùn)行結(jié)果進(jìn)行緩存，并在所有任務(wù)完成后進(jìn)行合并。所有事例模擬完成后，主節(jié)點(diǎn)發(fā)送消息關(guān)閉從節(jié)點(diǎn)以及自身進(jìn)程。圖6 中間框圖為從節(jié)點(diǎn)，這些節(jié)點(diǎn)除了啟動GPU 并行計(jì)算線程，還會啟動通信線程用于與主節(jié)點(diǎn)的通信。當(dāng)從節(jié)點(diǎn)任務(wù)為空時(shí)，通信線程會與主節(jié)點(diǎn)通信，拉取模擬任務(wù)。獲得步長子集后，從節(jié)點(diǎn)根據(jù)配置執(zhí)行單GPU 或多GPU 的模擬，完成模擬后將模擬結(jié)果發(fā)送回主節(jié)點(diǎn)。接著從節(jié)點(diǎn)又進(jìn)入空閑狀態(tài)，等待主節(jié)點(diǎn)的任務(wù)分配。該流程不斷循環(huán)，直至所有任務(wù)完成。

圖6 MPI-GPU 快速模擬示意圖Fig.6 Schematic diagram of MPI-GPU fast simulation

MPI-GPU 繆子快速模擬的主從節(jié)點(diǎn)以獨(dú)立的MPI程序啟動，并且任務(wù)的分配與節(jié)點(diǎn)之間解耦合，從而實(shí)現(xiàn)通信獨(dú)立性、系統(tǒng)容錯性。實(shí)現(xiàn)的MPI 通信服務(wù)在主從模式的基礎(chǔ)上采用C/S 模式維護(hù)主從之間的通信連接。主節(jié)點(diǎn)作為服務(wù)器，阻塞等待接收各個(gè)作為客戶端的從節(jié)點(diǎn)的連接請求，任意客戶端斷開連接或者無響應(yīng)都不會使服務(wù)器停止運(yùn)行。主從模式與C/S 模式結(jié)合無法保證數(shù)據(jù)一致性，但具備可用性和容錯性，這與證明了數(shù)據(jù)一致性、可用性和分區(qū)容錯性為不可能三角的CAP 定理［24］相契合。通過哈希算法分配任務(wù)到節(jié)點(diǎn)時(shí)，需將任務(wù)與節(jié)點(diǎn)解耦，避免因從節(jié)點(diǎn)故障導(dǎo)致模擬失敗。在實(shí)現(xiàn)時(shí)采用了帶虛擬節(jié)點(diǎn)的一致性哈希算法［25］，在節(jié)點(diǎn)發(fā)送變化時(shí)能進(jìn)行任務(wù)的自動分配：主節(jié)點(diǎn)在任務(wù)執(zhí)行過程中，將所有未返回結(jié)果的任務(wù)暫時(shí)持久化；若在計(jì)算過程中突發(fā)從節(jié)點(diǎn)計(jì)算速度緩慢或者崩潰，該節(jié)點(diǎn)的任務(wù)將被重新分配到其他臨近節(jié)點(diǎn)，從而實(shí)現(xiàn)輕量級的錯誤恢復(fù)與負(fù)載均衡。

MPI通信服務(wù)的設(shè)計(jì)如圖7 所示，主要由接口類和3個(gè)組件構(gòu)成：對外接口MPISvc，管理MPI通信的Comm組件，負(fù)責(zé)節(jié)點(diǎn)內(nèi)部事務(wù)的Context組件以及輔助工具。圖7 中GPUFastSimTool調(diào)用MPISvc，完成步長數(shù)據(jù)在不同節(jié)點(diǎn)之間的傳輸以及GPU快速模擬的調(diào)用。Comm組件底層的MPI 通信由ClientSend、ServerSend 和MsgRecv 構(gòu)成。ClientSend 和ServerSend 類為C/S 模式的核心實(shí)現(xiàn)，負(fù)責(zé)主從節(jié)點(diǎn)MPI 消息發(fā)送；MsgRecv 用于接收處理MPI 消息。Comm 組件中CommManager類是控制節(jié)點(diǎn)消息分發(fā)的接口，主從節(jié)點(diǎn)分別派生實(shí)現(xiàn)相應(yīng)的功能。Context 組件以BaseContex 為基類，以派生類MasterContext 和SlaveContext 類為核心，負(fù)責(zé)處理SNiPER 框架異步通信事務(wù)與錯誤恢復(fù)功能。Util組件提供線程池（ThreadPool 類）、進(jìn)程內(nèi)部的線程間通信（ThreadComm類）以及哈希值計(jì)算（MFSHash類）。MFSHash 類對序列化后的步長數(shù)據(jù)計(jì)算哈希值從而得到目標(biāo)GPU 節(jié)點(diǎn)。

圖7 MPI 通信服務(wù)UML 類圖Fig.7 UML class diagram of MPI communication service

4 性能測試與分析

針對本文完成的GPU 繆子快速模擬進(jìn)行如下的性能測試：首先對模擬的結(jié)果進(jìn)行驗(yàn)證，確保方法的正確性；然后分別對單GPU、多GPU、MPI-GPU 繆子快速模擬的3 種模式進(jìn)行性能測量。對于單GPU，測量了總光子數(shù)與加速比的關(guān)系，優(yōu)化得出最佳的光子數(shù)模擬區(qū)間。對于多GPU，測量了GPU 數(shù)目與加速比之間的關(guān)系，優(yōu)化多GPU 場景下程序的可擴(kuò)展性和線性加速比。對于MPI-GPU，更多的開銷將來自數(shù)據(jù)在網(wǎng)絡(luò)中的傳輸，通過測量跨節(jié)點(diǎn)情況下GPU 節(jié)點(diǎn)的利用率能夠優(yōu)化跨節(jié)點(diǎn)快速模擬的參數(shù)。這些參數(shù)的優(yōu)化為大規(guī)模數(shù)據(jù)產(chǎn)生提供重要的參考依據(jù)。

繆子快速模擬測試是在中國科學(xué)院高能物理研究所計(jì)算中心的GPU 計(jì)算集群上完成的。GPU 計(jì)算集群中的每個(gè)節(jié)點(diǎn)都配備了2 顆主頻為2.30 GHz的Intel Xeon Gold 6140 CPU 處理器，每個(gè)CPU 處理器包含18 個(gè)物理內(nèi)核。每個(gè)節(jié)點(diǎn)還搭載了8 塊NVIDIA Volta 構(gòu)架Tesla V100 SXM2 GPU 計(jì)算卡，同一個(gè)節(jié)點(diǎn)上的GPU 卡之間通過NVLink 互聯(lián)，單個(gè)GPU 擁有32 GB 內(nèi)存。InfiniBand 網(wǎng)絡(luò)交換機(jī)Mellanox MCS7520/MSB7800 實(shí)現(xiàn)集群計(jì)算節(jié)點(diǎn)間的互聯(lián)，節(jié)點(diǎn)間的數(shù)據(jù)傳輸速率為100 Gb/s，延時(shí)小于1 μs。整個(gè)GPU 集群由Slurm［26］進(jìn)行管理和作業(yè)的調(diào)度，集群節(jié)點(diǎn)軟硬件環(huán)境如表1 所示。由于該集群被多個(gè)用戶共享使用，在測試過程中申請的計(jì)算資源都以整個(gè)節(jié)點(diǎn)為單位，避免其他用戶的計(jì)算作業(yè)對測試結(jié)果的干擾。通過在專用節(jié)點(diǎn)進(jìn)行測試，可以實(shí)現(xiàn)測試環(huán)境和生產(chǎn)環(huán)境的隔離，避免對測試環(huán)境的干擾。

表1 測試節(jié)點(diǎn)軟硬件環(huán)境Table 1 Software and hardware environment of test node

CPU 繆子快速模擬和GPU 繆子快速模擬使用了相同的事例進(jìn)行測試，以確保運(yùn)行時(shí)加載相同的Step信息，驗(yàn)證GPU 繆子快速模擬算法結(jié)果的正確性。在對結(jié)果進(jìn)行校驗(yàn)時(shí)，選擇PMT 擊中時(shí)間以及每個(gè)PMT上nPE 的分布進(jìn)行比較。圖8 為GPU 繆子快速模擬和CPU 繆子快速模擬的模擬結(jié)果比較，縱軸使用對數(shù)坐標(biāo)展示事例數(shù)。圖8（a）為繆子事例中光子擊中PMT的時(shí)間分布，直方圖為GPU 繆子快速模擬結(jié)果，圓點(diǎn)為CPU 繆子快速模擬的結(jié)果，結(jié)果顯示擊中時(shí)間的分布一致。圖8（b）為繆子事例中擊中PMT 的光子數(shù)量分布，同樣地，直方圖為GPU 繆子快速模擬結(jié)果，圓點(diǎn)為CPU 繆子快速模擬的結(jié)果。

圖8 GPU 和CPU 繆子快速模擬的結(jié)果驗(yàn)證Fig.8 Result verification of GPU and CPU Muon fast simulation

從圖8 可以看出，在光子數(shù)量（nPE）小于150 時(shí)分布一致，而nPE 大于150 時(shí)也在統(tǒng)計(jì)誤差范圍內(nèi)。圖8（a）和圖8（b）中GPU 繆子快速模擬和CPU 繆子快速模擬的結(jié)果有細(xì)微的差別，這是由于兩者在直方圖抽樣過程中使用了不同的隨機(jī)數(shù)引擎造成的，CPU 繆子快速模擬使用線性同余隨機(jī)數(shù)（LCG）生成器，GPU 繆子快速模擬使用線性反饋移位寄存器的變種XORWOW 生成器。nPE>150 時(shí)兩者結(jié)果差距變大，是因?yàn)閳D8 縱坐標(biāo)取對數(shù)后，微小的差異在數(shù)據(jù)量較大時(shí)不明顯，nPE>150 時(shí)數(shù)據(jù)量小，抽樣次數(shù)少，結(jié)果波動較大，差異更明顯，但從整體上看，兩者結(jié)果基本一致。

在測試單GPU 加速模式中GPU 算法對CPU 算法的加速性能時(shí)，2 個(gè)算法利用相同的參數(shù)生成了數(shù)量一致的模擬結(jié)果，CPU 算法為單線程串行執(zhí)行，GPU 算法為單GPU 內(nèi)部并行。測量過程中的運(yùn)行時(shí)間僅包含快速模擬核心算法部分，以減少初始化加載參數(shù)直方圖和數(shù)據(jù)傳輸?shù)葘y量結(jié)果的影響。每個(gè)測量點(diǎn)都重復(fù)100 次后取平均值和標(biāo)準(zhǔn)差。在測試中需要改變產(chǎn)生的光子總數(shù)，測試過程中以215GeV 繆子事例為基準(zhǔn)，對Step 進(jìn)行采樣和重復(fù)2 種策略。例如，產(chǎn)生10%的總光子數(shù)時(shí)只需要按照10%的概率對每個(gè)Step 進(jìn)行篩選，產(chǎn)生2 倍的總光子數(shù)時(shí)，則對相同的Step 調(diào)用2 次。圖9 給出單GPU加速模式下加速比與總光子數(shù)的關(guān)系。圖中直方圖為典型的215GeV 繆子事例光子分布，加速比的定義為CPU 快速模擬的運(yùn)行時(shí)間除以GPU 快速模擬的運(yùn)行時(shí)間。

圖9 單GPU 加速模式加速比與光子數(shù)的關(guān)系Fig.9 Relationship of single-GPU mode acceleration ratio and number of photons

從圖9 可知，在光子數(shù)較少的情況下，GPU 加速效果不佳，未能充分利用GPU 硬件資源；隨著光子數(shù)增大，GPU 處理的數(shù)據(jù)量增大，GPU 算法加速比穩(wěn)步提升，波動趨緩，最終穩(wěn)定在33 倍附近，達(dá)到單GPU 性能極限。結(jié)合CPU 繆子快速模擬對全模擬的加速約為6 倍，單GPU 繆子快速模擬對全模擬加速約為198 倍。圖中的直方圖顯示了能量為215GeV 繆子事例的總光子數(shù)分布情況，可以看出大部分繆子事例產(chǎn)生的光子數(shù)處于GPU 加速比穩(wěn)定的區(qū)域。利用單GPU 加速繆子快速模擬可取得良好效果，模擬生成1 000 個(gè)事例數(shù)據(jù)從3 天縮短至2 h，使短時(shí)間內(nèi)生成大量數(shù)據(jù)成為可能。

在測量多GPU 加速模式提升GPU 算法的性能方面，本文比較在不同GPU 數(shù)量的條件下GPU 算法與CPU 算法的時(shí)間，并以加速比的形式展現(xiàn)。在多GPU 加速模式下，申請的配置為同一節(jié)點(diǎn)的單核CPU 與多塊GPU。圖10 顯示了在不同GPU 數(shù)量下，以同樣的隨機(jī)數(shù)種子生成相同的100 個(gè)事例的GPU加速比。圖中虛線由單GPU 加速比的線形倍增得到，實(shí)線是多GPU 加速比的測量值。從圖10 中可以看出，GPU 數(shù)量越多，加速比越高，趨近于線形增長。當(dāng)GPU 數(shù)量大于3 時(shí)，加速比增長率開始偏離理想情況，說明多GPU 并行時(shí)系統(tǒng)運(yùn)行效率下降，最大性能損耗為9.13%，仍保持較高效率。存在性能損耗的原因是：一方面，更多GPU 導(dǎo)致數(shù)據(jù)傳輸時(shí)間有所增加；另一方面，GPU 切換開銷極大，結(jié)果合并時(shí)也將多次啟動核函數(shù)計(jì)算。圖中8 塊GPU 并行實(shí)測加速比最高為253.7倍，是單卡性能的7.6倍，加速效果顯著。

圖10 多GPU 模式加速比Fig.10 Acceleration ratio of multi-GPU mode

GPU 繆子快速模擬利用分布式計(jì)算可進(jìn)一步擴(kuò)大并行計(jì)算規(guī)模，但當(dāng)需要并行計(jì)算的數(shù)據(jù)量不足時(shí)，分布式計(jì)算帶來的速度提升不足以抵消分布式通信引入的額外開銷。通過不斷提高光子數(shù)并測量節(jié)點(diǎn)的GPU 計(jì)算用時(shí)和系統(tǒng)總用時(shí)，來測試GPU 繆子快速模擬的數(shù)據(jù)量能否使MPI 分布式框架真實(shí)有效地利用GPU 節(jié)點(diǎn)，兩者的用時(shí)之比即為GPU 節(jié)點(diǎn)的利用率。

MPI-GPU 繆子快速模擬在集群中以集群作業(yè)的方式提交運(yùn)行，在通常情況下，主節(jié)點(diǎn)不必等待所有從節(jié)點(diǎn)開始執(zhí)行。在測試時(shí)主節(jié)點(diǎn)需等待從節(jié)點(diǎn)全部成功運(yùn)行方可解除阻塞等待，作業(yè)排隊(duì)的時(shí)間不包含在測試運(yùn)行時(shí)間內(nèi)，使結(jié)果更為準(zhǔn)確。相比單節(jié)點(diǎn)的GPU 并行計(jì)算，MPI-GPU 繆子快速模擬在分布式通信時(shí)有額外的時(shí)間開銷：MPI 初始化的時(shí)間，計(jì)算時(shí)Step 信息分發(fā)的時(shí)間，結(jié)果合并的時(shí)間。其中MPI 初始化時(shí)間在集群作業(yè)無需排隊(duì)的情況下可忽略不計(jì)。

在測試不同節(jié)點(diǎn)數(shù)量對MPI-GPU 版本的影響時(shí)，使用2～4 個(gè)節(jié)點(diǎn)進(jìn)行測試，每個(gè)節(jié)點(diǎn)僅使用單核單GPU，其中包含1 個(gè)主節(jié)點(diǎn)，主節(jié)點(diǎn)GPU 參與計(jì)算，并多次測試發(fā)送總量相同的光子數(shù)。圖11 所示為MPI-GPU 繆子快速模擬在跨節(jié)點(diǎn)分布式通信時(shí)，每次發(fā)送Step 信息對應(yīng)的光子數(shù)與運(yùn)行時(shí)GPU 節(jié)點(diǎn)利用率的關(guān)系。GPU 節(jié)點(diǎn)利用率是GPU 計(jì)算時(shí)間與從節(jié)點(diǎn)光子模擬總運(yùn)行時(shí)間的比值。

圖11 分布式計(jì)算GPU 節(jié)點(diǎn)利用率Fig.11 GPU node utilization rate of distributed computing

從圖11 可以看出，在發(fā)送總光子數(shù)固定的前提下，單次發(fā)送光子數(shù)量較少時(shí)，GPU 節(jié)點(diǎn)利用率低，因?yàn)閱未伟l(fā)送光子數(shù)量越少，MPI 通信次數(shù)越多，通信開銷增長而計(jì)算開銷不變。隨著發(fā)送光子數(shù)增長至500×106后，GPU 節(jié)點(diǎn)利用率持續(xù)增加直至平穩(wěn)，GPU 繆子快速模擬生成的光子數(shù)越多，分布式計(jì)算越能充分發(fā)揮多節(jié)點(diǎn)性能優(yōu)勢，且應(yīng)盡可能減少通信次數(shù)，提高單次通信數(shù)據(jù)量，單次發(fā)送Step 信息對應(yīng)的光子數(shù)應(yīng)保持在500×106，以提高節(jié)點(diǎn)利用率。隨著節(jié)點(diǎn)數(shù)量增加，GPU 利用率從90%下降至85%，因?yàn)楣?jié)點(diǎn)數(shù)增多主節(jié)點(diǎn)通信開銷增大，利用率稍有下降。當(dāng)系統(tǒng)節(jié)點(diǎn)數(shù)量小于5 個(gè)時(shí)，系統(tǒng)具有較好的橫向擴(kuò)展性，若配合縱向擴(kuò)展，啟用節(jié)點(diǎn)內(nèi)部多GPU 并行計(jì)算，可將加速比進(jìn)一步倍增。

在國際上，盡管不同高能物理實(shí)驗(yàn)的探測器復(fù)雜度和物理過程不盡相同，但快速模擬方法相比于GEANT4 全模擬的加速比一般在數(shù)十倍到數(shù)百倍之間。本文提出的GPU 繆子快速模擬方法相比于GEANT4 全模擬最終可以達(dá)到1 500 倍，加速比相對于傳統(tǒng)的快速模擬方法提高了一個(gè)量級。

5 結(jié)束語

繆子快速模擬對于中微子實(shí)驗(yàn)至關(guān)重要，但生成大量數(shù)據(jù)時(shí)存在性能瓶頸。為此，本文提出一種基于GPU 的分布式繆子快速模擬方法。通過改進(jìn)快速模擬算法，實(shí)現(xiàn)在GPU 并行閃爍光的產(chǎn)生與模擬，減少CPU 內(nèi)存開銷，并通過單GPU、多GPU、MPI-GPU 3 種模式對快速模擬進(jìn)行加速，利用GPU的高數(shù)據(jù)吞吐量，在單GPU 加速模式中加速比高達(dá)33 倍，多GPU 加速模式最高可達(dá)約250 倍；采用MPI分布式通信實(shí)現(xiàn)跨節(jié)點(diǎn)并行，4 個(gè)節(jié)點(diǎn)以內(nèi)GPU 節(jié)點(diǎn)利用率高于85%，橫向擴(kuò)展性良好，此外MPI-GPU框架的任務(wù)分配、容錯處理保證了快速模擬在集群中穩(wěn)定運(yùn)行。測試結(jié)果驗(yàn)證了GPU 繆子快速模擬系統(tǒng)結(jié)果的正確性。但多GPU 并行、MPI 多節(jié)點(diǎn)并行依然存在一定損耗，GPU 加速比還存在一定提升空間，下一步將提升MPI 通信和多GPU 并行效率，提供更通用的并行API 接口。