熱傳導(dǎo)方程的一種新的具有并行本性的區(qū)域分解算法及其并行程序設(shè)計(jì)

1、熱傳導(dǎo)方程的一種新的具有并行本性的區(qū)域分解算法及其并行程序設(shè)計(jì)1左風(fēng)麗*1 袁光偉*2北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所應(yīng)用科學(xué)計(jì)算部 高性能計(jì)算中心，北京 8009信箱 10088ABSTRACT A domain decomposition method based on Du Fort-Frankel scheme and the fully implicit scheme for heat conduction is presented and some characters are proved by numerical experiments in this paper. Two kin

2、ds of parallel strategies are designed and some programming techniques are discussed, which is less memory storage than conventional programming method leading to obtaining good result. By performance analysis of parallel numerical experiments, the method with intrinsic parallelism has the same accu

3、rate as the fully implicit scheme, which suit large scale scientific computing and engineering computing.KEYWORDS Heat conduction; Domain decomposition method; High accuracy; Highparallel efficiency1. 引言許多物理和力學(xué)問題的求解都可歸結(jié)為熱傳導(dǎo)方程的求解，其數(shù)值模擬常采用有限差分方法，在實(shí)際應(yīng)用問題中通常采用隱式差分格式，這是因?yàn)殡[式差分格式具有無條件穩(wěn)定的特征，對任意的網(wǎng)格比都能保證對計(jì)算誤差

4、的控制。隨著超大規(guī)模并行計(jì)算機(jī)的發(fā)展，對所要求解的問題的精度要求也越來越高，因而求解問題的規(guī)模變得越來越大。盡管單機(jī)的計(jì)算速度和內(nèi)存容量變得很大，但是仍難以完全滿足需求，其有效的解決方法就是構(gòu)造良好的具有并行本性的區(qū)域分解方法。本文討論求解一維熱傳導(dǎo)方程的區(qū)域分解方法，將求解的區(qū)域剖分成多個(gè)子區(qū)域，在子區(qū)域的內(nèi)邊界點(diǎn)使用Du FortFrankel格式，在子區(qū)域的內(nèi)點(diǎn)使用全隱式格式。它不但具有并行本性的特征，而且具有與純?nèi)[式格式一樣的精度和穩(wěn)定性。在構(gòu)造并行算法的具體實(shí)現(xiàn)時(shí)，提出一種動(dòng)態(tài)局部數(shù)據(jù)管理的方法，這樣做的好處不僅可以提高數(shù)據(jù)的局部性，而且可以任意改變內(nèi)存的大小，而不改變通信的次數(shù)和

5、通信量的大小，因此可以獲得良好的結(jié)果。我們知道，對熱傳導(dǎo)方程的數(shù)值求解已有很多的方法，通常采用的格式為隱式格式。隱式格式之所以被廣泛的應(yīng)用是因?yàn)樗菬o條件穩(wěn)定和絕對收斂的，但是，求解隱式格式需要1 本課題得到國家973項(xiàng)目（G1999032801），國家自然科學(xué)基金（19932010）和中物院重點(diǎn)基金資助。 左風(fēng)麗，碩士，助理研究員，大規(guī)?？茖W(xué)與工程計(jì)算并行研究方向。*2袁光偉，博士，研究員，大規(guī)模科學(xué)與工程計(jì)算并行研究方向。 *1解一個(gè)代數(shù)方程組，并且它所求解的方程組的階數(shù)與劃分節(jié)點(diǎn)個(gè)數(shù)是一樣的，這樣，隨著計(jì)算節(jié)點(diǎn)的增加計(jì)算量變得越來越大，甚至難以計(jì)算下去。由于對求解問題的要求越來越高，網(wǎng)格

