吳錦濤, 呂一品, 董 剛, 王 帥, 唐 科

(1.北京宇航系統(tǒng)工程研究所, 北京 100076; 2.北京電子工程總體研究所, 北京 100854; 3.南京理工大學(xué)瞬態(tài)物理重點(diǎn)實(shí)驗(yàn)室, 江蘇南京 210094)

引 言

高超聲速的可燃?xì)饬髟谡系K物作用下可形成駐定于障礙物的爆轟, 即斜爆轟(oblique detona-tion wave, ODW).斜爆轟的本質(zhì)是激波和對流預(yù)混火焰組合成的復(fù)合波, 這種由激波主導(dǎo)的燃燒, 具有很高的燃燒效率和傳播速度(達(dá)幾千米每秒).作為一種新型的能量轉(zhuǎn)換模式, 研究斜爆轟為發(fā)展高效的熱力推進(jìn)(如斜爆轟發(fā)動(dòng)機(jī),(oblique detona-tion wave engine, ODWE))提供了技術(shù)支持[1-3], 因此開展相關(guān)研究工作具有重要的科學(xué)價(jià)值和應(yīng)用前景.

為精確捕捉斜爆轟波的流場特性, 數(shù)值模擬時(shí)通常采用詳細(xì)的化學(xué)反應(yīng)機(jī)理對反應(yīng)源項(xiàng)進(jìn)行計(jì)算, 因此, 對源項(xiàng)的計(jì)算過程轉(zhuǎn)化為求解多個(gè)由反應(yīng)速率方程組成的強(qiáng)非線性常微分方程組.由于不同組分的生成/消耗的速率(反應(yīng)時(shí)間尺度)相差很大(實(shí)際燃燒和爆轟問題中尺度相差可達(dá)9個(gè)數(shù)量級(jí)[4]), 對該方程組的直接積分(direct integration, DI)求解呈現(xiàn)出強(qiáng)烈的剛性現(xiàn)象, 并花費(fèi)大量的計(jì)算用時(shí), 因而需要尋找減少化學(xué)反應(yīng)計(jì)算時(shí)間的方法.在此基礎(chǔ)上, 通過合理的約束條件對詳細(xì)化學(xué)反應(yīng)機(jī)理進(jìn)行簡化的方法應(yīng)運(yùn)而生[5-6], 然而, 簡化機(jī)理的約束條件并不適用于燃燒和爆轟問題中所有可能的化學(xué)熱力學(xué)狀態(tài), 可移植性較差.此外, 流場中的部分痕量物質(zhì)(如某些代表化學(xué)反應(yīng)活性的自由基)的特征時(shí)間尺度通常較小, 采用簡化機(jī)理方法時(shí), 這些組分會(huì)被認(rèn)為滿足約束條件而刪減, 導(dǎo)致求解過程的失真.

基于此, 通過改進(jìn)計(jì)算方法以在采用詳細(xì)機(jī)理的同時(shí)提升化學(xué)反應(yīng)過程的計(jì)算效率, 成為避免直接積分求解剛性方程組的又一途徑, 動(dòng)態(tài)自適應(yīng)建表(in situ adaptive tabulation, ISAT)方法[7]就是其中之一.與簡化化學(xué)反應(yīng)機(jī)理的方法不同, ISAT方法通過在計(jì)算機(jī)內(nèi)存中建立化學(xué)反應(yīng)數(shù)據(jù)表的方式, 將流場內(nèi)出現(xiàn)的化學(xué)動(dòng)/熱力學(xué)狀態(tài)以數(shù)據(jù)節(jié)點(diǎn)的形式進(jìn)行插入、問詢和取回等操作, 取代原直接積分求解.由于對數(shù)據(jù)節(jié)點(diǎn)的存/取操作耗時(shí)遠(yuǎn)遠(yuǎn)小于求解剛性方程組的耗時(shí), ISAT方法可以很大程度地減少化學(xué)反應(yīng)過程的計(jì)算時(shí)間.同時(shí), 由于動(dòng)態(tài)建表技術(shù)對內(nèi)存的占用需求不大, 該方法具有很好的可移植性.在低Mach數(shù)條件下的定常燃燒數(shù)值模擬中, 考慮到流場中一個(gè)時(shí)間步長內(nèi)化學(xué)反應(yīng)熱力學(xué)狀態(tài)的變化較小, 數(shù)據(jù)表內(nèi)的絕大多數(shù)節(jié)點(diǎn)可以被成功問詢并取回利用, 因此, 在這類問題中ISAT方法的計(jì)算效率很高.在均質(zhì)壓燃(homogeneous charge compression ignition, HCCI)發(fā)動(dòng)機(jī)的計(jì)算中, ISAT方法的化學(xué)反應(yīng)加速比可達(dá)100倍; 而在配對攪拌發(fā)生器(pairwise mixing stirred reactor, PMSR)的計(jì)算中, 采用ISAT方法計(jì)算化學(xué)反應(yīng)的時(shí)間僅為DI的1/42.

然而, 原始的ISAT方法主要局限于低Mach數(shù)下的穩(wěn)態(tài)燃燒問題和均相問題的計(jì)算中, 其在高M(jìn)ach數(shù)條件下可壓縮流的非定常燃燒/爆轟問題中的性能表現(xiàn)仍有待檢驗(yàn).事實(shí)上, 可壓縮流場中一個(gè)時(shí)間步長前后化學(xué)熱力學(xué)狀態(tài)的變化較大, 數(shù)據(jù)表內(nèi)的大量節(jié)點(diǎn)無法被成功問詢和取回, 導(dǎo)致ISAT方法的計(jì)算效率下降.同時(shí), 數(shù)據(jù)表內(nèi)不斷插入的冗余節(jié)點(diǎn)(未被取回利用)勢必導(dǎo)致表容量的過度增長, 進(jìn)而突破計(jì)算機(jī)內(nèi)存上限, 造成計(jì)算過程的中斷.為了避免這種情況, Dong等在ISAT方法的基礎(chǔ)上作了進(jìn)一步改進(jìn), 即, 添加了整體刪除和節(jié)點(diǎn)刪除兩種刪除機(jī)制, 并基于此提出一種動(dòng)態(tài)存儲(chǔ)/刪除建立熱力學(xué)數(shù)據(jù)表算法[8-9].在激波引發(fā)CH4/O2預(yù)混氣體爆轟問題的二維數(shù)值計(jì)算中, 動(dòng)態(tài)存儲(chǔ)/刪除算法計(jì)算化學(xué)反應(yīng)源項(xiàng)的加速比達(dá)17.88, 總加速比達(dá)4.75.

