WO2009079088A1 - Modeling subsurface processes on unstructured grid - Google Patents

Abstract

Embodiments of the invention involve forming a prismatic grid and solving a convection-diffusion problem using the prismatic grid and mixed finite element analysis. The prismatic grid may be formed by providing a triangular mesh on a plane of a model. The mesh is then coarsened to make cells that are less desirable larger. The coarsened grid is then projected to form the prismatic grid. Each cell of the grid is then assigned a plurality of degrees of freedom. Mixed finite element analysis of the grid produces a matrix, which is then solved to yield a solution to the convention-diffusion problem.

Description

MODELING SUBSURFACE PROCESSES ON UNSTRUCTURED GRID

CROSS-REFERENCE TO RELATED APPLICATION

[0001] This application claims the benefit of U.S. Provisional Patent Application 61/007,761 filed December 14, 2007 entitled MODELING-SUBSURFACE PROCESSES ON UNSTRUCTURED GRID, the entirety of which is incorporated by reference herein.

TECHNICAL FIELD

[0002] This application relates in general to computer modeling, and in specific to modeling subsurface processes with an unstructured grid.

BACKGROUND OF THE INVENTION

[0003] In geological exploration it is desirable to obtain information regarding the various formations and structures that exist beneath the Earth's surface. Such information may include determine geological strata, density, porosity, composition, etc. This information is then used to model the subsurface basin using the obtained data to predict the location of hydrocarbon reserves and aid in the extraction of hydrocarbon.

[0004] Unstructured grids have many appealing characteristics for modeling physical processes in complex geologic structures such as subsurface basins. Such grids may also be used in other industries, for example in the airspace industry or the auto industry. The basin or domain of interest may be modeled or represented as a set of layers of different thickness stacked together. The geological layers may be fractured along vertical or slanted surfaces and degenerate creating so-called pinch-outs. Pinch-outs are defined as parts of geological layers with near zero thickness. This complexity should be taken into account by the grid to produce a good model of the geological layers. An unstructured grid provides a better model than a structured grid. An unstructured grid may comprise a set of polyhedral elements or cells defined by their vertices and have a completely arbitrary topology. For example, a vertex of the grid can belong to a number of cells and each cell can have any number of edges or faces.

[0005] Many physical subsurface processes may be described by mathematical equations of convection-diffusion type. Examples of such processes can be fluid flow in porous media, temperature distribution, and/or pressure distribution. An important process for oil exploration is the temperature distribution or thermal modeling. Thermal modeling involves the heat moving from the magma below the crust and through the sedimentary layers and source rock. Source rocks are rocks that are involved in the formation of oil and other hydrocarbon materials. The oil and/or other hydrocarbon materials would be expelled from the source rocks and migrate elsewhere. The quality of hydrocarbon is determined by the temperature and pressure conditions inflicted on the source rocks and their surrounding area. The quality is also affected by the temperature and pressure conditions of the migration path between the source rocks and its current location. Thus, the pressure and temperature conditions of the basis throughout its history is important.

[0006] To more accurately model the processes, it is important to model not only the primary variables, such as pressure or temperature, but also model their fluxes, or the rates of flow of energy, fluids, etc. over any given surface There is a variety of known approaches for modeling these processes, such as finite difference, finite volume, or finite element methods. In these approaches, where a physical process is considered, the domain is covered by a grid. Then, the domain is approximated on a grid by introducing a set of unknowns called the degrees of freedom at specified locations of the grid cells and deriving algebraic equations for each location that connect the degree of freedom in that location with other degrees of freedom. The way of deriving such equations, as well as the locations of degrees of freedom, is different for different approaches mentioned above, but all these methods have a common feature, namely, that they only involve primary variables, such as temperature or pressure.

[0007] To compute fluxes, an interested person would first compute the desired primary variable using one of the above-described approaches, and then use numerical differentiation to compute the flux of the primary variable. All existing methods of numerical differentiation being accurate on regular grids, e.g. rectangular or parallelepiped grids, are inaccurate and very computationally expensive on unstructured grids, especially if the domain where the physical process is considered is highly heterogeneous. Moreover, the approaches for solving convection-diffusion problems using finite difference methods require Cartesian grids, and thus are not applicable in many subsurface applications, which have to employ unstructured grids. The finite element methods, being able to model complex geometries do not have local conservation property and can not be applied in many subsurface processes. Conversely, finite volume approaches are locally conservative and can be applied on a subset of unstructured grids which are locally orthogonal. However, when the unstructured grid does not posses local orthogonality property, finite volume method provides inaccurate solution. Thus, from all three classes of the approaches mentioned above, none is applicable for description of subsurface convection-diffusion processes in a basis modeled with an unstructured grid.

[0008] There is another mathematical approach to simultaneously approximate primary unknowns and their fluxes, called mixed finite element method, which is described in F. Brezzi and M. Fortin, "Mixed and hybrid finite element methods", Springer Verlag, Berlin 1991. Such method is proven to be locally mass conservative, accurate in the presence of heterogeneous medium, and provide accurate approximations to both, primary unknowns and fluxes. Until recently, the mixed finite element methods could not be directly applied to the domains covered by unstructured polyhedral grids, which are very common for the subsurface applications. A new version of mixed finite element method for diffusion-type equations on arbitrary polyhedral grids is proposed in Yu. Kuznetsov and S. Repin, "New mixed finite element method on polygonal and polyhedral meshes" , Russian Journal of Numerical Analysis and Mathematical Modeling, v.18, pp. 261-278, 2003.

BRIEF SUMMARY OF THE INVENTION

[0009] The present invention is directed to systems and methods which provide an accurate model of energy transfer and/or pressure distribution in basins being developed through geological times. At any given time, a basin is represented as a set of layers of different thicknesses stacked together. In some locations in the basin, the thickness of a layer degenerates to zero, forming a pinch-out. Embodiments of the invention use a prismatic mesh and mixed finite element analysis to model various processes in the basin, including energy transfer, e.g. thermal energy, and pressure. Thus, embodiments of the invention solve for both primary unknowns, e.g. temperature or pressure and secondary unknowns, e.g. temperature flux or pressure flux. One or more of the following aspects may be used to provide an accurate model of energy transfer and/or pressure distribution, e.g., a physical process, in basins being developed through geologic time. The model may be used to interpret a modern day reservoir, and in turn may be relied upon to control hydrocarbon production activities based on simulated results of the model. The production of hydrocarbons may be controlled, e.g., production rates from surface facilities may be controlled, wells may be strategically placed, and/or a reservoir generally characterized based on results interpreted from simulated basin model(s) generated by one or more of the following aspects. [0010] In one general aspect, a method for modeling on a computer a physical region, wherein the physical region includes a plurality of strata, the method includes receiving data that defines at least one physical characteristic of the physical region; providing a triangular mesh on a plane of a model of the physical region, wherein the mesh comprises a plurality of cells; coarsening the triangular mesh in a non-uniform manner such that cells that are less desirable are larger; and projecting the coarsened triangular mesh in a direction orthogonal to the plane in the physical region to form a prismatic grid, wherein each of the cells of the coarsened triangular mesh is separated into sub-cells according to the strata.

