A numerical simulation method and system for fractured-vug type reservoirs

Abstract

The application discloses a numerical simulation method and system for a fractured-vug type oil reservoir. The method comprises the following steps: integrating a plurality of simulation sub-zones with distributed input into a target simulation region; performing grid division on the target simulation region to obtain a calculation task; selecting a calculation scheme from a distributed parallel calculation scheme, a shared memory parallel calculation scheme and a hybrid parallel calculation scheme; dividing the calculation task according to the selected calculation scheme to obtain a corresponding number of sub-tasks; performing all the sub-tasks in parallel; solving a nonlinear equation by using a predetermined fully implicit method, an implicit pressure explicit saturation method or an adaptive implicit method, calling a third-party linear equation solving library and solving the linear equation based on a selected precondition method; and distributing all the simulation results obtained. The system: the processor realizes the method when executing the computer program stored in the memory. According to the application, large-scale fine simulation of the fractured-vug type oil reservoir can be realized.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of reservoir numerical simulation, and more particularly, relates to a fracture-cave type reservoir numerical simulation method and system. BACKGROUND

[0002] Numerical simulation is an important analysis and research method for oilfield development evaluation and production plan design. The core of reservoir numerical simulation is to solve a set of nonlinear partial differential equations that can reflect the three-phase flow of oil, water and gas by numerical methods, so as to achieve the purpose of simulating the oilfield development history, predicting the future dynamic of the reservoir, adjusting and formulating the best oil production mode, etc. How to accurately and efficiently solve these nonlinear partial differential equations and shorten the total time required for model operation is one of the key problems that should be solved in the development of reservoir simulators.

[0003] Parallel computing technology is considered to be an effective way to accelerate the solution process of large and complex computing problems, and can meet the requirements of large-scale fine simulation of reservoirs. At present, the main commercial reservoir simulation software improves the equation solving method and adopts domain decomposition and distributed parallel computing to improve the computing efficiency. In recent years, many commercial reservoir simulation software has been introduced and related research has been reported internationally. For example, Saudi Aramco first realized the simulation of a million grid scale in 2000, broke through the simulation of a billion grid scale in 2010, and realized the unimaginable simulation of a trillion grid scale (TeraPOWER, 2016) using 150,000 computing cores in 2016. At present, the reservoir simulators that have broken through the simulation of a billion grid scale internationally mainly include INTERSECT of Schlumberger, GigaPOWERS and TeraPOWERS of Saudi Aramco, ECHELON of Stone Ridge Technology Company and the self-use simulator of ExxonMobil. Among them, except for ECHELON which is a complete GPU simulator, the others are distributed memory parallel simulators based on MPI. However, the above simulators do not take into account the particularity of fracture-cave type reservoirs and cannot solve the problem of large-scale fine simulation of fracture-cave type reservoirs. SUMMARY

[0004] The present application aims to solve the problem that the existing reservoir simulators cannot perform large-scale fine simulation of fracture-cave type reservoirs.

[0005] In order to achieve the above-mentioned purpose, the present application provides a fracture-cave type reservoir numerical simulation method and system.

[0006] According to a first aspect of the present application, a fracture-cave type reservoir numerical simulation method is provided, which comprises the following steps:

[0007] Integrate a plurality of simulation sub-zones distributedly inputted into a target simulation region;

[0008] Mesh the target simulation region to obtain a calculation task;

[0009] Select a calculation scheme from a predetermined distributed parallel calculation scheme, a shared memory parallel calculation scheme and a hybrid parallel calculation scheme;

[0010] Divide the calculation task according to the selected calculation scheme to obtain a corresponding number of sub-tasks;

[0011] Parallelly execute all the sub-tasks: in this process, a predetermined fully implicit method, an implicit pressure explicit saturation method or an adaptive implicit method is used to solve a nonlinear equation, a third-party linear equation solving library is called and a selected precondition method is used to solve a linear equation;

[0012] Distribute all simulation results obtained by executing the sub-tasks.

[0013] Preferably, the distributed parallel calculation scheme specifically comprises:

[0014] Mesh the target simulation region after meshing to obtain a plurality of sub-domains;

[0015] Distribute the plurality of sub-domains to a corresponding plurality of calculation units respectively for parallel calculation;

[0016] Interaction between the plurality of calculation units is realized based on MPI communication.

[0017] Preferably, the shared memory parallel calculation scheme specifically comprises:

[0018] Distribute the calculation task to a calculation unit;

[0019] The calculation unit parallelly calculates the distributed calculation task based on a plurality of minimum calculation sub-units therein;

[0020] The plurality of minimum calculation sub-units share a memory space of the calculation unit.

[0021] Preferably, the hybrid parallel calculation scheme specifically comprises:

[0022] Mesh the target simulation region after meshing to obtain a plurality of sub-domains;

[0023] Distribute the plurality of sub-domains to a corresponding plurality of calculation units respectively for parallel calculation;

[0024] Interaction between the plurality of calculation units is realized based on MPI communication.

[0025] The computing unit performs parallel computation on the assigned sub-domain based on a plurality of minimum computing sub-units therein;

[0026] The plurality of minimum computing sub-units share the memory space of the computing unit.

[0027] Preferably, the constraint condition of the domain decomposition comprises:

[0028] The difference of the number of grids in any two sub-domains is less than a predetermined difference reference value;

[0029] The total number of grid connections on all sub-domain boundaries is minimum.

[0030] Preferably, in the process of the mesh partitioning of the target simulation region, further comprises:

[0031] Firstly globally numbering the grids in the target simulation region;

[0032] After the domain decomposition of the mesh-partitioned target simulation region to obtain a plurality of sub-domains, further comprises:

[0033] Locally numbering the grids in each sub-domain;

[0034] The numerical simulation method of the fractured-vuggy reservoir further comprises:

[0035] Secondly globally numbering the grids in the target simulation region for realizing the interaction between the computing units and the calling of the third-party linear equation solving library.

[0036] Preferably, the self-adaptive implicit method specifically comprises:

[0037] Judging whether the variation of the saturation of the target component in the target grid in a time step is greater than a predetermined variation reference value, if yes, using the fully implicit method to solve the nonlinear equation corresponding to the target grid, if no, using the implicit pressure explicit saturation method to solve the nonlinear equation corresponding to the target grid;

[0038] Alternatively, the self-adaptive implicit method specifically comprises:

[0039] Obtaining a CFL condition coefficient;