通過上述研究發(fā)現(xiàn), ISAT方法在使用詳細(xì)化學(xué)反應(yīng)機(jī)理的數(shù)值模擬中, 具有較好的可移植性, 但其計(jì)算效率因物理問題而異.隨著近年來大型計(jì)算設(shè)備的發(fā)展, 大規(guī)模并行計(jì)算日益盛行, 而以ISAT為代表的建表方法也不斷被應(yīng)用于并行數(shù)值模擬中.不同于串行執(zhí)行方式, 在并行執(zhí)行方式(以MPI并行通訊協(xié)議為例)的基礎(chǔ)上, 數(shù)值計(jì)算的每個(gè)分區(qū)都相當(dāng)于一個(gè)獨(dú)立的計(jì)算區(qū)域, 如何在并行環(huán)境中生成和構(gòu)建數(shù)據(jù)表、并完成各數(shù)據(jù)表之間的協(xié)同合作, 是建表技術(shù)在并行計(jì)算中應(yīng)用的一個(gè)關(guān)鍵.Lu等[10]首先開展了基于ISAT建表技術(shù)的并行計(jì)算策略的研究, 在其提出的所有策略中, 最直觀的策略是在每個(gè)MPI分區(qū)內(nèi)各自獨(dú)立地生成數(shù)據(jù)表, 計(jì)算過程中, 每個(gè)分區(qū)的化學(xué)熱力學(xué)狀態(tài)均在各自對應(yīng)的數(shù)據(jù)表內(nèi)進(jìn)行插入、問詢和取回, 彼此之間沒有數(shù)據(jù)傳遞.這種并行策略被稱之為純粹本地建表(purely local processing, PLP).此外, Lu等還提出了快速問詢(quick try, QT)、均勻隨機(jī)分布(uniform random distribution, URAN)和優(yōu)先分配(preferential distribution, PREF)等策略, 并將其應(yīng)用在CH4/O2預(yù)混氣體多重混合攪拌發(fā)生器的計(jì)算中.然而, 針對不同的計(jì)算工況, 尤其是流場內(nèi)熱力學(xué)狀態(tài)分布不均勻且隨時(shí)空變化較大的計(jì)算工況, 這些并行策略在計(jì)算效率上的表現(xiàn)并不令人滿意.這進(jìn)一步說明, 在并行計(jì)算中采用基于建表的化學(xué)加速方法, 其計(jì)算效率同樣取決于具體研究的物理問題, 針對二維斜爆轟這類流場內(nèi)熱力學(xué)狀態(tài)分布嚴(yán)重不均勻且隨時(shí)空劇烈變化的反應(yīng)流計(jì)算, Lu等[10]提出的策略顯然有待進(jìn)一步的檢驗(yàn).

不同于穩(wěn)態(tài)燃燒問題, 并行計(jì)算中爆轟波流場的各個(gè)分區(qū)的計(jì)算負(fù)擔(dān)往往差異很大(如由波陣面組成的間斷區(qū)和由波后燃燒產(chǎn)物組成的光滑區(qū)), 因而在提出并行策略時(shí), 應(yīng)考慮如何建立與并行分區(qū)相應(yīng)的數(shù)據(jù)表, 以保證計(jì)算過程中不同數(shù)據(jù)表操作之間的負(fù)載平衡性.此外, 由于爆轟問題的特殊性, 計(jì)算時(shí)數(shù)據(jù)表中不斷有新節(jié)點(diǎn)插入, 如何將數(shù)據(jù)表刪除機(jī)制[8-9]與并行策略相結(jié)合, 以在提升計(jì)算效率的同時(shí)避免因某個(gè)數(shù)據(jù)表容量突破內(nèi)存上限導(dǎo)致計(jì)算中斷.

本文針對二維2H2+O2斜爆轟數(shù)值模擬, 提出了幾種不同類型的并行策略, 并將這些并行策略與帶節(jié)點(diǎn)刪除機(jī)制的動(dòng)態(tài)存儲(chǔ)/刪除建立熱力學(xué)數(shù)據(jù)表方法相結(jié)合, 從而形成了適用于瞬態(tài)可壓縮反應(yīng)流問題的并行建表化學(xué)加速算法, 并從計(jì)算精度和計(jì)算效率兩方面分析了并行策略選取對算法性能的影響.

1 計(jì)算方法

1.1 基于ISAT的動(dòng)態(tài)存儲(chǔ)/刪除建表

考慮帶K種不同組分的可壓縮反應(yīng)流流場, 其化學(xué)熱力學(xué)狀態(tài)向量表達(dá)如下

Φ=(Y1,…,Yk,…,YK,T,ρ)

=(φ1, …,φk, …,φK,φK+1,φK+2)

(1)

式中,Yk(k=1, …,K)是組分k的質(zhì)量分?jǐn)?shù),T和ρ分別是反應(yīng)體系的溫度和密度.在燃燒過程中, 流場內(nèi)某單元的Φ0在固定步長Δt內(nèi)因化學(xué)反應(yīng)帶來的狀態(tài)變化為R(Φ), 具體表達(dá)形式為

(2)

式中,Φ0和Φ1分別代表更新前后的熱力學(xué)狀態(tài),S為化學(xué)反應(yīng)源項(xiàng).式(2)給出了狀態(tài)方程的DI求解過程.在爆轟等可壓縮反應(yīng)流的計(jì)算過程中, 由于不同組分生成/消耗的速率相差很大, 其參與反應(yīng)進(jìn)程的時(shí)間尺度分布范圍在10-9～1 s之間, 這使得S具有強(qiáng)剛性.在對化學(xué)反應(yīng)過程進(jìn)行DI求解時(shí), 耗時(shí)巨大.

不同于直接積分求解, 原始的ISAT方法通過對表中數(shù)據(jù)的問詢/取回操作替代ODE的直接積分求解過程, 從而獲取化學(xué)變化的近似解.考慮流場中一個(gè)待問詢的化學(xué)熱力學(xué)狀態(tài)矢量Φq, 如果數(shù)據(jù)表中存在某一狀態(tài)Φ(j)的節(jié)點(diǎn)與其足夠接近, 那么, 在實(shí)際ISAT計(jì)算過程中,Φq所對應(yīng)的化學(xué)變化R(Φq)可由存儲(chǔ)在數(shù)據(jù)表中的Φ(j)和對應(yīng)的R(Φ(j))通過Taylor級(jí)數(shù)展開獲得, 形式如下

R(Φq)=

(3)

其中, ΔΦ=Φq-Φ(j).式(3)對?R(Φ(j))/(?Φ(j))的計(jì)算和存儲(chǔ)產(chǎn)生了節(jié)點(diǎn)j狀態(tài)的梯度映射矩陣A, 通過對矩陣A進(jìn)行容差判據(jù)(εtol)的設(shè)定進(jìn)而可以構(gòu)造一個(gè)以Φ(j)為中心的超橢誤差域(ellipsoid of accuracy, EOA).如果狀態(tài)矢量Φq落在Φ(j)所對應(yīng)的EOA范圍內(nèi), 則R(Φq)將通過式(3)獲得的R(Φ(j))線性近似值, 這一過程稱為取回.由于對數(shù)據(jù)表的問詢/取回操作所耗費(fèi)的時(shí)間明顯小于直接積分式(2)求解的計(jì)算用時(shí), 因此, ISAT能夠起到加速化學(xué)反應(yīng)計(jì)算的效果.

