李睿涵 侯葉凡 李玉輝 劉召遠 張海紅 梁金剛

1（清華大學核能與新能源技術研究院 北京 100084）

2（齊魯工業(yè)大學（山東省科學院） 濟南 250353）

（li-rh20@mails.tsinghua.edu.cn）

隨著能源和環(huán)境問題日益嚴峻，社會亟需發(fā)展區(qū)別于傳統(tǒng)化石能源的新能源.其中，核能由于其穩(wěn)定性好、能量密度大、碳排放量低的優(yōu)點而成為解決能源與環(huán)境問題的有力手段.此外，由于以上優(yōu)點，核能還是航天器、潛艇、水面艦艇等特種設備的良好供能方式.因此，核能存在著較大的發(fā)展機遇.

由于核能實驗研究成本較高，故計算機模擬計算成為了核能領域的重要研究方式.其中，以反應堆堆芯核反應率分布為主要研究對象的“反應堆物理”計算有助于獲得反應堆有效增殖因數(shù)（keff）、功率分布等關鍵參數(shù)，并為堆芯熱工計算、核電廠系統(tǒng)計算等其他模擬研究提供必要的參數(shù)，是核能模擬計算研究的重要組成部分.蒙特卡羅（Monte Carlo，MC）方法[1]是使用計算機生成隨機數(shù)來模擬中子等微觀粒子在堆芯內(nèi)的運動、散射、吸收等隨機行為的計算方法，其數(shù)學模型較為貼近堆芯實際的物理過程，故計算結果較為精確.但同時MC 方法需要對堆芯的大量粒子分別進行模擬，因此也存在著計算量大的缺點.而由于堆芯中各粒子的行為是相對獨立的，故MC 方法可以將不同粒子直接分配給不同線程（進程）進行計算，這樣的并行模式天然具有非常好的并行可擴展性，與之相比，基于區(qū)域分解的并行模式的并行效率隨并行規(guī)模的增加而顯著降低[2]，因此是一種非常適合并行計算的方法.目前，MC 方法已經(jīng)在反應堆物理計算領域得到了廣泛地應用[3-6].

當前，隨著核技術的發(fā)展，一種新型的核燃料：顆粒燃料逐漸得到了越來越廣泛的應用.諸如球床式高溫氣冷堆[7]、空間反應堆[8]、氟鹽冷卻高溫堆[9]、耐事故燃料[10]等多種先進核能裝置中均使用了這種燃料.顆粒燃料是將核燃料活性成分，例如UO2制成小顆粒，然后彌散在基體中，而非傳統(tǒng)棒狀燃料使用整塊均勻的UO2.可見，顆粒燃料的使用必將導致堆芯幾何中存在大量隨機球體.以一種高溫氣冷堆HTRPM（high-temperature gas-cooled reactor pebble-bed module）為例，其堆芯存在約42 萬個隨機堆積的球形燃料元件，每個燃料元件中又包含約1 萬個隨機排列的燃料顆粒[7]，如圖1 所示，因此其全堆總共存在約40 億個隨機排列的球體.如此復雜的幾何結構會為反應堆物理計算，尤其是MC 計算帶來極大的困難，大大增加其計算耗時.例如經(jīng)本文驗證，在相同的MC 計算規(guī)模（即模擬的中子數(shù)量，根據(jù)統(tǒng)計學大數(shù)定律，計算的中子數(shù)越多，結果越精確，標準差越?。┮约安⑿幸?guī)模下，如果不使用任何加速方法，一個HTR-PM 的球形燃料元件的計算耗時將達到一個傳統(tǒng)壓水堆棒狀燃料元件的400 倍以上.因此有必要開展適用于顆粒燃料幾何的加速方法研究.

本文以其中幾何較為復雜的球床式高溫氣冷堆HTR-PM[7]為例，在開源MC 程序OpenMC[11]的基礎上進行加速方法的研究.

HTR-PM 是由清華大學核能與新能源技術研究院主導研發(fā)的球床式高溫氣冷堆，具有第4 代先進核能系統(tǒng)[12]特征.因其具有良好的固有安全性，HTR-PM受到了學界的廣泛認可.其示范堆山東石島灣核電站屬于我國16 個國家科技重大專項之一，目前已經(jīng)完成了建設并實現(xiàn)了并網(wǎng)發(fā)電.因此，針對HTR-PM的研究將具有很高的工程應用價值.

