多重介质油藏数值模拟异构并行算法研究

陈元科，张冬梅，崔书岳，张宇洋

(1.中国地质大学(武汉)计算机学院，湖北 武汉 430074；2.中国石化石油勘探开发研究院，北京100081)

1 引言

塔河油田是古岩溶与多期构造运动共同作用下生成的缝洞型碳酸盐岩油藏，储集体类型多样，具有非均质性较强等特点[1]。油藏数值模拟是定量评价地层流体流动规律的主要方法[2]，缝洞型油藏的储渗空间主要由溶蚀孔洞和裂缝组成，单元间连通性较差、油气水三相关系复杂[3 - 7]。针对该类复杂油藏，严侠等[8]提出基于离散缝洞网络模型和模拟有限差分法的刻画方法；张冬丽等[9]提出刻画不同尺度裂缝的数值模拟方法；康志江等[10]针对复杂油藏介质建立渗流耦合模型，形成多重介质数模的体系基础；邸元等[11]基于有限体积法采用非结构化网格建立多孔介质多相流数值模拟方法。KarstSim是中国石油化工股份有限公司研发的软件，能针对缝洞型碳酸盐岩油藏进行有效的数值模拟[11]，本文基于KarstSim开展并行多重介质油藏数值模拟算法研究。

缝洞型油藏复杂的多孔介质和油气水三相关系使数值模拟计算量巨大，为提高模拟效率必须开展并行计算研究[12]。目前主流的并行技术有3种：基于中央处理器CPU(Central Processing Unit)的分布式并行、基于图形处理器GPU(Graphic Processing Unit)的统一计算设备架构CUDA(Compute Unified Device Architecture)并行和基于CPU-GPU的异构并行。异构并行结合CUDA和分布式并行的技术特点，利用CPU实现粗粒度进程并行同时采用GPU实现细粒度线程并行。并行油藏数值模拟研究多集中在雅各比矩阵组装和线性方程组求解，雅各比矩阵组装分别计算各网格的流动项，线性方程求解多采用Krylov子空间求解方法中的共轭梯度CG(Conjugate Gradient)及改进算法，计算量大。Yang等[13]提出基于混合有限元法、隐式向后欧拉法和Newton-Krylov法的可拓展全隐式求解器，实现两相流油藏数值模拟在多核集群上的分布式并行；李毅等[14]通过网格分区、多CPU高性能计算库Aztec和网格分解软件METIS实现分布式并行缝洞型油藏数值模拟；Werneck等[15]基于分布式并行实现了并行CG法及改进算法并应用于单相流油藏模拟；李政等[16,17]针对油藏模拟中雅各比矩阵的稀疏模式，设计了一种块混合结构稀疏存储格式BHYB和线程组并行策略，并提出了基于CPU-GPU异构体系的并行BCPRP(Block Correction Proceduce via Reconstruction)预处理方法，在黑油模型上得到了较好的求解速度。基于分布式并行的油藏数值模拟的研究多通过网格划分实现并行方程求解，但网格划分通常会丢失一些流动信息；基于GPU的线性方程组求解应用非常广泛，结合矩阵预处理方法对矩阵信息的影响较少，比直接进行网格划分更为准确。

传统黑油模型难以用于缝洞型油藏数值模拟，本文针对多重介质模型数值模拟计算耗时的问题，提出一种基于CPU-GPU异构并行架构的数值模拟算法。基于GPU实现并行雅各比矩阵组装和线性方程组求解法，进一步研究多CPU下具有可扩展性的分布式并行数值模拟算法，提高缝洞型油藏数值模拟效率。

2 多重介质模型

多重介质缝洞型油藏中的介质主要由孔隙、裂缝和溶洞组成，裂缝是流动的主要通道，基质的孔隙度和渗透率较低，流体有油、气、水三相，其中气体也能溶解于油相，各相流体流动服从达西定律。等效模型划分为基质和介质，流动发生在基质之间的介质上，等效模型如图1所示[18]。

Figure 1 Equivalent model图1 等效模型

孔隙中油、气、水各组分的连续性方程如式(1)[10]所示：

(1)

裂缝、溶洞中油、气、水各组分的连续性方程如式(2)[10]所示：

(2)

孔隙、裂缝与溶洞中油、气、水各组分的运动方程如式(3)[10]所示：

(3)

其中,K表示绝对渗透率；Krβ代表各相具有的相对渗透率；μβ表示各相的黏度；pβ表示各相的压力；g表示重力加速度；D表示介质所在深度。

由连续性方程和运动方程可得孔隙、裂缝和溶洞中油、气、水各组分的质量平衡方程如式(4)所示：

(4)

其中,γog=(ρo+ρdg)g,γo=ρog,γw=ρwg。

各相饱和度和毛细管压力辅助方程如式(5)[11]所示：

(5)

其中,pcow、pcog分别表示油水、气油间的毛细管压力。