6、節(jié)點(diǎn)不斷加密，即使有高性能大存儲(chǔ)量的并行機(jī)，存儲(chǔ)空間也難以完全滿足需求。而基于MPI編程的區(qū)域分解技術(shù)是減少單機(jī)上計(jì)算量和存儲(chǔ)量的一條有效途徑。例如，周毓鱗沈隆鈞、袁光偉、張寶琳53 41,2、等人對具有并行本性的差分算法已開展多年的研究，諸如拋物型方程的實(shí)用的隱式區(qū)域分解方法和交替分段顯隱方法等。Kelly Black的利用配置法的區(qū)域分解方法是：在時(shí)間上采用有限差分，在多個(gè)空間子區(qū)域上關(guān)于x的二階偏導(dǎo)數(shù)利用正交多相式配置法進(jìn)行近似，在子區(qū)域的內(nèi)邊界點(diǎn)利用懲罰方法并采用Du Fort-Frankel格式，在子區(qū)域的內(nèi)點(diǎn)采用三點(diǎn)中心差分全隱式格式。文5中給出數(shù)值實(shí)驗(yàn)的節(jié)點(diǎn)數(shù)最多為256，劃分的

7、子區(qū)域數(shù)最多為32個(gè)。并且，劃分子區(qū)域數(shù)越多，計(jì)算的步數(shù)越少，在劃分的子區(qū)域數(shù)為32時(shí)，計(jì)算步數(shù)多數(shù)是幾十步。本文基于對以上方法的研究，從而提出一種新的高精度的區(qū)域分解方法，在時(shí)間上采用有限差分，在空間上，對求解區(qū)域進(jìn)行剖分，從而形成多個(gè)子區(qū)域，在子區(qū)域的內(nèi)邊界點(diǎn)上關(guān)于x的二階偏導(dǎo)數(shù)利用Du Fort-Frankel格式，在子區(qū)域的內(nèi)點(diǎn)上采用三點(diǎn)中心差分全隱式格式，并對此方法進(jìn)行并行化實(shí)現(xiàn)和分析。我們對所構(gòu)造的格式給出滿足無條件穩(wěn)定性必要條件的證明；在進(jìn)行數(shù)值模擬時(shí)，可以不考慮機(jī)器的存儲(chǔ)量，可以任給節(jié)點(diǎn)個(gè)數(shù)（數(shù)十萬），劃分的子區(qū)域數(shù)（數(shù)千個(gè)）也可以是任意的；我們可以不管劃分多少個(gè)子區(qū)域，計(jì)算（

8、迭代）的步數(shù)均可以達(dá)到上萬步，雖然格式帶來的累積誤差對計(jì)算精度有一定的影響，但考查的絕對誤差仍然可以達(dá)到10，這一結(jié)果是可以接受的；我們通過數(shù)值實(shí)驗(yàn)考察了計(jì)算格式的穩(wěn)定性、收斂性和可擴(kuò)展性等特性。本文是這樣組織的，在第二節(jié)，構(gòu)造新的具有并行本性的區(qū)域分解方法，這種方法可描述為：在時(shí)間方向上進(jìn)行有限差分離散，在空間方向上，對求解區(qū)域進(jìn)行剖分，從而形成多個(gè)子區(qū)域，在子區(qū)域的內(nèi)邊界點(diǎn)使用三層顯式 Du Fort Frankel 格式，在子區(qū)域的內(nèi)點(diǎn)使用三點(diǎn)中心差分純隱式格式；第三節(jié)，我們給出了該方法滿足無條件穩(wěn)定性的必要條件的證明；第四節(jié)，我們給出具有如下特點(diǎn)的并行算法：（）計(jì)算的總節(jié)點(diǎn)個(gè)數(shù)（數(shù)十萬

