堆芯計算的框架方法及其原型軟件NAC4R的設計1)

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

(中國核動力研究設計院，四川 成都 610213)

隨著并行計算技術發展以及電子計算機硬件性能的提升，反應堆工程領域具備研制較多大規模并行軟件的條件。并行計算聯合反應堆計算程序在高性能計算平臺上的高效應用成為科研設計和工程應用的研究熱點之一，它們從多個角度來提升核反應堆工程計算的效率和應用更多精細的模型[1,2]。伯克利并行實驗室提出的Dwarf導則將并行軟件的編寫按照行業進行抽象，總結出并行模式設計，并在研究中針對蒙特卡羅中子輸運計算形成了底層的框架軟件[3]。核能領域的CASL等項目也將并行工具和底層算法進行整合，形成適用于大規模集群計算機的軟件框架，并應用于工程耦合軟件的研發[4,5]。框架方法的實質便是基于復用構件和軟件結構進行軟件的開發。框架方法的共性研究在著名的FLASH研究文獻[6]中早期被總結為兩種技術路線：

1)先設計和開發框架，然后逐步添加求解器和其他計算能力；

2)圍繞小的求解器或專業應用進行持續研發，最終發展成為大型專業軟件。

第一種技術路線方面，桑迪亞國家實驗室研發的Trilinos是用于大規模工程和科學問題求解的軟件框架。該計算框架集成了較多知名軟件包，提供了豐富的數值計算功能;PETSc由阿貢國家實驗室維護，提供多個適用于并行計算的線性和非線性數值求解器和功能模塊;日本通過建立大規模并行數值計算技術研究小組，在京超級計算機上研制大規模的高性能數值軟件KMATHLIB，與底層AICS硬件緊密結合[7];Charm++是由Illinois大學基于C++編程語言研制的并行面向對象編程語言，它將集群計算機的處理器核心虛擬化[8]。第二種技術路線方面，PPM圍繞著物理學科的數值計算特點，研制了位于底層并行計算環境與學科應用之間的抽象中間件，減輕并行軟件的研制難度[9]。以上基于框架方法的研究思路在堆芯計算領域，尤其在國內尚處于無統一標準的局部應用階段(例如中國工程物理研究院的JASMIN軟件)，且缺乏細致的程序設計和對應的理論模型。

圍繞著堆芯計算軟件的應用重點，論文描述了適用于堆芯計算的框架原型軟件NAC4R的設計和初步編程實現，并列舉抽象程序的驗證和實際屏蔽輸運計算軟件Hydra-SN3D的重制工作來說明框架方法的有效性。

1 理論和設計

數值計算與并行計算環境可以按照不同層次進行融合[10]。例如在計算科學領域，中間件的概念來源于分布式計算技術，它基于網絡計算平臺上部署的分布式計算環境，用于提供開發工具和公共服務，從而支持分布式通信應用。本文提及的框架方法則側重于數值計算的功能設計，同時包含部分中間件的內容以便集群計算機的數據通信，重點關注數值計算內數據的傳遞和代數運算。NAC4R作為框架方法的原型設計，其研制過程將專業應用范圍限制在堆芯計算，且運行平臺僅考慮內部封閉的集群計算機系統，所以與硬件的配合程度應單獨予以考慮。

1.1 理論說明

框架方法映射至軟件系統的不同層次，NAC4R的基本結構框圖如圖1所示。

圖1 NAC4R的結構框架Fig.1 The structure of NAC4R

該設計提供多種基礎計算功能，然后利用Python接口層來進行自動化處理和提供外部應用函數接口。現階段在原型階段對基礎計算功能的模塊進行如下分類：

1)基礎數據模塊，對控制和計算所用的數據進行合理表示，兼顧串行情況和分布式內存并行情況；

2)常用功能模塊，提供基本的非數值計算操作，例如數據封裝、快速排序、查找和圖分區等功能；

