ARM計算環境下堆芯程序的移植*

明平洲，李治剛，劉 婷，蘆 韡，劉 東，曾 輝，余紅星

(1.中國核動力研究設計院，四川 成都 610213；2.中國核動力研究設計院核反應堆系統設計技術重點實驗室，四川 成都 610213)

1 引言

堆芯計算涉及到的學科主要有反應堆物理、熱工水力、力學和燃料性能等，它們的單學科計算或者多學科耦合計算均用于論證和分析堆芯內部存在的大量異質現象[1]。由于真實堆芯的解析解極難獲得，采用不同的數值離散方法對各種物理數學方程進行工程求解成為計算相關研究人員的核心工作[2]。現階段中國核動力研究設計院引進的基礎集群計算設施采用異構設計，國產飛騰處理器組成的計算分區能夠通過SLURM集群資源管理系統進行計算作業的分區投遞。目前國內高端處理器市場主要有龍芯、飛騰和兆芯等品牌，其中天津飛騰信息有限公司研制的CPU已加入到開放的ARM生態系統，飛騰64位通用處理器產品兼容ARMv8指令集[3]，采用28 nm工藝流片，具備高性能和低功耗等特點，定位為中高端Intel至強服務器芯片的替代品。區別于Intel商用處理器使用的X86復雜指令集體系，ARM是典型的精簡指令集體系，這使得集群計算機有不同特點的計算資源可使用。不同指令集體系下的處理器存在著不同特點[4]，在商用處理器受限或計算資源緊張的情況下，ARM計算資源有可能成為一種補充的數值計算資源。

本文研究移植的4個軟件描述為：

(1)NACK-R是使用粗網節塊法求解穩態中子擴散方程的C/C++程序[5]，常用于充當堆芯燃料管理軟件的核心求解器，提供堆芯各個工況點下的通量分布和功率分布。

(2)中國核動力研究設計院堆芯熱工水力分析軟件CORTH(CORe Thermal-Hydraulics)是用于堆芯熱工水力分析的子通道程序，可用于分析研究反應堆及其回路系統中冷卻劑的流動特性、熱量傳輸特性和燃料元件的傳熱特性，且正在持續進行新功能開發和細節改進[6]。

(3)特征線中子輸運計算方法具有計算精度高、幾何適應性強和形式簡單等優點。美國麻省理工學院研制的開源特征線輸運程序OpenMOC(Open Method Of Characteristics neutral particle transport code)使用Python和C++進行混合編程，內部的數值計算核心采用C++編程語言，早期版本用于任意二維幾何區域的輸運計算[7]。

(4)全堆芯的中子物理參數求解在多個層面上存在著均勻化計算。KYLIN2程序現階段在組件級別為堆芯程序提供組件庫的計算，其核心為使用特征線方法的中子輸運求解器，具備多種參數的計算功能[8]。

2 計算環境說明

中國核動力研究設計院引進的集群計算機具備多個異構計算結點，使用資源管理軟件SLURM劃分為多個分布式內存計算區域，其中2類分區分別對應于商業處理器環境和國產飛騰處理器環境。前者在運行穩定性和配套軟件上較為成熟；后者面向通用計算，但計算效率和運行穩定性略為不同，這源于ARM處理器結構主要應用于嵌入式數字設備等領域，通用ARM芯片在數值分析計算領域尚處于局部應用狀態。

硬件描述：飛騰FT-1500A是基于ARM64位架構的服務器芯片，面向計算服務器市場。在本單位購置的集群計算機上，每個邏輯上的飛騰計算結點使用1個飛騰FT-1500A處理器，具備16個處理器核心，每個處理器核心運行單個線程。每個飛騰計算結點上的內存為32 GB，具備16 MB的緩存。根據供應商提供的文檔可以獲知該通用計算處理器的主頻為1.6 GHz。

軟件描述：圍繞麒麟操作系統建立系統環境，關聯集群計算機的各個軟件層次。其中開發環境主要為GCC 4.9.3，安裝Glibc 2.19、make 4.0和開發工具鏈Binutils 2.24.90[9]。