OpenMC 是由美國麻省理工學院主導開發(fā)的開源MC 中子輸運程序.除麻省理工學院外，阿貢國家實驗室、密歇根大學等多個研究機構也對其開發(fā)工作做出了貢獻.OpenMC 程序具有較高的幾何建模和數(shù)據(jù)統(tǒng)計靈活性，并具備高性能并行能力，支持MP?（message passing interface）和OMP（open multi-processing）混合并行.作為開源程序，OpenMC 自2011 年問世以來已經(jīng)在世界范圍內(nèi)進行了大量的反應堆物理計算，其各項功能得到了學界的充分驗證，具有較高的可靠性.因此OpenMC 適合用于開展本項研究.

除HTR-PM 外，本文還建立了一種使用顆粒燃料的空間反應堆模型，以對優(yōu)化后的程序進行輔助驗證.

1 模型建立

1.1 隨機球體位置的生成

為建立顆粒燃料的高精度模型，必須準確設定每個燃料顆粒的位置.對于HTR-PM 中的數(shù)十億燃料顆粒來說，通過實驗方法測量其位置顯然是不現(xiàn)實的，因此顆粒的位置也必須通過計算機模擬計算得出.

OpenMC 程序本身具備生成隨機顆粒位置的功能.其使用的是RSP（random sequential packing）和CRP（close random packing）結合的方法[13].RSP 是在規(guī)定的空間內(nèi)隨機選取球體位置并去除重疊的球體的方法.這一方法在填充率較低時擁有較高的計算效率，但其無法生成38%以上填充率的隨機顆粒，且其生成的顆粒是均勻分布在空間內(nèi)的，無法模擬重力作用下的球堆.CRP 方法則可以生成填充率更高的隨機顆粒，但其計算效率較低，同時也無法模擬重力的影響（由于本文不涉及CRP 方法，故對其原理不再贅述）.OpenMC 在生成隨機顆粒時，會根據(jù)目標填充率進行判斷，如果填充率較低，則使用RSP 方法，否則將使用CRP 方法.除此之外，OpenMC 還支持用戶從外部輸入顆粒的位置，為其與其他程序的耦合提供了可能.

HTR-PM 堆芯中包含2 重隨機球體，首先是堆芯容器中堆積的燃料球，其次是每個燃料球中彌散的燃料顆粒，如圖1 所示.在建立OpenMC 模型時，可以首先建立單個燃料球及其內(nèi)部的燃料顆粒的模型，然后使用“重復結構”功能將該燃料球重復地堆積在堆芯容器內(nèi).如此一來，HTR-PM 模型中的隨機球體就分為了“堆芯容器中的燃料球”和“每個燃料球中的燃料顆?！? 部分.因此建模時無需生成一個40 億量級的球堆，而只需生成一個40 萬量級和一個1 萬量級的球堆即可.對于單個球內(nèi)的燃料顆粒，其填充率約為7%，且均勻分布在燃料球基體中，因此非常適合使用OpenMC 自帶的RSP 方法.而燃料球則是在重力的作用下堆積在堆芯容器內(nèi)，且填充率約為61%，因此不適合直接使用OpenMC 進行計算，而需要單獨進行處理.

本文使用了DEM（discrete element method）程序LAMMPS[14]來生成燃料球的位置信息.DEM 方法是使用經(jīng)典力學模擬每個顆粒運動過程的方法.在LAMMPS 模型中，首先設置好堆芯容器的幾何，然后在高處生成燃料球，再令其在重力作用下自然下落并堆積在容器內(nèi)即可.由于這一物理模型較為接近HTR-PM 的實際裝料過程，因此可以得到較為準確的結果.LAMMPS 模型的具體參數(shù)與文獻[15]中的一致.

值得一提的是，HTR-PM 模型并非包含42 萬個燃料球的平衡堆芯模型，而是包含約5 萬個燃料球和28 萬個石墨球的初裝堆模型[16]，這是因為當前HTRPM 示范堆山東石島灣核電站投產(chǎn)不久，并未達到平衡堆芯狀態(tài)，而僅有初裝堆的運行數(shù)據(jù).因此，建立初裝堆模型可以與石島灣核電站的實驗數(shù)據(jù)進行對比，以更加有力地驗證程序的準確性.

