基于OpenACC的GRAPES_GLOBAL模式長波輻射異構并行化研究

孫晨 王彬 顧文靜 魏敏

（國家氣象信息中心，北京 100081）

0　引言

GRAPES（Global/Regional Assimilation and Prediction System）是中國氣象局自主研發的靜力/非靜力多尺度通用數值預報系統，該系統是氣象與氣候研究的基礎和核心。

GRAPES系統的核心部分包括模式動力框架和物理過程，動力框架中計算時間最長的為半拉格朗日計算和求解Helmholtz大型線性代數方程組，物理過程中耗時較大的為輻射和微物理過程，其中輻射過程計算耗時占到GRAPES模式的40%，所以提升輻射過程的性能和計算效率，對整個GRAPES模式的性能提升具有重要科研價值和實際意義[1-2]。

近年來，并行計算技術日趨成熟，GPU巨大的計算能力已經吸引了越來越多的科研人員將其作為一個高效益低成本的高性能計算平臺，然而，要在GPU硬件上實現加速需要通過對底層的API進行編程來實現，程序編寫復雜、難度大，且容易形成高度依賴特定設備的代碼。因此，一種基于指令制導計算的編程模OpenACC應運而生[3]，OpenACC的編程機制是在源程序中添加少量編譯標識，編譯器根據作者的意圖自動產生低級語言代碼，無需學習新的編程語言和加速器硬件知識，變更迅速掌握。只添加少量編譯標識不破壞原代碼，開發速度快，即可并行執行又可恢復串行執行；硬件更新時，重新編譯一次代碼即可，不必手工修改代碼。OpenACC標準制定時就考慮了目前及將來多種加速器產品，同一份代碼可以在多種加速器設備上編譯、運行、無成本切換硬件平臺。

本文以GRAPES模式物理過程中的長波輻射為加速處理對象，依托國家超級計算無錫中心的“神威·太湖之光”高性能計算機系統，設計一種基于OpenACC的加速算法，最大程度繼承原始代碼的MPI 級并行，同時利用數量眾多的協處理器加速熱點函數。

1　“神威·太湖之光”簡介

2016年6月20日，新一期全球超級計算機500強榜單在德國法蘭克福舉辦的國際超算大會（ISC）上公布，“神威·太湖之光”超級計算機榮登榜首。

“神威·太湖之光”由國家并行計算機工程技術研究中心研制，是全球第一臺運行速度超過10億億次/s的超級計算機，峰值性能高達12.54億億次/s，持續性能達到9.3億億次/s。該系統部署在國家超級計算無錫中心，系統由40個運算機柜和8個網絡機柜組成。每個機柜有4個超級節點，每個超級節點包括32個節點板，每個節點板上有4個節點卡，每個節點卡有兩個節點，每個節點上裝有1個“申威26010”眾核處理器和32GB的DDR3內存，一臺機柜就有1024塊處理器。節點之間通過基于PCI-E3.0（外設部件互聯標準擴展3.0版）的神威高速網絡進行互聯。

“神威·太湖之光”系統全部采用“申威26010”眾核處理器，架構如圖1所示。

該芯片主頻為145 GHz，由4個核組組成，每個核組包含一個主核和一個從核簇，每個從核簇由64個從核組成，共260個處理器核心，同一個核組內的主核和64個從核共享主存，每個從核有獨立的、容量為64 KB的高速緩存（Local Data Memory，LDM），支持DMA（Direct Memory Access）的方式在主存和LDM之間傳輸數據。主核作為處理器核心，可以進行通信、IO、計算等操作，同一核組內的64個從核作為加速計算部件，用來加速主核代碼中計算密集部分。

圖1　“申威26010”處理器架構圖Fig. 1 Structure of the SW-26010 processor

2　OpenACC＊總述

