Level set函數快速步進重構的隱式分區并行

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

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

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

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

1 顯式分區并行重構

1.1 快速步進分區并行重構

Level set函數的重構(又稱為重新初始化)被認為是求解如下問題：

(1)

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

基于區域分解的并行重構預先將整個計算區域劃分成p個子區域，并分配給對應的線程，每個線程獨立完成各自子區域內的重構。這時，需要求解的問題形式為：

(2)

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

1.2 自適應分區策略

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

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

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

2 隱式分區并行重構

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

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

(3)

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

圖2 圖1重構的隱式分區并行策略Fig.2 Implicit domain decomposition method for the problem shown in Fig.1

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

與顯式分區并行重構比較，隱式分區并行重構沒有子區域間節點值傳遞環節和回滾環節。隱式分區并行重構方法的特點是允許線程重構其他線程重構過的節點，其中，節點被重新重構稱為重演。

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

為了保證線程間重構開銷平衡，減小節點重復重構的次數，隱式分區并行重構同樣采用自適應策略劃分交界面。定義交界面為Γ，線程總數為NP，具體步驟如下：

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

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

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

4) 當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)執行遞歸操作。

3 算例及討論

3.1 圓球

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

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

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

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

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

實際上，隱式分區并行重構中的節點重演和顯式分區并行算法中的節點回滾都是為了保證線程間計算的同步。在顯式分區并行算法中，定義了回滾參數：

(4)

式中：Nt節點回滾操作總次數;M計算區域內重構節點總數。

本文同樣定義一個節點重演參數，也用Fr表示，即:

(5)

式中Mr節點重演操作總次數。

圖6給出了采用隱式和顯式分區在不同線程數量下的Fr值。可以看出，采用隱式分區算法，節點重演的比例明顯低于顯式分區算法。

圖6 隱式與顯式分區下圓球并行重構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計算區域的(40,60,60)位置上，放置半徑為25、缺口寬度為10，深度為25, 缺口在y平面上的投影與x軸夾角呈-45°的Zalesak球，如圖7所示。本算例的整個重構范圍仍取φ∈[-12,12]。圖7給出了8線程下重構等值面f=-10, -8, -6, -4, -2, 2, 4, 6, 8,10與平面y=60以及平面z=60的交線。

圖7 8線程下不同等值面的重構結果與平面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平面上的控制分區如圖8所示，可以看出，與圓球算例同一平面上線程控制區域分布(圖4)相比，本算例中不同線程控制區域分布存在差異。φ=5等值面所包圍的體積以及重構值的計算誤差見表2。

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

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

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

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

圖9 隱式與顯式分區下Zalesak球并行重構加速比與耗時比較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計算區域內的啞鈴由半徑相同的2個球體和1個圓柱體構成。2個球體分別位于(40,60,60)和(70,30,40)處，半徑為20，連接2個球心的圓柱的橫截面半徑為10。啞鈴等值面同樣采用8個線程重構，其重構范圍依然限制在φ∈[-12,12]。圖11給出8線程下啞鈴重構等值面f=-10, -8, -6, -4, -2, 2, 4, 6, 8,10與平面y=45以及平面z=50的交線。圖12顯示平面y=50上不同線程控制的計算區域。

圖11 8線程下不同等值面的重構結果與與平面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上不同線程控制的區域(啞鈴)Fig.12 Regions governed by different threads at the plane of y=50(dumbbell)

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

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

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

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

4 結論

1)基于共享內存架構的隱式分區并行算法能夠有效節約計算時間。

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

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

4)隱式分區并行重構的節點重演次數要少于顯式分區并行重構的節點回滾次數。