CN104115035B - 在巨型储层中的多孔介质仿真中大规模线性系统的多级求解 - Google Patents

在巨型储层中的多孔介质仿真中大规模线性系统的多级求解 Download PDF

Info

Publication number
CN104115035B
CN104115035B CN201380008972.9A CN201380008972A CN104115035B CN 104115035 B CN104115035 B CN 104115035B CN 201380008972 A CN201380008972 A CN 201380008972A CN 104115035 B CN104115035 B CN 104115035B
Authority
CN
China
Prior art keywords
reservoir
unit
grid
computer
data handling
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201380008972.9A
Other languages
English (en)
Other versions
CN104115035A (zh
Inventor
乔治·A·皮塔
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Saudi Arabian Oil Co
Original Assignee
Saudi Arabian Oil Co
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Saudi Arabian Oil Co filed Critical Saudi Arabian Oil Co
Publication of CN104115035A publication Critical patent/CN104115035A/zh
Application granted granted Critical
Publication of CN104115035B publication Critical patent/CN104115035B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V99/00—Subject matter not provided for in other groups of this subclass
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]

Abstract

在对从被组织为细网格模型的巨型储层的储层仿真器压力方程得出的数十亿个方程的线性系统进行求解的过程中，多级方法获取要使用的快速和粗网格解，以作为细网格问题的更加精确的初始估计。在所述粗网格上执行的求解迭代在降低了处理时间的同时产生了更好的初始估计，因此在所述细网格级需要更少的较昂贵迭代。

Description