9、甚至更大）可以大規(guī)模的（即，并行算法具有良好的可擴(kuò)展性）；（）劃分的子區(qū)域包含的節(jié)點(diǎn)數(shù)也是任意選取的；（）我們采用動(dòng)態(tài)局部數(shù)據(jù)管理的方法：任意單機(jī)的存儲(chǔ)量僅僅是對任意一個(gè)子區(qū)域生成的方程組的系數(shù)的存儲(chǔ)，當(dāng)計(jì)算下一個(gè)子區(qū)域的值之前，重新生成方程組的系數(shù)矩陣，這樣，子區(qū)域的方程組的系數(shù)矩陣是一個(gè)臨時(shí)數(shù)組，每次計(jì)算子區(qū)域的值都要更新一次系數(shù)矩陣臨時(shí)數(shù)組的值）。第五節(jié)，我們給出數(shù)值結(jié)果的性能分析，包括求解精度的分析，并給出穩(wěn)定性、收斂性和可擴(kuò)展性的數(shù)值驗(yàn)證。6-42 具有并行本性差分格式的構(gòu)造我們考察如下的一維熱傳導(dǎo)方程：Ut=Uxx,(x,t)QT （1）U(0,t)=U(l,t)=0,t0,T （

10、2） U(x,0)=U0(x),x0,l （3） 其中，QT=(0,l)(0,T),T0,U=U0(x)是初始條件，滿足U0(0)=U0(l)=0?？臻g步長為h，時(shí)間步長，用平行線x=xj=jh(j=0,1,.,J)，t=t=n(n=0,1,.,N)劃分區(qū)域QT，其中，Jh=l,N=T,J和N是正整數(shù)。不失一般性，我們考慮將0,x分成n 2兩個(gè)子區(qū)域0,x1和0,x2的情形。分成多個(gè)子區(qū)域并沒有實(shí)質(zhì)上的區(qū)別。令k滿足1kJ1的正整數(shù)。對問題（1）（3）建立如下的具有并行本性的區(qū)域分解格式：在子區(qū)域的內(nèi)邊界x=xk處使用Du Fort-Frankel三層格式+11ununjjn+11nun+un

11、j+1(ujj)+uj12=h2 （4）在子區(qū)域的內(nèi)部使用純隱式格式 +1ununjj+1n+1+1un+unj+12ujj1 n+1=h2 （5） 離散初邊界條件為： U0 =UJn+1=0， 0nN （6）n+1u0 0jJ （7） j=U0(xj)，由格式（4）可首先顯式地計(jì)算出Uk。當(dāng)Uk計(jì)算出后，可由（5）、（6）分別并行求出+1+1un(1jk)和un(kjJ1)。 jjn+1而格式（4）的截?cái)嗾`差為O(+h+ 眾所周知，格式（5）的截?cái)嗾`差為O(+h)，222h2)。如果將空間區(qū)域僅分為有限多個(gè)子區(qū)域，即內(nèi)邊界點(diǎn)的個(gè)數(shù)為O(1)（該數(shù)為與h無關(guān)的任一固定數(shù)），那么，我們可以期望所

12、得到的解的收斂階為O(+h+/h)，數(shù)值實(shí)驗(yàn)也將證實(shí)這一點(diǎn)，數(shù)值實(shí)驗(yàn)表明，這里構(gòu)造的格式（4）、（5）的精度比純Du Fort-Frankel的精度要高。由于Du Fort-Frankel格式和純隱式格式均為無條件穩(wěn)定的，所以我們可以期望關(guān)于格式（4）（7）也是無條件穩(wěn)定的。我們可以用矩陣的方法證明，此方法滿足無條件穩(wěn)定的必要條件。223 并行算法設(shè)計(jì)因?yàn)樵趦?nèi)邊界點(diǎn)采用三層Du Fort-Frankel格式，所以，第一時(shí)間層uj,1jJ1 的計(jì)算需采用其他格式。首先，我們設(shè)計(jì)如下的并行方法：+11ununjjn+11nun+unj+1(ujj)+uj112+1ununjj=h2，在子區(qū)域的內(nèi)邊