[0011] Implementations of this aspect may include one or more of the following features. For example, coarsening may include projecting the data onto a plane; and using the data to determine which cells are less desirable. The model may include modeled features that model physical features in the physical region, and wherein using may include assigning a priority value to each cell, wherein the value is determined based on whether each cell is proximate to a modeled feature and a type of the modeled feature. Providing a triangular mesh may include providing a rectangular mesh on the plane; and splitting each cell of the rectangular mesh along at least one diagonal. Coarsening may include merging two adjacent triangles by eliminating a side common to the two adjacent triangles. The prismatic grid may include a plurality of prism cells, a plurality of pyramid cells, and a plurality of tetrahedron cells. The method may be used to model at least one flux of a physical process in the physical region, the method further including assigning a plurality of degrees of freedom for the flux in each sub-cell; applying mixed finite element analysis to each of the sub-cells to produce a matrix; and solving the matrix to determine the flux of the physical process in the region.

[0012] Assigning may include for each cell assigning one degree of freedom for the physical process; and assigning another degree of freedom for each face of the cell. Applying may include using a div-constant approach to form the finite element space. The physical process may be a convection-diffusion process. The physical process may be one of temperature and pressure and the physical region is a subsurface geological basin. The physical process may involve the formation of hydrocarbon material. The physical process may involve the movement of hydrocarbon material. The data may be derived from information from a sensor that measured the at least one physical characteristic of the physical region. [0013] In another general aspect, a method for modeling a physical process and a flux of the physical process on a computer includes forming an unstructured, prismatic grid that models a physical region, wherein the physical process operates within the physical region and the prismatic grid comprises a plurality of cells; assigning a plurality of degrees of freedom for the physical process and the flux for each cell; applying mixed finite element analysis to each of the cells to produce a matrix; and solving the matrix to determine the physical process and the flux in the region.

[0014] Implementations of this aspect may include one or more of the following features. For example, forming may include providing a triangular mesh on a plane of a model of the physical region, wherein the mesh comprises a plurality of cells; coarsening the triangular mesh in a non-uniform manner such that cells that are less desirable are larger; and projecting the coarsened triangular mesh in a direction orthogonal to the plane in the physical region to form the prismatic grid. The prismatic grid may include a plurality of prism cells, a plurality of pyramid cells, and a plurality of tetrahedron cells. Assigning may include assigning one degree of freedom for the physical process for each cell; and assigning another degree of freedom for each face of the cell for each cell. Applying may include using a div-constant approach to form the finite element space. The determined physical process and flux may be used to affect a change in the physical region. The physical process may be one of temperature and pressure and the physical region is a subsurface geological basin.

[0015] In another general aspect, a computer program product having a computer readable medium having computer program logic recorded thereon for modeling at a physical process and a flux of the physical process in a physical region, the computer program product including code for forming an unstructured, prismatic grid that models the physical region; code for applying mixed finite element analysis to the prismatic grid to produce a matrix; and code for solving the matrix thereby determining the physical process and the flux in the region.

[0016] Implementations of this aspect may include one or more of the following features. For example, code for forming may include code for providing a triangular mesh on a plane of a model of the physical region, wherein the mesh comprises a plurality of cells; coarsening the triangular mesh in a non-uniform manner such that cells that are less desirable are larger; and projecting the coarsened triangular mesh in a direction orthogonal to the plane in the physical region to form the prismatic grid. The prismatic grid may include a plurality of cells, and the code for applying may include assigning one degree of freedom for the physical process to each cell; assigning another degree of freedom for each face of the cell to each cell; and using a div-constant approach to form the finite element space. The code for solving may include using preconditioned conjugate gradient analysis to solve the matrix.

[0017] Embodiments of the invention operate by projecting some or most geological and geometrical features, such as pinch-out boundaries into horizontal plane. Note that projection can be non-orthogonal or slanted. Embodiments of the invention then create an unstructured grid resolving all the desired features on that plane. Note that the grid can be comprised of polygons, quadrilaterals, triangles, or combinations thereof. Embodiments of the invention then project the obtained grid back onto all boundary surfaces of all layers, thereby constructing a prismatic grid. The prismatic grid may comprise a plurality of cells, which can be prisms, tetrahedral shapes, pyramids, or combinations thereof. Note that the unstructured prismatic grid approximates boundary surfaces of all layers.

[0018] Embodiments of the invention may then operate by associating one degree of freedom per cell at the cell center for primary unknown and one degree of freedom per each face of the cells at the face center for normal components of flux. Embodiments of the invention then discretize the problem using a mixed finite element approach, for example the approach of Yu. Kuznetsov and S. Repin. The spatial discretization produces a sparse matrix equation. Embodiments of the invention may then solve the matrix equation to get both, primary unknowns and normal components of the flux at the faces of the cells. Thus, embodiments of the invention provide more accurate modeling without greatly expanding the number of unknowns that are required to be solved.

[0019] Embodiments of the invention may form the prismatic grid by providing a triangular mesh that covers a horizontal plane of the physical region. The mesh may then be coarsened in a non-uniform manner such that cells that are less desirable are larger, while leaving the desirable cells with a more fine format. The coarsened mesh is then projected in a vertical direction in the physical region to form the prismatic mesh

[0020] The foregoing has outlined rather broadly the features and technical advantages of the present invention in order that the detailed description of the invention that follows may be better understood. Additional features and advantages of the invention will be described hereinafter which form the subject of the claims of the invention. It should be appreciated by those skilled in the art that the conception and specific embodiment disclosed may be readily utilized as a basis for modifying or designing other structures for carrying out the same purposes of the present invention. It should also be realized by those skilled in the art that such equivalent constructions do not depart from the spirit and scope of the invention as set forth in the appended claims. The novel features which are believed to be characteristic of the invention, both as to its organization and method of operation, together with further objects and advantages will be better understood from the following description when considered in connection with the accompanying FIGURES. It is to be expressly understood, however, that each of the FIGURES is provided for the purpose of illustration and description only and is not intended as a definition of the limits of the present invention.

BRIEF DESCRIPTION OF THE DRAWINGS

