CN104136942B

CN104136942B - 用于大规模并行储层仿真的千兆单元的线性求解方法和装置

Info

Publication number: CN104136942B
Application number: CN201380008956.XA
Authority: CN
Inventors: 拉里·S·冯
Original assignee: Saudi Arabian Oil Co
Current assignee: Saudi Arabian Oil Co
Priority date: 2012-02-14
Filing date: 2013-02-14
Publication date: 2017-02-22
Anticipated expiration: 2033-02-14
Also published as: EP2815257A2; EP2815257B1; WO2013123108A2; CA2862954C; US20130211800A1; EP2815257B8; WO2013123108A3; CA2862954A1; US9208268B2; CN104136942A

Abstract

对储层数据应用线性求解方法，以采用结构化网格或非结构化网格以最低限度粗化的方式对由巨型油田的高分辨率储层仿真得到的大方程组进行求解。可以将精确仿真结果所需的分辨率下的全地质结构复杂度和不连续性考虑在内。提供了一种通用的非结构化方法，使得可以对多分支井附近的非常复杂的流动几何结构进行建模。

Description

用于大规模并行储层仿真的千兆单元的线性求解方法和装置

相关申请的交叉引用

本申请要求于2012年2月14日提交的美国临时申请No.61/598,633的优先权。为了美国专利实践的目的，该临时申请的内容通过引用整体合并于本申请。

技术领域

本发明涉及被称为巨型储层的计算机仿真。

背景技术

在石油和天然气产业中，大的地下油气储层的开发通常需要建立计算机仿真模型。这些地下油气储层通常为复杂的岩层，其包含以两种或多种流体相存在的石油流体混合物和水的形式的流体。在生产过程中，通过被钻入这些岩层并且在这些岩层中完成的井来生产石油混合物。有时，诸如水和/或气体之类的流体也被注入进这些岩层以提高石油流体的采收率。石油和天然气公司开始依靠储层仿真来作为提高它们开采石油储层的能力的重要工具。

储层仿真属于一般的多孔介质中的流体仿真的领域。但是，储层仿真通常涉及处于高压和高温下的地下地质构造中的多种油气成分和多个流体相。这些仿真器必须将这些油气流体和所包含的地下水的化学相行为考虑在内。

仿真模型包含描述了岩层和井的具体几何结构的数据、流体和岩石性质数据、以及与所关注的油气田中的特定储层相关的生产和注入历史。通过仿真器(被称为储层仿真器)形成仿真模型，该仿真器为一系列在数据处理系统上运行的计算机程序。运行这些模型的储层仿真器为计算机实施的数值方法、或基础数学模型的编码算法和数据构建。反映这些油气储层中流体运动的物理特性的数学模型为非线性偏微分方程组，该非线性偏微分方程组描述了：在这些储层中由流体生产和/或注入引起的瞬态多相、多成分流体流动和物料平衡行为，以及储层流体的压力-体积-温度(PVT)关系。

储层仿真器通过将体积细分为被称为网格块的连续单元以对地下储层以及所包括的周围的多孔岩层内的多相多成分流体流动和物料平衡进行仿真。网格块是应用基础数学模型的基本有限体积。基于仿真所需分辨率和所关注的储层的尺寸来改变网格块的数量。对于大储层，例如在业内被称为巨型储层的类型，其可以具有数十亿桶原始原油地质储量(original oil-in-place，OOIP)，因此其网格块的数量可以为数亿到超过十亿的数量级，从而具有足够的分辨率来反映储层内的流体动态、地层岩石孔隙度和渗透率非均质性、以及许多其它地质和沉积的复杂性。这种尺寸的储层的仿真可以被称为“千兆单元储层仿真”。

根据储层分层的沉积历史和随后对其塑形的地质或侵蚀过程，岩石分层的几何结构和由此对网格单元的描述将会非常复杂。图1为典型的采用角点几何(corner-point-geometry，CPG)描述的结构化网格G的示意图。对网格单元的几何描述通常包含诸如断层和裂缝、尖灭(pinch-out)和页岩壁垒等不连续性。这些不连续性导致即使对于采用结构化网格的储层仿真模型也需要非结构化单元连接，例如如图1中所示。这些不规则被称为非相邻连接。

地质复杂性是使储层仿真采用可以更精确地表示流动几何特性的非结构化单元连接的一个促进因素。图2A为具有多分支井20(分别在储层中的各自位置处被指出)的关注储层的被称为非结构化PEBI网格U的示例的一部分的显示。图2B为图2A中由参考数字22指出的部分的放大视图。由图2B中可以看出，非结构化网格U中的多分支井20周围的网格单元的离散化遵循局部流动几何特性，而不是如图1的结构化网格G的情况那样结构化的和直线的。

由于如今所钻的典型的井为战略性地布置为使储层接触面和油气采收率最大化的水平井或多分支井，因此这种复杂的井几何结构也是另一个促进因素。诸如图2A和2B中的网格之类的在井附近遵循局部流动方向的网格提供更好的数值精度。现代储层仿真器要求求解方法能够解决不规则的单元间连接。从而导致对解决包含一般非结构化矩阵的稀疏线性系统的需求。

精确的井周边流动建模和交叉流动的建模是现代储层仿真的另一个重要方面。所采用的求解器必须包括储层系统的全隐全耦合(fully-implicit fully-coupled)井建模以使得仿真是精确的和鲁棒的(robust)。全耦合式求解在大量并行计算的环境下还必须是快速和可扩展的。

多相多成分方程组的瞬态解涉及从储层的初始状态开始按照时间步长序列进行的质量和能量守恒的演变。对于每个时间步长，每个有限体积的非线性离散方程组均采用公知的广义牛顿法来进行线性化。成分i的一般组分守恒方程为：

其中

如果忽略弥散、化学反应和吸收，该组分方程简单为：

由于多孔介质的孔隙必定填充有流体，因此孔体积必定等于总的流体体积。这可以表达为：

其中孔体积V_φ仅为压力的函数，描述为：

压力和总摩尔数为主要变量。为了方程组的封闭性(closure)，其他采用的方程如下：

对每个相的摩尔分数的约束条件：

对每种成分的总摩尔数的约束条件：

对流体饱和度的约束条件：

其中：

并且

通过达西定律(Darcy’s Law)描述相速度：

其中K为渗透率张量，定义为：

上述方程中的符号的含义如下：