[0040] Judging whether the CFL condition coefficient is less than a predetermined maximum value of the CFL condition coefficient, if yes, using the implicit pressure explicit saturation method to solve the nonlinear equation corresponding to the target grid, if no, using the fully implicit method to solve the nonlinear equation corresponding to the target grid;

[0041] Or, when solving a nonlinear equation, the full implicit method is used to solve the grid corresponding nonlinear equation of the well or source-sink term; at the first time step of simulation calculation, the implicit pressure explicit saturation method is used to solve the nonlinear equation corresponding to the grid of other types of grid except the well or source-sink grid.

[0042] Preferably, the pre-condition method is selected in the following manner:

[0043] When the full implicit method is used to solve a linear equation, the CPR method is selected.

[0044] When the adaptive implicit method is used to solve a linear equation, the domain decomposition method or the ILU(0) method is selected, and the Jacob method is selected when the domain decomposition method fails.

[0045] Preferably, the distributed input of multiple simulation sub-zones and the distributed output of simulation results are realized based on a predetermined distributed input and output scheme.

[0046] The distributed input and output scheme comprises:

[0047] The target simulation region is divided into multiple simulation sub-zones.

[0048] Each simulation sub-zone is meshed to obtain a corresponding sub-model: the meshing manner on both sides of the boundary of two adjacent sub-models is consistent, and a corresponding virtual grid layer is arranged on both sides of the interface between the two adjacent sub-models and the corresponding entity grid, which serves as a marker layer in the subsequent sub-model merging process.

[0049] The obtained multiple sub-models are distributedly input into a numerical simulator.

[0050] The numerical simulator merges the input multiple sub-models to obtain the target simulation region: two adjacent sub-models are coordinate-matched through the pre-set virtual grid layer to realize merging, and all virtual grid layers are deleted after all sub-models are merged.

[0051] In response to an input simulation result output request, the numerical simulator outputs all simulation results in a corresponding distributed output manner.

[0052] According to a second aspect of the present application, a fractured-vuggy reservoir numerical simulation system is provided, which comprises a processor and a memory, and the processor realizes any of the above-mentioned fractured-vuggy reservoir numerical simulation methods when executing a computer program stored in the memory.

[0053] The present application has the following advantages:

[0054] The fracture-cave type oil reservoir numerical simulation method of the present application firstly integrates multiple simulation sub-zones distributedly input into a target simulation region; secondly, grid partition is performed on the target simulation region to obtain a calculation task; thirdly, a calculation scheme is selected from among a predetermined distributed parallel calculation scheme, a shared memory type parallel calculation scheme and a hybrid parallel calculation scheme; fourthly, the calculation task is divided according to the selected calculation scheme to obtain a corresponding number of sub-tasks; fifthly, all the sub-tasks are executed in parallel; and finally, all simulation results obtained by executing the sub-tasks are distributedly output.

[0055] The fracture-cave type oil reservoir numerical simulation method of the present application is directed to complex fracture-cave type carbonate rock reservoirs and comprehensively considers numerical simulation of multiple media of Darcy seepage and non-Darcy seepage in pores, fractures and fracture-caves. Due to the particularity of the reservoir of the fracture-cave type oil reservoir and the complexity of the seepage process, the fracture-cave type oil reservoir numerical simulation method of the present application uses domain decomposition technology and realizes simulation of a model of a grid scale of ten million through a distributed parallel calculation scheme, a shared memory type parallel calculation scheme and a hybrid parallel calculation scheme. The fracture-cave type oil reservoir numerical simulation method of the present application significantly improves the calculation efficiency of complex models by improving Newton iteration and solution of linear equations. In addition, the fracture-cave type oil reservoir numerical simulation method of the present application uses a distributed input / output scheme to solve the problem of difficult pre-processing and post-processing of massive data involved in large-scale models and has great significance for further development of such complex oil and gas fields.

[0056] The fracture-cave type oil reservoir numerical simulation system of the present application belongs to the same general inventive concept as the fracture-cave type oil reservoir numerical simulation method described above and therefore has the same beneficial effects as the fracture-cave type oil reservoir numerical simulation method described above, which will not be described again here.

[0057] Other features and advantages of the present application will be described in detail in the following specific embodiments. BRIEF DESCRIPTION OF DRAWINGS

[0058] The above and other objects, features and advantages of the present application will become more apparent from the following detailed description when taken in conjunction with the accompanying drawings in which like reference characters refer to like parts throughout the several views, in which:

[0059] Figure 1 An implementation flowchart of the fracture-cave type oil reservoir numerical simulation method according to Embodiment 1 of the present application is shown;

[0060] Figure 2 A domain decomposition schematic diagram of a two-dimensional simulation region according to Embodiment 1 of the present application is shown;

[0061] Figure 3 A domain decomposition schematic diagram of a three-dimensional simulation region according to Embodiment 1 of the present application is shown;

[0062] Figure 4 A structural schematic diagram of the hybrid memory structure according to Embodiment 1 of the present application is shown;

[0063] Figure 5 A structural schematic diagram of the large-scale model grid according to Embodiment 1 of the present application is shown;

[0064] Figure 6 A model diagram before S80 unit coarsening according to Embodiment 1 of the present application is shown;

[0065] Figure 7 A model block oil pressure distribution diagram according to Embodiment 1 of the present application is shown;

[0066] Figure 8 A comparative analysis diagram of the percentage of time consumed by each calculation link of the 50 million grid model and the 100 million grid model to the total time according to Embodiment 1 of the present application is shown. DETAILED DESCRIPTION

[0067] Preferred embodiments of the present application will be described in more detail below. Although the preferred embodiments of the present application are described below, it is to be understood that the present application can be carried out in various forms and should not be limited by the embodiments set forth herein. Rather, these embodiments are provided so that this application will be thorough and complete, and fully convey the scope of the present application to those skilled in the art.

[0068] Embodiment 1: Figure 1 An implementation flowchart of the fractured-vuggy reservoir numerical simulation method of the embodiments of the present application is shown. Referring to Figure 1 , the fractured-vuggy reservoir numerical simulation method of the embodiments of the present application includes the following steps:

[0069] Step S100, integrating a plurality of simulation partitions distributedly input into a target simulation region;

[0070] Step S200, performing grid partitioning on the target simulation region to obtain a calculation task;

