,2

1.中國科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心,北京 100190

2.中國科學(xué)院大學(xué)，北京 100049

引言

隨著核工業(yè)的持續(xù)發(fā)展與核電技術(shù)的不斷進(jìn)步，反應(yīng)堆物理的研究與發(fā)展日益重要。反應(yīng)堆物理主要研究反應(yīng)堆內(nèi)大量中子的統(tǒng)計(jì)行為，在反應(yīng)堆物理中用來描述中子在介質(zhì)中相互作用的，就是中子輸運(yùn)方程。

作為反應(yīng)堆系統(tǒng)設(shè)計(jì)與計(jì)算的基礎(chǔ)，中子輸運(yùn)方程求解方法可以分為確定論方法和統(tǒng)計(jì)方法。由于中子在堆內(nèi)的行為與中子截面，中子運(yùn)動(dòng)方向，中子能量等多種因素有關(guān)，中子輸運(yùn)方程的求解非常復(fù)雜。傳統(tǒng)反應(yīng)堆物理以確定論方法為主。在確定論方法中，特征線方法有著高效的并行特性，且具有適用于任意維度、任意幾何構(gòu)造求解的優(yōu)點(diǎn)，這使得特征線方法成為現(xiàn)在研究堆芯計(jì)算的主要方法。

在過去很長一段時(shí)間內(nèi)，國內(nèi)外對(duì)特征線方法的研究主要是在二維特征線方法上，反應(yīng)堆設(shè)計(jì)和工程中也一直采用二維特征線計(jì)算作為中等保真度模擬方法，二維特征線方法發(fā)展日漸成熟。然而對(duì)于高保真建模，第三維是正確預(yù)測(cè)異質(zhì)反應(yīng)堆中的中子泄漏以及軸向功率分布所必需的，擴(kuò)展到三維可以更精確地計(jì)算軸向泄漏和反應(yīng)速率。因此科研與生產(chǎn)迫切需要將二維特征線方法擴(kuò)展到三維。

針對(duì)將二維特征線方法擴(kuò)展到三維特征方法帶來的耗時(shí)、耗存儲(chǔ)等問題，三維特征線的加速方法成為國內(nèi)外研究的重要方向。然而，單純地使用數(shù)值方法來加速特征線方法在某些情況下不收斂，隨著計(jì)算規(guī)模的變化加速效果也會(huì)受到影響，因此還需進(jìn)行其他加速方法的研究，并行計(jì)算就是其中最為重要的一個(gè)方法。隨著計(jì)算機(jī)技術(shù)的不斷進(jìn)步，高性能計(jì)算成為計(jì)算加速的主要手段，并行計(jì)算也成為近十年來反應(yīng)堆堆芯數(shù)值計(jì)算領(lǐng)域內(nèi)的一個(gè)重點(diǎn)研究方向。因此，將數(shù)值加速方法與并行計(jì)算相結(jié)合解決三維特征線方法計(jì)算的加速問題，在學(xué)術(shù)與實(shí)際工程應(yīng)用中皆具有非常重要的意義。

基于以上分析，本文使用了一種簡化的幾何建模方式來完成特征線的布置，即在布置特征線時(shí)只存儲(chǔ)二維平面上特征線的信息，三維特征線的信息通過二維平面上的特征線信息來實(shí)時(shí)生成，從而大大減少了在增加第三維的情況下對(duì)內(nèi)存的需求，并結(jié)合區(qū)域分解的并行方式較好的模擬了中子在介質(zhì)運(yùn)動(dòng)過程中的能源分布情況，為實(shí)現(xiàn)工業(yè)級(jí)別的三維特征線方法求解程序提供了一個(gè)良好的開端。

1 中子輸運(yùn)方程的特征線理論方法

特征線法是求解中子輸運(yùn)方程的玻爾茲曼形式的一種主流方法。通過沿特征方向?qū)Σ柶澛斶\(yùn)方程將極角和方位角進(jìn)行離散，并對(duì)特定方位角和極角求積的方程特征形式積分，可以得到玻爾茲曼方程的多組能群的特征形式如下

