數值堆結構力學高精細模擬的并行實現和優(yōu)化技術

邢龍岳，王昭順，蔣章程，高付海，高 岳，馮仰德，王 玨

(1.北京科技大學，北京 100083；2.中國原子能科學研究院，北京 102413；3.中國科學院 計算機網絡信息中心，北京 100190)

核能作為一種清潔、安全和經濟的新型能源，是解決能源危機與環(huán)境問題的有效途徑，制約核電發(fā)展的關鍵是核反應堆的安全問題。在過去的20年中，歐洲、美國和中國的研究為虛擬反應堆的發(fā)展提供了重要信息[1-4]。由于反應堆堆芯組件具有多樣性和復雜性，實現全堆高精細結構力學模擬所需的計算資源和存儲資源遠非服務器和小型集群計算系統(tǒng)所能滿足。因此依托高性能計算技術，實現高精細的數值模擬就成了研究反應堆安全性的一條不可或缺的途徑。而我國異構體系結構的高性能計算系統(tǒng)目前居世界領先地位，且E級規(guī)模(每秒10×1018次數學計算)先進計算系統(tǒng)也即將來臨。基于超級計算機和先進建模技術的數值堆(CVR)結構力學模擬是在此背景下誕生的一個挑戰(zhàn)性的研究領域。CVR結構力學模擬是研究堆芯組件的靜、動態(tài)結構變形及應力應變。

異構計算系統(tǒng)中GPU處理器主要可分為NVIDIA GPU和AMD GPU。AMD GPU處理器實現了一種并行微體系結構，不僅為計算機圖形應用程序還為通用數據并行應用程序提供了一個極好的平臺。大多數需要高帶寬和計算密集型或數據密集型應用程序都可在AMD GPU處理器上高效運行[5]。如何對AMD GPU在異構計算系統(tǒng)進行核心優(yōu)化以獲得高性能對程序員來說是一個挑戰(zhàn)。高精細模擬采用三維空間有限元分析方法[6]，網格劃分是有限元分析的重要環(huán)節(jié)，是實現高精細數值模擬的基礎。對非結構化網格生成的研究，相關學者已進行了數十年，盡管已提出了多種方法[7-8]來提高網格生成的效率，但是對于前沿應用程序的模擬而言，生成大規(guī)模網格的時間仍難以接受。本文利用堆芯組件的同構性，通過使用單個組件或局部組件的網格文件自動生成全堆芯網格，從而極大地提高大規(guī)模(百億、千億級別)網格的劃分效率。通過數值模型推導將有限元網格轉換為稀疏線性剛度矩陣方程組Ku=r。由于剛度矩陣K具有對稱正定的特征，因此可使用大類算法Krylov子空間迭代法[9-10]中的共軛梯度(CG)[11]小類算法，并采用Jacobi[12]預處理進行高效數值求解。

1 CVR全堆芯大規(guī)模網格自動生成

CVR全堆芯大規(guī)模網格自動生成方法支持有限元分析前處理的完整功能，包括CVR幾何建模、有限元分析預處理、大規(guī)模并行網格生成3個階段。圖1為CVR全堆芯大規(guī)模網格自動生成技術流程圖。

圖1 CVR全堆芯大規(guī)模網格自動生成技術流程圖Fig.1 Flow chart of CVR core large-scale mesh automatic generation technology

1.1 數值堆幾何建模

采用反應堆全堆芯結構的標準尺寸，對各類型組件分別進行一次幾何建模，并將幾何模型的底面中心點坐標設置為(0，0，0)。以中國實驗快堆(CEFR)為例，采用常用的CAD模型文件格式(STEP、IGES、Brep、GEO等[13])，分別對其中的燃料組件、乏燃料組件、鋼屏蔽組件、硼屏蔽組件、安全棒組件、調節(jié)棒組件、補償棒組件、中子源組件進行幾何建模。圖2為CEFR全堆芯組件分布，表1列出CEFR全堆芯組件說明。

表1 CEFR全堆芯組件說明Table 1 CEFR whole core assembly description

圖2 CEFR全堆芯組件分布Fig.2 CEFR whole core assembly distribution