3)基礎線性代數操作，例如對BLAS函數接口規范、LAPACK函數接口規范提供包裹函數，便于直接使用FORTRAN函數。后續在接口規范一致的前提下也可以編程實現適用于特定并行硬件的基礎線性代數功能；

4)并行數據結構模塊，它由1)進行派生，基于著名的Epetra模塊，在分布式內存環境下對數據進行合理表示和存儲[11]；

5)線性方程組求解，現階段僅需提供Jacobi、Bi-CG、GMRES和LU分解的求解器，實現求解功能和考慮初步的健壯性，基于著名的Aztec模塊進行編程實現[12]；

6)輸入輸出模塊，對框架內部的輸入輸出功能進行編程實現，關注文本文件和二進制文件，同時提供框架內部固定的計算數據存儲格式；

7)事件追蹤模塊，提供對計算過程的時間線統計和計算時間分析，便于引入TAU等工具；

8)基礎數值計算模塊，提供基本的數值計算操作和數值函數，例如伽馬函數，勒讓德多項式展開系數，線性求積組等功能。

因此，NAC4R的基礎計算功能現階段由8個模塊組成，并強調使用C/C++語言和Fortran語言進行混合編程實現計算核心，一方面便于C/C++語言實現控制和封裝等管理特性，另一方面便于Fortran語言實現高效數值計算和利用已有的底層計算核心，例如BLAS、LAPACK等第三方內容。通過以上設計可以歸納出，NAC4R與框架方法的聯系體現在多項計算內容的抽象上。

(1)數據的表達

在集群計算機內部，數據結構主要分為三類：原語類型、復合類型和抽象數據類型。各種堆芯計算相關的具有物理意義的數據便由這些數據結構進行承載和存儲。NAC4R一方面提供專用的數據結構，諸如雙向鏈表、B-tree、紅黑樹以及圖，另一方面也根據堆芯計算的特點，強調數據有效的組織邏輯，表述為區域拓撲、粒子、網格等對象。這些數據的表達在NAC4R內部由類部件或者函數接口這種模塊形式來提供使用。

(2)數據通信

在分布式內存情況下實現高效的數據通信。現有的集群計算機在處理器核心之間主要使用MPI消息傳遞模型進行數據通信，為了開發用戶層面上函數接口的統一性和適應扁平化設計，并不直接使用MPI提供的函數，而是進行包裹函數庫(Wrapped Function Library)的研發設計，一方面減少輸入形參和實現功能的組合和統一管理，另一方面固定函數接口，使得包裹函數內部可以靈活更改，以適應不同的硬件平臺和更新或者替換底層的MPI庫(見圖2)。

圖2 包裹層示意圖Fig.2 Schematic of the wrapping layer

這種設計思路在PETSc、Trilinos、KMATHLIB和較多的框架軟件中均有體現，按照分布式通信集中管理的原則，NAC4R內的數據通信通過管理功能、數據轉換功能、通信區域分割等內容進行實現，它涉及MPI分布式內存并行和OpenMP或多線程共享式內存并行。

(3)文件I/O

在MPI-2規范中引入的并行I/O，也稱為MPI-IO，它是NAC4R數據通信模型中關于分布式文件讀寫的基礎。為了減小問題的復雜度，NAC4R對文件的讀寫僅僅針對二進制文件，并不針對文本文件(文本文件可以通過轉換接口來生成二進制文件)，且NAC4R也使用和擴展了Fisher所提出的SDF文件格式來減小數據存儲和傳遞的復雜性[13]。

(4)數值求解器

NAC4R的數值求解部分針對具體的堆芯計算特點，現階段已實現線性方程組Ax=b和非線性方程組的部分求解方法，例如圍繞Jacobi、Bi-CG、GMRES和LU分解的理論和算法進行設計和實現。為了提升適應性，當矩陣規模較小或精度要求較高時，使用直接法的接口更為有效。

1.2 編程考慮

