張子彬，李志輝*，白智勇，彭傲平

（1.中國(guó)空氣動(dòng)力研究與發(fā)展中心超高速空氣動(dòng)力研究所，綿陽(yáng)621000；2.北京航空航天大學(xué)北京前沿創(chuàng)新中心國(guó)家計(jì)算流力學(xué)實(shí)驗(yàn)室，北京100191）

1 引言

航天器軌道計(jì)算與飛行航跡數(shù)值預(yù)報(bào)的準(zhǔn)確可靠直接依賴(lài)于軌道動(dòng)力學(xué)方程空氣動(dòng)力項(xiàng)求解的快速高效［1］，研究建立適于大型航天器在軌與服役期滿軌道衰降過(guò)程全飛行流域高超聲速空氣動(dòng)力特性一體化模擬平臺(tái)［2］與融合修正沿軌道空氣動(dòng)力一體化快速計(jì)算方法是關(guān)鍵基礎(chǔ)［3］。當(dāng)航天器運(yùn)動(dòng)在200 km以上時(shí)，飛行繞流已屬于高稀薄自由分子流，此時(shí)氣體分子平均自由程將超過(guò)200 m，Knudsen數(shù)將大于10。對(duì)于這種情況，傳統(tǒng)觀點(diǎn)認(rèn)為航天器飛行阻力、升力等空氣動(dòng)力將不再隨飛行高度發(fā)生變化，因此應(yīng)該采用固定的氣動(dòng)力系數(shù)參與飛行動(dòng)力學(xué)軌道控制的［3］。但是，中國(guó)天宮飛行器在采用固定阻力系數(shù)2.2開(kāi)展低軌道飛行試驗(yàn)時(shí)，發(fā)現(xiàn)有時(shí)飛行偏差越來(lái)越大，須使用RCS姿控系統(tǒng)強(qiáng)制控回原標(biāo)稱(chēng)軌道，造成燃料消耗，甚至難以準(zhǔn)確有效控制［4］。因此，提出這樣的問(wèn)題：天宮飛行器低軌道飛行所受大氣阻尼系數(shù)是否會(huì)變化？為了解決這一問(wèn)題，李志輝等［4］使用修正Boettcher/Legge非對(duì)稱(chēng)橋函數(shù)，發(fā)展了適于高超聲速稀薄過(guò)渡流區(qū)氣動(dòng)力計(jì)算的當(dāng)?shù)鼗瘶蚬椒?，建立了?zhuān)門(mén)針對(duì)復(fù)雜巨大的天宮飛行器跨越大氣層不同高度、不同馬赫數(shù)、不同迎角與側(cè)滑角空氣動(dòng)力特性快速計(jì)算方法與程序軟件，以減小數(shù)值計(jì)算工作量，并融合彈（軌）道飛行動(dòng)力學(xué)計(jì)算推進(jìn)，得到了具有實(shí)用價(jià)值的結(jié)果。實(shí)際應(yīng)用發(fā)現(xiàn)，在空氣動(dòng)力融合軌道飛行航跡數(shù)值預(yù)報(bào)過(guò)程中，該類(lèi)算法程序的運(yùn)算速度與能力效率仍有待進(jìn)一步提升［3］。為此，本文研究基于CUDA架構(gòu)的GPU異構(gòu)高性能并行計(jì)算技術(shù)，開(kāi)展大型航天器軌道衰降過(guò)程空氣動(dòng)力特性一體化建模并行優(yōu)化設(shè)計(jì)與應(yīng)用嘗試。

隨著科學(xué)技術(shù)的迅猛發(fā)展，高性能計(jì)算已成為科學(xué)技術(shù)發(fā)展和重大工程設(shè)計(jì)中具有戰(zhàn)略意義的研究手段［5］。早期計(jì)算機(jī)性能的提升主要依賴(lài)CPU（中央處理器）主頻的提升，但由于CPU的功耗與頻率的三次方近似成正比，無(wú)限制地提升頻率已不可能。到2003年，CPU頻率的提升接近停止［6］。為了應(yīng)對(duì)這一情況，人們往往采用基于CPU的并行計(jì)算方法來(lái)解決這一問(wèn)題（例如消息傳遞接口，簡(jiǎn)稱(chēng)MPI），但由于單個(gè)的CPU最多只能集成十幾、幾十核，導(dǎo)致基于CPU的并行計(jì)算能力受到很大限制。

2007 年，隨著基于英偉達(dá)CUDA架構(gòu)圖形處理單元（GPU）硬件的出現(xiàn)，集成核心的能力比CPU要強(qiáng)大得多（例如，英偉達(dá)Tesla架構(gòu)的K20X芯片包含2688個(gè)核心），這使得GPU的運(yùn)算能力比CPU更勝一籌［7］。除此之外，CUDA架構(gòu)不但包含英偉達(dá)GPU硬件，而且包含一個(gè)軟件編程環(huán)境，在這樣的編程模型下，CPU負(fù)責(zé)整個(gè)任務(wù)的流程控制和任務(wù)的串行計(jì)算，而GPU采用并行算法完成任務(wù)中計(jì)算密集部分的計(jì)算，人性化的編程環(huán)境有力助推了GPU的廣泛應(yīng)用。

