數(shù)值堆結(jié)構(gòu)力學高精細模擬的并行實現(xiàn)和優(yōu)化技術(shù)

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

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

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

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

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

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

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

1.1 數(shù)值堆幾何建模

采用反應(yīng)堆全堆芯結(jié)構(gòu)的標準尺寸，對各類型組件分別進行一次幾何建模，并將幾何模型的底面中心點坐標設(shè)置為(0，0，0)。以中國實驗快堆(CEFR)為例，采用常用的CAD模型文件格式(STEP、IGES、Brep、GEO等[13])，分別對其中的燃料組件、乏燃料組件、鋼屏蔽組件、硼屏蔽組件、安全棒組件、調(diào)節(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。對每根組件進行編號后，計算求得每根組件的中心點坐標，并建立組件編號與組件中心點坐標的映射關(guān)系。采用上述方法確定CEFR全堆芯712根組件中心點坐標后，根據(jù)組件類型對每根組件進行編號。通過輸入所需組件的編號集合，即可導出所需的反應(yīng)堆各類型組件的幾何模型文件及幾何位置數(shù)據(jù)。圖3為CVR幾何建模實現(xiàn)流程。

圖3 CVR幾何建模實現(xiàn)流程Fig.3 CVR geometric modeling flow chart

1.2 有限元分析預(yù)處理

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

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

圖4 有限元分析預(yù)處理實現(xiàn)流程Fig.4 Finite element analysis preprocessing flow chart

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

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

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

表2 常用網(wǎng)格生成算法的時空復(fù)雜度Table 2 Time and space complexity of common mesh generation algorithm

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

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

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

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

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

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

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

2.1 異構(gòu)平臺上的并行計算

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

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

圖8 CPU與GPU結(jié)構(gòu)的差別Fig.8 Difference between CPU and GPU structures

2.2 并行方案

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

1) 密集內(nèi)存訪問

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

2) 同步要求

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

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

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

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

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

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

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

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

2) 數(shù)據(jù)預(yù)取

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

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

圖10 流水線并行中的數(shù)據(jù)預(yù)取Fig.10 Data prefetching in pipeline parallel

3) 數(shù)據(jù)對齊

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

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

3 性能與驗證

3.1 測試算例與測試目標

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

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

3.2 AMD GPU的加速性能

將CEFR單燃料組件劃分為1 947 584個有限元網(wǎng)格，通過有限元力學數(shù)值模型將其轉(zhuǎn)化為稀疏剛度矩陣線性方程組的形式Ku=r，然后使用Krylov子空間迭代法中的CG方法和Jacobi預(yù)處理，結(jié)合在AMD GPU中的并行方法對其進行計算。CPU計算核心從1個CPU Core擴展到32個，并使用1～4個AMD GPU與之結(jié)合。實驗結(jié)果展示了在其并行算法下使用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全堆芯百億網(wǎng)格劃分

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

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

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

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

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

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

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

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

圖15 六角形套管單組件彎曲變形量(a)和軸向伸長量(b)模擬結(jié)果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

結(jié)合異構(gòu)計算系統(tǒng)的AMD GPU，對整個CEFR全堆芯組件劃分為0.54～102億網(wǎng)格的計算模型進行靜力分析模擬，算例列于表8，靜力模擬結(jié)果的形變位移分布如圖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億網(wǎng)格量的CEFR全堆芯組件計算模型進行靜力模擬，計算節(jié)點從64個節(jié)點擴展到512個節(jié)點，每個節(jié)點使用258個計算核心，最高使用到132 096個。CPU+AMD GPU加速比與效率對比和并行弱可擴展性如圖17所示。在異構(gòu)AMD GPU計算系統(tǒng)上與CPU的版本相比，提高了30%以上的計算性能。另外，以9 504計算核心(CPU+AMD GPU)對應(yīng)1.584億網(wǎng)格量為基準擴展至615 168計算核心對應(yīng)102.528億網(wǎng)格量，測試出并行弱可擴展性達到53%。

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

4 結(jié)論

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

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