數值堆材料輻照損傷并行動力學蒙特卡羅軟件MISA-AKMC的設計與實現

儲根深，何遠杰，白 鶴，陳丹丹，任 帥，賈麗霞，吳 石，賀新福，楊 文，胡長軍

(1.智能超算融合應用技術教育部工程研究中心，北京 100083；2.北京科技大學，北京 100083；3.中國原子能科學研究院，北京 102413)

材料輻照損傷模擬是數值堆軟件的重要內容，采用多尺度模擬方法是研究材料輻照損傷問題的一個重要手段[1-6]。其中，動力學蒙特卡羅(KMC)方法是多尺度模擬中最為重要的方法之一，被廣泛用于研究微觀尺度的材料演化行為。不同于同樣原子尺度的分子動力學(MD)方法，KMC方法將原子運動軌跡粗化為體系組態躍遷，使該方法可在保持原子級別精度下有效將模擬的時間跨度從原子振動的尺度提高到組態躍遷的尺度，即實現秒甚至年量級的材料輻照行為模擬[7-9]。

KMC方法目前已有不少研究[10-19]，包括原子動力學蒙特卡羅(AKMC)方法、實體動力學蒙特卡羅(OKMC)方法、事件蒙特卡羅(EKMC)方法等。其中，AKMC方法以原子為基本對象，考慮原子和缺陷之間的躍遷導致的微結構演化行為，其在空間尺度上更精確且方便與同樣以原子為操作對象的MD方法進行耦合。但在面向實際應用需求時，其仍面臨著內存限制和復雜的計算量等挑戰。隨著模擬體系的增大，模擬所需內存需求的矛盾也開始凸顯出來。同時，雖然KMC模擬在時間尺度上相對于MD提升了數個數量級，但計算過程仍十分耗時，特別是在大體系模擬時。解決計算和內存的雙重矛盾最有效的方式就是并行化，通過計算任務分配以并行計算，并通過分布式內存方案減少單個節點上的內存占用。基于此，本文將采用AKMC方法設計數值堆材料輻照損傷動力學并行模擬軟件MISA-AKMC。

1 KMC方法與KMC并行化研究

和通用的蒙特卡羅思想類似，KMC方法的核心思想是隨機抽取事件。若可確定體系當前狀態下所有可能發生的事件及其發生的概率，即可將以這些事件發生的概率作為權重隨機抽取事件執行以改變體系的狀態，以此推進體系向前演化。圖1示出了一種通用的KMC執行流程，其在時間步循環內，一次計算所有可能的事件的躍遷速率(也即事件概率)，隨后通過隨機數選擇其中的1個事件，并行執行該事件；如此進行時間步的循環，即可推進系統向前演化，直到模擬到預期的時間。在材料微觀演化模擬過程中，事件可以是缺陷的位置躍遷或新缺陷的生成等。

圖1 KMC執行流程Fig.1 KMC execution process

由于串行的KMC中，體系完成的連續兩次演化是獨立、無記憶的，這個過程是一種典型的馬爾可夫過程(Markov process)。且計算兩個相鄰的事件之間具有依賴關系，即事件的發生是串行的，位于當前狀態的體系只發生1個事件，該事件發生且執行完成后才會進入到下一個事件的發生過程。這種依賴性給KMC并行化帶來了很大阻礙。

目前大多數已有的并行KMC方法基本是基于區域分解的，將模擬區域劃分為多個子區域，并將其分配到不同的計算單元上。同時，這些方法近似認為這些子區域之間的距離足夠遠，之間的相互影響不大，因此不同子區域上的事件可同時發生。這就是現有的基于空間分解的并行KMC算法的基本思想。典型的并行KMC算法有基于空間分解方案的SL(sub-lattice)KMC[20-21]、SP(synchronous parallel)KMC[22-24]、SR(synchronous relaxation)KMC[25-26]等。其中，SP KMC算法和SR KMC算法在執行過程中，每個時間步均需引入昂貴的全局通信，影響計算性能；而SL KMC算法執行過程不必全局通信，這也成為本文的并行KMC算法的首選。