原型軟件NAC4R內存在著面向結構和面向對象的軟件設計，考慮到應用需求，它自身存在著配置文件，可以指定具體的計算任務和集群計算機的資源等信息實現獨立運行。在NAC4R的編程實現過程中應考慮三種應用情況。

1)直接在上面集成專業應用的求解器，構成完整的專業軟件；

2)使用NAC4R提供的部分模塊，以函數庫的形式進行重用；

3)復用NAC4R的程序代碼，并利用提供的自動化功能進行并行性能預測和收集計算信息。

1.3 原型研制步驟

基礎模塊的開發為第一階段，體現為NAC4R在研制過程中實施單元測試，然后開發相應的示例程序和腳本語言包裹層。

第二階段對實際的堆芯樣本程序進行應用和分析，本文選取結構化網格條件下的三維離散縱標法的并行算法進行應用研究。針對其中的徑向掃描并行算法部分，樣本程序被解讀和分析，并在NAC4R原型軟件中快速進行重新開發，理解核心并行求解算法的編寫流程和效率提升模式(見表1)。

表1 樣本程序Hydra-SN3D的信息

第三階段是持續改進原型軟件NAC4R，并擴大應用范圍，對新的堆芯計算特點進行特征識別，即在NAC4R內部集成新的并行算法和共性計算內容。

1.4 離散縱標法

美國LANL實驗室研制的一系列離散縱標程序是SN方法的起源，北京應用物理與計算數學研究所前期圍繞中子輸運SN算法開展了多方面的研究，并重點介紹二維中子輸運方程SN算法的研究與應用情況[14]。2010年美國LANL實驗室的技術報告中宣講了走鵑超級計算機上運行的SN輸運程序Sweep3D，它作為PARTISN中子輸運程序的縮減版本，可以計算XYZ三維笛卡爾坐標系的幾何系統，使用二維處理器網絡結構來匹配幾何的分塊網格，從多個層次上實施并行來加速SN算法的求解[15]。離散縱標法的形式通常由波爾茲曼輸運方程的離散縱標形式進行建模，可以反映中子或光子的輻射效應。

Ωm·Img+(σA+σS)Img

(1)

輻射通量Img(x,y,z)按照能群、角方向進行離散化，公式表明沿著一個特定方向Ωm輻射通量將發生的變化，求解過程便是定量獲得此變化趨勢。計算時只對選定的若干個離散方向Ωm對中子輸運方程進行求解。從中子輸運方程求出φ(r,E,Ωm)后，關于方向Ω的有關積分則用數值積分來近似表示。

(2)

式中,求積系數ωm、離散方向及其數目取決于計算精度的要求(SN方法，又寫作SN)。此處下標N表示方向向量在某個坐標方向上(例如XYZ三維笛卡爾坐標系的某個坐標軸方向)的離散點數目。近期中國核動力研究設計院和西安交大聯合研制Hydra-SN3D是類似于DORT核心求解器的離散縱標法程序。程序將整個方程抽象為矩陣系統進行求解，這種完整性更貼近于復雜系統的模擬和仿真。假定輸運算符為Η，裂變算符為Γ。

Ηφ=Ω·

(3)

(4)

根據公式(3)和(4)的算符定義，離散縱標法的源迭代求解算法為：

Input:?0,v∑f,∑s,∑tOutput:keff,?S=Γ?;while|k-kold|kold≥εkor||S-Sold||||Sold||≥ε?do H?=S Sold=S S=Γ? k= S=1kS

對離散縱標法的框架方法實踐將基于上面描述內容實施抽象和編程實現。

2 抽象程序和數值實驗

框架方法同時強調抽象，為了有效在NAC4R內部實現對樣本程序的應用或研究，可以參照代理程序SNAP的研究思路[16]，根據問題域本身的特征進行抽象(本文強調堆芯計算)。一方面抽象其中的核心計算部分，按照NAC4R提供的程序部件和結構進行重寫;另一方面使用數學模型來模擬真實軟件的行為，研究新的模型、算法等。

