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

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

0 概述

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

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

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

1 繆子快速模擬算法

JUNO 模擬軟件由物理產(chǎn)生子、探測(cè)器模擬、電子學(xué)模擬3 個(gè)部分構(gòu)成，如圖1 所示。物理產(chǎn)生子用于產(chǎn)生原初粒子信息，之后再由探測(cè)器模擬對(duì)原初粒子進(jìn)行運(yùn)輸和跟蹤。模擬軟件基于Geant4 和SNiPER［22］開發(fā)，用于模擬粒子和探測(cè)器之間的相互作用，并記錄探測(cè)器中的響應(yīng)信息，最后由電子學(xué)模擬對(duì)擊中信息數(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

快速模擬方法將繆子徑跡分割成很多份短的步長(zhǎng)（Step），如圖2 所示。對(duì)于每個(gè)Step，需要根據(jù)位置和沉積能量產(chǎn)生探測(cè)器的響應(yīng)。Step 所在位置可以確定徑向位置R以及該位置與每個(gè)PMT 之間的夾角θ，可以唯一確定需要加載PMT 的光電子（nPE）直方圖和擊中時(shí)間（Hit time）直方圖數(shù)據(jù)?？焖倌M的第1 步是根據(jù)每個(gè)Step 的可見(jiàn)能量對(duì)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 繆子快速模擬

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

GPU 繆子快速模擬方法在初始化階段將預(yù)生成的直方圖數(shù)據(jù)經(jīng)過(guò)序列化存儲(chǔ)于GPU 全局內(nèi)存中，供后續(xù)內(nèi)核算法訪問(wèn)。在事例循環(huán)階段，該方法首先從Geant4 收集步長(zhǎng)信息，并進(jìn)行實(shí)時(shí)序列化，用于后續(xù)傳輸至GPU 內(nèi)存。該方法使用多GPU 及MPI進(jìn)行并行模擬，對(duì)步長(zhǎng)信息自動(dòng)劃分，并將劃分后的數(shù)據(jù)傳輸至對(duì)應(yīng)的GPU 全局內(nèi)存，解決了因單塊GPU 顯存的限制無(wú)法處理所有步長(zhǎng)信息的問(wèn)題，從而實(shí)現(xiàn)大規(guī)模并行。該方法實(shí)現(xiàn)的CUDA Kernel 函數(shù)將一個(gè)PMT 的模擬映射到一個(gè)CUDA 線程，通過(guò)啟動(dòng)大量的線程完成大規(guī)模的模擬計(jì)算，并且將中間結(jié)果存儲(chǔ)于GPU 全局內(nèi)存。在每個(gè)事例終止時(shí)，合并所有GPU 的數(shù)據(jù)，并將結(jié)果傳回CPU，從而完成一個(gè)事例的模擬。GPU 繆子快速模擬的Kernel 函數(shù)主要負(fù)責(zé)光子的模擬，用于替換串行的CPU 算法。CPU 算法中對(duì)于每個(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)用。對(duì)于每個(gè)物理事例，收集Geant4產(chǎn)生的步長(zhǎng)信息，然后拷貝至GPU，對(duì)應(yīng)于第3 行、第4行。隨后，調(diào)用CUDA Kernel函數(shù)，通過(guò)block和thread 2 個(gè)參數(shù)設(shè)置并發(fā)線程數(shù)，且總的并發(fā)數(shù)不少于PMT數(shù)目，對(duì)應(yīng)于第5 行。模擬完成后，再將結(jié)果合并拷貝回CPU，對(duì)應(yīng)于第6 行、第7 行。

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

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

對(duì)于偽代碼中SimulationHits 算法所描述的直方圖抽樣的過(guò)程，即代碼第26 行和第32 行中的getRandom，實(shí)現(xiàn)過(guò)程如下：首先在初始化加載直方圖數(shù)據(jù)時(shí)，將概率分布直方圖轉(zhuǎn)化為累計(jì)分布函數(shù)直方圖，直方圖在全局內(nèi)存中以定長(zhǎng)數(shù)組的形式保存。GPU 線程通過(guò)Step 信息中的R和θ找到對(duì)應(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 級(jí)索引的方式，保證Warp 中的線程能夠訪問(wèn)連續(xù)內(nèi)存空間的直方圖：第1 級(jí)索引為距離R，長(zhǎng)度為200；第2 級(jí)索引為角度θ，長(zhǎng)度為180。在圖3中，32 個(gè)線程訪問(wèn)的直方圖具有相同的距離R和相近的角度θ，能夠合并訪問(wèn)GPU 全局內(nèi)存空間，從而提高數(shù)據(jù)訪問(wèn)效率。

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

除訪問(wèn)直方圖數(shù)據(jù)外，GPU 線程還要訪問(wèn)步長(zhǎng)信息。預(yù)先產(chǎn)生的直方圖數(shù)據(jù)數(shù)目固定，而步長(zhǎng)的數(shù)目則與模擬事例的能量相關(guān)，能量越高產(chǎn)生的步長(zhǎng)數(shù)目越多，所以對(duì)步長(zhǎng)信息的實(shí)現(xiàn)與直方圖數(shù)據(jù)不同。針對(duì)CUDA 多次啟動(dòng)核函數(shù)開銷高的問(wèn)題，GPU 算法并不對(duì)每個(gè)Step 單獨(dú)調(diào)用快速模擬的方法，而是創(chuàng)建步長(zhǎng)信息的集合，將Step 信息集中收集，待Step 收集完成后，再一次性或分批調(diào)用并行化快速模擬算法，從而提高算法效率。通過(guò)引入步長(zhǎng)集合的設(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)用。通過(guò)配置參數(shù)，軟件在運(yùn)行時(shí)將動(dòng)態(tài)創(chuàng)建相應(yīng)的對(duì)象，從而支持單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è)備列表，并檢測(cè)是否支持NVLink。該功能可在任意2 塊GPU 之間互傳數(shù)據(jù)，速度遠(yuǎn)高于PCIe 總線傳輸。在多GPU 模式下，初始化和結(jié)果合并階段使用GPU 間點(diǎn)對(duì)點(diǎn)傳輸?shù)姆绞?，以提高?shù)據(jù)傳輸效率。CPU 將直方圖數(shù)據(jù)傳遞給編號(hào)為0 的GPU，然后其他GPU 通過(guò)GPUDirect 拷貝編號(hào)為0 的GPU 中的數(shù)據(jù)。在事例循環(huán)階段，采取CPU 多線程來(lái)控制和調(diào)度多GPU：每個(gè)CPU 線程控制各自的GPU 設(shè)備，將Step 信息均勻分配至GPU，并異步啟動(dòng)Kernel函數(shù)。CPU 主線程負(fù)責(zé)等待全部控制線程執(zhí)行結(jié)束，實(shí)現(xiàn)多GPU 設(shè)備同步。GPU Kernel 函數(shù)的結(jié)果由線程保存在各自開辟的定長(zhǎng)數(shù)組中，以避免線程間資源競(jìng)爭(zhēng)。對(duì)結(jié)果進(jìn)行合并的方式為：?jiǎn)蜧PU 模式的結(jié)果直接由GPU 傳回CPU；多GPU 模式的結(jié)果先匯總至編號(hào)0 的GPU 內(nèi)，再由編號(hào)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ì)算的核心思想是將步長(zhǎng)集合進(jìn)行切割，通過(guò)MPI 消息將步長(zhǎng)子集發(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)部分的流程控制、步長(zhǎng)子集的分配以及模擬結(jié)果的合并，從節(jié)點(diǎn)完成相應(yīng)步長(zhǎng)子集的快速模擬。由于網(wǎng)絡(luò)存在延遲，上述節(jié)點(diǎn)都以多線程的方式啟動(dòng)，并由專用線程負(fù)責(zé)節(jié)點(diǎn)之間的通信。在圖6 中，左右2 個(gè)框圖中的計(jì)算運(yùn)行于主節(jié)點(diǎn)。主節(jié)點(diǎn)收集步長(zhǎng)信息并暫存于數(shù)組中，收集完畢后將其拆分形成任務(wù)，并放入隊(duì)列中。主節(jié)點(diǎn)中的專用線程會(huì)監(jiān)聽(tīng)從節(jié)點(diǎn)的請(qǐng)求，并將任務(wù)分配給相應(yīng)的從節(jié)點(diǎn)。在任務(wù)調(diào)度時(shí)，會(huì)出現(xiàn)拆分后的任務(wù)數(shù)目多于子節(jié)點(diǎn)數(shù)目或者子節(jié)點(diǎn)運(yùn)行失敗的情況，主節(jié)點(diǎn)會(huì)將未執(zhí)行或執(zhí)行失敗的任務(wù)分發(fā)給空閑的從節(jié)點(diǎn)。完成相應(yīng)任務(wù)后的從節(jié)點(diǎn)會(huì)發(fā)送消息通知主節(jié)點(diǎn)，由主節(jié)點(diǎn)對(duì)從節(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)除了啟動(dòng)GPU 并行計(jì)算線程，還會(huì)啟動(dòng)通信線程用于與主節(jié)點(diǎn)的通信。當(dāng)從節(jié)點(diǎn)任務(wù)為空時(shí)，通信線程會(huì)與主節(jié)點(diǎn)通信，拉取模擬任務(wù)。獲得步長(zhǎng)子集后，從節(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程序啟動(dòng)，并且任務(wù)的分配與節(jié)點(diǎn)之間解耦合，從而實(shí)現(xiàn)通信獨(dú)立性、系統(tǒng)容錯(cuò)性。實(shí)現(xiàn)的MPI 通信服務(wù)在主從模式的基礎(chǔ)上采用C/S 模式維護(hù)主從之間的通信連接。主節(jié)點(diǎn)作為服務(wù)器，阻塞等待接收各個(gè)作為客戶端的從節(jié)點(diǎn)的連接請(qǐng)求，任意客戶端斷開連接或者無(wú)響應(yīng)都不會(huì)使服務(wù)器停止運(yùn)行。主從模式與C/S 模式結(jié)合無(wú)法保證數(shù)據(jù)一致性，但具備可用性和容錯(cuò)性，這與證明了數(shù)據(jù)一致性、可用性和分區(qū)容錯(cuò)性為不可能三角的CAP 定理［24］相契合。通過(guò)哈希算法分配任務(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ù)的自動(dòng)分配：主節(jié)點(diǎn)在任務(wù)執(zhí)行過(guò)程中，將所有未返回結(jié)果的任務(wù)暫時(shí)持久化；若在計(jì)算過(guò)程中突發(fā)從節(jié)點(diǎn)計(jì)算速度緩慢或者崩潰，該節(jié)點(diǎn)的任務(wù)將被重新分配到其他臨近節(jié)點(diǎn)，從而實(shí)現(xiàn)輕量級(jí)的錯(cuò)誤恢復(fù)與負(fù)載均衡。

MPI通信服務(wù)的設(shè)計(jì)如圖7 所示，主要由接口類和3個(gè)組件構(gòu)成：對(duì)外接口MPISvc，管理MPI通信的Comm組件，負(fù)責(zé)節(jié)點(diǎn)內(nèi)部事務(wù)的Context組件以及輔助工具。圖7 中GPUFastSimTool調(diào)用MPISvc，完成步長(zhǎng)數(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ù)與錯(cuò)誤恢復(fù)功能。Util組件提供線程池（ThreadPool 類）、進(jìn)程內(nèi)部的線程間通信（ThreadComm類）以及哈希值計(jì)算（MFSHash類）。MFSHash 類對(duì)序列化后的步長(zhǎng)數(shù)據(jù)計(jì)算哈希值從而得到目標(biāo)GPU 節(jié)點(diǎn)。

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

4 性能測(cè)試與分析

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

繆子快速模擬測(cè)試是在中國(guó)科學(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 卡之間通過(guò)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è)用戶共享使用，在測(cè)試過(guò)程中申請(qǐng)的計(jì)算資源都以整個(gè)節(jié)點(diǎn)為單位，避免其他用戶的計(jì)算作業(yè)對(duì)測(cè)試結(jié)果的干擾。通過(guò)在專用節(jié)點(diǎn)進(jìn)行測(cè)試，可以實(shí)現(xiàn)測(cè)試環(huán)境和生產(chǎn)環(huán)境的隔離，避免對(duì)測(cè)試環(huán)境的干擾。

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

CPU 繆子快速模擬和GPU 繆子快速模擬使用了相同的事例進(jìn)行測(cè)試，以確保運(yùn)行時(shí)加載相同的Step信息，驗(yàn)證GPU 繆子快速模擬算法結(jié)果的正確性。在對(duì)結(jié)果進(jìn)行校驗(yàn)時(shí)，選擇PMT 擊中時(shí)間以及每個(gè)PMT上nPE 的分布進(jìn)行比較。圖8 為GPU 繆子快速模擬和CPU 繆子快速模擬的模擬結(jié)果比較，縱軸使用對(duì)數(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ì)微的差別，這是由于兩者在直方圖抽樣過(guò)程中使用了不同的隨機(jī)數(shù)引擎造成的，CPU 繆子快速模擬使用線性同余隨機(jī)數(shù)（LCG）生成器，GPU 繆子快速模擬使用線性反饋移位寄存器的變種XORWOW 生成器。nPE>150 時(shí)兩者結(jié)果差距變大，是因?yàn)閳D8 縱坐標(biāo)取對(duì)數(shù)后，微小的差異在數(shù)據(jù)量較大時(shí)不明顯，nPE>150 時(shí)數(shù)據(jù)量小，抽樣次數(shù)少，結(jié)果波動(dòng)較大，差異更明顯，但從整體上看，兩者結(jié)果基本一致。

