陶蘇林， 李雨鴻， 周 寅， 劉 垚

(1.南京信息工程大學(xué)應(yīng)用氣象學(xué)院,南京 210044；2.南京大橋機(jī)器有限公司技術(shù)開發(fā)中心,南京 211101；3.遼寧省氣象科學(xué)研究所,沈陽 110166；4.94857部隊(duì)61分隊(duì),蕪湖 241007；5.寧夏氣象防災(zāi)減災(zāi)重點(diǎn)實(shí)驗(yàn)室,銀川 750002)

描述實(shí)際大氣問題建立的動力學(xué)模型通常是高維或無限維的非線性系統(tǒng),模型有效性受限于離散格式較多的自由度[1]。尋求既可保證模擬精度又可降低計(jì)算量的模型降階技術(shù)是區(qū)別于并行計(jì)算的有效思路。構(gòu)造降階模型基本思路是以較精確的數(shù)學(xué)描述刻畫原數(shù)值模型的主要動力學(xué)特征,并且這種數(shù)學(xué)描述的階數(shù)和計(jì)算代價遠(yuǎn)低于原數(shù)值模型[2]。

特征正交分解方法(proper orthogonal decomposition method, POD)提供已知大量數(shù)據(jù)在最小二乘意義下的一組正交基,以實(shí)現(xiàn)給定數(shù)據(jù)的最優(yōu)低維逼近,其降階效能在諸多研究中得到證明[3—6]。POD降階模型可較好模擬流場運(yùn)動細(xì)節(jié)、刻畫較強(qiáng)的非線性特征,但對非線性項(xiàng)的計(jì)算代價依然依賴于原始系統(tǒng)的維度m。離散經(jīng)驗(yàn)插值方法(discrete empirical interpolation method, DEIM)則進(jìn)一步給出了非線性項(xiàng)降階處理思路,即在精選的n個(n<

1 淺水波模式

引入Parrett和Cullen測試使用的一維淺水波模式[9]。該模式形式簡單,但其運(yùn)行依然可以達(dá)到極為復(fù)雜的程度,在測試氣象問題方面十分有效,且有典型性。該模式刻畫均質(zhì)不可壓縮無黏旋轉(zhuǎn)淺流體的基本運(yùn)動,忽略層結(jié),描述大氣與海洋運(yùn)動的重要特征:①流體是旋轉(zhuǎn)的,科氏加速度是重要的；②運(yùn)動的形態(tài)比是小量。Lorenc[10]基于該模式討論了有關(guān)非線性系統(tǒng)的資料同化問題。此外,Griffith等[11]也利用該模式構(gòu)建討論變分同化方案在估算系統(tǒng)誤差分量與時間相關(guān)的模式誤差分量研究中的有效性。

該模式包含科氏項(xiàng)和剛性底部地形邊界以描述淺層大氣運(yùn)動,其控制方程為[9]

{?u?t+u?u?x+?φ?x-fv+g?H?x=0

?v?t+u?v?x+fu=0

?φ?t+u?φ?x+φ?u?x=0

式(1)中,x∈Ω=(0,2πL)且t∈(0,T)；u=u(x,t)且v=v(x,t) 分別表示緯向和經(jīng)向風(fēng)速；f為科氏參數(shù)；φ=gη(x,t)為重力場,其中g(shù)為重力加速度,且η(x,t)>0為底部以上大氣流體的平均深度。剛性底邊界由H=H(x)定義。

式(1)雖然形式簡單,但為非線性,并且通過該框架可探究忽略層結(jié)情況下旋轉(zhuǎn)作用和底部地形對大氣運(yùn)動的影響。淺水波模式[式(1)]的離散格式參見文獻(xiàn)[12],該格式通過引入擴(kuò)散項(xiàng)來抑制二階有限差分格式在地形躍變(如山峰)附近產(chǎn)生的虛假震蕩。Parrett和Cullen[9]通過與Glimm法得到的參考解比較,證明使用的離散格式可以提供足夠精確的解。Griffith[13]也根據(jù)Parrett和Cullen提出的包含地形條件的案例測試了該離散格式,并得到了與Parrett和Cullen研究一致的結(jié)果。

若定義tk時刻的風(fēng)場和位勢場：

{uk=(uk1,uk2,…,ukj)T

vk=(vk1,vk2,…,vkj)T

φk=(φk1,φk2,…,φkj)T

并且：

{mk=(uk○φk)

nk=(vk○φk)

于是淺水波模式[式(1)]離散格式可寫成:

(4)

式(4)中,f1(uk-1,mk-1)、f2(φk-1)和f3(vk-1,mk-1)為非線性項(xiàng),具體形式為

(5)

其他矩陣形式如下:

(6)

(7)

(8)

(9)

(10)

(11)

(12)

取科氏參數(shù)f=7.292×10-5s-1,擴(kuò)散系數(shù)K=5.0×102m2·s-1,重力加速度g=9.8 m·s-2,空間尺度坐標(biāo)L=9.55×105m,時間尺度坐標(biāo)T=1.2×105s。利用這些參數(shù)設(shè)置使得淺水波模式[式(4)]可以模擬30°N的地球大氣運(yùn)動[11]。另外,模式剛性底邊界定義為

(13)

式(13)中,a=10Δx。

令空間格點(diǎn)數(shù)J=100,積分時間步數(shù)N=100；在[0,2πL]×[0,T]范圍內(nèi),淺水波模式[式(4)]的初始值條件定義為

(14)

式(14)中，大氣流體平均深度為ηm=1 m,并且j=1,2,…,J。因此,數(shù)值試驗(yàn)中大氣流體運(yùn)動是由中心位置的地形開始發(fā)展的。此外,為保證離散格式[式(4)]的穩(wěn)定性,時間步長Δt應(yīng)遵循CFL(courant-friedrichs-lewy)條件。

2 POD/DEIM降階解法

給定Hilbert空間L2內(nèi)一個自由度為m的非線性系統(tǒng):

(15)

式(15)中，空間x∈Ω,時間t∈[0,T],解y(x,t)=[y(x1,t),y(x2,t),…,y(xm,t)]T∈Rm×1,系數(shù)矩陣M,K∈Rm×m以及關(guān)于y(x,t)的非線性函數(shù)F={F[y(x1,t)],F[y(x2,t)],…,F[y(xm,t)]}T(l ?R時F:l→R)。定義該模型的緊湊形式為

(16)

則構(gòu)建該模型的POD/DEIM降階模型步驟如下。

(1)根據(jù)一定離散方法構(gòu)建原始數(shù)學(xué)物理模型的數(shù)值格式,將原始模型的連續(xù)形式轉(zhuǎn)變?yōu)橛善⒎址匠虡?gòu)成的離散形式,即確定模式。

(2)積分確定模式獲得約定時段的模型數(shù)值解,按照約定取樣方式選取瞬像(snapshot)[14]構(gòu)成瞬像集合Y。對瞬像集合Y進(jìn)行奇異值分解,獲得降序排列的特征值及其對應(yīng)的特征向量。這些特征向量構(gòu)成正交基函數(shù),進(jìn)一步用于構(gòu)建POD基函數(shù)。

(3)選取前r個正交基函數(shù)構(gòu)成一個子空間,重構(gòu)模式變量yPOD。

(4)將重構(gòu)解代入第1步的確定模式,并根據(jù)Galerkin投影變換進(jìn)一步將確定模式投影至由POD基函數(shù)張成的降維空間,獲得原始模式的降階形式

(17)

(5)進(jìn)一步利用DEIM方法處理POD降階模型[式(17)]中的非線性項(xiàng)。

(6)求解常微分方程,得到降階系數(shù)。返回第3步,利用r個正交基函數(shù)重構(gòu)模式變量yPOD。

2.1 降階基向量選取

2.1.1 獲得瞬像樣本

進(jìn)而基于數(shù)據(jù)集①～⑤構(gòu)建新的樣本數(shù)據(jù)集:{mk=0,m1,…,mN}, {nk=0,n1,…,nN}, {φk=0,φ1,…,φN}。其中:k=0,1,6,11,…,N-4時刻的狀態(tài)變量取自數(shù)據(jù)集①相應(yīng)時刻的狀態(tài)變量；k=2,7,12,…,N-3時刻的狀態(tài)變量取自數(shù)據(jù)集②；k=3,8,13,…,N-2時刻的狀態(tài)變量取自數(shù)據(jù)集③；k=4,9,14,…,N-1時刻的狀態(tài)變量取自數(shù)據(jù)集④；k=5,10,15,…,N時刻的狀態(tài)變量取自數(shù)據(jù)集⑤。

(18)

(19)

2.1.2 特征值分解

對矩陣A、B和C分別進(jìn)行奇異值分解,可以獲得按降序排列的特征值：

(20)

并定義：

(21)

2.1.3 確定POD基向量

選擇r滿足：r=argmin{I(r):I(r)≥γ0}且γ0∈[0,1]則可獲得與這r個特征值對應(yīng)的一組正交基向量,并構(gòu)成降階空間:

(22)

2.2 淺水波模式的POD降階模型

利用式(22)提供的POD基函數(shù)近似狀態(tài)變量mk、nk和φk:

(23)

式(23)中,αm,k、αn,k和αφ,k∈Rr×1為待求解系數(shù)。因此,只要獲得所有時刻的αm,k、αn,k和αφ,k,就可以得到所有時間層的狀態(tài)變量。

將式(23)代入淺水波模式[式(4)]，可得：

(24)

利用(Φm)T、(Φn)T和(Φφ)T分別左乘式(24),則可得到淺水波模式的POD降階模型:

(25)

該降階模型的初始條件為

(26)

周期邊界條件通過矩陣Mi(i=1,2,…,5)設(shè)定,具體見文獻(xiàn)[12]。因此,利用式(25)求解淺水波模式的POD降階模型,獲得系數(shù)αm,k、αn,k和αφ,k,從而重構(gòu)狀態(tài)變量。

2.3 淺水波模式的POD/DEIM降階模式

(27)

式(27)中，Pf1、Pf2和Pf3根據(jù)DEIM算法確定。并且可以得到：

(28)

因此,淺水波模式的POD/DEIM降階形式為

{αm,k=(Φm)TM1MPOD1∈Rr×J(+Φmαm,k-2)+

(Φm)TM2MPOD2∈Rr×J(+Φnαn,k-1)+

(Φm)TM3MPOD3∈Rr×J(+Φφαφ,k-1)+

αn,k=(Φn)TM2POD2∈Rr×J(+Φnαn,k-2)+

αφ,k=(Φφ)TM2POD2∈Rr×J(+Φφαφ,k-2)+

(Φφ)TM5MPOD5∈Rr×J(+Φmαm,k-1)-(Φφ)T

式(26)為該降階模型[式(29)]的初始條件,周期邊界條件通過矩陣Mi(i=1,2,…,5)設(shè)定。并且,由式(29)可知,POD/DEIM降階模型非線性項(xiàng)的計(jì)算僅利用s(s≤J)個空間格點(diǎn)。因此,POD/DEIM降階模型的計(jì)算代價較POD降階模型有了進(jìn)一步減小。

2.4 降階模型精度驗(yàn)證

(30)

3 算例與分析

令FULL-MODEL表示全階模型、POD-ROM表示POD降階模型、POD/DEIM-ROM表示POD/DEIM降階模型。算例結(jié)果基于CPU型號為Intel酷睿i3—2100(主頻3.10 GHz)、內(nèi)存2 GB以及安裝Windows XP Service Pack 3的計(jì)算機(jī)獲得。

3.1 降階模型計(jì)算效率

設(shè)計(jì)8種數(shù)值試驗(yàn)討論不同空間格點(diǎn)數(shù)時降階模型計(jì)算效率。設(shè)定降階模型需捕獲全階模型狀態(tài)變量99.8%以上的能量,即對于狀態(tài)變量m、n和φ均需滿足I(r)≥99.8%；設(shè)定[0,2πL]區(qū)間內(nèi)格點(diǎn)數(shù)分別為150、200、250、300、350、400、450和500,積分步數(shù)為格點(diǎn)數(shù)2倍,瞬像樣本數(shù)為積分步數(shù)1/2。

隨著格點(diǎn)數(shù)增加,FULL-MODEL、POD-ROM 及POD/DEIM-ROM這三種模型CPU耗時均呈增加趨勢,且FULL-MODEL的CPU耗時增幅最大(表1)?？臻g格點(diǎn)數(shù)相同時,三種模型耗費(fèi)CPU時間排序?yàn)?FULL-MODEL > POD-ROM > POD/DEIM-ROM耗時。由此可見,經(jīng)過非線性項(xiàng)處理的POD/DEIM-ROM模型計(jì)算效率最高。同時,隨著格點(diǎn)數(shù)及積分步數(shù)的增加,POD-ROM以及POD/DEIM-ROM所需POD基函數(shù)及非線性項(xiàng)基函數(shù)均呈增加趨勢。并且,為了捕獲全階模型超過99.8%的能量,狀態(tài)變量φ較m和n需選擇更多基函數(shù),非線性項(xiàng)f1比f2和f3也需使用更多基函數(shù)。

利用DEIM方法處理非線性項(xiàng)后,降階模型可捕獲的能量有所下降。POD/DEIM-ROM數(shù)值解(狀態(tài)變量u、v和φ)與FULL-MODEL數(shù)值解的均方根誤差均大于POD-ROM與FULL-MODEL數(shù)值解的均方根誤差(表2)。但隨格點(diǎn)數(shù)增加,兩種降階模型數(shù)值解誤差未明顯增加。POD-ROM模擬的各個狀態(tài)變量(u、v和φ)與FULL-MODEL數(shù)值解的均方根誤差變化在空間格點(diǎn)數(shù)增至300后趨于平緩。空間格點(diǎn)數(shù)小于300時,狀態(tài)變量u和φ的均方根誤差先減后增,狀態(tài)變量v的均方根誤差則呈逐漸增加趨勢。POD/DEIM-ROM模擬的狀態(tài)變量u和φ與FULL-MODEL的均方根誤差在格點(diǎn)數(shù)少于400時逐漸增加,之后有所減少；而狀態(tài)變量v與FULL-MODEL的均方根誤差在格點(diǎn)數(shù)為150的數(shù)值試驗(yàn)中達(dá)最小,其他趨勢與u和φ基本一致。

POD-ROM數(shù)值解(u、v和φ)與FULL-MODEL的相關(guān)系數(shù)總體達(dá)0.995以上,其中u和v相關(guān)系數(shù)達(dá)0.999(表3)。POD/DEIM-ROM數(shù)值解相關(guān)系數(shù)均低于POD-ROM與FULL-MODEL數(shù)值解相關(guān)系數(shù),表明非線性項(xiàng)的進(jìn)一步近似處理將丟失更多的模式能量。因此,使用DEIM方法進(jìn)一步處理模式〗非線性項(xiàng)后,需增加非線性項(xiàng)基函數(shù)維數(shù)以提高POD/DEIM-ROM捕獲全階模型的能量。

表1 格點(diǎn)數(shù)不同時POD-ROM和POD/DEIM-ROM所需CPU時間和POD基函數(shù)

表2 格點(diǎn)數(shù)不同時POD-ROM和POD/DEIM-ROM與FULL-MODEL的均方根誤差

表3 格點(diǎn)數(shù)不同時POD-ROM、POD/DEIM-ROM與FULL-MODEL的相關(guān)系數(shù)

此外,從POD/DEIM-ROM模擬的各個狀態(tài)變量與FULL-MODEL的相關(guān)系數(shù)來看,水平風(fēng)場v分量數(shù)值解效果最好,其次是u分量,位勢場φ模擬效果最差。隨格點(diǎn)數(shù)增加,POD-ROM和POD/DEIM-ROM數(shù)值解(u、v)與FULL-MODEL數(shù)值解相關(guān)系數(shù)無明顯變化,但POD/DEIM-ROM數(shù)值解φ相關(guān)系數(shù)明顯下降,且當(dāng)格點(diǎn)數(shù)增至500時,相關(guān)系數(shù)僅為0.508。

3.2 瞬像選擇與影響

瞬像選擇是降階模型應(yīng)用的重要步驟。瞬像樣本維數(shù)是降階模型的重要參數(shù),特別是使用DEIM方法處理非線性項(xiàng)后,還需考慮非線性項(xiàng)的瞬像選擇,相當(dāng)于增加了降階模型參數(shù)。討論降階模型對瞬像樣本選擇的敏感性,對于指導(dǎo)降階模型的應(yīng)用有很重要的意義。

根據(jù)狀態(tài)變量和非線性項(xiàng)瞬像維數(shù)的不同組合,設(shè)計(jì)了36種敏感試驗(yàn)。設(shè)定狀態(tài)變量m、n和φ均滿足I(r)≥99.8%；設(shè)定[0,2πL]區(qū)間內(nèi)格點(diǎn)數(shù)為100,積分步數(shù)為格點(diǎn)數(shù)的3倍；狀態(tài)變量m、n和φ以及非線性項(xiàng)f1、f2和f3的瞬像樣本數(shù)分別取10、20、50、100、200和300,并設(shè)定狀態(tài)變量m、n和φ的瞬像維數(shù)相等,非線性項(xiàng)f1、f2和f3的瞬像樣本數(shù)相等。

圖1 狀態(tài)變量和非線性項(xiàng)瞬像維數(shù)不同時POD-ROM和POD/DEIM-ROM與FULL-MODEL數(shù)值解的均方根誤差Fig.1 The RMSE of the POD-ROM and POD/DEIM-ROM in different snapshot dimension of the state variable and nonlinear term with FULL-MODEL

圖2 狀態(tài)變量和非線性項(xiàng)瞬像維數(shù)不同時POD-ROM和POD/DEIM-ROM與FULL-MODEL數(shù)值解的相關(guān)系數(shù)Fig.2 The correlation coefficient of the POD-ROM and POD/DEIM-ROM in different snapshot dimension of the state variable and nonlinear term with FULL-MODEL

僅從CPU耗時來看,狀態(tài)變量m、n和φ與非線性項(xiàng)f1、f2和f3的瞬像維數(shù)不同組合對POD-ROM和POD/DEIM-ROM計(jì)算效率沒有明顯影響。但從降階模型模擬的狀態(tài)變量與全階模型狀態(tài)變量的均方根誤差以及相關(guān)系數(shù)來看,不同組合的瞬像維數(shù)對降階模型數(shù)值解影響很明顯(圖1、圖2)。POD/DEIM-ROM數(shù)值解與FULL-MODEL狀態(tài)變量的均方根誤差較大,其對應(yīng)的相關(guān)系數(shù)也較低,POD-ROM結(jié)果卻與之相反。當(dāng)狀態(tài)變量瞬像維數(shù)為200時,POD-ROM和POD/DEIM-ROM數(shù)值解與FULL-MODEL數(shù)值解均方根誤差均較大。并且,POD-ROM數(shù)值解均方根誤差較其他類型的瞬像維數(shù)組合明顯增大,對應(yīng)的相關(guān)系數(shù)也明顯減小,POD/DEIM-ROM對瞬像選擇也有類似敏感性。因此,本試驗(yàn)當(dāng)區(qū)間[0,2πL]內(nèi)空間格點(diǎn)數(shù)為100、積分步數(shù)為300時,POD-ROM和POD/DEIM-ROM狀態(tài)變量瞬像維數(shù)應(yīng)避免設(shè)置為200。

3.3 DEIM插值點(diǎn)選擇與影響

利用POD/DEIM降階理論處理非線性項(xiàng)的關(guān)鍵是確定合理的DEIM插值點(diǎn)。為了討論不同數(shù)量DEIM插值點(diǎn)對POD/DEIM-ROM數(shù)值解的計(jì)算精度影響,對數(shù)值試驗(yàn)作如下設(shè)置:①設(shè)定[0,2πL]區(qū)間內(nèi)格點(diǎn)數(shù)為100,積分步數(shù)為300；②設(shè)定狀態(tài)變量m、n和φ以及非線性項(xiàng)f1、f2和f3的瞬像樣本數(shù)為100；③根據(jù)狀態(tài)變量(m、n和φ)瞬像矩陣POD譜(圖3)確定選擇狀態(tài)變量前10個POD基函數(shù),即令r=10；④設(shè)定非線性項(xiàng)f1、f2和f3的DEIM插值點(diǎn)分別為90、80、70、60、50、40、30、20和10,以測試POD/DEIM-ROM模擬效率。

隨著DEIM插值點(diǎn)維數(shù)減少,POD/DEIM-ROM降階模型耗費(fèi)的CPU時間呈明顯減少趨勢,但該降階模型數(shù)值解u、v和φ與FULL-MODEL數(shù)值解均方根誤差和相關(guān)系數(shù)無明顯變化(表4)。因此,在數(shù)值試驗(yàn)條件下,淺水波模式的POD/DEIM-ROM數(shù)值解受DEIM插值點(diǎn)維數(shù)變化影響不明顯,但計(jì)算效率得到了很大提高。此外,從各個狀態(tài)變量模擬結(jié)果來看,選擇相同數(shù)量DEIM插值點(diǎn)時,POD/DEIM-ROM對位勢場φ的模擬結(jié)果最差,均方根誤差最大,相關(guān)系數(shù)最小,水平風(fēng)場u和v分量模擬結(jié)果較好。

圖3 狀態(tài)變量(m、n和φ)與非線性項(xiàng)(f1、f2和f3)瞬像矩陣POD譜Fig.3 The POD spectrum corresponding to the solution snapshots for the state variables (m, n and φ) and nonlinear terms (f1, f2 and f3)

表4 DEIM插值點(diǎn)維數(shù)不同時POD/DEIM-ROM所需CPU時間及其與FULL-MODEL的均方根誤差和相關(guān)系數(shù)

全階模型FULL-MODEL預(yù)報(bào)的狀態(tài)變量u、v和φ如圖4。狀態(tài)變量瞬像維數(shù)為10,非線性項(xiàng)DEIM插值點(diǎn)維數(shù)為10時的POD/DEIM-ROM數(shù)值解u、v和φ則如圖5。

狀態(tài)變量瞬像維數(shù)為10,非線性項(xiàng)DEIM插值點(diǎn)維數(shù)為10時的POD/DEIM-ROM數(shù)值解u、v和φ與FULL-DODEL數(shù)值解均方根誤差如圖6所示。

POD/DEIM-ROM在前向積分中期和后期模擬結(jié)果最大值逐漸偏離全階模型結(jié)果。并且,隨著POD/DEIM-ROM積分時刻增加,狀態(tài)變量u、v和φ均方根誤差均呈明顯上升趨勢。而非線性項(xiàng)DEIM插值點(diǎn)維數(shù)在10～90間變化時,POD/DEIM-ROM模擬結(jié)果未發(fā)生明顯變化,表明不同DEIM插值點(diǎn)維數(shù)對POD/DEIM-ROM數(shù)值解結(jié)果影響相對較小。

圖4 全階模型預(yù)報(bào)的狀態(tài)變量Fig.4 The forecasted state variable by full model

圖5 降階模型預(yù)報(bào)的狀態(tài)變量Fig.5 The forecasted state variable by reduced-order model

圖6 降階模型與全階模型預(yù)報(bào)的狀態(tài)變量均方根誤差Fig.6 RMSE of the state variable between the reduced-order model and full model

4 討論

隨著對大氣運(yùn)動過程的深入理解,現(xiàn)有數(shù)值天氣預(yù)報(bào)模式不斷得到更新完善。比如,在原先的中尺度模式中嵌套小尺度模式,以提高模擬分辨率(特別是空間分辨率)。但這些復(fù)雜的數(shù)值模式需要消耗大量計(jì)算資源。應(yīng)用并行計(jì)算技術(shù)是解決這個問題的重要手段,如朱小謙[15]研究了高分辨率譜模式的可擴(kuò)展并行算法,提高了模式運(yùn)行效率。但如果可以直接對原始數(shù)值模式進(jìn)行降階處理,則可簡化很多額外算法研究。

采取直接對原始數(shù)值模式進(jìn)行降階處理的思路,對旋轉(zhuǎn)大氣中有限區(qū)域淺水波模式進(jìn)行降階處理。通過POD方法獲得淺水波模式的POD-ROM,并利用DEIM方法進(jìn)一步對非線性項(xiàng)進(jìn)行處理,獲得淺水波模式的POD/DEIM-ROM,試驗(yàn)結(jié)果表明POD/DEIM降階算法是有效的。取狀態(tài)變量瞬像矩陣前10個POD基函數(shù),POD-ROM就可以捕獲全階模式99.8%以上的能量,并且CPU耗時明顯減少,該優(yōu)勢在空間格點(diǎn)變得稠密時尤其明顯。這些結(jié)果與Cao等[16]、Du等[17]和Dimitriu等[18]等研究結(jié)論相同。而通過DEIM方法對淺水波模式的3個非線性項(xiàng)進(jìn)行近似處理后,POD/DEIM-ROM的CPU耗時進(jìn)一步減少。

但是,降階模型數(shù)值解與全階模式數(shù)值解仍然存在一定偏差。特別是經(jīng)過DEIM處理后的降階模型因?yàn)閬G失了全階模式的部分非線性(但可能是非常重要的)信息,其數(shù)值解與全階模式數(shù)值解均方根誤差明顯大于POD-ROM。因此,為了實(shí)現(xiàn)模型降階,同時又保持較高的計(jì)算精度,需考慮對降階模型進(jìn)行訂正。目前,有相關(guān)研究取得了初步成果。如,Du等[17]通過在瞬像矩陣創(chuàng)建過程應(yīng)用SobolevH1范數(shù)來訂正POD降階模型；Wang等[8]通過利用正則化技術(shù)對二維Burgers方程降階模型進(jìn)行訂正,提高了POD/DEIM降階模型精度。本文還未涉及降階模型訂正,將在未來研究中展開討論。

此外,還需注意POD/DEIM方法雖然進(jìn)一步提高了降階模型的計(jì)算效率,同時也引入了更多的可變參數(shù)。這些參數(shù)的設(shè)置對于降階模型的模擬精度影響不容忽視。3.2與3.3節(jié)試驗(yàn)表明淺水波模式的POD-ROM和POD/DEIM-ROM對瞬像和DEIM插值點(diǎn)維數(shù)等可變參數(shù)的敏感性。因此,如何選擇可變參數(shù)的最優(yōu)取值將是POD/DEIM降階模型應(yīng)用中的重要研究問題。

5 結(jié)論

利用特征正交分解方法(POD)和離散經(jīng)驗(yàn)插值方法(DEIM)建立了旋轉(zhuǎn)大氣中有限區(qū)域淺水波模式的降階模型,并研究了降階模型(ROM)的計(jì)算效率以及瞬像和DEIM插值點(diǎn)的選擇與影響。研究得出如下結(jié)論。

(1) POD/DEIM-ROM計(jì)算效率明顯高于POD-ROM和FULL-MODEL,并且可以捕獲FULL-MODEL超過99.8%的能量,特別當(dāng)空間格點(diǎn)數(shù)量明顯增加時,POD/DEIM-ROM計(jì)算效率最高。

(2) 狀態(tài)變量和非線性項(xiàng)瞬像樣本數(shù)變化對降階模型CPU耗時影響不大,但對數(shù)值解質(zhì)量影響比較明顯。POD/DEIM-ROM數(shù)值解質(zhì)量受DEIM插值點(diǎn)維數(shù)變化影響不明顯,但隨DEIM插值點(diǎn)數(shù)量減小,POD/DEIM-ROM的CPU耗時明顯減少。

實(shí)際應(yīng)用中,數(shù)值模式的POD/DEIM-ROM對瞬像和DEIM插值點(diǎn)維數(shù)有較強(qiáng)的敏感性。這些可變參數(shù)的最優(yōu)取值問題將是未來研究難點(diǎn)。

致謝:感謝南京信息工程大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院王曰朋教授和應(yīng)用氣象學(xué)院申雙和教授對本研究的支持和指導(dǎo)。