值得注意的是, 原始ISAT方法采用線性近似的目的在于讓數(shù)據(jù)表內(nèi)節(jié)點(diǎn)所存儲(chǔ)的熱力學(xué)狀態(tài)盡可能接近直接積分的結(jié)果, 然而, 當(dāng)流場內(nèi)化學(xué)反應(yīng)過程十分劇烈(如二維斜爆轟流場)時(shí), 設(shè)定EOA判據(jù)會(huì)導(dǎo)致表中節(jié)點(diǎn)難以在后續(xù)計(jì)算時(shí)被成功問詢并取回, 進(jìn)而影響其計(jì)算效率.Veljkovic[11]對比了線性近似和常數(shù)近似的計(jì)算結(jié)果發(fā)現(xiàn), 采用零階精度的常數(shù)近似方法, 可以在保證計(jì)算精度的基礎(chǔ)上改善建表算法在瞬態(tài)問題中的計(jì)算效率.因此, 本文所采用的動(dòng)態(tài)存儲(chǔ)/刪除建立熱力學(xué)數(shù)據(jù)表算法使用了常數(shù)近似取代原ISAT方法中的線性近似, 表達(dá)如下

R(Φq)=R(Φ(j))+O(|ΔΦ|)

(4)

在此基礎(chǔ)上, 本文重新定義了關(guān)于節(jié)點(diǎn)Φ(j)的誤差范圍, 如下

(5)

式中,e+(Φ(j))和e-(Φ(j))分別代指熱力學(xué)狀態(tài)Φ(j)的誤差上限和誤差下限, 等式右邊為狀態(tài)空間內(nèi)各分量的誤差上下限.上述定義實(shí)際上構(gòu)造了一個(gè)超矩形誤差域(hyper-rectangle of accuracy, ROA), 并且避免了式(3)中梯度映射矩陣A的計(jì)算和存儲(chǔ), 計(jì)算效率因此得到進(jìn)一步的提升.此外, 在包含間斷的瞬態(tài)可壓縮反應(yīng)流問題中, 隨著時(shí)間的進(jìn)行, 流場內(nèi)的熱力學(xué)狀態(tài)變化巨大, 表中的大量數(shù)據(jù)無法被成功取回, 這勢必導(dǎo)致數(shù)據(jù)表內(nèi)節(jié)點(diǎn)數(shù)量的急劇上升, 甚至可能會(huì)達(dá)到計(jì)算機(jī)內(nèi)存的使用上限并導(dǎo)致計(jì)算中斷.針對這一問題, Dong等[8-9]在改進(jìn)后的ISAT算法中引進(jìn)了節(jié)點(diǎn)刪除機(jī)制, 有效地控制了數(shù)據(jù)表的生長.基于以上對ISAT方法的改進(jìn), 本文動(dòng)態(tài)存儲(chǔ)/刪除建表算法的執(zhí)行流程可以描述如下:

(1)計(jì)算開始之前, 首先在計(jì)算機(jī)內(nèi)存中建立一個(gè)空的數(shù)據(jù)表.表中第1個(gè)節(jié)點(diǎn)j的熱力學(xué)狀態(tài)矢量由式(2)直接積分而來, 隨后, 節(jié)點(diǎn)j將記錄并保存該狀態(tài)矢量Φ(j)及其對應(yīng)的化學(xué)變化R(Φ(j))和初始誤差數(shù)組上下限e±(Φ(j)).

(2)計(jì)算開始后, 對流場中每一個(gè)待問詢的狀態(tài)矢量Φq, 都首先在數(shù)據(jù)表中進(jìn)行查找和比較, 如果表中存在某個(gè)節(jié)點(diǎn)Φ(j), 使得Φq中的所有分量都能同時(shí)滿足

(k=1,…,K,K+1,K+2)

(6)

則執(zhí)行式(4)中的常數(shù)近似操作, 將數(shù)據(jù)表中存儲(chǔ)的R(Φ(j))取回, 作為Φq化學(xué)反應(yīng)更新后的變化量, 即R(Φq).這種情況被稱為取回成功.

(3)如果式(6)中的條件不能得到滿足, 需要通過式(2)中的直接積分求解獲取Φq的新的化學(xué)變化量R(Φq).同時(shí), 記錄并選取在Φq查找路徑上的與之最接近的狀態(tài)矢量Φ(n), 將積分求得的R(Φq)和最接近節(jié)點(diǎn)中存儲(chǔ)的R(Φ(n))所包含的所有狀態(tài)分量進(jìn)行比較, 形式如下

(7)

對節(jié)點(diǎn)數(shù)據(jù)中的ROA誤差域執(zhí)行放大操作是動(dòng)態(tài)存儲(chǔ)/刪除建表算法中的關(guān)鍵步驟.雖然在當(dāng)前時(shí)間步內(nèi), 式(7)中的放大操作無法保證待問詢的狀態(tài)矢量Φq從表中成功取回, 但它能夠提高之后的時(shí)間步中狀態(tài)矢量Φ(n)在問詢過程中被成功取回的概率.假設(shè)下一步計(jì)算中某個(gè)待問詢的狀態(tài)矢量為Φq′, 如果它能夠與在上一步中經(jīng)過放大ROA范圍得到的狀態(tài)矢量(更新后的Φ(n))滿足式(6)中的條件, 這意味著

Φq′=Φ(n)+R(Φq)

(8)

另一方面, 對節(jié)點(diǎn)n而言, 經(jīng)過化學(xué)反應(yīng)更新后的狀態(tài)矢量Φn′又滿足

Φn′=Φ(n)+R(Φ(n))

(9)

因此, 結(jié)合式(8),(9)可以發(fā)現(xiàn), 組成Φq′和Φn′的所有狀態(tài)分量共同構(gòu)成放大后的ROA誤差域的邊界, 即

(10)

式(5)中的誤差數(shù)組上下限可以進(jìn)一步表示成

(11)

其中, 相對誤差和絕對誤差的初始值在本文中分別給定為εr=1×10-3和εa=1×10-8.當(dāng)問詢矢量滿足式(5)的條件時(shí), 節(jié)點(diǎn)數(shù)據(jù)在ROA誤差域內(nèi)的取回精度由εr的值決定.否則, ROA誤差范圍將被一定程度地放大, 因而在隨后的計(jì)算中, 經(jīng)過放大操作的ROA誤差域所對應(yīng)的新的取回精度將由式(10)中εtol而非εr決定.換句話說, 數(shù)據(jù)表中節(jié)點(diǎn)的初始εr越小, 計(jì)算過程中可能發(fā)生的放大操作越多, 因而該節(jié)點(diǎn)實(shí)際的取回精度最終將取決于εtol的設(shè)定; 相反, 數(shù)據(jù)表中節(jié)點(diǎn)的初始εr越大, 對ROA的放大操作發(fā)生的次數(shù)就越少, 因而該節(jié)點(diǎn)內(nèi)部所存儲(chǔ)的熱力學(xué)狀態(tài)將主要取決于初始εr的設(shè)定.

