第三章 數值方法
3.7 CUDA高速運算
由於此研究三維低速可壓縮流之計算量非常龐大，因此在模擬過程相當耗時。
過去的研究中，多使用平行化方法做為提高計算效能的方式，最常見的兩種平行 化方法為MPI (Message Passing Interface)及OpenMP (Open Multi-Processing)。前者 將程式平行化於並聯的多台電腦上運行，但效率隨電腦數增加而下降，且價格與 佔用空間亦十分龐大。後者則是使用具有多核心中央處理器的超級電腦，佔用空 間較小，但受限於目前的製程及散熱問題，計算核心數有限，效率提昇有限。
近年來，處理動畫(Animation)功能在研究與娛樂的殷切需求下，顯示卡 (Video Card)內圖形處理器(Graphics processing unit, GPU)的速度越來越快，如圖 3-3所示。在浮點運算數及記憶體頻寬上，GPU由於本身硬體架構的優勢，整體 的運算效能遠高於一般的多核心中央處理器(Central processing unit, CPU)。除此 之外，顯示卡之價格相對於購買高階多核心中央處理器較為低廉，因此於三維圖 形 顯 示 方 面 以 外 的 應 用 日 益 增 加 。 CUDA 平 台 （ Compute Unified Device Architecture）是由Nvidia公司所推出的整合技術，此技術是一種將GPU作為數據 平行化運算設備的統一處理架構。利用顯示卡內GPU具有數量眾多的流處理器 (Stream Processor，SP)，更適合用於平行化程式的計算特性，以達到提升數值計 算速度的目的。本研究發展利用CUDA於計算流體力學上，將CUDA平台成功運 用於平行化運算，以下將就CUDA平台的架構及執行模式做詳細的說明。
CUDA 平台的主要架構如圖 3-4(a)，右方可見主機(Host)端，內含中央處理 器(CPU)和主機記憶體。而設備(Device)端，即主機板插槽上的顯示卡。平行計 算過程中，依程式語言撰寫順序執行，部分在 CPU 上運算的為串行處理程式 (serial code)，部分藉自行開發在 CUDA 平台上的連結編輯程式，在 GPU 上執行 的為核心函數(kernel)。每個 kernel 的執行單位為網格(Grid)，網格中包含數個區 塊(Block)，區塊又由數個執行緒(Thread)所組成，當 kernel 開始執行時，這些執 行緒將分別由顯示卡內的流處理器(Stream Processor, SP)來作平行化運算。待此
48
kernel 結束，則接續下一個串行函數進行在 CPU 運算，依此順序執行至整個程 式完成運算工作。
在 CUDA 平台上作平行化運算，為得到最佳化的計算效能，在記憶體的配 置方式上極為重要，如圖 3-4(b)所示。由於在 GPU 的計算單元中，包含數個多 流處理器(Stream Multiprocessor，SM)，而每個 SM 中又包含 8 個單(或雙)精確浮 點數的流處理器(Stream Processor，SP)，在同 SM 中的 8 個 SP 共用同一共享記 憶體(shared memory)。如前述在 CUDA 平台中的 kernel 實質上是以區塊為單位 執行的，在同一區塊的所有執行緒需共享數據，因此會在同一個 SM 中執行，但 同一 SM 中可以存在多個區塊，如此可在某一區塊執行同步或存取記憶體時，讓 另一個區塊占用執行資源，用以隱藏延遲，更好的利用執行單元的資源。然而在 實際運行時，區塊內的執行緒將以執行緒束(warp)為單位執行，一個 warp 由 32 個連續執行緒所組成，且一個 SM 可以保存 32 個 warp(共 1024 個執行緒)。由 於各個 warp 在 CUDA 中的執行順序並無規律，因此需使用同步的方式，確保同 一區塊中的執行進度相同，一個 SM 可以完成 512 個執行緒同步的動作，且 warp 中的執行緒只與執行緒 ID 有關，排列如 ID 0~31 為一束，32~63 為第二束，以 此類推，因而在下述設定的區塊中執性緒數則可是 1~512 的任意值，並以 32 的 倍數為最佳，以避免多餘的執行緒仍占掉一個 warp。依此原則對區塊及執行緒 數作配置，以期達到平行化運算的最高性能。
在撰寫CUDA程式時，主要可分為下列三個步驟；第一步驟為配置顯示卡內 GPU的記憶體位址及設定區塊(Block)和執行緒(Thread)的大小，接著將本來利用 中央處理器CPU計算的資料由主機記憶體搬移至顯示卡上的記憶體。一般使用有 限差分法於三維的情形中，多使用三維陣列(i,j,k)將格點編號，然而若使用三維 陣列時，要將主機的資料搬移至顯示卡上的記憶體時，將非常複雜。因此需將三 維陣列轉換成一維陣列以利於資料搬移。其轉換公式如下所示：
Device ( assign numbers ) = Host ( i ny nz
× × + × +
j nz k )，其中nx、ny與nz分別49
為x、y與z方向的計算網格數目。
第二步驟為將程式平行運算。在步驟一將三維陣列轉換成一維陣列，且把主 機端的資料搬移至顯示卡的記憶體後，接著將一維陣列的迴圈平行化，利用流處 理器做計算並儲存於顯示卡上的記憶體。第三步驟為將 GPU 計算結果從顯示卡 的記憶體搬移至主機端的記憶體，接著將前述的一維陣列轉回原來的三維陣列，
並輸出結果。
為了使程式能編譯一次後即能在擁有不同流處理器(Stream Processor，SP)數 的顯示卡上執行，CUDA 平台使用了單指令多執行緒(Single Instruction, Multiple Thread，SMIT)的執行模型，將計算任務轉成可以大量平行執行的執行緒。實際 運行時，kernel 是以區塊為單位執行的，而網格則是用以表示一系列可以被平行 執行的區塊的集合。一個 kernel 會在一個網格中運行，但由於同一區塊中的所有 執行緒在同一時刻執行的指令並不一定相同，故在 CUDA 平台中，同一區塊中 的所有執行緒皆能存取共享記憶體(shared memory)，因此可以互通交換數據並能 快速進行同步的動作，不同區塊中的執行緒則否，在此模式的設定下雖然能同時 執行同一指令的執行緒數量有限，但通過共享記憶體和執行緒柵欄同步(Barrier) 的方式，確保同一區塊內的執行緒不僅能獨立平行運算，在需要資料同步時，亦 能待全部執行緒皆執行到相同指令時，程式才繼續往下運行，而不致計算錯誤。
此外，在記憶體最佳化的配置下，由於不同種類的記憶體將對應到不同階層 的運算單元，因此記憶體最佳化的程度也是影響 GPU 運算的重要因子。在最低 階的運算單元，執行緒中，每個執行緒有其可存取的 local memory 及 registers 記 憶體，為專屬於此執行緒的記憶單元，速度也最快。而在同一個區塊中，各個執 行緒將可存取一共用的共享記憶體，但相對的速度也較慢。最後則是在同一個網 格中，所有的執行緒皆可存取共用的 global memory，constant memory 及 texture memory，速度則最慢。此外，由於只有 global memory，constant memory 及 texture memory 可與主機端的記憶體進行互相存取的動作，因此善用各種不同的記憶體
50
在此為了解決此問題，以 Candler 和 Wright[38]所提出的 data-parallel LU relaxation method 來修改原始的 LUSGS 法，進一步達到全平行化的程式運算。
為了改寫 LUSGS 法，首先以一預設的R(ki,j,k)得到第零次運算的∆U(p0,()i−1,j,k)，如以 接著結合式(3-87)及(3-90)並改寫如下：
}
m為次疊代(subiterations)的次數，根據文獻[38]的測試，以mmax =4之效率為最 佳。因此最後可得：
51
表 3-1 精度係數值
β θc θd Order
1/3 0 0 2
1/3 -1/6 0 3
1/3 0 -1/6 4
1/3 -1/10 -1/15 5
1/3 -1/10 -1/15 6
52
圖 3 - 1 黎曼問題特徵值結構圖 y
x< 0 x= 0 x> x 0 λp
λi
λ2
λ1
1
λp+
λi
1
λm−
λm
53
(a)
(b)
圖 3 - 2 L 、1 L 、2 L 、3 L 與4 L 於管道兩端的方向性示意圖 5
1 2 3 4 5
L L L L L
0 u >
aperture
y
x
channel
surroundings fluid velocity →
1 2 3 4 5
L L L L L
0 u <
aperture
y
x channel
surroundings
fluid velocity ←
54
(a)
(b)
圖 3 - 3 GPU 及 CPU 效能比較[39](a)浮點運算數 (b)記憶體頻寬
55
(a)
(b)
圖 3 - 4(a)CUDA 程式架構 (b)CUDA 記憶體配置[39]
56