張留瑩，王鵬飛，張峰，劉海龍，林鵬飛，王濤，韋俊林，田少博， 姜金榮*，遲學(xué)斌

1. 中國科學(xué)院計算機網(wǎng)絡(luò)信息中心，北京 100190

2. 中國科學(xué)院大學(xué)，北京 100049

3. 中國科學(xué)院大氣物理研究所，北京 100029

引 言

海洋環(huán)流模式是一個包含遙感、計算機和動力學(xué)等多學(xué)科研究課題。通過遙感進行監(jiān)測和分析，利用動力學(xué)方法在理論上模擬海洋環(huán)流的變化并建立偏微分方程，通過計算機求解，預(yù)測在未來某個時刻海洋的高度、溫度、鹽度等重要因素。LICOM（LASG/IAP Climate system Ocean Model）是由中國科學(xué)院大氣物理研究所（IAP）大氣科學(xué)和地球流體力學(xué)數(shù)值模擬國家重點實驗室（LASG）自主研發(fā)的全球海洋環(huán)流模式，它是中國科學(xué)院地球系統(tǒng)模式CAS-ESM 的重要組成部分[1,23]，而CAS-ESM 又是開展全球變化研究的主要工具。

IAP 海洋環(huán)流模式經(jīng)歷了長足的發(fā)展，分辨率從20 世紀(jì) 80年代后期第一代模式[2]的5°（赤道上經(jīng)緯度1°約等于 111km）提高至2001年第四代模式LICOM1[3]的0.5°。目前，LICOM 已經(jīng)發(fā)展了三個版本，LICOM3[4]是一個分辨率為7km 的全球版本。模式的分辨率不斷提高使得模式的計算量大幅增加，這導(dǎo)致的一個嚴(yán)重問題是運行代價變的很大，一是運行時間很大，二是在大型機器運行的費用很高[5-6]。高分辨率海洋模式是科學(xué)研究的重要手段，也是提高海洋預(yù)報精度的重要技術(shù)。面向更大規(guī)模計算的全球超高分辨率海洋模式是未來幾十年的研究方向，對海洋環(huán)流模式的并行加速研究已經(jīng)成為了模式發(fā)展的必要條件。

目前國內(nèi)已經(jīng)有很多針對LICOM 進行的并行加速研究。2015年，王文浩等人基于LICOM2 應(yīng)用OpenMP 和指令優(yōu)化等手段進行MIC 移植優(yōu)化[7]，整體加速了2.09 倍。MIC 是一個眾核協(xié)處理器，每個處理器有幾十個核，每個核都是基于X86 架構(gòu)，控制性能與CPU 相似，雖然降低了程序?qū)崿F(xiàn)的難度，但是增加了優(yōu)化的代價。2017年，趙曉溪等人針對LICOM2 進行并行框架SC_Tangram 移植[8]，有效的增加了程序的可維護性和開發(fā)效率，但是程序性能與移植前基本一致。 2019年，Jinrong Jiang 等人基于LICOM2 使用OpenACC 從通信掩藏和循環(huán)優(yōu)化等方面進行GPU 移植優(yōu)化[9]，整體加速了6.6倍。GPU 與MIC 一樣都是眾核協(xié)處理器，不同的是GPU 有成百上千個核，每個核只有計算單元，多個核共用相同的控制單元。相比于MIC，GPU 的并行優(yōu)化更加可控和高效。OpenACC 是基于指令的GPU 并行編程標(biāo)準(zhǔn)，編譯器會自動的管理數(shù)據(jù)移動、啟動并行執(zhí)行和優(yōu)化循環(huán)、將代碼映射到GPU 中。OpenACC 為編程人員提供了一個相對學(xué)習(xí)成本較低的編程模型，但是在性能優(yōu)化方面存在較大的局限性，應(yīng)用程序不能充分的利用硬件資源，與完全優(yōu)化的CUDA 程序相比存在明顯的性能差距。Hoshino T 等人使用兩個微基準(zhǔn)測試和一個計算流體動力學(xué)應(yīng)用程序?qū)Ρ攘薕penACC 和CUDA 的性能差異[10]。兩項評估均表明OpenACC 性能比CUDA 約低50%，對于某些使用CUDA 精心優(yōu)化的應(yīng)用程序，這種性能差異可達98%。CUDA 相對于OpenACC 屬于更底層的編程語言，為編程人員提供了更多的硬件資源控制接口，性能優(yōu)化空間更大，但是程序?qū)崿F(xiàn)難度也更大、周期更長。目前已有的針對LICOM 加速研究成果并沒有充分發(fā)揮GPU 平臺的優(yōu)勢，而且程序的擴展性受到了限制。

