黃筱云，夏波，程永舟，江詩群

(1.長沙理工大學(xué) 水利工程學(xué)院，湖南 長沙 410114； 2.長沙理工大學(xué) 水沙科學(xué)與水災(zāi)害防治湖南省重點實驗室，湖南 長沙 116024)

Level set方法[1]、VOF方法[2]、波前追蹤法(front-tracking method)[3]和相場方法(phase-field method)[4]是追蹤自由表面的4種常用方法。Level set方法的最大優(yōu)勢是能夠容易處理自由表面的拓?fù)渥兓?，但在相同網(wǎng)格分辨率下其計算精度不如顯式的波前追蹤法，其質(zhì)量守恒性弱于VOF法[5]。為彌補level set方法的上述不足，主要的措施包括采用高精度的離散格式(如空間上5階WENO格式、時間上3階TVD Runge-Kutta格式)[6-8]、重構(gòu)level set方程的解[9-11]等。重構(gòu)將強制level set方程的解符合符號距離函數(shù)條件。為實現(xiàn)重新初始化，一種方式是迭代求解偽時間瞬態(tài)重構(gòu)方程[9-12]，另一種是直接求解Eikonal方程?，F(xiàn)階段，Eikonal方程的快速求解方法有快速步進法(fast marching method)[13-15]、快速掃描法(fast sweeping method)[16-18]2類。對于幾何特征簡單的交界面，快速掃描法的計算速度要快于快速步進法[19]。不過，快速掃描法需要對計算全域進行遍歷，而快速步進法只需遍歷自由表面附近的區(qū)域。為進一步提高快速步進重構(gòu)的效率，Herrmann[20]提出了基于區(qū)域分解的快速步進并行算法。Tugurlan[21]和Yang等[22]又分別改進了信息傳遞方式以提高分布式內(nèi)存架構(gòu)下并行計算的效率。黃筱云等[23-24]則提出了基于共享內(nèi)存架構(gòu)的快速步進重構(gòu)算法，解決了相互重疊內(nèi)邊界節(jié)點重構(gòu)值傳遞問題。在此基礎(chǔ)上，夏波等[25]提出了一種自適應(yīng)分區(qū)策略，確保各分區(qū)計算荷載平衡，進一步提高并行重構(gòu)效率。

在level set方法中，重新初始化會定期實施，甚至是每一步都需進行l(wèi)evel set函數(shù)的重構(gòu)。若采用文獻[25]的自適應(yīng)分區(qū)技術(shù)，區(qū)域分解將會帶來額外的計算開銷。因此，本文將充分利用共享內(nèi)存架構(gòu)的優(yōu)勢，提出一種隱式分區(qū)并行策略，避免區(qū)域分解帶來額外開銷。

1 顯式分區(qū)并行重構(gòu)

1.1 快速步進分區(qū)并行重構(gòu)

Level set函數(shù)的重構(gòu)(又稱為重新初始化)被認(rèn)為是求解如下問題：

(1)

式中：Ω表示整個計算區(qū)域；Γ為交界面。快速步進重構(gòu)采用1階迎風(fēng)格式離散方程(1)，即節(jié)點(i,j)上φi,j的空間導(dǎo)數(shù)由該節(jié)點與φ<φi,j的相鄰節(jié)點所決定[26]。離散后的方程(1)則變成一個關(guān)于φi,j的二次方程。若給定那些φ<φi,j相鄰節(jié)點的值，即可求出φi,j的值。因此，選擇從交界面開始，向外依次求解每個節(jié)點上的二次方程，便能求出整個域內(nèi)所有節(jié)點的值。具體操作時，先采用文獻[27]中的方法重構(gòu)交界面周圍網(wǎng)格節(jié)點，將其作為起始邊界；再試算起始邊界外側(cè)所有網(wǎng)格節(jié)點的值，并將這些節(jié)點根據(jù)其值放置到最小堆中，其中，處于根節(jié)點的網(wǎng)格節(jié)點被認(rèn)為滿足方程(1)，或者說重構(gòu)成功；然后，將根節(jié)點上的節(jié)點取出，試算其相鄰未重構(gòu)的節(jié)點，并將這些節(jié)點放置到最小堆中，或調(diào)整它們在最小堆中的位置；最后，將根節(jié)點的子孫節(jié)點中的網(wǎng)格節(jié)點依次向上遞補。