OpenACC是由OpenACC組織（www.openaccstandard.org）于2011年推出的眾核加速編程語言，以編譯指示的方式提供眾核編程所需的語言功能，其主要目的是降低眾核編程的難度。用于C/C++和Fortran的OpenACC API和指令負責把底層的GPU任務交給編譯的同時，又提供跨系統、跨主機CPU和加速設備的支持[3]。

“神威·太湖之光”計算機系統中的OpenACC*語言是在OpenACC2.0文本的基礎上，針對申威26010眾核處理器結構特點進行適當的精簡和擴充而來的。

2.1　執行模型

OpenACC*程序的執行模型是在host（主處理器）的指導下，host和device（加速設備）協作的加速執行模型，如圖2所示。

圖2　OpenACC＊執行模型Fig. 2 Execution model of the OpenACC＊

程序首先在host上啟動運行，以一個主線程串行執行，或者通過使用OpenMP或MPI等編程接口以多個主線程并行執行，計算密集區域則在主線程的控制下作為kernel或parallel（加速構件）被加載到加速設備上執行。加速構件的執行過程包括：在設備內存上分配所需數據空間；加載構件代碼（包含相關參數）至device；構件將所需數據從主存傳輸至設備內存；等待數據傳輸完成；device進行計算并將結果傳送回主存；釋放設備上的數據空間等。大多數情況下，host可以加載一系列構件，并在加速設備上逐個執行[4]。

OpenACC*支持三級并行機制：gang、worker、vector。gang是粗粒度并行，在加速設備上可以啟動一定數量的gang。worker 是細粒度并行，每個gang內包含有一定數量的worker。vector并行是在worker內通過SIMD或向量操作的指令級并行。

在“申威26010”中，一個運算控制核心（簡稱主核）僅控制一個運算核心陣列(加速設備，簡稱從核陣列)的運行，每個運算核心陣列內有64個運算核心（簡稱從核），每個運算核心可以運行一個加速線程。默認情況下，64個加速線程被組織成64個gang、每個gang內一個worker、worker內可以vector并行的邏輯視圖。通常情況下，盡量用滿64個加速線程可以獲得較好的運行性能。

2.2　存儲模型

在CUDA和OpenCL等低層次加速器編程語言中，主機和加速器存儲器分離的概念非常明確，內存間移動數據的語句占據用戶大部分代碼。

“申威26010”中從核和主核共享內存，從核可直接訪問主存空間，并在從核內提供加速線程私有的LDM，加速計算需要存放到LDM的數據由從核控制傳輸。存儲模型如圖3所示。

圖3　OpenACC＊存儲模型Fig. 3 Memory model of the OpenACC＊

OpenACC標準中的數據管理功能將不產生實際的作用，但直接訪問主存往往帶來性能損失，因此需要充分利用從核中的LDM，提升數據訪問效率。在“神威·太湖之光”中，OpenACC*對OpenACC標準所做的主要功能延伸和語法擴展就是為了解決共享內存架構下片內高速存儲空間的使用問題。

“申威26010”中，從一個加速線程的角度來看，其可見的數據空間有三種[4]：

1）主線程數據空間：位于主存，主線程的內存空間對其創建的加速線程直接可見，加速線程可以直接訪問相應的數據；對于主線程創建的多個加速線程而言，這部分空間是共享的；程序中在加速區外定義的變量，均位于該空間內。

2）加速線程私有空間：位于主存，每個加速線程有獨立的私有空間，程序中使用private、firstprivate子句修飾的變量將存放于該空間內。

3）加速線程本地空間：位于設備內存，每個加速線程有獨立的本地空間LDM，本地空間的訪問性能是三種空間中最高的。程序中使用local、copy等數據子句修飾的變量將由編譯系統控制全部或局部存放于LDM內。主存與本地空間的數據交互由加速線程控制。

3　長波輻射過程OpenACC優化方案