在巨型储层中的多孔介质仿真中大规模线性系统的多级求解
发明人：JORGE A.PITA
技术领域
本发明涉及被称为巨型储层的计算机仿真，特别是在油储层和可压缩单相干气储层中的单相流体的仿真过程中的大规模线性系统的多级求解。
背景技术
申请人是发明人且是本申请共同所有的美国专利第7,526,418号是一种储层组分模拟器，该模拟器在共享存储器的超级计算机、分布式存储器超级计算机或被配置为计算机数据处理单元(CPU’s)的个人计算机(PC)的集群中执行仿真。其他使用CPU’s的储层仿真成果是美国专利第7,516,056号和7,684,967号。
对由细网格模型(fine-grid model)(地震级(seismic-scale)12.5米至25米)组成的巨型方程组的压力方程和温度方程的线性求解可占据地下地质构造的储层和盆地仿真的总仿真时间的50％或以上。线性求解器确定牛顿迭代所需的“校正”以收敛于基本非线性偏微分方程组的解，该偏微分方程组以一系列时间步(time step)定义了储层流体的流体流量、物质平衡和压力-体积-温度条件。对于公知巨型储层而言，单元(cell)数量可以达到百万甚至十亿或更多。因此，该线性求解过程计算强度大且耗费时间。
即使当前可千兆级计算的系统，即，由成千上万处理器组成的能够每秒计算10的15次方的系统，在对细网格模型的仿真模型使用了不准确的初始猜测的情况下，也无法有效地求解这些方程组。计算机系统很多昂贵的线性迭代浪费在基于不准确的初始猜测寻找合理求解的方向上。
发明内容
简略地，本发明提供了在多个数据处理器的数据处理系统中对巨型地下储层的计算机化仿真的新的和改进的方法。所述计算机化仿真是对巨型地下储层的储层参数的方程的迭代线性求解，该巨型地下储层被仿真为这样的模型，该模型被划分为按有序的单元系统排列的多个单元。所述仿真还基于储层单元的地质和流体特性信息。根据本发明的方法包括如下计算机处理步骤：将源自原始细网格级的储层的有序单元系统的信息映射到从该原始网格减少了数量的粗单元网格中，并在计算机系统中针对该粗单元网格的储层参数对公设(postulated)系统解矩阵的进行初始化。根据本发明的方法还包括如下计算机处理步骤：针对所接收的粗单元网格的储层参数的初始化的公设系统解矩阵在计算机系统中执行预条件处理共轭梯度外推、在计算机系统中将针对粗单元网格的预条件处理共轭梯度外推的结果转换为原始单元网格、以及在计算机系统中对原始单元网格的储层参数执行迭代线性求解。
本发明还提供了在计算机系统中用于巨型地下储层的计算机化仿真的新的和改进的数据处理系统。所述计算机化仿真是对巨型地下储层的储层参数的方程的迭代线性求解，该巨型仿真为这样的模型，该模型被划分为按有序的单元系统排列的多个单元。所述仿真还基于储层单元的地质和流体特性信息。根据本发明的数据处理系统包括多个数据处理器，每个数据处理器并行执行如下步骤：将源自储层的有序单元系统的信息从原始细网格级映射到从该原始网格减少了数量的粗单元网格；针对该粗单元网格的储层参数在计算机系统中对公设系统解矩阵进行初始化；在计算机系统中针对所接收的粗单元网格的储层参数的初始化的公设系统解矩阵执行预条件处理共轭梯度外推；将针对粗单元网格的预条件处理共轭梯度外推的结果转换到原始单元网格；以及在计算机系统中对原始单元网格的储层参数执行迭代线性求解。所述数据处理系统还包括用于存储原始单元网格的确定的储层参数的存储器。
本发明还提供了一种具有其中存储了计算机可操作指令的计算机可读介质的新的和改进的数据存储装置，所述指令使得包括多个数据处理器的数据处理系统在对巨型地下储层进行计算机化仿真时根据存储在所述数据存储装置中的计算机可操作指令集进行操作。该巨型地下储层被仿真为这样的模型，该模型被划分为按有序的单元系统排列的多个单元。该仿真通过对储层参数的方程进行迭代线性求解来做出并且还基于储层的单元的地质和流体特性信息。存储在所述数据存储装置中的计算机可操作指令使得所述数据处理系统执行下列步骤：将源自储层的有序单元系统的信息从原始细网格级映射到从该原始网格减少了数量的粗单元网格；针对该粗单元网格的储层参数在计算机系统中对公设的系统解矩阵进行初始化；针对所接收的粗单元网格的储层参数的初始化的公设系统解矩阵在计算机系统中执行预条件处理共轭梯度外推；在计算机系统中将针对粗单元网格的预条件处理共轭梯度外推的结果转换到原始单元网格；以及在计算机系统中对原始单元网格的储层参数执行迭代线性求解。
附图说明
图1是组织成许多单元的巨型地下油气储层的组分模型的等距视图。
图2A和图2B是源自图1的油气储层模型的不同尺寸的相邻单元组的放大示意图。
图3A和图3B是将从图1所示的储层的粗网格模型获得的初始解插值到细网格模型中的示意图。
图4A是为根据本发明的储层仿真所安排的异构计算机系统的示意图。
图4B是为根据本发明的储层仿真所安排的千兆级计算机系统的示意图。
图5是根据本发明的在地下储层的计算机化储层仿真期间在图4A和图4B的计算机系统中执行的一组数据处理步骤的功能框图。
图5A是图5的数据处理步骤的一部分的功能框图。
图6是图1的储层的部分相邻单元网格的二维排序的示意图。
图7是图1的储层的部分相邻单元网格的根据本发明的重构的二维排序的代数矩阵的结构的示意图。
图8是图1的储层的部分相邻单元网格的根据本发明的三维排序的示意二维图(代数矩阵的结构)。
图9是根据本发明的用于应用多级求解插值的示例储层模型网格的示意图。
图10是传统单网格求解处理和根据本发明的多级求解针对诸如25米地震级的第一尺寸的细网格的每时间步的处理分钟的对比绘图。
图11是传统单网格求解处理和根据本发明的多级求解针对另一个较小的尺寸(15米地震级)的细网格的每时间步的处理分钟的对比绘图。
图12是根据本发明的其数据受到处理的实际储层的模型的图像的等距视图。
具体实施方式
在附图中，字母M表示一部分地下油气储层的简化模型，基于从储层单元获得的地质和流体特性信息根据本发明在估计生产寿命上对基于操作条件和参数的所述部分地下油气储层的生产结果进行了仿真。因此，所获得的结果是可用的并用于仿真历史绩效和用于预测储层的产出。基于该仿真的结果，诸如在美国专利申请第7,526,418中描述和示出的模型之类的模型可以随后被形成并可用于评估和分析。美国专利申请第7,526,418是本申请的受让方所拥有的并且通过引用合并于此。
如模型M所示出的在期望的储层寿命上对其生产数据进行了仿真的该类型的示例储层通常是本领域技术人员所公知的巨型储层。巨型储层在地下可以具有几英里的长度、宽度和深度，例如，可以具有三千亿立方英尺的体积或规模。
模型M被划分成许多合适尺寸的单元C，其中的若干个相邻单元在图2A和图2B中以放大图1的形式和比例示出。在为了分析而将储层体积以公知的地震级划分成单元时，如图2A中的粗网格20那样的典型单元每个在160英尺左右。应当理解的是在储层中的共有基准面中沿着单元的横向或平面尺寸是200米和100米的单元也是常见的“粗”网格模型。在如图2B中22处所示的公知的细网格模型类型中，典型的单元可以是沿着其横向或平面尺寸15到20英尺或以下的单元。例如，当粗网格是100米而细网格目标模型是25米时，可以使用两级多网格仿真：第一级使用100米而第二级使用50米，方法稍后描述。
因此，图1的模型M表示由具有上述尺寸的百万或更多单元组成的储层。如上所述，在地震级数据的情况下，单元的数量可以比图1的情况多几百倍或以上。将认识到的是，在附图中所示的形成模型M的单元C比起模型M而言被显著地放大以用于说明。关于该尺寸的模型的信息和复杂性如前所述，通过引用美国专利申请第7,526,418合并。
本发明提供求解多孔介质的单相流模型的数十亿方程的快速多级方法。在本文描述的实施例中，单相流可以是略微可压缩单向流体油储层的，或者是可压缩单向气储层的。本发明的方法可以在如图4A中示出为数据处理系统D的具有很多CPU’s和加速装置的异构计算系统中实施。本申请的方法也可以在如图4B中所示的公知的超级计算机S中实施，例如IBM(Blue)超级计算机，目前已知版本是蓝色基因/P。本申请已经在具有最高65,536个计算机核的大型计算机系统上测试，并且相比现有方法，提供了3倍以上的求解时间的加速。由于单相流系统的线性求解占用了总仿真时间的极大部分，因此整体仿真时间的节约是显著的。
如将要阐述的，根据本发明的多级方法通过在粗网格(例如图2A中20所示)上求解系统来获得解空间的便宜但准确的初始估计和猜测。从粗网格获得的初始估计或猜测如箭头21示意性示出般地细化为如图2B中22的更精细的地质级网格，并且，成为更加准确的估计，因此细网格上所需的迭代更少。这种方法避免了由于细网格模型的每次迭代几倍于较粗模型的迭代的时间而导致的处理时间的增加。
本发明针对多级系统提供了强健的预条件子(preconditioner)和共轭梯度加速器。本发明还提供了能够快速将粗网格细化到细网格中的快速可并行插值方法，即使计算硬件由成千上万的计算核组成的情况下。如将要阐述的，在65,536核的IBM蓝色基因/P超级计算机上计算的结果表明，通过本发明实现了三倍或以上的更快速的求解时间。
迭代线性解
本发明提供了用于求解一般线性系统Ax＝R的千兆级并行处理，其中A是系统矩阵(通常被求解的是非线性方程组的雅可比矩阵)，x是系统解的矢量，R是线性系统右手侧的矢量。在储层仿真和盆地建模中，该系统在采用牛顿方法进行的守恒方程的非线性迭代期间反复出现。
在牛顿方法中，雅可比矩阵(J)、非线性残差(F)和变量的非线性迭代更新(s)通过等式Js＝-F关联。因此与上述正则线性系统Ax＝R相比，A表示雅可比“J”，x表示解更新“s”而R表示右手侧的非线性残差“-F”(采用负数标记使得牛顿更新直接添加到当前变量状态)。
本发明的方法是基于这样的发现，即，在研究更廉价的方法来产生准确的初始估计或猜测以求解从雅可比线性化产生的线性方程组中，细网格的数十亿方程组所需的迭代数量可极大减少，大大节约了计算机时间。本发明通过将细网格直接抽取或插值成为较粗的网格来建立问题的较粗的表示法。随后使用该较小的线性系统来向线性系统的解迭代。在得到粗网格初始解之后，本发明将该解插值到细网格中，由此细网格具有更佳的初始猜测，从而收敛所需的迭代更少(在细网格中的迭代比在粗网格上的迭代更加耗费)。
初始粗网格解成为细网格系统的相关或关联数量的网格的初始估计猜测，因此针对细网格的线性系统的解可以在仅几次迭代后收敛，这是由于该解开始点预期与真实细网格解接近。该策略在图3A和图3B中概念地示出，其中从粗网格模型20的粗网格单元26延伸的箭头组24示意性地表示从根据本发明的初始解到细网格模型22的相关或关联数量的细网格单元28、30、32和34的转换。
为了得到粗网格模型20中的初始解，执行将细网格属性从细网格模型22映射到粗网格20中的初始映射。在该实施例中，细网格属性是在泡点压力以上的基本遵守单相行为的储层的压力和组分。在本发明的不同实施例中，初始映射可以通过不同技术来执行。已发现双线性插值是处理时间适中的最有效方法。
一旦获得粗网格解，也可以通过双线性插值实现在图3B中示意性示出的到细网格22的再映射，或者通过成本更低但不太准确的直接单射(direct injection)法来实现该再映射。在直接单射中，如图3A和图3B中所示，诸如24之类的单元的单个的粗单元值被直接复制到细网格模型22的四个单元28、30、32和34中。双线性插值推荐用于完全异构模型，而直接单射在具有平滑连续属性变化的模型中是有利的(即，更快)。
现在考虑根据本发明的异构计算机或数据处理系统，如图4A中所示，数据处理系统D提供用于储层的模型M中网格块之间每相关时间步的每次牛顿迭代的流体运动的计算机化仿真。数据处理系统D包括一个或多个中央处理单元或CPU’s 40。CPU或CPU’s 40与储层单元地质和流体特性信息的储层存储器或数据库42和用户界面44相连。用户界面44包括图形化显示46(用于显示图形化图像)、打印机或其他合适的图像形成机构以及用户输入装置48，以用于让用户操作、访问和提供处理结果、数据库记录和其他信息的输出形式。
储层存储器或数据库42通常是在外部数据存储计算机52的存储器50中。插入数据库42包含：包括模型M中单元的结构、位置和组织的数据；以及关于井、处理设施的数据，包括测量的静态井下压力数据在内的时间相关的井生产数据的数据，包括测量的井口压力和注入速度数据在内的时间相关的注入井数据、地质信息和流体特性信息以及其他储层生产记录和参数，以用于储层仿真，下文将进行描述。
数据处理系统D的CPU计算机40包括处理器54和与处理器54耦接的内部存储器56，该内部存储器56用于存储操作指令、控制信息并按需用作存储缓存器或传输缓存器。数据处理系统D包括存储在CPU或CPU’s 40的存储器56中的程序代码58。根据本发明，程序代码58是计算机可操作指令的形式，如将要阐述的，该指令可以使CPU’s 40来回地传输数据用于让许多图形处理单元或GPU’s 60进行处理以仿真储层中的流体运动。
应当注意到的是，程序代码58可以是微码、程序、例程或符号化计算机可操作语言的形式，它们可以提供特定的有序操作集，该有序操作集可以控制数据处理系统D的功能并指引其操作。程序代码58的指令可以以非瞬时性有形计算机可读的形式存储在存储器56中或者以上述形式存储在计算机盘、磁带、传统硬盘驱动器、电只读存储器、光存储装置或者其上存储了计算机可使用介质的其他适当的数据存储装置上。程序代码58也可以作为计算机可读介质包含在数据存储装置上。
图形单元或GPU’s 60是通用可编程图形处理单元，同样称作GPU’s。如将阐述的，GPU’s被编程以使用单独单元的线性化的方程组来确定未知量。在一些情形下，还应当理解的是，如果期望，在数据处理系统D中还可以使用除GPU’s以外的处理器节点作为处理器节点。
虽然本发明与使用的特定计算机硬件无关，但是本发明的一个实施例是基于适当数量的四核CPU’s和多核CPU’s的。在示例系统实施例中使用的CPU’s 40是Intel四核Nehalem处理器或Intel六核Westmere处理器的形式，在示例实施例中，优选地GPU’s 60是440-核的NVidia Fermi M2070Q或512-核的NVidia Fermi QuadroPlex7000GPU’s。然而，应当理解的是，如下将要阐述的，还可以使用其他计算机硬件。
本发明经由GPU’s 60采用双级(dual-tier)方法来加速，相比早期方法，可以获得接近数量级的速度改进。本发明在包括CPU’s40和GPU’s 60两者的异构(混合式)计算机环境中实现储层仿真。因此，本发明提供了用于巨型地下储层的细网格储层仿真的线性系统的基于计算机的系统。
在本发明中得到的储层仿真处理的线性化的方程组处理序列的加速可以大量节约计算机时间，降低成本并允许在给定的时间限制内进行更多的储层研究。如将要阐述的，在一些情况下处理时间加速了三倍。在本发明中，计算的加速使得比从前更快的确定成为可能，因此储层仿真器可以与现场测量值的实时数据采集同步。
在图4B中，示意性地示出了诸如IBM蓝色基因/P之类的能够千兆级计算的大规模并行计算机或数据处理系统S。该计算机系统的进一步细节在例如美国专利第7,680,048中阐述。
计算机系统S包括具有以规则阵列或格子方式逻辑排列了用于节点间通信的大量计算节点72的计算核70，其集合地执行批量处理。计算机系统S和计算核70的操作通常由控制子系统74控制。包括在前端节点76中的各种附加处理器执行一定的辅助数据处理功能，而文件服务器78向数据存储装置提供诸如可旋转磁盘驱动器80和82之类的接口或者其他输入/输出(I/O)源。功能网络84提供计算核70和其他系统组件之间的主要数据通信路径。例如，存储在连接到文件服务器78的存储装置中的数据被加载并通过功能网络84存储到其他系统组件。
计算核70包括许多输入/输出(I/O)节点86和计算节点72。计算节点72执行根据本发明的计算加强的储层仿真处理，该处理需要大量的并行进行的处理。每个I/O节点86对N个计算节点72的各个集合处理I/O操作。计算核70包括M个I/O节点86的集合以及它们的相关的计算节点集合，因此计算核70包括总共M乘N个计算节点72。在计算核70中的计算节点的数量可以非常庞大。例如，在一个M＝1024(1K)且N＝64的实施方式中，总计有64K个计算节点。
通常，应用程序编程代码和计算核用于执行用户应用程序处理所需的其他数据输入以及计算核执行用户应用程序处理而产生的数据输出在功能网络84上的计算核外部通信。集合内的计算节点72与相应的I/O节点86在相应的本地I/O树状网络88上通信。I/O节点86依次连接到功能网络84，在该功能网络84上，它们与连接到文件服务器78的I/O装置通信，或者与其他系统组件通信。因此，功能网络84为计算节点处理所有I/O，需要非常大的带宽。功能网络84可以是到多个以太网交换机的千兆比特以太网接口的集合。本地I/O树状网络88逻辑上可以看成是功能网络84的扩展，这是因为尽管它们是与功能网络84物理分隔开的并遵守不同的协议，但是I/O操作通过两者进行。
控制子系统78管理计算核70中的计算节点72的操作。控制子系统74优选地是微机系统，该微机系统包括其自有的处理器或处理器群、内部存储器和本地存储装置，并且具有用于与系统管理员或类似人员交互的随附的控制台90。
控制子系统74还包括内部数据库，其维护计算核70中的计算节点的一定状态信息以及各种控制和/或维护应用程序，这些应用程序在控制子系统74的处理器上执行，并且控制计算核70中硬件分配、指引到计算节点的数据预加载以及执行一定的诊断和维护功能。控制子系统74在控制系统网络92上与计算核70的节点交流控制和状态信息。网络92与一组硬件控制器94耦接，该组硬件控制器94与相关的计算节点72的集合的节点以及它们各自的I/O节点86在相应的本地硬件控制网络上通信。硬件控制器94和本地硬件控制网络可以逻辑上认为是控制系统网络92的扩展，尽管它们在物理上分隔开。控制系统网络和本地硬件控制网络以比功能网络84低得多的数据速率操作。
除了控制子系统74之外，前端节点76也在计算核外部最优执行，该前端节点76包括处理器和存储器的集合，其为了效率或其他原因执行一定的附属功能。涉及大量I/O操作的功能通常在前端节点中执行。例如，交互数据输入、应用程序代码编辑或其他用户接口功能通常由前端节点76处理，如同应用程序编译那样。前端节点76与功能网络84耦接以和文件服务器78通信，并可以包括或耦接到交互工作站。如上所述，处理器P的进一步细节在美国专利第7,680,048中阐述，其通过引用合并于此。
由于千兆级计算以及未来的百万兆级计算倾向于采用具有很多CPU的异构架构，而且也可能采用很多加速器装置(GPU或多核集成芯片)，根据本发明优化了粗网格求解和细网格求解两者所采取的基本线性求解方法。如下将要阐述的，该优化在预条件处理级和共轭梯度加速级都提供了不同硬件组件之间最小的通信和数据传输。因此，根据本发明的多级方法得益于无论施加何种硬件配置都高效率的整体单处理框架。
流程F(图5)示出了采用本发明的方法的储层仿真的基本计算机处理序列和在应用本发明的典型实施例期间采取的计算序列。
读取地质模型(步骤100)：根据本发明的仿真开始于将地质模型读取为输入和时不变数据。在步骤100期间读入的地质模型采用二进制数据形式，其包含每个储层模型属性的每网格单元的一个数值。这些属性包括下列：岩石渗透率张量；岩石孔隙度，x、y和z方向上的单个单元尺寸；每个单元的顶深；每个现有流体界面(油气界面、气水界面和油水界面，如适用)的x-y-z位置。
在步骤100期间读入的时不变数据包括每个分量的流体特性组分和热动力属性(临界温度、临界压力、临界体积、离心因子、分子量、等张比容、偏移参数和二元相互作用系数)。时不变数据还包括流体相关渗透率表格，该表格提供了考虑中的储层岩石的给定流体饱和度的相对渗透率的值，以及储层温度，因为本模型是等温的。
将模型离散化(步骤102)：基于连接的渗透率和单元几何形状，对每个单元执行岩石传递率的计算并存储到存储器中。取决于输入数据(例如块-面渗透率或块-中心渗透率)，存在本领域人员熟悉的若干用于传递率计算的模型。此外，每个单元的孔隙体积被计算并存储在存储器中。
初始化储层(步骤104)：在进行仿真之前，必须计算储层中流体的初始分布。该处理涉及在每个单元迭代压力。在每个点的压力等于“基准(datum)”压力加上该点上流体的静压头。因为单元处的静压头取决于单元上的流体的柱状密度，根据状态方程(或EOS,下面描述)，密度本身取决于压力和流体组分，因此求解本质上是迭代的。在每个单元处，计算出的压力被用于计算新的密度，从该新的密度重新计算新的静压头和单元压力。当以这种方式迭代的压力不再进一步变化时，系统被平衡，储层被称作“已初始化”。
读取和组织数据(步骤106)：在步骤106期间读入的周期性数据是时变数据，同样，其必须在仿真期间的每个时间步处读取。该数据包括已经在仿真的“历史”时期(用来校准仿真器的已知现场生产数据的时期)期间被观测的每口井的油、气和水速率。该数据还包括将在“预测”阶段(期望仿真器预测的现场生产时期)期间规定的生产方针。生产方针数据包括诸如每口井或井组需要的速率和应当施加到仿真上的约束(诸如最大气-油比、每口井被允许的最小井底-井口压力等)之类的数据。该数据可基于“历史”阶段期间实际现场测量值或基于“预测”阶段期间期望的生产能力而随时间阶段发生变化。
通过本发明，节约时间的主要策略是生成细单元参数值的成本低的良好的初始估计或猜测。此外，当异构系统用作计算架构时，从基本迭代处理的降低的通信需求得到额外的加速获益。
考虑图6，其示出了网格中的储层的简单的自然排列的相邻单元C，为简化起见，只以二维方式绘制。对储层单元的这种单元分组示意性地在图5的步骤106处示出。单元C被配置到棋盘格图案的两个区别组中(画影线的表示组1而空白的表示组2)。组2中空白单元的未知量将采用偶数下标识别而组1的画影线的单元的未知量将采用奇数下标识别。以图6所示的方式示出的单元组或单元排序在步骤106期间产生由5种连接(connection)(东“e”、西“w”，北“n”、南“s”和中-对角“d”)组成的有限差分模板。用这种方式将储层单元的交替的单元分组或组织到两个区别组中，产生矩阵结构，该矩阵结构之后在步骤106(图5)期间被重新排序或配置为如图7示意性示出的矩阵结构。
应当注意的是，图7的矩阵结构，以重新排序的形式配置在四个不同的象限中，奇数未知量组在矩阵的上部而偶数未知量组在矩阵的下部。通过重新排列混合计算机硬件的计算顺序，新型变体被用来使得CPU-GPU环境的求解处理最优化，不仅意在最小化CPU和GPU之间的通信量，而且留意到对GPU’s的当前存储限制。该过程在混合GPU-CPU计算中将该流量和GPU存储需求减半，这是本发明的一个重要特征。
矩阵填充是稀疏的，所以只有标记的条目需要被存储到存储器中。为了便于标记，图7中奇数对角线项的左上象限被标记为C1，偶数的西-东-北-南项的右上象限被标记为M1，偶数的西-东-北-南项的左下象限被标记为V2而偶数对角线项的右下象限被标记为C2。
扩展为三维网格是直接的，其中两条新的条目出现在矩阵中：“上”和“下”连接，其与C1和C2矩阵中的对角线项相邻。换句话说，三维网格将这些对角线矩阵变换为如图8示出的三对角线形式。
但是，应当注意到的是，以地震级成网格的地质模型在很多情况下可包括故障和尖灭(pinch-out)，其不符合概括于此的规则连接集合，它们将根据地质解释而插入到模型中，即使它们在地震数据中不可见。这些附加的连接有必要受到相应处理，在本实施例中，有必要在线性系统的右手侧中进行(以便维持同样的疏松结构)。取决于连接项的量级，如果连接脆弱，本发明按需执行。但是，对于连接强健，需要更多迭代或更强预条件子的情况，收敛可能是个问题。
在后面的说明中将看到，通过以如步骤106(图5)的示意性地示出的方式排序和重新排序未知量而进行的配置成功提供了线性解，该线性解可以每次以大约一半未知量的形式分成两部分处理。这对GPU是有利的，有两个原因：从CPU到GPU的数据传输量被减半且残留在GPU存储器(其小于CPU的存储器)中的数据量也被减半。在同构计算机系统中，这可以转化为对这些操作期间的工作阵列中的CPU存储器的节约，因此可以处理更大的模型。
在将单元重新排序为两个子组(如所述的，标记为1和2)之后，形成矩阵A：
因此通过处理将解决的原始矩阵方程组是：
初始化
粗网格初始化(步骤108)：采用本方法对步骤106产生的有组织的数据进行处理在初始化步骤108开始进行，在该步骤期间图6的矩阵的组1中的所有单元的初始估计被设置为0：
x1＝0 (3)
以使得组2残差始终为0的方式在步骤108期间确定组2中所有单元的初始估计。这排除了在GPU中更新组2残差的需要并将GPU负担的幅度减少了大约一半。数学意义上，这通过下面的原始系统的操作来实现。线性方程残差给出如下：
将组2的残差设置为0：
ρ2＝R2-V2x1-C2x2＝0 (5)
给定x1＝0，针对x2在步骤108期间可产生初始估计为：
这意味着组2的单元的三对角线系统的解。应当注意到的是组1的单元均不需要被传输到CPU’s进行该处理。
在步骤108期间从以上方程(4)计算组1的单元的初始残差(注意x1＝0)。
ρ1＝R1-M1x2 (7)
共轭梯度外推
粗网格外推(步骤110)：在步骤110期间，对步骤108得到的初始数据执行预条件处理的共轭梯度外推。在步骤110期间，预条件子Q被选择并在共轭梯度外推期间使用以改善原始矩阵A的条件数。
一旦预条件子被选定，对于任意共轭梯度法而言需要矩阵-向量乘积AQ-1v＝Pv。根据本发明的处理技术中使用的选定的预条件子Q是本领域公知的“右”预条件子，其“从右边”对矩阵A进行预条件处理，因此不会改变原始系统残差或右手侧的值，如下所示：
A x＝R (8)
AQ-1Qx＝R (9)
因此，让P＝AQ-1而y＝Qx，矩阵方程组变成：
P y＝R (10)
这意味着在通过共轭梯度加速针对未知量“y”对系统求解之后，真解“x”通过逆变换获得：
x＝Q-1y (11)
针对根据本发明的处理，Q被选为高斯赛德尔型的Z-线预条件子。术语“Z-线”出现是因为主对角线子矩阵C1和C2(图6和图7)包括3D储层网格的垂直柱(即，储层网格的Z-方向)：
应当注意到,如果方程式(12)与上面的方程式(1)比较，此处A＝Q+M1。因此，下列方程成立：
P v＝AQ-1v＝(Q+M1)Q-1v＝v+M1Q-1v
或者，以全矩阵表示：
根据本发明，由于组2的单元的残差一直是0：
矩阵-向量相乘，简化为：
其中：
通过展开方程(16)中的项，很容易针对δ1和δ2通过两步对该系统求解：
执行方程(18)使得矩阵-向量乘积V2δ1首先被计算，之后对三对角线系统求解δ2。
得到预条件处理的矩阵-向量乘积Pv，可以表示为：
一旦确定了预条件处理的矩阵-向量乘积Pv，使用下面描述的三种加速算法(ORTHOMIN、BiCGSTAB和TFQMR)的任意一种可以在步骤110(图5)期间执行共轭梯度加速。如上所述，加速器的选择影响GPU的性能且双共轭兰素斯法比起传统的Krylov法具有明显的存储器优点，因为Krylov法必须存储多个正交方向向量以便提供强健的性能。考虑到异构数据处理系统中GPU存储器的限制，该额外的存储是不切实际的，但是很容易通过基于CPU的同构系统S接纳，该系统中存储器充裕且不存在主机到装置的通信成本。
正交极小化共轭梯度法(ORTHOMIN(K))
该方法的计算步骤如下：
1)计算初始残差：
r0＝b-Ax0
2)设置初始方向向量：
p0＝r0
3)迭代k＝0,1,2…,直到收敛，进行：
4)
5)xk+1＝xk+αkpk
6)rk+1＝rk-αkApk
7)i＝j-K+1,…,k
8)
9)结束
需注意K是正交方向的数量，其存储在存储器中用于在迭代期间重复使用。该存储器通常在CPU上实施而很少在GPU架构中负担得起，这是因为GPU当前的总存储受限。由于这个原因，对于在同构架构(仅CPU)上的本发明的实施例而言，ORTHOMIN是优选的，但是对本发明的异构(CPU+GPU)实施例而言并非优选。
通过残差(r)的L2-范数的值来测量收敛标准。通常，当该范数小于规定公差(0.01或更小)时，处理已收敛。在本发明的模型中，相对规定公差来核对相对残差值(当前残差除以初始残差)的范数以避免初始残差的过大或过小值使收敛所需的迭代数量产生偏差，即：
双共轭梯度稳定法(BiCGSTAB)
该方法的计算步骤如下：
1)计算初始残差：
r0＝b-Ax0
2)设置初始方向矢量：
p0＝r0
3)迭代k＝0,1,2,…,直到收敛，执行：
4)
5)sk＝rk-αkApk
6)
7)xk+1＝xk+αkpk+wksk
8)rk+1＝sk-wkAsk
9)
10)pk+1＝rk+1+βk(pk-wkApk)
11)结束
通过残差(r)的L2-范数的值来测量收敛标准。通常，当该范数小于规定公差(0.01或更小)时，处理已收敛。在本发明的模型中，相对规定公差来核对相对残差值(当前残差除以初始残差)的范数以避免初始残差的过大或过小值使收敛所需的迭代数量产生偏差，即：
无转置准最小残差法(TFQMR)
该共轭梯度加速技术的计算步骤如下：
1)计算初始残差和初始变量：
r0＝b-Ax0
w0＝u0＝r0
v0＝Au0
d0＝θ0＝σ0＝0
2)计算初始误差范数：
τ0＝||r0||
3)迭代k＝0,1,2,…,直到收敛，执行：
4)如果k是偶数，计算：
uk+1＝uk-αkvk
5)wk+1＝wk-αkAuk
6)
7)
8)
9)τk+1＝τkθk+1ck+1
10)
11)xk+1＝xk+σk+1dk+1
12)如果k是奇数，计算：
uk+1＝wk+1+βk-1uk
vk+1＝Auk+1+βk-1(Auk+βk-1vk-1)
13)结束
同样地，通过残差(r)的L2-范数的值来测量收敛标准。通常，当该范数小于规定公差(0.01或更小)时，处理已收敛。在本发明的模型中，相对规定公差来核对相对残差值(当前残差除以初始残差)的范数以避免初始残差的过大或过小值使收敛所需的迭代数量产生偏差，即：
共轭梯度加速器的选择在如图4A所示的类型的本发明的采用CPU’s和GPU’s的异构实施例中非常重要，该类型在于2011年2月8日提交的共同拥有的美国专利申请第13/023,076中描述，其申请人是发明人。在图4B示出的同构计算机(例如IBM蓝色基因/P)类型中，选择的共轭梯度加速器的类型不那么重要。
通常，Orthomin加速器为同构(基于CPU)环境提供非常合适且强健的加速，但是在需要在CPU和GPU之间传输多个正交方向的异构计算机环境中受到限制。由于这样的原因，本发明通常将其他的共轭梯度变体用于异构环境，例如上述的BiCGSTAB(双共轭梯度稳定法)和TFQMR(无转置准最小残差法)。
转换解到原始系统中
粗网格转换(步骤112)：如前所述，在步骤110期间的共轭梯度外推处理产生了系统Py＝R(方程10)的解“y”，其中P＝AQ-1而y＝Qx。之后，根据图5，在步骤112期间由CPU执行处理以将解转换到原始系统中。在该处理中，如系统Qx＝y的三对角线解或如下矩阵所表示的，获得x＝Q-1y：
矩阵方程(20)的第一行产生：
目前为止，所有计算使用组1的未知量，因此，减轻了GPU传输、存储和计算的整个方程组的一半，因此节约了50％的时间和所需的计算机存储器。但是，应当注意到的是，这意味着y2不可用，因为针对组2的单元不存在残差更新，这是由于在处理被初始化时将它们设置为0。为了求解x2，原始系统第二行(2)可以展开为：
V2x1+C2x2＝R2 (22)
如下求解x2：
这只涉及对组2单元的三对角线解，因此，同样地，只有整个系统数据的一半需要被数据处理系统D的GPU’s 60传输和处理。对于基于CPU的同构系统而言，因为不需要主机-装置之间的通信，因此其优点在于工作阵列存储装置中的存储器节约。
多级插值
图9示意性地示出了根据本发明的将粗网格解转换到细单元网格的多级插值法的典型情形。在图9中，点P(x，y)处的值将被粗网格上更大网格间隔处(例如点Q、R、S和T)定义的值插值。本发明可以通过三种插值方法来执行，阐述如下。
方法1-双线性插值
点R和T之间的线性插值给出：
类似地，点Q和S之间的线性插值给出：
对期望位置P(x，y)插值
其中f(x,y1)和f(x,y2)由上面的前两个表达式给出。应当注意到的是，该插值方法在沿着与x或y方向平行的线上是线性的，但是沿着其他直线是二次的。
方法2-直接单射(Direct Injection)
在一些情况下，使用最近点替代上述方法1中的执行插值操作是足够准确的。在这种情况下，直接选择f(x,y)＝f(x1,y1)。尽管这样节约了一些计算机时间，但是通常结果不准确。但是，直接单射对于轻微异构的地质模型而言是有利的。
方法3-四点平均
该方法是成本妥协的方法，准确性介于前两者之间。采取最近的相邻点的直接算数平均。在这种情况下，选择f(x,y)：
对这些方法和其他插值法的进一步讨论可以在其他文献中得到(例如，1999年11月出版的SCAI GMD Report第69期由K.Stueben撰写的文章“A Review of AlgebraicMultigrid”)。尽管如上所述方法1沿着非平行线条是二次的(否则是线性的)，其并行化在多处理器环境中是极其高效的。三次法(例如，样条插值(splines))虽然更加准确，但是通常无法弥补带来的额外计算成本，因此在此不做考虑。
细网格仿真(步骤114)：在步骤114期间，基于在步骤110期间为粗网格获得并在步骤112期间转换为原始细网格模型的单相流的初始估计结果，使用更加准确的初始估计来执行在当前关注的时间步期间对储层的单相流的仿真。储层仿真过程于是按这样的方式进行。在步骤114期间的细网格求解使用与已描述的粗网格求解相同的预条件子和加速器，这在申请人为发明人的标题为“Seismic-scale Reservoir Simulation of GiantSubsurface Reservoirs using GPU-Accelerated Linear Equation SystemsCalculations”的专利申请第13/023,076号中呈现。如图5A所示，在步骤114期间的细网格仿真对于二级求解而言包括细网格初始化步骤114a、细网格外推步骤114b和细网格转换步骤114c。应当认识到的是三级求解将具有其中每个均由初始化、外推和转换运动三个步骤组成的两个序列的粗网格仿真以及上述三个中的一个细网格仿真。
解更新(步骤116)：在步骤116期间，在步骤114期间对线性方程组求解而获得的解矢量δx表示非线性迭代循环中上述牛顿方法的更新的解矢量。虽然在很大程度上这是本领域公知的“牛顿迭代”，但是在本发明中进行一些核对以抑制解矢量以便改善仿真的数值稳定性。因此，不总是采用完整的“牛顿步骤”。用户控制的参数可以提供给关注的流体运动变量。显然，这是取决于储层模型的，并且可以由熟悉被仿真的储层或现场的仿真工程师超越这些限制。
收敛测试(步骤118)：按照用户规定的公差来核对从步骤116产生的线性方程的各个残差。如果满足这些公差，则退出非线性迭代循环，解输出在步骤120期间针对目前时间步被写入到文件并且时间步在步骤122期间前进到下一级。
如果不满足这些公差，根据非线性迭代循环的处理返回到步骤108并继续。但是如果非线性迭代的数量变得过多，则做出削减时间步大小的决定并从步骤108开始针对相同时间级再次重复整个非线性迭代循环。迭代的过多数量是解已偏离并且变化可能太多以至于无法采用之前选择的时间步充分建模的迹象。期望的是时间步的削减不仅减少这些变化的幅度，而且还增加雅可比矩阵的对角线优势，其总是对线性求解器的收敛具有有利效果。
写输出(步骤120)：在三维网格的形式中确定的流体运动变量的测量值在每个时间步的结尾处以二进制格式作为磁盘I/O写出。同时关于其他所需数据或测量值的井信息(诸如每口井的气、油和水生产速度；每口井的气/油比(GOR)以及每口井的静态井口压力(SWP)之类)被写出。
使时间步前进(步骤122)：在步骤120期间针对当前时间步将解输出数据写入到文件后，时间步在步骤122期间前进到下一级并且处理返回到步骤106以继续仿真器的处理。
并行化
代码并行化遵循针对多核架构的当前标准的高性能计算(HPC)方法：将OpenMP共享存储器模型用于每个计算机节点内的并行化/多线程，而将MPI信息传递编程用于各个节点之间的并行化。
在同构架构(仅CPU)的环境中，这是当今大多数PC集群和超级计算机的典型设置。例如，IBM蓝色基因/P超级计算机由每机架(rack)1024个节点组成，每个节点由4个计算核组成。OpenMP被应用到这4个核中多线程执行代码，而MPI被用于在各个节点之间并行执行代码。因此，由1024个节点(4096个核)组成的机架具有1024个MPI进程(或任务)且每个进程4个OpenMP线程。
但是，应当注意到的是，通过本发明实施的计算机代码可以使用MPI/OpenMP资源的任意分配。例如，前述情况也可以运行为4096个MPI进程(或任务)且每进程1个OpenMP线程(有些人称之为“虚拟节点”并行化或VN并行化，这是因为节点中的每个核起到单独的MPI进程的作用)。也可以运行为2048个MPI进程(或任务)，每进程2个OpenMP线程(有些人称之为“双节点”并行化或DUAL并行化，其中每个节点执行2个任务，每个任务两个线程)。
已发现本发明的最有效的模式是上述第一模式(每节点1进程和4个OpenMP线程，有些人称之为“SMP节点”)。这对于非常大的需要数千个核来运行的情况是典型的：SMP模式比VN模式和DUAL模式有利的是其降低了同一时间运行数千个MPI进程的开销。通常认为在千兆级(以及一直到百万兆级)系统中，将多核作SMP使用(即，使用OpenMP)是有效利用这些系统的关键点。
图10示出了在使用24,576个核的IBM蓝色基因/P超级计算机上针对在储层的泡点压力以上包含单相流体的现有储层的60亿单元模型储层仿真器的结果(每个单元3个未知量，即，黑油型)。网格尺寸是地震级(25米)。直接在细网格上求解该系统的传统单网格求解每时间步平均运行3.4分钟。
求解15亿单元的线性系统(dx＝dy＝50米)的2级方法首先产生初始猜测，将该初始猜测插值到60亿单元的网格(dx＝dy＝25米)上并使用该进行了网格细分的初始猜测来求解细网格线性系统。该2级方法每个时间步花费0.96分钟，因此总体加速3.54倍。
求解15亿单元的线性系统(dx＝dy＝50米)的3级方法首先将该解作为初始猜测提供给30亿单元的线性系统，随后，该30亿单元的线性系统被提供给60亿单元的线性系统(dx＝dy＝25米)。该3级方法每时间步花费0.74分钟，因此总体加速4.59倍。从图10中可以看到对3级方法的总时间的单项贡献，表明更加昂贵的细网格迭代所花费的时间实际上相比其他两个粗级(15亿和30亿)被最小化了。
图11示出了在IBM蓝色基因/P超级计算机(使用49,152个核)上运行的实际储层(每口井3个未知量，即，黑油型)的100亿单元的储层模型的结果。网格尺寸是dx＝dy＝15米。直接在细网格上求解该系统的传统的单网格求解，每时间步平均运行2.6分钟。
求解25亿单元的线性系统(dy＝dy＝30米)的2级方法首先产生初始猜测，将该初始猜测插值到100亿单元的网格(dx＝dy＝15米)上并且使用该进行了网格细分的初始猜测来求解细网格线性系统。该2级方法每时间步花费0.9分钟，因此总体加速2.9倍。图12是该模型的黑白图像。在实际应用中，由色调和饱和度变化来表示在盆地模型中的温度分布。
图12示出了在网格尺寸dx＝43米、dy＝25米而dz＝13米的采用2级求解的情况下。在65,636核IBM蓝色基因/P超级计算机上运行的采用原型盆地模型的动力学的110亿个单元的温度方程的当今温度分布。图12是该模型的黑白图像。在实际应用中，由色调和饱和度变化来表示盆地模型中的温度分布。
虽然本发明的不同实施例可针对讨论中的特定情况进行调整，但是通用方法是首先基于求解粗网格问题以成本高效地使用2级插值获得近似解。该近似解之后用作实际的细网格问题的初始猜测或估计。随后应当验证细网格级上需要若干次(例如，小于5)迭代。还期望的是用两级或三级网格粗化来测试。如果使用三级比两级所减少的时间是很多的，则可以测试额外级。
应当注意到的是在该实施例中所主张的模型方程式使用压力和摩尔量而不是流体饱和度作为主要未知量。摩尔方程式的示例在标题为“Highly-Parallel ImplicitCompositional Reservoir Simulator for Multi-Million Cell Models”的共同拥有的美国专利No.7,526,418中阐述，其通过引用合并于此。流体饱和度对于从细到粗网格的多级插值而言是特别不利的，这是因为它们的值从一组单元到下一组单元会剧烈地变化。虽然上述引用的专利主题是针对储层仿真，但是相同的优点和考虑应用到盆地建模的温度方程。
从上述内容可以看出本发明提供了通过降低花费在线性方程组的昂贵的求解的时间(由单向流体导致，其能够占据总仿真时间的极大部分)来加速储层仿真器计算的特定部分的方法。该新方法的优势的增大直接与除了模型的断层、断裂和地层尖灭之外的地质非均质性和网格尺寸细化相关。具体地，地质的非均质性越高以及模型网格越精细，则(在计算时间和方法强健度上)多级方法相比单级求解的优势就越大。
本发明提供的多级方法形成廉价但精确的初始估计，该初始估计可以无需浪费对解空间搜索而预先解决细网格问题。本发明的多级方法提供了在可能由几万个计算机核组成的并行计算机上的多级之间的求解的快速插值/细化。
已经充分描述本发明，使得本领域的普通技术人员能够复现并获得在本文中所提及的结果。但是，本发明所述的技术领域中的任意技术人员可以进行于此未做描述的修改，将这些修改应用至确定的处理方法或使用其结果，要求在权利要求中的请求保护的主题；这些修改应当涵盖在本发明的范围内。
应当注意和理解的是可以在不脱离权利要求中阐述的本发明的精神和范围的情况下对上述详细描述的本发明做出改进和修改。