基于區(qū)域分解的并行重構(gòu)預(yù)先將整個計算區(qū)域劃分成p個子區(qū)域，并分配給對應(yīng)的線程，每個線程獨立完成各自子區(qū)域內(nèi)的重構(gòu)。這時，需要求解的問題形式為：

(2)

在文獻[23-25]中，計算區(qū)域?qū)⒀鼐W(wǎng)格線進行分割，故子區(qū)域邊界節(jié)點將相互重疊。每個子區(qū)域獨自實施重構(gòu)無法確保這些重疊的內(nèi)邊界節(jié)點的重構(gòu)值一致，故需要采取同步措施。在顯式分區(qū)并行算法中，當(dāng)重構(gòu)到邊界節(jié)點時，該節(jié)點重構(gòu)值將與相鄰子區(qū)域重疊節(jié)點比較，若其值小于重疊節(jié)點，則將該值傳遞給相鄰子區(qū)域。每個子區(qū)域在接收到傳遞值后，除了更新重疊節(jié)點的值外，還將其作為基準(zhǔn)，回滾域內(nèi)重構(gòu)過的網(wǎng)格節(jié)點。當(dāng)所有子區(qū)域內(nèi)所有節(jié)點都重構(gòu)完畢后，整個區(qū)域的重構(gòu)過程才算結(jié)束。具體過程見文獻[24]。

1.2 自適應(yīng)分區(qū)策略

Level set方法的重新初始化只需在交界面周圍一定范圍內(nèi)執(zhí)行，傳統(tǒng)均分計算區(qū)域的策略可能導(dǎo)致子區(qū)域間重構(gòu)計算開銷不一致，如圖1(a)所示。自適應(yīng)分區(qū)策略則通過平分交界面長度(或面積)來劃分計算區(qū)域，保證子區(qū)域間開銷基本一致，如圖1(b)所示。

圖1 2種顯式分區(qū)策略的比較Fig.1 Comparison between two explicit domain decompositions

自適應(yīng)分區(qū)是一種遞歸剖分過程。即先將計算區(qū)域分成2個子區(qū)域，再分別將子區(qū)域分成2個更小的子區(qū)域，然后繼續(xù)將這些更小的子區(qū)域分成2塊，直至劃分子區(qū)域的數(shù)量與線程數(shù)量相同。此外，文獻[25]中的自適應(yīng)分區(qū)策略還根據(jù)沿不同坐標(biāo)軸方向劃分所產(chǎn)生需要重構(gòu)內(nèi)邊界節(jié)點的數(shù)量來確定子區(qū)域一分為二的操作。這是因為需要重構(gòu)的內(nèi)邊界節(jié)點數(shù)量越少，子區(qū)域間信息傳遞次數(shù)以及域內(nèi)節(jié)點回滾的次數(shù)也越少。

2 隱式分區(qū)并行重構(gòu)

上面所述的并行算法需提前將整個計算區(qū)域沿網(wǎng)格線顯式地劃分成多個子區(qū)域。為保存重疊節(jié)點在不同線程下的計算值，實現(xiàn)相互重疊節(jié)點間的信息比較和傳遞，需要為這些節(jié)點分配額外的存儲空間。當(dāng)頻繁地實施自適應(yīng)分區(qū)并行重構(gòu)時，存儲空間的分配和回收操作會造成不小的計算開銷。

本文提出的隱式分區(qū)并行重構(gòu)不直接劃分計算區(qū)域，而是劃分交界面，將問題(1)將變成：

(3)