在SL KMC算法中，為實現正確的并行計算，需解決兩個問題：1) 邊界處的事件沖突；2) 不同計算單元上的時間同步(實際上，在其他兩種并行SP KMC和SR KMC算法中，也存在這兩個問題)。為解決邊界處的事件沖突，一般采用子區域著色方法。其將每個計算單元上的子區域再分為多個扇區，不同位置的扇區標記為不同顏色，且不同計算單元中處于同一位置的扇區標記為相同顏色。通過著色使相同顏色的扇區在空間上不相鄰。各計算單元均依次選擇相同顏色的扇區進行計算。KMC并行中的事件沖突問題如圖2a所示，當兩個計算單元上的事件發生在邊界處的鄰近位置時，可能會產生事件沖突。為解決這種沖突，可針對計算單元對應的子區域劃分扇區并著色。著色算法解決事件沖突如圖2b所示，圖中的一維著色示例，黃色與橙色的扇區分別被間隔開來，依次選擇黃色扇區和橙色扇區進行計算。當然，在實際程序中，考慮到模擬的三維空間和計算單元的三維邏輯分布，將每個計算單元內對應的子區域劃分為8個扇區，分別著色為S0、S1、…、S7，如圖2c所示。針對不同計算單元上的時間同步問題，SL并行KMC方法采用時間閾值的方式。給定時間閾值T，各計算單元上進行各自的KMC演化模擬，直到計算單元中當前扇區的時間增量達到閾值。在時間增量達到閾值后，即可與鄰居進程進行通信以交換數據(如同步躍遷對鄰居進程的更改，從鄰居進程處同步ghost區域等)。通信完成后，各計算單元切換到下一扇區，進行下一扇區的計算模擬，直到計算單元上的總時間增量達到閾值2T，如此往復。通過該時間閾值，計算單元在其時間增量到達nT時進行通信(其中n為正整數)，通信完成后進行下一扇區的模擬(同一計算單元上的各扇區的模擬依次進行)。

圖2 KMC并行中的事件沖突(a)、著色算法解決事件沖突(b)和三維空間上的扇區劃分與著色(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)

2 MISA-AKMC并行KMC框架

圖3示出了MISA-AKMC程序及并行KMC框架架構，為便于KMC模型的集成與可擴展，并方便軟件的核心優化，將MISA-AKMC并行部分分離開，形成KMC并行框架。該框架采用空間分解的并行思想實現了SL并行KMC方法，并提供額外的接口以供實現其他的并行KMC算法?？蚣茚槍诵牡牟⑿胁糠?，提供了多種加速優化方法，包括局部躍遷速率更新方法和轉發通信算法優化等。此外，框架還內置了晶格點和缺陷存儲等數據結構的表示，并向上層模型提供了KMC模型集成接口。

圖3 MISA-AKMC程序及并行KMC框架架構Fig.3 Architecture of MISA-AKMC program and parallel KMC framework

2.1 并行sub-lattice 方法實現

作為并行KMC框架的核心部分，本文實現的sub-lattice并行KMC算法計算流程如圖4所示。在本文的sub-lattice算法實現中，先依據通信閾值T計算出通信次數n，以此進行外層的時間步循環。在時間步循環內部，依次進行計算單元內的8個扇區的迭代：1) 扇區內的事件循環，對于每個扇區，重復在該扇區選擇并執行時間，直到該扇區的時間增量達到預設的閾值T；2) 通信，在一扇區s的計算完成后，還需進行通信以準備下一個扇區的計算，包括將扇區s對鄰居計算單元上的更改通信給鄰居計算單元，同時為計算下一扇區，還需從鄰居計算單元上通信獲取下一扇區的ghost區域；3) 進入下一扇區的計算。

圖4 MISA-AKMC并行KMC框架中的并行算法流程Fig.4 Parallel algorithm flow in parallel KMC framework of MISA-AKMC

2.2 并行框架優化

