Patent Publication Number: US-10769326-B2

Title: Parallel solution for fully-coupled fully-implicit wellbore modeling in reservoir simulation

Description:
CROSS REFERENCE TO RELATED APPLICATIONS 
     This application is a continuation of U.S. patent application Ser. No. 14/948,605 filed Nov. 23, 2015 and which claims priority to U.S. Provisional Patent Application No. 62/164,083 filed May 20, 2015, which are hereby incorporated by reference. 
    
    
     BACKGROUND OF THE INVENTION 
     1. Field of the Invention 
     The present invention relates to computerized simulation of hydrocarbon reservoirs, and in particular to simulation of wellbore flow in reservoirs with complex multi-lateral wells tightly coupled with thousands of reservoir grid cells in a large high resolution reservoir simulation model. 
     2. Description of the Related Art 
     Reservoir simulation is widely used in the petroleum industry to analyze by computerized processing the performance of subsurface hydrocarbon reservoirs and to manage and optimize production from such reservoirs. One type of reservoir simulation is what is known as fully-coupled full-implicit well-reservoir simulation. Presently, most wells drilled are multilateral horizontal wells with long reach to increase reservoir contact. At the same time, finer grid reservoir simulation model are used to enhance the fidelity of recovery process analysis and better optimize and plan future reservoir management operations. As a result, it is not uncommon to have a multi-lateral well penetrating several thousand grid cells. 
     To properly model the physics of flow in and around the wellbore, the well may also be segmented to represent in detail the flow physics inside the wellbore. This in turns leads to more accurate boundary condition for the inflow performance calculation for fluids flow into and out of a well cell. 
     So far as is known, current art computerized fully-coupled fully-implicit reservoir simulators have used what are known as the Rowsum or Colsum approximations in the construction of the preconditioning step of the iterative solver of the simulator. This is because earlier methods than the Rowsum or Colsum approximations used either pre-elimination of the well equations using the reservoir equations, or direct application of a preconditioning method on the composite matrix with both the well and reservoir equations. However, these earlier methods were only suitable if the numbers of grid cells penetrated by the wells were small. 
     Otherwise, the numbers of fill terms in the solution matrices of the earlier methods became too large and impractical. The earlier solver methods also in a number of cases were too complex to implement because the well equations and reservoir equations sets had different characteristics and difficulties. At the same time, the number of algebraic equations per well segment was typically different from the number of equations per reservoir grid cell. This complicated the solver book-keeping of processor node assignment and reservoir and well cell data distribution, and also reduced code complexity, impact on solution algorithm, as well as code efficiency. 
     The Rowsum or Colsum methods which came into use to replace earlier methods to account for the well influence matrix were simple to implement and were, as noted, typically used in current reservoir simulators. Unfortunately, this method was weak, and lacked the required robustness for the complex problems where a long-reach multilateral well might in normal cases penetrate thousands of grid cells. This was particularly true for a highly heterogeneous reservoir modeled with fine grids, and when the well might be cross-flowing, meaning that some of the well perforations may have fluid inflow into the wellbore while some other perforations may have backflow from the wellbore into the reservoir. 
     SUMMARY OF THE INVENTION 
     Briefly, the present invention provides a new and improved computer implemented method of simulating the flow of the multiphase flow of fluids in computerized reservoir simulation of flow of such multiphase fluids in wellbores of multilateral wells which have flow interchange at a plurality of locations along the extent of their lengths in the reservoir with cells of a subterranean hydrocarbon reservoir organized into a grid of reservoir cells based on input reservoir data, with the reservoir cells having multiphase fluid flow taking place therein. The computer implemented method partitions the reservoir into a plurality of reservoir grid cells in a grid formed of a number of contiguous grid cells; and partitions the multilateral wellbores into a plurality of wellbore cells in a grid formed of a number of contiguous wellbore cells along the lengths of the wellbores. Representations of state changes and pressures within the reservoir grid cells and flow interchange with the wellbore cells at locations of flow interchange with the wellbore cells are established. Representations of the flow interchange with the wellbore cells and flow interchange with the reservoir grid cells at locations of flow interchange with the reservoir grid cells are established for the wellbore cells. A series preconditioner is formed as a series based expansion by matrix-vector multiplication operations of well influence terms based on the established representations of flow interchange with the wellbore cells for the reservoir grid cells and the established representations of flow interchange with the reservoir grid cells for wellbore cells having flow interchange. A conjugate residual interactive matrix solution is applied in the computer to solve the representations of flow of the multiphase fluids in the reservoir grid cells and wellbore cells and obtain a residual. The steps of forming a series preconditioner by matrix-vector multiplication operations and applying a conjugate residual interactive matrix solution in the computer to solve for the representations of flow of the multiphase fluids of the wellbore cells are repeated until the obtained residual is within an established limit of accuracy. The computerized simulation of representations of flow of the multiphase fluids of the wellbore cells when the obtained residual is within the established limit of accuracy are stored, and an output display is formed of the stored computerized simulation of representations of flow of the multiphase fluids of the wellbore in the reservoir. 
     The present invention also provides a new and improved data processing system for computerized simulation of flow of multiphase fluids in wellbores of multilateral wells which have flow interchange at a plurality of locations along the extent of their lengths in the reservoir with cells of a subterranean hydrocarbon reservoir organized into a grid of reservoir cells based on input reservoir data, the reservoir cells having multiphase fluid flow taking place therein. The data processing system includes a processor which partitions the reservoir into a plurality of reservoir grid cells in a grid formed of a number of contiguous grid cells, and partitions the multilateral wellbores into a plurality of wellbore cells in a grid formed of a number of contiguous wellbore cells along the lengths of the wellbores. The processor also establishes for the reservoir grid cells representations of state changes and pressures within the reservoir grid cells and flow interchange with the wellbore cells at locations of flow interchange with the wellbore cells. The processor also establishes for the wellbore cells representations of the flow interchange with the wellbore cells and flow interchange with the reservoir grid cells at locations of flow interchange with the reservoir grid cells. The processor then forms a series preconditioner as a series based expansion by matrix-vector multiplication operations of well influence terms based on the established representations of flow interchange with the wellbore cells for the reservoir grid cells and the established representations of flow interchange with the reservoir grid cells for wellbore cells having flow interchange. The processor applies a conjugate residual interactive matrix solution in the computer to solve the representations of flow of the multiphase fluids in the reservoir grid cells and wellbore cells to obtain a residual. The processor then repeats the steps of forming a series preconditioner by matrix-vector multiplication operations and applying a conjugate residual interactive matrix solution in the computer to solve for the representations of flow of the multiphase fluids of the wellbore cells until the obtained residual is within an established limit of accuracy. The data processing system also includes memory storing the computerized simulation of representations of flow of the multiphase fluids of the wellbore cells when the obtained residual is within the established limit of accuracy, and a display which displays the stored computerized simulation of representations of flow of the multiphase fluids of the wellbore cells of the reservoir. 
     The present invention also provides a new and improved data storage device which has stored in a non-transitory computer readable medium computer operable instructions for causing a data processor to simulate flow of multiphase fluids in wellbores of multilateral wells which have flow interchange at a plurality of locations along the extent of their lengths in the reservoir with cells of a subterranean hydrocarbon reservoir organized into a grid of reservoir cells based on input reservoir data. The stored instructions cause the processor to partition the reservoir into a plurality of reservoir grid cells in a grid formed of a number of contiguous grid cells; and to partition the multilateral wellbores into a plurality of wellbore cells in a grid formed of a number of contiguous wellbore cells along the lengths of the wellbores. The stored instructions cause the processor to establish for the reservoir grid cells representations of state changes and pressures within the reservoir grid cells and flow interchange with the wellbore cells at locations of flow interchange with the wellbore cells, and to establish for the wellbore cells representations of the flow interchange with the wellbore cells and flow interchange with the reservoir grid cells at locations of flow interchange with the reservoir grid cells. The stored instructions also cause the processor to form a series preconditioner as a series based expansion by matrix-vector multiplication operations of well influence terms based on the established representations of flow interchange with the wellbore cells for the reservoir grid cells and the established representations of flow interchange with the reservoir grid cells for wellbore cells having flow interchange. The stored instructions cause the processor to apply a conjugate residual interactive matrix solution in the computer to solve the representations of flow of the multiphase fluids in the reservoir grid cells and wellbore cells to obtain a residual. The stored instructions cause the processor to repeat the steps of forming a series preconditioner by matrix-vector multiplication operations and applying a conjugate residual interactive matrix solution in the computer to solve for the representations of flow of the multiphase fluids of the wellbore cells until the obtained residual is within an established limit of accuracy. The stored instructions the cause the processor to store the computerized simulation of representations of flow of the multiphase fluids of the wellbore cells when the obtained residual is within the established limit of accuracy, and to form an output display of the stored computerized simulation of representations of flow of the multiphase fluids of the wellbore in the reservoir. 
     The present invention also provides a new and improved computer implemented method of simulation of flow of multiphase fluid in wellbores of multilateral wells in a subsurface hydrocarbon reservoir, the multilateral wells having interchange of fluids at a plurality of locations along the extent of their lengths in the reservoir with cells of a subterranean reservoir organized into a grid of reservoir cells based on input reservoir data, the reservoir cells having multiphase fluid flow taking place therein. The computer implemented method organizes in the computer a fully coupled nonlinear set of implicit reservoir equations with reservoir data, pressure equations and flow equations for reservoir cells having flow interchange with certain ones of the multilateral wells into a reservoir computation matrix, a vector of reservoir and fluid flow unknowns and a vector of reservoir residuals; and also organizes in the computer a fully coupled nonlinear set of implicit well equations with well data and flow equations for wellbore cells having flow interchange with certain ones of the reservoir cells into a wellbore computation matrix, a vector of fluid flow unknowns and a vector of wellbore residuals. The computer implemented method also organizes in the computer a well influence matrix based on the flow interchange of the wellbore cells and the reservoir cells, as well as a full system computation matrix comprising the reservoir computation matrix and the wellbore computation matrix, a vector of full system unknowns and a vector of full system residuals. The computer implemented method extracts pressure coefficients of the reservoir computation matrix and the wellbore computation matrix, and extracts pressure residuals from the full system residuals. The computer implemented then solves an approximate pressure solution for pressures within the reservoir and wellbore cells of the full system computation matrix by minimizing the extracted pressure residuals, and updates fluid pressures and the residuals for the reservoir cells of the full system computation matrix based on the approximate pressure solution. The computer implemented method then calculates an approximate full system update for the full system computation matrix, the well influence matrix and the updated pressures and residuals. The computer implemented method then combines the approximate full system update with the updated fluid pressures, updates the full system residuals and determines the multiphase fluid flow by solving the full system computation matrix using the fully coupled nonlinear set of conservation equations and the updated system residuals. 
     The present invention also provides a new and improved data processing system for simulation of simulation of flow of multiphase fluid in wellbores of multilateral wells in a subsurface hydrocarbon reservoir, the multilateral wells having interchange of fluids at a plurality of locations along the extent of their lengths in the reservoir with cells of a subterranean reservoir organized into a grid of reservoir cells based on input reservoir data, the reservoir cells having multiphase fluid flow taking place therein. The data processing system includes a processor which organizes in the computer a fully coupled nonlinear set of implicit reservoir equations with reservoir data, pressure equations and flow equations for reservoir cells having flow interchange with certain ones of the multilateral wells into a reservoir computation matrix, a vector of reservoir and fluid flow unknowns and a vector of reservoir residuals. The processor also organizes in the computer a fully coupled nonlinear set of implicit well equations with well data and flow equations for wellbore cells having flow interchange with certain ones of the reservoir cells into a wellbore computation matrix, a vector of fluid flow unknowns and a vector of wellbore residuals. The processor then organizes in the computer a well influence matrix based on the flow interchange of the wellbore cells and the reservoir cells, and a full system computation matrix formed of the reservoir computation matrix and the wellbore computation matrix, a vector of full system unknowns and a vector of full system residual. The processor extracts pressure coefficients of the reservoir computation matrix and the wellbore computation matrix, and also extracts pressure residuals from the full system residuals. The processor then solves an approximate pressure solution for pressures within the reservoir and wellbore cells of the full system computation matrix by minimizing the extracted pressure residuals. The processor then updates the fluid pressures and the residuals for the reservoir cells of the full system computation matrix based on the approximate pressure solution. The processor calculates an approximate full system update for the full system computation matrix, the well influence matrix and the updated pressures and residuals, and combines the approximate full system update with the updated fluid pressures, and updates the full system residuals. The processor then determines the multiphase fluid flow by solving the full system computation matrix using the fully coupled nonlinear set of conservation equations and the updated system residuals. 
     The present invention further provides a new and improved data processing system for simulation of simulation of flow of multiphase fluid in wellbores of multilateral wells in a subsurface hydrocarbon reservoir, the multilateral wells having interchange of fluids at a plurality of locations along the extent of their lengths in the reservoir with cells of a subterranean reservoir organized into a grid of reservoir cells based on input reservoir data, the reservoir cells having multiphase fluid flow taking place therein. The data processing system includes a processor which organizes in the computer a fully coupled nonlinear set of implicit reservoir equations with reservoir data, pressure equations and flow equations for reservoir cells having flow interchange with certain ones of the multilateral wells into a reservoir computation matrix, a vector of reservoir and fluid flow unknowns and a vector of reservoir residuals. The processor also organizes in the computer a fully coupled nonlinear set of implicit well equations with well data and flow equations for wellbore cells having flow interchange with certain ones of the reservoir cells into a wellbore computation matrix, a vector of fluid flow unknowns and a vector of wellbore residuals. The processor also organizes in the computer a well influence matrix based on the flow interchange of the wellbore cells and the reservoir cells, and organizes in the computer a full system computation matrix composed of the reservoir computation matrix and the wellbore computation matrix, a vector of full system unknowns and a vector of full system residuals. The processor extracts pressure coefficients of the reservoir computation matrix and the wellbore computation matrix, and then extracts pressure residuals from the full system residuals. The processor then solves for an approximate pressure solution for pressures within the reservoir and wellbore cells of the full system computation matrix by minimizing the extracted pressure residuals. The processor then updates the fluid pressures and the residuals for the reservoir cells of the full system computation matrix based on the approximate pressure solution. The processor then calculates an approximate full system update for the full system computation matrix, the well influence matrix and the updated pressures and residuals. The processor combines the approximate full system update with the updated fluid pressures, updates the full system residuals and determines the multiphase fluid flow by solving the full system computation matrix using the fully coupled nonlinear set of conservation equations and the updated system residuals. 
     The present invention further provides a new and improved data storage device having stored in a non-transitory computer readable medium computer operable instructions for causing a processor to simulate flow of multiphase fluid in wellbores of multilateral wells in a subsurface hydrocarbon reservoir, the multilateral wells having interchange of fluids at a plurality of locations along the extent of their lengths in the reservoir with cells of a subterranean reservoir organized into a grid of reservoir cells based on input reservoir data, the reservoir cells having multiphase fluid flow taking place therein. The instructions stored in the data storage device causing the processor to organize in the computer a fully coupled nonlinear set of implicit reservoir equations with reservoir data, pressure equations and flow equations for reservoir cells having flow interchange with certain ones of the multilateral wells into a reservoir computation matrix, a vector of reservoir and fluid flow unknowns and a vector of reservoir residuals. The stored instructions also cause the processor to organize in the computer a fully coupled nonlinear set of implicit well equations with well data and flow equations for wellbore cells having flow interchange with certain ones of the reservoir cells into a wellbore computation matrix, a vector of fluid flow unknowns and a vector of wellbore residuals. The stored instructions further cause the processor to organize in the computer a well influence matrix based on the flow interchange of the wellbore cells and the reservoir cells, and then organize in the computer a full system computation matrix comprising the reservoir computation matrix and the wellbore computation matrix, a vector of full system unknowns and a vector of full system residuals. The instructions cause the processor to extract pressure coefficients of the reservoir computation matrix and the wellbore computation matrix, and extract pressure residuals from the full system residuals. The instructions then cause the processor to solve an approximate pressure solution for pressures within the reservoir and wellbore cells of the full system computation matrix by minimizing the extracted pressure residuals. The instructions cause the processor to then update fluid pressures and the residuals for the reservoir cells of the full system computation matrix based on the approximate pressure solution, and calculate an approximate full system update for the full system computation matrix, the well influence matrix and the updated pressures and residuals. The instructions cause the processor to combine the approximate full system update with the updated fluid pressures, and updating the full system residuals. The instructions then cause the processor to determine the multiphase fluid flow by solving the full system computation matrix using the fully coupled nonlinear set of conservation equations and the updated system residuals. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  is an isometric view of a computerized model of a subsurface reservoir structured grid. 
         FIG. 2  is an image of a three dimensional view of a several multi-lateral wells and the penetrated grid blocks in a reservoir simulation model. 
         FIG. 3  is a plan view of several multi-lateral wellbores traversing finite-volume grid blocks in a reservoir model. 
         FIGS. 4A and 4B  are schematic diagrams of an example of a fully-coupled fully-implicit well-reservoir Jacobian matrix for an example small model. 
         FIGS. 5A and 5B  are schematic diagrams of an example of a reduced system matrix with additional non-zero dense block derivatives induced by fully-implicit well coupling terms. 
         FIG. 6  is a functional block diagram or flow chart of data processing steps for fully-coupled fully-implicit wellbore modeling in reservoir simulation according to certain embodiments of the present invention. 
         FIG. 7  is a functional block diagram or flow chart of data processing steps for fully-coupled fully-implicit wellbore modeling in reservoir simulation according to other embodiments of the present invention. 
         FIG. 8  is a schematic diagram of small model example of a reservoir subdivided into subdomains to be processed by computer processing cores according to the present invention. 
         FIG. 9  is a schematic diagram of partitioned matrix and vector data according to the present invention for the model of  FIG. 8 . 
         FIG. 10  is an isometric view of a computerized model of a subsurface reservoir structured grid reservoir simulation model for a test example of processing according to the present invention. 
         FIGS. 11, 12, 13 and 14  are comparison data plots of processing results according to the present invention against processing results according to the prior art for the model of  FIG. 8 . 
         FIG. 15  is a schematic diagram of a computer network for fully-coupled fully-implicit wellbore modeling according to the present invention. 
         FIG. 16  is a schematic diagram of an application server or computer node of the computer network of  FIG. 15 . 
     
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
     Reservoir Simulator 
     As an introduction, a more detailed description of known computerized reservoir simulation of the types identified above is provided. For the recovery of oil and gas from subterranean reservoirs, wellbores are drilled into these formations for the recovery of hydrocarbon fluid. During the recovery process, fluids such as water and/or gas are injected into the injector wells and the fluid mixture in the pore space is produced from the producer wells. In order to predict the future performance of these reservoirs and to evaluate alternative development plan, reservoir simulators are used to run simulation models. These models are first calibrated with a history matching step using existing production data. The calibrated models are then used to evaluate future operation scenarios. For example, the history-matched models may be used to determine when and where to drill additional wells in order to recover more of the remaining hydrocarbon in place. 
     The reservoir simulator is a computer-implemented software methodology which solves a system of discrete balance equations for each grid block. The discrete equations are typically formed from a finite-volume discretization of the governing system of non-nonlinear partial differential equations describing the mass, momentum, and energy conservation equations within the reservoir. 
       FIG. 1  illustrates a typical reservoir simulation domain D discretized into millions of finite volumes. Wells are drilled into the reservoir formation to inject or produce fluids into or from its pore space.  FIG. 2  illustrates a 3D-view of several multi-branch complex wells  20  penetrating many grid blocks within a reservoir model M and  FIG. 3  illustrates a 2D plan view of a reservoir model M- 1  with many multilateral wells  22 . 
     The well equations which represent the imposed constraints, either the injecting or producing fluid phase rates, or the well bottom-hole-pressure, and sometimes also include the mass, momentum, and energy balance within the wellbore itself of wells such as shown at  20  and  22 . The system, including the wells and the reservoir, represents a tightly coupled highly nonlinear system where the derivatives can be discontinuous locally making it difficult to solve. For robustness and numerical stability, implicit methods are generally applied to solve the coupled system of equations. Since the well flow terms represent the dominant boundary conditions, a robust implicit scheme for the strongly coupled system between wells and reservoir is usually highly desirable. Thus, what is known as the fully-coupled fully-implicit method is the predominant current methodology for the treatment of wells in reservoir simulation. However, due to its complexity, the coupled well-reservoir solution method is handled in some approximate way during the solution process. The present invention provides methodology to address these shortcomings. 
     The background description of Saudi-Aramco&#39;s GigaPOWERS reservoir simulator can be found in Dogru et al. (SPE119272, “A Next-Generation Parallel Reservoir Simulator for Giant Reservoirs,” Proceedings of the SPE Reservoir Simulation Symposium, The Woodlands, Tex., USA, 2-4 Feb. 2009, 29 pp.) The transient solution of the multiphase multicomponent system involves the evolution of mass and energy conservation in a sequence of time steps from the initial condition of the reservoir. For each time step, the system of nonlinear discrete equations for each finite volume is linearized using what is known as the generalized Newton&#39;s method. 
     Nomenclature 
     The symbols in the equations to be set forth below have the following meanings: 
     p=pressure 
     q=production rate 
     x i =Mole fraction 
     V j =Phase Volume 
     S j =Phase Saturation 
     c i =Overall Concentration of species i 
     ϕ=porosity 
     ρ=density 
     μ=viscosity 
     ω=mass fraction 
     R=Homogeneous reaction rate 
     D=Dispersion Coefficient 
     u=velocity 
     V i =Rock pore volume 
     n i   t =Overall number of mole
         Superscripts:       

     ref=reference 
     p=a fluid phase 
     t=total
         Subscripts:       

     i=component index 
     j=phase index 
     Species Conservation 
     A general species conservation equation for a species component i is given by: 
     
       
         
           
             
               
                 
                   
                     
                       
                         ∂ 
                         
                           c 
                           i 
                         
                       
                       
                         ∂ 
                         t 
                       
                     
                     + 
                     
                       ∇ 
                       
                         · 
                         
                           
                             ∑ 
                             
                               j 
                               = 
                               1 
                             
                             
                               n 
                               p 
                             
                           
                           ⁢ 
                           
                             ( 
                             
                               
                                 
                                   ρ 
                                   j 
                                 
                                 ⁢ 
                                 
                                   ω 
                                   ij 
                                 
                                 ⁢ 
                                 
                                   u 
                                   j 
                                 
                               
                               - 
                               
                                 
                                   ϕρ 
                                   j 
                                 
                                 ⁢ 
                                 
                                   S 
                                   j 
                                 
                                 ⁢ 
                                 
                                   
                                     
                                       D 
                                       -&gt; 
                                     
                                     ij 
                                   
                                   · 
                                   
                                     ∇ 
                                     
                                       ω 
                                       ij 
                                     
                                   
                                 
                               
                             
                             ) 
                           
                         
                       
                     
                   
                   = 
                   
                     
                       ϕ 
                       ⁢ 
                       
                         
                           ∑ 
                           
                             j 
                             = 
                             1 
                           
                           
                             n 
                             p 
                           
                         
                         ⁢ 
                         
                           ( 
                           
                             
                               S 
                               j 
                             
                             ⁢ 
                             
                               R 
                               ij 
                             
                           
                           ) 
                         
                       
                     
                     + 
                     
                       
                         ( 
                         
                           1 
                           - 
                           ϕ 
                         
                         ) 
                       
                       ⁢ 
                       
                         R 
                         is 
                       
                     
                     + 
                     
                       
                         q 
                         . 
                       
                       i 
                     
                   
                 
               
               
                 
                   ( 
                   1 
                   ) 
                 
               
             
             
               
                 
                   
                       
                   
                   ⁢ 
                   Where 
                 
               
               
                 
                     
                 
               
             
             
               
                 
                   
                       
                   
                   ⁢ 
                   
                     
                       c 
                       i 
                     
                     = 
                     
                       
                         ϕ 
                         ⁢ 
                         
                           
                             ∑ 
                             
                               j 
                               = 
                               1 
                             
                             
                               n 
                               p 
                             
                           
                           ⁢ 
                           
                             
                               ρ 
                               j 
                             
                             ⁢ 
                             
                               ω 
                               ij 
                             
                             ⁢ 
                             
                               S 
                               j 
                             
                           
                         
                       
                       + 
                       
                         
                           ( 
                           
                             1 
                             - 
                             ϕ 
                           
                           ) 
                         
                         ⁢ 
                         
                           ρ 
                           s 
                         
                         ⁢ 
                         
                           ω 
                           is 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   2 
                   ) 
                 
               
             
           
         
       
     
     If dispersion, chemical reaction and absorption are ignored, the species equation simplifies to 
     
       
         
           
             
               
                 
                   
                     
                       
                         1 
                         
                           V 
                           b 
                         
                       
                       ⁢ 
                       
                         
                           ∂ 
                           
                             n 
                             i 
                             t 
                           
                         
                         
                           ∂ 
                           t 
                         
                       
                     
                     + 
                     
                       ∇ 
                       
                         · 
                         
                           
                             ∑ 
                             
                               j 
                               = 
                               1 
                             
                             
                               n 
                               p 
                             
                           
                           ⁢ 
                           
                             ( 
                             
                               
                                 ρ 
                                 j 
                               
                               ⁢ 
                               
                                 x 
                                 ij 
                               
                               ⁢ 
                               
                                 u 
                                 j 
                               
                             
                             ) 
                           
                         
                       
                     
                   
                   = 
                   
                     
                       q 
                       . 
                     
                     i 
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
     Since the pore space of porous medium must be filled with fluids present, the pore volume must be equal to the total fluid volume. This can be expressed as: 
     
       
         
           
             
               
                 
                   
                     
                       ∑ 
                       
                         i 
                         = 
                         1 
                       
                       
                         n 
                         p 
                       
                     
                     ⁢ 
                     
                       V 
                       j 
                     
                   
                   = 
                   
                     V 
                     ϕ 
                   
                 
               
               
                 
                   ( 
                   4 
                   ) 
                 
               
             
           
         
       
     
     Where the pore volume, V ϕ , is a function of pressure alone and described as:
 
 V   ϕ   =V   ϕ   ref   e   C     r     (P     ref     -P)   (5)
 
     Pressure and the overall number of moles are the primary variables. For closure, the other equations used are: 
     
       
         
           
             
               
                 
                   
                     
                       ∑ 
                       
                         i 
                         = 
                         1 
                       
                       
                         n 
                         c 
                       
                     
                     ⁢ 
                     
                       x 
                       ij 
                     
                   
                   = 
                   1 
                 
               
               
                 
                   ( 
                   6 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       
                         ∑ 
                         
                           j 
                           = 
                           1 
                         
                         
                           n 
                           p 
                         
                       
                       ⁢ 
                       
                         
                           n 
                           j 
                           p 
                         
                         ⁢ 
                         
                           x 
                           ij 
                         
                       
                     
                     = 
                     
                       n 
                       i 
                       t 
                     
                   
                   ⁢ 
                   
                       
                   
                 
               
               
                 
                   ( 
                   7 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       ∑ 
                       
                         i 
                         = 
                         1 
                       
                       
                         n 
                         p 
                       
                     
                     ⁢ 
                     
                       S 
                       j 
                     
                   
                   = 
                   1 
                 
               
               
                 
                   ( 
                   8 
                   ) 
                 
               
             
             
               
                 
                   
                     S 
                     j 
                   
                   = 
                   
                     
                       V 
                       j 
                     
                     
                       
                         ∑ 
                         
                           j 
                           = 
                           1 
                         
                         
                           n 
                           p 
                         
                       
                       ⁢ 
                       
                         V 
                         j 
                       
                     
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
             
               
                 
                   
                     V 
                     j 
                   
                   = 
                   
                     
                       n 
                       j 
                       p 
                     
                     
                       ρ 
                       j 
                     
                   
                 
               
               
                 
                   ( 
                   10 
                   ) 
                 
               
             
             
               
                 
                   
                     u 
                     j 
                   
                   = 
                   
                     
                       - 
                       K 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         λ 
                         j 
                       
                       ⁡ 
                       
                         ( 
                         
                           
                             ∇ 
                             
                               P 
                               j 
                             
                           
                           - 
                           
                             
                               γ 
                               j 
                             
                             ⁢ 
                             
                               ∇ 
                               D 
                             
                           
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   11 
                   ) 
                 
               
             
             
               
                 
                   K 
                   = 
                   
                     [ 
                     
                       
                         
                           
                             k 
                             xx 
                           
                         
                         
                           
                             k 
                             xy 
                           
                         
                         
                           
                             k 
                             xz 
                           
                         
                       
                       
                         
                           
                             k 
                             yx 
                           
                         
                         
                           
                             k 
                             yy 
                           
                         
                         
                           
                             k 
                             yz 
                           
                         
                       
                       
                         
                           
                             k 
                             zx 
                           
                         
                         
                           
                             k 
                             zy 
                           
                         
                         
                           
                             k 
                             zz 
                           
                         
                       
                     
                     ] 
                   
                 
               
               
                 
                   ( 
                   12 
                   ) 
                 
               
             
           
         
       
     
     A typical well rate relationship for a hydrocarbon component i in moles/day used in a reservoir simulator has the form: 
     
       
         
           
             
               
                 
                   
                     q 
                     i 
                   
                   = 
                   
                     
                       ∑ 
                       
                         k 
                         = 
                         1 
                       
                       
                         n 
                         / 
                         y 
                       
                     
                     ⁢ 
                     
                       
                         
                           WI 
                           k 
                         
                         ⁡ 
                         
                           ( 
                           
                             
                               
                                 λ 
                                 o 
                               
                               ⁢ 
                               
                                 ρ 
                                 o 
                               
                               ⁢ 
                               
                                 x 
                                 i 
                               
                             
                             + 
                             
                               
                                 λ 
                                 g 
                               
                               ⁢ 
                               
                                 ρ 
                                 g 
                               
                               ⁢ 
                               
                                 y 
                                 i 
                               
                             
                           
                           ) 
                         
                       
                       ⁢ 
                       
                         ( 
                         
                           
                             p 
                             k 
                           
                           - 
                           
                             P 
                             bh 
                           
                           - 
                           
                             
                               γ 
                               well 
                             
                             ⁢ 
                             Δ 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               z 
                               k 
                             
                           
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   13 
                   ) 
                 
               
             
           
         
       
     
     Layer Well Index 
     WI is commonly referred to as the layer well index. It relates the well-cell finite-difference grid pressure to the well bottom-hole-pressure through the steady state pressure profile near the well and the determination of an equivalent radius r o  for the well cell. What are known as Peaceman formulae are frequently used for determining layer well index in reservoir simulators and this has been later extended to anisotropic rectangular grids and other more complex situations. 
     In Equation (13), λ represents the fluid phase mobility and γ well  is the wellbore gravity gradient which depends on the fluid mixture density in the wellbore. P bh  is the bottom-hole-pressure at well reference depth and p k  is the well cell pressure. A generalized form of reservoir simulation with the Peaceman well index is contained as Appendix B in Fung et al.: “Unconstrained Voronoi Grids for Densely Spaced Complex Wells in Full-Field Reservoir Simulation,” SPE Journal, October 2014, pp 803-815. 
     The wellbore flow and constraints can be represented at varying level of details. In the simplest case, the entire well is considered as a single storage at steady-state condition and the fluid content is exactly the mixture of the total inflow from all the contributing layers. In this case, the fully-implicit well equation is simply Equation (13) where all variables are taken at the new time level. The next step is to include a mass balance equation of each fluid component for the well, but the entire well is still considered as a single storage. In this case, wellbore transients can be accounted for and the well equation set for the entire wellbore network will have nwvar=(nvar+1) equations where nvar is the number of fluid components. However, a more detailed well model can have the wellbores divided into several segments where the well-segment fluid component and total mass balances, plus the well constraint equations are solved simultaneously with the reservoir system of equations. For thermal simulation, an energy balance equation for each well segment is also needed. 
     Fully-Coupled Systems 
     The fully-coupled system with implicit well equations and implicit reservoir equations is as follows: 
                       [           A   WW           A   WR               A   RW           A   RR           ]     ⁢     {           X   W               X   R           }       =     {           R   W               R   R           }             (   14   )               
Where the Jacobian derivatives of the nonlinear well and reservoir equations are:
 
                 A   WW     =       ∂     F   w         ∂     X   w           ;       A   WR     =       ∂     F   w         ∂     X   R           ;       A   RW     =       ∂     F   R         ∂     X   W           ;       A   RR     =       ∂     F   R         ∂     X   R                 
with X W  and X R  being the well and reservoir variables, R W  and R R  being the initial well and reservoir residuals. Since each grid cell can have multiple equations, and the number of equations for the well grids and the reservoir grids can be different, nvar is used as the number of reservoir equations per grid and nwvar as the number of well equations per grid. It is to be noted that nvar in general is not equal to nwvar. Therefore, each non-zero in A WW  is a dense block of (nwvar*nwvar); each non-zero in A RR  is a dense block of (nvar*nvar); each non-zero in A WR  is a dense block of (nwvar*nvar), and in A RW  is a dense block of (nvar*nwvar).
 
     The solution method for an implicit reservoir system is typically a preconditioned iterative method which is suitable for asymmetric ill-conditioned large sparse matrices. For a serial to limited parallel system, the nested factorization, NF, method or the defined-level incomplete lower-upper, ILU(k), triangular factorization with a domain decomposition method, such as the additive Schwarz method, are commonly used. A pressure predictor-corrector method known as the constraint pressure residual, CPR, method as a two-stage preconditioner is also well known in current art. The Krylov subspace methods for an asymmetric system matrix frequently used in the current art are the ORTHOMIN algorithm and the GMRES algorithm. Alternatively, the BICGSTAB algorithm is also used but is less popular due to the higher work counts per iteration. 
     As pointed out above, two earlier methods prior to the Rowsum or Colsum approximations were inefficient and complicated to apply for complex well and large coupled system of equations. These were solution of the fully coupled system directly with a preconditioned iterative method and solving the reduced system by preconditioning the full reduced system matrix. They can be useful for serial computation of a small coupled well-reservoir simulation system, but impractical for large scale parallel application of a robust simulator. 
     Solving the Fully Coupled System Directly Using a Preconditioned Iterative Method 
     In this method, the composite matrix was solved directly. That is, the matrix equations were set: 
     
       
         
           
             
               
                 
                   
                     
                       A 
                       · 
                       x 
                     
                     = 
                     R 
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     Where 
                     ⁢ 
                     
                       : 
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       A 
                       = 
                       
                         [ 
                         
                           
                             
                               
                                 A 
                                 WW 
                               
                             
                             
                               
                                 A 
                                 WR 
                               
                             
                           
                           
                             
                               
                                 A 
                                 RW 
                               
                             
                             
                               
                                 A 
                                 RR 
                               
                             
                           
                         
                         ] 
                       
                     
                     ; 
                     
                       x 
                       = 
                       
                         { 
                         
                           
                             
                               
                                 X 
                                 W 
                               
                             
                           
                           
                             
                               
                                 X 
                                 R 
                               
                             
                           
                         
                         } 
                       
                     
                     ; 
                     
                       R 
                       = 
                       
                         { 
                         
                           
                             
                               
                                 R 
                                 W 
                               
                             
                           
                           
                             
                               
                                 R 
                                 R 
                               
                             
                           
                         
                         } 
                       
                     
                   
                 
               
               
                 
                   ( 
                   15 
                   ) 
                 
               
             
           
         
       
     
     This approach solves the complete system together and obtains the solution vector [X W , X R ] simultaneously.  FIG. 4B  is small example model of the fully-coupled fully-implicit well-reservoir Jacobian matrix. An example of a small reservoir model M- 2  with fifteen grid cells  40  and a wellbore  42  penetrating five grid cells at the center is shown in  FIG. 4A . A matrix  44  in  FIG. 4B  depicts schematically the connectivity of the cells  40  at the reservoir (A RR ) and  42  of the wellbore (A RW ). The letters w, x, wr, and rw in the matrix of  FIG. 4  represent the dense sub-matrices non-zeros for the well segment derivatives, reservoir grid derivatives, well-reservoir coupling derivatives, and reservoir-well coupling derivatives. This method may be used for small problems with a few wells in a small to moderate size model. It is a method which was popular in serial computation many years ago. 
     This method is more amenable to implementation in algebraic preconditioning method such as the level-of-fill incomplete lower upper ILU(k) triangular factorization method. It has been found not to be suitable for a structured-grid based preconditioning method such as the nested factorization (NF) method. The bookkeeping of cell assignments and matrix organization is more complex due to the irregular dense blocks to be factored as well as the negative implication on parallel domain decomposition methods. This was because wells can traverse multiple subdomains requiring more network communication to factorize and solve. 
     Solving the Reduced System by Preconditioning the Full Reduced System Matrix 
     This solution method first decouples the well equation from the reservoir equation to form a reduced system:
 
 Ā·X   R =[ A   RR   −A   RW   A   WW   −1   A   WR ]· X   R =[ R   R   −A   RW   A   WW   −1   R   W ]  (16)
 
     The reduced system is solved to find the reservoir variables [X R ]. The solution for the well equations is obtained by the following back solve step after the reduced system has been solved to obtain the reservoir variables:
 
 X   W   =A   WW   −1 ( R   W   −A   WR   X   R )  (17)
 
 FIG. 5B  shows the same example of reservoir model M- 2  of  FIGS. 4A and 5A  with a reduced system matrix  50  according to Equation (16).  FIG. 5B  is a small example of the reduced system matrix with the additional non-zero dense block derivatives induced by the fully-implicit well coupling terms. When the wells in the model penetrate very few grid cells, this method could be acceptable. If the wells in reservoir, however, penetrate many grid cells, the number of additional non-zeros in the reduced system matrix  50  became excessive which render this method unattractive. This was because the reduced system matrix [A RR −A RW A WW   −1 A wr ] had many well-induced non-zero terms, making it far more expensive to solve as compared to the reservoir matrix without the well terms.
 
     Solving the Reduced System Using Rowsum or Colsum Preconditioning 
     The primary method of the current start-of-the-art for fully-implicit coupled system solution used in present simulators is either the Rowsum or Colsum preconditioning method. While the method is simple and suitable for wells with a moderate number of well perforated grid cells, it may start to break down with slow convergence, or fail to converge for the overall solver, when the simulation model contains many complex multi-lateral long-reach wells which may perforate thousands of grid cells with complex heterogeneity and with very high perforation multiphase flow rates. In these complex situations, the coupled solver can stall or fail to converge all together due to excessive well-influence-coefficient couplings which are only weakly preconditioned relative to the inter-cell flow terms. 
     Because the reduced system matrix [A RR −A RW A WW   −1 A WR ] can have significantly more non-zeros than the original A RR  matrix, it is not explicitly computed. For preconditioning, the current art do either a column sum (Colsum) or a row sum (Rowsum) is performed to diagonalize the well influence matrix. In the Nested Factorization (NF) preconditioning algorithm, NF is applied to the matrix:
 
A RR −Colsum(A RW A WW   −1 A WR )  (18)
 
     The Colsum matrix is a block diagonal matrix with block (nvar*nvar) where nvar is the number of equations to be solved per reservoir grid cells. In the incomplete lower upper (ILU) triangular factorization preconditioning algorithm, ILU is applied to the matrix:
 
A RR −Rowsum(A RW A WW   −1 A WR )  (19)
 
     The Rowsum matrix is a block diagonal matrix with block (nvar*nvar) where nvar is the number of equations to be solved per reservoir grid cells. The └A RW A WW   −1 A WR ┘{V} matrix-vector product is computed as a series of matrix-vector multiplies whenever the product is needed. 
     Rowsum and Colsum preconditioning each require far fewer operations than preconditioning the full system matrix used earlier. However, if the wells penetrate many heterogeneous layers with high flow rates, the convergence of the solver is affected. In modern simulation, a full-field model may have thousands of wells, as can be seen in  FIGS. 2 and 3 , and each well can penetrate thousands of grid cells. In such field models the current Rowsum and Colsum preconditioning may have poor convergence or fail to converge for some problems. 
     Earlier preconditioning of the full system matrix as described above becomes prohibitively expensive while the other earlier method of direct solution of the fully coupled system is complex, also expensive, and difficult to parallelize for multiple processors. Thus, so far as is known, none of the prior methods have been sufficiently robust and efficient for large scale parallel reservoir simulation. 
     The Present Invention 
     To overcome the aforementioned difficulties, the present invention provides a computer implemented methodology which is robust, yet efficient to compute, based on and highly parallelizable preconditioner for the reduced system matrix. The present invention provides improvements to the existing technological processes of reservoir simulation used to analyze the performance of hydrocarbon reservoirs and evaluate their development. The present invention also improves the functioning of computers in performance of reservoir simulation in reducing the processing time lost due to poor convergence or failure to converge in the simulator processing, and also in reducing the complexity of communication within the network of the data processing system. In Equation (16):
 
 Ā =[ A   RR   −A   RW   A   WW   −1   A   WR ]  (20)
 
     The present invention takes the form of three embodiments. The first is based on a homogeneous one-level preconditioner. The second embodiment is a homogeneous two-level CPR-type preconditioner. The third embodiment is a heterogeneous combinative preconditioning method. The preconditioners used in the methodology of the three embodiments act as an accelerator for one of the Krylov subspace iterative algorithms such as the parallelized GCR(k), GMRES, or BICGSTAB methods. The present invention provides speedup in convergence of the coupled solution of the reservoir simulation. The present invention also provides an ability to generate good approximate update vectors to reduce the residuals of the coupled well-reservoir problem. Thus, the methodology of the present invention is fast, robust, and efficient to compute, highly parallelizable, and generate close approximate solutions. 
     Nomenclature 
     A RR =Reservoir Grid Jacobian Matrix 
     A WR =Well to Reservoir Coupling Jacobian Matrix 
     A RW =Reservoir to Well Coupling Jacobian Matrix 
     A WW =Multi-Segment Well Jacobian Matrix 
     A P =Pressure Coefficient Jacobian Matrix 
     Ā=Reduced System Jacobian Matrix 
     C=Pressure Coefficient Selector 
     C T =Transpose of the Pressure Coefficient Selector 
     E R =Remainder Component of Reservoir Grid Jacobian Matrix 
     E W =Well Influence Coefficient Jacobian Matrix 
     M −1 =Full-System Preconditioner 
     M CPR   −1 =Constraint Pressure Residual Preconditioner 
     M N   −1 =N-Term Line-Solve Power-Series Preconditioner 
     P R =Partitioned Component of Reservoir Grid Jacobian Matrix 
     R, r=Residual Vector 
     W=Pressure decoupling operator matrix 
     X=Solution Vector 
     Subscripts: 
     C=correction 
     P=pressure 
     R=reservoir 
     W=well 
     T=total 
     Embodiment 1 
     Parallel One-Level Line-Solve Power-Series (LSPS) Preconditioning with Variable Orders for Well-Coupling and Reservoir-Grid Non-Zeros Fill Terms 
     The homogeneous one-level preconditioning method is described as follows: In this method, the following substitutions are made:
 
 A   RR =[ P   R   +E   R ]  (21)
 
 E   W   =−A   RW   A   WW   −1   A   WR   (22)
 
     Therefore,
 
 Ā=A   RR   +E   W   (23)
 
     And the reduced-system residual:
 
 R =[ R   R   −A   R   A   WW   −1   R   W ]  (24)
 
     This gives the simplified form of the reduced system matrix equation:
 
 Ā·X   R =[ P   R   +E   R   +E   W ] X   R   =R   (25)
 
     Further:
 
 E   T   =E   R   +E   W   =E   R   −A   RW   A   WW   −1   A   WR   (26)
 
     And then:
 
 Ā =[ P   R   +E   T ]=[ P   R   +E   R   +E   W ]  (27)
 
     To obtain the familiar form of a linear system:
 
 Ā·X   R   =R   (28)
 
     In the above reduced system equations, P R  is the part of A RR  whose non-zeroes are used to generate the LU factors so that the result of P R   −1 ·V can be easily computed. For example, one suitable choice of P R  is the block tridiagonal matrix for a Z-Line ordered reservoir system of equations. 
     Another possible choice would be the maximum connection factor ordered block tridiagonal portion of the reservoir matrix. In a fractured dual-porosity dual-permeability system, or a fractured multi-modal porosity system, all the connection factors of all the pore space partitions at the same spatial location and the Z-line or the direction with the maximum connection factor are all part of P R . 
     An approximate inverse preconditioner using N-term power series, M N   −1 , can be written as: 
     
       
         
           
             
               
                 
                   
                     
                       
                         A 
                         _ 
                       
                       
                         - 
                         1 
                       
                     
                     ≈ 
                     
                       M 
                       N 
                       
                         - 
                         1 
                       
                     
                   
                   = 
                   
                     
                       [ 
                       
                         I 
                         + 
                         
                           
                             ∑ 
                             
                               K 
                               = 
                               1 
                             
                             N 
                           
                           ⁢ 
                           
                             
                               
                                 ( 
                                 
                                   - 
                                   1 
                                 
                                 ) 
                               
                               K 
                             
                             ⁢ 
                             
                               
                                 ( 
                                 
                                   
                                     P 
                                     R 
                                     
                                       - 
                                       1 
                                     
                                   
                                   ⁢ 
                                   
                                     E 
                                     T 
                                   
                                 
                                 ) 
                               
                               K 
                             
                           
                         
                       
                       ] 
                     
                     ⁢ 
                     
                       P 
                       R 
                       
                         - 
                         1 
                       
                     
                   
                 
               
               
                 
                   ( 
                   29 
                   ) 
                 
               
             
           
         
       
     
     Further, due to the different significance of the E T  components, if we keep terms for E R  and E W  up to different orders to save work: 
     
       
         
           
             
               
                 
                   
                     M 
                     N 
                     
                       - 
                       1 
                     
                   
                   = 
                   
                     
                       ( 
                       
                         
                           ∏ 
                           
                             K 
                             = 
                             
                               M 
                               + 
                               1 
                             
                           
                           N 
                         
                         ⁢ 
                         
                           ( 
                           
                             I 
                             + 
                             
                               
                                 ( 
                                 
                                   
                                     P 
                                     R 
                                     
                                       - 
                                       1 
                                     
                                   
                                   ⁢ 
                                   
                                     E 
                                     R 
                                   
                                 
                                 ) 
                               
                               
                                 2 
                                 
                                   K 
                                   - 
                                   1 
                                 
                               
                             
                           
                           ) 
                         
                       
                       ) 
                     
                     · 
                     
                       ( 
                       
                         
                           ∏ 
                           
                             K 
                             = 
                             2 
                           
                           M 
                         
                         ⁢ 
                         
                           ( 
                           
                             I 
                             + 
                             
                               
                                 ( 
                                 
                                   
                                     P 
                                     R 
                                     
                                       - 
                                       1 
                                     
                                   
                                   ⁢ 
                                   
                                     E 
                                     T 
                                   
                                 
                                 ) 
                               
                               
                                 2 
                                 
                                   K 
                                   - 
                                   1 
                                 
                               
                             
                           
                           ) 
                         
                       
                       ) 
                     
                     · 
                     
                       ( 
                       
                         I 
                         - 
                         
                           
                             P 
                             R 
                             
                               - 
                               1 
                             
                           
                           ⁢ 
                           
                             E 
                             T 
                           
                         
                       
                       ) 
                     
                     · 
                     
                       P 
                       R 
                       
                         - 
                         1 
                       
                     
                   
                 
               
               
                 
                   ( 
                   30 
                   ) 
                 
               
             
           
         
       
     
     Satisfying M≥2; N&gt;M; otherwise, the respective terms in the above matrix equation is the identity matrix [I]. Some special cases can be illustrated as follows: 
     For example, by keeping only the first order term for E W , then 
     
       
         
           
             
               
                 
                   
                     
                       M 
                       N 
                       
                         - 
                         1 
                       
                     
                     = 
                     
                       
                         ( 
                         
                           
                             ∏ 
                             
                               K 
                               = 
                               2 
                             
                             N 
                           
                           ⁢ 
                           
                             ( 
                             
                               I 
                               + 
                               
                                 
                                   ( 
                                   
                                     
                                       P 
                                       R 
                                       
                                         - 
                                         1 
                                       
                                     
                                     ⁢ 
                                     
                                       E 
                                       R 
                                     
                                   
                                   ) 
                                 
                                 
                                   2 
                                   
                                     K 
                                     - 
                                     1 
                                   
                                 
                               
                             
                             ) 
                           
                         
                         ) 
                       
                       · 
                       
                         ( 
                         
                           I 
                           - 
                           
                             
                               P 
                               R 
                               
                                 - 
                                 1 
                               
                             
                             ⁢ 
                             
                               E 
                               T 
                             
                           
                         
                         ) 
                       
                       · 
                       
                         P 
                         R 
                         
                           - 
                           1 
                         
                       
                     
                   
                   ; 
                   
                     N 
                     ≥ 
                     2 
                   
                 
               
               
                 
                   ( 
                   31 
                   ) 
                 
               
             
           
         
       
     
     If the same order for E R  and E W  are kept, then the preconditioner becomes 
     
       
         
           
             
               
                 
                   
                     
                       M 
                       N 
                       
                         - 
                         1 
                       
                     
                     = 
                     
                       
                         ( 
                         
                           
                             ∏ 
                             
                               K 
                               = 
                               2 
                             
                             N 
                           
                           ⁢ 
                           
                             ( 
                             
                               I 
                               + 
                               
                                 
                                   ( 
                                   
                                     
                                       P 
                                       R 
                                       
                                         - 
                                         1 
                                       
                                     
                                     ⁢ 
                                     
                                       E 
                                       T 
                                     
                                   
                                   ) 
                                 
                                 
                                   2 
                                   
                                     K 
                                     - 
                                     1 
                                   
                                 
                               
                             
                             ) 
                           
                         
                         ) 
                       
                       · 
                       
                         ( 
                         
                           I 
                           - 
                           
                             
                               P 
                               R 
                               
                                 - 
                                 1 
                               
                             
                             ⁢ 
                             
                               E 
                               T 
                             
                           
                         
                         ) 
                       
                       · 
                       
                         P 
                         R 
                         
                           - 
                           1 
                         
                       
                     
                   
                   ; 
                   
                     N 
                     ≥ 
                     2 
                   
                 
               
               
                 
                   ( 
                   32 
                   ) 
                 
               
             
           
         
       
     
     If N is set to 2 and only the first order terms for E W , then
 
 M   N   −1 =( I +( P   R   −1   E   R ) 2 )·( I−P   R   −1   E   T )· P   R   −1   (33)
 
     It is not necessary to build the approximate inverse matrix directly. The approximate inverse preconditioner for the reduced coupled fully-implicit system is applied as a series of matrix-vector multiply whenever an approximate solution update is needed. This embodiment of preconditioner addresses the smoothing of the error components arising from the well influence coefficients directly but at a moderate additional cost which can be control to the desired optimal level for speed and robustness of the coupled system solution. This is an advantage over techniques prior which can cost excessive solver iterations or convergence failure when there is a large complex simulation system with significant coupling between wells and reservoirs. As noted, the Rowsum or Colsum method has been too weak to be effective. 
     The preconditioned Krylov subspace algorithm for solving the reduced system Equation (25) is illustrated below using the generalized constraint residual GCR(k) method as follows:
 
Compute  r   o   =R−Ā·X   o   (34)
 
 Set p   o   =M   N   −1   ·r   o   (35)
 
     For j=0, 1 . . . , step until convergence DO: 
     
       
         
           
             
               
                 
                   
                     A 
                     _ 
                   
                   = 
                   
                     [ 
                     
                       
                         P 
                         R 
                       
                       + 
                       
                         E 
                         R 
                       
                       + 
                       
                         E 
                         W 
                       
                     
                     ] 
                   
                 
               
               
                 
                   ( 
                   36 
                   ) 
                 
               
             
             
               
                 
                   
                     α 
                     j 
                   
                   = 
                   
                     
                       ( 
                       
                         
                           r 
                           j 
                         
                         , 
                         
                           
                             A 
                             _ 
                           
                           ⁢ 
                           
                             p 
                             j 
                           
                         
                       
                       ) 
                     
                     
                       ( 
                       
                         
                           A 
                           _ 
                         
                         ⁢ 
                         
                           
                             p 
                             j 
                           
                           · 
                           
                             A 
                             _ 
                           
                         
                         ⁢ 
                         
                           p 
                           j 
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   37 
                   ) 
                 
               
             
             
               
                 
                   
                     x 
                     
                       j 
                       + 
                       1 
                     
                   
                   = 
                   
                     
                       x 
                       j 
                     
                     + 
                     
                       
                         α 
                         j 
                       
                       ⁢ 
                       
                         p 
                         j 
                       
                     
                   
                 
               
               
                 
                   ( 
                   38 
                   ) 
                 
               
             
             
               
                 
                   
                     r 
                     
                       j 
                       + 
                       1 
                     
                   
                   = 
                   
                     
                       r 
                       j 
                     
                     - 
                     
                       
                         α 
                         j 
                       
                       ⁢ 
                       
                         Ap 
                         j 
                       
                     
                   
                 
               
               
                 
                   ( 
                   39 
                   ) 
                 
               
             
             
               
                 
                   
                     Compute 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       β 
                       ij 
                     
                   
                   = 
                   
                     - 
                     
                       
                         ( 
                         
                           
                             
                               A 
                               _ 
                             
                             ⁢ 
                             
                               M 
                               
                                 - 
                                 1 
                               
                             
                             ⁢ 
                             
                               r 
                               
                                 j 
                                 + 
                                 1 
                               
                             
                           
                           , 
                           
                             
                               A 
                               _ 
                             
                             ⁢ 
                             
                               p 
                               i 
                             
                           
                         
                         ) 
                       
                       
                         ( 
                         
                           
                             
                               A 
                               _ 
                             
                             ⁢ 
                             
                               p 
                               i 
                             
                           
                           , 
                           
                             
                               A 
                               _ 
                             
                             ⁢ 
                             
                               p 
                               i 
                             
                           
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   40 
                   ) 
                 
               
             
           
         
       
     
     For i=s, . . . , j: 
     
       
         
           
             
               
                 
                   
                     p 
                     
                       j 
                       + 
                       1 
                     
                   
                   = 
                   
                     
                       
                         M 
                         N 
                         
                           - 
                           1 
                         
                       
                       · 
                       
                         r 
                         
                           j 
                           + 
                           1 
                         
                       
                     
                     + 
                     
                       
                         ∑ 
                         
                           i 
                           = 
                           k 
                         
                         j 
                       
                       ⁢ 
                       
                         
                           β 
                           ij 
                         
                         ⁢ 
                         
                           p 
                           i 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   41 
                   ) 
                 
               
             
           
         
       
     
     Enddo 
     For s=0, the above algorithm is a generalized constraint or GCR based. For s=max (0, j−k+1), the algorithm is ORTHOMIN based. GCR(k) is simply a restart version of GCR. The GMRES and BICGSTAB algorithms can also be used alternatively and similarly. 
     Flow chart  60  of  FIG. 6  illustrates the structure of the logic of the present invention for the parallel-level line-solve power series preconditioning embodiment as implemented in a data processing system or computer such as that shown at S in  FIG. 15 . Those skilled in the art will appreciate that the flow charts illustrate the structures of computer program code elements including logic circuits on an integrated circuit that function according to this invention. Manifestly, the invention is practiced in its essential embodiment by a machine component that renders the program code elements in a form that instructs a digital processing apparatus (that is, a computer) to perform a sequence of function steps corresponding to those shown. 
     It is important to note that, while the present invention has been, and will continue to be, described in the context of a fully functional computer system, those skilled in the art will appreciate that the present invention is capable of being distributed as a program product in a variety of forms, and that the present invention applies equally regardless of the particular type of non-transitory signal-bearing media utilized to actually carry out the distribution. Examples of non-transitory signal-bearing media include: recordable-type media, such as floppy disks, hard disk drives, and CD ROMs. 
     It should be understood that the processing described herein can be implemented in a variety of other types of reservoir simulators. It can be run on a variety of computer platforms, such as single CPU, a shared memory parallel or massively parallel processing computer, a distributed memory super-computer, and a variety of PC clusters, such as a self-made PC cluster, or a production PC cluster. 
     In the sequence of  FIG. 6 , a flow chart  60  indicates a sequence of processing steps according to the present invention. The processing sequence is entered at step  61  from a main reservoir simulator processing sequence of the type described herein. During a first processing step  61  a composite Jacobian matrix according to Equation (15) is constructed for the nonlinear well and reservoir equations. 
     Then, during a step  63 , an iteration counter i is set to 0, and an estimate xi and is selected so that a matrix calculation
 
 r   i   =b −[ A ] x   i  
 
may be made.
 
     In a step  65 , parallel one-level line-solve power-series (LSPS) with variable orders for well-coupling and reservoir-grid non-zeros fill terms is applied as described above in connection with Equations (29) through (33) where equations (30) is the general form of the preconditioner described in equation (29) and equations (31), (32), (33) are specific examples of the forms of the preconditioner. 
     In step  65 , the approximate solution vector X R =Ā −1 ·R≈M N   −1 ·R is according to Equation (28), the components of which are detailed by equations (21) through (27) which involves the fully-coupled fully-implicit well influence Jacobian matrix. The approximate inverse preconditioner is applied within the Krylov subspace iterative method, as will be described. 
     This computation is performed using the truncated Neumann series expansion in the manner described above. If a parallel computation array is being used, communication between adjacent processors is done during step  65 , in a manner to be described below. 
     Next, during step  67 , a truncated general conjugate residual or GCR(k) method is applied to solve the system of well and reservoir equations, as described above in connection with the methodology of Equations (34) through (41). Again if a parallel computation array is in use, communication between adjacent processors is performed. Besides GCR(k), alternative Krylov subspace iterative methods, such as parallel versions of GMRES or BICGSTAB methods, may also be used in this step. 
     In step  69 , the residual result obtained during step  67  is compared with a user specified tolerance of residual. If the obtained residual result is not within the user specified tolerance, the iteration counter is incremented and processing returns to step  65  for another processing cycle in the manner described above for subsequent iterations. The well solution vector is computed using Matrix equation (17) after the reservoir solution vector has been computed. 
     Processing proceeds in the foregoing manner until a solution value for the residual is found within the user specified tolerance range during step  67  of an iteration cycle. At this time, the results obtained for the well and reservoir equations of the fluids yielding the satisfactory residual in the grid cells are stored in memory of the processor or processors proceeding according to the instruction steps of  FIG. 6 . Control of the processor or processors then reverts to the main reservoir simulator sequence. Displays of the results of satisfactory processing according to  FIG. 6  are made on user request. 
     The procedure is the same for both the serial application and the parallel application of the method. In the parallel application of the method, an interprocessor communication step is performed prior to the matrix-vector multiplication where the elements of the intermediate solution vector situated at the boundary of a data partition need to be exchanged between the processors sharing an internal domain boundary. 
     Embodiment 2 
     Parallel CPR-Type Line-Solve Power-Series Preconditioning 
     What is known as the constraint pressure residual, CPR, preconditioning method was first described in Wallis, J. R. et al.: “Constrained Residual Acceleration of Conjugate Residual Methods,” SPE 13563, Proceedings of the 8th SPE Reservoir Simulation Symposium, Dallas, USA Feb. 10-13, 1985. A variant of the CPR method was further discussed in Fung &amp; Dogru “Parallel Unstructured Solver Methods for Simulation of Complex Giant Reservoirs,” SPE (December-2008), pp. 440-446. 
     The constraint pressure residual procedure is adapted and improved to solve the coupled well-reservoir system, with the matrix terms used being defined according to the Nomenclature definitions. The CPR preconditioning according to the present invention involves a pressure predictor-corrector step and can be written as:
 
 M   CPR   −1   =M   −1 [ I−ÃCA   p   −1   C   T ]+( CA   p   −1   C   T )  (42)
 
     Where:
 
 Ã=W*A   (43)
 
     And the pressure matrix is
 
 A   p   =C   T   ÃC   (44)
 
     C is given by 
     
       
         
           
             
               
                 
                   C 
                   = 
                   
                     [ 
                     
                       
                         
                           
                             e 
                             p 
                           
                         
                         
                           
                               
                           
                         
                         
                           
                               
                           
                         
                         
                           
                               
                           
                         
                       
                       
                         
                           
                               
                           
                         
                         
                           
                             e 
                             p 
                           
                         
                         
                           
                               
                           
                         
                         
                           
                               
                           
                         
                       
                       
                         
                           
                               
                           
                         
                         
                           
                               
                           
                         
                         
                           ⋱ 
                         
                         
                           
                               
                           
                         
                       
                       
                         
                           
                               
                           
                         
                         
                           
                               
                           
                         
                         
                           
                               
                           
                         
                         
                           
                             e 
                             p 
                           
                         
                       
                     
                     ] 
                   
                 
               
               
                 
                   ( 
                   45 
                   ) 
                 
               
             
           
         
       
     
     Assuming pressure is a first unknown for each cell, then, let e p  be the nvar×1 vector where nvar is the number of equations per grid cell: 
     
       
         
           
             
               
                 
                   
                     e 
                     p 
                   
                   = 
                   
                     [ 
                     
                       
                         
                           1 
                         
                       
                       
                         
                           0 
                         
                       
                       
                         
                           ⋮ 
                         
                       
                       
                         
                           0 
                         
                       
                     
                     ] 
                   
                 
               
               
                 
                   ( 
                   46 
                   ) 
                 
               
             
           
         
       
     
     The aim of W is to perform an IMPES-like reduction step on A. For example, W can be computed as:
 
 W =DIAG −1 ( A )  (47)
 
     The notation DIAG(A) denotes the main diagonal nvar*nvar sub-matrix blocks of A. Note that W needs to be chosen such that Ã p  is close to being positive definite for PAMG to work well as the pressure solver. This is not a requirement for LSPS or any of the ILU variants. 
     Thus, the CPR preconditioning step x=M CPR   −1 r can be outlined as follows:
         1. Restrict the full system residual to the pressure system
 
 r   P   =C   T   Wr   (48)
   2. Iteratively solve the pressure system
 
 A   p   x   p   =r   p   (49)
   3. Expand pressure solution to full system
 
 s=Cx   p   (50)
   4. Correct the full system residual from pressure residual
 
 r   c   =r−A·s   (51)
   5. Solve the full system using second stage preconditioning
 
 Mx=r   c   (52)
   6. Correct full system solution from pressure solution
 
 x   c   =x+s   (53)
   7. The pressure system preconditioning can be accomplished in a number of ways, including the use of a parallel algebraic multi-grid or PAMG method, or the ILU(k) or the LSPS method as described in Embodiment 1 above.       

     The preferred full system preconditioning method is the LSPS method. The pressure system solver is normally limited to a few Krylov subspace iterations of the GCR(k) or the GMRES, or the BICGSTAB method. For the PAMG method, it is usually limited to a single local vector V multi-grid cycle. The full-system preconditioned solution in Equation (52) also uses one of the Krylov subspace method as outlined in Equations (34) through (41). 
       FIG. 7  shows a schematic flow diagram  70  for a solver method for parallel constraint pressure residual or CPR-type line-solver power series preconditioning according to the second embodiment of the present invention. The approximate inverse preconditioning method used in determining an approximate solution vector involves the fully-coupled well influence matrix E W =−A RW A WW   −1 A WR  which is applied within one of Krylov subspace iterative algorithms during both the pressure solution and full-system solution. The well solution vector is computed using Matrix equation (17) after the reservoir solution vector has been computed. The processing shown in  FIG. 7  is performed in the data processing systems of  FIG. 15 . 
     Step  71  as indicated in  FIG. 7  scales the block rows of the system matrix by the block diagonal inverse. Since data store for the entire block rows of each subdomain is local, this step is fully parallel and requires no synchronization. Step  72  as indicated, computes upper and lower decomposition [L] [U] matrices for [P]. The [P] matrix is the partitioned Jacobian matrix [P R ] described in equation (25) and the inverse [P R   −1 ] as in equations (29) to (33) is computed as [U] −1 [L] −1 . 
     Step  73  extracts the pressure coefficients matrix from the full-system coefficient matrix. [C] and [C] T  are the restriction and prolongation operators in the algebraic extraction process. Step  74  simply initializes the full-system iteration process. 
     Step  75  extracts the pressure residual from the full system residuals. Step  76  performs an approximate pressure solution. The solution method is the preconditioned Krylov subspace method. The preconditioner can be the parallel line-solve power-series method as shown in equation (29) or a parallel algebraic multigrid method PAMG. The pressure solution only requires a rough approximate solution. Step  76  includes the MPI communication method at each matrix-vector multiplication. This method includes communication hiding strategy to enhance parallel scalability. The pressure solution is used to constrain the full system residual to accelerate the full-system solution. 
     Step  77  expands the pressure solution vector to a full-system solution vector. This is used to reduce the full-system residual in step  78 . In step  79 , the reduced full-system residual is then used to compute an approximate full-system update using the matrix equation as indicated. The combined update from both the pressure system and full system is computed in step  80 . Step  81  is the full-system version of the parallel GCR(k) or GMRES(k) algorithm which optimizes on the full-system update such that the residuals are minimized. Step  82  is the full-system convergence check. The i-indexed iterative loop from steps  74  to  82  is the full-system iterative solution loop. 
     Steps  76 ,  78 ,  79 , and  81  use embedded matrix-vector multiplication which involves MPI communication with communication hiding algorithm. The matrix-vector multiplications for interior cells are done simultaneously with the MPI communication. When communication is complete, local vector V now contains all ghost cell information and the matrix-vector multiplication can be done for the remaining boundary cell. Steps  76  and  81  involve distributed vector dot product which requires MPI reduction operations. The other steps are fully parallel steps. 
     Embodiment 3 
     Parallel Combinatory Preconditioning 
     This heterogeneous preconditioning method uses the preconditioning according to Embodiment 1 above for the well influence coefficient matrix, but includes a different method for the reservoir coefficient matrix. This is helpful if it is desired to apply different preconditioners for the reservoir terms and the well-influence terms. The combinatory preconditioner is written as:
 
 A   −1   ≈M   C   −1   =M   R   −1 [ I−A ·[ M   W   −1 ]]+ M   W   −1   (54)
 
     The general form of the well-influence-term preconditioner can be written as: 
     
       
         
           
             
               
                 
                   
                     
                       
                         M 
                         N 
                         
                           - 
                           1 
                         
                       
                       = 
                       
                         
                           ( 
                           
                             
                               ∏ 
                               
                                 K 
                                 = 
                                 2 
                               
                               N 
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               ( 
                               
                                 I 
                                 + 
                                 
                                   
                                     ( 
                                     
                                       
                                         P 
                                         R 
                                         
                                           - 
                                           1 
                                         
                                       
                                       ⁢ 
                                       
                                         E 
                                         W 
                                       
                                     
                                     ) 
                                   
                                   
                                     2 
                                     
                                       K 
                                       - 
                                       1 
                                     
                                   
                                 
                               
                               ) 
                             
                           
                           ) 
                         
                         · 
                         
                           ( 
                           
                             I 
                             - 
                             
                               
                                 P 
                                 R 
                                 
                                   - 
                                   1 
                                 
                               
                               ⁢ 
                               
                                 E 
                                 W 
                               
                             
                           
                           ) 
                         
                         · 
                         
                           P 
                           R 
                           
                             - 
                             1 
                           
                         
                       
                     
                     ; 
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     N 
                     ≥ 
                     2 
                   
                 
               
               
                 
                   ( 
                   55 
                   ) 
                 
               
             
           
         
       
     
     Frequently, only a low order approximation is needed to improve convergence. The M R   −1  may be another effective preconditioner for the reservoir terms. For example, the ILU(k) or ILUT for the reservoir part of the A matrix with a suitable domain decomposition method for parallel application may be chosen here. If M R   −1  is the ILU(k) preconditioner, then:
 
M R   −1 =M ILU   −1   (56)
 
     Processing according to Embodiment 3 of the present invention in data processing system S is performed according to the methodology of  FIG. 6 , with the exception that for the preconditioner applied according to step  65  is according to Equations (54) through (56) above. 
     Again, a two-stage CPR preconditioner can be equally applied in the parallel combinatory preconditioning. In this application, the pressure system preconditioning method can be parallel algebraic multigrid or PAMG, while the full-system preconditioning method can be the parallel combinatory preconditioning method. This method is used as the preconditioning method in a Krylov subspace iterative algorithm. 
     One such algorithm GCR(k) was discussed previously in Equations (34) through (41). The Krylov subspace method using the GCR(k) algorithm includes an approximate inverse preconditioning method for calculating the approximate solution vector which involves the fully-coupled fully-implicit well influence Jacobian matrix E W =−A RW A WW   −1 A WR . The approximate inverse preconditioner is applied within the Krylov subspace iterative algorithm. 
     Distributed Parallel Method for Fully-Coupled Fully-Implicit Well-Reservoir Solution 
       FIGS. 8 and 9  illustrate schematically organization of grid data for distributed parallel processing of the unstructured subdomains of the reservoir grid of Domain D of  FIG. 1  in order to process the data for the coupled implicit well bores such as shown at  20  and  22  in  FIGS. 2 and 3 , respectively. The data for distributed parallel processing for the reservoir domain is partitioned for the structured-grid and unstructured-grid into multiple subdomains, with each subdomain assigned to a different compute process is described in U.S. Pat. Nos. 8,433,551 and 8,386,227, respectively, in each of which Applicant is named as an inventor. The distributed parallel processing for the fully-coupled fully-implicit well-reservoir solution is accomplished in the following manner. 
     For each of the wells  20  and  22 , the full set of all the active wells in the model are distributed to the participating compute processes of the data processing system S. For example, as a simplified model, where there are four wells and four compute processes, each well is run on a different compute core of the data processing system S, and the compute cores jointly solves the fully-coupled simulation model of the reservoir grid blocks and the wells. Thus, each processing core and its associated data storage works on and stores data for each one of the wells. 
     In another example, if there are two wells but four processing cores for as shown in  FIG. 8 , then, processor cores for two processors receive data for the two wells, respectively, but the remaining cores do not have the ownership of any processing for the wells.  FIG. 8  shows four subdomains  88  into which the reservoir grid blocks are partitioned, with two complex wells  20 , penetrating certain ones of the subdomains  88  as indicated. The reservoir grid of domain D has two complex wells  20 . Well  1  intersects grid blocks in subdomains  1 ,  3 , and  4 . Well  2  intersects grid blocks in subdomain  2 .  FIG. 9  illustrates schematically the allocation and distribution of matrix data and vector data for the solution of the well-reservoir coupled model of  FIG. 8 . The matrix shown in  FIG. 9  is comparable to those shown in  FIGS. 4 and 5 , wherein  FIGS. 4 and 5  illustrate the matrix dense blocks in detail for a serial instance, but  FIG. 9  shows the high-level schematic layout where each named block represents a matrix of many internal dense blocks. The labels on the block indicate the types of dense blocks contained within. For example, the label W 1 .D 1  indicates the dense blocks contained inside are to do with the connections between the well  1  and subdomain  1 , and so on. Data for each of the subdomains  1 ,  2 ,  3 , and  4  are allocated to processing cores P 1 , P 2 , P 3 , and P 4 . Well  1  is allocated to P 1  and well  2  is allocated to P 2 . 
     The vector shown in  FIG. 9  represents either the residual vector, the solution vector, or other block-length vectors in the solution process. In  FIG. 9 , each long horizontal box  90  schematically represents the collection data belonging to a compute process P 1 , P 2 , P 3  or P 4 . The well data groups W 1 , W 1 .D 1 , W 1 .D 3 , and W 1 .D 4  of  FIG. 9  represent the well-segment internal data and well-segment to well-segment data, as well as well-to-grid data stored in process P 1 . Well data group W 1 .D 1  represents the matrix data of the interaction terms between well segments belonging to well  1  and the perforated grid cell data belonging to subdomain  1 , and so on. Grid data groups D 1 .W 1 , D 1 , D 1 .D 2 , and D 1 .D 3  are the grid-to-well data, grid cell internal data, grid-to-grid interaction data within subdomain  1  and between subdomain  1  and another subdomain, in this example subdomains  2  and  3 . The grid data group D 1 .D 2  shown in  FIG. 9  represents the inter-block flow-term derivatives in the Jacobian matrix of the grid cells on the boundary of subdomain  1  adjacent to the boundary of subdomain  2 . The inter-block flow-term derivatives of grid data group D 1 .D 2  are also stored in process P 1 . 
     Inter-process communications are provided according to the present invention to construct part of the Jacobian matrix representing boundaries of subdomains or well-to-grid grid-to-well interaction terms which do not belong to the same compute process. Inter-process communications are also provided to solve the resulting system of algebraic equations. For example, the well-to-grid Jacobian term W 1 .D 3  shown in  FIG. 9  indicates the transfer of grid block data in P 3  which are penetrated by well  1  to process P 1  to calculate them and data for well  1  from process P 1  is need by process P 3  to do processing for D 3 .W 1  and so on. The coupling terms involving information partially owned by another process thus has those data transferred locally to complete the construction of the coupling term. 
     Similarly, inter-process communications for grid-to-grid flow term for grid blocks located at the boundaries of subdomains are also provided. For the construction of the grid data group D 2 .D 1  and D 1 .D 2  for grid cells adjacent the respective subdomain boundaries, data exchange between processes P 1  and P 2  is provided. Similarly, for construction of grid data groups D 1 .D 3  and D 3 .D 1 , data exchange between processes P 1  and P 3  is provided, and so on. 
     In the solution process, data elements in the vector which are required by the matrix-vector operation are inter-process communicated to the processor P which has processing ownership of the matrix components. For example, vector components from VD 3  and VD 4  are communicated to process P 1  for matrix-vector operations with the matrix components W 1 .D 3  and W 1 .D 4 . The required vector components from VW 1  are also communicated to P 3  and P 4  for matrix-vector operation with D 3 .W 1  and D 4 .W 1 . The vector can be an intermediate solution vector, or a residual vector, or an intermediate product resulting from a partially completed matrix-vector multiplication sequence. 
     With the present invention, as shown in  FIG. 9 , communication hiding is achieved by organizing the data so that the grid blocks and connections which do not share boundaries with other subdomains are ordered first. They are followed by boundary grid cells belonging to this subdomain called the inner halo and finally boundary grid cells belonging to neighboring subdomain called the outer halo. Similarly, internal cell-to-cell connections and their connection factor are ordered first and the cell-to-cell connections between a cell in this subdomain and a cell belonging to a neighboring subdomain are organized at the end in consecutive banks of each neighboring subdomains. Computation for internal grid blocks and connections and communication can be carried out simultaneously while the computation for the inter-subdomain connection can only start after the inter-process communication has been completed. 
     The processing and method steps of  FIGS. 6 and 7  are suitable for deployment on a variety of today&#39;s HPC hardware. These are typically rack mounted hardware with several compute nodes which contains multiple CPUs with multi-core architecture. Nodes are interconnected with conventional low latency high bandwidth networks, switches, and routers. 
     The typical HPC environment for use with this simulation system is today&#39;s multi-node, multi-CPU, multi-core compute clusters. An example such cluster is illustrated at C in the data processing system S of  FIGS. 15 and 16 . The cluster C is formed of a plurality of computer nodes  150  ( FIGS. 15 and 16 ) which are provided with data in parallel as indicated by and arrow  152  by a router server or servers  154 . If desired, several such router servers may be used for this purpose. Original simulation or input data of the types described above is stored in a suitable number of data storage/file servers  156 . The router servers  154  under control of computer code  155  stored in memory transfer input simulation data in parallel from the storage servers  156 , as well as simulation processing results as indicated by an arrow  158 , to and from the computer nodes  150  of the cluster C. The program code  155  according to the present invention is in the form of non-transitory computer operable instructions causing the server or servers  154  to index, order and transfer the data. Typically, the data processing system D includes a set of work stations  157  of suitable, conventional type which are connected to the system by a network  159 . 
     The computer nodes  150  of the cluster C include a plurality of processors or cores  160  of the type illustrated in  FIG. 6  operating in parallel under instructions of a computer code or program product  162  stored in stored in memory  164  of the computer nodes  150 . The program code  162  according to the present invention is in the form of non-transitory computer operable instructions causing the data processors  160  to simulate fluid flow in a fractured subterranean reservoir where the pore space has been characterized as a multi-continuum. 
     It should be noted that program codes  155  and  162  may be in the form of microcode, programs, routines, or symbolic computer operable languages that provide a specific set of ordered operations that control the functioning of the data processing system D and direct its operation. The instructions of program codes  155  and  162  may be stored in memory of the servers  154  or processor nodes  150 , or on computer diskette, magnetic tape, conventional hard disk drive, electronic read-only memory, optical storage device, or other appropriate data storage device having a non-transitory computer usable medium stored thereon. Program code  160  may also be contained on a data storage device such as server  156  as a computer readable medium, as shown. 
     RAM and cache memory are distributed and local to each compute node and are shared by the processing cores on each the node. The physics simulated by the system of the present invention is a tightly coupled, global multiphase flow problem which is both convective and diffusive in nature. A high bandwidth, low latency network is thus preferred to minimize inter-process communication overhead. The message passing interface (MPI) standard is used for inter-process communication operations while MPI- 2  is used for parallel I/O operations. Disk storage for simulation or model data and processing output results are typically on centralized NAS, SAN, GPFS, or other parallel file systems. For smaller scale parallelism, local hard disk storage which resides on the cluster can also be used. Parallel distributed I/O methods are used to minimize read/write time from/to disk during simulation. 
       FIGS. 11, 12, 13 and 14  are comparative data plots of processing of the reservoir simulation model R shown in  FIG. 10  according to the present invention against processing according to the prior Rowsum/Colsum preconditioning described above. The reservoir model was an 11.9 million grid-cell (399*607*49), three-phase, black-oil reservoir model with 2,095 wells.  FIG. 11  illustrates comparative results of static well pressure over time,  FIG. 12  illustrates comparative results of borehole pressure or BHP over time, and  FIGS. 13 and 14  illustrate comparative results over time of oil rates and water cuts, respectively. 
     In  FIGS. 11 through 14  it is shown that the fully-implicit well solution (FIW) produces expected static well pressure and borehole pressure value, whereas a prior-art method (SIW) did not give the correct values. Further, the fully-implicit well solution (FIW) and a prior art (SIW) produces comparable results for oil rates and water cuts for this particular well. 
     With the present invention, it is preferable to use the solver preconditioner by generating the [L][U] decomposition of [P] in the initial stage of the solve invocation. Whenever the inverse is needed during the solver preconditioning operations expressed in Equations (29), (42), or (54), a forward-backward substitution is used to generate the solution. Equations (30) through (33) are instances of the preconditioner Equation (29) where the number of series terms may be different for the matrices E R  and E W . The matrix-vector operations are done in a sequence from right to left as the computation and inter-process communication required to generate the required approximate solution vector. 
     From the foregoing, it can be seen that the present invention provides a methodology which handles the well influence coefficient matrix at comparable accuracy and robustness to those for the grid-to-grid flow terms. The method is highly parallelizable and runs faster than prior methods, as it requires fewer solver iterations to converge to the same tolerances. 
     The present invention provides parallel solution for fully-coupled fully-implicit wellbore modeling in reservoir simulation which is robust, and particularly applicable for complex wells with thousands of interaction terms with the reservoir grid cells. It is highly parallelizable and well suited for implementation in modern HPC hardware. 
     The invention has been sufficiently described so that a person with average knowledge in the field of reservoir modeling and simulation may reproduce and obtain the results mentioned in the invention herein. Nonetheless, any skilled person in the field of technique, subject of the invention herein, may carry out modifications not described in the request herein, to apply these modifications to a determined structure and methodology, or in the use and practice thereof, requires the claimed matter in the following claims; such structures and processes shall be covered within the scope of the invention. 
     It should be noted and understood that there can be improvements and modifications made of the present invention described in detail above without departing from the spirit or scope of the invention as set forth in the accompanying claims.