待移植程序或軟件的部分數值求解內容依賴于第三方數值函數庫，因此移植過程需要圍繞2類軟件進行移植、修訂和測試。一類是底層系統軟件，用于為上層的堆芯數值計算提供通用計算功能，在ARM計算環境中擺脫對商業處理器環境中的函數庫依賴；另一類是堆芯計算本身的求解內容，需要進行程序代碼改寫和修訂，從而以同一套程序代碼適應不同硬件的計算環境。表1統計了真實移植過程所涉及的軟件。

Table 1 Softwares to be ported

本文首先驗證飛騰處理器在ARM環境中進行移植的可行性和正確性，記錄上述軟件的移植過程和數值實驗。然后利用OpenMP對KYLIN2程序并行進行優化，論證單結點多個飛騰處理器核心的并行能力。這里選取的對比對象為Intel商用處理器，其中Intel商用處理器的主頻約為飛騰處理器主頻的2倍，理論上飛騰處理器的串行運行效率為Intel商用處理器串行運行效率的一半。分析KYLIN2程序并行優化時，按照最多使用12個處理器核心進行討論。

3 移植和測試

整個移植工作采用自底向上的流程以適應由飛騰處理器所構建的集群計算機單結點ARM計算環境。測試過程將對比Intel商用處理器所在計算環境和ARM計算環境的效率、數值穩定性和正確性。

3.1 系統軟件的移植

參與數據前后處理的DATA函數庫和HDF5讀寫庫直接使用ARM計算環境提供的C/C++編譯器進行編譯，然后使用軟件提供的基準程序進行測試，確保文本文件或二進制文件的正確讀寫。

參與數值計算的PETSc、SparseLib++、Lapack和SuperLU在移植過程中重點用于求解線性方程組和參與特征值計算等數值問題，這些函數庫由多個模塊組成，底層依賴于向量、矩陣和索引集等數據結構，內部采用C/C++語言和Fortran語言混合編程。本文將其移植到飛騰計算結點后使用GNU編譯器(gcc，g++和gfortran)以及基于GNU編譯器生成的并行編譯器(mpicc，mpicxx和mpifort)進行編譯，能夠正常生成靜態庫和動態庫。與此同時，基于已獲取的熱工分析相關線性方程組的研究工作[10]，本文使用移植的庫進行線性方程組求解例題的測試，以證明飛騰處理器構建的ARM通用計算環境能夠正常完成基本的數值計算，且收斂精度均可達到預期。

3.2 NACK-R

NACK-R是用于堆芯燃料管理分析的C/C++版本的粗網節塊法擴散計算程序(CORCA-3D擴散計算模塊的參考原型)，本文將其移植到ARM環境中，并參照NACK-F的測試內容，使用4個例題進行對比測試，計算結果主要統計有效增殖因子keff和例題的運行時間(各個粗網節塊的通量值也進行了對比，計算結果保持一致)。Intel商用處理器的計算環境和飛騰處理器的計算環境所使用的編譯器均選擇g++。

編譯的優化選項首先選取為O0，保證不引入編譯器的特定優化行為，統計結果如表2所示。

Table 2 Test Comparison of NACK-R using O0 optimization in Intel environment and Phytium environment

編譯的優化選項為O2，相應的統計結果如表3所示。

Table 3 Test comparison of NACK-R using O2 optimization in Intel environment and Phytium environment

編譯的優化選項為O3，相應的統計結果如表4所示。

Table 4 Test comparison of NACK-R using O3 optimization in Intel environment and Phytium environment