這意味著將一個以交界面為起始邊界的重構(gòu)問題變成以交界面一部分為起始邊界的重構(gòu)問題的集合。每個線程將獨立求解一個交界面部分的重構(gòu)問題。圖1所描述的橢圓形交界面邊界重構(gòu)問題若按照隱式分區(qū)并行策略處理，將變成圖2所示的情況。即將交界面分成長度一致的4個部分，每個線程完成以交界面的一部分為起始邊界的重構(gòu)。

圖2 圖1重構(gòu)的隱式分區(qū)并行策略Fig.2 Implicit domain decomposition method for the problem shown in Fig.1

在具體的操作時，無需為每個線程復(fù)制整個計算區(qū)域，而是允許每個線程均能獲取所有節(jié)點。在線程更新網(wǎng)格節(jié)點值之前，應(yīng)預(yù)先比較寫入值與節(jié)點既有值的大小，若寫入值更小，則更新該節(jié)點，當(dāng)多個線程更新同一節(jié)點出現(xiàn)沖突時，可通過OpenMP中的原子操作解決；若節(jié)點值更小，寫入值將不被采納，而且線程將不能通過該節(jié)點向外繼續(xù)求解其他未重構(gòu)的節(jié)點。此時，線程自動將該節(jié)點作為其子區(qū)域的邊界節(jié)點。

與顯式分區(qū)并行重構(gòu)比較，隱式分區(qū)并行重構(gòu)沒有子區(qū)域間節(jié)點值傳遞環(huán)節(jié)和回滾環(huán)節(jié)。隱式分區(qū)并行重構(gòu)方法的特點是允許線程重構(gòu)其他線程重構(gòu)過的節(jié)點，其中，節(jié)點被重新重構(gòu)稱為重演。

在快速步進重構(gòu)法中，節(jié)點根據(jù)重構(gòu)狀態(tài)分為3類：Far表示未重構(gòu)、Close表示重構(gòu)中、Alive表示已重構(gòu)。而在新的隱式分區(qū)并行重構(gòu)算法中，節(jié)點狀態(tài)將被設(shè)定為線程私有。這樣，被一個線程重構(gòu)過的節(jié)點可以被另一個線程繼續(xù)重構(gòu)。當(dāng)一個線程重構(gòu)結(jié)束時，所有節(jié)點會被標(biāo)記成Alive或保留Far標(biāo)記。節(jié)點保留Far標(biāo)記是因為該節(jié)點在其他線程下的重構(gòu)值更小。當(dāng)重構(gòu)結(jié)束時，重構(gòu)的節(jié)點勢必屬于某個線程，即以某個線程的計算結(jié)果作為其重構(gòu)值，而屬于某個線程重構(gòu)節(jié)點的集合構(gòu)成該線程控制的區(qū)域，這樣，整個重構(gòu)區(qū)域分解成多個子區(qū)域。

為了保證線程間重構(gòu)開銷平衡，減小節(jié)點重復(fù)重構(gòu)的次數(shù)，隱式分區(qū)并行重構(gòu)同樣采用自適應(yīng)策略劃分交界面。定義交界面為Γ，線程總數(shù)為NP，具體步驟如下：

1) 令o=1,n=Np,Γo=Γ。

2) 確定交界面Γo的坐標(biāo)中值(xc,yc)。

3) 確定交界面Γ0分別與x=xc及y=yc相交的交點個數(shù)mx和my。

4) 當(dāng)mxxc部分移到Γ0+n/2，其余部分保留在Γ0中；否則，將交界面中y>yc的部分移給Γ0+n/2。

5) 若n>2，則o←o,n←n/2,Γo←Γo和o←o+n/2,n←n/2,Γo←Γo+n/2，

6) 返回步驟2)執(zhí)行遞歸操作。

3 算例及討論

3.1 圓球

