串行油藏模拟器并行化的一种便利途径.doc

串行油藏模拟器并行化的一种便利途径*赵国忠* 本文是作者中国科学院软件研究所并行软件研究开发中心访问期间完成的；并得到国家高技术研究发展计划（八六三）课题“大规模整体油田的精细油藏数值模拟”的资助，课题编号为“863-306-ZD11-03-1”大庆油田有限责任公司勘探开发研究院）摘 要 本文给出了一种借助于PETSc[1]来实现串行油藏模拟器并行化的便利途径文中详细论述了并行策略及并行化过程中一些重要技术环节的实践考虑黑油IMPES情况下的实现实例表明，按此途径可在短时间内搞出具有一定实用程度的并行模拟器该途径对于其它类型的油藏模拟器如组分、热采和化学驱等也有一定的适用价值主题词 油藏 模拟器 并行化 途径1 引 言油藏数值模拟自从诞生之日起就一直跟随计算机硬件的发展，来满足石油开采业的技术需求，其求解规模越来越大，模拟的主要目的也已从过去的注重油藏整体动态趋势发展到研究地下流体的精细分布状况这无疑使每次模拟的计算开销越来越大十几年前，发达国家出现了并行计算环境（共享或分布式的），与此同时一些油藏模拟器开发者就开始研究串行模拟器的并行化近五年来，我国石油行业先后引进了一些并行计算机，也开始了对80年代中期所引进的源码油藏模拟器的并行化探讨。

然而由于并行实现对特定并行计算环境的依赖很强，加之我们对原串行模拟器的消化和理解程度也有限，实用的国产并行模拟器仍未推出迫于油田实际应用的需要，我国只好再一次引进数值模拟技术大庆油田、新疆油田、北京石油勘探开发研究院和胜利油田等先后引进了PVIP但这次软件商并不出售其源码了，且售价居高不下，对用户的使用也有种种限制，如在机型、处理器个数等方面我们将来如对硬件扩容、更新，还要不断在软件上投入在这种情况下，实现原有串行模拟器的并行化，是摆在油藏模拟专业人员面前的一个不可回避的问题油藏数值模拟过程中所发生的计算可以分成耦合运算和非耦合运算两类非耦合运算的并行只涉及程序实现技术，主要体现在数据的区域分解上，并行程序可大量沿用串行源码而耦合运算主要发生性系统的求解过程中，必须使用不同于串行时的并行求解方法因此，大型、稀疏、非对称线性系统的并行解法开发与实现是串行模拟器的并行化的技术关键所在中国科学院并行软件研究与开发中心于1993年开始在这样的算法和实现方面都做了大量的基础工作，并以SimBestⅡ的ESPIDO串行求解器为基础研制出了并行求解器，现已在几种并行环境下成功实现了对若干大规模实例的运行测试，这在串行模拟器的并行化方面走出了一条成功之路。

然而，这条路对于油藏模拟器普通开发者来说是很困难的美国Argonne国家实验室推出的微分方程并行求解软件包PETSc（Portable，Parallel，Extended Toolkit for Scientific Computation）的出现，使油藏模拟器开发者快速实现原有串行模拟器的并行化成为可能对此，本文通过对PETSc部分构件的使用，以黑油模拟器之IMPES选件为例，探讨了这种可能性愿本文的工作对有志于把串行油藏模拟器并行化的同行们有借鉴意义2 总的并行策略油藏数值模拟计算一般分为初始化和时间循环两个阶段初始化主要完成初始条件的设置，其发生的全部运算均是一次性的，但程序代码却不少，若搞并行也会付出很多复杂的编码工作，所以一般不对初始化部分作并行如果并行环境是分布式的，由于单个节点内存的有限，有时也可考虑这部分的并行本文讨论的并行化途径对于共享式并行环境或某一节点的内存能够满足实际问题初始化要求的分布式并行环境都是实用的另外，虚拟共享技术也能弥补分布式环境下单一节点内存的不足2.1 串行求解的一般步骤在谋划具体的并行策略之前，把握串行求解的一般步骤是必要的串行油藏数值模拟计算的时间循环阶段一般遵循以下主要步骤： ⑴ 若非本次时间循环之第一时间步，则变量更新，按用户控制打印输出一些信息；⑵ 若已达下一控制时刻，则输入控制数据，若遇STOP键词则转向第⑽步；⑶ 时间步长的设置；⑷ 雅可比矩阵的计算；⑸ 井的处理；⑹ 若走IMPES，消去非压力方程，形成单一未知量线性系统；⑺ 网格方程与井方程的耦合，解线性系统；⑻ 若未达最小牛顿迭代次数或不收敛，则转到第⑷步；⑼ 若收敛或已达最大牛顿迭代次数，则转到第⑴步；⑽ 退出时间循环，打印作业统计信息，中止程序的运行。