初裝堆狀態(tài)下，堆芯的燃料球和石墨球并非均勻混合，而是先在堆芯底部堆積約23 萬個石墨球，然后在其上方混合堆積約5 萬個燃料球和5 萬個石墨球.由于LAMMPS 程序只關注球的經(jīng)典力學特征，不關注其核物理性質(zhì)，而燃料球和石墨球的經(jīng)典力學特征，如質(zhì)量、摩擦系數(shù)、彈性系數(shù)等相近，因此這2 種球在LAMMPS 模型中不作區(qū)分.在LAMMPS給出球的位置后，再人工確定其中哪些是石墨球、哪些是燃料球，最后將數(shù)據(jù)輸入到OpenMC 中.

LAMMPS 程序得出的球床堆積率為61%，與文獻[16]中的符合較好.

與HTR-PM 不同，本文模擬的空間反應堆[17]堆芯為包含1 519 個圓形冷卻劑孔道的六棱柱，其幾何中僅存在1 種隨機球體，即彌散在堆芯固體區(qū)域內(nèi)的大量燃料顆粒.該堆芯中顆粒填充率約為30%，且均勻分布在基體中，因此使用RSP 方法計算即可.但是，由于全堆燃料顆粒數(shù)量極多（約3 億個），故直接使用OpenMC 在整個堆芯內(nèi)部隨機生成顆粒位置是不實際的，此外OpenMC 還不支持在帶有冷卻劑孔道的復雜幾何中生成隨機顆粒位置.綜上所述，本文在建立空間堆模型時同樣使用了“重復結構”功能：將一個冷卻劑孔道周圍的六棱柱區(qū)域視為一個單元，使用RSP 方法在此實心六棱柱內(nèi)部隨機生成燃料顆粒，并人工去除生成在冷卻劑孔道中的顆粒位置，然后將1 519 個這樣的單元重復填充在堆芯中，即可實現(xiàn)空間堆隨機顆粒的建模.

1.2 HTR-PM 與空間堆全堆模型

在1.1 節(jié)所述的生成隨機顆粒位置方法的基礎上，本文根據(jù)石島灣核電站初裝堆相關參數(shù)[16]，建立了HTR-PM 堆芯的OpenMC 模型，如圖2 所示.并根據(jù)文獻[17]建立了空間堆堆芯的OpenMC 模型，如圖3 所示.

Fig.2 OpenMC model of the HTR-PM core圖2 HTR-PM 堆芯的OpenMC 模型

Fig.3 OpenMC model of the space reactor core圖3 空間堆堆芯的OpenMC 模型

2 加速方法

2.1 網(wǎng)格加速的原理

前文已經(jīng)提到，大量隨機球體的幾何結構會大大增加MC 方法的計算量.究其原因，是程序無法對隨機排列的球體進行有效地檢索.例如在圖4（a）所示的情形中，大量隨機排列的燃料顆?？臻g中存在一個中子，而程序需要確認中子具體位于哪個顆粒內(nèi).為此，程序只能將中子的坐標一一帶入顆粒的幾何信息中進行驗證，如式（1）所示：

Fig.4 Basic principle of lattice speed up圖4 網(wǎng)格加速的基本原理

其中r0是中子位置，ri是第i個燃料顆粒的球心位置，Ri是第i個燃料顆粒的半徑.如果對于某一燃料顆粒i，式（1）得到滿足，則中子位于該顆粒內(nèi)或其表面，而如果中子不位于任何燃料顆粒內(nèi)，則中子位于基體中.

顯然，由于模型中顆粒的數(shù)量極多，以上“遍歷所有顆?！钡倪^程必然耗費大量時間，因此有必要進行優(yōu)化.

針對這一問題，“網(wǎng)格加速”方法是一種有效的手段.如圖4（b）所示，在空間中劃分規(guī)則的網(wǎng)格，使得每個單元格的范圍為

其中，(x0,y0,z0) 為網(wǎng)格空間角點坐標；Px，Py，Pz分別為x，y，z方向的網(wǎng)格柵距；i，j，k分別為x，y，z方向的網(wǎng)格編號.

完成網(wǎng)格劃分后，即可建立單元格與燃料顆粒的對應關系，將與某單元格相交的所有燃料顆粒的幾何信息與該單元格編號 (i,j,k)關聯(lián)（該過程僅在程序初始化時進行，計算過程中無需重復）.如此，則當進行中子位置的判斷時，程序即可首先根據(jù)式（2）求出中子所在的單元格 (i,j,k)，然后僅遍歷該單元格內(nèi)的所有燃料顆粒即可.例如圖4（b）所示的情形中，只需對3 個實線的燃料顆粒進行驗證，而無需驗證虛線的顆粒，這就大大節(jié)省了計算時間.以上即為網(wǎng)格加速方法的基本原理.

