王 民

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

0 引 言

有限差分法用于解決難以或無(wú)法分析解決的微分方程，其中包括拉普拉斯方程。求解拉普拉斯方程通常是為了找到平衡解，因此這是一個(gè)可以從物理問(wèn)題中產(chǎn)生的方程。借助于有限差分實(shí)現(xiàn)二維拉普拉斯方程[1-3]在穩(wěn)態(tài)熱傳導(dǎo)中的應(yīng)用有助于解決其溫度方面的問(wèn)題。本文在基于有限差分的拉普拉斯方程程序編寫(xiě)基礎(chǔ)之上，對(duì)其進(jìn)行CUDA 并行化[4-5]，提升數(shù)據(jù)計(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)分析可以用來(lái)分析恒定中心熱負(fù)荷對(duì)系統(tǒng)或部件的影響[6-7]。通常情況下，在進(jìn)行臨時(shí)熱分析之前，要進(jìn)行穩(wěn)態(tài)熱分析以確定初始溫度分布。穩(wěn)態(tài)熱分析可用于分析溫度等參數(shù)，這些參數(shù)是由恒定熱負(fù)荷引起的，通過(guò)有限元素的計(jì)算，將連續(xù)的傳熱問(wèn)題離散化為一系列離散的代數(shù)方程組，然后利用計(jì)算機(jī)進(jìn)行求解。

1.2 有限差分的拉普拉斯

1.2.1 有限差分

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

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

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

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

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

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

1.2.2 二維的拉普拉斯方程

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

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

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

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

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

（1）串行程序算法

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

圖3 串行核心代碼

在計(jì)算過(guò)程中使用固定次數(shù)的迭代。請(qǐng)注意，圖3 中第二個(gè)數(shù)組g[][]是用來(lái)保存從舊值中計(jì)算出的點(diǎn)的新值；數(shù)組h[][]被更新為g[][]中的新值；乘以0.25 而不是除以4 是為了使計(jì)算更有效率，因?yàn)槌朔ㄍǔ１瘸ǜ行?。提高順序代碼效率的方法可以延續(xù)到GPU 代碼，并且在所有情況下都應(yīng)該這樣做。

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

圖4 初始化薄板方法

（2）CUDA 程序算法

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

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

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

2 正確性驗(yàn)證以及性能分析

2.1 正確性驗(yàn)證

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

2.2 性能分析

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

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

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

圖5 加速比和并行效率

3 結(jié) 語(yǔ)

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