直接三維中子輸運特征線法并行計算軟件的實現

汪 岸，胡長軍，王 根，曹 敏，單浩棟，胡 赟，楊 文，劉天才

(1.北京科技大學，北京 100083；2.中國原子能科學研究院，北京 102413)

數值反應堆是在超級計算機上實現的核反應堆各種物理過程及其耦合高保真數值模擬、預測的軟件系統。它是實際反應堆“外在”和“內在”的鏡像，是反應堆設計、建安、運行、退役、全周期從微觀機理到宏觀現象的研究平臺。近二十年來，歐美發達國家均啟動了數值反應堆的相關研究項目[1]，在各國的超級計算機上完成了大量對現役核反應堆的模擬實驗。為促進我國核能的安全高效發展，以國產超算為依托實現核反應堆的高精細模擬，國家重點研發計劃項目“數值反應堆原型系統開發及示范應用”的研究團隊設計并開發了面向E級超算的數值反應堆原型系統CVR1.0[2]。

中子輸運計算是CVR1.0的核心功能之一，它為數值堆中的其他模擬過程提供耦合計算所需的精細中子通量密度分布。在中子輸運的數值方法中，特征線法(MOC)[3]以其優秀的幾何適應性、較高的計算精度、簡潔易實現的計算模型及強大的并行潛力[4]逐漸成為研究熱點。三維特征線法同樣帶來了對計算資源的龐大需求，其中大量的幾何網格、特征線和通量數據給存儲帶來壓力，相應的海量源迭代浮點運算次數使計算相當耗時。隨著近年來計算機硬件水平的提升和高性能計算技術的發展[5-6]，大規模并行特征線法已具備了實現的基礎，國內外也陸續出現可大規模并行的三維特征線法程序，如OpenMOC[7]、NECP-X[8]、APOLLO3[9]等。這些程序在各自的應用領域內均取得顯著的成果，但是或使用結構化建模只能描述單一堆型、未能充分利用特征線的幾何適應性，或使用2D/1D耦合的間接三維特征線算法[10]，目前面向多種堆芯的直接三維特征線并行程序仍在發展階段。

為進一步開拓特征線方法的幾何適應性，借助反應堆模擬與高性能計算的學科交叉與專業互補，探究直接三維特征線方法的應用潛力，CVR1.0項目團隊基于流行的OpenMOC方法[11]，開發了直接三維中子輸運并行計算程序ANT-MOC，以期實現多堆型pin-by-pin精細計算。

本文建立三維特征線法中子輸運計算程序ANT-MOC的架構，給出ANT-MOC功能和性能兩方面的測試結果，驗證ANT-MOC中關鍵實現技術的有效性。

1 理論基礎

1.1 直接三維特征線法

中子輸運方程(玻爾茲曼方程)是一個用于研究介質內的中子密度分布的守恒方程，表示了因輸運、吸收、散射及裂變反應產生的中子在反應過程中的守恒關系。ANT-MOC能夠求解穩態輸運方程：

ΣT(r,E)Ψ(r,Ω,E)=Q(r,Ω,E)

(1)

式中：Ω為角方向向量；Ψ為角中子通量；r為空間位置向量；E為中子能量；Q為角中子源；ΣT為中子總截面。

輸運方程經過空間、角度、能群等多種離散化之后，其求解轉化為沿不同軌跡(空間離散后假想的中子運動軌跡)的特征線微分方程：

(2)

(3)

式中，j為FSR的編號，FSR總數為J。

將式(3)沿特征線進行積分求平均，得到該條特征線上的平均中子角通量：

(4)

式中，sk為特征線在該網格下的截斷長度。將網格內沿所有方向軌跡的平均中子角通量加權求和，即可解出該網格下的中子通量密度。沿著特征線方向循序求解，便得到中子通量密度在全計算空間下的分布情況。

以上是一次典型的特征線輸運計算流程。特征線計算需要提前獲取特征線在空間網格下的截斷長度，因而射線追蹤成為特征線法特有的幾何預處理步驟。不同的射線追蹤方式也決定了后續特征線輸運計算的具體流程。由式(2)～(4)可看出，特征線輸運計算以單條特征線為最小計算單元，計算過程中不同軌跡上的特征方程沒有強耦合關系，所以可同時進行多條特征線的射線追蹤計算，因而使得特征線方法十分有利于并行計算。

1.2 幾何建模方法