2.2 虛擬網(wǎng)格與真實網(wǎng)格

當前OpenMC 中已經(jīng)支持網(wǎng)格加速方法，但其使用的是“真實網(wǎng)格”，而本文為OpenMC 添加的是“虛擬網(wǎng)格”的方法，本節(jié)將介紹二者的不同.

在真實網(wǎng)格模型中，網(wǎng)格分界面是模型中真實存在的曲面，模型將被網(wǎng)格面劃分為若干長方體.如圖5（a）所示，使用真實網(wǎng)格方法將模型劃分為左右2 個單元格，則最終的模型是由左右2 個長方體拼接而成的.

Fig.5 Division patterns of different lattices圖5 不同網(wǎng)格的劃分方式

而虛擬網(wǎng)格則是一套假想的網(wǎng)格，其模型幾何中不存在網(wǎng)格面，不會對模型進行分割.如2.1 節(jié)所述，網(wǎng)格加速方法本質(zhì)上是建立一種各球體與其所在區(qū)域的對應關系，以便根據(jù)粒子位置縮小所需遍歷的球體數(shù)量.因此，網(wǎng)格面的存在是沒有必要的，只需按照網(wǎng)格位置對各球體進行編號即可.例如圖5（b）所示的情形下，使用虛擬網(wǎng)格劃分模型，則只需對左側(cè)4 個球體添加編號“1”，代表其屬于左側(cè)單元格；對右側(cè)4 個球體添加編號“2”，代表其屬于右側(cè)單元格；對中間球體同時添加編號“1，2”，代表其同時屬于2 個單元格.模型的幾何不存在變化.

虛擬網(wǎng)格方法與真實網(wǎng)格相比，具有4 個優(yōu)點：

1）虛擬網(wǎng)格對存儲空間的占用更小.這是因為真實網(wǎng)格的模型幾何由大量的單元格幾何體拼接而成，其復雜度遠超虛擬網(wǎng)格模型的單一幾何體.以圖2所示的HTR-PM 堆芯模型為例，為達到良好的加速效果，需要做出較細的網(wǎng)格劃分，例如將堆芯容器內(nèi)堆積的燃料球以及石墨球劃分為33×33×154 的網(wǎng)格，同時將單個燃料球內(nèi)的燃料顆粒劃分為24×24×24 的網(wǎng)格.顯然，如果使用真實網(wǎng)格，則如此大規(guī)模的單元格必將占據(jù)大量存儲空間.表1 展示了該網(wǎng)格劃分下虛擬網(wǎng)格與真實網(wǎng)格模型占用的存儲空間對比，可見使用虛擬網(wǎng)格方法將大大節(jié)約存儲空間.

Table 1 Storage Space Occupation of HTR-PM Model表1 HTR-PM 模型存儲空間占用

2）虛擬網(wǎng)格的計算速度更快.具體表現(xiàn)在2 方面.①在程序初始化階段，雖然虛擬網(wǎng)格模型需要額外對球體進行編號，但由于其模型幾何大大簡化，使得初始化的整體耗時遠小于真實網(wǎng)格模型.②在程序計算階段，由于虛擬網(wǎng)格模型不存在網(wǎng)格面，故程序無需對中子穿過這些網(wǎng)格面的過程進行模擬，因此虛擬網(wǎng)格模型也將在計算階段節(jié)約時間.表2 展示了相同計算規(guī)模和并行規(guī)模下，HTR-PM 堆芯模型的計算時間對比.可見虛擬網(wǎng)格模型的初始化時間和計算時間分別約為真實網(wǎng)格模型的13%和82%，擁有更高的計算效率.

Table 2 Time-Consumption of HTR-PM Model表2 HTR-PM 模型計算耗時 s

表3 則展示了空間堆模型的計算時間對比，其虛擬網(wǎng)格的初始化時間和計算時間分別為真實網(wǎng)格的13%和46%，與HTR-PM 模型相比效率的提升更為明顯.

Table 3 Time-Consumption of the Space Reactor Model表3 空間堆模型計算耗時 s