GPU（Graphics Processing Unit） 已經(jīng)成為高性能計算系統(tǒng)和科學(xué)應(yīng)用中一個重要的加速解決方案[11-13,23,24]，它是一個通用的、可編程的、數(shù)據(jù)和任務(wù)并行的、低成本的、高性能處理器，CUDA為GPU 提供了更加高效、通用的編程模型。2013年，Huadong Xiao 等人使用CUDA 將天氣預(yù)報模式GRAPES 中最耗時的WSM6 模塊進行GPU 并行加速優(yōu)化[14]，獲得了近140 倍的加速。2012年，Mielikainen 等人通過使用CUDA 將天氣預(yù)報模式WRF 中的Goddard 太陽輻射傳輸模塊進行GPU 移植優(yōu)化[15]，獲得了116 倍的加速。2012年，郭松等人采用CUDA 編程模型將POP 海洋環(huán)流模式移植到GPU 平臺上[16]，獲得了8.47 倍的加速。2013年，王春暉等人使用CUDA 將非靜壓海洋數(shù)值模式進行GPU 加速優(yōu)化[17]，獲得了142 倍的加速。GPU 為大型應(yīng)用的加速提供了可行方案，為LICOM 更加便捷高效的研究提供可能。LICOM 網(wǎng)格點與GPU 線程一一對應(yīng)的關(guān)系使其具有天然的可并行性。與此同時也存在著很多困難，例如LICOM 海洋模式具有變量多、計算分散、分支計算多、計算訪存比低、上下文依賴嚴(yán)重、通信同步多、代碼量大等特點。

本文使用CUDA C 對LICOM3 進行GPU 移植，將網(wǎng)格點映射為CUDA 線程進行加速，并從主機與設(shè)備的數(shù)據(jù)傳輸、GPU 內(nèi)存使用方面對并行算法進行了優(yōu)化。第1 節(jié)介紹LICOM 程序的基本計算流程和特點；第2 節(jié)介紹GPU 架構(gòu)、CUDA 編程模型和LICOM 在GPU 上實現(xiàn)；第3 節(jié)介紹LICOM 在GPU 上的優(yōu)化手段；第4 節(jié)介紹實驗環(huán)境，進行正確性驗證和性能分析。

1 LICOM 算法基本原理

LICOM 數(shù)值模式可以歸結(jié)為利用離散的時間和空間差分求解一組偏微分方程[18]，包括動量方程、熱力學(xué)方程、鹽度方程、海水狀態(tài)方程和太陽輻射、長波輻射等一系列參數(shù)化過程。海洋環(huán)流主要受海洋表面的風(fēng)力作用影響和海水溫鹽度變化的影響，同時還受地球重力、地球自轉(zhuǎn)和不同的波動等多方面的影響，其中重力外波速度快尺度小、重力內(nèi)波速度慢尺度大等特點，使得模式的動量方程分為正壓模態(tài)（Barotropic）和斜壓模態(tài)（Baroclinic）兩個獨立的部分，正壓模態(tài)負責(zé)二維外波積分，斜壓模態(tài)負責(zé)三維內(nèi)波積分。

1.1 LICOM 程序結(jié)構(gòu)

如圖1 所示，LICOM 模式運行的流程總體上分為三個階段：

（1）預(yù)處理階段。主要進行處理器設(shè)置和數(shù)據(jù)劃分等；

（2）模式初始化階段。主要進行初始數(shù)據(jù)的導(dǎo)入和迭代變量的設(shè)置[19]，通過Netcdf 將nc 格式的源數(shù)據(jù)讀入模式；

圖1 LICOM 計算流程Fig.1 LICOM calculation flowchart