構造實體幾何(CSG)方法是一種直觀便捷的計算機建模方法，在主流蒙特卡羅方法[12]和部分特征線法[11]計算程序中得到廣泛應用。CSG建模的基本思路是，對若干由曲面圍成的不可分幾何實體(圓柱、球體等)進行幾何變換(平移、旋轉等)和布爾運算，以生成復雜幾何。ANT-MOC的幾何建模支持棱柱面、多種陣列分布形式，并且實現了交、并、補等布爾運算，能夠很好地建模壓水堆、快堆等重復結構較多的堆芯幾何。

圖1示出一假想的多柵元單組件結構，其可由組件邊緣鈉和柵元陣列兩部分合并得到。組件邊緣部分可視為外套管內部(正六邊形)和柵元陣列內部(棱柱面內部)的差集；柵元陣列可視為單個柵元結構在六邊形陣列上的重復排列。該建模過程中只有底層的基本幾何實體作為數據存儲，實體間的運算只在射線追蹤時才會進行，由此可在略微增加計算量的情況下，大大減少后續計算所需的存儲量。

圖1 ANT-MOC中的單組件CSG建模Fig.1 Assembly modeling with CSG implementation of ANT-MOC

1.3 基于軌跡鏈的射線追蹤

ANT-MOC支持3種幾何邊界：周期、真空、全反射。程序中生成的二維軌跡滿足循環軌跡的要求，即各軌跡可經由平移和對稱變換最終在幾何邊界首尾相連，構成軌跡鏈(周期性軌跡和反射性軌跡)。如圖2所示，在邊長為a×b的幾何上分布了5條方位角為φ的軌跡，這些軌跡通過簡單的平移變換可在邊界相連。六邊形幾何的情形下，平移變換沿y軸和α軸(x軸逆時針旋轉30°)進行。

圖2 四邊形幾何上的周期性軌跡示意圖Fig.2 Periodic track scheme on rectangular geometry

基于3DGT方法[13]和OTF追蹤方法[14]，ANT-MOC實現了針對四邊形和六邊形幾何的三維射線追蹤算法，該算法以軌跡鏈為最小單位生成源迭代所需的軌跡。根據循環軌跡的理論[15]，軌跡鏈的角度φ、長度L和間距δ⊥等關鍵參數可由公式給出。

周期性軌跡鏈在矩形幾何上的關鍵公式為：

(5)

式中：a、b分別為矩形的長和寬；n和m分別為沿x軸和y軸平移變換的次數(正負表示方向)。

周期性軌跡鏈在六邊形幾何上的關鍵公式為：

(6)

式中：a為六邊形幾何邊長；n和m分別是沿α軸和y軸平移變換的次數(正負表示方向)。

若將同一方位角下的周期性軌跡數量記為k，則可得到該角度下任意兩段相鄰軌跡的間距：

δ⊥(k)=δ⊥/k

(7)

基于軌跡鏈的算法在處理邊界通量時較為簡單，對于壓水堆和快堆兩種幾何配置能夠統一處理。一方面，在同一軌跡鏈中邊界通量總是沿計算方向傳遞；另一方面，不同軌跡鏈之間不存在邊界通量交換，這為后續的并行算法實現創造了條件。

2 程序實現

2.1 功能模塊

ANT-MOC在計算流程上可分為輸入輸出、幾何建模、射線追蹤、迭代求解4大部分。整個程序使用C++編寫，運用面向對象設計的思想分解各功能模塊，形成了易于擴展和維護的代碼結構。

圖3示出ANT-MOC的關鍵功能模塊和輔助工具。ANT-MOC主要使用MPI+OpenMP混合編程實現針對軌跡鏈、軌跡、能群的多級并行算法。基于OpenMP編程接口，預處理、軌跡生成、通量計算等多處計算密集的代碼，均采用多線程和SIMD向量化加速計算，提高處理器利用率。在具有多節點的集群上，MPI非阻塞通信接口用于實現通量數據的跨節點交換。在核心計算功能中，CSG建模和射線追蹤算法能夠適用于壓水堆和快堆幾何，是輸運掃描和軌跡鏈并行算法的基礎。ANT-MOC支持多種數據格式的輸入和輸出，在運行時能對幾何輸入和宏觀截面數據進行必要檢查，在運行后可輸出通量和軌跡的可視化數據，有利于對結果的觀察分析。除基本功能外，為保證軟件質量、便于代碼迭代更新，開發團隊吸取數值反應堆的開發經驗，建立了完善的代碼版本控制和自動化測試。為便于用戶分析輸出結果，ANT-MOC主體程序配套了兩個輔助工具：用于分析日志、生成表格的日志分析工具，以及用于通量數據文件格式轉換、截面數據生成、截面數據操縱的數據處理工具。