該程序移植使用C/C++編程語言，從統計數據來看，串行運行時飛騰處理器隨著編譯優化等級的變化能夠持續對效率進行改善，O3優化等級比起O2優化等級仍然有10%左右的性能提升。Intel商用處理器在O3優化等級和O2優化等級的計算效率基本相同。整體上Intel商用處理器性能高于飛騰計算結點的飛騰處理器，O0優化等級情況下數據量較小時兩者相差約7倍，在數據量較大時相差約9倍，高于頻率之間的2倍差異。隨著編譯優化等級的提升，飛騰處理器本身的運行效率有了較大提升，2種處理器之間的性能差異在縮小，這表明C/C++數值計算程序在使用精簡指令集的情況下可以通過編譯優化等級的改變來明顯提升計算效率，此時2種硬件之間性能差距也接近頻率差距的2倍。移植和測試結果也表明2種硬件環境中的數值計算結果保持一致，遵循IEEE的浮點數標準規范。

3.3 CORTH

CORTH程序采用矩陣形式求解具有滑速比的四方程均勻流模型，可以實現按照全堆芯柵元級別劃分子通道的數值計算，求解整個堆芯的焓升和壓降。本文移植的CORTH版本為0.8.4，其主要編程語言為C/C++。為了保證數值穩定性以及Intel處理器環境和飛騰處理器環境能夠針對同一套程序代碼進行編譯，向ARM環境的移植過程需要對CORTH 0.8.4版本完成以下2方面的代碼修訂工作：

(1)去除Intel 商業函數庫MKL(Math Kernel Library)的依賴。CORTH程序使用MKL函數庫內的mkl_dcsrgemv函數對能量守恒和動量守恒形成的線性方程組進行迭代求解，這里將求解函數替換為SparseLib++提供的GMRES函數。該函數同樣使用廣義最小殘差的Krylov子空間方法來迭代求解線性方程組，能夠獲得與原函數mkl_dcsrgemv相同的計算精度，且串行計算效率更高。

(2)編譯過程需要在ARM環境中重新編譯CORTH程序使用的第三方內容，例如2.1節列出的DATA類和熱工物性函數庫。該移植過程涉及到Fortran和C++混合編程，為了獲得不同硬件環境中程序代碼的一致性，需要按照GNU編譯器推薦的標準混合編程規范重寫計算函數的聲明，例如熱工物性函數庫中的func函數，重寫的頭文件聲明如下所示：

#ifdef__cplusplus

extern"C"

{

#endif

externdoublefunc_(double*,double*,double*);

#ifdef__cplusplus

}

#endif

它在Intel處理器環境和ARM處理器環境中均能夠被正確識別和混合編程應用。

當C/C++編程語言內調用Fortran語言形式的熱工物性函數接口時，需要在函數名稱func后添加下劃線“_”，且形參均為引用傳遞。

完成CORTH移植后,本文針對ACP1000的堆芯構造4種子通道劃分的例題進行串行計算的正確性和效率測試。2種硬件環境下均使用GNU編譯器，根據3.2節的移植分析結論，編譯優化等級統一設置為O2。

Table 5 Results and running time of CORTH

從表5的統計數據來看，引入SparseLib++的迭代法后，其數值穩定性和正確性在以飛騰結點代表的ARM環境和Intel處理器環境中保持一致，表明了飛騰處理器具備通用計算的特點。計算性能方面，Intel處理器明顯高于飛騰處理器，但不同的問題規模下飛騰處理器與Intel處理器的效率差距始終保持在3～4倍，略高于兩者之間的處理器頻率差異。這是由于ARM環境中飛騰處理器具有更小的緩存，間接影響CORTH程序在執行數值計算時的性能。

3.4 OpenMOC

本文將中子輸運的偏微分方程轉換為常微分方程，然后使用特征線形式進行數值求解。相比其它中子輸運理論數值計算方法，特征線方法具有計算精度高、幾何適應性強和形式簡單等優點。麻省理工學院研制的二維特征線輸運程序OpenMOC使用Python和C++進行混合編程，內部的數值計算核心采用C++語言編程，用于任意二維幾何區域的輸運計算。本文將OpenMOC 0.1.4b版本移植到ARM環境之后使用C5G7 1/4堆芯問題進行串行運行測試，以探究特征線方法的運行特點。在Intel處理器環境中使用相同程序對比驗證，記錄keff的數值和計算效率。各個網格柵元的通量經過可視化顯示后可以快速觀察2種環境中運行結果是否相同。圖1所示移植結果表明在Intel處理器與飛騰處理器上的計算保持一致。

