热传导方程的一种新的具有并行本性的

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格式精度高且具有与纯隐式格式相同的二阶精度，而且当网格比很大时，格式仍是稳定的。在对程序优化前，我们首先考虑格式是否具有并行本性（数据的