EP2791858A1

EP2791858A1 - System and method for simulation of gas desorption in a reservoir using a multi-porosity approach

Info

Publication number: EP2791858A1
Application number: EP11877491.8A
Authority: EP
Inventors: John Edwin KILLOUGH
Original assignee: Landmark Graphics Corp
Current assignee: Landmark Graphics Corp
Priority date: 2011-12-16
Filing date: 2011-12-16
Publication date: 2014-10-22
Also published as: AR089218A1; AU2011383286A1; AU2011383286B2; MX2014006667A; BR112014014677A2; CA2858319C; CA2858319A1; US20140350906A1; MX341255B; CN103999093A; WO2013089784A1; EP2791858A4

Abstract

A hydrocarbon shale reservoir simulation system and method represented by a model having at least four different types of porosity nodes is described. The method includes the computer-implemented steps of characterizing porosity nodes within the model as one of natural fracture pore systems, matrix pore systems, induced fracture pore systems or vug pore systems. Following characterization, transfer terms between nodes are identified. Transfer terms may include transfer terms between vug nodes, matrix nodes, natural fracture nodes and induced fracture nodes. Once transfer terms have been assigned, the linear system for the model can be solved utilizing a linear solver. The method further includes the steps of utilizing the characterized pore nodes to define one or more subgrids that represent a zone within the reservoir, wherein the zone includes at least one node of each porosity type; and wherein the linear solver is applied by subgrid or associated subgrids.

Description

SYSTEM AND METHOD FOR SIMULATION OF GAS DESORPTION IN A RESERVOIR USING A MULTI-POROSITY APPROACH

Cross-Reference to Related Application(s)

This application is a continuation-in-part, and claims priority to, PCT Application No. PCT/US201 1/44178, filed July 15, 201 1 , assigned to the assignee of the present application, and hereby incorporated by reference in its entirety.

Background

Reservoir simulation is an area of reservoir engineering that employs computer models to predict the transport of fluids, such as petroleum, water, and gas, within a reservoir. Reservoir simulators are used by petroleum producers in determining how best to develop new fields, as well as in generating production forecasts on which investment decisions are based in connection with developed fields.

Fractured reservoirs present special challenges for simulation because of the multiple porosity systems or structures that may be present in these types of reservoirs. Fractured reservoirs are traditionally modeled by representing the porous media using two co-exiting pore systems or structures interconnected by flow networks, in what is referred to as dual porosity analysis. One type of pore system used in the prior art is the rock matrix, defined with matrix nodes, is characterized by high pore volume and low conductivity. The other type of pore system used in the prior art are induced fractures, and defined with fracture nodes, is characterized by low pore volume and high conductivity. These prior art reservoir simulation methods and systems typically treat absorbed gas within the reservoir as residing in the rock matrix pores of the reservoir. For example, in one simulated representation, referred to as dual-porosity, single-permeability ("DPSP"), matrix simulation nodes communicate only with fracture simulation nodes, and the analysis focuses on mass transfer and fluid flow of hydrocarbons between matrix nodes and fracture nodes. In DPSP, fracture nodes can also communicate with other fracures , which communicate with both matrix simulation nodes as well as other fracture simulation nodes. In another simulated representation, referred to as dual-porosity, dual-permeability

("DPDP"), matrix simulation nodes communicate with both fracture simulation nodes and as well as other matrix simulation nodes, and the analysis focuses on mass transfer and fluid flow of hydrocarbons between matrix nodes and fracture nodes as well as between matrix nodes and other matrix nodes. Those of ordinary skilled in the art will appreciate that "nodes" as used herein refer to the an elemental representation of pore structures within a simulated reservoir., while "zones" refer to a collection nodes within the simulated reservoir. Unknowns such as pressures and composition are solved for, typically on a node by node basis, at desired time and/or depth increments.

One particular type of reservoir encountered in oil and gas reservoir simulation is a shale reservoir. Shale reservoirs typically include large pores or vugs. Vugs are pore spaces that are comparatively larger than pore spaces of the rock matrix. Kerogen resides in this system of vugs within the porous rock matrix. Vugs may or may not be connected to one another. "Separate vugs" are vugs that are interconnected only through the interparticle porosity, i.e., the rock matrix porosity, and are not interconnected to one another (as are matrix pore volumes and fracture pore volumes). "Touching vugs" are vugs that are interconnected to one another. Because of their separate physical and mechanical characteristics, the fluid retention and transport properties of vug pore systems are different from those of both the matrix and fracture systems, and have not heretofore been adequately addressed with analysis utilizing only matrix porosity systems and induced porosity systems. In other words, because of the geologic complexities of shale reservoirs, traditional dual porosity reservoir modeling techniques do not adequately predict mass transfer and fluid flow characteristics of shale reservoirs.