将介质划分为不同区块，各区块再划分为溶洞、裂缝、基质3个单元。各单元之间、相连区块之间各单元的流体流动有连续性方程：

(6)

其中qβ是β相的源汇项。

用有限体积法对式(6)进行空间离散得到非线性方程[11]：

(7)

在满足达西定律的流动中，流动项可表示为式(8)[11]所示：

(8)

(9)

γi,j定义为单元i和单元j之间的连通度：

(10)

Ψ是单元中流体的流动势：

(11)

式(10)中,Aij表示单元i和单元j交界面的面积；di,dj分别表示单元i和单元j的中心点到2个单元交界面的距离；Kij+1/2表示单元i和单元j的平均渗透率；式(11)中，Li表示单元i的中心点到高度参考面的距离。

3 异构并行多重介质油藏数值模拟算法

多重介质模型的迭代过程包括雅各比矩阵组装和线性方程组求解2个主要耗时步骤。本文基于消息传递接口MPI(Message Passing Interface)在多个GPU上实现雅各比矩阵组装，各分布式结点分别调用GPU计算分配任务；方程求解部分通过CUDA并行求解库实现，该库包含多种在GPU上实现的并行矩阵预处理方法和并行基础线性方程求解方法，如不完全LU分解法(incomplete Lower-Upper factorization)和稳定双共轭梯度法BiCGStab(BiConjugate Gradient Stabilized method)。

3.1 基于CPU-GPU异构的多重介质油藏数值模拟算法

本文对雅可比矩阵组装和线性方程组求解展开并行工作，进而降低数值模拟总耗时，基于分布式并行雅可比矩阵组装算法和GPU并行BiCGStab法的线性方程组求解方法实现异构并行多重介质油藏数值模拟算法。通过MPI为各CPU分配计算任务，各CPU通过GPU实现雅可比矩阵组装。GPU并行时将数据从主机内存传输到设备内存，针对数据交换延迟，利用GPU的异步发送功能将数据划分为多个数据段，发送的同时执行计算任务，以减少时间开销。调用GPU完成计算后各CPU通过MPI通信函数传输结果。CUDA线性方程求解库不支持多GPU并行且存在数据依赖关系，故将各自组装的矩阵数据规约传输到主CPU，由主CPU执行GPU并行求解算法。根据收敛情况，将解向量广播到其他各CPU用于更新流动项。多重介质油藏数值模拟异构并行算法流程如图2所示。

Figure 2 Flow chart of heterogeneous parallel algorithm of multi-media oil numerical simulation图2 多重介质油藏数值模拟异构并行算法流程图

3.2 基于GPU并行的矩阵组装和方程求解

缝洞型油藏非均质性强，多重介质模型建立的线性方程组系数矩阵非常稀疏，为提高求解效率，使用块压缩稀疏行BSR(Block Compressed Sparse Row)格式对稀疏矩阵进行压缩。各行组的非零元素表示为：

ABlock=Ai,j

(12)

其中,ABlock表示i行j列非零位置的3×3子矩阵。

本文基于CUDA实现GPU并行雅可比矩阵组装和GPU并行预处理BiCGStab法。在组装雅可比矩阵的过程中，要计算各网格通过介质与相连网格产生的流动信息，单元中共有n个网格对应矩阵中的n个行组。矩阵中的行组被划分为油气水三相，子矩阵结构如图3所示。对矩阵行组的各非零子矩阵调用3×3个GPU线程实现并行计算。

Figure 3 Jacobi matrix structure图3 雅可比矩阵结构

预处理BiCGStab法通过调整矩阵结构提高迭代求解的收敛性能。本文通过cuSparse库对矩阵进行不完全ILU分解，得到预处理后的上三角矩阵和下三角矩阵，再对线性方程组的解空间进行搜索，其中矩阵向量乘是主要的计算步骤。GPU实现并行矩阵向量乘的过程如图4所示，GPU分配矩阵向量乘任务对应的数据，各线程调用一个计算单元执行一次乘操作，具体通过CUDA库函数实现迭代求解，收敛得到当前时刻三相流动项的变化值。

Figure 4 Parallel matrix vector multiplication图4 并行矩阵向量乘

3.3 基于CPU并行的雅可比矩阵组装

为进一步提高矩阵组装并行效率，基于异构并行架构实现多CPU调用多GPU。统计所有非零位置子矩阵的数量M，为N个CPU核心分配相近数量的子矩阵计算任务，再由N个CPU核心分别调用对应GPU实现负载均衡，如图5所示。

Figure 5 Calculation task partition图5 计算任务分配

