Method and device for determining hydrocarbon production for a reservoir

Abstract

The present invention relates to a method for determining hydrocarbon production for a reservoir. The method comprises determining a projector matrix based on a Jacobian matrix function of the gridded model, then splitting the Jacobian matrix into subsets of consecutive lines. For each subset of consecutive lines, creating a respective square matrix based on said subset. A determining eigenvectors and respective eigenvalues associated with the respective square matrix and then determining relevant eigenvectors having respective eigenvalues below a predetermined threshold. The projector is determined as a concatenation of the extended eigenvectors ordered according to multiple criteria: the respective order value of the subset; and the respective eigenvalue of the relevant eigenvector.

Claims

1. A method for determining hydrocarbon production for a reservoir, wherein the method comprises: modeling the reservoir with a gridded model, said model having partitions represented by a Jacobian matrix function of unknowns for computing a flow rate, the unknowns comprising: a pressure, hydrocarbon component compositions, and a temperature; permuting the Jacobian matrix function of the gridded model according to an ordering of the unknowns to obtain a first matrix the unknowns each having a respective index in said permutation, the permutation is performed so that a set of all unknowns related to a same partition in the partitions have a contiguous set of indexes in the ordering and so that inside each contiguous set of indexes the unknowns of the set of unknowns that are connected to unknowns in another set of unknowns have indexes greater than the indexes of unknowns that have no connection outside the set of unknowns; determining a projector matrix function of the first matrix, said determining of the projector matrix comprising: splitting the first matrix into subsets of consecutive lines, each subset of consecutive lines having a respective order value function of an index of a line in said subset; for each subset of consecutive lines: creating a respective square matrix based on said subset; determining extended eigenvectors and respective eigenvalues associated with the eigenvalues of the respective square matrix; determining relevant eigenvectors, the relevant eigenvector being the determined eigenvectors having respective eigenvalues below a predetermined threshold, determining the projector matrixes a concatenation of the relevant eigenvectors ordered according to multiple criteria as follows: firstly, the respective order value of the subset for which the relevant eigenvector is determined; secondly, the respective eigenvalue of the relevant eigenvector; determining a preconditioner operator based on the projector matrix, wherein the preconditioner operator is a function of Z(Z.sup.1AZ).sup.−1Z.sup.t, Z being the determined projector matrix and A being the Jacobian matrix or a matrix obtained after a transformation of the Jacobian matrix; and determining hydrocarbon production for the reservoir based on the preconditioner operator.

2. The method according to claim 1, wherein a first subset in the subsets of consecutive lines has a line corresponding to a first line of the matrix; the other subset in the subsets of consecutive lines has a line immediately subsequent to a line of another subset; a concatenation of the subsets of consecutive lines being the Jacobian matrix.

3. The method according to claim 1, wherein: the respective order value of a subset having a line subsequent to a line of a second subset is greater that the order value of the second subset.

4. The method according to claim 1, wherein: line of the respective square matrix having index equal to the index of a second line in the subset of consecutive lines is identical to said second line; line of the respective square matrix having no index equal to the index of a second line in the subset of consecutive lines is a line with 0-values.

5. A non-transitory computer readable storage medium, having stored thereon a computer program comprising program instructions, the computer program being loadable into a data-processing unit and adapted to cause the data-processing unit to carry out the steps of any of claim 1 when the computer program is run by the data-processing device.

6. A device for determining hydrocarbon production for a reservoir, wherein the device comprises an interface to receive information of the reservoir: wherein the device comprises a processor adapted to: model the reservoir with a gridded model based on the received information of the reservoir said model having partitions represented by a Jacobian matrix function of unknowns for computing a flow rate, the unknowns comprising: a pressure, a hydrocarbon component compositions, and a temperature; permuting the Jacobian matrix function of the e gridded model according to an ordering of the unknowns to obtain a first matrix, the unknowns each having a respective index in said permutation, the permutation is performed so that a set of all unknowns related to a same partition in the partitions have a contiguous set of indexes in the ordering and so that inside each contiguous set of indexes the unknowns of the set of unknowns that are connected to unknowns in another set of unknowns have indexes greater than the indexes of unknowns that have no connection outside the set of unknowns; determine a projector matrix function of the first matrix, said determining of the projector matrix comprising: splitting the Jacobian matrix into subsets of consecutive lines, each subset of consecutive lines having a respective order value function of an index of a line in said subset; for each subset of consecutive lines: creating a respective square matrix based on said subset; determining eigenvectors and respective eigen values associated with the respective square matrix; determining relevant eigenvectors, the relevant eigenvector being the determined eigenvectors having respective eigen values below a predetermined threshold, extending the relevant eigenvectors by adding zero values to the indexes corresponding to line outside the slice; determining the projector matrix as a concatenation of the extended eigenvectors ordered according to a multiple criteria as follows: firstly, the respective order value of the subset for which the relevant eigenvector is determined; secondly, the respective eigenvalue of the relevant eigenvector; determine a preconditioner operator based on the projector matrix, wherein the preconditioner operator is a function of Z(Z.sup.1AZ).sup.−1Z.sup.t, Z being the determined projector matrix and A being the Jacobian matrix or a matrix obtained after a transformation of the Jacobian matrix; and determine hydrocarbon production for the reservoir based on the preconditioner operator.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The present invention is illustrated by way of example, and not by way of limitations, in the figures of the accompanying drawings, in which like reference numerals refer to similar elements and in which:

(2) FIG. 1 is a diagram describing a prior art process of solving a flow problem;

(3) FIG. 2a is a representation of the selection of a matrix for the individual processing by a plurality of processors;

(4) FIG. 2b is a representation of the creation of a square matrix based on the previous selection;

(5) FIG. 2c is a representation of a determination of eigenvalues and eigenvectors for the created square matrix;

(6) FIG. 2d is a representation of the reduction of the number of eigenvectors based on the associated eigenvalues;

(7) FIG. 2e is a representation of a determination of the projector Z based on the determined eigenvectors for all processors;

(8) FIG. 2f is a representation of the projector Z obtained after concatenation of the extended eigenvector obtained in each subset of the global matrix. For each extended eigenvector the dash part corresponds to a eigenvector of a square matrix and the remaining corresponds to the zeros added to obtain a vector of the dimension of the global matrix A.

(9) FIG. 3 is a representation of multilevel approach for the projector Z;

(10) FIG. 4 is a representation of the determination of R.sub.i for the determination of the local solvers (part of the preconditioner);

(11) FIG. 5 is a possible embodiment for a device that enables the present invention.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

(12) To explain a possible embodiment for determining a proper projector Z in a parallelized way, FIGS. 2a to 2f provides detailed steps.

(13) For the purpose of simplification, Jacobian J.sub.X1 or a transformation of the Jacobian matrix (see above) is noted in the following description A.

(14) The purpose of determining the projector Z is to determine a proper preconditioner matrix M.sup.−1, for instance M.sup.−1=Σ.sub.Ω.sub.i R.sub.i.sup.tA.sub.i.sup.−1R.sub.i+Z(Z.sup.tAZ).sup.−1Z.sup.t and thus to ease the determination of hydrocarbon production in reservoirs for oil/gas industry.

(15) The matrix 201 is an example of A. Said matrix is a sparse matrix n×n, in which many values (blank zones in the representation of matrix 201) are null.

(16) This is mainly due to the fact that the interactions between cells are limited to the neighboring cells.

(17) In order to parallelize the task, a slice of the matrix is determined (202) for each processor p: the first processor receives the first k.sub.1 lines (e.g. the order number of said subset comprising the first k.sub.1 lines may be 1), the second processor receives the next k.sub.2 subsequent lines (e.g. the order number of said subset comprising the next k.sub.2 subsequent lines may be 2), the third processor receives the next k.sub.3 subsequent lines (e.g. the order number of said subset comprising the next k.sub.3 subsequent lines may be 3), etc.

(18) This slicing is a splitting of the Jacobian matrix (201) into subsets of consecutive lines (202).

(19) The number k.sub.1, k.sub.2, k.sub.3, etc. may be identical but it is not mandatory: it is possible to have a different number of lines if the load/capacity/performance/number of available flops/etc. are different for each processor.

(20) In FIG. 2b, each processor p then determines A.sub.p a square matrix 203 with the received lines (at the same position that in the original matrix), the square matrix may be regarded: as the matrix A where the non-received lines are set to zero. a null matrix where the received lines are copied at a position identical to their original position in the matrix A a matrix where: line having index equal to the index of a second line in the received lines is identical to said second line; line having no index equal to the index of a second line in the received lines is a line with 0-values.

(21) To detail the eigenvector calculations we need to distinguish the case where the matrix A is symmetric positive definite (SPD) or not.

(22) Let's first consider the case where the matrix A is a SPD matrix, in this case the computation of the square matrix and the eigenvalues is done as follows (we choose to denote the subset of line A.sub.p* which is the restriction of the A.sub.p matrix to the received lines as illustrated in FIG. 2c). we construct the A.sub.pp matrix which is a square matrix obtained by considering the restriction of A.sub.p* to the columns which indexes are in the p.sup.th subset of indexes; we compute the square matrix B.sub.pp which is obtained by adding all the entries A.sub.p*(i,j) that are outside the square matrix A.sub.pp in the diagonal coefficient A.sub.pp(i,i) of A.sub.pp (the indexing X(i,j) refers to a coefficient in the i.sup.th line and j.sup.th column the considered matrix X); we compute the eigenvalues λ and eigenvectors v of the general eigenvalue problem: B.sub.p.Math.v=λ.Math.A.sub.pp.Math.v such that λ<θ with θ ∈ [0, 1] (it can be proven that all eigenvalues λ are in [0, 1]).

