/* SOApprox.h -- tools for function approximation. */ /* Last edited on 2003-09-15 17:20:27 by ra014852 */ #ifndef SOApprox_H #define SOApprox_H #include #include #include #include #include #include #include #include #include #include /* The procedures in this interface are mostly wrappers for procedures in other interfaces, such as {SOMatrix.h}, {SOFunction.h}, with some additional facilities for measuring, printing, and plotting the approximation errors. */ /* READING THINGS FROM FILES */ /* The procedures in this section read some object from the file {fileName}. The file name must include the extension; see {addext} in {js.h}. If {fileName == "-"}, these procedures will read from {stdin}. They fail with error if the file does not exist. Procedures that return pointers will return NULL if {fileName == ""}. */ typedef struct Resv_Func { char *name; /* Name of the reservoir function. */ char *descr; /* Descriptor of the reservoir function. */ int degree; /* For polynomial approximated functions, its degree. */ double *coeff; /* For polynomial approximated functions, its coefficients. */ SOFunction *F; /* Pressure, Saturation or their differentials. */ FuncMap *FMap; /* The mapping function. May be null. */ int plot; /* Binary verification for plotting: **1: 1D, *1*: 2D, 1**: 3D. */ double fRange; /* Plotting range, from -fRange to fRange. */ double fStep; /* Plotting step. */ int changrng; /* Changes automatically fRange and fStep. */ int ps_eps; /* O: Ps, 1: Eps plot file format. */ int plot_times; /* Number of images to be plotted. */ SOPlot_PSStream *ps; /* Plotting stream descriptor. */ } Resv_Func; typedef struct Resv_Well { int type; /* Water Influx (0) or Oil/Water Outflux (1). */ int ind; /* Index of the basis function that represents the well. */ double *pos; /* Point coordenates of the well. */ double tflux; /* Total well flux. */ double wflux; /* Water flux. */ double oflux; /* Oil flux. */ double qflux; /* Well flux by volume. */ double vol; /* Well volume. */ } Resv_Well; typedef struct Resv_Func_vec { nat nel; Resv_Func *el; } Resv_Func_vec; typedef struct Resv_Well_vec { nat nel; Resv_Well *el; } Resv_Well_vec; SOGrid_Tree *SOApprox_ReadTree(char *fileName); /* Reads a dyadic grid. Returns NULL if {fileName == ""}. */ SOFunction *SOApprox_ReadFunction(char *fileName); /* Reads a function. Returns NULL if {fileName == ""}. */ Basis SOApprox_ReadBasis(char *fileName); /* Reads a function basis.*/ SOMatrix SOApprox_ReadMatrix(char *fileName); /* Reads an {SOMatrix}. */ SOListMatrix SOApprox_ReadListMatrix(char *fileName); /* Reads an {SOListMatrix}. */ /* OBTAINING RIGIDITY MATRICES */ SOMatrix SOApprox_GetBasisMatrix(char *matName, bool cholesky); /* Obtains the rigidity matrix {E[i,j] == } for some function basis {B}. The matrix is read from the file {M-ev.mat}, where {M} is the given {matName}. If {cholesky} is TRUE, instead of {E} returns the Cholesky lower triangular factor {EL} of {E}, which is read from file {M-ev-ch.mat}. */ SOMatrix SOApprox_GetTransferMatrix(char *matName); /* Obtains the basis transfer matrix {F[i,j] == } for two function bases {B1,B2}. The matrix is read from the file {M-ev.mat}, where {M} is the given {matName}. */ void SOApprox_SimoGetRightHandSide ( Basis bas, /* Basis for the approximation space. */ SOListMatrix Hev, /* Evaluation dot product matrix. */ double reserv_box, Resv_Well_vec wells, SOFunction *saturation, /* Current Water saturation approximation function. */ SOFunction *pressure, /* Current Oil pressure approximation function. */ FuncMap **Frelperm, FuncMap **Fvisc, FuncMap **Fform, FuncMap **Fgravden, FuncMap *Fabperm, double_vec b ); void SOApprox_SimoGetLeftMatrix ( Basis bas, /* Basis for the approximation space. */ SOMatrix *H, /* IN/OUT Ridigity matrix. */ SOListMatrix Hev, /* Evaluation dot product matrix. */ SOListMatrix Hgr, /* Gradient dot product matrix. */ double res_side, SOFunction *saturation, /* Current Water saturation approximation function. */ SOFunction *pressure, /* Current Oil pressure approximation function. */ SOFunction *porosity, /* Porosity function. */ FuncMap *Fpor, FuncMap **Frelperm, FuncMap **Fvisc, FuncMap **Fform, FuncMap *Fabperm, int_vec *intersec ); void SOApprox_GetBasisIndeces ( Basis bas, SOGrid_Dim pDim, int numpos, int *indices, double *basispos ); /* Given a certain number of points within a SOGrid, finds the basis {bas} functions that are closer to them and return their indices. */ SOMatrix SOApprox_GetHelmholtzMatrix(char *matName, double c, bool cholesky); /* Obtains the differential operator matrix {H[i,j] = }, where {B} is some function basis, for the Helmholtz operator {D(f)(p) = SLap(f)(p) + c*f(p)} with given coefficient {c}. The matrix {H} is actually computed according to the formula {H[i,j] = + c*}. where {SGrd} is the gradient operator. The matrices {G[i,j] = } and {E[i,j] = } are read from files {M-gr.mat} and {M-ev.mat}, where {M} is the given {matName}. If {cholesky} is true, computes and returns the Cholesky lower triangular factor {HL} of {H}, instead of {H} itself. */ /* SOLVING LINEAR SYSTEMS The procedures in this section are used to solve the linear system {A * x == y}, which occurs in the approximation of {F} by {bas} (where {A[i,j] = }), or in the solution of the differential equation {D(f)(p) == F(p)} (where {A[i,j] == }. The basis elements are assumed to be scalars (range dimension = 1); the target function {F} may have any range dimension {n}. The matrix {A} (and its triangular factors, where aplicable) are assumed to be square of size {m × m}, where {m} is the number of elements in the basis. The matrices {x} and {y} have {m} rows and {n} columns. They are stored as vectors of size {n*m}, in row-by-row order. */ void SOApprox_ComputeRightHandSide ( SOFunction *f, /* Function to approximate. */ FuncMap *FMap, /* Right-hand side functional. */ Basis bas, /* Basis for the approximation space. */ SOFunction *w, /* Weight for integrals. */ SOGrid_Tree *tree, /* SOGrid to use for integration. */ double_vec y, /* IN/OUT: right-hand-side matrix, linearized. */ bool verbose /* TRUE mumbles and prints the change in {y}. */ ); /* (Re)computes the right-hand-side matrix {y[i,j] == }, where {F[j](p) = FMap(f(p), p)[j]}. In particular, if {FMap == NULL}, reduces to {y[i,j] == }. If {verbose} is true, also computes and prints the distance between the newly computed value of {y} and its input value. */ void SOApprox_GaussElim ( double *TotH, /* (WORK) temporary storage area. */ SOMatrix H, /* Left side of the system. */ double_vec b, /* Right side of the system. */ double_vec x, /* (OUT) Coordinates of new solution relative to {bas}. */ bool verbose ); void SOApprox_ListGaussElim ( double *TotH, /* (WORK) temporary storage area. */ SOListMatrix H, /* Left side of the system. */ double_vec b, /* Right side of the system. */ double_vec x, /* (OUT) Coordinates of new solution relative to {bas}. */ bool verbose ); void SOApprox_CholeskySolve ( SOMatrix L, /* Cholesky factor of system's matrix. */ int n, /* Number of columns in soln and rhs matrices. */ double_vec y, /* Right-hand-side matrix of system. */ double_vec x, /* (OUT) Coordinates of new solution relative to {bas}. */ double_vec r, /* (WORK) temporary storage area. */ double_vec s, /* (WORK) temporary storage area. */ double_vec t, /* (WORK) temporary storage area. */ bool verbose /* TRUE mumbles along the way. */ ); /* Solves the system {A x == y} on the unknown matrix {x}, for a square positive-definite symmetric matrix {A}, given the matrix {y} and the lower triangular factor {L} of the Cholesky decomposition {L * L^t} of the matrix {A}. The solution is returned in the argument {x}. The vectors {r}, {s}, and {t} must have size {m}, and are used as temporary storage areas. */ void SOApprox_GaussSeidelIteration( SOMatrix A, /* System's matrix. */ int n, /* Number of columns in soln and rhs matrices. */ double_vec y, /* Right-hand-side of system. */ double omega, /* Overrelaxation factor. */ double_vec xOld, /* (IN) Current guess at the solution. */ double_vec xNew, /* (OUT) New solution. */ double_vec r, /* (WORK) temporary storage area. */ double_vec s, /* (WORK) temporary storage area. */ bool verbose /* TRUE mumbles along the way. */ ); /* Performs one iteration of the Gauss-Seidel algorithm for solving the system {A x == y}. Upon entry, {xOld} should be the coordinates of a guess for the solution, relative to some basis. Upon exit, {xNew} will have been updated by performing one iteration of the Gauss-Seidel algorithm. The vectors {r} and {s} must have size {m}, and are used as temporary storage areas. */ void SOApprox_GaussSeidelSolve ( SOMatrix A, /* System's matrix. */ int n, /* Number of columns in soln and rhs matrices. */ double_vec y, /* Right-hand-side of system. */ double omega, /* Relaxation factor. */ nat maxIter, /* Maximum number of iterations. */ double absTol, /* Stopping criterion: absolute change in solution. */ double relTol, /* Stopping criterion: relative change in solution. */ double_vec xOld, /* (IN) Initial guess at solution. */ double_vec xNew, /* (OUT) New solution. */ double_vec r, /* (WORK) temporary storage area. */ double_vec s, /* (WORK) temporary storage area. */ bool verbose /* TRUE mumbles along the way. */ ); /* Attempts to solve the linear system {A x == y} by multiple iterations of the Gauss-Seidel method. Upon entry, {xOld} should be an initial guess for the solution. Upon exit, {xNew} will be the tentative solution for the system. The procedure will stop when the L2 difference between successive trials is at most {absTol}, or the relative difference is less than {relTol}, or after at most {maxIter} iterations. The vectors {r} and {s} must have size {m}, and are used as temporary storage areas. */ void SOApprox_ConjugateGradientSolve ( SOMatrix A, /* System's matrix. */ int n, /* Number of columns in soln and rhs matrices. */ double_vec y, /* Right-hand-side of system. */ double_vec xOld, /* (IN) Current solution. */ double_vec xNew, /* (OUT) New solution. */ double_vec r, /* (WORK) temporary storage area. */ double_vec s, /* (WORK) temporary storage area. */ bool verbose /* TRUE mumbles along the way. */ ); /* Attempts to solve the linear system {A x == y} by conjugate gradient minimization of the squared norm of the residual matrix {y - A x}. Upon entry, {xOld} should be the a guess for the solution. Upon exit, {xNew} will be a tentative solution for the system. The vectors {r} and {s} must have size {m}, and are used as temporary storage areas. */ double SOApprox_UpdateSolution ( double_vec xNew, /* (IN) New solution vector. */ double_vec x /* (I/O)Current solution vector. */ ); /* This procedure is handy when using iterative methods like GaussSeidel. It replaces the previous solution matrix {x} by the newly computed solution matrix {xNew}, comparing the two. Returns the L2 norm of the difference. */ void SOApprox_GuessSol(double_vec x); /* Fills {x} with random values. */ SOFunction *SOApprox_BuildIterFunction ( Basis bas, /* Basis for space {V[r]}. */ double_vec u, /* Coordinates of {g} relative to {bas} */ char *basName, /* Name of the (Basis bas) description file */ SOGrid_Dim pDim, /* Domain dimension of the new function */ SOGrid_Dim fDim /* Range dimension of the new function */ ); SOFunction *SOApprox_BuildFunction ( Basis bas, /* Basis for space {V[r]}. */ double_vec u, /* Coordinates of {g} relative to {basis} */ char *solName, /* Name prefix for output file. */ char *basName, /* Name of the (Basis bas) description file */ SOFunction *f /* True solution (for comparisons) */ ); /* Returns the function {g} with coefficient matrix {u} relative to the given {basis}. Also writes {g} to disk (with the given {solName} extended with "-g.fun") and compares it to the given function {f}, with {SOApprox_PrintSomeValues}. */ void SOApprox_PlotSingleFunction ( SOFunction *f, FuncMap *FMap, SOGrid_Tree *tree, char *fName, PlotOptions *plt, double *fObsMin, /* (IN/OUT) Min value seen during plot. */ double *fObsMax /* (IN/OUT) Max value seen during plot. */ ); /* Plots the function {f} alone, to a file called "{fName}.ps" or "{fName}.eps", as specified by the plotting parameters in {plt}. Uses {SO2DPlot_FunctionFigure}. The parameters {} and {} must be initialized by the client, and are updated with the maximum and minimum values of {f} found during the plot. */ void SOApprox_PrintSomeValues(SOFunction *g, SOFunction *f); /* Prints values of {f(p)}, {g(p)}, and their difference at a few points in the root cell of the grid. */ void SOApprox_PrintMaxErrorValues ( SOFunction *g, /* Approximation. */ SOFunction *f, /* Target function. */ double *gMax, /* OUT: max abs function value, per coordinate. */ double *eMax /* OUT: max abs error value, per coordinate. */ ); /* Computes and prints the maximum absolute value {gMax[j]} of the function {g}, and the maximum absolute error {eMax[j]} when compared to {f}, for each range coordinate {j}. */ void SOApprox_GetLimitValues ( SOFunction *f, /* General function. */ double *maxP, /* OUT: point with maximum value. */ double *fMax, /* OUT: max function values. */ double *fMin /* OUT: min function values. */ ); /* Simply returns on fMax and fMin the maximum and minimum values of function {f}. */ #endif