CEFR全堆芯712根組件的分布方式為多層正六邊形陣列，且有1個中心組件，其二維中心點坐標為坐標軸原點(0，0)，相鄰組件的中心距為61 mm。對每根組件進行編號后，計算求得每根組件的中心點坐標，并建立組件編號與組件中心點坐標的映射關系。采用上述方法確定CEFR全堆芯712根組件中心點坐標后，根據組件類型對每根組件進行編號。通過輸入所需組件的編號集合，即可導出所需的反應堆各類型組件的幾何模型文件及幾何位置數據。圖3為CVR幾何建模實現流程。

圖3 CVR幾何建模實現流程Fig.3 CVR geometric modeling flow chart

1.2 有限元分析預處理

調用三維有限元網格生成器Gmsh[14]，對1.1節(jié)中導出的幾何模型進行并行網格劃分，輸出不同類型組件的網格文件。通過自編寫的接口程序調用ParMetis/Metis庫[15]和Zoltan庫[16]，對網格數據進行區(qū)域分解與負載均衡處理，劃分得到各類型單個組件的多個分區(qū)的網格數據，并將所得的網格數據寫入不同的網格文件中，用于實現后續(xù)有限元求解程序的并行求解。

對各類型組件的網格子文件集添加有限元求解所必需的控制參數包括：邊界條件、初始條件、材料屬性。其中，邊界條件包括位移邊界條件和力邊界條件；初始條件包括初始位移、初始速度、初始加速度、初始溫度和初始狀態(tài)變量；材料屬性包括密度、彈性模量、泊松比和熱膨脹系數。圖4為有限元分析預處理實現流程。對于相同類型的組件，其材料屬性是相同的，但是邊界條件和初始條件可能不同，如組件上溫度的分布情況。針對這種情況，需進行兩次物理參數的添加：第1次是對基準組件網格添加材料屬性等恒定不變的物理參數(當組件底面都被固定時，不同組件的位移邊界條件也是相同的，因此相關物理參數也可在第1次就進行添加)；第2次是在進行1.3節(jié)所述的大規(guī)模并行網格生成時，實現不同組件間存在差異的物理參數的自動分配和添加。

圖4 有限元分析預處理實現流程Fig.4 Finite element analysis preprocessing flow chart

1.3 大規(guī)模并行網格生成

使用傳統(tǒng)的網格劃分軟件對全堆芯組件進行大規(guī)模網格劃分，需通過網格生成算法進行節(jié)點生成、網格構造、網格排序、網格質量優(yōu)化等數據操作，常用的網格生成算法時空復雜度列于表2，其中，n為單根組件的網格數量。以百億網格文件為例，相關網格數據的大小接近1 TB，對計算機的計算能力和存儲能力要求極高。目前常用的網格劃分工具，最大網格劃分速率基本維持在8 000～9 000 s-1。假設劃分速率為10 000 s-1，完成百億網格數量的劃分任務需要277.8 h。為了解決這個問題，本文設計并實現了大規(guī)模并行網格生成程序。該程序基于同類型組件間具有同構性的特點，不需對所有的組件進行整體網格生成、構造、排序和質量優(yōu)化，只需對于其中一個基準組件進行網格生成，再通過對基準組件的網格進行網格的復制和平移，即可生成多組件乃至全堆芯組件的網格，大幅減少與網格生成相關的大量運算。同類型組件間網格的坐標空間變換方式有平移、旋轉、鏡像等(圖5)。

圖5 同類型組件網格幾何變換示意圖Fig.5 Geometric transformation of the same type assembly mesh

表2 常用網格生成算法的時空復雜度Table 2 Time and space complexity of common mesh generation algorithm

針對全堆芯組件進行大規(guī)模網格生成，分別采用組件間網格串行生成方法、組件間網格并行生成方法及本文提出的大規(guī)模網格并行生成(組件間網格并行復制)算法，各自所需的時空復雜度列于表3。其中，m為堆芯中所需生成網格的組件數量。可看出，組件間網格并行復制方法的時空復雜度是最優(yōu)的。

表3 組件間網格生成算法的時空復雜度Table 3 Time and space complexity of mesh generation algorithm between components