為進一步提升MISA-AKMC程序性能，本文還針對其中核心并行部分進行了優化。需強調的是，由于MISA-AKMC采用并行算法和計算模型分離的設計，針對并行框架部分的優化和改動不會影響到KMC模型部分的實現，這是本文設計的KMC程序的一大優勢。

1) 局部更新優化

為模擬的精確性，將通信閾值T取得較小，一般近似的等于單次事件產生的時間增量。這必然導致扇區內的事件循環中事件發生的數量很少，這些事件影響的區域也有限。故不必在每個事件發生后“全局地”更新該扇區的躍遷速率，而只需更新上一個事件影響的區域內的躍遷速率。計算躍遷速率是KMC模擬中最耗時的，通過局部躍遷速率更新的優化方法，可極大減少更新躍遷速率的計算開銷。

考慮到單次躍遷影響區域可能覆蓋到其他扇區，甚至其他計算單元對應的區域，如圖5所示。所以，事件發生后，不僅需更新當前扇區的總躍遷速率，還需更新受影響的其他扇區的躍遷速率和鄰居計算單元上對應扇區的躍遷速率。為此，將事件的影響區域劃分為3個部分。其中，第1個部分為影響區域和當前計算扇區的交集，記為S；第2部分為影響區域與當前計算單元上但除當前計算扇區以外的區域的交集，記為C；第3部分為影響區域與ghost區域的交集(影響區域不會超過ghost區域的邊界)，ghost區域是計算單元模擬區域外面的一層額外區域，但該區域是來自鄰居進程上的副本，將該區域記為G。在更新躍遷速率時，針對S區域，直接更新；針對C區域，找到該區域所屬的扇區，再更新該扇區的總躍遷速率；對于區域G，先更改ghost區域的狀態，后續通信時，將更改同步到鄰居進程，再在鄰居進程上重新計算這部分受影響區域的總躍遷速率。

圖5 KMC事件的可能影響區域Fig.5 Possible effect region of KMC event

2) 進程通信優化

在三維KMC模擬中，每個扇區與7個鄰居計算單元相鄰。在圖4所示的算法中，針對兩個通信過程，每個扇區均與7個鄰居計算單元進行通信。直接的通信方式，完成1輪8個扇區的計算，需進行224次通信(7鄰居×8扇區×2個通信操作×2(發送/接收操作))，如SPPARKS程序[27]。在Wu等[28]開發的Crystal-KMC程序中，提出了通信合并的方案，在1個扇區計算完成后的同步模擬區域的通信可與下一扇區計算前的同步ghost區域的通信進行合并，這樣可將通信次數降為182。

實際上，該方法還可進一步優化：(1) 1輪計算(1輪指8個扇區依次完成各自的計算)的最后1個扇區和下一輪的第1個扇區也可通信合并；(2) 調整扇區的計算順序，進一步降低通信次數。通過進一步優化，通信次數可降為128。

本文的通信優化方法，采用轉發通信的通信策略，實現了更少的通信次數。該方法僅與面相鄰的鄰居計算單元直接通信，并通過面相鄰的進程轉發來自角相鄰與邊相鄰的計算單元之間的通信。采用通信轉發算法，在1個扇區和7個鄰居計算單元通信時，僅需進行3個面相鄰的計算單元上的3次通信，即可完成和7個鄰居進程的通信。該方法單獨使用，可將通信次數從224降為96。進一步地，將轉發通信策略配合[28]中的通信合并方法，可將一輪的通信次數降為56，極大地壓縮通信次數。本文所采用的通信優化方法是性能和通信策略層面的優化，因此該通信優化方法不會對模擬精度造成影響。

在通信實現層面，并行KMC框架提供了統一的通信抽象層，KMC模型的實現上只需實現給定數據的打包和解包即可。在該通信抽象層中，主要實現了：1) 進程與鄰居進程通信，獲取其ghost區域數據，并給鄰居進程發送對應的ghost區域數據，以支持計算進程邊界處缺陷的躍遷速率；2) 進程與鄰居進程通信，將運動到進程模擬區域外的缺陷數據發送給鄰居進程，并接收從鄰居進程發送的運動到該進程的缺陷數據。通過這兩個通信步驟，結合高效的通信優化策略，可支持并行KMC模擬中各進程間的數據交換與高效計算。

