劉 旭,張曦煌,劉 釗,呂小敬,朱光輝

(1.江南大學(xué) 物聯(lián)網(wǎng)工程學(xué)院,江蘇 無錫 214122; 2.國家超級計算無錫中心,江蘇 無錫 214072; 3.清華大學(xué),北京 100084; 4.無錫航天江南數(shù)據(jù)系統(tǒng)科技有限公司,江蘇 無錫 214000)

0 概述

近年來,科學(xué)家通過超級計算機模擬研究揭示了早期宇宙中星系與暗物質(zhì)的形成過程。據(jù)估計,暗物質(zhì)約占宇宙物質(zhì)的85%,占總能量密度的1/4。由于暗物質(zhì)對整個電磁光譜而言都是不可見的,因此需要在超級計算機上對暗物質(zhì)演化模型的建立過程進行數(shù)值模擬。宇宙學(xué)模擬對于科學(xué)家研究非線性結(jié)構(gòu)的形成以及暗物質(zhì)、暗能量等假想形式具有重要作用,并且宇宙N體模擬已成為超級計算機的主要應(yīng)用。

超級計算機作為大規(guī)模高分辨率模擬和應(yīng)用的重要工具,在過去幾十年發(fā)展迅速。在20年內(nèi),頂級超級計算機的峰值性能已經(jīng)從每秒1012次浮點運算(TFLOPS)增長到每秒1015次浮點運算(PFLOPS),現(xiàn)在正朝著每秒1018次浮點運算(EFLOPS)邁進。這種快速的增長趨勢很大程度上是源于使用GPU或多核芯片等加速器設(shè)備的異構(gòu)架構(gòu)興起。我國的神威太湖之光超級計算機由40 960塊自主研發(fā)的多核芯片SW26010組成,整個系統(tǒng)具有1 000多萬個內(nèi)核,峰值性能為125 PFLOPS。神威太湖之光強大的計算能力使其成為解決超大規(guī)模問題的理想平臺。因此,自2016年神威太湖之光發(fā)布以來,已經(jīng)在這臺超級計算機上進行了大量的前沿研究和多個領(lǐng)域的應(yīng)用。然而,至今沒有任何軟件能夠在神威太湖之光上模擬具有數(shù)萬億個粒子的超大N體系統(tǒng)。

本文在神威太湖之光上基于PHotoNs[1]設(shè)計并研發(fā)一款高效N體模擬軟件SwPHoToNs。針對粒子網(wǎng)格(Particle Mesh,PM)和快速多極子方法(Fast Multipole Method,FMM)的混合算法,提出多層次分解和負載均衡方案、執(zhí)行樹遍歷和引力計算的流水線策略以及向量化引力計算和超越函數(shù)優(yōu)化算法。通過以上并行算法設(shè)計和優(yōu)化策略,實施多個大規(guī)模并行模擬實驗,其中最大算例包含6 400億個粒子,并擴展到神威太湖之光超級計算系統(tǒng)的5 200 000個計算核心上。

1 背景介紹

1.1 N體問題相關(guān)算法

在一個粒子動力學(xué)系統(tǒng)中,每個粒子在物理力(如引力)的影響下與其他所有粒子進行相互作用。因此,如果直接按照定義對系統(tǒng)進行模擬,對于每個粒子需要進行O(N)個操作,總計算復(fù)雜度為O(N2),其中N為系統(tǒng)中的粒子總數(shù)。直接引力模擬算法又稱為粒子-粒子(Particle Particle,P2P)算法,具有很高的計算精度,但其隨著計算規(guī)模的增大,計算效率會降低。