（3）模式迭代計算階段。是整個模式的核心計算部分，是本文進行GPU 移植目標(biāo)程序。核心的7 個積分計算模塊分別為正壓斜模態(tài)積分的預(yù)處理（readyc 模塊，readyt 模塊）、正壓模態(tài)積分（barotr 模 塊）、斜壓模態(tài)積分（bclinc 模塊）、熱鹽積分（tracer模塊）、海冰積分（icesnow 模塊）、GFDL 全對流調(diào)整（convadj 模塊），包含在月數(shù)循環(huán)、月內(nèi)循環(huán)、日內(nèi)時間步循環(huán)這三層循環(huán)內(nèi)。正壓模態(tài)積、斜壓模態(tài)積分和熱鹽積分采用不同時間步長的蛙跳形式進行迭代[18]，每個模式天的積分迭代中，正壓模態(tài)每積分Nbarotr 時間步（步長為Tbarotr），斜壓模態(tài)就積分一個時間步（步長為Tbclinc=Nbarotr×Tbarotr）；在正斜壓模態(tài)耦合積分Nbclinc 個時間步后，熱鹽積分一個時間步（步長為Ttracer = Nbclinc × Tbclinc）。

1.2 LICOM 海洋網(wǎng)格塊并行方案

海洋在水平方向以交錯的經(jīng)緯度為坐標(biāo)被離散為的二維網(wǎng)格點，緯度方向分為imt_global個網(wǎng)格點，經(jīng)度方向分為jmt_global個網(wǎng)格點。垂直方向以海洋深度為坐標(biāo)分為km層，海洋被離散成一個三維結(jié)構(gòu)，在這基礎(chǔ)上將海洋網(wǎng)格點組織為網(wǎng)格塊，每個網(wǎng)格塊包含了BLCKX×BLCKY×km個網(wǎng)格點。網(wǎng)格塊與網(wǎng)格塊之間存在正壓模態(tài)通信，層與層之間會存在相互的壓力和垂直速度的計算。如圖2 所示，整個海洋被分為了ceil(imt_global/BLCKX)×ceil(jmt_global/BLCKY)個海洋網(wǎng)格塊，其中ceil 為向上取整函數(shù)。為了提高差分求解的效率，LICOM 采用MPI進行網(wǎng)格塊并行加速，為每個進程分配NBLOCKS_CLINIC個網(wǎng)格塊，網(wǎng)格塊內(nèi)部的網(wǎng)格點串行計算。

圖2 網(wǎng)格塊劃分Fig.2 Ocean grid

每個網(wǎng)格有自己獨立的計算數(shù)組，數(shù)組的規(guī)模隨著并行度的增大而減小。因為差分計算需要鄰接網(wǎng)格的邊界值，所以需要頻繁的進行網(wǎng)格邊界數(shù)據(jù)通信，每個網(wǎng)格的計算數(shù)組根據(jù)二維邏輯上相鄰進程的數(shù)組邊界來更新。所以每個網(wǎng)格的計算數(shù)組除了存儲本網(wǎng)格內(nèi)的數(shù)據(jù)外還存儲著鄰接網(wǎng)格的邊界數(shù)據(jù)（又稱為“偽邊界”，本地數(shù)組的邊界稱為“實邊界”）。LICOM 采用“可伸縮偽邊界”的策略[18]，即偽邊界的條數(shù)可根據(jù)具體問題來確定。如圖2 右的深色部分顯示了海洋網(wǎng)格塊的“偽邊界”，在程序中聲明為變量ghost并賦值為2。相比于一條偽邊界，將偽邊界設(shè)置為兩條，可以進行兩次計算數(shù)組更新只進行通信一次，能夠有效的降低通信 計算比。

圖3 顯示了在一次迭代中網(wǎng)格點的計算過程[18]，首先進行網(wǎng)格塊兩條偽邊界的更新：網(wǎng)格塊1 中的實邊界c 和d 分別傳遞給網(wǎng)格塊2 的偽邊界b 和a，網(wǎng)格塊2 中的實邊界c 和d 傳遞給網(wǎng)格塊1 的偽邊界b 和a；然后進行實邊界更新：網(wǎng)格塊1 內(nèi)部向右更新至b、網(wǎng)格塊2 內(nèi)部向左更新至b；最后進行內(nèi)部網(wǎng)格點的更新：網(wǎng)格塊1 內(nèi)部向右更新至c、網(wǎng)格塊2 內(nèi)部向左更新至c。

