邢龍?jiān)?，王昭順，蔣章程，高付海，高 岳，馮仰德，王 玨

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 異構(gòu)平臺(tái)上的并行計(jì)算

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

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

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

2.2 并行方案

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

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

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

2) 同步要求

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

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

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

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

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

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

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

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

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

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

提素質(zhì)、強(qiáng)能力，創(chuàng)新人才工作新機(jī)制。一直以來西南銅業(yè)公司黨委匡正選人用人風(fēng)氣，完善選人用人機(jī)制，堅(jiān)持嚴(yán)管和厚愛結(jié)合、激勵(lì)和約束并重，優(yōu)化干部考核評價(jià)機(jī)制。抓好各級干部的領(lǐng)導(dǎo)、監(jiān)管、屬地和執(zhí)行責(zé)任，提升黨政共融，黨務(wù)干部要抓安全、懂生產(chǎn)、謀經(jīng)營，行政干部要講政治、會(huì)黨建、抓文化。完善激勵(lì)機(jī)制和容錯(cuò)糾錯(cuò)機(jī)制。開展好選人用人制度修訂、培訓(xùn)及考核、中層管理人員畫像、后備干部選拔等工作。堅(jiān)持黨管人才原則，推進(jìn)員工人才成長體系建設(shè)，明晰各條人才通道的晉升、考核、退出、流動(dòng)機(jī)制，為培養(yǎng)公司“素質(zhì)優(yōu)良、結(jié)構(gòu)合理、隊(duì)伍穩(wěn)定”的員工隊(duì)伍打下扎實(shí)基礎(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ì)算，從而導(dǎo)致性能損失。因此，需根據(jù)向量寄存器上緩存線的大小來設(shè)置動(dòng)態(tài)分配的數(shù)據(jù)對齊字節(jié)數(shù)，圖11為512位向量寄存器上的數(shù)據(jù)未對齊。

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

3 性能與驗(yàn)證

3.1 測試算例與測試目標(biāo)

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

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

3.2 AMD GPU的加速性能

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

圖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)格文件；再將各種類型組件之間的相對位置信息寫入配置文件；最后自動(dòng)生成全堆芯大規(guī)模網(wǎng)格文件。CEFR全堆芯百億網(wǎng)格生成實(shí)現(xiàn)流程如圖13所示。

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

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

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

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

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

表5 IAEA標(biāo)準(zhǔn)例題靜力學(xué)驗(yàn)證算例Table 5 Static verification example of IAEA standard example

通過對程序計(jì)算結(jié)果和理論值的對比，分別給出了六角形套管單組件在給定線性溫度梯度下彎曲變形量(表6)和軸向伸長量(表7)在高度為4 000、3 000和2 500 mm處的偏差情況。組件彎曲變形量和軸向伸長量的數(shù)值解與理論值偏差較小，驗(yàn)證了計(jì)算結(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)計(jì)算系統(tǒng)的AMD GPU，對整個(gè)CEFR全堆芯組件劃分為0.54～102億網(wǎng)格的計(jì)算模型進(jìn)行靜力分析模擬，算例列于表8，靜力模擬結(jié)果的形變位移分布如圖16所示。

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

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

4 結(jié)論

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

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