Brief Description of the Drawings

A more complete understanding of the present disclosure and advantages thereof may be acquired by referring to the following description taken in conjunction with the accompanying figures, wherein: Fig. 1 illustrates an example of a reservoir simulation model comprising multiple wells.

Fig. 2 illustrates a representation of an example formation comprising a complex network of artificially-induced fractures. Fig. 3 illustrates a simulation grid of a formation comprising a highly deviated wellbore surrounded by natural fractures and a complex network of artificially-induced fractures.

Detailed Description To overcome the above-noted and other limitations of the current approaches, one or more embodiments described herein comprise a reservoir simulator including a unique manner of handling gas desorption in shale gas reservoir simulations by rigorously simulating the flow mechanism that occurs therein.

It has been found that the mechanism for desorption of gas in a shale gas reservoir is based on the existence of four separate porosity systems, each of which is incorporated in the method and system of the invention. In the method and system of the invention, each of these four porosity systems is separately characterized and incorporated into the model. The four porosity systems are the matrix porosity system, the induced fracture porosity systems, the natural fracture porosity system and the vug porosity system. As explained above, heretofore, only the matrix porosity system and the induced fracture porosity systems have been used in reservoir modeling in the past. The method and system of the invention incorporate natural fracture porosity systems and vug porosity systems.

The innermost of these porosity systems is the kerogen vugs, which contain the gas saturation as wetting fluid. The other porosity systems, which are the rock matrix, the induced fracture network and the natural fracture network, function as conduits for the gas contained in the kerogen of the shale. Rather than residing in pores throughout the porous rock matrix, the adsorbed gas is generally found only in the kerogen vugs. Natural fractures exist near the vugs, which natural fractures may or may not be open. The framework rock matrix of the porous medium connects the complex natural fractures to the hydraulic induced fractures near the well. The only accurate way to treat this system in accordance with the flow characteristics thereof, and in particular to account for the vugs throughout the matrix, requires a multiple porosity simulation system in which the vuggy portions of the formation containing the kerogen are connected to the rock matrix and the natural fracture system. The matrix and natural fractures are connected to the induced fractures near the wellbore, which fractures will in turn have different properties than the natural fractures due to the presence of fracking fluid and perhaps proppant, and hence one reason why natural fractures pore systems and induced fracture pore systems are characterized and separately analyzed in the method and system of the invention.

The simulation system described in the aforementioned PCT Application No. PCT/US201 1/44178 provides a unique tool for simulating general multi-porosity systems in which fluid flow through several porosity systems is modeled using various equations and connectivities in accordance with characteristics of the porosity systems. The embodiments described herein utilize this feature to simulate in a unique fashion the desorption of gases from shale gas reservoirs. The equations provided hereinbelow are transfer functions derived from field observations and laboratory measurements of the desorption process from the kerogen vugs, matrix and fractures of a reservoir. The transfer functions are then utilized by the simulation system to simulate the complex fracture system for the shale as reservoir coupled with the kerogen desorption. In one embodiment of the invention, vuggy porosity is used to model the kerogen desorption from within a complex fracture system comprised of both induced and natural fractures. Fig. 1 is a block diagram of an exemplary computer system 100 adapted for implementing a reservoir simulation system as described herein. In one embodiment, the computer system 100 includes at least one processor 102, storage 104, I/O devices 106, and a display 108 interconnected via a system bus 109. Software instructions executable by the processor 102 for implementing a reservoir simulation system 1 10 in accordance with the embodiments described herein, may be stored in storage 104. Although not explicitly shown in Fig. 1, it will be recognized that the computer system 100 may be connected to one or more public and/or private networks via appropriate network connections. It will also be recognized that the software instructions comprising the reservoir simulation system 1 10 may be loaded into storage 104 from a CD-ROM or other appropriate storage media.

In one embodiment, a portion of the reservoir simulation system 1 10 is