圓球位于尺寸為100×100×100的計算區(qū)域內(nèi)，其球心位置為(40,60,60)，半徑為25。重構(gòu)區(qū)域的范圍設(shè)定為φ∈[-12,12]。圖3為8線程下重構(gòu)等值面f=-10, -8, -6, -4, -2, 2, 4, 6, 8,10與平面x=40以及平面z=60的交線。圖4給出了y=50平面上不同線程控制的分區(qū)，可以看出，分區(qū)大小基本相等，且其邊界垂直于交界面邊界，這是因為信息沿交界面法向方向傳播。表1給出φ=10等值面所包圍的體積以及重構(gòu)的計算誤差，其中，體積計算和誤差估計采用文獻[25]中的方法，整體計算精度為1階。圖5比較了隱式分區(qū)和顯式分區(qū)下不同線程數(shù)量帶來的加速比。8線程下顯式分區(qū)并行的加速比達到5，而隱式分區(qū)并行的加速比接近4。雖然相同線程數(shù)量下隱式分區(qū)并行重構(gòu)加速比略小于顯式分區(qū)計算，但前者實際計算時間要小于后者。

圖5 隱式與顯式分區(qū)下圓球并行重構(gòu)加速比與耗時比較Fig.5 Comparison of speedup & runtime between the implicit and explicit domain decomposition methods in parallel reconstruction of the sphere

表1 圓球等值面φ=10并行重構(gòu)的計算損失和誤差Table 1 Computational loss and error of the parallel reconstructed sphere isosurface of φ=10

圖4 平面y=50上不同線程控制的區(qū)域(圓球)Fig.4 Regions governed by different threads at the plane of y=50(sphere)

圖3 8線程下不同等值面的重構(gòu)結(jié)果與平面x=40和平面z=60的交線 (圓球)Fig.3 Intersections between sphere′s reconstructed isosurfaces and two planes of x=40 & z=60 under 8 threads (sphere)

實際上，隱式分區(qū)并行重構(gòu)中的節(jié)點重演和顯式分區(qū)并行算法中的節(jié)點回滾都是為了保證線程間計算的同步。在顯式分區(qū)并行算法中，定義了回滾參數(shù)：

(4)

式中：Nt節(jié)點回滾操作總次數(shù);M計算區(qū)域內(nèi)重構(gòu)節(jié)點總數(shù)。

本文同樣定義一個節(jié)點重演參數(shù)，也用Fr表示，即:

(5)

式中Mr節(jié)點重演操作總次數(shù)。