[0071] Step S300, selecting a calculation scheme from among a predetermined distributed parallel computing scheme, a shared memory parallel computing scheme and a hybrid parallel computing scheme;

[0072] Step S400, dividing the calculation task according to the selected calculation scheme to obtain a corresponding number of subtasks;

[0073] Step S500, all subtasks are executed in parallel: in this process, a predetermined fully implicit method, an implicit pressure explicit saturation method or an adaptive implicit method is used to solve the nonlinear equation, a third-party linear equation solving library is called and a linear equation is solved based on a selected precondition method;

[0074] Step S600, all simulation result obtained by executing the subtasks is distributedly output.

[0075] Further, in step S300 of the embodiment of the present application, the distributed parallel computing scheme is specifically:

[0076] The target simulation region after grid partitioning is domain-decomposed to obtain a plurality of subdomains;

[0077] The plurality of subdomains are respectively distributed to a plurality of corresponding computing units for parallel computation;

[0078] Interaction between the plurality of computing units is realized based on MPI communication.

[0079] Further, in step S300 of the embodiment of the present application, the shared memory parallel computing scheme is specifically:

[0080] The computing task is distributed to a computing unit;

[0081] The computing unit performs parallel computation on the distributed computing task based on a plurality of minimum computing subunits therein;

[0082] The plurality of minimum computing subunits share the memory space of the computing unit.

[0083] Further, in step S300 of the embodiment of the present application, the hybrid parallel computing scheme is specifically:

[0084] The target simulation region after grid partitioning is domain-decomposed to obtain a plurality of subdomains;

[0085] The plurality of subdomains are respectively distributed to a plurality of corresponding computing units for parallel computation;

[0086] Interaction between the plurality of computing units is realized based on MPI communication;

[0087] The computing unit performs parallel computation on the distributed subdomain based on a plurality of minimum computing subunits therein;

[0088] The plurality of minimum computing subunits share the memory space of the computing unit.

[0089] Further, in the embodiment of the present application, the constraint condition of the domain decomposition includes:

[0090] The difference of the number of grids in any two sub-domains is less than a predetermined difference reference value;

[0091] The total number of grid connections on all sub-domain boundaries is minimized.

[0092] Further, in the process of the grid division of the target simulation region, the method further comprises:

[0093] The grids in the target simulation region are globally numbered for the first time;

[0094] After the domain decomposition of the grid-divided target simulation region, the method further comprises:

[0095] The grids in each sub-domain are locally numbered;

[0096] The numerical simulation method of the fractured-vuggy reservoir further comprises:

[0097] The grids in the target simulation region are globally numbered for the second time, so as to realize the interaction between the calculation units and the calling of the third-party linear equation solving library.

[0098] Further, in the step S500 of the embodiment of the method, the adaptive implicit method specifically comprises:

[0099] It is judged whether the variation of the saturation of the target component in the target grid in a time step is greater than a predetermined variation reference value, if yes, the full-implicit method is used to solve the nonlinear equation corresponding to the target grid, if no, the implicit pressure explicit saturation method is used to solve the nonlinear equation corresponding to the target grid;

[0100] Alternatively, the adaptive implicit method specifically comprises:

[0101] The CFL condition coefficient is obtained;

[0102] It is judged whether the CFL condition coefficient is less than a predetermined maximum value of the CFL condition coefficient, if yes, the implicit pressure explicit saturation method is used to solve the nonlinear equation corresponding to the target grid, if no, the full-implicit method is used to solve the nonlinear equation corresponding to the target grid;

[0103] Alternatively, when solving the nonlinear equation, the full-implicit method is used to solve the nonlinear equation corresponding to the grid of the well or source-sink term, and the implicit pressure explicit saturation method is used to solve the nonlinear equation corresponding to the grid of other types except the grid of the well or source-sink term in the first time step of the simulation calculation.

[0104] Further, in the step S500 of the embodiment of the method, the selection manner of the preconditioning method specifically comprises:

[0105] When the linear equations are solved by using a fully implicit method, the CPR method is selected;

[0106] When the linear equations are solved by using an adaptive implicit method, a domain decomposition method or an ILU(0) method is selected, and the Jacob method is selected when the domain decomposition method fails.

[0107] Further, in the embodiment of the present application, the distributed input of multiple simulation sub-zones and the distributed output of simulation results are based on a predetermined distributed input and output scheme;

[0108] The distributed input and output scheme comprises:

[0109] The target simulation region is divided into multiple simulation sub-zones;

[0110] Each simulation sub-zone is meshed to obtain a corresponding sub-model: the meshing methods on both sides of the boundary of two adjacent sub-models are consistent, and a corresponding virtual grid layer is arranged on both sides of the interface between the two adjacent sub-models and the corresponding entity grid, which serves as a marker layer in the subsequent sub-model merging process;

[0111] The obtained multiple sub-models are distributedly input into a numerical simulator;

[0112] The numerical simulator merges the input multiple sub-models to obtain the target simulation region: two adjacent sub-models are matched in coordinates through the pre-set virtual grid layer to realize merging, and all virtual grid layers are deleted after all sub-models are merged;

[0113] In response to an input simulation result output request, the numerical simulator outputs all simulation results by using a corresponding distributed output method.

[0114] The following describes the method for numerically simulating a fractured-vuggy reservoir in more detail:

[0115] A hybrid parallel computing scheme:

[0116] A.1 Domain decomposition technology

[0117] For a successful parallelization system, the effective partitioning method of the unstructured grid region is very important. Firstly, in order to achieve better numerical simulation performance, the parallel simulator needs to distribute the grid to different processors evenly, that is, the number of grid blocks distributed to each processor should be substantially the same; secondly, the number of connections on the region boundary should be minimized. The purpose of the first condition is to balance the computing work between different processors, and the purpose of the second condition is to minimize the time consumed by the exchange of information between different processors. In the simulator, the grid of the model is composed of a set of one-dimensional, two-dimensional or three-dimensional grid blocks, and the contact surface between different grids is represented by connection information. The overall grid system is an unstructured grid system. In parallel computing simulation, the domain decomposition technique is used to make the CPU participating in the calculation balanced in terms of computing tasks, memory requirements and communication. The software first constructs the link data in the CSR format according to the link condition of the model grid, and if the grid is input through a distributed multi-grid file, the distributed CSR format will be used to construct the link data. According to these data, the software uses the ParMETIS graph partitioning parallel library (Karypis and Schloegel, 2013) to decompose the simulation region, and decomposes the simulation region into a corresponding number of sub-regions according to the number of CPU or core (computing unit) participating in the calculation. Each computing unit will be responsible for the calculation task of a sub-region.

