二維穩(wěn)態(tài)熱傳導(dǎo)CUDA 并行

王 民

（煙臺黃金職業(yè)學(xué)院，山東 煙臺 265401）

0 引 言

有限差分法用于解決難以或無法分析解決的微分方程，其中包括拉普拉斯方程。求解拉普拉斯方程通常是為了找到平衡解，因此這是一個可以從物理問題中產(chǎn)生的方程。借助于有限差分實現(xiàn)二維拉普拉斯方程[1-3]在穩(wěn)態(tài)熱傳導(dǎo)中的應(yīng)用有助于解決其溫度方面的問題。本文在基于有限差分的拉普拉斯方程程序編寫基礎(chǔ)之上，對其進(jìn)行CUDA 并行化[4-5]，提升數(shù)據(jù)計算的效率，從而得到更好的擴(kuò)展性。

1 二維穩(wěn)態(tài)熱傳導(dǎo)應(yīng)用程序

1.1 二維穩(wěn)態(tài)熱傳導(dǎo)理論

穩(wěn)態(tài)熱傳導(dǎo)分析可以用來分析恒定中心熱負(fù)荷對系統(tǒng)或部件的影響[6-7]。通常情況下，在進(jìn)行臨時熱分析之前，要進(jìn)行穩(wěn)態(tài)熱分析以確定初始溫度分布。穩(wěn)態(tài)熱分析可用于分析溫度等參數(shù)，這些參數(shù)是由恒定熱負(fù)荷引起的，通過有限元素的計算，將連續(xù)的傳熱問題離散化為一系列離散的代數(shù)方程組，然后利用計算機(jī)進(jìn)行求解。

1.2 有限差分的拉普拉斯

1.2.1 有限差分

有限差分法是一種數(shù)值方法，用于解決偏微分方程和積分方程。在有限差分法中，微分方程中的導(dǎo)數(shù)是用有限差分公式近似計算的[8-10]，可以將[a，b]的區(qū)間劃分為n個長度為h的相等子區(qū)間，如圖1 所示。

圖1 有限差分法的區(qū)間劃分

有限差分法基本上是一種近似導(dǎo)數(shù)的數(shù)值方法，所以先要分析如何取導(dǎo)數(shù)。一個函數(shù)f(x)的導(dǎo)數(shù)的定義如下：

通常情況下，在有限差分方法中使用中心差分公式，因為它們能產(chǎn)生更好的精度。微分方程只在網(wǎng)格點上執(zhí)行，而第一和第二導(dǎo)數(shù)分別為：

這些有限差分表達(dá)式被用來取代微分方程中y的導(dǎo)數(shù)，如果微分方程是線性的，就會產(chǎn)生一個n+1 的線性代數(shù)方程系統(tǒng)。如果微分方程是非線性的，代數(shù)方程也將是非線性的。

在數(shù)值計算過程中，將解決方案區(qū)域劃分為不同的網(wǎng)格，最終網(wǎng)格點用于替換不同網(wǎng)格中的連續(xù)解決方案區(qū)域。將要求解的流動變量存儲在每個網(wǎng)格點，并用相應(yīng)的微分商替換偏微分方程的微分元素，從而將偏微分方程轉(zhuǎn)換為代數(shù)微分方程，獲得離散點處具有有限未知變量的微分方程。通過求解微分方程，獲得了網(wǎng)格點處流動變量的數(shù)值解。

1.2.2 二維的拉普拉斯方程

拉普拉斯方程是一個二階偏微分方程，出現(xiàn)在許多科學(xué)和工程領(lǐng)域，如電力、流體流動和穩(wěn)定熱傳導(dǎo)等。在一個領(lǐng)域中解決該方程，需要指定某些條件，即未知函數(shù)在域的邊界必須滿足的條件。在空間上二維拉普拉斯方程滿足如下等式：

對x、y方向進(jìn)行有限差分，并基于4 點格式執(zhí)行計算。考慮到一塊薄板的頂部和底部是完全絕緣的，它的每個邊緣都有對應(yīng)的已知溫度。因此，目標(biāo)是找到穩(wěn)態(tài)溫度的分布。區(qū)域內(nèi)部的溫度取決于它周圍的溫度。我們可以將該區(qū)域劃分為細(xì)網(wǎng)狀的點h(i,j)。一個內(nèi)部點的溫度可以被看作是四個相鄰點溫度的平均值，如圖2 所示。