圖3 網(wǎng)格塊邊界通信Fig.3 Ocean grid boundary communication

1.3 LICOM3 程序性能測試與分析

本文對LICOM3 原CPU 程序進行性能測試并分析熱點模塊。如圖4 所示，分別使用1（20）、2（40）、4（80）、8（160）個計算節(jié)點（MPI 進程數(shù)量）測試了單個模式天內(nèi)各個模塊迭代計算的耗時占比。我們可以看出隨著計算節(jié)點個數(shù)的增加、程序的并行度在增加、海洋網(wǎng)格塊的大小在減小，整個模式的熱點模塊也隨之改變。對于readyc 這樣計算密集的模塊，在計算節(jié)點個數(shù)較少的情況下成為程序中影響最大的因素。對于barotr、bclinc、tracer等通信密集的模塊，隨著并行度的增加，單個進程通信量的占比在增加，在計算節(jié)點個數(shù)較多的情況下成為程序中影響最大的因素，限制著程序的可擴展性。對于icesnow、convadj 這兩個模塊，運行時間占比一直小于1%，對整個LICOM 程序性能的影響比較小。

圖4 各個模塊的運行時間占比Fig.4 Time proportion of each module

2 LICOM 在GPU 上的實現(xiàn)

2.1 GPU 架構(gòu)

GPU 是一種眾核架構(gòu)處理器[20]，如圖5 所示，每個GPU 卡包含多個流式多處理器SM（Streaming Multiprocessors），每個SM 包含多個流處理器SP（Streaming Processors）、共享內(nèi)存等部件。以Tesla K20 為例，每個處理器包含15 組SM，每組SM 包含192 個SP。一個SM 上可以同時承載著多個線程塊（Block），可同時并發(fā)執(zhí)行數(shù)百個線程（Thread）。每個Block 內(nèi)的Thread 被劃分為每32 個一組的線程束（wrap），以wrap 為單位進行調(diào)度。

GPU 是CPU 的一種協(xié)處理器，并不是一個獨立運行的平臺，必須通過PCIe 總線與CPU 相連，所以我們將CPU 稱為主機端（Host），GPU 稱為設(shè)備端（Device）。CPU+GPU 異構(gòu)系統(tǒng)上執(zhí)行的應(yīng)用通常由CPU 進行數(shù)據(jù)初始化、設(shè)備端內(nèi)存分配、數(shù)據(jù)拷貝和核函數(shù)的啟動，GPU 負責(zé)計算密集的程序部分的加速。這樣的異構(gòu)系統(tǒng)能夠與充分發(fā)揮各自的特長，極大程度上加速了應(yīng)用程序 的計算。

2.2 CUDA 編程模型

CUDA 編程模型[20]使得編程人員能夠使用熟悉的C 語言實現(xiàn)核函數(shù)（Kernel ），多個Kernel 可以重疊運行在同一個GPU 上。由于CUDA 編程模型是異步的，Kernel 也可以與CPU 計算重疊執(zhí)行。為了保持同步，可以顯式的通過調(diào)用cudaDeviceSynchronize函數(shù)來實現(xiàn)，也可以隱式的通過調(diào)用變量拷貝cudaMemcpy 函數(shù)來實現(xiàn)。

CUDA 程序是靈活、可擴展的，線程按照線程、線程塊、線程格的多層次模型組織起來，通過各個層次的編號可以唯一的標(biāo)識某個線程。一個Kernel開啟一個線程格（Grid），Grid 會被組織成三維的線程塊，相同Grid 內(nèi)的Block 可以通過blockIdx(x, y, z)唯一的標(biāo)識。Block 被組織成三維的線程，相同Block 內(nèi)的Thread 通過threadIdx(x, y, z)唯一標(biāo)識。

CUDA 存儲模型按照訪問速度從快到慢分別為寄存器、共享內(nèi)存、常量內(nèi)存、全局內(nèi)存和局部內(nèi)存。如圖6 所示，寄存器和共享內(nèi)存都是片上存儲空間，供SM 上的所有活躍線程使用，其中寄存器是不可編程的且屬于線程私有，共享內(nèi)存由同一個Block 內(nèi)的線程共享。全局內(nèi)存、局部內(nèi)存和常量內(nèi)存都位于片外顯存中，其中全局內(nèi)存和常量內(nèi)存由同一個Grid 內(nèi)的線程共享，而局部內(nèi)存屬于線程私有且不能合并訪存。