[0021] For a more complete understanding of the present invention, reference is now made to the following descriptions taken in conjunction with the accompanying drawing, in which:

[0022] FIGURE 1 depicts a domain being partitioned into layers, according to embodiments of the invention;

[0023] FIGURES 2A and 2B depict a domain and the domain covered with a rectangular mesh, according to embodiments of the invention;

[0024] FIGURE 3 depicts a non-uniformly coarsened triangular grid, according to embodiments of the invention;

[0025] FIGURE 4 depicts different types of cells formed according to embodiments of the invention;

[0026] FIGURE 5 depicts a 3D prismatic grid formed according to embodiments of the invention;

[0027] FIGURE 6 depicts a tetrahedral cell used by embodiments of the invention;

[0028] FIGURE 7 depicts a pyramidal cell used by embodiments of the invention;

[0029] FIGURES 8A-8D depict a prismatic cell used by embodiments of the invention, and the separated into three tetrahedrons, according to embodiments of the invention; [0030] FIGURE 9 depicts independent face splitting of neighboring cells, according to embodiments of the invention;

[0031] FIGURE 10 depicts a method of forming a prismatic grid, according to embodiments of the invention;

[0032] FIGURE 11 depicts a method of solving a matrix, according to embodiments of the invention; and

[0033] FIGURE 12 depicts a block diagram of a computer system which is adapted to use the present invention.

DETAILED DESCRIPTION OF THE INVENTION

[0034] Embodiments of the invention are useful for modeling subsurface oil fields. The examples of the embodiments described herein may reference such oil fields. However, the embodiments may be used to model other domains involving other materials and/or processes. For example, other hydrocarbon materials may be involved, such as coal. Embodiments of the invention may be useful for mining or tunneling. Embodiments of the invention may be used for other domain types, e.g. the atmosphere, and may be useful for modeling the weather, temperature, and/or pollution. Another domain may be the oceans, and embodiments may be used to measure sound, temperature, saliently, and/or pollution. Any type of stratified domain may be modeled using embodiments of the invention. Any type of material that moves through a convection-diffusion process may be modeled using embodiments of the invention. Any type of flux that is present in the domain or material may be modeled using embodiments of the invention.

[0035] As stated earlier, embodiments of the invention can be applied to any convection- diffusion process. The following is an example of a 3D convection-diffusion type equation

- V - (KVp) + c - p = f in Ω (1.1)

where p is an unknown function (called as pressure), K=K(x) is a diffusion tensor, c is a nonnegative function, /is a source function, and Ω c i?³ is a bounded computational domain. It is assumed that K is a uniformly positive definite matrix and the boundary δΩ of the domain Ω is partitioned into two non-overlapping sets T_D and T^. [0036] Equation (1.1) is complemented with the boundary conditions

P = g_D on τ_{D (i 2)}

(KVp)-n + σ-p = g_N on T_N

where n is the outward unit normal vector to T_N, σ is a nonnegative function, and g_# and g_N are given functions. It is assumed that equations (1.1)-(1.2) have a unique solution.

[0037] The partial differential equations (1.1)-(1.2) may be replaced by the equivalent first order system:

u + KVp = 0 in Ω

(1.3) V-u + c ^■ p = f in Ω

P = g_D on r, _{(i 4)}

-u-n + σ-p = g_N on T_N

[0038] Equations (1.3)-(1.4) are the mixed formulation of equations (1.1)-(1.2). Note that in this way the primary unknown/? and its flux u may be approximated simultaneously.

[0039] As stated earlier, embodiments of the invention may operate with different domains. Thus, let G be a domain in R² with a regularly shaped boundary dG, e.g. piecewise smooth and angles between pieces are greater than 0. Let computational domain Ω be defined as follows

Ω = {(x,y,z)<ΞR :(x,y)≡G,Z_am(x,y)≤z≤Z_aia(x,y

where Z_m1n(X₅J;) and Z_max(x,j) are smooth surfaces.

[0040] Let mbea positive integer and z = Z₁(Xy), i = 0,..., m be single-valued continuous functions defined on G such that

Z₀(x, y) ≡ Z_^(x,y) in G

Z₁ __j (x, y) < Z₁ (x, y) in G i = \,m zΛ*,y) ^≡ z_^(^χ, y) m G

[0041] These functions define the interfaces between geological layers. In other words, the computational domain Ω can be split into m sub-domains (strips or layers) which are defined as follows for all i=\,...,m. Ω, = {(x, y, z) e Ω : (x, y) e G, Z_M (x, y) < z < Z_t (x, y)}.

[0042] FIGURE 1 depicts an example of the partitioning of computational domain Ω 100 into a plurality of sub-domains or layers 101-107. Note that FIGURE 1 depicts the different layers in an exploded view, however for computation, the layer need not be separated. Further note that it is assumed that sub-domains Ω_z, satisfy the cone condition, where the boundaries of the sub-domains do not have singular points (zero angles, etc) and, in addition, that all the sets

P₁ = ^y) G G : Z_i__l(x,y) = Z_i(x,y)}

comprise of a finite number of polygons. The boundary of corresponding set P₁ is denoted by 3P₁.

[0043] FIGURE 1 depicts the different strata of a basin. The data used to form the different layers of FIGURE 1 may be determined by various techniques, such as stratigraphic analysis and/or seismic inversion, using sensors to measure various characteristics of the basin.

[0044] The boundaries of the sets P₁ may be projected to the flat plane in the following manner. For any given point (x,y,z) from 3P₁, the projected point has coordinates (x,y,0). All such points organize the set of closed line like those of FIGURE 2A, which are used to create plane triangulation, examples of which are shown in FIGURE 2B and FIGURE 3.

[0045] The variational mixed formulation of differential equations (1.3)-(1.4) can be written as follows: find u e H_Λv(Ω.) , p e L₂(Q) , and λ e L₂ (T_N ) such that

for all , where

[0046] It is noted that λ is the restriction of the pressure function p = p(x) onto T_N. In this formulation, all the boundary conditions are natural.

[0047] In case σ = 0 on T_N the variational mixed formulation can be written in the different form as follows: find u e H_Λv (Ω) , u • n = -g_N on T_N and p e L₂(Q) such that

for all v e H _dlv (Ω) , v • n = 0 on T_N and q ^ L₂ (Ω) .

[0048] In the following analysis, the equation (1.5) is considered, although the conclusions can be also be applied to equation (1.6) without loss of generality .