3）虛擬網(wǎng)格的數(shù)據(jù)統(tǒng)計更加便捷.由于真實網(wǎng)格模型存在網(wǎng)格面對模型進行分割，其難免會切斷模型中隨機分布的球體，使得一個完整的球在模型中實際由2 個，甚至多個“半球”拼接而成.例如圖5（a）中，中間的球體被網(wǎng)格面分割為了2 部分.因此如果用戶需要以單個球為單位統(tǒng)計某參數(shù)（例如功率分布），其必須人工將多個“半球”指定為一個整體進行統(tǒng)計，這為用戶操作帶來了困難.而虛擬網(wǎng)格模型中不存在網(wǎng)格面，其所有球均為一個完整的球體，因此用戶可以便捷地進行數(shù)據(jù)統(tǒng)計.

4）使用虛擬網(wǎng)格加速的模型擁有更好的并行可擴展性.這是因為虛擬網(wǎng)格模型幾何更簡單，使得單個中子輸運過程的耗時更短、隨機性更小.由于在并行計算中，程序只有在每個進程各自完成其“負責”的全部中子的計算后，才會進入下一階段，故更小的隨機性就意味著進程間相互等待的時間更短，即程序擁有更好的負載均衡.圖6（a）展示了HTR-PM 虛擬網(wǎng)格和真實網(wǎng)格的并行強可擴展性對比，其中每代中子數(shù)保持10 752 000 不變，總代數(shù)為20，并行規(guī)模為112～10 752 核；圖6（b）為弱可擴展性對比，其中每代中子數(shù)為核數(shù)的1 000 倍，總代數(shù)為20，并行規(guī)模為112～131 600 核.本文針對圖6 的每個數(shù)據(jù)點均進行了10 次重復測量，并標出了其平均值以及10 次測量結果的變化范圍.其中某次測量的并行效率是用測量得到的計算時間與基點（112 核）下10 次測量的平均時間計算得出的.由于數(shù)據(jù)的波動性，某些“并行效率”可能超過了100%.

Fig.6 The parallel scalability comparison of the virtual lattice and physical lattice圖6 虛擬網(wǎng)格與真實網(wǎng)格并行可擴展性對比

可見，在強擴展問題中，當并行規(guī)模達到約3 000核時，虛擬網(wǎng)格的并行效率顯著高于真實網(wǎng)格；而對于弱擴展問題，當并行規(guī)模達到約3 000 核時虛擬網(wǎng)格的并行效率略高于真實網(wǎng)格；而達到約10 000 核時該趨勢則更為顯著.因此可以認為虛擬網(wǎng)格擁有更好的并行可擴展性.此外，真實網(wǎng)格測量結果的波動范圍明顯遠大于虛擬網(wǎng)格，這可能是因為真實網(wǎng)格數(shù)據(jù)量更大、中子歷史更復雜，其在緩存讀寫等方面具有更大的不確定性.

2.3 虛擬網(wǎng)格加速在OpenMC 程序的實現(xiàn)

在OpenMC 程序中，首先需要對初始化部分進行修改，以便生成虛擬網(wǎng)格.具體的方法為：根據(jù)式（2）的規(guī)則建立單元格編號與其位置的對應關系，并建立對應的3 維數(shù)組，使得數(shù)組編號與單元格編號一一對應.然后對于每一個單元格，遍歷所有隨機球體，如果某一球體與單元格存在交集，則程序?qū)⒃撉虻男畔?，包括其編號、球心位置、半徑儲存在單元格對應的?shù)組元素內(nèi).在對所有的單元格完成遍歷后，就成功建立了燃料顆粒的幾何信息與單元格編號的關聯(lián)，即生成了虛擬網(wǎng)格.

除生成虛擬網(wǎng)格外，還涉及計算過程中對虛擬網(wǎng)格的利用問題，需要對OpenMC 進行修改.

OpenMC 計算中的幾何處理過程包含3 項主要任務：

1）已知中子的位置，求其位于哪個幾何體內(nèi)；

2）已知中子的位置及其運動方向和運動距離，求其即將與哪個曲面相撞；

3）已知中子正在穿過某一曲面，求其即將進入哪個幾何體.

其中，任務1 已經(jīng)在2.1 節(jié)中進行了討論，此處不再贅述.