(23) In addition, since A is a symmetric positive definite (SPD) matrix then A.sub.pp is also SPD.

(24) Thus, we can write the Cholesky decomposition of A.sub.pp=L.Math.L.sup.t where L is a lower triangular matrix (see https://en.wikipedia.org/wiki/Cholesky_decomposition) to transform the generalized eigenvalue problem into:
B.sub.ppv=λL.Math.L.sup.tv

(25) Thus it is noticeable that the eigenvectors v can be equivalently computed as v=(L.sup.1).sup.−1w with w being the eigenvector solution to the eigenvalue problem:
L.sup.−1B.sub.pp(L.sup.t).sup.−1w=λw

(26) The adaptation to non strictly SPD matrix of the projector Z construction relies upon this reformulated eigenvalue problem.

(27) When A is not a symmetric positive definite matrix (which is the case if we want to apply our method as a more parallel replacement for AMG in the CPR-AMG method where the matrix is not perfectly symmetric) we adapt of the method as follows: instead of computing the Cholesky decomposition, we compute the L.D.U decomposition of the square matrix A.sub.pp (as defined in https://en.wikipedia.org/wiki/LU_decomposition. In practice, for reservoir simulation, we assume that this decomposition is always possible without considering a permutation of A.sub.pp as in the general formula. But if such a case was encountered a simple line permutation of the A.sub.p* matrix would also allow to write the L.D.U decomposition without having to use a permutation matrix). In this decomposition L is a lower triangular matrix with unitary coefficients on the diagonal, D is a diagonal matrix and U is an upper triangular matrix with unitary coefficients on the diagonal; we consider the matrix S that is the diagonal matrix obtained from D such that the diagonal coefficient S(i,i) of S are computed as S(i,i)=sqrt(abs(D(i,i))) where S(i,i) and D(i,i) are respectively the coefficient in line i and column i of S and D, abs(x) is the absolute value of x and sqrt(x) is the square root of x. we compute the square matrix M.sub.p=S.sup.−1 L.sup.−1 B.sub.pp U.sup.−1 S.sup.−1; we build the square symmetric matrix K.sub.p=M.sub.p(M.sub.p).sup.t.-given the chosen threshold θ>0 we select all the eigenvectors w solution of K.sub.p.Math.w=λ w such that λ<θ.sup.2. finally we obtain the set of eigenvectors v needed by computing v=U.sup.−1S.sup.−1.Math.w for any w computed at the previous step.

(28) One can notice that if we apply these more general computations in the case where A is SPD, we end up with the same eigenvectors as with the computations for the previous method dedicated to the purely SPD case (but in a less optimized manner since we do at least some extra computations to obtain K.sub.p).

(29) Once the eigenvectors have been computed each vector v is prolongated into a vector x of the same dimension than the A matrix by completing with zero values for the indexes that do not correspond to a received line.

(30) In other words, in the vector x, the values that are not in zone 204 (i.e. zone 204 corresponded to the received lines) may be set to zero. When the n eigenvalues {λ.sub.k}.sub.k∈n and the n eigenvectors {x.sub.k}.sub.k∈n are determined, it is possible to order (205) the eigenvectors {x.sub.k}.sub.k∈n according to the order of the corresponding λ.sub.k (ascending order). In the example of FIG. 2d, λ.sub.1≤λ.sub.2≤λ.sub.3≤ . . . λ.sub.n.

(31) Then, all the non-identified eigenvectors, for each processor having received lines from the A matrix, may be concatenated, in the order of the sequence of lines (i.e. order value of the subset of lines) that has been received by the processor: if the processor Proc1 has received the first lines, the identified eigenvectors for said processor is retrieved and forms the first columns of a new matrix Z (in the example of FIG. 2e, three eigenvectors have been identified for Proc1, said processor Proc1 having received the first lines of the matrix A). if the processor Proc2 has received a number of immediate subsequent lines (i.e. after the lines of Proc1), the identified eigenvectors for said processor is retrieved and forms the immediate subsequent columns of the new matrix Z (in the example of FIG. 2e, two eigenvectors have been identified for Proc2, said processor Proc2 having received the immediate subsequent lines of the matrix A, after the lines of Proc1). if the processor Proc3 has received a number of immediate subsequent lines (i.e. after the lines of Proc2), the identified eigenvectors for said processor is retrieved and forms the immediate subsequent columns of the new matrix Z (in the example of FIG. 2e, four eigenvectors have been identified for Proc3, said processor Proc3 having received the immediate subsequent lines of the matrix A, after the lines of Proc2).

(32) This is equivalent to determining eigenvectors having respective eigenvalues below a predetermined threshold (θ) as relevant eigenvectors and concatenating said relevant eigenvectors according: firstly, the respective order value of the subset for which the relevant eigenvector is determined; secondly, the respective eigenvalue of the relevant eigenvector.

(33) It is noted that there are two criteria for the ordering: it means that, first, the groups of eigenvectors (i.e. one group per processor) are ordered according respective order value of the subset for which the relevant eigenvectors are determined, and then each group of eigenvectors are ordered according to the respective eigenvalue.

(34) As noted in FIG. 2c, the eigenvectors have null values in zones that do not correspond to the received lines. Therefore, and as shown in FIG. 2f, the Z matrix is formed of diagonal blocks, outside said block the value of Z is zero (see zones 207 and 208). The Z matrix has a height of n, the width of matrix Z is m, m being lesser or equal to n (depending of the value of the threshold θ: the lesser θ is, the lesser m is).

(35) Said Z matrix allows projecting the matrix A (LEV1) on a coarser grid, Z.sup.tAZ=A.sub.c (see LEV2, FIG. 3). If the same process of determining a projector is reapplied to A.sub.c, it is possible to determine a coarser level LEV3.

(36) Thanks to this approach, it is possible to have a multilevel approach: if the matrix A.sub.c is too big, the process of determination of the projector Z may be recursively applied on A.sub.c until the size of A.sub.c is correct (i.e. adequate for the use of the users computing said matrix).

(37) In order to determine the preconditioner M, the following formula is used:

(38) $M^{-1}=\underset{{\Omega }_i}{.Math.}{R}_i^t{A}_i^{-1}{R}_i+{Z\left(Z^{t}AZ\right)}^{-1}{Z}^t$
with {Ω.sub.i}.sub.i a set of local domain/space.

(39) It is noted that other formula may exist: the core of the invention is to provide coarsening method and thus the computation of Z.

(40) For each space Ω.sub.i, it is defined matrix R.sub.i. R.sub.i being a matrix has shown in FIG. 4: R.sub.i is a matrix with n column (width) and j lines (height); R.sub.i has a j×j block being an identity block 401 (1 for the diagonal, 0 for the other values), the other values of the matrix R.sub.i (except the identity block 401) being 0. The identity block 401 for R.sub.1 corresponds to columns 1 (i.e. the first) to j of the matrix R.sub.1, the identity block 401 for R.sub.2 corresponds to columns j+1 to 2j of the matrix R.sub.2, . . . , the identity block 401 for R.sub.k corresponds to columns (k−1)j+1 to kj of the matrix R.sub.2, etc.

(41) In addition, the matrix A.sub.i is defined as A.sub.i=(R.sub.iAR.sub.i.sup.t).

(42) These choices are adapted to the reservoir flow determination, digital rock physics, PNM (Pore Network model), etc. and to any fluid flow in porous media.

(43) Part of the steps described above can represent steps of an example of a computer program which may be executed by the device of FIG. 5.

(44) FIG. 5 is a possible embodiment for a device that enables the present invention.

(45) In this embodiment, the device 500 comprise a computer, this computer comprising a memory 505 to store program instructions loadable into a circuit and adapted to cause circuit 504 to carry out the steps of the present invention when the program instructions are run by the circuit 504.

(46) The memory 505 may also store data and useful information for carrying the steps of the present invention as described above.

(47) The circuit 504 may be for instance: a processor or a processing unit adapted to interpret instructions in a computer language, the processor or the processing unit may comprise, may be associated with or be attached to a memory comprising the instructions, or the association of a processor/processing unit and a memory, the processor or the processing unit adapted to interpret instructions in a computer language, the memory comprising said instructions, or an electronic card wherein the steps of the invention are described within silicon, or a programmable electronic chip such as a FPGA chip (for «Field-Programmable Gate Array»).

(48) This computer comprises an input interface 503 for the reception of data used for the above method (e.g. the reservoir model) according to the invention and an output interface 506 for providing the projector Z or the preconditioner M to an external device 507.

(49) To ease the interaction with the computer, a screen 501 and a keyboard 502 may be provided and connected to the computer circuit 504.

(50) Expressions such as “comprise”, “include”, “incorporate”, “contain”, “is” and “have” are to be construed in a non-exclusive manner when interpreting the description and its associated claims, namely construed to allow for other items or components which are not explicitly defined also to be present. Reference to the singular is also to be construed in be a reference to the plural and vice versa.

(51) A person skilled in the art will readily appreciate that various parameters disclosed in the description may be modified and that various embodiments disclosed may be combined without departing from the scope of the invention.

(52) It is noted that the previous examples are 2D examples for the sake of simplicity, but 3D examples may be derived from these 2D\examples without difficulties.