網格數據主要包括網格節(jié)點數據、網格單元數據和控制參數。網格節(jié)點數據用于說明每個網格節(jié)點的坐標；網格單元數據用于說明每個網格單元由哪些網格節(jié)點構成；控制參數用于說明影響問題求解的各參數選項。得到添加了控制參數(包括物理參數、與求解相關的參數、計算結果輸出路徑等)的網格子文件集后，通過大規(guī)模并行網格生成程序讀取單個組件的網格子文件集及幾何位置數據，并進行網格節(jié)點坐標的移動、轉換，網格節(jié)點及網格單元編號的變換，由于不同組件的控制參數可能存在差異，因此還需進行相關控制參數的變換，輸出所需的CEFR全堆芯或多組件的網格數據[17]。其中，在進行控制參數的變換前，需在腳本文件中按不同組件的編號寫入不同組件間具有差異性的控制參數，在進行網格復制轉換的同時，程序會將不同參數分配給指定的組件，即進行第2次物理參數的添加。

圖6 大規(guī)模并行網格生成程序流程Fig.6 Large-scale parallel mesh generator flow chart

2 大規(guī)模稀疏剛度矩陣線性方程組并行求解算法設計

本文采用ANSYS對簡單立方體靜力算例的用時估算(圖7)，可得出百億網格量的求解用時將近116 d。另外，由于在大規(guī)模計算中系數矩陣過于龐大，采用Jacobi預處理的方法進行加速。但隨著每個并行區(qū)域中剛度矩陣的規(guī)模逐漸增加，預處理的迭代次數會明顯增加，收斂速度變得十分緩慢。這不但增加了每個計算核心所處理的浮點計算量(FLOPs)，而且增加了更多的計算核心資源使用，導致產生大量的通信時間。

圖7 ANSYS立方體靜力算例用時Fig.7 Time for ANSYS to solve cube static example

2.1 異構平臺上的并行計算

GPU提供了一種擁有強大浮點計算能力和極高存儲器帶寬的低成本、低功耗通用的CPU+GPU異構計算平臺[18-19]。其特點可支持高度并行、多線程并擁有眾核處理器。以天河2號超級計算機Intel(R) Xeon(R) CPU E5-2692v2為例，單CPU雙精度浮點計算能力為211.2 GFLOPs，數據傳輸能力為50～60 GB/s，而單AMD GPU(Vega series)雙精度浮點計算能力則可達到6.25 TFLOPs，是前者的28倍以上，且與顯存之間的傳輸最高可達1 TB/s。

CPU和GPU浮點計算能力差異顯著的原因是因為GPU是專門為圖形渲染而設計的，而圖形渲染計算特點為計算密集型、高并行度計算。因此GPU更致力于數據處理而不是數據緩存和流量控制，為了達到高密度計算的要求，GPU芯片上有更多的算術邏輯單元。CPU與GPU結構的差別如圖8所示。具體來說，GPU更適合解決大規(guī)模并行計算的問題，即適用于同一程序同時對多組具有較高計算密度的數據進行運算。因為同樣的程序代碼在每個數據元素上執(zhí)行時，對精密流量控制的要求較低，尤其適合大規(guī)模稀疏剛度矩陣向量乘。

圖8 CPU與GPU結構的差別Fig.8 Difference between CPU and GPU structures

2.2 并行方案

在包含AMD GPU的異構計算系統(tǒng)中，并行計算過程中允許使用大量Thread。可根據并行編程模型為每個數據或數據集設置一個Thread，并以AMD GPU中每個計算單元中可執(zhí)行的最大Thread數量來劃分，從而獲得最優(yōu)的并行計算性能。為此需對實際的應用特征進行分析并給出方案。

1) 密集內存訪問

CEFR靜力數值模型構造的稀疏剛度矩陣線性方程組Ku=r在多個計算節(jié)點間并行計算時，各計算單元需要將內存中的數據讀入到顯存，并將計算的結果寫入到內存。并行規(guī)模越大，這種訪存操作所消耗的時間將會大大超過計算所消耗的時間，因此在大規(guī)模計算時，訪存的時間開銷將成為提高計算效率的最大瓶頸。解決此類問題，可采用內存與顯存之間的數據分塊和雙緩沖技術，以減少訪存的時間。