[0118] The coarse-grained parallel computing scheme is realized by domain decomposition, which can realize the parallelization of the assembly of Jacobian matrix, the solution of linear equations, Newton iteration, the calculation of EOS parameters and I / O. Figure 2 The domain decomposition schematic diagram of the two-dimensional simulation region of the embodiment of the application is shown, Figure 2 The basic idea of domain decomposition is given, and the figure shows that a small grid region (composed of 12 units) is decomposed into three regions. The grid blocks are distributed to different processors by the partitioning method, and all the grid blocks are rearranged by the labels in the respective assigned processors. The units corresponding to these grid blocks are explicitly stored in the processors and are defined by a set of indicators. This set of indicators is called the update group. Specifically, the update group is further divided into two subgroups: the internal group and the boundary group. The units in the internal group only involve the information of the current processor when updating. The units in the boundary group are connected to at least one unit in other processors, and the units in the boundary group need information from other processors when updating the numerical value. The units in other processors connected to the units in the boundary group are called the external group. The units in the external group need the information of the units in the boundary group connected thereto when updating the information. Figure 3 The domain decomposition schematic diagram of the three-dimensional simulation region of the embodiment of the application is shown. For referenceFigure 3 A three-dimensional simulation region is divided into 16 sub-domains.

[0119] A.2 Three-level parallel computing scheme

[0120] The three-level parallel computing scheme is adopted in the embodiment of the present application. The first-level parallel computing scheme is a task-driven parallelization with coarse granularity through domain decomposition (based on each computing node or each CPU). The second-level parallel computing scheme is parallelization of the computing cores sharing the memory. This part of parallel computing does not need domain decomposition, and the workload distribution among the computing cores is realized by the system through multi-thread parallelism. The third-level parallel computing scheme is data-driven parallelism, i.e. GPU parallel computing. The embodiment of the present application only designs the third-level parallel computing scheme, and reserves the relevant interface but does not make specific development.

[0121] In the embodiment of the present application, the three-level parallel computing scheme allows only the first-level parallel computing, taking each computing core as a coarse-grained parallel computing unit, or only the second-level parallel computing without domain decomposition. The scheme of only using the second-level parallel computing is more suitable for small-scale simulation computing using a single CPU (multi-core, OPENMP multi-thread parallelism). In the embodiment of the present application, the first-level parallel computing between distributed or shared memory CPUs is realized by MPI, and the hybrid parallel computing scheme based on distribution and shared memory is realized by combining MPI and OPENMP. The corresponding hybrid memory structure is shown in Figure 4 The hybrid memory structure shown in Figure 4 can realize the distributed parallel computing scheme, the shared memory parallel computing scheme or the hybrid parallel computing scheme. The hybrid parallel computing scheme refers to that each computing unit (including several computing cores, CPUs or GPUs, or a node of a supercomputer) will be responsible for the calculation of a sub-domain after domain decomposition. Although the calculation of each sub-domain is independent, it needs to coordinate each computing unit through MPI communication, and the whole model calculation is simulated as a whole. All the computing cores, CPUs or GPUs in each computing unit share the memory space, and they jointly complete the calculation in a sub-domain. According to experience, in the case without GPU participation, pure MPI first-level parallel computing can achieve the best effect if the communication speed between computing units is good. In actual application, the memory structure characteristics of multi-core CPUs and supercomputer nodes should be considered, and one of the three parallel computing schemes, i.e. hybrid, distributed or shared memory, can be selected to achieve the best parallel computing effect.

[0122] In the design of parallel computing scheme, an important task is to design a convenient communication scheme. The grid numbering is important for convenient communication and calling of third-party linear equation solver. After domain decomposition, the global numbering of the grids in each subdomain may not be continuous. In order to facilitate, the following numbering methods are adopted in the software. Firstly, the global numbering, that is, the order numbering of the grid when the user inputs. After domain decomposition, each grid in the subdomain has its own local numbering. This numbering can be based on the original natural order of the grid in the subdomain, or the minimum bandwidth reordering numbering. In addition, the software also reorders the global numbering to obtain a new global numbering. The new global numbering ensures that the grid numbering in each subdomain is continuous and increasing. According to the new global numbering of a certain grid, the subdomain and position of the grid can be easily determined. Therefore, the embodiments of the present application have the following three grid numbering schemes:

[0123] 1) Global numbering, mainly used for input and output control and early grid preprocessing.

[0124] 2) Local numbering, used for calculation in the subdomain.

[0125] 3) New global numbering, used for communication and calling of third-party parallel linear equation solver.

[0126] A key task of communication between computing units (or CPUs) is the exchange of parameters on the boundary. In simulation calculation, the equation corresponding to a grid is related to the calculation results of each iteration step of the adjacent grid. If the grid is located on the boundary and the adjacent grid is allocated to other computing units, communication between computing units is required to exchange parameters. The embodiments of the present application design a set of functions to realize efficient and accurate exchange of parameters between computing units.

[0127] A.3 Shared memory multi-thread parallel processing

[0128] OPENMP is used to make the software multi-threaded and parallel. Multi-threaded parallelism is usually only effective for loops, so the time-consuming processes are designed as large loops and the effective range of the parameters in the loop (public or private) is carefully considered to avoid interference between parameters. The three most time-consuming parts in the simulation process are the linear equation solving, the Jacobian matrix assembly and the state equation calculation. The linear equation solving can use third-party linear solvers that support multi-threaded parallel computing, such as AMGCL, PETSC and TRILINOS, etc. The Jacobian matrix assembly and the state equation calculation are redesigned as large loop bodies based on the number of model grids to achieve effective multi-threaded parallel computing. Tests show that the multi-threaded parallel effect of these two parts is very good, and the parallel speedup is very obvious when the grid number has a certain scale. For the multi-threaded parallel effect of linear equation solving, different linear equation solvers will have a big difference, and in general, the parallel effect of the Jacobian matrix assembly and the state equation calculation is better.

[0129] A.4 Information exchange between processors