GRAPES模式的輻射模型采用歐洲中期天氣預報中心（ECMWF）的長短波輻射方案，該方案將整個大氣層劃分成水平方向和垂直方向的三維網格，水平方向為橫跨地球表面的二維網格，垂直方向則表示大氣的層數[5]。數值天氣預報模式是一種非線性化離散計算模式，其計算量巨大[6]，最初的都是通過CPU計算實現，GPU出現之后，雖然可以把該算法移植到GPU上去并行執行，但編程難度大且易出錯。因此，文中設計一種基于OPENACC加速方法，以實現通過簡化的編程算法和簡潔的Fortran語言代碼完成對長波輻射過程并行處理過程。

3.1　輻射過程的構成

輻射過程是GRAPES系統中一個重要的物理過程，包含云結構的描述和輻射算法（RRTM模塊）兩部分。其中云對調解輻射平衡起著至關重要的作用[7]，本系統在云產生器的基礎上利用McICA進行輻射計算。輻射過程首先調用mcica_subcol_lw函數實現云結構的描述，之后調用RRTM模塊進行輻射計算。RRTM模塊初始化函數rrtmg_lw_init讀取數據文件中的輸入數據；然后調用rrtmg_lw_rad函數計算輻射傳輸的過程；最后調用輻射模型子函數rrtm_lw；rrtm_lw子函數一次計算一個垂直列的輻射傳輸值[6]。它主要包含以下幾個基本步驟：inatm，準備大氣剖面；cldprmc，計算云層光學深度；setcoef，計算這種大氣剖面下輻射傳輸算法所需的各種量；taumol，計算氣體光學深度和普朗克分數；rtrnmc，計算晴空和云層的輻射傳輸值。其中，taumol中的16個子函數taugb，計算了全部16個長波波域譜帶的氣態光學厚度和普朗克（Planck）函數。輻射過程的函數包結構如圖4所示。

圖4　長波輻射過程的函數結構Fig. 4 Structure diagram of the long-wave radiation function

3.2　并行優化方案

OPENACC并行化的唯一目的是利用從核充足的硬件資源來提高程序運行速度，縮短運行時間。整個長波輻射過程中最耗時間的是RRTMG_LW主過程，計算并行化的目標是將循環迭代步分散到多個不同的線程上執行，這些線程運行在多個加速器核心上，從而將計算任務由CPU轉移到加速器上，減輕CPU的負擔。

結合程序代碼的分析，確定需要眾核加速的代碼段，并根據源代碼數據結構進行預處理，對相關數組和循環進行調整，以適合OpenACC加速。優化的基本思想是確定將程序中的標量和關鍵數組盡可能多的放到局存LDM中，并設法提高程序中數據的傳輸效率，將需要使用眾核加速的循環的前后使用OpenACC加速編譯指示進行標注，其中PARALLEL表明該部分代碼是需要加速執行的并行代碼；LOOP表示下面緊跟著的循環需要在加速線程間進行并行劃分；DATA表示主從核之間的數據拷貝。

rrtm_lw函數包含一層循環，涵蓋inatm、cldprmc、setcoef、taumol和rtrnmc五個子函數調用，每個子函數內部有多個do循環分塊，最外層的iplon循環閾值為[1，90]，根據OPENACC以及申威“26010”處理器硬件特點，可以采取的并行化方案有如下幾種：

3.2.1從核并行化方案1：最外層直接并行

結合輻射過程的計算模型為單柱計算模型，核心迭代針對每個單柱串行進行多個物理過程的計算，每個單柱之間不存在數據依賴關系。可以在該層循環外層添加OpenACC加速指示語句，單個從核將獨立運行循環內部的所有函數，完成全部計算后將結果傳回主核，代碼結構示意如圖5所示，此時從核使用率為70.3%。

圖5　iplon層直接并行代碼結構示意圖Fig. 5 Structure of direct parallel methods on the iplon layer