[0049] Embodiments of the invention use prismatic grids, which provide many appealing characteristics in modeling convection-diffusion subsurface processes. In many cases, a domain can be represented as a set of layers of different thickness stacked together. Slice- wise geological structure and unstructured geometrical features in stacked layers are represented by prismatic grids satisfactorily. The 2D geometrical data is provided by post processing of geostatistical information. Usually, the data is material properties associated with the nodes or cells of a 2D fine rectangular grid. Note that there may be millions of nodes. The presence of material data in a node implies a computational node, whereas their absence implies a node-outlier. The set of the computational nodes defines the computational domain.

[0050] The interfaces between geological layers are geostatistical data and may intersect each other resulting in topologically incorrect situations. The bottom interface is the bottom boundary of the lowest geological layer and is represented by geostatistical data as well. The geological layers may be fractured along vertical surfaces and degenerate. Pinch-outs are defined as parts of geological layers with thickness modulo not greater than a user defined threshold δ ≥ 0. Fault polylines are defined as intersections of the bottom geological interfaces and faults.

[0051] To simplify the description of the algorithms, the grids to be used in embodiments of the invention should satisfy some very natural requirements. The objective prismatic grid should be a logical product of a ID grid and a 2D triangulation, whose nodes form a subset of the geostatistical rectangular grid with no nodes-outliers, with the node density equal to that in the original grid in user-specified regions of interest. Besides, the triangulation must be refined in the vicinity of user-defined faults and wells, as well as automatically detected pinch-outs. Also, the number of triangles in the 2D grid should not be greater than a user- defined number. Regarding the prismatic grid, lateral faces of prisms have to approximate the geological interfaces and form shape regular 2D triangulations.

[0052] The process described below is an example of a process that may be used to construct prismatic grid. Note that other processes may be used. Furthermore, the mixed finite element method is not bound by that type of prismatic grid. First, is the generation of 2D regular triangulation refined towards projections of pinch-outs onto the bottom geological interface, fault polylines, and points representing wells. Next, the 2D triangulation is projected onto surfaces defined by functions Z_}(x,y), i=\,..., m, to form the resulting 3D prismatic grid. This process is further described in the following paragraphs.

[0053] To form the 2D regular triangular grid, the an exemplary process may begin with a rectangular grid. Given coordinates of the nodes in x- and j-direction, a rectangular conforming mesh G^h covering domain G is generated, which is composed of cells having at least one of the four nodes with material data. For example, FIGURE 2A depicts a domain 200, and FIGURE 2B depicts a rectangular grid 201 covering the domain 200. Without loss of generality, it is assumed that the mesh G^h is composed of squares, i.e. mesh sizes in x- and j-direction are equal, h_x=h_y=h.

[0054] Next, according to an embodiment, each rectangular cell is split by its diagonal into two triangles. One processes that may be used to form the triangles is described below, note that other processes may be used. The choice between two possible diagonals may be made according to the following rule. Let each rectangular cell be assigned an integer equal to the sum of minimal x- and j-indices of its nodes. For the cells with even numbers, the splitting diagonal has the node of the cell with minimal x- and j-indices. For the cells with odd number, the other diagonal is chosen. The above process specifies the triangulation uniquely for a given set of nodes with a given material data distribution. Alternating the directions of the diagonals reduces issues of grid orientation. As an alternative process to form triangles, each rectangular cell may be split into four triangles by using both diagonals. [0055] The generated triangulation is projected onto the bottom geological interface as described in paragraph [0036]. Note that there can be defined regions of interest ω₁ in the domain G, where the modification of the grid is not necessary or not desired.

[0056] Let P,^h denote the maximal subset of rectangular elements of G^h, which belong to P₁. If there is no element of G^h which belongs to P₁ but there is the vertex of G^h which belongs to P₁ then this vertex is said to belong to P* . Then the set is defined

which is known as a "pinch-out" projection. Based on that definition, "pinch-out" projections are the subset of edges and vertices belonging to G^h.

[0057] Next, according to an embodiment, different priorities are assigned to the triangles. A priority, or an integer marker, is assigned to each triangle of the fine grid. Values of the priority control coarsening process.

[0058] At the beginning, zero priority is assigned to all the triangles. For triangles whose closures intersect faults, their priorities are changed to 1. To find triangles intersecting faults, the following method may be used. First, the triangles are extracted from the fault triangulations, which intersect the most bottom geological interface. Second, each extracted triangle is checked for the intersection with the fine grid triangles. For triangles whose closures intersect pinch-outs, their priorities are changed to 2. These triangles are defined as those where the following condition is violated: in all triangle nodes, the thickness of a geological layer is either greater than δ or less than -δ. For triangles whose closures contain well points, their priorities are changed to 3. For triangles that belongs to a user defined region of interest ω_z, their priorities are changed to 4. For triangles that satisfy several conditions, their priorities are changed to the maximum priority value.

[0059] After the priorities have been assigned, the triangular grid may be non-uniformly coarsened. The fine portion of the grid may have a large number of equally small triangles. These areas are more desirable because they contain more information, include interesting geologic features, e.g. wells, faults, pinch-outs, and/or are indicated as desirable by a user. The coarsened portions of the grid are not as desirable as the finer portions of the grid. The grid may have a range of coarsening, where the most coarsened indicates the portions having little or no desirable qualities, and the areas with no coarsening indicates the most desirable areas. Coarsened areas between the most coarsened and no coarsened indicate area with some desirable aspects.

[0060] Coarsening is a sequence of triangle-merging procedures. For example, two triangles may be coupled into one by elimination of their common side. This procedure comprises two stages. First, certain triangles are marked for coarsening. Second, they are coarsened. It should be noted that the grid conformity may cause coarsening of unmarked triangles. Each coarse triangle inherits the maximal priority of the two merging triangles. In addition to the priority, each triangle is assigned another integer denoted as level. Any triangle of the initial fine grid has level 1. Coarsening may be applied to a pair of triangles of the same level j, and result in a coarser triangle with level j+1.

[0061] Below is one example of a coarsening procedure according to embodiments of the invention. Note that other procedures may be used.

[0062] The coarsening procedure can be described as the loop:

1) Set k=\.

2) Form the set M from triangles with zero priority, which coarsening will not cause coarsening other triangles with non-zero priority.

3) IfM is empty then go to 6.

4) Coarsen triangles from M.

5) Go to 2.

6) If the number of triangles in the new grid is not greater than the user defined threshold Nt_usr, then Stop.

7) Form the set M from triangles with non-zero priority not greater than k, which coarsening will not cause coarsening other triangles with priority greater than k.

8) Coarsen triangles from M.

9) If k<3 , then set k=k+ 1 , otherwise set A= 1.