以上步骤仅是并行化时所必须考虑的，此外还有许多细节，如物质平衡分析和时间步截断等，这里就不触及了2.2 计算任务的并行分配我们知道，计算量将主要发生在上述第⑷、⑹和⑺步大量的实际算例表明，串行时消耗在这三步之外的计算时间不足全部时间的2%所以要得到好的并行效率，这三步的计算及有关数据必须分配到各个并行节点上其他部分若搞并行并不是不能实现，但为此而作出的通讯开销是相当多的，反而得不偿失………………………CALL MPI_INIT();CALL MPI_COMM_SIZE();CALL MPI_COMM_RANK();IF(RANK.EQ.0)THEN… ! CODES RUN IN MASTER PROCESSORELSE… ! CODES RUN IN SLAVE PROCESSORSENDIFCALL MPI_FINALIZE();………………………图1 用MPI实现主从方式代码示意为了尽可能减少并行编码的工作量，采取主从式（Master/Slave）的并行实现策略是便利的具体来说可以由从进程承担上述第⑷、⑹和⑺步的计算任务，而由主进程来承担剩下的输入、输出、井的处理和求解过程的控制等任务。

并行平台采用流行的MPI标准信息传递界面（PETSc必须有MPI的支持），进程间通讯效率高，且可在不同的系统环境下源码移植主从式的并行实现对于PVM是便利的，因PVM提供了主进程生成子进程的机制，主进程和子进程运行不同的可执行码不过在MPI下，我们可用同一个可执行码通过逻辑控制实现主从方式图1给出了简略的FORTRAN分支程序示意，其中略去了MPI调用的形参，rank表示MPI通讯子内进程的编号2.3 数据通讯的实际考虑主进程与每一个从进程均要发生信息传递从进程启动后的第一件事是把所有不依赖于时间的数据一次性从主进程接受过来；每一时间步结束之前，主进程应自从进程接收网格压力、饱和度和油气比变量用于物质平衡分析和可能的打印输出；每一牛顿迭代步各从进程要把差分方程余量和未知量变化的局域最大绝对值传递给主进程，由主进程汇总后形成迭代控制数据再传回各从进程相邻区域所在进程之间也要进行必要的通信每一牛顿迭代步相邻子域应传递边界网格未知量的变化，用于子域内的变量更新和便于下一牛顿步的雅可比矩阵计算3 数据的分配与内存的管理一个油藏模拟器的数据可以分为标量数据和网格数组数据两类大量的存储消耗在后者上，原则讲后者应该分配给各从进程局部存放，但由于输出的需要，仍然要在主进程中全局存放象压力、饱和度和油气比等几个主要网格数组。