3 并行框架下的材料演化模型

金屬材料的許多性能主要決定于缺陷的類型和濃度，缺陷會改變完美晶體的結構。缺陷中最簡單的是點缺陷，包括間隙和空位。目前很多KMC程序只能考慮空位的作用，而缺少間隙的引入。因此，針對并行KMC框架，開發了面向材料輻照微觀演化的空位-間隙模型。

針對BCC的晶格結構，定義如下基本KMC反應事件(圖6)：1) 空位躍遷，空位可向8個1nn近鄰處的原子位置躍遷；2) 間隙躍遷，雙原子間隙dumbbell對中的1個原子向8個1nn近鄰中的位于同一平面上的4個晶格點躍遷；3) 空位-間隙復合，在空位或間隙躍遷事件發生后，若間隙和空位位于一定范圍內，會發生復合反應。

a——空位可向1nn附近的8個鄰居晶格點躍遷；b——間隙可向1nn附近的且位于同一平面的4個晶格點躍遷圖6 KMC空位-間隙躍遷計算模型示意圖Fig.6 KMC model for vacancy-interstitial transition

對于躍遷速率的計算，本文通過Arrhenius方程進行計算：

其中：v為嘗試頻率；Ef和Ei分別為躍遷后、躍遷前體系的能量；kB為玻爾茲曼常量；T為體系溫度；e0與空位近鄰處原子類型有關。為計算躍遷前后體系的能量，本文引入勢函數的計算方法，包括vicent_2007勢函數[10]和EAM勢函數[29]。

3.1 模型表示

為高效地基于并行KMC框架進行模型的實現，首先需針對模型中的空位和間隙缺陷進行表示(圖7)，并提供事件列表、事件類型的表示。為此，在本模型的實現層面，從框架提供的晶格點索引出發，通過晶格點的id和類型確定對應的晶格點位置是否為缺陷或單原子。針對不同類型的缺陷，如空位或間隙，分別采用列表進行存儲，并可通過晶格點id進行索引。其中，列表中的每項包括了該缺陷向周圍各不同方向躍遷的速率，對于間隙缺陷，還包括其取向信息。這種表示方式建立起了晶格點與缺陷之間的一一映射，并將缺陷的躍遷速率和躍遷事件統一起來管理與更新，以及通過不同的缺陷列表實現了不同躍遷事件的分類。

圖7 MISA-AKMC中針對間隙和空位缺陷的表示Fig.7 Representation of interstitial and vacancy defects in MISA-AKMC

3.2 模型針對框架接口的實現

為實現KMC模型，需遵循并行KMC框架提供的接口，其中主要包括以下部分的實現：1) 躍遷速率接口，依據給定區域更新該區域的缺陷對應的躍遷速率；2) 通信接口，通信時的數據打包和解包接口。其中，第1個接口的實現為KMC模型的核心，通過勢函數計算給定區域內所有缺陷的躍遷速率；第2個接口是通信的必要依賴，數據打包和解包串聯起了并行KMC框架內的高效通信算法。由于模型中，晶格點類型和缺陷類型多樣，這對高效的通信接口實現提出了挑戰。本文采用元數據加字節數據的方式來進行數據打包和解包。在發送的數據包的頭部包含通信的元信息，包括需通信的單原子數、空位數和間隙數。在元數據后依次添加需通信的單原子數據、空位數據(如空位id、躍遷速率等數據)和間隙數據。通過這種數據通信方式，可實現將復雜的數據類型進行打包和解包，依賴并行框架獲得高效的通信性能。

4 MISA-AKMC測試與分析

4.1 測試環境

對MISA-AKMC的測試均在曙光超算平臺上進行，僅使用了其中的CPU核，系統配置和環境列于表1。

表1 曙光超算平臺的相關參數Table 1 Parameter of Sugon supercomputer