13、界點(diǎn)，n=1,2,.,N （8）=+1n+1+1un+unj+12ujj1h2，在子區(qū)域的內(nèi)點(diǎn)，n=0,1,.,N （9）以及如下的純顯式格式+1ununjjnnunj+12uj+uj11=h2，在子區(qū)域的內(nèi)邊界點(diǎn)，n=0 （10）注1：uj,1jJ1的計(jì)算也可完全采用純隱式格式。經(jīng)數(shù)值實(shí)驗(yàn)知，這兩種方法的精度是非常接近的，所以下文討論根據(jù)格式（8）（10）。注2：本文討論的問題是針對一維情形，所以，這里所指的子區(qū)域是由某條子線段上的離散點(diǎn)組成，在進(jìn)行并行算法設(shè)計(jì)時(shí)，開辟的存儲(chǔ)空間大小僅為子區(qū)域的大小，對任一個(gè)子區(qū)域所求解方程組的系數(shù)都要?jiǎng)討B(tài)生成。因此，不管求解的節(jié)點(diǎn)數(shù)有多大，當(dāng)劃分的子區(qū)域足

14、夠多的情形下，存儲(chǔ)量可以任意小，這與一般的MPI編程時(shí)對區(qū)域的劃分（每個(gè)處理機(jī)的存儲(chǔ)量至少是總存儲(chǔ)量的1/P，P為處理機(jī)的個(gè)數(shù)）是不同的。4 數(shù)值實(shí)驗(yàn)及結(jié)果分析41 數(shù)值實(shí)驗(yàn)環(huán)境我們對所求解的偏微分方程的處邊值問題（1）、（2）和（3）進(jìn)行數(shù)值實(shí)驗(yàn)，其中，l的取值為，初始值為U0(x)=sin(x)，邊界值為U(0,t)=U(,t)=0，其精確解為U(x,t)=etsinx。根據(jù)算法二設(shè)計(jì)并行程序。將劃分的多個(gè)子區(qū)域（子區(qū)域數(shù)大于處理機(jī)的個(gè)數(shù)）按子區(qū)域的順序均分到多個(gè)處理機(jī)上。測試環(huán)境：聯(lián)想奔和某一并行機(jī)。編程環(huán)境：MPI （Message Passing Interface）。42 精度分析

15、 我們從兩個(gè)方面考察格式數(shù)值解的精度，一方面，不同網(wǎng)格比對精度的影響，另一方面，不同分段數(shù)對精度的影響。首先，我們來看不同網(wǎng)格比對精度的影響，節(jié)點(diǎn)數(shù) 127999子區(qū)域的隱式點(diǎn)數(shù)99模型1 計(jì)算步數(shù) 1000分段數(shù) 12802表1(某并行機(jī)) Error: 絕對誤差 Aerror :Error/(+h)Error Aerror 時(shí)間步長 網(wǎng)格比（/h2）1E-84E-8 1E-6 1E-516.6 66.4 1660.05 16600.54.629E-13 1.381E-12 1.585E-8 1.372E-54.366E-5 3.402E-5 1.584E-2 1.3717我們可從表1的結(jié)果

16、可以看出，隨著網(wǎng)格比的增加，精度有所下降，但是，其最大絕對 誤差仍控制在1E-5，精度仍能滿足要求。其次，我們看不同分段數(shù)對精度的影響，計(jì)算的模型如下：模型2節(jié)點(diǎn)個(gè)數(shù) 1279計(jì)算時(shí)間步數(shù)1000分段數(shù)（Nseg）Error Aerror我們可從表2的結(jié)果可以看出，隨著分段數(shù)的增加，其數(shù)值解的精度隨之下降。其 原因是在內(nèi)邊界點(diǎn)Du Fort-Frankel格式的引入，我們知道，Du FortFrankel格式與微 分方程相容的充分必要條件是與h的比值趨向于0，反之，如果/h是一常數(shù)，此時(shí)， 差分格式就不再與方程（1）相容，而與雙曲方程網(wǎng)格比 1664 1.742E-5 1.731E2表2 8

