基于天河二號超算的網格無關性及并行研究

徐藝峰，李 健，王 杰，范冰豐，王 鋼,3+

(1.中山大學 電子與信息工程學院，廣東 廣州 510000；2.中山大學 先進技術研究院，廣東 廣州 510275；3.中山大學 光電材料與技術國家重點實驗室，廣東 廣州 510275)

0 引 言

隨著計算機軟硬件的不斷優化和升級，求解問題的規模也日益增大。由于單個計算機的計算和存儲能力有限，而大型并行計算機可以實現多個處理部件與設備間的高效互聯，因此并行計算逐漸成為求解CFD問題的重要方法[3,4]。“天河二號”是我國自主研發的超級計算機平臺，坐落于中山大學的國家超級計算廣州中心[5,6]，計算和存儲能力世界領先[7,8]，是研究大規模并行計算的最佳選擇。

MOCVD(metal-organic chemical vapor deposition)生長過程復雜[9]，包括湍流層流模型、質量與熱傳輸、內部化學反應、三大守恒、工藝參數的設置與薄膜生長條件等[10]。實際生長要對每個因素進行全面而系統的研究因而成本巨大，利用計算流體力學(CFD)通過建立仿真模型來進行數值計算就體現出強大的優越性，并成為了國內外的研究熱點[11,12]。通過計算機進行各參數數值的計算，能夠得到比較詳細的資料，省時省力。

本文通過在“天河二號”上測試Fluent 15.0對MOCVD腔體模型的并行計算能力，探究了網格無關性，得出了最佳并行規模組合，使計算效率大大提高。

1 計算環境及并行計算原理

1.1 計算環境

2017年6月19日，最新的全球超級計算機500強榜單公布，來自中國的超級計算機“神威·太湖之光”和“天河二號”第三次攜手奪得前二。隨著行業和應用的需要，國家相關部門的重視，中國超級計算機正在飛速發展，在性能與應用上不斷邁步[13]。

位于無錫的“神威·太湖之光”由40個機柜、共160個超級節點組成，每個超級節點包含256個計算節點，每個節點裝有1個1.5 GHz、260核的SW26010眾核處理器和32GB DDR3內存，全系統總Linpack峰值浮點計算能力為為125.4359PFlops[14]。

本次測試中的“天河二號”的硬件系統包含有計算陣列、服務陣列、存儲子系統、互聯通信子系統、監控診斷子系統等，如圖1所示。計算陣列全系統包含125個計算機柜，每個計算機柜包含4個計算插框，每個計算插框包含16個計算刀片，每個計算刀片包含2個計算節點，因此計算陣列一共包含16 000個計算節點。每個計算節點包含：2個Intel(R)Xeon(R)CPU E5-2692 v2 @2.2 GHz 12核心處理器、3個Intel Xeon Phi 31S1P 57核心協處理計算卡、64 GB內存、高速互聯接口、2個以太網接口。一個E5-2692處理器Linpack峰值浮點計算能力為0.2112TFlops，一個Phi 31S1P協處理計算卡Linpack峰值浮點計算能力為1.003TFlops，總的峰值性能為3.4314 TFlops。因此計算陣列Linpack峰值浮點計算能力為54902.4TFlops，也就是每秒鐘5.49億億次[15]。

圖1 天河二號結構

軟件系統包含有系統操作、運行支撐等4種環境。當前操作系統版本為Red Hat Enterprise Linux Server release 6.2。應用開發環境包括串行編程語言、并行開發工具和并行編程模型。并行編程模型定義請參見文獻[16]。

1.2 Fluent并行計算

在天河二號超級計算機上的測試采用的是由安世亞太公司(Ansys公司的中國代理商)提供的測試版Ansys15.0。Fluent軟件是一個應用于模擬和分析復雜幾何區域內的流體流動與傳熱現象的專業軟件，該軟件是當今世界CFD仿真領域最為全面的軟件包之一[17]。

Fluent并行計算是用多處理器來計算大規模問題，計算可以在一臺機器上執行，也可以同時在多個不同機器上執行。并行處理主要目的是為了減少仿真時間，可以使用速度更快的機器，例如更快的CPU，內存，緩存以及CPU和內存之間的通信帶寬，也可以使用更快的互聯，例如較小的延遲和更高的帶寬，還可以使用更好的負載均衡，例如載荷均勻分布以及CPU的運算過程中不會空轉。Fluent并行計算的基本原理，如圖2所示。Fluent并行處理包括一臺主機進程，一系列的計算進程(處理器)以及計算進程之間的相互作用。Cortex功能是處理Fluent的用戶界面和基本的圖形功能。主機進程不存儲任何網格和計算數據，主要解釋來自cortex的命令，然后通過socket把這些命令發給計算進程0，計算進程0再將獲取的命令傳遞給其它計算進程。每個計算進程上都有一個MPI(message passing interface)，計算進程之間相互交換數據信息通過MP庫實現[18]。

