賈志強

(中國石油大學(xué)勝利學(xué)院 網(wǎng)絡(luò)信息中心,山東 東營 257000)

矢量場可視化是科學(xué)可視化的重要分支,在油氣田油層分析、氣象預(yù)報、海洋數(shù)值模擬等領(lǐng)域都有著廣泛的應(yīng)用[1]。流線可視化是目前矢量場可視化中應(yīng)用最廣泛的方式之一,是業(yè)務(wù)人員分析矢量場最重要的工具。此外,流線還是很多其他可視化方式(如LIC、FTLE等)的基礎(chǔ),受到了可視化領(lǐng)域廣泛的關(guān)注。隨著數(shù)據(jù)規(guī)模的不斷增大,流線可視化的效率逐漸成為大規(guī)模矢量場數(shù)據(jù)快速可視化處理與分析的瓶頸。為此,許多學(xué)者致力于高效流線可視化方法的研究。目前通常的方法是采用并行加速技術(shù),利用計算機集群或者超級計算機,對流場數(shù)據(jù)和可視化任務(wù)進行劃分后分派給各個節(jié)點,多個節(jié)點同時進行計算。流線并行主要有數(shù)據(jù)塊并行和種子點并行兩種方式。數(shù)據(jù)塊并行的典型方法有:Nouanesengsy等[2]提出的負(fù)載敏感數(shù)據(jù)劃分策略、Sujudi等[3]提出的空間分解方法等;種子點并行的典型方法有:Chen等[4]提出的最優(yōu)化動態(tài)數(shù)據(jù)劃分策略、Chen等[5]提出的空間一致性方法等。此外,彭寶云等[6]以及Zhang等[7]利用數(shù)據(jù)預(yù)取技術(shù)進行流線高效可視計算的方法。以上方法雖然從各種角度對流線可視化過程進行加速,但是在流線生成過程中都是利用數(shù)值積分進行計算。針對該問題,筆者首先利用代數(shù)運算將流線表示成矩陣形式,進而將流線在網(wǎng)格單元內(nèi)的積分運算轉(zhuǎn)化成矩陣表示的方程與網(wǎng)格單元各表面求交點的問題,并利用牛頓迭代法進行求解,從而避免繁瑣的數(shù)值積分運算,有效提高大規(guī)模矢量場流線可視化的效率。

1 非結(jié)構(gòu)網(wǎng)格的流線矩陣表示

對于矢量場中的流線,可以看成是如下微分方程的解:

對于大規(guī)?？茖W(xué)與工程計算中最常用的線性矢量場,可以將其表示成如下的形式:

v(x)=Ax+b.

非結(jié)構(gòu)網(wǎng)格由于其表達的靈活性,目前得到越來越廣泛的應(yīng)用。四面體網(wǎng)格是應(yīng)用最廣泛的非結(jié)構(gòu)網(wǎng)格,且任意的結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格都可剖分為四面體網(wǎng)格,故本文中處理四面體網(wǎng)格線性流場。

由常微分方程的理論可知,四面體網(wǎng)格線性流場條件下,微分方程(1)的解可表示為

x(t)=e(t-t0)Ax0+(e(t-t0)A-I)A-1b.

(2)

為了得到流線方程(2)的表達式,將矩陣A進行特征值分解。利用線性代數(shù)知識可以證明,當(dāng)A有3個不同的非零實特征值,或者兩個復(fù)數(shù)特征值、一個非零實特征值時,流線方程(2)都可以表示成如下的形式:

x(t)=SP(t)S-1x0+SQ(t)S-1b.

(3)

式中,矩陣S是A的3個線性無關(guān)的特征向量組成的矩陣(由于A有3個不同的特征值,故必有3個線性無關(guān)的特征向量);當(dāng)A有3個不同的非零實數(shù)特征值時(分別用λ1、λ2和λ3表示),矩陣

當(dāng)A有兩個復(fù)數(shù)特征值、一個非零實特征值時,利用線性代數(shù)知識可推導(dǎo)出類似的P(t)和Q(t)的表達式。

需要說明的是,雖然推導(dǎo)過程看似繁瑣,但由于上述矩陣都是三階矩陣,其計算過程并不耗時,且可在預(yù)處理階段計算,得到的各矩陣當(dāng)作每個網(wǎng)格單元的屬性進行存儲。后續(xù)進行流線可視化處理與分析時,各矩陣無需再次計算。

2 非結(jié)構(gòu)網(wǎng)格的流線可視化方法

2.1 網(wǎng)格單元分類

根據(jù)各網(wǎng)格單元矩陣的特征值,將非結(jié)構(gòu)四面體網(wǎng)格單元分為3類:

(1)平行單元。4個網(wǎng)格頂點的矢量值完全相同;(2)普通單元。網(wǎng)格單元矩陣具有3個非零的實特征值,或有兩個復(fù)數(shù)特征值和一個非零的實特征值;(3)其他單元。除以上兩種情況外的其他網(wǎng)格單元。

2.2 網(wǎng)格單元處理形式

對應(yīng)上述3種情況,分別進行如下處理,得到網(wǎng)格單元內(nèi)部的流線段:

(1)對于平行單元,由于各網(wǎng)格頂點矢量值完全相同,單元內(nèi)部矢量場無變化,此時流線在網(wǎng)格單元內(nèi)部流向與各網(wǎng)格頂點矢量相同,直接計算流向所在射線與網(wǎng)格單元各表面的交點即可,得到的最小正數(shù)t所代表的交點就是流線流經(jīng)該網(wǎng)格單元的出點,連接入點和出點即可得到該網(wǎng)格單元內(nèi)部的流線段;(2)對于普通單元,通過單元矩陣特征值分解,可得到形如式(3)的矩陣表達式,然后利用牛頓迭代法求解式(3)與網(wǎng)格單元各表面的交點,得到的最小正數(shù)t所代表的交點就是流線流經(jīng)該網(wǎng)格單元的出點,連接入點和出點即可得到該網(wǎng)格單元內(nèi)部的流線段;(3)對于其他單元,按照四階隆格-庫塔進行數(shù)值積分,計算網(wǎng)格單元內(nèi)部的流線段。雖然這種情況需要較為耗時的數(shù)值積分,但有統(tǒng)計數(shù)據(jù)表明,對于典型的矢量場數(shù)據(jù),此類網(wǎng)格單元數(shù)目只占總網(wǎng)格單元數(shù)目的1%左右,并不影響流線可視化的整體效率。

2.3 普通單元流線出點的處理

對于網(wǎng)格單元的任一表面,分別用v1、v2和v3表示該表面的3個頂點,則流線與表面的交點計算等價于求解下述方程:

x(t)·((v1-v3)·(v2-v3))=v3·((v1-v3)×(v2-v3)).

(4)

式中,×表示向量的叉積、·表示向量的點積。

聯(lián)合式(3)與式(4),利用牛頓迭代法可求出t的值。求解時,初值設(shè)為0即可。經(jīng)過測試發(fā)現(xiàn),對于絕大部分的求解,牛頓迭代法經(jīng)過4次便可收斂,因而可以顯著節(jié)省數(shù)值積分的時間,提高流線可視化的效率。對于四面體網(wǎng)格單元的每個表面,都可利用牛頓迭代法求出一個t值;4個t值中最小的正數(shù)t代表的交點就是流線流經(jīng)該網(wǎng)格單元的出點。連接入點與出點,就可得到該網(wǎng)格單元內(nèi)部的流線段。

整體算法描述如下:給定流線種子點,首先判斷其在哪個網(wǎng)格單元內(nèi)部,然后利用上述方法求解該網(wǎng)格單元的流線出點;此出點也是與該網(wǎng)格單元相鄰的網(wǎng)格單元的流線入點,然后再利用上述方法求解出點。此過程一直進行下去,直到流線到達矢量場邊界/臨界點或形成封閉軌道。該算法可用偽碼描述如下:

輸入:流線種子點p

輸出:由p生成的流線

數(shù)組E,記錄形成流線的一系列點,初始化為只有一個元素p

找到p所在的網(wǎng)格單元cell

while(cell不為空)

switch(cell.type)//根據(jù)網(wǎng)格單元的不同類型分別處理

caseparallel//平行單元

v=cell.vectorAt(p)//得到p點處的矢量

t=INFINITY//將t設(shè)為無窮大

for(i=1to4)//對cell的4個表面依次處理

tmp=cell.plane[i].distanceTo(p)//得到p到各平面的距離