implemented using reservoir simulation software. In this embodiment, a "subgrid" data type is used to offer a generalized formulation design. In one embodiment, this data type may be Fortran. The subgrid defines the grid domain and interconnectivity properties of the nodes of the various porosity structures. It also tracks various node variables, such as pressure, composition, fluid saturation, etc. Subgrids are designated as being of a particular porosity type, e.g., natural fracture, matrix, induced fracture and vug. The nodes that constitute these grids are correspondingly referred to as natural fracture nodes, matrix nodes, induced fracture nodes and vug nodes. Subgrids of different porosity types occupying the same physical space are said to be "associated". Connections between porosity types, and in particular, the nodes of the porosity types, are represented as external connections, subgrid to associated subgrids. Internal (or intragrid) connections, and in particular, the nodes of a subgrid, represent flow connections within a porosity type.

The modeling of a shale gas reservoir generally involves defining one or more elongated, highly deviated production wellbore, typically thousands of feet in length, with multiple hydraulic fracture zones disposed substantially perpendicular to the wellbore, depending on the stress field in the formation. For certain formations, the stress field is such that a complex fracture system is induced between the large fractures emanating from the well. One representation of such fractures for an example formation is presented in Fig. 2 and designated by a reference numeral 200. The representation 200 has been derived from a finite element model of the porous media of the formation following high pressure injection of fracture fluid and proppant. Heavier lines, such as those designated by reference numeral 202, represent fractures induced by hydraulic fracturing, as described above, and which have been modeled in the prior art, i.e., induced fracture porosity systems. Narrower lines and triangular features, such as those designated by reference numerals 204 and 206, respectively, represent a possible finite volume grid with which to model the flow of fluid (primarily gas and water) in the complex fracture network and eventually to a horizontal production wellbore via the induced fracture, as illustrated in Fig. 3. Specifically, Fig. 3 illustrates a simulation grid 300 for an elongated, substantially horizontal wellbore 302, surrounded by induced fractures 304 and a complex natural fracture network 306. It has been found that there are two features of the physics in reservoir simulation that are required to properly model shale gas flow to a wellbore in a reservoir: non-Darcy flow and gas desorption. In the method and system of the invention, these two features are considered when modeling mass transfer and fluid flow between the matrix nodes, natural fracture nodes, induced fracture nodes and vug nodes. Non-Darcy flow is fluid flow that deviates from Darcy's assumption that fluid flow in the formation will be laminar. Non-Darcy flow is typically observed in high-velocity gas flow induced pressure differentials between the formation and the wellbore. Specifically, when the flow at the wellbore reaches velocities in excess of the Reynolds number for Darcy (or laminar) flow, turbulent flow results and non-Darcy analysis must be utilized. The effect of non-Darcy flow is a rate-dependent skin effect. That is, as the velocity within the wellbore increases, there is an increase in the pressure drop between the wellbore and the fracture.

Thus, the typical equation for flow in the reservoir is modified to account for the effects of non-Darcy flow using the Forchheimer parameter β as shown in equation (1) below:

Where:

P=pressure dP , . .

— = pressure drop in a direction x

dx

μ = viscosity

K = permeability

k_r = relative permeability

A = cross sectional area to flow

β = Forchheimer parameter

p = density

q = flow rate For high velocity flow occurring in the fractures and near the wellbore, non-Darcy flow results in a significant increase in the pressure drop and therefore plays an important role in properly modeling shale gas production. Because prior art techniques did not tend to model natural fracture pore systems, to the extent non-Darcy flow analysis has been used in the past for reservoir modeling, it has only been utilized to model flow in matrix pore systems and induced fracture systems.

Unfortunately, inclusion of the effect of equation (1) in reservoir fluid flow represents significantly more effort than was required for the skin factor. Since velocity depends not only on pressure drop but also on viscosity and relative permeability, a highly non-linear dependence is added to the flow equations for gridblock-to-gridblock or fracture-to-fracture non-Darcy flow treatment. The skin factor only requires a minor modification to the coefficient for the pressure loss between the wellbore and the reservoir or fractures. Inclusion of the non-Darcy effect adds a significant non-linear term to the pressure equations and requires that this term be included in the linearization for the Newton-Raphson iteration to solve for the flow in the wellbore and reservoir. In turn, this may increase the number of non-linear iterations and therefore increase overall

