劉壯，黃小猛

（清華大學(xué) 地球系統(tǒng)科學(xué)系，北京 100084）

0 引言

近年來，隨著計算機計算能力的提升與數(shù)值方法的發(fā)展，非結(jié)構(gòu)網(wǎng)格在地球系統(tǒng)模式中得到了越來越多的應(yīng)用[1-2]。不同于結(jié)構(gòu)網(wǎng)格，非結(jié)構(gòu)網(wǎng)格沒有極點附近網(wǎng)格收縮而出現(xiàn)奇異的問題；同時，非結(jié)構(gòu)網(wǎng)格可以更靈活地處理各種邊界條件、實現(xiàn)局部網(wǎng)格加密等。

在各種非結(jié)構(gòu)網(wǎng)格中，球面質(zhì)心Voronoi 網(wǎng)格[3]（Spherical Centroid Voronoi Tessellations，SCVT）因具備良好的性質(zhì)而得到廣泛應(yīng)用。生成SCVT 一般采用的Lloyd 算法[3]是一個不動點迭代算法，其計算量與生成的網(wǎng)格點數(shù)成正比。為應(yīng)對快速生成高分辨率SCVT 的挑戰(zhàn)，Jacobsen等[4]提出了基于區(qū)域分解和球極投影的并行算法，并在GitHub 開源了其代碼（MPI-SCVT：https：//github.com/douglasjacobsen/MPI-SCVT）。數(shù)值實驗顯示了MPI-SCVT 相對于串行算法的加速效果。

盡管MPI-SCVT 可以很好地適用于一般分辨率下SCVT 的生成，當(dāng)所需網(wǎng)格的分辨率非常高（例如2 km 以下）時，MPI-SCVT 的使用會出現(xiàn)一系列問題。首先，MPI-SCVT 內(nèi)置的幾種網(wǎng)格初始化方法生成的網(wǎng)格點質(zhì)量較差，與最終的收斂結(jié)果距離較遠，這使得迭代步數(shù)過多，程序整體運行時間長。其次，高分辨率導(dǎo)致迭代過程中通信的數(shù)據(jù)量大，MPI-SCVT 使用的boost MPI（http：//www.boost.org） 單次通信數(shù)據(jù)量超過一定大?。s2 GB）時會出錯。再者，MPI-SCVT 調(diào)用外部庫Triangle[5]實現(xiàn)平面Delaunay 三角網(wǎng)格[4]的構(gòu)建，當(dāng)單次輸入網(wǎng)格點數(shù)超過約42 600 000時，Triangle 無法正常運行。最后，高分辨率導(dǎo)致輸出變量較大，使用外部庫NetCDF-CXX 無法完成所需NetCDF 文件的輸出。

本文針對以上問題，提出一套整體解決方案。首先，通過采用正二十面體二分加密的方法，提高初始網(wǎng)格質(zhì)量，減少收斂所需的迭代步數(shù)。其次，通過數(shù)據(jù)分段、多次通信，實現(xiàn)迭代過程中大數(shù)據(jù)的傳輸。再者，通過增加分片處理區(qū)域的數(shù)量，減少每次調(diào)用Triangle 時輸入的網(wǎng)格點數(shù)，順利完成Delaunay 三角網(wǎng)格的構(gòu)建。最后，通過更換Paralell NetCDF 庫（后文簡稱PnetCDF），將大變量數(shù)據(jù)分塊，多次調(diào)用輸出接口進行輸出，消除了原來單變量不能超過4 GB 大小的限制。

1 SCVT 生成方法

這里給出SCVT 的基本定義及生成算法，其詳細介紹可以參考文獻[3-6]。

1.1 網(wǎng)格定義

首先介紹Voronoi 單元。記d（x，y）為球面Ω 上兩點x和y 的測地線距離，給定球面上的一組點，那么Vi={x ∈Ω:d（x，zi）＜d（x，zj） for j=1，…，k，j ≠i}稱為點zi對應(yīng)的Voronoi 單元。

（2）三角形內(nèi)部都不包含其他頂點；

（3）三角形兩兩不相交；

Delaunay 三角剖分是滿足所有三角形的外接圓內(nèi)部都不包含其他頂點的三角剖分。

下面給出Voronoi 單元與Delaunay 三角形的關(guān)系，其詳細證明可以參見文獻[6]。Voronoi 單元與Delaunay 三角形互為對偶：Voronoi 單元的邊是Delaunay 三角形邊的中垂線，Voronoi 單元的頂點是Delaunay 三角形的外心。因此，給定Delaunay 三角剖分，即可構(gòu)建出Voronoi 網(wǎng)格剖分。圖1 給出了Delaunay 三角剖分及其對應(yīng)的Voronoi剖分示意圖。

圖1 Delaunay 三角剖分和Voronoi 剖分

最后，給出SCVT 網(wǎng)格定義。給定Ω 上的密度函數(shù)ρ（x），Voronoi 單元Vi的質(zhì)心定義為：