時(shí)至今日，GPU和并行技術(shù)的使用已成為當(dāng)今數(shù)值計(jì)算的潮流和趨勢(shì)，并被廣泛應(yīng)用于包括計(jì)算流體力學(xué)在內(nèi)的許多領(lǐng)域。同樣地，它們也為求解稀薄氣體流動(dòng)問(wèn)題提供了一條行之有效的途徑。特別是使用NVIDIA的Tesla GPU在計(jì)算流體動(dòng)力學(xué)中的應(yīng)用加速明顯，針對(duì)不可壓縮的納維—斯托克斯（Navier-Stokes）方程數(shù)值求解模型的并行實(shí)現(xiàn)，在1024×32×1024的規(guī)模下使用4個(gè)Tesla C870 s與使用AMD Opteron 2.4 GHz的性能相比約提高100倍；針對(duì)格子Boltzmann方法在格子大小128×128規(guī)模下，使用Tesla 8-series GPU相比于使用Intel Xeon的性能約提高了130倍［8］?？梢?jiàn)，GPU和CPU融合的并行架構(gòu)是未來(lái)加速通用計(jì)算必然趨勢(shì)。

本文基于英偉達(dá)CUDA架構(gòu)的GPU，以微型計(jì)算機(jī)為平臺(tái)，根據(jù)近地軌道類(lèi)天宮飛行器空氣動(dòng)力特性當(dāng)?shù)鼗瘶蚝瘮?shù)關(guān)聯(lián)快速算法的運(yùn)算流程及對(duì)應(yīng)程序軟件的編寫(xiě)特點(diǎn)，提出相應(yīng)并行優(yōu)化方案，并付諸實(shí)施，以提升程序計(jì)算效率，減少程序運(yùn)行時(shí)間。并以類(lèi)天宮一號(hào)目標(biāo)飛行器無(wú)控飛行軌道衰降過(guò)程空氣動(dòng)力特性一體化建模作為算例，給出其在低軌道飛行過(guò)程相應(yīng)高度飛行繞流狀態(tài)的氣動(dòng)特性，檢驗(yàn)檢測(cè)本文并行計(jì)算模型方案的并行加速有效性與高效率，以此為基礎(chǔ)，逐步開(kāi)拓GPU并行計(jì)算在融合彈（軌）道數(shù)值預(yù)報(bào)關(guān)于空氣動(dòng)力計(jì)算的快速高效與準(zhǔn)確性。

2 空氣動(dòng)力特性當(dāng)?shù)鼗?jì)算方法

融合航天器在軌飛行與離軌再入飛行航跡彈（軌）道數(shù)值預(yù)報(bào)的跨流域高超聲速空氣動(dòng)力計(jì)算所用的當(dāng)?shù)鼗?jì)算方法，是基于自由分子流與連續(xù)流當(dāng)?shù)貥蚝瘮?shù)理論關(guān)聯(lián)的橋公式法［9-10］。自由分子流和牛頓連續(xù)流作為航天器再入飛行過(guò)程2個(gè)對(duì)應(yīng)極限流態(tài)，其中稀薄過(guò)渡流區(qū)域的氣動(dòng)力系數(shù)可以用當(dāng)?shù)貧鈩?dòng)力的解析來(lái)表述，當(dāng)?shù)孛嬖獨(dú)鈩?dòng)力系數(shù)采用基于自由分子流與連續(xù)流的當(dāng)?shù)鼗瘶蚝瘮?shù)光滑搭接［11］。這是一種加權(quán)函數(shù)，依賴(lài)于獨(dú)立的關(guān)聯(lián)或相似參數(shù)（如克努森數(shù)Kn），通常由實(shí)驗(yàn)結(jié)果與數(shù)值計(jì)算分析確定。橋公式法根據(jù)自由分子流和牛頓連續(xù)流不同的物理定律，采用兩個(gè)不同的數(shù)學(xué)函數(shù)FCont和FFM來(lái)描述，而在中間區(qū)域給出一個(gè)新的函數(shù)Fb，即橋函數(shù)。由此，對(duì)于給定的入射角θ，當(dāng)?shù)孛嬖辖?jīng)歸一化壓力和摩擦力系數(shù)表達(dá)為式（1）［1］：

式中，F(xiàn)b，p和Fb，τ分別是壓力和摩擦力橋函數(shù)，依賴(lài)于獨(dú)立參數(shù)Kn、TW/T0和θ。

