多線程并行在CHF一體化開發系統中的應用

明平洲 李治剛 潘俊杰 蘆韡

摘 要

CHF（Critical Heat Flux）一體化開發系統根據臨界熱流密度的實驗數據與子通道計算進行相互對比和統計學分析，最終對CHF經驗關系式進行修正和優化。由于CHF一體化開發系統的運行效率偏低，且反復調用子通道計算模塊，所以引入多線程并行編程技術來加速計算。通過OpenMP和POSIX線程兩種編程方案取得不同的加速效果，區別在于后者編程難度更大。19組實驗數據問題的CHF一體化開發系統運行效率被提升5倍左右，接近加速的理論上限。

關鍵詞

CHF經驗關系式;子通道計算;并行計算;多線程

中圖分類號： TM75 ? ? ? ? ? ? ? ? ?文獻標識碼： A

DOI：10.19694/j.cnki.issn2095-2457.2020.15.015

堆芯的熱工水力設計準則之一是在正常、運行瞬態及中等頻率事故工況下，堆芯不發生偏離泡核沸騰。當前計算偏離泡核沸騰比率（DNBR，Departure From Nucleate Boiling Ratio）的主流是采用子通道程序。目前采用子通道模型進行堆芯熱工水力分析的程序很多，世界上主要的核電設計公司均擁有各自的子通道程序，如法國阿海瑪公司開發的FLICA程序[1]，美國西北太平洋國家實驗室開發的COBRA程序[2]是典型代表，以上程序都已在工程上得到應用。中國核動力研究設計院設計所近期也研制了具備全堆芯計算能力的子通道程序CORTH。

對于堆芯的穩態工況和預計事故工況都要分別計算出最小的DNBR值。子通道計算中要準確地得到DNBR值，首先必須準確地計算出臨界熱流密度值CHF，記為qCHF。預測堆芯的qCHF較為復雜，堆芯軸向和徑向熱流密度的分布、棒的間隙以及流體、定位格架等多種因素都會對臨界熱流密度產生一定影響。經驗關系式在預測堆芯qCHF方面使用較多，例如著名的W-2關系式、W-3關系式和WRB-1公式等，它們對應于不同的工況和沸騰條件。查表方法近年來也是重要預測熱流密度的重要方案，它應用范圍更廣，便于對特殊情況進行修正[3]。經驗關系式涉及的大量實驗數據和計算數據往往通過多個步驟，多種程序工具進行處理，耗時較長。CHF一體化開發系統是一種偏向于經驗關系式方法的軟件系統，它對堆芯的qCHF進行預測，確定和優化CHF的經驗關系式。系統聯合實驗數據和子通道計算數據，縮短了數據處理周期。本文是使用多線程并行技術對研制的CHF一體化開發系統進行并行優化的研究，并歸納此類應用模式下的并行優化特點。

1 CHF一體化開發系統描述

1.1 計算原理

1.1.1 堆芯熱工水力學參數

壓水堆中燃料元件表面的最大熱流密度應小于臨界熱流密度，DNBR實質上是指利用合適的熱流量預測方法得到的冷卻劑通道中燃料元件表面某一點的臨界熱流量與該位置的實際熱流量比值[4]。

其中z為堆芯軸向的位置，子通道模型中徑向則被劃分為多個冷卻劑通道。DNBR的值沿軸向是變化的，其最小值（DNBR）min稱為該子通道的最小DNBR。在壓水堆熱工設計中，把整個堆芯（所有子通道）最小DNBR不小于某個值作為熱工設計準則之一。

研制的CHF一體化開發系統構造了類似于已公開發表內容[5]的關系式。關系式所需的熱工水力本地參數可以由子通道程序計算得出。現階段關系式適用于壓水堆的熱工水力分析。

1.1.2 非線性回歸

CHF經驗關系式的準確性主要與關系式的各項系數相關。根據實驗數據和子通道計算產生的中間結果，可以采用非線性回歸方法進行系數的確定、修正或者優化。

常用的非線性回歸方法有：共軛梯度法、Quasi-Newton法等。Quasi-Newton法是求解非線性優化問題最有效的方法之一，本文選用Quasi-Newton法擬合CHF經驗關系式，該方法在系統使用的第三方工具Scilab[6]內的調用格式為：

[x，err]=datafit（g，z，x0，qn）（3）