圖6 CUDA 編程模型Fig.6 CUDA programming model

2.3 LICOM 在GPU 上的實現(xiàn)

本文將LICOM 迭代計算過程中的7 個模塊使用CUDA C 進行GPU 移植，其他的預(yù)處理部分、模塊初始化部分以及MPI 通信部分仍保留原來的Fortran 程序。在MPI 并行計算海洋網(wǎng)格塊的基礎(chǔ)之上，結(jié)合使用CUDA 線程并行計算海洋網(wǎng)格點。如圖7所示，每個MPI進程負責(zé)一個海洋網(wǎng)格塊的計算，海洋網(wǎng)格塊內(nèi)的網(wǎng)格點由GPU 線程并行執(zhí)行。

圖7 海洋網(wǎng)格點與CUDA 線程的映射Fig.7 Mapping of ocean grid points to CUDA threads

表1 和表2 顯示了LICOM 程序片段的Fortran實現(xiàn)與CUDA 實現(xiàn)的對應(yīng)關(guān)系，CPU 中的全局變量映射到GPU 上的全局內(nèi)存中，例如CPU 全局變量VIV，在GPU 申請為全局內(nèi)存變量d_viv。

表2 CUDA 版本的程序片段Table 2 CUDA version of vinteg function

主機和設(shè)備的數(shù)據(jù)傳輸是整個應(yīng)用程序的瓶頸，為了減少CPU-GPU 頻繁的內(nèi)存拷貝，本文將GPU的全局內(nèi)存申請和初始數(shù)據(jù)的拷貝前置到模式初始化階段。在模式迭代計算的中間，遇到MPI 通信（如barotr 模塊）的時候，將GPU 變量顯式的拷貝到CPU，進行MPI 通信后再將數(shù)據(jù)拷貝到GPU 上繼續(xù)進行迭代計算。當(dāng)模式計算完成后再將數(shù)據(jù)拷貝回CPU 輸出計算結(jié)果。但是中間的MPI 通信是必要的，同時也影響著整個程序的效率。CUDA 版本的LICOM 程序?qū)崿F(xiàn)過程大致如下：

（1）CPU 進行預(yù)處理和模式初始化；

（2）GPU 進行內(nèi)存分配，并且將模式變量初始值拷貝至GPU；

（3）GPU 進行模式迭代計算；

（4）GPU 把計算結(jié)果拷貝回CPU，輸出結(jié)果。

3 LICOM 在GPU 上的優(yōu)化

在GPU 上進行程序的性能優(yōu)化的關(guān)鍵是合理的減少內(nèi)存訪問延遲和設(shè)備主機的數(shù)據(jù)傳輸延遲[21]。通過nvprof 性能分析工具測試核函數(shù)和數(shù)據(jù)傳輸?shù)臅r間占比，如圖8 所示，主機和設(shè)備數(shù)據(jù)傳輸（memcpyHtoD 和memcpyDtoH）占總時間的42%，核函數(shù)的執(zhí)行占總時間的58%。

圖8 核函數(shù)和數(shù)據(jù)傳輸?shù)臅r間比例Fig.8 Proportion of kernel and data transmission

3.1 數(shù)據(jù)傳輸優(yōu)化

本文通過使用固定內(nèi)存和CUDA 流進行數(shù)據(jù)傳輸優(yōu)化，減少程序中因MPI 通信帶來的數(shù)據(jù)傳輸?shù)拈_銷。主機上的變量一般存儲在可分頁內(nèi)存（Pageable Host Memory）中，GPU 設(shè)備不可以直接訪問可分頁內(nèi)存。如圖9 左所示，在進行主機和設(shè)備的數(shù)據(jù)傳輸之前，需要開辟一個臨時的固定內(nèi)存（Pinned Host Memory），將需要傳輸?shù)臄?shù)據(jù)從分頁內(nèi)存中拷貝到固定內(nèi)存，最后才將數(shù)據(jù)從固定內(nèi)存?zhèn)鬏數(shù)皆O(shè)備內(nèi)存中。