式中ProjΩ表示通過乘以一個常數(shù)因子將點投影到球面Ω 上。如果Ω 為單位球面，那么：

1.2 Lloyd 算法及并行化

生成SCVT 一般采用Lloyd 算法，如算法1[3]所述。

算法1 生成單位球面SCVT 的Lloyd 算法

輸入：密度函數(shù)ρ（x），網(wǎng)格點數(shù)k。

Jacobsen 等[4]提出了基于區(qū)域分解和球極投影的并行算法，稱為MPI-SCVT，如算法2 所述。

算法2 MPI-SCVT 并行算法

輸入：密度函數(shù)ρ（x），網(wǎng)格點數(shù)k，MPI 進程數(shù)N。

（1）在球面均勻灑N 個點，分別以這N 個點為圓心、以圓心到其鄰居點（與其同屬于一個Delaunay 三角形的點）的距離的最大值為半徑做圓，得到N 個球面區(qū)域，將其分配給N 個進程。

（3）各進程掃描落在本區(qū)域內(nèi)的網(wǎng)格點，之后以本區(qū)域中心的球心對稱點為極點、以此區(qū)域中心處的切平面為投影平面，將本區(qū)域內(nèi)的網(wǎng)格點投影到切平面上（如圖2 所示）。這樣，區(qū)域內(nèi)每個球面網(wǎng)格點即可對應(yīng)于一個平面網(wǎng)格點。然后調(diào)用Triangle 生成平面Delaunay 三角網(wǎng)格。

圖2 球極投影示意圖（t 為區(qū)域中心，f 為其球心對稱點，p 為區(qū)域內(nèi)一點，q 為p 的球極投影）

（4）各進程根據(jù)平面Delaunay 三角網(wǎng)格生成球面Delaunay 三角網(wǎng)格，利用球面Delaunay 三角網(wǎng)格計算Voronoi單元積分。計算方法是用三角形外心與三邊中點的連線將其分為三個子部分，之后計算各子部分的積分，將其累加到對應(yīng)頂點的Voronoi 單元積分，如圖3 所示。為避免重復(fù)計算，每個進程只計算最靠近中心點為本進程負責(zé)區(qū)域中心點的頂點的Voronoi 單元積分，如果某個頂點到兩個區(qū)域中心距離相等，那么選擇進程號小的進程計算。

圖3 Voronoi 單元積分計算方法

2 高分辨率網(wǎng)格生成解決方案

2.1 初始點集生成

如算法1 和算法2 所述，Lloyd 算法是一個不動點迭代算法，因此其收斂速度與初始點集密切相關(guān)。MPISCVT 內(nèi)置了三種網(wǎng)格點初始化方法：Monte Carlo 法、Generalized Spiral 法和Fibbonaci 法。然而，幾種方法生成的初始網(wǎng)格質(zhì)量并不好，與最終收斂結(jié)果距離較遠。

在此，采用球面正二十面體二分加密的方法，如圖4所示。首先給定球面正二十面體網(wǎng)格點，形成球面Delaunay三角網(wǎng)格（G0 網(wǎng)格），然后連接各三角形三邊中點，將原來一個三角形一分為四，得到右上網(wǎng)格（G1 網(wǎng)格）。接著將右上網(wǎng)格的每個三角形用同樣的方法一分為四，得到左下網(wǎng)格（G2 網(wǎng)格），依此進行下去，即可逐步得到高分辨率網(wǎng)格。

圖4 球面正二十面體二分加密方法

表1 中展示了使用幾種初始化方法后達到相同收斂精度所需的迭代步數(shù)。此處判定收斂的變量為相鄰兩次迭代網(wǎng)格點變化量的l1范數(shù)的平均值，收斂參數(shù)分別為1×10-8、1×10-9和1×10-10?？梢钥吹?，對于不同分辨率的網(wǎng)格，二十面體二分加密法（Icosahedron）所需的迭代步數(shù)都遠少于其他幾種方法。值得注意的是，隨著收斂準則變嚴格（即收斂參數(shù)變?。?，其他幾種初始化方法不能收斂的情況變得更多。當(dāng)網(wǎng)格分辨率較高時，只有二十面體二分加密法可以收斂。此測試結(jié)果表明使用二十面體二分加密的初始化方法可以獲得更好的初始網(wǎng)格。

表1 相鄰兩次迭代格點變化量l1 范數(shù)達到收斂條件所需的迭代步數(shù)

2.2 解決大變量問題

2.2.1 MPI 通信問題

MPI-SCVT 初始化及迭代過程中需要使用MPI broadcast 進行數(shù)據(jù)廣播，當(dāng)數(shù)據(jù)大小超過一定量（約2 GB）時會出錯。因此，對broadcast 代碼進行修改，設(shè)置傳輸數(shù)據(jù)數(shù)量上限。對于未超過上限的數(shù)據(jù)像原來一樣進行一次通信，對于超過上限的數(shù)據(jù)，將數(shù)據(jù)分段，進行多次通信，即可解決此問題。