[0130] For MPI parallel, information exchange between processors is an important part of parallel computing. After the grid is partitioned, if the grids between different processors are connected, the processors where these connected grid blocks are located need to exchange information in real time when performing parallel computing, such as checking the convergence of the Newton iteration process, which requires global exchange between processors. In addition to the exchange between processors when solving linear equations, exchange between processors of adjacent partitions is also required when updating the values of main variables. The management of information exchange between different processors is usually implemented by a function. When this function is activated by a processor, the vector unit corresponding to the external group grid block will exchange.

[0131] B Nonlinear equation and linear equation solving method:

[0132] B.1 Solution of nonlinear equations

[0133] Nonlinear equation solving is a key part of numerical simulation. For the time term in the nonlinear equation, explicit difference and implicit difference can be used. Explicit difference is forward difference in the time domain, so the value at the new time step can be obtained, that is, X(n+1) = Q[X(n)], where X(n) is the value at time T, and X(n+1) is the value at time T+1, that is, explicit difference can completely calculate the value at the next time step from the known or calculated value.

[0134] The implicit difference needs not only the known or calculated value but also the value of the new time step to obtain the value of the new time step, i.e. X(n+1)=Q[X(n), X(n+1)], which means that the implicit difference method must be realized by iteration and solving the simultaneous equations, which may encounter two problems: one is that the iteration process does not necessarily converge, and the other is that the simultaneous equations may be ill-conditioned and have no definite solution. The biggest advantage of the implicit solution method is that it has unconditional stability.

[0135] In order to improve the calculation efficiency, three solving modes of full implicit method (FIM), implicit pressure explicit saturation method (IMPES) and adaptive implicit method (AIM) are adopted in the embodiments of the present application. Different methods can be adopted in the program for different grids. The main differences among the three methods are as follows: (1) the full implicit method (FIM) is to simultaneously solve the pressure and saturation at the current time according to the parameters at the current time and the pressure and saturation at the last time by using Newton-Raphson iteration; (2) the implicit pressure explicit saturation method (IMPES) is to solve the oil pressure at the current time by using the implicit method, to obtain the water pressure at the current time according to the oil pressure at the current time and the interfacial capillary pressure at the last time, to obtain the gas pressure at the current time by using a similar method, and to obtain the oil saturation and water saturation at the current time by substituting the oil pressure and the water pressure into the respective mass conservation equations, and finally to calculate the gas saturation according to the saturation constraint condition; and (3) the adaptive implicit method (AIM) is to determine which method to use to calculate the saturation and pressure according to the change of the saturation at the last time step or other methods, to use the full implicit method to solve the saturation and pressure in the area where the saturation changes greatly, and to use the IMPES method in other areas.

[0136] The above three methods have respective characteristics and advantages. The FIM has the advantage of unconditional stability, but the efficiency is slightly poor; the IMPES needs less memory and calculation workload, and thus is more efficient, but is conditionally stable; and the AIM has the advantages of both, and has good numerical stability and high calculation efficiency.

[0137] The actual application test shows that the AIM has good calculation effect, however, the AIM needs to switch between the FIM and the IMPES according to the fluid conditions of the grid when solving a certain grid, and the appropriate switching is the key to the successful application of the AIM, and the embodiments of the present application realize two switching conditions.

[0138] The first condition is to switch according to the change of the saturation, i.e. if the change amount s of the saturation of a component in a grid is greater than a given standard in a time step, the nonlinear equation corresponding to the grid is solved by using the FIM, otherwise, the nonlinear equation corresponding to the grid is solved by using the IMPES.

[0139] The second one is based on the Courant-Friedrichs-Lewy (CFL) condition (Joachim Moortgat, 2017). The basic method is to calculate the CFL condition coefficient C according to the following equation:

[0140]

[0141] The calculation of the above formula C can also be written as:

[0142] C = (f x +f y +f z ) / (Vd) (2)

[0143] In the above formula, f x , f y and f z are the flow rates of the fluid in the x, y and z directions, d is the density of the fluid, and V is the grid volume. When C < C max , use IMPES calculation, otherwise use FIM calculation, C max is usually taken as 1.0.

[0144] In the embodiments of the present application, the grid of the well or source-sink term is always solved by FIM. In the first time step of simulation calculation, in addition to the well and source-sink grid using FIM processing, all others use IMPES processing.

[0145] B.2 Solution of linear equations

[0146] The solution of linear equations is the most time-consuming task in numerical simulation, usually accounting for more than 60% of the entire simulation time, sometimes as high as 90%. Therefore, the efficiency of the linear equation solver is the key factor determining the efficiency of the simulator. In the embodiments of the present application, the standard distributed CSR (Compressed Sparse Row) matrix storage method is adopted, which can easily call third-party linear equation solving libraries, including some of the most famous parallel solving libraries in the world. The following libraries have been implemented:

[0147] AMGCL (Demidov, 2018);

[0148] FASP (FASP Developers, 2017);

[0149] PETSC (Balay et al., 2017);

[0150] TRILINOS (Trilinos project, 2018).

[0151] AMGCL (Demidov and Shevchenko 2012; Kindratenko 2014) is an open-source C++ library for solving linear equations, which is mainly used to solve large sparse linear systems arising from the discretization of partial differential equations on unstructured grids. The library supports shared memory and distributed memory parallel computing, supports modern multi-core architectures, allows the use of modern large parallel processors through OpenMP, OpenCL or CUDA technologies, and has high independence and flexible scalability. AMGCL provides preconditioning methods such as AMG, CPR and ILU, and the selection of preconditioning iteration methods is the key to solving large problems. Using Krylov subspace methods and AMG to generate preconditioning operators is one of the most effective methods for solving large-scale computing problems.

[0152] PETSc is the acronym for Portable, Extensible Toolkit for Scientific Computation, which is a suite of data structures and routines that provide building blocks for implementing large-scale applications on parallel and serial computers. PETSc uses the MPI standard for message passing communication. PETSc includes a suite of parallel linear / nonlinear solvers and time integrators written in Fortran, C, C++, or Python. PETSc provides parallel matrix and vector set routines and many other mechanisms needed by parallel applications. PETSc library uses a hierarchical organization pattern, and users can choose the most appropriate level of abstraction according to specific problems. Object-oriented programming technology makes the library very flexible for users. PETSc is a high-level software toolset that enables many applications to run efficiently.