if(tmp

t=tmp

id=i

E.add(p+t*v)//得到出點,并加入到數(shù)組E中

break

casenormal//普通單元

t=INFINITY//將t設(shè)為無窮大

for(i=1to4)//依次求流線與4個表面的交點

tmp=Newton(p,cell)//利用牛頓迭代法求交點對應(yīng)的t

if(tmp>0 &&tmp

t=tmp

id=i

E.add(getExitPoint(t))//得到出點,并加入到數(shù)組E中

break

caseother//上述兩種情況之外的其他情況

p=RungeKutta-4(p,cell)//利用4階隆格-庫塔求下一積分點

while(cell.contains(p))//如果p還在cell內(nèi),則一直積分下去

p=RungeKutta-4(p,cell)

E.add(p)//將p加入到數(shù)組E中

id=findFaceID(p)//根據(jù)p找到其對應(yīng)的網(wǎng)格單元表面

break

cell=cell.getNeighbourCell(id)//根據(jù)id找到相鄰的下一個網(wǎng)格單元

依次連接E中的各點,得到流線并輸出

當(dāng)網(wǎng)格分辨率較低時(即單個網(wǎng)格單元尺寸較大),上述折線段連線形成的流線可視化效果可能不能令人滿意。此時,可以利用4階貝塞爾曲線表示網(wǎng)格單元內(nèi)部的流線段。若用P、Q分別表示流線在此網(wǎng)格單元的流入點和流出點,用vP和vQ分別表示P、Q處的矢量值,用tin和tout分別表示流線進入和流出網(wǎng)格單元的時間,則貝塞爾曲線的4個控制點分別為:B0=P、B1=P+vP(tout-tin)/3、B2=Q-vQ(tout-tin)/3、B3=Q。

隨著高性能計算機的不斷發(fā)展,數(shù)值模擬的精度越來越細(xì),網(wǎng)格分辨率越來越高,甚至已遠(yuǎn)高于顯示屏幕的分辨率(即1個網(wǎng)格單元顯示在屏幕上不到1個像素)。因此,在通常情況下,直接連接各網(wǎng)格單元的流線入點和出點形成的流線已經(jīng)可以滿足矢量場可視化的精度需要。

3 算法驗證

在i7-3.6GHz處理器、8GB內(nèi)存、Window7 64位操作系統(tǒng)上用C++實現(xiàn)了本文算法。為了驗證本文方法的有效性,分別用二維矢量場數(shù)據(jù)和三維矢量場數(shù)據(jù)進行測試。第一個數(shù)據(jù)為油氣田油層數(shù)據(jù),該數(shù)據(jù)為二維數(shù)據(jù),共有1 920 000個網(wǎng)格單元,生成200條流線。為測試本文中流線可視化效率,與最新的文獻[8]的流線可視化時間進行對比。對于本例,文獻[8]算法耗時為0.206 2 s,本文耗時為0.141 7 s,節(jié)省時間為31.28%。第二個數(shù)據(jù)是三維湍流數(shù)據(jù),共有82 906 875個網(wǎng)格單元,生成1 000條流線。對于該例,文獻[8]算法耗時50.68 s,本文算法耗時35.83 s,節(jié)省時間為29.30%。由以上兩個測試案例可以看出,本文方法可有效提高矢量場流線可視化的效率。

4 結(jié)束語

針對大規(guī)模矢量場流線可視化效率較低的問題,本文提出一種非結(jié)構(gòu)網(wǎng)格流線的快速可視化方法。該方法首先將流線表示成矩陣形式,從而將流線在網(wǎng)格單元內(nèi)的積分運算轉(zhuǎn)化成矩陣表示的方程與網(wǎng)格單元各表面求交點的問題,并利用牛頓迭代法進行求解。在實驗中發(fā)現(xiàn),對于該方程求解,牛頓迭代法通常4次左右便可收斂,從而可以避免繁瑣的數(shù)值積分運算。實驗結(jié)果表明,相比目前已有的流線可視化方法,本文方法可節(jié)省大約30%的時間。

[1] 李思昆.大規(guī)模流場科學(xué)計算可視化[M].北京:國防工業(yè)出版社,2013.

[2] NOUANESENGSY B, SHENG H W. Load-balanced parallel stre-

amline generation on large scale vector fields[J]. IEEE Tran-

saction on Visualization and Computer Graphics, 2011,17(12):85-94.

[3] SUJUDI D,HAIMISY.Integration of particle paths and streamlines in a spatially-decomposed computation[J]. Parallel Computation-

al Fluid Dynamics,1996,47(3):315-322.

[4] CHEN L,FUJISHIROI L.Optimizing parallel performance of str-

eamline visualization for large distributed flow datasets[J]. IEEE Pacific Visualization Symposium, 2008,21(4):87-94.

[5] CHEN M C, SHA D, HART J D. Fast coherent particle advection through time-varying unstructured flow datasets[J]. IEEE Transactions on Visualization and Computer Graphics, 2016,22(8):1959-1972.

[6] 彭寶云,王文珂,李思昆.大規(guī)模流場矢量線可視化的數(shù)據(jù)預(yù)取方法[J].計算機輔助設(shè)計與圖形學(xué)學(xué)報,2016,28(3):464- 470.

[7] ZHANG J, GUO H,YUAN X. Efficient unsteady flow visualization with high-order access dependencies[J].IEEE Pacific Visualization Symposium (PacificVis 2016), 2016,12(2):56-60.

[8] WANG W,WANG Z K, LI S. Batch advection for the piecewise linear vector field on simplicial grids[J]. Computers & Graphics, 2016,54(2):75-83.