Add MEX wrapper for computing pressure and fluxes given interface
pressures and output from 'mex_schur_comp_symm'.
This commit is contained in:
Bård Skaflestad
parent
de74b2565c
commit
98fcce29f6
218
mex_compute_press_flux.c
Normal file
218
mex_compute_press_flux.c
Normal file
@ -0,0 +1,218 @@
|
||||
#include <string.h>
|
||||
|
||||
#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);
|
||||
}
|
||||
}
|
||||
75
mex_compute_press_flux.m
Normal file
75
mex_compute_press_flux.m
Normal file
@ -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
|
||||
@ -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);
|
||||
|
||||
Loading…
Reference in New Issue
Block a user