實(shí)際計(jì)算過程中, 通常會(huì)設(shè)定相同的εr和εtol, 以保證節(jié)點(diǎn)參數(shù)的取回精度不受放大操作的影響.不同的是,εr的值在節(jié)點(diǎn)插入數(shù)據(jù)表的時(shí)刻已經(jīng)給定, 因而只能局部地影響該節(jié)點(diǎn)的數(shù)據(jù)取回精度; 而εtol可以在整個(gè)計(jì)算過程中動(dòng)態(tài)地控制數(shù)據(jù)節(jié)點(diǎn)的ROA誤差范圍, 因而具有全局的影響.考慮到εtol對動(dòng)態(tài)存儲(chǔ)/刪除建表算法精度有明顯影響, Dong等[8]在激波聚焦引發(fā)爆轟波起爆的數(shù)值模擬中, 研究了不同容差判據(jù)(εtol)下建表加速算法的表現(xiàn).結(jié)果表明, 一味地增大εtol雖然能夠明顯地提升計(jì)算效率, 但勢必會(huì)導(dǎo)致計(jì)算精度的下降.為同時(shí)兼顧建表加速算法在計(jì)算精度和計(jì)算效率上的表現(xiàn), 本文將容差判據(jù)設(shè)定為εtol=1×10-3, 該大小的容差判據(jù)在瞬態(tài)可壓縮反應(yīng)流問題[8-9]中被證明是可靠的.另外, 通過設(shè)定絕對誤差(εa=1×10-8)保證數(shù)據(jù)節(jié)點(diǎn)中狀態(tài)矢量各個(gè)分量φ(j)的非負(fù)性.在實(shí)際計(jì)算過程中, 對任意分量k, 需要滿足φk(j)≥εa/(1-εr)≈1.001×10-8, 因此, 本文設(shè)定一個(gè)足夠小的下限指標(biāo)1.005×10-8, 當(dāng)計(jì)算過程中φk(j)小于這個(gè)值時(shí), 令φk(j)=1.005×10-8, 以保證反應(yīng)系統(tǒng)中所有化學(xué)熱力學(xué)狀態(tài)分量在給定誤差數(shù)組上下限時(shí)的物理真實(shí)性.需要注意的是,εa的量綱隨著狀態(tài)分量的量綱而發(fā)生改變.當(dāng)φi(j)代表組分濃度時(shí)φa是一個(gè)無量綱量(φa=1×10-8); 當(dāng)φk(j)代表溫度時(shí),εa等于1×10-8K; 當(dāng)φk(j)代表密度時(shí),εa等于1×10-8kg/m3.

(4)針對包含復(fù)雜流場結(jié)構(gòu)的瞬態(tài)可壓縮反應(yīng)流問題(如本文中的氣相爆轟波傳播), 流場內(nèi)化學(xué)熱力學(xué)狀態(tài)隨時(shí)空的變化十分劇烈, 這會(huì)導(dǎo)致數(shù)據(jù)表中已插入的大量數(shù)據(jù)節(jié)點(diǎn)無法被成功取回利用.因此, 隨著計(jì)算的進(jìn)行, 不斷有新的節(jié)點(diǎn)會(huì)插入到數(shù)據(jù)表中并導(dǎo)致其過度增長.對于有限的計(jì)算機(jī)內(nèi)存空間而言, 持續(xù)增長的數(shù)據(jù)表尺寸最終將達(dá)到計(jì)算機(jī)內(nèi)存的使用上限, 從而導(dǎo)致整個(gè)計(jì)算過程被迫中斷.為了避免這一現(xiàn)象, 本文提出一種數(shù)據(jù)表尺寸控制策略, 該控制策略由單表尺寸(single table size,Msin)控制方法和總表尺寸(total table size,Mtot)控制方法組成.其中, 單表尺寸控制方法意味著通過單個(gè)數(shù)據(jù)表內(nèi)節(jié)點(diǎn)數(shù)量的刪除上限來控制表的尺寸, 計(jì)算過程中, 當(dāng)某個(gè)數(shù)據(jù)表中節(jié)點(diǎn)數(shù)量到達(dá)Msin時(shí), 節(jié)點(diǎn)刪除操作將被激活, 并檢查表內(nèi)已插入的所有數(shù)據(jù)節(jié)點(diǎn), 將之前的計(jì)算過程中未被重復(fù)取回利用的節(jié)點(diǎn)刪除.總表尺寸控制方法則是針對并行計(jì)算多數(shù)據(jù)表專門提出的控制方法, 其通過所有數(shù)據(jù)表內(nèi)節(jié)點(diǎn)數(shù)量總和的刪除上限來控制表的尺寸.在計(jì)算的每一個(gè)時(shí)間步內(nèi), 首先統(tǒng)計(jì)所有數(shù)據(jù)表中的節(jié)點(diǎn)數(shù)量總和, 當(dāng)總和大于Mtot時(shí), 節(jié)點(diǎn)刪除操作將在全局范圍內(nèi)激活, 并清除所有數(shù)據(jù)表內(nèi)未被重復(fù)取回的節(jié)點(diǎn).此外, 本文進(jìn)一步優(yōu)化了原節(jié)點(diǎn)刪除機(jī)制[12], 當(dāng)數(shù)據(jù)表內(nèi)節(jié)點(diǎn)數(shù)量到達(dá)上限(Msin或Mtot)后, 整個(gè)之前的計(jì)算過程(而非上一時(shí)間步)中被重復(fù)利用過的節(jié)點(diǎn)數(shù)據(jù), 都將被保留.優(yōu)化后的節(jié)點(diǎn)刪除機(jī)制雖然會(huì)導(dǎo)致數(shù)據(jù)表內(nèi)節(jié)點(diǎn)的累積, 但卻能夠進(jìn)一步減少重復(fù)的直接積分運(yùn)算.

1.2 并行策略

在基于MPI通訊協(xié)議的并行化學(xué)反應(yīng)流數(shù)值模擬中, 原本完整的計(jì)算區(qū)域被分成若干個(gè)計(jì)算分區(qū), 因此, 在數(shù)值模擬的過程中往往存在多個(gè)數(shù)據(jù)表, 每個(gè)數(shù)據(jù)表對其中的一個(gè)分區(qū)進(jìn)行操作, 故實(shí)現(xiàn)并行加速算法的重要手段是需要考慮并行計(jì)算中各個(gè)數(shù)據(jù)表的操作及其之間相互作用的方法或策略.其中, 最直觀和最易實(shí)現(xiàn)的一種并行策略是Lu等[10]提出的純粹本地建表(purely local process-ing, PLP)策略, 見圖1(a).在PLP策略中, 所有ISAT數(shù)據(jù)表均建立在其本身所對應(yīng)的計(jì)算分區(qū)中, 在實(shí)際計(jì)算過程中, 所有數(shù)據(jù)表操作只能在對應(yīng)的(即本地的)分區(qū)內(nèi)執(zhí)行, 彼此之間沒有數(shù)據(jù)傳遞.在有些燃燒問題中, 由于流場內(nèi)不同并行計(jì)算分區(qū)之間的熱力學(xué)狀態(tài)分布相似, 因此, 各分區(qū)內(nèi)的數(shù)據(jù)表所承擔(dān)的計(jì)算任務(wù)相當(dāng), 各個(gè)數(shù)據(jù)表的表操作效率基本相同, 從而使得PLP策略能夠滿足化學(xué)加速計(jì)算的需求.然而, 在包含間斷的瞬態(tài)可壓縮反應(yīng)流問題中, 不同分區(qū)之間的化學(xué)熱力學(xué)狀態(tài)分布可能差異很大, 并且隨時(shí)間會(huì)發(fā)生劇烈變化, 因此, 傳統(tǒng)的PLP策略[10]并不總能保證各分區(qū)所對應(yīng)的數(shù)據(jù)表在化學(xué)反應(yīng)處理(操作效率)上的均衡, 從而導(dǎo)致計(jì)算效率的下降.因此, 本文針對二維爆轟波數(shù)值模擬提出了一種全新的并行策略, 轉(zhuǎn)置建表策略(transposed processing, TP), 如圖1(b)所示.在該策略中, 數(shù)據(jù)表(T1-Tn)按序建立在其對應(yīng)分區(qū)(S1-Sn)轉(zhuǎn)置后的位置上.

