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

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

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

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

引 言

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

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

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

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

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

1 LICOM 算法基本原理

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

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

如圖1 所示，LICOM 模式運(yùn)行的流程總體上分為三個(gè)階段：

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

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

圖1 LICOM 計(jì)算流程Fig.1 LICOM calculation flowchart

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

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

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

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

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

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

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

1.3 LICOM3 程序性能測(cè)試與分析

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

圖4 各個(gè)模塊的運(yùn)行時(shí)間占比Fig.4 Time proportion of each module

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

2.1 GPU 架構(gòu)

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

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

2.2 CUDA 編程模型

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

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

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

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

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

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

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

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

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

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

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

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

（3）GPU 進(jìn)行模式迭代計(jì)算；

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

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

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

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

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

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

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

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

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

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

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

4 實(shí)驗(yàn)結(jié)果與分析

4.1 實(shí)驗(yàn)平臺(tái)

本文實(shí)驗(yàn)所用的硬件環(huán)境為中科院計(jì)算機(jī)網(wǎng)絡(luò)信息中心新一代超級(jí)計(jì)算機(jī)“元”。本文主要使用“元”上的GPGPU 異構(gòu)計(jì)算系統(tǒng)，該系統(tǒng)共有30 臺(tái)曙光I620-G15，其中每個(gè)GPGPU 計(jì)算節(jié)點(diǎn)配置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 顯示本文所用實(shí)驗(yàn)環(huán)境上單個(gè)節(jié)點(diǎn)配置情況。

4.2 正確性驗(yàn)證

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

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

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

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

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

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

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

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

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

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

4.3 性能分析

本文對(duì)LICOM3 的原CPU 版本和GPU 版本分別通過gettimeofday、cudaEvent 來(lái)獲取運(yùn)行時(shí)間，分別進(jìn)行獨(dú)立的性能測(cè)試。

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

圖16 GPU 版本與原CPU 版本的執(zhí)行時(shí)間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 顯示了對(duì)比原CPU 版本，GPU 版本化前后的加速倍數(shù)，使用1 至8 個(gè)計(jì)算節(jié)點(diǎn)加速了9.31至1.27 倍。圖18 顯示了在使用8 個(gè)計(jì)算節(jié)點(diǎn)的情況下，GPU 版本中各個(gè)模塊耗時(shí)占比，從圖中可以看出tracer、barotr、readyc 模塊是占比最大的三個(gè)模塊，隨著海洋網(wǎng)格塊規(guī)模的減少，它們會(huì)變得更加突出。一是因?yàn)镚PU 適合大規(guī)模加速，當(dāng)海洋網(wǎng)格塊的規(guī)模減小，計(jì)算密集型模塊的計(jì)算密度隨之減?。欢且?yàn)橥ㄐ琶芗K中的通信比增大。這造成了在海洋網(wǎng)格塊規(guī)模小的時(shí)候，數(shù)據(jù)傳輸時(shí)間占比較大、數(shù)據(jù)訪存延遲未被合適的隱藏。圖19 顯示了GPU版本的可擴(kuò)展性，使用2、4、8 個(gè)計(jì)算節(jié)點(diǎn)的加速比分別1.62、2.12、2.5。

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

圖18 GPU 版本各個(gè)模塊耗時(shí)占比（8 個(gè)計(jì)算節(jié)點(diǎn)）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)化方法來(lái)提高數(shù)據(jù)傳輸效率，測(cè)試結(jié)果表明移植優(yōu)化工作取得了一定加速效果。由于LICOM3 計(jì)算的邊界同步通信比較多，造成其可擴(kuò)展性較差，后續(xù)將通過邊界通信優(yōu)化和算法優(yōu)化來(lái)提高模式的可擴(kuò)展性。

利益沖突聲明

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