為此,學(xué)者們提出了多種以輕微降低計算精度為代價來提高計算效率的方法,其中樹形法[2]得到廣泛應(yīng)用,而由BARNES和 HUT在1986年提出的Barnes-Hut分級樹算法就是一種典型的樹形法。在Barnes-Hut分級樹算法中,根據(jù)牛頓萬有引力定律,兩個粒子之間的相互作用力會隨著距離的增加而減小,在考慮遠程一組粒子對單個粒子的作用時,不是簡單地計算組內(nèi)每一個粒子的作用,而是用它們的總質(zhì)量和質(zhì)量中心進行近似計算。樹形法的基本思想是使用樹形數(shù)據(jù)結(jié)構(gòu)將所有粒子分為粒子集。在計算受力時,進行分級樹的遍歷,對節(jié)點和粒子的距離進行判定,若距離足夠遠則利用節(jié)點信息進行計算,否則繼續(xù)細分進行比較,直至遍歷到單個粒子。因此,可以極大地減少粒子之間的相互作用,將粒子間的計算復(fù)雜度降低至O(NlogaN)。

PM算法[3]通過對空間進行離散化進一步提高計算效率。但由于網(wǎng)格分辨率的限制,PM算法在建模近程力時存在不足,因此研究人員提出多種混合方法,如粒子-粒子-粒子網(wǎng)格(Particle-Particle/Particle-Mesh,P3M)算法和TreePM(Tree+Particle-Mesh)算法等,使用PM算法計算長程力,P2P算法或樹形算法計算短程力。這些算法的時間復(fù)雜度一般為O(NlogaN)。

由于模擬系統(tǒng)中粒子總數(shù)已增長至數(shù)千億甚至達到萬億,上述算法的時間復(fù)雜度已無法滿足性能要求,因此文獻[4]設(shè)計了在給定精度下時間復(fù)雜度為O(N)的FMM算法,其在某些方面與樹形法相似,但其使用的是勢而不是力。FMM算法通過層次劃分和位勢函數(shù)的多極子計算各點的位勢,其計算精度和劃分的層次有關(guān),因此可以達到任意的計算精度。

1.2 N體模擬

盡管隨著算法的改進,N體模擬的時間復(fù)雜度已降至O(NlogaN)甚至O(N),但從本質(zhì)上而言,N體模擬具有統(tǒng)計性,這意味著N值越大,相空間的采樣越精確和完整,從而可以得到更準確的結(jié)果。隨著粒子規(guī)模的爆炸型增長,傳統(tǒng)PC機的計算能力已經(jīng)遠遠無法滿足N體模擬的性能要求,因此需要大型并行計算機來解決該問題,而超級計算機正因為具有出色的計算和存儲能力,所以大規(guī)模的N體仿真被移植到這些計算設(shè)備上。早期的超級計算機多數(shù)為同構(gòu)架構(gòu),這意味著系統(tǒng)中只使用一種處理器或CPU核心。在當時的Caltech/JPL Mark III、Intel Touchstone Delta和Intel ASCI Red超級計算機上,WARREN和SALMON是第一批進行N體仿真的宇宙學(xué)家。1992年,他們用樹形法模擬了一個包含1 715萬個粒子的宇宙暗物質(zhì)模型,獲得了N體宇宙模擬領(lǐng)域的第一個戈登·貝爾獎[5]。在當時的Intel Touchstone Delta超級計算機上,模擬的持續(xù)性能達到5.1 GFLOPS～5.4 GFLOPS。為在Intel ASCI Red超級計算機上保持更好的負載平衡,WARREN和SALMON等人進一步改進了樹代碼,并進行一個包含3.22億多個粒子的宇宙學(xué)模擬,模擬了大爆炸之后宇宙中物質(zhì)的演化過程。仿真結(jié)果達到431 GFLOPS的持續(xù)性能,相當于超級計算機峰值性能的33%。在美洲豹超級計算機上,WARREN利用樹形法進行暗物質(zhì)暈的質(zhì)量函數(shù)模擬[6],仿真包含690億個粒子,單精度持續(xù)性能達到1.5 PFLOPS,相當于整個系統(tǒng)峰值性能的43%。ISHIYAMA等人報道了使用TreePM方法進行萬億個粒子的模擬,在K計算機上實現(xiàn)雙精度4.45 PFLOPS的性能,達到理論峰值性能的42%[7]。

