diff --git a/mex_compute_press_flux.c b/mex_compute_press_flux.c new file mode 100644 index 00000000..12ffd356 --- /dev/null +++ b/mex_compute_press_flux.c @@ -0,0 +1,218 @@ +#include + +#include "mex.h" +#include "matrix.h" + +#include "hybsys.h" + + +#define MAX(a,b) (((a) > (b)) ? (a) : (b)) + + +/* ---------------------------------------------------------------------- */ +static int +verify_args(int nlhs, int nrhs, const mxArray *prhs[]) +/* ---------------------------------------------------------------------- */ +{ + int ok; + + ok = (nlhs == 2) && (nrhs == 6); + ok = ok && mxIsDouble(prhs[0]); + ok = ok && mxIsDouble(prhs[1]); + ok = ok && (mxIsDouble(prhs[2]) || mxIsInt32(prhs[2])); + ok = ok && (mxIsDouble(prhs[3]) || mxIsInt32(prhs[3])); + ok = ok && mxIsDouble(prhs[4]); + ok = ok && mxIsDouble(prhs[5]); + + return ok; +} + + +/* ---------------------------------------------------------------------- */ +static void +count_cf(const mxArray *nconn, int *nc, int *max_ncf, int *ncf_tot) +/* ---------------------------------------------------------------------- */ +{ + int c, *pi; + double *pd; + + *nc = mxGetNumberOfElements(nconn); + + *max_ncf = *ncf_tot = 0; + + if (mxIsDouble(nconn)) { + pd = mxGetPr(nconn); + + for (c = 0; c < *nc; c++) { + *max_ncf = MAX(*max_ncf, pd[c]); + *ncf_tot += pd[c]; + } + } else { + pi = mxGetData(nconn); + + for (c = 0; c < *nc; c++) { + *max_ncf = MAX(*max_ncf, pi[c]); + *ncf_tot += pi[c]; + } + } +} + + +/* ---------------------------------------------------------------------- */ +static void +deallocate_aux_arrays(int *nconn, int *conn, + double *src, double *gflux, double *work) +/* ---------------------------------------------------------------------- */ +{ + /* Apparently mxFree() makes no guarantee regarding NULL arguments. + * Institute a belt-and-suspenders approach to releasing resources. + */ + if (work != NULL) { mxFree(work); } + if (gflux != NULL) { mxFree(gflux); } + if (src != NULL) { mxFree(src); } + if (conn != NULL) { mxFree(conn); } + if (nconn != NULL) { mxFree(nconn); } +} + + +/* ---------------------------------------------------------------------- */ +static int +allocate_aux_arrays(int max_ncf, int nc, int ncf_tot, + int **nconn, int **conn, + double **src, double **gflux, + double **work) +/* ---------------------------------------------------------------------- */ +{ + int ret, *n, *c; + double *s, *g, *w; + + n = mxMalloc(nc * sizeof *n); + c = mxMalloc(ncf_tot * sizeof *c); + s = mxMalloc(nc * sizeof *s); + g = mxMalloc(ncf_tot * sizeof *g); + w = mxMalloc(max_ncf * sizeof *w); + + if ((n == NULL) || (c == NULL) || + (s == NULL) || (g == NULL) || (w == NULL)) { + deallocate_aux_arrays(n, c, s, g, w); + + *nconn = NULL; + *conn = NULL; + *src = NULL; + *gflux = NULL; + *work = NULL; + + ret = 0; + } else { + *nconn = n; + *conn = c; + *src = s; + *gflux = g; + *work = w; + + ret = 1; + } + + return ret; +} + + +/* ---------------------------------------------------------------------- */ +static void +get_nconn(const mxArray *M_nconn, int *nconn) +/* ---------------------------------------------------------------------- */ +{ + size_t c, nc; + + int *pi; + double *pd; + + nc = mxGetNumberOfElements(M_nconn); + + if (mxIsDouble(M_nconn)) { + pd = mxGetPr(M_nconn); + + for (c = 0; c < nc; c++) { nconn[c] = pd[c]; } + } else { + pi = mxGetData(M_nconn); + + for (c = 0; c < nc; c++) { nconn[c] = pi[c]; }; + } +} + + +/* ---------------------------------------------------------------------- */ +static void +get_conn(const mxArray *M_conn, int *conn) +/* ---------------------------------------------------------------------- */ +{ + size_t nel, i; + + int *pi; + double *pd; + + nel = mxGetNumberOfElements(M_conn); + + if (mxIsDouble(M_conn)) { + pd = mxGetPr(M_conn); + + for (i = 0; i < nel; i++) { conn[i] = pd[i] - 1; } + } else { + pi = mxGetData(M_conn); + + for (i = 0; i < nel; i++) { conn[i] = pi[i] - 1; } + } +} + + +/* + * [v, p] = mex_compute_press_flux(BI, pi, nconn, conn, F, L) + */ + +/* ---------------------------------------------------------------------- */ +void +mexFunction(int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) +/* ---------------------------------------------------------------------- */ +{ + int ok, i, nc, max_ncf, ncf_tot, *nconn, *conn; + double *src, *gflux, *work, *Binv, *pi, *ptr; + struct hybsys *sys; + + ok = verify_args(nlhs, nrhs, prhs); + + if (ok) { + count_cf(prhs[2], &nc, &max_ncf, &ncf_tot); + + allocate_aux_arrays(max_ncf, nc, ncf_tot, + &nconn, &conn, &src, &gflux, &work); + + plhs[0] = mxCreateDoubleMatrix(ncf_tot, 1, mxREAL); + plhs[1] = mxCreateDoubleMatrix(nc, 1, mxREAL); + + sys = hybsys_allocate(max_ncf, nc, ncf_tot); + hybsys_init(max_ncf, ncf_tot, sys); + + ptr = mxGetPr(prhs[4]); + memcpy(sys->F, ptr, ncf_tot * sizeof *sys->F); + + ptr = mxGetPr(prhs[5]); + memcpy(sys->L, ptr, nc * sizeof *sys->L); + + get_nconn(prhs[2], nconn); + get_conn (prhs[3], conn); + + Binv = mxGetPr(prhs[0]); + pi = mxGetPr(prhs[1]); + + for (i = 0; i < nc; i++) { src[i] = 0.0; } /* No sources */ + for (i = 0; i < ncf_tot; i++) { gflux[i] = 0.0; } /* No gravity */ + + hybsys_compute_press_flux(nc, nconn, conn, gflux, src, Binv, sys, + pi, mxGetPr(plhs[1]), mxGetPr(plhs[0]), + work, max_ncf); + + hybsys_free(sys); + deallocate_aux_arrays(nconn, conn, src, gflux, work); + } +} diff --git a/mex_compute_press_flux.m b/mex_compute_press_flux.m new file mode 100644 index 00000000..e0546c73 --- /dev/null +++ b/mex_compute_press_flux.m @@ -0,0 +1,75 @@ +function varargout = mex_compute_press_flux(varargin) +%Derive pressure and flux from hybrid system using compiled C code. +% +% SYNOPSIS: +% [v, p] = mex_compute_press_flux(BI, lam, nconn, conn, F, L) +% +% PARAMETERS: +% BI - Inner product values. +% +% nconn - diff(G.cells.facePos) +% +% RETURNS: +% S - A SUM(nconn .^ 2)-by-1 array of unassembled system matrix values, +% ordered by cells. +% +% r - A SUM(nconn)-by-1 array of unassemble system rhs values, ordered by +% cells. +% +% F - A SUM(nconn)-by-1 array of C'*inv(B) values, ordered by cells. +% +% L - A G.cells.num-by-1 array of C'*inv(B)*C values, ordered by cells. +% +% NOTE: +% As the return value 'BI' is but a simple data array value, it must be +% subsequently assembled into the 'S.BI' sparse matrix before being used +% to solve a flow problem using, e.g., the 'solveIncompFlow' function. +% +% Moreover, the 'solveIncompFlow' function expects its 'S' parameter to +% specify a 'type' field which is consistent with the kind of matrix +% stored within 'S'. In the case of 'ip_simple', the 'type' must be the +% string value 'hybrid'. +% +% EXAMPLE: +% G = computeGeometry(processGRDECL(makeModel3([100, 60, 15]))); +% K = logNormLayers(G.cartDims, [10, 300, 40, 0.1, 100]); +% rock.perm = bsxfun(@times, [1, 100, 0.1], K(:)); +% rock.perm = convertFrom(rock.perm(G.cells.indexMap, :), ... +% milli*darcy); +% +% t0 = tic; +% BI = mex_ip_simple(G, rock); +% toc(t0) +% +% nconn = diff(G.cells.facePos); +% +% [S, r, F, L] = mex_schur_comp(BI, nconn); +% +% [i, j] = blockDiagIndex(diff(G.cells.facePos), ... +% diff(G.cells.facePos)); +% +% S = struct('BI', sparse(i, j, BI), 'type', 'hybrid', 'ip', 'ip_simple') +% +% t0 = tic; +% S2 = computeMimeticIP(G, rock) +% toc(t0) +% +% norm(S.BI - S2.BI, inf) / norm(S2.BI, inf) +% +% SEE ALSO: +% computeMimeticIP, solveIncompFlow, blockDiagIndex. + +%{ +#COPYRIGHT# +%} + +% $Date: 2010-08-04 13:44:48 +0200 (Wed, 04 Aug 2010) $ +% $Revision: 4999 $ + + buildmex -O -v -largeArrayDims -DCOMPILING_FOR_MATLAB=1 ... + mex_compute_press_flux.c hybsys.c ... + -lmwlapack -lmwblas + + % Call MEX'ed edition. + [varargout{1:nargout}] = mex_compute_press_flux(varargin{:}); +end diff --git a/test_mex_schur_comp_symm.m b/test_mex_schur_comp_symm.m index 7b44c47a..ef44ad79 100644 --- a/test_mex_schur_comp_symm.m +++ b/test_mex_schur_comp_symm.m @@ -1,15 +1,21 @@ run ../../startup -G = computeGeometry(cartGrid([2,1])); +G = computeGeometry(cartGrid([200, 1])); rock.perm = ones(G.cells.num, 1); BI = mex_ip_simple(G, rock); nconn = diff(G.cells.facePos); +conn = G.cells.faces(:, 1); [S, r, F, L] = mex_schur_comp_symm(BI, nconn); [i, j] = blockDiagIndex(nconn, nconn); -SS = sparse(double(G.cells.faces(i, 1)), ... - double(G.cells.faces(j, 1)), S ); +SS = sparse(double(conn(i)), double(conn(j)), S); R = accumarray(G.cells.faces(:, 1), r); -x = SS \ R %#ok +SS(1) = SS(1) * 2; +R([1, G.cells.num+1]) = [1, -1]; + +x = SS \ R; +[v, p] = mex_compute_press_flux(BI, x, nconn, conn, F, L); + +plotCellData(G, p);