基于神威太湖之光的宇宙學(xué)多體模擬

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

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

0 概述

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

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

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

1 背景介紹

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

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

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

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

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

1.2 N體模擬

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

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

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

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

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

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

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

神威太湖之光采用基于高密度運(yùn)算超級(jí)節(jié)點(diǎn)的網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)。1個(gè)超級(jí)節(jié)點(diǎn)由256個(gè)處理器組成,1個(gè)超級(jí)節(jié)點(diǎn)內(nèi)的所有處理器通過(guò)定制的網(wǎng)絡(luò)交換機(jī)完全連接,實(shí)現(xiàn)高效的all-to-all通信,連接所有超級(jí)節(jié)點(diǎn)的網(wǎng)絡(luò)拓?fù)涫且粋€(gè)使用自定義高速互連網(wǎng)絡(luò)的兩級(jí)胖樹,網(wǎng)絡(luò)鏈路帶寬為16 GB/s,對(duì)分帶寬達(dá)到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以及主流的并行編程標(biāo)準(zhǔn)(如MPI和OpenACC),提供了一個(gè)輕量級(jí)線程庫(kù)Athread,并允許程序員執(zhí)行更細(xì)粒度的優(yōu)化。

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

2 SwPHoToNs算法描述及平臺(tái)實(shí)現(xiàn)

2.1 SwPHoToNs框架

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

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

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

(1)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 優(yōu)化策略

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

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

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

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

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

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

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

由于Sunway處理器在CPE上沒(méi)有對(duì)超越函數(shù)進(jìn)行向量化的指令,因此只能將向量存儲(chǔ)到數(shù)組中,運(yùn)用查表+多項(xiàng)式逼近的方法對(duì)exp函數(shù)和erf函數(shù)逐一進(jìn)行求解。由于向量無(wú)法直接查表,因此針對(duì)exp函數(shù)的多項(xiàng)式(2)～式(4)和erf函數(shù)的定義式(5)、多項(xiàng)式(6)和多項(xiàng)式(7)展開,只采用多項(xiàng)式逼近的方法實(shí)現(xiàn)這些函數(shù)的SIMD版本,雖然會(huì)在erf函數(shù)的部分輸入域時(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)

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

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

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

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

4 實(shí)驗(yàn)結(jié)果與分析

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

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

圖8 z為0且邊長(zhǎng)為2 Gpc時(shí)的宇宙密度場(chǎng)Fig.8 Cosmic density field when the z is 0 andthe side length is 2 Gpc

4.2 計(jì)算性能與可擴(kuò)展性分析

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

程序總浮點(diǎn)運(yùn)算次數(shù)的計(jì)算公式為:

FFLOPS=I×C×SSteps

(8)

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

程序持續(xù)浮點(diǎn)運(yùn)算性能的計(jì)算公式具體如下:

PPerformance=FFLOPS/T

(9)

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

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

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

4.3 計(jì)算性能對(duì)比

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

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

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

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