儲根深,何遠(yuǎn)杰,白 鶴,陳丹丹,任 帥,賈麗霞,吳 石,賀新福,楊 文,胡長軍
(1.智能超算融合應(yīng)用技術(shù)教育部工程研究中心,北京 100083;2.北京科技大學(xué),北京 100083;3.中國原子能科學(xué)研究院,北京 102413)
材料輻照損傷模擬是數(shù)值堆軟件的重要內(nèi)容,采用多尺度模擬方法是研究材料輻照損傷問題的一個重要手段[1-6]。其中,動力學(xué)蒙特卡羅(KMC)方法是多尺度模擬中最為重要的方法之一,被廣泛用于研究微觀尺度的材料演化行為。不同于同樣原子尺度的分子動力學(xué)(MD)方法,KMC方法將原子運動軌跡粗化為體系組態(tài)躍遷,使該方法可在保持原子級別精度下有效將模擬的時間跨度從原子振動的尺度提高到組態(tài)躍遷的尺度,即實現(xiàn)秒甚至年量級的材料輻照行為模擬[7-9]。
KMC方法目前已有不少研究[10-19],包括原子動力學(xué)蒙特卡羅(AKMC)方法、實體動力學(xué)蒙特卡羅(OKMC)方法、事件蒙特卡羅(EKMC)方法等。其中,AKMC方法以原子為基本對象,考慮原子和缺陷之間的躍遷導(dǎo)致的微結(jié)構(gòu)演化行為,其在空間尺度上更精確且方便與同樣以原子為操作對象的MD方法進(jìn)行耦合。但在面向?qū)嶋H應(yīng)用需求時,其仍面臨著內(nèi)存限制和復(fù)雜的計算量等挑戰(zhàn)。隨著模擬體系的增大,模擬所需內(nèi)存需求的矛盾也開始凸顯出來。同時,雖然KMC模擬在時間尺度上相對于MD提升了數(shù)個數(shù)量級,但計算過程仍十分耗時,特別是在大體系模擬時。解決計算和內(nèi)存的雙重矛盾最有效的方式就是并行化,通過計算任務(wù)分配以并行計算,并通過分布式內(nèi)存方案減少單個節(jié)點上的內(nèi)存占用?;诖?,本文將采用AKMC方法設(shè)計數(shù)值堆材料輻照損傷動力學(xué)并行模擬軟件MISA-AKMC。
和通用的蒙特卡羅思想類似,KMC方法的核心思想是隨機抽取事件。若可確定體系當(dāng)前狀態(tài)下所有可能發(fā)生的事件及其發(fā)生的概率,即可將以這些事件發(fā)生的概率作為權(quán)重隨機抽取事件執(zhí)行以改變體系的狀態(tài),以此推進(jìn)體系向前演化。圖1示出了一種通用的KMC執(zhí)行流程,其在時間步循環(huán)內(nèi),一次計算所有可能的事件的躍遷速率(也即事件概率),隨后通過隨機數(shù)選擇其中的1個事件,并行執(zhí)行該事件;如此進(jìn)行時間步的循環(huán),即可推進(jìn)系統(tǒng)向前演化,直到模擬到預(yù)期的時間。在材料微觀演化模擬過程中,事件可以是缺陷的位置躍遷或新缺陷的生成等。
圖1 KMC執(zhí)行流程Fig.1 KMC execution process
由于串行的KMC中,體系完成的連續(xù)兩次演化是獨立、無記憶的,這個過程是一種典型的馬爾可夫過程(Markov process)。且計算兩個相鄰的事件之間具有依賴關(guān)系,即事件的發(fā)生是串行的,位于當(dāng)前狀態(tài)的體系只發(fā)生1個事件,該事件發(fā)生且執(zhí)行完成后才會進(jìn)入到下一個事件的發(fā)生過程。這種依賴性給KMC并行化帶來了很大阻礙。
目前大多數(shù)已有的并行KMC方法基本是基于區(qū)域分解的,將模擬區(qū)域劃分為多個子區(qū)域,并將其分配到不同的計算單元上。同時,這些方法近似認(rèn)為這些子區(qū)域之間的距離足夠遠(yuǎn),之間的相互影響不大,因此不同子區(qū)域上的事件可同時發(fā)生。這就是現(xiàn)有的基于空間分解的并行KMC算法的基本思想。典型的并行KMC算法有基于空間分解方案的SL(sub-lattice)KMC[20-21]、SP(synchronous parallel)KMC[22-24]、SR(synchronous relaxation)KMC[25-26]等。其中,SP KMC算法和SR KMC算法在執(zhí)行過程中,每個時間步均需引入昂貴的全局通信,影響計算性能;而SL KMC算法執(zhí)行過程不必全局通信,這也成為本文的并行KMC算法的首選。
在SL KMC算法中,為實現(xiàn)正確的并行計算,需解決兩個問題:1) 邊界處的事件沖突;2) 不同計算單元上的時間同步(實際上,在其他兩種并行SP KMC和SR KMC算法中,也存在這兩個問題)。為解決邊界處的事件沖突,一般采用子區(qū)域著色方法。其將每個計算單元上的子區(qū)域再分為多個扇區(qū),不同位置的扇區(qū)標(biāo)記為不同顏色,且不同計算單元中處于同一位置的扇區(qū)標(biāo)記為相同顏色。通過著色使相同顏色的扇區(qū)在空間上不相鄰。各計算單元均依次選擇相同顏色的扇區(qū)進(jìn)行計算。KMC并行中的事件沖突問題如圖2a所示,當(dāng)兩個計算單元上的事件發(fā)生在邊界處的鄰近位置時,可能會產(chǎn)生事件沖突。為解決這種沖突,可針對計算單元對應(yīng)的子區(qū)域劃分扇區(qū)并著色。著色算法解決事件沖突如圖2b所示,圖中的一維著色示例,黃色與橙色的扇區(qū)分別被間隔開來,依次選擇黃色扇區(qū)和橙色扇區(qū)進(jìn)行計算。當(dāng)然,在實際程序中,考慮到模擬的三維空間和計算單元的三維邏輯分布,將每個計算單元內(nèi)對應(yīng)的子區(qū)域劃分為8個扇區(qū),分別著色為S0、S1、…、S7,如圖2c所示。針對不同計算單元上的時間同步問題,SL并行KMC方法采用時間閾值的方式。給定時間閾值T,各計算單元上進(jìn)行各自的KMC演化模擬,直到計算單元中當(dāng)前扇區(qū)的時間增量達(dá)到閾值。在時間增量達(dá)到閾值后,即可與鄰居進(jìn)程進(jìn)行通信以交換數(shù)據(jù)(如同步躍遷對鄰居進(jìn)程的更改,從鄰居進(jìn)程處同步ghost區(qū)域等)。通信完成后,各計算單元切換到下一扇區(qū),進(jìn)行下一扇區(qū)的計算模擬,直到計算單元上的總時間增量達(dá)到閾值2T,如此往復(fù)。通過該時間閾值,計算單元在其時間增量到達(dá)nT時進(jìn)行通信(其中n為正整數(shù)),通信完成后進(jìn)行下一扇區(qū)的模擬(同一計算單元上的各扇區(qū)的模擬依次進(jìn)行)。
圖2 KMC并行中的事件沖突(a)、著色算法解決事件沖突(b)和三維空間上的扇區(qū)劃分與著色(c)Fig.2 Event conflict in KMC parallelism (a), solving even conflict by coloring method (b), and sector dividing and coloring in three-dimensional space (c)
圖3示出了MISA-AKMC程序及并行KMC框架架構(gòu),為便于KMC模型的集成與可擴展,并方便軟件的核心優(yōu)化,將MISA-AKMC并行部分分離開,形成KMC并行框架。該框架采用空間分解的并行思想實現(xiàn)了SL并行KMC方法,并提供額外的接口以供實現(xiàn)其他的并行KMC算法??蚣茚槍诵牡牟⑿胁糠?,提供了多種加速優(yōu)化方法,包括局部躍遷速率更新方法和轉(zhuǎn)發(fā)通信算法優(yōu)化等。此外,框架還內(nèi)置了晶格點和缺陷存儲等數(shù)據(jù)結(jié)構(gòu)的表示,并向上層模型提供了KMC模型集成接口。
圖3 MISA-AKMC程序及并行KMC框架架構(gòu)Fig.3 Architecture of MISA-AKMC program and parallel KMC framework
作為并行KMC框架的核心部分,本文實現(xiàn)的sub-lattice并行KMC算法計算流程如圖4所示。在本文的sub-lattice算法實現(xiàn)中,先依據(jù)通信閾值T計算出通信次數(shù)n,以此進(jìn)行外層的時間步循環(huán)。在時間步循環(huán)內(nèi)部,依次進(jìn)行計算單元內(nèi)的8個扇區(qū)的迭代:1) 扇區(qū)內(nèi)的事件循環(huán),對于每個扇區(qū),重復(fù)在該扇區(qū)選擇并執(zhí)行時間,直到該扇區(qū)的時間增量達(dá)到預(yù)設(shè)的閾值T;2) 通信,在一扇區(qū)s的計算完成后,還需進(jìn)行通信以準(zhǔn)備下一個扇區(qū)的計算,包括將扇區(qū)s對鄰居計算單元上的更改通信給鄰居計算單元,同時為計算下一扇區(qū),還需從鄰居計算單元上通信獲取下一扇區(qū)的ghost區(qū)域;3) 進(jìn)入下一扇區(qū)的計算。
圖4 MISA-AKMC并行KMC框架中的并行算法流程Fig.4 Parallel algorithm flow in parallel KMC framework of MISA-AKMC
為進(jìn)一步提升MISA-AKMC程序性能,本文還針對其中核心并行部分進(jìn)行了優(yōu)化。需強調(diào)的是,由于MISA-AKMC采用并行算法和計算模型分離的設(shè)計,針對并行框架部分的優(yōu)化和改動不會影響到KMC模型部分的實現(xiàn),這是本文設(shè)計的KMC程序的一大優(yōu)勢。
1) 局部更新優(yōu)化
為模擬的精確性,將通信閾值T取得較小,一般近似的等于單次事件產(chǎn)生的時間增量。這必然導(dǎo)致扇區(qū)內(nèi)的事件循環(huán)中事件發(fā)生的數(shù)量很少,這些事件影響的區(qū)域也有限。故不必在每個事件發(fā)生后“全局地”更新該扇區(qū)的躍遷速率,而只需更新上一個事件影響的區(qū)域內(nèi)的躍遷速率。計算躍遷速率是KMC模擬中最耗時的,通過局部躍遷速率更新的優(yōu)化方法,可極大減少更新躍遷速率的計算開銷。
考慮到單次躍遷影響區(qū)域可能覆蓋到其他扇區(qū),甚至其他計算單元對應(yīng)的區(qū)域,如圖5所示。所以,事件發(fā)生后,不僅需更新當(dāng)前扇區(qū)的總躍遷速率,還需更新受影響的其他扇區(qū)的躍遷速率和鄰居計算單元上對應(yīng)扇區(qū)的躍遷速率。為此,將事件的影響區(qū)域劃分為3個部分。其中,第1個部分為影響區(qū)域和當(dāng)前計算扇區(qū)的交集,記為S;第2部分為影響區(qū)域與當(dāng)前計算單元上但除當(dāng)前計算扇區(qū)以外的區(qū)域的交集,記為C;第3部分為影響區(qū)域與ghost區(qū)域的交集(影響區(qū)域不會超過ghost區(qū)域的邊界),ghost區(qū)域是計算單元模擬區(qū)域外面的一層額外區(qū)域,但該區(qū)域是來自鄰居進(jìn)程上的副本,將該區(qū)域記為G。在更新躍遷速率時,針對S區(qū)域,直接更新;針對C區(qū)域,找到該區(qū)域所屬的扇區(qū),再更新該扇區(qū)的總躍遷速率;對于區(qū)域G,先更改ghost區(qū)域的狀態(tài),后續(xù)通信時,將更改同步到鄰居進(jìn)程,再在鄰居進(jìn)程上重新計算這部分受影響區(qū)域的總躍遷速率。
圖5 KMC事件的可能影響區(qū)域Fig.5 Possible effect region of KMC event
2) 進(jìn)程通信優(yōu)化
在三維KMC模擬中,每個扇區(qū)與7個鄰居計算單元相鄰。在圖4所示的算法中,針對兩個通信過程,每個扇區(qū)均與7個鄰居計算單元進(jìn)行通信。直接的通信方式,完成1輪8個扇區(qū)的計算,需進(jìn)行224次通信(7鄰居×8扇區(qū)×2個通信操作×2(發(fā)送/接收操作)),如SPPARKS程序[27]。在Wu等[28]開發(fā)的Crystal-KMC程序中,提出了通信合并的方案,在1個扇區(qū)計算完成后的同步模擬區(qū)域的通信可與下一扇區(qū)計算前的同步ghost區(qū)域的通信進(jìn)行合并,這樣可將通信次數(shù)降為182。
實際上,該方法還可進(jìn)一步優(yōu)化:(1) 1輪計算(1輪指8個扇區(qū)依次完成各自的計算)的最后1個扇區(qū)和下一輪的第1個扇區(qū)也可通信合并;(2) 調(diào)整扇區(qū)的計算順序,進(jìn)一步降低通信次數(shù)。通過進(jìn)一步優(yōu)化,通信次數(shù)可降為128。
本文的通信優(yōu)化方法,采用轉(zhuǎn)發(fā)通信的通信策略,實現(xiàn)了更少的通信次數(shù)。該方法僅與面相鄰的鄰居計算單元直接通信,并通過面相鄰的進(jìn)程轉(zhuǎn)發(fā)來自角相鄰與邊相鄰的計算單元之間的通信。采用通信轉(zhuǎn)發(fā)算法,在1個扇區(qū)和7個鄰居計算單元通信時,僅需進(jìn)行3個面相鄰的計算單元上的3次通信,即可完成和7個鄰居進(jìn)程的通信。該方法單獨使用,可將通信次數(shù)從224降為96。進(jìn)一步地,將轉(zhuǎn)發(fā)通信策略配合[28]中的通信合并方法,可將一輪的通信次數(shù)降為56,極大地壓縮通信次數(shù)。本文所采用的通信優(yōu)化方法是性能和通信策略層面的優(yōu)化,因此該通信優(yōu)化方法不會對模擬精度造成影響。
在通信實現(xiàn)層面,并行KMC框架提供了統(tǒng)一的通信抽象層,KMC模型的實現(xiàn)上只需實現(xiàn)給定數(shù)據(jù)的打包和解包即可。在該通信抽象層中,主要實現(xiàn)了:1) 進(jìn)程與鄰居進(jìn)程通信,獲取其ghost區(qū)域數(shù)據(jù),并給鄰居進(jìn)程發(fā)送對應(yīng)的ghost區(qū)域數(shù)據(jù),以支持計算進(jìn)程邊界處缺陷的躍遷速率;2) 進(jìn)程與鄰居進(jìn)程通信,將運動到進(jìn)程模擬區(qū)域外的缺陷數(shù)據(jù)發(fā)送給鄰居進(jìn)程,并接收從鄰居進(jìn)程發(fā)送的運動到該進(jìn)程的缺陷數(shù)據(jù)。通過這兩個通信步驟,結(jié)合高效的通信優(yōu)化策略,可支持并行KMC模擬中各進(jìn)程間的數(shù)據(jù)交換與高效計算。
金屬材料的許多性能主要決定于缺陷的類型和濃度,缺陷會改變完美晶體的結(jié)構(gòu)。缺陷中最簡單的是點缺陷,包括間隙和空位。目前很多KMC程序只能考慮空位的作用,而缺少間隙的引入。因此,針對并行KMC框架,開發(fā)了面向材料輻照微觀演化的空位-間隙模型。
針對BCC的晶格結(jié)構(gòu),定義如下基本KMC反應(yīng)事件(圖6):1) 空位躍遷,空位可向8個1nn近鄰處的原子位置躍遷;2) 間隙躍遷,雙原子間隙dumbbell對中的1個原子向8個1nn近鄰中的位于同一平面上的4個晶格點躍遷;3) 空位-間隙復(fù)合,在空位或間隙躍遷事件發(fā)生后,若間隙和空位位于一定范圍內(nèi),會發(fā)生復(fù)合反應(yīng)。
a——空位可向1nn附近的8個鄰居晶格點躍遷;b——間隙可向1nn附近的且位于同一平面的4個晶格點躍遷圖6 KMC空位-間隙躍遷計算模型示意圖Fig.6 KMC model for vacancy-interstitial transition
對于躍遷速率的計算,本文通過Arrhenius方程進(jìn)行計算:
其中:v為嘗試頻率;Ef和Ei分別為躍遷后、躍遷前體系的能量;kB為玻爾茲曼常量;T為體系溫度;e0與空位近鄰處原子類型有關(guān)。為計算躍遷前后體系的能量,本文引入勢函數(shù)的計算方法,包括vicent_2007勢函數(shù)[10]和EAM勢函數(shù)[29]。
為高效地基于并行KMC框架進(jìn)行模型的實現(xiàn),首先需針對模型中的空位和間隙缺陷進(jìn)行表示(圖7),并提供事件列表、事件類型的表示。為此,在本模型的實現(xiàn)層面,從框架提供的晶格點索引出發(fā),通過晶格點的id和類型確定對應(yīng)的晶格點位置是否為缺陷或單原子。針對不同類型的缺陷,如空位或間隙,分別采用列表進(jìn)行存儲,并可通過晶格點id進(jìn)行索引。其中,列表中的每項包括了該缺陷向周圍各不同方向躍遷的速率,對于間隙缺陷,還包括其取向信息。這種表示方式建立起了晶格點與缺陷之間的一一映射,并將缺陷的躍遷速率和躍遷事件統(tǒng)一起來管理與更新,以及通過不同的缺陷列表實現(xiàn)了不同躍遷事件的分類。
圖7 MISA-AKMC中針對間隙和空位缺陷的表示Fig.7 Representation of interstitial and vacancy defects in MISA-AKMC
為實現(xiàn)KMC模型,需遵循并行KMC框架提供的接口,其中主要包括以下部分的實現(xiàn):1) 躍遷速率接口,依據(jù)給定區(qū)域更新該區(qū)域的缺陷對應(yīng)的躍遷速率;2) 通信接口,通信時的數(shù)據(jù)打包和解包接口。其中,第1個接口的實現(xiàn)為KMC模型的核心,通過勢函數(shù)計算給定區(qū)域內(nèi)所有缺陷的躍遷速率;第2個接口是通信的必要依賴,數(shù)據(jù)打包和解包串聯(lián)起了并行KMC框架內(nèi)的高效通信算法。由于模型中,晶格點類型和缺陷類型多樣,這對高效的通信接口實現(xiàn)提出了挑戰(zhàn)。本文采用元數(shù)據(jù)加字節(jié)數(shù)據(jù)的方式來進(jìn)行數(shù)據(jù)打包和解包。在發(fā)送的數(shù)據(jù)包的頭部包含通信的元信息,包括需通信的單原子數(shù)、空位數(shù)和間隙數(shù)。在元數(shù)據(jù)后依次添加需通信的單原子數(shù)據(jù)、空位數(shù)據(jù)(如空位id、躍遷速率等數(shù)據(jù))和間隙數(shù)據(jù)。通過這種數(shù)據(jù)通信方式,可實現(xiàn)將復(fù)雜的數(shù)據(jù)類型進(jìn)行打包和解包,依賴并行框架獲得高效的通信性能。
對MISA-AKMC的測試均在曙光超算平臺上進(jìn)行,僅使用了其中的CPU核,系統(tǒng)配置和環(huán)境列于表1。
表1 曙光超算平臺的相關(guān)參數(shù)Table 1 Parameter of Sugon supercomputer
為驗證開發(fā)的并行MISA-AKMC程序的正確性和有效性,選用表2中數(shù)據(jù)對Cu團簇析出過程進(jìn)行驗證。模擬采用MPI單進(jìn)程進(jìn)行,依據(jù)SL KMC算法將每個進(jìn)程劃分為8個扇區(qū)。
表2 Fe-Cu-V(鐵-銅-空位)體系模擬算例的輸入?yún)?shù)Table 2 Input parameter of simulation case for Fe-Cu-V (iron-copper-vacancy) system
在模擬過程中,統(tǒng)計不同時刻孤立的溶質(zhì)Cu原子濃度和數(shù)量,并與相同輸入條件下的其他文獻(xiàn)[30]中的KMC模擬結(jié)果進(jìn)行對比,如圖8所示,其中的藍(lán)色曲線為MISA-AKMC模擬的Cu析出過程的孤立溶質(zhì)Cu原子濃度的演化曲線,黑色曲線為文獻(xiàn)[30]中對應(yīng)的模擬結(jié)果。模擬結(jié)果表明,MISA-AKMC模擬結(jié)果的曲線和文獻(xiàn)[30]中的結(jié)果十分相近,兩曲線誤差在允許范圍內(nèi)。
圖8 Cu溶質(zhì)析出過程中的孤立溶質(zhì)Cu原子比例的演化曲線Fig.8 Evolution of number of isolated solute atoms as a function of time
將模擬過程中的各階段進(jìn)行可視化,示于圖9。從可視化結(jié)果能觀察到Cu原子聚集現(xiàn)象,驗證了模擬結(jié)果的正確性。
圖9 MISA-AKMC模擬的Cu溶質(zhì)析出過程可視化Fig.9 Visualization of Cu solute precipitation simulation
為驗證并行KMC框架和基于該框架開發(fā)的計算模型的性能,進(jìn)行強擴展性測試。在曙光平臺上,采用80×80×80的模擬盒子,其中合金比例為95.98at%Fe-1.34at%Cu。模擬的通信步數(shù)為10,通信時間閾值為3.0×10-7s,即總模擬的蒙特卡羅時間固定為3.0×10-6s。
由于MISA-AKMC中采用的事件選擇算法時間復(fù)雜度為O(N),N為MPI進(jìn)程中的缺陷數(shù)量。為維持強擴展性測試中的問題總負(fù)載不變的要求,需給不同規(guī)模的輸入算例設(shè)置不同的缺陷數(shù)。這是因為每個進(jìn)程上,執(zhí)行躍遷后的時間增量取決于該進(jìn)程上的總躍遷速率,而總躍遷速率與進(jìn)程中缺陷數(shù)量和類型等因素有關(guān)。假設(shè)每個缺陷對應(yīng)的事件的躍遷速率大致相等,那么容易得到,在保持固定負(fù)載的前提下,運行kp個MPI進(jìn)程時,體系中引入的總?cè)毕輸?shù)量應(yīng)為p進(jìn)程時的k倍。這里,對于進(jìn)程數(shù)分別為1、2、4、8、16、32時,初始的空位缺陷數(shù)量分別為128、256、512、1 024、2 048、4 096。當(dāng)MPI進(jìn)程數(shù)從1增加到32時,相同負(fù)載下的加速比約為8,且運行時間會持續(xù)減少,并行具有持續(xù)的加速效果,這說明,該算例可強擴展到32倍或更高倍數(shù)的進(jìn)程數(shù)量。從MISA-AKMC的并行效率結(jié)果(圖10)來看,當(dāng)MPI進(jìn)程數(shù)擴展到基準(zhǔn)代32倍時,并行效率從100%降到24%左右,軟件并行效率良好。
圖10 MISA-AKMC并行效率測試結(jié)果Fig.10 Result of parallel efficiency test of MISA-AKMC
面向材料輻照微觀結(jié)構(gòu)演化模擬計算需求,本文基于KMC方法開發(fā)了并行KMC框架和用于模擬核材料輻照損傷的并行程序MISA-AKMC,并介紹了其實現(xiàn)細(xì)節(jié),包括sub-lattice并行方法、通信策略以及程序優(yōu)化。通過將MISA-AKMC用于模擬Fe-Cu合金的富Cu團簇析出過程,驗證了程序的正確性。在此基礎(chǔ)上,進(jìn)行了程序的強可擴展性測試,從單核擴展到32倍核心,強擴展的并行效率保持在24%以上,性能測試表明MISA-AKMC具有良好的擴展性。
本文涉及的KMC實現(xiàn)代碼均在github上開源,可通過http:∥github.com/misa-kmc/misa-akmc獲取。針對MISA-AKMC未來的開發(fā),將進(jìn)一步完善勢函數(shù)計算,結(jié)合機器學(xué)習(xí)方法引入更高效和更精確的勢函數(shù);同時會考慮off-lattice機制,進(jìn)一步豐富原子尺度材料模擬社區(qū)。