圖2 四點差分區(qū)域劃分

對于這種計算，采用與內(nèi)部點相鄰的點來描述邊緣是很方便的。點h(i,j)的內(nèi)部點是指0＜i＜n，0＜j＜n范圍內(nèi)的點，并且有(n-1)×(n-1)個點。邊緣點是當(dāng)i=0，i=n，j=0或j=n時，有對應(yīng)于固定溫度值的網(wǎng)格邊緣點。因此，h(i,j)的全部范圍是0 ≤x≤n，0 ≤y≤n，并且有(n+1)×(n+1)個點，可以通過迭代來計算每個點的溫度方程。

1.2.3 串行和CUDA 并行化設(shè)計

（1）串行程序算法

假設(shè)每個點的溫度被保存在一個數(shù)組h[i][j]中，邊界點h[0][x]、h[x][0]、h[n][x]、h[x][n]（0 ≤x≤n）已經(jīng)被初始化為邊緣溫度。串行核心代碼如圖3 所示。

圖3 串行核心代碼

在計算過程中使用固定次數(shù)的迭代。請注意，圖3 中第二個數(shù)組g[][]是用來保存從舊值中計算出的點的新值；數(shù)組h[][]被更新為g[][]中的新值；乘以0.25 而不是除以4 是為了使計算更有效率，因為乘法通常比除法更有效。提高順序代碼效率的方法可以延續(xù)到GPU 代碼，并且在所有情況下都應(yīng)該這樣做。

將所要計算的薄板的兩側(cè)溫度保持為20 ℃，其中一側(cè)有較短的一段保持在100 ℃，如圖4 所示。

圖4 初始化薄板方法

（2）CUDA 程序算法

采用非共享方式實現(xiàn)對程序的并行。分別對x、y方向進(jìn)行分塊計算，每個塊所啟動的線程組數(shù)為x=(n+block.x-1)/block.x,y=(n+block.y-1)/block.y。

將核心計算部分封裝到global 函數(shù)中，將每次的迭代計算結(jié)果拷貝回從設(shè)備重新參與計算工作。完成整個迭代后計算結(jié)果將被拷貝到主存中。

在從設(shè)備上執(zhí)行計算，基于線程以及設(shè)備號獲取每個數(shù)據(jù)塊的坐標(biāo)點。

2 正確性驗證以及性能分析

2.1 正確性驗證

對于程序正確性的驗證，主要通過比較串行數(shù)據(jù)結(jié)果與并行數(shù)據(jù)結(jié)果來完成，程序上將計算結(jié)果保存到csv 文件中，借助于編寫Python 腳本來讀取兩個文件的數(shù)值，對比每個數(shù)值點的值，當(dāng)誤差值在10 ～14 以內(nèi)時，說明計算結(jié)果正確。通過對點的比較完成程序正確性驗證。

2.2 性能分析

在測試中，本文選擇的數(shù)組長度的計算規(guī)模分別是1 024、2 048，所使用的CUDA 方向上塊的大小分別是（1，1）、（2，2）、（4，4）、（8，8）。迭代次數(shù)為固定值100。

通過測試記錄串行時間以及并行時間，然后根據(jù)記錄時間計算出程序的加速比S和并行效率E，通過這兩個值分析程序的性能，如圖5 所示。

從圖5 中可以看出，在不同數(shù)據(jù)規(guī)模下程序的加速比和并行效率都得到改善。在相同規(guī)模下，隨著分塊數(shù)的增加，程序的加速比自增，且其并行效率也符合了小于1 的要求。在數(shù)據(jù)規(guī)模分別為1 024 和2 048 兩種情況下，最大加速比值將是最小加速比值的N倍。

圖5 加速比和并行效率

3 結(jié) 語

本文采用CUDA 實現(xiàn)了二維穩(wěn)態(tài)熱傳導(dǎo)程序的并行化。分析測試數(shù)據(jù)得出，在執(zhí)行并行化后，隨著分塊數(shù)的增加，其程序性能要比原來提升幾倍。這對其實際應(yīng)用有很大的幫助。