圖3 ANT-MOC架構Fig.3 ANT-MOC architecture

2.2 面向快堆的預處理

為計算軌跡在堆芯幾何穿行時被網格截斷后的線段長度，在生成三維軌跡時，需要獲得軌跡與幾何內所有網格的交點，ANT-MOC利用兩個核心算法完成：網格尋找算法和最近交點算法。網格尋找算法用于沿CSG樹形結構搜索給定坐標點所在的網格，最近交點算法用于計算軌跡從某點出發與幾何網格的最近交點。對于每一條軌跡，已知其起點、終點坐標和CSG幾何數據結構時，從起點出發反復執行網格尋找算法和最近交點算法，即可依次得到該軌跡與整個幾何體的所有交點。

在程序實現中，針對壓水堆幾何的算法易于實現，只要建立三維坐標點在CSG幾何不同深度(如陣列和陣列中的格點)的坐標變換，以及軌跡與不同幾何體的交點計算公式[16]。引入快堆相關的幾何體(如六邊形陣列、棱柱面)后，為表達圖1所示的快堆組件排布、組件與柵元陣列之間的鈉填充等結構，使軌跡生成、交點計算能夠正常執行，需要改進算法。

六邊形陣列和棱柱面(陣列的包圍面)兩種幾何對象可用于構建快堆幾何，其中，六邊形陣列的實現采用了兩組經典的斜坐標軸實現相關算法[12]。以圖4所示的幾何為例，第1組坐標軸是β-γ-δ軸，第2組為β-δ軸，沿各軸的單位向量分別為：

圖4 六邊形幾何預處理示意Fig.4 Illustration of hexagonal geometry preprocessing

δ=(0,1)

(8)

β-γ-δ軸恰好分別穿過正六邊形幾何的3對對邊，利用單位向量、軌跡方向向量等向量的內積運算可迅速求出軌跡到最近交點的距離。β-δ軸則相當于原x-y軸進行剪切變換后得到的新坐標軸，利用簡單的基變換公式，容易知道任何一個三維坐標點所在的陣列格點。

上述最近交點算法能夠很好地處理坐標點在陣列內的情況，但無法直接用于坐標點在陣列外的情況。圖4演示了軌跡從陣列和外套管之間的鈉區域進入陣列的情形。假設內部柵元陣列的柵距為wr，柵元排列為nr圈(即沿β方向有nr個柵元)，為了計算關鍵點，構造1個外接圓，其半徑router為：

(9)

同時構造1個包含該外接圓的外接陣列，陣列的柵距wouter和柵元排列圈數nouter分別為：

wouter=wr

(10)

借助上述兩個輔助結構，軌跡與柵元陣列邊緣的交點轉化為軌跡與外接陣列內部格點邊界的交點，此交點的計算可使用前述最近交點算法。此外，當軌跡起點遠離柵元陣列邊緣時，為了減少計算量，需計算軌跡與輔助外接圓的交點，將起點移置該交點后再執行最近交點算法。

2.3 軌跡鏈并行算法

軌跡鏈并行算法的核心，是利用軌跡鏈的分布規律實現相應的負載平衡策略以減少進程等待時間。如圖5所示，軌跡鏈被分配到各MPI進程，每條軌跡鏈中的軌跡由OpenMP線程并行處理，單條軌跡涉及的所有能群的計算由SIMD加速。

圖5 ANT-MOC的多級并行策略Fig.5 Multi-level parallel strategy of ANT-MOC

軌跡鏈并行算法的基本思想是從二維軌跡鏈的數量、長度規律出發，以“均勻分配長度”為目標建立1個組合優化問題，求解這個問題得到每個進程應有的二維軌跡鏈數量、編號等信息，再結合3DGT方法生成三維軌跡。假定以I={1，2，…，a}表示a個方位角的索引集，P={1，2，…，p}表示p個進程的索引集，xi,r表示第r個進程應該在第i個方位角下分配的軌跡鏈數量。用軌跡的長度之和作為計算負載的估計值，得到負載的表達式為：

(11)

式中：Li和ni分別為第i個方位角下軌跡鏈的長度和數量，根據循環軌跡原理，同一方位角的軌跡鏈有相同長度；W為總負載；Wr為第r個進程的負載。根據負載的表達式，以均勻分配軌跡總長度為目標構造混合二次規劃問題：

(12)

ANT-MOC采用動態規劃和貪心算法聯合求解上述問題，得到軌跡鏈的具體分配方案。特別地，若將同一方位角下的軌跡鏈長度求和、數量設為1，得到的便是軌跡以方位角為單位的分配方案。