[0153] Trilinos is an open-source software library designed to be used as a building block for developing scientific applications. "Trilinos" is Greek for "a string of pearls", and the software library uses this to indicate that many software packages are connected together by a common infrastructure, which can be used to construct and use sparse matrices, dense matrices and vectors, iterative and direct solutions of linear systems, parallel multilevel algebraic preconditioning, solutions of nonlinear, eigenvalue and time-dependent problems, PDE-constrained optimization problems, partitioning and load balancing of distributed data structures, automatic differentiation, and discretization of partial differential equations. Trilinos is developed by the Sandia National Laboratories in the United States based on a set of core algorithms that already exist, and it uses the functions of software interfaces such as BLAS, LAPACK and MPI. Trilinos supports distributed memory parallel computing through the Message Passing Interface (MPI). In addition, some Trilinos packages gradually support shared memory parallel computing through the KOKOS package in them. The KOKOS package provides a common C++ interface on various parallel programming models, including OpenMP, POSIX threads and CUDA.

[0154] The third-party linear equation solver libraries mentioned above have advantages and disadvantages. According to the different solving problems and the characteristics of the computer platform, a suitable solver can be selected, and more libraries can be called through PETSc.

[0155] The preconditioning method is the key to effectively solve linear equations. The selection of the preconditioning method is summarized as follows:

[0156] (1) FIM preconditioning method:

[0157] The CPR preconditioning method can be selected. It is generally believed that in the fully implicit case, CPR is one of the best preconditioning methods. CPR is a two-step preconditioning method, which first solves the pressure equation from the equation system, and then corrects the solution of the entire system according to the pressure solution. For large-scale models and models solving parabolic partial differential equations, the AMG (Algebraic Multigrid) preconditioning method or the linear equation solving method also has high efficiency, but it requires that the equation matrix must be diagonally dominant. The CPR-AMG preconditioning method is also a good choice, that is, the first step of the CPR preconditioning method is solved by AMG, and the second step is solved by ILU or Domain Decomposition method.

[0158] (2) AIM preconditioning method:

[0159] Test results show that the Domain Decomposition (additive Schwarz) or ILU (0) method has good effect. For large-scale models (above ten million grids), the Domain Decomposition method may fail, at which time the JACOB preconditioning method can be selected. For the solution of linear equations, test results show that BICGSTAB or GMRES has good effect.

[0160] Massive data preprocessing and post-processing method for large-scale C model:

[0161] In large-scale models, the number of grids can reach tens of millions or even hundreds of millions. Such models not only require a large amount of computing work and memory for simulation, but also are very challenging for pre-processing and post-processing. The computing and memory requirements of the model can be distributed to each computing unit participating in the calculation through domain decomposition parallel computing. Due to the development of computer hardware, large-scale simulation is not a problem. However, for large-scale simulation, the pre-processing and post-processing of massive data are also very challenging. The embodiments of the present application realize the pre-processing and post-processing of large-scale simulation through block modeling. The basic idea of block modeling is to divide a large model into several sub-models according to certain rules. Users can perform separate pre-processing and post-processing on these sub-models. The computer automatically combines all sub-models together as a large model for simulation calculation. These sub-models are only independent in grid partitioning, and the boundary condition processing is still considered globally. The block modeling process is as follows:

[0162] According to the simulation area, design and divide into several simulation partitions. The boundary of the divided simulation partition requires division along the vertical (same X) or horizontal direction (same Y), or along a certain horizontal layer (same Z).

[0163] Respectively, grid partitioning is performed on all simulation partitions to establish sub-models. The number of grids of each sub-model should not exceed five million. On both sides of the boundary of adjacent sub-models, the grid partitioning method is required to be basically consistent, that is, the number of grids is consistent, the grid shape is similar, and the spatial position is continuous. On the partitioning boundary of the sub-model on both sides of the adjacent interface, one layer of virtual grid corresponding to the boundary grid is added, which is only a thin layer (adopting a thin layer can ensure that the coordinates of the virtual grid on both sides of the boundary are basically consistent) and is used as a marker layer for the later grid merging of the sub-model. It will not participate in simulation calculation.

[0164] Merging sub-model grids. The software uses parallel reading of sub-model grids to improve file reading efficiency. After reading, the adjacent sub-model grids are matched through the coordinates of the set virtual grid, all sub-models are merged, the virtual grid is deleted, and finally the complete large-scale model grid is assembled. The structure of the large-scale model grid is shown in Figure 5 In the grid name processing, a grid naming method with a maximum of 13 characters is adopted. The sub-area, layer, row and column information after grid merging can be saved through the grid name.

[0165] Domain decomposition. According to the number of CPUs or GPUs participating in model simulation, the complete model grid after merging is re-decomposed for parallel simulation calculation.

[0166] • Output. Simulation results can be output in different ways according to user requirements, such as a single file or multiple files, to facilitate subsequent processing.

[0167] Regarding output, the software implements multi-file output functionality. For large-scale models with massive output data volumes, the software allows users to select the output variables, output time, or output frequency. Furthermore, the multi-file output scheme can be determined in three ways: by grid name, by X, Y, and Z coordinates, or automatically divided into multiple files based on the domain decomposition results.

[0168] The following specific example illustrates the practical effect of the numerical simulation method for fractured-vuggy reservoirs according to embodiments of the present invention:

[0169] Specific applications:

[0170] S80 is a production model with 4,926,352 active grids and 12,223,952 valid links. This model is based on a real geological model, which had over 26 million grids. Grids with zero or near-zero porosity and permeability were removed. The porosity and permeability distribution in the model is based on 3D seismic results (such as...). Figure 6 (As shown) The results were coarsened using the arithmetic mean with the grid volume as the weight and the directional averaging method, respectively.

[0171] The model simulates 13 wells operating under a hybrid working system, simulating continuous oil production for 10 years. The oilfield block has a horizontal distribution area of ​​approximately 20 km². 2 The oil reservoir is located at a depth of 3600-5000m. The initial conditions of the model assume that the oil saturation of all grids is 1.0, the pressure is at gravity equilibrium, and there is no gas. Figure 7 The model presents the oil pressure distribution after 10 years of continuous oil production from 13 oil wells, utilizing the block output display function of large-scale simulation results.

[0172] For ultra-large-scale oilfield development models, to address the difficulty in displaying simulation results due to excessive grid numbers, the following approach can be adopted: Figure 7 The block display method shown divides the simulation results of different locations in the oilfield into blocks for plotting. This embodiment of the invention divides the blocks according to the X-coordinate, but it also allows automatic block output along the Y-coordinate, Z-coordinate, or grid name. Theoretically, any size grid model can be post-processed and displayed in this way.