在測(cè)試單GPU 加速模式中GPU 算法對(duì)CPU 算法的加速性能時(shí)，2 個(gè)算法利用相同的參數(shù)生成了數(shù)量一致的模擬結(jié)果，CPU 算法為單線程串行執(zhí)行，GPU 算法為單GPU 內(nèi)部并行。測(cè)量過(guò)程中的運(yùn)行時(shí)間僅包含快速模擬核心算法部分，以減少初始化加載參數(shù)直方圖和數(shù)據(jù)傳輸?shù)葘?duì)測(cè)量結(jié)果的影響。每個(gè)測(cè)量點(diǎn)都重復(fù)100 次后取平均值和標(biāo)準(zhǔn)差。在測(cè)試中需要改變產(chǎn)生的光子總數(shù)，測(cè)試過(guò)程中以215GeV 繆子事例為基準(zhǔn)，對(duì)Step 進(jìn)行采樣和重復(fù)2 種策略。例如，產(chǎn)生10%的總光子數(shù)時(shí)只需要按照10%的概率對(duì)每個(gè)Step 進(jìn)行篩選，產(chǎn)生2 倍的總光子數(shù)時(shí)，則對(duì)相同的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)步提升，波動(dòng)趨緩，最終穩(wěn)定在33 倍附近，達(dá)到單GPU 性能極限。結(jié)合CPU 繆子快速模擬對(duì)全模擬的加速約為6 倍，單GPU 繆子快速模擬對(duì)全模擬加速約為198 倍。圖中的直方圖顯示了能量為215GeV 繆子事例的總光子數(shù)分布情況，可以看出大部分繆子事例產(chǎn)生的光子數(shù)處于GPU 加速比穩(wěn)定的區(qū)域。利用單GPU 加速繆子快速模擬可取得良好效果，模擬生成1 000 個(gè)事例數(shù)據(jù)從3 天縮短至2 h，使短時(shí)間內(nèi)生成大量數(shù)據(jù)成為可能。