(a)Purely local processing strategy

在二維數(shù)值模擬中, 當(dāng)爆轟波傳播方向與分區(qū)邊界相垂直, TP策略可以均勻劃分包含間斷的劇烈反應(yīng)區(qū), 并將它們分配到各數(shù)據(jù)表中以保持表操作的負(fù)載均衡.相應(yīng)地, 當(dāng)爆轟波傳播方向與分區(qū)邊界平行時(shí), PLP策略同樣具有上述計(jì)算優(yōu)勢.因此, PLP策略和TP策略被認(rèn)為是兩個(gè)等價(jià)策略.針對流場中不同的爆轟波構(gòu)型, 選取合適的并行策略可以有效地提升數(shù)據(jù)表操作之間的負(fù)載均衡性.因?yàn)檫@兩種策略決定了并行計(jì)算中各數(shù)據(jù)表的初始結(jié)構(gòu), 本文將PLP策略[10]和TP策略合稱為底層建表策略.

在實(shí)際計(jì)算過程中, 二維爆轟波的傳播方向并不總與分區(qū)邊界正交(如本文所研究的二維斜爆轟計(jì)算), 這意味著, 針對某種爆轟波構(gòu)型, 即使選取了適用的底層建表策略, 仍不能總保證每個(gè)數(shù)據(jù)表在操作效率上的一致.另外, 考慮到瞬態(tài)問題中化學(xué)熱力學(xué)狀態(tài)分布的隨機(jī)性和不穩(wěn)定性, 在某個(gè)具體的時(shí)間步內(nèi), 各數(shù)據(jù)表所執(zhí)行的表操作并不完全一致.當(dāng)隨機(jī)性和不穩(wěn)定性發(fā)生在流場內(nèi)化學(xué)反應(yīng)比較劇烈的區(qū)域時(shí), 這種不一致性會(huì)愈發(fā)明顯.在此基礎(chǔ)上, 本文提出了另一類并行策略, 即, 負(fù)載均衡策略, 以進(jìn)一步改善數(shù)據(jù)表操作上的負(fù)載均衡性.所提出的負(fù)載均衡策略包括數(shù)據(jù)分?jǐn)偛呗?data apportion processing, DAP)和數(shù)據(jù)交換策略(data exchange processing, DEP), 兩種策略的原理見圖2.

(a)Data apportion processing strategy

DAP策略的基本思想可以理解為: 對計(jì)算過程中的每個(gè)時(shí)間步, 通過比較不同數(shù)據(jù)表的表操作時(shí)間, 將操作時(shí)間長(操作效率低)的數(shù)據(jù)表中的一部分計(jì)算任務(wù)通過MPI協(xié)議函數(shù)分?jǐn)偟讲僮鲿r(shí)間短(操作效率高)的數(shù)據(jù)表中, 讓后者幫助前者處理.當(dāng)針對這部分任務(wù)的數(shù)據(jù)表操作完成之后, 再將其返回到原數(shù)據(jù)表中, 以盡可能提高不同數(shù)據(jù)表操作效率之間的一致性.DAP策略的具體執(zhí)行方式如圖2(a)所示.在每個(gè)時(shí)間步(大于等于第2步)計(jì)算開始之前, 首先對所有數(shù)據(jù)表操作耗費(fèi)的累計(jì)CPU時(shí)間(墻上時(shí)間)做統(tǒng)計(jì)和排序.依次地, 將耗時(shí)最長的數(shù)據(jù)表(記為Tmax)中一半的節(jié)點(diǎn)通過分?jǐn)偛僮鱾鬟f給耗時(shí)最短的數(shù)據(jù)表(Tmin), 將耗時(shí)第二長的數(shù)據(jù)表的一半節(jié)點(diǎn)分?jǐn)偨o耗時(shí)第二短的數(shù)據(jù)表, 直到所有數(shù)據(jù)表都完成配對為止.需要注意的是, 當(dāng)并行計(jì)算中分區(qū)/數(shù)據(jù)表的個(gè)數(shù)n為偶數(shù)時(shí), 配對數(shù)量為n/2; 當(dāng)n為奇數(shù)時(shí), 配對數(shù)量為(n-1)/2, 而剩余的單個(gè)數(shù)據(jù)表因其計(jì)算效率恰好居中, 不對其進(jìn)行數(shù)據(jù)分配的操作.以上數(shù)據(jù)傳遞過程結(jié)束之后, 各個(gè)數(shù)據(jù)表所需處理的計(jì)算任務(wù)將不完全相同, 從數(shù)量上來看, 作為接收方的數(shù)據(jù)表需要承擔(dān)的計(jì)算任務(wù)是作為發(fā)送方的數(shù)據(jù)表任務(wù)的3倍; 但是, 從計(jì)算時(shí)間上來看, 配對數(shù)據(jù)表雙方在該時(shí)間步內(nèi)進(jìn)行表操作所耗費(fèi)的CPU時(shí)間將比分?jǐn)傊案咏咏? 這是由各數(shù)據(jù)表在操作效率上的差異導(dǎo)致的.當(dāng)此時(shí)間步內(nèi)所有數(shù)據(jù)表的化學(xué)反應(yīng)計(jì)算完成之后, 被分?jǐn)偟囊话牍?jié)點(diǎn)將返回原數(shù)據(jù)表, 該時(shí)間步內(nèi)的計(jì)算結(jié)束.接下來, 在下一個(gè)時(shí)間步的化學(xué)反應(yīng)計(jì)算開始之前, 再對各數(shù)據(jù)表操作的累計(jì)CPU時(shí)間重新排序, 然后再進(jìn)行分?jǐn)偛僮?可以看出, DAP策略的根本目的在于改善數(shù)據(jù)表操作效率上的負(fù)載均衡性, 而其中每個(gè)表需要處理的計(jì)算任務(wù)多少相互之間并不相同.

相應(yīng)地, DEP策略同樣依賴于MPI通訊協(xié)議對各數(shù)據(jù)表的計(jì)算任務(wù)進(jìn)行重新分配, 但該策略是在數(shù)據(jù)交換操作的基礎(chǔ)上完成的.DEP策略的具體執(zhí)行方式如圖2(b)所示, 在每一計(jì)算時(shí)間步開始之前, 數(shù)據(jù)表根據(jù)操作時(shí)間進(jìn)行排序和配對的方式和前面的DAP策略相同, 但不同的是配對后的兩個(gè)數(shù)據(jù)表的數(shù)據(jù)傳遞將從DAP中的單向傳遞關(guān)系變化為DEP中的雙向傳遞關(guān)系.因此, DEP策略能夠在提升數(shù)據(jù)表操作負(fù)載均衡性的同時(shí), 保證各個(gè)數(shù)據(jù)表所需處理的任務(wù)量一致.需要注意的是, 數(shù)據(jù)表在化學(xué)反應(yīng)操作中耗費(fèi)的CPU時(shí)間主要依賴于熱力學(xué)狀態(tài)的分布而非計(jì)算任務(wù)的多少, 因此, 單獨(dú)比較DAP和DEP策略, 并不能直接得出誰更好的結(jié)論.通常認(rèn)為, 當(dāng)各個(gè)數(shù)據(jù)表需要處理的操作任務(wù)(流場內(nèi)的化學(xué)熱力學(xué)狀態(tài)分布)差異較小時(shí), DEP策略的計(jì)算效率要高于DAP策略, 反之, 則DAP策略的效率更高.

