OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Update remaining infrastructure (and documentation) to account for

Browse Source 
'mex_ip_simple' now returning an indirection array/data array pair
  (akin to G.cells.facePos and G.cells.faces(:,1)) to account for
  wells.
This commit is contained in:
Bård Skaflestad 2010-08-12 17:14:57 +00:00
parent f2f971d8ed
commit d4bde033a4
9 changed files with 277 additions and 253 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

97
hybsys.c
View File

@ -20,19 +20,19 @@
/* ---------------------------------------------------------------------- */
struct hybsys *
hybsys_allocate(int max_ncf, int nc, int ncf_tot)
hybsys_allocate(int max_nconn, int nc, int nconn_tot)
/* ---------------------------------------------------------------------- */
{
struct hybsys *new;
new = malloc(1 * sizeof *new);
if (new != NULL) {
new->work = malloc(max_ncf * sizeof *new->work);
new->one = malloc(max_ncf * sizeof *new->one );
new->S = malloc(max_ncf * max_ncf * sizeof *new->S );
new->L = malloc(nc * sizeof *new->L );
new->F = malloc(ncf_tot * sizeof *new->F );
new->r = malloc(ncf_tot * sizeof *new->r );
new->work = malloc(max_nconn * sizeof *new->work);
new->one = malloc(max_nconn * sizeof *new->one );
new->S = malloc(max_nconn * max_nconn * sizeof *new->S );
new->L = malloc(nc * sizeof *new->L );
new->F = malloc(nconn_tot * sizeof *new->F );
new->r = malloc(nconn_tot * sizeof *new->r );
if ((new->work == NULL) || (new->one == NULL) || (new->S == NULL) ||
(new->L == NULL) || (new->F == NULL) || (new->r == NULL)) {
@ -66,16 +66,16 @@ hybsys_free(struct hybsys *sys)
/* ---------------------------------------------------------------------- */
void
hybsys_init(int max_ncf, int ncf_tot, struct hybsys *sys)
hybsys_init(int max_nconn, int nconn_tot, struct hybsys *sys)
/* ---------------------------------------------------------------------- */
{
int i;
for (i = 0; i < max_ncf; i++) {
for (i = 0; i < max_nconn; i++) {
sys->one[i] = 1.0;
}
for (i = 0; i < ncf_tot; i++) {
for (i = 0; i < nconn_tot; i++) {
sys->r[i] = 0.0;
}
}
@ -83,12 +83,12 @@ hybsys_init(int max_ncf, int ncf_tot, struct hybsys *sys)
/* ---------------------------------------------------------------------- */
void
hybsys_compute_components(int nc, const int *nconn,
hybsys_compute_components(int nc, const int *pconn,
const double *gflux, const double *src,
const double *Binv, struct hybsys *sys)
/* ---------------------------------------------------------------------- */
{
int c, i, p1, p2;
int c, i, p1, p2, nconn;
double csrc, a1, a2;
MAT_SIZE_T incx, incy;
@ -98,7 +98,9 @@ hybsys_compute_components(int nc, const int *nconn,
p1 = p2 = 0;
for (c = 0; c < nc; c++) {
nrows = ncols = lda = nconn[c];
p1 = pconn[c];
nconn = pconn[c + 1] - pconn[c];
nrows = ncols = lda = nconn;
/* F <- C' * inv(B) == (inv(B) * ones(n,1))' in single cell */
a1 = 1.0; a2 = 0.0;
@ -113,7 +115,7 @@ hybsys_compute_components(int nc, const int *nconn,
csrc = src[c] - ddot_(&nrows, sys->one, &incx, &gflux[p1], &incy);
/* r <- v_g */
for (i = 0; i < nconn[c]; i++) {
for (i = 0; i < nconn; i++) {
sys->r[p1 + i] = gflux[p1 + i];
}
@ -121,8 +123,7 @@ hybsys_compute_components(int nc, const int *nconn,
a1 = csrc / sys->L[c];
daxpy_(&nrows, &a1, &sys->F[p1], &incx, &sys->r[p1], &incy);
p1 += nconn[c];
p2 += nconn[c] * nconn[c];
p2 += nconn * nconn;
}
}
@ -170,31 +171,34 @@ hybsys_compute_cellmatrix(int c, int nconn, int p1, int p2,
/* ---------------------------------------------------------------------- */
void
hybsys_compute_press_flux(int nc, const int *nconn, const int *conn,
hybsys_compute_press_flux(int nc, const int *pconn, const int *conn,
const double *gflux, const double *src,
const double *Binv, const struct hybsys *sys,
const double *pi, double *press, double *flux,
double *work, const int lwork)
/* ---------------------------------------------------------------------- */
{
int c, i, p1, p2;
int c, i, nconn, p1, p2;
double a1, a2;
MAT_SIZE_T incx, incy, nrows, ncols, lda;
incx = incy = 1;
p1 = p2 = 0;
p2 = 0;
a1 = a2 = 1.0;
for (c = 0; c < nc; c++) {
assert (lwork >= nconn[c]);
p1 = pconn[c];
nconn = pconn[c + 1] - p1;
assert (lwork >= nconn);
/* Serialise interface pressures for cell */
for (i = 0; i < nconn[c]; i++) {
for (i = 0; i < nconn; i++) {
work[i] = pi[conn[p1 + i]];
}
nrows = ncols = lda = nconn[c];
nrows = ncols = lda = nconn;
/* Solve Lp = g - C'*v_g + F'*pi (for cell pressure) */
press[c] = src[c];
@ -203,7 +207,7 @@ hybsys_compute_press_flux(int nc, const int *nconn, const int *conn,
press[c] /= sys->L[c];
/* Form rhs of system B*v = B*v_g + C*p - D*pi */
for (i = 0; i < nconn[c]; i++) {
for (i = 0; i < nconn; i++) {
flux[p1 + i] = gflux[p1 + i];
work[i] = press[c] - work[i];
}
@ -213,8 +217,7 @@ hybsys_compute_press_flux(int nc, const int *nconn, const int *conn,
&a1, &Binv[p2], &lda, work, &incx,
&a2, &flux[p1], &incy);
p1 += nconn[c];
p2 += nconn[c] * nconn[c];
p2 += nconn * nconn;
}
}
@ -226,7 +229,8 @@ hybsys_compute_press_flux(int nc, const int *nconn, const int *conn,
/* ---------------------------------------------------------------------- */
static MAT_SIZE_T *
hybsys_build_ia(int nc, int nf, int *nconn, int *conn)
hybsys_build_ia(int nc, int nf,
const int *pconn, const int *conn)
/* ---------------------------------------------------------------------- */
{
MAT_SIZE_T *ia = malloc((nf+1) * sizeof *ia);
@ -243,7 +247,7 @@ hybsys_build_ia(int nc, int nf, int *nconn, int *conn)
int c, pos = 0;
for(c=0; c<nc; ++c)
{
int n = nconn[c];
int n = pconn[c + 1] - pconn[c];
for (i=pos; i<pos+n; ++i)
{
mxAssert(conn[i]<nf, "conn out of bounds");
@ -264,7 +268,8 @@ hybsys_build_ia(int nc, int nf, int *nconn, int *conn)
/* ---------------------------------------------------------------------- */
static MAT_SIZE_T*
hybsys_build_ja(int nc, int nf, int *nconn, int *conn,
hybsys_build_ja(int nc, int nf,
const int *pconn, const int *conn,
MAT_SIZE_T *ia, int *work)
/* ---------------------------------------------------------------------- */
{
@ -283,7 +288,7 @@ hybsys_build_ja(int nc, int nf, int *nconn, int *conn,
int c, pos = 0;
for(c=0; c<nc; ++c)
{
int n = nconn[c];
int n = pconn[c + 1] - pconn[c];
for (i=pos; i<pos+n; ++i)
{
int fi = conn[i];
@ -308,8 +313,10 @@ hybsys_build_ja(int nc, int nf, int *nconn, int *conn,
}
/* ---------------------------------------------------------------------- */
static void
hybsys_build_sa_and_b(int nc, int nf, int *nconn, int *conn, MAT_SIZE_T *ia,
double *S, double *R, int *work, double **sa, double **b)
hybsys_build_sa_and_b(int nc, int nf,
const int *pconn, const int *conn, MAT_SIZE_T *ia,
const double *S, const double *R,
int *work, double **sa, double **b)
/* ---------------------------------------------------------------------- */
{
*sa = malloc(ia[nf] * sizeof **sa);
@ -326,13 +333,13 @@ hybsys_build_sa_and_b(int nc, int nf, int *nconn, int *conn, MAT_SIZE_T *ia,
(*b) [i] = 0;
}
double *s = S;
double *r = R;
const double *s = S;
const double *r = R;
int c, pos = 0;
for(c=0; c<nc; ++c)
{
int n = nconn[c];
int n = pconn[c + 1] - pconn[c];
for (i=pos; i<pos+n; ++i)
{
int fi = conn[i];
@ -367,14 +374,15 @@ hybsys_build_sa_and_b(int nc, int nf, int *nconn, int *conn, MAT_SIZE_T *ia,
/* ---------------------------------------------------------------------- */
static void
hybsys_build_matrix_structure(int nc, int nf, int *nconn, int *conn,
hybsys_build_matrix_structure(int nc, int nf,
const int *pconn, const int *conn,
MAT_SIZE_T **ia, MAT_SIZE_T **ja)
/* ---------------------------------------------------------------------- */
{
int *work = malloc(nf * sizeof *work);
*ia = hybsys_build_ia(nc, nf, nconn, conn);
*ja = hybsys_build_ja(nc, nf, nconn, conn, *ia, work);
*ia = hybsys_build_ia(nc, nf, pconn, conn);
*ja = hybsys_build_ja(nc, nf, pconn, conn, *ia, work);
free(work);
}
@ -382,14 +390,15 @@ hybsys_build_matrix_structure(int nc, int nf, int *nconn, int *conn,
/* ---------------------------------------------------------------------- */
static void
hybsys_assemble_global_system(int nc, int nf, int *nconn, int *conn,
double *S, double *R,
hybsys_assemble_global_system(int nc, int nf,
const int *pconn, const int *conn,
const double *S, const double *R,
double **sa, double **b, MAT_SIZE_T *ia)
/* ---------------------------------------------------------------------- */
{
int *work = malloc(nf * sizeof *work);
hybsys_build_sa_and_b(nc, nf, nconn, conn, ia, S, R, work, sa, b);
hybsys_build_sa_and_b(nc, nf, pconn, conn, ia, S, R, work, sa, b);
free(work);
}
@ -397,11 +406,13 @@ hybsys_assemble_global_system(int nc, int nf, int *nconn, int *conn,
/* ---------------------------------------------------------------------- */
void
hybsys_assemble(int nc, int nf, int *nconn, int *conn, double *S, double *R,
hybsys_assemble(int nc, int nf,
const int *pconn, const int *conn,
const double *S, const double *R,
struct Sparse*A, double **b)
/* ---------------------------------------------------------------------- */
{
A->m = A->n = nf;
hybsys_build_matrix_structure(nc, nf, nconn, conn, &A->ia, &A->ja);
hybsys_assemble_global_system(nc, nf, nconn, conn, S, R, &A->sa, b, A->ia);
hybsys_build_matrix_structure(nc, nf, pconn, conn, &A->ia, &A->ja);
hybsys_assemble_global_system(nc, nf, pconn, conn, S, R, &A->sa, b, A->ia);
}

27
hybsys.h
View File

@ -10,34 +10,34 @@ struct hybsys {
double *F; /* C' * inv(B) */
double *r; /* system rhs per half face */
double *S; /* system matrix in single cell */
double *one; /* ones(max_ncf, 1) */
double *work; /* work array (SIZE [max_ncf, 1]) */
double *one; /* ones(max_nconn, 1) */
double *work; /* work array (SIZE [max_nconn, 1]) */
};
struct Sparse
{
int m;
int n;
MAT_SIZE_T *ia;
MAT_SIZE_T *ja;
double *sa;
int m;
int n;
MAT_SIZE_T *ia;
MAT_SIZE_T *ja;
double *sa;
};
struct hybsys *
hybsys_allocate(int max_ncf, int nc, int ncf_tot);
hybsys_allocate(int max_nconn, int nc, int nconn_tot);
void
hybsys_free(struct hybsys *sys);
void
hybsys_init(int max_ncf, int ncf_tot, struct hybsys *sys);
hybsys_init(int max_nconn, int nconn_tot, struct hybsys *sys);
void
hybsys_compute_components(int nc, const int *nconn,
hybsys_compute_components(int nc, const int *pconn,
const double *gflux, const double *src,
const double *Binv, struct hybsys *sys);
@ -50,13 +50,16 @@ hybsys_compute_cellmatrix(int c, int nconn, int p1, int p2,
const double *Binv, struct hybsys *sys);
void
hybsys_compute_press_flux(int nc, const int *nconn, const int *conn,
hybsys_compute_press_flux(int nc, const int *pconn, const int *conn,
const double *gflux, const double *src,
const double *Binv, const struct hybsys *sys,
const double *pi, double *press, double *flux,
double *work, const int lwork);
void
hybsys_assemble(int nc, int nf, int *nconn, int *conn, double *S, double *R,
hybsys_assemble(int nc, int nf,
const int *pconn, const int *conn,
const double *S, const double *R,
struct Sparse *A, double **b);
#endif /* HYBSYS_H_INCLUDED */

130
mex_compute_press_flux.c
View File

@ -30,29 +30,31 @@ verify_args(int nlhs, int nrhs, const mxArray *prhs[])
/* ---------------------------------------------------------------------- */
static void
count_cf(const mxArray *nconn, int *nc, int *max_ncf, int *ncf_tot)
count_conns(const mxArray *M_pconn, int *nc, int *max_nconn, int *nconn_tot)
/* ---------------------------------------------------------------------- */
{
int c, *pi;
double *pd;
int c, nconn, *pi;
double *pd;
*nc = mxGetNumberOfElements(nconn);
*nc = mxGetNumberOfElements(M_pconn) - 1;
*max_ncf = *ncf_tot = 0;
*max_nconn = *nconn_tot = 0;
if (mxIsDouble(nconn)) {
pd = mxGetPr(nconn);
if (mxIsDouble(M_pconn)) {
pd = mxGetPr(M_pconn);
for (c = 0; c < *nc; c++) {
*max_ncf = MAX(*max_ncf, pd[c]);
*ncf_tot += pd[c];
nconn = pd[c + 1] - pd[c];
*max_nconn = MAX(*max_nconn, nconn);
*nconn_tot += nconn;
}
} else {
pi = mxGetData(nconn);
pi = mxGetData(M_pconn);
for (c = 0; c < *nc; c++) {
*max_ncf = MAX(*max_ncf, pi[c]);
*ncf_tot += pi[c];
nconn = pi[c + 1] - pi[c];
*max_nconn = MAX(*max_nconn, nconn);
*nconn_tot += nconn;
}
}
}
@ -60,7 +62,7 @@ count_cf(const mxArray *nconn, int *nc, int *max_ncf, int *ncf_tot)
/* ---------------------------------------------------------------------- */
static void
deallocate_aux_arrays(int *nconn, int *conn,
deallocate_aux_arrays(int *pconn, int *conn,
double *src, double *gflux, double *work)
/* ---------------------------------------------------------------------- */
{
@ -71,32 +73,32 @@ deallocate_aux_arrays(int *nconn, int *conn,
if (gflux != NULL) { mxFree(gflux); }
if (src != NULL) { mxFree(src); }
if (conn != NULL) { mxFree(conn); }
if (nconn != NULL) { mxFree(nconn); }
if (pconn != NULL) { mxFree(pconn); }
}
/* ---------------------------------------------------------------------- */
static int
allocate_aux_arrays(int max_ncf, int nc, int ncf_tot,
int **nconn, int **conn,
allocate_aux_arrays(int max_nconn, int nc, int nconn_tot,
int **pconn, int **conn,
double **src, double **gflux,
double **work)
/* ---------------------------------------------------------------------- */
{
int ret, *n, *c;
int ret, *p, *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);
p = mxMalloc((nc + 1) * sizeof *p);
c = mxMalloc(nconn_tot * sizeof *c);
s = mxMalloc(nc * sizeof *s);
g = mxMalloc(nconn_tot * sizeof *g);
w = mxMalloc(max_nconn * sizeof *w);
if ((n == NULL) || (c == NULL) ||
if ((p == NULL) || (c == NULL) ||
(s == NULL) || (g == NULL) || (w == NULL)) {
deallocate_aux_arrays(n, c, s, g, w);
deallocate_aux_arrays(p, c, s, g, w);
*nconn = NULL;
*pconn = NULL;
*conn = NULL;
*src = NULL;
*gflux = NULL;
@ -104,7 +106,7 @@ allocate_aux_arrays(int max_ncf, int nc, int ncf_tot,
ret = 0;
} else {
*nconn = n;
*pconn = p;
*conn = c;
*src = s;
*gflux = g;
@ -119,54 +121,48 @@ allocate_aux_arrays(int max_ncf, int nc, int ncf_tot,
/* ---------------------------------------------------------------------- */
static void
get_nconn(const mxArray *M_nconn, int *nconn)
copy_M_int_vector(const mxArray *M_a, int *a)
/* ---------------------------------------------------------------------- */
{
size_t c, nc;
size_t nel, i;
int *pi;
double *pd;
nc = mxGetNumberOfElements(M_nconn);
nel = mxGetNumberOfElements(M_a);
if (mxIsDouble(M_nconn)) {
pd = mxGetPr(M_nconn);
if (mxIsDouble(M_a)) {
pd = mxGetPr(M_a);
for (c = 0; c < nc; c++) { nconn[c] = pd[c]; }
for (i = 0; i < nel; i++) { a[i] = pd[i] - 1; }
} else {
pi = mxGetData(M_nconn);
pi = mxGetData(M_a);
for (c = 0; c < nc; c++) { nconn[c] = pi[c]; };
for (i = 0; i < nel; i++) { a[i] = pi[i] - 1; }
}
}
/* ---------------------------------------------------------------------- */
static void
get_pconn(const mxArray *M_pconn, int *pconn)
/* ---------------------------------------------------------------------- */
{
copy_M_int_vector(M_pconn, pconn);
}
/* ---------------------------------------------------------------------- */
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; }
}
copy_M_int_vector(M_conn, conn);
}
/*
* [v, p] = mex_compute_press_flux(BI, pi, nconn, conn, F, L)
* [v, p] = mex_compute_press_flux(BI, pi, connPos, conns, F, L)
*/
/* ---------------------------------------------------------------------- */
@ -175,44 +171,44 @@ mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
/* ---------------------------------------------------------------------- */
{
int ok, i, nc, max_ncf, ncf_tot, *nconn, *conn;
int ok, i, nc, max_nconn, nconn_tot, *pconn, *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);
count_conns(prhs[2], &nc, &max_nconn, &nconn_tot);
allocate_aux_arrays(max_ncf, nc, ncf_tot,
&nconn, &conn, &src, &gflux, &work);
allocate_aux_arrays(max_nconn, nc, nconn_tot,
&pconn, &conn, &src, &gflux, &work);
plhs[0] = mxCreateDoubleMatrix(ncf_tot, 1, mxREAL);
plhs[1] = mxCreateDoubleMatrix(nc, 1, mxREAL);
plhs[0] = mxCreateDoubleMatrix(nconn_tot, 1, mxREAL);
plhs[1] = mxCreateDoubleMatrix(nc, 1, mxREAL);
sys = hybsys_allocate(max_ncf, nc, ncf_tot);
hybsys_init(max_ncf, ncf_tot, sys);
sys = hybsys_allocate(max_nconn, nc, nconn_tot);
hybsys_init(max_nconn, nconn_tot, sys);
ptr = mxGetPr(prhs[4]);
memcpy(sys->F, ptr, ncf_tot * sizeof *sys->F);
memcpy(sys->F, ptr, nconn_tot * sizeof *sys->F);
ptr = mxGetPr(prhs[5]);
memcpy(sys->L, ptr, nc * sizeof *sys->L);
memcpy(sys->L, ptr, nc * sizeof *sys->L);
get_nconn(prhs[2], nconn);
get_pconn(prhs[2], pconn);
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 */
for (i = 0; i < nc; i++) { src[i] = 0.0; } /* No sources */
for (i = 0; i < nconn_tot; i++) { gflux[i] = 0.0; } /* No gravity */
hybsys_compute_press_flux(nc, nconn, conn, gflux, src, Binv, sys,
hybsys_compute_press_flux(nc, pconn, conn, gflux, src, Binv, sys,
pi, mxGetPr(plhs[1]), mxGetPr(plhs[0]),
work, max_ncf);
work, max_nconn);
hybsys_free(sys);
deallocate_aux_arrays(nconn, conn, src, gflux, work);
deallocate_aux_arrays(pconn, conn, src, gflux, work);
}
}

38
mex_compute_press_flux.m
View File

@ -2,26 +2,31 @@ 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)
% [v, p] = mex_compute_press_flux(BI, lam, connPos, conns, F, L)
%
% PARAMETERS:
% BI - Inner product values. Typically computed using function
% 'mex_ip_simple'.
% BI - Inner product values. Typically computed using function
% 'mex_ip_simple'.
%
% lam - Interface pressure values. One scalar value for each face in
% the discretised reservoir model.
% lam - Interface pressure values. One scalar value for each face in
% the discretised reservoir model.
%
% nconn - Number of connections per cell. Typically equals the result of
% connPos - Indirection map of size [G.cells.num,1] into 'conns' table
% (i.e., the connections or DOFs). Specifically, the DOFs
% connected to cell 'i' are found in the submatrix
%
% diff(G.cells.facePos)
% conns(connPos(i) : connPos(i + 1) - 1)
%
% conn - Actual connections per cell. Typically equals the result of
% conns - A (connPos(end)-1)-by-1 array of cell connections
% (local-to-global DOF mapping in FEM parlance).
%
% G.cells.faces(:,1)
% F - Second-to-last return value from 'mex_schur_comp_symm'.
%
% F - Second-to-last return value from 'mex_schur_comp_symm'.
% L - Last return value from 'mex_schur_comp_symm'.
%
% L - Last return value from 'mex_schur_comp_symm'.
% NOTE:
% The (connPos,conns) array pair is expected to be the output of function
% 'mex_ip_simple'.
%
% RETURNS:
% v - A SUM(nconn)-by-1 array of half-contact fluxes, ordered by cells.
@ -40,21 +45,20 @@ function varargout = mex_compute_press_flux(varargin)
% rock.perm = convertFrom(rock.perm(G.cells.indexMap, :), ...
% milli*darcy);
%
% nconn = diff(G.cells.facePos);
% conn = double(G.cells.faces(:,1));
% [BI, connPos, conns] = mex_ip_simple(G, rock);
%
% BI = mex_ip_simple(G, rock);
% nconn = diff(connPos);
%
% [i, j] = blockDiagIndex(nconn, nconn);
% [S, r, F, L] = mex_schur_comp_symm(BI, nconn);
% [S, r, F, L] = mex_schur_comp_symm(BI, connPos, conns);
%
% SS = sparse(conn(i), conn(j), S);
% SS = sparse(double(conns(i)), double(conns(j)), S);
% R = accumarray(conn, r);
%
% lam = SS \ R;
%
% t0 = tic;
% [v, p] = mex_compute_press_flux(BI, lam, nconn, conn, F, L);
% [v, p] = mex_compute_press_flux(BI, lam, connPos, conns, F, L);
% toc(t0)
%
% SEE ALSO:

39
mex_ip_simple.m
View File

@ -2,24 +2,33 @@ function varargout = mex_ip_simple(varargin)
%Compute 'ip_simple' inner product values using compiled C code.
%
% SYNOPSIS:
% BI = mex_ip_simple(G, rock, nconn, conn)
% [BI, connPos, conns] = mex_ip_simple(G, rock)
% [BI, connPos, conns] = mex_ip_simple(G, rock, W)
%
% PARAMETERS:
% G - Grid data structure.
% G - Grid data structure.
%
% rock - Rock data structure. Must contain a valid field 'perm'.
% rock - Rock data structure. Must contain a valid field 'perm'.
%
% nconn - Number of connections per cell. Often coincides with
% DIFF(G.cells.facePos), but may be larger if any cells are
% perforated by one or more wells.
%
% conn - Connection data per cell. Often coincides with
% G.cells.faces(:,1) but will contain additional data if a cell is
% perforated by one or more wells.
% W - Well data structure as defined by function 'addWell'.
% OPTIONAL. Only processed if present and valid.
%
% RETURNS:
% BI - A SUM(nconn .^ 2)-by-1 array of inner product values, ordered by
% the cells of the input grid.
% BI - An array of inner product values, ordered by the cells of the
% input grid. The array contains SUM(DIFF(connPos) .^ 2)
% elements.
%
% connPos - Indirection map of size [G.cells.num,1] into 'conns' table
% (i.e., the connections or DOFs). Specifically, the DOFs
% connected to cell 'i' are found in the submatrix
%
% conns(connPos(i) : connPos(i + 1) - 1)
%
% In the absence of wells, connPos == G.cells.facePos.
%
% conns - A (connPos(end)-1)-by-1 array of cell connections
% (local-to-global DOF mapping in FEM parlance). In the
% absence of wells, conns == G.cells.faces(:,1).
%
% NOTE:
% As the return value 'BI' is but a simple data array value, it must be
@ -38,13 +47,11 @@ function varargout = mex_ip_simple(varargin)
% rock.perm = convertFrom(rock.perm(G.cells.indexMap, :), ...
% milli*darcy);
%
% nconn = diff(G.cells.facePos);
% conn = G.cells.faces(:,1);
%
% t0 = tic;
% BI = mex_ip_simple(G, rock, nconn, conn);
% [BI, connPos, conns] = mex_ip_simple(G, rock);
% toc(t0)
%
% nconn = diff(connPos);
% [i, j] = blockDiagIndex(nconn, nconn);
%
% S = struct('BI', sparse(i, j, BI), 'type', 'hybrid', 'ip', 'ip_simple')

130
mex_schur_comp_symm.c
View File

@ -26,32 +26,21 @@ verify_args(int nlhs, int nrhs, const mxArray *prhs[])
/* ---------------------------------------------------------------------- */
static void
count_cf(const mxArray *nconn, int *nc, int *max_ncf, int *ncf_tot)
static int
count_cellconn(int nc, const int *pconn)
/* ---------------------------------------------------------------------- */
{
int c, *pi;
double *pd;
int c, nconn, max_nconn;
*nc = mxGetNumberOfElements(nconn);
max_nconn = 0;
*max_ncf = *ncf_tot = 0;
for (c = 0; c < nc; c++) {
nconn = pconn[c + 1] - pconn[c];
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];
}
max_nconn = MAX(max_nconn, nconn);
}
return max_nconn;
}
@ -71,27 +60,24 @@ deallocate_aux_arrays(int *nconn, double *src, double *gflux)
/* ---------------------------------------------------------------------- */
static int
allocate_aux_arrays(int nc, int ncf_tot,
int **nconn, double **src, double **gflux)
allocate_aux_arrays(int nc, int nconn_tot,
double **src, double **gflux)
/* ---------------------------------------------------------------------- */
{
int ret, *n;
int ret;
double *s, *g;
n = mxMalloc(nc * sizeof *n);
s = mxMalloc(nc * sizeof *s);
g = mxMalloc(ncf_tot * sizeof *g);
s = mxMalloc(nc * sizeof *s);
g = mxMalloc(nconn_tot * sizeof *g);
if ((n == NULL) || (s == NULL) || (g == NULL)) {
deallocate_aux_arrays(n, s, g);
if ((s == NULL) || (g == NULL)) {
deallocate_aux_arrays(NULL, s, g);
*nconn = NULL;
*src = NULL;
*gflux = NULL;
ret = 0;
} else {
*nconn = n;
*src = s;
*gflux = g;
@ -103,26 +89,38 @@ allocate_aux_arrays(int nc, int ncf_tot,
/* ---------------------------------------------------------------------- */
static void
get_nconn(const mxArray *M_nconn, int *nconn)
static int
get_pconn(const mxArray *M_pconn, int **pconn)
/* ---------------------------------------------------------------------- */
{
size_t c, nc;
int ret;
size_t e, ne;
int *pi;
double *pd;
nc = mxGetNumberOfElements(M_nconn);
ne = mxGetNumberOfElements(M_pconn);
if (mxIsDouble(M_nconn)) {
pd = mxGetPr(M_nconn);
*pconn = mxMalloc(ne * sizeof **pconn);
for (c = 0; c < nc; c++) { nconn[c] = pd[c]; }
if (*pconn != NULL) {
if (mxIsDouble(M_pconn)) {
pd = mxGetPr(M_pconn);
for (e = 0; e < ne; e++) { (*pconn)[e] = pd[e] - 1; }
} else {
pi = mxGetData(M_pconn);
for (e = 0; e < ne; e++) { (*pconn)[e] = pi[e] - 1; };
}
ret = ne - 1;
} else {
pi = mxGetData(M_nconn);
for (c = 0; c < nc; c++) { nconn[c] = pi[c]; };
ret = -1;
}
return ret;
}
@ -175,7 +173,7 @@ get_number_of_faces(int nc, int *nconn, int *conn)
/*
* [S, r, F, L] = mex_schur_comp_symm(BI, nconn)
* [S, r, F, L] = mex_schur_comp_symm(BI, connPos, conns)
*/
/* ---------------------------------------------------------------------- */
@ -184,7 +182,8 @@ mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
/* ---------------------------------------------------------------------- */
{
int ok, nc, ncf_tot, max_ncf, sum_ncf2, p1, p2, c, i, *nconn;
int ok, nc, nconn_tot, max_nconn, sum_nconn2;
int p2, c, i, nconn, *pconn;
double *Binv, *ptr, *src, *gflux;
struct hybsys *sys;
@ -197,39 +196,42 @@ mexFunction(int nlhs, mxArray *plhs[],
ok = verify_args(nlhs, nrhs, prhs);
if (ok) {
count_cf(prhs[1], &nc, &max_ncf, &ncf_tot);
nc = get_pconn(prhs[1], &pconn);
allocate_aux_arrays(nc, ncf_tot, &nconn, &src, &gflux);
nconn_tot = pconn[nc];
max_nconn = count_cellconn(nc, pconn);
sum_ncf2 = mxGetNumberOfElements(prhs[0]);
plhs[0] = mxCreateDoubleMatrix(sum_ncf2, 1, mxREAL);
plhs[1] = mxCreateDoubleMatrix(ncf_tot, 1, mxREAL);
plhs[2] = mxCreateDoubleMatrix(ncf_tot, 1, mxREAL);
plhs[3] = mxCreateDoubleMatrix(nc, 1, mxREAL);
allocate_aux_arrays(nc, nconn_tot, &src, &gflux);
sys = hybsys_allocate(max_ncf, nc, ncf_tot);
hybsys_init(max_ncf, ncf_tot, sys);
sum_nconn2 = mxGetNumberOfElements(prhs[0]);
plhs[0] = mxCreateDoubleMatrix(sum_nconn2, 1, mxREAL);
plhs[1] = mxCreateDoubleMatrix(nconn_tot, 1, mxREAL);
plhs[2] = mxCreateDoubleMatrix(nconn_tot, 1, mxREAL);
plhs[3] = mxCreateDoubleMatrix(nc, 1, mxREAL);
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 */
sys = hybsys_allocate(max_nconn, nc, nconn_tot);
hybsys_init(max_nconn, nconn_tot, sys);
for (i = 0; i < nc; i++) { src[i] = 0.0; } /* No sources */
for (i = 0; i < nconn_tot; i++) { gflux[i] = 0.0; } /* No gravity */
#if 0
src[0] = 1;
src[nc-1]=-1;
#endif
Binv = mxGetPr(prhs[0]);
get_nconn(prhs[1], nconn);
hybsys_compute_components(nc, nconn, gflux, src, Binv, sys);
hybsys_compute_components(nc, pconn, gflux, src, Binv, sys);
ptr = mxGetPr(plhs[0]);
p1 = p2 = 0;
p2 = 0;
for (c = 0; c < nc; c++) {
hybsys_compute_cellmatrix(c, nconn[c], p1, p2, Binv, sys);
nconn = pconn[c + 1] - pconn[c];
memcpy(ptr + p2, sys->S, nconn[c] * nconn[c] * sizeof *ptr);
hybsys_compute_cellmatrix(c, nconn, pconn[c], p2, Binv, sys);
p1 += nconn[c];
p2 += nconn[c] * nconn[c];
memcpy(ptr + p2, sys->S, nconn * nconn * sizeof *ptr);
p2 += nconn * nconn;
}
#if defined(ASSEMBLE_AND_SOLVE_UMFPACK) && ASSEMBLE_AND_SOLVE_UMFPACK
@ -249,15 +251,15 @@ mexFunction(int nlhs, mxArray *plhs[],
#endif
ptr = mxGetPr(plhs[1]);
memcpy(ptr, sys->r, ncf_tot * sizeof *ptr);
memcpy(ptr, sys->r, nconn_tot * sizeof *ptr);
ptr = mxGetPr(plhs[2]);
memcpy(ptr, sys->F, ncf_tot * sizeof *ptr);
memcpy(ptr, sys->F, nconn_tot * sizeof *ptr);
ptr = mxGetPr(plhs[3]);
memcpy(ptr, sys->L, nc * sizeof *ptr);
memcpy(ptr, sys->L, nc * sizeof *ptr);
hybsys_free(sys);
deallocate_aux_arrays(nconn, src, gflux);
deallocate_aux_arrays(pconn, src, gflux);
}
}

30
mex_schur_comp_symm.m
View File

@ -2,14 +2,23 @@ function varargout = mex_schur_comp_symm(varargin)
%Compute hybrid system component matrices using compiled C code.
%
% SYNOPSIS:
% [S, r, F, L] = mex_schur_comp_symm(BI, nconn)
% [S, r, F, L] = mex_schur_comp_symm(BI, connPos, conns)
%
% PARAMETERS:
% BI - Inner product values.
% BI - Inner product values.
%
% nconn - Number of connections per cell. Often coincides with
% DIFF(G.cells.facePos), but may be larger if any cells are
% perforated by one or more wells.
% connPos - Indirection map of size [G.cells.num,1] into 'conns' table
% (i.e., the connections or DOFs). Specifically, the DOFs
% connected to cell 'i' are found in the submatrix
%
% conns(connPos(i) : connPos(i + 1) - 1)
%
% conns - A (connPos(end)-1)-by-1 array of cell connections
% (local-to-global DOF mapping in FEM parlance).
%
% NOTE:
% The (connPos,conns) array pair is expected to be the output of function
% 'mex_ip_simple'.
%
% RETURNS:
% S - A SUM(nconn .^ 2)-by-1 array of unassembled system matrix values,
@ -29,19 +38,18 @@ function varargout = mex_schur_comp_symm(varargin)
% rock.perm = convertFrom(rock.perm(G.cells.indexMap, :), ...
% milli*darcy);
%
% nconn = diff(G.cells.facePos);
% conn = G.cells.faces(:,1);
% [BI, connPos, conns] = mex_ip_simple(G, rock);
%
% BI = mex_ip_simple(G, rock, nconn, conn);
% nconn = diff(connPos);
%
% t0 = tic;
% [S, r, F, L] = mex_schur_comp_symm(BI, nconn);
% [S, r, F, L] = mex_schur_comp_symm(BI, connPos, conns);
% toc(t0)
%
% [i, j] = blockDiagIndex(nconn, nconn);
%
% SS = sparse(double(conn(i)), double(conn(j)), S);
% R = accumarray(conn, r);
% SS = sparse(double(conns(i)), double(conns(j)), S);
% R = accumarray(conns, r);
%
% lam = SS \ R;
%

25
qfs_mex.m
View File

@ -4,30 +4,25 @@ g = computeGeometry(cartGrid(cartDims, physDims));
rock.perm = ones(g.cells.num, 1);
nconn = diff(g.cells.facePos);
nconn([1,end]) = nconn([1, end]) + 1;
conn = zeros(sum(nconn), 1);
W = addWell([], g, rock, 1);
W = addWell(W , g, rock, g.cells.num);
pos = cumsum([1; double(nconn)]);
ii = mcolon(pos(1 : end-1), ...
pos(1 : end-1) - 1 + double(diff(g.cells.facePos)));
[BI, connPos, conns] = mex_ip_simple(g, rock, W);
conn(ii) = g.cells.faces(:,1);
conn([nconn(1), end]) = g.faces.num + (1:2);
nconn = diff(connPos);
BI = mex_ip_simple(g, rock, nconn, conn);
BI([nconn(1)^2, end]) = 1; % Fake production indices...
BI([nconn(1)^2, end]) = [ W.WI ];
[S, r, F, L] = mex_schur_comp_symm(BI, nconn, conn);
[S, r, F, L] = mex_schur_comp_symm(BI, connPos, conns);
[i, j] = blockDiagIndex(nconn, nconn);
SS = sparse(double(conn(i)), double(conn(j)), S);
R = accumarray(conn, r);
SS = sparse(double(conns(i)), double(conns(j)), S);
R = accumarray(conns, r);
R(end-1 : end) = [1, -1];
R(end-1 : end) = 1000 * [1, -1];
lam = SS \ R;
[flux, press] = mex_compute_press_flux(BI, lam, nconn, conn, F, L);
[flux, press] = mex_compute_press_flux(BI, lam, connPos, conns, F, L);
plotCellData(g, press);

14
test_mex_schur_comp_symm.m
View File

@ -2,22 +2,20 @@ run ../../startup
G = computeGeometry(cartGrid([200, 1], [1, 1]));
rock.perm = ones(G.cells.num, 1);
nconn = diff(G.cells.facePos);
conn = G.cells.faces(:, 1);
[BI, connPos, conns] = mex_ip_simple(G, rock);
BI = mex_ip_simple(G, rock, nconn, conn);
[S, r, F, L] = mex_schur_comp_symm(BI, nconn, conn);
[S, r, F, L] = mex_schur_comp_symm(BI, connPos, conns);
nconn = diff(connPos);
[i, j] = blockDiagIndex(nconn, nconn);
SS = sparse(double(conn(i)), double(conn(j)), S);
R = accumarray(conn, r);
SS = sparse(double(conns(i)), double(conns(j)), S);
R = accumarray(conns, r);
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);
[v, p] = mex_compute_press_flux(BI, x, connPos, conns, F, L);
plotCellData(G, p);