數(shù)值求解該方程近似化方法是將角度離散。首先，將方位角離散成 M 等分，同時(shí)給定不同方位角Ωm相應(yīng)的權(quán)重系數(shù)wm，其中 m?{1,.....M}，不同方位角的角通量相加得到所求區(qū)域的標(biāo)通量。其次，將極角離散成 P 等分，不同極角的權(quán)重系數(shù)為wp，其中

得到：

數(shù)值求解該方程的第二個(gè)近似是假設(shè)散射源是各向同性的，這使得總源可以使用標(biāo)通量來表示：

數(shù)值求解該方程的第三個(gè)近似是假定每個(gè)平源區(qū)的源為常量，這意味著源沿著特征方向 k 在平源區(qū)FSGi的進(jìn)入點(diǎn)s′與出射點(diǎn)s′出的源是不變的，即：

數(shù)值求解該方程的第四個(gè)近似是假設(shè)材料屬性在每個(gè)平源區(qū) FSR 上都是恒定的。

可得入射點(diǎn)和出射點(diǎn)之間通量的改變量為：

數(shù)值求解該方程的第五個(gè)近似是假設(shè)一條射線代表其附近的一小塊區(qū)域。

標(biāo)量通量的最終形式可以根據(jù)沿每個(gè)軌道段的角通量變化來簡化得到，有：

以上就是使用特征線方法所解決的中子輸運(yùn)方程的解的形式。

特征線方法的迭代求解過程可以概括為以下4點(diǎn)：

(1)首先將整個(gè)求解區(qū)域進(jìn)行網(wǎng)格劃分。

(2)使用經(jīng)過離散角度后的不同方向的特征線對(duì)區(qū)域進(jìn)行掃描，計(jì)算出射線與網(wǎng)格的交點(diǎn)，從而得到每一個(gè)射線線段的長度。

(3)根據(jù)邊界處已知的入射通量，逐段求出每一條射線線段的出射通量。

(4)通過散射源的內(nèi)迭代和裂變?cè)吹耐獾蠼庵钡绞諗俊?/p>

求解過程的流程圖如圖1所示。

2 特征線并行方法與實(shí)現(xiàn)

2.1 三維特征線生成算法

本文使用一種簡化的特征線布置方式，可以大大的減少特征線布置所需的內(nèi)存。首先生成徑向平面內(nèi)的特征線，軸向上的特征線在計(jì)算時(shí)再根據(jù)徑向上的特征線實(shí)時(shí)生成。

2.1.1 徑向平面特征線生成

為了保證 2D 軌道的循環(huán)纏繞，對(duì)于一個(gè)特定角度 ?i，在x和y軸上必須有整數(shù)個(gè)軌道數(shù)。

通過使用δ?的輸入值，計(jì)算沿特定角度 ?i的軸的軌道整數(shù)個(gè)數(shù)：

最后校正方位角和特征線之間的距離，并用公式(1)和公式(2)校正后的值進(jìn)行徑向平面的特征線布置。

2.1.2 軸向平面的特征線生成

為了生成 3D 軌跡，除了產(chǎn)生 2D 軌跡所需的參數(shù)之外，還需要軸向上軌跡間隔δθ和0＜θ＜π的極角數(shù)目nθ。利用這種方法創(chuàng)建極角來均勻地細(xì)分角度空間，使得 0＜θ＜π/2 中的每個(gè)極角與互補(bǔ)的極角配對(duì)。

其中θi,nθ-j-1是角度θi,j的補(bǔ)角。

為了保證三維軌道的循環(huán)軌道纏繞，必須滿足兩個(gè)條件：

(1)對(duì)于每個(gè)方位角 ?i，極角θi,j和2D 軌道周期沿著 2D 軌道周期的3D 軌道投影的開始和結(jié)束之間的距離必須是軌道間隔整數(shù)倍。

(2)對(duì)于每個(gè)方位角 ?i和極角θi,j，必須在幾何深度Δz 上沿著 z 軸具有整數(shù)倍的軌跡間隔。

