類天宮飛行器軌道衰降過程空氣動力特性一體化建模并行優(yōu)化設(shè)計(jì)

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

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

發(fā)展可分段描述的非對稱壓力與摩擦力系數(shù)關(guān)聯(lián)橋函數(shù)快速計(jì)算模型，需要根據(jù)飛行器外形及繞流特點(diǎn)調(diào)試確定Knn，p、ΔKnp1、ΔKnp2及Knn，τ、ΔKnτ1、ΔKnτ26個關(guān)聯(lián)參數(shù)，由于天宮飛行器的大尺度復(fù)雜結(jié)構(gòu)與在軌飛行跨流域多尺度非平衡繞流特點(diǎn)，需要依靠高性能大規(guī)模并行計(jì)算機(jī)開展求解Boltzmann模型方程氣體動理論統(tǒ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］，將飛行器物面劃分為若干個面元，利用四邊形或三角形來逼近復(fù)雜外形，圖1繪出類似表征天宮飛行器兩艙體的大尺度圓柱體計(jì)算模型表面三角非結(jié)構(gòu)網(wǎng)格布局，可計(jì)算出每個面元的中心點(diǎn)坐標(biāo)、中心點(diǎn)處的法向及面元面積［13］。由于類天宮飛行器屬于復(fù)雜非規(guī)則巨大板艙組合體，為了解決類天宮飛行器及附件復(fù)雜物形表征的困難，引入分區(qū)網(wǎng)格與非結(jié)構(gòu)網(wǎng)格生成技術(shù)，進(jìn)行復(fù)雜物形表征處理，以便計(jì)算面元所受氣動力，將所有面元上的氣動力加起來就可以得到整個飛行器的氣動力。

圖1 類天宮飛行器兩艙體圓柱體外形表面非結(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)?shù)鼗?jì)算程序沒有采用并行方案，因此只能使用其中1個核心和1個線程，且每次只能執(zhí)行1次循環(huán)，造成了計(jì)算資源的閑置與浪費(fèi)。為了滿足氣動融合軌道飛行航跡數(shù)值預(yù)報對空氣動力特性實(shí)時計(jì)算快速高效要求，需要發(fā)展多核處理單元的并行優(yōu)化方法［14］。

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

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

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

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

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

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

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

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

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

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

圖3 基于CUDA架構(gòu)的GPU內(nèi)存對齊訪問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到另一臺GPU設(shè)備的數(shù)據(jù)傳輸，通過測試，可選用cudaMemcpy函數(shù)作為傳輸接口，具體格式為：cudaMemcpy（dst，src，count，kdir）。其中，dst、src分別為目標(biāo)變量和源變量，count為拷貝的數(shù)據(jù)長度，kdir為可省略參數(shù)。

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

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

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

因此，本文將算法求解部分的最內(nèi)層循環(huán)經(jīng)過展開與合并處理后，整理為一個循環(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)過優(yōu)化的程序代碼為：

attributes subroutine CUDA（input，output）

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

do i=1，ntr！ntr=121822

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

end do

end subroutine

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

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

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

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

integer：：t，tr！ 使用私有存儲器

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

real，parameter：：PI=3.14159！ 使用常量存儲器

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

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

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

實(shí)施對于天宮一號目標(biāo)飛行器低地球軌道飛行過程氣動力系數(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)航天器氣動特性計(jì)算，參考面積選取可分為2部分考慮：第一部分是本體迎風(fēng)面積，即飛行器本體的底面積；第二部分是帆板的迎風(fēng)面積，由于太陽電池帆板需要圍繞中心軸轉(zhuǎn)動，于是帆板迎風(fēng)面積采用等效面積計(jì)算；總的參考面積為上述兩部分面積之和［4］。

本算例采用多核CPU+GPU架構(gòu)計(jì)算機(jī)進(jìn)行運(yùn)算，CPU及GPU性能參數(shù)見表1所列。對比試驗(yàn)中，串行程序僅由CPU獨(dú)立完成，共調(diào)用1個核心和一個線程。并行程序的算法求解運(yùn)算由GPU完成，調(diào)用640個核心和131 072個線程，其余步驟如數(shù)據(jù)載入、結(jié)果收集及試驗(yàn)監(jiān)測等過程仍由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è)滑角β變化趨勢其中藍(lán)色曲線表示CPU串行程序計(jì)算結(jié)果，而紅色方格代表并行程序計(jì)算結(jié)果。由圖看出，兩者完全吻合，本文經(jīng)過對156次計(jì)算結(jié)果進(jìn)行統(tǒng)計(jì)，除部分個別數(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并行程序在求解過程中，單步求解各模塊程序段的平均消耗時間，其中GPU運(yùn)算是在GPU中發(fā)生的過程，數(shù)據(jù)歸約是在CPU中各節(jié)點(diǎn)進(jìn)行。由表2看出：①GPU并行程序雖然增加了數(shù)據(jù)輸入、輸出的過程，但總的運(yùn)行時間仍比CPU串行程序計(jì)算時間有了明顯減小；②經(jīng)過并行優(yōu)化后，數(shù)據(jù)在CPU與GPU之間傳輸?shù)臅r間也有了明顯減少。

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

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

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

4.3 性能分析

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

比較分析還可看出：

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

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

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

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

4.4 類天宮飛行器軌道衰降過程空氣動力特性計(jì)算分析

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

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

5 結(jié)論

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

圖8 天宮飛行器軌道衰降過程氣動升力系數(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 天宮飛行器軌道衰降過程氣動系數(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）針對氣動融合軌道數(shù)值預(yù)報快速高效特點(diǎn)，發(fā)展了系統(tǒng)、算法、語句級別的3個層次并行優(yōu)化方法。建立了基于GPU內(nèi)存對齊訪問格式與數(shù)據(jù)傳輸函數(shù)方案、循環(huán)展開與合并、優(yōu)化存儲器使用等CPU+GPU并行移植優(yōu)化技術(shù)，快速提升了GPU數(shù)據(jù)傳輸速度和運(yùn)算效率。

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

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

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