對(duì)過(guò)渡流區(qū)域，當(dāng)?shù)貥蚬剑?2］可以提供依賴(lài)于獨(dú)立參數(shù)的歸一化壓力系數(shù)和摩擦力系數(shù)。當(dāng)?shù)貥蚝瘮?shù)的表達(dá)式有許多種，本文研究發(fā)展修正的Boettcher/Legge非對(duì)稱(chēng)橋函數(shù)表達(dá)式，其中，非對(duì)稱(chēng)的壓力橋函數(shù)可表示為式（2）：

式中，Knn，p、ΔKnp1、ΔKnp2根據(jù)飛行器外形及繞流特點(diǎn)為可調(diào)關(guān)聯(lián)參數(shù)。

非對(duì)稱(chēng)摩擦力橋函數(shù)可表示為式（3）：

發(fā)展可分段描述的非對(duì)稱(chēng)壓力與摩擦力系數(shù)關(guān)聯(lián)橋函數(shù)快速計(jì)算模型，需要根據(jù)飛行器外形及繞流特點(diǎn)調(diào)試確定Knn，p、ΔKnp1、ΔKnp2及Knn，τ、ΔKnτ1、ΔKnτ26個(gè)關(guān)聯(lián)參數(shù)，由于天宮飛行器的大尺度復(fù)雜結(jié)構(gòu)與在軌飛行跨流域多尺度非平衡繞流特點(diǎn)，需要依靠高性能大規(guī)模并行計(jì)算機(jī)開(kāi)展求解Boltzmann模型方程氣體動(dòng)理論統(tǒng)一算法等跨流域空氣動(dòng)力學(xué)數(shù)值計(jì)算典型飛行狀態(tài)模擬結(jié)果［2］，最優(yōu)擬合計(jì)算修正確定，使當(dāng)?shù)鼗こ逃?jì)算結(jié)果最大程度地滿足實(shí)際計(jì)算精度［9］。

關(guān)于飛行器物形處理，選取一種較為通用和近似的面元非結(jié)構(gòu)網(wǎng)格物形表征處理方法［4］，將飛行器物面劃分為若干個(gè)面元，利用四邊形或三角形來(lái)逼近復(fù)雜外形，圖1繪出類(lèi)似表征天宮飛行器兩艙體的大尺度圓柱體計(jì)算模型表面三角非結(jié)構(gòu)網(wǎng)格布局，可計(jì)算出每個(gè)面元的中心點(diǎn)坐標(biāo)、中心點(diǎn)處的法向及面元面積［13］。由于類(lèi)天宮飛行器屬于復(fù)雜非規(guī)則巨大板艙組合體，為了解決類(lèi)天宮飛行器及附件復(fù)雜物形表征的困難，引入分區(qū)網(wǎng)格與非結(jié)構(gòu)網(wǎng)格生成技術(shù)，進(jìn)行復(fù)雜物形表征處理，以便計(jì)算面元所受氣動(dòng)力，將所有面元上的氣動(dòng)力加起來(lái)就可以得到整個(gè)飛行器的氣動(dòng)力。

圖1 類(lèi)天宮飛行器兩艙體圓柱體外形表面非結(jié)構(gòu)三角形面元網(wǎng)格Fig.1 Unstructured triangular surface elementmesh on cylinder surface for two-capsule body of Tiangong type spacecraft

3 并行優(yōu)化方法

傳統(tǒng)的空氣動(dòng)力特性當(dāng)?shù)鼗?jì)算程序沒(méi)有采用并行方案，因此只能使用其中1個(gè)核心和1個(gè)線程，且每次只能執(zhí)行1次循環(huán)，造成了計(jì)算資源的閑置與浪費(fèi)。為了滿足氣動(dòng)融合軌道飛行航跡數(shù)值預(yù)報(bào)對(duì)空氣動(dòng)力特性實(shí)時(shí)計(jì)算快速高效要求，需要發(fā)展多核處理單元的并行優(yōu)化方法［14］。

為此，對(duì)類(lèi)天宮飛行器空氣動(dòng)力特性當(dāng)?shù)鼗惴ǖ倪\(yùn)算流程及對(duì)應(yīng)程序代碼進(jìn)行整體分析發(fā)現(xiàn)：

1）程序熱點(diǎn)代碼較為集中。原有程序主要執(zhí)行的代碼集中在算法的求解部分，具體表現(xiàn)為多個(gè)連續(xù)的3層嵌套循環(huán)。其中最內(nèi)層的循環(huán)次數(shù)最多，達(dá)到十萬(wàn)級(jí)別，適合數(shù)百個(gè)以上核心同時(shí)計(jì)算，有利于發(fā)揮GPU設(shè)備的性能優(yōu)勢(shì)。

2）數(shù)據(jù)依賴(lài)性較弱。數(shù)據(jù)依賴(lài)是指下一條對(duì)數(shù)據(jù)操作的指令必須等待上一條操作該數(shù)據(jù)的指令完成。原程序各個(gè)循環(huán)間保持相互獨(dú)立，不存在數(shù)據(jù)依賴(lài)關(guān)系，有利于發(fā)揮GPU設(shè)備的性能優(yōu)勢(shì)。