其中，g表示擬合目標函數，z為存放測量數據的矩陣，err為誤差返回值，x0為系數初始值，x為擬合所得系數。qn則制定使用的數據擬合方法，這里選定為Quasi-Newton法。

1.2 系統的程序結構

CHF一體化開發系統目前采用結構化程序設計，采用C++編程語言書寫控制流程和數據存儲，同時預留接口，對多個第三方程序進行集成。整個程序具備模塊化編程的特點，各個重要功能可以進行合理替換。

上圖是整個CHF一體化開發系統的控制流程框圖，其中子通道程序計算和統計學檢驗依賴于第三方程序。圖中子通道計算模塊需要對多組數據調用外部的子通道程序進行計算處理，也是本文中需要進行并行編程應用的計算區域。

1.2.1 第三方程序調用

子通道程序的建模和計算采用中國核動力研究設計院自主研制的CORTH子通道程序。為了與CHF一體化開發系統較好的耦合，對CORTH程序進行代碼級別的集成;CHF一體化開發系統調用解釋型腳本來耦合SCILAB的接口，實現非線性回歸等功能;調用第三方畫圖軟件Gnuplot實現對關鍵參數分布、數據趨勢進行圖像化輸出

1.2.2 模塊功能

CHF實驗數據預處理：對CHF實驗獲得的大量數據點進行預處理，如去除實驗壞點等，形成開發數據庫和驗證數據庫（對比驗證）。

實驗數據子通道建模和計算：批量調用子通道程序，對CHF實驗數據進行計算和預測。

CHF經驗關系式系數求解：根據CHF開發數據庫對應的子通道計算結果確定并優化CHF關系式的形式，并求解得到CHF關系式的各項系數。

建立M/P數據庫：將前一階段計算得到的CHF關系式嵌入至子通道程序中，再次對CHF開發數據庫調用子通道程序計算，提取臨界位置的DNBR等參數，形成M/P數據庫。

確定DNBR限值：對M/P數據庫進行統計學分析和檢驗，計算DNBR限值。

CHF經驗關系式的驗證和評價：采用CHF驗證數據庫，調用已嵌入新CHF關系式的子通道程序進行建模和計算，獲得驗證數據庫的M/P數據庫，結合開發數據庫的M/P數據庫，評估新CHF關系式的預測率并對熱工參數的變化趨勢等進行評價。

2 并行編程方案描述

根據程序結構的描述，多組實驗數據涉及的子通道程序建模和計算具備并行度。它們調用第三方程序CORTH完成子通道的計算相互獨立，可以使用多線程并行編程技術來對計算加速。

CHF一體化開發系統占用的內存在2G以下，運行于工作站服務器，多線程并行編程屬于共享式內存并行策略，可以按照如下兩類算法形式來實施：

2.1 OpenMP并行

基于編譯器提供的OpenMP并行編程[7]方法，對CHF一體化開發系統中的全局變量進行轉換，降低數據相關性，然后添加OpenMP指導語句。

OpenMP并行編程使用的是編譯器指導語句，開發人員一般不需要進行顯式確定調度策略，由于實驗數據CHFExpData與數據組數i相關，不均勻性會帶來一定負載均衡問題。

2.2 POSIX線程并行

類UNIX計算服務器提供了POSIX標準線程庫[8]，相比于OpenMP并行編程更貼近操作系統底層。CHF一體化開發系統在初始階段生成多個線程，構成池式控制結構，通過互斥開關來使用這些線程完成并行計算。

POSIX多線程并行充分利用了線程細粒度運行任務的特點，以單次調用子通道計算為基本任務，避免了負載均衡等問題，但線程池方案存在資源互斥，編程難度大，且輸出結果的次序也會影響后續的統計分析。

3 實驗結果與分析

CHF一體化開發系統及其并行優化的初始版本目前已經研制完畢，這里選取HP Z800計算服務器作為運行環境。

整個實驗涉及的CHF實驗數據分為開發實驗數據和驗證實驗數據兩類，其中開發實驗數據共19組，下表是每組實驗數據對應的CHF實驗樣本點個數。

從輸入來看，19組實驗數據進行建模和計算存在著不均勻性，多線程并行的負載均衡需要進行考慮。按照整個CHF一體化開發系統的運行過程，完成整個計算需要對這19組實驗數據前后迭代四次。實驗結果先統計單次迭代運行效率，即多線程并行一次處理19組實驗數據的性能。整個CHF一體化開發系統的運行時間（迭代四次）對比將在最后給出。