10) Go to l. [0063] The output triangular grid G_h has nodes coinciding with certain nodes of the projected input rectangular mesh and fine triangles in the region of interest, as well as triangles refining towards well points, fault polylines, and pinch-outs. FIGURE 3 depicts an example of the above grid coarsening procedure. The resulting grid 300 depicts non-uniform areas of most coarse triangles 301, coarsened triangles 302, and fine triangles 303. Note that the fine triangles may have some coarsening or no coarsening. Note that in FIGURE 3 the triangles have been formed by using the two diagonal method described above.

[0064] After coarsening, the 3D prismatic grid may be formed. Let βh be a triangle in the triangulation , £=1,2,3, be the vertices of en. Consider three

vertical lines (x,y) = d^k\ £=1,2,3, in R³ and denote by d^k'^h ', £=1,2,3, their intersections with the surfaces z = Z₁(X, y), i=0,...,m. Then, any polyhedron with the vertices located on the neighboring surfaces is either a vertical prism (all six vertices are different), or a pyramid (two vertices coincide, e.g. the corresponding vertex

or a tetrahedron (two pairs of the vertices coincide, i.e. two vertices of the set a®, £=1,2,3, belong to P¹¹). FIGURE 4 depicts a portion of a 3D prismatic grid with the three types of prisms, namely a vertical prism 401, a pyramid 402, and a tetrahedron 403. In other words, a pyramid is a prism with one edge disappearing, and a tetrahedron is a prism with two edges disappearing.

[0065] Performing the above operation for all the triangles in G_h will provide the partitioning Q_f1 of the domain Ω. In particular cell, each surface z = Z_t(x,y) is approximated by piecewise triangular surface z =

which comprises of the top (bottom) triangles of the prisms as well as the particular faces the pyramids and the tetrahedrons. FIGURE 5 depicts an example of a 3D prismatic grid 500.

[0066] After completion of the 3D prismatic grid, the grid may be subject to mixed finite element analysis. In the previous section it was indicated that the grid Qf₁ comprises elements

[Ek) which are either vertical prisms, or pyramids, or tetrahedrons. To formulate the mixed finite element (MFE) method for equation (1.5), the finite element subspaces of the spaces

H_Λv(Ω.) , L₂(Q), and L₂(T_N) should be defined.

[0067] The finite element space Lf₁ a L₂(Q) comprises functions pu which are constants on each grid cell Ek a Qf₁. The finite element space A_h a L₂(T_N) comprises functions Xf₁ which are constants on each intersection of a grid cell E_k in Ω_/, with boundary part T^. These intersections may be either quadrilaterals or triangles.

[0068] One problem in mixed finite element methods is the design of finite element subspaces Vh of the space H_Λv(Ω.) . For computational efficiency only those finite element vector-functions should be considered which have constant normal components on the interfaces Tu between neighboring cells Ek and Ei, k>l, as well as on the intersections T₁^ of a cell Ek with boundary T^. The dimension of the finite element subspace of the space H_ώv(Ω) is equal to the total number of different interfaces {T_M} and { f_w }. This finite element space may be constructed based on "div-constant" approach described in Yu. Kuznetsov and S. Repin, "New mixed finite element method on polygonal and polyhedral meshes", Russian Journal of Numerical Analysis and Mathematical Modeling, v.18, pp. 261- 278, 2003.

[0069] For a tetrahedral cell, T, the finite element space V^I_T coincides with a classical lowest order Raviart-Thomas finite element space RTo(T) (see F. Brezzi and M. Fortin, "Mixed and hybrid finite element methods", Springer Verlag, Berlin 1991). A finite element vector- valued function w_/, e R T₀(T) has four degrees of freedom (DOF), i.e.

where φ_z(x) are the basis vector- functions associated with the faces γ_z of the tetrahedron T, /=1,2,3,4.

[0070] Denote by γ₇ the face of tetrahedron T opposite to the vertex A₁, namely, γi is the face A2A3A4, γ₂ is the face A3A4A1, 73 is the face A4A1A2, and γ₄ is the face A1A2A3. Let n_} be the outward unit vector on the face γ₇, and h be the length of the perpendicular from the vertex A₁ onto the face γ_;,_/=l, 2,3,4. Such a tetrahedral cell 600 is shown in FIGURE 6.

[0071] The space of the lowest order Raviart-Thomas elements on a tetrahedron T can be defined as

Rr₀(T) = spanfø, φ₂, φ₃, φ_Λ), where the basis vector functions φ_z satisfy the conditions is the

Kronecker symbol, which is equal to 1 if / =j, and 0 otherwise.

[0072] Straightforward calculations show that the basis functions can be defined explicitly by

where x® are the coordinates of the vertices A₁, /=1,2,3,4.

[0073] For a Pyramidal cell, P, when a cell P e Ω_/, is a quadrilateral pyramid the finite element space Vh\? is constructed using the "div-constant" approach.

[0074] The faces of pyramid P are denoted by γ₇,y=l,2,3,4,5, namely, γi is the face A₂A₃A₅, γ₂ is the face A\A₃A/\, γ₃ is the face A1A2A4A5, γ₄ is the face A1A2A3, and 75 is the face A3A4A5. Such a pyramid cell 700 is shown in FIGURE 7.

[0075] To describe the "div-const" approach, the pyramid P is partitioned into two tetrahedrons Ti and T₂. It can be done in two different ways and both ways provide workable algorithm described below. Thus, without loss of generality, it is assumed that the pyramid P is divided into two tetrahedrons Ti = A1A2A3A4 and T₂ = A2A3A4A5. Let H₁ denote the unit outer normal vectors to the faces J₁, /=1,2,3,4,5, and by «6 the unit normal vector to the interface y₆ = A2A1A4 between tetrahedrons Ti and T₂ directed from Ti to T₂.

[0076] On each tetrahedron, the classical lowest order Raviart-Thomas space of vector- functions is construct, the sets of basis functions are determined, and vector

field W_h in P is defined as follows:

[0077] This representation shows that the vector-function \Xh is linear on each of the tetrahedrons Ti and T₂, belongs to the space Hώ_v(P), and satisfies the required conditions

To determine the unknown value M₆ of the normal component of flux on the interface γ₆ the following condition should be operative

V- W_h ≡ const in P. (2.2)

[0078] From definition (2.1) of u^ the expressions for divergence of u^ in Ti and T₂ are obtained. Application of the divergence operator for both sides of (2.1) yields

and

[0079] Since it follows that

[0080] The value of M₆ under assumption (2.2) can be found in a different way. First, the Stokes' formula is applied

and the value of V • U_f1 is determined by

where |γ_z| is the area of corresponding face γ_z and |P| is the volume of the pyramid. Second, the value of u₆ is determined from (2.4):

