賈志強(qiáng)

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

矢量場可視化是科學(xué)可視化的重要分支,在油氣田油層分析、氣象預(yù)報(bào)、海洋數(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ù),利用計(jì)算機(jī)集群或者超級計(jì)算機(jī),對流場數(shù)據(jù)和可視化任務(wù)進(jìn)行劃分后分派給各個(gè)節(jié)點(diǎn),多個(gè)節(jié)點(diǎn)同時(shí)進(jìn)行計(jì)算。流線并行主要有數(shù)據(jù)塊并行和種子點(diǎn)并行兩種方式。數(shù)據(jù)塊并行的典型方法有:Nouanesengsy等[2]提出的負(fù)載敏感數(shù)據(jù)劃分策略、Sujudi等[3]提出的空間分解方法等;種子點(diǎn)并行的典型方法有:Chen等[4]提出的最優(yōu)化動(dòng)態(tài)數(shù)據(jù)劃分策略、Chen等[5]提出的空間一致性方法等。此外,彭寶云等[6]以及Zhang等[7]利用數(shù)據(jù)預(yù)取技術(shù)進(jìn)行流線高效可視計(jì)算的方法。以上方法雖然從各種角度對流線可視化過程進(jìn)行加速,但是在流線生成過程中都是利用數(shù)值積分進(jìn)行計(jì)算。針對該問題,筆者首先利用代數(shù)運(yùn)算將流線表示成矩陣形式,進(jìn)而將流線在網(wǎng)格單元內(nèi)的積分運(yùn)算轉(zhuǎn)化成矩陣表示的方程與網(wǎng)格單元各表面求交點(diǎn)的問題,并利用牛頓迭代法進(jìn)行求解,從而避免繁瑣的數(shù)值積分運(yùn)算,有效提高大規(guī)模矢量場流線可視化的效率。

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

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

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

v(x)=Ax+b.

非結(jié)構(gòu)網(wǎng)格由于其表達(dá)的靈活性,目前得到越來越廣泛的應(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)的表達(dá)式,將矩陣A進(jìn)行特征值分解。利用線性代數(shù)知識可以證明,當(dāng)A有3個(gè)不同的非零實(shí)特征值,或者兩個(gè)復(fù)數(shù)特征值、一個(gè)非零實(shí)特征值時(shí),流線方程(2)都可以表示成如下的形式:

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

(3)

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

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

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

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

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

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

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

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

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

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

2.3 普通單元流線出點(diǎn)的處理

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

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

(4)

式中,×表示向量的叉積、·表示向量的點(diǎn)積。

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

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

輸入:流線種子點(diǎn)p

輸出:由p生成的流線

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

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

while(cell不為空)

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

caseparallel//平行單元

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

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

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

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

if(tmp

t=tmp

id=i

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

break

casenormal//普通單元

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

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

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

if(tmp>0 &&tmp

t=tmp

id=i

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

break

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

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

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找到相鄰的下一個(gè)網(wǎng)格單元

依次連接E中的各點(diǎn),得到流線并輸出

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

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

3 算法驗(yàn)證

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

4 結(jié)束語

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

[1] 李思昆.大規(guī)模流場科學(xué)計(jì)算可視化[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].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),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.