这里所说的标量数据还包括一些与网格数和井数无关的小数组，如相对渗透率和PVT关系表数据等，其占用的内存是微不足道的，所以可在各进程中重复存放3.1 数据分区 由于直井一般会占大多数，为了不必要地增加编码的复杂度，只要允许用户在X、Y两个方向上任意划分区域，就足以避开区域分界线与斜井、水平井或断层相交的情况得到用户的分区信息后，作少量的编码易得满足上述条件的分区方案，当然应尽可能使各区域内的局部网格数接近，这有利于加载平衡区域的编号可任意考虑到区域间自动信息交换的要求，区域边界线的编号应先一个方向，再另一个方向分区后应准确存储续用信息，如区域、区域边界线的位置、尺度及编号等另外，区域内的井也应再进行编号，且要保留域内编号与原编号的关系3.2 内存的管理主进程启动后，读重启文件，按记录要求申请足够的内存，读入全部初始信息待把应该由各子域存储的数据分别传过去后，备份必要数据，释放并重新申请动态内存，再导入备份数据这样可使主进程占用更少的内存从进程则按各自需要分别申请动态内存4 雅可比矩阵计算的并行雅可比矩阵的计算是可精确并行化的考虑到边界网格方程的雅可比元素与相邻子域对应的边界网格数据有关，此阶段我们不妨把各子域的网格系统在内边界处向外扩一排，使当前的虚拟网格系统包括与该子域雅可比矩阵计算有关的所有网格。

若用NXD（I）和NYD（I）表示子域I的方向实际网格数，用NXV（I）和NYV（I）表示子域I的方向虚拟网格数，则对于图2示例中的区域2（以斜线为背景的网格是外扩的虚拟网格），分别有NXD（2）=4、NYD（2）=5、NXV（2）=6和NYV（2）=6这样，我们可以把每一个子域当作一个独立的模型套用原来串行程序的雅可比计算源码在虚拟网格系统的数组数据按自然顺序正确设置后，唯一需要改变的是对于第I子域，用NXV（I）和NYV（I）分别替换原来的NX和NY在图2示例中NX=11、NY=9这样安排避免了雅可比计算过程中信息传递的频繁发生，但比串行增加了总的计算量，不过只要在代码中稍加逻辑判断即可消除此处设置一个局域整型数组是必要的，便于雅可比计算完成后把结果从虚拟网格向实际网格映射，用于后续对线性系统的求解123456图2 子域内雅可比计算的虚拟网格系统5 用PETSc并行求解线性系统5.1 PETSc简介PETSc是一个可在多种操作系统平台下运行、可扩展的大规模科学计算并行软件包，适合于对常／偏微分方程的并行求解它主要由数据对象、数据与网格管理、线性系统求解器（SLES）、非线性系统求解器（SNES）和常微分方程求解器（TS）等五个层次内的若干构件以及一些系统实用工具组成。

数据对象主要包括向量（VEC）和矩阵（MAT），是PETSc最基础构件所在的层次；SLES主要包括子空间方法（KSP）和预条件子（PC），是PETSc目前最成熟和构件被使用最多的层次用户可以根据自己的需要，部分或全部地使用这些构件来开发求解某一实际问题的并行应用软件另外，PETSc还为高级用户提供了自行扩展构件的良好机制PETSc2.0.28同时支持在MPI环境下的C和FORTRAN两种调用5.2 线性系统的分区装配 用PETSc的SLES并行求解油藏数值模拟中产生的线性系统，要获得较好的并行效率关键在于矩阵的装配SLES的存储方案是，把矩阵中属于某子域的行都交由对应进程来存储按3.1节所述的数据分区方案，由某一子域内节点方程和井方程导出的行的元素恰好在对应的进程内，这就保证了矩阵的局部装配过程中不发生子域间数据的相互传递PETSc的开发和释放者强烈地建议用户在矩阵的局部装配之前用D_NNZ和O_NNZ两个整型数组准确地告诉装配函数该矩阵的对角块和非对角块上非零元素的位置，使装配过程中尽量减少内存申请次数和便于按行压缩存储，以此获得并行求解的高效率解决这一问题的办法是先装配，在确认装配无误后再逻辑跟踪装配过程，准确记录非零元素的位置，形成须预先运行的子程序来设置D_NNZ和O_NNZ数组。

实际求解时不必在每一次之前都运行此程序，而只需在井层生产／注入状态发生变化的某一时间步（含一次模拟运算的第一时间步）的第一次牛顿迭代之前运行就可以了5.3 对PETSc并行构件的选用 SLES层次内提供了很多子空间迭代方法和串行预条件子，而并行预条件子只有块雅可比（BJACOBI）和加法Schwarz（ASM）两种实现这些构件的不同组合能用于求解不同的问题数值实验表明，GMRES[2]和BCGS[3。