17、2.036E-4 0.202416 5.843E-4 0.5808 32 1.351E-3 1.34272uu2u2+2=0 txt相容，如圖1表示引入Du Fort-Frankel格式的誤差曲線圖。所以，我們在選取空 間步長和時(shí)間步長時(shí)，不僅要考慮網(wǎng)格比的大小，而且還應(yīng)考慮 與h的比值（即， Du Fort-Frankel格式的影響），這樣才能達(dá)到所要求的精度。4.3 結(jié)果的收斂性分析我們選取的模型如下表，在表4列出，網(wǎng)格比保持不變，時(shí)間步長和空間步長變小 對誤差影響子區(qū)域的隱式點(diǎn)數(shù)模型3 網(wǎng)格比計(jì)算的步數(shù)6319 166表3 節(jié)點(diǎn)數(shù) Error 1279 1.7415E-5 2559 2.

18、2666E-6 5119 1.7021E-7 51199 2.458E-111000時(shí)間步長1.E-3 2.5E-4 6.25E-5 6.25E-7空間步長 2.454E-3 1.227E-3 6.136E-4 6.136E-5 Aerror 1.7311E-2 9.012E-3 2.7070E-3 3.9100E-5我們從表3的結(jié)果可以看出，固定網(wǎng)格比，隨著時(shí)間步長和空間步長變小，不管絕對 誤差還是相對誤差都是減小的，因此從數(shù)值的實(shí)驗(yàn)結(jié)果來看，此差分格式是收斂的。4.4結(jié)果的穩(wěn)定性分析我們從兩個(gè)方面來考察格式的穩(wěn)定性，其一， 不同網(wǎng)格比對誤差的影響，其二，長時(shí)間計(jì)算對誤差的影響。首先，我們考

19、察不同網(wǎng)格比對誤差的影響，我們選取如下的模型模型4計(jì)算的節(jié)點(diǎn)數(shù) 127999分段數(shù) 12800網(wǎng)格比 時(shí)間步長 S-maxerror我們從表4的數(shù)值結(jié)果可以看出，此區(qū)域分解方法是無條件穩(wěn)定的。其次，我們從長時(shí)間的計(jì)算能否保證精度，來考察此區(qū)域分解方法的穩(wěn)定性，計(jì)算的節(jié)點(diǎn)數(shù)12799計(jì)算步數(shù)5E+3 2E+4 1E+5表5 網(wǎng)格比66.4 4.8475E-7 1.8278E-6 6.6378E-6664 4.2840E-4 9.4916E-4 1.9314E-5模型5 分段數(shù) 1280空間步長 2.454E-4表4（S-maxerror：分子區(qū)域后的最大誤差）1660 16600 66401 1

20、E-6 1E-5 4E-5 1.4E-8 2E-6 3E-51660041E-4 2E-4空間步長 2.454E-5計(jì)算步數(shù) 200我們從表5的結(jié)果可以看出，在長時(shí)間計(jì)算情形下，此區(qū)域分解方法仍能夠保證精度， 這比較適合大型科學(xué)工程計(jì)算的要求。4.5 并行算法的可擴(kuò)展性分析我們?nèi)∽訁^(qū)域的隱式節(jié)點(diǎn)數(shù)為39，表6 列出不同節(jié)點(diǎn)數(shù)的在不同的處理機(jī)上計(jì)算 10000個(gè)時(shí)間步的模擬結(jié)果。 表6（Nodes：節(jié)點(diǎn)個(gè)數(shù)，Nseg：分段數(shù)，NCPU：處理機(jī)個(gè)數(shù)，Wtime：花費(fèi)的CPU時(shí)間）Nodes 1279 2559 5119 10239Nseg 32 64 128 256NCPU 4 8 16 32Wt