Figure 1 Core layout and flux distribution of 2D C5G7 quarter core圖1 二維C5G7 1/4堆芯的布置和標通量分布

2種硬件環境下的編譯器均選擇GNU編譯器，編譯過程支持C++0x特性。固定C5G7例題的數值計算分辨率參數(能群數7，方位角個數16，特征線間距0.1 cm)，編譯優化選項首先取為O0等級，OpenMOC性能對比如表6所示。

Table 6 Test comparison of OpenMOC using O0 optimization in Intel environment and Phytium environment

編譯優化選項為O2等級，OpenMOC性能對比如表7所示。

Table 7 Test comparison of OpenMOC using O2 optimization in Intel environment and Phytium environment

編譯優化選項為O3等級，OpenMOC性能對比如表8所示。

Table 8 Test comparison of OpenMOC using O3 optimization in Intel environment and Phytium environment

上述測試完成之后可以明確，OpenMOC在飛騰結點和Intel處理器環境中的計算結果保持一致，但在每種編譯優化等級條件下飛騰結點的計算效率均低于Intel處理器的計算效率，比值相差約為3～4倍，高于兩者之間的頻率差異。隨著編譯優化等級的提升，2種計算環境中的效率差距有所縮小，不同于3.2節移植程序的特點，O2選項和O3選項在2種硬件環境下對于計算效率的改善相差不大，說明OpenMOC程序編程實現的特征線方法使用O3的優化等級并不會對生成的可執行程序的效率有更好的改善。

3.5 KYLIN2

KYLIN2程序使用二維特征線輸運算法集成燃耗計算用于生成組件截面庫，本文選取0.8.8版本KYLIN2移植到ARM環境中，源代碼按照多個目錄進行組織，每個目錄代表一個功能模塊，其整體計算步驟如下所示：

(1)計算參數的初始化和內存空間分配；

(2)讀取輸入文件；

(3)獲取計算基干的數據庫文件；

(4)判斷計算類型，執行組件計算或燃料分支計算；

(5)判斷燃耗步是否結束，結束則整個計算完成，否則轉回步驟(2)。

移植過程的主要工作分為2個方面，一方面是重寫編譯腳本來對這種多目錄組織的源代碼進行定制編譯和鏈接，規范化整個程序的軟件結構和執行順序；另一方面是移植和修訂程序使用的第三方內容。在飛騰計算環境中不能使用MKL函數庫，因此需要對KYLIN2的程序代碼進行調整，去除依賴性。為了不破壞已有的程序代碼，在移植修訂過程中通過宏開關來對改動代碼進行處理，宏的判斷標準為是否是ARM64硬件平臺。

#ifdef__aarch64__

(新增的修訂代碼)

#else

(原有代碼)

#endif

調用的MKL求解函數LAPACKE_dgesv用于求解Ax=B，其中A和B為矩陣。實質為使用LAPACKE中的求解函數，安裝Netlib提供的LAPACK和LAPACKE，替換MKL求解函數。調用MKL求解函數PARDISO，該函數使用直接消去法求解稀疏線性方程組Ax=b，其中A為矩陣，b為右手向量。使用SuperLU軟件(SuperLU軟件較為知名，是由勞倫斯伯克利國家實驗室研制的稀疏矩陣直接法求解函數庫)來應用直接消去法進行求解，替換MKL的求解函數。

移植完成之后，可在Intel結點和飛騰結點上使用同一套KYLIN2程序代碼。本文使用表9列舉的2個組件例題進行數值測試和驗證，統計相應的計算結果來論證其正確性和計算效率。

Table 9 Assembly testing examples for KYLIN2 program

