1、熱傳導方程的一種新的具有并行本性的區域分解算法及其并行程序設計1左風麗*1 袁光偉*2北京應用物理與計算數學研究所應用科學計算部 高性能計算中心，北京 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. 引言許多物理和力學問題的求解都可歸結為熱傳導方程的求解，其數值模擬常采用有限差分方法，在實際應用問題中通常采用隱式差分格式，這是因為隱式差分格式具有無條件穩定的特征，對任意的網格比都能保證對計算誤差

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

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

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

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

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

9、甚至更大）可以大規模的（即，并行算法具有良好的可擴展性）；（）劃分的子區域包含的節點數也是任意選取的；（）我們采用動態局部數據管理的方法：任意單機的存儲量僅僅是對任意一個子區域生成的方程組的系數的存儲，當計算下一個子區域的值之前，重新生成方程組的系數矩陣，這樣，子區域的方程組的系數矩陣是一個臨時數組，每次計算子區域的值都要更新一次系數矩陣臨時數組的值）。第五節，我們給出數值結果的性能分析，包括求解精度的分析，并給出穩定性、收斂性和可擴展性的數值驗證。6-42 具有并行本性差分格式的構造我們考察如下的一維熱傳導方程：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),T>0,U=U0(x)是初始條件，滿足U0(0)=U0(l)=0。空間步長為h，時間步長，用平行線x=xj=jh(j=0,1,.,J)，t=t=n(n=0,1,.,N)劃分區域QT，其中，Jh=l,N=T,J和N是正整數。不失一般性，我們考慮將0,x分成n 2兩個子區域0,x1和0,x2的情形。分成多個子區域并沒有實質上的區別。令k滿足1<k<J1的正整數。對問題（1）（3）建立如下的具有并行本性的區域分解格式：在子區域的內邊界x=xk處使用Du Fort-Frankel三層格式+

11、11ununjjn+11nun+unj+1(ujj)+uj12=h2 （4）在子區域的內部使用純隱式格式 +1ununjj+1n+1+1un+unj+12ujj1 n+1=h2 （5） 離散初邊界條件為： U0 =UJn+1=0， 0nN （6）n+1u0 0jJ （7） j=U0(xj)，由格式（4）可首先顯式地計算出Uk。當Uk計算出后，可由（5）、（6）分別并行求出+1+1un(1jk)和un(kjJ1)。 jjn+1而格式（4）的截斷誤差為O(+h+ 眾所周知，格式（5）的截斷誤差為O(+h)，222h2)。如果將空間區域僅分為有限多個子區域，即內邊界點的個數為O(1)（該數為與h無關

12、的任一固定數），那么，我們可以期望所得到的解的收斂階為O(+h+/h)，數值實驗也將證實這一點，數值實驗表明，這里構造的格式（4）、（5）的精度比純Du Fort-Frankel的精度要高。由于Du Fort-Frankel格式和純隱式格式均為無條件穩定的，所以我們可以期望關于格式（4）（7）也是無條件穩定的。我們可以用矩陣的方法證明，此方法滿足無條件穩定的必要條件。223 并行算法設計因為在內邊界點采用三層Du Fort-Frankel格式，所以，第一時間層uj,1jJ1 的計算需采用其他格式。首先，我們設計如下的并行方法：+11ununjjn+11nun+unj+1(ujj)+uj112+

13、1ununjj=h2，在子區域的內邊界點，n=1,2,.,N （8）=+1n+1+1un+unj+12ujj1h2，在子區域的內點，n=0,1,.,N （9）以及如下的純顯式格式+1ununjjnnunj+12uj+uj11=h2，在子區域的內邊界點，n=0 （10）注1：uj,1jJ1的計算也可完全采用純隱式格式。經數值實驗知，這兩種方法的精度是非常接近的，所以下文討論根據格式（8）（10）。注2：本文討論的問題是針對一維情形，所以，這里所指的子區域是由某條子線段上的離散點組成，在進行并行算法設計時，開辟的存儲空間大小僅為子區域的大小，對任一個子區域所求解方程組的系數都要動態生成。因此，不管

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