本段核心代碼中，5個子函數中包含眾多變量，多為二維數組，少數為三維數組（標量及一維數組產生局存壓力較小）。其中輸入數組22個，包含溫度，壓力，痕量氣體濃度等多個物理量，總大小約為22×140×60(or 61)×8=1478400 byte；計算過程中的中間變量超過40個，總大小約為2688000 byte；輸出數組多為一維數組或標量。即使通過調整代碼使計算過程中的大多數中間變量空間能及時釋放，棧容量需求也遠超局存大小（65536 byte），以0.5°分辨率計算過程為例，無論是添加routine函數將部分從核函數的棧數組放到主存，還是使用分塊語句指導數據依次拷入拷出，都會產生極大的通訊開銷，影響運行速度。

另外，子函數cldprmc與setcoef只有簡單的賦值與計算過程，在原本的運行過程中耗時很少，但與前后子函數存在數據依賴，所以此方案下，同樣需要拷貝程序與數據在從核運行，運行時間大量增加。

3.2.2從核并行化方案2：循環內移

對于氣象模式，運算過程中的中間變量數組通常個數很多并且數據量極大，由于局存壓力需要多次多步拷貝，數據通信會造成較大的時間開銷影響加速效果，例如長波輻射過程的rrtmg_lw函數。

此時可以通過循環內移，對大函數進行分塊，內移工作從rrtmg_lw函數調用的函數inatm開始，依次對cldprmc、setcoef、taumol和rtrnmc子函數進行循環內移達到子函數分塊的目的，之后，對于適合從核并行的inatm、taumol和rtrnmc子函數按OPENACC規則添加導語及子句使其在從核運行，不適合從核運行的cldprmc、setcoef子函數仍在主核運行，代碼結構示意如圖6所示。此時，數據拷貝通訊量縮短至方案1的20%左右，與此同時，可根據每個子函數的最外層循環數，合并外層循環，增加從核運行指示語句，最大化擴大子函數運行的并行度。

圖6　循環內移代碼結構示意圖Fig. 6 Diagram of the loop inside displacement

循環內移使程序原本的計算順序發生了改變，由從核單線程串行計算子函數1-5變為先計算所有子函數1，保存計算結果后計算子函數2，以此類推。所以，代碼中大部分中間變量以及輸出結果需要升維用于解決數據依賴，保證正確性。

循環內移后，局存壓力顯著降低，特別是rtrnmc函數，只有6個二維輸入變量和2個二維中間變量，經過計算和試驗，只需最低12個循環分塊即可正確運行，通訊次數降低，運算速度顯著提高。同時，inatm函數作為大氣剖面準備模塊，輸入數據較多，中間變量較少，循環內移后，從核運行時所需的循環分塊個數以及數據拷貝次數也明顯降低，運行速度有所提高。

3.2.3從核并行化方案3：循環分列

為最大程度利用異構眾核處理器的結構優勢，同時進一步降低局存壓力，減少通訊時間，可以對循環內移后的子函數進行循環改寫，擴大并行度使每個從核承擔更少的計算任務，從而降低計算所需通訊量。通過分析子函數循環結構，對于迭代層數較多，耗時較長的計算熱點進行進一步從核并行化改寫。同時，通過collapse循環合并語句，可以使代碼運行時的并行度顯著增加。以inatm函數為例，改寫循環后k由子函數內層轉至函數外，為保證數據正確性，需改寫原函數，部分常量數據需要升維以符合最外層k循環，部分變量數據需要轉置以適合最高效的訪存方式，新函數名稱以inatm_trans表示，循環分列前與分列后的代碼結構示意如圖7所示，此時從核使用率為99.3%。

圖7　inatm函數循環分列前后結構示意圖Fig. 7 Diagrams before and behind the loop splitting