正確性方面的對比按照燃耗深度與有效增殖因子(keff)的關系圖給出，計算效率則給出同一例題在不同計算環境中的總運行時間。2個例題均需要迭代多個燃耗步，有限增殖因子隨燃耗深度的變化趨勢如圖2所示，其中參考解選取的是原始KYLIN2程序使用MKL函數庫在Intel環境中的運行結果，根據3.4節的分析結果，統一選取的編譯優化選項為O2等級。圖2中橫坐標均為“燃耗步數”，縱坐標均為“有效增殖因子”，二者均無單位。

Figure 2 Burn-up change trend graph of two examples圖2 2個例題的燃耗變化趨勢圖

從圖2中可以看到，2種硬件環境中增殖因子隨燃耗深度的變化趨勢互相吻合，經過較長時間的迭代計算之后數值穩定性均未發生變化。2個例題在2種硬件環境中的串行運行時間如表10所示。

根據以上的移植和測試，可以得出KYLIN2程序運行這2個例題在飛騰計算環境與Intel計算環境中的效率差距在6～8倍，遠高于處理器頻率之間的差別。這一方面是由于組件計算過程中使用的數據量較大；另一方面是由于ARM精簡指令的運行方式在數值計算的細節上相比Intel處理器略有不足，KYLIN2程序可供優化的數值操作較多。盡管如此，KYLIN2程序底層的數值求解方法的替換對數值穩定性沒有影響，飛騰計算結點上數值計算的準確性與Intel計算環境沒有明顯差異。

Table 10 Serial running time of two examples

4 KYLIN2的OpenMP并行優化

實際工程應用需要反復調用KYLIN2，生成不同工況點下的大量組件截面庫，計算效率成為了應用的瓶頸。由飛騰處理器構建的集群計算機結點通常具備多個處理器核心，具備共享式內存并行的能力。對KYLIN2的輸運計算模塊編寫相應的OpenMP共享式內存并行版本KYLIN2-OpenMP，一方面在飛騰處理器為代表的ARM計算環境中提升KYLIN2程序的計算效率；另一方面也用于論證單個飛騰計算結點的并行性能，增強國產芯片的實用性。本文使用2種常見于學術論文的組件基準例題(壓水堆組件AFA3G和C5G7基準例題UO2組件)進行測試，給出對正確性和計算效率的對比結果。

先對KYLIN2程序不使用粗網加速算法的情況進行統計，根據3.4節和3.5節的討論，編譯優化選項選為O2等級，相應的計算結果和串行計算效率統計如表11所示。飛騰計算結點在求解AFA3G組件問題時耗時較長，與Intel環境中的耗時相差1個數量級(10倍)。

Table 11 Results and running time of KYLIN2 without GCMFD表11 無粗網加速的KYLIN2運行結果和時間統計

粗網加速算法將減少數值迭代次數，縮短串行運行時間。表12統計相應的運行時間，編譯優化選項仍然選為O2等級。

Table 12 Results and running time of KYLIN2 with GCMFD表12 粗網加速的KYLIN2運行結果和時間統計

特征線方法內數值方法的局部改變并不會改變程序內數據的內存布局，處理器之間的性能差異仍然保持在1個數量級左右?，F階段實際應用KYLIN2程序時往往在Intel處理器環境中使用單個處理器核心來運行例題，表12中數據表明飛騰處理器不能直接應用于堆芯程序的研發階段，因其有可能降低程序的運行效率?？紤]到KYLIN2程序在特征線非模塊化布置情況下，對組件問題的能群循環參數引入OpenMP共享式內存并行是論證ARM環境中使用飛騰處理器多個處理器核心的較好選擇。特征線方法內能群并行的OpenMP并行算法偽代碼列舉如下：

setcurrent←0;

setangular_flux←0;

#pragma omp parallel for default(shared) private(current) schedule(guided)

forg=1:GROUPS

forray=1:RAYS

forp=1:POLARS

current←FC(g,ray,p);

#pragma omp parallel reduction(+:sum)

angularflux←sweepFC(g,ray,p,current);

current←FB(g,ray,p);

#pragma omp parallel reduction(+:sum)