在业内，这被称为牛顿迭代或非线性迭代。在每次牛顿迭代中，线性方程组被构造为采用被称为雅可比矩阵的矩阵和被称为残差的右侧向量来对系统主要变量的变化进行求解。

现今用于求解线性方程组的行业惯例是通过预条件迭代(preconditionediterative)方法。该迭代算法可以为以下三种形式中的一种：一种称为ORTHOMIN；第二种称为GMRES；或第三种称为BICGSTAB方法。对于鲁棒性，所使用的预条件子(preconditioner)(而不是迭代器的选择)已更加重要。在储层仿真行业中采用的典型的最新技术的预条件子可以为一级预条件子，例如嵌套分解(nested factorization，NF)方法、或不完全LU分解(ILU)方法的变形。这些均广泛应用在现今的商用仿真器中。

较新的实践包括采用多级预条件处理(multi-stage preconditioning)。储层仿真中的控制方程的本质为抛物线和双曲面的混合。对于方程组的抛物线部分的处理，如果原始问题的压力成分首先通过被称为约束压力残差(CPR)的算法从整个方程组中解耦，那么多级代数多重网格(AMG)或几何多重网格(GMG)可以非常有效。近似压力问题可以有效地通过多级方法(例如AMG或GMG)来处理。残差的剩余双曲面成分可以采用合适并且成本更低的一级预条件子来容易地处理。这实际上是一种在本领域内公知的被称为“分治”预条件处理策略。但是，到目前所知，GMG多级方法不能处理非相邻连接。因此，AMG方法是通常用于解决现实仿真问题的通用多级方法。

在并行计算中计算处理算法的效率的度量为其可扩展性(scalability)。如果一种方法在具有单个CPU的计算机上解决一个计算问题需要一小时，而如果将该工作准确地划分并分配给两个CPU，解决同一问题的时间是0.5小时，并且使用四个CPU解决同一问题的时间是0.25小时，并可以此类推，那么该方法是完全可扩展的、或是具有100％的并行效率。即，没有更好的并行化了。完全可扩展方法会是理想的情况。现实中，许多原因可以导致实际系统的求解时间远长于该理想的情况。

以上描述的实践或策略很好地适用于利用单个CPU的单核进行串行计算、以及采用数十个CPU核进行小规模并行计算。但是，当处理数量增加时，它们会遭受并行可扩展性的严重损失。这是因为上述嵌套分解(NF)或不完全LU分解(ILU)具有大量的不容易实现并行化的递归成分。此外，诸如代数多重网格(AMG)之类的多级方法会在粗化的网格级的聚合和平滑中引发显著的通信开销，从而降低并行可扩展性。

储层和油/气田的仿真模型变得日益庞大和复杂。与此同时，快速发展的计算机硬件也变得便宜并且快速。现今硬件具有向大规模并行化发展的趋势：从过去的具有单核CPU的高性能计算或HPC系统，到现今采用的由数千个多核CPU组成的PC集群，到最近出现的异构性CPU-GPU系统。仿真实施者现在可以采用便宜的具有数万个计算机核心共同工作以解决大规模仿真问题的硬件。

就目前所知，现有的求解方法由于较差的可扩展性问题和/或鲁棒性问题而不能处理千兆的储层仿真。当前可用的求解方法限至于约千万网格单元和数十个处理的并行可扩展性。如果没有大规模并行可扩展性，那么针对历史匹配仿真或预测仿真运行的仿真将会变得很慢且难以实施。这限制了可以在流体仿真中采用的模型的大小。对于大储层，这意味着过度粗化。

关于储层的一些其他因素还进一步提高了千兆级储层仿真的复杂度。对于常见的诸如断层、裂缝、尖灭和页岩壁垒之类的地质不连续性的处理对形成鲁棒的可扩展求解方法带来了困难。采用结构化的网格难以对多分支井附近的复杂流体几何结构进行建模。非结构化的网格需要专门的非结构化方法，而且使得计算负荷难以在处理器之间平衡和难以鲁棒并且有效地求解。

而且，全耦合全隐的井解决方案是精确的储层仿真所需要的方法。这是可以降低可扩展性的另一个因素，特别是当储层包含数千个复杂的分支井时。此外，许多石油储层岩石具有多模式的多孔隙-渗透(porosity-permeability)和多尺度断层网络。

发明内容

简而言之，本发明提供一种新的和改进的计算机实施方法以利用储层仿真器对地下储层的流体特性进行仿真。流体特性取决于流体流动、物料平衡和压力-体积-温度关系。储层被划分为许多个子域，每个子域代表了储层的包含多相多成分流体的有限体积并且被组织成一组具有流体相关数据的网格单元。储层仿真在由至少一个主节点和多个处理器节点形成的集群计算机中执行，其中对每个处理器节点分配给所划分的储层的一个子域。根据本发明的计算机实施方法通过牛顿迭代对网格单元的关注流体特性参数进行仿真，该计算机实施方法由一系列计算机实施步骤组成。每个储层子域的数据被排列成牛顿矩阵A，该牛顿矩阵A由以下部分组成：矩阵P，其由子域的块对角线和与块对角线相邻的非零数据块组成；矩阵E，其由未包含在块对角矩阵P中的非零数据块组成。矩阵E由子域内的内部网格块的矩阵E’和具有储层中的相邻子域的网格单元的边界网格单元块的矩阵E″组成。针对每个子域的网格块在所分配的处理器节点中采用为子域选择的局部预条件子来执行流体特性的并行近似求解。基于求得的并行近似解来更新储层的牛顿矩阵的全系统残差。根据全局约化空间矩阵将该全系统残差限制在全局约化空间中。采用全局约化空间预条件子来执行对全系统残差的并行近似求解。随后，在分配的处理器节点中，由每个子域的并行近似解和全系统残差的并行近似解形成组合近似解的更新，以形成解向量。采用解向量更新全系统残差。对全系统应用迭代求解器以获得储层的网格单元的流体特性的测量值。对通过应用迭代求解器获得的流体特性的测量值进行检查，以确定是否实现收敛。如果未实现收敛，则返回至执行并行近似求解的步骤。如果实现了收敛，则处理返回至储层仿真器。