3）循環(huán)內(nèi)外數(shù)據(jù)傳輸需求較少。長(zhǎng)期以來(lái)，GPU與CPU之間數(shù)據(jù)傳輸較慢，一直是GPU在應(yīng)用過(guò)程中的詬病。具體來(lái)說(shuō)，由于GPU與CPU之間的數(shù)據(jù)傳輸速度遠(yuǎn)低于GPU、CPU芯片與各自寄存器的傳輸速度，甚至不及CPU與內(nèi)存、GPU與顯存之間的傳輸速度，因此造成引入GPU設(shè)備被迫增加的數(shù)據(jù)傳輸時(shí)間，會(huì)小于引入GPU設(shè)備所節(jié)省時(shí)間的情況，即損耗大于收益。而本文所涉及的程序由于循環(huán)內(nèi)外數(shù)據(jù)傳輸需求較少，可以避免這一情況，因此有利于發(fā)揮GPU設(shè)備的性能優(yōu)勢(shì)。

由上述特點(diǎn)可見(jiàn)，原有程序擁有改寫(xiě)為并行程序的潛力，并可以較好地發(fā)揮GPU設(shè)備的性能優(yōu)勢(shì)。因此，針對(duì)原有程序的編寫(xiě)特點(diǎn)，我們引入CUDA架構(gòu)的GPU設(shè)備，根據(jù)調(diào)試過(guò)程中出現(xiàn)的問(wèn)題，分別采用系統(tǒng)級(jí)別、算法級(jí)別以及語(yǔ)句級(jí)別3個(gè)層次的優(yōu)化手段對(duì)其進(jìn)行優(yōu)化。

3.1 系統(tǒng)級(jí)別優(yōu)化

系統(tǒng)級(jí)別手段優(yōu)化主要考慮程序性能的控制因素。對(duì)本文要解決的問(wèn)題，首先分析原程序是否存在以下3個(gè)限制因素：①處理器利用率過(guò)高還是過(guò)低；②存儲(chǔ)器帶寬利用是否不足；③阻塞處理器運(yùn)算的其他因素。

CPU、GPU及存儲(chǔ)器連接關(guān)系如圖2所示，測(cè)試發(fā)現(xiàn)，GPU與CPU以及GPU與顯存的數(shù)據(jù)傳輸速度比較緩慢，因此對(duì)這兩部分著重進(jìn)行了優(yōu)化設(shè)計(jì)如下：

1）對(duì)齊訪問(wèn)。由于在GPU顯存數(shù)據(jù)傳輸時(shí)，可以通過(guò)32、64或128 Byte的作業(yè)來(lái)訪問(wèn)，這些作業(yè)都對(duì)齊到它們的尺寸。因此未對(duì)齊的訪問(wèn)會(huì)增加訪問(wèn)次數(shù)，浪費(fèi)傳輸資源。就本程序而言，由于算法使用GPU運(yùn)算時(shí)有很多零散的參數(shù)需要傳入，因此將其整合為一個(gè)數(shù)組，進(jìn)行一次性傳輸，方便訪問(wèn)數(shù)據(jù)的對(duì)齊。圖3為本文設(shè)計(jì)使用的GPU內(nèi)存對(duì)齊訪問(wèn)格式。

圖2 CPU、GPU及存儲(chǔ)器連接關(guān)系示意圖Fig.2 Connection diagram of CPU，GPU and correspondingmemory

圖3 基于CUDA架構(gòu)的GPU內(nèi)存對(duì)齊訪問(wèn)Fig.3 Aligned access of GPU memory based on CUDA

2）使用數(shù)據(jù)傳輸函數(shù)。由于GPU設(shè)備傳統(tǒng)的數(shù)據(jù)傳輸方式效率較低，因此高版本CUDA架構(gòu)的GPU設(shè)備提供了豐富的數(shù)據(jù)傳輸函數(shù)，可以完成CPU到GPU或者GPU到另一臺(tái)GPU設(shè)備的數(shù)據(jù)傳輸，通過(guò)測(cè)試，可選用cudaMemcpy函數(shù)作為傳輸接口，具體格式為：cudaMemcpy（dst，src，count，kdir）。其中，dst、src分別為目標(biāo)變量和源變量，count為拷貝的數(shù)據(jù)長(zhǎng)度，kdir為可省略參數(shù)。

3.2 算法級(jí)別優(yōu)化

算法級(jí)別的優(yōu)化通常涉及程序部分，而這部分包含一個(gè)或多個(gè)函數(shù)，甚至只有一段語(yǔ)句。算法級(jí)別優(yōu)化主要涉及算法實(shí)現(xiàn)時(shí)要考慮的問(wèn)題，通常算法要考慮數(shù)據(jù)的組織、算法實(shí)現(xiàn)的策略。