15、Interface）。42 精度分析 我們從兩個方面考察格式數值解的精度，一方面，不同網格比對精度的影響，另一方面，不同分段數對精度的影響。首先，我們來看不同網格比對精度的影響，節點數 127999子區域的隱式點數99模型1 計算步數 1000分段數 12802表1(某并行機) Error: 絕對誤差 Aerror :Error/(+h)Error Aerror 時間步長 網格比（/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

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

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

18、 1.7415E-5 2559 2.2666E-6 5119 1.7021E-7 51199 2.458E-111000時間步長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的結果可以看出，固定網格比，隨著時間步長和空間步長變小，不管絕對 誤差還是相對誤差都是減小的，因此從數值的實驗結果來看，此差分格式是收斂的。4.4結果的穩定性分析我們從兩個方面來考察格式的穩定性，其一， 不同網格比對誤差的影響，其二，

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

20、1660 16600 66401 1E-6 1E-5 4E-5 1.4E-8 2E-6 3E-51660041E-4 2E-4空間步長 2.454E-5計算步數 200我們從表5的結果可以看出，在長時間計算情形下，此區域分解方法仍能夠保證精度， 這比較適合大型科學工程計算的要求。4.5 并行算法的可擴展性分析我們取子區域的隱式節點數為39，表6 列出不同節點數的在不同的處理機上計算 10000個時間步的模擬結果。 表6（Nodes：節點個數，Nseg：分段數，NCPU：處理機個數，Wtime：花費的CPU時間）Nodes 1279 2559 5119 10239Nseg 32 64 128 2

21、56NCPU 4 8 16 32Wtime 2.6333 2.6398 2.6613 2.71224.6加速比分析 從表6和圖2的結果可以看出，并行算法具有良好的可擴展性。我們從兩個方面考察加速比，一方面，串行優化加速比：劃分多個子區域后單機上運行與全隱式格式在單機上運行相比獲得的加速比；另一方面，并行優化加速比：劃分多個子區域后在多個處理機上并行運行與全隱式格式在單機上運行相比獲得的加速比。我們首先考察串行優化獲得的加速比，考察以下模型模型6節點個數1279表7（Nseg：劃分的子區域數，Wtime：并行執行的CPU時間，Sp：并行獲得的加速）Nseg 1 2 4 8 16 32Wtime

22、18.86 8.95 3.93 1.11 0.58 0.32Sp 1 2.107 4.799 16.99 32.517 58.938從表7的結果不難分析出，縮小單機的存儲量，提高數據的局部性，Cache的性能大大空間步長 2.45437E-03 時間步長 4.000E-04 網格比 66.40185 7地提高，從而能夠獲得良好的優化效果。其次，我們考察并行優化獲得的加速比模型7時間步長4.000E-07 節點個數 空間步長 網格比 分段數目 127999 2.45437E-05 664.0185 1280表8（NCPU：處理機的個數，Wtime：并行執行的時間，Sp：并行獲得的加速）NCPU

23、1 2 4 8 16 32Wtime 45.777 22.910 11.460 5.737 2.996 1.475Sp 1 1.9976 3.9375 7.7077 15.403 31.064節點個數127999表9（NCPU：處理機的個數，Wtime：并行執行的時間，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測得的結果可以看出，對具有并行本性的格式進行并行化優化，獲得的并行加速比幾乎是線性的，這就提示我

24、們，在并行算法設計時應考慮數據的局部性（并行本性），這樣，通信量就很小，因此，可以獲得良好的加速比。從表8與表9的結果相比較來看，在相同數目處理機的情形下，通信量是相同的，但是，劃分的子區域數不一樣（即，在程序設計時，體現在內存容量是不同的），表8計算的時間是表9的的5倍（表8中隱式點的個數是表9中隱式點的個數的10倍）。空間步長 2.45437E-05 模型8 時間步長 4.000E-07 網格比 664.0185 分段數目 128005 結論本文考察熱傳導方程的一類具有并行本性的差分格式。證明了該格式滿足絕對穩定的的必要條件，數值驗證了格式比純Du Fort-Frankel格式精度高且具有與純隱式格式相同的二階精度，而且當網格比很大時，格式仍是穩定的。在對程序優化前，我們首先考慮格式是否具有并行本性（數據的