本发明还提供一种新的和改进的数据处理系统，用于利用储层仿真器对地下储层的流体特性进行仿真。该数据处理系统包括由至少一个主节点和多个处理器节点形成的集群计算机。储层的流体特性取决于流体流动、物料平衡和压力-体积-温度关系，并且该储层被划分为多个子域。每个子域代表了储层的包含多相多成分流体的有限体积并且被组织成一组具有流体相关数据的网格单元。对数据集群计算机的每个处理器节点分配所划分的储层的一个子域。集群计算机通过牛顿迭代对网格单元的关注流体特性参数进行仿真。数据处理系统中的处理器通过以计算机实施步骤的处理序列进行操作。每个储层子域的数据被排列成牛顿矩阵A，该牛顿矩阵A由以下部分组成：矩阵P，其由子域的块对角线和与块对角线相邻的非零数据块组成；矩阵E，其由未包含在块对角矩阵P中的非零数据块组成。矩阵E由子域内的内部网格块的矩阵E’和具有储层中的相邻子域的网格单元的边界网格单元块的矩阵E″组成。针对每个子域的网格块在所分配的处理器节点中采用为子域选择的局部预条件子来执行流体特性的并行近似求解。基于求得的并行近似解来更新储层的牛顿矩阵的全系统残差。根据全局约化空间矩阵将该全系统残差限制在全局约化空间中。采用全局约化空间预条件子来执行全系统残差的并行近似求解。随后，在分配的处理器节点中，由每个子域的并行近似解和全系统残差的并行近似解形成组合近似解的更新，以形成解向量。采用解向量更新全系统残差。对全系统应用迭代求解器以获得储层的网格单元的流体特性的测量值。对通过应用迭代求解器获得的流体特性的测量值进行检查以确定是否实现收敛。如果未实现收敛，则返回至执行并行近似求解的步骤。如果实现了收敛，则返回至储层仿真器。

本发明还提供一种新的和改进的数据存储装置，该数据存储装置具有存储在计算机可读介质中的计算机可操作指令，该指令用于使储层仿真器对地下储层的流体特性进行仿真。流体特性取决于流体流动、物料平衡和压力-体积-温度关系。储层被划分为许多个子域，每个子域代表储层的包含多相多成分流体的有限体积并且被组织成一组具有流体相关数据的网格单元。该数据处理系统包括由至少一个主节点和多个处理器节点形成的集群计算机，对每个处理器节点分配所划分的储层的一个子域。该数据处理系统还通过牛顿迭代对网格单元的关注流体特性参数进行仿真。数据存储装置包括用于使数据处理系统按以下处理序列进行操作的指令。将每个储层子域的数据排列成牛顿矩阵A，该牛顿矩阵A由以下部分组成：矩阵P，其由子域的块对角线和与块对角线相邻的非零数据块组成；矩阵E，其由未包含在块对角线矩阵P中的非零数据块组成。矩阵E由子域内的内部网格块的矩阵E’和具有储层中的相邻子域的网格单元的边界网格单元块的矩阵E″组成。针对每个子域的网格块在所分配的处理器节点中采用为子域选择的局部预条件子来执行流体特性的并行近似求解。基于并行近似求解来更新储层的牛顿矩阵的全系统残差。根据全局约化空间矩阵将该全系统残差限制在全局约化空间中。采用全局约化空间预条件子来执行对全系统残差的并行近似求解。随后，在所分配的处理器节点中，由每个子域的并行近似解和全系统残差的并行近似解形成组合近似解的更新，以形成解向量。采用解向量更新全系统残差。对全系统应用迭代求解器以获得储层的网格单元的流体特性的测量值。对通过应用迭代求解器获得的流体特性的测量值进行检查以确定是否实现收敛。如果未实现收敛，则返回至执行并行近似求解的步骤。如果实现了收敛，则返回至储层仿真器。

附图说明

图1为地下储层的典型结构化网格模型的等距视图。

图2A和图2B为地下储层的典型非结构化网格的示例的显示。

图3为根据本发明的用于处理的矩阵数据的示例划分的示意图。

图4A、图4B和图4C为根据本发明的矩阵数据域的划分和组织的示意图。

图5为根据本发明的用于储层仿真的数据处理系统中执行的一组数据处理步骤的功能框图。

图6为根据本发明的用于储层仿真的数据处理系统的示意框图。

图7为根据本发明的储层仿真的求解器可扩展性的曲线图。

图8为根据本发明的作为对五种仿真模型运行的仿真的总运行时间的百分比的仿真求解时间的示意图。

具体实施方式

利用本发明，提供了高度可扩展的大规模并行储层仿真中的多级预条件处理。通过本发明的储层仿真而求解的控制非线性偏微分方程是针对互连的有限体积(称为网格单元)的瞬态多相多成分流体流动而写出的离散质量、动量和能量的平衡。对于被称为巨型储层的储层，网格单元的数量可以为数百万直至超过十亿。此外，一个仿真模型可以包括数个储层。这些储层可以通过周围地质构造和地层中的连通的含水层和/或地下断层系统以液压连接。这些储层还可以通过钻入不同地层以进行多层合采或注入处理的多分支井来连接。

尽管针对只具有数百万个网格单元的模型的处理方法是有竞争力的、鲁棒的、并且提供出色的可扩展性，但本发明特别以需要鲁棒并且高度可扩展的预条件子来有效地求解问题的大型系统求解方法为目标。就目前所知，现有技术并未充分解决并行可扩展性需求的问题，由此对于大规模应用而言遭受了并行性能的严重损失。

一种被称为通用牛顿-拉普森(Newton-Raphson)方法的技术被用于产生形式为[A]{v}＝{R}的系统方程，其中[A]代表系统雅可比矩阵，{v}代表更新的解向量并且{R}代表残差向量。本发明涉及在系统矩阵[A]和残差向量{R}产生之后的求解方法。具有数年历史匹配和/或未来性能预测的单个仿真通常需要数百个至数千个时间步长，每个时间步长需要若干次牛顿-拉普森迭代，并且每次牛顿-拉普森迭代将需要若干次求解迭代以将系统收敛到可接受的容差范围之内。

为了求解在由数百个计算节点(每个计算节点都由多个CPU形成，并且每个CPU由多个核构成)组成的高性能计算(HPC)集群上的大规模仿真模型，需要将这些模型细分为被称为子域网格单元的组或块。每个子域都由解域的连续网格单元的子集组成，并且每个子域执行几乎等量的计算工作。对于较大的模型，可能需要将域划分为数千个或数万个子域。对子域数据的并行计算即是分治策略。