一般串行算法求解的程序部分表現(xiàn)為多個(gè)連續(xù)的三層嵌套循環(huán)，且循環(huán)主要集中在最內(nèi)層，達(dá)到十萬(wàn)級(jí)別。以NVIDIA公司的GTX750Ti圖形處理單元為例，它擁有640個(gè)計(jì)算核心，理論上可支持131 072個(gè)線程同時(shí)計(jì)算［14］。從線程數(shù)來(lái)說(shuō)，每個(gè)循環(huán)都可以交給一個(gè)指定的線程進(jìn)行計(jì)算，具有較好的并行性；從核心數(shù)來(lái)講，十萬(wàn)余個(gè)循環(huán)可以拆分給多個(gè)GPU芯片共同完成，具有較好的可擴(kuò)展性。

因此，本文將算法求解部分的最內(nèi)層循環(huán)經(jīng)過(guò)展開(kāi)與合并處理后，整理為一個(gè)循環(huán)，再整體移植入核函數(shù)，從而利用GPU強(qiáng)大的并行計(jì)算能力提升運(yùn)算效率。具體代碼如下所述。

原有串行程序代碼為：

do i=1，ntr！ntr=121822

……！數(shù)組的循環(huán)運(yùn)算A

end do

……！數(shù)組的共同運(yùn)算B

do i=1，ntr！

……！數(shù)組的循環(huán)運(yùn)算C

end do

……！數(shù)組的共同運(yùn)算D

經(jīng)過(guò)優(yōu)化的程序代碼為：

attributes subroutine CUDA（input，output）

！GPU設(shè)備核函數(shù)開(kāi)頭

do i=1，ntr！ntr=121822

……！運(yùn)算A+B+C+D

end do

end subroutine

3.3 語(yǔ)句級(jí)別優(yōu)化

語(yǔ)句級(jí)別的優(yōu)化包括對(duì)函數(shù)、循環(huán)、指令等語(yǔ)句細(xì)節(jié)的優(yōu)化。

CUDA架構(gòu)的GPU設(shè)備顯式地將存儲(chǔ)器分成4種不同的類(lèi)型：全局存儲(chǔ)器（globalmemory）、常量存儲(chǔ)器（constantmemory）、共享存儲(chǔ)器（share memory）或稱(chēng)局部存儲(chǔ)器（localmemory）、私有存儲(chǔ)器（private memory）。雖然其中除全局存儲(chǔ)器有較大容量之外（如tesla K20全局內(nèi)存為5 GB），其余3種存儲(chǔ)器的存儲(chǔ)容量只有64KB左右［14］，但是這4種存儲(chǔ)器的讀取速度卻有很大差異：全局存儲(chǔ)器是最慢的，其延遲為幾百個(gè)時(shí)鐘；常量存儲(chǔ)器次之，延遲為幾十個(gè)周期；共享存儲(chǔ)器的延遲在幾個(gè)至幾十個(gè)時(shí)鐘周期之內(nèi)；而私有存儲(chǔ)器通常就是硬件的寄存器，它是GPU的存儲(chǔ)器中讀取最快的，幾乎不用耗費(fèi)時(shí)間［14］。而相比較全局變量耗費(fèi)的幾百個(gè)周期，1個(gè)時(shí)鐘周期可以做4次整形加法、2次浮點(diǎn)乘加，3個(gè)周期可以做一次乘法，十幾個(gè)周期可以做一次乘法［15］。因此如果能夠合理利用常量存儲(chǔ)器、共享存儲(chǔ)器和私有存儲(chǔ)器，盡量減少讀取全局存儲(chǔ)器所產(chǎn)生的時(shí)間耗費(fèi)，那么所產(chǎn)生的收益將是相當(dāng)可觀的。

以下列出了私有存儲(chǔ)器、共享存儲(chǔ)器和常量存儲(chǔ)器的聲明使用方式：

integer：：t，tr！ 使用私有存儲(chǔ)器

real，shared：：s（64） ！ 使用共享存儲(chǔ)器

real，parameter：：PI=3.14159！ 使用常量存儲(chǔ)器

4 算法驗(yàn)證與結(jié)果分析

將所建立類(lèi)天宮飛行器空氣動(dòng)力特性當(dāng)?shù)鼗兴惴皩?duì)應(yīng)程序按提出的并行計(jì)算方案進(jìn)行CPU+GPU移植優(yōu)化后，對(duì)類(lèi)天宮一號(hào)目標(biāo)飛行器外形進(jìn)行運(yùn)算，分析優(yōu)化前后程序運(yùn)行性能。

4.1 計(jì)算模型與條件

實(shí)施對(duì)于天宮一號(hào)目標(biāo)飛行器低地球軌道飛行過(guò)程氣動(dòng)力系數(shù)一體化計(jì)算，采用基于三角非結(jié)構(gòu)網(wǎng)格表征物面，繪出計(jì)算模型表面三角非結(jié)構(gòu)網(wǎng)格布局如圖4所示。