對于任務2，由于任意曲面都可以表示為f(x,y,z)=0，因此在已知中子當前位置r0和其運動方向rdir（單位向量）后，只需求解關于運動距離l的方程：

即可求得中子運動到該曲面上所需的距離，如果方程沒有正實數(shù)解，則中子不會與該曲面相撞，此時程序會將碰撞距離l記為正無窮.在模型空間的所有曲面中，l值最小的曲面即為即將與中子相撞的曲面.此處值得一提的是，假如最小相撞距離大于中子本次運動的距離，則中子不會與任何曲面碰撞，而是會在當前所處的幾何體內(nèi)發(fā)生核反應.

可見，這項任務的計算量同樣會受到遍歷的曲面數(shù)量的影響.在大多數(shù)情形下，所需遍歷的曲面數(shù)量不多，無需進行優(yōu)化，而當“中子位于基體內(nèi)，求其即將與哪個球體相撞”時，需要對所有球面進行遍歷，其計算量過大，需要借助虛擬網(wǎng)格加以簡化.簡化的方法與任務1 類似，程序首先可以根據(jù)中子的位置找出中子所在的單元格，然后僅遍歷該單元格內(nèi)的所有球面即可，如果這些球面均不滿足碰撞要求，則只需對下一個單元格內(nèi)的球面進行遍歷，如此循環(huán)直至找到即將碰撞的球面，這樣就有效減少了需要檢驗的曲面數(shù)量.

任務3 的邏輯則較為簡單.在中子穿過曲面離開某幾何體時，程序只需遍歷與該幾何體相鄰的所有幾何體，判斷中子是否在其內(nèi)部即可.但這一過程同樣存在計算量問題：當中子穿過球體表面進入基體時，程序需要判斷中子是否在基體空間內(nèi)，這需要對組成基體空間的大量曲面（這些曲面包括基體內(nèi)的所有球面）一一進行驗證，極為耗時.而實際上這一過程是沒有必要的，中子穿過球體表面離開球體時必然會進入基體.因此本文針對這種特殊情況進行了優(yōu)化，取消了其判斷過程，直接認定中子進入了基體，即可顯著減小計算量.

表4 展示了在56 個進程的純MP? 并行下，對于高溫氣冷堆單個燃料球模型（中子歷史為100 代，每代1 萬個）使用不同優(yōu)化水平時的計算耗時對比.可見對3 項任務的優(yōu)化均可使計算效率得到明顯提升.值得一提的是，如果使用真實網(wǎng)格方法對該模型進行加速，其計算耗時為65 s，可見其計算耗時高于虛擬網(wǎng)格方法.

Table 4 Comparison of the Time-Consumption of Single Pebble Model in HTR-PM Under Different Virtual Lattice Optimizations表4 不同虛擬網(wǎng)格優(yōu)化下HTR-PM 單個燃料球模型計算耗時對比

而對于HTR-PM 全堆模型（中子歷史為1 000 代，每代10 萬個），在相同并行規(guī)模下（10 個MP? 進程，每個進程56 個OMP 線程）與不進行優(yōu)化相比，優(yōu)化全部任務后其計算耗時可以從估算耗時約10 700 h 減小到539 s，計算效率大幅度提升，加速效果更為顯著.

3 結果分析

3.1 程序運行環(huán)境

本文所做的計算均在國家超級計算濟南中心的山河超級計算平臺開展，其采用?ntel Xeon 6258R 型號28 核CPU 處理器，主頻2.7 GHz，每個計算節(jié)點包含2 個CPU，其內(nèi)存為192 GB.使用的操作系統(tǒng)為Linux version 3.10.0-957.el7.x86_64，編譯所用的軟件為cmake/3.21.1，gcc/9.4.0，hdf5/1.8.13，mpich/3.4.2.

需要特別指出的是，表1 所示真實網(wǎng)格算例的內(nèi)存占用（213 GB）超過了山河超算平臺單個節(jié)點的內(nèi)存，該算例實際是在北京超算平臺測試的.而本文其他真實網(wǎng)格模型則略微縮小了網(wǎng)格規(guī)模，使其能夠在山河平臺運行.

3.2 計算結果分析

本節(jié)將介紹HTR-PM 和空間堆模型的定量化計算結果，以驗證優(yōu)化后程序的性能及其正確性.