隨著摩爾定律的終結(jié),將配備了GPU、FPGA和GRAPE等加速器的超級計算機用于為科學(xué)和工程應(yīng)用提供計算能力,這些計算系統(tǒng)被稱為異構(gòu)超級計算機。雖然面對仿真軟件的復(fù)雜性、加速器與主機CPU之間數(shù)據(jù)傳輸?shù)睦щy性以及在異構(gòu)超級計算機上編程的挑戰(zhàn)性等問題,但是應(yīng)用程序可以從這些加速器中獲得更好的高性能計算效果。GRAPE就是一種專門用于引力計算的加速器,MAKINO和FUKUSHIGE等人成功地進行了一系列基于GRAPE的天體物理N體模擬[8-11]。在過去的十幾年中,由于GPU一直是頂級超級計算系統(tǒng)的主要組成部分,因此產(chǎn)生了大量面向GPU的N體問題解算代碼。HAMADA等人在2009年利用樹形法和FMM對256個GPU集群進行42 TFLOPS、包含16億個粒子的N體模擬[12],仿真結(jié)果表明GPU集群在成本/性能、功耗/性能以及物理尺寸/性能方面都優(yōu)于PC集群。2010年,HAMADA和NITADORI將代碼移植到DEGIMA集群中[13],DEGIMA集群是由288個GTX295 GPU組成的PC機集群,可以用于更大規(guī)模的天體物理模擬。仿真包含30多億個粒子,性能達到190 TFLOPS,計算效率為35%。BéDORF等人于2014年使用樹形法對包含2 420億個粒子的銀河系進行了單精度24.77 PFLOPS的N體模擬[14]。泰坦超級計算機在單精度下的理論峰值性能為73.2 PFLOPS,N體模擬的計算效率達到33.8%。

神威太湖之光超級計算機是第一個峰值性能超過100 PFLOPS的超級計算機,其將不同種類的內(nèi)核集成到一個處理器中,具有處理進程管理、計算加速和內(nèi)存訪問操作的異構(gòu)架構(gòu)。雖然目前已在神威太湖之光上進行了許多大氣、生物學(xué)、計算流體動力學(xué)等前沿模擬,但大規(guī)模的宇宙學(xué)模擬一直未見報道。IWASAWA等人[15]使用粒子模擬平臺FDPS對行星環(huán)進行N體模擬,算例規(guī)模達到80億個粒子,共擴展到4 096個Sunway節(jié)點。實測性能達到理論峰值性能的35%。本文設(shè)計的SwPHoToNs是為了利用整個超級計算系統(tǒng)來執(zhí)行超大規(guī)模天體物理N體模擬而設(shè)計。在宇宙學(xué)模擬中,本文成功地將代碼擴展到20 000多個Sunway節(jié)點,相當于可用計算資源的一半。

1.3 神威太湖之光超級系統(tǒng)

神威太湖之光是我國自主研發(fā)的超級計算系統(tǒng),具有超過1 000萬個計算核心,最高性能為125 PFLOPS,持續(xù)性能為93 PFLOPS。整個系統(tǒng)由40 960個SW26010處理器組成,這是上海高性能集成電路設(shè)計中心設(shè)計的一款具有260個異構(gòu)內(nèi)核的多核處理器,最高性能為3.06 TFLOPS,性能功率比為10 GFLOPS/W。整個處理器分為4個核心組(CG),每個核心組包含1個管理處理單元(MPE)、1個智能內(nèi)存處理單元(IMPE)和64個計算處理單元(CPE)。SW26010體系結(jié)構(gòu)如圖1所示。

圖1 SW26010體系結(jié)構(gòu)Fig.1 Architecture of SW26010