圖4 天宮飛行器計(jì)算模型表面非結(jié)構(gòu)網(wǎng)格布局Fig.4 Unstructured triangular surface elementmesh on the com puting model of Tiangong spacecraft

關(guān)于極不規(guī)則大型復(fù)雜結(jié)構(gòu)航天器氣動(dòng)特性計(jì)算，參考面積選取可分為2部分考慮：第一部分是本體迎風(fēng)面積，即飛行器本體的底面積；第二部分是帆板的迎風(fēng)面積，由于太陽(yáng)電池帆板需要圍繞中心軸轉(zhuǎn)動(dòng)，于是帆板迎風(fēng)面積采用等效面積計(jì)算；總的參考面積為上述兩部分面積之和［4］。

本算例采用多核CPU+GPU架構(gòu)計(jì)算機(jī)進(jìn)行運(yùn)算，CPU及GPU性能參數(shù)見(jiàn)表1所列。對(duì)比試驗(yàn)中，串行程序僅由CPU獨(dú)立完成，共調(diào)用1個(gè)核心和一個(gè)線程。并行程序的算法求解運(yùn)算由GPU完成，調(diào)用640個(gè)核心和131 072個(gè)線程，其余步驟如數(shù)據(jù)載入、結(jié)果收集及試驗(yàn)監(jiān)測(cè)等過(guò)程仍由CPU完成。

表1 試驗(yàn)設(shè)備性能列表Table 1 Performance list of test equipments

4.2 GPU并行算法與CPU串行計(jì)算結(jié)果驗(yàn)證

圖5與圖6分別繪出天宮飛行器角α為22°，阻力系數(shù)與升力系數(shù)隨側(cè)滑角β變化趨勢(shì)其中藍(lán)色曲線表示CPU串行程序計(jì)算結(jié)果，而紅色方格代表并行程序計(jì)算結(jié)果。由圖看出，兩者完全吻合，本文經(jīng)過(guò)對(duì)156次計(jì)算結(jié)果進(jìn)行統(tǒng)計(jì)，除部分個(gè)別數(shù)值末位略有不同且結(jié)果僅差0.000 01，誤差緣由可能是CPU和GPU芯片取舍原則不同所致，其余數(shù)據(jù)完全一致，這種GPU與CPU計(jì)算的高精度一致性，進(jìn)一步證實(shí)本文GPU并行算法實(shí)現(xiàn)的正確性與高效率。

圖5 GPU加速前后阻力系數(shù)隨側(cè)滑角變化計(jì)算驗(yàn)證Fig.5 Computing verification of drag variation w ith sideslip angle before and after GPU acceleration

圖6 GPU加速前后升力系數(shù)隨側(cè)滑角變化計(jì)算驗(yàn)證Fig.6 Com puting verification of lift variation w ith sideslip angle before and after GPU acceleration

表2列出原串行程序CPU計(jì)算與本文GPU并行算法程序及優(yōu)化后的GPU并行程序在求解過(guò)程中，單步求解各模塊程序段的平均消耗時(shí)間，其中GPU運(yùn)算是在GPU中發(fā)生的過(guò)程，數(shù)據(jù)歸約是在CPU中各節(jié)點(diǎn)進(jìn)行。由表2看出：①GPU并行程序雖然增加了數(shù)據(jù)輸入、輸出的過(guò)程，但總的運(yùn)行時(shí)間仍比CPU串行程序計(jì)算時(shí)間有了明顯減?。虎诮?jīng)過(guò)并行優(yōu)化后，數(shù)據(jù)在CPU與GPU之間傳輸?shù)臅r(shí)間也有了明顯減少。

表3給出了CPU串行、并行及引入GPU設(shè)備并優(yōu)化后的程序平均消耗時(shí)間，從表中可以看出，受制于CPU線程與計(jì)算核心數(shù)影響（CPU i7-4770 K只能使用8個(gè)線程與4個(gè)處理單元），CPU并行程序在單步求解時(shí)間方面的加速比只能達(dá)到7.0左右，而無(wú)法超過(guò)8，而GPU設(shè)備的相應(yīng)加速比可以達(dá)到13.0，受此影響，CPU并行程序與優(yōu)化后的GPU并行程序的整體加速比仍有一定差距，雖然近年來(lái)CPU的最大線程數(shù)和性能有所提升，但GPU設(shè)備的并行能力也在大幅度增長(zhǎng)，因此仍然沒(méi)有實(shí)質(zhì)性的區(qū)別，這樣的案例進(jìn)一步說(shuō)明了采用GPU并行的重要性和研究?jī)r(jià)值。

表2 算法單步求解各階段平均消耗時(shí)間Table 2 Average runtime per step in solving process/ms

表3 CPU串行、并行及GPU優(yōu)化后程序平均消耗時(shí)間Table 3 The average runtime time of CPU serial program，parallel program and GPU optim ized program

4.3 性能分析

