软件学报
JOURNAL OF SOFTWARE
1999年 第2期 No.2 1999



非均匀区域油藏模拟负载平衡的分区并行算法*
舒继武1,2　赵金熙1,2 周维四3 张德富1,2
　　摘要　基于分布式并行计算机系统,对一类非均匀区域的油藏数值模拟问题,采用了区域分解方法并行求解,给出了并行求解的负载平衡模型及区域负载平衡的一种有效分区算法,从而将这类油藏数值模拟问题均衡映射到并行环境中,高效地并行求解.在黑油油藏数值模拟并行软件的研究中,实验结果表明,该算法有利于提高加速比.
　　关键词　油藏数值模拟,区域分解,负载平衡,并行算法.
　　中图法分类号　TP301
Zonal Parallel Algorithm of Load-balancing on Nonuniform Domain Simulation of Reservoir 
SHU Ji-wu1,2 ZHAO Jin-xi1,2 ZHOU Wei-si3 ZHANG De-fu1,2 
　　Abstract　 Based on the distributed parallel environment, for one of the nonuniform domain numerical simulation problems of reservoir, the authors give a load-balance model of parallel computing and a zonal parallel algorithm of load-balance using domain decomposition method, so some well and effective solutions which balancly map the nonuniform domain numerical simaulation of reservoir to parallel environment are found. In the study of block oil reservoir numerical simulation parallel software, the practical results show that the zonal parallel algorithm can improve the performance in speedup.
　　Key words　Numerical simulation of reservoir, domain decomposition, load-balancing, parallel algorithm.
　　油藏数值模拟问题,无不归结为求解大型偏微分方程的问题.这类问题的求解需要大量的计算时间,因此,如何提高求解这类问题的计算速度成为人们非常关心的问题.区域分解法是应运而生的、适应并行计算机工作原理的偏微分方程数值算法.它体现了“分而治之”的思想,受到计算数学界和工程界的热切关注.该方法从问题的“模型级”发掘并行性,将计算区域分割成若干独立的规模较小的子区域,将原问题的求解转化为各子区域上子问题的求解.一个或若干个子问题用一个处理机来处理,当所有的子问题均衡地映射到多个处理机上时,通过并行求解子问题而获得整个区域上大问题的解,它有利于问题的并行求解,提高了解题速度,具有广阔的应用前景.
　　区域分解并行计算,在结合并行机的特有体系结构时,要考虑区域的划分及负载平衡.文献［1］基于二维网格连接的分布式并行处理机结构,在给出几种区域分解方法的计算复杂性和通信耗费的基础上对其加速比和效率进行了分析.文献［2］对区域的网格划分与并行计算的效率进行了讨论,并寻找其最优的网格分裂.它们都假设区域均匀,区域划分为等大小的子区域,认为子区域有完全相同的计算量和通信要求.但在实际的问题中,如油藏模拟,油区经常是不均匀的,而在选定区域时,为了便于处理,常选定一规则区域,如长方形区域.若将区域划分成与处理机数目相同的等大小的子区域,各子区域的收敛速度就不同,对应分配到处理机上将引起处理机的负载不均衡.因此,要维持一种相对稳定的高效并行处理,需采用负载平衡的分区并行计算方法.本文以在油藏数值模拟中,运用区域分解方法并行求解非均匀区域的黑油压力问题为背景,给出了并行求解的负载平衡模型及区域负载平衡的一种分区算法.
1 区域分解并行算法
　　区域分解法主要有3类:迭代子结构类型,子区域间无重叠;Schwarz类型,子区域间有重叠;迭代子结构和Schwarz的结合类型.
　　下面以Schwarz类型的区域分解算法来说明其并行求解问题的基本步骤.设有一个模型为:
　　　　　　　　　　　　　(1)
对区域Ω进行分解:令Ω有分解Ω=Ωi,Ωi是开集,且对每个i都可找到j≠i,有Ωi∩Ωj≠.对
问题（1）用平均法构造的并行Schwarz方法如下［3］:
　　步骤1: 选择初始近似P0,n＝0.
　　步骤2: 并行计算子区域上边值问题
　　　　　　　　　　　(2)
　　步骤3: 延拓Pn+1i的定义到式(3),即定义
　　　　　　　　　　　　(3)
并取平均
　　　　　　　　　　　　　　(4)
置n=n＋1,若达到问题精度,则终止;否则,转步骤2.
　　因此,对于一般的数值计算问题,应用区域分解法并行求解问题的基本过程为:
　　.将给定的实际问题的求解区域,分割为几个子区域.
　　.将各个子区域的求解,归结于解线性代数方程组.
　　.给各个处理机分配数据和任务.
　　.整个区域迭代循环:
　　do{
　　　处理机Pi上子区域的迭代循环（与其他处理机并行）:
　　　do{
　　　　　子区域迭代求解
　　　}while（未达到对应子区域的收敛条件）
　　　求子区域间边界节点（无重叠区域分解）或公共区域（重叠区域分解）
　　}while（整个区域未达到其收敛条件）
　　.整个区域收敛,得到问题的解. 
2 负载平衡的分区算法
　　用区域分解方法求解,任务的完成是等待最慢结点机上的消息,这样,导致整体并行性能退化到与最慢结点机相同.为了解决这个问题,需要做到:① 要提高网络传输性能;② 要采用高效的算法,使得各子区域收敛速度最快;③ 要采用负载平衡的分区方法,使得每个结点机上完成的时间大致相等.本文主要考虑区域分解方法并行求解非均匀问题的负载平衡.
2.1 非均匀区域的划分
　　在非均匀区域的油藏数值模拟中,有的子区域由于包含更多的“死节点（即不需迭代计算的节点）”,达到收敛条件时迭代次数很少,而有的子区域由于包含很少的“死节点”,达到收敛条件时迭代次数很多.若划分的子区域数目大于处理机的个数,由于在开始迭代时,不能估计出其迭代次数,为了能较好地使负载平衡,子区域的数目必须远大于处理机的个数,这样,由于划分的子区域太多,整个区域收敛速度减慢.这里划分的子区域数目等于处理机的个数,并在区域迭代求解过程中,根据本处理机上的局部状态信息和整个区域的全局状态来动态地改变区域的划分,从而达到任务最终能均衡地运行的目的.
　　设处理机集为:P={P1,P2,P3,...,Pm},Pi为处理机i,处理机数为m.Ω为一油藏模拟区域,分解为m个子区域,令Ω有分解Ω=Ωi,即子区域集Ω={Ω1,Ω2,...,Ωm },Ωi为第i个子域.
　　并行求解的每一迭代涉及到区域的并行迭代和负载平衡算法,在P个处理机上完成max-iter次迭代的总的执行时间［4］为
　　　　　　　　　　（5）
其中k为迭代步,i为处理机Pi（i=1,2,...,P）,,i表示处理机Pi在第k次迭代步的执行时间,,i表示处理机Pi在第k次迭代后负载平衡算法的时间.并行迭代过程中,计算和同步通信是交替进行的,而处理机Pi一次迭代的时间受3个因素的影响.
Tapp,i(k)=f(tcomp,tcomm,tid),　　　　　　　　　　　（6）
其中tcomp是一次迭代的计算时间,tcomm是迭代通信时间,tid是同步等待时间.在迭代收敛的情况下,区域的划分应以总的执行时间尽量短为目标.
2.2 并行计算的几点说明
　　根据上述分析,对油藏数值模拟中区域分解并行计算作如下几点说明:
　　（1） 油区是不均匀的,而选定的计算域是一规则的区域,如长方形区域.子区域网格离散化相同,各子域的求解方法可以不同,各子区域的收敛条件也可以不同.
　　（2） 在并行计算环境中,网络的带宽会由于环境中其他用户的介入或因本身任务而发生动态的变化,结果是消息通过网络延迟发生变化,这里不考虑网络传输性能及结点性能的变化.同一处理机上的子区域之间的通信代价为0,不在同一处理机上的子区域之间的通信代价只与其相邻区域的边界大小或公共区域大小有关.
　　（3） 处理机上计算和通信不重叠.
　　（4） 区域分解并行迭代时,整个区域的每一迭代步包括所有子区域的多次迭代,直至分别达到各子区域的收敛条件.整个区域的迭代步相对于各子区域的迭代步少得多.这里将每一个子区域的每一次收敛到该子区域的收敛条件的迭代过程视为一个任务,任务迁移使得相邻子区域间的边界移动,引起子区域大小的变化,从而将一部分迭代过程迁移到其他处理机上.本文将不加区分地使用子区域的每一次迭代收敛过程和任务.
　　（5） 在实际应用中,处理机的数目一般较少,而油藏模拟的区域一般很大,划分后的子区域的计算量和收敛速度与该子区域需计算的节点数、每次的迭代初值、所用算法等有关.
　　（6） 处理机间的连接是线性连接,计算区域是长方形（体）区域,区域的划分与处理机连接一一对应,子区域的数目与处理机的数目相等.每个处理机只与相邻的处理机交换数据.
2.3 负载平衡模型及分区算法
　　设整个区域的迭代步为k=1,2,...,在每一迭代步后,可根据处理机的信息来启动负载平衡算法,负载平衡算法涉及的时刻可以用k表示.
　　定义1. w（k）:并行系统在k时刻的状态-任务的分布向量.wi（k）:处理机Pi在区域第k步迭代后的任务量,则时刻k+1的任务分布包括:k时刻的任务分布、负载平衡引起的任务传递.处理机Pi的负载平衡模型可表述为
wi(k+1)=wi(k)-wi,l.d,s(k)+wi,l.d,r(k),　　　　　　　　　（7）
其中wi,l.d,s（k）,wi,l.d,r(k)分别表示负载平衡引起的发送和接收的任务.
　　定义2. A（k）:处理机参与矩阵,对于P个处理机的并行系统,它是一个p×p的矩阵,每个处理机只能与其相邻的处理机有任务的重新分配.如果处理机Pi和Pj参与任务迁移,则对应元素置为“1”,否则, 对应元素置为“0”.
　　定义3. M（k）:描述任务的迁移,它是一个p×p的矩阵.每个处理机只能与其相邻的处理机有任务的迁移.Mij:处理机Pi向Pj迁移任务的因子,如Pi送给Pj的任务可表示为:Mij（k）.wi（k）.
　　定义4. H:允许任务迁移的路径矩阵.Hij:如果处理机Pi允许向Pj迁移任务,则对应元素置为“1”,否则, 对应元素置为“0”.因此,Mij≠0仅当Hij=1.
　　根据上述的定义,式（7）表述为
wi(k+1)=wi(k)-Mij(k)Aii(k)wi(k)+Mji(k)Ajj(k)wj(k).　　　　（8）
上式中,第1个和项表示从处理机Pi到其相邻处理机的任务迁移;第2个和项表示从Pi的相邻处理机到Pi的任务迁移.定义:Mii(k)=-Mij(k),则负载平衡模型可以表述为
w(k+1)=w(k)+［M(k)］TA(k)w(k).　　　　　　　　　　（9）
　　定义5. :处理机的理想任务分布.
= 
　　定义6. :处理机Pi的相邻处理机任务分布的平均信息.
= 
其中d为Pi的相邻处理机数,Ri为Pi的相邻处理机集.
　　 和可分别作为处理机任务分布的全局信息和局部信息.
　　根据上述定义,将处理机Pi第k次迭代完成的时间,i记为其在k时刻的任务量wi(k),采用负载平衡的分区并行算法可表述为:
　　（1） 将区域均分为子区域Ω1,...,Ωm,并对应地分配到处理机P1,...,Pm上.令初始集为S={1,2,...,m}.
　　（2） 应用区域分解并行算法,各子区域在处理机上并行执行,直至达到各子区域的收敛条件;分别统计其迭代完成时间,i,记为wi(k).
　　（3） 求解子区域间的边界节点或公共子区域.
　　（4） 如果整个区域达到其收敛条件,则终止;否则,转（5）.
　　（5） 从S中选择满足i(k)=max｛p(k),p∈S｝的处理机Pi及处理机Pl,Pl满足l(k)=min
｛p(k),p∈S｝.如果i-l≤α(σ为处理机间容许的负载差）,转（2）;否则,转（6）.
　　（6） 在Pi的相邻处理机集中选择Pj,满足j(k)=min｛wp(k),p∈Ri｝（Ωi的相邻区域集记为ΩRi）,如果i-j＞α,处理机Pi的部分任务需迁移到Pj上,即需在域Ωi与其邻域Ωj之间重新划分区域.Pi迁移到Pj的任务量约为:wi(k)-k(k);否则,转（7）.
　　（7） S=S-{i},如果S≠,转（5）;否则,重置S={1,2,...,m},转（2）.
3 计算实例
　　在黑油数值模拟并行化软件开发的研究中,使用了区域分解并行计算,获得了较好的效果.下面是一油藏模拟中黑油压力问题模型,
(10)
其中
Ct=Cr+CoSo+CwSw+CgSg, Pcow=Po-Pw, Pcgo=Pg-Po, So+Sw+Sg=1 
这里,下标o,w,g分别为油、水、气相,下标sc为地面标准条件,B为体积系数,Rso为溶解气油比,q为注入或采出的体积流量,S为饱和度,P为压力,Pc为毛管压力,Φ为孔隙度,ρ为密度,L1为边界条件算子,t为时间,f是关于时间的压力函数.
　　对于上述模型的这种大型边值问题,如应用迭代子结构类型的区域分解方法,其离散模型从总体上说都可化为一个系数矩阵加边对角块的线性方程组.将区域Ω分成互不相交的m个子区域Ωi,i＝1,2,...,m,若用步长为h的正方形网格,用7点差分格式离散方程,可得到如下线性方程组［5］.
　　　　　　（11）
式中A1,1,A2,2,...,Am,m分别对应于子区域Ω1,Ω2,...,Ωm上的矩阵.式（10）的求解可运用迭代法,以前一阶段的近似解作为下一阶段的初始量,使整个问题收敛.对每一迭代步中的每个子区域问题Aj,jpj=fj-Ａj,m+1pm+1（j=1,2,...,m）的求解,可以对应地分配到m个处理机上并行求解,在每一迭代步后,可以根据负载平衡算法,重新划分区域,形成新的方程组
　　　　　　（12）
　　对新形成的子区域问题:A′j,jp′j=f′j-A′j,m+1p′m+1（j=1,2,...,m）,仍对应在m个处理机上并行求解.这样,经过负载平衡算法,各处理机的任务较先前的均衡.计算实例为:区域x,y,z方向网格数分别为170、50、5,节点规模为42500的上述黑油模型,沿平行于z方向剖分区域,该区域中非均匀地分布着许多“死节点”.运用区域分解法并行计算30个阶段（时间段）.区域收敛判别为:‖Pi+1-Pi‖∞≤10-4,Pi+1,Pi分别为整个区域的第i＋1迭代步和第i迭代步的压力.各子区域采用LU预处理广义极小剩余法（GMRES（m））求解,取GMRES（m）算法中的截断参数m为10,各子区域收敛判别为‖pm-p1‖∞≤m*10-4,pm,p1分别为解子区域的GMRES（m）算法的第m迭代步和第1迭代步压力.计算环境是用PVM通过以太网分别连接2个、3个、4个、5个SGI indy构成几种分布式并行计算系统.图1给出了该黑油数值模型采用负载平衡方法和无负载平衡时的加速比的比较.实验结果表明,采用负载平衡方法能使系统并行加速比提高10％～40％.如对4个处理机的并行系统在无负载平衡时,上述模型计算30个阶段的加速比为2.51,采用本文的负载平衡分区算法后其加速比为3.25,加速比提高29.5%.

图1　加速比比较
本文研究得到国家863高科技项目基金资助.
作者介绍：舒继武,1969年生,博士生,主要研究领域为并行处理和分布式计算.
　　　　　赵金熙,1950年生,博士,副教授，主要研究领域为计算数学,数值计算.
　　　　　周维四,1942年生,教授级高工,主要研究领域为计算技术和油藏模拟.
　　　　　张德富,1937年生,教授,博士生导师,主要研究领域为计算机软件,并行处理,分布式计算.
本文通讯联系人:舒继武，北京 100084,清华大学计算机科学与技术系
作者单位：舒继武　赵金熙　张德富　南京大学计算机软件新技术国家重点实验室　南京 210093　南京大学计算机科学与技术系　南京　210093
　　　　　周维四　胜利油田地质科学研究院 东营 257000
参考文献
　［1］ Hoffmann K H, Zou Jun. Parallel efficiency of domain decomposition methods. Parallel Computing, 1993,19（12）:1375～1391
　［2］Chan T F, Shao J P. Parallel complexity of domain decomposition methods and optimal coarse grid size. Parallel Computing, 1995,21（7）:1033～1049
　［3］吕涛,石济民,林振宝等.区域分解法——偏微分方程数值解新技术.北京:北京科学技术出版社,1992(L? Tao, Shi Ji-min, Lin Zhen-bao et al. Domain Decomposition Methods——New Technology for Numerical Solver of Partial Differential Equations. Beijing: Beijing Science and Technology Press, 1992)
　［4］Franklin M A, Govindan V. A general matrix iterative model for dynamic load balancing. Parallel Computing, 1996,22（7）:969～989
　［5］赵金熙,舒继武,张德富等.大型对称不定箭形线性方程组的分解方法.高等学校计算数学学报,1997,19（1）:83～92(Zhao Jin-xi, Shu Ji-wu, Zhang De-fu et al. A direct method for solving symmetric indefinite arrow linear system. Numerical Mathematics, A Journal of Chinese Universities, 1997,19（1）:83～92)
本文1997-11-13收到原稿,1998-03-03收到修改稿