computation time for the reservoir simulation.

Gas desorption for shale development is an important, heretofore underutilized parameter in shale formation modeling. It is estimated that in some shale formations, more than 50% of the gas production will be due to desorption. To the extent gas desorption has been modeled in shale reservoirs, its use has been limited to desorption from the shale matrix. It has not heretofore been applied to desorption analysis from kerogen vugs.

Because economics are highly dependent on ultimate recovery from the formation, gas desorption must be treated for a shale gas reservoir simulation system to have any credibility. Moreover, it must be applied in a way that accounts for the existence of kerogen vugs within the reservoir. Desorption is described by the Langmuir equation (equation (2) below) for isothermal desorption characteristics:

Where:

Vg = volume of gas contained in the porous medium

VL = asymptotic adsorption volume

PL = pressure at which the adsorbed volume reaches VL

P = reservoir pressure

Use of equation (2) in the simulator results in a modification similar to that for dual porosity, single permeability ("DPSP"), in which a source of gas, i.e., vug nodes, are included in each grid, the volume of which depends on the change in matrix pressure over a timestep in the simulator.

For more rigorous treatment of the physics, consideration of sorption time and desorption effects on formation permeability may also be necessary. Sorption time is the time it takes for 63.2% of the gas to be desorbed as calculated using equation (2). In the case of shale gas, this time is generally extremely short and can be ignored. Similarly, the effect of desorption on matrix permeability is generally very small for shale gas and can also be easily ignored. In practice, when pressure is lowered in the horizontal production wellbore 302 (Fig. 3), the pressure is almost instantaneously lowered in all of the fracture system

(including both induced fracture and natural fractures) to which the wellbore is connected. For those fractures connected directly to a kerogen vug, the pressure is also reduced from the initial pressure. From equation (2) above, the kerogen must release gas into the surrounding reservoir fractures and matrix based on the Langmuir equation. Although the VL and PL parameters can be experimentally determined, often these are estimated by analogy and rules of thumb. Flow through the multiple fractures and matrix represents a significant difference with conventional treatment in which only the matrix contains the adsorbed gas which is directly in contact with fractures. The more complex treatment using vugs should allow a more realistic simulation of desorption than is currently achieved since the complex geometry of the porous medium is correctly characterized.

With reference to Figure 4, a flowchart is shown illustrating the steps of the process of the invention. The process is utilized to model flow characteristics to a wellbore of a shale reservoir having kerogen vugs and is preferably performed in conjunction with a three dimensional model of a reservoir. In step 400, reservoir characterization is undertake in which at least three, and preferably four different pore types are described based on fractured shale characteristics. In one embodiment, at least three different pore types are identified, selected from the group consisting of natural fracture pore systems, matrix pore systems, induced fracture pore systems and vug pore systems. In one embodiment, four different pore types are identified, namely natural fracture pore systems, matrix pore systems, induced fracture pore systems and vug pore systems. In any event, in step 402, the pore types are utilized to create one or more subgrids that represent a zone within the reservoir. Each zone includes a plurality of nodes of at least one of the pore types. In one embodiment, a subgrid for at least three different pore types is created for a zone. In one embodiment, a subgrid for each of the four pore types is created for a zone.

In step 404, once the shale reservoir has been described sufficiently, connectivities or transfer terms, if any, between the nodes are identified and assigned. This may include connectivity between similar nodes within the same subgrid, such as between matrix nodes within a subgrid, or may include connectivity between the nodes of one subgrid and the nodes of another associated subgrid, such as between vug nodes and natural fracture nodes or between matrix nodes and vug nodes. These transfer terms are the parameters that effect flow rates among the various porosity types, such as, for example, initial pore pressures, basis to the nodes of the subgrids. In one embodiment, the model consists of at least three different pore types and associated volumes which contain fluids which are to be modeled. In one embodiment, the model consists of at least four different pore types and associated volumes which contain fluids which are to be modeled. In step 406, know magnitudes for the transfer terms may be assigned, such as, for example, densities, volumes, flow rates and compressibilities.

In step 408, source terms are now incorporated as boundary conditions to the model in such a way that extraction of the gas is consistent with the wellbore's induced fractures. Put another way, to initiate flow in the simulation, a wellbore pressure is selected and incorporated into the model. This pressure affects the flow in the induced fractures, which, in turn by virtue of the transfer terms, affects flow between the other porosity types.