Claims (22)

1.一种在多个数据处理器的数据处理系统中对巨型地下储层进行计算机化仿真的方法，该计算机化仿真是针对巨型地下储层的储层参数的方程进行迭代线性求解，所述巨型地下储层被仿真为划分成以有序单元系统排列的多个单元的模型，所述仿真进一步基于储层的单元的地质和流体特性信息，所述方法包括计算机处理步骤：
(a)将来自所述储层的有序单元系统的信息从原始网格的原始细网格级映射到从该原始网格降低了数目的粗单元网格中；
(b)在计算机系统中针对所述粗单元网格的储层参数来对公设系统解矩阵进行初始化；
(c)通过确定预条件处理的矩阵-矢量乘积在所述计算机系统中对接收的所述粗单元网格的储层参数的初始化的公设系统解矩阵执行预条件处理的共轭梯度外推；
(d)在所述计算机系统中将所述粗单元网格的预条件处理的共轭梯度外推的结果转换到所述原始网格；以及
(e)在所述计算机系统中针对所述原始网格的单元的储层参数执行迭代线性求解。
2.根据权利要求1所述的方法，其中所述数据处理系统包括显示器，并且所述方法进一步包括形成针对所述原始网格的单元的储层参数执行迭代线性求解的结果的输出显示的步骤。
3.根据权利要求1所述的方法，其中所述储层参数包括在所述储层的单元中的流体流动。
4.根据权利要求1所述的方法，其中所述储层参数包括在所述储层的单元中的单相流体流动。
5.根据权利要求1所述的方法，其中在所述计算机系统中的针对所述粗单元网格的储层参数的预条件处理的共轭梯度外推的处理的所述计算机处理步骤包括步骤：形成储层中的第一组相互间隔开的单元的解的矩阵值的初始估计。
6.根据权利要求5所述的方法，进一步包括计算机处理步骤：
确定与所述第一组相互间隔开的单元相邻的第二组相互间隔开的单元的解的矩阵值的初始估计。
7.根据权利要求5所述的方法，进一步包括计算机处理步骤：
确定所述第一组相互间隔开的单元的残差矩阵值。
8.根据权利要求1所述的方法，其中所述执行预条件处理的共轭梯度外推的计算机处理步骤包括步骤：针对粗单元网格的储层参数在公设系统解矩阵上形成矩阵预条件子。
9.根据权利要求6所述的方法，其中所述在计算机系统中将预条件处理的共轭梯度外推的结果进行转换的计算机处理步骤包括步骤：
对矩阵值执行三对角线求解以求解储层中的所述第一组相互间隔开的单元；以及
对矩阵值执行三对角线求解以求解储层中的所述第二组相互间隔开的单元。
10.根据权利要求1所述的方法，其中所述对原始网格的单元执行迭代线性求解的步骤包括计算机处理步骤：
在所述计算机系统中针对所述原始网格的单元对公设系统解矩阵进行初始化；以及
在所述计算机系统中对接收到的所述原始网格的单元的储层参数的初始化的公设系统解矩阵执行预条件处理的共轭梯度外推。
11.一种用于在计算机系统中对巨型地下储层进行计算机化仿真的数据处理系统，该计算机化仿真是针对巨型地下储层的储层参数的方程的迭代线性求解，所述巨型地下储层被仿真为划分成以有序单元系统排列的多个单元的模型，所述仿真进一步基于储层的单元的地质和流体特性信息，所述数据处理系统包括：
(a)多个数据处理器，每个以并行方式执行如下步骤：
(1)将来自所述储层的有序单元系统的信息从原始网格的原始细网格级映射到从该原始网格降低了数目的粗单元网格中；
(2)在所述计算机系统中针对所述粗单元网格的储层参数来对公设系统解矩阵进行初始化；
(3)通过确定预条件处理的矩阵-矢量乘积在所述计算机系统中对接收的所述粗单元网格的储层参数的初始化的公设系统解矩阵执行预条件处理的共轭梯度外推；
(4)将所述粗单元网格的预条件处理的共轭梯度外推的结果转换到所述原始网格；以及
(5)在所述计算机系统中针对所述原始网格的单元的储层参数执行迭代线性求解；和
(b)存储器，其用于存储所述原始网格的单元的确定的储层参数。
12.根据权利要求11所述的数据处理系统，其中所述多个处理器包括异构的计算机节点组。
13.根据权利要求11所述的数据处理系统，其中所述多个处理器包括同构的计算机节点组。
14.根据权利要求11所述的数据处理系统，其中所述数据处理系统进一步包括形成所述原始网格的单元的储层参数的输出图像的显示器。
15.根据权利要求11所述的数据处理系统，其中所述储层参数包括在该所述储层的单元中的流体流动。
16.根据权利要求11所述的数据处理系统，其中所述储层参数包括在所述储层的单元中的单相流体流动。
17.根据权利要求11所述的数据处理系统，其中所述处理器在执行针对粗单元网格的储层参数的预条件处理的共轭梯度外推的步骤时执行形成储层中的第一组相互间隔开的网格的解的矩阵值的初始估计的步骤。
18.根据权利要求17所述的数据处理系统，其中所述处理器进一步执行计算机处理步骤：
确定与所述第一组相互间隔开的单元相邻的第二组相互间隔开的单元的解的矩阵值的初始估计。
19.根据权利要求17所述的数据处理系统，其中所述处理器进一步执行计算机处理步骤：
确定所述第一组相互间隔开的单元的残差矩阵值。
20.根据权利要求11所述的数据处理系统，其中所述处理器在执行预条件处理的共轭梯度外推的步骤时执行针对粗单元网格的储层参数在公设系统解矩阵上形成矩阵预条件子的步骤。
21.根据权利要求18所述的数据处理系统，其中所述处理器在执行将预条件处理的共轭梯度外推的结果进行转换的步骤时执行如下步骤：
对矩阵值执行三对角线求解以求解储层中的所述第一组相互间隔开的单元；以及
对矩阵值执行三对角线求解以求解储层中的所述第二组相互间隔开的单元。
22.根据权利要求11所述的数据处理系统，其中所述处理器在执行原始网格的单元的储层参数的迭代线性求解时执行计算机处理步骤：
在所述计算机系统中针对所述原始网格的单元对公设系统解矩阵进行初始化；以及
在所述计算机系统中对接收到的所述原始网格的单元的储层参数的初始化的公设系统解矩阵执行预条件处理的共轭梯度外推。
CN201380008972.9A 2012-02-09 2013-02-08 在巨型储层中的多孔介质仿真中大规模线性系统的多级求解 Expired - Fee Related CN104115035B (zh)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
US201261596948P 2012-02-09 2012-02-09
US61/596,948 2012-02-09
PCT/US2013/025280 WO2013119906A1 (en) 2012-02-09 2013-02-08 Multi-level solution of large-scale linear systems in simulation of porous media in giant reservoirs