MPE是一個完整的64位RISC內(nèi)核,頻率為1.45 GHz,主要用于處理程序的流量控制、I/O和通信功能。CPE采用更簡化的微體系結(jié)構(gòu),能最大限度地聚合計算能力。每個CPE都有16 KB的L1指令緩存和64 KB的本地數(shù)據(jù)內(nèi)存(Local Data Memory,LDM),可以配置為用戶控制的快速緩沖區(qū)。基于三層結(jié)構(gòu)(REG-LDM-MEM)的內(nèi)存層次結(jié)構(gòu)由FANG等人[16]提出。CPE可以直接訪問帶寬為8 GB/s的全局內(nèi)存,也可通過REG-LDM-MEM內(nèi)存層次結(jié)構(gòu)獲得更高的帶寬。每個CPE有兩個管道來處理指令的解碼和執(zhí)行。指令重排方案可以降低指令序列的依賴性,從而提高指令執(zhí)行效率。在每個CPE集群中,一個包含8行通信總線和8列通信總線的網(wǎng)狀網(wǎng)絡(luò)能夠在寄存器級上實現(xiàn)快速的數(shù)據(jù)通信和共享,從而允許CPE之間高效的數(shù)據(jù)重用和通信。

神威太湖之光采用基于高密度運算超級節(jié)點的網(wǎng)絡(luò)拓撲結(jié)構(gòu)。1個超級節(jié)點由256個處理器組成,1個超級節(jié)點內(nèi)的所有處理器通過定制的網(wǎng)絡(luò)交換機完全連接,實現(xiàn)高效的all-to-all通信,連接所有超級節(jié)點的網(wǎng)絡(luò)拓撲是一個使用自定義高速互連網(wǎng)絡(luò)的兩級胖樹,網(wǎng)絡(luò)鏈路帶寬為16 GB/s,對分帶寬達到70 TB/s。

為適應(yīng)SW26010處理器,本文開發(fā)神威太湖之光軟件系統(tǒng),包括定制的64位Linux操作系統(tǒng)內(nèi)核、全局文件系統(tǒng)、編譯器、作業(yè)管理系統(tǒng)等。Sunway編譯系統(tǒng)支持C/C++、Fortran以及主流的并行編程標準(如MPI和OpenACC),提供了一個輕量級線程庫Athread,并允許程序員執(zhí)行更細粒度的優(yōu)化。

為神威太湖之光設(shè)計的應(yīng)用程序通常采用兩級方法將不同算法映射到管理和計算核心,類似于Summit、Piz Daint、Titan等異構(gòu)超級計算機,提供了兩種使用CPE執(zhí)行更細粒度優(yōu)化的應(yīng)用程序:1)Sunway OpenACC,其是從OpenACC 2.0標準擴展而來,可以在CPE集群上管理并行任務(wù);2)輕量級線程庫Athread,可以用于控制每個CPE。Sunway OpenACC利用CPE網(wǎng)格的快捷方式,但其缺少一些用于細粒度調(diào)優(yōu)的重要特性,如寄存器通信、向量化等。相反地,使用Athread接口可以利用這些特性及多核處理器的計算能力。因此,本文采用“MPI+Athread”作為并行化解決方案。

2 SwPHoToNs算法描述及平臺實現(xiàn)

2.1 SwPHoToNs框架

PHoToNs是針對大規(guī)模并行宇宙模擬[1]而設(shè)計,SwPHoToNs在此基礎(chǔ)上針對Sunway平臺進行設(shè)計和優(yōu)化。PHoToNs軟件將PM和FMM兩種常用算法相結(jié)合,如圖2所示。

圖2 PM+FMM算法Fig.2 PM+FMM algorithm

本文使用PM算法用于遠程引力計算,FMM算法用于短程引力計算,因此將P2P的方程修改為式(1):

(1)

其中,s為分離度,Rs為特征尺度。因此,短程力約在5Rs的范圍外消失,修正后的方程引入了PM上的長距離高斯平滑,但卻引入了兩個計算非常耗時的誤差函數(shù)(erf)和指數(shù)函數(shù)(exp)。