并行处理中的问题在于被建模的物理现象具有全局影响，而细分为子域的处理并不将该储层现象考虑在内。就目前所知，过去的细分会导致大量的计算或存储负担。因此，致使分治策略不能有效地减少处理时间。

成功的迭代求解器需要鲁棒的预条件子，许多鲁棒的预条件子具有大量串行或递归成分，这些串行和递归成分不容易修正为并行实施。另一方面，已知的高度并行的预条件子不会是非常鲁棒的。如果不能适当地解决模型的全局性质，则由于将解空间并行划分为子域，导致计算开销显著增加，从而使得并行计算不能有效地加速求解。

如上文所述，本发明提供了多级求解方法，该方法具有串行方法的鲁棒性，但同时还达到了大规模并行效率。为了方便参考，本发明提供的方法被称为全局约化空间预条件方法。该全局约化空间预条件方法可以采用许多已知的鲁棒的串行方法作为局部成分，以构建可修正为用于在数千个处理节点上运行的实施方式的鲁棒的可扩展并行方法。

为了实现可移植性，本发明可以实施为采用诸如FORTRAN、C或C++之类的标准高级计算语言的计算机代码以及采用用于分布式存储器的MPI标准、用于共享存储器的OpenMP标准或混合范式实施方法的并行化。这允许实施的系统容易地在多种HPC硬件平台上运行。既然本发明在存在数百万个困难的储层仿真要解决时展示出了高效率，涉及数十亿个困难方程的系统也会通过本发明而尤其受益。

对于由大量有限体积(其中每个有限体积可能具有一大组控制方程)组成的仿真系统，需要将这些有限体积分配给HPC集群的多个计算节点，从而拥有足够的资源(计算核和存储器)来有效地求解系统。如果意图采用N个进程(采用索引i来确定进程的等级)来求解全局问题，那么有限体积的全局数量和相关方程组分布至N个进程并且每个进程i处理全局体积的子集和方程组：

[A]_i{υ}_i＝{R}_i (1)

i＝1，2，3…N

这种划分为子集的方法被称为域划分(domain partitioning)。通过划分得到的每个子域包含大致等量的计算工作，使得工作负荷在每个进程之间大致平衡以得到好的并行可扩展性。对超大仿真系统的划分同样需要每个子域的工作量大致相等，并且同时最小化与相邻的域相通和从相邻的域获得数据的要求。对于涉及单个物理特性并针对每个有限体积具有相同的守恒方程组的仿真模型，域的划分需要每个子域包含相同数量的网格单元并且最小化在子域边界处的网格单元的数量。

由于储层仿真系统紧密耦合，为了具有并行可扩展性，需要全局解，并且为了调和解空间的全局性质，有效解通常要求计算处理之间的大量数据通信。但是，对于大规模的并行化，先前认为有必要增加计算处理的数量以减少仿真时间。通信成本是并行化开销的组成部分，其最终由于子域尺寸的降低而导致可扩展性的损失。

由此发生这样一种情况：计算处理的进一步增加并不会导致仿真时间的降低。一些已知的用于储层仿真的鲁棒的预条件方法通过依赖于顺序的执行来递归，这导致并行可扩展性的范围受到限制。嵌套分解(NF)是属于该类别的一个示例。另一方面，高度并行预条件方法并不是鲁棒的并且需要过多的迭代才能收敛，或者对于困难的问题甚至可能不能收敛。多色块高斯-赛德尔(multi-color block Gauss-Seidel，MCBGS)或块雅可比预条件子是属于该类别的示例。本发明的方法提供了能克服这些困难的多级方法。

图5的流程图F结合图6所示的数据处理系统P示出了用于根据本发明的用于大规模并行储层仿真的本发明的基本计算机处理序列。在步骤50期间，针对储层仿真系统的独立子域选择和构建局部预条件子。根据特定子域的仿真问题的本质，针对一个子域构建或选择的预条件子可能不同于其他子域的预条件子。

接下来在步骤52(图5)期间构建或建立本发明的基础并行预条件子。如图3所示，步骤52的预条件子涉及将储层的每个子域的网格单元矩阵子结构化(matrix sub-structuring)为核心子矩阵[P]和补充子矩阵[E],其中[A]＝[P]+[E]。矩阵[P]由块对角线和该对角线附近的非零块组成。重新排列子域内的网格单元使得最大的连接因子与对角线最接近。这种排序策略被称为最大传输率排序(maximum transmissibility ordering)。

每个块行的宽度根据该行上的包含在子域的矩阵[p]中的非零块的数量而改变。子域的[E]矩阵由未包含在[P]矩阵中的非零块组成。[E]矩阵包括子域内的内部网格单元连接的部分[E′]和具有相邻的子域中的网格单元的所有边界网格单元连接[E″]。

由于系统矩阵[A]为并行存储器分布，因此在图5的步骤52期间可将矩阵数据存储为针对[P]、[E′]、[E″]部分中的每一个部分的块压缩稀疏行(CSR)的3个列表。步骤52的预条件子可以被称为变量块行-线求解幂级数(Variable-Block-Row-Line-Solve-Power -Series，VBR-LSPS)预条件子，并且具有以下数学形式：

步骤52的VBR-LSPS预条件子为基于方程2中的截尾级数的近似逆预条件子，其中N是为预条件子保留的级数项，并且无论何时需要{y}＝[P]^-1{x}的结果，通过首先在设置步骤和向前/向后代入步骤中生成LU分解来计算[P]^-1。

图3中示出了四个子域的子结构化的系统矩阵[A]的作为步骤50和52的结果的数据结构。并行数据被看作分布式存储器，每个域数据存储在独立的存储器分区中。

图4A为示出了本发明的全局约化空间求解方法的基本概念的示意图。示出了被划分为九个子域40的解域S，其中每个子域i包括几乎相同的计算工作量。如果所有子域均采用相同的局部方法和相同的方程组，意味着将有限体积等量细分至每个子域。

根据本发明的方法基于局部化步骤54(图5)，该局部化步骤54用于对采用由步骤50选择的局部预条件子指示的每个子域I进行并行近似求解。子域I的矩阵方程(1)可以展开为：

其中，A_I代表对应于子域i的内部网格单元的矩阵系数；A_B代表对应于子域i的边界网格单元的矩阵系数；B和C为内部网格单元和边界网格单元之间的连接系数；以及E为子域i内的边界有限体积和子域i的所有相邻子域的边界有限体积之间的连接系数。局部化步骤54从矩阵方程(1)中移除E和v_E以形成独立集，对于子域i该独立集可以被展开为：