2 二維斜爆轟數(shù)值模擬

2.1 控制方程和計(jì)算格式

針對二維2H2+O2斜爆轟問題的數(shù)值模擬, 本文采用基于理想氣體的帶二維詳細(xì)化學(xué)反應(yīng)機(jī)理的多組分Euler方程求解, 控制方程表達(dá)如下

(12)

其中

計(jì)算時(shí)采用分裂算法將控制方程(12)中的對流項(xiàng)和化學(xué)反應(yīng)源項(xiàng)解耦處理, 在規(guī)定的每個(gè)計(jì)算時(shí)間步內(nèi), 首先在不考慮化學(xué)反應(yīng)的情況下計(jì)算流動(dòng)過程, 然后利用流動(dòng)過程更新后的流場參數(shù)計(jì)算化學(xué)反應(yīng)過程.針對反應(yīng)源項(xiàng)S采用經(jīng)典2H2+O2詳細(xì)機(jī)理[14], 其中包含9種化學(xué)反應(yīng)組分和19個(gè)基元反應(yīng), 分別通過并行建表化學(xué)加速算法和直接積分方法求解.流動(dòng)過程(F和G)采用Lax-Friedrichs通量分裂格式結(jié)合高精度9階WENO格式[15]計(jì)算.最后, 采用3階顯式Runge-Kutta方法處理時(shí)間項(xiàng)U, 計(jì)算的定時(shí)間步長為1.0×10-9s.

2.2 計(jì)算模型

針對二維斜爆轟波流場建立的數(shù)值計(jì)算模型如圖3所示.圖中矩形區(qū)為計(jì)算區(qū)域,x方向上的長度(Lx)為15 mm,y方向上的長度(Ly)為 7.5 mm.初始來流為當(dāng)量比H2∶O2∶N2=2∶1∶4的超聲速可燃預(yù)混氣體, 初始溫度和初始?jí)毫Ψ謩e設(shè)置為T0=298.15 K和p0=1 atm(1 atm=1.013 25×105Pa).圖3中楔形區(qū)域(灰色部分)的半錐角θ代表了來流與楔平面之間的角度, 本文計(jì)算中設(shè)定來流Mach數(shù)和半錐角分別為Ma=7.0和θ=30°.計(jì)算初始時(shí)刻, 斜激波和斜爆轟波在流場中的位置和角度可以根據(jù)Ma和θ由斜激波關(guān)系式[16]求出.在整個(gè)計(jì)算區(qū)域中, 楔形表面對應(yīng)的邊界采用絕熱滑移剛性壁面條件, 其他邊界均采用零梯度邊界條件.為精確捕捉二維斜爆轟流場內(nèi)的波系結(jié)構(gòu), 使用尺寸為Δx=Δy=0.01 mm的均勻正交網(wǎng)格覆蓋圖3所示的整個(gè)計(jì)算區(qū)域范圍.基于該網(wǎng)格尺寸, 流場內(nèi)半反應(yīng)區(qū)的長度包含了約24個(gè)網(wǎng)格, 足以描述斜爆轟波流場的精細(xì)結(jié)構(gòu).

圖3 二維斜爆轟流場示意圖(OSW—斜激波, ODW—斜爆轟波, CW/DW—燃燒/爆燃波)

計(jì)算開始后, 由于可燃預(yù)混氣流的超聲速傳播, 首先會(huì)在圖3所示的楔面上方形成以楔形頂點(diǎn)為起點(diǎn)的OSW, 楔面與斜激波之間的角度由來流的Mach數(shù)決定.此后, 來流中的可燃預(yù)混氣體被斜激波誘導(dǎo)燃燒, 并通過化學(xué)反應(yīng)釋放能量支持斜激波傳播, 進(jìn)而在流場(見圖3)下游形成ODW.從OSW到ODW的發(fā)展過程, 通常伴隨著流場內(nèi)復(fù)雜的波系結(jié)構(gòu)變化, 包括DW和CW.此外, 計(jì)算初期斜爆轟波陣面在下游所處的位置并不駐定.為了避免上述這些非物理過程的影響, 本文選取穩(wěn)定階段的10 000個(gè)時(shí)間步(8 000～18 000)作為斜爆轟波數(shù)值模擬的一個(gè)完整的計(jì)算過程.圖4給出了計(jì)算結(jié)束時(shí)刻(時(shí)間步18 000), 采用直接方法計(jì)算化學(xué)反應(yīng)得到的斜爆轟波流場密度分布圖, 圖中清晰地顯示了包含斜激波、斜爆轟波、爆燃波以及它們匯聚形成的三波點(diǎn)等波系結(jié)構(gòu).定義從楔頂點(diǎn)到三波點(diǎn)的水平距離(x方向)為斜爆轟波的誘導(dǎo)區(qū)長度, 經(jīng)統(tǒng)計(jì), 本文數(shù)值模擬的誘導(dǎo)區(qū)長度為5.0 mm, 這與文獻(xiàn)[16]中計(jì)算得到的誘導(dǎo)區(qū)長度一致, 表明了本文計(jì)算模型和數(shù)值方法的合理性及可靠性.

圖4 直接積分求解的計(jì)算結(jié)束時(shí)刻(時(shí)間步18 000)2H2+O2斜爆轟波流場的密度分布圖

3 結(jié)果與討論

為考察基于不同并行策略的化學(xué)加速算法在二維斜爆轟波數(shù)值模擬中的表現(xiàn), 本文分別采用了基于PLP策略的3種算法與基于TP策略的3種算法, 對化學(xué)反應(yīng)過程進(jìn)行加速計(jì)算, 并與DI求解化學(xué)反應(yīng)過程的結(jié)果對比, 考察并行建表化學(xué)加速算法在計(jì)算精度和計(jì)算效率兩個(gè)方面的表現(xiàn), 見表1.

表1 并行建表化學(xué)加速算法在二維斜爆轟波數(shù)值模擬中的計(jì)算性能

本文所有算例(包括直接積分算例)均采用15個(gè)MPI并行分區(qū), 每個(gè)分區(qū)對應(yīng)一個(gè)數(shù)據(jù)表, 且將數(shù)據(jù)表的容差判據(jù)統(tǒng)一設(shè)定為εtol=1×10-3.根據(jù)計(jì)算機(jī)內(nèi)存空間上限, 將總表尺寸和單表尺寸分別設(shè)置為M5=12.0×107和Msin=1.8×107.

3.1 計(jì)算精度