循環分列的目的是通過減少單個從核承擔計算任務量的方式降低數據通信量，同時增大并行度且簡化程序結構，需要在一定程度上改寫代碼結構使其重組，合并，擴大最外層循環大小。整個rrtmg_lw函數改寫合并循環后，同時并行的線程增大至5400個，每個從核的輸入數據量為22×140×8=24640 byte，中間變量大小為7×140×8=7840 byte，總數據量已經低于局存大小，不再需要數據分塊拷貝，重新調整其通訊策略后，運行速度得到進一步提升。

對于某些原本不適合進行從核并行的計算核心段需要循環改寫后再進行從核并行，否則只能通過循環內移分塊跳過，使其在主核完成計算過程，影響模式運行效率。例如計算氣態光學厚度和普朗克量的taumol函數包含16個子函數對應lay變量的16種分割方式，不具備從核并行條件，通過合并成為一個最外層以lay為變量的大循環后，順利在從核運行且使用率超過90%。

在實際情況下，方案1，2，3實現數據并行計算的不同思路通常需要結合使用。針對不同核心函數不同的計算結構與數據規模，分析設計出各自最高效的通訊策略。本次工作中主要選取了長波輻射部分的2個主函數和7個子函數作為核心段，最終實現 “粗粒度MPI并行+細粒度眾核并行”多級異構并行方案。

3.2.4其他并行優化措施

使用ldm數學函數庫：在OpenACC加速指導語句的基礎上，使用ldm數學函數庫，經測試，ldm數學函數精度與默認庫保持一致，計算性能提升30%～80%不等，進一步提升了GRAPES系統整體運算效率。

3.2.5結果與分析

試驗采用分辨率為0.5°×0.5°的GRAPES全球模式長波輻射方案，積分步長為36，優化選項使用-O2。由于cldprmc、setcoef耗時較少，沒有進行優化。加速后，對比純主核運行輸出的全球長波輻射通量與模式輸出值，經驗證誤差，相對誤差的量級在10－5至l0－5之間，在允許范圍內。優化結果對比如表1所示。

表1　眾核加速結果Table 1 The result of many-core acceleration

4　展望與不足

隨著編譯器的進一步優化和硬件技術的發展，OpenACC加速與CUDA等在底層技術實現上的差距將會越來越小，而且支持OpenACC加速指令轉換將不僅僅針對CUDA設備，還包括其他更多廠商的硬件加速設備，從而能極大地提高OpenACC加速指令的普適性，增大程序可移植性。

文中以GRAPES全球模式RRTMG輻射模塊長波過程為加速對象，目前只應用到為數不多的熱點函數，就能得到較好的加速效果，優化還有很大的提升空間。

下一步將繼續深入研究源代碼數據結構，加強數據預處理，減少從核與主存之間的通訊頻率，充分利用從核中高速緩存LDM，提升數據訪問效率；同時嘗試多種OpenACC加速指令組和使用，進一步提升加速效果。

[1]伍湘君, 金之雁, 黃麗萍, 等．GRAPES模式軟件框架與實現. 應用氣象學報, 2005, 16(4): 539-546．

[2]陳德輝, 沈學順．新一代數值預報系統GRAPES研究進展. 應用氣象學報, 2006, 17(6): 773-777．

[3]曾文權, 胡玉貴, 何擁軍, 等．一種基于OPENACC的GPU加速實現高斯模糊算法. 計算機技術與發展, 2013, 23(7): 147-150．

[4]何滄平. OpenACC并行編程實戰. 北京: 機械工業出版社, 2016．

[5]郭妙, 金之雁, 周斌. 基于通用圖形處理器的GRAPES長波輻射并行方案. 應用氣象學報, 2012, 23(3): 348-354．

[6]鄭芳, 許先斌, 向冬冬, 等. 基于GPU的GRAPES數值預報系統中RRTM模塊的并行化研究. 計算機科學, 2012, 39(6A): 370-374.

[7]荊現文, 張華. McICA 云—輻射方案在國家氣候中心全球氣候模式中的應用與評估. 大氣科學, 2012, 36(5): 945-958.