有效增殖因數(shù)（keff）可以描述核反應堆的臨界程度，是衡量反應堆物理計算結果的重要指標.山東石島灣核電站于2021 年成功臨界，其臨界時keff=1.而本文按照其參數(shù)建立的HTR-PM 虛擬網(wǎng)格模型的keff計算結果為0.998 54±0.000 10，二者相差146 pcm.同時，如果使用OpenMC 程序本身的真實網(wǎng)格方法計算，則其keff=0.998 51±0.000 10，與虛擬網(wǎng)格僅相差3 pcm.可見，無論是與實驗值對比，還是與OpenMC 程序本身對比，HTR-PM 虛擬網(wǎng)格模型的結果都符合得很好.

而本文模擬的空間堆還處于理論階段，不存在相應實驗裝置.同時，文獻[17]并未對空間堆燃料顆粒進行顯式建模，而是將材料混合后建立了均勻燃料模型.研究表明，這樣的處理方法會使計算結果存在很大偏差[18-20]，不具有參考價值.綜上所述，本文的空間堆虛擬網(wǎng)格模型缺少實驗或計算結果作為驗證基準，因此其只能與OpenMC 自身的真實網(wǎng)格方法進行對比驗證.經(jīng)計算得到：虛擬網(wǎng)格模型的計算結果為keff=1.105 47±0.000 07，真實網(wǎng)格的結果為keff=1.105 51±0.00007，二者符合較好.

綜上所述，可以認為使用虛擬網(wǎng)格優(yōu)化后的程序具有良好的準確性.圖7 展示了HTR-PM 的堆芯功率分布，其中白球為石墨球，無功率.得益于2.2 節(jié)所述的虛擬網(wǎng)格數(shù)據(jù)統(tǒng)計方面的便利性，該功率分布結果可以精確到單個燃料球.

Fig.7 Power distribution of the HTR-PM core圖7 HTR-PM 堆芯功率分布

圖8 展示了空間堆模型的堆芯徑向功率分布，結果可精確到全堆尺寸的1/8000.

Fig.8 Power distribution of the space reactor core圖8 空間堆堆芯功率分布

3.3 并行性能分析

3.3.1 并行策略優(yōu)化

OpenMC 支持MP? 和OMP 混合并行，而在并行總核數(shù)相同的情況下，不同的MP? 進程數(shù)和OMP 線程數(shù)搭配得到的計算效率不同.因此，有必要尋找一種最優(yōu)的并行策略.表5 展示了同樣使用單節(jié)點56核并行的情況下，HTR-PM 虛擬網(wǎng)格模型在不同線程、進程數(shù)搭配時的計算時間，其中每種搭配分別測量了10 次結果，計算了其平均值、標準差.

Table 5 Comparison of the Time-Consumption of HTRPM Model with Different Parallel Strategies表5 不同并行策略下HTR-PM 模型計算耗時比較

由表5 可見，當進程數(shù)由1 變?yōu)? 時，計算效率存在明顯提升.這可能和山河超算的硬件條件有關，其單個56 核的節(jié)點實際包含了2 個28 核的CPU，因此相比單個進程的設置，2 個進程的設置可以簡化CPU 之間的通信，進而節(jié)省計算時間.此后隨著進程數(shù)的增加，計算效率存在先升后降的趨勢，這是因為MP? 并行雖然計算效率更高，但其通信也更為耗時，因此過少或過多的MP? 進程數(shù)設置都不利于獲得更高的計算效率.

考慮到相比于“2 個進程、28 個線程”的并行策略，其他設置對效率的提升不顯著，且更多的MP? 進程設置會大量占用節(jié)點內(nèi)存，不利于大規(guī)模數(shù)據(jù)統(tǒng)計，因此本文最終選取的并行策略是：每個節(jié)點2 個進程、每個進程28 個線程.

3.3.2 大規(guī)模并行效率分析

本文借助山河超算集群對HTR-PM 模型進行了大規(guī)模并行計算，測試了其在10 萬核規(guī)模下的并行性能.圖9 展示了HTR-PM 模型的強、弱可擴展性，其計算規(guī)模與2.2 節(jié)相同：對于強可擴展問題，每代中子數(shù)保持10 752 000 不變，總代數(shù)為20，并行規(guī)模為112～10 752 核（2～192 節(jié)點）；對于弱可擴展問題，每代中子數(shù)為核數(shù)的1 000 倍，總代數(shù)為20，并行規(guī)模為112～131 600 核（2～2 350 節(jié)點）.圖9 中統(tǒng)計的“計算時間”是指排除程序初始化以及寫輸出所用時間后的運行時間.對于每個數(shù)據(jù)點，本文均進行了10次測量，并選取了其平均值進行分析.