2) 同步要求

計算規(guī)模越大，迭代次數就會越多，這是Jacobi迭代法特質決定的。在Jacobi迭代中，除第一次計算外，每次迭代需各計算單元同步并更新上一次迭代的數據，這需要AMD GPU結束內部所有線程才能實現全局同步。為了減少同步開銷，首先要減少迭代次數，可采用分塊Jacobi方法，這種方法對剛度矩陣的構造是有要求的，因此需在預處理組裝剛度矩陣時進行額外的處理。該方法能大幅減少計算的迭代次數；其次只為每個MPI進程分配1個塊，這樣可有效減少主機端與設備端額外的數據傳輸，并均衡各進程間的負載，降低對同步的時間開銷。

在求解Ku=r時，每個計算單元之間數據是相互依賴的，同一個單指令多數據(SIMD)單元執(zhí)行的線程在所有的計算節(jié)點中會尋址連續(xù)的內存塊，因此將每個進程中的向量計算數據連續(xù)存儲，這樣相鄰線程執(zhí)行時只需對上一個線程的數據索引加1就可得到當前線程的向量計算數據，從而提高計算性能。

2.3 AMD GPU中的并行優(yōu)化設計

有限元方法得出的結果是近似值，因此要實現反應堆堆芯結構力學高精細模擬，則網格量需達到十分龐大的規(guī)模，這也將生成巨大的Ku=r方程組，所以面向AMD GPU，本文基于Krylov子空間迭代法，提出了幾種并行優(yōu)化設計。

1) HIP編程模型中的線程部署

在AMD GPU上編程的語言稱為HIP，它構建了可遷移的運行時體系結構，這個體系結構允許HIP程序可動態(tài)地在多個平臺進行切換。HIP是顯式編程模型，需要在程序中寫出詳細的并行控制語句，包括數據傳輸、核函數的啟動等。所謂核函數Kernel就是運行在AMD GPU上的函數。在CPU上運行部分稱為主機端(host，執(zhí)行管理+啟動)，在AMD GPU上運行部分為設備端(device，執(zhí)行計算)。程序流程是在主機端串行執(zhí)行，當到設備端執(zhí)行時則首先將需計算的數據通過hipMemcpy傳遞給AMD GPU(即將CPU存儲中的內容傳遞給AMD GPU的顯存)，然后調用核函數啟動函數hipLaunchKernelGGL使AMD GPU開始計算，Kernel端計算完成后再用hipMemcpy將計算的數據從AMD GPU拷貝回CPU。其中hipMemcpy是阻塞的，數據拷貝后才可進行后面的程序，hipLaunchKernelGGL是非阻塞的，執(zhí)行后程序繼續(xù)往后執(zhí)行，但在Kernel未計算完成之前，最后一個hipMemcpy不會開始。

在HIP中對設備端Kernel使用三級并行編程模型對其并行化進行線程部署，如圖9所示。按照3D布局部署Grid、Block和Thread，以平衡每個矩陣維度上計算域的分解。相鄰的Block共享位于邊界上的數據，且Block中的每個Thread負責在每次迭代中更新數據，而Threads實際是以Wavefront執(zhí)行的，為實現最優(yōu)并行化，映射時將Block中的Thread數量設置為Wavefront執(zhí)行數的倍數，這樣可達到更好的計算效果。

圖9 并行化策略的線程部署Fig.9 Thread deployment of parallelization strategy

2) 數據預取

在執(zhí)行核函數時，主機端計算是阻塞的，數據拷貝完成后才可進行后面的程序。若多個核函數同時執(zhí)行則需很大的時間開銷。若將一個核函數放入一個邏輯隊列中，多個向量計算分別建立多個核函數，就可在AMD GPU中采用流水線并行實現數據預取，如圖10所示。這種方案可充分利用CPU與AMD GPU的計算資源，提高計算效率。