FMM算法簡單而言是將大系統(tǒng)切成小立方體,只計算小立方體之間的作用力。因此,1萬個粒子可能只切成10個立方體,從而大幅降低計算復(fù)雜度,并且采用樹結(jié)構(gòu)代碼有利于快速計算和并行劃分[17]。在FMM算法中,粒子被組織成K-D樹,樹的葉子是一個粒子包。樹的每個節(jié)點存儲2個多極矩M和L,將粒子間的力轉(zhuǎn)化為多極矩之間的相互作用,其中計算引力的算符有P2M、M2M、M2L、L2L、L2P和P2P。

P2P算符是兩個包中粒子之間的交互,是N體直接求解的核心部分。其他操作符則用于不同樹節(jié)點的矩之間的交互。本文使用對偶樹遍歷來確定包和樹節(jié)點之間的交互關(guān)系。對偶樹遍歷的優(yōu)點是可以減少節(jié)點對的數(shù)量[19]。節(jié)點或包之間的打開角度(作為控制參數(shù),由兩者距離和包的寬度確定)為0.4,如果遞歸檢查其子節(jié)點的交互關(guān)系,則此節(jié)點需滿足該打開角度。

為提高編程的并行效果,本文將引力計算分解為一系列任務(wù),對在CPE上處理的6個FMM算符中計算量最大的P2P和M2L構(gòu)造一個任務(wù)池。

2.2 區(qū)域分解與并行優(yōu)化

如圖3所示,計算框由K-D樹分解。每一個計算節(jié)點作為頂層樹的一個節(jié)點,對應(yīng)一個連續(xù)的空間區(qū)域和一個長方體,也是由局部粒子構(gòu)成的局部樹的根節(jié)點。因此,所有粒子都鏈接到一個完整的全局K-D樹中。在圖3(c)中,點劃線圈內(nèi)粒子間的相互作用由P2P處理,虛線圈內(nèi)粒子間的相互作用由M2L處理。

圖3 區(qū)域分解 Fig.3 Domain decomposition

在所有進程中記錄頂層樹,邊界信息需要在進程之間進行通信,在此使用一個額外的遍歷根據(jù)打開角度估計需要發(fā)送給其他進程的本地樹節(jié)點數(shù)量。實際上,引力的長-短分裂保證了邊界交換只發(fā)生在截止半徑內(nèi)。

在引力計算中,最大的工作量集中在兩個部分,粒子間的引力作用P2P和M2L。由于這兩個算符包含的計算量大致相同,因此算符的出現(xiàn)次數(shù)是負載計算單位。通過測量P2P和M2L的數(shù)量(M2L的數(shù)量相比P2P的數(shù)量少)來估計域之間的工作負載平衡性。

在將交互計算任務(wù)分配給每個MPE控制的CPE集群前,通過執(zhí)行樹遍歷獲取交互任務(wù)列表。交互任務(wù)在CPE之間均勻分布,以確保細粒度的負載平衡。因此,每個CPE都可以通過DMA在主存和LDM之間傳輸數(shù)據(jù),并獨立處理交互任務(wù)。

在樹形法中,如果需要對質(zhì)點進行引力計算,則需要遍歷樹獲取交互信息。大多數(shù)GPU代碼[13-14]是通過GPU進行遍歷和引力計算過程,本文提出一種流水線策略,允許在MPE和CPE集群上同時執(zhí)行遍歷和引力計算過程。使用數(shù)組存儲粒子包的ID,并在遍歷過程中判斷是否需要計算這些粒子包。一旦具有足夠的交互任務(wù)(如8 192個包),就啟動CPE集群進行引力計算,同時MPE會一直執(zhí)行遍歷過程,通過調(diào)整參數(shù)使遍歷過程和引力計算處于平衡狀態(tài)并且相互重疊。