在測(cè)量多GPU 加速模式提升GPU 算法的性能方面，本文比較在不同GPU 數(shù)量的條件下GPU 算法與CPU 算法的時(shí)間，并以加速比的形式展現(xiàn)。在多GPU 加速模式下，申請(qǐng)的配置為同一節(jié)點(diǎn)的單核CPU 與多塊GPU。圖10 顯示了在不同GPU 數(shù)量下，以同樣的隨機(jī)數(shù)種子生成相同的100 個(gè)事例的GPU加速比。圖中虛線由單GPU 加速比的線形倍增得到，實(shí)線是多GPU 加速比的測(cè)量值。從圖10 中可以看出，GPU 數(shù)量越多，加速比越高，趨近于線形增長(zhǎng)。當(dāng)GPU 數(shù)量大于3 時(shí)，加速比增長(zhǎng)率開始偏離理想情況，說(shuō)明多GPU 并行時(shí)系統(tǒng)運(yùn)行效率下降，最大性能損耗為9.13%，仍保持較高效率。存在性能損耗的原因是：一方面，更多GPU 導(dǎo)致數(shù)據(jù)傳輸時(shí)間有所增加；另一方面，GPU 切換開銷極大，結(jié)果合并時(shí)也將多次啟動(dòng)核函數(shù)計(jì)算。圖中8 塊GPU 并行實(shí)測(cè)加速比最高為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ì)算帶來(lái)的速度提升不足以抵消分布式通信引入的額外開銷。通過(guò)不斷提高光子數(shù)并測(cè)量節(jié)點(diǎn)的GPU 計(jì)算用時(shí)和系統(tǒng)總用時(shí)，來(lái)測(cè)試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í)行。在測(cè)試時(shí)主節(jié)點(diǎn)需等待從節(jié)點(diǎn)全部成功運(yùn)行方可解除阻塞等待，作業(yè)排隊(duì)的時(shí)間不包含在測(cè)試運(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è)無(wú)需排隊(duì)的情況下可忽略不計(jì)。