提素質、強能力，創(chuàng)新人才工作新機制。一直以來西南銅業(yè)公司黨委匡正選人用人風氣，完善選人用人機制，堅持嚴管和厚愛結合、激勵和約束并重，優(yōu)化干部考核評價機制。抓好各級干部的領導、監(jiān)管、屬地和執(zhí)行責任，提升黨政共融，黨務干部要抓安全、懂生產、謀經營，行政干部要講政治、會黨建、抓文化。完善激勵機制和容錯糾錯機制。開展好選人用人制度修訂、培訓及考核、中層管理人員畫像、后備干部選拔等工作。堅持黨管人才原則，推進員工人才成長體系建設，明晰各條人才通道的晉升、考核、退出、流動機制，為培養(yǎng)公司“素質優(yōu)良、結構合理、隊伍穩(wěn)定”的員工隊伍打下扎實基礎。

圖10 流水線并行中的數據預取Fig.10 Data prefetching in pipeline parallel

3) 數據對齊

在對齊的邊界上創(chuàng)建數據可提高CPU上的SIMD性能。若向量數據未與緩存線大小對齊，編譯器必須額外處理未從緩存線邊界開始的數據部分的計算，從而導致性能損失。因此，需根據向量寄存器上緩存線的大小來設置動態(tài)分配的數據對齊字節(jié)數，圖11為512位向量寄存器上的數據未對齊。

圖11 512位向量寄存器上的數據未對齊Fig.11 Data alignment on 512 bites vector register

3 性能與驗證

3.1 測試算例與測試目標

對CEFR堆芯組件與全堆芯組件進行百億網格劃分和形變位移的測試，測試算例邊界載荷為下過度接口處鉸支，約束所有平動自由度，燃料組件添加溫度載荷的變化區(qū)間為300～500 ℃。表4列出了CEFR堆芯燃料組件物理參數設置。測試目標包括：1) 采用這種方法的可行性和正確性，通過國際原子能機構(IAEA)標準例題的解析解進行對比；2) 采用這種方法后對反應堆結構力學大規(guī)模并行模擬的加速效果的提升；3) 多節(jié)點集群上優(yōu)化程序的并行可擴展性。

表4 CEFR堆芯燃料組件物理參數設置Table 4 Physical parameter setting of CEFR core fuel assembly

3.2 AMD GPU的加速性能

將CEFR單燃料組件劃分為1 947 584個有限元網格，通過有限元力學數值模型將其轉化為稀疏剛度矩陣線性方程組的形式Ku=r，然后使用Krylov子空間迭代法中的CG方法和Jacobi預處理，結合在AMD GPU中的并行方法對其進行計算。CPU計算核心從1個CPU Core擴展到32個，并使用1～4個AMD GPU與之結合。實驗結果展示了在其并行算法下使用AMD GPU對計算性能的加速效果。圖12為單燃料組件形變位移分布和CPU+AMD GPU加速效果，可看出，當1個節(jié)點中使用4個CPU Core和4個AMD GPU時，并行計算效率最高，而在32個CPU Core和4個AMD GPU的計算中執(zhí)行時間最短，但并行強可擴展性不佳。

圖12 單燃料組件形變位移分布(a)和CPU+AMD GPU加速效果(b)Fig.12 Deformation displacement distribution of single fuel assembly (a) and CPU+ADM GPU acceleration (b)

3.3 功能與規(guī)模測試

1) CEFR全堆芯百億網格劃分

首先分別為燃料組件、乏燃料組件、硼屏蔽組件、事故棒組件建立CAD幾何模型；然后針對不同類型組件劃分并生成單組件網格文件；再將各種類型組件之間的相對位置信息寫入配置文件；最后自動生成全堆芯大規(guī)模網格文件。CEFR全堆芯百億網格生成實現流程如圖13所示。

圖13 CEFR全堆芯百億網格生成Fig.13 CEFR whole core 10 billion meshes generation

反應堆全堆芯大規(guī)模網格自動生成方法還可將生成后的網格均勻地劃分為多個網格區(qū)域，以8組件為1組的網格區(qū)域分解與負載均衡效果如圖14所示，其中不同的色塊代表不同的劃分區(qū)域。將單組件劃分到1 432萬網格，每根組件可劃分96個分區(qū)，可實現102億的全堆芯網格量，且最高可劃分137億，數據存儲量為1.067 TB，網格數據由節(jié)點位置數據和點、線、面、體幾何模型數據組合而成，包含坐標(double)、編號ID(int)、類型ID(int)、物理ID(int)、點ID(int)。全堆芯網格文件以二進制文件存儲，網格文件數量等同于分區(qū)數量，并可轉換為vtk格式使用Paraview可視化。