2.1 性能預測

堆芯計算的一次完整計算可以按照功能和先后次序進行分解，總可以形成圖3所示的多個計算任務構成的網狀結構。

圖3 計算任務的網狀結構Fig.3 Mesh structure for computing tasks

此時每個計算任務可以使用一個進程或線程代表的計算資源完成計算，然后發送數據到它的后繼，其后繼在接收到所有所需數據之后開始計算。在數學上，這種互聯網狀結構可以通過有向圖或者Petri網進行表達和研究。有向圖是由一組定點和邊(V,E)進行定義的數學模型，當有向圖中不存在循環結構時，有向圖又被稱為有向無循環圖，它仍然遵循有向圖的各種定理，例如反向后序序列等價于圖模型中的拓撲排序，反映了各個計算任務的一種調度關系。當有向無循環圖的連接邊E不存在權重值時，可以直接由有向圖對應的矩陣M進行調度算法的設計。

1)Mi,j在初始情況下為0，當單元被調度完成計算后，可以置Mi,j為1；

2)Mi,j(i≠j)表明了單元j依賴的前驅單元，如果Mi,j≠0，則單元j的調度需要等待單元i完成計算。

離散縱標輸運計算的徑向掃描計算具備圖4所示的流水線結構，它在SN輸運計算過程中可能由多進程對較大的分塊區域按照此互聯關系進行并行求解，也可能由多線程對分塊區域內部的網格按照此互聯關系進行并行求解，前者會產生進程之間的數據通信開銷，后者則共享數據內存，只存在著計算的先后次序和多線程的管理開銷。

圖4 SN方法的流水線結構Fig.4 Pipeline of the SN method

假定每個網格區域的體積一致時，此時每個網格區域的計算任務的計算量相同，當數據量較小且不考慮數據通信時，對應于多線程并行情況，此時理想并行性能(加速比)為：

(5)

理想性能的取得需要計算資源(進程或線程)滿足≥min{i,j}。當計算資源并不充足時，可以根據上述抽象和NAC4R提供的并行性能預測工具生成不同概率分布的輸入，模擬資源數的增長與實際的并行性能，獲得指導性建議來輔助并行編程。以I=50,J=50為例，可以使用工具模擬，按照時間片輪轉的方式得到以下的運行時間趨勢。

圖5顯示的趨勢表明每個網格區域的計算量相同，且不考慮數據通信時，當計算資源在10(十個進程或線程)以內時，該例題規模下并行所帶來的增益最為明顯。考慮到預測的復雜性，還可以采用Petri網數學模型對更為復雜的情況進行性能預測，例如當每個網格區域的計算量不相同，考慮數據通信之后的多進程并行性能。2.1節解釋的框架方法反映的是一種程序抽象的能力，并行算法的性能預測是一種用途。

圖5 SN徑向掃描的并行性能預測Fig.5 Parallel performanceprediction of radial scanning

2.2 離散縱標法的徑向并行掃描

與2.1描述的抽象保持一致，真實的離散縱標法的掃描算法在每個卦限內每個離散方向均需要遍歷掃描各個網格，它在結構化網格的條件下應用KBA并行掃描算法。此時求解過程對應于通量矩的線性系統的代數運算。

(6)

L——差分輸運算符;

M——通量矩離散化算符矩陣;

S——散射矩陣。

待求解對象ψ代表角通量，φ為通量矩。結構化網格情況下KBA掃描計算的實質是求解L-1，它將三維網格進行區域劃分，假設總的網格在每個方向上為(I,J,K)，各個區域(Ia,Jb,K)被映射至每個進程計算任務進行并行求解。此時每個卦限內的各區域XY方向上的掃描計算由多層循環結構表示，它的算法描述為：

I,J,K←0,1,2for k in d_begin[K]:d_end[K] ?up,j←?in for j in d_begin[J]:d_end[J] for i in d_begin[I]:d_end[I] ?psi←?up callT(i,j,k,?psi) end for ?up←?psi end for ?out←?up,j end for ?out←?up,k