计算任务所需数据从异构框架中的主机端发送到设备端，从设备端开始GPU并行。按求解执行顺序将数据划分为多个部分，先发送第1部分，后续发送的同时执行求解任务，以减少时间开销。调用MPI数据交换函数时进程需等待，调用次数越多总等待时间越长。本文将通信数据统一存放到临时数组，只需调用一次数据交换函数，降低了时间开销。各结点完成计算任务后，通过MPI异步数据传输函数MPI_REDUCE将各CPU计算的矩阵数据传输回主CPU；用REDUCE函数同步数据以减少各CPU在数据传输后的缓冲区读取次数和时间开销；数据同步到主CPU后，再由主CPU开启CUDA核函数进行GPU线性方程组求解。

4 数值实验

实验将本文提出的并行算法和串行油藏数值模拟软件的运行结果进行比较。对塔河油田3个不同规模的问题进行模拟计算，验证并行算法在模拟结果上的准确性和可扩展性。并行平台由12个主频为3.4 GHz的4核心CPU(Intel core i5-7500)和12块NVIDIA GeForceGTX1060组成，搭载Linux操作系统，各结点安装MPICH2和CUDA9.1。

4.1 并行模拟结果对比

塔河油田S04单元模型网格如图6所示，该模型网格数是399 158，连接数1 064 405，稀疏矩阵非零元素数量2 527 968。除裂缝和溶洞介质外，油藏基质部分渗透率低，各单元有不同孔隙度和绝对渗透率，相对渗透率和毛细管压力曲线相同。S04单元的基本参数如表1所示。

Figure 6 Grid of unit S04图6 S04单元网格

Table 1 Basic parameters of unit S04 表1 S04单元基本参数

本文采用3个的总计12个CPU核心对S04单元开展实验，对比模拟结果的累积产油量、单元中TK455井的含水率和油藏压力。图7所示为累积产油量曲线，图8所示为TK455井的含水率曲线，图9所示为油藏压力曲线。

Figure 7 Comparison of simulated remaining oil quality of unit S04图7 S04单元模拟累积产油量对比

Figure 8 Comparison of the moisture ratio of well TK455图8 TK455井模拟含水率对比

Figure 9 Comparison of bottom pressure data of unit S04图9 S04单元模拟油藏压力数据对比

对比并行算法与串行模拟曲线可知，模拟结果几乎相同，本文提出的并行算法能正确地模拟计算S04单元。

4.2 并行加速效果对比

并行加速效果由加速比进行衡量，包括矩阵组装加速比、求解加速比和总体加速比。对3个不同规模的单元数据应用异构并行算法开展数值模拟。615B规模相对较小，单元数约为4万，非零元素数量约22万，线性方程组系数矩阵的稀疏度较低;S80规模和线性方程组规模都比615B单元更大，单元数约9万，非零元素数量约48万；S04单元规模最大，单元数接近40万，矩阵非零元素数量达到252万。615B单元与S80单元的基本参数如表2所示，其中稀疏程度定义为：非零元素数量/(单元数)2。

Table 2 Basic parameters of unit 615B and S80 表2 615B单元与S80单元的基本参数

本文调用2～12个CPU核心，各CPU核心分别对应MPI环境中的一个进程，并分别调用一块GPU。分别对615B、S80和S04单元开展异构并行数值模拟实验，结果如图10所示。随着CPU核心数增加，各单元线性方程组组装加速比接近线性；当CPU核心数超过8时，算法对615B单元和S04单元仍有较好的可扩展性，但对S80单元的可扩展性明显减弱。

Figure 10 Matrix assembly parallel scalability curve of each unit图10 各单元矩阵组装并行可扩展性曲线

对比线性方程组求解部分的预处理耗时、迭代求解耗时、平均各迭代耗时和求解总耗时。串行软件仅调用一个CPU核心，并行算法调用12个CPU核心和12个GPU。实现超过2倍的预处理加速比和3倍以上的迭代平均加速比，总体加速效果达到预期，结果如表3所示。

Table 3 Equation system solving speedup of each unit表3 各单元方程组求解加速比

配置12个CPU核心3个单元数值模拟的总体并行加速比如表4所示，各单元数据都实现了2倍以上的总体加速比。

综上所述，本文设计的异构并行算法能够充分实现加速效果，模拟结果也与串行油藏数值模拟软件结果基本一致。

Table 4 Overall parallel speedup of each unit表4 各单元总体并行加速比

5 结束语

本文采用多重介质模型刻画缝洞型油藏，但数值模拟耗时巨大。为提高多重介质模型计算速度，本文基于异构并行架构，针对多重介质模型耗时最高的线性方程组组装和求解部分，设计多重介质油藏数值模拟异构并行算法。

本文针对塔河油田缝洞型油藏3个不同规模的单元615B、S80、S04进行并行数值模拟实验，615B和S04单元的实验结果表明，多重介质油藏数值模拟异构并行算法具有较好的准确性，能有效提高缝洞型油藏数值计算速度。S80单元的模拟结果表明，并行矩阵组装的可扩展性有优化的空间。未来将对超大规模网格的油藏数据开展实验，优化多重介质油藏数值模拟异构并行算法，提高计算效率。