圖2 Fluent并行計算結構原理

各計算進程計算之前，需要將整個三維流場區域分成多個計算區域。對網格進行分區域，需要選擇分區的方法，設置分區的數量，以及對網格分割的質量進行最優化。分區域的主要目的是為了平衡每個計算區域的網格單元數量、盡量減少各個計算區域邊界的表面積以及盡量減少計算區域的數量，確保每個處理器都有相同的負載以及各個計算區域之間數據交換能同時進行。本此測試采用METIS算法進行分區[19]，這種分區方法對計算區域的數量沒有限制，該算法將根據實際處理器的數量自動生成相同數量的計算區域。然后分配給相應的各計算進程去計算，由主機進程調度各個處理器的計算，在每一次全區的掃描過程中，由各個處理器完成計算并在邊界完成數據交換。最后由主機進程根據收斂準則進行判別，若收斂則計算結束，不收斂則繼續迭代求解。

2 MOCVD模型簡介及邊界條件設定

MOCVD的整個反應腔體為立式旋轉基座結構，如圖3所示。MOCVD的基本原理是以V族、Ⅵ族元素的氫化物和Ⅱ族、Ⅲ族元素的金屬有機化合物作為外延生長的源材料，經過熱分解反應與化合反應等一系列化學反應后，在襯底表面沉積出各種Ⅲ-V族、Ⅱ-Ⅵ族化合物半導體薄膜材料[20]。MOCVD腔體原幾何模型非常復雜，對原幾何模型在不影響流場流動的情況下進行簡化，如圖4所示，原上蓋的兩層進氣結構簡化為一層進氣入口，去掉了觀察窗等區域。

圖3 MOCVD反應腔體

圖4 MOCVD反應腔體簡化

本文選用SIMPLEC算法來求解控制方程組[12]，根據實際流動狀態設置CFD邊界條件如下：

(1)MO源入口和O源入口均設置為Velocity Inlet(速度入口)，O源混合氣體流量為1560 sccm，O2摩爾百分數為4.6%，MO摩爾百分數為0.22%，MFC1入口流量為147 sccm，MFC2入口流量為632 sccm，MFC3入口流量為747 sccm，MFC4入口流量為575 sccm，MFC5入口流量為173 sccm。

(2)出口設置為Pressure Outlet(壓力出口)為0 Pa，腔體內部壓強為10torr。

(3)襯底的轉速和溫度分別為750 rpm和723 k，反應堆內壁和外壁設置為絕熱和無滑移的。

3 網格無關性

從理論上來說，當網格數量越密時，模型越準確，但計算量也越大。由于目前的超級計算機是在Linux系統下操作，同時計算所保存的數據特別大，從超級計算機上傳和下載數據不便，因此有必要探索網格數量對結果的影響。網格無關性是指，當網格數量達到一定數量時，隨著網格數量再次增加，沉積率基本上變化不大的情況[21]。實驗中求沉積率的方法，如圖5所示，將轉盤表面放有硅片的區域沿著徑向分成9片圓環，標記為9個點，對每片圓環沉積率求得的平均值即每個點上的沉積率[22]。

本文模型中通過對網格數量作均勻加密處理，計算得出每組Case的9個點的沉積率數值，并比較每組沉積率的相關系數和平均沉積率的誤差，以探索網格無關性，具體情況如表1及圖6所示。

圖5 轉盤表面分區

CaseCells對Case01加密倍數Case0117萬1倍Case0268萬4倍Case03136萬8倍Case04408萬24倍Case05680萬40倍Case06901萬53倍Case071020萬60倍

圖6 不同網格數量的網格