mesh_angular_flux←sweepFB(g,ray,p,current);

endfor

endfor

endfor

該算法中函數FC和FB不存在數據爭用的情況，遍歷特征線更新角中子通量密度的操作sweepFC和sweepFB需要使用OpenMP的互斥操作。KYLIN2-OpenMP多線程版本在2種硬件環境中進行了較多測試，限于篇幅這里僅給出C5G7 UO2組件例題不使用粗網加速算法的統計結果(AFA3G組件例題有類似結果)。實驗時考慮使用的最多處理器核心數量為12個，線程數取值參數取為1，2，4，8或12。2種硬件環境中的計算結果保持一致，這里重點分析計算所消耗的時間。

如圖3所示，2種硬件環境下C5G7組件例題的計算效率在O2優化等級下均通過多線程并行取得效率的提升(處理器核數超過8個之后并不能帶來計算效率的改善，這是由于C5G7例題的能群循環結構只有7個，超過數量7之后OpenMP并行增益達到上限)。使用單個線程運行時飛騰計算結點耗時11 123 s，在Intel環境中耗時1 205 s。使用8個線程時飛騰計算結點運行時間為1 598 s，Intel環境下的運行時間為120 s。整體趨勢上飛騰計算結點的計算效率與在Intel環境中的計算效率差距在縮小，且使用8個線程后飛騰處理器運行C5G7組件例題的效率與參照對象Intel處理器的單核運行時間接近。當集群計算機的商用計算資源緊張時，由于國產計算結點使用人數較少，此多線程并行版本的臨時解決方案可達成使用ARM計算環境臨時替換商用計算環境來執行工程計算任務的需求(現階段商用處理器計算結點經常使用單線程來運行KYLIN2程序)。

Figure 3 Running time of assembly example after OpenMP programming圖3 引入OpenMP并行后的組件例題運行時間

5 結束語

在內部集群計算機的飛騰計算結點上進行堆芯程序的移植和測試，探討國產芯片構建的計算環境應用于工程計算的可行性，形成的結論為：

(1)所有的程序均移植成功，由飛騰處理器構建的ARM計算環境具備通用計算的特點，適用于堆芯數值計算程序的研制和運行。移植的程序在飛騰計算結點上的串行計算效率均低于在參照對象Intel處理器上的串行計算效率，且不管在何種編譯優化等級條件下，計算時間的比值均超過處理器頻率之間的2倍差異，小規模堆芯例題的差異保持在3～4倍，較大規模堆芯例題的串行計算效率差異更大。

(2)移植過程的堆芯程序能夠在不同硬件環境中獲得一致的計算結果(偏差范圍內)，確保了數值計算的穩定性。

(3)飛騰計算結點上的數值計算運行效率較低，一方面是由于底層函數庫為了適應精簡指令集特性和通用計算的特性，需要執行更多的指令來完成某些數值計算；另一方面也由于飛騰處理器本身與參照對象Intel處理器之間的硬件架構差異。Intel商用處理器的復雜指令集存在多種復合指令，完成部分計算需要更少的指令周期，但通過對KYLIN2程序進行局部的OpenMP多線程優化，可以使得國產飛騰處理器對應的ARM計算環境能夠臨時替換商用處理器的串行運行環境，成為工程計算的補充資源。

現階段各個樣本程序在集群計算機的飛騰計算結點上逐步開展移植和運行，以期充分利用這部分計算資源。國產處理器及其生態環境持續發展拓寬了核反應堆工程計算與計算機科學的融合領域。橫向借鑒日本、美國等經驗，核反應堆工程計算在后續發展過程可以覆蓋更多的計算硬件，提升適用性和交叉驗證，更快和更準確地為反應堆設計和應用進行服務。本文形成的結論和記錄的技術細節將隨著后續在國產計算結點上開展程序研制和移植的工作進行深化和擴展。

致 謝

本研究工作感謝中國核動力研究設計院設計所九室提供的樣本程序和集群運行環境，同時也感謝相關人員對論文的建議。