[0173] In order to test the 10 million grid model, the embodiment of the present application uses an existing 5 million grid model to copy a set, by adding a character in front of each grid name, to ensure that the grid name is not repeated (8 characters are used as the grid name for the 5 million grid, and 9 characters are used as the grid name for the 10 million grid), and the original model and the copied model are combined to form a 10 million grid model with doubled grid number, link number and well number. Through such operation, the formed model has a total of 9,852,704 active grids and 24,447,904 effective links. The model simulates 26 wells to produce in a mixed mode, and simulates continuous oil production for 10 years. The horizontal distribution area of the oilfield block is about 40km 2 , and the reservoir burial depth range is 3600-5000m. This 10 million grid model retains the characteristics of S80 strong non-uniformity and the complexity of the actual production process, and can be used as a typical actual 10 million grid case.

[0174] The simulation calculation of the 5 million grid model is run on a supercomputer, which uses 6 Intel Xeon Broadwell 2.6GHz CPUs, 84 cores (3 nodes are involved in calculation, each node has 2 CPUs, and each CPU has 14 cores), to calculate. The software automatically divides the simulation area into 84 calculation sub-areas for calculation by 84 cores. The absolute time used is 04:45:49, and the memory usage is 19.79GB. The TRILINOS linear equation solver library is selected to solve the linear equation, the JACOB preconditioning method is selected, the BICGSTAB method is used for solving, and pure MPI parallel calculation is adopted.

[0175] The simulation calculation of the 10 million grid model uses 16 nodes, 32 CPUs and 256 cores (each node has 2 CPUs, and each CPU has 8 cores) to calculate, uses Intel Xeon Sandy Bridge 2.6GHz CPU, and the calculation time is 08:25:10. The memory usage is 108.98GB. The TRILINOS linear equation solver library is selected to solve the linear equation, the JACOB preconditioning method is selected, and pure MPI parallel calculation is adopted. Figure 8 The time length of each calculation link of the two models of different calculation scales is analyzed, which accounts for the percentage of the total time length. It can be seen that the calculation time of the model is mostly spent on solving the linear equation, followed by constructing the Jacobian matrix, which accounts for more than 84% and more than 6%, respectively. The input reading of the model accounts for only 1% or even less. By comparing the time consumption of the two models, it can be found that the proportion of the linear equation solving time of the 10 million grid model is about 4% less than that of the 5 million grid model, and the proportion of the time consumed by the Jacobian matrix construction and model iteration operations is 1-2% higher.

[0176] 2600 million grid model application case:

[0177] The in-situ geological model of the above discussed model has more than 2600 million grids, and the embodiment of the present application directly uses the grids to perform simulation calculation on S80. Since the grids with zero porosity and permeability are not removed, a very small value is assigned to these grids, and for those grids with a volume close to 0, a volume of 0.01m3is assigned. Through these treatments, the 2600 million grids of S80 can be normally simulated and calculated, but since the grids with ultra-low porosity and permeability are not removed, the Newton iteration process of the model converges worse. 3

[0178] The model test is performed on the same calculation platform as the 1000 million grid model, and the calculation uses Intel Xeon Sandy Bridge 2.6GHz CPU, and the model is simulated and tested through 64 computing cores, 128 computing cores and 256 computing cores respectively. Due to the limitation of computing resources, more computing core tests are not performed. The main input files of the model include a 7.8GB grid file (MESH), a 2.4GB initial condition file (INCON) and a 2.2GB heterogeneous spatial parameter distribution file. The effective grid number of the model is 26,826,264, the effective link number is 80,184,113, and the simulation is performed on 13 wells for production. The test shows that the simulation calculation of the model can be successfully performed using 64 computing cores, 128 computing cores and 256 computing cores. Since it is a two-component oil-water model, the total number of equations to be solved is 53,652,528.

[0179] ​The input of data and preparation of simulation calculation before the time step Newton iteration calculation need to be carried out, including domain decomposition, reading in data and distributing data, grid sorting, source and sink item and initial condition corresponding grid search, etc. For such a large scale model, it is also very challenging, and the retrieval speed is greatly improved by optimizing algorithm, such as using dichotomy to retrieve maximum, and in addition, distributed parallel computing is used to ensure the feasibility of simulation calculation of such a large scale. The time spent in this part of 64-core, 128-core and 256-core calculation is 467 seconds, 620 seconds and 915 seconds respectively. Since data distribution and parallel computing coordination need more time, the more cores used, the more time consumed in this part. Most of the time is spent on reading and publishing data, such as in the case of 128-core: reading and distributing initial condition file takes 245 seconds; reading and distributing grid data takes 290 seconds; reading and distributing non-uniform spatial parameter distribution file (svparam.dat) takes 56 seconds, which accounts for more than 95% of the total model preparation time. This part is only a one-time time-consuming, which is not significant compared to the overall simulation running time. The average time spent in simulation calculation per Newton iteration step is: 2.28 seconds for solving linear equations, 0.54 seconds for assembling Jacobian matrix and 0.68 seconds for other operations in 128-core; 1.08 seconds for solving linear equations, 0.19 seconds for assembling Jacobian matrix and 0.65 seconds for other operations in 256-core. The results show that the most time-consuming linear equation solving and Jacobian matrix assembling have very obvious parallel effect, reaching or exceeding the ideal speedup ratio.

[0180] The tests of different number of computing cores are different in total memory usage, the peak memory occupation of 64-core, 128-core and 256-core tests is 82.06GB, 140.13GB and 275.62GB respectively; in distributed memory parallel computing, part of the data needs to be repeatedly occupied in memory in each core, so the more cores used, the more total memory occupied, but the average memory occupied by each core is 1.28GB, 1.10GB and 1.09GB respectively. When the number of computing cores increases, the memory occupied by each core will be less. The average memory of a single core of a supercomputer usually ranges from several GB to several tens of GB, so from the memory requirement aspect, the technical solution of the embodiment of the application can run a model with more than one hundred million grids without any problem.

[0181] The fracture-cave type oil reservoir numerical simulation method of the embodiment realizes super large scale fine simulation of the fracture-cave type oil reservoir, improves the efficiency and simulation scale of parallel simulation by considering the characteristics of different computer hardware platforms, such as multi-core CPU and node structure of supercomputer, and adopting a hybrid three-level parallel computing scheme combining shared memory and distributed memory. The embodiment develops a good efficient software interface, successfully calls international well-known linear equation solving libraries such as PETSc, TRILINOS, FASP and AMGCL, adds efficient preconditioning methods such as CPR and AMG, ensures that there is a suitable linear equation solver for selection under many different hardware conditions, and realizes the practical application of super large scale model simulation through a distributed input output scheme.