根據(jù)常用并行算法性能衡量手段，采用加速比對(duì)計(jì)算試驗(yàn)進(jìn)行分析。加速比的定義為順序執(zhí)行時(shí)間與并行執(zhí)行時(shí)間比值。由表3可知，CPU串行程序消耗時(shí)間0.747 h，GPU并行程序消耗時(shí)間0.140 h，因此加速比為5.3，這說(shuō)明，將GPU并行加速引入類(lèi)天宮飛行器無(wú)控飛行軌降再入解體過(guò)程氣動(dòng)力系數(shù)一體化計(jì)算相當(dāng)快速高效。

比較分析還可看出：

1）在單次求解過(guò)程中，CPU循環(huán)運(yùn)算平均消耗75.5 ms，而優(yōu)化后的GPU運(yùn)算平均消耗5.8 ms，GPU并行計(jì)算加速比約為13?？梢?jiàn)，第3章引入GPU設(shè)備、發(fā)展的GPU并行優(yōu)化策略是行之有效加速手段；

2）通過(guò)優(yōu)化，單次求解過(guò)程中的數(shù)據(jù)傳輸耗時(shí)由原來(lái)的9.2 ms降為2.1 ms，減少了77%；數(shù)據(jù)歸約耗時(shí)由原來(lái)的2.0ms降為0.7 ms，減少了65%；GPU運(yùn)算過(guò)程也由原來(lái)的4 ms降為3 ms；合計(jì)時(shí)間由優(yōu)化前的GPU運(yùn)行時(shí)間15.2 ms降為5.8 ms，減少了61%。可見(jiàn)，本文介紹的GPU優(yōu)化手段成效顯著；

3）在試驗(yàn)中共使用線程512×256=131 072個(gè)，而GPU計(jì)算核心僅調(diào)用640個(gè)，這是因?yàn)槊總€(gè)核心上被分配了多個(gè)線程，這增加了線程等待時(shí)間，因此，當(dāng)調(diào)用核心數(shù)增加時(shí)，程序仍有進(jìn)一步擴(kuò)展的空間；

4）在單次求解過(guò)程中，GPU并行程序數(shù)據(jù)輸入輸出的平均耗時(shí)為2.6+6.6=9.2 ms，是單次求解過(guò)程平均耗時(shí)（15.2 ms）的近60%，即使在GPU并行優(yōu)化后，輸入輸出的平均耗時(shí)為0.9+1.2=2.1ms，是單次串行求解過(guò)程平均耗時(shí)（5.8 ms）的近36.2%?？梢?jiàn)，由于GPU與CPU之間的傳輸速度過(guò)慢，輸入輸出占用了GPU調(diào)用的相當(dāng)多時(shí)間，因此，如何發(fā)展高效快速的輸入輸出仍是限制GPU推廣應(yīng)用需要進(jìn)一步研究解決的瓶頸之一。

4.4 類(lèi)天宮飛行器軌道衰降過(guò)程空氣動(dòng)力特性計(jì)算分析

使用類(lèi)天宮飛行器軌道衰降過(guò)程空氣動(dòng)力特性一體化建模GPU并行算法程序，對(duì)天宮飛行器無(wú)控飛行軌道衰降過(guò)程340～120 km氣動(dòng)特性詳細(xì)計(jì)算，圖7與圖8分別繪出太陽(yáng)電池帆板面與天宮飛行器本體主軸夾角10°、迎角α=10°進(jìn)行340～120 km不同高度340 km、260 km、180 km、140 km、120 km的氣動(dòng)阻力、升力隨不同側(cè)滑角變化輪廓。可看出對(duì)于所計(jì)算的不同飛行繞流狀態(tài)，均在側(cè)滑角β=0°阻力系數(shù)最小、升力系數(shù)最大，說(shuō)明天宮飛行器存在嚴(yán)重的板艙非規(guī)則構(gòu)型，一旦無(wú)控飛行姿態(tài)出現(xiàn)側(cè)滑角飛行繞流現(xiàn)象，迎風(fēng)面積、阻力系數(shù)都會(huì)迅速增大，升力系數(shù)快速減小，致失穩(wěn)加速飛行繞流狀態(tài)，印證了天宮一號(hào)目標(biāo)飛行器無(wú)控飛行軌降過(guò)程數(shù)值預(yù)報(bào)發(fā)現(xiàn)的，其自2017年9月初出現(xiàn)失穩(wěn)后三軸旋轉(zhuǎn)加速直至最終再入大氣層解體的無(wú)控飛行繞流現(xiàn)象。

圖7 天宮飛行器軌道衰降過(guò)程氣動(dòng)阻力系數(shù)GPU優(yōu)化并行計(jì)算結(jié)果與CPU串行計(jì)算驗(yàn)證Fig.7 Verification of GPU parallel program and CPU serial computation for aerodynam ic drag variation of Tiangong spacecraft during orbital declining at 340～120 km