第一個(gè)條件保證當(dāng) 2D 反射周期完成時(shí) 3D 軌道周期回繞到另一個(gè)3D 軌道上。第二個(gè)條件保證當(dāng)周期完成時(shí)，包含頂面或底面反射的3D 軌跡周期。

圖1 迭代求解算法流程Fig.1 Iterative solution algorithm flow

2.1.3 三維全局射線追蹤方法

使用每個(gè)方位角的循環(huán)軌道周期長度，可以計(jì)算在一個(gè)完整的2D 周期之后在一組 3D 軌道的開始和結(jié)束之間的軌道間隔的整數(shù)倍:

然后需要將沿著 z 軸的軌道間隔的數(shù)量設(shè)置為整數(shù)倍。通過考慮 z 方向上的軌道間隔的數(shù)量與沿著2D 軌道長度間隔之間的關(guān)系可以找到z 軸上的軌道的數(shù)量：

從而得到z 軸上特征線之間的距離：

最后校正極角和特征線之間的距離，并用使用公式(2.10)和公式(2.11)校正后的值布置軸向的特征線。

2.2 基于區(qū)域分解的并行設(shè)計(jì)

本文實(shí)現(xiàn)的程序采用 MPI 編程模型，并結(jié)合區(qū)域分解的并行方式。

空間分解采用分布式存儲(chǔ)模型,這對(duì)緩解完整核心問題所需的內(nèi)存負(fù)擔(dān)是必要的?？臻g分解背后的基本思想是將完整的核心問題分成足夠的空間域，以便每個(gè)域的大小都可以在單個(gè)進(jìn)程中管理。為了實(shí)現(xiàn)空間域分解，主要問題是跨空間子域邊界傳遞邊界條件信息的開銷。圖2 展示了在二維空間中，特征線在邊界處的數(shù)據(jù)交換，除非引入近似值，否則特征射線必須跨越這些邊界直接連接，以便傳輸界面角通量。

從特征線方法可知特征線方法的并行方式可以從能群、角度、射線以及幾何區(qū)域選取，結(jié)合本文所使用的特征線布置方式，本文采用區(qū)域并行的方式，使用 MPI 提供的切割函數(shù)，將求解的幾何區(qū)域進(jìn)行切分，將分解后的各個(gè)子區(qū)域分配至不同的進(jìn)程，從而減少每個(gè)進(jìn)程上的計(jì)算量。如圖3所示，這里假設(shè)將求解區(qū)域劃分為4個(gè)子區(qū)域。

上圖的4個(gè)子區(qū)域分配到對(duì)應(yīng)的4個(gè)進(jìn)程上，進(jìn)程間在反射邊界進(jìn)行數(shù)據(jù)通信，每個(gè)進(jìn)程可以接受來自上下、左右、前后 6個(gè)方向進(jìn)程的數(shù)據(jù)，進(jìn)程間數(shù)據(jù)交換的示意圖如圖4所示。

其中 scr_id為源進(jìn)程號(hào)，myid為當(dāng)前進(jìn)程號(hào)，dest_id為目的進(jìn)程號(hào)，f_psi、b_psi 分別代表消息傳遞的前后方向。

圖2 邊界數(shù)據(jù)交換Fig.2 Boundary data exchange

圖3 區(qū)域分解示意圖Fig.3 Area decomposition diagram

2.3 基于三維特征線方法的并行求解器

基于三維特征線方法的并行求解器的核心算法主要包括特征線的生成、幾何區(qū)域分解、傳輸掃描算法、進(jìn)程間數(shù)據(jù)通信。

2.3.1 特征線生成

特征線的生成是求解的第一步，它涉及徑向平面的特征線生成和軸向平面的特征線生成兩個(gè)方面。首先是徑向平面的特征線生成過程：根據(jù)輸入的方位角的個(gè)數(shù)以及特征線之間的間隔來生成特征線，為不同方位角的特征線分配不同的權(quán)重，表1給出了二維平面特征線生成的偽代碼算法。