具體地，MPI-SCVT 使用broadcast 廣播的大變量是points，其類型為vector＜pnt＞，pnt 為自定義的class 類型，單個實例所占用內(nèi)存空間為48 B。因此，設(shè)置傳輸數(shù)據(jù)數(shù)量上限參數(shù)num_pts_limit_bcast 為40 000 000。在通信前，主進程將全部網(wǎng)格點數(shù)num_pts 廣播給其他進程。之后在傳輸數(shù)據(jù)前進行判斷，如果num_pts＜=num_pts_limit_bcast，那么直接進行一次通信；否則，將數(shù)據(jù)分段，進行多次通信。此處需要注意的是，對于非源進程，boost MPI的broadcast 函數(shù)對于未分配內(nèi)存的vector，可以根據(jù)源進程vector 大小自動為其分配內(nèi)存并使用接收到的數(shù)據(jù)賦值；對于已分配內(nèi)存的vecotor，將自動覆蓋vector 里面的內(nèi)容。另一方面，對于多次通信的情況，因為每次通信一部分數(shù)據(jù)，所以需要指定points 向量的起始位置，此時非源進程的points 必須是已經(jīng)分配過內(nèi)存的。因此，在進行通信前為未分配內(nèi)存的points 分配內(nèi)存，即可保持代碼統(tǒng)一且簡潔。

2.2.2 Triangle 調(diào)用問題

Triangle 是生成平面Delaunay 三角網(wǎng)格的程序，當(dāng)輸入的網(wǎng)格點數(shù)超過約42 600 000 時無法正常運行。為此，將球面分成更多的區(qū)域（即用更多進程運行），減少每個區(qū)域包含的網(wǎng)格點數(shù)量，即可解決此問題。

2.2.3 NetCDF 格式問題輸出問題

大變量NetCDF 格式文件輸出存在同樣的問題。MPI-SCVT 使用NetCDF-CXX 進行輸出，除了最后一個變量，其余變量大小均不能超過4 GB。為支持更多大變量的輸出，將NetCDF-CXX 庫替換為PnetCDF庫，即可消除變量大小的限制。此處需要注意的是，因為PnetCDF單次輸出數(shù)據(jù)大小限制在2 GB 以內(nèi)，所以對于超過2 GB大小的變量，需要將其分塊，進行多次輸出。

在進行分塊輸出時，對于高維數(shù)據(jù)，要確定在哪個維度方向進行劃分。一般需要在與網(wǎng)格點數(shù)（或者三角形個數(shù)）相關(guān)的維度進行。如果此維度不是變化最慢的維度，那么需要分配一個臨時buffer 數(shù)組將每次要輸出的數(shù)據(jù)拷貝到相應(yīng)的位置上。這樣做的原因是，PnetCDF每次輸出的數(shù)據(jù)必須是一段連續(xù)內(nèi)存的數(shù)據(jù)，如果不進行拷貝，那么無法在對應(yīng)文件的位置寫入正確的數(shù)據(jù)。

現(xiàn)以二維變量分塊輸出為例進行說明。如圖5 所示，假設(shè)num_pts_limit_bcast=10，要輸出的二維數(shù)據(jù)大小為2×10，那么需要將大小為10 的維度劃分，分兩次進行輸出，即每次輸出2×5 的數(shù)據(jù)塊。圖中給出了兩次調(diào)用PnetCDF 輸出接口時輸入的start 和count 數(shù)組。如果變量存儲順序符合C 語言習(xí)慣，即一行一行順序存儲，那么兩次輸出時，分別需要將數(shù)組的白色部分和灰色部分拷貝到一個連續(xù)數(shù)組，將此連續(xù)數(shù)組作為輸入?yún)?shù)；如果變量符合Fortran 語言習(xí)慣，即一列一列順序存儲，則不需要拷貝到臨時數(shù)組，只需分別將兩次輸出的變量起始地址傳遞給PnetCDF 輸出接口。

圖5 PnetCDF 二維數(shù)據(jù)分塊輸出示意圖

3 結(jié)論

本文介紹了生成高分辨率SCVT 的解決方案，重點包括初始點集生成方法以及應(yīng)對大變量問題的幾種方法。通過采用正二十面體二分加密的方法，顯著減少了收斂所需的迭代步數(shù)；將大變量分塊處理，解決了MPI通信、Triangle 調(diào)用以及NetCDF 格式文件輸出的問題。通過采用這一整套方案，可順利生成全球1.9 km 分辨率網(wǎng)格。

最后需要指出的是，本文介紹的正二十面體二分加密初始化方法應(yīng)用領(lǐng)域為生成全球準均勻網(wǎng)格，即密度函數(shù)ρ（x）≡1 的情況。對于變分辨率網(wǎng)格（ρ（x）?1），此方法生成的網(wǎng)格與最終迭代收斂的網(wǎng)格之間的差距并不會比其他幾種方法小，因此不再有優(yōu)勢。對于ρ（x）?1時，如何生成相應(yīng)的初始網(wǎng)格，使其與迭代收斂網(wǎng)格差距盡可能小，是下一步的研究方向。