在二維斜爆轟波流場中, 由斜激波和斜爆轟波連接組成的波陣面與計(jì)算區(qū)域的兩個(gè)方向(x和y方向)存在明顯的傾斜角(見圖4).在此基礎(chǔ)上, PLP策略和TP策略在本文流場中的應(yīng)用效果存在“先天的”相似性, 故基于兩種策略的各算法在計(jì)算精度上的表現(xiàn)也比較相似.因此, 本文只考察基于TP策略的3種并行建表化學(xué)加速算法(TP，TP/DAP和TP/DEP)的計(jì)算精度.圖5首先給出了計(jì)算結(jié)束時(shí)刻(時(shí)間步18 000)3種算法計(jì)算的2H2+O2斜爆轟波二維流場圖, 并將其與對應(yīng)的直接積分結(jié)果作對比.結(jié)果表明, 基于TP策略的3種算法計(jì)算得到的壓力、密度和溫度的分布結(jié)果均與DI的結(jié)果高度一致, 反映了本文提出的并行建表化學(xué)加速算法在斜爆轟數(shù)值模擬中的精確性.

(a)Density

為定量分析帶數(shù)據(jù)表尺寸控制策略的并行建表化學(xué)加速算法在二維斜爆轟波數(shù)值模擬中的計(jì)算精度, 圖6給出了結(jié)束時(shí)刻(時(shí)間步18 000), 基于TP策略的3種算法和DI計(jì)算的緊靠楔形表面的一層網(wǎng)格內(nèi)(y=0 mm), 壓力、溫度以及組分OH和H2O2的質(zhì)量分?jǐn)?shù)沿x方向的分布曲線.結(jié)果表明, 基于TP的3種算法對不同物理量的計(jì)算結(jié)果均與DI結(jié)果保持高度一致, 再次反映了并行建表化學(xué)加速算法的精確性.圖6中的結(jié)果還指出, 各算例(包括DI)對應(yīng)的積分量分布曲線在流場下游位置均存在一定程度的波動(dòng).這是由瞬態(tài)反應(yīng)流問題本身的數(shù)值不穩(wěn)定導(dǎo)致的, 與并行建表化學(xué)加速算法的使用無關(guān).綜上, 通過與DI計(jì)算結(jié)果的對比發(fā)現(xiàn), 并行建表化學(xué)加速算法在二維斜爆轟數(shù)值模擬中的計(jì)算精度可靠, 因而可以進(jìn)一步考察算法在斜爆轟波數(shù)值模擬中的計(jì)算效率.

(a)Pressure

3.2 計(jì)算效率

為考察帶數(shù)據(jù)表尺寸控制策略的并行建表化學(xué)加速算法在二維斜爆轟波數(shù)值模擬中的計(jì)算效率, 圖7首先給出了各算法在整個(gè)計(jì)算過程(時(shí)間步8 000～18 000)中的加速比Rs曲線, 加速比定義為采用直接積分求解化學(xué)反應(yīng)過程的時(shí)間與采用加速算法求解化學(xué)反應(yīng)過程的時(shí)間之比.總體來說, 各算法的加速比曲線都經(jīng)歷了初始階段的持續(xù)增長過程, 且計(jì)算初期的增長率更高.對比兩種由底層建表策略(PLP和TP)直接構(gòu)成的算法可以發(fā)現(xiàn), 純粹TP算法的加速比明顯高于純粹PLP算法, 這表明TP算法通過轉(zhuǎn)置操作帶來的數(shù)據(jù)傳遞有利于提升不同數(shù)據(jù)表操作之間的負(fù)載均衡性.

圖7 2H2+O2斜爆轟波數(shù)值模擬中各并行建表化學(xué)加速算法在整個(gè)計(jì)算過程(時(shí)間步8 000～18 000)中的加速比曲線

此外, 比較純粹算法和復(fù)合算法在加速比上的差異可以發(fā)現(xiàn), 當(dāng)選擇不同的底層建表策略時(shí), 負(fù)載均衡策略(包括DAP和DEP)對計(jì)算效率的作用也不同.在基于PLP策略的算法中, 復(fù)合算法PLP/DAP和PLP/DEP的加速效果均高于純粹PLP算法; 而在基于TP策略的算法中, TP/DAP和TP/DEP兩種復(fù)合算法的加速比反而不如純粹TP算法.計(jì)算結(jié)束時(shí)刻(穩(wěn)定階段), 各算法的化學(xué)反應(yīng)加速比取值見表1, 所有算例中, 純粹TP算法的計(jì)算效率最高(Rs=3.71).

為了分析不同算法在計(jì)算效率上的差異, 圖8選取了加速比最高和最低的兩種算法(純粹TP算法和TP/DEP算法), 繪制了幾個(gè)典型數(shù)據(jù)表(T4,T6,T13和T14)在整個(gè)計(jì)算過程中的建立/刪除成長過程曲線.在選取的4個(gè)典型數(shù)據(jù)表中,T4和T6對應(yīng)的流場區(qū)域內(nèi)均包含有部分爆轟波陣面; T13對應(yīng)的區(qū)域?yàn)榱鲌鰞?nèi)的已燃?xì)怏w; 而T14對應(yīng)的區(qū)域內(nèi)包含了三波點(diǎn)和反應(yīng)誘導(dǎo)區(qū).圖8的結(jié)果表明, 純粹TP算法中不同數(shù)據(jù)表執(zhí)行節(jié)點(diǎn)刪除操作的次數(shù)具有一定的差異; 而TP/DEP算法中4個(gè)典型數(shù)據(jù)表執(zhí)行節(jié)點(diǎn)刪除操作的次數(shù)幾乎相同, 說明疊加DEP策略可以有效改善數(shù)據(jù)表計(jì)算負(fù)載上的不均衡.此外, 通過觀察圖8中不同數(shù)據(jù)表執(zhí)行節(jié)點(diǎn)刪除操作的時(shí)間步可以發(fā)現(xiàn), 引入數(shù)據(jù)表尺寸控制策略, 表操作尤其是節(jié)點(diǎn)刪除操作發(fā)生的時(shí)機(jī)具有很高的同步性, 避免了計(jì)算過程中不同數(shù)據(jù)表操作之間相互等待的時(shí)間耗費(fèi).

(a)Single TP algorithm

表1給出了不同算法中所有15個(gè)數(shù)據(jù)表執(zhí)行節(jié)點(diǎn)刪除次數(shù)的標(biāo)準(zhǔn)差(σN)取值, 用于表征各算法在數(shù)據(jù)表操作負(fù)載均衡性上的差異.結(jié)果顯示TP/DEP算法的取值要低于純粹TP算法, 再次驗(yàn)證了DEP策略對負(fù)載均衡性的提升.然而, 通過表1中數(shù)據(jù)表執(zhí)行節(jié)點(diǎn)刪除操作的平均次數(shù)(表1中aN)可以發(fā)現(xiàn), 純粹TP算法在整個(gè)斜爆轟數(shù)值模擬過程中所有15個(gè)數(shù)據(jù)表節(jié)點(diǎn)刪除操作的平均次數(shù)只有2.2次, 明顯小于TP/DAP和TP/DEP兩種復(fù)合算法.節(jié)點(diǎn)刪除操作的次數(shù)代表了計(jì)算過程中數(shù)據(jù)表內(nèi)插入節(jié)點(diǎn)的頻次, 通常由數(shù)據(jù)表內(nèi)節(jié)點(diǎn)的重復(fù)取回比例決定, 更小的節(jié)點(diǎn)刪除操作次數(shù)意味著表內(nèi)更多的節(jié)點(diǎn)被取回, 也即直接積分插入節(jié)點(diǎn)的數(shù)量更少.表1進(jìn)一步給出了各算法中所有15個(gè)數(shù)據(jù)表內(nèi)節(jié)點(diǎn)取回率(定義為數(shù)據(jù)節(jié)點(diǎn)成功取回次數(shù)與數(shù)據(jù)查詢次數(shù)的比值)的平均值, 可以發(fā)現(xiàn)不同算法在平均取回率上的表現(xiàn)相近, 而純粹TP算法中數(shù)據(jù)節(jié)點(diǎn)的插入次數(shù)明顯更少, 意味著其中數(shù)據(jù)表節(jié)點(diǎn)的取回次數(shù)更高.綜合來看, 純粹TP算法中各數(shù)據(jù)表的計(jì)算負(fù)載比較均衡且同步性較高, 在此基礎(chǔ)上, 通過DEP策略只能一定程度地提升計(jì)算效率.同時(shí), 由于純粹TP算法執(zhí)行節(jié)點(diǎn)刪除操作的平均次數(shù)更少, 意味著需要重新積分插入新節(jié)點(diǎn)的次數(shù)也相對更少, 在數(shù)據(jù)節(jié)點(diǎn)取回率相近的基礎(chǔ)上, 控制節(jié)點(diǎn)插入的數(shù)量可以有效避免大量直接積分求解剛性方程的時(shí)間, 提升化學(xué)反應(yīng)過程的計(jì)算效率.