现在可以在步骤54期间针对每个子域在单个处理上采用没有任何进程间通信的单线程或多线程来独立地求解属于计算处理的子域40的子域(40)的独立集。

图4B中示意性地示出了单独子域40的独立并行处理，其中示出的子域40彼此分隔。对单独子域的独立处理在涉及多物理特性的非均质系统中特别有用，在这种非均质系统中每个子域可以在求解的方程中具有不同的刚度等级，并且针对一个子域选择的局部化方法可能需要比其他子域选择的方法更强。

在本发明中，通过基于步骤50的合适的局部预条件方法，在步骤54期间针对每个单独子域的独立集来确定并行独立近似解。在分组中的每个处理可以彼此独立地实行该步骤并且在有需要时可以采用不同的预条件方法。有许多可以在步骤54期间采用的具有不同优点和计算成本的鲁棒的预条件方法。应当理解，这种局部预条件方法也可以为多线程的，从而采用共享存储器并行性以提高性能。

举例而言，如果用于仿真的控制方程为椭圆型方程，那么代数多重网格(AMG)预条件子是合适的。如果储层仿真的问题为双曲型方程的形式，那么诸如ILU(k)方法或ILUT方法之类的不完全LU分解是合适的。对于直接7点结构化矩阵(straightforward7-pointstructured matrix)，NF是快速和鲁棒的。如果子域包含非常刚性的方程(即，其数值解通常不稳定的方程)，那么可能需要对这些方程采用具有局部旋转的稀疏完全LU分解。

对于额外的预条件子的研究，例如参见：Yousef Saad,“Iterative Methods forSparse Linear Systems”,2003,SIAM,ISBN0-89871-534-2。由于递归算法和存储器通信要求，许多这些方法对于良好的并行化性能具有有限的范围。但是，对于在单个计算处理上针对每个单独子域执行步骤54的局部预条件处理，这些方法非常适合。

如在多孔介质中的很多多相多成分流动和传输问题一样，对于涉及本质上为混合的抛物线-双曲面的方程的问题，可以利用约束压力残差(CPR)方法。CPR方法将抛物线部分从全系统中解耦，并且采用适合的预条件子来产生椭圆变量的近似解。

随后在步骤54期间采用该近似解来预设用于求解全系统耦合的双曲抛物问题的全系统预条件子。

CPR算法可以写成以下方程形式：

在步骤54期间采用CPR算法对压力矩阵的解耦允许在步骤54期间对压力变量使用AMG多级求解器。在产生近似压力解时采用全局约化空间方法以达到近理想的可扩展性。全系统预条件子可以为高度可扩展并且鲁棒的VBR-LSPS预条件子。或者，该预条件子可以为与用于可扩展全系统预条件处理的全局约化空间方法相结合的其他局部预条件子(例如ILU(k)方法)。

通常，局部预条件处理步骤54可以表达为：

其中选择使其可以容易地根据A和条件数求解出来。选择预条件处理以正确地保存初始残差并且方便应用通常基于相对残差范数的收敛性判据。相对残差范数为当前迭代的残差向量的L2范数与初始的残差向量的L2范数的比值。

在步骤54期间局部化问题的近似解可以写为：

对于子域40，在方程7中表达的近似解在子域边界(如图4C中的40b所示)处具有显著的剩余误差。随后采用在步骤54期间获得的方程7的子域解执行步骤56以计算在步骤54期间获得的残差的全局更新。这可以写为：

步骤56涉及矩阵向量相乘的步骤，该矩阵为全矩阵而并非局部子集的约化矩阵。向量分量υ′_E属于相邻子域，并且必须由拥有这些子域的相邻进程传递。如果这些采用消息传递接口(MPI)标准来实施，那么其对应于点对点的非同步通信步骤。此外，可以容易地实施计算和通信的重叠以达到可扩展的并行性能。

在(多个)步骤56期间可以采用若干矩阵格式。本发明并未排除对矩阵存储方案的选择。在一个实施例中，分布的块稀疏矩阵被子结构化为三个块CSR存储。对于全局约化空间方法而言矩阵格式的选择是灵活的，并且可以进行修改以满足存储器访问的特定需要以及满足HPC硬件的存储需要。分布式并行矩阵数据存储对于大规模并行实施具有有效的结果。在2010年9月7日提交的题为“Machine,Computer Program Product and Method toGenerate Unstructured Grids and Carry Out Parallel Reservoir Simulation”的共同拥有的待审美国专利申请No.12/876,727中描述了该并行数据系统。

由于局部子域的解忽视了由所有子域的解空间之间的全局相互作用导致的某些贡献，因此方程8中的全局更新残差向量可以在子域边界附近集聚显著的残差。在步骤52中构造的并行全局约化空间预条件子现在可以被用于解决这些集聚的边界残差。在步骤58中的全局约化空间矩阵和约化空间残差可以构造为：

A_IR＝S^t·A_I·S (9a)

F＝S^t·C (9b)

D＝B·S (9c)

R′_IR＝S^t·R′_I (9d)其中S为约束算子(restriction operator)，S^t为插值算子。S和S^t的行为是将图4a中的全局精细空间中的问题转换为如图4C中的全局约化空间中的粗化问题。由于目的在于解决集聚的边界残差，因此仅在域的内部粗化空间。因此，这些行为是对两级之间的矩阵和残差的多网格精细网格至粗化网格的转移操作。

在最简单的形式中，约化系统为代表子域的内部有限体积40c的单个粗化单元。如图4C中所示，内部有限体积40c为与邻域中的有限体积40b没有显著连接的有限体积的子集。许多其他集结法可以用于产生粗化空间，并且适合本发明采用。聚合度较低的内部粗化的使用成本可能更高，但是其可以产生更好的近似。另一种策略为采用从子域边界的距离加权集结。这些都是根据本发明可以考虑采用的可行的选择。

根据步骤52、58和60的全局约化空间矩阵方程为以下形式：

在步骤60期间，采用高度并行全局预条件方法根据方程10来确定有效且鲁棒的并行近似解向量。其涉及计算处理之间的MPI通信，这类似于计算方程(8)中的全局残差的要求。在步骤60期间，全局约化空间预条件子通过将全局约化空间系统矩阵划分为以寻求近似解，其中包括对角块和中的其他大系数并且是可以有效地找到逆矩阵的部分，或者可以计算稀疏LU分解矩阵。