[0182] In the embodiment 2, a fracture-cave type oil reservoir numerical simulation system is correspondingly provided based on the fracture-cave type oil reservoir numerical simulation method provided in the embodiment 1, and the fracture-cave type oil reservoir numerical simulation system comprises a processor and a memory, and the processor realizes the fracture-cave type oil reservoir numerical simulation method as described in the embodiment 1 when executing a computer program stored in the memory.

[0183] The fracture-cave type oil reservoir numerical simulation system of the embodiment has the same beneficial effects as the fracture-cave type oil reservoir numerical simulation method described in the embodiment 1, and thus will not be described here again to avoid repetition.

[0184] The above has described the embodiments of the present application, and the above description is exemplary, not exhaustive, and is not limited to the disclosed embodiments. Many modifications and changes are obvious to those skilled in the art without departing from the scope and spirit of the described embodiments.

Claims

1. A numerical simulation method for fractured-vuggy reservoirs, characterized in that, include: Integrate multiple simulation partitions of the distributed input into a target simulation region; The target simulation region is meshed to obtain the computational task; Choose one computing scheme from the predetermined distributed parallel computing scheme, shared memory parallel computing scheme, and hybrid parallel computing scheme; Based on the selected computing scheme, the computing task is divided into a corresponding number of subtasks; All subtasks are executed in parallel: In this process, the nonlinear equations are solved using a predefined fully implicit method, implicit pressure explicit saturation method, or adaptive implicit method, and a third-party linear equation solver library is called to solve the linear equations based on the selected precondition method. All simulation results obtained by executing the subtasks will be output in a distributed manner. Specifically, the adaptive implicit method is as follows: Determine whether the change in saturation of the target component in the target grid within a time step is greater than a predetermined change reference value. If yes, use the fully implicit method to solve the nonlinear equation corresponding to the target grid. If no, use the implicit pressure explicit saturation method to solve the nonlinear equation corresponding to the target grid. Alternatively, the adaptive implicit method is specifically as follows: Obtain the CFL condition coefficients; Determine whether the CFL condition coefficient is less than the predetermined maximum value of the CFL condition coefficient. If yes, use the implicit pressure explicit saturation method to solve the nonlinear equation corresponding to the target mesh. If no, use the fully implicit method to solve the nonlinear equation corresponding to the target mesh. Alternatively, when solving nonlinear equations, the fully implicit method is used to solve the nonlinear equations corresponding to the well or source / sink term grids; in the first time step of the simulation calculation, the nonlinear equations corresponding to other types of grids, except for the well or source / sink term grids, are solved using the implicit pressure explicit saturation method. Among them, the distributed input of multiple simulation partitions and the distributed output of simulation results are implemented based on a predetermined distributed input-output scheme; The distributed input / output scheme includes: The target simulation region is divided into multiple simulation partitions; Each simulation partition is meshed to obtain the corresponding sub-model: the meshing method on both sides of the boundary of two adjacent sub-models is the same, and a corresponding virtual mesh layer is set on both sides of the interface between two adjacent sub-models and the corresponding solid mesh. The virtual mesh layer is used as a marker layer in the subsequent sub-model merging process. The resulting multiple sub-models are then distributed and input into the numerical simulator. The numerical simulator merges multiple input sub-models to obtain the target simulation region: adjacent sub-models are merged by coordinate matching through a pre-set virtual mesh layer, and all virtual mesh layers are deleted after merging all sub-models. In response to the input simulation result output request, the numerical simulator outputs all simulation results using a corresponding distributed output method.

2. The numerical simulation method for fractured-vuggy reservoirs according to claim 1, characterized in that, The distributed parallel computing scheme is specifically as follows: The target simulation region after mesh partitioning is decomposed into multiple subdomains; The multiple subdomains are respectively assigned to multiple corresponding computing units for parallel computing; The interaction between the multiple computing units is achieved through MPI communication.

3. The numerical simulation method for fractured-vuggy reservoirs according to claim 1, characterized in that, The shared memory parallel computing scheme is specifically as follows: The computational task is assigned to the computational unit; The computing unit performs parallel computing on the assigned computing tasks based on multiple smallest computing sub-units therein; The plurality of smallest computational subunits share the memory space of the computational unit.

4. The numerical simulation method for fractured-vuggy reservoirs according to claim 1, characterized in that, The hybrid parallel computing scheme is specifically as follows: The target simulation region after mesh partitioning is decomposed into multiple subdomains; The multiple subdomains are respectively assigned to multiple corresponding computing units for parallel computing; The interaction between the multiple computing units is achieved based on MPI communication; The computing unit performs parallel computation on the allocated subdomain based on multiple smallest computing sub-units therein; The plurality of smallest computational subunits share the memory space of the computational unit.

5. The numerical simulation method for fractured-vuggy reservoirs according to claim 2 or 4, characterized in that, The constraints for the domain decomposition include: The difference in the number of grid cells in any two subdomains is less than a predetermined difference reference value; The total number of mesh connections on all subdomain boundaries is minimized.

6. The numerical simulation method for fractured-vuggy reservoirs according to claim 5, characterized in that, The process of meshing the target simulation region also includes: The grid within the target simulation area is first globally numbered; After performing domain decomposition on the target simulation region after mesh partitioning to obtain multiple subdomains, the method further includes: Locally number the grid within each subdomain; The numerical simulation method for fractured-vuggy reservoirs also includes: The grid within the target simulation region is assigned a second global number to facilitate interaction between the computing units and the invocation of the third-party linear equation solution library.

7. The numerical simulation method for fractured-vuggy reservoirs according to claim 1, characterized in that, The selection method for the preconditioning method is specifically as follows: When using a fully implicit method to solve linear equations, choose the CPR method; When using an adaptive implicit method to solve linear equations, choose either the domain decomposition method or the ILU(0) method. If the domain decomposition method fails, choose the Jacob method.

8. A numerical simulation system for fractured-vuggy reservoirs, characterized in that, It includes a processor and a memory, wherein the processor executes a computer program stored in the memory to implement the numerical simulation method for fractured-vuggy reservoirs as described in any one of claims 1-7.

Images