根據(jù)徑向平面某一方向的特征線作為軸向特征線在徑向上的投影，軸向上根據(jù)不同極角來生成特征線。表2給出了軸向平面特征線生成的偽代碼。

圖4 進(jìn)程間數(shù)據(jù)交換示意圖Fig.4 Schematic diagram of data exchange between processes

表1 徑向平面特征線生成算法偽代碼Table1 Radial plane feature line generation algorithm pseudo code

2.3.2 幾何區(qū)域分解

本文中使用 MPI_Cart_create 函數(shù)來分割需要求解的區(qū)域，然后通過 MPI_Cart_shift 函數(shù)來獲取當(dāng)前進(jìn)程的相鄰進(jìn)程號(hào)，表3給出了相應(yīng)的偽代碼算法。

上述偽代碼中 ndims 表示劃分后的維度，ndims表示不同維度分配的進(jìn)程數(shù)，period表示進(jìn)程通信是否首尾相鄰。

2.3.3 掃描迭代算法

掃描迭代算法是整個(gè)程序的核心部分，其主要是計(jì)算對(duì)于某個(gè)網(wǎng)格內(nèi)某一方向上的特征線對(duì)該網(wǎng)格的貢獻(xiàn)度，即平均角通量，然后通過對(duì)不同方向的貢獻(xiàn)度進(jìn)行規(guī)約計(jì)算得到該網(wǎng)格的標(biāo)通量的一個(gè)過程，最后通過平均角通量來更新散射源的內(nèi)迭代和通過標(biāo)通量來更新裂變?cè)吹耐獾鷣砬蠼鈫栴}知道問題收斂。更具體而言，對(duì)于某一條具體的特征線，其主要是計(jì)算該條特征線在網(wǎng)格內(nèi)的通量改變量，從而得到不同方向的平均角通量。表4 列出了單次特征線掃描迭代算法。

表2 軸向平面特征線生成算法偽代碼Table2 Axial plane feature line generation algorithm pseudo code

上述傳輸掃描算法涉及方位角、軌跡、不同F(xiàn)SR中的段，能量組和極角上的五個(gè)嵌套環(huán)路偽代碼中符號(hào) m、k、s、g、p 分別表示方位角號(hào)、特征線號(hào)、線段號(hào)、能群號(hào)、極角號(hào)；相對(duì)應(yīng)的大些字母表示其集合。

在計(jì)算完某網(wǎng)格內(nèi)的平均角通量后，我們需要更新散射源和裂變?cè)矗源藖砼卸▎栴}是否收斂。表5中給出了如何更新源的偽代碼。

2.3.4 邊界數(shù)據(jù)通信

每個(gè)進(jìn)程內(nèi)都單獨(dú)的進(jìn)行傳輸掃描算法，當(dāng)傳輸掃描算法結(jié)束后在邊界上進(jìn)程間需要進(jìn)行數(shù)據(jù)通信，從而在進(jìn)行下一輪的傳輸掃描算法，表6給出了進(jìn)程間數(shù)據(jù)通信的偽代碼。

上述偽代碼中 Send Buffer、Recv Buffer 分別代表發(fā)送和接收消息的隊(duì)列，Number of Elements 代表發(fā)送的數(shù)據(jù)大小，send_dest[i]、rec_sources [j]分別代表目的進(jìn)程號(hào)和源進(jìn)程號(hào)。

表3 幾何區(qū)域分解算法偽代碼Table3 Set region decomposition algorithm pseudo code

表4 特征線掃描算法偽代碼Table.4 Feature line scan algorithm pseudo code

表5 源更新算法偽代碼Table.5 Source update algorithm pseudo code

表6 進(jìn)程間數(shù)據(jù)通信算法偽代碼Table6 Interprocess data communication algorithm pseudo code

3 程序性能測(cè)試與分析

3.1 測(cè)試平臺(tái)

測(cè)試平臺(tái)為中科院計(jì)算機(jī)網(wǎng)絡(luò)信息中心的超級(jí)計(jì)算機(jī)“元”上的部分節(jié)點(diǎn)，節(jié)點(diǎn)的具體配置表如 7所示。