Publications (2)

Publication Number Publication Date
CN104115035A CN104115035A (zh) 2014-10-22
CN104115035B true CN104115035B (zh) 2017-12-01

Family

ID=47843376

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201380008972.9A Expired - Fee Related CN104115035B (zh) 2012-02-09 2013-02-08 在巨型储层中的多孔介质仿真中大规模线性系统的多级求解

Country Status (5)

Country Link
US (1) US9378311B2 (zh)
EP (1) EP2812737B1 (zh)
CN (1) CN104115035B (zh)
CA (1) CA2862900C (zh)
WO (1) WO2013119906A1 (zh)

Families Citing this family (20)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011059537A1 (en) * 2009-11-12 2011-05-19 Exxonmobil Upstream Research Company Method and system for rapid model evaluation using multilevel surrogates
US9626466B2 (en) 2010-11-23 2017-04-18 Exxonmobil Upstream Research Company Variable discretization method for flow simulation on complex geological models
WO2013148021A1 (en) 2012-03-28 2013-10-03 Exxonmobil Upstream Research Company Method for mutiphase flow upscaling
GB2512372B (en) * 2013-03-28 2020-07-29 Total Sa Method of modelling a subsurface volume
WO2015031749A1 (en) * 2013-08-30 2015-03-05 Schlumberger Canada Limited Stratigraphic function
US10208577B2 (en) * 2013-10-09 2019-02-19 Chevron U.S.A. Inc. Method for efficient dynamic gridding
US9910173B2 (en) * 2013-11-15 2018-03-06 Schlumberger Technology Corporation Saturation end-point adjustment
US10634814B2 (en) * 2014-01-17 2020-04-28 Conocophillips Company Advanced parallel “many-core” framework for reservoir simulation
JP5894645B1 (ja) * 2014-08-29 2016-03-30 株式会社日立製作所 半導体装置及びその品質管理方法
JP5901712B2 (ja) * 2014-08-29 2016-04-13 株式会社日立製作所 半導体装置および情報処理装置
JP6021864B2 (ja) * 2014-08-29 2016-11-09 株式会社日立製作所 半導体装置および情報処理装置
JP5903471B2 (ja) * 2014-08-29 2016-04-13 株式会社日立製作所 半導体装置および情報処理装置
US10534877B2 (en) * 2015-05-19 2020-01-14 Schlumberger Technology Corporation Adaptive multiscale multi-fidelity reservoir simulation
US10520643B2 (en) * 2015-10-20 2019-12-31 Pgs Geophysical As Geophysical inversion using sparse modeling
US11126762B2 (en) * 2018-02-28 2021-09-21 Saudi Arabian Oil Company Locating new hydrocarbon fields and predicting reservoir performance from hydrocarbon migration
US10822925B2 (en) 2018-04-26 2020-11-03 Saudi Arabian Oil Company Determining pressure distribution in heterogeneous rock formations for reservoir simulation
US11461514B2 (en) * 2018-09-24 2022-10-04 Saudi Arabian Oil Company Reservoir simulation with pressure solver for non-diagonally dominant indefinite coefficient matrices
CN111565159B (zh) * 2020-04-13 2022-08-23 重庆邮电大学 基于无转置极小残差的迭代大规模mimo信号检测方法
US20230332490A1 (en) * 2022-04-14 2023-10-19 Saudi Arabian Oil Company Method and system for performing reservoir simulations using look-ahead models
CN117436370B (zh) * 2023-12-06 2024-03-19 山东省计算中心（国家超级计算济南中心） 面向流体力学网格生成的超定矩阵方程并行方法及系统

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20080167849A1 (en) * 2004-06-07 2008-07-10 Brigham Young University Reservoir Simulation
US7526418B2 (en) 2004-08-12 2009-04-28 Saudi Arabian Oil Company Highly-parallel, implicit compositional reservoir simulator for multi-million-cell models
US7349833B2 (en) 2004-10-01 2008-03-25 Seiko Epson Corporation 2D central difference level set projection method for ink-jet simulations
US7596480B2 (en) 2005-04-14 2009-09-29 Saudi Arabian Oil Company Solution method and apparatus for large-scale simulation of layered formations
US7516056B2 (en) 2005-04-26 2009-04-07 Schlumberger Technology Corporation Apparatus, method and system for improved reservoir simulation using a multiplicative overlapping Schwarz preconditioning for adaptive implicit linear systems
EP1893845A1 (en) 2005-06-14 2008-03-05 Services Pétroliers Schlumberger Apparatus, method and system for improved reservoir simulation using an algebraic cascading class linear solver
US7680048B2 (en) 2006-10-06 2010-03-16 International Business Machiens Corporation Method and apparatus for routing data in an inter-nodal communications lattice of a massively parallel computer system by dynamically adjusting local routing strategies
CN101661514B (zh) * 2008-05-21 2014-04-16 中国石化股份胜利油田分公司地质科学研究院 一种油藏黑油模型数值模拟系统
CA2759199A1 (en) * 2008-10-09 2010-04-15 Chevron U.S.A. Inc. Iterative multi-scale method for flow in porous media
WO2011097055A2 (en) 2010-02-02 2011-08-11 Conocophillips Company Multilevel percolation aggregation solver for petroleum reservoir simulations