[0081] At this point the basis functions

on pyramid P are reconstructed which satisfy the conditions

[0082] Since the values of normal components of the basis functions on the faces of pyramid P are known, the values of its normal components , i=\, 2,3,4,5, on the internal face γ₆ can be found from the formulas (2.5) and (2.6). Therefore, the explicit representation of basis functions is given by

[0083] For a prismatic cell with a triangular base, let Il e Ω_/,. The finite element space V^n is constructed using the "div-const" approach.

[0084] Denote by γ₇ the faces of prism IT, namely, γi is the face A₂A₃A₅A₆, y₂ is the face AiA₃A₄A₆, γ₃ is the face AiA₂A₄A₅ , γ₄ is the bottom face AiA₂A₃, and γ₅ is the top face A₄A₅A₆. Such a prism cell 800 is shown in FIGURE 8A.

[0085] In order to apply the "div-const" approach according to an embodiment, prism Il is partitioned into three tetrahedrons T₁, T₂, and T₃. The partitioned tetrahedrons 801, 802, 803 are shown in FIGURES 8B-8D. Without loss of generality, it can be assumed that prism FI is divided into tetrahedrons Ti = AiA₂A₃A₆, T₂ = AiA₄A₅A₆, and T₃ = AiA₂A₅A₆. The partitioned tetrahedrons 801, 802, 803 are shown in FIGURES 8B-8D.

[0086] Let H₁, /=1,2,3,4,5, denote the unit outward normal vectors to faces J₁, n₆ denote the unit normal vector to the interface y₆ = AiA₂A₆ between tetrahedrons Ti and T₃ directed from Ti to T₃, and, finally, by nη the unit normal vector to the interface 7₇ = AiA₅A₆ between tetrahedrons T₂ and T₃ directed from T₂ to T₃. On each tetrahedron the classical lowest order Raviart-Thomas finite element space is constructed, the sets of basis functions

A= 1,2,3, are determined, and vector field u^ in Il is defined as follows

[0087] With this representation, the vector- function \Xh is linear on each of the tetrahedrons T₁, T₂, and T3, belongs to the space H^v(FI), and satisfies the required conditions

[0088] The natural choice for the unknown normal components Ue and uη on subsidiary interfaces y₆ and 77, respectively, comes from the condition where the divergence of u^ is a constant on IT, i.e.

V- W_h ≡ const in II. (2.8)

[0089] The unknown values Ue and uη are obtained in a similar way as discussed with respect to pyramidal cells. First, the value of V-u* on Tl is determined using Stokes' formula

therefore

where |γ_z| is the area of corresponding face γ_z and |II| is the volume of the prism. Second, the values of uβ and uη are computed as

and

[0090] At this point, the basis functions

on prism Il are reconstructed which satisfy the following conditions:

[0091] Since the values of normal components of the basis functions on the faces of the prism Il are known, the values of normal components

and

, z-1,2,3,4,5, on the internal faces y₆ and 77 can be determined, respectively, by formulas (2.10) and (2.11). Therefore, the explicit representation of basis functions is given by

[0092] It is reasonable to emphasize the following distinction of the "div-const" approach. The partition of a grid cell into tetrahedrons is independent from the partition of its neighboring cells, e.g. a quadrilateral face which is a common for two neighboring cells can be divided in different ways from different sides. For example, FIGURE 9 depicts two neighboring cells 901, 902 that have adjoining surfaces 903, 904 split in different manners. This feature simplifies generation of prismatic grids and discretization of boundary value problems on them comparing with tetrahedral grids.

[0093] For general prismatic cells, the approach described above for a prismatic cell with a triangular base may be extended to cover a generic prism II. Every prism may be independently partitioned into tetrahedral and the finite element space

is constructed using the "div-const" approach described in the previous sections.

[0094] The above paragraphs describe having one degree of freedom per each face of the cells at the face center for normal components of flux. Each cell is also assigned one degree of freedom associated at the cell center for each primary unknown. For example, primary unknown(s) may include ambient temperature and/or fluid pressure.

[0095] The following paragraphs describe using a hybrid version of mixed finite element (MFE) analysis on the cells of the 3D prismatic grid with the introduction of an additional degrees of freedom associated with the faces of the prisms. These degrees of freedom allow the separation of the flux of one cell from the flux of another cell. These degrees of freedom are known in mathematical literature as LaGrange multipliers. The introduction of the additional degrees of freedom allows for a simplification of the structure of the numerical problem, as the original variables, such as cell centered temperature (or pressure) and fluxes on the boundaries of one prismatic cell, become disconnected from the temperature and fluxes of any other cells. Thus, those unknowns can be eliminated, and that allows the reduction of the number of unknowns.

[0096] With the definitions provided above for the degrees of freedom, the MFE method can be introduced as follows: find u^ e \h,Ph <≡ Lh, and λ_h e A_h such that

for all v e V_f1, q e L_h, and μ e A_h . The finite element problem (3.1) results in the system of linear algebraic equations

with the saddle point matrix

where M = M^τ is a positive definite matrix, and D = D^τ and Σ = ∑^τ are either positive definite or positive semi-definite matrices. It can be shown that system (3.2) has the unique solution.

[0097] Iterative methods for algebraic systems with symmetric saddle point matrices are well developed. However, the efficient preconditioning technique for the saddle point matrices which arise from MFE discretizations on polyhedral grids is still a concern. Symmetric positive definite matrices are much better objects for efficient preconditioning. The system (3.2) can be transformed to the equivalent system with a symmetric positive definite matrix by using the hybridization of the mixed finite element problem (3.1). In the next paragraphs, this method is described as preferable way of solving problem (3.1).

[0098] For hybridization of mixed finite element analysis, let Ek be a grid cell in Ω/, and V^ and Lf₁^ be the restrictions of the finite element spaces \u and Lh onto Ek, respectively. Also, new finite element space A)₁ is created, which is the space of functions Xh = λ/,(x) that are defined on the interfaces Tu between grid cells as well as on the intersections of grid cells with boundary parts of T_N. On each of the interfaces a function λ_h e A)₁ equals to a constant.

[0099] To introduce the mixed hybrid finite element (MHFE) problem two additional finite element spaces and a number of the bilinear forms and linear functionals have to be defined. Two new finite element spaces are

where n is the number of cells Ek in Qh- Note that the dimension of any of the spaces V^ is at most five, and the dimension of each L^ is equal to one. [0100] For elements u,v e \h,p, q <≡ Lh, and A, μ e Ah the following bilinear forms are introduced

where

are the values of λ, μ e A_h, ti_ki is the unit normal to

directed from

and

[0101] Also, the following linear functionals are defined