2.4 數據輸入輸出

ANT-MOC的輸入數據分為3部分：幾何模型(XML格式)、多群截面(HDF5格式)和運行時參數(TOML或XML格式)。其中，運行時參數包括幾何輸入文件路徑、多群截面文件路徑以及控制程序行為的控制參數。輸出數據也分為3部分：程序日志(純文本格式)、原始通量數據(HDF5格式)和可視化數據(VTK格式)。原始通量數據記錄了每個FSR中的反應率、中子通量分布，可視化數據則是用戶提取的感興趣的通量分布數據。

程序輸出的反應率、通量和軌跡分布數據文件可由ANT-MOC數據處理工具轉換為VTK可視化格式，圖6示出多組件輸運計算的裂變率和軌跡分布經ParaView[17]可視化后的結果。

a——裂變率分布；b——三維軌跡分布圖6 127組件的ParaView可視化示意圖Fig.6 127 assemblies visualized scheme by ParaView

3 驗證與確認

ANT-MOC的驗證工作遵循科學計算軟件驗證與確認(V&V)的一般方法，包括計算模型的驗證、實驗結果對比、參數敏感性分析、不確定性量化等多項任務。國內外關于計算流體力學、結構力學等學科計算程序的V&V技術已有系統研究，關于輸運計算程序的V&V工作在CASL、NEAMS等美國數值反應堆項目中也有比較成熟的經驗。在文獻[18]中，作者使用Takeda、C5G7等壓水堆國際基準題對ANT-MOC的輸運計算能力做了初步驗證，表明ANT-MOC在小型柵元均勻化和pin-by-pin問題中能得到接近參考值的結果。

本文基于國際基準題Takeda和中國原子能科學研究院制作的中國實驗快堆(CEFR)算例，首先驗證ANT-MOC并行程序的串行等價性，隨后使用一般的敏感性分析方法研究ANT-MOC的參數范圍，最后給出CEFR全堆計算的并行性能測試結果。

3.1 串行等價性測試

數值模擬并行程序由于受到多種不確定性和誤差的影響，在相同參數條件下的串行和并行計算結果可能不同。本文基于Takeda Model 1例題完成4類測試：單進程單線程(串行計算)、單進程多線程、多進程單線程、多進程多線程。所有測試均在同一集群上完成，每計算節點包含2個Intel Xeon E5-2620 v3處理器，節點內共享64 GB內存，節點間使用千兆以太網互連。測試結果中的keff和中子通量分布作為判斷計算結果相等與否的依據。比較中子通量分布時，使用了均方根誤差(RMSE)：

(13)

式中：i為FSR編號；n為FSR總數量；φi為串行測試第i個FSR的標通量；φ′i為某一測試結果中第i個FSR的標通量。

測試參數中，FSR取邊長為1.25 cm的立方體，方位角數量取12，徑向軌跡間距取0.1 cm，極角數量取6，軸向軌跡間距取0.1 cm。線程數量與OpenMP相關，取值分別為1、12、24；進程數量與MPI相關，取值分別為1、4、6、12、24。

表1列出4類測試的結果。其中，每一類并行測試都包含兩個進程數或線程數不同的測試結果，7組測試在收斂時有相同的迭代次數，keff誤差在10-6以內，中子通量的均方根誤差在10-8以內。結果表明，ANT-MOC并行計算結果與原串行算法的計算結果幾乎相同，細微差別可能是由并行算法中各浮點運算順序不同導致的。

表1 串行、并行計算比較Table 1 Comparison between serial and parallel computing

3.2 參數敏感性分析

ANT-MOC輸入中的關鍵控制參數有兩類，一是控制FSR數量的參數，二是控制軌跡數量的參數。本文基于CEFR單個燃料組件幾何配置和514群宏觀截面輸入，選取5個對特征線法較為重要的控制參數依次測試。對于每個參數，分析參數值變化對keff的影響，隨后固定該參數的值并進行下一參數的測試與分析。如圖7所示，整個CEFR燃料組件幾何從上到下分為12層不同材料。在柵元陣列部分，柵元排列為5圈，相鄰柵元中心距離為0.695 cm。

圖7 CEFR燃料組件幾何示意圖Fig.7 Geometry scheme of CEFR fuel assembly

表2列出敏感性分析選取的輸入參數。由表2可見：扇區劃分數影響FSR數量，相應的測試數量較少；其余參數影響軌跡數量，相應的測試數量較多。除表2中列出的參數外，程序運行時使用的進程數量和線程數量也可能不相同，根據串行等價性測試結論，認為計算資源的數量對計算結果的影響可以忽略。