矩阵不包括所有在E中的系数。的构建可以包括子域中的网格单元的重新排列以产生局部最小带宽系统。随后，近似为：

方程(11)与方程(2)类似，但是其在计算机处理步骤60期间应用于全局约化空间系统。矩阵包括A_IR，D，F和A_B在的块行宽度内的部分。矩阵包括矩阵E和矩阵A_B中未包括在中的分量。对A_B进行划分的方法包括基于最大传播能力和随后[P]的宽度的行指定进行的局部重排序。

根据本发明，在首先通过计算将或分解为LU矩阵之后，本质上通过将一系列矩阵向量相乘的步骤来计算出近似全局约化空间解。在设定步骤52中执行分解步骤。在步骤50中的全局约化空间预条件有效地构造了子域边界解，由于步骤54中的局部子域解不能充分地求解子域边界残差，因此这种子域边结界是对局部子域解的补充。随后约化空间解插入回步骤62中的子域：

步骤64随后将来自步骤54的由方程(7)表达的局部近似解以及来自方程(12)的子域i的全局约化空间解组合在一起成为如下所示的全局近似解：

如方程(13)所示的在步骤64中执行的组合近似解的更新是在步骤66中用于以更新迭代求解器的全残差的解向量。在步骤66之后，步骤68应用预条件处理后的Krylov子空间迭代求解方法。其可以为已知的传统ORTHOMIN、GMRES、或BICGSTAB算法中的一个的并行实施。将在步骤54中执行的局部预条件处理和在步骤60中执行的全局预条件处理进行结合显著地加快了在步骤68中的迭代求解器的收敛并且解决了在一些Krylov子空间迭代内的收敛困难问题。

在步骤70中，对于特定的求解器迭代i(i仅为求解器迭代计数器)，检查由执行步骤68所得的残差以确定是否实现预定容差范围内的收敛。如果没有，处理返回至步骤54以继续按上述方式进行处理。如果在步骤70期间确定实现了预定容差范围内的收敛，那么如步骤72所示，数据处理的处理结果和控制均返回至储层仿真器。储层仿真器的非线性求解器采用解向量来更新其主要变量。随后对所有的辅助变量进行非线性更新并且计算非线性残差向量。如果主要变量的解未收敛到预定的物料平衡容差和变化之内，则调用下一个牛顿迭代。如果解已经收敛，则仿真器将前进至下一个时间步长。

本发明特别适合很大的紧密耦合的方程组，该方程组通常需要在现有的HPC集群上的许多计算节点来进行求解。超大的仿真模型可以形成高度不对称并且呈现非正定块稀疏矩阵的复杂方程组。不连续的系数可以横跨十六个数量级。

通常，迭代求解器是低效的。本发明的上述具有并行全局约化空间预条件子的方法产生可以一致地减少全局空间中的系统残差的全局近似解。本发明允许Krylov求解器在几乎相同次数的迭代内收敛，好像该系统没有以并行分布的方式细分到多个计算处理。由此该系统可以针对超大、难以解决的问题达到很高的并行效率。

现在考虑根据本发明的数据处理系统，如图6所示，针对根据本发明的计算机仿真提供数字处理系统P以用于大规模并行储层仿真。数据处理系统P包括一个或多个中央处理单元或CPU82。作为集群节点的(多个)CPU82具有与其相关联的用于以下数据的储层存储器或数据库84：通用输入参数、来自井的核心样本数据、单元组织数据和信息、以及数据处理结果。用户界面86可操作地连接至CPU82，该用户界面86包括用于显示图形图像的图形显示88、打印机或其他合适的图像形成机构以及用户输入装置90，从而为用户提供访问和操作的权限，并且提供处理结果、数据库记录或其他信息的输出形式。

储层存储器和数据库84通常在外部数据存储计算机94的存储器92中。数据库84包含包括模型中的单元的结构、位置和组织结构在内的数据以及在根据图5的处理方法的储层仿真中采用的数据通用输入参数、来自油的核心样本数据、单元组织数据和信息、以及数据处理结果。

数据处理系统P的CPU或计算机82包括主节点96和连接至主节点96的内部存储器98，该内部存储器98用于存储操作指令和控制信息、并且在需要时作为存储或转移缓冲区。数据处理系统P包括存储在存储器98中的程序代码100。根据本发明的程序代码100为计算机可操作指令的形式，使得主节点30来回传递数据和指令，该数据和指令用于通过处理器节点执行逐个单元地仿真储层中的单独单元的储层性质或属性的处理，将在下文中详细描述。

要注意，程序代码100可以为微代码、程序、例程、或是提供一组特定的有序操作以控制数据处理系统P的功能并指导其操作的符号化计算机可操作语言的形式。程序代码100的指令可以永久性地存储在存储器98中或在计算机磁盘、磁带、传统硬盘驱动器、电子只读存储器、光学存储设备、或其他合适的存储了永久性计算机可用介质的数据存储装置上。程序代码100也可以包含在作为永久性计算机可读介质的数据存储装置中。

处理器节点102为通用可编程数据处理单元，其被编程为执行上述数据岩石物理算法的处理并且逐个单元地仿真储层中的单独单元的储层性质或属性。处理器节点102在(多个)主节点82的控制下操作。

在512节点(1024个CPU，每个CPU6核)的PC集群上运行的仿真产生了计算结果。报告的仿真结果用于针对分布式存储器并行化的采用FORTRAN90/95和MPI标准的实施。这些CPU为在2.93GHz下运行的英特尔Westmere CPU(Xeon-X5670)。互联结构为QLoigic公司的改进的QDR无限带宽(Infiniband)交换机。但是，可以理解的是，还可以采用其他计算机硬件。采用FORTRAN90/95+MPI标准通信接口，可以采用许多支持这些广泛接受的高级语言和标准的HPC硬件。

对于大规模并行应用，方法的可扩展性对于根据计算节点的数量的增加而成比例地加速仿真时间非常重要。对具有28,664,550个网格单元(20,355,022个活动单元)的全油田模型的仿真运行被用于示出本发明的并行可扩展性。该模型具有3000口油井，并且储层类型为断裂的三倍空隙三倍渗透率(triple-porosity triple-permeability，TPTP)储层。