3.1 OpenMP并行編程

考慮到共有19組實驗數據，而OpenMP并行僅僅針對外層循環實驗組數進行，所以這里選取4個，8個和10個線程三種情況進行并行實驗，觀察OpenMP并行編程情況下的收益。

其中4個線程對應并行總共調用CORTH子通道次數為五次（=5），8個線程對應并行調用三次，10個線程對應并行調用兩次。在OpenMP并行編程的情況下，計算時間的波動較為明顯，這是由于CHF開發實驗數據的不均勻性造成，OpenMP的默認調度方式無法做到負載均衡。

考慮10個線程運行情況，在OpenMP的編譯指導語句中引入動態調度。圖2是10個線程情況下OpenMP默認調度方式與動態調度方式的效率對比。由圖中所示，動態調度方式比起默認調度方式有了一定改進，各個線程的負載相對均衡，但運行時間比起默認調度方式反倒有所惡化。這也從側面說明動態調度方式在本計算例題情況下不是最優選擇。

3.2 POSIX線程并行編程

利用POSIX標準線程庫開發，對于每組實驗數據內的單個實驗樣本點均從線程池中獲取線程進行運行，當線程池中沒有線程存在時則阻塞掛起，這種并行編程方式粒度較細。

為了兼顧通用性，OpenMP并行編程往往針對外層循環。CHF一體化開發系統中內層循環的任務數與外層循環的序號相關，具有一定隨機性，所以直接采用POSIX線程的并行編程方式能夠帶來更多的收益。

對比表6和表7可以發現，OpenMP并行編程的編程難度較小，但從實驗結果來看存在負載均衡問題;另一種POSIX線程的編程難度較大，但從實驗結果來看并行收益更高。

此外，完整運行CHF一體化開發系統對19組實驗數據進行處理。串行計算時除開調用子通道計算的內容，每次完成19組實驗數據的處理耗時約320秒左右，整個計算時間為1687秒左右，則可并行區域占總計算的百分比為80%左右。根據阿姆達爾定律，并行優化之后帶來的加速比上限約為5倍左右。實際工程實現之后接近這一理論值，證明可并行區域具有天然并行的特性，但各個任務之間存在著負載不均衡的缺點，工程編程實現時通過POSIX線程的策略能避免此缺點。

4 結束語

本文對中國核動力研究設計院的CHF一體化開發系統進行了描述。引入了兩種并行編程方案，在調用子通道計算的過程中改善其計算效率。形成如下結論：

1）此問題域中，OpenMP并行編程實現簡便，但存在負載均衡問題，POSIX線程池的方法實現困難，但負載較為平均。

2）CHF一體化開發系統的各個子任務與實驗數據的規模相關，所使用的19組實驗數據可并行區域帶來的加速比上限約為5倍左右，目前的并行編程工程實現接近這一理論值，改善了CHF一體化開發系統的運行效率，為多次執行CHF實驗數據的分析和處理奠定了一定基礎。

后續該系統將應用于CHF實驗數據的評價和篩選，且數據之間的相關性以及CHF關系式優化等內容將被進一步研究。

參考文獻

[1]Fillion P，Chanoine A，Dellacherie S，et al.FLICA-OVAP：A new platform for core thermal–hydraulic studies[J].Nuclear Engineering and Design，2011，241（11）：4348- 4358.

[2]Kelly J E，Loomis J N，Wolf L，et al.LWR core thermal-hydraulic analysis： assessment and comparison of the range of applicability of the codes COBRA IIIC/MIT and COBRA IV-I[R].OSTI-5493728，1978.

[3]于平安，朱瑞安，喻真烷.核反應堆熱工分析，上海交通大學出版社[M].2003.

[4]蘇光輝，秋穗正，等.核動力系統熱工水力計算方法，清華大學出版社[M].2013.

[5]FRAMATOME.FC-FRAMATOME Critical Heat Flux Correlation for AFA-2G and AFA-3G Fuel Assemblies-Topical Report （EPD/DC-454，Rev.A），December， 2000.

[6]胡寶鋼，趙星.科學計算自有軟件—SCILAB教程[M].清華大學出版社，2003.

[7]Salvini S.Unlocking the Power of OpenMP[C].Proc. Fifth European Workshop OpenMP （EWOMP03）， invited lecture，2003.

[8]March V，Zhang R，See S，et al.Survey on Parallel Programming Model[C].Network and Parallel Computing，2008.