劉壯，黃小猛

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

0 引言

近年來，隨著計(jì)算機(jī)計(jì)算能力的提升與數(shù)值方法的發(fā)展，非結(jié)構(gòu)網(wǎng)格在地球系統(tǒng)模式中得到了越來越多的應(yīng)用[1-2]。不同于結(jié)構(gòu)網(wǎng)格，非結(jié)構(gòu)網(wǎng)格沒有極點(diǎn)附近網(wǎng)格收縮而出現(xiàn)奇異的問題；同時(shí)，非結(jié)構(gòu)網(wǎng)格可以更靈活地處理各種邊界條件、實(shí)現(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]是一個(gè)不動(dòng)點(diǎn)迭代算法，其計(jì)算量與生成的網(wǎng)格點(diǎn)數(shù)成正比。為應(yīng)對快速生成高分辨率SCVT 的挑戰(zhàn)，Jacobsen等[4]提出了基于區(qū)域分解和球極投影的并行算法，并在GitHub 開源了其代碼（MPI-SCVT：https：//github.com/douglasjacobsen/MPI-SCVT）。數(shù)值實(shí)驗(yàn)顯示了MPI-SCVT 相對于串行算法的加速效果。

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

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

1 SCVT 生成方法

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

1.1 網(wǎng)格定義

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

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

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

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

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

圖1 Delaunay 三角剖分和Voronoi 剖分

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

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

1.2 Lloyd 算法及并行化

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

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

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

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

算法2 MPI-SCVT 并行算法

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

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

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

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

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

圖3 Voronoi 單元積分計(jì)算方法

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

2.1 初始點(diǎn)集生成

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

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

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

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

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

2.2 解決大變量問題

2.2.1 MPI 通信問題

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

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

2.2.2 Triangle 調(diào)用問題

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

2.2.3 NetCDF 格式問題輸出問題

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

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

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

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

3 結(jié)論

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

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