圖14 網格區(qū)域分解(a)與負載均衡(b)Fig.14 Mesh domain decomposition (a) and load balancing (b)

2) CEFR全堆芯大規(guī)模靜力分析

首先采用六角形套管單組件的靜力學熱形變位移計算模型，與IAEA標準例題的解析解對比，驗證靜力學計算結果的正確性(表5)。

表5 IAEA標準例題靜力學驗證算例Table 5 Static verification example of IAEA standard example

通過對程序計算結果和理論值的對比，分別給出了六角形套管單組件在給定線性溫度梯度下彎曲變形量(表6)和軸向伸長量(表7)在高度為4 000、3 000和2 500 mm處的偏差情況。組件彎曲變形量和軸向伸長量的數值解與理論值偏差較小，驗證了計算結果的正確性，可視化效果如圖15所示。

圖15 六角形套管單組件彎曲變形量(a)和軸向伸長量(b)模擬結果Fig.15 Simulation result of bending deformation displacement (a) and axial elongation (b) of single hexagon casing assembly

表6 六角形套管單組件在給定線性溫度梯度下的彎曲變形量Table 6 Bending deformation of single hexagon casing assembly under given linear temperature gradient

表7 六角形套管單組件在給定線性溫度梯度下的軸向伸長量Table 7 Axial elongation of single hexagon casing assembly under given linear temperature gradient

結合異構計算系統(tǒng)的AMD GPU，對整個CEFR全堆芯組件劃分為0.54～102億網格的計算模型進行靜力分析模擬，算例列于表8，靜力模擬結果的形變位移分布如圖16所示。

表8 CEFR全堆芯組件的靜力學算例Table 8 Static example of CEFR whole core assembly

圖16 CEFR全堆芯形變位移分布Fig.16 Deformation displacement distribution of CEFR whole core

對具有0.54億網格量的CEFR全堆芯組件計算模型進行靜力模擬，計算節(jié)點從64個節(jié)點擴展到512個節(jié)點，每個節(jié)點使用258個計算核心，最高使用到132 096個。CPU+AMD GPU加速比與效率對比和并行弱可擴展性如圖17所示。在異構AMD GPU計算系統(tǒng)上與CPU的版本相比，提高了30%以上的計算性能。另外，以9 504計算核心(CPU+AMD GPU)對應1.584億網格量為基準擴展至615 168計算核心對應102.528億網格量，測試出并行弱可擴展性達到53%。

圖17 CPU+AMD GPU加速比與效率對比和并行弱可擴展性Fig.17 CPU+AMD GPU speedup, efficiency and parallel weak scalability

4 結論

本文利用CEFR堆芯組件的同構性，通過使用單個組件或局部組件的網格文件自動生成全堆芯大規(guī)模網格，從而解決了傳統(tǒng)方法受計算機內存資源、網格文件大小的限制而導致的大規(guī)模網格劃分難以實現的問題。包括網格劃分功能在內，本文提供了有限元分析前處理的完整功能，使得生成的全堆芯網格文件能直接作為有限元求解軟件的輸入文件。

基于AMD GPU體系結構，研究了CEFR堆芯數值模型構造的稀疏剛度矩陣線性方程組Ku=r的并行計算。采用適合于對稱正定剛度矩陣K的Krylov子空間迭代法的CG法求解并使用Jacobi預處理加速。結合AMD GPU體系結構，采用數據預取和數據對齊的并行優(yōu)化方法，在AMD GPU體系結構計算系統(tǒng)上，通過對CEFR單組件和全堆芯組件的模擬計算，大規(guī)模稀疏剛度線性方程組Ku=r的計算效率得到了顯著的提高。另外，在數據并行性方面，根據AMD GPU體系結構的特點優(yōu)化了內存訪問和數據布局。