本发明涉及化学驱精细油藏数值模拟，特别是涉及到油藏数值模拟中压力方程的并行直接求解方法。
背景技术：
1、油藏数值模拟是利用计算机研究油气藏中多相流体渗流规律的数值计算方法，用以解决油藏开发过程中难以求解的复杂渗流及工程问题，是剩余油精细预测和高效挖潜的关键技术之一。目前，我国大部分可开采油田都已进入开发后期，油田的综合含水率居高不下。普通的油藏数值模拟已经不能满足当前开采需求，为了精确了解剩余油分布，提高石油开采率，需要采用精细油藏数值模拟进一步优化，然而精细油藏数值模拟的网格数十分庞大，达到几百万甚至数千万个网格，随着模型规模的增大，模拟运算时间越来越长，模拟速度成为数值模拟技术广泛应用的瓶颈。用普通的计算机或工作站进行模拟需要耗费几天甚至几个月时间，或者因为内存不足根本无法模拟。
2、油藏数值模拟一般包括采集数据、源量分配、压力预处理、求解压力方程、求解饱和度和组分浓度、分析剩余油和综合含水率几个模块。其中，求解压力方程是目前油藏数值模拟中最耗时的模块，占比高达55％以上甚至更高。为了解决精细油藏数值模拟因网格数多而导致耗时长的问题，本发明着重研究利用mpi实现对压力方程模块的并行求解，提高精细油藏数值模拟的计算效率。
3、在申请号：cn201510135020.4的中国专利申请中，涉及到一种化学驱油藏数值模拟方法，该化学驱油藏数值模拟方法包括:步骤1，进行网格模型优化设计，建立静态模型；步骤2，在建立的静态模型的基础上，建立动态模型；以及步骤3，在建立好动态模型的基础上，分时段对参数进行调整，实现分时段历史拟合，不同时段重点调整的参数不同，不同层次拟合精度要求不同。该化学驱油藏数值模拟方法改进数值模拟软件应用水平，提高油藏描述精度，为油藏开发后期提供更加精确的技术支持。
4、在申请号：cn201510026314.3的中国专利申请中，涉及到一种基于并行自调整差分进化的三元复合驱优化方法。该方法首先以净现值为性能指标，以油、水两相的渗流方程和碱、表面活性剂、聚合物的对流扩散吸附方程为支配方程，以驱替剂协同方程、物化参数代数方程以及驱替剂边界条件作为约束方程，以各驱替剂注入浓度为优化变量，建立优化模型。其次，为了改善传统差分进化算法的全局和局部搜索能力，提高收敛速度，避免进入局部最优，通过对进化过程引入并行策略，对缩放因子进行自调整，提出一种新的并行自调整差分进化算法。最后利用该算法优化三元复合驱模型的注入浓度，得到最优注采策略。
5、在申请号：cn201710702432.0的中国专利申请中，涉及到一种对缝洞型油藏剩余油进行快速分析的方法及系统，其包括建立缝洞型油藏的数学模型；对缝洞型油藏进行网格单元的划分，基于数学模型建立数值模型；对数值模型进行迭代求解，并根据用于描述缝洞型油藏的流体压力和饱和度的分布对缝洞型油藏的剩余油进行分析；其中，在对数值模型进行迭代求解时的每一次迭代步骤中，包括：基于网格单元的特性，从数值模型中提取油相压力方程；基于代数多重网格算法对所述油相压力方程进行处理，以得到第一压力值；以第一压力值作为对所述数值模型进行迭代时的初始解。采用该方法，可对缝洞型油藏剩余油进行快速有效分析。
6、在申请号：cn201710253502.9的中国专利申请中，涉及到一种油藏渗流数值模拟方法及装置，其中方法包括：获取非均质砂岩油藏各网格区域的初始状态信息、初始状态信息中各个参数对应的变化函数以及非均质砂岩油藏各个渗透率区间的属性信息；将上述参数代入预设的油藏渗流数值模型，得到包含各个渗透率区间对应的渗流控制方程的渗流控制方程组；结合渗流控制方程组的定解条件对渗流控制方程组进行求解，得到非均质砂岩油藏生产一定时间后各个参数的参数值，进而准确描述非均质砂岩油藏生产预设时间后的剩余油分布特征并制定合理的挖潜方式；该发明充分考虑了各储层的渗流特征的差异，减小了评估结果与实际情况之间的误差。
7、以上现有技术均与本发明有较大区别，未能解决我们想要解决的技术问题，为此我们发明了一种新的油藏数值模拟中压力方程的并行直接求解方法。
技术实现思路
1、本发明的目的是提供一种可以解决随着油藏模拟网格规模增大，模拟耗时过长的问题的油藏数值模拟中压力方程的并行直接求解方法。
2、本发明的目的可通过如下技术措施来实现：油藏数值模拟中压力方程的并行直接求解方法，该油藏数值模拟中压力方程的并行直接求解方法包括：
3、步骤1，根据油藏数学模型得到压力方程，采集相关数据，对数据做预处理；
4、步骤2，根据井模型进行源量分配；
5、步骤3，进行压力预处理；
6、步骤4，求解压力方程；
7、步骤5，计算油藏数值模拟的组分浓度、相饱和度；
8、步骤6，分析剩余油和综合含水率，确定矿场实施方案。
9、本发明的目的还可通过如下技术措施来实现：
10、在步骤1中，油藏数值模拟基本的数学模型基于如下假设：油藏中局部热力学平衡、固相不流动、岩石和流体微可压缩、费克弥散、理想混合、流体渗流满足达西定律；在基本假设下，压力求解的精度与速度对油藏数值模拟的整体求解具有先决性的，压力方程为如下形式：
11、
12、其中，φ是参考压力下的孔隙度，ct为总压缩系数，t是时间，p为压力，k为岩石的渗透率，λrw、λro是水相、油相的流度，γw、γo是水相、油相的相对密度，pc是毛管压力，d是深度，div是散度，▽是梯度，qw和qo分别为水、油的源汇流量。
13、在步骤1中，采集初始压力，初始压力下的密度，粘度，孔隙度，对数据做预处理，计算传导率和渗透率。
14、在步骤2中，注入井按注入速度分配到各层，对压力方程的隐式求解添加井约束。
15、在步骤3中，进行压力预处理，计算压力方程系数和右端项；压力方程按照七点中心差分格式离散，离散之后，用油水两相传输系数之和确定六个方向的离散系数，同时利用综合压缩系数计算时间导数相系数，形成主对角元素；毛管力项、重力项显式处理置于右端；系数矩阵为七对角对称性镶边矩阵a、右端项b。
16、在步骤4中，基于lu分解的直接求解方法并行求解压力线性代数方程组ap＝b；首先把系数矩阵a分解为下三角矩阵l和上三角矩阵u，即并行求解矩阵l和矩阵u使得a＝lu；然后分别并行三角求解uy＝b和lp＝y求得压力p。
17、步骤4包括：
18、步骤41，存储代数方程组；
19、步骤42，分发代数方程组；
20、步骤43，缩放代数方程组；
21、步骤44，对方程组进行行排序；
22、步骤45，对方程组进行列排序；
23、步骤46，进行代数方程组符号分解；
24、步骤47，进行代数方程组并行数值分解及三角求解；
25、步骤48，收集求解得到的压力p。
26、在步骤41中，为了节省存储空间，采用行压缩形式存储系数矩阵a；在行压缩形式中，按行存储矩阵a中所有非零元素，使用数组nzval[]存储非零元素的数值，使用数组colind[]存储非零元素所在的列号，数组rowptr[]存储非零元素所在行的信息；由于非零元素是按行存储的，矩阵中的一行元素是连续的，数组rowptr[]只需存储每一行的第一个非零元素的位置；此外，还需要额外存储矩阵a的维数n和非零元个数nnz。
27、在步骤42中，为了使各处理器并行lu分解，把代数方程组按行分块并且分发给各个处理器；若使用的处理器个数为p个，则把矩阵a划分成a1,a2,...,ap共p块子矩阵，右端项b以相同形式划分成b1,b2,...,bp共p块子向量，然后分发到各个处理器，第i个处理器存储行压缩形式的子矩阵ai和数组形式的bi。
28、在步骤43中，为了避免lu分解过程出现浮点异常错误，按行和列对系数矩阵ai和右端项bi进行平衡缩放使得系数矩阵的各行各列元素的数量级相同；首先计算行缩放对角矩阵dr与列缩放对角矩阵dc，然后使用行缩放矩阵dr与列缩放矩阵dc缩放代数方程组的系数矩阵ai和右端项bi，得到ai←dr·ai·dc，bi←dr·bi·dc。
29、在步骤44中，为避免lu分解过程因主对角线元素过小导致的不稳定性，对系数矩阵ai和右端项bi进行行排序使得系数矩阵的各个主对角线元素最大化；通过列选主元的方式，依次比较矩阵每行的对角元素与对角线下各行同列元素的大小，并通过行交换的方式使得每行的对角元素大于对角线下各行同列元素；通过行排序对角矩阵pr记录矩阵进行行交换后各行的顺序，然后使用该矩阵重新排序代数方程组的系数矩阵ai和右端项bi，得到a←pr·a，b←pr·b。
30、在步骤45中，为了提高lu分解过程的效率，对系数矩阵ai和右端项bi进行列排序使得系数矩阵分解时产生的填充元素的个数最小化；使用多层嵌套分割排序算法，重新排序以系数矩阵为邻接矩阵的图g(a)的每个顶点，通过列排序对角矩阵pc记录排序后各顶点的顺序，然后使用该矩阵以对称的方式重新排序代数方程组的系数矩阵ai和右端项bi，得到
31、在步骤46中，为了得到分解ai＝liui中li和ui矩阵的非零元素的个数和位置，并为其分配内存，对系数矩阵进行符号分解；使用深度优先搜索遍历以系数矩阵为邻接矩阵的图g(a)的每个顶点，得到li和ui矩阵中非零元素的个数与位置。
32、在步骤47中，使用高斯消元法并行计算ai＝liui中li和ui矩阵非零元素的数值，然后使用回代法分别进行并行三角求解uiyi＝bi和lipi＝yi求得压力pi；该步骤的并行性在于步骤42中代数方程组的分发，由于代数方程组的各块存储在不同的处理器上面，对于分解与求解过程的各个计算，各个处理器仅计算自己存储的部分所以，该步骤的计算是由多个处理器共同计算完成的。
33、在步骤48中，收集各处理器求解得到的压力p1,p2,...,pp共p块子向量，并按下标顺序合并得到压力p。
34、在步骤5，利用油藏数学模型计算油藏数值模拟的组分浓度、相饱和度。
35、在步骤6，油藏开发研究人员通过分析剩余油和综合含水率等数据，确定矿场实施方案。
36、本发明中的油藏数值模拟中压力方程的并行直接求解方法，针对油藏数值模拟向大规模、精细化的方向发展而带来的油藏网格规模增加、模拟时间过长的问题，油藏数值模拟可以在高性能计算机上进行，对其中的压力求解并行化，可以缩短压力求解的时间，提高压力方程的求解精度与速度，压力方程求解越快，计算越准确，接下来的组分浓度、相饱和度等的计算才能更加快速准确，进而减少油藏数值模拟的时间，消除油藏数值模拟过程的性能瓶颈，为油藏进一步提高采油率提供技术支持。
1.油藏数值模拟中压力方程的并行直接求解方法，其特征在于，该油藏数值模拟中压力方程的并行直接求解方法包括：
2.根据权利要求1所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，在步骤1中，油藏数值模拟基本的数学模型基于如下假设：油藏中局部热力学平衡、固相不流动、岩石和流体微可压缩、费克弥散、理想混合、流体渗流满足达西定律；在基本假设下，压力求解的精度与速度对油藏数值模拟的整体求解具有先决性的，压力方程为如下形式：
3.根据权利要求2所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，在步骤1中，采集初始压力，初始压力下的密度，粘度，孔隙度，对数据做预处理，计算传导率和渗透率。
4.根据权利要求1所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，在步骤2中，注入井按注入速度分配到各层，对压力方程的隐式求解添加井约束。
5.根据权利要求2所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，在步骤3中，进行压力预处理，计算压力方程系数和右端项；压力方程按照七点中心差分格式离散，离散之后，用油水两相传输系数之和确定六个方向的离散系数，同时利用综合压缩系数计算时间导数相系数，形成主对角元素；毛管力项、重力项显式处理置于右端；系数矩阵为七对角对称性镶边矩阵a、右端项b。
6.根据权利要求5所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，在步骤4中，基于lu分解的直接求解方法并行求解压力线性代数方程组ap＝b；首先把系数矩阵a分解为下三角矩阵l和上三角矩阵u，即并行求解矩阵l和矩阵u使得a＝lu；然后分别并行三角求解uy＝b和lp＝y求得压力p。
7.根据权利要求6所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，步骤4包括：
8.根据权利要求7所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，在步骤41中，为了节省存储空间，采用行压缩形式存储系数矩阵a；在行压缩形式中，按行存储矩阵a中所有非零元素，使用数组nzval[]存储非零元素的数值，使用数组colind[]存储非零元素所在的列号，数组rowptr[]存储非零元素所在行的信息；由于非零元素是按行存储的，矩阵中的一行元素是连续的，数组rowptr[]只需存储每一行的第一个非零元素的位置；此外，还需要额外存储矩阵a的维数n和非零元个数nnz。
9.根据权利要求8所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，在步骤42中，为了使各处理器并行lu分解，把代数方程组按行分块并且分发给各个处理器；若使用的处理器个数为p个，则把矩阵a划分成a1,a2,...,ap共p块子矩阵，右端项b以相同形式划分成b1,b2,...,bp共p块子向量，然后分发到各个处理器，第i个处理器存储行压缩形式的子矩阵ai和数组形式的bi。
10.根据权利要求9所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，在步骤43中，为了避免lu分解过程出现浮点异常错误，按行和列对系数矩阵ai和右端项bi进行平衡缩放使得系数矩阵的各行各列元素的数量级相同；首先计算行缩放对角矩阵dr与列缩放对角矩阵dc，然后使用行缩放矩阵dr与列缩放矩阵dc缩放代数方程组的系数矩阵ai和右端项bi，得到ai←dr·ai·dc，bi←dr·bi·dc。
11.根据权利要求10所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，在步骤44中，为避免lu分解过程因主对角线元素过小导致的不稳定性，对系数矩阵ai和右端项bi进行行排序使得系数矩阵的各个主对角线元素最大化；通过列选主元的方式，依次比较矩阵每行的对角元素与对角线下各行同列元素的大小，并通过行交换的方式使得每行的对角元素大于对角线下各行同列元素；通过行排序对角矩阵pr记录矩阵进行行交换后各行的顺序，然后使用该矩阵重新排序代数方程组的系数矩阵ai和右端项bi，得到a←pr·a，b←pr·b。
12.根据权利要求11所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，在步骤45中，为了提高lu分解过程的效率，对系数矩阵ai和右端项bi进行列排序使得系数矩阵分解时产生的填充元素的个数最小化；使用多层嵌套分割排序算法，重新排序以系数矩阵为邻接矩阵的图g(a)的每个顶点，通过列排序对角矩阵pc记录排序后各顶点的顺序，然后使用该矩阵以对称的方式重新排序代数方程组的系数矩阵ai和右端项bi，得到
13.根据权利要求12所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，在步骤46中，为了得到分解ai＝liui中li和ui矩阵的非零元素的个数和位置，并为其分配内存，对系数矩阵进行符号分解；使用深度优先搜索遍历以系数矩阵为邻接矩阵的图g(a)的每个顶点，得到li和ui矩阵中非零元素的个数与位置。
14.根据权利要求13所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，在步骤47中，使用高斯消元法并行计算ai＝liui中li和ui矩阵非零元素的数值，然后使用回代法分别进行并行三角求解uiyi＝bi和lipi＝yi求得压力pi；该步骤的并行性在于步骤42中代数方程组的分发，由于代数方程组的各块存储在不同的处理器上面，对于分解与求解过程的各个计算，各个处理器仅计算自己存储的部分所以，该步骤的计算是由多个处理器共同计算完成的。
15.根据权利要求14所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，在步骤48中，收集各处理器求解得到的压力p1,p2,...,pp共p块子向量，并按下标顺序合并得到压力p。
16.根据权利要求1所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，在步骤5，利用油藏数学模型计算油藏数值模拟的组分浓度、相饱和度。
17.根据权利要求1所述的油藏数值模拟中压力方程的并行直接求解方法，其特征在于，在步骤6，油藏开发研究人员通过分析剩余油和综合含水率这些数据，确定矿场实施方案。