我們分別對(duì)程序進(jìn)行了弱擴(kuò)展、強(qiáng)擴(kuò)展二個(gè)方面的性能測(cè)試，一個(gè)進(jìn)程運(yùn)行在一個(gè)節(jié)點(diǎn)上，具體測(cè)試結(jié)果見后續(xù)小節(jié)。

3.2 弱擴(kuò)展性測(cè)試結(jié)果

弱擴(kuò)展性的測(cè)試分別在1個(gè)進(jìn)程、2個(gè)進(jìn)程、4個(gè)進(jìn)程、8個(gè)進(jìn)程、16個(gè)進(jìn)程下進(jìn)行測(cè)試，相對(duì)應(yīng)的測(cè)試規(guī)模如表8所示。

測(cè)試結(jié)果結(jié)果如圖5所示。

由圖5 可知，隨著求解規(guī)模的增大，程序運(yùn)行時(shí)間略有上升，但整體的求解時(shí)間并沒有大幅度的波動(dòng)，這一定程度上表明我們程序的弱擴(kuò)展性較好。

表7 測(cè)試平臺(tái)參數(shù)Table7 Test platform parameters

表8 弱擴(kuò)展測(cè)試程序規(guī)模參數(shù)表Table8 Weak extension test program size parameter table

圖5 弱擴(kuò)展測(cè)試效果圖Fig.5 Weakly extended test renderings

表9 強(qiáng)擴(kuò)展測(cè)試程序規(guī)模參數(shù)表Table9 Strong extension test program size parameter table

表10 相對(duì)單個(gè)進(jìn)程運(yùn)行時(shí)間的時(shí)間加速比Table10 Time-to-acceleration ratio relative to the running time of a single process

圖6 強(qiáng)擴(kuò)展測(cè)試效果圖Fig.6 Strong expansion test renderings

3.3 強(qiáng)擴(kuò)展性測(cè)試結(jié)果

強(qiáng)擴(kuò)展性的測(cè)試分別在1個(gè)進(jìn)程、2個(gè)進(jìn)程、4個(gè)進(jìn)程、8個(gè)進(jìn)程、16個(gè)進(jìn)程下進(jìn)行測(cè)試，測(cè)試規(guī)模參數(shù)如表9所示。

測(cè)試結(jié)果如圖6所示。

表10給出了隨著進(jìn)程數(shù)增加實(shí)際的加速效果比。

由圖6 與表10 可知，在求解問題規(guī)模不變的情況下，求解的時(shí)間隨著進(jìn)程數(shù)的增加而減少。在4個(gè)進(jìn)程以內(nèi)時(shí)，程序表現(xiàn)出很好的擴(kuò)展性，但是在擴(kuò)展到8個(gè)進(jìn)程以上時(shí)，由于每個(gè)進(jìn)程分到的計(jì)算量太少而無法重復(fù)發(fā)揮每個(gè)節(jié)點(diǎn)的性能，導(dǎo)致程序的并行效果受到影響。

4 總結(jié)

作為核反應(yīng)堆系統(tǒng)設(shè)計(jì)與計(jì)算的基礎(chǔ)，堆芯計(jì)算的方式在一直不斷地得發(fā)展。在反應(yīng)堆計(jì)算的諸多數(shù)值方法當(dāng)中，特征線方法以其良好的并行性、幾何適應(yīng)性強(qiáng)等優(yōu)點(diǎn)，成為國內(nèi)外研究的重點(diǎn)。

本文實(shí)現(xiàn)了一個(gè)基于三維特征線方法的求解中子輸運(yùn)方程的計(jì)算程序，并 將該程序在多核 CPU 節(jié)點(diǎn)上進(jìn)行測(cè)試。實(shí)驗(yàn)結(jié)果表明程序的強(qiáng)弱擴(kuò)展性、加速比都達(dá)到了預(yù)期的效果，為后續(xù)設(shè)計(jì)更加復(fù)雜的利用三維特征線方法計(jì)算中子輸運(yùn)問題打下良好的基礎(chǔ)。