In step 410, a linear solver is utilized to solve for any unknown magnitudes of the transfer terms associated with the nodes. In one embodiment, non-linear equations are selected to model the reservoir and the subgrids and nodes thereof. In one embodiment, the linear solver methodology is applied by subgrid or associated subgrids. The Newton- Raphson method is then applied to linearize these non-linear equations. The linear solver then can be applied to the linear equations to solve for the unknowns. In one embodiment, this step may be iterated utilizing the resultant magnitudes until a desired degree of convergence is achieved between the linear and non-linear equations. In step 412, optionally, once a desired degree of convergence is obtained and the magnitudes of the unknowns are identified, time may be incremented and/or the wellbore parameters, such as the boundary conditions of pressure, may be altered to achieve a desired level of mass transfer and fluid flow for the modeled reservoir.

The foregoing methods and systems described herein are particularly useful in drilling wellbores in shale reservoirs. First a shale reservoir is modeled as described herein to design a well completion plan for a well. In an embodiment, the drilling well completion plan includes the selection of a fracturing plan, which may include the selection of fracture zones and their positioning, fracturing fluids, proppants and fracturing pressures. In other embodiments, the drilling well completion plan may include selecting a particular trajeocry of the wellbore or selecting a desired wellbore pressure to facilitate mass transfer and fluid flow to the wellbore. Based on the model, a drilling plan may be implemented and a wellbore drilled in accordance with the plan. Thereafter, in one embodiment, fracturing may be carried out in accordance with the model to enhance flow from the reservoir to the wellbore. In another embodiment, wellbore pressure may be adjusted in accordance with the model to achieve a desired degree of mass transfer and fluid flow. Those of ordinary skilled in the art will appreciate that while the method of the invention has been described statically as part of implementation of a drilling plan, the method can also be implemented dynamically. Thus, a drilling plan may be implemented and data from the drilling process, and in particular, the actual flow characteristics of the reservoir, may be used to update the model for the drilling of additional wellbores within the reservoir. After implementing the drilling plan, the system of the invention may be utilized during the drilling process on the fly or iteratively to calculate and re-calculate connectivity characteristics of the reservoir over a period of time as parameters change or are clarified or adjusted. In either case, the results of the dynamic calculations may be utilized to alter a previously implemented drilling plan. For example, the dynamic calculations may result in the utilization of a heavier or lighter fracturing fluids.

While certain features and embodiments of the invention have been described in detail herein, it will be readily understood that the invention encompasses all modifications and enhancements within the scope and spirit of the following claims. Furthermore, no limitations are intended in the details of construction or design herein shown, other than as described in the claims below. Moreover, those skilled in the art will appreciate that description of various components as being oriented vertically or horizontally are not intended as limitations, but are provided for the convenience of describing the invention.

It is therefore evident that the particular illustrative embodiments disclosed above may be altered or modified and all such variations are considered within the scope and spirit of the present invention. Also, the terms in the claims have their plain, ordinary meaning unless otherwise explicitly and clearly defined by the patentee.

Claims

What is claimed is: 1. A method for performing simulation of a shale reservoir represented by a model, the method comprising computer-implemented steps of:

characterizing at least three different porosity types for the modeled reservoir, the at least three different porosity types selected from the group consisting of natural fracture pore systems, matrix pore systems, induced fracture pore systems and vug pore systems; identifying transfer terms between the at least three porosity types; and

solving a linear system for the model using a linear solver.

2. The method of claim 1, wherein the formation is characterized with at least one of each of the porosity types of natural fracture pore systems, matrix pore systems, induced fracture pore systems and vug pore systems.

3. The method of claim 1 , wherein the characterized pore types are utilized to create one or more subgrids that represent a zone within the reservoir.

4. The method of claim 3, wherein each zone includes a plurality of nodes of at least one of the porosity types.

5. The method of claim 3, wherein a subgrid for at least three different porosity types is created for a zone.

6. The method of claim 3, wherein a subgrid for each of the four porosity types is created for a zone.

7. The method of claim 5, wherein each zone includes a plurality of nodes of at least one of the porosity types.

8. The method of claim 7, wherein transfer terms between the same node types within the same subgrid are identified.

9. The method of claim 7, wherein transfer terms between different node types in different subgrids are identified.