4.2 正確性驗證

為驗證開發的并行MISA-AKMC程序的正確性和有效性，選用表2中數據對Cu團簇析出過程進行驗證。模擬采用MPI單進程進行，依據SL KMC算法將每個進程劃分為8個扇區。

表2 Fe-Cu-V(鐵-銅-空位)體系模擬算例的輸入參數Table 2 Input parameter of simulation case for Fe-Cu-V (iron-copper-vacancy) system

在模擬過程中，統計不同時刻孤立的溶質Cu原子濃度和數量，并與相同輸入條件下的其他文獻[30]中的KMC模擬結果進行對比，如圖8所示，其中的藍色曲線為MISA-AKMC模擬的Cu析出過程的孤立溶質Cu原子濃度的演化曲線，黑色曲線為文獻[30]中對應的模擬結果。模擬結果表明，MISA-AKMC模擬結果的曲線和文獻[30]中的結果十分相近，兩曲線誤差在允許范圍內。

圖8 Cu溶質析出過程中的孤立溶質Cu原子比例的演化曲線Fig.8 Evolution of number of isolated solute atoms as a function of time

將模擬過程中的各階段進行可視化，示于圖9。從可視化結果能觀察到Cu原子聚集現象，驗證了模擬結果的正確性。

圖9 MISA-AKMC模擬的Cu溶質析出過程可視化Fig.9 Visualization of Cu solute precipitation simulation

4.3 軟件性能測試

為驗證并行KMC框架和基于該框架開發的計算模型的性能，進行強擴展性測試。在曙光平臺上，采用80×80×80的模擬盒子，其中合金比例為95.98at%Fe-1.34at%Cu。模擬的通信步數為10，通信時間閾值為3.0×10-7s，即總模擬的蒙特卡羅時間固定為3.0×10-6s。

由于MISA-AKMC中采用的事件選擇算法時間復雜度為O(N)，N為MPI進程中的缺陷數量。為維持強擴展性測試中的問題總負載不變的要求，需給不同規模的輸入算例設置不同的缺陷數。這是因為每個進程上，執行躍遷后的時間增量取決于該進程上的總躍遷速率，而總躍遷速率與進程中缺陷數量和類型等因素有關。假設每個缺陷對應的事件的躍遷速率大致相等，那么容易得到，在保持固定負載的前提下，運行kp個MPI進程時，體系中引入的總缺陷數量應為p進程時的k倍。這里，對于進程數分別為1、2、4、8、16、32時，初始的空位缺陷數量分別為128、256、512、1 024、2 048、4 096。當MPI進程數從1增加到32時，相同負載下的加速比約為8，且運行時間會持續減少，并行具有持續的加速效果，這說明，該算例可強擴展到32倍或更高倍數的進程數量。從MISA-AKMC的并行效率結果(圖10)來看，當MPI進程數擴展到基準代32倍時，并行效率從100%降到24%左右，軟件并行效率良好。

圖10 MISA-AKMC并行效率測試結果Fig.10 Result of parallel efficiency test of MISA-AKMC

5 結論和展望

面向材料輻照微觀結構演化模擬計算需求，本文基于KMC方法開發了并行KMC框架和用于模擬核材料輻照損傷的并行程序MISA-AKMC，并介紹了其實現細節，包括sub-lattice并行方法、通信策略以及程序優化。通過將MISA-AKMC用于模擬Fe-Cu合金的富Cu團簇析出過程，驗證了程序的正確性。在此基礎上，進行了程序的強可擴展性測試，從單核擴展到32倍核心，強擴展的并行效率保持在24%以上，性能測試表明MISA-AKMC具有良好的擴展性。

本文涉及的KMC實現代碼均在github上開源，可通過http:∥github.com/misa-kmc/misa-akmc獲取。針對MISA-AKMC未來的開發，將進一步完善勢函數計算，結合機器學習方法引入更高效和更精確的勢函數；同時會考慮off-lattice機制，進一步豐富原子尺度材料模擬社區。