图7示出了240个CPU核至1920个CPU核之间的可扩展性曲线图。图7为用于确定为以下表1中的模型3的2870万单元模型(2040万个活动单元)的求解器可扩展性曲线图，该模型为具有3相黑油(BO)流体描述并且采用全隐式公式的断裂的三倍空隙三倍渗透率的模型。图7中的检验结果显示出采用了多达1920个CPU核，本发明保持了超线性的可扩展性。在采用1920个核时，仿真的总挂钟时间(wall-clock time)为1.324小时，该时间对于实用目的而言是完全可接受的。

现有技术的仿真处理通常会遭受可扩展性损失并且不能实际地求解超过千万单元的储层仿真问题。举例而言，如果仿真处理方法可以达到90％的合适的并行效率，那么在不考虑还需添加多少并行硬件时，期望最高提速为十倍。由于求解器可扩展性是最主要的难点，因此在仿真领域中所熟知的是随着问题大小的增加，求解器所使用的总运行时间的百分比可达95+％。即，CPU时间完全由求解时间所决定。

下文的表1示出了根据本发明的方法的五个仿真运行的统计以及这些运行的求解时间百分比。在表1中，BO表示黑油模型，COMP 表示成分(compositional)，SPSP表示单倍孔隙单倍渗透率，DPDP表示双倍孔隙双倍渗透率，TPTP表示三倍孔隙三倍渗透率，并且Nc为油气中成分的数量。

表1

名称	模型1	模型2	模型3	模型4	模型5
						#网格单元	1,031,923,800	171,987,300	28,664,550	12,293,424	12,293,424
#活动单元	1,031,918,940	171,987,300	20,355,022	5,941,741	5,940,515
						#井	2,959	2959	2959	557	623
模型类型	SPSP	SPSP	断裂TPTP	断裂DPDP	断裂DPDP
						流体类型	3相BO	3相BO	3相BO	COMP,Nc＝9	3相BO
Duration	60	60	60	52	42
						#方程/单元	3	3	3	11	3
#CPU核	4,800	3000	2160	1,200	1,200
						总挂钟时间	14.464hrs	6.646hrs	1.324hrs	6.736hrs	2.935hrs.
求解器％	58.78％	51.99％	50.90％	54.15％	53.25％
						#时间步长	3,074	3,461	2423	6,981	5,375
#求解器迭代	16,778	23,443	21824	31,715	78,381
						时间/迭代/单元-核	8.484μsec	9.255μsec	11.79μsec	83.61μsec	14.5μsec

可以看到模型之间的差异极大：尺寸(从模型4和5的六百万个单元到模型1的超过十亿个单元)；储层类型；以及流体类型，但求解器百分比保持在50％至60％的范围中。求解时间并未增加至超过现有技术中的总运行时间。图8中示出了表1中的各情况的求解器百分比。

由上文可以看出本发明的处理方法提供了大规模的并行化，还提供了与现有技术中的强串行方法相同的鲁棒特性。本发明可以通过数年的储层历史和预测仿真来求解非常困难的十亿单元的多相多成分流体流动模型。本发明可以采用数千个CPU核鲁棒地以近完美的可扩展性来求解大范围的仿真问题。

被编程用于分布式存储器的储层仿真器可以采用本发明作为求解引擎。仿真器可以为采用结构化网格(Catesian或CPG)、非结构化网格(PEBI)、或其两者的结合的类型。仿真器可以包括诸如废弃单元(dead cell)、尖灭、断层和裂缝之类的地质复杂性。储层描述可以是单倍孔隙单倍渗透率(SPSP)至多倍孔隙多倍渗透率 (MPMP)。MPMP系统的逐单元的指定是最常见的并且包括作为系统的具体示例的已知的双倍孔隙(DP)和双倍孔隙双倍渗透率(DPDP)系统。储层流体可以为多相多成分的。

根据本发明，全局约化空间求解的目标在于子域边界修正。该程序在大规模并行实施中是可扩展的。全局约化空间方法为将解空间集中在子域边界的网格粗化方法。全局约化空间的近似解可以采用高度可扩展并且依然鲁棒的VBR-PSLS预条件子来获得。

通过选择鲁棒的局部预条件子来实现域内部预条件处理，该局部预条件子例如是用于压力的AMG、或用于通用预条件处理的ILU变型。此处可以基于解空间的难度和性质来对一些预条件子进行选择。当独立应用时，这些局部预条件子能够有效地减少子域内部的残差。如果不采用全局约化空间方法，则采用这些局部预条件子会留下显著的边界残差，导致过多的求解迭代。

本发明有利地解决了在鲁棒的串行预条件子(其可以包含限制并行化范围的算法)的大规模并行应用中的可扩展性问题。本发明有效地构建了子域边界解。随后将全局约化解和局部解组合。组合的解和对应更新的残差在加速储层仿真中通常使用的Krylov迭代求解器的收敛方面成为有效的估计。

本发明提供了一种针对多孔介质储层仿真问题中的多相多成分流体流动实现鲁棒且高度可扩展的解的独特且灵活的方法。该仿真器可以具有数百万至数十亿的范围内的联立方程组。对于现有的求解方法在方程组达到约千万时即不可用，其原因在于可扩展性的损失和/或鲁棒性的损耗导致过多迭代并且由此导致计算成本太高。本发明采用数千个计算核来实现大规模并行可扩展性，从而有效地求解这种耦合系统。该系统保留了近完美的可扩展性并且迭代次数并未增加。

从上述说明可以看出本发明具有近完美的可扩展性，并且可以通过简单地增加用于求解问题的CPU核的数量来缩短中等大小的仿真模型的运行时间。由于该求解方法不存在对仿真模型大小的理论限制，因此可以构建千兆单元或十亿单元模型来获得具有最小粗化率的高分辨率模型。该方法可以通过具有足够多的CPU核和存储器的商品化的HPC硬件来求解以获得好的计算运转时间。由此本发明可以用于克服模型的尺寸限制和速度限制。

本发明已经被充分描述以使得在本领域中的普通技术人员可以重现和获得在本文中提到的发明结果。但是，本发明的技术、主题所属的领域中的技术人员可以实施本文中未描述的修改以将这些修改应用至确定的处理方法，或者在使用本发明的结果时需要以下权利要求中要求保护的主题；这些修改均涵盖在本发明的保护范围之内。

应注意和理解，可以在不脱离所附权利要求书中阐述的本发明的精神和范围的情况下，对以上详细描述的本发明做出改进和修改。