在測(cè)試不同節(jié)點(diǎn)數(shù)量對(duì)MPI-GPU 版本的影響時(shí)，使用2～4 個(gè)節(jié)點(diǎn)進(jìn)行測(cè)試，每個(gè)節(jié)點(diǎn)僅使用單核單GPU，其中包含1 個(gè)主節(jié)點(diǎn)，主節(jié)點(diǎn)GPU 參與計(jì)算，并多次測(cè)試發(fā)送總量相同的光子數(shù)。圖11 所示為MPI-GPU 繆子快速模擬在跨節(jié)點(diǎn)分布式通信時(shí)，每次發(fā)送Step 信息對(duì)應(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ù)越多，通信開銷增長(zhǎng)而計(jì)算開銷不變。隨著發(fā)送光子數(shù)增長(zhǎng)至500×106后，GPU 節(jié)點(diǎn)利用率持續(xù)增加直至平穩(wěn)，GPU 繆子快速模擬生成的光子數(shù)越多，分布式計(jì)算越能充分發(fā)揮多節(jié)點(diǎn)性能優(yōu)勢(shì)，且應(yīng)盡可能減少通信次數(shù)，提高單次通信數(shù)據(jù)量，單次發(fā)送Step 信息對(duì)應(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)一步倍增。

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

5 結(jié)束語(yǔ)

繆子快速模擬對(duì)于中微子實(shí)驗(yàn)至關(guān)重要，但生成大量數(shù)據(jù)時(shí)存在性能瓶頸。為此，本文提出一種基于GPU 的分布式繆子快速模擬方法。通過(guò)改進(jìn)快速模擬算法，實(shí)現(xiàn)在GPU 并行閃爍光的產(chǎn)生與模擬，減少CPU 內(nèi)存開銷，并通過(guò)單GPU、多GPU、MPI-GPU 3 種模式對(duì)快速模擬進(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ù)分配、容錯(cuò)處理保證了快速模擬在集群中穩(wěn)定運(yùn)行。測(cè)試結(jié)果驗(yàn)證了GPU 繆子快速模擬系統(tǒng)結(jié)果的正確性。但多GPU 并行、MPI 多節(jié)點(diǎn)并行依然存在一定損耗，GPU 加速比還存在一定提升空間，下一步將提升MPI 通信和多GPU 并行效率，提供更通用的并行API 接口。