圖6給出了采用隱式和顯式分區(qū)在不同線程數(shù)量下的Fr值?？梢钥闯觯捎秒[式分區(qū)算法，節(jié)點重演的比例明顯低于顯式分區(qū)算法。

圖6 隱式與顯式分區(qū)下圓球并行重構(gòu)Fr值的對比Fig.6 Comparison of the Fr value between the implicit and explicit domain decomposition methods in parallel reconstruction of the sphere

3.2 Zalesak球

同樣在100×100×100計算區(qū)域的(40,60,60)位置上，放置半徑為25、缺口寬度為10，深度為25, 缺口在y平面上的投影與x軸夾角呈-45°的Zalesak球，如圖7所示。本算例的整個重構(gòu)范圍仍取φ∈[-12,12]。圖7給出了8線程下重構(gòu)等值面f=-10, -8, -6, -4, -2, 2, 4, 6, 8,10與平面y=60以及平面z=60的交線。

圖7 8線程下不同等值面的重構(gòu)結(jié)果與平面y=60和平面z=60的交線(Zalesak圓球)Fig.7 Intersections between Zalesak sphere′s reconstructed isosurfaces and two planes of y=60 & z=60 under 8 threads (Zalesak sphere)

不同線程在y=50平面上的控制分區(qū)如圖8所示，可以看出，與圓球算例同一平面上線程控制區(qū)域分布(圖4)相比，本算例中不同線程控制區(qū)域分布存在差異。φ=5等值面所包圍的體積以及重構(gòu)值的計算誤差見表2。

圖8 平面y=50上不同線程控制的區(qū)域(Zalesak球)Fig.8 Regions governed by different threads at the plane of y=50(Zalesak sphere)

表2 Zalesak球等值面φ=10并行重構(gòu)的計算損失和誤差Table 2 Computational loss and error of parallel reconstruction of the Zalesak′s sphere isosurface of φ=5

Zalesak球隱式與顯式分區(qū)并行重構(gòu)的加速比和計算時間如圖9所示，容易看出，2種并行算法的加速比相近，8線程下的加速比均接近4，但隱式分區(qū)并行重構(gòu)算法的計算時間明顯要少于顯式分區(qū)，幾乎只有前者耗時的1/2。從圖10中同樣可以看出，在Zalesak球等值面的并行重構(gòu)中，隱式分區(qū)的節(jié)點重演次數(shù)也少于顯式算法的節(jié)點回滾次數(shù)，前者的Fr值不到后者的1/3。

圖10 隱式與顯式分區(qū)下Zalesak球并行重構(gòu)Fr值的對比Fig.10 Comparison of the Fr value between the implicit and explicit domain decomposition methods in parallel reconstruction of the Zalesak sphere

圖9 隱式與顯式分區(qū)下Zalesak球并行重構(gòu)加速比與耗時比較Fig.9 Comparison of speedup & runtime between the implicit and explicit domain decomposition methods in parallel reconstruction of the Zalesak sphere

3.3 啞鈴

放置在100×100×100計算區(qū)域內(nèi)的啞鈴由半徑相同的2個球體和1個圓柱體構(gòu)成。2個球體分別位于(40,60,60)和(70,30,40)處，半徑為20，連接2個球心的圓柱的橫截面半徑為10。啞鈴等值面同樣采用8個線程重構(gòu)，其重構(gòu)范圍依然限制在φ∈[-12,12]。圖11給出8線程下啞鈴重構(gòu)等值面f=-10, -8, -6, -4, -2, 2, 4, 6, 8,10與平面y=45以及平面z=50的交線。圖12顯示平面y=50上不同線程控制的計算區(qū)域。

圖11 8線程下不同等值面的重構(gòu)結(jié)果與與平面y=45和平面z=50的交線(啞鈴)Fig.11 Intersections between dumbbell′s reconstructed isosurfaces and two planes of y=45 & z=50 under 8 threads (dumbbell)

圖12 平面y=50上不同線程控制的區(qū)域(啞鈴)Fig.12 Regions governed by different threads at the plane of y=50(dumbbell)

表3給出了φ=5等值面所包圍體積以及重構(gòu)的計算誤差。從圖13可以看出，相同線程數(shù)量下隱式分區(qū)并行計算的加速比略大于顯式分區(qū)，且隱式分區(qū)并行計算耗時幾乎只有顯式分區(qū)的1/2。如圖14所示，在啞鈴等值面的并行重構(gòu)中，隱式分區(qū)的節(jié)點重演次數(shù)也少于顯式算法的節(jié)點回滾次數(shù)。

表3 啞鈴等值面φ=5并行重構(gòu)的計算損失和誤差Table 3 Computational loss and error of the parallel isosurface reconstruction of the dumbbell isosurface of φ=5

圖13 隱式與顯式分區(qū)下啞鈴并行重構(gòu)加速比與耗時比較Fig.13 Comparison of speedup & runtime between the implicit and explicit domain decomposition methods in parallel reconstruction of the dumbbell

圖14 隱式與顯式分區(qū)下啞鈴并行重構(gòu)Fr值的對比Fig.14 Comparison of the Fr value between the implicit and explicit domain decomposition methods in parallel reconstruction of the dumbbell

4 結(jié)論

1)基于共享內(nèi)存架構(gòu)的隱式分區(qū)并行算法能夠有效節(jié)約計算時間。

2)相同線程數(shù)量下，隱式分區(qū)并行算法的加速比與顯式分區(qū)并行算法基本一致，8線程下并行算法的加速比接近4。

3)在相同線程數(shù)量下，隱式分區(qū)并行算法的實際耗時要小于顯式分區(qū)并行算法。

4)隱式分區(qū)并行重構(gòu)的節(jié)點重演次數(shù)要少于顯式分區(qū)并行重構(gòu)的節(jié)點回滾次數(shù)。