零拷貝內(nèi)存其實是一塊主機上的固定內(nèi)存，通過使用零拷貝內(nèi)存，避免了從分頁內(nèi)存到固定內(nèi)存的拷貝過程。同時零拷貝內(nèi)存映射到了設(shè)備地址中，允許設(shè)備直接訪問，避免了顯式的數(shù)據(jù)傳輸，而且CUDA 自動的利用核函數(shù)的執(zhí)行隱藏了數(shù)據(jù)傳輸。但是零拷貝內(nèi)存適合少量的讀寫，每一次的讀寫都會產(chǎn)生一次數(shù)據(jù)傳輸，過多的使用或者不合理的使用零拷貝內(nèi)存可能會帶來性能的下降。

圖9 數(shù)據(jù)傳輸Fig.9 Data transmission

CUDA 流是一系列在設(shè)備端異步執(zhí)行操作，同一個流內(nèi)操作按照啟動順序執(zhí)行，不同的流之間相互獨立。結(jié)合cudaMemcpyAsync 函數(shù)，通過CUDA流的控制也可以達到隱藏數(shù)據(jù)傳輸。

3.2 內(nèi)存訪問優(yōu)化

本文將核函數(shù)中頻繁訪問的全局變量優(yōu)化為寄存器內(nèi)存，用寄存器內(nèi)存暫存中間計算量。對于通過cudaMemset 將全局變量初始化為0 的部分，優(yōu)化為寄存器內(nèi)存并在Kernel 內(nèi)部進行初始化。通過CUDA 提供的cudaDeviceProp 可以查詢到本文實驗所用的Tesla K20 允許每個Block（或SM）最多可以使用65536 個寄存器（每個寄存器內(nèi)存為32-bit）。在核函數(shù)中定義局部變量，在少量使用局部變量的情況下會分配寄存器內(nèi)存，寄存器內(nèi)存的訪問延遲較低。但是過多的定義局部變量會適得其反，因為當(dāng)寄存器的數(shù)量不夠用的時候，將分配訪問速度較慢且不能合并訪存的局部內(nèi)存。所以適當(dāng)?shù)膶⒑撕瘮?shù)局部數(shù)組變量改變?yōu)槿謨?nèi)存變量，程序能夠有效的合并訪存，也可以帶來性能提升。

4 實驗結(jié)果與分析

4.1 實驗平臺

本文實驗所用的硬件環(huán)境為中科院計算機網(wǎng)絡(luò)信息中心新一代超級計算機“元”。本文主要使用“元”上的GPGPU 異構(gòu)計算系統(tǒng)，該系統(tǒng)共有30 臺曙光I620-G15，其中每個GPGPU 計算節(jié)點配置2 顆 Nvidia Tesla K20、2 顆 Intel Xeon E5-2680 V2(10 core |2.8 GHz)。本文使用Intel MPI 4.1.3.049搭建MPI環(huán)境、CUDA 6.0.37 搭建CUDA 環(huán)境。通過Intel MPI 提供的mpiifort 編譯器編譯Fortran 程序、CUDA 提供的nvcc 編譯器編譯CUDA 程序。表3 顯示本文所用實驗環(huán)境上單個節(jié)點配置情況。

4.2 正確性驗證

本文對LICOM3 的原CPU 版本和GPU 版本進行獨立的測試，對兩個版本的模式運行至第1年2月1日瞬時的計算結(jié)果進行做差，分別從海表高度（SSH，Sea Surface Height）、海表溫度（SST，Sea Surface Temperatures）、海表鹽度（SSS，Sea Surface Salinity）三個方面進行驗證。從圖10-15 可以看出，GPU 版本與原CPU 版本的模擬過程基本一致，計算偏差均在可接受的范圍內(nèi)。

圖10 顯示了LICOM3 原CPU 版本運行至第1年2月1日的瞬時SSH。圖11 顯示了GPU 版本在第1年2月1日 的瞬時SSH 差。GPU 版本的SSH計算結(jié)果絕對誤差平均值為-0.000080145。

圖12 顯示了LICOM3 原CPU 版本運行至第1年2月1日的瞬時SST。圖13 顯示了GPU 版本在第1年2月1日的瞬時SST 差。 GPU 版本的SST計算結(jié)果絕對誤差平均值為-0.00027874。

圖10 CPU 版本的瞬時SSH（米）Fig.10 Instantaneous SSH of CPU version (meter)