21、ime 2.6333 2.6398 2.6613 2.71224.6加速比分析 從表6和圖2的結(jié)果可以看出，并行算法具有良好的可擴(kuò)展性。我們從兩個(gè)方面考察加速比，一方面，串行優(yōu)化加速比：劃分多個(gè)子區(qū)域后單機(jī)上運(yùn)行與全隱式格式在單機(jī)上運(yùn)行相比獲得的加速比；另一方面，并行優(yōu)化加速比：劃分多個(gè)子區(qū)域后在多個(gè)處理機(jī)上并行運(yùn)行與全隱式格式在單機(jī)上運(yùn)行相比獲得的加速比。我們首先考察串行優(yōu)化獲得的加速比，考察以下模型模型6節(jié)點(diǎn)個(gè)數(shù)1279表7（Nseg：劃分的子區(qū)域數(shù)，Wtime：并行執(zhí)行的CPU時(shí)間，Sp：并行獲得的加速）Nseg 1 2 4 8 16 32Wtime 18.86 8.95 3.93 1.

22、11 0.58 0.32Sp 1 2.107 4.799 16.99 32.517 58.938從表7的結(jié)果不難分析出，縮小單機(jī)的存儲(chǔ)量，提高數(shù)據(jù)的局部性，Cache的性能大大空間步長 2.45437E-03 時(shí)間步長 4.000E-04 網(wǎng)格比 66.40185 7地提高，從而能夠獲得良好的優(yōu)化效果。其次，我們考察并行優(yōu)化獲得的加速比模型7時(shí)間步長4.000E-07 節(jié)點(diǎn)個(gè)數(shù) 空間步長 網(wǎng)格比 分段數(shù)目 127999 2.45437E-05 664.0185 1280表8（NCPU：處理機(jī)的個(gè)數(shù)，Wtime：并行執(zhí)行的時(shí)間，Sp：并行獲得的加速）NCPU 1 2 4 8 16 32Wtime

23、 45.777 22.910 11.460 5.737 2.996 1.475Sp 1 1.9976 3.9375 7.7077 15.403 31.064節(jié)點(diǎn)個(gè)數(shù)127999表9（NCPU：處理機(jī)的個(gè)數(shù)，Wtime：并行執(zhí)行的時(shí)間，Sp：并行獲得的加速） NCPU 1 2 4 8 16 32Wtime 9.257 4.634 2.351 1.201 0.601 0.298 Sp 1 1.9981 3.9995 7.9792 15.279 31.035從表8和表9測得的結(jié)果可以看出，對具有并行本性的格式進(jìn)行并行化優(yōu)化，獲得的并行加速比幾乎是線性的，這就提示我們，在并行算法設(shè)計(jì)時(shí)應(yīng)考慮數(shù)據(jù)的局部

24、性（并行本性），這樣，通信量就很小，因此，可以獲得良好的加速比。從表8與表9的結(jié)果相比較來看，在相同數(shù)目處理機(jī)的情形下，通信量是相同的，但是，劃分的子區(qū)域數(shù)不一樣（即，在程序設(shè)計(jì)時(shí)，體現(xiàn)在內(nèi)存容量是不同的），表8計(jì)算的時(shí)間是表9的的5倍（表8中隱式點(diǎn)的個(gè)數(shù)是表9中隱式點(diǎn)的個(gè)數(shù)的10倍）?？臻g步長 2.45437E-05 模型8 時(shí)間步長 4.000E-07 網(wǎng)格比 664.0185 分段數(shù)目 128005 結(jié)論本文考察熱傳導(dǎo)方程的一類具有并行本性的差分格式。證明了該格式滿足絕對穩(wěn)定的的必要條件，數(shù)值驗(yàn)證了格式比純Du Fort-Frankel格式精度高且具有與純隱式格式相同的二階精度，而且當(dāng)網(wǎng)格比很大時(shí)，格式仍是穩(wěn)定的。在對程序優(yōu)化前，我們首先考慮格式是否具有并行本性