where

[0102] With the above definitions, the equivalent mixed hybrid formulation of the finite element problem (3.1) reads as follows: find such that

a(u_A,v) + b(\,p_h) + φ,λ_h) = I_D(Ύ) b(μ_h,q) - d(p_h,q) = l_f(q) (3.4) c(yy_h,μ) - σ(λ_h,μ) = l_N(μ)

for all

[0103] Finite element problem (3.4) results in the system of linear algebraic equations

with the saddle point matrix

is the block diagonal matrix with the symmetric positive definite submatrices Mk, k=l,...,n, D is a diagonal positive definite or semi-definite matrix, and Σ is a diagonal positive semi- definite matrix. The components of the right-hand side subvectors g_D , f , g_N are defined by the linear functionals in (3.4).

[0104] The matrix A has a very useful representation where

is the local saddle point matrix for cell E₁ and N₁ is the corresponding assemble matrix.

[0105] Respectively, the right hand side of system (3.5) can be written as follows:

[0106] It is important to observe that matrix A₁ and subvectors g_{D l} , f _t , and g_{N 1} can be obtained by applying the local mixed finite element discretization for the following problem:

where if , k= 1 , ... ,S₁ are the faces of cell E₁.

[0107] The important properties of matrices M₁, B₁, C₁, D₁, and ∑_z in (2.6), can be demonstrated on an internal cell E₁, i.e. 3E₁ n 3Ω._h = 0 . Let cell E₁ have S₁ faces then

M₁ e R^s'^*s' is a symmetric positive definite matrix, then

where V_E is the volume of cell E₁ and . For an internal cell matrix ∑_z= 0.

[0108] Let e, = (l 1 ... l)^r e i?"' . Then 5_; = -CIe₁ . This property holds for any A₁ from (3.6).

[0109] It is pertinent to note that primary variables U₁ and P₁ , i=\,...,n, are only connected within a single cell. So, these unknowns can easily be excluded:

u = M ^ι(g_D - Cλ - Bp) (3.7)

and

p = (B^TM ¹B + DY(B^TM ^lg_D - f - B^τM ^ιCλ ). (3.8) [0110] Due to the structure of matrices M, B, and D, the matrix B¹M ^"~! B + D is diagonal, therefore, it is invertible.

[0111] Using relationships (3.7) and (3.8) system (3.5) is transformed to the system:

Sl = ξ (3.9)

where

S = C¹M ¹C - C^TM ¹CB(B^TM ¹B + Z)) ' B¹M ¹C + Σ

and

[0112] The matrix S is called "condensed matrix". This matrix is symmetric and positive definite except the case of Neumann boundary conditions when S is semi-positive definite, but has simple kernel - constant vector. This matrix is global in nature that connects all of the nodes or cells together. The large system of linear equations may be solved simultaneously.

[0113] Any iterative method can be applied to solve the system of linear equations (3.9) with that matrix. One method is Preconditioned Conjugate Gradient method (PCG), however other methods may be used. Note that in case of semi-positive defmiteness PCG should be performed in the subspace orthogonal to the kernel. After solving system (3.9), primary unknowns p and u can be recovered locally element-by-element using equations (3.8) and (3.7), respectively.

[0114] The matrix S can be also presented as where

and N₁ are the corresponding assembling matrices. The right hand side of (3.9) has a similar representation. [0115] Note that the embodiments of the invention may operate with a single primary unknown, e.g. temperature or pressure, and its associated flux. Other embodiments may operate with more than one primary unknown.

[0116] The various processes and methods outlined above may be combined in one or more different methods, used in one or more different systems, used in one or more different computer program products, according to various embodiments of the invention.

[0117] For example, one exemplary method 1000 may be to form a prismatic grid as shown in FIGURE 10. The geological and geometrical features of interest, such as pinch-out boundaries, fault lines, or well locations are projected into horizontal plane using orthogonal projection 1001. A fine rectangular conforming mesh is generated that covers all features of the projected domain of the same size as the fine grid on which the material data is provided 1002. The rectangular grid is separated in triangles 1003. The various lines and points on the grid can represent, for example, the fault lines and well locations. The triangles are coarsened in non-uniform manner 1004. It is desirable to keep fine triangulation near some geologic or geometric features, but having coarser resolution away from these features will allow for easier analysis. Such a grid is comprised only of triangles. The coarsened grid is projected vertically onto all boundary surfaces of all layers to form the prismatic grid 1005. Such a grid will contain cells, which can be triangular prisms, tetrahedrons, or pyramids. The unstructured prismatic grid built in such a way approximates boundary surfaces of all layers. Note that in convection-diffusion subsurface problems the input data is associated with millions of nodes.

[0118] Another method may be to solve a convection-diffusion problem for a geologic basin 1100 as shown in FIGURE 11. A grid is formed that models the basin, 1101. Note that the grid may be formed by the method 1000 shown in FIGURE 10 or another method may be used. The method associates one degree of freedom per cell of the grid at the cell center for primary unknown and one degree of freedom per each face of the cells at the face center for normal components of flux 1102. The grid, with associated degrees of freedom is analyzed using a mixed finite element approach 1103. This analysis produces a sparse matrix equation. The method may then solve the matrix equation to get both, the primary unknown(s) and normal components of the flux of the unknown(s) at the faces of the cells 1104. [0119] Note that any of the functions described herein may be implemented in hardware, software, and/or firmware, and/or any combination thereof. When implemented in software, the elements of the present invention are essentially the code segments to perform the necessary tasks. The program or code segments can be stored in a computer readable medium or transmitted by a computer data signal. The "computer readable medium" may include any medium that can store or transfer information. Examples of the computer readable medium include an electronic circuit, a semiconductor memory device, a ROM, a flash memory, an erasable ROM (EROM), a floppy diskette, a compact disk CD-ROM, an optical disk, a hard disk, a fiber optic medium, etc. The computer data signal may include any signal that can propagate over a transmission medium such as electronic network channels, optical fibers, air, electromagnetic, RF links, etc. The code segments may be downloaded via computer networks such as the Internet, Intranet, etc.

[0120] FIGURE 12 illustrates computer system 1200 adapted to use the present invention. Central processing unit (CPU) 1201 is coupled to system bus 1202. The CPU 1201 may be any general purpose CPU, such as an Intel Pentium processor. However, the present invention is not restricted by the architecture of CPU 1201 as long as CPU 1201 supports the inventive operations as described herein. Bus 1202 is coupled to random access memory (RAM) 1203, which may be SRAM, DRAM, or SDRAM. ROM 1204 is also coupled to bus 1202, which may be PROM, EPROM, or EEPROM. RAM 1203 and ROM 1204 hold user and system data and programs as is well known in the art.