根據(jù)牛頓萬有引力定律,粒子之間的相互作用為all-to-all。因此,在力的計算過程中,目標進程需要與其他所有進程進行通信,這是影響N體求解器可擴展性的主要因素。為有效增強代碼的可擴展性,將通信和FFT過程與P2P過程重疊。在MPE上處理通信和FFT時,驅(qū)動CPE集群進行P2P計算,此時不對遍歷過程和引力計算進行相互掩蓋,而是在任務(wù)池中存儲所有的交互任務(wù)后驅(qū)動CPE集群,該策略具有較好的可擴展性。如圖4所示,此時的掩蓋策略與在本地進程進行引力計算時的掩蓋遍歷方式不同,因此針對兩種計算情況構(gòu)建了不同的任務(wù)池。

圖4 針對兩種情況的掩蓋方式Fig.4 Cover-up modes of two situations

3 核心函數(shù)的眾核優(yōu)化

3.1 優(yōu)化策略

如表1所示,P2P算符在整個仿真過程中所消耗的時間是最多的。P2P和P2P_ex函數(shù)的作用都是通過粒子之間的距離及粒子質(zhì)量計算出粒子的加速情況,是天體模擬的核心引力計算函數(shù)。區(qū)別是兩者使用的數(shù)據(jù)不同,P2P負責計算自身MPE中粒子的影響,P2P_ex負責計算其他MPE中的粒子的影響。因此,通過這兩部分的從核計算進行性能優(yōu)化。

表1 函數(shù)占用時間比較Table 1 Comparison of function occupancy time

在計算粒子間的相互作用時,由于力的作用是相互的,因此計算結(jié)果應(yīng)累加到兩者之上,但是該做法會導(dǎo)致所有粒子需記錄已受過哪些粒子的作用,進行的判斷與記錄會耗費大量的時間和空間。因此,本文將計算結(jié)果只累加到一個粒子上,將計算中累加計算結(jié)果的一種粒子稱為計算粒子,而另一種稱為影響粒子。使用如圖5所示的偽代碼描述計算粒子包與影響粒子包之間的P2P計算和從核化過程。

圖5 P2P函數(shù)偽代碼Fig.5 Pseudo-code of P2P function

對于眾核計算程序的優(yōu)化,通過提高DMA帶寬以降低DMA通信時間是一個關(guān)鍵策略,而另一個關(guān)鍵策略是將DMA通信時間隱藏在計算時間下,通過測試發(fā)現(xiàn)DMA通信時間比計算時間短,因此建立同等大小的數(shù)據(jù)緩沖池,將DMA通信時間進行隱藏,異步DMA方案的具體實現(xiàn)過程如圖6所示。

圖6 異步DMA方案Fig.6 Asynchronous DMA scheme

3.2 超越函數(shù)的優(yōu)化

由于Sunway處理器在CPE上沒有對超越函數(shù)進行向量化的指令,因此只能將向量存儲到數(shù)組中,運用查表+多項式逼近的方法對exp函數(shù)和erf函數(shù)逐一進行求解。由于向量無法直接查表,因此針對exp函數(shù)的多項式(2)～式(4)和erf函數(shù)的定義式(5)、多項式(6)和多項式(7)展開,只采用多項式逼近的方法實現(xiàn)這些函數(shù)的SIMD版本,雖然會在erf函數(shù)的部分輸入域時產(chǎn)生誤差(最大誤差為1.5×10-7),但加速效果較好。

exp(x)=2k×exp(r)

(2)

(3)

c(r)=r-(P1×r+P2×r+…+P5×r)

(4)

(5)

erf(x)≈1-(a1t+a2t2+…+a5t5)e-x2

(6)

(7)

其中,計算式(2)中的2k時,本文針對浮點型數(shù)據(jù)的存儲形式,對于整型數(shù)據(jù)的指數(shù)k加偏差后位移,再保存為浮點型數(shù)據(jù),即為2k的結(jié)果,以此避免了冪函數(shù)計算。

3.3 眾核優(yōu)化的加速效果