上述掃描計算過程中，外推模型的計算子函數T具備獨立性，僅僅與網格序號(i,j,k)相關，通量矩變量φup、φpsi等存在著先后依賴關系。對于該流水線并行的并行算法，可引入全局標志位來輔助流水線結構的執行，此時的循環結構使用OpenMP多線程并行的算法可以設計為：

#pragma omp parallel forfor j in d_begin[J]:d_end[J] for i in d_begin[I]:d_end[I] wait(done[j][i]) ?psi←?up call T(i,j,k,?psi) done[next(j)][i]=1 end forend for

框架方法在原型軟件NAC4R中單獨提供一類機制來統一這種流水線結構的執行。一方面由于循環結構并行算法的不恰當編程實現將造成性能的惡化，另一方面以上算法結構較難高效和正確編寫成功。NAC4R統一使用OpenMP提供的omp_set_lock和omp_unset_lock操作提前分配和管理互斥數據結構，便于OpenMP多線程在流水線結構下并發運行。計算用戶只需要編寫循環結構本身，將計算和數據表示相互獨立，也便于后續保留效率較高的程序編寫模式。表2統計了編譯優化等級為-O0時的數值實驗計算效率，兩個預設例題的計算結果在樣本程序、NAC4R改寫之后的程序及其并行版本均保持一致。

表2 預設例題的計算時間統計

兩個例題的單進程計算效率得到了較大提升，且對于每個基準例題，按照3.1節提供的預測信息，對每個進程使用8個線程實施的多線程并行均起到了接近圖5所示的加速效果。這給后續引入進程和線程多級并行的策略帶來信心。

表3則統計和對比了不同編譯優化等級下NAC4R內重制的離散縱標法模塊的計算時間，如果簡單按照前面列出的全局標志位來編寫OpenMP并行算法，O2優化等級將無法得到一致的數值結果。NAC4R內部統一提供的多線程鎖機制的管理方式確保了計算結果保持不變，且表中數據表明高編譯優化等級下流水線掃描結構仍然可以使用OpenMP多線程來提升計算性能。

表3 不同編譯優化等級的預設例題計算時間統計

樣本程序Hydra-SN3D在NAC4R內部的重制工作前后共花費2個月，共3個人力，其中主要程序開發人員僅花費1個月便完成Hydra-SN3D軟件的編制和測試。NAC4R這種統一結構的框架方法及其編程理念有利于程序的研制、改寫或重組，可復用的工程應用核心就在于提供了較多共性計算功能和使用單個模塊部件及其組合這種形式來組織數值計算內容。原型軟件NAC4R內部存在較多可復用的程序片段，它們是以往并行軟件研制經驗和模式的總結，使得開發用戶花費更多時間對問題進行抽象，匹配已有模式，從而縮減編寫程序的時間。

3 結論

隨著應用的深入和計算硬件的改善，超大規模集群計算機下堆芯計算軟件的呈現將產生較大變化。NAC4R原型軟件作為框架方法的研究思路在本文進行了介紹，通過討論，得到如下結論。

1)框架方法可以取得軟件復用的目標，堆芯計算經過抽象來適應已有的模式或基礎計算功能可以簡化數值程序的編寫流程，且便于并行算法相關內容的分析。提出的NAC4R原型軟件對框架方法的實例化初步可行，從基礎功能上給出了較多計算的理論來源和編程實現考慮。

2)離散縱標法的屏蔽輸運程序Hydra-SN3D在NAC4R內被重新抽象和編程實現，通過數值實驗和統計數據論證了框架方法應用于堆芯計算的可行性和有效性。

后續將進一步對原型軟件進行發展，以適應更多堆芯計算內容，例如非結構化網格的離散縱標法程序和兩相流仿真程序等，擴大框架方法的適用范圍和提煉其中的并行模式。