圖9進(jìn)一步給出了純粹TP算法中, 幾個(gè)典型數(shù)據(jù)表(T4,T6,T13和T14)在狀態(tài)空間中的可達(dá)區(qū)域和可重復(fù)區(qū)域, 以表征2H2+O2斜爆轟波數(shù)值模擬中流場不同位置對應(yīng)的數(shù)據(jù)表內(nèi)的熱力學(xué)狀態(tài)分布.可達(dá)區(qū)域(圖9中灰色散點(diǎn))可以理解為整個(gè)計(jì)算過程中斜爆轟波流場內(nèi)實(shí)際出現(xiàn)的所有化學(xué)熱力學(xué)狀態(tài)的集合; 可重復(fù)區(qū)域(圖9中黑色散點(diǎn))則可以理解為整個(gè)計(jì)算過程中流場內(nèi)重復(fù)出現(xiàn)的化學(xué)熱力學(xué)狀態(tài)的集合.

圖9 純粹TP算法中典型數(shù)據(jù)表(T4, T6, T13和T14)在化學(xué)熱力學(xué)狀態(tài)空間中的散點(diǎn)分布圖

從可達(dá)區(qū)域的范圍來看, 不同數(shù)據(jù)表內(nèi)出現(xiàn)的化學(xué)熱力學(xué)狀態(tài)范圍明顯不同, 包含三波點(diǎn)和反應(yīng)誘導(dǎo)區(qū)在內(nèi)的數(shù)據(jù)表T14的可達(dá)區(qū)域范圍明顯比其他數(shù)據(jù)表更為復(fù)雜.從可重復(fù)區(qū)域?qū)蛇_(dá)區(qū)域的覆蓋情況來看, 數(shù)據(jù)表中能夠重復(fù)取回利用的熱力學(xué)狀態(tài)范圍較大, 尤其在數(shù)據(jù)表T14對應(yīng)的計(jì)算區(qū)域內(nèi), 大量數(shù)據(jù)節(jié)點(diǎn)可以被重復(fù)取回利用(經(jīng)統(tǒng)計(jì)純粹TP算法中數(shù)據(jù)表T14的取回率達(dá)到91.35%).考慮到該區(qū)域內(nèi)的化學(xué)熱力學(xué)狀態(tài)十分復(fù)雜且變化劇烈, 對應(yīng)數(shù)據(jù)表內(nèi)可重復(fù)區(qū)域的覆蓋率高意味著直接積分操作的次數(shù)少, 進(jìn)而解釋了純粹TP算法更高的化學(xué)反應(yīng)加速比.

結(jié)合表1的結(jié)果可以發(fā)現(xiàn), 更大的可重復(fù)區(qū)域范圍使得數(shù)據(jù)表中能夠取回重復(fù)利用的熱力學(xué)狀態(tài)更多, 計(jì)算過程中節(jié)點(diǎn)刪除操作的平均次數(shù)更少.在此基礎(chǔ)上, 每個(gè)數(shù)據(jù)表內(nèi)包含的有效節(jié)點(diǎn)數(shù)量已經(jīng)基本涵蓋了對應(yīng)分區(qū)內(nèi)的熱力學(xué)狀態(tài), 而通過DAP和DEP進(jìn)行數(shù)據(jù)傳遞后, 表內(nèi)被引入大量其他計(jì)算區(qū)域?qū)?yīng)的化學(xué)熱力學(xué)狀態(tài)節(jié)點(diǎn), 導(dǎo)致了數(shù)據(jù)表內(nèi)節(jié)點(diǎn)插入和刪除的次數(shù)增加, 影響了化學(xué)反應(yīng)過程的計(jì)算效率.因此, 純粹TP算法在2H2+O2斜爆轟數(shù)值模擬中的計(jì)算效率更高.

4 結(jié)論

本文以帶節(jié)點(diǎn)刪除機(jī)制的ISAT方法為基礎(chǔ), 針對可壓縮反應(yīng)流的大規(guī)模并行數(shù)值模擬手段, 構(gòu)造了幾種基于數(shù)據(jù)表操作的并行策略, 并由此提出了一系列并行建表化學(xué)加速算法, 將其運(yùn)用到二維2H2+O2斜爆轟數(shù)值模擬中, 考察了其在計(jì)算精度和計(jì)算效率上的表現(xiàn).研究結(jié)果表明, 本文提出的基于不同策略的所有并行建表化學(xué)加速算法, 均能在不損失計(jì)算精度的基礎(chǔ)上有效提升化學(xué)反應(yīng)過程的計(jì)算效率.對比不同算法的化學(xué)反應(yīng)加速比發(fā)現(xiàn), 并行策略的選取對計(jì)算效率具有明顯的影響.在包含15個(gè)分區(qū)的所有算例中, 基于PLP策略的復(fù)合算法(PLP/DAP和PLP/DEP)比純粹PLP算法的計(jì)算效率更高; 而基于TP策略的復(fù)合算法(TP/DAP和TP/DEP)則比純粹TP算法的計(jì)算效率更低.在所有算例中, 純粹TP算法在計(jì)算結(jié)束時(shí)刻的化學(xué)反應(yīng)加速比最高(Rs=3.71).通過分析斜爆轟流場化學(xué)熱力學(xué)狀態(tài)空間的分布規(guī)律發(fā)現(xiàn), 狀態(tài)空間內(nèi)可重復(fù)區(qū)域?qū)蛇_(dá)區(qū)域的覆蓋率越高, 數(shù)據(jù)表內(nèi)節(jié)點(diǎn)被重復(fù)取回利用的次數(shù)就越多, 整個(gè)計(jì)算過程中數(shù)據(jù)表執(zhí)行節(jié)點(diǎn)刪除操作的次數(shù)也就越少, 并最終體現(xiàn)為更高的化學(xué)反應(yīng)加速比.

致謝感謝國家自然科學(xué)基金(11702032)和航天一院CALT聯(lián)合基金(CALT2020-ZD01)對本文的大力支持.