43串行油藏模拟器并行化的一种便利途径*赵国忠（大庆油田有限责任公司勘探开发研究院）摘 要 本文给出了一种借助于 PETSc1来实现串行油藏模拟器并行化的便利途径。文中详细论述了并行策略及并行化过程中一些重要技术环节的实践考虑。黑油 IMPES 情况下的实现实例表明，按此途径可在短时间内搞出具有一定实用程度的并行模拟器。该途径对于其它类型的油藏模拟器如组分、热采和化学驱等也有一定的适用价值。主题词 油藏 模拟器 并行化 途径1引引 言言油藏数值模拟自从诞生之日起就一直跟随计算机硬件的发展，来满足石油开采业的技 术需求，其求解规模越来越大，模拟的主要目的也已从过去的注重油藏整体动态趋势发展 到研究地下流体的精细分布状况。这无疑使每次模拟的计算开销越来越大。十几年前，发 达国家出现了并行计算环境（共享或分布式的） ，与此同时一些油藏模拟器开发者就开始 研究串行模拟器的并行化。近五年来，我国石油行业先后引进了一些并行计算机，也开始 了对 80 年代中期所引进的源码油藏模拟器的并行化探讨。然而由于并行实现对特定并行 计算环境的依赖很强，加之我们对原串行模拟器的消化和理解程度也有限，实用的国产并 行模拟器仍未推出。迫于油田实际应用的需要，我国只好再一次引进数值模拟技术。大庆 油田、新疆油田、北京石油勘探开发研究院和胜利油田等先后引进了 PVIP。但这次软件商 并不出售其源码了，且售价居高不下，对用户的使用也有种种限制，如在机型、处理器个 数等方面。我们将来如对硬件扩容、更新，还要不断在软件上投入。在这种情况下，实现 原有串行模拟器的并行化，是摆在油藏模拟专业人员面前的一个不可回避的问题。 油藏数值模拟过程中所发生的计算可以分成耦合运算和非耦合运算两类。非耦合运算 的并行只涉及程序实现技术，主要体现在数据的区域分解上，并行程序可大量沿用串行源 码。而耦合运算主要发生在线性系统的求解过程中，必须使用不同于串行时的并行求解方 法。因此，大型、稀疏、非对称线性系统的并行解法开发与实现是串行模拟器的并行化的 技术关键所在。中国科学院并行软件研究与开发中心于 1993 年开始在这样的算法和实现 方面都做了大量的基础工作，并以 SimBest的 ESPIDO 串行求解器为基础研制出了并行 求解器，现已在几种并行环境下成功实现了对若干大规模实例的运行测试，这在串行模拟 器的并行化方面走出了一条成功之路。然而，这条路对于油藏模拟器普通开发者来说是很 困难的。 美国 Argonne 国家实验室推出的微分方程并行求解软件包 PETSc（Portable，Parallel，Extended Toolkit for Scientific Computation）的出现，使油藏模 拟器开发者快速实现原有串行模拟器的并行化成为可能。对此，本文通过对 PETSc 部分构 本文是作者中国科学院软件研究所并行软件研究开发中心访问期间完成的；并得到国家高技术研究发展计划（八六三）课题“大规模整体油田的精细油藏数值模拟”的资助，课题编号为“863-306-ZD11-03- 1” 。44 CALL MPI_INIT(); CALL MPI_COMM_SIZE(); CALL MPI_COMM_RANK(); IF(RANK.EQ.0)THEN ! CODES RUN IN MASTER PROCESSOR ELSE ! CODES RUN IN SLAVE PROCESSORS ENDIF CALL MPI_FINALIZE();图1 用MPI实现主从方式代码示意件的使用，以黑油模拟器之 IMPES 选件为例，探讨了这种可能性。愿本文的工作对有志 于把串行油藏模拟器并行化的同行们有借鉴意义。2总的并行策略总的并行策略油藏数值模拟计算一般分为初始化和时间循环两个阶段。初始化主要完成初始条件的 设置，其发生的全部运算均是一次性的，但程序代码却不少，若搞并行也会付出很多复杂 的编码工作，所以一般不对初始化部分作并行。如果并行环境是分布式的，由于单个节点 内存的有限，有时也可考虑这部分的并行。本文讨论的并行化途径对于共享式并行环境或 某一节点的内存能够满足实际问题初始化要求的分布式并行环境都是实用的。另外，虚拟 共享技术也能弥补分布式环境下单一节点内存的不足。 2.1串行求解的一般步骤串行求解的一般步骤 在谋划具体的并行策略之前，把握串行求解的一般步骤是必要的。串行油藏数值模拟 计算的时间循环阶段一般遵循以下主要步骤： 若非本次时间循环之第一时间步，则变量更新，按用户控制打印输出一些信 息； 若已达下一控制时刻，则输入控制数据，若遇 STOP 键词则转向第步； 时间步长的设置； 雅可比矩阵的计算； 井的处理； 若走 IMPES，消去非压力方程，形成单一未知量线性系统； 网格方程与井方程的耦合，解线性系统； 若未达最小牛顿迭代次数或不收敛，则转到第步； 若收敛或已达最大牛顿迭代次数，则转到第步； 退出时间循环，打印作业统计信息，中止程序的运行。以上步骤仅是并行化时所必须考虑的，此外还有许多细节，如物质平衡分析和时间步 截断等，这里就不触及了。 2.2计算任务的并行分配计算任务的并行分配 我们知道，计算量将主要发生在上述第、和步。大量的实际算例表明，串行时 消耗在这三步之外的计算时间不足全部时间的 2%。所以要得到好的并行效率，这三步的 计算及有关数据必须分配到各个并行节点上。其他部分若搞并行并不是不能实现，但为此 而作出的通讯开销是相当多的，反而得不偿失。45为了尽可能减少并行编码的工作量，采取主从式（Master/Slave）的并行实现策略是便 利的。具体来说可以由从进程承担上述第、和步的计算任务，而由主进程来承担剩 下的输入、输出、井的处理和求解过程的控制等任务。并行平台采用流行的 MPI 标准信息 传递界面（PETSc 必须有 MPI 的支持） ，进程间通讯效率高，且可在不同的系统环境下源 码移植。主从式的并行实现对于 PVM 是便利的，因 PVM 提供了主进程生成子进程的机制， 主进程和子进程运行不同的可执行码。不过在 MPI 下，我们可用同一个可执行码通过逻辑 控制实现主从方式。图 1 给出了简略的 FORTRAN 分支程序示意，其中略去了 MPI 调用 的形参，rank 表示 MPI 通讯子内进程的编号。 2.3 数据通讯的实际考虑数据通讯的实际考虑 主进程与每一个从进程均要发生信息传递。从进程启动后的第一件事是把所有不依赖 于时间的数据一次性从主进程接受过来；每一时间步结束之前，主进程应自从进程接收网 格压力、饱和度和油气比变量用于物质平衡分析和可能的打印输出；每一牛顿迭代步各从 进程要把差分方程余量和未知量变化的局域最大绝对值传递给主进程，由主进程汇总后形 成迭代控制数据再传回各从进程。 相邻区域所在进程之间也要进行必要的通信。每一牛顿迭代步相邻子域应传递边界网 格未知量的变化，用于子域内的变量更新和便于下一牛顿步的雅可比矩阵计算。3数据的分配与内存的管理数据的分配与内存的管理一个油藏模拟器的数据可以分为标量数据和网格数组数据两类。大量的存储消耗在后 者上，原则讲后者应该分配给各从进程局部存放，但由于输出的需要，仍然要在主进程中 全局存放象压力、饱和度和油气比等几个主要网格数组。这里所说的标量数据还包括一些 与网格数和井数无关的小数组，如相对渗透率和 PVT 关系表数据等，其占用的内存是微不 足道的，所以可在各进程中重复存放。 3.1 数据分区数据分区 由于直井一般会占大多数，为了不必要地增加编码的复杂度，只要允许用户在 X、Y 两个方向上任意划分区域，就足以避开区域分界线与斜井、水平井或断层相交的情况。得 到用户的分区信息后，作少量的编码易得满足上述条件的分区方案，当然应尽可能使各区 域内的局部网格数接近，这有利于加载平衡。区域的编号可任意。考虑到区域间自动信息 交换的要求，区域边界线的编号应先一个方向，再另一个方向。分区后应准确存储续用信 息，如区域、区域边界线的位置、尺度及编号等。另外，区域内的井也应再进行编号，且 要保留域内编号与原编号的关系。 3.2 内存的管理内存的管理 主进程启动后，读重启文件，按记录要求申请足够的内存，读入全部初始信息。待把 应该由各子域存储的数据分别传过去后，备份必要数据，释放并重新申请动态内存，再导 入备份数据。这样可使主进程占用更少的内存。从进程则按各自需要分别申请动态内存。4雅可比矩阵计算的并行雅可比矩阵计算的并行雅可比矩阵的计算是可精确并行化的。考虑到边界网格方程的雅可比元素与相邻子域46对应的边界网格数据有关，此阶段我们不妨把各子域的网格系统在内边界处向外扩一排， 使当前的虚拟网格系统包括与该子域雅可比矩阵计算有关的所有网格。若用 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。 这样安排避免了雅可比计算过程中信息传递的频繁发生，但比串行增加了总的计算量， 不过只要在代码中稍加逻辑判断即可消除。此处设置一个局域整型数组是必要的，便于雅 可比计算完成后把结果从虚拟网格向实际网格映射，用于后续对线性系统的求解。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 的存储方案是，把矩阵中属于某子域的行都交由对应进程来123456图2 子域内雅可比计算的虚拟网格系统47存储。按 3.1 节所述的数据分区方案，由某一子域内节点方程和井方程导出的行的元素恰 好在对应的进程内，这就保证了矩阵的局部装配过程中不发生子域间数据的相互传递。 PETSc 的开发和释放者强烈地建议用户在矩阵的局部装配之前用 D_NNZ 和 O_NNZ 两个整型数组准确地告诉装配函数该矩阵的对角块和非对角块上非零元素的位置，使装配 过程中尽量减少内存申请次数和便于按行压缩存储，以此获得并行求解的高效率。解决这 一问题的办法是先装配，在确认装配无误后再逻辑跟踪装配过程，准确记录非零元素的位 置，形成须预先运行的子程序来设置 D_NNZ 和 O_NNZ 数组。实际求解时不必在每一次 之前都运行此程序，而只需在井层生产注入状态发生变化的某一时间步（含一次模拟运 算的第一时间步）的第一次牛顿迭代之前运行就可以了。 5.3 对对 PE