胡斌星, 李新國(guó), 常武權(quán)
(1.西北工業(yè)大學(xué) 航天學(xué)院,西安 710072; 2.中國(guó)運(yùn)載火箭技術(shù)研究院,北京 100076)
隨著對(duì)飛行器性能要求的不斷提高,飛行器自重比例逐步減小,彈性結(jié)構(gòu)與其氣動(dòng)、控制系統(tǒng)的耦合問(wèn)題愈發(fā)嚴(yán)重,箭體在飛行過(guò)程中振動(dòng)特性對(duì)其動(dòng)態(tài)性能影響較大,對(duì)其準(zhǔn)確表征逐步成為新一代飛行器設(shè)計(jì)的核心技術(shù)內(nèi)容之一。由于火箭和導(dǎo)彈的長(zhǎng)細(xì)比較大,以往常用方法是將箭體簡(jiǎn)化為柔性梁分析;即使在有助推器的情況下,仍將助推器也視為一維梁與芯級(jí)組成空間梁系,基于特定的小偏差增量模型分析彈性振動(dòng)對(duì)控制系統(tǒng)的影響[1];在新一代高超聲速飛行器的分析中也常以自由梁或懸臂梁建模分析與氣動(dòng)、推進(jìn)系統(tǒng)的耦合關(guān)系[2]。在設(shè)計(jì)驗(yàn)證階段,通常將彈體離散成若干個(gè)梁?jiǎn)卧驈椈?質(zhì)點(diǎn)(站點(diǎn))模型,以此為基礎(chǔ)建立飛行器模態(tài)信息并提前裝訂至仿真計(jì)算機(jī);但在實(shí)際工程應(yīng)用中,由于材料、飛行器幾何外形的不規(guī)則以及簡(jiǎn)化過(guò)程模型精度損失等問(wèn)題的影響,過(guò)少的站點(diǎn)數(shù)已無(wú)法準(zhǔn)確表達(dá)飛行器自身的彈性特性,而過(guò)多的單元數(shù)已無(wú)法滿足仿真過(guò)程的實(shí)時(shí)性要求。為快速的在仿真過(guò)程中獲取較精確的振動(dòng)特性以便為姿態(tài)控制設(shè)計(jì)、響應(yīng)分析提供參考,勢(shì)必采用并行計(jì)算的方法解決[3]。一些文獻(xiàn)提出利用OpenMP(Open Multiple Processing)或MPI(Multi Point Interface)技術(shù)提高計(jì)算速度,但OpenMP技術(shù)從理論上僅能提供當(dāng)前計(jì)算機(jī)CPU核心數(shù)的最大加速性能,其受線程創(chuàng)建及任務(wù)分發(fā)的影響較大;而MPI在仿真中由于通信的延遲其實(shí)時(shí)性并不良好,且兩者維護(hù)上存在著較大的不便[4-5]。較之單純的CPU集群運(yùn)算,GPU(Graphics Processing Unit)的帶寬更高、延遲更低、單位浮點(diǎn)運(yùn)算的功耗更低,故現(xiàn)有的超級(jí)計(jì)算機(jī)將GPU作為協(xié)處理器以提升浮點(diǎn)運(yùn)算能力。而自2007年NVIDIA公司正式發(fā)布統(tǒng)一計(jì)算設(shè)備架構(gòu)(CUDA, Compute Unifed Device Architecture)并行計(jì)算技術(shù)以來(lái),由于其采用標(biāo)準(zhǔn)C語(yǔ)言擴(kuò)展的API形式大大降低可開(kāi)發(fā)門檻,通過(guò)PCIE總線與CPU連接可以保證一定的實(shí)時(shí)性,故近年來(lái)在有限元的顯式動(dòng)力學(xué)[6]、隱式動(dòng)力學(xué)[7]、結(jié)構(gòu)力學(xué)[8]、分子動(dòng)力學(xué)、方面有著廣泛應(yīng)用,并已成為高性能計(jì)算領(lǐng)域的重要分支。Fleischmann等[9]就MPI和OpenMP 及CUDA結(jié)合,在矩陣不完全LU分解和子空間迭代法求解特征值問(wèn)題的加速效果做了非常詳細(xì)的對(duì)比,驗(yàn)證了GPU在數(shù)值計(jì)算方面能取得明顯的性能提升;Cheng等[10]對(duì)克里金插值的GPU加速算法做了詳盡闡述,說(shuō)明了插值計(jì)算GPU加速的可能性;李濤等[11]實(shí)現(xiàn)了基于線程池的GPU任務(wù)并行計(jì)算模式,就矩陣運(yùn)算等一般性的運(yùn)算做了驗(yàn)證分析,但并未針對(duì)GPU特性實(shí)現(xiàn)細(xì)粒度的性能優(yōu)化。在多GPU運(yùn)算方面,賴劍奇等實(shí)現(xiàn)了多GPU的并行可壓縮流求解器,驗(yàn)證了多GPU良好的可擴(kuò)展性;Gao等[12-13]也發(fā)表了大型稀疏線性系統(tǒng)求解的多GPU解算方案,可上述文獻(xiàn)皆未關(guān)注除解算器以外的粗粒度并行問(wèn)題,僅在有限元法解算中實(shí)現(xiàn)了加速,設(shè)計(jì)適宜實(shí)時(shí)仿真的多GPU架構(gòu)并未闡述。因此本文根據(jù)以上研究?jī)?nèi)容,從GPU硬件特點(diǎn)著手,在單臺(tái)多GPU的工作站上設(shè)計(jì)了細(xì)長(zhǎng)體彈性飛行器彈性模塊異步并發(fā)并行計(jì)算的架構(gòu)設(shè)計(jì),實(shí)現(xiàn)了基于算法的GPU的動(dòng)態(tài)性能優(yōu)化,首次提出了動(dòng)態(tài)分配線程塊的八叉樹(shù)構(gòu)建及索引方法用以實(shí)現(xiàn)多線程快速索引氣動(dòng)數(shù)據(jù),對(duì)不同計(jì)算規(guī)模的性能進(jìn)行系統(tǒng)分析。對(duì)并行程序設(shè)計(jì)及優(yōu)化,三維數(shù)據(jù)快速索引插值方面具有重要指導(dǎo)意義。
為便于說(shuō)明細(xì)長(zhǎng)體彈性飛行器振動(dòng)響應(yīng)算法的并行實(shí)現(xiàn),在此結(jié)合算法分析計(jì)算流程,并分析各步驟的計(jì)算量及比例,以便決定并行算法設(shè)計(jì)的重點(diǎn)。如前所述,將全彈體簡(jiǎn)化為非均質(zhì)梁模型,以EJ(x)表示彎曲剛度,m(x)表示線質(zhì)量,l表示參考長(zhǎng)度,q(x,t)表示作用在梁元上的橫向分布外力,忽略梁元的轉(zhuǎn)動(dòng)慣性,且在無(wú)阻尼的情況下,依據(jù)梁的彎曲理論和力矩平衡方程,得到梁的彎曲振動(dòng)微分方程
(1)
當(dāng)q(x,t)=0時(shí)得到該振動(dòng)微分方程的通解部分;通過(guò)自由梁兩端的邊界條件可得特解;把梁的受迫振動(dòng)展開(kāi)成固有振型的級(jí)數(shù)形式如下
y(x,t)=yc(t)+θ(t)(x-xc)+
(2)
鑒于單位長(zhǎng)度分布外力可表示為乘積q(x,t)=q(t)P(x),把式(2)代入式(1)對(duì)x積分變換后得
(3)
式(3)中等式右端廣義力為
式(3)中第一項(xiàng)表示質(zhì)心運(yùn)動(dòng)方程,第二項(xiàng)表示繞橫軸的旋轉(zhuǎn)運(yùn)動(dòng)方程。當(dāng)給定外力q(t)P(x)時(shí),可以由最后一組方程求得梁彎曲振動(dòng)的廣義坐標(biāo)。經(jīng)推導(dǎo),可將細(xì)長(zhǎng)體飛行器的法向和橫向彈性振動(dòng)方程數(shù)值解算步驟歸納為
(4)
(5)
Qfayi=Ya(xj)×fiy(xj)
(6)
Qfayi+Py+Feiy+Fpy+Fdy
(7)
(8)
(9)
圖1 飛行動(dòng)力學(xué)仿真中彈性模塊計(jì)算流程
表1 不同站點(diǎn)數(shù)CPU版插值計(jì)算耗時(shí)及占計(jì)算總時(shí)長(zhǎng)百分比
Tab.1 Interpolation calculation accounts for the total length of calculation using different sites
站點(diǎn)數(shù)20080032001280051200CPU耗時(shí)/ms1.382.4718.676.9269.4占計(jì)算總時(shí)長(zhǎng)百分比/%6177879195
如上節(jié)所述,廣義質(zhì)量的計(jì)算包含兩次長(zhǎng)度為N的向量元素相乘和長(zhǎng)度為N的向量歸約問(wèn)題[14]。在本文中每個(gè)線程完成元素相乘的工作后,選用交錯(cuò)配對(duì)的歸約算法[15]。交錯(cuò)配對(duì)歸約算法與鄰近配對(duì)歸約算法相比,其計(jì)算的跨越步長(zhǎng)不再是相鄰線程的數(shù)據(jù),而是為線程塊大小的一半起始,每次迭代減少一半,因此與鄰近配對(duì)歸約算法相比,交錯(cuò)配對(duì)算法每次迭代減少的工作線程不發(fā)生變化,使得盡可能的減少了線程束的失速。經(jīng)測(cè)試表明,交錯(cuò)配對(duì)歸約算法較鄰近配對(duì)歸約算法,在不同計(jì)算規(guī)模下有1.34倍~1.69倍的加速比。
針對(duì)此類問(wèn)題,通常用八叉樹(shù)描述三維空間的樹(shù)狀數(shù)據(jù)結(jié)構(gòu):樹(shù)中任意節(jié)點(diǎn)的子節(jié)點(diǎn)數(shù)要么為零要么為八,通過(guò)增加每個(gè)節(jié)點(diǎn)指針的數(shù)量(由兩個(gè)變?yōu)榘嘶蚓艂€(gè))和每次索引后的計(jì)算量(由左右邊界的判斷變?yōu)槿较蛏线吔缗c中點(diǎn)的判斷)減少索引深度,以達(dá)到提高速度的目的。而NVIDIA公司提供了動(dòng)態(tài)并行API(CDP,CUDA Dynamic Parallelism API)并演示了通過(guò)遞歸的方式構(gòu)建四叉樹(shù)的示例,故本文以此為模板構(gòu)建氣動(dòng)數(shù)據(jù)表的八叉樹(shù)。其GPU內(nèi)動(dòng)態(tài)創(chuàng)建線程塊的流程圖如圖2所示:在獲取數(shù)據(jù)表數(shù)據(jù)大小后確認(rèn)八叉樹(shù)的深度;待滿足建表要求后,將父線程塊一分為八并以線程束為基本單位統(tǒng)計(jì)每個(gè)子塊內(nèi)的點(diǎn)數(shù);待線程束完成其計(jì)算后將相關(guān)數(shù)據(jù)傳與八個(gè)子線程塊,以此遞歸直至不滿足條件為止。
圖2 動(dòng)態(tài)構(gòu)建線程塊示意圖
國(guó)內(nèi)外針對(duì)多機(jī)并行計(jì)算時(shí)多選用MPI實(shí)現(xiàn)數(shù)據(jù)的交互及同步[18],但在單機(jī)多GPU算例中未詳細(xì)闡述針對(duì)所研究對(duì)象的程序架構(gòu)優(yōu)化問(wèn)題。在以上文獻(xiàn)中均使用的同步函數(shù),使得設(shè)備在完成任務(wù)前不會(huì)將控制權(quán)交還主機(jī)線程,導(dǎo)致主機(jī)線程的阻塞。故在此可通過(guò)粗粒度的CPU和細(xì)粒度的GPU兩方面著手實(shí)現(xiàn)并行程序的性能優(yōu)化:前者實(shí)則為程序主架構(gòu)的優(yōu)化,通過(guò)開(kāi)啟一個(gè)有著若干線程的線程池,初始化階段構(gòu)建一個(gè)任務(wù)結(jié)構(gòu)體,內(nèi)含CUDA計(jì)算的識(shí)別號(hào)、所使用的設(shè)備號(hào)、流指針、主機(jī)端數(shù)據(jù)指針及設(shè)備端數(shù)據(jù)指針,在程序初始化時(shí)根據(jù)當(dāng)前所擁有的設(shè)備數(shù)將計(jì)算任務(wù)分發(fā)給不同的線程,運(yùn)用異步方式始終保證CPU的控制權(quán),直至GPU運(yùn)算完成觸發(fā)回調(diào)函數(shù)即可按邏輯進(jìn)行下一步的計(jì)算;而后者則可通過(guò)某個(gè)計(jì)算能力2.x以上的GPU在執(zhí)行計(jì)算時(shí)同時(shí)進(jìn)行數(shù)據(jù)的傳輸,使GPU得執(zhí)行引擎與存儲(chǔ)引擎可同時(shí)工作[19]。
鑒于在本文數(shù)值解算過(guò)程中彈性振動(dòng)方程解算模塊解耦為軸法橫三個(gè)方向的獨(dú)立計(jì)算,且沒(méi)有相互耦合的中間變量,程序主架構(gòu)的設(shè)計(jì)如圖3所示。主線程維護(hù)一個(gè)先進(jìn)先出的隊(duì)列完成任務(wù)分配并追蹤線程池內(nèi)從線程的任務(wù)進(jìn)度和狀態(tài),而線程池內(nèi)從線程負(fù)責(zé)完成實(shí)際的GPU函數(shù)調(diào)用并產(chǎn)生局部解算結(jié)果。
圖3 程序多線程異步架構(gòu)示意圖
而圖4演示了某個(gè)線程內(nèi)GPU運(yùn)算的時(shí)序?qū)Ρ?。GPU采用同步運(yùn)行時(shí),一旦傳輸數(shù)據(jù)量較大可能導(dǎo)致GPU設(shè)備的空置。若像圖中在初始化階段進(jìn)行數(shù)據(jù)分割,即將某一方向中的原始數(shù)據(jù)劃分為若干個(gè)小塊,以流的形式可以實(shí)現(xiàn)數(shù)據(jù)傳輸?shù)碾[藏以提高計(jì)算效率。同時(shí)自Kepler架構(gòu)后HyperQ技術(shù)使得運(yùn)行多個(gè)CPU核或線程在某一GPU上啟動(dòng)任務(wù),以多個(gè)工作隊(duì)列的形式提高了GPU的利用率,避免了偽依賴。其偽代碼如下
圖4 串行與異步流形式運(yùn)算時(shí)序?qū)Ρ葓D
1) CreateStream
2) While i 3) cudaHostAlloc && cudamalloc 4) cudaMemcpyAsync(H2D) 5) interpolate<< 6) cudaMemcpyAsync(D2H) 7) cudaStreamAddCallback(stream,callbackfunc, data,flag) 8)i++, return 3) 上述第1)步首先依照任務(wù)結(jié)構(gòu)體內(nèi)的流指針創(chuàng)建流對(duì)象,并得到計(jì)算設(shè)備其內(nèi)部的各種硬件參數(shù);第2)步~第8)步依照給定的流對(duì)象完成內(nèi)存拷貝、核函數(shù)計(jì)算、函數(shù)回調(diào)等待數(shù)據(jù)回傳的過(guò)程。其中第5)步的第三個(gè)模板參數(shù)即低訪問(wèn)延遲的共享內(nèi)存的大小,在核函數(shù)內(nèi)通過(guò)extern __shared__關(guān)鍵詞動(dòng)態(tài)分配共享內(nèi)存空間,以起到節(jié)省每個(gè)流多處理器(SM,Streaming Multiprocessors)上寶貴的一級(jí)緩存空間,防止核函數(shù)調(diào)用失敗。鑒于本文所用最大的三維氣動(dòng)插值表不超過(guò)1 200個(gè)數(shù)據(jù)點(diǎn),因此可通過(guò)將線性化的氣動(dòng)數(shù)據(jù)表放入共享內(nèi)存中以實(shí)現(xiàn)讀取氣動(dòng)數(shù)據(jù)的加速。 為使程序在任意復(fù)雜度的數(shù)值算法中保證高性能的通適性,需考慮第5)步內(nèi)前兩個(gè)參數(shù)對(duì)硬件的利用率有影響:一個(gè)流多處理器在線程數(shù)、線程塊數(shù)及寄存器數(shù)方面皆有硬件限制,具體參數(shù)因GPU計(jì)算能力不同而異。在計(jì)算模式和算法已定的情況下,每個(gè)線程內(nèi)的寄存器使用數(shù)量即已確定,而所需分配共享內(nèi)存的大小直接取決于氣動(dòng)數(shù)據(jù)表的大小,故最為影響流多處理器理論效率的即為核函數(shù)的第二個(gè)參數(shù),即每個(gè)線程塊內(nèi)的線程數(shù)。依照參考文獻(xiàn)[20]建立流多處理器利用率和共享內(nèi)存大小、線程塊內(nèi)線程束的函數(shù),在程序初始化階段依照該函數(shù)關(guān)系建立效率索引表,如圖5所示為單線程使用25個(gè)寄存器的情況下每個(gè)流多處理器的理論效率索引圖。故在程序中可依照算法和氣動(dòng)數(shù)據(jù)表的大小決定參數(shù)2,以期實(shí)現(xiàn)利用率的最大化及核函數(shù)的動(dòng)態(tài)參數(shù)優(yōu)化。 本文建立單臺(tái)計(jì)算機(jī)、多GPU并行求解器,模擬細(xì)長(zhǎng)體彈性飛行器飛行過(guò)程進(jìn)行仿真,并對(duì)單GPU及多GPU并行的性能及拓展性進(jìn)行對(duì)比。因本文重點(diǎn)為加速計(jì)算的方法研究,暫選用某型飛行器的氣動(dòng)數(shù)據(jù),但鑒于數(shù)據(jù)敏感且不失一般性,計(jì)算時(shí)每個(gè)時(shí)步較上一個(gè)時(shí)步在氣動(dòng)系數(shù)的各個(gè)維度上均添加一個(gè)隨機(jī)量。受計(jì)算資源限制,算例平臺(tái)硬件采用主頻3.4的Intel E3 1230V5,兩片英偉達(dá)GTX1060 3G顯卡;編程環(huán)境為VS2017下的MSVC14.0編譯器及CUDA 9.0的NVCC編譯器。時(shí)間統(tǒng)計(jì)時(shí)采用CUDA自帶的cudaEvenet_t計(jì)時(shí),其精度為0.5 μs;考慮彈性模塊與其他模塊的數(shù)據(jù)耦合關(guān)系,以計(jì)算最頻繁的控制周期為實(shí)時(shí)性判定準(zhǔn)則,在此選用10 ms為限。 圖5 流多處理器利用率索引示意圖 程序的正確性及快速性是衡量程序性能的重要指標(biāo)。選用某型細(xì)長(zhǎng)體飛行器的攻角作為標(biāo)準(zhǔn),模擬起飛時(shí)刻至第一級(jí)分離過(guò)程彈性模塊的計(jì)算流程,如圖6所示,剛體飛行器與彈性體飛行器的彈道傾角變化總體一致且相差不大。其原因是在傳統(tǒng)剛體飛行器的彈道仿真中,由于控制系統(tǒng)的魯棒性,在有一定偏差條件下仍能夠滿足飛行姿態(tài)的穩(wěn)定并保證對(duì)標(biāo)準(zhǔn)彈道的跟蹤,故彈性變形引起的各附加量輸出至指導(dǎo)控制計(jì)算機(jī)后對(duì)整體彈道的影響不大。而如圖7所示,通過(guò)傳統(tǒng)CPU與GPU運(yùn)算結(jié)果對(duì)比,可認(rèn)為CPU與GPU的運(yùn)算結(jié)果并無(wú)二致。 由圖8可知,單GPU情況下在站點(diǎn)數(shù)較少時(shí)加速效果并不明顯,因?yàn)橛?jì)算均需通過(guò)低速的PCI-E總線傳輸計(jì)算數(shù)據(jù),且啟動(dòng)的線程數(shù)不足以隱藏?cái)?shù)據(jù)傳輸和核函數(shù)啟動(dòng)的時(shí)間開(kāi)銷,即異構(gòu)架構(gòu)帶來(lái)額外的通信及函數(shù)啟動(dòng)時(shí)間無(wú)法抵消多核運(yùn)算帶來(lái)的性能提升。但盡管選用常規(guī)的Geforce卡,其沒(méi)有ECC顯存且每個(gè)流多處理器內(nèi)的雙精度計(jì)算單元比例更少,采用CUDA技術(shù)的雙精度計(jì)算的求解時(shí)間在單元個(gè)數(shù)較多時(shí)仍能保證較好的加速比,經(jīng)優(yōu)化后在站點(diǎn)數(shù)25 000個(gè)左右時(shí)能保證仿真計(jì)算的實(shí)時(shí)性;隨站點(diǎn)數(shù)的增加,加速比達(dá)到13左右且斜率下降趨于穩(wěn)定。而當(dāng)氣動(dòng)數(shù)據(jù)表較小足以放入共享內(nèi)存中時(shí),同理在站點(diǎn)數(shù)較少加速效果不明顯,隨著站點(diǎn)數(shù)的增加,計(jì)算由I/O密集型轉(zhuǎn)向計(jì)算密集型,加速效果逐漸明顯;與常規(guī)從全局顯存讀取數(shù)據(jù)的加速比不同,在站點(diǎn)數(shù)3 000左右有一個(gè)下跌的趨勢(shì),其原因可能是多線程訪問(wèn)共享內(nèi)存導(dǎo)致的訪存沖突進(jìn)而影響了運(yùn)算效率。 圖6 彈性體與剛形體彈道傾角對(duì)比 Fig.6 Comparison of trajectory inclination angle between rigid and elastic vehicle 圖7 CPU與GPU運(yùn)算誤差結(jié)果對(duì)比 Fig.7 Comparison of error results between CPU and GPU 圖8 不同站點(diǎn)數(shù)CPU與GPU計(jì)算結(jié)果對(duì)比 Fig.8 Comparison of the calculate results between CPU and GPU with different scales 圖9 不同站點(diǎn)數(shù)多GPU同/異步架構(gòu)耗時(shí)對(duì)比 Fig.9 Time consuming comparison of multi GPU synchronous/asynchronous architecture with different scales 而異步較之于同步設(shè)計(jì)、多GPU較之于單GPU在數(shù)據(jù)規(guī)模較小時(shí)差異不大,甚至效率更低,如圖9所示。其原因是多GPU情況需要進(jìn)行數(shù)據(jù)分割和開(kāi)啟額外的顯卡資源。此時(shí)異步操作更需要線程同步等待,加速比更為不理想;隨站點(diǎn)數(shù)的增加,多GPU加速比提升更快,加速效果更明顯,直至50 000個(gè)站點(diǎn)時(shí)也未出現(xiàn)加速比趨于平緩的趨勢(shì),其原因是數(shù)據(jù)分割時(shí)間所占總時(shí)長(zhǎng)百分比愈來(lái)愈小,并行部分的傳輸較計(jì)算操作在程序中所占比例逐步減少得到充分的延遲隱藏,使得流多處理器效率更高,也從側(cè)面映證該方案適宜求解大規(guī)模的數(shù)值運(yùn)算。 (1) 基于CUDA技術(shù)實(shí)現(xiàn)了細(xì)長(zhǎng)體飛行器彈性模塊的快速求解方法設(shè)計(jì),在多GPU環(huán)境下取截?cái)嗄B(tài)40階,可保證至少1 200個(gè)站點(diǎn)的彈性飛行器的實(shí)時(shí)仿真計(jì)算。 (2) 針對(duì)GPU架構(gòu)提出了動(dòng)態(tài)并行分配線程塊的八叉樹(shù)的氣動(dòng)系數(shù)索引方法,優(yōu)化了雙線性插值的索引步驟,效果較好。 (3) 設(shè)計(jì)了單機(jī)多GPU的異步并行架構(gòu),并仿真說(shuō)明多GPU并行能夠進(jìn)一步顯著提高計(jì)算效率,證明了所建并行架構(gòu)的擴(kuò)展性與快速性。在此研究基礎(chǔ)之上,開(kāi)展針對(duì)超過(guò)三維的插值表優(yōu)化索引模式研究,對(duì)多機(jī)協(xié)同計(jì)算做進(jìn)一步的可拓展性分析。3 數(shù)值算例分析
3.1 測(cè)試平臺(tái)
3.2 計(jì)算精度及效率分析
4 結(jié) 論