圖9（a）、（b）分別繪出太陽(yáng)電池帆板面與天宮飛行器本體主軸夾角10°、α=22°、β=0°時(shí)天宮一號(hào)目標(biāo)飛行器軌降進(jìn)入臨近再入段240～120 km不同高度240 km、220 km、200 km、180 km、160 km、140 km、120 km的氣動(dòng)阻力、升力系數(shù)隨飛行高度變化趨勢(shì)?？煽闯鰧?duì)于所計(jì)算的不同飛行繞流狀態(tài)，隨著高度降低，阻力系數(shù)、升力系數(shù)呈現(xiàn)快速非線性下降趨勢(shì)，印證了天宮一號(hào)飛行航跡在此高度軌降過(guò)程中，將于20 d內(nèi)到達(dá)120 km再入大氣層解體墜毀的無(wú)控飛行變化現(xiàn)象。

5 結(jié)論

1）通過(guò)分析近地軌道類(lèi)天宮飛行器空氣動(dòng)力特性當(dāng)?shù)鼗瘶蚝瘮?shù)關(guān)聯(lián)算法及其對(duì)應(yīng)程序具有多連續(xù)三層嵌套循環(huán)熱點(diǎn)代碼集中、循環(huán)間彼此獨(dú)立不存在數(shù)據(jù)依賴(lài)關(guān)系及循環(huán)內(nèi)外數(shù)據(jù)傳輸需求少等特點(diǎn)，引入了CUDA架構(gòu)的GPU并行加速設(shè)備，提出了基于多核處理單元的并行優(yōu)化設(shè)計(jì)方案。

圖8 天宮飛行器軌道衰降過(guò)程氣動(dòng)升力系數(shù)GPU優(yōu)化并行計(jì)算結(jié)果與CPU串行計(jì)算驗(yàn)證Fig.8 Verification of GPU parallel program and CPU serial computation for aerodynam ic lift variation of Tiangong spacecraft during orbital declining at 340～120 km

圖9 天宮飛行器軌道衰降過(guò)程氣動(dòng)系數(shù)隨飛行高度GPU并行優(yōu)化計(jì)算與CPU串行計(jì)算驗(yàn)證Fig.9 Verification of GPU parallel program and CPU serial com putation for aerodynam ic drag and lift variation of Tiangong spacecraft during orbital declining at 340～120 km.

2）針對(duì)氣動(dòng)融合軌道數(shù)值預(yù)報(bào)快速高效特點(diǎn)，發(fā)展了系統(tǒng)、算法、語(yǔ)句級(jí)別的3個(gè)層次并行優(yōu)化方法。建立了基于GPU內(nèi)存對(duì)齊訪問(wèn)格式與數(shù)據(jù)傳輸函數(shù)方案、循環(huán)展開(kāi)與合并、優(yōu)化存儲(chǔ)器使用等CPU+GPU并行移植優(yōu)化技術(shù)，快速提升了GPU數(shù)據(jù)傳輸速度和運(yùn)算效率。

3）經(jīng)本文GPU并行算法與CPU串行計(jì)算試驗(yàn)驗(yàn)證表明，類(lèi)天宮飛行器空氣動(dòng)力特性當(dāng)?shù)鼗こ趟惴ǔ绦蚪?jīng)GPU并行移植優(yōu)化后，運(yùn)行程序消耗時(shí)間由CPU串行計(jì)算0.747 h，減小至GPU并行計(jì)算0.140 h，并行加速比為5.3，證實(shí)首次將GPU并行加速引入類(lèi)天宮飛行器無(wú)控飛行軌降再入解體過(guò)程氣動(dòng)力系數(shù)一體化計(jì)算相當(dāng)快速高效。

4）GPU作為圖形處理工具，在架構(gòu)上為了集成更多的計(jì)算核心而犧牲了存儲(chǔ)性能和單核計(jì)算能力，GPU雖然在并行加速能力上遠(yuǎn)勝于CPU，但是在數(shù)據(jù)傳輸及單核處理能力方面較弱。計(jì)算程序只有在被改寫(xiě)成數(shù)據(jù)依賴(lài)性好，傳輸需求少、計(jì)算量大的運(yùn)算模式，才能發(fā)揮GPU設(shè)備優(yōu)勢(shì)。本文提出的并行優(yōu)化設(shè)計(jì)方案，可以為其他算法程序在大規(guī)模并行計(jì)算實(shí)現(xiàn)方面提供經(jīng)驗(yàn)和借鑒參考。

5）本文工作屬初步階段性，將GPU并行移植優(yōu)化算法應(yīng)用于其他跨流域空氣動(dòng)力學(xué)數(shù)值模擬方法，特別是GPU與CPU之間傳輸速度慢，輸入輸出占用GPU調(diào)用較多時(shí)間，如何發(fā)展高效快速的輸入輸出是限制GPU推廣應(yīng)用瓶頸，均有待進(jìn)一步研究解決。