Claims

1.一种计算机实施方法，通过牛顿迭代对网格单元的关注流体特性参数进行仿真，其中在计算机实施的利用储层仿真器对地下储层的流体特性的仿真中，所述流体特性取决于流体流动、物料平衡和压力-体积-温度关系，所述储层被划分为多个子域，每个子域代表了所述储层的包含多相多成分流体的有限体积并且被组织成一组具有流体相关数据的网格单元，所述储层仿真在由至少一个主节点和多个处理器节点形成的集群计算机中执行，其中对每个处理器节点分配所划分的储层的一个子域，所述计算机实施方法包括步骤：

将每个储层子域的数据排列成雅可比矩阵A，所述雅可比矩阵A由以下部分组成：矩阵P，其由所述子域的块对角线和与所述块对角线相邻的非零数据块组成；矩阵E，其由未包含在所述矩阵P中的非零数据块组成，所述矩阵E由所述子域内的内部网格块的矩阵E’和具有所述储层中的相邻子域的网格单元的边界网格单元块的矩阵E”组成；

针对每个子域的所述网格块在所分配的处理器节点中采用为所述子域选择的局部预条件子来执行所述流体特性的并行近似求解；

基于求得的并行近似解来更新所述储层的所述雅可比矩阵的全系统残差；

根据全局约化空间矩阵将所述全系统残差限制在全局约化空间中；

采用全局约化空间预条件子来执行对所述全系统残差的并行近似求解；

在所分配的处理器节点中，由每个子域的并行近似解和全系统残差的并行近似解来形成组合近似解的更新，以形成解向量；

采用所述解向量更新所述全系统残差；

对所述全系统应用迭代求解器以获得所述储层的所述网格单元的每个成分的流体物料平衡的测量值；

对通过应用所述迭代求解器获得的每个成分的流体物料平衡的测量值进行检验，以确定是否实现收敛；以及

如果未实现收敛，则返回至执行并行近似求解的步骤；或者

如果实现了收敛，则返回至所述储层仿真器。

2.如权利要求1所述的计算机实施方法，其中执行并行近似求解的步骤得到子域边界残差，并且还包括步骤：

将所述子域边界残差限制在全局约化空间中。

3.如权利要求1所述的计算机实施方法，其中对所述雅可比矩阵A中的数据执行并行近似求解的步骤包括移除所述矩阵E的数据的步骤。

4.如权利要求1所述的计算机实施方法，其中更新所述全系统残差的步骤包括针对所述子域执行所述雅可比矩阵A和局部解向量v的矩阵向量相乘的步骤。

5.如权利要求4所述的计算机实施方法，其中针对子域执行矩阵向量相乘的步骤产生所述储层中的相邻子域的边界网格单元块的残差。

6.如权利要求5所述的计算机实施方法，还包括将边界网格单元块的残差从所述子域的处理器传递至相邻子域的处理器的步骤。

7.如权利要求1所述的计算机实施方法，其中用于对全系统残差执行并行近似求解的步骤的所述全局约化空间预条件子为全局约化空间矩阵所述全局约化空间矩阵由包括对角块的矩阵和由矩阵的其他块组成的矩阵所组成。

8.一种利用储层仿真器对地下储层的流体特性进行仿真的数据处理系统，所述数据处理系统包括由至少一个主节点和多个处理器节点形成的集群计算机，所述流体特性取决于流体流动、物料平衡和压力-体积-温度关系，所述储层被划分为多个子域，每个子域代表了所述储层的包含多相多成分流体的有限体积并且被组织成一组具有流体相关数据的网格单元，对每个处理器节点分配所划分的储层的一个子域，所述集群计算机通过牛顿迭代对网格单元的关注流体特性参数进行仿真，所述数据处理系统还包括：

主处理器，其用于执行将每个储层子域的数据排列成雅可比矩阵A的步骤，所述雅可比矩阵A由以下部分组成：矩阵P，其由所述子域的块对角线和与所述块对角线相邻的非零数据块组成；矩阵E，其由未包含在所述矩阵P中的非零数据块组成，所述矩阵E由所述子域内的内部网格块的矩阵E’和具有所述储层中的相邻子域的网格单元的边界网格单元块的矩阵E”组成；

所述处理器节点包括以下装置：

针对每个子域的所述网格块在所分配的处理器节点中采用为所述子域选择的局部预条件子来执行对所述流体特性的并行近似求解的装置；

基于求得的并行近似解来更新所述储层的所述雅可比矩阵的全系统残差的装置；

根据全局约化空间矩阵将所述全系统残差限制在全局约化空间中的装置；

采用全局约化空间预条件子执行对所述全系统残差的并行近似求解的装置；

在所述分配的处理器节点中由每个子域的并行近似解和全系统残差的并行近似解来形成组合近似解的更新以形成解向量的装置；

采用所述解向量来更新所述全系统残差的装置；

对所述全系统应用迭代求解器以获得所述储层的所述网格单元的流体特性的测量值的装置；

对通过应用所述迭代求解器获得的流体特性的测量值进行检验以确定是否实现收敛的装置；以及

如果未实现收敛成，则返回至执行并行近似求解的步骤的装置；或

如果实现了收敛，则返回至所述储层仿真器的装置。

9.如权利要求8所述的数据处理系统，其中执行并行近似求解的装置得到子域边界残差，所述数据处理系统还包括：

将所述子域边界残差限制在全局约化空间中的装置。

10.如权利要求8所述的数据处理系统，其中对所述雅可比矩阵A中的数据执行并行近似求解的装置包括移除所述矩阵E的数据的装置。

11.如权利要求8所述的数据处理系统，其中更新所述全系统残差的装置包括对所述子域执行所述雅可比矩阵A的矩阵向量相乘的装置。

12.如权利要求11所述的数据处理系统，其中对所述子域执行矩阵向量相乘的装置产生所述储层中的相邻子域的边界网格单元块的残差。

13.如权利要求12所述的数据处理系统，所述处理器节点还包括将边界网格单元块的残差从所述子域的处理器节点传递至相邻子域的处理器节点的装置。

14.如权利要求8所述的数据处理系统，其中用于所述处理器节点执行全系统残差的并行近似求解的所述全局约化空间预条件子为全局约化空间矩阵所述全局约化空间矩阵由包括对角块的矩阵和由矩阵的其他块组成的矩阵所组成。