Also Published As

Publication number Publication date
WO2013119906A1 (en) 2013-08-15
CA2862900C (en) 2017-07-04
US9378311B2 (en) 2016-06-28
US20130226540A1 (en) 2013-08-29
EP2812737B1 (en) 2018-04-11
CN104115035A (zh) 2014-10-22
EP2812737A1 (en) 2014-12-17
CA2862900A1 (en) 2013-08-15

Similar Documents

Publication Publication Date Title
CN104115035B (zh) 在巨型储层中的多孔介质仿真中大规模线性系统的多级求解
Zhou Parallel general-purpose reservoir simulation with coupled reservoir models and multisegment wells
US10762258B2 (en) Parallel solution for fully-coupled fully-implicit wellbore modeling in reservoir simulation
Ciaurri et al. Application of derivative-free methodologies to generally constrained oil production optimization problems
CN102648425B (zh) 用于储集层仿真的多尺度有限体积方法
CN104136942B (zh) 用于大规模并行储层仿真的千兆单元的线性求解方法和装置
CN101310272B (zh) 地下水流仿真中使用的多尺度有限体积法
DK2673465T3 (en) RESERVOIR SIMULATION In seismic SCALE OF GIGANTIC underground reservoirs USING GPU ACCELERATED LINEAR EQUATION SYSTEMS
CN109964151A (zh) 利用加速的含水层计算的并行储层模拟
Gyulassy et al. The parallel computation of Morse-Smale complexes
Tang et al. Use of low-fidelity models with machine-learning error correction for well placement optimization
Schornbaum Block-structured adaptive mesh refinement for simulations on extreme-scale supercomputers
Farahi et al. Model-based production optimization under geological and economic uncertainties using multi-objective particle swarm method
Fuchs et al. Finite volume methods for wave propagation in stratified magneto-atmospheres
Young et al. High-Performance Black Oil Computations
Yang et al. Accelerating linear solvers for reservoir simulation on GPU workstations
Deng A Parallel-In-Time Gradient-Type Method For Optimal Control Problems
Kompantsev Viability of Physics-Aware Deep-Learning based Proxy Reservoir Simulator for use in Well Control Optimization
SOUZA A MsCV framework using a non-orthodox MPFA-D for the simulation of two-phase flows on truly unstructured grids
la Cour Christensen Multilevel techniques for Reservoir Simulation
Niu Implementation of Linear Constraint Handling Techniques for Generating Set Search Method for Well Placement Optimization Problem on FieldOpt Platform
Hardeman A geometric modeling framework for the analysis of fractured rock volumes
Eigestad Reservoir simulation with imposed flux continuity conditions on heterogeneous and anisotropic media for general geometries, and the inclusion of hysteresis in forward modeling
Garrett Special requirements for computer simulation of injection projects
Wang et al. Solver Preconditioning Using the Combinatorial Multilevel Method

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20171201

Termination date: 20200208

CF01 Termination of patent right due to non-payment of annual fee