圖11 GPU 版本的瞬時SSH 差（米）Fig.11 Instantaneous SSH deviation of GPU version (meter)

圖12 CPU 版本的瞬時SST（℃）Fig.12 Instantaneous SST of CPU version (°C)

圖13 GPU 版本的瞬時SST 差（℃）Fig.13 Instantaneous SST deviation of GPU version (°C)

圖14 顯示了LICOM3 原CPU 版本運行至第1年2月1日的瞬時SSS，圖15 顯示了GPU 版本在第1年2月1日的瞬時SSS 差， GPU 版本的SSS 計算結(jié)果絕對誤差平均值為-0.0021483。

圖14 CPU 版本的瞬時 SSS（PSU）Fig.14 Instantaneous SSS of CPU version (PSU)

圖15 GPU 版本的瞬時SSS 差（PSU）Fig.15 Instantaneous SSS deviation of GPU version (PSU)

4.3 性能分析

本文對LICOM3 的原CPU 版本和GPU 版本分別通過gettimeofday、cudaEvent 來獲取運行時間，分別進行獨立的性能測試。

圖16 顯示了在使用1（20）、2（40）、4（80）、 8（160）個計算節(jié)點（MPI 進程數(shù)量）的情況下，原CPU 版本和GPU 版本單個模式天的運行時間。在使用1 個計算節(jié)點的情況下（開辟20 個進程），原CPU 版本的單個模式天運行時間為102.21 秒，CUDA 版本的運行時間下降到了10.98 秒。

圖16 GPU 版本與原CPU 版本的執(zhí)行時間Fig.16 Runtime of the GPU version and the CPU version

圖17 優(yōu)化前后GPU 版本比原CPU 版本的加速倍數(shù)Fig.17 Speedup of the before and after optimization GPU version compared to the CPU version

圖17 顯示了對比原CPU 版本，GPU 版本化前后的加速倍數(shù)，使用1 至8 個計算節(jié)點加速了9.31至1.27 倍。圖18 顯示了在使用8 個計算節(jié)點的情況下，GPU 版本中各個模塊耗時占比，從圖中可以看出tracer、barotr、readyc 模塊是占比最大的三個模塊，隨著海洋網(wǎng)格塊規(guī)模的減少，它們會變得更加突出。一是因為GPU 適合大規(guī)模加速，當(dāng)海洋網(wǎng)格塊的規(guī)模減小，計算密集型模塊的計算密度隨之減??；二是因為通信密集模塊中的通信比增大。這造成了在海洋網(wǎng)格塊規(guī)模小的時候，數(shù)據(jù)傳輸時間占比較大、數(shù)據(jù)訪存延遲未被合適的隱藏。圖19 顯示了GPU版本的可擴展性，使用2、4、8 個計算節(jié)點的加速比分別1.62、2.12、2.5。

實驗結(jié)果表明，GPU 版本相對于原CPU 版本取得一些加速效果，并具有一定的可擴展性。今后的工作將進一步的擴大數(shù)據(jù)規(guī)模，在更高分辨率的模式上進行測試。為了進一步的增加程序的可擴展性，可以從模式通信算法方面進一步的優(yōu)化，減少通信次數(shù)；還可以通過使用支持CUDA-aware MPI的GPU 平臺，使得通信函數(shù)可以直接訪問GPU 內(nèi)存，避免了數(shù)據(jù)傳輸?shù)幕ㄤN。

圖18 GPU 版本各個模塊耗時占比（8 個計算節(jié)點）Fig.18 Time proportion of each module in GPU version (8 nodes)

圖19 GPU 版本加速比Fig.19 The speedup of GPU version

5 結(jié)論

本文初步完成了全球海洋環(huán)流模式LICOM3 的CUDA C 移植，并利用固定內(nèi)存和CUDA 流等優(yōu)化方法來提高數(shù)據(jù)傳輸效率，測試結(jié)果表明移植優(yōu)化工作取得了一定加速效果。由于LICOM3 計算的邊界同步通信比較多，造成其可擴展性較差，后續(xù)將通過邊界通信優(yōu)化和算法優(yōu)化來提高模式的可擴展性。

利益沖突聲明

所有作者聲明不存在利益沖突關(guān)系。