[0121] Bus 1202 is also coupled to input/output (I/O) controller card 1205, communications adapter card 1211, user interface card 1208, and display card 1209. The I/O adapter card 1205 connects to storage devices 1206, such as one or more of a hard drive, a CD drive, a floppy disk drive, a tape drive, to the computer system. The I/O adapter 1205 is may connected to printer, which would allow the system to print paper copies of information such as document, photographs, articles, etc. Note that the printer may a printer (e.g. inkjet, laser, etc.), a fax machine, or a copier machine. Communications card 1211 is adapted to couple the computer system 1200 to a network 1212, which may be one or more of a telephone network, a local (LAN) and/or a wide-area (WAN) network, an Ethernet network, and/or the Internet network. User interface card 1208 couples user input devices, such as keyboard 1213 and pointing device 1207, to the computer system 1200. User interface card 1208 may also provides sound output to a user via speaker(s). The display card 1209 is driven by CPU 1201 to control the display on display device 1210. [0122] Although the present invention and its advantages have been described in detail, it should be understood that various changes, substitutions and alterations can be made herein without departing from the spirit and scope of the invention as defined by the appended claims. Moreover, the scope of the present application is not intended to be limited to the particular embodiments of the process, machine, manufacture, composition of matter, means, methods and steps described in the specification. As one of ordinary skill in the art will readily appreciate from the disclosure of the present invention, processes, machines, manufacture, compositions of matter, means, methods, or steps, presently existing or later to be developed that perform substantially the same function or achieve substantially the same result as the corresponding embodiments described herein may be utilized according to the present invention. Accordingly, the appended claims are intended to include within their scope such processes, machines, manufacture, compositions of matter, means, methods, or steps.

Claims

CLAIMSWhat is claimed is:

1. A method for modeling on a computer a physical region, wherein the physical region comprises a plurality of strata, the method comprising: receiving data that defines at least one physical characteristic of the physical region; providing a triangular mesh on a plane of a model of the physical region, wherein the mesh comprises a plurality of cells; coarsening the triangular mesh in a non-uniform manner such that cells that are less desirable are larger; and projecting the coarsened triangular mesh in a direction orthogonal to the plane in the physical region to form a prismatic grid, wherein each of the cells of the coarsened triangular mesh is separated into sub-cells according to the strata.

2. The method of claim 1, wherein coarsening comprises: projecting the data onto a plane; and using the data to determine which cells are less desirable.

3. The method of claim 2, wherein the model includes modeled features that model physical features in the physical region, and wherein using comprises: assigning a priority value to each cell, wherein the value is determined based on whether each cell is proximate to a modeled feature and a type of the modeled feature.

4. The method of claim 1 , wherein providing a triangular mesh comprises: providing a rectangular mesh on the plane; and splitting each cell of the rectangular mesh along at least one diagonal.

5. The method of claim 1, wherein coarsening comprises: merging two adjacent triangles by eliminating a side common to the two adjacent triangles.

6. The method of claim 1, wherein the prismatic grid comprises a plurality of prism cells, a plurality of pyramid cells, and a plurality of tetrahedron cells.

7. The method of claim 1 , wherein the method is used to model at least one flux of a physical process in the physical region, the method further comprising: assigning a plurality of degrees of freedom for the flux in each sub-cell; applying mixed finite element analysis to each of the sub-cells to produce a matrix; and solving the matrix to determine the flux of the physical process in the region.

8. The method of claim 7, wherein assigning comprises for each cell: assigning one degree of freedom for the physical process; and assigning another degree of freedom for each face of the cell.

9. The method of claim 7, wherein applying comprises: using a div-constant approach to form the finite element space.

10. The method of claim 7, wherein the physical process is a convection-diffusion process.

11. The method of claim 7, wherein the physical process is one of temperature and pressure and the physical region is a subsurface geological basin.

12. The method of claim 7, wherein the physical process involves the formation of hydrocarbon material.

13. The method of claim 7, wherein the physical process involves the movement of hydrocarbon material.

14. The method of claim 1, further comprising: deriving the data from information from a sensor that measured the at least one physical characteristic of the physical region.

15. A method for modeling a physical process and a flux of the physical process on a computer, the method comprising: forming an unstructured, prismatic grid that models a physical region, wherein the physical process operates within the physical region and the prismatic grid comprises a plurality of cells; assigning a plurality of degrees of freedom for the physical process and the flux for each cell; applying mixed finite element analysis to each of the cells to produce a matrix; and solving the matrix to determine the physical process and the flux in the region.

16. The method of claim 15, wherein forming comprises: providing a triangular mesh on a plane of a model of the physical region, wherein the mesh comprises a plurality of cells; coarsening the triangular mesh in a non-uniform manner such that cells that are less desirable are larger; and projecting the coarsened triangular mesh in a direction orthogonal to the plane in the physical region to form the prismatic grid.

17. The method of claim 15, wherein the prismatic grid comprises a plurality of prism cells, a plurality of pyramid cells, and a plurality of tetrahedron cells.

18. The method of claim 15, wherein assigning comprises: assigning one degree of freedom for the physical process for each cell; and assigning another degree of freedom for each face of the cell for each cell.

19. The method of claim 15, wherein applying comprises: using a div-constant approach to form the finite element space.

20. The method of claim 15, further comprising: using the determined physical process and flux to affect a change in the physical region.

21. The method of claim 15, wherein the physical process is one of temperature and pressure and the physical region is a subsurface geological basin.

22. A computer program product having a computer readable medium having computer program logic recorded thereon for modeling at a physical process and a flux of the physical process in a physical region, the computer program product comprising: code for forming an unstructured, prismatic grid that models the physical region; code for applying mixed finite element analysis to the prismatic grid to produce a matrix; and code for solving the matrix thereby determining the physical process and the flux in the region.

23. The computer program product of claim 22, wherein the code for forming comprises: code for providing a triangular mesh on a plane of a model of the physical region, wherein the mesh comprises a plurality of cells; coarsening the triangular mesh in a non-uniform manner such that cells that are less desirable are larger; and projecting the coarsened triangular mesh in a direction orthogonal to the plane in the physical region to form the prismatic grid.

24. The computer program product of claim 22, wherein the prismatic grid comprises a plurality of cells, and the code for applying comprises: assigning one degree of freedom for the physical process to each cell; assigning another degree of freedom for each face of the cell to each cell; and using a div-constant approach to form the finite element space.

25. The computer program product of claim 22, wherein the code for solving comprises: using preconditioned conjugate gradient analysis to solve the matrix.