Fig.9 Parallel scalability of HTR-PM virtual lattice model圖9 HTR-PM 虛擬網(wǎng)格模型的并行可擴展性

當計算規(guī)模不隨并行規(guī)模變化時，程序在10 752核并行下與112 核并行相比具有83.4%的并行效率，其理想計算時間與實際計算時間分別為46 s 和55 s.而當計算規(guī)模與并行規(guī)模成正比時，程序在131 600核并行下與112 核并行相比具有83.1%的并行效率，其理想計算時間與實際計算時間分別為46 s 和56 s.可見程序在大規(guī)模并行時仍然具有較高的并行效率，其并行可擴展性良好.

此外注意到在弱擴展問題中，當并行規(guī)模由14 336核（256 節(jié)點）增大到28 672 核（512 節(jié)點）時，程序的并行效率突然下降，由95%降為84%，而此后并行效率隨并行規(guī)模的下降反而不明顯.目前分析這一現(xiàn)象可能由MP? 通信協(xié)議的變化引起.OpenMC 程序的通信耗時主要花費在進程間的中子源同步（synchronizing fission bank）上，在弱可擴展問題中，雖然每個進程的中子數(shù)是固定的，但進程間所需傳遞的中子數(shù)量仍會隨進程總數(shù)的增加而不斷增大[21].而MP? 在通信的信息量較大時，通信協(xié)議會隨數(shù)據(jù)量有所調(diào)整，從而導致計算耗時的階躍式增加[22].針對這一現(xiàn)象，在后續(xù)研究工作中將開展更加深入的分析和測試.

4 總結與展望

本文在開源MC 程序OpenMC 的基礎上，開發(fā)了適用于顆粒燃料的虛擬網(wǎng)格加速方法.通過構建虛擬網(wǎng)格并對中子輸運過程的3 項主要任務進行優(yōu)化，有效提高了程序的計算效率.該方法的計算耗時不但遠低于不使用網(wǎng)格的算例，且與OpenMC 此前的真實網(wǎng)格方法相比，計算耗時也更小.針對HTR-PM模型和空間堆模型，虛擬網(wǎng)格的計算時間分別減小到了真實網(wǎng)格的82%和46%.同時虛擬網(wǎng)格相較真實網(wǎng)格還具有內(nèi)存占用小、數(shù)據(jù)統(tǒng)計方便、并行可擴展性好的優(yōu)點.

在定量化結果方面，HTR-PM 模型的keff計算結果與石島灣核電站實驗值符合較好，因此可以認為程序具有良好的準確性.除此之外，本文還計算了HTR-PM 和空間堆的功率分布，并得益于虛擬網(wǎng)格方法數(shù)據(jù)統(tǒng)計的便利性，實現(xiàn)了HTR-PM 精確到單個燃料球的堆芯功率分布計算.

最后，本文對優(yōu)化后的程序進行了大規(guī)模并行測試，結果表明：對于強擴展問題，使用虛擬網(wǎng)格加速的HTR-PM 模型在10 752 核并行下與112 核并行相比，并行效率為83.4%；而對于弱擴展問題，程序在131 600 核并行下與112 核并行相比并行效率為83.1%，因此可以認為程序具有良好的并行可擴展性.

本文所建模型及測試更側(cè)重于對程序并行性能的探討，模型本身的計算并不需要10 萬核規(guī)模的并行.在此后的研究中，本文優(yōu)化的程序?qū)⒖梢杂糜谡嬲枰笠?guī)模并行的、更加復雜的模型中，例如帶有燃耗的HTR 模型.

同時，本文程序的并行性能仍存在優(yōu)化空間.例如進行燃耗計算時，OpenMC 在數(shù)據(jù)讀取和輸出部分的并行可擴展性較差，導致其占用了大量時間.這些問題均有待進一步的研究.

作者貢獻聲明：李睿涵負責算法優(yōu)化的實現(xiàn)與部分模型的建立以及論文的撰寫；侯葉凡負責部分模型的建立；李玉輝負責開展大規(guī)模并行計算；劉召遠指導并行計算工作以及修改論文；張海紅指導建模工作；梁金剛指導算法優(yōu)化以及修改論文.