通過計算加密后的網格，得出每個Case對應每個點的沉積率數值，其中沉積率的單位是μm/h，結果見表2。比較每個Case與Case07之間的相關系數(correl)，即可得出它們間變化趨勢的相似程度，而比較每個Case與Case07之間的誤差(error)，即可得出它們間平均沉積率數值上的差異程度。從表2中可以看出，對于相關系數，Case01到Case03，相關系數并沒有呈增加的趨勢，反而下降，表現出了一定的不穩定性，而Case04到Case06，相關系數不斷穩定接近于1，因此Case04到Case07的變化趨勢相似程度較高。對于誤差，Case01到Case03誤差較大，約為2%到3%，而Case04到Case06誤差縮小近一個數量級，達到約0.5%到0.7%，因此Case04到Case07的平均沉積率數值差異較小。綜合以上兩點說明，從Case04開始沉積率曲線已經趨于穩定，與Case07之間差異不大，可以認為此時已經達到網格無關性，為了能快速的進行大規模計算，沒必要繼續做加密處理，因此在進行數值計算過程中選用case04進行計算，而對于沒有達到網格無關性的Case，如Case01，由于其誤差也小于5%，因此工程上也可以應用它來進行調機，以快速得到結果。

表2 加密網格的計算結果

4 并行計算結果與分析

測試選取的模型是上面提到的兩個不同數量級網格單元，并且均是實驗常用的MOCVD反應腔體的三維網格模型，具體情況見表3。

表3 MOCVD模型的基本情況

并行計算時間：算法在并行機上求解問題時，從第一個任務進程開始執行到最后一個任務進程執行完畢所需的時間。包含CPU計算、并行開銷和算法輸入輸出所需時間總和。本次測試中，對于問題規模大的情況，如單核計算網格數量大的情況，計算時間是估算值。所有的計算時間均是指從開始迭代到迭代結束為止所用的時間，不包括啟動軟件以及讀入Case的時間，因為啟動軟件以及讀入Case的時間相對比較小，可以忽略不計。

并行加速比： Sn=T1/Tn， 即問題規模不變的情況下并行規模n時的并行加速比，并行規模通常為處理器核數，本次并行規模的核數選取依次為：1、2、4、8、12、16、32、64、128、160、240。

并行計算效率： En=Sn/n， 即問題規模不變的情況下，并行規模核數為n時的并行計算效率，假設單核迭代時，計算效率為1。

首先單獨考慮Case A的情況，即問題規模保持不變，如圖7所示，發現隨著并行規模核數的增加，加速比呈現先增加后減小的趨勢，在處理器核數為16時達到峰值，Case B也呈現出相同的趨勢，并且在處理器核數為128時達到峰值。加速比增加的原因是多核的并行計算，大大縮短了計算時間。加速比減小的原因是由于并行規模相對于問題規模較大，各個計算區域的信息數據交換非常消耗時間，進程之間數據交換的時間比進程計算的時間要長，從而導致整個并行計算時間相對增加，因此減少分界面以及計算區域的數量可以減少數據交換所需的時間。并且各個計算區域的信息數據交換方式還有待進一步的提高。

圖7 Case A的加速比曲線

對比Case A和Case B的加速比曲線，如圖8所示。當并行規模不變時，隨著問題規模的增加，加速比逐漸增加，因此對于問題規模比較大的情況，使用多核并行計算會大大縮短計算時間，從而降低數值模擬的周期。問題規模較小時，讀入Case的時間可以忽略不計，問題規模較大時，隨著并行規模加大，讀入Case的時間也相對的減少了。

圖8 Case B的加速比曲線

實驗結論：

(1)由于網格數量大小影響沉積率結果以及計算速度，因此在進行迭代計算之前，對于網格劃分時，要充分考慮到網格數量問題，實驗結果表明，對于MOCVD腔體模型，當網格數量達到408萬時，可以認為已經達到網格無關性。

(2)當問題規模保持不變時，隨著并行規模的增加，并行計算的加速比呈現先增加后減小的趨勢，計算效率也逐漸降低。對于17萬個六面體單元的MOCVD腔體模型，處理器核數為16時計算效率達到最高，加速比達到9，對于408萬個六面體單元的MOCVD腔體模型，處理器核數為128時計算效率達到最高，加速比達到45。

5 結束語

本次測試在國家超級計算廣州中心的平臺上進行，利用天河二號超級計算機的計算資源，由安世亞太公司提供的測試版軟件Ansys 15.0。測試驗證了網格無關性，并選取了17萬個六面體單元和408萬個六面體單元的兩個MOCVD腔體模型，對MOCVD反應腔體的溫度場及流場進行數值模擬，結果表明，前者在處理器核數為16時，加速比達到最高值9，后者在處理器核數為128時，加速比達到最高值45。

通過本次測試可以看出，對于MOCVD腔體數值模擬的大規模計算問題，由于需要計算的工況很多，而且網格數目巨大，利用普通的計算機很難完成，因此超級計算機并行計算成為了理想選擇。Fluent在天河二號上的應用嘗試，大大提高了仿真實驗的計算效率，適用與大規模并行計算，結果可為后續計算提供參考。