10. The method of claim 7, wherein transfer terms between the different node types within the same subgrid are identified.

1 1. The method of claim 7, wherein transfer terms between nodes are assigned on a nodal basis to the nodes of the subgrids.

12. The method of claim 11, wherein the transfer terms include initial pore pressures, fluid distributions and volumes.

13. The method of claim 11 , wherein the characterized pore types include at least one of each of the porosity types of natural fracture pore systems, matrix pore systems, induced fracture pore systems and vug pore systems, and wherein the zone includes at least one node of each porosity type.

14. The method of claim 1, wherein the step of solving the linear system comprises selecting non-linear equations to represent the modeled reservoir; and linearizing the nonlinear equations for subsequent solving using the linear solver.

15. The method of claim 1 , wherein the characterized pore types are utilized to create one or more subgrids that represent a zone within the reservoir, and wherein the linear solver is applied by subgrid or associated subgrids.

16. The method of claim 1, wherein the step of solving a linear system is iterated utilizing the resultant magnitudes until a desired degree of convergence is achieved between the linear and non-linear equations.

17. The method of 1 , further comprising the step of altering the wellbore pressure of the model to achieve a desired level of mass transfer and fluid flow for the modeled reservoir.

18. A computer program product comprising non-transitory computer-readable medium having stored thereon instructions executable by a computer for causing the computer perform simulation of a reservoir represented by a model having a plurality of porosity nodes, the instructions for causing the computer to:

characterizing at least three different porosity types for the modeled reservoir, the at least three different porosity types selected from the group consisting of natural fracture pore systems, matrix pore systems, induced fracture pore systems and vug pore systems; identifying transfer terms between the nodes of at least three porosity types; and solving a linear system for the model using a linear solver.

19. The computer program product of claim 18, wherein the characterized pore types include at least one of each of the porosity types of natural fracture pore systems, matrix pore systems, induced fracture pore systems and vug pore systems, wherein the characterized pore types are utilized to create one or more subgrids that represent a zone within the reservoir, and wherein the zone includes at least one node of each porosity type.

20. The computer program product of claim 19, wherein transfer terms between nodes are assigned on a nodal basis to the nodes of the subgrids.

21. The computer program product of claim 21 , wherein the linear solver is applied by subgrid or associated subgrids.

22. A method for drilling one or more wellbores in shale reservoir, which method comprises:

modeling an oil and gas shale reservoir having natural fracture pore systems, matrix pore systems, induced fracture pore systems and vug pore systems;

characterizing at least three different porosity types in the modeled reservoir, wherein the characterized pore types are utilized to create one or more subgrids that represent a zone within the reservoir, and wherein the zone includes at least one node of each porosity type;

assigning transfer terms between the at least three porosity types, wherein transfer terms between nodes are assigned on a nodal basis to the nodes of the subgrids; and solving a linear system for the model using a linear solver; preparing equipment to construct a portion of said wellbore; based on the modeled reservoir, selecting a characteristic for the wellbore; and

drilling a wellbore in accordance with the selected characteristic.

23. The method of claim 22, wherein the selected characteristic is the trajectory of the wellbore.

24. The method of claim 22, wherein the selected characteristic is the pressure of the wellbore.

25. The method of claim 24, further comprising the step of iteratively altering the wellbore pressure of the model to identify a wellbore pressure at which a desired level of mass transfer and fluid flow for the modeled reservoir is achieved; and utilizing the identified wellbore pressure as the selected characteristic.

26. The method of claim 22, further comprising the steps of drilling a first wellbore in the reservoir; recording values associated with mass transfer and fluid flow around the first wellbore; and utilizing the recorded vales as the values associated with a portion of the assigned transfer terms between the at least three porosity types; and drilling a second wellbore in the reservoir, wherein the second wellbore is the wellbore drilled in accordance with the selected characteristic.

27. The method of claim 22, wherein the characterized pore types include at least one of each of the porosity types of natural fracture pore systems, matrix pore systems, induced fracture pore systems and vug pore systems, wherein the characterized pore types are utilized to create one or more subgrids that represent a zone within the reservoir, and wherein the zone includes at least one node of each porosity type.

28. The method of claim 27, wherein transfer terms between nodes are assigned on a nodal basis to the nodes of the subgrids.

29. The method of claim 28, wherein the linear solver is applied by subgrid or associated subgrids.