參數敏感性分析如圖8所示，以表2所列的參數編號順序進行，首先分析FSR扇區劃分數，此時固定軸向軌跡間距為0.05 cm，徑向軌跡間距為0.05 cm，方位角數量為16，極角數量為6。由圖8a可見，扇區劃分數增多時，keff波動幅度變小。為節省計算資源，后續參數的分析都給定扇區劃分數為0。

表2 敏感性分析選取的輸入參數Table 2 Input parameter for sensitivity analysis

圖8b～e表明，在影響軌跡數量的4個關鍵參數逐漸精細時，keff的波動幅度也逐漸變小，有收斂的趨勢。根據該結果，可以為同類型CEFR幾何配置(如多組件、全堆計算)給出1組使ANT-MOC計算結果較為可靠的推薦參數范圍：扇區劃分數不小于4，軸向軌跡間距不超過0.5 cm，徑向軌跡間距不超過0.05 cm，方位角數量不小于64，極角數量不小于12。

圖8 參數敏感性分析Fig.8 Parameter sensitivity analysis

3.3 并行性能測試

由于程序中的每條軌跡與幾何相交后被截斷為數條線段，ANT-MOC單次迭代的計算實際上與與總線段數量呈正相關，在固定負載的條件下增加計算資源完成的測試就是強可擴展性測試。測試均在中國科學院網絡信息中心的國產集群上完成，該集群每節點包含1個國產32核x86處理器，擁有8個DDR4內存通道，節點間使用200 Gbps的Infiniband HDR胖樹網絡互聯。

根據文獻[19]中給出的CEFR首爐堆芯布置編寫幾何輸入卡，以6群截面為輸入，方位角數量取64，徑向軌跡間距取0.005 cm，極角數量取4，軸向軌跡間距取0.25 cm。該條件下，各測試均在55次迭代后收斂。以4 096處理器核為基準擴展到98 304核，共使用6種不同核數規模，每種規模測試3次，記錄總時間、求解時間、射線追蹤時間、通信時間并取平均值作為最終結果。相應的并行效率計算公式[20]為：

(14)

式中：T1為程序在4 096核上執行的時間；N為核數相對于4 096的倍數，使用98 304核時，N=24；TN為程序在4 096×N核上執行的時間。

如圖9所示，以4 096核為基準，ANT-MOC在98 304核上的并行效率較高，達到52.60%，其中通信時間占比從2.54%增長到17.50%。圖10計算了射線追蹤部分的并行效率，在98 304核上僅為24.39%，且效率隨核數增加迅速下降。結合ANT-MOC的并行算法來分析，迭代計算耗時受軌跡和線段數量影響較大；通信時間主要消耗在MPI全局歸約，因此隨核數增加沒有減小；源迭代之前的射線追蹤過程以軌跡鏈為單位，其開銷相比于簡單的追蹤算法有所增加。綜上，性能測試結果表明，ANT-MOC并行算法有較好的強可擴展性，但通信算法和并行射線追蹤仍有較大改進余地。

圖9 ANT-MOC并行效率Fig.9 Parallel efficiency of ANT-MOC

圖10 射線追蹤并行效率Fig.10 Parallel efficiency of ray tracing

4 結論

本文闡述了三維特征線法程序ANT-MOC的理論基礎、程序實現與相關測試結果。CSG建模算法中對圓柱、棱柱等多種幾何體的實現使快堆堆芯建模更加方便、準確，幾何體的交、并、補運算實現同時簡化了用戶輸入和編程復雜度。利用循環軌跡原理，以軌跡鏈為基礎的射線追蹤算法能適用于壓水堆和快堆兩種堆芯幾何。在此基礎上建立的MPI+OpenMP多級并行算法，能夠利用軌跡鏈的分布規律實現負載平衡，在一定程度上提高并行效率。串行等價性測試表明，ANT-MOC的并行算法能夠得到與串行特征線法基本相同的計算結果。CEFR參數敏感性分析初步為快堆高精細模擬確定了合理的參數范圍。在以4 096核為基準擴展到98 304核的全堆計算性能測試中，軌跡鏈并行算法的強可擴展性超過50%，證明ANT-MOC具備大規模并行計算的能力。

綜上，ANT-MOC作為面向數值反應堆的直接三維中子輸運并行計算程序，具備對多種堆型高精細真實幾何下的大規模并行計算能力。后續研究中應以更多全堆例題為對象，通過對ANT-MOC參數敏感性的研究，確定合理的參數選取原則，并繼續優化程序的計算模型和代碼結構，進一步完善程序的適用性和穩定性。