優(yōu)化策略產(chǎn)生的加速效果與P2P算符的運行時間如圖7所示,其中柱狀表示一次迭代的運行時間,基準是一個MPE上對于執(zhí)行P2P原始代碼的運行時間。將計算任務(wù)分配給其控制的CPE集群,在解決了部分訪址非連續(xù)問題后獲得了38.27倍的加速。雖然理論上SIMD優(yōu)化可將計算性能提升4倍,但是實際上除了加載和存儲的開銷外,exp函數(shù)和erf函數(shù)是計算中主要耗時的部分,約占總計算量的90%,且不存在相應(yīng)的SIMD函數(shù),只能將向量封裝入數(shù)組中分別進行計算。因此,將函數(shù)循環(huán)展開并進行向量化后,只獲得1.30倍的加速。在使用多項式逼近方法實現(xiàn)向量化超越函數(shù)后,既避免了向量的封裝與加載,又加快了兩個函數(shù)的計算。雖然產(chǎn)生了誤差,但是誤差大小是可接受的,且能通過泰勒公式更高階的展開進行控制。經(jīng)過函數(shù)計算方法的改進,獲得了2.20倍的加速。將DMA通信時間隱藏后,獲得了1.51倍的加速。最終實現(xiàn)了165.29倍的加速,在該規(guī)模的模擬中,對于P2P函數(shù)單次計算的時間從6 804.800 s減少到41.170 s。

圖7 優(yōu)化策略產(chǎn)生的加速效果與P2P算符的計算時間Fig.7 The acceleration effect produced by the optimizationstrageties and calculation time of P2P operator

4 實驗結(jié)果與分析

4.1 宇宙學(xué)模擬效果分析

宇宙大規(guī)模結(jié)構(gòu)是由初始的微小隨機密度波動演化而來。根據(jù)Zel’dovich近似方法[20],通過粒子遷移建立仿真的初始條件。該方法是由初始功率譜編譯卷積,早期可以用微擾理論進行估算。本文主要研究一個具有臨界密度ΩM=0.25、Ω∧=0.75、σ8=0.8和哈勃常數(shù)參數(shù)H0=70 km/s的平面模型。在紅移時用一個巨大的周期盒(z=49)進行模擬,該尺度可以滿足宇宙學(xué)天空測量體積的統(tǒng)計量要求。具體而言,體積接近于8h-3Gpc3(其中,pc為秒差距,是天文學(xué)上的一種長度單位,h是以100 km/s/Mpc為單位的哈勃常數(shù)),從z=49模擬演變到z=0,形成的宇宙密度場如圖8所示。本文共使用109個粒子參與到宇宙演變研究的實際模擬過程中,最大的基準測試算例包含6 400億個粒子,總質(zhì)量接近于太陽質(zhì)量的108倍,這是首次在Sunway平臺上實現(xiàn)同類宇宙學(xué)模擬。

圖8 z為0且邊長為2 Gpc時的宇宙密度場Fig.8 Cosmic density field when the z is 0 andthe side length is 2 Gpc

4.2 計算性能與可擴展性分析

由于程序執(zhí)行過程中的絕大部分浮點操作都集中在P2P算符中,因此忽略程序中其他的浮點計算,僅統(tǒng)計P2P核心中的浮點運算次數(shù)作為宇宙模擬的總浮點計算次數(shù)。使用該方法計算出的總浮點計算次數(shù)略小于程序的真實浮點計算次數(shù)。在本文實現(xiàn)中,每次P2P過程中包含3次減法運算、12次乘法運算、6次加法運算、1次指數(shù)函數(shù)(exp)運算、1次誤差函數(shù)(erf)運算和1次倒數(shù)平方根函數(shù)(rsqrt)運算。為便于與其他工作進行性能比較,本文對rsqrt函數(shù)使用文獻[13-14]中提到的4次浮點運算,對exp和erf函數(shù)分別使用21次和15次浮點運算,因此1次P2P交互共產(chǎn)生61次浮點運算。

程序總浮點運算次數(shù)的計算公式為:

FFLOPS=I×C×SSteps

(8)

其中,I表示每次迭代計算中粒子間的相互作用(P2P)總次數(shù),C表示每次P2P過程中包含的浮點運算次數(shù),SSteps表示程序運行的總迭代步數(shù)。

程序持續(xù)浮點運算性能的計算公式具體如下:

PPerformance=FFLOPS/T

(9)

其中,T是程序運行的總時間。本文采用的最大算例中包含6 400億個粒子,每次迭代計算中的粒子相互作用總數(shù)約為3.2×1016,每次P2P中包含61次浮點操作,總迭代步數(shù)為64,程序運行總時間約為4 243.2 s。將以上數(shù)據(jù)代入式(8)和式(9),計算得到程序的持續(xù)性能達到29.44 PFLOPS,計算效率為48.3%,并行效率約為85%,說明大部分通信時間與計算時間相重疊。

圖9(a)、圖9(b)為弱可擴展性性能與并行效率測試結(jié)果。對宇宙學(xué)天體演變測量中的統(tǒng)計量進行一系列以盒徑為2h-1Gpc標定的宇宙學(xué)模擬。本文從1 024個節(jié)點開始,逐漸擴展到20 000個節(jié)點,每個節(jié)點包含4 000萬個粒子(每個CG有1 000萬個粒子),其中最大的粒子總數(shù)為6 400億,實現(xiàn)了29.44 PFLOPS的持續(xù)計算速度和84.6%的并行效率。圖9(c)為強可擴展性加速比測試結(jié)果,粒子數(shù)量為549億,在20 000個節(jié)點上的并行效率達到85.6%。

圖9 大規(guī)模實驗性能測試結(jié)果Fig.9 Performance test results of large-scale experiment

4.3 計算性能對比

表2給出了HAMADA & NITADORI[16]、ISHIYAMA[7]、 BéDORF[14]、FDPS[15]和SwPHoToNs大規(guī)模宇宙學(xué)N體模擬平臺的性能比較,同時介紹了所有模擬平臺使用的計算機名稱、計算節(jié)點體系結(jié)構(gòu)以及內(nèi)存帶寬浮點比率(B/F)。雖然SwPHoToNs硬件平臺在B/F方面處于劣勢,但相比其他平臺具有更高的計算效率和浮點性能,并且本文對SwPHoToNs擴展到神威太湖之光全系統(tǒng)的計算結(jié)果進行了預(yù)估。

表2 大規(guī)模宇宙學(xué)N體模擬平臺的性能比較Table 2 Performance comparison of large-scale cosmological N-body simulation platforms

5 結(jié)束語

本文設(shè)計了在神威太湖之光上進行大規(guī)模宇宙學(xué)N體模擬的軟件SwPHoToNs。為將SwPHoToNs算法映射到神威太湖之光計算系統(tǒng)的5 200 000個計算核心中,提出多層次分解和負載均衡方案、樹遍歷和引力計算的流水線策略以及向量化引力計算算法等多種實現(xiàn)并行可擴展性和提高計算效率的技術(shù),并在神威太湖之光上的5 200 000個計算核心上進行了宇宙學(xué)模擬,模擬算例包含6 400億個粒子,弱可擴展性并行效率為84.6%,持續(xù)計算速度為29.44 PFLOPS。與BéDORF[14]宇宙學(xué)N體模擬平臺相比,SwPHoToNs將最大模擬算例中的粒子數(shù)提高了2.6倍,浮點計算速度由單精度24.77 PFLOPS提高到雙精度29.44 PFLOPS,浮點計算效率高達48.3%。由于本文模擬工作建立在半機神威太湖之光系統(tǒng)上,因此還沒有實現(xiàn)真正意義上的萬億級天體物理模擬,預(yù)估整機神威太湖之光系統(tǒng)能夠支持包含多達1.2萬億個粒子的宇宙學(xué)模擬,持續(xù)計算速度將達到56.00 PFLOPS。